X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=d6e4f9a0b0a317b9cc770bdcc51a5c4a2214e02b;hb=86121c917e091118d633145c4b37e687c7628cc6;hp=3a5449b47eb4444616b64eb93940f73a1ab14d84;hpb=e94530da2f51a4f0384cbdce6cfaf34e87c85c54;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index 3a5449b47eb..d6e4f9a0b0a 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -13,142 +13,102 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.24 2002/10/22 15:02:15 alibrary -Introducing Riostream.h - -Revision 1.23 2002/10/14 14:57:32 hristov -Merging the VirtualMC branch to the main development branch (HEAD) - -Revision 1.19.4.2 2002/08/28 15:06:52 alibrary -Updating to v3-09-01 - -Revision 1.19.4.1 2002/06/10 14:43:06 hristov -Merged with v3-08-02 - -Revision 1.21 2002/05/28 14:24:57 hristov -Correct warning messages - -Revision 1.20 2002/04/30 11:47:30 morsch -KeepPhysics method called by PurifyKine added (N. Carrer, A.M.) - -Revision 1.19 2002/03/12 11:06:03 morsch -Add particle status code to argument list of SetTrack(..). - -Revision 1.18 2002/02/20 16:14:41 hristov -fParticleBuffer points to object which doesn't belong to AliStack, do not delete it (J.Chudoba) - -Revision 1.17 2001/12/05 08:51:56 hristov -The default constructor now creates no objects (thanks to r.Brun). Some corrections required by the previous changes. - -Revision 1.16 2001/11/20 09:27:55 hristov -Possibility to investigate a primary of not yet loaded particle (I.Hrivnacova) - -Revision 1.15 2001/09/04 15:10:37 hristov -Additional protection is included to avoid some problems using Hijing - -Revision 1.14 2001/08/30 09:44:06 hristov -VertexSource_t added to avoid the warnings - -Revision 1.13 2001/08/29 13:31:42 morsch -Protection against (fTreeK == 0) in destructor. - -Revision 1.12 2001/07/27 13:03:13 hristov -Default Branch split level set to 99 - -Revision 1.11 2001/07/27 12:34:20 jchudoba -remove the dummy argument in GetEvent method - -Revision 1.10 2001/07/20 10:13:54 morsch -In Particle(Int_t) use GetEntriesFast to speed up the procedure. - -Revision 1.9 2001/07/03 08:10:57 hristov -J.Chudoba's changes merged correctly with the HEAD - -Revision 1.6 2001/05/31 06:59:06 fca -Clean setting and deleting of fParticleBuffer - -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 +#include "AliLog.h" #include "AliStack.h" -#include "AliRun.h" -#include "AliModule.h" -#include "AliHit.h" ClassImp(AliStack) -//_____________________________________________________________________________ -AliStack::AliStack(Int_t size) - : TVirtualMCStack() +//_______________________________________________________________________ +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), + fTrackLabelMap(0) { // - // 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() - : TVirtualMCStack() +//_______________________________________________________________________ +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), + fTrackLabelMap(0) { // - // Default constructor + // Constructor // - - // Create the particles arrays - fParticles = 0; - fParticleMap = 0; - 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), + fTrackLabelMap(0) +{ + // Copy constructor +} -//_____________________________________________________________________________ + +//_______________________________________________________________________ +void AliStack::Copy(TObject&) const +{ + AliFatal("Not implemented!"); +} + +//_______________________________________________________________________ AliStack::~AliStack() { // @@ -160,7 +120,6 @@ AliStack::~AliStack() delete fParticles; } delete fParticleMap; - if (fTreeK) delete fTreeK; } // @@ -168,9 +127,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, Int_t is) +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 @@ -196,25 +155,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 // - 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]); - + 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); - SetTrack(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); + 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, Double_t weight, Int_t is) + TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is) { // // Load a track on the stack @@ -237,10 +202,10 @@ 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++]) @@ -251,37 +216,46 @@ void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, 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!! if(parent>=0) { - particle = (TParticle*) 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 = 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 { + // + // 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 @@ -294,64 +268,15 @@ TParticle* AliStack::GetNextTrack(Int_t& itrack) itrack = fCurrent; track->SetBit(kDoneBit); } - else + else itrack = -1; - - 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(); + fCurrentTrack = track; + return track; } -*/ //_____________________________________________________________________________ -TParticle* AliStack::GetPrimaryForTracking(Int_t i) +TParticle* AliStack::PopPrimaryForTracking(Int_t i) { // // Returns i-th primary particle if it is flagged to be tracked, @@ -366,7 +291,6 @@ TParticle* AliStack::GetPrimaryForTracking(Int_t i) return 0; } - //_____________________________________________________________________________ void AliStack::PurifyKine() { @@ -378,118 +302,191 @@ void AliStack::PurifyKine() TObjArray &particles = *fParticleMap; int nkeep=fHgwmk+1, parent, i; TParticle *part, *father; - TArrayI map(particles.GetLast()+1); + fTrackLabelMap.Set(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); + 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 = (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; 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(particles.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) 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(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 = dynamic_cast(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); + } } - } - } - - // 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()); - 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(particles.At(i)); + fParticleFileMap[i]=static_cast(TreeK()->GetEntries()); + TreeK()->Fill(); + particles[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; iGetEntries(); - fTreeK->Fill(); - particles[i]=fParticleBuffer=0; - } - - for (i=nkeep; iRemoveAt(i); + fNtrack=nkeep; + fHgwmk=nkeep-1; +} - Int_t toshrink = fNtrack-fHgwmk-1; - fLoadPoint-=toshrink; - for(i=fLoadPoint; iRemoveAt(i); +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; - 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; + 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 + // + 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 } Bool_t AliStack::KeepPhysics(TParticle* part) @@ -504,7 +501,7 @@ Bool_t AliStack::KeepPhysics(TParticle* part) // Int_t parent = part->GetFirstMother(); if (parent >= 0 && parent <= fHgwmk) { - TParticle* father = (TParticle*) Particles()->At(parent); + TParticle* father = dynamic_cast(Particles()->At(parent)); Int_t kf = father->GetPdgCode(); kf = TMath::Abs(kf); Int_t kfl = kf; @@ -523,22 +520,20 @@ 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); } @@ -547,30 +542,28 @@ void AliStack::FinishEvent() TObject *part; for(Int_t i=0; iAt(i))) { - fParticleBuffer = (TParticle*) part; - fParticleFileMap[i]= (Int_t) fTreeK->GetEntries(); - fTreeK->Fill(); - //PH (*fParticleMap)[i]=fParticleBuffer=0; + fParticleBuffer = dynamic_cast(part); + fParticleFileMap[i]= static_cast(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) { // @@ -581,7 +574,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; @@ -605,12 +598,12 @@ void AliStack::KeepTrack(Int_t track) } //_____________________________________________________________________________ -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; @@ -619,6 +612,18 @@ void AliStack::Reset(Int_t size) ResetArrays(size); } +//_____________________________________________________________________________ +void AliStack::Reset(Int_t size) +{ + // + // Reset stack data including fTreeK + // + + Clean(size); + + fTreeK = 0x0; +} + //_____________________________________________________________________________ void AliStack::ResetArrays(Int_t size) { @@ -638,18 +643,16 @@ void AliStack::ResetArrays(Int_t 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; } //_____________________________________________________________________________ @@ -657,8 +660,7 @@ TParticle* AliStack::Particle(Int_t i) { // // Return particle with specified ID - - //PH if(!(*fParticleMap)[i]) { + if(!fParticleMap->At(i)) { Int_t nentries = fParticles->GetEntriesFast(); // algorithmic way of getting entry index @@ -668,17 +670,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); + // 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 (TParticle *) (*fParticleMap)[i]; - return (TParticle *) fParticleMap->At(i); + return dynamic_cast(fParticleMap->At(i)); } //_____________________________________________________________________________ @@ -697,9 +700,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); + + if (current) + return current->GetFirstMother(); + else { + AliWarning("Current track not found in the stack"); + return -1; + } +} + //_____________________________________________________________________________ Int_t AliStack::GetPrimary(Int_t id) { @@ -731,9 +760,7 @@ void AliStack::DumpPart (Int_t i) const // // Dumps particle i in the stack // - - //PH ((TParticle*) (*fParticleMap)[i])->Print(); - ((TParticle*) fParticleMap->At(i))->Print(); + dynamic_cast(fParticleMap->At(i))->Print(); } //_____________________________________________________________________________ @@ -745,8 +772,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"); @@ -811,7 +836,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); @@ -831,11 +856,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; } } @@ -843,89 +866,185 @@ 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; } - //__________________________________________________________________________________________ -void AliStack::MakeTree(Int_t event, const char *file) + +TTree* AliStack::TreeK() { -// -// 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); - } +//returns TreeK + return fTreeK; } +//__________________________________________________________________________________________ -//_____________________________________________________________________________ -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 + } + + // 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", "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::FinishRun() -{ -// Clean TreeK information - if (fTreeK) { - delete fTreeK; fTreeK = 0; - } -} -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 = (TTree*)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 = (Int_t)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 +// + + 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 ? +} +