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