X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=747481180ce3d4ddf9fb18433b60ca86fbb65a58;hb=07d955deedabdc3187198fc7fc09e02b1f1f7064;hp=ac9e0bd015c5cae20c2db09ddcfef9a0c2f71b17;hpb=17f4a3d927d84b025145c276c41d8fdb8d28d1c5;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index ac9e0bd015c..747481180ce 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -13,100 +13,119 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.5 2001/05/30 12:18:46 hristov -Loop variables declared once - -Revision 1.4 2001/05/25 07:25:20 hristov -AliStack destructor corrected (I.Hrivnacova) - -Revision 1.3 2001/05/22 14:33:16 hristov -Minor changes - -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 "AliStack.h" -#include "AliRun.h" -#include "AliModule.h" +#include "AliLog.h" #include "AliHit.h" -//#include "ETrackBits.h" +#include "AliModule.h" +#include "AliRun.h" +#include "AliMC.h" +#include "AliRunLoader.h" +#include "AliStack.h" ClassImp(AliStack) -//_____________________________________________________________________________ -AliStack::AliStack(Int_t size) +//_______________________________________________________________________ +AliStack::AliStack(): + fParticles(0), + fParticleMap(0), + fParticleFileMap(0), + fParticleBuffer(0), + fCurrentTrack(0), + fTreeK(0), + fNtrack(0), + fNprimary(0), + fCurrent(-1), + fCurrentPrimary(-1), + fHgwmk(0), + fLoadPoint(0), + fEventFolderName(AliConfig::GetDefaultEventFolderName()) { // - // Constructor + // Default constructor // - - // Create the particles arrays - fParticles = new TClonesArray("TParticle",1000); - fParticleMap = new TObjArray(size); - fParticleBuffer = 0; - fNtrack = 0; - fNprimary = 0; - fCurrent = -1; - fCurrentPrimary = -1; - fTreeK = 0; } - -//_____________________________________________________________________________ -AliStack::AliStack() +//_______________________________________________________________________ +AliStack::AliStack(Int_t size, const char* evfoldname): + fParticles(new TClonesArray("TParticle",1000)), + fParticleMap(new TObjArray(size)), + fParticleFileMap(0), + fParticleBuffer(0), + fCurrentTrack(0), + fTreeK(0), + fNtrack(0), + fNprimary(0), + fCurrent(-1), + fCurrentPrimary(-1), + fHgwmk(0), + fLoadPoint(0), + fEventFolderName(evfoldname) { // - // Default constructor + // Constructor // - - // Create the particles arrays - fParticles = new TClonesArray("TParticle",1000); - fParticleMap = new TObjArray(10000); - fParticleBuffer = 0; - fNtrack = 0; - fCurrent = -1; - fNprimary = 0; - fCurrentPrimary = -1; - fTreeK = 0; } +//_______________________________________________________________________ +AliStack::AliStack(const AliStack& st): + TVirtualMCStack(st), + fParticles(new TClonesArray("TParticle",1000)), + fParticleMap(new TObjArray(*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), + fEventFolderName(0) +{ + // Copy constructor + ConnectTree(); +} -//_____________________________________________________________________________ + +//_______________________________________________________________________ +void AliStack::Copy(TObject&) const +{ + AliFatal("Not implemented!"); +} + +//_______________________________________________________________________ AliStack::~AliStack() { // // Destructor // - delete fParticleBuffer; if (fParticles) { fParticles->Delete(); delete fParticles; } delete fParticleMap; - delete fTreeK; + //PH??? delete fTreeK; } // @@ -114,9 +133,9 @@ 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, Float_t *pmom, + Float_t *vpos, Float_t *polar, Float_t tof, + TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) { // // Load a track on the stack @@ -133,10 +152,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; // @@ -146,52 +161,31 @@ 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 @@ -215,14 +209,14 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, // 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, + TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter, px, py, pz, e, vx, vy, vz, tof); particle->SetPolarisation(polx, poly, polz); @@ -234,12 +228,19 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, // Declare that the daughter information is valid particle->SetBit(kDaughtersBit); // Add the particle to the stack + + fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!! if(parent>=0) { - particle = (TParticle*) fParticleMap->At(parent); - particle->SetLastDaughter(fNtrack); - if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); + particle = dynamic_cast(fParticleMap->At(parent)); + if (particle) { + particle->SetLastDaughter(fNtrack); + if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); + } + else { + AliError(Form("Parent %d does not exist",parent)); + } } else { // @@ -254,52 +255,41 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, } //_____________________________________________________________________________ -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(); - track->SetBit(kDoneBit); - //cout << "Filled params" << endl; - } - else - mtrack=-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()); + if (track) { + itrack = fCurrent; + track->SetBit(kDoneBit); } - fTimer.Start(); + 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 + // + + TParticle* particle = Particle(i); + + if (!particle->TestBit(kDoneBit)) + return particle; + else + return 0; +} //_____________________________________________________________________________ void AliStack::PurifyKine() @@ -315,17 +305,20 @@ void AliStack::PurifyKine() TArrayI map(particles.GetLast()+1); // Save in Header total number of tracks before compression - // If no tracks generated return now if(fHgwmk+1 == fNtrack) return; // First pass, invalid Daughter information - for(i=0; i(particles.At(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); @@ -335,7 +328,7 @@ void AliStack::PurifyKine() // 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); + part = dynamic_cast(particles.At(fHgwmk+1)); particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit); // Second pass, build map between old and new numbering for(i=fHgwmk+1; i(particles.At(nkeep)); // as the parent is always *before*, it must be already // in place. This is what we are checking anyway! @@ -359,10 +352,10 @@ void AliStack::PurifyKine() // Fix daughters information for (i=fHgwmk+1; i(particles.At(i)); parent = part->GetFirstMother(); if(parent>=0) { - father = (TParticle *)particles.At(parent); + father = dynamic_cast(particles.At(parent)); if(father->TestBit(kDaughtersBit)) { if(iGetFirstDaughter()) father->SetFirstDaughter(i); @@ -375,16 +368,19 @@ void AliStack::PurifyKine() } } } + // Now loop on all registered hit lists - TList* hitLists = gAlice->GetHitLists(); + TList* hitLists = gAlice->GetMCApp()->GetHitLists(); TIter next(hitLists); TCollection *hitList; - while((hitList = (TCollection*)next())) { + + while((hitList = dynamic_cast(next()))) { TIter nexthit(hitList); AliHit *hit; - while((hit = (AliHit*)nexthit())) { - hit->SetTrack(map[hit->GetTrack()]); + + while((hit = dynamic_cast(nexthit()))) { + hit->SetTrack(map[hit->GetTrack()]); } } @@ -392,54 +388,213 @@ void AliStack::PurifyKine() // 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())) { + 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; iGetEntries(); - fTreeK->Fill(); + fParticleBuffer = dynamic_cast(particles.At(i)); + fParticleFileMap[i]=static_cast(TreeK()->GetEntries()); + TreeK()->Fill(); particles[i]=fParticleBuffer=0; - } + } for (i=nkeep; iRemoveAt(i); fNtrack=nkeep; fHgwmk=nkeep-1; - // delete [] map; } -//_____________________________________________________________________________ -void AliStack::FinishEvent() +void 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; + // - // Write out the kinematics that was not yet filled + // Howmany secondaries have been produced ? + Int_t nNew = fNtrack - fHgwmk - 1; + + if (nNew > 0) { + Int_t i, j; + TObjArray &particles = *fParticleMap; + TArrayI map1(nNew); + // + // Copy pointers to temporary array + TParticle** tmp = new TParticle*[nNew]; + + for (i = 0; i < nNew; i++) { + if (particles.At(fHgwmk + 1 + i)) { + tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i)); + } else { + tmp[i] = 0x0; + } + map1[i] = -99; + } + + + // + // Reset LoadPoint + // + fLoadPoint = 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 =dynamic_cast(particles.At(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) { + particles[fLoadPoint] = tmp[j]; + // Re-establish daughter information + parP->SetLastDaughter(fLoadPoint); + if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(fLoadPoint); + // Set Mother information + if (i != -1) { + tmp[j]->SetFirstMother(map1[i]); + } + // Build the map + map1[j] = fLoadPoint; + // Increase load point + fLoadPoint++; + } + } // children + } // parents + + delete[] tmp; + + // + // Build map for remapping of hits + // + TArrayI map(fNtrack); + for (i = 0; i < fNtrack; i ++) { + if (i <= fHgwmk) { + map[i] = i; + } else{ + map[i] = map1[i - fHgwmk -1]; + } + } + + // 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 // - // Update event header + 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()); + } // new particles poduced +} + +Bool_t AliStack::KeepPhysics(TParticle* part) +{ + // + // Some particles have to kept on the stack for reasons motivated + // by physics analysis. Decision is put here. + // + Bool_t keep = kFALSE; + // + // 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 >= 4) { + keep = kTRUE; + } + } + return keep; +} + +//_____________________________________________________________________________ +void AliStack::FinishEvent() +{ +// +// Write out the kinematics that was not yet filled +// + +// 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); } @@ -448,28 +603,29 @@ void AliStack::FinishEvent() TObject *part; for(Int_t i=0; iAt(i))) { - fParticleBuffer = (TParticle*) part; - fParticleFileMap[i]= (Int_t) fTreeK->GetEntries(); - fTreeK->Fill(); - (*fParticleMap)[i]=fParticleBuffer=0; + fParticleBuffer = dynamic_cast(part); + fParticleFileMap[i]= static_cast(TreeK()->GetEntries()); + TreeK()->Fill(); + //PH (*fParticleMap)[i]=fParticleBuffer=0; + 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) { // @@ -480,7 +636,7 @@ void AliStack::FlagTrack(Int_t track) Int_t curr=track; while(1) { - particle=(TParticle*)fParticleMap->At(curr); + particle=dynamic_cast(fParticleMap->At(curr)); // If the particle is flagged the three from here upward is saved already if(particle->TestBit(kKeepBit)) return; @@ -509,12 +665,13 @@ void AliStack::Reset(Int_t size) // // Resets stack // - + fNtrack=0; fNprimary=0; fHgwmk=0; fLoadPoint=0; fCurrent = -1; + fTreeK = 0x0; ResetArrays(size); } @@ -525,21 +682,27 @@ void AliStack::ResetArrays(Int_t size) // Resets stack arrays // - fParticles->Clear(); - fParticleMap->Clear(); - if (size>0) fParticleMap->Expand(size); + 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); } //_____________________________________________________________________________ -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; @@ -551,46 +714,86 @@ TParticle* AliStack::Particle(Int_t i) // // Return particle with specified ID - if(!(*fParticleMap)[i]) { - Int_t nentries = fParticles->GetEntries(); + //PH if(!(*fParticleMap)[i]) { + 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); + TreeK()->GetEntry(entry); new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); fParticleMap->AddAt((*fParticles)[nentries],i); } - return (TParticle *) (*fParticleMap)[i]; + //PH return dynamic_cast((*fParticleMap)[i]); + return dynamic_cast(fParticleMap->At(i)); +} + +//_____________________________________________________________________________ +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 +// + Int_t entry; + if (idAt(fCurrent); + + if (current) + return current->GetFirstMother(); + 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; } } @@ -602,7 +805,8 @@ void AliStack::DumpPart (Int_t i) const // Dumps particle i in the stack // - ((TParticle*) (*fParticleMap)[i])->Print(); + //PH dynamic_cast((*fParticleMap)[i])->Print(); + dynamic_cast(fParticleMap->At(i))->Print(); } //_____________________________________________________________________________ @@ -614,8 +818,7 @@ void AliStack::DumpPStack () Int_t i; - printf( - "\n\n=======================================================================\n"); + printf("\n\n=======================================================================\n"); for (i=0;i(particles[i]); if (particle) { printf("-> %d ",i); particle->Print(); printf("--------------------------------------------------------------\n"); @@ -680,7 +882,7 @@ void AliStack::CleanParents() TParticle *part; int i; for(i=0; i(particles.At(i)); if(part) if(!part->TestBit(kDaughtersBit)) { part->SetFirstDaughter(-1); part->SetLastDaughter(-1); @@ -700,11 +902,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 = dynamic_cast(fParticleMap->At(i)); if ((particle) && (!particle->TestBit(kDoneBit))) { fCurrent=i; - //cout << "GetNextParticle() - secondary " - // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; return particle; } } @@ -712,87 +912,232 @@ TParticle* AliStack::GetNextParticle() // take next primary if all secondaries were done while (fCurrentPrimary>=0) { fCurrent = fCurrentPrimary; - particle = (TParticle*) fParticleMap->At(fCurrentPrimary--); + particle = dynamic_cast(fParticleMap->At(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; } +//__________________________________________________________________________________________ +TTree* AliStack::TreeK() +{ +//returns TreeK + if (fTreeK) + { + return fTreeK; + } + else + { + AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName); + if (rl == 0x0) + { + AliFatal(Form("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 + AliWarning(Form("Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data())); + } + } + return fTreeK;//never reached +} //__________________________________________________________________________________________ -void AliStack::MakeTree(Int_t event, const char *file) + +void AliStack::ConnectTree() { // -// 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); - } +// Creates branch for writing particles +// + 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 + } + + // 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(AliRunLoader::GetKineBranchName()); + if(branch == 0x0) + { + branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &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::BeginEvent(Int_t event) + +void AliStack::BeginEvent() { // start a new event -// -// - fNprimary = 0; - fNtrack = 0; - - char hname[30]; - if(fTreeK) { - fTreeK->Reset(); - sprintf(hname,"TreeK%d",event); - fTreeK->SetName(hname); - } + Reset(); } //_____________________________________________________________________________ void AliStack::FinishRun() { // Clean TreeK information - if (fTreeK) { - delete fTreeK; fTreeK = 0; - } } - //_____________________________________________________________________________ -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); + fTreeK = 0x0; - if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticleBuffer); + if (TreeK() == 0x0) //forces connecting + { + AliError("cannot find Kine Tree for current event"); + return kFALSE; + } + + Int_t size = (Int_t)TreeK()->GetEntries(); + ResetArrays(size); + return kTRUE; +} +//_____________________________________________________________________________ - else - Error("GetEvent","cannot find Kine Tree for event:%d\n",event); +void AliStack::SetEventFolderName(const char* foldname) +{ + //Sets event folder name + fEventFolderName = foldname; +} + +Bool_t AliStack::IsStable(Int_t pdg) const +{ // -// printf("\n primaries %d", fNprimary); -// printf("\n tracks %d", fNtrack); +// Decide whether particle (pdg) is stable // - Int_t size = (Int_t)fTreeK->GetEntries(); - ResetArrays(size); + + const Int_t kNstable = 14; + Int_t i; + + Int_t pdgStable[kNstable] = { + kGamma, // Photon + kElectron, // Electron + kMuonPlus, // Muon + kPiPlus, // Pion + kKPlus, // Kaon + kProton, // Proton + kNeutron, // Neutron + kLambda0, // Lambda_0 + kSigmaMinus, // Sigma Minus + kSigma0, // Sigma_0 + 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 > 20) 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 +// +// Check if this is a heavy flavor decay product + Int_t imo = p->GetFirstMother(); + TParticle* pm = Particle(imo); + Int_t mpdg = TMath::Abs(pm->GetPdgCode()); + 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 = p->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 ? +} +