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