X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=bcedbac77186e2b802ace93a255330d5615cc5d6;hb=fb3430ac68b6928f9e0d35e74bf71f65785c98a5;hp=4d1040f7df396564f92723db7003bf19a4f5f2cd;hpb=fe046ade3136e16900c79722861d96d4594459e0;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index 4d1040f7df3..bcedbac7718 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -17,40 +17,45 @@ /////////////////////////////////////////////////////////////////////////////// // // -// 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 -#include +#include +#include "AliLog.h" #include "AliStack.h" -#include "AliRun.h" -#include "AliModule.h" -#include "AliHit.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), fCurrent(-1), fCurrentPrimary(-1), fHgwmk(0), - fLoadPoint(0) + fLoadPoint(0), + fTrackLabelMap(0) { // // Default constructor @@ -58,18 +63,20 @@ AliStack::AliStack(): } //_______________________________________________________________________ -AliStack::AliStack(Int_t size): - 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), fCurrent(-1), fCurrentPrimary(-1), fHgwmk(0), - fLoadPoint(0) + fLoadPoint(0), + fTrackLabelMap(0) { // // Constructor @@ -78,29 +85,29 @@ AliStack::AliStack(Int_t size): //_______________________________________________________________________ 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(AliStack&) const +void AliStack::Copy(TObject&) const { - Fatal("Copy","Not implemented!\n"); + AliFatal("Not implemented!"); } //_______________________________________________________________________ @@ -110,12 +117,7 @@ AliStack::~AliStack() // Destructor // - if (fParticles) { - fParticles->Delete(); - delete fParticles; - } - delete fParticleMap; - delete fTreeK; + fParticles.Clear(); } // @@ -123,15 +125,15 @@ AliStack::~AliStack() // //_____________________________________________________________________________ -void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, - Float_t *vpos, Float_t *polar, Float_t tof, - TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) +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 @@ -161,17 +163,17 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, // mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); - SetTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e, + PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e, vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2], mech, ntr, weight, is); } else { - Warning("SetTrack", "Particle type %d not defined in PDG Database !\n", pdg); - Warning("SetTrack", "Particle skipped !\n"); + AliWarning(Form("Particle type %d not defined in PDG Database !", pdg)); + AliWarning("Particle skipped !"); } } //_____________________________________________________________________________ -void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, +void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t e, Double_t vx, Double_t vy, Double_t vz, Double_t tof, Double_t polx, Double_t poly, Double_t polz, @@ -180,8 +182,8 @@ void AliStack::SetTrack(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 @@ -198,51 +200,58 @@ void AliStack::SetTrack(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::SetTrack: 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++; } //_____________________________________________________________________________ -TParticle* AliStack::GetNextTrack(Int_t& itrack) +TParticle* AliStack::PopNextTrack(Int_t& itrack) { // // Returns next track from stack of particles @@ -262,58 +271,8 @@ TParticle* AliStack::GetNextTrack(Int_t& itrack) return track; } -/* //_____________________________________________________________________________ -void AliStack::GetNextTrack(Int_t& itrack, Int_t& pdg, - Double_t& px, Double_t& py, Double_t& pz, Double_t& e, - Double_t& vx, Double_t& vy, Double_t& vz, Double_t& tof, - Double_t& polx, Double_t& poly, Double_t& polz) -{ - // - // Return next track from stack of particles - // - - - TParticle* track = GetNextParticle(); -// cout << "GetNextTrack():" << fCurrent << fNprimary << endl; - - if (track) { - itrack = fCurrent; - pdg = track->GetPdgCode(); - px = track->Px(); - py = track->Py(); - pz = track->Pz(); - e = track->Energy(); - vx = track->Vx(); - vy = track->Vy(); - vz = track->Vz(); - tof = track->T(); - TVector3 pol; - track->GetPolarisation(pol); - polx = pol.X(); - poly = pol.Y(); - polz = pol.Z(); - track->SetBit(kDoneBit); -// cout << "Filled params" << endl; - } - else - itrack = -1; - - // - // stop and start timer when we start a primary track - Int_t nprimaries = fNprimary; - if (fCurrent >= nprimaries) return; - if (fCurrent < nprimaries-1) { - fTimer.Stop(); - track=(TParticle*) fParticleMap->At(fCurrent+1); - // track->SetProcessTime(fTimer.CpuTime()); - } - fTimer.Start(); -} - -*/ -//_____________________________________________________________________________ -TParticle* AliStack::GetPrimaryForTracking(Int_t i) +TParticle* AliStack::PopPrimaryForTracking(Int_t i) { // // Returns i-th primary particle if it is flagged to be tracked, @@ -322,163 +281,252 @@ TParticle* AliStack::GetPrimaryForTracking(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->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 // - - TObjArray* modules = gAlice->Modules(); - TIter nextmod(modules); - AliModule *detector; - while((detector = dynamic_cast(nextmod()))) { - detector->RemapTrackHitIDs(map.GetArray()); - detector->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(TreeK()->GetEntries()); + TreeK()->Fill(); + fParticleMap[i]=fParticleBuffer=0; + } - // 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(fTreeK->GetEntries()); - fTreeK->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 // by physics analysis. Decision is put here. // Bool_t keep = kFALSE; + + Int_t parent = part->GetFirstMother(); + if (parent >= 0 && parent <= fHgwmk) { + TParticle* father = GetParticleMapEntry(parent); // // Keep first-generation daughter from primaries with heavy flavor // - Int_t parent = part->GetFirstMother(); - if (parent >= 0 && parent <= fHgwmk) { - TParticle* father = dynamic_cast(Particles()->At(parent)); Int_t kf = father->GetPdgCode(); kf = TMath::Abs(kf); Int_t kfl = kf; // meson ? if (kfl > 10) kfl/=100; // baryon - if (kfl > 10) kfl/=10; - if (kfl > 10) kfl/=10; + if (kfl > 10) kfl/=10; + if (kfl > 10) kfl/=10; if (kfl >= 4) { keep = kTRUE; } + // + // e+e- from pair production of primary gammas + // + if ((part->GetUniqueID()) == kPPair) keep = kTRUE; + } + // + // Decay(cascade) from primaries + // + if ((part->GetUniqueID() == kPDecay) && (parent >= 0)) { + // Particles from decay + TParticle* father = GetParticleMapEntry(parent); + Int_t imo = parent; + while((imo > fHgwmk) && (father->GetUniqueID() == kPDecay)) { + imo = father->GetFirstMother(); + father = GetParticleMapEntry(imo); + } + if ((imo <= fHgwmk)) keep = kTRUE; } return keep; } @@ -486,54 +534,51 @@ Bool_t AliStack::KeepPhysics(TParticle* part) //_____________________________________________________________________________ void AliStack::FinishEvent() { - // - // Write out the kinematics that was not yet filled - // +// +// Write out the kinematics that was not yet filled +// - // Update event header +// Update event header - - if (!fTreeK) { + if (!TreeK()) { // Fatal("FinishEvent", "No kinematics tree is defined."); // Don't panic this is a probably a lego run return; - } CleanParents(); - if(fTreeK->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); - fParticleFileMap[i]= static_cast(fTreeK->GetEntries()); - fTreeK->Fill(); - //PH (*fParticleMap)[i]=fParticleBuffer=0; + TParticle *part; + for(Int_t i=0; i(TreeK()->GetEntries()); + TreeK()->Fill(); 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 { + 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; - } -// cout << "Nof particles: " << fNtrack << endl; - //Reset(); + // 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; + } + } } - //_____________________________________________________________________________ + void AliStack::FlagTrack(Int_t track) { // @@ -544,7 +589,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; @@ -564,16 +609,16 @@ 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; fNprimary=0; fHgwmk=0; @@ -582,22 +627,27 @@ void AliStack::Reset(Int_t size) 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); } //_____________________________________________________________________________ @@ -606,13 +656,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; } //_____________________________________________________________________________ @@ -620,10 +668,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); @@ -631,17 +678,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])); } - - fTreeK->GetEntry(entry); - new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); - fParticleMap->AddAt((*fParticles)[nentries],i); + // Load particle at entry into fParticleBuffer + TreeK()->GetEntry(entry); + // 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); } //_____________________________________________________________________________ @@ -660,9 +708,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("CurrentTrackParent", "Current track not found in the stack"); + AliWarning("Current track not found in the stack"); return -1; } } @@ -711,9 +768,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(); } //_____________________________________________________________________________ @@ -725,8 +780,7 @@ void AliStack::DumpPStack () Int_t i; - printf( - "\n\n=======================================================================\n"); + printf("\n\n=======================================================================\n"); for (i=0;i(particles[i]); + TParticle* particle = GetParticleMapEntry(i); if (particle) { printf("-> %d ",i); particle->Print(); printf("--------------------------------------------------------------\n"); @@ -775,6 +827,15 @@ void AliStack::DumpLoadedStack() const "\n=======================================================================\n\n"); } +//_____________________________________________________________________________ +void AliStack::SetCurrentTrack(Int_t track) +{ + fCurrent = track; + if (fCurrent < fNprimary) fCurrentTrack = Particle(track); +} + + +//_____________________________________________________________________________ // // protected methods // @@ -787,11 +848,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); @@ -811,11 +871,9 @@ 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; - //cout << "GetNextParticle() - secondary " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; return particle; } } @@ -823,90 +881,195 @@ 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))) { - //cout << "GetNextParticle() - primary " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; return particle; } } // nothing to be tracked fCurrent = -1; - //cout << "GetNextParticle() - none " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; + + return particle; } - //__________________________________________________________________________________________ -void AliStack::MakeTree(Int_t event, const char * /*file*/) -{ -// -// Make Kine tree and creates branch for writing particles -// - TBranch *branch=0; - // Make Kinematics Tree - char hname[30]; - if (!fTreeK) { - sprintf(hname,"TreeK%d",event); - fTreeK = new TTree(hname,"Kinematics"); - // Create a branch for particles - branch = fTreeK->Branch("Particles", "TParticle", &fParticleBuffer, 4000); - fTreeK->Write(0,TObject::kOverwrite); - } -} -//_____________________________________________________________________________ -void AliStack::BeginEvent(Int_t event) +void AliStack::ConnectTree(TTree* tree) { -// start a new event // +// Creates branch for writing particles // - fNprimary = 0; - fNtrack = 0; + + fTreeK = tree; - char hname[30]; - if(fTreeK) { - fTreeK->Reset(); - sprintf(hname,"TreeK%d",event); - fTreeK->SetName(hname); - } -} + AliDebug(1, "Connecting TreeK"); + if (fTreeK == 0x0) + { + if (TreeK() == 0x0) + { + AliFatal("Parameter is NULL");//we don't like such a jokes + return; + } + return;//in this case TreeK() calls back this method (ConnectTree) + //tree after setting fTreeK, the rest was already executed + //it is safe to return now + } -//_____________________________________________________________________________ -void AliStack::FinishRun() -{ -// Clean TreeK information - if (fTreeK) { - delete fTreeK; fTreeK = 0; - } + // Create a branch for particles + + AliDebug(2, Form("Tree name is %s",fTreeK->GetName())); + + if (fTreeK->GetDirectory()) + { + AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName())); + } + else + AliWarning("DIR IS NOT SET !!!"); + + TBranch *branch=fTreeK->GetBranch("Particles"); + if(branch == 0x0) + { + branch = fTreeK->Branch("Particles", &fParticleBuffer, 4000); + AliDebug(2, "Creating Branch in Tree"); + } + else + { + AliDebug(2, "Branch Found in Tree"); + branch->SetAddress(&fParticleBuffer); + } + if (branch->GetDirectory()) + { + AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName())); + } + else + AliWarning("Branch Dir is NOT SET"); } //_____________________________________________________________________________ -Bool_t AliStack::GetEvent(Int_t event) + +Bool_t AliStack::GetEvent() { // // Get new event from TreeK // Reset/Create the particle stack - if (fTreeK) delete fTreeK; - - // Get Kine Tree from file - char treeName[20]; - sprintf(treeName,"TreeK%d",event); - fTreeK = dynamic_cast(gDirectory->Get(treeName)); - - if (fTreeK) - fTreeK->SetBranchAddress("Particles", &fParticleBuffer); - else { - // Error("GetEvent","cannot find Kine Tree for event:%d\n",event); - Warning("GetEvent","cannot find Kine Tree for event:%d\n",event); - return kFALSE; - } -// printf("\n primaries %d", fNprimary); -// printf("\n tracks %d", fNtrack); - - Int_t size = static_cast(fTreeK->GetEntries()); + Int_t size = (Int_t)TreeK()->GetEntries(); ResetArrays(size); return kTRUE; } +//_____________________________________________________________________________ + +Bool_t AliStack::IsStable(Int_t pdg) const +{ + // + // 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 = 18; + 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 + kNuE, // Electron Neutrino + kNuMu, // Muon Neutrino + kNuTau // Tau Neutrino + }; + + 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 for Sigma0 + if ((mpdg == 3212) && (imo < GetNprimary())) return kTRUE; +// +// Check if it comes from a pi0 decay +// + 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 ? +}