X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=c86db151e8f227ba132b89f87a8c55d3b2cdaa22;hb=26acf84d05f6b99b3a8bb9d15c1597f9e88a7332;hp=03bfdf84e4d5497446de648222bae9c2f082007f;hpb=e191bb57c96a66bf4f60381b9a3f8c511bc6188e;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index 03bfdf84e4d..c86db151e8f 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -17,33 +17,36 @@ /////////////////////////////////////////////////////////////////////////////// // // -// Particles stack class +// Particles stack class // +// Implements the TMCVirtualStack of the Virtual Monte Carlo // +// Holds the particles transported during simulation // +// Is used to compare results of reconstruction with simulation // +// Author A.Morsch // // // /////////////////////////////////////////////////////////////////////////////// -#include -#include +#include #include +#include #include #include +#include #include +#include -#include "AliHit.h" -#include "AliModule.h" -#include "AliRun.h" -#include "AliMC.h" -#include "AliRunLoader.h" +#include "AliLog.h" #include "AliStack.h" ClassImp(AliStack) //_______________________________________________________________________ AliStack::AliStack(): - fParticles(0), - fParticleMap(0), + fParticles("TParticle", 1000), + fParticleMap(), fParticleFileMap(0), fParticleBuffer(0), + fCurrentTrack(0), fTreeK(0), fNtrack(0), fNprimary(0), @@ -51,7 +54,7 @@ AliStack::AliStack(): fCurrentPrimary(-1), fHgwmk(0), fLoadPoint(0), - fEventFolderName(AliConfig::GetDefaultEventFolderName()) + fTrackLabelMap(0) { // // Default constructor @@ -59,11 +62,12 @@ AliStack::AliStack(): } //_______________________________________________________________________ -AliStack::AliStack(Int_t size, const char* evfoldname): - fParticles(new TClonesArray("TParticle",1000)), - fParticleMap(new TObjArray(size)), +AliStack::AliStack(Int_t size, const char* /*evfoldname*/): + fParticles("TParticle",1000), + fParticleMap(size), fParticleFileMap(0), fParticleBuffer(0), + fCurrentTrack(0), fTreeK(0), fNtrack(0), fNprimary(0), @@ -71,7 +75,7 @@ AliStack::AliStack(Int_t size, const char* evfoldname): fCurrentPrimary(-1), fHgwmk(0), fLoadPoint(0), - fEventFolderName(evfoldname) + fTrackLabelMap(0) { // // Constructor @@ -80,29 +84,29 @@ AliStack::AliStack(Int_t size, const char* evfoldname): //_______________________________________________________________________ AliStack::AliStack(const AliStack& st): - TVirtualMCStack(st), - fParticles(0), - fParticleMap(0), - fParticleFileMap(0), - fParticleBuffer(0), - fTreeK(0), - fNtrack(0), - fNprimary(0), - fCurrent(-1), - fCurrentPrimary(-1), - fHgwmk(0), - fLoadPoint(0) + TVirtualMCStack(st), + fParticles("TParticle",1000), + fParticleMap(*(st.Particles())), + fParticleFileMap(st.fParticleFileMap), + fParticleBuffer(0), + fCurrentTrack(0), + fTreeK((TTree*)(st.fTreeK->Clone())), + fNtrack(st.GetNtrack()), + fNprimary(st.GetNprimary()), + fCurrent(-1), + fCurrentPrimary(-1), + fHgwmk(0), + fLoadPoint(0), + fTrackLabelMap(0) { - // - // Copy constructor - // - st.Copy(*this); + // Copy constructor } + //_______________________________________________________________________ void AliStack::Copy(TObject&) const { - Fatal("Copy","Not implemented!\n"); + AliFatal("Not implemented!"); } //_______________________________________________________________________ @@ -112,12 +116,7 @@ AliStack::~AliStack() // Destructor // - if (fParticles) { - fParticles->Delete(); - delete fParticles; - } - delete fParticleMap; - //PH??? delete fTreeK; + fParticles.Clear(); } // @@ -125,15 +124,15 @@ AliStack::~AliStack() // //_____________________________________________________________________________ -void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, - Float_t *vpos, Float_t *polar, Float_t tof, +void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom, + const Float_t *vpos, const Float_t *polar, Float_t tof, TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) { // // Load a track on the stack // - // done 0 if the track has to be transported - // 1 if not + // done 1 if the track has to be transported + // 0 if not // parent identifier of the parent track. -1 for a primary // pdg particle code // pmom momentum GeV/c @@ -167,8 +166,8 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2], mech, ntr, weight, is); } else { - Warning("PushTrack", "Particle type %d not defined in PDG Database !\n", pdg); - Warning("PushTrack", "Particle skipped !\n"); + AliWarning(Form("Particle type %d not defined in PDG Database !", pdg)); + AliWarning("Particle skipped !"); } } @@ -182,8 +181,8 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, // // Load a track on the stack // - // done 0 if the track has to be transported - // 1 if not + // done 1 if the track has to be transported + // 0 if not // parent identifier of the parent track. -1 for a primary // pdg particle code // kS generation status code @@ -200,45 +199,52 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, // Note: the energy is not calculated from the static mass but // it is passed by argument e. - const Int_t kFirstDaughter=-1; const Int_t kLastDaughter=-1; - - TClonesArray &particles = *fParticles; + + TParticle* particle - = new(particles[fLoadPoint++]) + = new(fParticles[fLoadPoint++]) TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter, px, py, pz, e, vx, vy, vz, tof); - + particle->SetPolarisation(polx, poly, polz); particle->SetWeight(weight); particle->SetUniqueID(mech); - if(!done) particle->SetBit(kDoneBit); + + + if(!done) { + particle->SetBit(kDoneBit); + } else { + particle->SetBit(kTransportBit); + } + + // Declare that the daughter information is valid particle->SetBit(kDaughtersBit); // Add the particle to the stack - fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!! + + fParticleMap.AddAtAndExpand(particle, fNtrack);//CHECK!! if(parent>=0) { - particle = dynamic_cast(fParticleMap->At(parent)); - if (particle) { - particle->SetLastDaughter(fNtrack); - if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); - } - else { - printf("Error in AliStack::PushTrack: Parent %d does not exist\n",parent); - } - } - else { - // - // This is a primary track. Set high water mark for this event - fHgwmk = fNtrack; - // - // Set also number if primary tracks - fNprimary = fHgwmk+1; - fCurrentPrimary++; + particle = GetParticleMapEntry(parent); + if (particle) { + particle->SetLastDaughter(fNtrack); + if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); + } + else { + AliError(Form("Parent %d does not exist",parent)); + } + } else { + // + // This is a primary track. Set high water mark for this event + fHgwmk = fNtrack; + // + // Set also number if primary tracks + fNprimary = fHgwmk+1; + fCurrentPrimary++; } ntr = fNtrack++; } @@ -274,140 +280,211 @@ TParticle* AliStack::PopPrimaryForTracking(Int_t i) TParticle* particle = Particle(i); - if (!particle->TestBit(kDoneBit)) + if (!particle->TestBit(kDoneBit)) { + fCurrentTrack = particle; return particle; + } else return 0; } //_____________________________________________________________________________ -void AliStack::PurifyKine() +Bool_t AliStack::PurifyKine() { // // Compress kinematic tree keeping only flagged particles // and renaming the particle id's in all the hits // - TObjArray &particles = *fParticleMap; - int nkeep=fHgwmk+1, parent, i; + int nkeep = fHgwmk + 1, parent, i; TParticle *part, *father; - TArrayI map(particles.GetLast()+1); + fTrackLabelMap.Set(fParticleMap.GetLast()+1); // Save in Header total number of tracks before compression - // If no tracks generated return now - if(fHgwmk+1 == fNtrack) return; + if(fHgwmk+1 == fNtrack) return kFALSE; // First pass, invalid Daughter information for(i=0; i(particles.At(i)))) { + // Preset map, to be removed later + if(i<=fHgwmk) fTrackLabelMap[i]=i ; + else { + fTrackLabelMap[i] = -99; + if((part=GetParticleMapEntry(i))) { // // Check of this track should be kept for physics reasons - if (KeepPhysics(part)) KeepTrack(i); + if (KeepPhysics(part)) KeepTrack(i); // - part->ResetBit(kDaughtersBit); - part->SetFirstDaughter(-1); - part->SetLastDaughter(-1); + part->ResetBit(kDaughtersBit); + part->SetFirstDaughter(-1); + part->SetLastDaughter(-1); + } } - } } // Invalid daughter information for the parent of the first particle // generated. This may or may not be the current primary according to // whether decays have been recorded among the primaries - part = dynamic_cast(particles.At(fHgwmk+1)); - particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit); + part = GetParticleMapEntry(fHgwmk+1); + fParticleMap.At(part->GetFirstMother())->ResetBit(kDaughtersBit); // Second pass, build map between old and new numbering for(i=fHgwmk+1; iTestBit(kKeepBit)) { - - // This particle has to be kept - map[i]=nkeep; - // If old and new are different, have to move the pointer - if(i!=nkeep) particles[nkeep]=particles.At(i); - part = dynamic_cast(particles.At(nkeep)); - - // as the parent is always *before*, it must be already - // in place. This is what we are checking anyway! - if((parent=part->GetFirstMother())>fHgwmk) - if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent); - else part->SetFirstMother(map[parent]); - - nkeep++; - } + if(fParticleMap.At(i)->TestBit(kKeepBit)) { + // This particle has to be kept + fTrackLabelMap[i]=nkeep; + // If old and new are different, have to move the pointer + if(i!=nkeep) fParticleMap[nkeep]=fParticleMap.At(i); + part = GetParticleMapEntry(nkeep); + // as the parent is always *before*, it must be already + // in place. This is what we are checking anyway! + if((parent=part->GetFirstMother())>fHgwmk) { + if(fTrackLabelMap[parent]==-99) Fatal("PurifyKine","fTrackLabelMap[%d] = -99!\n",parent); + else part->SetFirstMother(fTrackLabelMap[parent]);} + nkeep++; + } } // Fix daughters information for (i=fHgwmk+1; i(particles.At(i)); - parent = part->GetFirstMother(); - if(parent>=0) { - father = dynamic_cast(particles.At(parent)); - if(father->TestBit(kDaughtersBit)) { - - if(iGetFirstDaughter()) father->SetFirstDaughter(i); - if(i>father->GetLastDaughter()) father->SetLastDaughter(i); - } else { - // Initialise daughters info for first pass - father->SetFirstDaughter(i); - father->SetLastDaughter(i); - father->SetBit(kDaughtersBit); + part = GetParticleMapEntry(i); + parent = part->GetFirstMother(); + if(parent>=0) { + father = GetParticleMapEntry(parent); + if(father->TestBit(kDaughtersBit)) { + + if(iGetFirstDaughter()) father->SetFirstDaughter(i); + if(i>father->GetLastDaughter()) father->SetLastDaughter(i); + } else { + // Initialise daughters info for first pass + father->SetFirstDaughter(i); + father->SetLastDaughter(i); + father->SetBit(kDaughtersBit); + } } - } } - - // Now loop on all registered hit lists - TList* hitLists = gAlice->GetMCApp()->GetHitLists(); - TIter next(hitLists); - TCollection *hitList; - while((hitList = dynamic_cast(next()))) { - TIter nexthit(hitList); - AliHit *hit; - while((hit = dynamic_cast(nexthit()))) { - hit->SetTrack(map[hit->GetTrack()]); - } - } - - // - // This for detectors which have a special mapping mechanism - // for hits, such as TPC and TRD // + // Now the output bit, from fHgwmk to nkeep we write everything and we erase + if(nkeep > fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5)); + for (i=fHgwmk+1; i(TreeK()->GetEntries()); + TreeK()->Fill(); + fParticleMap[i]=fParticleBuffer=0; + } - TObjArray* modules = gAlice->Modules(); - TIter nextmod(modules); - AliModule *detector; - while((detector = dynamic_cast(nextmod()))) { - detector->RemapTrackHitIDs(map.GetArray()); - detector->RemapTrackReferencesIDs(map.GetArray()); - } - // - gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray()); - - // Now the output bit, from fHgwmk to nkeep we write everything and we erase - if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5)); - - for (i=fHgwmk+1; i(particles.At(i)); - fParticleFileMap[i]=static_cast(TreeK()->GetEntries()); - TreeK()->Fill(); - particles[i]=fParticleBuffer=0; - } + for (i = nkeep; i < fNtrack; ++i) fParticleMap[i]=0; + + Int_t toshrink = fNtrack-fHgwmk-1; + fLoadPoint-=toshrink; + + for(i=fLoadPoint; iRemoveAt(i); +Bool_t AliStack::ReorderKine() +{ +// +// In some transport code children might not come in a continuous sequence. +// In this case the stack has to be reordered in order to establish the +// mother daughter relation using index ranges. +// + if(fHgwmk+1 == fNtrack) return kFALSE; - fNtrack=nkeep; - fHgwmk=nkeep-1; - // delete [] map; + // + // Howmany secondaries have been produced ? + Int_t nNew = fNtrack - fHgwmk - 1; + + if (nNew > 0) { + Int_t i, j; + TArrayI map1(nNew); + // + // Copy pointers to temporary array + TParticle** tmp = new TParticle*[nNew]; + + for (i = 0; i < nNew; i++) { + if (fParticleMap.At(fHgwmk + 1 + i)) { + tmp[i] = GetParticleMapEntry(fHgwmk + 1 + i); + } else { + tmp[i] = 0x0; + } + map1[i] = -99; + } + + + // + // Reset LoadPoint + // + Int_t loadPoint = fHgwmk + 1; + // + // Re-Push particles into stack + // The outer loop is over parents, the inner over children. + // -1 refers to the primary particle + // + for (i = -1; i < nNew-1; i++) { + Int_t ipa; + TParticle* parP; + if (i == -1) { + ipa = tmp[0]->GetFirstMother(); + parP = GetParticleMapEntry(ipa); + } else { + ipa = (fHgwmk + 1 + i); + // Skip deleted particles + if (!tmp[i]) continue; + // Skip particles without children + if (tmp[i]->GetFirstDaughter() == -1) continue; + parP = tmp[i]; + } + // Reset daughter information + + Int_t idaumin = parP->GetFirstDaughter() - fHgwmk - 1; + Int_t idaumax = parP->GetLastDaughter() - fHgwmk - 1; + parP->SetFirstDaughter(-1); + parP->SetLastDaughter(-1); + for (j = idaumin; j <= idaumax; j++) { + // Skip deleted particles + if (!tmp[j]) continue; + // Skip particles already handled + if (map1[j] != -99) continue; + Int_t jpa = tmp[j]->GetFirstMother(); + // Check if daughter of current parent + if (jpa == ipa) { + fParticleMap[loadPoint] = tmp[j]; + // Re-establish daughter information + parP->SetLastDaughter(loadPoint); + if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(loadPoint); + // Set Mother information + if (i != -1) { + tmp[j]->SetFirstMother(map1[i]); + } + // Build the map + map1[j] = loadPoint; + // Increase load point + loadPoint++; + } + } // children + } // parents + + delete[] tmp; + + // + // Build map for remapping of hits + // + fTrackLabelMap.Set(fNtrack); + for (i = 0; i < fNtrack; i ++) { + if (i <= fHgwmk) { + fTrackLabelMap[i] = i; + } else{ + fTrackLabelMap[i] = map1[i - fHgwmk -1]; + } + } + } // new particles poduced + + return kTRUE; } -Bool_t AliStack::KeepPhysics(TParticle* part) +Bool_t AliStack::KeepPhysics(const TParticle* part) { // // Some particles have to kept on the stack for reasons motivated @@ -419,7 +496,7 @@ Bool_t AliStack::KeepPhysics(TParticle* part) // Int_t parent = part->GetFirstMother(); if (parent >= 0 && parent <= fHgwmk) { - TParticle* father = dynamic_cast(Particles()->At(parent)); + TParticle* father = GetParticleMapEntry(parent); Int_t kf = father->GetPdgCode(); kf = TMath::Abs(kf); Int_t kfl = kf; @@ -444,7 +521,6 @@ void AliStack::FinishEvent() // Update event header - if (!TreeK()) { // Fatal("FinishEvent", "No kinematics tree is defined."); // Don't panic this is a probably a lego run @@ -452,35 +528,35 @@ void AliStack::FinishEvent() } CleanParents(); - if(TreeK()->GetEntries() ==0) { + if(TreeK()->GetEntries() ==0) { // set the fParticleFileMap size for the first time fParticleFileMap.Set(fHgwmk+1); } Bool_t allFilled = kFALSE; - TObject *part; - for(Int_t i=0; iAt(i))) { - fParticleBuffer = dynamic_cast(part); + TParticle *part; + for(Int_t i=0; i(TreeK()->GetEntries()); TreeK()->Fill(); - //PH (*fParticleMap)[i]=fParticleBuffer=0; fParticleBuffer=0; - fParticleMap->AddAt(0,i); + fParticleMap.AddAt(0,i); // When all primaries were filled no particle!=0 // should be left => to be removed later. - if (allFilled) printf("Why != 0 part # %d?\n",i); + if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i)); } else - { + { // // printf("Why = 0 part # %d?\n",i); => We know. // break; // we don't break now in order to be sure there is no // particle !=0 left. // To be removed later and replaced with break. if(!allFilled) allFilled = kTRUE; - } + } + } } //_____________________________________________________________________________ @@ -494,7 +570,7 @@ void AliStack::FlagTrack(Int_t track) Int_t curr=track; while(1) { - particle=dynamic_cast(fParticleMap->At(curr)); + particle = GetParticleMapEntry(curr); // If the particle is flagged the three from here upward is saved already if(particle->TestBit(kKeepBit)) return; @@ -514,14 +590,14 @@ void AliStack::KeepTrack(Int_t track) // Flags a track to be kept // - fParticleMap->At(track)->SetBit(kKeepBit); + fParticleMap.At(track)->SetBit(kKeepBit); } //_____________________________________________________________________________ -void AliStack::Reset(Int_t size) +void AliStack::Clean(Int_t size) { // - // Resets stack + // Reset stack data except for fTreeK // fNtrack=0; @@ -529,26 +605,30 @@ void AliStack::Reset(Int_t size) fHgwmk=0; fLoadPoint=0; fCurrent = -1; - fTreeK = 0x0; ResetArrays(size); } +//_____________________________________________________________________________ +void AliStack::Reset(Int_t size) +{ + // + // Reset stack data including fTreeK + // + + Clean(size); + delete fParticleBuffer; fParticleBuffer = 0; + fTreeK = 0x0; +} + //_____________________________________________________________________________ void AliStack::ResetArrays(Int_t size) { // // Resets stack arrays // - - if (fParticles) - fParticles->Clear(); - else - fParticles = new TClonesArray("TParticle",1000); - if (fParticleMap) { - fParticleMap->Clear(); - if (size>0) fParticleMap->Expand(size);} - else - fParticleMap = new TObjArray(size); + fParticles.Clear(); + fParticleMap.Clear(); + if (size>0) fParticleMap.Expand(size); } //_____________________________________________________________________________ @@ -557,13 +637,11 @@ void AliStack::SetHighWaterMark(Int_t) // // Set high water mark for last track in event // - - fHgwmk = fNtrack-1; - fCurrentPrimary=fHgwmk; - - // Set also number of primary tracks - fNprimary = fHgwmk+1; - fNtrack = fHgwmk+1; + + fHgwmk = fNtrack-1; + fCurrentPrimary=fHgwmk; + // Set also number of primary tracks + fNprimary = fHgwmk+1; } //_____________________________________________________________________________ @@ -571,10 +649,9 @@ TParticle* AliStack::Particle(Int_t i) { // // Return particle with specified ID - - //PH if(!(*fParticleMap)[i]) { - if(!fParticleMap->At(i)) { - Int_t nentries = fParticles->GetEntriesFast(); + + if(!fParticleMap.At(i)) { + Int_t nentries = fParticles.GetEntriesFast(); // algorithmic way of getting entry index // (primary particles are filled after secondaries) Int_t entry = TreeKEntry(i); @@ -582,17 +659,18 @@ TParticle* AliStack::Particle(Int_t i) // and the fParticleFileMap[i] give the same; // give the fatal error if not if (entry != fParticleFileMap[i]) { - Fatal("Particle", + AliFatal(Form( "!! The algorithmic way and map are different: !!\n entry: %d map: %d", - entry, fParticleFileMap[i]); + entry, fParticleFileMap[i])); } - + // Load particle at entry into fParticleBuffer TreeK()->GetEntry(entry); - new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); - fParticleMap->AddAt((*fParticles)[nentries],i); + // Add to the TClonesarray + new (fParticles[nentries]) TParticle(*fParticleBuffer); + // Store a pointer in the TObjArray + fParticleMap.AddAt(fParticles[nentries],i); } - //PH return dynamic_cast((*fParticleMap)[i]); - return dynamic_cast(fParticleMap->At(i)); + return GetParticleMapEntry(i); } //_____________________________________________________________________________ @@ -611,9 +689,18 @@ TParticle* AliStack::ParticleFromTreeK(Int_t id) const Int_t AliStack::TreeKEntry(Int_t id) const { // -// return entry number in the TreeK for particle with label id -// return negative number if label>fNtrack +// Return entry number in the TreeK for particle with label id +// Return negative number if label>fNtrack +// +// The order of particles in TreeK reflects the order of the transport of primaries and production of secondaries: // +// Before transport there are fNprimary particles on the stack. +// They are transported one by one and secondaries (fNtrack - fNprimary) are produced. +// After the transport of each particles secondaries are written to the TreeK +// They occupy the entries 0 ... fNtrack - fNprimary - 1 +// The primaries are written after they have been transported and occupy +// fNtrack - fNprimary .. fNtrack - 1 + Int_t entry; if (idAt(fCurrent); + TParticle* current = GetParticleMapEntry(fCurrent); if (current) return current->GetFirstMother(); else { - Warning("GetCurrentParentTrackNumber", "Current track not found in the stack"); + AliWarning("Current track not found in the stack"); return -1; } } @@ -662,9 +749,7 @@ void AliStack::DumpPart (Int_t i) const // // Dumps particle i in the stack // - - //PH dynamic_cast((*fParticleMap)[i])->Print(); - dynamic_cast(fParticleMap->At(i))->Print(); + GetParticleMapEntry(i)->Print(); } //_____________________________________________________________________________ @@ -691,9 +776,9 @@ void AliStack::DumpPStack () printf("\n=======================================================================\n\n"); // print particle file map - printf("\nParticle file map: \n"); - for (i=0; i(particles[i]); + TParticle* particle = GetParticleMapEntry(i); if (particle) { printf("-> %d ",i); particle->Print(); printf("--------------------------------------------------------------\n"); @@ -724,6 +808,15 @@ void AliStack::DumpLoadedStack() const "\n=======================================================================\n\n"); } +//_____________________________________________________________________________ +void AliStack::SetCurrentTrack(Int_t track) +{ + fCurrent = track; + if (fCurrent < fNprimary) fCurrentTrack = Particle(track); +} + + +//_____________________________________________________________________________ // // protected methods // @@ -736,11 +829,10 @@ void AliStack::CleanParents() // Set parent/daughter relations // - TObjArray &particles = *fParticleMap; TParticle *part; int i; for(i=0; i(particles.At(i)); + part = GetParticleMapEntry(i); if(part) if(!part->TestBit(kDaughtersBit)) { part->SetFirstDaughter(-1); part->SetLastDaughter(-1); @@ -760,7 +852,7 @@ TParticle* AliStack::GetNextParticle() // search secondaries //for(Int_t i=fNtrack-1; i>=0; i--) { for(Int_t i=fNtrack-1; i>fHgwmk; i--) { - particle = dynamic_cast(fParticleMap->At(i)); + particle = GetParticleMapEntry(i); if ((particle) && (!particle->TestBit(kDoneBit))) { fCurrent=i; return particle; @@ -770,7 +862,7 @@ TParticle* AliStack::GetNextParticle() // take next primary if all secondaries were done while (fCurrentPrimary>=0) { fCurrent = fCurrentPrimary; - particle = dynamic_cast(fParticleMap->At(fCurrentPrimary--)); + particle = GetParticleMapEntry(fCurrentPrimary--); if ((particle) && (!particle->TestBit(kDoneBit))) { return particle; } @@ -778,54 +870,26 @@ TParticle* AliStack::GetNextParticle() // nothing to be tracked fCurrent = -1; + + return particle; } //__________________________________________________________________________________________ -TTree* AliStack::TreeK() -{ -//returns TreeK - if (fTreeK) - { - return fTreeK; - } - else - { - AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName); - if (rl == 0x0) - { - Fatal("TreeK","Can not get RunLoader from event folder named %s",fEventFolderName.Data()); - return 0x0;//pro forma - } - fTreeK = rl->TreeK(); - if ( fTreeK ) - { - ConnectTree(); - } - else - { - //don't panic - could be Lego - if (AliLoader::GetDebug()) - { - Warning("TreeK","Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data()); - } - } - } - return fTreeK;//never reached -} -//__________________________________________________________________________________________ - -void AliStack::ConnectTree() +void AliStack::ConnectTree(TTree* tree) { // // Creates branch for writing particles // - if (AliLoader::GetDebug()) Info("ConnectTree","Connecting TreeK"); + + fTreeK = tree; + + AliDebug(1, "Connecting TreeK"); if (fTreeK == 0x0) { if (TreeK() == 0x0) { - Fatal("ConnectTree","Parameter is NULL");//we don't like such a jokes + AliFatal("Parameter is NULL");//we don't like such a jokes return; } return;//in this case TreeK() calls back this method (ConnectTree) @@ -835,50 +899,34 @@ void AliStack::ConnectTree() // Create a branch for particles - if (AliLoader::GetDebug()) - Info("ConnectTree","Tree name is %s",fTreeK->GetName()); + AliDebug(2, Form("Tree name is %s",fTreeK->GetName())); if (fTreeK->GetDirectory()) { - if (AliLoader::GetDebug()) - Info("ConnectTree","and dir is %s",fTreeK->GetDirectory()->GetName()); + AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName())); } else - Warning("ConnectTree","DIR IS NOT SET !!!"); + AliWarning("DIR IS NOT SET !!!"); - TBranch *branch=fTreeK->GetBranch(AliRunLoader::fgkKineBranchName); + TBranch *branch=fTreeK->GetBranch("Particles"); if(branch == 0x0) { - branch = fTreeK->Branch(AliRunLoader::fgkKineBranchName, "TParticle", &fParticleBuffer, 4000); - if (AliLoader::GetDebug()) Info("ConnectTree","Creating Branch in Tree"); + branch = fTreeK->Branch("Particles", "TParticle", &fParticleBuffer, 4000); + AliDebug(2, "Creating Branch in Tree"); } else { - if (AliLoader::GetDebug()) Info("ConnectTree","Branch Found in Tree"); + AliDebug(2, "Branch Found in Tree"); branch->SetAddress(&fParticleBuffer); } if (branch->GetDirectory()) { - if (AliLoader::GetDebug()) - Info("ConnectTree","Branch Dir Name is %s",branch->GetDirectory()->GetName()); + AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName())); } else - Warning("ConnectTree","Branch Dir is NOT SET"); -} -//__________________________________________________________________________________________ - - -void AliStack::BeginEvent() -{ -// start a new event - Reset(); + AliWarning("Branch Dir is NOT SET"); } -//_____________________________________________________________________________ -void AliStack::FinishRun() -{ -// Clean TreeK information -} //_____________________________________________________________________________ Bool_t AliStack::GetEvent() @@ -887,22 +935,117 @@ Bool_t AliStack::GetEvent() // Get new event from TreeK // Reset/Create the particle stack - fTreeK = 0x0; - - if (TreeK() == 0x0) //forces connecting - { - Error("GetEvent","cannot find Kine Tree for current event\n"); - return kFALSE; - } - Int_t size = (Int_t)TreeK()->GetEntries(); ResetArrays(size); return kTRUE; } //_____________________________________________________________________________ -void AliStack::SetEventFolderName(const char* foldname) +Bool_t AliStack::IsStable(Int_t pdg) const { - //Sets event folder name - fEventFolderName = foldname; + // + // Decide whether particle (pdg) is stable + // + + + // All ions/nucleons are considered as stable + // Nuclear code is 10LZZZAAAI + if(pdg>1000000000)return kTRUE; + + const Int_t kNstable = 15; + Int_t i; + + Int_t pdgStable[kNstable] = { + kGamma, // Photon + kElectron, // Electron + kMuonPlus, // Muon + kPiPlus, // Pion + kKPlus, // Kaon + kK0Short, // K0s + kK0Long, // K0l + kProton, // Proton + kNeutron, // Neutron + kLambda0, // Lambda_0 + kSigmaMinus, // Sigma Minus + kSigmaPlus, // Sigma Plus + 3312, // Xsi Minus + 3322, // Xsi + 3334, // Omega + }; + + Bool_t isStable = kFALSE; + for (i = 0; i < kNstable; i++) { + if (pdg == TMath::Abs(pdgStable[i])) { + isStable = kTRUE; + break; + } + } + + return isStable; } + +//_____________________________________________________________________________ +Bool_t AliStack::IsPhysicalPrimary(Int_t index) +{ + // + // Test if a particle is a physical primary according to the following definition: + // Particles produced in the collision including products of strong and + // electromagnetic decay and excluding feed-down from weak decays of strange + // particles. + // + TParticle* p = Particle(index); + Int_t ist = p->GetStatusCode(); + + // + // Initial state particle + if (ist > 1) return kFALSE; + + Int_t pdg = TMath::Abs(p->GetPdgCode()); + + if (!IsStable(pdg)) return kFALSE; + + if (index < GetNprimary()) { +// +// Particle produced by generator + return kTRUE; + } else { +// +// Particle produced during transport +// + + Int_t imo = p->GetFirstMother(); + TParticle* pm = Particle(imo); + Int_t mpdg = TMath::Abs(pm->GetPdgCode()); +// Check if it comes from a pi0 decay +// +// What about the pi0 Dalitz ?? +// if ((mpdg == kPi0) && (imo < GetNprimary())) return kTRUE; + +// Check if this is a heavy flavor decay product + Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg)))); + // + // Light hadron + if (mfl < 4) return kFALSE; + + // + // Heavy flavor hadron produced by generator + if (imo < GetNprimary()) { + return kTRUE; + } + + // To be sure that heavy flavor has not been produced in a secondary interaction + // Loop back to the generated mother + while (imo >= GetNprimary()) { + imo = pm->GetFirstMother(); + pm = Particle(imo); + } + mpdg = TMath::Abs(pm->GetPdgCode()); + mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg)))); + + if (mfl < 4) { + return kFALSE; + } else { + return kTRUE; + } + } // produced by generator ? +}