123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691 |
- // -------------------------------------------------------------------------
- // ----- FairStack source file -----
- // ----- Created 10/08/04 by D. Bertini / V. Friese -----
- // ----- adopted for NICA/MPD 29/03/10 (litvin) -----
- // ----- adopted for NICA/MPD 20/12/19 (ABychkov) -----
- // ----- added external decayer 30/04/20 (AZinchenko) -----
- // -------------------------------------------------------------------------
- #include "MpdStack.h"
- #include "MpdDecayer.h"
- #include "MpdDecayerPyt8.h"
- #include "MpdMCTrack.h"
- #include "FairDetector.h"
- #include "FairMCPoint.h"
- #include "FairRootManager.h"
- #include "TError.h"
- #include "TLorentzVector.h"
- #include "TParticle.h"
- #include "TRefArray.h"
- #include "TParticlePDG.h"
- #include "TVirtualMC.h"
- #include <list>
- #include <iostream>
- using std::cout;
- using std::endl;
- using std::pair;
- using std::vector;
- // ----- Default constructor -------------------------------------------
- MpdStack::MpdStack(Int_t size) // {
- // fParticles = new TClonesArray("TParticle", size);
- // fTracks = new TClonesArray("MpdMCTrack", size);
- // fCurrentTrack = -1;
- // fNPrimaries = fNParticles = fNTracks = 0;
- // fIndex = 0;
- // fStoreSecondaries = kTRUE;
- // fMinPoints = 1;
- // fEnergyCut = 0.;
- // fStoreMothers = kTRUE;
- //}
- : FairGenericStack(),
- fStack(),
- fParticles(new TClonesArray("TParticle", size)),
- fTracks(new TClonesArray("MpdMCTrack", size)),
- fStoreMap(),
- fStoreIter(),
- fIndexMap(),
- fIndexIter(),
- fPointsMap(),
- fCurrentTrack(-1),
- fNPrimaries(0),
- fNParticles(0),
- fNTracks(0),
- fIndex(0),
- fDecayFlag(0),
- fStoreSecondaries(kTRUE),
- fMinPoints(1),
- //fEnergyCut(0.),
- fEnergyCut(-0.001), //AZ
- fStoreMothers(kTRUE),
- fMinMotherMass(0.4),
- fMaxMotherMass(6.1),
- fRadiusCut(172), //175 NG ecal outer XYradius, some mm smaller
- fVzCut(311), //296 NG barrel half length, TODO does not include EndCaps
- fNoZDC(kTRUE) //NG
- {
- fDecays = new TClonesArray("TParticle");
- }
- // -------------------------------------------------------------------------
- // ----- Destructor ----------------------------------------------------
- MpdStack::~MpdStack() {
- if (fParticles) {
- fParticles->Delete();
- delete fParticles;
- }
- if (fTracks) {
- fTracks->Delete();
- delete fTracks;
- }
- if (fDecays) {
- fDecays->Delete();
- delete fDecays;
- }
- }
- // -------------------------------------------------------------------------
-
- // ----- Virtual public method PushTrack -------------------------------
- void MpdStack::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) {
- // --> 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);
- //AZ - global polarization transfer to decay Lambda
- if (TMath::Abs(pdgCode) == 3122 && proc == kPDecay) {
- TParticle *moth = (TParticle*) partArray.UncheckedAt(parentId);
- TVector3 polar;
- moth->GetPolarisation(polar);
- Double_t value = moth->GetWeight();
- if (TMath::Abs(moth->GetPdgCode()) == 3312) value *= 0.927; // Xi+-
- else if (TMath::Abs(moth->GetPdgCode()) == 3322) value *= 0.900; // Xi0
- else if (TMath::Abs(moth->GetPdgCode()) == 3212) {
- // Sigma0
- value *= (1./3);
- polar *= -1.0;
- }
- particle->SetPolarisation(polar);
- particle->SetWeight(value);
- }
-
- // --> Increment counter
- if (parentId < 0) fNPrimaries++;
- // --> Set argument variable
- ntr = trackId;
-
- //AZ
- if (fDecayFlag && (!particle->GetPDG()->Stable() || (particle->GetPdgCode() == 221) || (particle->GetPdgCode() == 111)) && particle->GetPDG()->DecayList()) {
- // Decay unstable particle
- //MpdDecayer::Instance()->Decay(particle);
- //Int_t npart = MpdDecayer::Instance()->ImportParticles(fDecays);
- MpdDecayerPyt8::Instance()->Decay(particle);
- Int_t npart = MpdDecayerPyt8::Instance()->ImportParticles(fDecays);
- // Copy decay products to vector to avoid TClonesArray overwriting
- // for recursive decays
- vector<TParticle> vDecays;
-
- for (Int_t j = 0; j < npart; ++j) {
- TParticle tpart = *((TParticle*) fDecays->UncheckedAt(j));
- tpart.Print();
- vDecays.push_back(tpart);
- }
- for (Int_t j = 0; j < npart; ++j) { // skip mother particle
- if (j == 0) {
- particle->SetStatusCode(11); // decayed particle
- continue; // skip mother particle
- }
- //TParticle *part = (TParticle*) fDecays->UncheckedAt(j);
- TParticle *part = &vDecays[j];
- Int_t toDo = 1;
- if (TString(gMC->GetName()) == "TGeant4") toDo = 0;
- //PushTrack(toBeDone, trackId, part->GetPdgCode(),
- PushTrack(toDo, trackId, part->GetPdgCode(),
- part->Px(), part->Py(), part->Pz(),
- part->Energy(), part->Vx(), part->Vy(), part->Vz(),
- time, polx, poly, polz, proc, ntr,
- weight, is);
- }
- if (parentId < 0) fNPrimaries += (npart-1); // treat decay products as primaries (dirty trick)
- } else {
- // --> Push particle on the stack if toBeDone is set
- if (toBeDone == 1) fStack.push(particle);
- }
- // --> Push particle on the stack if toBeDone is set
- //AZ if (toBeDone == 1) fStack.push(particle);
- }
- void MpdStack::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,
- new(partArray[fNParticles++]) TParticle(pdgCode, trackId, secondMotherID,
- nPoints, daughter1Id,
- daughter2Id, px, py, pz, e,
- vx, vy, vz, time);
- particle->SetPolarisation(polx, poly, polz);
- particle->SetWeight(weight);
- particle->SetUniqueID(proc);
- //AZ - global polarization transfer to decay Lambda
- if (TMath::Abs(pdgCode) == 3122 && proc == kPDecay) {
- TParticle *moth = (TParticle*) partArray.UncheckedAt(parentId);
- TVector3 polar;
- moth->GetPolarisation(polar);
- Double_t value = moth->GetWeight();
- if (TMath::Abs(moth->GetPdgCode()) == 3312) value *= 0.927; // Xi+-
- else if (TMath::Abs(moth->GetPdgCode()) == 3322) value *= 0.900; // Xi0
- else if (TMath::Abs(moth->GetPdgCode()) == 3212) {
- // Sigma0
- value *= (1./3);
- polar *= -1.0;
- }
- particle->SetPolarisation(polar);
- particle->SetWeight(value);
- }
- // --> Increment counter
- if (parentId < 0) fNPrimaries++;
- // --> Set argument variable
- ntr = trackId;
- //AZ
- if (fDecayFlag && (!particle->GetPDG()->Stable() || (particle->GetPdgCode() == 221) || (particle->GetPdgCode() == 111)) && particle->GetPDG()->DecayList()) {
- // Decay unstable particle
- //MpdDecayer::Instance()->Decay(particle);
- //Int_t npart = MpdDecayer::Instance()->ImportParticles(fDecays);
- MpdDecayerPyt8::Instance()->Decay(particle);
- Int_t npart = MpdDecayerPyt8::Instance()->ImportParticles(fDecays);
- // Copy decay products to vector to avoid TClonesArray overwriting
- // for recursive decays
- vector<TParticle> vDecays;
- for (Int_t j = 0; j < npart; ++j) {
- TParticle tpart = *((TParticle*) fDecays->UncheckedAt(j));
- tpart.Print();
- vDecays.push_back(tpart);
- }
-
- for (Int_t j = 0; j < npart; ++j) {
- if (j == 0) {
- particle->SetStatusCode(11); // decayed particle
- continue; // skip mother particle
- }
- //TParticle *part = (TParticle*) fDecays->UncheckedAt(j);
- TParticle *part = &vDecays[j];
- Int_t toDo = 1;
- if (TString(gMC->GetName()) == "TGeant4") toDo = 0;
- //PushTrack(toBeDone, trackId, part->GetPdgCode(),
- PushTrack(toDo, trackId, part->GetPdgCode(),
- part->Px(), part->Py(), part->Pz(),
- part->Energy(), part->Vx(), part->Vy(), part->Vz(),
- time, polx, poly, polz, proc, ntr,
- weight, is);
- }
- if (parentId < 0) fNPrimaries += (npart-1); // treat decay products as primaries (dirty trick)
- } else {
- // --> Push particle on the stack if toBeDone is set
- if (toBeDone == 1) fStack.push(particle);
- }
- // --> Push particle on the stack if toBeDone is set
- //AZ if (toBeDone == 1) fStack.push(particle);
- }
- // -------------------------------------------------------------------------
-
- // ----- Virtual method PopNextTrack -----------------------------------
- TParticle* MpdStack::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* MpdStack::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
- //AZ-15.05.20 if (iPrim < 0 || iPrim >= fNPrimaries) {
- if (iPrim < 0 || iPrim >= fParticles->GetEntriesFast()) {
- cout << fParticles->GetEntriesFast() << endl;
- cout << "-E- MpdStack: Primary index out of range! " << iPrim << endl;
- Fatal("MpdStack::PopPrimaryForTracking", "Index out of range");
- }
- // Return the iPrim-th TParticle from the fParticle array. This should be
- // a primary.
- TParticle* part = (TParticle*)fParticles->At(iPrim);
- /*AZ-30.04.20
- if ( ! (part->GetMother(0) < 0) ) {
- cout << "-E- MpdStack:: Not a primary track! " << iPrim << endl;
- Fatal("MpdStack::PopPrimaryForTracking", "Not a primary track");
- }
- */
- if (part->GetStatusCode() == 11) return (TParticle*)NULL; //AZ-30.04.20
- return part;
- }
- // -------------------------------------------------------------------------
- // ----- Virtual public method GetCurrentTrack -------------------------
- TParticle* MpdStack::GetCurrentTrack() const {
- TParticle* currentPart = GetParticle(fCurrentTrack);
- if ( ! currentPart) {
- cout << "-W- MpdStack: Current track not found in stack!" << endl;
- Warning("MpdStack::GetCurrentTrack", "Track not found in stack");
- }
- return currentPart;
- }
- // -------------------------------------------------------------------------
-
- // ----- Public method AddParticle -------------------------------------
- void MpdStack::AddParticle(TParticle* oldPart) {
- TClonesArray& array = *fParticles;
- TParticle* newPart = new(array[fIndex]) TParticle(*oldPart);
- newPart->SetWeight(oldPart->GetWeight());
- newPart->SetUniqueID(oldPart->GetUniqueID());
- fIndex++;
- }
- // -------------------------------------------------------------------------
- // ----- Public method FillTrackArray ----------------------------------
- void MpdStack::FillTrackArray() {
- cout << "-I- MpdStack: Filling MCTrack array..." << 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() ) {
- cout << "-E- MpdStack: Particle " << iPart
- << " not found in storage map!" << endl;
- Fatal("MpdStack::FillTrackArray",
- "Particle not found in storage map.");
- }
- Bool_t store = (*fStoreIter).second;
- if (store) {
- MpdMCTrack* track =
- new( (*fTracks)[fNTracks]) MpdMCTrack(GetParticle(iPart));
- fIndexMap[iPart] = fNTracks;
- // --> Set the number of points in the detectors for this track
- for (Int_t iDet=kSTS; iDet<=kMCORD; iDet++) {
- pair<Int_t, Int_t> a(iPart, iDet);
- track->SetNPoints(iDet, fPointsMap[a]);
- }
- fNTracks++;
- }
- else fIndexMap[iPart] = -2;
- }
- // --> Map index for primary mothers
- fIndexMap[-1] = -1;
- // --> Screen output
- Print(0);
- }
- // -------------------------------------------------------------------------
- // ----- Public method UpdateTrackIndex --------------------------------
- void MpdStack::UpdateTrackIndex(TRefArray* detList) {
- cout << "-I- MpdStack: Updating track indizes...";
- Int_t nColl = 0;
- // First update mother ID in MCTracks
- for (Int_t i=0; i<fNTracks; i++) {
- MpdMCTrack* track = (MpdMCTrack*)fTracks->At(i);
- Int_t iMotherOld = track->GetMotherId();
- if (iMotherOld < 0) continue; //AZ
- fIndexIter = fIndexMap.find(iMotherOld);
- if (fIndexIter == fIndexMap.end()) {
- cout << "-E- MpdStack: Particle index " << iMotherOld
- << " not found in dex map! " << endl;
- Fatal("MpdStack::UpdateTrackIndex",
- "Particle index not found in map");
- }
- 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()) {
- cout << "-E- MpdStack: Particle index " << iTrack
- << " not found in index map! " << endl;
- Fatal("MpdStack::UpdateTrackIndex",
- "Particle index not found in map");
- }
- point->SetTrackID((*fIndexIter).second);
- //point->SetLink(FairLink("MCTrack", (*fIndexIter).second));
- }
- } // Collections of this detector
- } // List of active detectors
- cout << "...stack and " << nColl << " collections updated." << endl;
- }
- // -------------------------------------------------------------------------
- // ----- Public method Reset -------------------------------------------
- void MpdStack::Reset() {
- fIndex = 0;
- fCurrentTrack = -1;
- fNPrimaries = fNParticles = fNTracks = 0;
- while (! fStack.empty() ) fStack.pop();
- fParticles->Clear();
- fTracks->Clear();
- fPointsMap.clear();
- }
- // -------------------------------------------------------------------------
- // ----- Public method Register ----------------------------------------
- void MpdStack::Register() {
- FairRootManager::Instance()->Register("MCTrack", "Stack", fTracks,kTRUE);
- }
- // -------------------------------------------------------------------------
- // ----- Public method Print --------------------------------------------
- void MpdStack::Print(Int_t iVerbose) const {
- cout << "-I- MpdStack: Number of primaries = "
- << fNPrimaries << endl;
- cout << " Total number of particles = "
- << fNParticles << endl;
- cout << " Number of tracks in output = "
- << fNTracks << endl;
- if (iVerbose) {
- for (Int_t iTrack=0; iTrack<fNTracks; iTrack++)
- ((MpdMCTrack*) fTracks->At(iTrack))->Print(iTrack);
- }
- }
- // -------------------------------------------------------------------------
- // ----- Public method AddPoint (for current track) --------------------
- void MpdStack::AddPoint(DetectorIdMPD 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 MpdStack::AddPoint(DetectorIdMPD 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 MpdStack::GetCurrentParentTrackNumber() const {
- TParticle* currentPart = GetCurrentTrack();
- if ( currentPart ) return currentPart->GetFirstMother();
- else return -1;
- }
- // -------------------------------------------------------------------------
- // ----- Public method GetParticle -------------------------------------
- TParticle* MpdStack::GetParticle(Int_t trackID) const {
- if (trackID < 0 || trackID >= fNParticles) {
- cout << "-E- MpdStack: Particle index " << trackID
- << " out of range." << endl;
- Fatal("MpdStack::GetParticle", "Index out of range");
- }
- return (TParticle*)fParticles->At(trackID);
- }
- // -------------------------------------------------------------------------
- // ----- Private method SelectTracks -----------------------------------
- void MpdStack::SelectTracks() {
- //#define EDEBUG
- #ifdef EDEBUG
- static Int_t lEDEBUGcounter=0;
- if (lEDEBUGcounter<=20)
- std::cout << "EDEBUG-- MpdStack::SelectTracks(): entered " << fNParticles << endl;
- lEDEBUGcounter++;
- #endif
- // --> 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;
-
- #ifdef EDEBUG
- if (lEDEBUGcounter<=20)
- thisPart->Print();
- lEDEBUGcounter++;
- #endif
- #undef EDEBUG
- // --> Get track parameters
- Int_t iMother = thisPart->GetMother(0);
- TLorentzVector p;
- thisPart->Momentum(p);
- Double_t energy = p.E();
- Double_t mass = thisPart->GetMass();
- Double_t eKin = energy - mass;
- // --> Calculate number of points
- Int_t nPoints = 0, nPointsZDC = 0, nPointsNoZDC = 0;
- for (Int_t iDet=kSTS; iDet<=kMCORD; iDet++) {
- pair<Int_t, Int_t> a(i, iDet);
- if ( fPointsMap.find(a) != fPointsMap.end() ) {
- nPoints += fPointsMap[a];
- if (iDet == kZDC) nPointsZDC += fPointsMap[a]; //NG save points in ZDC
- else nPointsNoZDC += fPointsMap[a]; // points not in ZDC
- }
- }
-
- // --> 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;
- if (fNoZDC && nPointsZDC && !nPointsNoZDC) store = kFALSE; //NG cut many secondary zdc tracks
- // !!!!!AZ - store products of potentially interesting decays
- if (nPoints == 0) {
- store = kFALSE;
- if (fStoreMothers) {
- TParticle* mother = GetParticle(iMother);
- if (mass < 0.1 && thisPart->P() < 0.01) store = kFALSE; // P < cut on e+- and gamma
- else if (mother->GetPDG()->Mass() > fMinMotherMass && mother->GetPDG()->Mass() < fMaxMotherMass && //0.4 and 6.1
- TMath::Abs(mother->GetPDG()->Mass()-0.940) > 0.004 && mother->R() < fRadiusCut && abs(mother->Vz()) < fVzCut)
- store = kTRUE; // except nucleons and heavy fragments
- }
- } //no points
- }
- // --> Set storage flag
- fStoreMap[i] = store;
- }
- //AZ If flag is set, flag recursively mothers of selected tracks and their sisters
- if (fStoreMothers) {
- // Fill mother map
- std::multimap<Int_t,Int_t> moths;
- std::multimap<Int_t,Int_t>::iterator it;
- pair<std::multimap<Int_t,Int_t>::iterator,std::multimap<Int_t,Int_t>::iterator> ret;
- std::map<Int_t,Bool_t> copyMap = fStoreMap;
- for (Int_t i=0; i<fNParticles; i++) {
- Int_t iMother = GetParticle(i)->GetMother(0);
- if (iMother >= 0) moths.insert(pair<Int_t,Int_t>(iMother,i));
- }
- for (Int_t i=0; i<fNParticles; i++) {
- if (copyMap[i]) { //if particle is to be stored
- Int_t iMother = GetParticle(i)->GetMother(0);
- while (iMother >= 0) { //and it's mother is not primary
- TParticle* mother = GetParticle(iMother);
- if (mother->GetPDG()->Mass() < 0.1 && mother->P() < 0.01) break;
- //if (mother->GetPDG()->Mass() < 0.1 && mother->P() < 0.001) break;
- if (abs(mother->Vz()) > fVzCut || mother->R() > fRadiusCut) break; //NG and doesn't originate within Rcut/Zcut
- fStoreMap[iMother] = kTRUE;//store the mother
- ret = moths.equal_range(iMother);
- for (it = ret.first; it != ret.second; ++it) {
- TParticle* daught = GetParticle(it->second);
- if (abs(daught->Vz()) > fVzCut || daught->R() > fRadiusCut) continue; //NG and doesn't originate within Rcut/Zcut
- if (daught->GetPDG()->Mass() < 0.1 && daught->P() < 0.01) continue;
- //if (daught->GetPDG()->Mass() < 0.1 && daught->P() < 0.001) continue;
- fStoreMap[it->second] = kTRUE; // sister
- }
- iMother = mother->GetMother(0); //mother of the mother, loop again
- } // while (iMother >= 0)
- } // if (copyMap[i])
- } // for (Int_t i=0; i<fNParticles;
- // AZ - If mother is not stored, find the nearest stored ancestor
- for (Int_t i = 0; i < fNParticles; ++i) {
- if (fStoreMap[i]) {
- Int_t iMother = GetParticle(i)->GetMother(0);
- while (iMother >= 0 && fStoreMap[iMother] == kFALSE) {
- TParticle* mother = GetParticle(iMother);
- iMother = mother->GetMother(0); //mother of the mother, loop again
- }
- if (iMother >= 0) GetParticle(i)->SetFirstMother(iMother);
- }
- }
- } // if (fStoreMothers)
-
- //*/
- }
- // -------------------------------------------------------------------------
- ClassImp(MpdStack)
|