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 | // // |
20 | // Particles stack class |
21 | // // |
22 | /////////////////////////////////////////////////////////////////////////////// |
23 | |
9e1a0ddb |
24 | |
88cb7938 |
25 | #include <TFile.h> |
26 | #include <TFolder.h> |
9e1a0ddb |
27 | #include <TObjArray.h> |
28 | #include <TParticle.h> |
fe046ade |
29 | #include <TParticlePDG.h> |
9e1a0ddb |
30 | #include <TTree.h> |
9e1a0ddb |
31 | |
9e1a0ddb |
32 | #include "AliHit.h" |
88cb7938 |
33 | #include "AliModule.h" |
34 | #include "AliRun.h" |
5d12ce38 |
35 | #include "AliMC.h" |
88cb7938 |
36 | #include "AliRunLoader.h" |
37 | #include "AliStack.h" |
9e1a0ddb |
38 | |
39 | ClassImp(AliStack) |
40 | |
e2afb3b6 |
41 | //_______________________________________________________________________ |
42 | AliStack::AliStack(): |
43 | fParticles(0), |
44 | fParticleMap(0), |
45 | fParticleFileMap(0), |
46 | fParticleBuffer(0), |
47 | fTreeK(0), |
48 | fNtrack(0), |
49 | fNprimary(0), |
50 | fCurrent(-1), |
51 | fCurrentPrimary(-1), |
52 | fHgwmk(0), |
88cb7938 |
53 | fLoadPoint(0), |
e191bb57 |
54 | fEventFolderName(AliConfig::GetDefaultEventFolderName()) |
9e1a0ddb |
55 | { |
56 | // |
e2afb3b6 |
57 | // Default constructor |
9e1a0ddb |
58 | // |
9e1a0ddb |
59 | } |
60 | |
e2afb3b6 |
61 | //_______________________________________________________________________ |
88cb7938 |
62 | AliStack::AliStack(Int_t size, const char* evfoldname): |
e2afb3b6 |
63 | fParticles(new TClonesArray("TParticle",1000)), |
64 | fParticleMap(new TObjArray(size)), |
65 | fParticleFileMap(0), |
66 | fParticleBuffer(0), |
67 | fTreeK(0), |
68 | fNtrack(0), |
69 | fNprimary(0), |
70 | fCurrent(-1), |
71 | fCurrentPrimary(-1), |
72 | fHgwmk(0), |
88cb7938 |
73 | fLoadPoint(0), |
74 | fEventFolderName(evfoldname) |
e2afb3b6 |
75 | { |
76 | // |
77 | // Constructor |
78 | // |
79 | } |
9e1a0ddb |
80 | |
e2afb3b6 |
81 | //_______________________________________________________________________ |
82 | AliStack::AliStack(const AliStack& st): |
83 | TVirtualMCStack(st), |
84 | fParticles(0), |
85 | fParticleMap(0), |
86 | fParticleFileMap(0), |
87 | fParticleBuffer(0), |
88 | fTreeK(0), |
89 | fNtrack(0), |
90 | fNprimary(0), |
91 | fCurrent(-1), |
92 | fCurrentPrimary(-1), |
93 | fHgwmk(0), |
94 | fLoadPoint(0) |
9e1a0ddb |
95 | { |
96 | // |
e2afb3b6 |
97 | // Copy constructor |
9e1a0ddb |
98 | // |
e2afb3b6 |
99 | st.Copy(*this); |
9e1a0ddb |
100 | } |
101 | |
e2afb3b6 |
102 | //_______________________________________________________________________ |
6c4904c2 |
103 | void AliStack::Copy(TObject&) const |
e2afb3b6 |
104 | { |
105 | Fatal("Copy","Not implemented!\n"); |
106 | } |
9e1a0ddb |
107 | |
e2afb3b6 |
108 | //_______________________________________________________________________ |
9e1a0ddb |
109 | AliStack::~AliStack() |
110 | { |
111 | // |
112 | // Destructor |
113 | // |
114 | |
115 | if (fParticles) { |
116 | fParticles->Delete(); |
117 | delete fParticles; |
118 | } |
119 | delete fParticleMap; |
88cb7938 |
120 | //PH??? delete fTreeK; |
9e1a0ddb |
121 | } |
122 | |
123 | // |
124 | // public methods |
125 | // |
126 | |
127 | //_____________________________________________________________________________ |
642f15cf |
128 | void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, |
88cb7938 |
129 | Float_t *vpos, Float_t *polar, Float_t tof, |
130 | TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) |
9e1a0ddb |
131 | { |
132 | // |
133 | // Load a track on the stack |
134 | // |
135 | // done 0 if the track has to be transported |
136 | // 1 if not |
137 | // parent identifier of the parent track. -1 for a primary |
138 | // pdg particle code |
139 | // pmom momentum GeV/c |
140 | // vpos position |
141 | // polar polarisation |
142 | // tof time of flight in seconds |
143 | // mecha production mechanism |
144 | // ntr on output the number of the track stored |
145 | // |
146 | |
9e1a0ddb |
147 | // const Float_t tlife=0; |
148 | |
149 | // |
150 | // Here we get the static mass |
151 | // For MC is ok, but a more sophisticated method could be necessary |
152 | // if the calculated mass is required |
153 | // also, this method is potentially dangerous if the mass |
154 | // used in the MC is not the same of the PDG database |
155 | // |
fe046ade |
156 | TParticlePDG* pmc = TDatabasePDG::Instance()->GetParticle(pdg); |
157 | if (pmc) { |
158 | Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); |
159 | Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+ |
160 | pmom[1]*pmom[1]+pmom[2]*pmom[2]); |
161 | |
ccf7a81f |
162 | // 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", |
163 | // mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); |
164 | |
9e1a0ddb |
165 | |
642f15cf |
166 | PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e, |
fe046ade |
167 | vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2], |
168 | mech, ntr, weight, is); |
169 | } else { |
642f15cf |
170 | Warning("PushTrack", "Particle type %d not defined in PDG Database !\n", pdg); |
171 | Warning("PushTrack", "Particle skipped !\n"); |
fe046ade |
172 | } |
9e1a0ddb |
173 | } |
174 | |
175 | //_____________________________________________________________________________ |
642f15cf |
176 | void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, |
9e1a0ddb |
177 | Double_t px, Double_t py, Double_t pz, Double_t e, |
178 | Double_t vx, Double_t vy, Double_t vz, Double_t tof, |
179 | Double_t polx, Double_t poly, Double_t polz, |
98490ea9 |
180 | TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is) |
9e1a0ddb |
181 | { |
182 | // |
183 | // Load a track on the stack |
184 | // |
185 | // done 0 if the track has to be transported |
186 | // 1 if not |
187 | // parent identifier of the parent track. -1 for a primary |
188 | // pdg particle code |
189 | // kS generation status code |
190 | // px, py, pz momentum GeV/c |
191 | // vx, vy, vz position |
192 | // polar polarisation |
193 | // tof time of flight in seconds |
194 | // mech production mechanism |
195 | // ntr on output the number of the track stored |
196 | // |
197 | // New method interface: |
198 | // arguments were changed to be in correspondence with TParticle |
199 | // constructor. |
200 | // Note: the energy is not calculated from the static mass but |
201 | // it is passed by argument e. |
202 | |
203 | |
9e1a0ddb |
204 | const Int_t kFirstDaughter=-1; |
205 | const Int_t kLastDaughter=-1; |
206 | |
207 | TClonesArray &particles = *fParticles; |
208 | TParticle* particle |
209 | = new(particles[fLoadPoint++]) |
47c8bcbe |
210 | TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter, |
9e1a0ddb |
211 | px, py, pz, e, vx, vy, vz, tof); |
212 | |
213 | particle->SetPolarisation(polx, poly, polz); |
214 | particle->SetWeight(weight); |
215 | particle->SetUniqueID(mech); |
216 | |
217 | if(!done) particle->SetBit(kDoneBit); |
218 | |
219 | // Declare that the daughter information is valid |
220 | particle->SetBit(kDaughtersBit); |
221 | // Add the particle to the stack |
222 | fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!! |
223 | |
224 | if(parent>=0) { |
e2afb3b6 |
225 | particle = dynamic_cast<TParticle*>(fParticleMap->At(parent)); |
1305ca2f |
226 | if (particle) { |
227 | particle->SetLastDaughter(fNtrack); |
228 | if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); |
229 | } |
230 | else { |
642f15cf |
231 | printf("Error in AliStack::PushTrack: Parent %d does not exist\n",parent); |
1305ca2f |
232 | } |
9e1a0ddb |
233 | } |
234 | else { |
235 | // |
236 | // This is a primary track. Set high water mark for this event |
237 | fHgwmk = fNtrack; |
238 | // |
239 | // Set also number if primary tracks |
240 | fNprimary = fHgwmk+1; |
241 | fCurrentPrimary++; |
242 | } |
243 | ntr = fNtrack++; |
244 | } |
245 | |
246 | //_____________________________________________________________________________ |
642f15cf |
247 | TParticle* AliStack::PopNextTrack(Int_t& itrack) |
b9d0a01d |
248 | { |
249 | // |
250 | // Returns next track from stack of particles |
251 | // |
252 | |
253 | |
254 | TParticle* track = GetNextParticle(); |
255 | |
256 | if (track) { |
257 | itrack = fCurrent; |
258 | track->SetBit(kDoneBit); |
259 | } |
fe046ade |
260 | else |
b9d0a01d |
261 | itrack = -1; |
fe046ade |
262 | |
263 | fCurrentTrack = track; |
b9d0a01d |
264 | return track; |
265 | } |
266 | |
b9d0a01d |
267 | //_____________________________________________________________________________ |
642f15cf |
268 | TParticle* AliStack::PopPrimaryForTracking(Int_t i) |
b9d0a01d |
269 | { |
270 | // |
271 | // Returns i-th primary particle if it is flagged to be tracked, |
272 | // 0 otherwise |
273 | // |
274 | |
275 | TParticle* particle = Particle(i); |
276 | |
277 | if (!particle->TestBit(kDoneBit)) |
278 | return particle; |
279 | else |
280 | return 0; |
281 | } |
282 | |
9e1a0ddb |
283 | //_____________________________________________________________________________ |
284 | void AliStack::PurifyKine() |
285 | { |
286 | // |
287 | // Compress kinematic tree keeping only flagged particles |
288 | // and renaming the particle id's in all the hits |
289 | // |
290 | |
291 | TObjArray &particles = *fParticleMap; |
292 | int nkeep=fHgwmk+1, parent, i; |
293 | TParticle *part, *father; |
294 | TArrayI map(particles.GetLast()+1); |
295 | |
296 | // Save in Header total number of tracks before compression |
297 | |
298 | // If no tracks generated return now |
299 | if(fHgwmk+1 == fNtrack) return; |
300 | |
9e1a0ddb |
301 | // First pass, invalid Daughter information |
302 | for(i=0; i<fNtrack; i++) { |
303 | // Preset map, to be removed later |
304 | if(i<=fHgwmk) map[i]=i ; |
305 | else { |
306 | map[i] = -99; |
e2afb3b6 |
307 | if((part=dynamic_cast<TParticle*>(particles.At(i)))) { |
e44c5340 |
308 | // |
309 | // Check of this track should be kept for physics reasons |
88cb7938 |
310 | if (KeepPhysics(part)) KeepTrack(i); |
e44c5340 |
311 | // |
76969085 |
312 | part->ResetBit(kDaughtersBit); |
313 | part->SetFirstDaughter(-1); |
314 | part->SetLastDaughter(-1); |
315 | } |
9e1a0ddb |
316 | } |
317 | } |
318 | // Invalid daughter information for the parent of the first particle |
319 | // generated. This may or may not be the current primary according to |
320 | // whether decays have been recorded among the primaries |
e2afb3b6 |
321 | part = dynamic_cast<TParticle*>(particles.At(fHgwmk+1)); |
9e1a0ddb |
322 | particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit); |
323 | // Second pass, build map between old and new numbering |
324 | for(i=fHgwmk+1; i<fNtrack; i++) { |
325 | if(particles.At(i)->TestBit(kKeepBit)) { |
326 | |
327 | // This particle has to be kept |
328 | map[i]=nkeep; |
329 | // If old and new are different, have to move the pointer |
330 | if(i!=nkeep) particles[nkeep]=particles.At(i); |
e2afb3b6 |
331 | part = dynamic_cast<TParticle*>(particles.At(nkeep)); |
9e1a0ddb |
332 | |
333 | // as the parent is always *before*, it must be already |
334 | // in place. This is what we are checking anyway! |
335 | if((parent=part->GetFirstMother())>fHgwmk) |
336 | if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent); |
337 | else part->SetFirstMother(map[parent]); |
338 | |
339 | nkeep++; |
340 | } |
341 | } |
342 | |
343 | // Fix daughters information |
344 | for (i=fHgwmk+1; i<nkeep; i++) { |
e2afb3b6 |
345 | part = dynamic_cast<TParticle*>(particles.At(i)); |
9e1a0ddb |
346 | parent = part->GetFirstMother(); |
347 | if(parent>=0) { |
e2afb3b6 |
348 | father = dynamic_cast<TParticle*>(particles.At(parent)); |
9e1a0ddb |
349 | if(father->TestBit(kDaughtersBit)) { |
350 | |
351 | if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i); |
352 | if(i>father->GetLastDaughter()) father->SetLastDaughter(i); |
353 | } else { |
354 | // Initialise daughters info for first pass |
355 | father->SetFirstDaughter(i); |
356 | father->SetLastDaughter(i); |
357 | father->SetBit(kDaughtersBit); |
358 | } |
359 | } |
360 | } |
361 | |
362 | // Now loop on all registered hit lists |
5d12ce38 |
363 | TList* hitLists = gAlice->GetMCApp()->GetHitLists(); |
9e1a0ddb |
364 | TIter next(hitLists); |
365 | TCollection *hitList; |
e2afb3b6 |
366 | while((hitList = dynamic_cast<TCollection*>(next()))) { |
9e1a0ddb |
367 | TIter nexthit(hitList); |
368 | AliHit *hit; |
e2afb3b6 |
369 | while((hit = dynamic_cast<AliHit*>(nexthit()))) { |
9e1a0ddb |
370 | hit->SetTrack(map[hit->GetTrack()]); |
371 | } |
372 | } |
373 | |
374 | // |
375 | // This for detectors which have a special mapping mechanism |
376 | // for hits, such as TPC and TRD |
377 | // |
2b22f272 |
378 | |
9e1a0ddb |
379 | TObjArray* modules = gAlice->Modules(); |
380 | TIter nextmod(modules); |
381 | AliModule *detector; |
e2afb3b6 |
382 | while((detector = dynamic_cast<AliModule*>(nextmod()))) { |
9e1a0ddb |
383 | detector->RemapTrackHitIDs(map.GetArray()); |
ae9676e4 |
384 | detector->RemapTrackReferencesIDs(map.GetArray()); |
9e1a0ddb |
385 | } |
2b22f272 |
386 | // |
5d12ce38 |
387 | gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray()); |
2b22f272 |
388 | |
9e1a0ddb |
389 | // Now the output bit, from fHgwmk to nkeep we write everything and we erase |
390 | if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5)); |
391 | |
392 | for (i=fHgwmk+1; i<nkeep; ++i) { |
e2afb3b6 |
393 | fParticleBuffer = dynamic_cast<TParticle*>(particles.At(i)); |
88cb7938 |
394 | fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries()); |
395 | TreeK()->Fill(); |
17f4a3d9 |
396 | particles[i]=fParticleBuffer=0; |
88cb7938 |
397 | } |
9e1a0ddb |
398 | |
399 | for (i=nkeep; i<fNtrack; ++i) particles[i]=0; |
400 | |
76969085 |
401 | Int_t toshrink = fNtrack-fHgwmk-1; |
9e1a0ddb |
402 | fLoadPoint-=toshrink; |
403 | for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i); |
404 | |
405 | fNtrack=nkeep; |
406 | fHgwmk=nkeep-1; |
407 | // delete [] map; |
408 | } |
409 | |
e44c5340 |
410 | Bool_t AliStack::KeepPhysics(TParticle* part) |
411 | { |
412 | // |
413 | // Some particles have to kept on the stack for reasons motivated |
414 | // by physics analysis. Decision is put here. |
415 | // |
416 | Bool_t keep = kFALSE; |
417 | // |
418 | // Keep first-generation daughter from primaries with heavy flavor |
419 | // |
420 | Int_t parent = part->GetFirstMother(); |
421 | if (parent >= 0 && parent <= fHgwmk) { |
e2afb3b6 |
422 | TParticle* father = dynamic_cast<TParticle*>(Particles()->At(parent)); |
e44c5340 |
423 | Int_t kf = father->GetPdgCode(); |
424 | kf = TMath::Abs(kf); |
425 | Int_t kfl = kf; |
426 | // meson ? |
427 | if (kfl > 10) kfl/=100; |
428 | // baryon |
429 | if (kfl > 10) kfl/=10; |
430 | if (kfl > 10) kfl/=10; |
431 | if (kfl >= 4) { |
432 | keep = kTRUE; |
433 | } |
434 | } |
435 | return keep; |
436 | } |
437 | |
9e1a0ddb |
438 | //_____________________________________________________________________________ |
439 | void AliStack::FinishEvent() |
440 | { |
88cb7938 |
441 | // |
442 | // Write out the kinematics that was not yet filled |
443 | // |
9e1a0ddb |
444 | |
88cb7938 |
445 | // Update event header |
9e1a0ddb |
446 | |
447 | |
88cb7938 |
448 | if (!TreeK()) { |
9e1a0ddb |
449 | // Fatal("FinishEvent", "No kinematics tree is defined."); |
450 | // Don't panic this is a probably a lego run |
451 | return; |
9e1a0ddb |
452 | } |
453 | |
454 | CleanParents(); |
88cb7938 |
455 | if(TreeK()->GetEntries() ==0) { |
9e1a0ddb |
456 | // set the fParticleFileMap size for the first time |
457 | fParticleFileMap.Set(fHgwmk+1); |
458 | } |
459 | |
460 | Bool_t allFilled = kFALSE; |
461 | TObject *part; |
462 | for(Int_t i=0; i<fHgwmk+1; ++i) |
463 | if((part=fParticleMap->At(i))) { |
e2afb3b6 |
464 | fParticleBuffer = dynamic_cast<TParticle*>(part); |
88cb7938 |
465 | fParticleFileMap[i]= static_cast<Int_t>(TreeK()->GetEntries()); |
466 | TreeK()->Fill(); |
da2af24c |
467 | //PH (*fParticleMap)[i]=fParticleBuffer=0; |
468 | fParticleBuffer=0; |
469 | fParticleMap->AddAt(0,i); |
9e1a0ddb |
470 | |
471 | // When all primaries were filled no particle!=0 |
472 | // should be left => to be removed later. |
473 | if (allFilled) printf("Why != 0 part # %d?\n",i); |
474 | } |
88cb7938 |
475 | else |
476 | { |
9e1a0ddb |
477 | // // printf("Why = 0 part # %d?\n",i); => We know. |
478 | // break; |
88cb7938 |
479 | // we don't break now in order to be sure there is no |
480 | // particle !=0 left. |
481 | // To be removed later and replaced with break. |
482 | if(!allFilled) allFilled = kTRUE; |
483 | } |
9e1a0ddb |
484 | } |
9e1a0ddb |
485 | //_____________________________________________________________________________ |
88cb7938 |
486 | |
9e1a0ddb |
487 | void AliStack::FlagTrack(Int_t track) |
488 | { |
489 | // |
490 | // Flags a track and all its family tree to be kept |
491 | // |
492 | |
493 | TParticle *particle; |
494 | |
495 | Int_t curr=track; |
496 | while(1) { |
e2afb3b6 |
497 | particle=dynamic_cast<TParticle*>(fParticleMap->At(curr)); |
9e1a0ddb |
498 | |
499 | // If the particle is flagged the three from here upward is saved already |
500 | if(particle->TestBit(kKeepBit)) return; |
501 | |
502 | // Save this particle |
503 | particle->SetBit(kKeepBit); |
504 | |
505 | // Move to father if any |
506 | if((curr=particle->GetFirstMother())==-1) return; |
507 | } |
508 | } |
509 | |
510 | //_____________________________________________________________________________ |
ed54bf31 |
511 | void AliStack::KeepTrack(Int_t track) |
9e1a0ddb |
512 | { |
513 | // |
514 | // Flags a track to be kept |
515 | // |
516 | |
517 | fParticleMap->At(track)->SetBit(kKeepBit); |
518 | } |
519 | |
520 | //_____________________________________________________________________________ |
521 | void AliStack::Reset(Int_t size) |
522 | { |
523 | // |
524 | // Resets stack |
525 | // |
88cb7938 |
526 | |
9e1a0ddb |
527 | fNtrack=0; |
528 | fNprimary=0; |
529 | fHgwmk=0; |
530 | fLoadPoint=0; |
531 | fCurrent = -1; |
88cb7938 |
532 | fTreeK = 0x0; |
9e1a0ddb |
533 | ResetArrays(size); |
534 | } |
535 | |
536 | //_____________________________________________________________________________ |
537 | void AliStack::ResetArrays(Int_t size) |
538 | { |
539 | // |
540 | // Resets stack arrays |
541 | // |
542 | |
be3380a3 |
543 | if (fParticles) |
544 | fParticles->Clear(); |
545 | else |
546 | fParticles = new TClonesArray("TParticle",1000); |
547 | if (fParticleMap) { |
548 | fParticleMap->Clear(); |
549 | if (size>0) fParticleMap->Expand(size);} |
550 | else |
551 | fParticleMap = new TObjArray(size); |
9e1a0ddb |
552 | } |
553 | |
554 | //_____________________________________________________________________________ |
e2afb3b6 |
555 | void AliStack::SetHighWaterMark(Int_t) |
9e1a0ddb |
556 | { |
557 | // |
558 | // Set high water mark for last track in event |
559 | // |
560 | |
561 | fHgwmk = fNtrack-1; |
562 | fCurrentPrimary=fHgwmk; |
563 | |
564 | // Set also number of primary tracks |
565 | fNprimary = fHgwmk+1; |
566 | fNtrack = fHgwmk+1; |
567 | } |
568 | |
569 | //_____________________________________________________________________________ |
570 | TParticle* AliStack::Particle(Int_t i) |
571 | { |
572 | // |
573 | // Return particle with specified ID |
574 | |
da2af24c |
575 | //PH if(!(*fParticleMap)[i]) { |
576 | if(!fParticleMap->At(i)) { |
911a3148 |
577 | Int_t nentries = fParticles->GetEntriesFast(); |
9e1a0ddb |
578 | // algorithmic way of getting entry index |
579 | // (primary particles are filled after secondaries) |
e94530da |
580 | Int_t entry = TreeKEntry(i); |
9e1a0ddb |
581 | // check whether algorithmic way and |
582 | // and the fParticleFileMap[i] give the same; |
583 | // give the fatal error if not |
584 | if (entry != fParticleFileMap[i]) { |
585 | Fatal("Particle", |
586 | "!! The algorithmic way and map are different: !!\n entry: %d map: %d", |
587 | entry, fParticleFileMap[i]); |
588 | } |
589 | |
88cb7938 |
590 | TreeK()->GetEntry(entry); |
9e1a0ddb |
591 | new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); |
592 | fParticleMap->AddAt((*fParticles)[nentries],i); |
593 | } |
e2afb3b6 |
594 | //PH return dynamic_cast<TParticle *>((*fParticleMap)[i]); |
595 | return dynamic_cast<TParticle*>(fParticleMap->At(i)); |
9e1a0ddb |
596 | } |
597 | |
e94530da |
598 | //_____________________________________________________________________________ |
599 | TParticle* AliStack::ParticleFromTreeK(Int_t id) const |
600 | { |
601 | // |
602 | // return pointer to TParticle with label id |
603 | // |
604 | Int_t entry; |
605 | if ((entry = TreeKEntry(id)) < 0) return 0; |
606 | if (fTreeK->GetEntry(entry)<=0) return 0; |
607 | return fParticleBuffer; |
608 | } |
609 | |
610 | //_____________________________________________________________________________ |
611 | Int_t AliStack::TreeKEntry(Int_t id) const |
612 | { |
613 | // |
614 | // return entry number in the TreeK for particle with label id |
615 | // return negative number if label>fNtrack |
616 | // |
617 | Int_t entry; |
618 | if (id<fNprimary) |
619 | entry = id+fNtrack-fNprimary; |
620 | else |
621 | entry = id-fNprimary; |
622 | return entry; |
623 | } |
624 | |
5d8718b8 |
625 | //_____________________________________________________________________________ |
642f15cf |
626 | Int_t AliStack::GetCurrentParentTrackNumber() const |
5d8718b8 |
627 | { |
628 | // |
629 | // Return number of the parent of the current track |
630 | // |
631 | |
632 | TParticle* current = (TParticle*)fParticleMap->At(fCurrent); |
633 | |
634 | if (current) |
635 | return current->GetFirstMother(); |
636 | else { |
642f15cf |
637 | Warning("GetCurrentParentTrackNumber", "Current track not found in the stack"); |
5d8718b8 |
638 | return -1; |
639 | } |
640 | } |
641 | |
9e1a0ddb |
642 | //_____________________________________________________________________________ |
f5f55563 |
643 | Int_t AliStack::GetPrimary(Int_t id) |
9e1a0ddb |
644 | { |
645 | // |
646 | // Return number of primary that has generated track |
647 | // |
648 | |
649 | int current, parent; |
9e1a0ddb |
650 | // |
651 | parent=id; |
652 | while (1) { |
653 | current=parent; |
f5f55563 |
654 | parent=Particle(current)->GetFirstMother(); |
9e1a0ddb |
655 | if(parent<0) return current; |
656 | } |
657 | } |
658 | |
659 | //_____________________________________________________________________________ |
660 | void AliStack::DumpPart (Int_t i) const |
661 | { |
662 | // |
663 | // Dumps particle i in the stack |
664 | // |
665 | |
e2afb3b6 |
666 | //PH dynamic_cast<TParticle*>((*fParticleMap)[i])->Print(); |
667 | dynamic_cast<TParticle*>(fParticleMap->At(i))->Print(); |
9e1a0ddb |
668 | } |
669 | |
670 | //_____________________________________________________________________________ |
671 | void AliStack::DumpPStack () |
672 | { |
673 | // |
674 | // Dumps the particle stack |
675 | // |
676 | |
3ab6f951 |
677 | Int_t i; |
678 | |
88cb7938 |
679 | printf("\n\n=======================================================================\n"); |
3ab6f951 |
680 | for (i=0;i<fNtrack;i++) |
9e1a0ddb |
681 | { |
682 | TParticle* particle = Particle(i); |
683 | if (particle) { |
684 | printf("-> %d ",i); particle->Print(); |
685 | printf("--------------------------------------------------------------\n"); |
686 | } |
687 | else |
688 | Warning("DumpPStack", "No particle with id %d.", i); |
689 | } |
690 | |
88cb7938 |
691 | printf("\n=======================================================================\n\n"); |
9e1a0ddb |
692 | |
693 | // print particle file map |
694 | printf("\nParticle file map: \n"); |
3ab6f951 |
695 | for (i=0; i<fNtrack; i++) |
9e1a0ddb |
696 | printf(" %d th entry: %d \n",i,fParticleFileMap[i]); |
697 | } |
698 | |
699 | |
700 | //_____________________________________________________________________________ |
701 | void AliStack::DumpLoadedStack() const |
702 | { |
703 | // |
704 | // Dumps the particle in the stack |
705 | // that are loaded in memory. |
706 | // |
707 | |
708 | TObjArray &particles = *fParticleMap; |
709 | printf( |
710 | "\n\n=======================================================================\n"); |
711 | for (Int_t i=0;i<fNtrack;i++) |
712 | { |
e2afb3b6 |
713 | TParticle* particle = dynamic_cast<TParticle*>(particles[i]); |
9e1a0ddb |
714 | if (particle) { |
715 | printf("-> %d ",i); particle->Print(); |
716 | printf("--------------------------------------------------------------\n"); |
717 | } |
718 | else { |
719 | printf("-> %d Particle not loaded.\n",i); |
720 | printf("--------------------------------------------------------------\n"); |
721 | } |
722 | } |
723 | printf( |
724 | "\n=======================================================================\n\n"); |
725 | } |
726 | |
727 | // |
728 | // protected methods |
729 | // |
730 | |
731 | //_____________________________________________________________________________ |
732 | void AliStack::CleanParents() |
733 | { |
734 | // |
735 | // Clean particles stack |
736 | // Set parent/daughter relations |
737 | // |
738 | |
739 | TObjArray &particles = *fParticleMap; |
740 | TParticle *part; |
741 | int i; |
742 | for(i=0; i<fHgwmk+1; i++) { |
e2afb3b6 |
743 | part = dynamic_cast<TParticle*>(particles.At(i)); |
9e1a0ddb |
744 | if(part) if(!part->TestBit(kDaughtersBit)) { |
745 | part->SetFirstDaughter(-1); |
746 | part->SetLastDaughter(-1); |
747 | } |
748 | } |
749 | } |
750 | |
751 | //_____________________________________________________________________________ |
752 | TParticle* AliStack::GetNextParticle() |
753 | { |
754 | // |
755 | // Return next particle from stack of particles |
756 | // |
757 | |
758 | TParticle* particle = 0; |
759 | |
760 | // search secondaries |
761 | //for(Int_t i=fNtrack-1; i>=0; i--) { |
762 | for(Int_t i=fNtrack-1; i>fHgwmk; i--) { |
e2afb3b6 |
763 | particle = dynamic_cast<TParticle*>(fParticleMap->At(i)); |
9e1a0ddb |
764 | if ((particle) && (!particle->TestBit(kDoneBit))) { |
765 | fCurrent=i; |
9e1a0ddb |
766 | return particle; |
767 | } |
768 | } |
769 | |
770 | // take next primary if all secondaries were done |
771 | while (fCurrentPrimary>=0) { |
772 | fCurrent = fCurrentPrimary; |
e2afb3b6 |
773 | particle = dynamic_cast<TParticle*>(fParticleMap->At(fCurrentPrimary--)); |
9e1a0ddb |
774 | if ((particle) && (!particle->TestBit(kDoneBit))) { |
9e1a0ddb |
775 | return particle; |
776 | } |
777 | } |
778 | |
779 | // nothing to be tracked |
780 | fCurrent = -1; |
9e1a0ddb |
781 | return particle; |
782 | } |
88cb7938 |
783 | //__________________________________________________________________________________________ |
9e1a0ddb |
784 | |
88cb7938 |
785 | TTree* AliStack::TreeK() |
786 | { |
787 | //returns TreeK |
788 | if (fTreeK) |
789 | { |
790 | return fTreeK; |
791 | } |
792 | else |
793 | { |
794 | AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName); |
795 | if (rl == 0x0) |
796 | { |
797 | Fatal("TreeK","Can not get RunLoader from event folder named %s",fEventFolderName.Data()); |
798 | return 0x0;//pro forma |
799 | } |
800 | fTreeK = rl->TreeK(); |
504b172d |
801 | if ( fTreeK ) |
88cb7938 |
802 | { |
504b172d |
803 | ConnectTree(); |
88cb7938 |
804 | } |
805 | else |
806 | { |
504b172d |
807 | //don't panic - could be Lego |
d0d4a6b3 |
808 | if (AliLoader::GetDebug()) |
504b172d |
809 | { |
810 | Warning("TreeK","Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data()); |
811 | } |
88cb7938 |
812 | } |
813 | } |
814 | return fTreeK;//never reached |
815 | } |
9e1a0ddb |
816 | //__________________________________________________________________________________________ |
88cb7938 |
817 | |
818 | void AliStack::ConnectTree() |
9e1a0ddb |
819 | { |
820 | // |
88cb7938 |
821 | // Creates branch for writing particles |
822 | // |
d0d4a6b3 |
823 | if (AliLoader::GetDebug()) Info("ConnectTree","Connecting TreeK"); |
88cb7938 |
824 | if (fTreeK == 0x0) |
825 | { |
826 | if (TreeK() == 0x0) |
827 | { |
828 | Fatal("ConnectTree","Parameter is NULL");//we don't like such a jokes |
829 | return; |
830 | } |
831 | return;//in this case TreeK() calls back this method (ConnectTree) |
832 | //tree after setting fTreeK, the rest was already executed |
833 | //it is safe to return now |
834 | } |
835 | |
836 | // Create a branch for particles |
837 | |
d0d4a6b3 |
838 | if (AliLoader::GetDebug()) |
88cb7938 |
839 | Info("ConnectTree","Tree name is %s",fTreeK->GetName()); |
840 | |
841 | if (fTreeK->GetDirectory()) |
842 | { |
d0d4a6b3 |
843 | if (AliLoader::GetDebug()) |
88cb7938 |
844 | Info("ConnectTree","and dir is %s",fTreeK->GetDirectory()->GetName()); |
845 | } |
846 | else |
847 | Warning("ConnectTree","DIR IS NOT SET !!!"); |
848 | |
024a7e64 |
849 | TBranch *branch=fTreeK->GetBranch(AliRunLoader::GetKineBranchName()); |
88cb7938 |
850 | if(branch == 0x0) |
851 | { |
024a7e64 |
852 | branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000); |
d0d4a6b3 |
853 | if (AliLoader::GetDebug()) Info("ConnectTree","Creating Branch in Tree"); |
88cb7938 |
854 | } |
855 | else |
856 | { |
d0d4a6b3 |
857 | if (AliLoader::GetDebug()) Info("ConnectTree","Branch Found in Tree"); |
88cb7938 |
858 | branch->SetAddress(&fParticleBuffer); |
859 | } |
860 | if (branch->GetDirectory()) |
861 | { |
d0d4a6b3 |
862 | if (AliLoader::GetDebug()) |
88cb7938 |
863 | Info("ConnectTree","Branch Dir Name is %s",branch->GetDirectory()->GetName()); |
864 | } |
865 | else |
866 | Warning("ConnectTree","Branch Dir is NOT SET"); |
9e1a0ddb |
867 | } |
88cb7938 |
868 | //__________________________________________________________________________________________ |
9e1a0ddb |
869 | |
88cb7938 |
870 | |
871 | void AliStack::BeginEvent() |
9e1a0ddb |
872 | { |
873 | // start a new event |
88cb7938 |
874 | Reset(); |
9e1a0ddb |
875 | } |
876 | |
81694e6a |
877 | //_____________________________________________________________________________ |
9e1a0ddb |
878 | void AliStack::FinishRun() |
879 | { |
880 | // Clean TreeK information |
9e1a0ddb |
881 | } |
5d8718b8 |
882 | //_____________________________________________________________________________ |
88cb7938 |
883 | |
884 | Bool_t AliStack::GetEvent() |
9e1a0ddb |
885 | { |
886 | // |
887 | // Get new event from TreeK |
ccf7a81f |
888 | |
9e1a0ddb |
889 | // Reset/Create the particle stack |
88cb7938 |
890 | fTreeK = 0x0; |
9e1a0ddb |
891 | |
88cb7938 |
892 | if (TreeK() == 0x0) //forces connecting |
893 | { |
894 | Error("GetEvent","cannot find Kine Tree for current event\n"); |
ccf7a81f |
895 | return kFALSE; |
88cb7938 |
896 | } |
ccf7a81f |
897 | |
88cb7938 |
898 | Int_t size = (Int_t)TreeK()->GetEntries(); |
9e1a0ddb |
899 | ResetArrays(size); |
ccf7a81f |
900 | return kTRUE; |
9e1a0ddb |
901 | } |
88cb7938 |
902 | //_____________________________________________________________________________ |
903 | |
904 | void AliStack::SetEventFolderName(const char* foldname) |
905 | { |
906 | //Sets event folder name |
907 | fEventFolderName = foldname; |
908 | } |