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