AliRoot Core  3dc7879 (3dc7879)
AliStack.cxx
Go to the documentation of this file.
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 /* $Id$ */
17 
19 // //
20 // Particles stack class //
21 // Implements the TMCVirtualStack of the Virtual Monte Carlo //
22 // Holds the particles transported during simulation //
23 // Is used to compare results of reconstruction with simulation //
24 // Author A.Morsch //
25 // //
27 
28 
29 #include <TClonesArray.h>
30 #include <TObjArray.h>
31 #include <TPDGCode.h>
32 #include <TMCProcess.h>
33 #include <TParticle.h>
34 #include <TParticlePDG.h>
35 #include <TDatabasePDG.h>
36 #include <TTree.h>
37 #include <TDirectory.h>
38 
39 #include "AliLog.h"
40 #include "AliStack.h"
41 #include "AliMiscConstants.h"
42 
43 
44 ClassImp(AliStack)
45 
46 
47 const Char_t* AliStack::fgkEmbedPathsKey = "embeddingBKGPaths";
48 
49 //_______________________________________________________________________
51  fParticles("TParticle", 1000),
52  fParticleMap(),
53  fParticleFileMap(0),
54  fParticleBuffer(0),
55  fCurrentTrack(0),
56  fTreeK(0),
57  fNtrack(0),
58  fNprimary(0),
59  fCurrent(-1),
60  fCurrentPrimary(-1),
61  fHgwmk(0),
62  fLoadPoint(0),
63  fTrackLabelMap(0),
64  fMCEmbeddingFlag(kFALSE)
65 {
66  //
67  // Default constructor
68  //
69 }
70 
71 //_______________________________________________________________________
72 AliStack::AliStack(Int_t size, const char* /*evfoldname*/):
73  fParticles("TParticle",1000),
74  fParticleMap(size),
75  fParticleFileMap(0),
76  fParticleBuffer(0),
77  fCurrentTrack(0),
78  fTreeK(0),
79  fNtrack(0),
80  fNprimary(0),
81  fNtransported(0),
82  fCurrent(-1),
83  fCurrentPrimary(-1),
84  fHgwmk(0),
85  fLoadPoint(0),
86  fTrackLabelMap(0),
87  fMCEmbeddingFlag(kFALSE)
88 {
89  //
90  // Constructor
91  //
92 }
93 
94 //_______________________________________________________________________
96  TVirtualMCStack(st),
97  fParticles("TParticle",1000),
98  fParticleMap(*(st.Particles())),
100  fParticleBuffer(0),
101  fCurrentTrack(0),
102  fTreeK((TTree*)(st.fTreeK->Clone())),
103  fNtrack(st.GetNtrack()),
104  fNprimary(st.GetNprimary()),
106  fCurrent(-1),
107  fCurrentPrimary(-1),
108  fHgwmk(0),
109  fLoadPoint(0),
110  fTrackLabelMap(0),
111  fMCEmbeddingFlag(kFALSE)
112 {
113  // Copy constructor
114 }
115 
116 
117 //_______________________________________________________________________
118 void AliStack::Copy(TObject&) const
119 {
120  AliFatal("Not implemented!");
121 }
122 
123 //_______________________________________________________________________
125 {
126  //
127  // Destructor
128  //
129 
130  fParticles.Clear();
131 }
132 
133 //
134 // public methods
135 //
136 
137 //_____________________________________________________________________________
138 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom,
139  const Float_t *vpos, const Float_t *polar, Float_t tof,
140  TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
141 {
142  //
143  // Load a track on the stack
144  //
145  // done 1 if the track has to be transported
146  // 0 if not
147  // parent identifier of the parent track. -1 for a primary
148  // pdg particle code
149  // pmom momentum GeV/c
150  // vpos position
151  // polar polarisation
152  // tof time of flight in seconds
153  // mecha production mechanism
154  // ntr on output the number of the track stored
155  //
156 
157  // const Float_t tlife=0;
158 
159  //
160  // Here we get the static mass
161  // For MC is ok, but a more sophisticated method could be necessary
162  // if the calculated mass is required
163  // also, this method is potentially dangerous if the mass
164  // used in the MC is not the same of the PDG database
165  //
166  TParticlePDG* pmc = TDatabasePDG::Instance()->GetParticle(pdg);
167  if (pmc) {
168  Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
169  Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
170  pmom[1]*pmom[1]+pmom[2]*pmom[2]);
171 
172 // 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",
173 // mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS);
174 
175 
176  PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e,
177  vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2],
178  mech, ntr, weight, is);
179  } else {
180  AliWarning(Form("Particle type %d not defined in PDG Database !", pdg));
181  AliWarning("Particle skipped !");
182  }
183 }
184 
185 //_____________________________________________________________________________
186 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg,
187  Double_t px, Double_t py, Double_t pz, Double_t e,
188  Double_t vx, Double_t vy, Double_t vz, Double_t tof,
189  Double_t polx, Double_t poly, Double_t polz,
190  TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is)
191 {
192  //
193  // Load a track on the stack
194  //
195  // done 1 if the track has to be transported
196  // 0 if not
197  // parent identifier of the parent track. -1 for a primary
198  // pdg particle code
199  // kS generation status code
200  // px, py, pz momentum GeV/c
201  // vx, vy, vz position
202  // polar polarisation
203  // tof time of flight in seconds
204  // mech production mechanism
205  // ntr on output the number of the track stored
206  //
207  // New method interface:
208  // arguments were changed to be in correspondence with TParticle
209  // constructor.
210  // Note: the energy is not calculated from the static mass but
211  // it is passed by argument e.
212 
213  const Int_t kFirstDaughter=-1;
214  const Int_t kLastDaughter=-1;
215 
216 
217  TParticle* particle
218  = new(fParticles[fLoadPoint++])
219  TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
220  px, py, pz, e, vx, vy, vz, tof);
221 
222  particle->SetPolarisation(polx, poly, polz);
223  particle->SetWeight(weight);
224  particle->SetUniqueID(mech);
225 
226 
227 
228  if(!done) {
229  particle->SetBit(kDoneBit);
230  } else {
231  particle->SetBit(kTransportBit);
232  fNtransported++;
233  }
234 
235 
236 
237  // Declare that the daughter information is valid
238  particle->SetBit(kDaughtersBit);
239  // Add the particle to the stack
240 
241  fParticleMap.AddAtAndExpand(particle, fNtrack);//CHECK!!
242 
243  if(parent>=0) {
244  particle = GetParticleMapEntry(parent);
245  if (particle) {
246  particle->SetLastDaughter(fNtrack);
247  if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
248  }
249  else {
250  AliError(Form("Parent %d does not exist",parent));
251  }
252  } else {
253  //
254  // This is a primary track. Set high water mark for this event
255  fHgwmk = fNtrack;
256  //
257  // Set also number if primary tracks
258  fNprimary = fHgwmk+1;
259  fCurrentPrimary++;
260  }
261  ntr = fNtrack++;
262 }
263 
264 //_____________________________________________________________________________
265 TParticle* AliStack::PopNextTrack(Int_t& itrack)
266 {
267  //
268  // Returns next track from stack of particles
269  //
270 
271 
272  TParticle* track = GetNextParticle();
273 
274  if (track) {
275  itrack = fCurrent;
276  track->SetBit(kDoneBit);
277  }
278  else
279  itrack = -1;
280 
282  return track;
283 }
284 
285 //_____________________________________________________________________________
287 {
288  //
289  // Returns i-th primary particle if it is flagged to be tracked,
290  // 0 otherwise
291  //
292 
293  TParticle* particle = Particle(i);
294 
295  if (!particle->TestBit(kDoneBit)) {
296  fCurrentTrack = particle;
297  return particle;
298  }
299  else
300  return 0;
301 }
302 
303 //_____________________________________________________________________________
304 Bool_t AliStack::PurifyKine(Float_t rmax, Float_t zmax)
305 {
306  //
307  // Compress kinematic tree keeping only flagged particles
308  // and renaming the particle id's in all the hits
309  //
310 
311  int nkeep = fHgwmk + 1, parent, i;
312  TParticle *part, *father;
313  fTrackLabelMap.Set(fParticleMap.GetLast()+1);
314 
315  // Save in Header total number of tracks before compression
316  // If no tracks generated return now
317  if(fHgwmk+1 == fNtrack) return kFALSE;
318 
319  // First pass, invalid Daughter information
320  for(i=0; i<fNtrack; i++) {
321  // Preset map, to be removed later
322  if(i<=fHgwmk) fTrackLabelMap[i]=i ;
323  else {
324  fTrackLabelMap[i] = -99;
325  if((part=GetParticleMapEntry(i))) {
326 //
327 // Check of this track should be kept for physics reasons
328  if (KeepPhysics(part, rmax, zmax)) KeepTrack(i);
329 //
330  part->ResetBit(kDaughtersBit);
331  part->SetFirstDaughter(-1);
332  part->SetLastDaughter(-1);
333  }
334  }
335  }
336  // Invalid daughter information for the parent of the first particle
337  // generated. This may or may not be the current primary according to
338  // whether decays have been recorded among the primaries
339  part = GetParticleMapEntry(fHgwmk+1);
340  fParticleMap.At(part->GetFirstMother())->ResetBit(kDaughtersBit);
341  // Second pass, build map between old and new numbering
342  for(i=fHgwmk+1; i<fNtrack; i++) {
343  if(fParticleMap.At(i)->TestBit(kKeepBit)) {
344  // This particle has to be kept
345  fTrackLabelMap[i]=nkeep;
346  // If old and new are different, have to move the pointer
347  if(i!=nkeep) fParticleMap[nkeep]=fParticleMap.At(i);
348  part = GetParticleMapEntry(nkeep);
349  // as the parent is always *before*, it must be already
350  // in place. This is what we are checking anyway!
351  if((parent=part->GetFirstMother())>fHgwmk) {
352  if(fTrackLabelMap[parent]==-99) Fatal("PurifyKine","fTrackLabelMap[%d] = -99!\n",parent);
353  else part->SetFirstMother(fTrackLabelMap[parent]);}
354  nkeep++;
355  }
356  }
357 
358  // Fix daughters information
359  for (i=fHgwmk+1; i<nkeep; i++) {
360  part = GetParticleMapEntry(i);
361  parent = part->GetFirstMother();
362  if(parent>=0) {
363  father = GetParticleMapEntry(parent);
364  if(father->TestBit(kDaughtersBit)) {
365 
366  if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
367  if(i>father->GetLastDaughter()) father->SetLastDaughter(i);
368  } else {
369  // Initialise daughters info for first pass
370  father->SetFirstDaughter(i);
371  father->SetLastDaughter(i);
372  father->SetBit(kDaughtersBit);
373  }
374  }
375  }
376  //
377  // Now the output bit, from fHgwmk to nkeep we write everything and we erase
378  if(nkeep > fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
379  for (i=fHgwmk+1; i<nkeep; ++i) {
381  fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries());
382  TreeK()->Fill();
384  }
385 
386  for (i = nkeep; i < fNtrack; ++i) fParticleMap[i]=0;
387 
388  Int_t toshrink = fNtrack-fHgwmk-1;
389  fLoadPoint-=toshrink;
390 
391  for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles.RemoveAt(i);
392  fNtrack=nkeep;
393  fHgwmk=nkeep-1;
394  return kTRUE;
395 }
396 
397 
399 {
400 //
401 // In some transport code children might not come in a continuous sequence.
402 // In this case the stack has to be reordered in order to establish the
403 // mother daughter relation using index ranges.
404 //
405  if(fHgwmk+1 == fNtrack) return kFALSE;
406 
407  //
408  // Howmany secondaries have been produced ?
409  Int_t nNew = fNtrack - fHgwmk - 1;
410 
411  if (nNew > 0) {
412  Int_t i, j;
413  TArrayI map1(nNew);
414  //
415  // Copy pointers to temporary array
416  TParticle** tmp = new TParticle*[nNew];
417 
418  for (i = 0; i < nNew; i++) {
419  if (fParticleMap.At(fHgwmk + 1 + i)) {
420  tmp[i] = GetParticleMapEntry(fHgwmk + 1 + i);
421  } else {
422  tmp[i] = 0x0;
423  }
424  map1[i] = -99;
425  }
426 
427 
428  //
429  // Reset LoadPoint
430  //
431  Int_t loadPoint = fHgwmk + 1;
432  //
433  // Re-Push particles into stack
434  // The outer loop is over parents, the inner over children.
435  // -1 refers to the primary particle
436  //
437  for (i = -1; i < nNew-1; i++) {
438  Int_t ipa;
439  TParticle* parP;
440  if (i == -1) {
441  ipa = tmp[0]->GetFirstMother();
442  parP = GetParticleMapEntry(ipa);
443  } else {
444  ipa = (fHgwmk + 1 + i);
445  // Skip deleted particles
446  if (!tmp[i]) continue;
447  // Skip particles without children
448  if (tmp[i]->GetFirstDaughter() == -1) continue;
449  parP = tmp[i];
450  }
451  // Reset daughter information
452 
453  Int_t idaumin = parP->GetFirstDaughter() - fHgwmk - 1;
454  Int_t idaumax = parP->GetLastDaughter() - fHgwmk - 1;
455  parP->SetFirstDaughter(-1);
456  parP->SetLastDaughter(-1);
457  for (j = idaumin; j <= idaumax; j++) {
458  // Skip deleted particles
459  if (!tmp[j]) continue;
460  // Skip particles already handled
461  if (map1[j] != -99) continue;
462  Int_t jpa = tmp[j]->GetFirstMother();
463  // Check if daughter of current parent
464  if (jpa == ipa) {
465  fParticleMap[loadPoint] = tmp[j];
466  // Re-establish daughter information
467  parP->SetLastDaughter(loadPoint);
468  if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(loadPoint);
469  // Set Mother information
470  if (i != -1) {
471  tmp[j]->SetFirstMother(map1[i]);
472  }
473  // Build the map
474  map1[j] = loadPoint;
475  // Increase load point
476  loadPoint++;
477  }
478  } // children
479  } // parents
480 
481  delete[] tmp;
482 
483  //
484  // Build map for remapping of hits
485  //
486  fTrackLabelMap.Set(fNtrack);
487  for (i = 0; i < fNtrack; i ++) {
488  if (i <= fHgwmk) {
489  fTrackLabelMap[i] = i;
490  } else{
491  fTrackLabelMap[i] = map1[i - fHgwmk -1];
492  }
493  }
494  } // new particles poduced
495 
496  return kTRUE;
497 }
498 
499 Bool_t AliStack::KeepPhysics(const TParticle* part, Float_t rmax, Float_t zmax)
500 {
501  //
502  // Some particles have to kept on the stack for reasons motivated
503  // by physics analysis. Decision is put here.
504  //
505  Bool_t keep = kFALSE;
506 
507 
508 
509  Int_t parent = part->GetFirstMother();
510  if (parent >= 0 && parent <= fHgwmk) {
511  // Keep 1st generation secondaries in a pre-defined r-z range
512  Float_t vx = part->Vx();
513  Float_t vy = part->Vy();
514  Float_t vz = part->Vz();
515  Float_t r = TMath::Sqrt(vx * vx + vy * vy);
516  if (r < rmax && TMath::Abs(vz) < zmax) return kTRUE;
517  //
518  TParticle* father = GetParticleMapEntry(parent);
519  //
520  // Keep first-generation daughter from primaries with heavy flavor
521  //
522  Int_t kf = father->GetPdgCode();
523  kf = TMath::Abs(kf);
524  Int_t kfl = kf;
525  // meson ?
526  if (kfl > 10) kfl/=100;
527  // baryon
528  if (kfl > 10) kfl/=10;
529  if (kfl > 10) kfl/=10;
530  if (kfl >= 4) {
531  keep = kTRUE;
532  }
533  //
534  // e+e- from pair production of primary gammas
535  //
536  if ((part->GetUniqueID()) == kPPair) keep = kTRUE;
537  }
538  //
539  // Decay(cascade) from primaries
540  //
541  if ((part->GetUniqueID() == kPDecay) && (parent >= 0)) {
542  // Particles from decay
543  TParticle* father = GetParticleMapEntry(parent);
544  Int_t imo = parent;
545  while((imo > fHgwmk) && (father->GetUniqueID() == kPDecay)) {
546  imo = father->GetFirstMother();
547  father = GetParticleMapEntry(imo);
548  }
549  if ((imo <= fHgwmk)) keep = kTRUE;
550  }
551  return keep;
552 }
553 
554 //_____________________________________________________________________________
556 {
557 //
558 // Write out the kinematics that was not yet filled
559 //
560 
561 // Update event header
562 
563  if (!TreeK()) {
564 // Fatal("FinishEvent", "No kinematics tree is defined.");
565 // Don't panic this is a probably a lego run
566  return;
567  }
568 
569  CleanParents();
570  if(TreeK()->GetEntries() ==0) {
571  // set the fParticleFileMap size for the first time
572  fParticleFileMap.Set(fHgwmk+1);
573  }
574 
575  Bool_t allFilled = kFALSE;
576  TParticle *part;
577  for(Int_t i=0; i<fHgwmk+1; ++i) {
578  if((part=GetParticleMapEntry(i))) {
579  fParticleBuffer = part;
580  fParticleFileMap[i]= static_cast<Int_t>(TreeK()->GetEntries());
581  TreeK()->Fill();
582  fParticleBuffer=0;
583  fParticleMap.AddAt(0,i);
584 
585  // When all primaries were filled no particle!=0
586  // should be left => to be removed later.
587  if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i));
588  }
589  else
590  {
591  // // printf("Why = 0 part # %d?\n",i); => We know.
592  // break;
593  // we don't break now in order to be sure there is no
594  // particle !=0 left.
595  // To be removed later and replaced with break.
596  if(!allFilled) allFilled = kTRUE;
597  }
598  }
599  AliInfoF("Ntrack=%d kept from %d transported\n",fNtrack,fNtransported);
600 }
601 //_____________________________________________________________________________
602 
604 {
605  //
606  // Flags a track and all its family tree to be kept
607  //
608 
609  TParticle *particle;
610 
611  Int_t curr=track;
612  while(1) {
613  particle = GetParticleMapEntry(curr);
614 
615  // If the particle is flagged the three from here upward is saved already
616  if(particle->TestBit(kKeepBit)) return;
617 
618  // Save this particle
619  particle->SetBit(kKeepBit);
620 
621  // Move to father if any
622  if((curr=particle->GetFirstMother())==-1) return;
623  }
624 }
625 
626 //_____________________________________________________________________________
628 {
629  //
630  // Flags a track to be kept
631  //
632 
633  fParticleMap.At(track)->SetBit(kKeepBit);
634 }
635 
636 //_____________________________________________________________________________
637 void AliStack::Clean(Int_t size)
638 {
639  //
640  // Reset stack data except for fTreeK
641  //
642 
643  fNtrack=0;
644  fNprimary=0;
645  fNtransported=0;
646  fHgwmk=0;
647  fLoadPoint=0;
648  fCurrent = -1;
649  ResetArrays(size);
650 }
651 
652 //_____________________________________________________________________________
653 void AliStack::Reset(Int_t size)
654 {
655  //
656  // Reset stack data including fTreeK
657  //
658 
659  Clean(size);
660  delete fParticleBuffer; fParticleBuffer = 0;
661  fTreeK = 0x0;
662 }
663 
664 //_____________________________________________________________________________
665 void AliStack::ResetArrays(Int_t size)
666 {
667  //
668  // Resets stack arrays
669  //
670  fParticles.Clear();
671  fParticleMap.Clear();
672  if (size>0) fParticleMap.Expand(size);
673 }
674 
675 //_____________________________________________________________________________
677 {
678  //
679  // Set high water mark for last track in event
680  //
681 
682  fHgwmk = fNtrack-1;
684  // Set also number of primary tracks
685  fNprimary = fHgwmk+1;
686 }
687 
688 //_____________________________________________________________________________
689 TParticle* AliStack::Particle(Int_t i, Bool_t useInEmbedding)
690 {
691  //
692  // Return particle with specified ID
693  if (GetMCEmbeddingFlag() && !useInEmbedding) {
694  AliError("Method should not be called by user in embedding mode, returning dummy particle");
695  return GetDummyParticle();
696  }
697  if (i==gkDummyLabel) return 0;
698 
699  if(!fParticleMap.At(i)) {
700  Int_t nentries = fParticles.GetEntriesFast();
701  // algorithmic way of getting entry index
702  // (primary particles are filled after secondaries)
703  Int_t entry = TreeKEntry(i,useInEmbedding);
704  // check whether algorithmic way and
705  // and the fParticleFileMap[i] give the same;
706  // give the fatal error if not
707  if (entry != fParticleFileMap[i]) {
708  AliFatal(Form(
709  "!! The algorithmic way and map are different: !!\n entry: %d map: %d",
710  entry, fParticleFileMap[i]));
711  }
712  // Load particle at entry into fParticleBuffer
713  TreeK()->GetEntry(entry);
714  // Add to the TClonesarray
715  new (fParticles[nentries]) TParticle(*fParticleBuffer);
716  // Store a pointer in the TObjArray
717  fParticleMap.AddAt(fParticles[nentries],i);
718  }
719  return GetParticleMapEntry(i);
720 }
721 
722 //_____________________________________________________________________________
723 TParticle* AliStack::ParticleFromTreeK(Int_t id, Bool_t useInEmbedding) const
724 {
725 //
726 // return pointer to TParticle with label id
727 //
728  if (GetMCEmbeddingFlag() && !useInEmbedding) {
729  AliError("Method should not be called by user in embedding mode, returning dummy particle");
730  return GetDummyParticle();
731  }
732  Int_t entry;
733  if ((entry = TreeKEntry(id,useInEmbedding)) < 0) return 0;
734  if (fTreeK->GetEntry(entry)<=0) return 0;
735  return fParticleBuffer;
736 }
737 
738 //_____________________________________________________________________________
739 Int_t AliStack::TreeKEntry(Int_t id, Bool_t useInEmbedding) const
740 {
741 //
742 // Return entry number in the TreeK for particle with label id
743 // Return negative number if label>fNtrack
744 //
745 // The order of particles in TreeK reflects the order of the transport of primaries and production of secondaries:
746 //
747 // Before transport there are fNprimary particles on the stack.
748 // They are transported one by one and secondaries (fNtrack - fNprimary) are produced.
749 // After the transport of each particles secondaries are written to the TreeK
750 // They occupy the entries 0 ... fNtrack - fNprimary - 1
751 // The primaries are written after they have been transported and occupy
752 // fNtrack - fNprimary .. fNtrack - 1
753 
754  if (GetMCEmbeddingFlag() && !useInEmbedding) {
755  AliError("Method should not be called by user in embedding mode, returning -1");
756  return -1;
757  }
758 
759 
760  Int_t entry;
761  if (id<fNprimary)
762  entry = id+fNtrack-fNprimary;
763  else
764  entry = id-fNprimary;
765  return entry;
766 }
767 
768 //_____________________________________________________________________________
770 {
771  //
772  // Return number of the parent of the current track
773  //
774 
775  TParticle* current = GetParticleMapEntry(fCurrent);
776 
777  if (current)
778  return current->GetFirstMother();
779  else {
780  AliWarning("Current track not found in the stack");
781  return -1;
782  }
783 }
784 
785 //_____________________________________________________________________________
786 Int_t AliStack::GetPrimary(Int_t id, Bool_t useInEmbedding)
787 {
788  //
789  // Return number of primary that has generated track
790  //
791 
792  int current, parent;
793  //
794  parent=id;
795  while (1) {
796  current=parent;
797  TParticle* part = Particle(current,useInEmbedding);
798  if (!part || ( parent=part->GetFirstMother() )<0 ) return current;
799  }
800 }
801 
802 //_____________________________________________________________________________
803 void AliStack::DumpPart (Int_t i) const
804 {
805  //
806  // Dumps particle i in the stack
807  //
808  GetParticleMapEntry(i)->Print();
809 }
810 
811 //_____________________________________________________________________________
813 {
814  //
815  // Dumps the particle stack
816  //
817 
818  Int_t i;
819 
820  printf("\n\n=======================================================================\n");
821  for (i=0;i<fNtrack;i++)
822  {
823  TParticle* particle = Particle(i);
824  if (particle) {
825  printf("-> %d ",i); particle->Print();
826  printf("--------------------------------------------------------------\n");
827  }
828  else
829  Warning("DumpPStack", "No particle with id %d.", i);
830  }
831 
832  printf("\n=======================================================================\n\n");
833 
834  // print particle file map
835  // printf("\nParticle file map: \n");
836  // for (i=0; i<fNtrack; i++)
837  // printf(" %d th entry: %d \n",i,fParticleFileMap[i]);
838 }
839 
840 
841 //_____________________________________________________________________________
843 {
844  //
845  // Dumps the particle in the stack
846  // that are loaded in memory.
847  //
848 
849  printf(
850  "\n\n=======================================================================\n");
851  for (Int_t i=0;i<fNtrack;i++)
852  {
853  TParticle* particle = GetParticleMapEntry(i);
854  if (particle) {
855  printf("-> %d ",i); particle->Print();
856  printf("--------------------------------------------------------------\n");
857  }
858  else {
859  printf("-> %d Particle not loaded.\n",i);
860  printf("--------------------------------------------------------------\n");
861  }
862  }
863  printf(
864  "\n=======================================================================\n\n");
865 }
866 
867 //_____________________________________________________________________________
869 {
870  fCurrent = track;
871  if (fCurrent < fNprimary) fCurrentTrack = Particle(track);
872 }
873 
874 
875 //_____________________________________________________________________________
876 //
877 // protected methods
878 //
879 
880 //_____________________________________________________________________________
882 {
883  //
884  // Clean particles stack
885  // Set parent/daughter relations
886  //
887 
888  TParticle *part;
889  int i;
890  for(i=0; i<fHgwmk+1; i++) {
891  part = GetParticleMapEntry(i);
892  if(part) if(!part->TestBit(kDaughtersBit)) {
893  part->SetFirstDaughter(-1);
894  part->SetLastDaughter(-1);
895  }
896  }
897 }
898 
899 //_____________________________________________________________________________
901 {
902  //
903  // Return next particle from stack of particles
904  //
905 
906  TParticle* particle = 0;
907 
908  // search secondaries
909  //for(Int_t i=fNtrack-1; i>=0; i--) {
910  for(Int_t i=fNtrack-1; i>fHgwmk; i--) {
911  particle = GetParticleMapEntry(i);
912  if ((particle) && (!particle->TestBit(kDoneBit))) {
913  fCurrent=i;
914  return particle;
915  }
916  }
917 
918  // take next primary if all secondaries were done
919  while (fCurrentPrimary>=0) {
921  particle = GetParticleMapEntry(fCurrentPrimary--);
922  if ((particle) && (!particle->TestBit(kDoneBit))) {
923  return particle;
924  }
925  }
926 
927  // nothing to be tracked
928  fCurrent = -1;
929 
930 
931  return particle;
932 }
933 //__________________________________________________________________________________________
934 
936 {
937 //
938 // Creates branch for writing particles
939 //
940 
941  fTreeK = tree;
942 
943  AliDebug(1, "Connecting TreeK");
944  if (fTreeK == 0x0)
945  {
946  if (TreeK() == 0x0)
947  {
948  AliFatal("Parameter is NULL");//we don't like such a jokes
949  return;
950  }
951  return;//in this case TreeK() calls back this method (ConnectTree)
952  //tree after setting fTreeK, the rest was already executed
953  //it is safe to return now
954  }
955 
956  // Create a branch for particles
957 
958  AliDebug(2, Form("Tree name is %s",fTreeK->GetName()));
959 
960  if (fTreeK->GetDirectory())
961  {
962  AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName()));
963  }
964  else
965  AliWarning("DIR IS NOT SET !!!");
966 
967  TBranch *branch=fTreeK->GetBranch("Particles");
968  if(branch == 0x0)
969  {
970  branch = fTreeK->Branch("Particles", &fParticleBuffer, 4000);
971  AliDebug(2, "Creating Branch in Tree");
972  }
973  else
974  {
975  AliDebug(2, "Branch Found in Tree");
976  branch->SetAddress(&fParticleBuffer);
977  }
978  if (branch->GetDirectory())
979  {
980  AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName()));
981  }
982  else
983  AliWarning("Branch Dir is NOT SET");
984 }
985 
986 //_____________________________________________________________________________
987 
989 {
990 //
991 // Get new event from TreeK
992 
993  // Reset/Create the particle stack
994  Int_t size = (Int_t)TreeK()->GetEntries();
995  ResetArrays(size);
996  return kTRUE;
997 }
998 //_____________________________________________________________________________
999 
1000 Bool_t AliStack::IsStable(Int_t pdg) const
1001 {
1002  //
1003  // Decide whether particle (pdg) is stable
1004  //
1005 
1006 
1007  // All ions/nucleons are considered as stable
1008  // Nuclear code is 10LZZZAAAI
1009  if(pdg>1000000000)return kTRUE;
1010 
1011  const Int_t kNstable = 18;
1012  Int_t i;
1013 
1014  Int_t pdgStable[kNstable] = {
1015  kGamma, // Photon
1016  kElectron, // Electron
1017  kMuonPlus, // Muon
1018  kPiPlus, // Pion
1019  kKPlus, // Kaon
1020  kK0Short, // K0s
1021  kK0Long, // K0l
1022  kProton, // Proton
1023  kNeutron, // Neutron
1024  kLambda0, // Lambda_0
1025  kSigmaMinus, // Sigma Minus
1026  kSigmaPlus, // Sigma Plus
1027  3312, // Xsi Minus
1028  3322, // Xsi
1029  3334, // Omega
1030  kNuE, // Electron Neutrino
1031  kNuMu, // Muon Neutrino
1032  kNuTau // Tau Neutrino
1033  };
1034 
1035  Bool_t isStable = kFALSE;
1036  for (i = 0; i < kNstable; i++) {
1037  if (pdg == TMath::Abs(pdgStable[i])) {
1038  isStable = kTRUE;
1039  break;
1040  }
1041  }
1042 
1043  return isStable;
1044 }
1045 
1046 //_____________________________________________________________________________
1047 Bool_t AliStack::IsPhysicalPrimary(Int_t index, Bool_t useInEmbedding)
1048 {
1049  //
1050  // Test if a particle is a physical primary according to the following definition:
1051  // Particles produced in the collision including products of strong and
1052  // electromagnetic decay and excluding feed-down from weak decays of strange
1053  // particles.
1054  //
1055  TParticle* p = Particle(index,useInEmbedding);
1056  if (!p) return kFALSE;
1057  Int_t ist = p->GetStatusCode();
1058  Int_t pdg = TMath::Abs(p->GetPdgCode());
1059  //
1060  // Initial state particle
1061  // Solution for K0L decayed by Pythia6
1062  // ->
1063  if ((ist > 1) && (pdg!=130) && index < GetNprimary()) return kFALSE;
1064  if ((ist > 1) && index >= GetNprimary()) return kFALSE;
1065  // <-
1066 
1067 
1068  if (!IsStable(pdg)) return kFALSE;
1069  if (index < GetNprimary()) {
1070 //
1071 // Particle produced by generator
1072  // Solution for K0L decayed by Pythia6
1073  // ->
1074  Int_t ipm = p->GetFirstMother();
1075  if (ipm > -1) {
1076  TParticle* ppm = Particle(ipm, useInEmbedding);
1077  if (TMath::Abs(ppm->GetPdgCode()) == 130) return kFALSE;
1078  }
1079  // <-
1080  return kTRUE;
1081  } else {
1082 //
1083 // Particle produced during transport
1084 //
1085 
1086  Int_t imo = p->GetFirstMother();
1087  TParticle* pm = Particle(imo,useInEmbedding);
1088  Int_t mpdg = TMath::Abs(pm->GetPdgCode());
1089 // Check for Sigma0
1090  if ((mpdg == 3212) && (imo < GetNprimary())) return kTRUE;
1091 //
1092 // Check if it comes from a pi0 decay
1093 //
1094  if ((mpdg == kPi0) && (imo < GetNprimary())) return kTRUE;
1095 
1096 // Check if this is a heavy flavor decay product
1097  Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
1098  //
1099  // Light hadron
1100  if (mfl < 4) return kFALSE;
1101 
1102  //
1103  // Heavy flavor hadron produced by generator
1104  if (imo < GetNprimary()) {
1105  return kTRUE;
1106  }
1107 
1108  // To be sure that heavy flavor has not been produced in a secondary interaction
1109  // Loop back to the generated mother
1110  while (imo >= GetNprimary()) {
1111  imo = pm->GetFirstMother();
1112  pm = Particle(imo,useInEmbedding);
1113  }
1114  mpdg = TMath::Abs(pm->GetPdgCode());
1115  mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
1116 
1117  if (mfl < 4) {
1118  return kFALSE;
1119  } else {
1120  return kTRUE;
1121  }
1122  } // produced by generator ?
1123 }
1124 
1125 Bool_t AliStack::IsSecondaryFromWeakDecay(Int_t index, Bool_t useInEmbedding) {
1126 
1127  // If a particle is not a physical primary, check if it comes from weak decay
1128 
1129  if(IsPhysicalPrimary(index,useInEmbedding)) return kFALSE;
1130 
1131  TParticle* particle = Particle(index, useInEmbedding);
1132  if (!particle) return kFALSE;
1133  Int_t uniqueID = particle->GetUniqueID();
1134 
1135  Int_t indexMoth = particle->GetFirstMother();
1136  if(indexMoth < 0) return kFALSE; // if index mother < 0 and not a physical primary, is a non-stable product or one of the beams
1137  TParticle* moth = Particle(indexMoth,useInEmbedding);
1138  Float_t codemoth = (Float_t)TMath::Abs(moth->GetPdgCode());
1139  // mass of the flavour
1140  Int_t mfl = 0;
1141  // Protect the "rootino" case when codemoth is 0
1142  if (TMath::Abs(codemoth)>0) mfl = Int_t (codemoth / TMath::Power(10, Int_t(TMath::Log10(codemoth))));
1143 
1144  if(mfl == 3 && uniqueID == kPDecay) return kTRUE;// The first mother is strange and it's a decay
1145  if(codemoth == 211 && uniqueID == kPDecay) return kTRUE;// pion+- decay products
1146  if(codemoth == 13 && uniqueID == kPDecay) return kTRUE;// muon decay products
1147 
1149  if (TMath::Abs(moth->GetPdgCode()) > 1000000000 && uniqueID == kPDecay) {
1150  if ((moth->GetPdgCode() / 10000000) % 10 != 0) return kTRUE;
1151  }
1152 
1153  return kFALSE;
1154 
1155 }
1156 Bool_t AliStack::IsSecondaryFromMaterial(Int_t index, Bool_t useInEmbedding) {
1157 
1158  // If a particle is not a physical primary, check if it comes from material
1159 
1160  if(IsPhysicalPrimary(index,useInEmbedding)) return kFALSE;
1161  if(IsSecondaryFromWeakDecay(index,useInEmbedding)) return kFALSE;
1162  TParticle* particle = Particle(index,useInEmbedding);
1163  if (!particle) return kFALSE;
1164  Int_t indexMoth = particle->GetFirstMother();
1165  if(indexMoth < 0) return kFALSE; // if index mother < 0 and not a physical primary, is a non-stable product or one of the beams
1166  return kTRUE;
1167 
1168 }
1169 
1170 //__________________________________________
1172 {
1173  static TParticle dummy(21,999,-1,-1,-1,-1,1,1,999,999,0,0,0,0);
1174  return &dummy;
1175 }
virtual Int_t GetCurrentParentTrackNumber() const
Definition: AliStack.cxx:769
void SetHighWaterMark(Int_t hgwmk)
Definition: AliStack.cxx:676
TArrayI fParticleFileMap
Map of particles in the supporting TClonesArray.
Definition: AliStack.h:110
Int_t fNtrack
Particle stack.
Definition: AliStack.h:114
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
Bool_t GetEvent()
Definition: AliStack.cxx:988
Bool_t KeepPhysics(const TParticle *part, Float_t rmax=-1, Float_t zmax=-1.)
Definition: AliStack.cxx:499
static TParticle * GetDummyParticle()
Definition: AliStack.cxx:1171
TArrayI fTrackLabelMap
Next free position in the particle buffer.
Definition: AliStack.h:121
TParticle * fCurrentTrack
Pointer to current particle for writing.
Definition: AliStack.h:112
AliStack()
Definition: AliStack.cxx:50
void CleanParents()
Definition: AliStack.cxx:881
Int_t fNtransported
Definition: AliStack.h:116
Int_t fCurrent
Definition: AliStack.h:117
Bool_t IsSecondaryFromWeakDecay(Int_t index, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:1125
virtual TParticle * PopPrimaryForTracking(Int_t i)
Definition: AliStack.cxx:286
Float_t p[]
Definition: kNNTest.C:133
Int_t fHgwmk
Last primary track returned from the stack.
Definition: AliStack.h:119
AliTPCfastTrack * track
#define AliInfoF(message,...)
Definition: AliLog.h:499
void DumpPStack()
Definition: AliStack.cxx:812
TVectorD vz
Definition: driftITSTPC.C:87
void ResetArrays(Int_t size)
Definition: AliStack.cxx:665
#define AliWarning(message)
Definition: AliLog.h:541
Bool_t ReorderKine()
Definition: AliStack.cxx:398
TTree * fTreeK
Pointer to particle currently transported.
Definition: AliStack.h:113
Int_t GetNprimary() const
Definition: AliStack.h:137
Bool_t fMCEmbeddingFlag
Map of track labels.
Definition: AliStack.h:122
TTree * tree
Bool_t IsPhysicalPrimary(Int_t i, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:1047
virtual Int_t GetNtrack() const
Definition: AliStack.h:134
virtual ~AliStack()
Definition: AliStack.cxx:124
Int_t TreeKEntry(Int_t id, Bool_t useInEmbedding=kFALSE) const
Definition: AliStack.cxx:739
void ConnectTree(TTree *tree)
Definition: AliStack.cxx:935
void Copy(TObject &st) const
Definition: AliStack.cxx:118
virtual void SetCurrentTrack(Int_t track)
Definition: AliStack.cxx:868
virtual TParticle * PopNextTrack(Int_t &track)
Definition: AliStack.cxx:265
void polar()
Definition: polar.C:2
Int_t GetPrimary(Int_t id, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:786
Int_t GetNtransported() const
Definition: AliStack.h:140
const TObjArray * Particles() const
Definition: AliStack.h:146
void KeepTrack(Int_t itrack)
Definition: AliStack.cxx:627
#define AliFatal(message)
Definition: AliLog.h:640
TParticle * fParticleBuffer
Definition: AliStack.h:111
Bool_t IsSecondaryFromMaterial(Int_t index, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:1156
void Clean(Int_t size=0)
Definition: AliStack.cxx:637
virtual void PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom, const Float_t *vpos, const Float_t *polar, Float_t tof, TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
Definition: AliStack.cxx:138
void DumpLoadedStack() const
Definition: AliStack.cxx:842
Int_t fNprimary
Definition: AliStack.h:115
TParticle * Particle(Int_t id, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:689
#define AliDebug(logLevel, message)
Definition: AliLog.h:300
Bool_t IsStable(Int_t pdg) const
Definition: AliStack.cxx:1000
void Reset(Int_t size=0)
Definition: AliStack.cxx:653
TParticle * GetParticleMapEntry(Int_t id) const
Definition: AliStack.h:151
void FinishEvent()
Definition: AliStack.cxx:555
Int_t fLoadPoint
Last track purified.
Definition: AliStack.h:120
Int_t fCurrentPrimary
Last track returned from the stack.
Definition: AliStack.h:118
#define AliError(message)
Definition: AliLog.h:591
Bool_t PurifyKine(Float_t rmax=-1., Float_t zmax=-1.)
Definition: AliStack.cxx:304
void FlagTrack(Int_t track)
Definition: AliStack.cxx:603
const int gkDummyLabel
TObjArray fParticleMap
Pointer to list of particles.
Definition: AliStack.h:109
void DumpPart(Int_t i) const
Definition: AliStack.cxx:803
TClonesArray fParticles
Definition: AliStack.h:108
TParticle * ParticleFromTreeK(Int_t id, Bool_t useInEmbedding=kFALSE) const
Definition: AliStack.cxx:723
TParticle * GetNextParticle()
Definition: AliStack.cxx:900
TTree * TreeK() const
Definition: AliStack.h:81
Bool_t GetMCEmbeddingFlag() const
Definition: AliStack.h:92