]>
Commit | Line | Data |
---|---|---|
9e1a0ddb | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
acd84897 | 16 | /* $Id$ */ |
9e1a0ddb | 17 | |
18 | /////////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
af7ba10c | 20 | // Particles stack class // |
21 | // Implements the TMCVirtualStack of the Virtual Monte Carlo // | |
22 | // Holds the particles transported during simulation // | |
23 | // Is used to compare results of reconstruction with simulation // | |
24 | // Author A.Morsch // | |
9e1a0ddb | 25 | // // |
26 | /////////////////////////////////////////////////////////////////////////////// | |
27 | ||
9e1a0ddb | 28 | |
7ca4655f | 29 | #include <TClonesArray.h> |
9e1a0ddb | 30 | #include <TObjArray.h> |
7ca4655f | 31 | #include <TPDGCode.h> |
9e1a0ddb | 32 | #include <TParticle.h> |
fe046ade | 33 | #include <TParticlePDG.h> |
9e1a0ddb | 34 | #include <TTree.h> |
9e1a0ddb | 35 | |
21bf7095 | 36 | #include "AliLog.h" |
9e1a0ddb | 37 | #include "AliHit.h" |
88cb7938 | 38 | #include "AliModule.h" |
39 | #include "AliRun.h" | |
5d12ce38 | 40 | #include "AliMC.h" |
88cb7938 | 41 | #include "AliRunLoader.h" |
42 | #include "AliStack.h" | |
9e1a0ddb | 43 | |
44 | ClassImp(AliStack) | |
45 | ||
e2afb3b6 | 46 | //_______________________________________________________________________ |
47 | AliStack::AliStack(): | |
48 | fParticles(0), | |
49 | fParticleMap(0), | |
50 | fParticleFileMap(0), | |
51 | fParticleBuffer(0), | |
90e48c0c | 52 | fCurrentTrack(0), |
e2afb3b6 | 53 | fTreeK(0), |
54 | fNtrack(0), | |
55 | fNprimary(0), | |
56 | fCurrent(-1), | |
57 | fCurrentPrimary(-1), | |
58 | fHgwmk(0), | |
88cb7938 | 59 | fLoadPoint(0), |
e191bb57 | 60 | fEventFolderName(AliConfig::GetDefaultEventFolderName()) |
9e1a0ddb | 61 | { |
62 | // | |
e2afb3b6 | 63 | // Default constructor |
9e1a0ddb | 64 | // |
9e1a0ddb | 65 | } |
66 | ||
e2afb3b6 | 67 | //_______________________________________________________________________ |
88cb7938 | 68 | AliStack::AliStack(Int_t size, const char* evfoldname): |
e2afb3b6 | 69 | fParticles(new TClonesArray("TParticle",1000)), |
70 | fParticleMap(new TObjArray(size)), | |
71 | fParticleFileMap(0), | |
72 | fParticleBuffer(0), | |
90e48c0c | 73 | fCurrentTrack(0), |
e2afb3b6 | 74 | fTreeK(0), |
75 | fNtrack(0), | |
76 | fNprimary(0), | |
77 | fCurrent(-1), | |
78 | fCurrentPrimary(-1), | |
79 | fHgwmk(0), | |
88cb7938 | 80 | fLoadPoint(0), |
81 | fEventFolderName(evfoldname) | |
e2afb3b6 | 82 | { |
83 | // | |
84 | // Constructor | |
85 | // | |
86 | } | |
9e1a0ddb | 87 | |
e2afb3b6 | 88 | //_______________________________________________________________________ |
89 | AliStack::AliStack(const AliStack& st): | |
90 | TVirtualMCStack(st), | |
91 | fParticles(0), | |
92 | fParticleMap(0), | |
93 | fParticleFileMap(0), | |
94 | fParticleBuffer(0), | |
90e48c0c | 95 | fCurrentTrack(0), |
e2afb3b6 | 96 | fTreeK(0), |
97 | fNtrack(0), | |
98 | fNprimary(0), | |
99 | fCurrent(-1), | |
100 | fCurrentPrimary(-1), | |
101 | fHgwmk(0), | |
90e48c0c | 102 | fLoadPoint(0), |
103 | fEventFolderName() | |
9e1a0ddb | 104 | { |
105 | // | |
e2afb3b6 | 106 | // Copy constructor |
9e1a0ddb | 107 | // |
e2afb3b6 | 108 | st.Copy(*this); |
9e1a0ddb | 109 | } |
110 | ||
e2afb3b6 | 111 | //_______________________________________________________________________ |
6c4904c2 | 112 | void AliStack::Copy(TObject&) const |
e2afb3b6 | 113 | { |
21bf7095 | 114 | AliFatal("Not implemented!"); |
e2afb3b6 | 115 | } |
9e1a0ddb | 116 | |
e2afb3b6 | 117 | //_______________________________________________________________________ |
9e1a0ddb | 118 | AliStack::~AliStack() |
119 | { | |
120 | // | |
121 | // Destructor | |
122 | // | |
123 | ||
124 | if (fParticles) { | |
125 | fParticles->Delete(); | |
126 | delete fParticles; | |
127 | } | |
128 | delete fParticleMap; | |
88cb7938 | 129 | //PH??? delete fTreeK; |
9e1a0ddb | 130 | } |
131 | ||
132 | // | |
133 | // public methods | |
134 | // | |
135 | ||
136 | //_____________________________________________________________________________ | |
642f15cf | 137 | void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, |
88cb7938 | 138 | Float_t *vpos, Float_t *polar, Float_t tof, |
139 | TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) | |
9e1a0ddb | 140 | { |
141 | // | |
142 | // Load a track on the stack | |
143 | // | |
144 | // done 0 if the track has to be transported | |
145 | // 1 if not | |
146 | // parent identifier of the parent track. -1 for a primary | |
147 | // pdg particle code | |
148 | // pmom momentum GeV/c | |
149 | // vpos position | |
150 | // polar polarisation | |
151 | // tof time of flight in seconds | |
152 | // mecha production mechanism | |
153 | // ntr on output the number of the track stored | |
154 | // | |
155 | ||
9e1a0ddb | 156 | // const Float_t tlife=0; |
157 | ||
158 | // | |
159 | // Here we get the static mass | |
160 | // For MC is ok, but a more sophisticated method could be necessary | |
161 | // if the calculated mass is required | |
162 | // also, this method is potentially dangerous if the mass | |
163 | // used in the MC is not the same of the PDG database | |
164 | // | |
fe046ade | 165 | TParticlePDG* pmc = TDatabasePDG::Instance()->GetParticle(pdg); |
166 | if (pmc) { | |
167 | Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
168 | Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+ | |
169 | pmom[1]*pmom[1]+pmom[2]*pmom[2]); | |
170 | ||
ccf7a81f | 171 | // 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", |
172 | // mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); | |
173 | ||
9e1a0ddb | 174 | |
642f15cf | 175 | PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e, |
fe046ade | 176 | vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2], |
177 | mech, ntr, weight, is); | |
178 | } else { | |
21bf7095 | 179 | AliWarning(Form("Particle type %d not defined in PDG Database !", pdg)); |
180 | AliWarning("Particle skipped !"); | |
fe046ade | 181 | } |
9e1a0ddb | 182 | } |
183 | ||
184 | //_____________________________________________________________________________ | |
642f15cf | 185 | void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, |
9e1a0ddb | 186 | Double_t px, Double_t py, Double_t pz, Double_t e, |
187 | Double_t vx, Double_t vy, Double_t vz, Double_t tof, | |
188 | Double_t polx, Double_t poly, Double_t polz, | |
98490ea9 | 189 | TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is) |
9e1a0ddb | 190 | { |
191 | // | |
192 | // Load a track on the stack | |
193 | // | |
194 | // done 0 if the track has to be transported | |
195 | // 1 if not | |
196 | // parent identifier of the parent track. -1 for a primary | |
197 | // pdg particle code | |
198 | // kS generation status code | |
199 | // px, py, pz momentum GeV/c | |
200 | // vx, vy, vz position | |
201 | // polar polarisation | |
202 | // tof time of flight in seconds | |
203 | // mech production mechanism | |
204 | // ntr on output the number of the track stored | |
205 | // | |
206 | // New method interface: | |
207 | // arguments were changed to be in correspondence with TParticle | |
208 | // constructor. | |
209 | // Note: the energy is not calculated from the static mass but | |
210 | // it is passed by argument e. | |
211 | ||
212 | ||
9e1a0ddb | 213 | const Int_t kFirstDaughter=-1; |
214 | const Int_t kLastDaughter=-1; | |
6308b375 | 215 | |
216 | ||
9e1a0ddb | 217 | TClonesArray &particles = *fParticles; |
218 | TParticle* particle | |
219 | = new(particles[fLoadPoint++]) | |
47c8bcbe | 220 | TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter, |
9e1a0ddb | 221 | px, py, pz, e, vx, vy, vz, tof); |
222 | ||
223 | particle->SetPolarisation(polx, poly, polz); | |
224 | particle->SetWeight(weight); | |
225 | particle->SetUniqueID(mech); | |
226 | ||
227 | if(!done) particle->SetBit(kDoneBit); | |
228 | ||
229 | // Declare that the daughter information is valid | |
230 | particle->SetBit(kDaughtersBit); | |
231 | // Add the particle to the stack | |
6308b375 | 232 | |
233 | ||
9e1a0ddb | 234 | fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!! |
235 | ||
236 | if(parent>=0) { | |
e2afb3b6 | 237 | particle = dynamic_cast<TParticle*>(fParticleMap->At(parent)); |
1305ca2f | 238 | if (particle) { |
239 | particle->SetLastDaughter(fNtrack); | |
240 | if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); | |
241 | } | |
242 | else { | |
21bf7095 | 243 | AliError(Form("Parent %d does not exist",parent)); |
1305ca2f | 244 | } |
9e1a0ddb | 245 | } |
246 | else { | |
247 | // | |
248 | // This is a primary track. Set high water mark for this event | |
249 | fHgwmk = fNtrack; | |
250 | // | |
251 | // Set also number if primary tracks | |
252 | fNprimary = fHgwmk+1; | |
253 | fCurrentPrimary++; | |
254 | } | |
255 | ntr = fNtrack++; | |
256 | } | |
257 | ||
258 | //_____________________________________________________________________________ | |
642f15cf | 259 | TParticle* AliStack::PopNextTrack(Int_t& itrack) |
b9d0a01d | 260 | { |
261 | // | |
262 | // Returns next track from stack of particles | |
263 | // | |
264 | ||
265 | ||
266 | TParticle* track = GetNextParticle(); | |
267 | ||
268 | if (track) { | |
269 | itrack = fCurrent; | |
270 | track->SetBit(kDoneBit); | |
271 | } | |
fe046ade | 272 | else |
b9d0a01d | 273 | itrack = -1; |
fe046ade | 274 | |
275 | fCurrentTrack = track; | |
b9d0a01d | 276 | return track; |
277 | } | |
278 | ||
b9d0a01d | 279 | //_____________________________________________________________________________ |
642f15cf | 280 | TParticle* AliStack::PopPrimaryForTracking(Int_t i) |
b9d0a01d | 281 | { |
282 | // | |
283 | // Returns i-th primary particle if it is flagged to be tracked, | |
284 | // 0 otherwise | |
285 | // | |
286 | ||
287 | TParticle* particle = Particle(i); | |
288 | ||
289 | if (!particle->TestBit(kDoneBit)) | |
290 | return particle; | |
291 | else | |
292 | return 0; | |
293 | } | |
294 | ||
9e1a0ddb | 295 | //_____________________________________________________________________________ |
296 | void AliStack::PurifyKine() | |
297 | { | |
298 | // | |
299 | // Compress kinematic tree keeping only flagged particles | |
300 | // and renaming the particle id's in all the hits | |
301 | // | |
302 | ||
303 | TObjArray &particles = *fParticleMap; | |
304 | int nkeep=fHgwmk+1, parent, i; | |
305 | TParticle *part, *father; | |
306 | TArrayI map(particles.GetLast()+1); | |
307 | ||
308 | // Save in Header total number of tracks before compression | |
9e1a0ddb | 309 | // If no tracks generated return now |
310 | if(fHgwmk+1 == fNtrack) return; | |
311 | ||
9e1a0ddb | 312 | // First pass, invalid Daughter information |
1fed0e8b | 313 | for(i=0; i<fNtrack; i++) { |
9e1a0ddb | 314 | // Preset map, to be removed later |
315 | if(i<=fHgwmk) map[i]=i ; | |
316 | else { | |
317 | map[i] = -99; | |
e2afb3b6 | 318 | if((part=dynamic_cast<TParticle*>(particles.At(i)))) { |
e44c5340 | 319 | // |
320 | // Check of this track should be kept for physics reasons | |
88cb7938 | 321 | if (KeepPhysics(part)) KeepTrack(i); |
e44c5340 | 322 | // |
76969085 | 323 | part->ResetBit(kDaughtersBit); |
324 | part->SetFirstDaughter(-1); | |
325 | part->SetLastDaughter(-1); | |
326 | } | |
9e1a0ddb | 327 | } |
328 | } | |
329 | // Invalid daughter information for the parent of the first particle | |
330 | // generated. This may or may not be the current primary according to | |
331 | // whether decays have been recorded among the primaries | |
e2afb3b6 | 332 | part = dynamic_cast<TParticle*>(particles.At(fHgwmk+1)); |
9e1a0ddb | 333 | particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit); |
334 | // Second pass, build map between old and new numbering | |
335 | for(i=fHgwmk+1; i<fNtrack; i++) { | |
336 | if(particles.At(i)->TestBit(kKeepBit)) { | |
337 | ||
338 | // This particle has to be kept | |
339 | map[i]=nkeep; | |
340 | // If old and new are different, have to move the pointer | |
341 | if(i!=nkeep) particles[nkeep]=particles.At(i); | |
e2afb3b6 | 342 | part = dynamic_cast<TParticle*>(particles.At(nkeep)); |
9e1a0ddb | 343 | |
344 | // as the parent is always *before*, it must be already | |
345 | // in place. This is what we are checking anyway! | |
346 | if((parent=part->GetFirstMother())>fHgwmk) | |
347 | if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent); | |
348 | else part->SetFirstMother(map[parent]); | |
349 | ||
350 | nkeep++; | |
351 | } | |
352 | } | |
353 | ||
354 | // Fix daughters information | |
355 | for (i=fHgwmk+1; i<nkeep; i++) { | |
e2afb3b6 | 356 | part = dynamic_cast<TParticle*>(particles.At(i)); |
9e1a0ddb | 357 | parent = part->GetFirstMother(); |
358 | if(parent>=0) { | |
e2afb3b6 | 359 | father = dynamic_cast<TParticle*>(particles.At(parent)); |
9e1a0ddb | 360 | if(father->TestBit(kDaughtersBit)) { |
361 | ||
362 | if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i); | |
363 | if(i>father->GetLastDaughter()) father->SetLastDaughter(i); | |
364 | } else { | |
365 | // Initialise daughters info for first pass | |
366 | father->SetFirstDaughter(i); | |
367 | father->SetLastDaughter(i); | |
368 | father->SetBit(kDaughtersBit); | |
369 | } | |
370 | } | |
371 | } | |
a1d92a3f | 372 | |
9e1a0ddb | 373 | |
374 | // Now loop on all registered hit lists | |
5d12ce38 | 375 | TList* hitLists = gAlice->GetMCApp()->GetHitLists(); |
9e1a0ddb | 376 | TIter next(hitLists); |
377 | TCollection *hitList; | |
6308b375 | 378 | |
e2afb3b6 | 379 | while((hitList = dynamic_cast<TCollection*>(next()))) { |
9e1a0ddb | 380 | TIter nexthit(hitList); |
381 | AliHit *hit; | |
1fed0e8b | 382 | |
e2afb3b6 | 383 | while((hit = dynamic_cast<AliHit*>(nexthit()))) { |
6308b375 | 384 | hit->SetTrack(map[hit->GetTrack()]); |
9e1a0ddb | 385 | } |
386 | } | |
387 | ||
388 | // | |
389 | // This for detectors which have a special mapping mechanism | |
390 | // for hits, such as TPC and TRD | |
391 | // | |
2b22f272 | 392 | |
9e1a0ddb | 393 | TObjArray* modules = gAlice->Modules(); |
394 | TIter nextmod(modules); | |
395 | AliModule *detector; | |
e2afb3b6 | 396 | while((detector = dynamic_cast<AliModule*>(nextmod()))) { |
9e1a0ddb | 397 | detector->RemapTrackHitIDs(map.GetArray()); |
ae9676e4 | 398 | detector->RemapTrackReferencesIDs(map.GetArray()); |
9e1a0ddb | 399 | } |
2b22f272 | 400 | // |
5d12ce38 | 401 | gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray()); |
2b22f272 | 402 | |
9e1a0ddb | 403 | // Now the output bit, from fHgwmk to nkeep we write everything and we erase |
404 | if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5)); | |
405 | ||
406 | for (i=fHgwmk+1; i<nkeep; ++i) { | |
e2afb3b6 | 407 | fParticleBuffer = dynamic_cast<TParticle*>(particles.At(i)); |
88cb7938 | 408 | fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries()); |
409 | TreeK()->Fill(); | |
17f4a3d9 | 410 | particles[i]=fParticleBuffer=0; |
88cb7938 | 411 | } |
9e1a0ddb | 412 | |
413 | for (i=nkeep; i<fNtrack; ++i) particles[i]=0; | |
414 | ||
76969085 | 415 | Int_t toshrink = fNtrack-fHgwmk-1; |
9e1a0ddb | 416 | fLoadPoint-=toshrink; |
6308b375 | 417 | |
418 | ||
9e1a0ddb | 419 | for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i); |
420 | ||
421 | fNtrack=nkeep; | |
422 | fHgwmk=nkeep-1; | |
6308b375 | 423 | } |
424 | ||
425 | void AliStack::ReorderKine() | |
426 | { | |
427 | // | |
428 | // In some transport code children might not come in a continuous sequence. | |
429 | // In this case the stack has to be reordered in order to establish the | |
430 | // mother daughter relation using index ranges. | |
431 | // | |
432 | if(fHgwmk+1 == fNtrack) return; | |
433 | ||
434 | // | |
435 | // Howmany secondaries have been produced ? | |
436 | Int_t nNew = fNtrack - fHgwmk - 1; | |
a1d92a3f | 437 | |
6308b375 | 438 | if (nNew > 0) { |
439 | Int_t i, j; | |
440 | TObjArray &particles = *fParticleMap; | |
441 | TArrayI map1(nNew); | |
442 | // | |
443 | // Copy pointers to temporary array | |
444 | TParticle** tmp = new TParticle*[nNew]; | |
445 | ||
446 | for (i = 0; i < nNew; i++) { | |
447 | if (particles.At(fHgwmk + 1 + i)) { | |
6308b375 | 448 | tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i)); |
6308b375 | 449 | } else { |
450 | tmp[i] = 0x0; | |
451 | } | |
452 | map1[i] = -99; | |
453 | } | |
454 | ||
455 | ||
456 | // | |
457 | // Reset LoadPoint | |
458 | // | |
459 | fLoadPoint = fHgwmk + 1; | |
460 | // | |
461 | // Re-Push particles into stack | |
462 | // The outer loop is over parents, the inner over children. | |
463 | // -1 refers to the primary particle | |
464 | // | |
a1d92a3f | 465 | for (i = -1; i < nNew-1; i++) { |
6308b375 | 466 | Int_t ipa; |
467 | TParticle* parP; | |
468 | if (i == -1) { | |
469 | ipa = tmp[0]->GetFirstMother(); | |
470 | parP =dynamic_cast<TParticle*>(particles.At(ipa)); | |
471 | } else { | |
472 | ipa = (fHgwmk + 1 + i); | |
473 | // Skip deleted particles | |
474 | if (!tmp[i]) continue; | |
475 | // Skip particles without children | |
476 | if (tmp[i]->GetFirstDaughter() == -1) continue; | |
477 | parP = tmp[i]; | |
478 | } | |
479 | // Reset daughter information | |
a1d92a3f | 480 | |
481 | Int_t idaumin = parP->GetFirstDaughter() - fHgwmk - 1; | |
482 | Int_t idaumax = parP->GetLastDaughter() - fHgwmk - 1; | |
6308b375 | 483 | parP->SetFirstDaughter(-1); |
484 | parP->SetLastDaughter(-1); | |
a1d92a3f | 485 | for (j = idaumin; j <= idaumax; j++) { |
6308b375 | 486 | // Skip deleted particles |
487 | if (!tmp[j]) continue; | |
488 | // Skip particles already handled | |
489 | if (map1[j] != -99) continue; | |
490 | Int_t jpa = tmp[j]->GetFirstMother(); | |
491 | // Check if daughter of current parent | |
492 | if (jpa == ipa) { | |
493 | particles[fLoadPoint] = tmp[j]; | |
494 | // Re-establish daughter information | |
495 | parP->SetLastDaughter(fLoadPoint); | |
496 | if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(fLoadPoint); | |
497 | // Set Mother information | |
498 | if (i != -1) { | |
499 | tmp[j]->SetFirstMother(map1[i]); | |
500 | } | |
501 | // Build the map | |
502 | map1[j] = fLoadPoint; | |
503 | // Increase load point | |
504 | fLoadPoint++; | |
505 | } | |
506 | } // children | |
507 | } // parents | |
508 | ||
509 | delete[] tmp; | |
510 | ||
511 | // | |
512 | // Build map for remapping of hits | |
513 | // | |
514 | TArrayI map(fNtrack); | |
515 | for (i = 0; i < fNtrack; i ++) { | |
516 | if (i <= fHgwmk) { | |
517 | map[i] = i; | |
518 | } else{ | |
519 | map[i] = map1[i - fHgwmk -1]; | |
520 | } | |
521 | } | |
522 | ||
523 | // Now loop on all registered hit lists | |
524 | ||
525 | TList* hitLists = gAlice->GetMCApp()->GetHitLists(); | |
526 | TIter next(hitLists); | |
527 | TCollection *hitList; | |
528 | ||
529 | while((hitList = dynamic_cast<TCollection*>(next()))) { | |
530 | TIter nexthit(hitList); | |
531 | AliHit *hit; | |
532 | while((hit = dynamic_cast<AliHit*>(nexthit()))) { | |
533 | hit->SetTrack(map[hit->GetTrack()]); | |
534 | } | |
535 | } | |
536 | ||
537 | // | |
538 | // This for detectors which have a special mapping mechanism | |
539 | // for hits, such as TPC and TRD | |
540 | // | |
541 | ||
542 | TObjArray* modules = gAlice->Modules(); | |
543 | TIter nextmod(modules); | |
544 | AliModule *detector; | |
545 | while((detector = dynamic_cast<AliModule*>(nextmod()))) { | |
546 | detector->RemapTrackHitIDs(map.GetArray()); | |
547 | detector->RemapTrackReferencesIDs(map.GetArray()); | |
548 | } | |
549 | gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray()); | |
550 | } // new particles poduced | |
9e1a0ddb | 551 | } |
552 | ||
e44c5340 | 553 | Bool_t AliStack::KeepPhysics(TParticle* part) |
554 | { | |
555 | // | |
556 | // Some particles have to kept on the stack for reasons motivated | |
557 | // by physics analysis. Decision is put here. | |
558 | // | |
559 | Bool_t keep = kFALSE; | |
560 | // | |
561 | // Keep first-generation daughter from primaries with heavy flavor | |
562 | // | |
563 | Int_t parent = part->GetFirstMother(); | |
564 | if (parent >= 0 && parent <= fHgwmk) { | |
e2afb3b6 | 565 | TParticle* father = dynamic_cast<TParticle*>(Particles()->At(parent)); |
e44c5340 | 566 | Int_t kf = father->GetPdgCode(); |
567 | kf = TMath::Abs(kf); | |
568 | Int_t kfl = kf; | |
569 | // meson ? | |
570 | if (kfl > 10) kfl/=100; | |
571 | // baryon | |
572 | if (kfl > 10) kfl/=10; | |
573 | if (kfl > 10) kfl/=10; | |
574 | if (kfl >= 4) { | |
575 | keep = kTRUE; | |
576 | } | |
577 | } | |
578 | return keep; | |
579 | } | |
580 | ||
9e1a0ddb | 581 | //_____________________________________________________________________________ |
582 | void AliStack::FinishEvent() | |
583 | { | |
88cb7938 | 584 | // |
585 | // Write out the kinematics that was not yet filled | |
586 | // | |
9e1a0ddb | 587 | |
88cb7938 | 588 | // Update event header |
9e1a0ddb | 589 | |
590 | ||
88cb7938 | 591 | if (!TreeK()) { |
9e1a0ddb | 592 | // Fatal("FinishEvent", "No kinematics tree is defined."); |
593 | // Don't panic this is a probably a lego run | |
594 | return; | |
9e1a0ddb | 595 | } |
596 | ||
597 | CleanParents(); | |
88cb7938 | 598 | if(TreeK()->GetEntries() ==0) { |
9e1a0ddb | 599 | // set the fParticleFileMap size for the first time |
600 | fParticleFileMap.Set(fHgwmk+1); | |
601 | } | |
602 | ||
603 | Bool_t allFilled = kFALSE; | |
604 | TObject *part; | |
605 | for(Int_t i=0; i<fHgwmk+1; ++i) | |
606 | if((part=fParticleMap->At(i))) { | |
e2afb3b6 | 607 | fParticleBuffer = dynamic_cast<TParticle*>(part); |
88cb7938 | 608 | fParticleFileMap[i]= static_cast<Int_t>(TreeK()->GetEntries()); |
609 | TreeK()->Fill(); | |
da2af24c | 610 | //PH (*fParticleMap)[i]=fParticleBuffer=0; |
611 | fParticleBuffer=0; | |
612 | fParticleMap->AddAt(0,i); | |
9e1a0ddb | 613 | |
614 | // When all primaries were filled no particle!=0 | |
615 | // should be left => to be removed later. | |
21bf7095 | 616 | if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i)); |
9e1a0ddb | 617 | } |
88cb7938 | 618 | else |
619 | { | |
9e1a0ddb | 620 | // // printf("Why = 0 part # %d?\n",i); => We know. |
621 | // break; | |
88cb7938 | 622 | // we don't break now in order to be sure there is no |
623 | // particle !=0 left. | |
624 | // To be removed later and replaced with break. | |
625 | if(!allFilled) allFilled = kTRUE; | |
626 | } | |
9e1a0ddb | 627 | } |
9e1a0ddb | 628 | //_____________________________________________________________________________ |
88cb7938 | 629 | |
9e1a0ddb | 630 | void AliStack::FlagTrack(Int_t track) |
631 | { | |
632 | // | |
633 | // Flags a track and all its family tree to be kept | |
634 | // | |
635 | ||
636 | TParticle *particle; | |
637 | ||
638 | Int_t curr=track; | |
639 | while(1) { | |
e2afb3b6 | 640 | particle=dynamic_cast<TParticle*>(fParticleMap->At(curr)); |
9e1a0ddb | 641 | |
642 | // If the particle is flagged the three from here upward is saved already | |
643 | if(particle->TestBit(kKeepBit)) return; | |
644 | ||
645 | // Save this particle | |
646 | particle->SetBit(kKeepBit); | |
647 | ||
648 | // Move to father if any | |
649 | if((curr=particle->GetFirstMother())==-1) return; | |
650 | } | |
651 | } | |
652 | ||
653 | //_____________________________________________________________________________ | |
ed54bf31 | 654 | void AliStack::KeepTrack(Int_t track) |
9e1a0ddb | 655 | { |
656 | // | |
657 | // Flags a track to be kept | |
658 | // | |
659 | ||
660 | fParticleMap->At(track)->SetBit(kKeepBit); | |
661 | } | |
662 | ||
663 | //_____________________________________________________________________________ | |
664 | void AliStack::Reset(Int_t size) | |
665 | { | |
666 | // | |
667 | // Resets stack | |
668 | // | |
88cb7938 | 669 | |
9e1a0ddb | 670 | fNtrack=0; |
671 | fNprimary=0; | |
672 | fHgwmk=0; | |
673 | fLoadPoint=0; | |
674 | fCurrent = -1; | |
88cb7938 | 675 | fTreeK = 0x0; |
9e1a0ddb | 676 | ResetArrays(size); |
677 | } | |
678 | ||
679 | //_____________________________________________________________________________ | |
680 | void AliStack::ResetArrays(Int_t size) | |
681 | { | |
682 | // | |
683 | // Resets stack arrays | |
684 | // | |
685 | ||
be3380a3 | 686 | if (fParticles) |
687 | fParticles->Clear(); | |
688 | else | |
689 | fParticles = new TClonesArray("TParticle",1000); | |
690 | if (fParticleMap) { | |
691 | fParticleMap->Clear(); | |
692 | if (size>0) fParticleMap->Expand(size);} | |
693 | else | |
694 | fParticleMap = new TObjArray(size); | |
9e1a0ddb | 695 | } |
696 | ||
697 | //_____________________________________________________________________________ | |
e2afb3b6 | 698 | void AliStack::SetHighWaterMark(Int_t) |
9e1a0ddb | 699 | { |
700 | // | |
701 | // Set high water mark for last track in event | |
702 | // | |
6308b375 | 703 | |
704 | ||
9e1a0ddb | 705 | fHgwmk = fNtrack-1; |
706 | fCurrentPrimary=fHgwmk; | |
9e1a0ddb | 707 | // Set also number of primary tracks |
708 | fNprimary = fHgwmk+1; | |
709 | fNtrack = fHgwmk+1; | |
710 | } | |
711 | ||
712 | //_____________________________________________________________________________ | |
713 | TParticle* AliStack::Particle(Int_t i) | |
714 | { | |
715 | // | |
716 | // Return particle with specified ID | |
717 | ||
da2af24c | 718 | //PH if(!(*fParticleMap)[i]) { |
719 | if(!fParticleMap->At(i)) { | |
911a3148 | 720 | Int_t nentries = fParticles->GetEntriesFast(); |
9e1a0ddb | 721 | // algorithmic way of getting entry index |
722 | // (primary particles are filled after secondaries) | |
e94530da | 723 | Int_t entry = TreeKEntry(i); |
9e1a0ddb | 724 | // check whether algorithmic way and |
725 | // and the fParticleFileMap[i] give the same; | |
726 | // give the fatal error if not | |
727 | if (entry != fParticleFileMap[i]) { | |
21bf7095 | 728 | AliFatal(Form( |
9e1a0ddb | 729 | "!! The algorithmic way and map are different: !!\n entry: %d map: %d", |
21bf7095 | 730 | entry, fParticleFileMap[i])); |
9e1a0ddb | 731 | } |
732 | ||
88cb7938 | 733 | TreeK()->GetEntry(entry); |
9e1a0ddb | 734 | new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); |
735 | fParticleMap->AddAt((*fParticles)[nentries],i); | |
736 | } | |
e2afb3b6 | 737 | //PH return dynamic_cast<TParticle *>((*fParticleMap)[i]); |
738 | return dynamic_cast<TParticle*>(fParticleMap->At(i)); | |
9e1a0ddb | 739 | } |
740 | ||
e94530da | 741 | //_____________________________________________________________________________ |
742 | TParticle* AliStack::ParticleFromTreeK(Int_t id) const | |
743 | { | |
744 | // | |
745 | // return pointer to TParticle with label id | |
746 | // | |
747 | Int_t entry; | |
748 | if ((entry = TreeKEntry(id)) < 0) return 0; | |
749 | if (fTreeK->GetEntry(entry)<=0) return 0; | |
750 | return fParticleBuffer; | |
751 | } | |
752 | ||
753 | //_____________________________________________________________________________ | |
754 | Int_t AliStack::TreeKEntry(Int_t id) const | |
755 | { | |
756 | // | |
757 | // return entry number in the TreeK for particle with label id | |
758 | // return negative number if label>fNtrack | |
759 | // | |
760 | Int_t entry; | |
761 | if (id<fNprimary) | |
762 | entry = id+fNtrack-fNprimary; | |
763 | else | |
764 | entry = id-fNprimary; | |
765 | return entry; | |
766 | } | |
767 | ||
5d8718b8 | 768 | //_____________________________________________________________________________ |
642f15cf | 769 | Int_t AliStack::GetCurrentParentTrackNumber() const |
5d8718b8 | 770 | { |
771 | // | |
772 | // Return number of the parent of the current track | |
773 | // | |
774 | ||
775 | TParticle* current = (TParticle*)fParticleMap->At(fCurrent); | |
776 | ||
777 | if (current) | |
778 | return current->GetFirstMother(); | |
779 | else { | |
21bf7095 | 780 | AliWarning("Current track not found in the stack"); |
5d8718b8 | 781 | return -1; |
782 | } | |
783 | } | |
784 | ||
9e1a0ddb | 785 | //_____________________________________________________________________________ |
f5f55563 | 786 | Int_t AliStack::GetPrimary(Int_t id) |
9e1a0ddb | 787 | { |
788 | // | |
789 | // Return number of primary that has generated track | |
790 | // | |
791 | ||
792 | int current, parent; | |
9e1a0ddb | 793 | // |
794 | parent=id; | |
795 | while (1) { | |
796 | current=parent; | |
f5f55563 | 797 | parent=Particle(current)->GetFirstMother(); |
9e1a0ddb | 798 | if(parent<0) return current; |
799 | } | |
800 | } | |
801 | ||
802 | //_____________________________________________________________________________ | |
803 | void AliStack::DumpPart (Int_t i) const | |
804 | { | |
805 | // | |
806 | // Dumps particle i in the stack | |
807 | // | |
808 | ||
e2afb3b6 | 809 | //PH dynamic_cast<TParticle*>((*fParticleMap)[i])->Print(); |
810 | dynamic_cast<TParticle*>(fParticleMap->At(i))->Print(); | |
9e1a0ddb | 811 | } |
812 | ||
813 | //_____________________________________________________________________________ | |
814 | void AliStack::DumpPStack () | |
815 | { | |
816 | // | |
817 | // Dumps the particle stack | |
818 | // | |
819 | ||
3ab6f951 | 820 | Int_t i; |
821 | ||
88cb7938 | 822 | printf("\n\n=======================================================================\n"); |
3ab6f951 | 823 | for (i=0;i<fNtrack;i++) |
9e1a0ddb | 824 | { |
825 | TParticle* particle = Particle(i); | |
826 | if (particle) { | |
827 | printf("-> %d ",i); particle->Print(); | |
828 | printf("--------------------------------------------------------------\n"); | |
829 | } | |
830 | else | |
831 | Warning("DumpPStack", "No particle with id %d.", i); | |
832 | } | |
833 | ||
88cb7938 | 834 | printf("\n=======================================================================\n\n"); |
9e1a0ddb | 835 | |
836 | // print particle file map | |
837 | printf("\nParticle file map: \n"); | |
3ab6f951 | 838 | for (i=0; i<fNtrack; i++) |
9e1a0ddb | 839 | printf(" %d th entry: %d \n",i,fParticleFileMap[i]); |
840 | } | |
841 | ||
842 | ||
843 | //_____________________________________________________________________________ | |
844 | void AliStack::DumpLoadedStack() const | |
845 | { | |
846 | // | |
847 | // Dumps the particle in the stack | |
848 | // that are loaded in memory. | |
849 | // | |
850 | ||
851 | TObjArray &particles = *fParticleMap; | |
852 | printf( | |
853 | "\n\n=======================================================================\n"); | |
854 | for (Int_t i=0;i<fNtrack;i++) | |
855 | { | |
e2afb3b6 | 856 | TParticle* particle = dynamic_cast<TParticle*>(particles[i]); |
9e1a0ddb | 857 | if (particle) { |
858 | printf("-> %d ",i); particle->Print(); | |
859 | printf("--------------------------------------------------------------\n"); | |
860 | } | |
861 | else { | |
862 | printf("-> %d Particle not loaded.\n",i); | |
863 | printf("--------------------------------------------------------------\n"); | |
864 | } | |
865 | } | |
866 | printf( | |
867 | "\n=======================================================================\n\n"); | |
868 | } | |
869 | ||
870 | // | |
871 | // protected methods | |
872 | // | |
873 | ||
874 | //_____________________________________________________________________________ | |
875 | void AliStack::CleanParents() | |
876 | { | |
877 | // | |
878 | // Clean particles stack | |
879 | // Set parent/daughter relations | |
880 | // | |
881 | ||
882 | TObjArray &particles = *fParticleMap; | |
883 | TParticle *part; | |
884 | int i; | |
885 | for(i=0; i<fHgwmk+1; i++) { | |
e2afb3b6 | 886 | part = dynamic_cast<TParticle*>(particles.At(i)); |
9e1a0ddb | 887 | if(part) if(!part->TestBit(kDaughtersBit)) { |
888 | part->SetFirstDaughter(-1); | |
889 | part->SetLastDaughter(-1); | |
890 | } | |
891 | } | |
892 | } | |
893 | ||
894 | //_____________________________________________________________________________ | |
895 | TParticle* AliStack::GetNextParticle() | |
896 | { | |
897 | // | |
898 | // Return next particle from stack of particles | |
899 | // | |
900 | ||
901 | TParticle* particle = 0; | |
902 | ||
903 | // search secondaries | |
904 | //for(Int_t i=fNtrack-1; i>=0; i--) { | |
905 | for(Int_t i=fNtrack-1; i>fHgwmk; i--) { | |
e2afb3b6 | 906 | particle = dynamic_cast<TParticle*>(fParticleMap->At(i)); |
9e1a0ddb | 907 | if ((particle) && (!particle->TestBit(kDoneBit))) { |
908 | fCurrent=i; | |
9e1a0ddb | 909 | return particle; |
910 | } | |
911 | } | |
912 | ||
913 | // take next primary if all secondaries were done | |
914 | while (fCurrentPrimary>=0) { | |
915 | fCurrent = fCurrentPrimary; | |
e2afb3b6 | 916 | particle = dynamic_cast<TParticle*>(fParticleMap->At(fCurrentPrimary--)); |
9e1a0ddb | 917 | if ((particle) && (!particle->TestBit(kDoneBit))) { |
9e1a0ddb | 918 | return particle; |
919 | } | |
920 | } | |
921 | ||
922 | // nothing to be tracked | |
923 | fCurrent = -1; | |
6308b375 | 924 | |
925 | ||
9e1a0ddb | 926 | return particle; |
927 | } | |
88cb7938 | 928 | //__________________________________________________________________________________________ |
9e1a0ddb | 929 | |
88cb7938 | 930 | TTree* AliStack::TreeK() |
931 | { | |
932 | //returns TreeK | |
933 | if (fTreeK) | |
934 | { | |
935 | return fTreeK; | |
936 | } | |
937 | else | |
938 | { | |
939 | AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName); | |
940 | if (rl == 0x0) | |
941 | { | |
21bf7095 | 942 | AliFatal(Form("Can not get RunLoader from event folder named %s",fEventFolderName.Data())); |
88cb7938 | 943 | return 0x0;//pro forma |
944 | } | |
945 | fTreeK = rl->TreeK(); | |
504b172d | 946 | if ( fTreeK ) |
88cb7938 | 947 | { |
504b172d | 948 | ConnectTree(); |
88cb7938 | 949 | } |
950 | else | |
951 | { | |
504b172d | 952 | //don't panic - could be Lego |
21bf7095 | 953 | AliWarning(Form("Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data())); |
88cb7938 | 954 | } |
955 | } | |
956 | return fTreeK;//never reached | |
957 | } | |
9e1a0ddb | 958 | //__________________________________________________________________________________________ |
88cb7938 | 959 | |
960 | void AliStack::ConnectTree() | |
9e1a0ddb | 961 | { |
962 | // | |
88cb7938 | 963 | // Creates branch for writing particles |
964 | // | |
21bf7095 | 965 | AliDebug(1, "Connecting TreeK"); |
88cb7938 | 966 | if (fTreeK == 0x0) |
967 | { | |
968 | if (TreeK() == 0x0) | |
969 | { | |
21bf7095 | 970 | AliFatal("Parameter is NULL");//we don't like such a jokes |
88cb7938 | 971 | return; |
972 | } | |
973 | return;//in this case TreeK() calls back this method (ConnectTree) | |
974 | //tree after setting fTreeK, the rest was already executed | |
975 | //it is safe to return now | |
976 | } | |
977 | ||
978 | // Create a branch for particles | |
979 | ||
21bf7095 | 980 | AliDebug(2, Form("Tree name is %s",fTreeK->GetName())); |
88cb7938 | 981 | |
982 | if (fTreeK->GetDirectory()) | |
983 | { | |
21bf7095 | 984 | AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName())); |
88cb7938 | 985 | } |
986 | else | |
21bf7095 | 987 | AliWarning("DIR IS NOT SET !!!"); |
88cb7938 | 988 | |
024a7e64 | 989 | TBranch *branch=fTreeK->GetBranch(AliRunLoader::GetKineBranchName()); |
88cb7938 | 990 | if(branch == 0x0) |
991 | { | |
024a7e64 | 992 | branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000); |
21bf7095 | 993 | AliDebug(2, "Creating Branch in Tree"); |
88cb7938 | 994 | } |
995 | else | |
996 | { | |
21bf7095 | 997 | AliDebug(2, "Branch Found in Tree"); |
88cb7938 | 998 | branch->SetAddress(&fParticleBuffer); |
999 | } | |
1000 | if (branch->GetDirectory()) | |
1001 | { | |
21bf7095 | 1002 | AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName())); |
88cb7938 | 1003 | } |
1004 | else | |
21bf7095 | 1005 | AliWarning("Branch Dir is NOT SET"); |
9e1a0ddb | 1006 | } |
88cb7938 | 1007 | //__________________________________________________________________________________________ |
9e1a0ddb | 1008 | |
88cb7938 | 1009 | |
1010 | void AliStack::BeginEvent() | |
9e1a0ddb | 1011 | { |
1012 | // start a new event | |
88cb7938 | 1013 | Reset(); |
9e1a0ddb | 1014 | } |
1015 | ||
81694e6a | 1016 | //_____________________________________________________________________________ |
9e1a0ddb | 1017 | void AliStack::FinishRun() |
1018 | { | |
1019 | // Clean TreeK information | |
9e1a0ddb | 1020 | } |
5d8718b8 | 1021 | //_____________________________________________________________________________ |
88cb7938 | 1022 | |
1023 | Bool_t AliStack::GetEvent() | |
9e1a0ddb | 1024 | { |
1025 | // | |
1026 | // Get new event from TreeK | |
ccf7a81f | 1027 | |
9e1a0ddb | 1028 | // Reset/Create the particle stack |
88cb7938 | 1029 | fTreeK = 0x0; |
9e1a0ddb | 1030 | |
88cb7938 | 1031 | if (TreeK() == 0x0) //forces connecting |
1032 | { | |
21bf7095 | 1033 | AliError("cannot find Kine Tree for current event"); |
ccf7a81f | 1034 | return kFALSE; |
88cb7938 | 1035 | } |
ccf7a81f | 1036 | |
88cb7938 | 1037 | Int_t size = (Int_t)TreeK()->GetEntries(); |
9e1a0ddb | 1038 | ResetArrays(size); |
ccf7a81f | 1039 | return kTRUE; |
9e1a0ddb | 1040 | } |
88cb7938 | 1041 | //_____________________________________________________________________________ |
1042 | ||
1043 | void AliStack::SetEventFolderName(const char* foldname) | |
1044 | { | |
1045 | //Sets event folder name | |
1046 | fEventFolderName = foldname; | |
1047 | } | |
1fed0e8b | 1048 | |
1049 | Bool_t AliStack::IsStable(Int_t pdg) const | |
1050 | { | |
1051 | // | |
1052 | // Decide whether particle (pdg) is stable | |
1053 | // | |
1054 | ||
1055 | const Int_t kNstable = 14; | |
1056 | Int_t i; | |
1057 | ||
1058 | Int_t pdgStable[kNstable] = { | |
1059 | kGamma, // Photon | |
1060 | kElectron, // Electron | |
1061 | kMuonPlus, // Muon | |
1062 | kPiPlus, // Pion | |
1063 | kKPlus, // Kaon | |
1064 | kProton, // Proton | |
1065 | kNeutron, // Neutron | |
1066 | kLambda0, // Lambda_0 | |
1067 | kSigmaMinus, // Sigma Minus | |
1068 | kSigma0, // Sigma_0 | |
1069 | kSigmaPlus, // Sigma Plus | |
1070 | 3312, // Xsi Minus | |
1071 | 3322, // Xsi | |
1072 | 3334, // Omega | |
1073 | }; | |
1074 | ||
1075 | Bool_t isStable = kFALSE; | |
1076 | for (i = 0; i < kNstable; i++) { | |
1077 | if (pdg == TMath::Abs(pdgStable[i])) { | |
1078 | isStable = kTRUE; | |
1079 | break; | |
1080 | } | |
1081 | } | |
1082 | ||
1083 | return isStable; | |
1084 | } | |
1085 | ||
1086 | Bool_t AliStack::IsPhysicalPrimary(Int_t index) | |
1087 | { | |
1088 | // | |
1089 | // Test if a particle is a physical primary according to the following definition: | |
1090 | // Particles produced in the collision including products of strong and | |
1091 | // electromagnetic decay and excluding feed-down from weak decays of strange | |
1092 | // particles. | |
1093 | // | |
1094 | TParticle* p = Particle(index); | |
1095 | Int_t ist = p->GetStatusCode(); | |
1096 | ||
1097 | // | |
1098 | // Initial state particle | |
1099 | if (ist > 20) return kFALSE; | |
1100 | ||
1101 | Int_t pdg = TMath::Abs(p->GetPdgCode()); | |
1102 | ||
1103 | if (!IsStable(pdg)) return kFALSE; | |
1104 | ||
1105 | if (index < GetNprimary()) { | |
1106 | // | |
1107 | // Particle produced by generator | |
1108 | return kTRUE; | |
1109 | } else { | |
1110 | // | |
1111 | // Particle produced during transport | |
1112 | // | |
1113 | // Check if this is a heavy flavor decay product | |
1114 | Int_t imo = p->GetFirstMother(); | |
1115 | TParticle* pm = Particle(imo); | |
1116 | Int_t mpdg = TMath::Abs(pm->GetPdgCode()); | |
1117 | Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg)))); | |
1118 | // | |
1119 | // Light hadron | |
1120 | if (mfl < 4) return kFALSE; | |
1121 | ||
1122 | // | |
1123 | // Heavy flavor hadron produced by generator | |
1124 | if (imo < GetNprimary()) { | |
1125 | return kTRUE; | |
1126 | } | |
1127 | ||
1128 | // To be sure that heavy flavor has not been produced in a secondary interaction | |
1129 | // Loop back to the generated mother | |
1130 | while (imo >= GetNprimary()) { | |
1131 | imo = p->GetFirstMother(); | |
1132 | pm = Particle(imo); | |
1133 | } | |
1134 | mpdg = TMath::Abs(pm->GetPdgCode()); | |
1135 | mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg)))); | |
1136 | ||
1137 | if (mfl < 4) { | |
1138 | return kFALSE; | |
1139 | } else { | |
1140 | return kTRUE; | |
1141 | } | |
1142 | } // produced by generator ? | |
1143 | } | |
1144 |