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