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