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