]>
Commit | Line | Data |
---|---|---|
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> |
415d9f5c | 30 | |
31 | #include "AliLog.h" | |
32 | #include "AliMCEvent.h" | |
33 | #include "AliStack.h" | |
34 | #include "AliTrackReference.h" | |
35 | #include "AliHeader.h" | |
83626b6b | 36 | #include "AliGenEventHeader.h" |
415d9f5c | 37 | |
38 | ||
39 | AliMCEvent::AliMCEvent(): | |
40 | AliVEvent(), | |
41 | fStack(0), | |
42 | fMCParticles(new TClonesArray("AliMCParticle",1000)), | |
43 | fMCParticleMap(0), | |
fb9aeda5 | 44 | fHeader(new AliHeader()), |
93df0e9b | 45 | fTRBuffer(0), |
46 | fTrackReferences(new TClonesArray("AliTrackReference", 1000)), | |
415d9f5c | 47 | fTreeTR(0), |
48 | fTmpTreeTR(0), | |
49 | fTmpFileTR(0), | |
50 | fNprimaries(-1), | |
51 | fNparticles(-1) | |
52 | { | |
53 | // Default constructor | |
54 | } | |
55 | ||
56 | AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) : | |
57 | AliVEvent(mcEvnt), | |
58 | fStack(0), | |
59 | fMCParticles(0), | |
60 | fMCParticleMap(0), | |
fb9aeda5 | 61 | fHeader(new AliHeader()), |
93df0e9b | 62 | fTRBuffer(0), |
415d9f5c | 63 | fTrackReferences(0), |
64 | fTreeTR(0), | |
65 | fTmpTreeTR(0), | |
66 | fTmpFileTR(0), | |
67 | fNprimaries(-1), | |
678a0069 | 68 | fNparticles(-1) |
415d9f5c | 69 | { |
70 | // Copy constructor | |
71 | } | |
72 | ||
73 | ||
74 | AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt) | |
75 | { | |
76 | // assignment operator | |
77 | if (this!=&mcEvnt) { | |
78 | AliVEvent::operator=(mcEvnt); | |
79 | } | |
80 | ||
81 | return *this; | |
82 | } | |
83 | ||
84 | void AliMCEvent::ConnectTreeE (TTree* tree) | |
85 | { | |
86 | // Connect the event header tree | |
87 | tree->SetBranchAddress("Header", &fHeader); | |
88 | } | |
89 | ||
90 | void AliMCEvent::ConnectTreeK (TTree* tree) | |
91 | { | |
92 | // Connect the kinematics tree to the stack | |
93 | fStack = fHeader->Stack(); | |
94 | fStack->ConnectTree(tree); | |
95 | // | |
96 | // Load the event | |
97 | fStack->GetEvent(); | |
98 | fNparticles = fStack->GetNtrack(); | |
99 | fNprimaries = fStack->GetNprimary(); | |
ed97dc98 | 100 | Int_t iev = fHeader->GetEvent(); |
101 | Int_t ievr = fHeader->GetEventNrInRun(); | |
102 | ||
103 | AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n", | |
104 | iev, ievr, fNparticles, fNprimaries)); | |
415d9f5c | 105 | |
106 | // This is a cache for the TParticles converted to MCParticles on user request | |
107 | if (fMCParticleMap) { | |
678a0069 | 108 | fMCParticleMap->Clear(); |
ee2774b7 | 109 | fMCParticles->Delete(); |
678a0069 | 110 | if (fNparticles>0) fMCParticleMap->Expand(fNparticles); |
111 | } | |
415d9f5c | 112 | else |
93df0e9b | 113 | fMCParticleMap = new TRefArray(fNparticles); |
415d9f5c | 114 | } |
115 | ||
116 | void AliMCEvent::ConnectTreeTR (TTree* tree) | |
117 | { | |
118 | // Connect the track reference tree | |
119 | fTreeTR = tree; | |
120 | ||
121 | if (fTreeTR->GetBranch("AliRun")) { | |
122 | if (fTmpFileTR) { | |
123 | fTmpFileTR->Close(); | |
124 | delete fTmpFileTR; | |
125 | } | |
126 | // This is an old format with one branch per detector not in synch with TreeK | |
127 | ReorderAndExpandTreeTR(); | |
128 | } else { | |
129 | // New format | |
93df0e9b | 130 | fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer); |
415d9f5c | 131 | } |
132 | } | |
133 | ||
134 | Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs) | |
135 | { | |
136 | // Retrieve entry i | |
137 | if (i < 0 || i >= fNparticles) { | |
138 | AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range")); | |
139 | particle = 0; | |
140 | trefs = 0; | |
141 | return (-1); | |
142 | } | |
143 | particle = fStack->Particle(i); | |
144 | if (fTreeTR) { | |
145 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); | |
93df0e9b | 146 | trefs = fTRBuffer; |
415d9f5c | 147 | return trefs->GetEntries(); |
148 | } else { | |
149 | trefs = 0; | |
150 | return -1; | |
151 | } | |
152 | } | |
153 | ||
154 | ||
155 | void AliMCEvent::Clean() | |
156 | { | |
157 | // Clean-up before new trees are connected | |
fb9aeda5 | 158 | // if (fHeader) { |
159 | // delete fHeader; | |
160 | // fHeader = 0; | |
161 | // } | |
c6360f33 | 162 | |
163 | delete fStack; fStack = 0; | |
415d9f5c | 164 | |
165 | // Clear TR | |
93df0e9b | 166 | if (fTRBuffer) { |
ee2774b7 | 167 | fTRBuffer->Delete(); |
93df0e9b | 168 | delete fTRBuffer; |
169 | fTRBuffer = 0; | |
415d9f5c | 170 | } |
171 | } | |
172 | ||
173 | void AliMCEvent::FinishEvent() | |
174 | { | |
f1ac8a97 | 175 | // Clean-up after event |
c6360f33 | 176 | // |
177 | fStack->Reset(0); | |
178 | fMCParticles->Delete(); | |
678a0069 | 179 | fMCParticleMap->Clear(); |
568d5aed | 180 | if (fTRBuffer) |
181 | fTRBuffer->Delete(); | |
c6360f33 | 182 | fTrackReferences->Delete(); |
2bb794ba | 183 | fNparticles = -1; |
184 | fNprimaries = -1; | |
185 | fStack = 0; | |
186 | ||
415d9f5c | 187 | } |
188 | ||
189 | ||
93df0e9b | 190 | |
415d9f5c | 191 | void AliMCEvent::DrawCheck(Int_t i, Int_t search) |
192 | { | |
193 | // | |
194 | // Simple event display for debugging | |
195 | if (!fTreeTR) { | |
196 | AliWarning("No Track Reference information available"); | |
197 | return; | |
198 | } | |
199 | ||
200 | if (i > -1 && i < fNparticles) { | |
201 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); | |
202 | } else { | |
203 | AliWarning("AliMCEvent::GetEntry: Index out of range"); | |
204 | } | |
205 | ||
93df0e9b | 206 | Int_t nh = fTRBuffer->GetEntries(); |
415d9f5c | 207 | |
208 | ||
209 | if (search) { | |
210 | while(nh <= search && i < fNparticles - 1) { | |
211 | i++; | |
212 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); | |
93df0e9b | 213 | nh = fTRBuffer->GetEntries(); |
415d9f5c | 214 | } |
215 | printf("Found Hits at %5d\n", i); | |
216 | } | |
217 | TParticle* particle = fStack->Particle(i); | |
218 | ||
219 | TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500); | |
220 | Float_t x0 = particle->Vx(); | |
221 | Float_t y0 = particle->Vy(); | |
222 | ||
223 | Float_t x1 = particle->Vx() + particle->Px() * 50.; | |
224 | Float_t y1 = particle->Vy() + particle->Py() * 50.; | |
225 | ||
226 | TArrow* a = new TArrow(x0, y0, x1, y1, 0.01); | |
227 | h->Draw(); | |
228 | a->SetLineColor(2); | |
229 | ||
230 | a->Draw(); | |
231 | ||
232 | for (Int_t ih = 0; ih < nh; ih++) { | |
93df0e9b | 233 | AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih); |
415d9f5c | 234 | TMarker* m = new TMarker(ref->X(), ref->Y(), 20); |
235 | m->Draw(); | |
236 | m->SetMarkerSize(0.4); | |
237 | ||
238 | } | |
239 | } | |
240 | ||
241 | ||
242 | void AliMCEvent::ReorderAndExpandTreeTR() | |
243 | { | |
244 | // | |
245 | // Reorder and expand the track reference tree in order to match the kinematics tree. | |
246 | // Copy the information from different branches into one | |
247 | // | |
248 | // TreeTR | |
249 | ||
250 | fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate"); | |
251 | fTmpTreeTR = new TTree("TreeTR", "TrackReferences"); | |
93df0e9b | 252 | if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100); |
253 | fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0); | |
415d9f5c | 254 | |
255 | ||
256 | // | |
257 | // Activate the used branches only. Otherwisw we get a bad memory leak. | |
258 | fTreeTR->SetBranchStatus("*", 0); | |
259 | fTreeTR->SetBranchStatus("AliRun.*", 1); | |
260 | fTreeTR->SetBranchStatus("ITS.*", 1); | |
261 | fTreeTR->SetBranchStatus("TPC.*", 1); | |
262 | fTreeTR->SetBranchStatus("TRD.*", 1); | |
263 | fTreeTR->SetBranchStatus("TOF.*", 1); | |
264 | fTreeTR->SetBranchStatus("FRAME.*", 1); | |
265 | fTreeTR->SetBranchStatus("MUON.*", 1); | |
266 | // | |
267 | // Connect the active branches | |
268 | TClonesArray* trefs[7]; | |
269 | for (Int_t i = 0; i < 7; i++) trefs[i] = 0; | |
270 | if (fTreeTR){ | |
271 | // make branch for central track references | |
272 | if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]); | |
273 | if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]); | |
274 | if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]); | |
275 | if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]); | |
276 | if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]); | |
277 | if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]); | |
278 | if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]); | |
279 | } | |
280 | ||
281 | Int_t np = fStack->GetNprimary(); | |
282 | Int_t nt = fTreeTR->GetEntries(); | |
283 | ||
284 | // | |
285 | // Loop over tracks and find the secondaries with the help of the kine tree | |
286 | Int_t ifills = 0; | |
287 | Int_t it = 0; | |
288 | Int_t itlast = 0; | |
289 | TParticle* part; | |
290 | ||
291 | for (Int_t ip = np - 1; ip > -1; ip--) { | |
292 | part = fStack->Particle(ip); | |
293 | // printf("Particle %5d %5d %5d %5d %5d %5d \n", | |
294 | // ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), | |
295 | // part->GetLastDaughter(), part->TestBit(kTransportBit)); | |
296 | ||
297 | // Determine range of secondaries produced by this primary during transport | |
298 | Int_t dau1 = part->GetFirstDaughter(); | |
299 | if (dau1 < np) continue; // This particle has no secondaries produced during transport | |
300 | Int_t dau2 = -1; | |
301 | if (dau1 > -1) { | |
302 | Int_t inext = ip - 1; | |
303 | while (dau2 < 0) { | |
304 | if (inext >= 0) { | |
305 | part = fStack->Particle(inext); | |
306 | dau2 = part->GetFirstDaughter(); | |
307 | if (dau2 == -1 || dau2 < np) { | |
308 | dau2 = -1; | |
309 | } else { | |
310 | dau2--; | |
311 | } | |
312 | } else { | |
313 | dau2 = fStack->GetNtrack() - 1; | |
314 | } | |
315 | inext--; | |
316 | } // find upper bound | |
317 | } // dau2 < 0 | |
318 | ||
319 | ||
320 | // printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2); | |
321 | // | |
322 | // Loop over reference hits and find secondary label | |
323 | // First the tricky part: find the entry in treeTR than contains the hits or | |
324 | // make sure that no hits exist. | |
325 | // | |
326 | Bool_t hasHits = kFALSE; | |
327 | Bool_t isOutside = kFALSE; | |
328 | ||
329 | it = itlast; | |
330 | while (!hasHits && !isOutside && it < nt) { | |
331 | fTreeTR->GetEntry(it++); | |
332 | for (Int_t ib = 0; ib < 7; ib++) { | |
333 | if (!trefs[ib]) continue; | |
334 | Int_t nh = trefs[ib]->GetEntries(); | |
335 | for (Int_t ih = 0; ih < nh; ih++) { | |
336 | AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih); | |
337 | Int_t label = tr->Label(); | |
338 | if (label >= dau1 && label <= dau2) { | |
339 | hasHits = kTRUE; | |
340 | itlast = it - 1; | |
341 | break; | |
342 | } | |
343 | if (label > dau2 || label < ip) { | |
344 | isOutside = kTRUE; | |
345 | itlast = it - 1; | |
346 | break; | |
347 | } | |
348 | } // hits | |
349 | if (hasHits || isOutside) break; | |
350 | } // branches | |
351 | } // entries | |
352 | ||
353 | if (!hasHits) { | |
354 | // Write empty entries | |
355 | for (Int_t id = dau1; (id <= dau2); id++) { | |
356 | fTmpTreeTR->Fill(); | |
357 | ifills++; | |
358 | } | |
359 | } else { | |
360 | // Collect all hits | |
361 | fTreeTR->GetEntry(itlast); | |
362 | for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) { | |
363 | for (Int_t ib = 0; ib < 7; ib++) { | |
364 | if (!trefs[ib]) continue; | |
365 | Int_t nh = trefs[ib]->GetEntries(); | |
366 | for (Int_t ih = 0; ih < nh; ih++) { | |
367 | AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih); | |
368 | Int_t label = tr->Label(); | |
369 | // Skip primaries | |
370 | if (label == ip) continue; | |
371 | if (label > dau2 || label < dau1) | |
372 | printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n", | |
373 | itlast, label, dau1, dau2); | |
374 | if (label == id) { | |
375 | // secondary found | |
376 | tr->SetDetectorId(ib-1); | |
93df0e9b | 377 | Int_t nref = fTRBuffer->GetEntriesFast(); |
378 | TClonesArray &lref = *fTRBuffer; | |
415d9f5c | 379 | new(lref[nref]) AliTrackReference(*tr); |
380 | } | |
381 | } // hits | |
382 | } // branches | |
383 | fTmpTreeTR->Fill(); | |
ee2774b7 | 384 | fTRBuffer->Delete(); |
415d9f5c | 385 | ifills++; |
386 | } // daughters | |
387 | } // has hits | |
388 | } // tracks | |
389 | ||
390 | // | |
391 | // Now loop again and write the primaries | |
392 | // | |
393 | it = nt - 1; | |
394 | for (Int_t ip = 0; ip < np; ip++) { | |
395 | Int_t labmax = -1; | |
396 | while (labmax < ip && it > -1) { | |
397 | fTreeTR->GetEntry(it--); | |
398 | for (Int_t ib = 0; ib < 7; ib++) { | |
399 | if (!trefs[ib]) continue; | |
400 | Int_t nh = trefs[ib]->GetEntries(); | |
401 | // | |
402 | // Loop over reference hits and find primary labels | |
403 | for (Int_t ih = 0; ih < nh; ih++) { | |
404 | AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih); | |
405 | Int_t label = tr->Label(); | |
406 | if (label < np && label > labmax) { | |
407 | labmax = label; | |
408 | } | |
409 | ||
410 | if (label == ip) { | |
411 | tr->SetDetectorId(ib-1); | |
93df0e9b | 412 | Int_t nref = fTRBuffer->GetEntriesFast(); |
413 | TClonesArray &lref = *fTRBuffer; | |
415d9f5c | 414 | new(lref[nref]) AliTrackReference(*tr); |
415 | } | |
416 | } // hits | |
417 | } // branches | |
418 | } // entries | |
419 | it++; | |
420 | fTmpTreeTR->Fill(); | |
ee2774b7 | 421 | fTRBuffer->Delete(); |
415d9f5c | 422 | ifills++; |
423 | } // tracks | |
424 | // Check | |
425 | ||
426 | ||
427 | // Clean-up | |
428 | delete fTreeTR; fTreeTR = 0; | |
429 | ||
430 | for (Int_t ib = 0; ib < 7; ib++) { | |
431 | if (trefs[ib]) { | |
432 | trefs[ib]->Clear(); | |
433 | delete trefs[ib]; | |
434 | trefs[ib] = 0; | |
435 | } | |
436 | } | |
437 | ||
438 | if (ifills != fStack->GetNtrack()) | |
439 | printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", | |
440 | ifills, fStack->GetNtrack()); | |
441 | ||
442 | fTmpTreeTR->Write(); | |
443 | fTreeTR = fTmpTreeTR; | |
444 | } | |
445 | ||
446 | AliMCParticle* AliMCEvent::GetTrack(Int_t i) const | |
447 | { | |
448 | // Get MC Particle i | |
93df0e9b | 449 | AliMCParticle *mcParticle = 0; |
450 | TParticle *particle = 0; | |
451 | TClonesArray *trefs = 0; | |
452 | Int_t ntref = 0; | |
453 | TRefArray *rarray = 0; | |
415d9f5c | 454 | // Out of range check |
455 | if (i < 0 || i >= fNparticles) { | |
456 | AliWarning(Form("AliMCEvent::GetEntry: Index out of range")); | |
457 | mcParticle = 0; | |
458 | return (mcParticle); | |
459 | } | |
460 | ||
461 | ||
415d9f5c | 462 | // |
463 | // First check of the MC Particle has been already cached | |
464 | if(!fMCParticleMap->At(i)) { | |
93df0e9b | 465 | // Get particle from the stack |
466 | particle = fStack->Particle(i); | |
467 | // Get track references from Tree TR | |
468 | if (fTreeTR) { | |
469 | fTreeTR->GetEntry(fStack->TreeKEntry(i)); | |
470 | trefs = fTRBuffer; | |
471 | ntref = trefs->GetEntriesFast(); | |
472 | rarray = new TRefArray(ntref); | |
473 | Int_t nen = fTrackReferences->GetEntriesFast(); | |
474 | for (Int_t j = 0; j < ntref; j++) { | |
475 | // Save the track references in a TClonesArray | |
476 | AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]); | |
477 | // Save the pointer in a TRefArray | |
478 | new ((*fTrackReferences)[nen]) AliTrackReference(*ref); | |
479 | rarray->AddAt((*fTrackReferences)[nen], j); | |
480 | nen++; | |
481 | } // loop over track references for entry i | |
482 | } // if TreeTR available | |
415d9f5c | 483 | Int_t nentries = fMCParticles->GetEntriesFast(); |
eee13e8d | 484 | new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i); |
415d9f5c | 485 | mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]); |
486 | fMCParticleMap->AddAt(mcParticle, i); | |
487 | } else { | |
488 | mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i)); | |
489 | } | |
415d9f5c | 490 | return mcParticle; |
491 | } | |
492 | ||
83626b6b | 493 | AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());} |
415d9f5c | 494 | |
495 | ClassImp(AliMCEvent) |