]>
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 | ||
16 | /* | |
17 | $Log$ | |
3ab6f951 | 18 | Revision 1.4 2001/05/25 07:25:20 hristov |
19 | AliStack destructor corrected (I.Hrivnacova) | |
20 | ||
81694e6a | 21 | Revision 1.3 2001/05/22 14:33:16 hristov |
22 | Minor changes | |
23 | ||
ed54bf31 | 24 | Revision 1.2 2001/05/17 05:49:39 fca |
25 | Reset pointers to daughters | |
26 | ||
76969085 | 27 | Revision 1.1 2001/05/16 14:57:22 alibrary |
28 | New files for folders and Stack | |
29 | ||
9e1a0ddb | 30 | */ |
31 | ||
32 | /////////////////////////////////////////////////////////////////////////////// | |
33 | // // | |
34 | // Particles stack class | |
35 | // // | |
36 | /////////////////////////////////////////////////////////////////////////////// | |
37 | ||
38 | #include <iostream.h> | |
39 | ||
40 | #include <TObjArray.h> | |
41 | #include <TParticle.h> | |
42 | #include <TTree.h> | |
43 | #include <TFile.h> | |
44 | #include <TFolder.h> | |
45 | #include <TROOT.h> | |
46 | ||
47 | #include "AliStack.h" | |
48 | #include "AliRun.h" | |
49 | #include "AliModule.h" | |
50 | #include "AliHit.h" | |
51 | //#include "ETrackBits.h" | |
52 | ||
53 | ClassImp(AliStack) | |
54 | ||
55 | //_____________________________________________________________________________ | |
56 | AliStack::AliStack(Int_t size) | |
57 | { | |
58 | // | |
59 | // Constructor | |
60 | // | |
61 | ||
62 | // Create the particles arrays | |
63 | fParticles = new TClonesArray("TParticle",1000); | |
64 | fParticleMap = new TObjArray(size); | |
65 | fParticleBuffer = new TParticle(); | |
66 | fNtrack = 0; | |
67 | fNprimary = 0; | |
68 | fCurrent = -1; | |
69 | fCurrentPrimary = -1; | |
70 | fTreeK = 0; | |
71 | } | |
72 | ||
73 | ||
74 | //_____________________________________________________________________________ | |
75 | AliStack::AliStack() | |
76 | { | |
77 | // | |
78 | // Default constructor | |
79 | // | |
80 | ||
81 | // Create the particles arrays | |
82 | fParticles = new TClonesArray("TParticle",1000); | |
83 | fParticleMap = new TObjArray(10000); | |
84 | fParticleBuffer = new TParticle(); | |
85 | fNtrack = 0; | |
86 | fCurrent = -1; | |
87 | fNprimary = 0; | |
88 | fCurrentPrimary = -1; | |
89 | fTreeK = 0; | |
90 | } | |
91 | ||
92 | ||
93 | //_____________________________________________________________________________ | |
94 | AliStack::~AliStack() | |
95 | { | |
96 | // | |
97 | // Destructor | |
98 | // | |
99 | ||
81694e6a | 100 | delete fParticleBuffer; |
9e1a0ddb | 101 | if (fParticles) { |
102 | fParticles->Delete(); | |
103 | delete fParticles; | |
104 | } | |
105 | delete fParticleMap; | |
9e1a0ddb | 106 | delete fTreeK; |
107 | } | |
108 | ||
109 | // | |
110 | // public methods | |
111 | // | |
112 | ||
113 | //_____________________________________________________________________________ | |
114 | void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, | |
115 | Float_t *vpos, Float_t *polar, Float_t tof, | |
116 | AliMCProcess mech, Int_t &ntr, Float_t weight) | |
117 | { | |
118 | // | |
119 | // Load a track on the stack | |
120 | // | |
121 | // done 0 if the track has to be transported | |
122 | // 1 if not | |
123 | // parent identifier of the parent track. -1 for a primary | |
124 | // pdg particle code | |
125 | // pmom momentum GeV/c | |
126 | // vpos position | |
127 | // polar polarisation | |
128 | // tof time of flight in seconds | |
129 | // mecha production mechanism | |
130 | // ntr on output the number of the track stored | |
131 | // | |
132 | ||
133 | Float_t mass; | |
134 | const Int_t kfirstdaughter=-1; | |
135 | const Int_t klastdaughter=-1; | |
136 | const Int_t kS=0; | |
137 | // const Float_t tlife=0; | |
138 | ||
139 | // | |
140 | // Here we get the static mass | |
141 | // For MC is ok, but a more sophisticated method could be necessary | |
142 | // if the calculated mass is required | |
143 | // also, this method is potentially dangerous if the mass | |
144 | // used in the MC is not the same of the PDG database | |
145 | // | |
146 | mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
147 | Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+ | |
148 | pmom[1]*pmom[1]+pmom[2]*pmom[2]); | |
149 | ||
150 | // 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", | |
151 | // mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); | |
152 | ||
153 | TClonesArray &particles = *fParticles; | |
154 | TParticle* particle | |
155 | = new(particles[fLoadPoint++]) | |
156 | TParticle(pdg, kS, parent, -1, kfirstdaughter, klastdaughter, | |
157 | pmom[0], pmom[1], pmom[2], e, vpos[0], vpos[1], vpos[2], tof); | |
158 | particle->SetPolarisation(TVector3(polar[0],polar[1],polar[2])); | |
159 | particle->SetWeight(weight); | |
160 | particle->SetUniqueID(mech); | |
161 | if(!done) particle->SetBit(kDoneBit); | |
162 | ||
163 | ||
164 | // Declare that the daughter information is valid | |
165 | particle->SetBit(kDaughtersBit); | |
166 | // Add the particle to the stack | |
167 | fParticleMap->AddAtAndExpand(particle, fNtrack); | |
168 | ||
169 | if(parent>=0) { | |
170 | particle = (TParticle*) fParticleMap->At(parent); | |
171 | particle->SetLastDaughter(fNtrack); | |
172 | if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); | |
173 | } | |
174 | else { | |
175 | // | |
176 | // This is a primary track. Set high water mark for this event | |
177 | fHgwmk = fNtrack; | |
178 | // | |
179 | // Set also number if primary tracks | |
180 | fNprimary = fHgwmk+1; | |
181 | fCurrentPrimary++; | |
182 | } | |
183 | ntr = fNtrack++; | |
184 | } | |
185 | ||
186 | //_____________________________________________________________________________ | |
187 | void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, | |
188 | Double_t px, Double_t py, Double_t pz, Double_t e, | |
189 | Double_t vx, Double_t vy, Double_t vz, Double_t tof, | |
190 | Double_t polx, Double_t poly, Double_t polz, | |
191 | AliMCProcess mech, Int_t &ntr, Float_t weight) | |
192 | { | |
193 | // | |
194 | // Load a track on the stack | |
195 | // | |
196 | // done 0 if the track has to be transported | |
197 | // 1 if not | |
198 | // parent identifier of the parent track. -1 for a primary | |
199 | // pdg particle code | |
200 | // kS generation status code | |
201 | // px, py, pz momentum GeV/c | |
202 | // vx, vy, vz position | |
203 | // polar polarisation | |
204 | // tof time of flight in seconds | |
205 | // mech production mechanism | |
206 | // ntr on output the number of the track stored | |
207 | // | |
208 | // New method interface: | |
209 | // arguments were changed to be in correspondence with TParticle | |
210 | // constructor. | |
211 | // Note: the energy is not calculated from the static mass but | |
212 | // it is passed by argument e. | |
213 | ||
214 | ||
215 | const Int_t kS=0; | |
216 | const Int_t kFirstDaughter=-1; | |
217 | const Int_t kLastDaughter=-1; | |
218 | ||
219 | TClonesArray &particles = *fParticles; | |
220 | TParticle* particle | |
221 | = new(particles[fLoadPoint++]) | |
222 | TParticle(pdg, kS, parent, -1, kFirstDaughter, kLastDaughter, | |
223 | px, py, pz, e, vx, vy, vz, tof); | |
224 | ||
225 | particle->SetPolarisation(polx, poly, polz); | |
226 | particle->SetWeight(weight); | |
227 | particle->SetUniqueID(mech); | |
228 | ||
229 | if(!done) particle->SetBit(kDoneBit); | |
230 | ||
231 | // Declare that the daughter information is valid | |
232 | particle->SetBit(kDaughtersBit); | |
233 | // Add the particle to the stack | |
234 | fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!! | |
235 | ||
236 | if(parent>=0) { | |
237 | particle = (TParticle*) fParticleMap->At(parent); | |
238 | particle->SetLastDaughter(fNtrack); | |
239 | if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); | |
240 | } | |
241 | else { | |
242 | // | |
243 | // This is a primary track. Set high water mark for this event | |
244 | fHgwmk = fNtrack; | |
245 | // | |
246 | // Set also number if primary tracks | |
247 | fNprimary = fHgwmk+1; | |
248 | fCurrentPrimary++; | |
249 | } | |
250 | ntr = fNtrack++; | |
251 | } | |
252 | ||
253 | //_____________________________________________________________________________ | |
254 | void AliStack::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom, | |
255 | Float_t &e, Float_t *vpos, Float_t *polar, | |
256 | Float_t &tof) | |
257 | { | |
258 | // | |
259 | // Return next track from stack of particles | |
260 | // | |
261 | ||
262 | ||
263 | TParticle* track = GetNextParticle(); | |
264 | // cout << "GetNextTrack():" << fCurrent << fNprimary << endl; | |
265 | ||
266 | if(track) { | |
267 | mtrack=fCurrent; | |
268 | ipart=track->GetPdgCode(); | |
269 | pmom[0]=track->Px(); | |
270 | pmom[1]=track->Py(); | |
271 | pmom[2]=track->Pz(); | |
272 | e =track->Energy(); | |
273 | vpos[0]=track->Vx(); | |
274 | vpos[1]=track->Vy(); | |
275 | vpos[2]=track->Vz(); | |
276 | TVector3 pol; | |
277 | track->GetPolarisation(pol); | |
278 | polar[0]=pol.X(); | |
279 | polar[1]=pol.Y(); | |
280 | polar[2]=pol.Z(); | |
281 | tof=track->T(); | |
282 | track->SetBit(kDoneBit); | |
283 | //cout << "Filled params" << endl; | |
284 | } | |
285 | else | |
286 | mtrack=-1; | |
287 | ||
288 | // | |
289 | // stop and start timer when we start a primary track | |
290 | Int_t nprimaries = fNprimary; | |
291 | if (fCurrent >= nprimaries) return; | |
292 | if (fCurrent < nprimaries-1) { | |
293 | fTimer.Stop(); | |
294 | track=(TParticle*) fParticleMap->At(fCurrent+1); | |
295 | // track->SetProcessTime(fTimer.CpuTime()); | |
296 | } | |
297 | fTimer.Start(); | |
298 | } | |
299 | ||
300 | ||
301 | //_____________________________________________________________________________ | |
302 | void AliStack::PurifyKine() | |
303 | { | |
304 | // | |
305 | // Compress kinematic tree keeping only flagged particles | |
306 | // and renaming the particle id's in all the hits | |
307 | // | |
308 | ||
309 | TObjArray &particles = *fParticleMap; | |
310 | int nkeep=fHgwmk+1, parent, i; | |
311 | TParticle *part, *father; | |
312 | TArrayI map(particles.GetLast()+1); | |
313 | ||
314 | // Save in Header total number of tracks before compression | |
315 | ||
316 | // If no tracks generated return now | |
317 | if(fHgwmk+1 == fNtrack) return; | |
318 | ||
9e1a0ddb | 319 | // First pass, invalid Daughter information |
320 | for(i=0; i<fNtrack; i++) { | |
321 | // Preset map, to be removed later | |
322 | if(i<=fHgwmk) map[i]=i ; | |
323 | else { | |
324 | map[i] = -99; | |
76969085 | 325 | if((part=(TParticle*) particles.At(i))) { |
326 | part->ResetBit(kDaughtersBit); | |
327 | part->SetFirstDaughter(-1); | |
328 | part->SetLastDaughter(-1); | |
329 | } | |
9e1a0ddb | 330 | } |
331 | } | |
332 | // Invalid daughter information for the parent of the first particle | |
333 | // generated. This may or may not be the current primary according to | |
334 | // whether decays have been recorded among the primaries | |
335 | part = (TParticle *)particles.At(fHgwmk+1); | |
336 | particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit); | |
337 | // Second pass, build map between old and new numbering | |
338 | for(i=fHgwmk+1; i<fNtrack; i++) { | |
339 | if(particles.At(i)->TestBit(kKeepBit)) { | |
340 | ||
341 | // This particle has to be kept | |
342 | map[i]=nkeep; | |
343 | // If old and new are different, have to move the pointer | |
344 | if(i!=nkeep) particles[nkeep]=particles.At(i); | |
345 | part = (TParticle*) particles.At(nkeep); | |
346 | ||
347 | // as the parent is always *before*, it must be already | |
348 | // in place. This is what we are checking anyway! | |
349 | if((parent=part->GetFirstMother())>fHgwmk) | |
350 | if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent); | |
351 | else part->SetFirstMother(map[parent]); | |
352 | ||
353 | nkeep++; | |
354 | } | |
355 | } | |
356 | ||
357 | // Fix daughters information | |
358 | for (i=fHgwmk+1; i<nkeep; i++) { | |
359 | part = (TParticle *)particles.At(i); | |
360 | parent = part->GetFirstMother(); | |
361 | if(parent>=0) { | |
362 | father = (TParticle *)particles.At(parent); | |
363 | if(father->TestBit(kDaughtersBit)) { | |
364 | ||
365 | if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i); | |
366 | if(i>father->GetLastDaughter()) father->SetLastDaughter(i); | |
367 | } else { | |
368 | // Initialise daughters info for first pass | |
369 | father->SetFirstDaughter(i); | |
370 | father->SetLastDaughter(i); | |
371 | father->SetBit(kDaughtersBit); | |
372 | } | |
373 | } | |
374 | } | |
375 | ||
376 | // Now loop on all registered hit lists | |
377 | TList* hitLists = gAlice->GetHitLists(); | |
378 | TIter next(hitLists); | |
379 | TCollection *hitList; | |
380 | while((hitList = (TCollection*)next())) { | |
381 | TIter nexthit(hitList); | |
382 | AliHit *hit; | |
383 | while((hit = (AliHit*)nexthit())) { | |
384 | hit->SetTrack(map[hit->GetTrack()]); | |
385 | } | |
386 | } | |
387 | ||
388 | // | |
389 | // This for detectors which have a special mapping mechanism | |
390 | // for hits, such as TPC and TRD | |
391 | // | |
392 | ||
393 | TObjArray* modules = gAlice->Modules(); | |
394 | TIter nextmod(modules); | |
395 | AliModule *detector; | |
396 | while((detector = (AliModule*)nextmod())) { | |
397 | detector->RemapTrackHitIDs(map.GetArray()); | |
398 | } | |
399 | ||
400 | // Now the output bit, from fHgwmk to nkeep we write everything and we erase | |
401 | if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5)); | |
402 | ||
403 | for (i=fHgwmk+1; i<nkeep; ++i) { | |
404 | fParticleBuffer = (TParticle*) particles.At(i); | |
405 | fParticleFileMap[i]=(Int_t) fTreeK->GetEntries(); | |
406 | fTreeK->Fill(); | |
407 | particles[i]=0; | |
408 | } | |
409 | ||
410 | for (i=nkeep; i<fNtrack; ++i) particles[i]=0; | |
411 | ||
76969085 | 412 | Int_t toshrink = fNtrack-fHgwmk-1; |
9e1a0ddb | 413 | fLoadPoint-=toshrink; |
414 | for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i); | |
415 | ||
416 | fNtrack=nkeep; | |
417 | fHgwmk=nkeep-1; | |
418 | // delete [] map; | |
419 | } | |
420 | ||
421 | //_____________________________________________________________________________ | |
422 | void AliStack::FinishEvent() | |
423 | { | |
424 | // | |
425 | // Write out the kinematics that was not yet filled | |
426 | // | |
427 | ||
428 | // Update event header | |
429 | ||
430 | ||
431 | if (!fTreeK) { | |
432 | // Fatal("FinishEvent", "No kinematics tree is defined."); | |
433 | // Don't panic this is a probably a lego run | |
434 | return; | |
435 | ||
436 | } | |
437 | ||
438 | CleanParents(); | |
439 | if(fTreeK->GetEntries() ==0) { | |
440 | // set the fParticleFileMap size for the first time | |
441 | fParticleFileMap.Set(fHgwmk+1); | |
442 | } | |
443 | ||
444 | Bool_t allFilled = kFALSE; | |
445 | TObject *part; | |
446 | for(Int_t i=0; i<fHgwmk+1; ++i) | |
447 | if((part=fParticleMap->At(i))) { | |
448 | fParticleBuffer = (TParticle*) part; | |
449 | fParticleFileMap[i]= (Int_t) fTreeK->GetEntries(); | |
450 | fTreeK->Fill(); | |
451 | (*fParticleMap)[i]=0; | |
452 | ||
453 | // When all primaries were filled no particle!=0 | |
454 | // should be left => to be removed later. | |
455 | if (allFilled) printf("Why != 0 part # %d?\n",i); | |
456 | } | |
457 | else { | |
458 | // // printf("Why = 0 part # %d?\n",i); => We know. | |
459 | // break; | |
460 | // we don't break now in order to be sure there is no | |
461 | // particle !=0 left. | |
462 | // To be removed later and replaced with break. | |
463 | if(!allFilled) allFilled = kTRUE; | |
464 | } | |
465 | //cout << "Nof particles: " << fNtrack << endl; | |
466 | //Reset(); | |
467 | } | |
468 | ||
469 | //_____________________________________________________________________________ | |
470 | void AliStack::FlagTrack(Int_t track) | |
471 | { | |
472 | // | |
473 | // Flags a track and all its family tree to be kept | |
474 | // | |
475 | ||
476 | TParticle *particle; | |
477 | ||
478 | Int_t curr=track; | |
479 | while(1) { | |
480 | particle=(TParticle*)fParticleMap->At(curr); | |
481 | ||
482 | // If the particle is flagged the three from here upward is saved already | |
483 | if(particle->TestBit(kKeepBit)) return; | |
484 | ||
485 | // Save this particle | |
486 | particle->SetBit(kKeepBit); | |
487 | ||
488 | // Move to father if any | |
489 | if((curr=particle->GetFirstMother())==-1) return; | |
490 | } | |
491 | } | |
492 | ||
493 | //_____________________________________________________________________________ | |
ed54bf31 | 494 | void AliStack::KeepTrack(Int_t track) |
9e1a0ddb | 495 | { |
496 | // | |
497 | // Flags a track to be kept | |
498 | // | |
499 | ||
500 | fParticleMap->At(track)->SetBit(kKeepBit); | |
501 | } | |
502 | ||
503 | //_____________________________________________________________________________ | |
504 | void AliStack::Reset(Int_t size) | |
505 | { | |
506 | // | |
507 | // Resets stack | |
508 | // | |
509 | ||
510 | fNtrack=0; | |
511 | fNprimary=0; | |
512 | fHgwmk=0; | |
513 | fLoadPoint=0; | |
514 | fCurrent = -1; | |
515 | ResetArrays(size); | |
516 | } | |
517 | ||
518 | //_____________________________________________________________________________ | |
519 | void AliStack::ResetArrays(Int_t size) | |
520 | { | |
521 | // | |
522 | // Resets stack arrays | |
523 | // | |
524 | ||
525 | fParticles->Clear(); | |
526 | fParticleMap->Clear(); | |
527 | if (size>0) fParticleMap->Expand(size); | |
528 | } | |
529 | ||
530 | //_____________________________________________________________________________ | |
ed54bf31 | 531 | void AliStack::SetHighWaterMark(Int_t nt) |
9e1a0ddb | 532 | { |
533 | // | |
534 | // Set high water mark for last track in event | |
535 | // | |
536 | ||
537 | fHgwmk = fNtrack-1; | |
538 | fCurrentPrimary=fHgwmk; | |
539 | ||
540 | // Set also number of primary tracks | |
541 | fNprimary = fHgwmk+1; | |
542 | fNtrack = fHgwmk+1; | |
543 | } | |
544 | ||
545 | //_____________________________________________________________________________ | |
546 | TParticle* AliStack::Particle(Int_t i) | |
547 | { | |
548 | // | |
549 | // Return particle with specified ID | |
550 | ||
551 | if(!(*fParticleMap)[i]) { | |
552 | Int_t nentries = fParticles->GetEntries(); | |
553 | // algorithmic way of getting entry index | |
554 | // (primary particles are filled after secondaries) | |
555 | Int_t entry; | |
556 | if (i<fNprimary) | |
557 | entry = i+fNtrack-fNprimary; | |
558 | else | |
559 | entry = i-fNprimary; | |
560 | // check whether algorithmic way and | |
561 | // and the fParticleFileMap[i] give the same; | |
562 | // give the fatal error if not | |
563 | if (entry != fParticleFileMap[i]) { | |
564 | Fatal("Particle", | |
565 | "!! The algorithmic way and map are different: !!\n entry: %d map: %d", | |
566 | entry, fParticleFileMap[i]); | |
567 | } | |
568 | ||
569 | fTreeK->GetEntry(entry); | |
570 | new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); | |
571 | fParticleMap->AddAt((*fParticles)[nentries],i); | |
572 | } | |
573 | return (TParticle *) (*fParticleMap)[i]; | |
574 | } | |
575 | ||
576 | //_____________________________________________________________________________ | |
577 | Int_t AliStack::GetPrimary(Int_t id) const | |
578 | { | |
579 | // | |
580 | // Return number of primary that has generated track | |
581 | // | |
582 | ||
583 | int current, parent; | |
584 | TParticle *part; | |
585 | // | |
586 | parent=id; | |
587 | while (1) { | |
588 | current=parent; | |
589 | part = (TParticle *)fParticleMap->At(current); | |
590 | parent=part->GetFirstMother(); | |
591 | if(parent<0) return current; | |
592 | } | |
593 | } | |
594 | ||
595 | //_____________________________________________________________________________ | |
596 | void AliStack::DumpPart (Int_t i) const | |
597 | { | |
598 | // | |
599 | // Dumps particle i in the stack | |
600 | // | |
601 | ||
602 | ((TParticle*) (*fParticleMap)[i])->Print(); | |
603 | } | |
604 | ||
605 | //_____________________________________________________________________________ | |
606 | void AliStack::DumpPStack () | |
607 | { | |
608 | // | |
609 | // Dumps the particle stack | |
610 | // | |
611 | ||
3ab6f951 | 612 | Int_t i; |
613 | ||
9e1a0ddb | 614 | printf( |
615 | "\n\n=======================================================================\n"); | |
3ab6f951 | 616 | for (i=0;i<fNtrack;i++) |
9e1a0ddb | 617 | { |
618 | TParticle* particle = Particle(i); | |
619 | if (particle) { | |
620 | printf("-> %d ",i); particle->Print(); | |
621 | printf("--------------------------------------------------------------\n"); | |
622 | } | |
623 | else | |
624 | Warning("DumpPStack", "No particle with id %d.", i); | |
625 | } | |
626 | ||
627 | printf( | |
628 | "\n=======================================================================\n\n"); | |
629 | ||
630 | // print particle file map | |
631 | printf("\nParticle file map: \n"); | |
3ab6f951 | 632 | for (i=0; i<fNtrack; i++) |
9e1a0ddb | 633 | printf(" %d th entry: %d \n",i,fParticleFileMap[i]); |
634 | } | |
635 | ||
636 | ||
637 | //_____________________________________________________________________________ | |
638 | void AliStack::DumpLoadedStack() const | |
639 | { | |
640 | // | |
641 | // Dumps the particle in the stack | |
642 | // that are loaded in memory. | |
643 | // | |
644 | ||
645 | TObjArray &particles = *fParticleMap; | |
646 | printf( | |
647 | "\n\n=======================================================================\n"); | |
648 | for (Int_t i=0;i<fNtrack;i++) | |
649 | { | |
650 | TParticle* particle = (TParticle*) particles[i]; | |
651 | if (particle) { | |
652 | printf("-> %d ",i); particle->Print(); | |
653 | printf("--------------------------------------------------------------\n"); | |
654 | } | |
655 | else { | |
656 | printf("-> %d Particle not loaded.\n",i); | |
657 | printf("--------------------------------------------------------------\n"); | |
658 | } | |
659 | } | |
660 | printf( | |
661 | "\n=======================================================================\n\n"); | |
662 | } | |
663 | ||
664 | // | |
665 | // protected methods | |
666 | // | |
667 | ||
668 | //_____________________________________________________________________________ | |
669 | void AliStack::CleanParents() | |
670 | { | |
671 | // | |
672 | // Clean particles stack | |
673 | // Set parent/daughter relations | |
674 | // | |
675 | ||
676 | TObjArray &particles = *fParticleMap; | |
677 | TParticle *part; | |
678 | int i; | |
679 | for(i=0; i<fHgwmk+1; i++) { | |
680 | part = (TParticle *)particles.At(i); | |
681 | if(part) if(!part->TestBit(kDaughtersBit)) { | |
682 | part->SetFirstDaughter(-1); | |
683 | part->SetLastDaughter(-1); | |
684 | } | |
685 | } | |
686 | } | |
687 | ||
688 | //_____________________________________________________________________________ | |
689 | TParticle* AliStack::GetNextParticle() | |
690 | { | |
691 | // | |
692 | // Return next particle from stack of particles | |
693 | // | |
694 | ||
695 | TParticle* particle = 0; | |
696 | ||
697 | // search secondaries | |
698 | //for(Int_t i=fNtrack-1; i>=0; i--) { | |
699 | for(Int_t i=fNtrack-1; i>fHgwmk; i--) { | |
700 | particle = (TParticle*) fParticleMap->At(i); | |
701 | if ((particle) && (!particle->TestBit(kDoneBit))) { | |
702 | fCurrent=i; | |
703 | //cout << "GetNextParticle() - secondary " | |
704 | // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; | |
705 | return particle; | |
706 | } | |
707 | } | |
708 | ||
709 | // take next primary if all secondaries were done | |
710 | while (fCurrentPrimary>=0) { | |
711 | fCurrent = fCurrentPrimary; | |
712 | particle = (TParticle*) fParticleMap->At(fCurrentPrimary--); | |
713 | if ((particle) && (!particle->TestBit(kDoneBit))) { | |
714 | //cout << "GetNextParticle() - primary " | |
715 | // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; | |
716 | return particle; | |
717 | } | |
718 | } | |
719 | ||
720 | // nothing to be tracked | |
721 | fCurrent = -1; | |
722 | //cout << "GetNextParticle() - none " | |
723 | // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; | |
724 | return particle; | |
725 | } | |
726 | ||
727 | //__________________________________________________________________________________________ | |
728 | void AliStack::MakeTree(Int_t event, const char *file) | |
729 | { | |
730 | // | |
731 | // Make Kine tree and creates branch for writing particles | |
732 | // | |
733 | TBranch *branch=0; | |
734 | // Make Kinematics Tree | |
735 | char hname[30]; | |
736 | // printf("\n MakeTree called %d", event); | |
737 | if (!fTreeK) { | |
738 | sprintf(hname,"TreeK%d",event); | |
739 | fTreeK = new TTree(hname,"Kinematics"); | |
740 | // Create a branch for particles | |
741 | branch = fTreeK->Branch("Particles", "TParticle", &fParticleBuffer, 4000, 1); | |
742 | fTreeK->Write(0,TObject::kOverwrite); | |
743 | } | |
744 | } | |
745 | ||
81694e6a | 746 | //_____________________________________________________________________________ |
9e1a0ddb | 747 | void AliStack::BeginEvent(Int_t event) |
748 | { | |
749 | // start a new event | |
750 | // | |
751 | // | |
752 | fNprimary = 0; | |
753 | fNtrack = 0; | |
754 | ||
755 | char hname[30]; | |
756 | if(fTreeK) { | |
757 | fTreeK->Reset(); | |
758 | sprintf(hname,"TreeK%d",event); | |
759 | fTreeK->SetName(hname); | |
760 | } | |
761 | } | |
762 | ||
81694e6a | 763 | //_____________________________________________________________________________ |
9e1a0ddb | 764 | void AliStack::FinishRun() |
765 | { | |
766 | // Clean TreeK information | |
767 | if (fTreeK) { | |
768 | delete fTreeK; fTreeK = 0; | |
769 | } | |
770 | } | |
771 | ||
81694e6a | 772 | //_____________________________________________________________________________ |
9e1a0ddb | 773 | void AliStack::GetEvent(Int_t event) |
774 | { | |
775 | // | |
776 | // Get new event from TreeK | |
777 | // Reset/Create the particle stack | |
778 | if (fTreeK) delete fTreeK; | |
779 | ||
780 | // Get Kine Tree from file | |
781 | char treeName[20]; | |
782 | sprintf(treeName,"TreeK%d",event); | |
783 | fTreeK = (TTree*)gDirectory->Get(treeName); | |
784 | ||
785 | if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticleBuffer); | |
786 | ||
787 | else | |
788 | Error("GetEvent","cannot find Kine Tree for event:%d\n",event); | |
789 | // | |
790 | // printf("\n primaries %d", fNprimary); | |
791 | // printf("\n tracks %d", fNtrack); | |
792 | // | |
793 | Int_t size = (Int_t)fTreeK->GetEntries(); | |
794 | ResetArrays(size); | |
795 | } |