ERStack.cxx

// -------------------------------------------------------------------------
// -----                       ERStack source file                     -----
// -----           Created 11/12/15  by V.Schetinin                    -----
// -------------------------------------------------------------------------
#include "ERStack.h"
#include "ERMCTrack.h"

#include "FairDetector.h"
#include "FairMCPoint.h"
#include "FairRootManager.h"

#include "TLorentzVector.h"
#include "TClonesArray.h"
#include "TParticle.h"
#include "TRefArray.h"

#include <iostream>
using namespace std;

#include <list>

using std::pair;


// -----   Default constructor   -------------------------------------------
ERStack::ERStack(Int_t size)
  : FairGenericStack(),
    fStack(),
    fParticles(new TClonesArray("TParticle", size)),
    fTracks(new TClonesArray("ERMCTrack", size)),
    fStoreMap(),
    fStoreIter(),
    fIndexMap(),
    fIndexIter(),
    fPointsMap(),
    fCurrentTrack(-1),
    fNPrimaries(0),
    fNParticles(0),
    fNTracks(0),
    fIndex(0),
    fStoreSecondaries(kTRUE),
    fMinPoints(1),
    fEnergyCut(0.),
    fStoreMothers(kTRUE)
{

}

// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
ERStack::~ERStack()
{
  if (fParticles) {
    fParticles->Delete();
    delete fParticles;
  }
  if (fTracks) {
    fTracks->Delete();
    delete fTracks;
  }
}
// -------------------------------------------------------------------------

void ERStack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode,
                         Double_t px, Double_t py, Double_t pz,
                         Double_t e, Double_t vx, Double_t vy, Double_t vz,
                         Double_t time, Double_t polx, Double_t poly,
                         Double_t polz, TMCProcess proc, Int_t& ntr,
                         Double_t weight, Int_t is)
{

  PushTrack( toBeDone, parentId, pdgCode,
             px,  py,  pz,
             e,  vx,  vy,  vz,
             time,  polx,  poly,
             polz, proc, ntr,
             weight, is, -1);
}


// -----   Virtual public method PushTrack   -------------------------------
void ERStack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode,
                         Double_t px, Double_t py, Double_t pz,
                         Double_t e, Double_t vx, Double_t vy, Double_t vz,
                         Double_t time, Double_t polx, Double_t poly,
                         Double_t polz, TMCProcess proc, Int_t& ntr,
                         Double_t weight, Int_t is,Int_t secondparentID)
{

  // --> Get TParticle array
  TClonesArray& partArray = *fParticles;

//cerr << "PushTrack " <<endl;
  // --> Create new TParticle and add it to the TParticle array
  Int_t trackId = fNParticles;
  Int_t nPoints = 0;
  Int_t daughter1Id = -1;
  Int_t daughter2Id = -1;
  TParticle* particle =
    new(partArray[fNParticles++]) TParticle(pdgCode, trackId, parentId,
        nPoints, daughter1Id,
        daughter2Id, px, py, pz, e,
        vx, vy, vz, time);
  particle->SetPolarisation(polx, poly, polz);
  particle->SetWeight(weight);
  particle->SetUniqueID(proc);

  // --> Increment counter
  if (parentId < 0) { fNPrimaries++; }

  // --> Set argument variable
  ntr = trackId;

  // --> Push particle on the stack if toBeDone is set

  if (toBeDone == 1) {
    fStack.push(particle);
    //cerr << "PushTrack " << particle <<endl;
  }

}
// -------------------------------------------------------------------------



// -----   Virtual method PopNextTrack   -----------------------------------
TParticle* ERStack::PopNextTrack(Int_t& iTrack)
{
  //cerr << "PopNextTrack " << iTrack << endl;
  // If end of stack: Return empty pointer
  if (fStack.empty()) {
    iTrack = -1;
    return NULL;
  }

  // If not, get next particle from stack
  TParticle* thisParticle = fStack.top();
  fStack.pop();

  if ( !thisParticle) {
    iTrack = 0;
    return NULL;
  }

  fCurrentTrack = thisParticle->GetStatusCode();
  iTrack = fCurrentTrack;
  
  return thisParticle;

}
// -------------------------------------------------------------------------



// -----   Virtual method PopPrimaryForTracking   --------------------------
TParticle* ERStack::PopPrimaryForTracking(Int_t iPrim)
{
  //cerr << "PopPrimaryForTracking " << iPrim << endl;
  // Get the iPrimth particle from the fStack TClonesArray. This
  // should be a primary (if the index is correct).

  // Test for index
  if (iPrim < 0 || iPrim >= fNPrimaries) {
    LOG(FATAL) << "Primary index out of range! " << iPrim << std::endl;
  }

  // Return the iPrim-th TParticle from the fParticle array. This should be
  // a primary.
  TParticle* part = (TParticle*)fParticles->At(iPrim);
  if ( ! (part->GetMother(0) < 0) ) {
    LOG(FATAL) << "Not a primary track! " << iPrim << std::endl;
  }

  return part;

}
// -------------------------------------------------------------------------



// -----   Virtual public method GetCurrentTrack   -------------------------
TParticle* ERStack::GetCurrentTrack() const
{
  
  TParticle* currentPart = GetParticle(fCurrentTrack);
  //cerr << "GetCurrentTrack : " <<  currentPart << endl;
  if ( ! currentPart) {
    LOG(WARNING) << "Current track not found in stack!" << std::endl;
  }
  return currentPart;
}
// -------------------------------------------------------------------------



// -----   Public method AddParticle   -------------------------------------
void ERStack::AddParticle(TParticle* oldPart)
{
  //cerr << "AddParticle " << endl;
  TClonesArray& array = *fParticles;
  TParticle* newPart = new(array[fIndex]) TParticle(*oldPart);
  newPart->SetWeight(oldPart->GetWeight());
  newPart->SetUniqueID(oldPart->GetUniqueID());
  fIndex++;
}
// -------------------------------------------------------------------------



// -----   Public method FillTrackArray   ----------------------------------
void ERStack::FillTrackArray()
{

  // std::cout << "Filling MCTrack array..." << std::endl;
  //cerr << fNParticles << endl;
  // --> Reset index map and number of output tracks
  fIndexMap.clear();
  fNTracks = 0;

  // --> Check tracks for selection criteria
  SelectTracks();

  // --> Loop over fParticles array and copy selected tracks
  for (Int_t iPart=0; iPart<fNParticles; iPart++) {
    fStoreIter = fStoreMap.find(iPart);
    if (fStoreIter == fStoreMap.end() ) {
      LOG(FATAL) << "Particle " << iPart
                 << " not found in storage map!" << std::endl;
    }
    Bool_t store = (*fStoreIter).second;
    //cerr << "store = " << store << endl;
    if (store) {
      ERMCTrack* track =
        new( (*fTracks)[fNTracks]) ERMCTrack(GetParticle(iPart));
      fIndexMap[iPart] = fNTracks;
      // --> Set the number of points in the detectors for this track
      /*for (Int_t iDet=kREF; iDet<=kPSD; iDet++) {
        pair<Int_t, Int_t> a(iPart, iDet);
        track->SetNPoints(iDet, fPointsMap[a]);
      }*/
      fNTracks++;
    } else { fIndexMap[iPart] = -2; }

  }
  //cerr << "fNTracks = " << fNTracks << endl;
  // --> Map index for primary mothers
  fIndexMap[-1] = -1;

  // --> Screen output
  Print(0);

}
// -------------------------------------------------------------------------



// -----   Public method UpdateTrackIndex   --------------------------------
void ERStack::UpdateTrackIndex(TRefArray* detList)
{

  // std::cout << "Updating track indizes..." << std::endl;
  Int_t nColl = 0;

  // First update mother ID in MCTracks
  for (Int_t i=0; i<fNTracks; i++) {
    ERMCTrack* track = (ERMCTrack*)fTracks->At(i);
    Int_t iMotherOld = track->GetMotherId();
    fIndexIter = fIndexMap.find(iMotherOld);
    if (fIndexIter == fIndexMap.end()) {
      LOG(FATAL) << "Particle index " << iMotherOld
                 << " not found in dex map! " << std::endl;
    }
    track->SetMotherId( (*fIndexIter).second );
  }

  // Now iterate through all active detectors
  TIterator* detIter = detList->MakeIterator();
  detIter->Reset();
  FairDetector* det = NULL;
  while( (det = (FairDetector*)detIter->Next() ) ) {


    // --> Get hit collections from detector
    Int_t iColl = 0;
    TClonesArray* hitArray;
    while ( (hitArray = det->GetCollection(iColl++)) ) {
      nColl++;
      Int_t nPoints = hitArray->GetEntriesFast();

      // --> Update track index for all MCPoints in the collection
      for (Int_t iPoint=0; iPoint<nPoints; iPoint++) {
        FairMCPoint* point = (FairMCPoint*)hitArray->At(iPoint);
        Int_t iTrack = point->GetTrackID();

        fIndexIter = fIndexMap.find(iTrack);
        if (fIndexIter == fIndexMap.end()) {
          LOG(FATAL) << "Particle index " << iTrack
                     << " not found in index map! " << std::endl;
        }
        point->SetTrackID((*fIndexIter).second);
        point->SetLink(FairLink("MCTrack", (*fIndexIter).second));
      }

    }   // Collections of this detector
  }     // List of active detectors

  delete detIter;
  // std::cout << "...stack and " << nColl << " collections updated." << std::endl;

}
// -------------------------------------------------------------------------



// -----   Public method Reset   -------------------------------------------
void ERStack::Reset()
{
  fIndex = 0;
  fCurrentTrack = -1;
  fNPrimaries = fNParticles = fNTracks = 0;
  while (! fStack.empty() ) { fStack.pop(); }
  fParticles->Clear();
  fTracks->Clear();
  fPointsMap.clear();
}
// -------------------------------------------------------------------------



// -----   Public method Register   ----------------------------------------
void ERStack::Register()
{
  FairRootManager::Instance()->Register("MCTrack", "Stack", fTracks,kTRUE);
}
// -------------------------------------------------------------------------



// -----   Public method Print  --------------------------------------------
void ERStack::Print(Int_t iVerbose) const
{
  // std::cout << "Number of primaries        = "
  //           << fNPrimaries << std::endl;
  // std::cout << "Total number of particles  = "
  //           << fNParticles << std::endl;
  // std::cout << "Number of tracks in output = "
  //           << fNTracks << std::endl;
  /*for (Int_t iTrack=0; iTrack<fNTracks; iTrack++) {
    ((ERMCTrack*) fTracks->At(iTrack))->Print(iTrack);
  }*/
}
// -------------------------------------------------------------------------



// -----   Public method AddPoint (for current track)   --------------------
void ERStack::AddPoint(DetectorId detId)
{
  Int_t iDet = detId;
  pair<Int_t, Int_t> a(fCurrentTrack, iDet);
  if ( fPointsMap.find(a) == fPointsMap.end() ) { fPointsMap[a] = 1; }
  else { fPointsMap[a]++; }
}
// -------------------------------------------------------------------------



// -----   Public method AddPoint (for arbitrary track)  -------------------
void ERStack::AddPoint(DetectorId detId, Int_t iTrack)
{
  if ( iTrack < 0 ) { return; }
  Int_t iDet = detId;
  pair<Int_t, Int_t> a(iTrack, iDet);
  if ( fPointsMap.find(a) == fPointsMap.end() ) { fPointsMap[a] = 1; }
  else { fPointsMap[a]++; }
}
// -------------------------------------------------------------------------




// -----   Virtual method GetCurrentParentTrackNumber   --------------------
Int_t ERStack::GetCurrentParentTrackNumber() const
{ 
  TParticle* currentPart = GetCurrentTrack();
  //cerr << "GetCurrentParentTrackNumber " <<currentPart->GetFirstMother() << endl;
  if ( currentPart ) { return currentPart->GetFirstMother(); }
  else { return -1; }
}
// -------------------------------------------------------------------------



// -----   Public method GetParticle   -------------------------------------
TParticle* ERStack::GetParticle(Int_t trackID) const
{

  if (trackID < 0 || trackID >= fNParticles) {
    LOG(FATAL) << "Particle index " << trackID
               << " out of range." << std::endl;
  }
  //cerr << "GetParticle: " << trackID << endl;
  return (TParticle*)fParticles->At(trackID);
}
// -------------------------------------------------------------------------



// -----   Private method SelectTracks   -----------------------------------
void ERStack::SelectTracks()
{

  // cerr << "SelectTracks" << endl;
  // --> Clear storage map
  fStoreMap.clear();

  // --> Check particles in the fParticle array
  for (Int_t i=0; i<fNParticles; i++) {

    TParticle* thisPart = GetParticle(i);
    Bool_t store = kTRUE;

    // --> Get track parameters
    Int_t iMother   = thisPart->GetMother(0);
    TLorentzVector p;
    thisPart->Momentum(p);
    Double_t energy = p.E();
    Double_t mass   = p.M();
//    Double_t mass   = thisPart->GetMass();
    Double_t eKin = energy - mass;
    if(eKin < 0.0) { eKin=0.0; }
    // --> Calculate number of points
    Int_t nPoints = 0;
    /*
    for (Int_t iDet=kMVD; iDet<=kPSD; iDet++) {
      pair<Int_t, Int_t> a(i, iDet);
      if ( fPointsMap.find(a) != fPointsMap.end() ) {
        nPoints += fPointsMap[a];
      }
    }
    */
    // --> Check for cuts (store primaries in any case)
    if (iMother < 0) { store = kTRUE; }
    else {
      //if (!fStoreSecondaries) { store = kFALSE; }
      if (nPoints < fMinPoints) { store = kFALSE; }
      if (eKin < fEnergyCut) { store = kFALSE; }
    }
    // --> Set storage flag
    fStoreMap[i] = store;
  }

  // --> If flag is set, flag recursively mothers of selected tracks
  if (fStoreMothers) {
    for (Int_t i=0; i<fNParticles; i++) {
      if (fStoreMap[i]) {
        Int_t iMother = GetParticle(i)->GetMother(0);
        while(iMother >= 0) {
          fStoreMap[iMother] = kTRUE;
          iMother = GetParticle(iMother)->GetMother(0);
        }
      }
    }
  }

}
// -------------------------------------------------------------------------



ClassImp(ERStack)