X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=0550ea8e8861f18c6a7845862e54b281c75088ef;hb=b9d1a7e2fd1fffbe9f0d37c43424b0943c915d52;hp=e9159ff3ba926902b9867e171bf4719190ef29e9;hpb=ed54bf31b1a24842d58cefb6e7b115846621aa5b;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index e9159ff3ba9..0550ea8e886 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -13,91 +13,111 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.2 2001/05/17 05:49:39 fca -Reset pointers to daughters - -Revision 1.1 2001/05/16 14:57:22 alibrary -New files for folders and Stack - -*/ +/* $Id$ */ /////////////////////////////////////////////////////////////////////////////// // // -// 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" -//#include "ETrackBits.h" ClassImp(AliStack) -//_____________________________________________________________________________ -AliStack::AliStack(Int_t size) +//_______________________________________________________________________ +AliStack::AliStack(): + fParticles("TParticle", 1000), + fParticleMap(), + fParticleFileMap(0), + fParticleBuffer(0), + fCurrentTrack(0), + fTreeK(0), + fNtrack(0), + fNprimary(0), + fCurrent(-1), + fCurrentPrimary(-1), + fHgwmk(0), + fLoadPoint(0), + fTrackLabelMap(0) { // - // Constructor + // Default constructor // - - // Create the particles arrays - fParticles = new TClonesArray("TParticle",1000); - fParticleMap = new TObjArray(size); - fParticleBuffer = new TParticle(); - fNtrack = 0; - fNprimary = 0; - fCurrent = -1; - fCurrentPrimary = -1; - fTreeK = 0; } - -//_____________________________________________________________________________ -AliStack::AliStack() +//_______________________________________________________________________ +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), + fTrackLabelMap(0) { // - // Default constructor + // Constructor // - - // Create the particles arrays - fParticles = new TClonesArray("TParticle",1000); - fParticleMap = new TObjArray(10000); - fParticleBuffer = new TParticle(); - fNtrack = 0; - fCurrent = -1; - fNprimary = 0; - fCurrentPrimary = -1; - fTreeK = 0; } +//_______________________________________________________________________ +AliStack::AliStack(const AliStack& st): + 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 +} -//_____________________________________________________________________________ + +//_______________________________________________________________________ +void AliStack::Copy(TObject&) const +{ + AliFatal("Not implemented!"); +} + +//_______________________________________________________________________ AliStack::~AliStack() { // // Destructor // - if (fParticles) { - fParticles->Delete(); - delete fParticles; - } - delete fParticleMap; - delete fParticleBuffer; - delete fTreeK; + fParticles.Clear(); } // @@ -105,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, - AliMCProcess mech, Int_t &ntr, Float_t weight) +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 @@ -124,10 +144,6 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, // ntr on output the number of the track stored // - Float_t mass; - const Int_t kfirstdaughter=-1; - const Int_t klastdaughter=-1; - const Int_t kS=0; // const Float_t tlife=0; // @@ -137,58 +153,37 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, // also, this method is potentially dangerous if the mass // used in the MC is not the same of the PDG database // - mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); - Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+ - pmom[1]*pmom[1]+pmom[2]*pmom[2]); - -// printf("Loading mass %f ene %f No %d ip %d parent %d done %d pos %f %f %f mom %f %f %f kS %d m \n", -// mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); - - TClonesArray &particles = *fParticles; - TParticle* particle - = new(particles[fLoadPoint++]) - TParticle(pdg, kS, parent, -1, kfirstdaughter, klastdaughter, - pmom[0], pmom[1], pmom[2], e, vpos[0], vpos[1], vpos[2], tof); - particle->SetPolarisation(TVector3(polar[0],polar[1],polar[2])); - particle->SetWeight(weight); - particle->SetUniqueID(mech); - if(!done) particle->SetBit(kDoneBit); + TParticlePDG* pmc = TDatabasePDG::Instance()->GetParticle(pdg); + if (pmc) { + Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); + Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+ + pmom[1]*pmom[1]+pmom[2]*pmom[2]); + +// printf("Loading mass %f ene %f No %d ip %d parent %d done %d pos %f %f %f mom %f %f %f kS %d m \n", +// mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); - - // Declare that the daughter information is valid - particle->SetBit(kDaughtersBit); - // Add the particle to the stack - fParticleMap->AddAtAndExpand(particle, fNtrack); - if(parent>=0) { - particle = (TParticle*) fParticleMap->At(parent); - particle->SetLastDaughter(fNtrack); - if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); - } - 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++; + 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 { + 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, - AliMCProcess mech, Int_t &ntr, Float_t weight) + TMCProcess mech, Int_t &ntr, Double_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 // kS generation status code @@ -205,262 +200,385 @@ 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 kS=0; const Int_t kFirstDaughter=-1; const Int_t kLastDaughter=-1; - - TClonesArray &particles = *fParticles; + + TParticle* particle - = new(particles[fLoadPoint++]) - TParticle(pdg, kS, parent, -1, kFirstDaughter, kLastDaughter, + = 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 = (TParticle*) fParticleMap->At(parent); - particle->SetLastDaughter(fNtrack); - if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); - } - 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++; } //_____________________________________________________________________________ -void AliStack::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom, - Float_t &e, Float_t *vpos, Float_t *polar, - Float_t &tof) +TParticle* AliStack::PopNextTrack(Int_t& itrack) { // - // Return next track from stack of particles + // Returns next track from stack of particles // TParticle* track = GetNextParticle(); -// cout << "GetNextTrack():" << fCurrent << fNprimary << endl; - - if(track) { - mtrack=fCurrent; - ipart=track->GetPdgCode(); - pmom[0]=track->Px(); - pmom[1]=track->Py(); - pmom[2]=track->Pz(); - e =track->Energy(); - vpos[0]=track->Vx(); - vpos[1]=track->Vy(); - vpos[2]=track->Vz(); - TVector3 pol; - track->GetPolarisation(pol); - polar[0]=pol.X(); - polar[1]=pol.Y(); - polar[2]=pol.Z(); - tof=track->T(); + + if (track) { + itrack = fCurrent; track->SetBit(kDoneBit); - //cout << "Filled params" << endl; } - else - mtrack=-1; + else + itrack = -1; + + fCurrentTrack = track; + return track; +} +//_____________________________________________________________________________ +TParticle* AliStack::PopPrimaryForTracking(Int_t i) +{ + // + // Returns i-th primary particle if it is flagged to be tracked, + // 0 otherwise // - // 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()); + + TParticle* particle = Particle(i); + + if (!particle->TestBit(kDoneBit)) { + fCurrentTrack = particle; + return particle; } - fTimer.Start(); -} - + 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; iResetBit(kDaughtersBit); - part->SetFirstDaughter(-1); - part->SetLastDaughter(-1); + // 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); +// + 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 = (TParticle *)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 = (TParticle*) 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; iGetFirstMother(); - if(parent>=0) { - father = (TParticle *)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 = (TCollection*)next())) { - TIter nexthit(hitList); - AliHit *hit; - while((hit = (AliHit*)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 = (AliModule*)nextmod())) { - detector->RemapTrackHitIDs(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; + } + + for (i = nkeep; i < fNtrack; ++i) fParticleMap[i]=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)); + Int_t toshrink = fNtrack-fHgwmk-1; + fLoadPoint-=toshrink; + + for(i=fLoadPoint; iGetEntries(); - fTreeK->Fill(); - particles[i]=0; - } - for (i=nkeep; i 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; +} - Int_t toshrink = fNtrack-fHgwmk-1; - fLoadPoint-=toshrink; - for(i=fLoadPoint; iRemoveAt(i); +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; - fNtrack=nkeep; - fHgwmk=nkeep-1; - // delete [] map; + 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 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 >= 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; } //_____________________________________________________________________________ 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 = (TParticle*) part; - fParticleFileMap[i]= (Int_t) fTreeK->GetEntries(); - fTreeK->Fill(); - (*fParticleMap)[i]=0; + TParticle *part; + for(Int_t i=0; i(TreeK()->GetEntries()); + TreeK()->Fill(); + fParticleBuffer=0; + 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) { // @@ -471,7 +589,7 @@ void AliStack::FlagTrack(Int_t track) Int_t curr=track; while(1) { - particle=(TParticle*)fParticleMap->At(curr); + particle = GetParticleMapEntry(curr); // If the particle is flagged the three from here upward is saved already if(particle->TestBit(kKeepBit)) return; @@ -491,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; @@ -509,31 +627,40 @@ 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 // - - fParticles->Clear(); - fParticleMap->Clear(); - if (size>0) fParticleMap->Expand(size); + fParticles.Clear(); + fParticleMap.Clear(); + if (size>0) fParticleMap.Expand(size); } //_____________________________________________________________________________ -void AliStack::SetHighWaterMark(Int_t nt) +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; } //_____________________________________________________________________________ @@ -541,47 +668,96 @@ TParticle* AliStack::Particle(Int_t i) { // // Return particle with specified ID - - if(!(*fParticleMap)[i]) { - Int_t nentries = fParticles->GetEntries(); + + if(!fParticleMap.At(i)) { + Int_t nentries = fParticles.GetEntriesFast(); // algorithmic way of getting entry index // (primary particles are filled after secondaries) - Int_t entry; - if (iGetEntry(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); } - return (TParticle *) (*fParticleMap)[i]; + return GetParticleMapEntry(i); } //_____________________________________________________________________________ -Int_t AliStack::GetPrimary(Int_t id) const +TParticle* AliStack::ParticleFromTreeK(Int_t id) const +{ +// +// return pointer to TParticle with label id +// + Int_t entry; + if ((entry = TreeKEntry(id)) < 0) return 0; + if (fTreeK->GetEntry(entry)<=0) return 0; + return fParticleBuffer; +} + +//_____________________________________________________________________________ +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 +// +// 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 (idGetFirstMother(); + else { + AliWarning("Current track not found in the stack"); + return -1; + } +} + +//_____________________________________________________________________________ +Int_t AliStack::GetPrimary(Int_t id) { // // Return number of primary that has generated track // int current, parent; - TParticle *part; // parent=id; while (1) { current=parent; - part = (TParticle *)fParticleMap->At(current); - parent=part->GetFirstMother(); + parent=Particle(current)->GetFirstMother(); if(parent<0) return current; } } @@ -592,8 +768,7 @@ void AliStack::DumpPart (Int_t i) const // // Dumps particle i in the stack // - - ((TParticle*) (*fParticleMap)[i])->Print(); + GetParticleMapEntry(i)->Print(); } //_____________________________________________________________________________ @@ -603,9 +778,10 @@ void AliStack::DumpPStack () // Dumps the particle stack // - printf( - "\n\n=======================================================================\n"); - for (Int_t i=0;i %d ",i); particle->Print(); printf("--------------------------------------------------------------\n"); @@ -653,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 // @@ -665,11 +848,10 @@ void AliStack::CleanParents() // Set parent/daughter relations // - TObjArray &particles = *fParticleMap; TParticle *part; int i; for(i=0; iTestBit(kDaughtersBit)) { part->SetFirstDaughter(-1); part->SetLastDaughter(-1); @@ -689,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 = (TParticle*) fParticleMap->At(i); + particle = GetParticleMapEntry(i); if ((particle) && (!particle->TestBit(kDoneBit))) { fCurrent=i; - //cout << "GetNextParticle() - secondary " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; return particle; } } @@ -701,84 +881,196 @@ TParticle* AliStack::GetNextParticle() // take next primary if all secondaries were done while (fCurrentPrimary>=0) { fCurrent = fCurrentPrimary; - particle = (TParticle*) 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]; - // printf("\n MakeTree called %d", event); - if (!fTreeK) { - sprintf(hname,"TreeK%d",event); - fTreeK = new TTree(hname,"Kinematics"); - // Create a branch for particles - branch = fTreeK->Branch("Particles", "TParticle", &fParticleBuffer, 4000, 1); - 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"); } -void 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 = (TTree*)gDirectory->Get(treeName); + 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; - if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticleBuffer); + 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; +} - else - Error("GetEvent","cannot find Kine Tree for event:%d\n",event); +//_____________________________________________________________________________ +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()) { // -// printf("\n primaries %d", fNprimary); -// printf("\n tracks %d", fNtrack); +// Particle produced by generator + return kTRUE; + } else { // - Int_t size = (Int_t)fTreeK->GetEntries(); - ResetArrays(size); -} +// 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 +// +// 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 ? +}