415d9f5c |
1 | /************************************************************************** |
2 | * Copyright(c) 1998-2007, 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 | // Class for Kinematic Events |
20 | // Author: Andreas Morsch, CERN |
21 | //------------------------------------------------------------------------- |
22 | #include <TArrow.h> |
23 | #include <TMarker.h> |
24 | #include <TH2F.h> |
25 | #include <TTree.h> |
26 | #include <TFile.h> |
27 | #include <TParticle.h> |
28 | #include <TClonesArray.h> |
93df0e9b |
29 | #include <TRefArray.h> |
93836e1b |
30 | #include <TList.h> |
5df3496b |
31 | #include <TArrayF.h> |
415d9f5c |
32 | |
33 | #include "AliLog.h" |
34 | #include "AliMCEvent.h" |
5df3496b |
35 | #include "AliMCVertex.h" |
415d9f5c |
36 | #include "AliStack.h" |
37 | #include "AliTrackReference.h" |
38 | #include "AliHeader.h" |
83626b6b |
39 | #include "AliGenEventHeader.h" |
415d9f5c |
40 | |
41 | |
93836e1b |
42 | Int_t AliMCEvent::fgkBgLabelOffset(10000000); |
43 | |
44 | |
415d9f5c |
45 | AliMCEvent::AliMCEvent(): |
46 | AliVEvent(), |
47 | fStack(0), |
7aad0c47 |
48 | fMCParticles(0), |
415d9f5c |
49 | fMCParticleMap(0), |
fb9aeda5 |
50 | fHeader(new AliHeader()), |
93df0e9b |
51 | fTRBuffer(0), |
52 | fTrackReferences(new TClonesArray("AliTrackReference", 1000)), |
415d9f5c |
53 | fTreeTR(0), |
54 | fTmpTreeTR(0), |
55 | fTmpFileTR(0), |
56 | fNprimaries(-1), |
93836e1b |
57 | fNparticles(-1), |
58 | fSubsidiaryEvents(0), |
59 | fPrimaryOffset(0), |
7aad0c47 |
60 | fSecondaryOffset(0), |
5df3496b |
61 | fExternal(0), |
62 | fVertex(0) |
415d9f5c |
63 | { |
64 | // Default constructor |
65 | } |
66 | |
67 | AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) : |
68 | AliVEvent(mcEvnt), |
93836e1b |
69 | fStack(mcEvnt.fStack), |
70 | fMCParticles(mcEvnt.fMCParticles), |
71 | fMCParticleMap(mcEvnt.fMCParticleMap), |
72 | fHeader(mcEvnt.fHeader), |
73 | fTRBuffer(mcEvnt.fTRBuffer), |
74 | fTrackReferences(mcEvnt.fTrackReferences), |
75 | fTreeTR(mcEvnt.fTreeTR), |
76 | fTmpTreeTR(mcEvnt.fTmpTreeTR), |
77 | fTmpFileTR(mcEvnt.fTmpFileTR), |
78 | fNprimaries(mcEvnt.fNprimaries), |
79 | fNparticles(mcEvnt.fNparticles), |
80 | fSubsidiaryEvents(0), |
81 | fPrimaryOffset(0), |
7aad0c47 |
82 | fSecondaryOffset(0), |
5df3496b |
83 | fExternal(0), |
84 | fVertex(mcEvnt.fVertex) |
415d9f5c |
85 | { |
86 | // Copy constructor |
87 | } |
88 | |
89 | |
90 | AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt) |
91 | { |
92 | // assignment operator |
93 | if (this!=&mcEvnt) { |
94 | AliVEvent::operator=(mcEvnt); |
95 | } |
96 | |
97 | return *this; |
98 | } |
99 | |
100 | void AliMCEvent::ConnectTreeE (TTree* tree) |
101 | { |
102 | // Connect the event header tree |
103 | tree->SetBranchAddress("Header", &fHeader); |
104 | } |
105 | |
106 | void AliMCEvent::ConnectTreeK (TTree* tree) |
107 | { |
108 | // Connect the kinematics tree to the stack |
7aad0c47 |
109 | if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000); |
110 | // |
415d9f5c |
111 | fStack = fHeader->Stack(); |
112 | fStack->ConnectTree(tree); |
113 | // |
114 | // Load the event |
115 | fStack->GetEvent(); |
116 | fNparticles = fStack->GetNtrack(); |
117 | fNprimaries = fStack->GetNprimary(); |
477fb828 |
118 | |
ed97dc98 |
119 | Int_t iev = fHeader->GetEvent(); |
120 | Int_t ievr = fHeader->GetEventNrInRun(); |
477fb828 |
121 | AliDebug(1, Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n", |
ed97dc98 |
122 | iev, ievr, fNparticles, fNprimaries)); |
415d9f5c |
123 | |
124 | // This is a cache for the TParticles converted to MCParticles on user request |
125 | if (fMCParticleMap) { |
678a0069 |
126 | fMCParticleMap->Clear(); |
ee2774b7 |
127 | fMCParticles->Delete(); |
678a0069 |
128 | if (fNparticles>0) fMCParticleMap->Expand(fNparticles); |
129 | } |
415d9f5c |
130 | else |
93df0e9b |
131 | fMCParticleMap = new TRefArray(fNparticles); |
415d9f5c |
132 | } |
133 | |
134 | void AliMCEvent::ConnectTreeTR (TTree* tree) |
135 | { |
136 | // Connect the track reference tree |
137 | fTreeTR = tree; |
138 | |
139 | if (fTreeTR->GetBranch("AliRun")) { |
140 | if (fTmpFileTR) { |
141 | fTmpFileTR->Close(); |
142 | delete fTmpFileTR; |
143 | } |
144 | // This is an old format with one branch per detector not in synch with TreeK |
145 | ReorderAndExpandTreeTR(); |
146 | } else { |
147 | // New format |
93df0e9b |
148 | fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer); |
415d9f5c |
149 | } |
150 | } |
151 | |
152 | Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs) |
153 | { |
154 | // Retrieve entry i |
155 | if (i < 0 || i >= fNparticles) { |
156 | AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range")); |
157 | particle = 0; |
158 | trefs = 0; |
159 | return (-1); |
160 | } |
161 | particle = fStack->Particle(i); |
162 | if (fTreeTR) { |
163 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); |
93df0e9b |
164 | trefs = fTRBuffer; |
415d9f5c |
165 | return trefs->GetEntries(); |
166 | } else { |
167 | trefs = 0; |
168 | return -1; |
169 | } |
170 | } |
171 | |
172 | |
173 | void AliMCEvent::Clean() |
174 | { |
175 | // Clean-up before new trees are connected |
c6360f33 |
176 | delete fStack; fStack = 0; |
415d9f5c |
177 | |
178 | // Clear TR |
93df0e9b |
179 | if (fTRBuffer) { |
ee2774b7 |
180 | fTRBuffer->Delete(); |
93df0e9b |
181 | delete fTRBuffer; |
182 | fTRBuffer = 0; |
415d9f5c |
183 | } |
184 | } |
185 | |
04cb11d4 |
186 | #include <iostream> |
187 | |
415d9f5c |
188 | void AliMCEvent::FinishEvent() |
189 | { |
f1ac8a97 |
190 | // Clean-up after event |
c6360f33 |
191 | // |
67d77d71 |
192 | if (fStack) fStack->Reset(0); |
c6360f33 |
193 | fMCParticles->Delete(); |
04cb11d4 |
194 | |
195 | if (fMCParticleMap) |
196 | fMCParticleMap->Clear(); |
04cb11d4 |
197 | if (fTRBuffer) { |
568d5aed |
198 | fTRBuffer->Delete(); |
04cb11d4 |
199 | } |
200 | // fTrackReferences->Delete(); |
04cb11d4 |
201 | fTrackReferences->Clear(); |
2bb794ba |
202 | fNparticles = -1; |
203 | fNprimaries = -1; |
204 | fStack = 0; |
93836e1b |
205 | // fSubsidiaryEvents->Clear(); |
206 | fSubsidiaryEvents = 0; |
415d9f5c |
207 | } |
208 | |
209 | |
93df0e9b |
210 | |
415d9f5c |
211 | void AliMCEvent::DrawCheck(Int_t i, Int_t search) |
212 | { |
213 | // |
214 | // Simple event display for debugging |
215 | if (!fTreeTR) { |
216 | AliWarning("No Track Reference information available"); |
217 | return; |
218 | } |
219 | |
220 | if (i > -1 && i < fNparticles) { |
221 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); |
222 | } else { |
223 | AliWarning("AliMCEvent::GetEntry: Index out of range"); |
224 | } |
225 | |
93df0e9b |
226 | Int_t nh = fTRBuffer->GetEntries(); |
415d9f5c |
227 | |
228 | |
229 | if (search) { |
230 | while(nh <= search && i < fNparticles - 1) { |
231 | i++; |
232 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); |
93df0e9b |
233 | nh = fTRBuffer->GetEntries(); |
415d9f5c |
234 | } |
235 | printf("Found Hits at %5d\n", i); |
236 | } |
237 | TParticle* particle = fStack->Particle(i); |
238 | |
239 | TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500); |
240 | Float_t x0 = particle->Vx(); |
241 | Float_t y0 = particle->Vy(); |
242 | |
243 | Float_t x1 = particle->Vx() + particle->Px() * 50.; |
244 | Float_t y1 = particle->Vy() + particle->Py() * 50.; |
245 | |
246 | TArrow* a = new TArrow(x0, y0, x1, y1, 0.01); |
247 | h->Draw(); |
248 | a->SetLineColor(2); |
249 | |
250 | a->Draw(); |
251 | |
252 | for (Int_t ih = 0; ih < nh; ih++) { |
93df0e9b |
253 | AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih); |
415d9f5c |
254 | TMarker* m = new TMarker(ref->X(), ref->Y(), 20); |
255 | m->Draw(); |
256 | m->SetMarkerSize(0.4); |
257 | |
258 | } |
259 | } |
260 | |
261 | |
262 | void AliMCEvent::ReorderAndExpandTreeTR() |
263 | { |
264 | // |
265 | // Reorder and expand the track reference tree in order to match the kinematics tree. |
266 | // Copy the information from different branches into one |
267 | // |
268 | // TreeTR |
269 | |
270 | fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate"); |
271 | fTmpTreeTR = new TTree("TreeTR", "TrackReferences"); |
93df0e9b |
272 | if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100); |
04cb11d4 |
273 | fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 64000, 0); |
415d9f5c |
274 | |
275 | |
276 | // |
277 | // Activate the used branches only. Otherwisw we get a bad memory leak. |
278 | fTreeTR->SetBranchStatus("*", 0); |
279 | fTreeTR->SetBranchStatus("AliRun.*", 1); |
280 | fTreeTR->SetBranchStatus("ITS.*", 1); |
281 | fTreeTR->SetBranchStatus("TPC.*", 1); |
282 | fTreeTR->SetBranchStatus("TRD.*", 1); |
283 | fTreeTR->SetBranchStatus("TOF.*", 1); |
284 | fTreeTR->SetBranchStatus("FRAME.*", 1); |
285 | fTreeTR->SetBranchStatus("MUON.*", 1); |
286 | // |
287 | // Connect the active branches |
288 | TClonesArray* trefs[7]; |
289 | for (Int_t i = 0; i < 7; i++) trefs[i] = 0; |
290 | if (fTreeTR){ |
291 | // make branch for central track references |
292 | if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]); |
293 | if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]); |
294 | if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]); |
295 | if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]); |
296 | if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]); |
297 | if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]); |
298 | if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]); |
299 | } |
300 | |
301 | Int_t np = fStack->GetNprimary(); |
302 | Int_t nt = fTreeTR->GetEntries(); |
303 | |
304 | // |
305 | // Loop over tracks and find the secondaries with the help of the kine tree |
306 | Int_t ifills = 0; |
307 | Int_t it = 0; |
308 | Int_t itlast = 0; |
309 | TParticle* part; |
310 | |
311 | for (Int_t ip = np - 1; ip > -1; ip--) { |
312 | part = fStack->Particle(ip); |
313 | // printf("Particle %5d %5d %5d %5d %5d %5d \n", |
314 | // ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), |
315 | // part->GetLastDaughter(), part->TestBit(kTransportBit)); |
316 | |
317 | // Determine range of secondaries produced by this primary during transport |
318 | Int_t dau1 = part->GetFirstDaughter(); |
319 | if (dau1 < np) continue; // This particle has no secondaries produced during transport |
320 | Int_t dau2 = -1; |
321 | if (dau1 > -1) { |
322 | Int_t inext = ip - 1; |
323 | while (dau2 < 0) { |
324 | if (inext >= 0) { |
325 | part = fStack->Particle(inext); |
326 | dau2 = part->GetFirstDaughter(); |
327 | if (dau2 == -1 || dau2 < np) { |
328 | dau2 = -1; |
329 | } else { |
330 | dau2--; |
331 | } |
332 | } else { |
333 | dau2 = fStack->GetNtrack() - 1; |
334 | } |
335 | inext--; |
336 | } // find upper bound |
337 | } // dau2 < 0 |
338 | |
339 | |
340 | // printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2); |
341 | // |
342 | // Loop over reference hits and find secondary label |
343 | // First the tricky part: find the entry in treeTR than contains the hits or |
344 | // make sure that no hits exist. |
345 | // |
346 | Bool_t hasHits = kFALSE; |
347 | Bool_t isOutside = kFALSE; |
348 | |
349 | it = itlast; |
350 | while (!hasHits && !isOutside && it < nt) { |
351 | fTreeTR->GetEntry(it++); |
352 | for (Int_t ib = 0; ib < 7; ib++) { |
353 | if (!trefs[ib]) continue; |
354 | Int_t nh = trefs[ib]->GetEntries(); |
355 | for (Int_t ih = 0; ih < nh; ih++) { |
356 | AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih); |
357 | Int_t label = tr->Label(); |
358 | if (label >= dau1 && label <= dau2) { |
359 | hasHits = kTRUE; |
360 | itlast = it - 1; |
361 | break; |
362 | } |
363 | if (label > dau2 || label < ip) { |
364 | isOutside = kTRUE; |
365 | itlast = it - 1; |
366 | break; |
367 | } |
368 | } // hits |
369 | if (hasHits || isOutside) break; |
370 | } // branches |
371 | } // entries |
372 | |
373 | if (!hasHits) { |
374 | // Write empty entries |
375 | for (Int_t id = dau1; (id <= dau2); id++) { |
376 | fTmpTreeTR->Fill(); |
377 | ifills++; |
378 | } |
379 | } else { |
380 | // Collect all hits |
381 | fTreeTR->GetEntry(itlast); |
382 | for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) { |
383 | for (Int_t ib = 0; ib < 7; ib++) { |
384 | if (!trefs[ib]) continue; |
385 | Int_t nh = trefs[ib]->GetEntries(); |
386 | for (Int_t ih = 0; ih < nh; ih++) { |
387 | AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih); |
388 | Int_t label = tr->Label(); |
389 | // Skip primaries |
390 | if (label == ip) continue; |
391 | if (label > dau2 || label < dau1) |
392 | printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n", |
393 | itlast, label, dau1, dau2); |
394 | if (label == id) { |
395 | // secondary found |
396 | tr->SetDetectorId(ib-1); |
93df0e9b |
397 | Int_t nref = fTRBuffer->GetEntriesFast(); |
398 | TClonesArray &lref = *fTRBuffer; |
415d9f5c |
399 | new(lref[nref]) AliTrackReference(*tr); |
400 | } |
401 | } // hits |
402 | } // branches |
403 | fTmpTreeTR->Fill(); |
ee2774b7 |
404 | fTRBuffer->Delete(); |
415d9f5c |
405 | ifills++; |
406 | } // daughters |
407 | } // has hits |
408 | } // tracks |
409 | |
410 | // |
411 | // Now loop again and write the primaries |
412 | // |
413 | it = nt - 1; |
414 | for (Int_t ip = 0; ip < np; ip++) { |
415 | Int_t labmax = -1; |
416 | while (labmax < ip && it > -1) { |
417 | fTreeTR->GetEntry(it--); |
418 | for (Int_t ib = 0; ib < 7; ib++) { |
419 | if (!trefs[ib]) continue; |
420 | Int_t nh = trefs[ib]->GetEntries(); |
421 | // |
422 | // Loop over reference hits and find primary labels |
423 | for (Int_t ih = 0; ih < nh; ih++) { |
424 | AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih); |
425 | Int_t label = tr->Label(); |
426 | if (label < np && label > labmax) { |
427 | labmax = label; |
428 | } |
429 | |
430 | if (label == ip) { |
431 | tr->SetDetectorId(ib-1); |
93df0e9b |
432 | Int_t nref = fTRBuffer->GetEntriesFast(); |
433 | TClonesArray &lref = *fTRBuffer; |
415d9f5c |
434 | new(lref[nref]) AliTrackReference(*tr); |
435 | } |
436 | } // hits |
437 | } // branches |
438 | } // entries |
439 | it++; |
440 | fTmpTreeTR->Fill(); |
ee2774b7 |
441 | fTRBuffer->Delete(); |
415d9f5c |
442 | ifills++; |
443 | } // tracks |
444 | // Check |
445 | |
446 | |
447 | // Clean-up |
448 | delete fTreeTR; fTreeTR = 0; |
449 | |
450 | for (Int_t ib = 0; ib < 7; ib++) { |
451 | if (trefs[ib]) { |
452 | trefs[ib]->Clear(); |
453 | delete trefs[ib]; |
454 | trefs[ib] = 0; |
455 | } |
456 | } |
457 | |
458 | if (ifills != fStack->GetNtrack()) |
459 | printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", |
460 | ifills, fStack->GetNtrack()); |
461 | |
462 | fTmpTreeTR->Write(); |
463 | fTreeTR = fTmpTreeTR; |
464 | } |
465 | |
7aad0c47 |
466 | AliVParticle* AliMCEvent::GetTrack(Int_t i) const |
415d9f5c |
467 | { |
468 | // Get MC Particle i |
93836e1b |
469 | // |
7aad0c47 |
470 | |
471 | if (fExternal) { |
472 | return ((AliVParticle*) (fMCParticles->At(i))); |
473 | } |
474 | |
93836e1b |
475 | // |
476 | // Check first if this explicitely accesses the subsidiary event |
7aad0c47 |
477 | |
93836e1b |
478 | if (i >= BgLabelOffset()) { |
479 | if (fSubsidiaryEvents) { |
480 | AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1)); |
481 | return (bgEvent->GetTrack(i - BgLabelOffset())); |
482 | } else { |
483 | return 0; |
484 | } |
485 | } |
486 | |
487 | // |
93df0e9b |
488 | AliMCParticle *mcParticle = 0; |
489 | TParticle *particle = 0; |
490 | TClonesArray *trefs = 0; |
491 | Int_t ntref = 0; |
492 | TRefArray *rarray = 0; |
93836e1b |
493 | |
494 | |
495 | |
415d9f5c |
496 | // Out of range check |
497 | if (i < 0 || i >= fNparticles) { |
498 | AliWarning(Form("AliMCEvent::GetEntry: Index out of range")); |
499 | mcParticle = 0; |
500 | return (mcParticle); |
501 | } |
93836e1b |
502 | |
415d9f5c |
503 | |
93836e1b |
504 | if (fSubsidiaryEvents) { |
505 | AliMCEvent* mc; |
506 | Int_t idx = FindIndexAndEvent(i, mc); |
507 | return (mc->GetTrack(idx)); |
508 | } |
415d9f5c |
509 | |
415d9f5c |
510 | // |
93836e1b |
511 | // First check If the MC Particle has been already cached |
415d9f5c |
512 | if(!fMCParticleMap->At(i)) { |
93df0e9b |
513 | // Get particle from the stack |
514 | particle = fStack->Particle(i); |
515 | // Get track references from Tree TR |
516 | if (fTreeTR) { |
517 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); |
518 | trefs = fTRBuffer; |
519 | ntref = trefs->GetEntriesFast(); |
520 | rarray = new TRefArray(ntref); |
521 | Int_t nen = fTrackReferences->GetEntriesFast(); |
522 | for (Int_t j = 0; j < ntref; j++) { |
523 | // Save the track references in a TClonesArray |
524 | AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]); |
525 | // Save the pointer in a TRefArray |
526 | new ((*fTrackReferences)[nen]) AliTrackReference(*ref); |
527 | rarray->AddAt((*fTrackReferences)[nen], j); |
528 | nen++; |
529 | } // loop over track references for entry i |
530 | } // if TreeTR available |
415d9f5c |
531 | Int_t nentries = fMCParticles->GetEntriesFast(); |
eee13e8d |
532 | new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i); |
415d9f5c |
533 | mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]); |
534 | fMCParticleMap->AddAt(mcParticle, i); |
93836e1b |
535 | |
536 | TParticle* part = mcParticle->Particle(); |
537 | Int_t imo = part->GetFirstMother(); |
538 | Int_t id1 = part->GetFirstDaughter(); |
539 | Int_t id2 = part->GetLastDaughter(); |
540 | if (fPrimaryOffset > 0 || fSecondaryOffset > 0) { |
541 | // Remapping of the mother and daughter indices |
542 | if (imo < fNprimaries) { |
543 | mcParticle->SetMother(imo + fPrimaryOffset); |
544 | } else { |
545 | mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries); |
546 | } |
547 | |
548 | if (id1 < fNprimaries) { |
549 | mcParticle->SetFirstDaughter(id1 + fPrimaryOffset); |
550 | mcParticle->SetLastDaughter (id2 + fPrimaryOffset); |
551 | } else { |
552 | mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries); |
553 | mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries); |
554 | } |
555 | |
556 | |
557 | if (i > fNprimaries) { |
558 | mcParticle->SetLabel(i + fPrimaryOffset); |
559 | } else { |
560 | mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries); |
561 | } |
562 | |
563 | |
564 | } else { |
565 | mcParticle->SetFirstDaughter(id1); |
566 | mcParticle->SetLastDaughter (id2); |
567 | mcParticle->SetMother (imo); |
568 | } |
569 | |
415d9f5c |
570 | } else { |
571 | mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i)); |
572 | } |
93836e1b |
573 | |
574 | |
415d9f5c |
575 | return mcParticle; |
576 | } |
577 | |
5df3496b |
578 | AliGenEventHeader* AliMCEvent::GenEventHeader() const {return (fHeader->GenEventHeader());} |
93836e1b |
579 | |
580 | |
581 | void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event) |
582 | { |
583 | // Add a subsidiary event to the list; for example merged background event. |
584 | if (!fSubsidiaryEvents) { |
585 | TList* events = new TList(); |
586 | events->Add(new AliMCEvent(*this)); |
587 | fSubsidiaryEvents = events; |
588 | } |
589 | |
590 | fSubsidiaryEvents->Add(event); |
591 | } |
592 | |
593 | Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const |
594 | { |
595 | // Find the index and event in case of composed events like signal + background |
596 | TIter next(fSubsidiaryEvents); |
597 | next.Reset(); |
598 | if (oldidx < fNprimaries) { |
599 | while((event = (AliMCEvent*)next())) { |
600 | if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break; |
601 | } |
602 | return (oldidx - event->GetPrimaryOffset()); |
603 | } else { |
604 | while((event = (AliMCEvent*)next())) { |
605 | if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break; |
606 | } |
607 | return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries()); |
608 | } |
609 | } |
610 | |
611 | Int_t AliMCEvent::BgLabelToIndex(Int_t label) |
612 | { |
613 | // Convert a background label to an absolute index |
614 | if (fSubsidiaryEvents) { |
615 | AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1)); |
616 | label -= BgLabelOffset(); |
617 | if (label < bgEvent->GetNumberOfPrimaries()) { |
618 | label += bgEvent->GetPrimaryOffset(); |
619 | } else { |
620 | label += (bgEvent->GetSecondaryOffset() - fNprimaries); |
621 | } |
622 | } |
623 | return (label); |
624 | } |
625 | |
626 | |
627 | Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) |
628 | { |
629 | // |
630 | // Delegate to subevent if necesarry |
631 | |
632 | |
633 | if (!fSubsidiaryEvents) { |
634 | return fStack->IsPhysicalPrimary(i); |
635 | } else { |
636 | AliMCEvent* evt = 0; |
637 | Int_t idx = FindIndexAndEvent(i, evt); |
638 | return (evt->IsPhysicalPrimary(idx)); |
639 | } |
640 | } |
641 | |
642 | void AliMCEvent::InitEvent() |
643 | { |
35152e4b |
644 | // |
645 | // Initialize the subsidiary event structure |
93836e1b |
646 | if (fSubsidiaryEvents) { |
647 | TIter next(fSubsidiaryEvents); |
648 | AliMCEvent* evt; |
649 | fNprimaries = 0; |
650 | fNparticles = 0; |
651 | |
652 | while((evt = (AliMCEvent*)next())) { |
653 | fNprimaries += evt->GetNumberOfPrimaries(); |
654 | fNparticles += evt->GetNumberOfTracks(); |
655 | } |
656 | |
657 | Int_t ioffp = 0; |
658 | Int_t ioffs = fNprimaries; |
659 | next.Reset(); |
660 | |
661 | while((evt = (AliMCEvent*)next())) { |
662 | evt->SetPrimaryOffset(ioffp); |
663 | evt->SetSecondaryOffset(ioffs); |
664 | ioffp += evt->GetNumberOfPrimaries(); |
665 | ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries()); |
666 | } |
667 | } |
668 | } |
35152e4b |
669 | void AliMCEvent::PreReadAll() |
670 | { |
671 | // Preread the MC information |
672 | Int_t i; |
673 | // secondaries |
674 | for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++) |
675 | { |
676 | GetTrack(i); |
677 | } |
678 | // primaries |
679 | for (i = 0; i < fStack->GetNprimary(); i++) |
680 | { |
681 | GetTrack(i); |
682 | } |
683 | |
684 | |
685 | } |
93836e1b |
686 | |
415d9f5c |
687 | |
5df3496b |
688 | const AliVVertex * AliMCEvent::GetPrimaryVertex() const |
689 | { |
690 | // Create a MCVertex object from the MCHeader information |
570308e4 |
691 | TArrayF v; |
692 | GenEventHeader()->PrimaryVertex(v) ; |
5df3496b |
693 | if (!fVertex) { |
5df3496b |
694 | fVertex = new AliMCVertex(v[0], v[1], v[2]); |
570308e4 |
695 | } else { |
696 | ((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]); |
5df3496b |
697 | } |
5df3496b |
698 | return fVertex; |
699 | } |
700 | |
701 | |
415d9f5c |
702 | ClassImp(AliMCEvent) |