// ------------------------------------------------------------------------- // ----- FairGenericStack source file ----- // ----- Created 10/08/04 by D. Bertini ----- // ------------------------------------------------------------------------- #include "FairGenericStack.h" #include "FairLogger.h" #include "TLorentzVector.h" #include "TParticle.h" // ----- Default constructor ------------------------------------------- FairGenericStack::FairGenericStack() : TVirtualMCStack(), fStack(), fParticles(new TClonesArray("TParticle", 100)), fTracks(NULL), fStoreMap(), fStoreIter(), fIndexMap(), fIndexIter(), fPointsMap(), fCurrentTrack(-1), fNPrimaries(0), fNParticles(0), fNTracks(0), fIndex(0), fStoreSecondaries(kTRUE), fMinPoints(1), fEnergyCut(0.), fStoreMothers(kTRUE) { } // ------------------------------------------------------------------------- // ----- Constructor with estimated array dimension -------------------- FairGenericStack::FairGenericStack(Int_t size) : TVirtualMCStack(), fStack(), fParticles(new TClonesArray("TParticle", size)), fTracks(NULL), fStoreMap(), fStoreIter(), fIndexMap(), fIndexIter(), fPointsMap(), fCurrentTrack(-1), fNPrimaries(0), fNParticles(0), fNTracks(0), fIndex(0), fStoreSecondaries(kTRUE), fMinPoints(1), fEnergyCut(0.), fStoreMothers(kTRUE) { } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- FairGenericStack::~FairGenericStack() { if (fParticles) { fParticles->Delete(); delete fParticles; } if (fTracks) { fTracks->Delete(); delete fTracks; } } // ------------------------------------------------------------------------- // ----- Virtual method PushTrack -------------------------------------- void FairGenericStack::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); } void FairGenericStack::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 secondMotherID) { // --> Get TParticle array TClonesArray& partArray = *fParticles; // --> 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); } } // ----- Virtual method PopNextTrack ----------------------------------- TParticle* FairGenericStack::PopNextTrack(Int_t& iTrack) { // 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* FairGenericStack::PopPrimaryForTracking(Int_t iPrim) { // 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 << FairLogger::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 << FairLogger::endl; } return part; } // ------------------------------------------------------------------------- // ----- Public method AddParticle ------------------------------------- void FairGenericStack::AddParticle(TParticle* oldPart) { TClonesArray& array = *fParticles; TParticle* newPart = new(array[fIndex]) TParticle(*oldPart); newPart->SetWeight(oldPart->GetWeight()); newPart->SetUniqueID(oldPart->GetUniqueID()); fIndex++; } TParticle* FairGenericStack::GetParticle(Int_t trackID) const { if (trackID < 0 || trackID >= fNParticles) { LOG(FATAL) << "Particle index " << trackID << " out of range." << FairLogger::endl; } return (TParticle*)fParticles->At(trackID); } // ------------------------------------------------------------------------- // ----- Public method FillTrackArray ---------------------------------- void FairGenericStack::FillTrackArray() { } // ------------------------------------------------------------------------- // ----- Public method UpdateTrackIndex -------------------------------- void FairGenericStack::UpdateTrackIndex(TRefArray* detList) { } // ------------------------------------------------------------------------- // ----- Public method Reset ------------------------------------------- void FairGenericStack::Reset() { fIndex = 0; fCurrentTrack = -1; fNPrimaries = fNParticles = fNTracks = 0; while (! fStack.empty() ) { fStack.pop(); } fParticles->Clear(); // fTracks->Clear(); fPointsMap.clear(); } // ------------------------------------------------------------------------- // ----- Public method Register ---------------------------------------- void FairGenericStack::Register() { } // ------------------------------------------------------------------------- // ----- Public method Print -------------------------------------------- void FairGenericStack::Print(Int_t iVerbose) const { } // ------------------------------------------------------------------------- // ----- Virtual method SetCurrentTrack -------------------------------- void FairGenericStack::SetCurrentTrack(Int_t iTrack) { fCurrentTrack = iTrack; } // ------------------------------------------------------------------------- // ----- Virtual method GetNtrack -------------------------------------- Int_t FairGenericStack::GetNtrack() const { return fNParticles; } // ------------------------------------------------------------------------- // ----- Virtual method GetNprimary ------------------------------------ Int_t FairGenericStack::GetNprimary() const { return fNPrimaries; } // ------------------------------------------------------------------------- // ----- Virtual method GetCurrentTrack -------------------------------- TParticle* FairGenericStack::GetCurrentTrack() const { TParticle* currentPart = GetParticle(fCurrentTrack); if ( ! currentPart) { LOG(WARNING) << "Current track not found in stack!" << FairLogger::endl; } return currentPart; } // ------------------------------------------------------------------------- // ----- Virtual method GetCurrentTrackNumber -------------------------- Int_t FairGenericStack::GetCurrentTrackNumber() const { return fCurrentTrack; } // ------------------------------------------------------------------------- // ----- Virtual method GetCurrentParentTrackNumber -------------------- Int_t FairGenericStack::GetCurrentParentTrackNumber() const { TParticle* currentPart = GetCurrentTrack(); if ( currentPart ) { return currentPart->GetFirstMother(); } else { return -1; } } void FairGenericStack::SelectTracks(Int_t begin, Int_t end) { // --> Clear storage map fStoreMap.clear(); // --> Check particles in the fParticle array for (Int_t i=0; i 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; // --> Calculate number of points Int_t nPoints = 0; for (Int_t iDet=begin; iDet 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; iGetMother(0); while(iMother >= 0) { fStoreMap[iMother] = kTRUE; iMother = GetParticle(iMother)->GetMother(0); } } } } } // ------------------------------------------------------------------------- ClassImp(FairGenericStack)