X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStack.cxx;h=a0405e41eda543558330718e75c8f3a70828c7cb;hb=817e1004756df1b947a1f02e3a3e907336274186;hp=049b29037f2d5cf384f35ef606c1c4093be5129b;hpb=f3069e25a66d07cec3702f384e979d4a5d6a0469;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliStack.cxx b/STEER/AliStack.cxx index 049b29037f2..a0405e41eda 100644 --- a/STEER/AliStack.cxx +++ b/STEER/AliStack.cxx @@ -29,8 +29,10 @@ #include #include #include +#include #include #include +#include #include #include @@ -123,8 +125,8 @@ AliStack::~AliStack() // //_____________________________________________________________________________ -void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, - Float_t *vpos, Float_t *polar, Float_t tof, +void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom, + const Float_t *vpos, const Float_t *polar, Float_t tof, TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) { // @@ -206,7 +208,7 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, = 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); @@ -225,11 +227,10 @@ void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, particle->SetBit(kDaughtersBit); // Add the particle to the stack - fParticleMap.AddAtAndExpand(particle, fNtrack);//CHECK!! if(parent>=0) { - particle = dynamic_cast(fParticleMap.At(parent)); + particle = GetParticleMapEntry(parent); if (particle) { particle->SetLastDaughter(fNtrack); if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); @@ -280,8 +281,10 @@ TParticle* AliStack::PopPrimaryForTracking(Int_t i) TParticle* particle = Particle(i); - if (!particle->TestBit(kDoneBit)) + if (!particle->TestBit(kDoneBit)) { + fCurrentTrack = particle; return particle; + } else return 0; } @@ -308,7 +311,7 @@ Bool_t AliStack::PurifyKine() if(i<=fHgwmk) fTrackLabelMap[i]=i ; else { fTrackLabelMap[i] = -99; - if((part=dynamic_cast(fParticleMap.At(i)))) { + if((part=GetParticleMapEntry(i))) { // // Check of this track should be kept for physics reasons if (KeepPhysics(part)) KeepTrack(i); @@ -322,7 +325,7 @@ Bool_t 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 = dynamic_cast(fParticleMap.At(fHgwmk+1)); + part = GetParticleMapEntry(fHgwmk+1); fParticleMap.At(part->GetFirstMother())->ResetBit(kDaughtersBit); // Second pass, build map between old and new numbering for(i=fHgwmk+1; i(fParticleMap.At(nkeep)); + 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((parent=part->GetFirstMother())>fHgwmk) { if(fTrackLabelMap[parent]==-99) Fatal("PurifyKine","fTrackLabelMap[%d] = -99!\n",parent); - else part->SetFirstMother(fTrackLabelMap[parent]); + else part->SetFirstMother(fTrackLabelMap[parent]);} nkeep++; } } // Fix daughters information for (i=fHgwmk+1; i(fParticleMap.At(i)); + part = GetParticleMapEntry(i); parent = part->GetFirstMother(); if(parent>=0) { - father = dynamic_cast(fParticleMap.At(parent)); + father = GetParticleMapEntry(parent); if(father->TestBit(kDaughtersBit)) { if(iGetFirstDaughter()) father->SetFirstDaughter(i); @@ -363,7 +366,7 @@ Bool_t AliStack::PurifyKine() // 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(fParticleMap.At(i)); + fParticleBuffer = GetParticleMapEntry(i); fParticleFileMap[i]=static_cast(TreeK()->GetEntries()); TreeK()->Fill(); fParticleMap[i]=fParticleBuffer=0; @@ -403,7 +406,7 @@ Bool_t AliStack::ReorderKine() for (i = 0; i < nNew; i++) { if (fParticleMap.At(fHgwmk + 1 + i)) { - tmp[i] = (TParticle*) (fParticleMap.At(fHgwmk + 1 + i)); + tmp[i] = GetParticleMapEntry(fHgwmk + 1 + i); } else { tmp[i] = 0x0; } @@ -425,7 +428,7 @@ Bool_t AliStack::ReorderKine() TParticle* parP; if (i == -1) { ipa = tmp[0]->GetFirstMother(); - parP =dynamic_cast(fParticleMap.At(ipa)); + parP = GetParticleMapEntry(ipa); } else { ipa = (fHgwmk + 1 + i); // Skip deleted particles @@ -482,30 +485,35 @@ Bool_t AliStack::ReorderKine() return kTRUE; } -Bool_t AliStack::KeepPhysics(TParticle* part) +Bool_t AliStack::KeepPhysics(const TParticle* part) { // // Some particles have to kept on the stack for reasons motivated // by physics analysis. Decision is put here. // Bool_t keep = kFALSE; + + Int_t parent = part->GetFirstMother(); + if (parent >= 0 && parent <= fHgwmk) { + TParticle* father = GetParticleMapEntry(parent); // // Keep first-generation daughter from primaries with heavy flavor // - Int_t parent = part->GetFirstMother(); - if (parent >= 0 && parent <= fHgwmk) { - TParticle* father = dynamic_cast(Particles()->At(parent)); Int_t kf = father->GetPdgCode(); kf = TMath::Abs(kf); Int_t kfl = kf; // meson ? if (kfl > 10) kfl/=100; // baryon - if (kfl > 10) kfl/=10; - if (kfl > 10) kfl/=10; + if (kfl > 10) kfl/=10; + if (kfl > 10) kfl/=10; if (kfl >= 4) { keep = kTRUE; } + // + // e+e- from pair production of primary gammas + // + if ((part->GetUniqueID()) == kPPair) keep = kTRUE; } return keep; } @@ -532,10 +540,10 @@ void AliStack::FinishEvent() } Bool_t allFilled = kFALSE; - TObject *part; - for(Int_t i=0; i(part); + TParticle *part; + for(Int_t i=0; i(TreeK()->GetEntries()); TreeK()->Fill(); fParticleBuffer=0; @@ -546,14 +554,15 @@ void AliStack::FinishEvent() if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i)); } else - { + { // // printf("Why = 0 part # %d?\n",i); => We know. // break; // we don't break now in order to be sure there is no // particle !=0 left. // To be removed later and replaced with break. if(!allFilled) allFilled = kTRUE; - } + } + } } //_____________________________________________________________________________ @@ -567,7 +576,7 @@ void AliStack::FlagTrack(Int_t track) Int_t curr=track; while(1) { - particle=dynamic_cast(fParticleMap.At(curr)); + particle = GetParticleMapEntry(curr); // If the particle is flagged the three from here upward is saved already if(particle->TestBit(kKeepBit)) return; @@ -667,7 +676,7 @@ TParticle* AliStack::Particle(Int_t i) // Store a pointer in the TObjArray fParticleMap.AddAt(fParticles[nentries],i); } - return dynamic_cast(fParticleMap.At(i)); + return GetParticleMapEntry(i); } //_____________________________________________________________________________ @@ -713,7 +722,7 @@ Int_t AliStack::GetCurrentParentTrackNumber() const // Return number of the parent of the current track // - TParticle* current = (TParticle*)fParticleMap.At(fCurrent); + TParticle* current = GetParticleMapEntry(fCurrent); if (current) return current->GetFirstMother(); @@ -746,7 +755,7 @@ void AliStack::DumpPart (Int_t i) const // // Dumps particle i in the stack // - dynamic_cast(fParticleMap.At(i))->Print(); + GetParticleMapEntry(i)->Print(); } //_____________________________________________________________________________ @@ -773,9 +782,9 @@ void AliStack::DumpPStack () printf("\n=======================================================================\n\n"); // print particle file map - printf("\nParticle file map: \n"); - for (i=0; i(fParticleMap[i]); + TParticle* particle = GetParticleMapEntry(i); if (particle) { printf("-> %d ",i); particle->Print(); printf("--------------------------------------------------------------\n"); @@ -805,6 +814,15 @@ void AliStack::DumpLoadedStack() const "\n=======================================================================\n\n"); } +//_____________________________________________________________________________ +void AliStack::SetCurrentTrack(Int_t track) +{ + fCurrent = track; + if (fCurrent < fNprimary) fCurrentTrack = Particle(track); +} + + +//_____________________________________________________________________________ // // protected methods // @@ -820,7 +838,7 @@ void AliStack::CleanParents() TParticle *part; int i; for(i=0; i(fParticleMap.At(i)); + part = GetParticleMapEntry(i); if(part) if(!part->TestBit(kDaughtersBit)) { part->SetFirstDaughter(-1); part->SetLastDaughter(-1); @@ -840,7 +858,7 @@ TParticle* AliStack::GetNextParticle() // search secondaries //for(Int_t i=fNtrack-1; i>=0; i--) { for(Int_t i=fNtrack-1; i>fHgwmk; i--) { - particle = dynamic_cast(fParticleMap.At(i)); + particle = GetParticleMapEntry(i); if ((particle) && (!particle->TestBit(kDoneBit))) { fCurrent=i; return particle; @@ -850,7 +868,7 @@ TParticle* AliStack::GetNextParticle() // take next primary if all secondaries were done while (fCurrentPrimary>=0) { fCurrent = fCurrentPrimary; - particle = dynamic_cast(fParticleMap.At(fCurrentPrimary--)); + particle = GetParticleMapEntry(fCurrentPrimary--); if ((particle) && (!particle->TestBit(kDoneBit))) { return particle; } @@ -864,13 +882,6 @@ TParticle* AliStack::GetNextParticle() } //__________________________________________________________________________________________ -TTree* AliStack::TreeK() -{ -//returns TreeK - return fTreeK; -} -//__________________________________________________________________________________________ - void AliStack::ConnectTree(TTree* tree) { // @@ -947,7 +958,7 @@ Bool_t AliStack::IsStable(Int_t pdg) const // Nuclear code is 10LZZZAAAI if(pdg>1000000000)return kTRUE; - const Int_t kNstable = 14; + const Int_t kNstable = 18; Int_t i; Int_t pdgStable[kNstable] = { @@ -956,15 +967,19 @@ Bool_t AliStack::IsStable(Int_t pdg) const kMuonPlus, // Muon kPiPlus, // Pion kKPlus, // Kaon + kK0Short, // K0s + kK0Long, // K0l kProton, // Proton kNeutron, // Neutron kLambda0, // Lambda_0 kSigmaMinus, // Sigma Minus - kSigma0, // Sigma_0 kSigmaPlus, // Sigma Plus 3312, // Xsi Minus 3322, // Xsi 3334, // Omega + kNuE, // Electron Neutrino + kNuMu, // Muon Neutrino + kNuTau // Tau Neutrino }; Bool_t isStable = kFALSE; @@ -978,6 +993,7 @@ Bool_t AliStack::IsStable(Int_t pdg) const return isStable; } +//_____________________________________________________________________________ Bool_t AliStack::IsPhysicalPrimary(Int_t index) { // @@ -991,7 +1007,7 @@ Bool_t AliStack::IsPhysicalPrimary(Int_t index) // // Initial state particle - if (ist > 20) return kFALSE; + if (ist > 1) return kFALSE; Int_t pdg = TMath::Abs(p->GetPdgCode());