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