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