Line data Source code
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>
29 : #include <TList.h>
30 : #include <TArrayF.h>
31 :
32 : #include "AliLog.h"
33 : #include "AliMCEvent.h"
34 : #include "AliMCVertex.h"
35 : #include "AliStack.h"
36 : #include "AliTrackReference.h"
37 : #include "AliHeader.h"
38 : #include "AliGenEventHeader.h"
39 : #include "AliGenHijingEventHeader.h"
40 : #include "AliGenCocktailEventHeader.h"
41 :
42 :
43 : Int_t AliMCEvent::fgkBgLabelOffset(10000000);
44 :
45 :
46 : AliMCEvent::AliMCEvent():
47 1 : AliVEvent(),
48 1 : fStack(0),
49 1 : fMCParticles(0),
50 1 : fMCParticleMap(0),
51 3 : fHeader(new AliHeader()),
52 1 : fAODMCHeader(0),
53 1 : fTRBuffer(0),
54 3 : fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
55 1 : fTreeTR(0),
56 1 : fTmpTreeTR(0),
57 1 : fTmpFileTR(0),
58 1 : fNprimaries(-1),
59 1 : fNparticles(-1),
60 1 : fSubsidiaryEvents(0),
61 1 : fPrimaryOffset(0),
62 1 : fSecondaryOffset(0),
63 1 : fExternal(0),
64 1 : fVertex(0),
65 1 : fNBG(-1)
66 5 : {
67 : // Default constructor
68 2 : }
69 :
70 : AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
71 0 : AliVEvent(mcEvnt),
72 0 : fStack(mcEvnt.fStack),
73 0 : fMCParticles(mcEvnt.fMCParticles),
74 0 : fMCParticleMap(mcEvnt.fMCParticleMap),
75 0 : fHeader(mcEvnt.fHeader),
76 0 : fAODMCHeader(mcEvnt.fAODMCHeader),
77 0 : fTRBuffer(mcEvnt.fTRBuffer),
78 0 : fTrackReferences(mcEvnt.fTrackReferences),
79 0 : fTreeTR(mcEvnt.fTreeTR),
80 0 : fTmpTreeTR(mcEvnt.fTmpTreeTR),
81 0 : fTmpFileTR(mcEvnt.fTmpFileTR),
82 0 : fNprimaries(mcEvnt.fNprimaries),
83 0 : fNparticles(mcEvnt.fNparticles),
84 0 : fSubsidiaryEvents(0),
85 0 : fPrimaryOffset(0),
86 0 : fSecondaryOffset(0),
87 0 : fExternal(0),
88 0 : fVertex(mcEvnt.fVertex),
89 0 : fNBG(mcEvnt.fNBG)
90 0 : {
91 : // Copy constructor
92 0 : }
93 :
94 :
95 : AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
96 : {
97 : // assignment operator
98 0 : if (this!=&mcEvnt) {
99 0 : AliVEvent::operator=(mcEvnt);
100 0 : }
101 :
102 0 : return *this;
103 : }
104 :
105 : void AliMCEvent::ConnectTreeE (TTree* tree)
106 : {
107 : // Connect the event header tree
108 2 : tree->SetBranchAddress("Header", &fHeader);
109 1 : }
110 :
111 : void AliMCEvent::ConnectTreeK (TTree* tree)
112 : {
113 : // Connect Kinematics tree
114 8 : fStack = fHeader->Stack();
115 4 : fStack->ConnectTree(tree);
116 : //
117 : // Load the event
118 4 : fStack->GetEvent();
119 :
120 4 : UpdateEventInformation();
121 4 : }
122 :
123 : void AliMCEvent::ConnectHeaderAndStack(AliHeader* header)
124 : {
125 : // fill MC event information from stack and header
126 :
127 0 : fHeader = header;
128 0 : fStack = fHeader->Stack();
129 :
130 0 : UpdateEventInformation();
131 0 : }
132 :
133 : void AliMCEvent::UpdateEventInformation()
134 : {
135 : // bookkeeping for next event
136 :
137 : // Connect the kinematics tree to the stack
138 10 : if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
139 :
140 : // Initialize members
141 4 : fNparticles = fStack->GetNtrack();
142 4 : fNprimaries = fStack->GetNprimary();
143 :
144 4 : Int_t iev = fHeader->GetEvent();
145 4 : Int_t ievr = fHeader->GetEventNrInRun();
146 12 : AliDebug(1, Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
147 : iev, ievr, fNparticles, fNprimaries));
148 :
149 : // This is a cache for the TParticles converted to MCParticles on user request
150 4 : if (fMCParticleMap) {
151 3 : fMCParticleMap->Clear();
152 3 : fMCParticles->Delete();
153 6 : if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
154 : }
155 : else
156 2 : fMCParticleMap = new TObjArray(fNparticles);
157 4 : }
158 :
159 : void AliMCEvent::ConnectTreeTR (TTree* tree)
160 : {
161 : // Connect the track reference tree
162 8 : fTreeTR = tree;
163 :
164 4 : if (fTreeTR->GetBranch("AliRun")) {
165 0 : if (fTmpFileTR) {
166 0 : fTmpFileTR->Close();
167 0 : delete fTmpFileTR;
168 : }
169 : // This is an old format with one branch per detector not in synch with TreeK
170 0 : ReorderAndExpandTreeTR();
171 0 : } else {
172 : // New format
173 4 : fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
174 : }
175 4 : }
176 :
177 : Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
178 : {
179 : // Retrieve entry i
180 0 : if (i < 0 || i >= fNparticles) {
181 0 : AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
182 0 : particle = 0;
183 0 : trefs = 0;
184 0 : return (-1);
185 : }
186 0 : particle = fStack->Particle(i);
187 0 : if (fTreeTR) {
188 0 : fTreeTR->GetEntry(fStack->TreeKEntry(i));
189 0 : trefs = fTRBuffer;
190 0 : return trefs->GetEntries();
191 : } else {
192 0 : trefs = 0;
193 0 : return -1;
194 : }
195 0 : }
196 :
197 :
198 : void AliMCEvent::Clean()
199 : {
200 : // Clean-up before new trees are connected
201 0 : delete fStack; fStack = 0;
202 :
203 : // Clear TR
204 0 : if (fTRBuffer) {
205 0 : fTRBuffer->Delete();
206 0 : delete fTRBuffer;
207 0 : fTRBuffer = 0;
208 0 : }
209 0 : }
210 :
211 : #include <iostream>
212 :
213 : void AliMCEvent::FinishEvent()
214 : {
215 : // Clean-up after event
216 : //
217 12 : if (fStack) fStack->Reset(0);
218 4 : fMCParticles->Delete();
219 :
220 4 : if (fMCParticleMap)
221 4 : fMCParticleMap->Clear();
222 4 : if (fTRBuffer) {
223 4 : fTRBuffer->Delete();
224 4 : }
225 : // fTrackReferences->Delete();
226 4 : fTrackReferences->Clear();
227 4 : fNparticles = -1;
228 4 : fNprimaries = -1;
229 4 : fStack = 0;
230 : // fSubsidiaryEvents->Clear();
231 4 : fSubsidiaryEvents = 0;
232 4 : fNBG = -1;
233 4 : }
234 :
235 :
236 :
237 : void AliMCEvent::DrawCheck(Int_t i, Int_t search)
238 : {
239 : //
240 : // Simple event display for debugging
241 0 : if (!fTreeTR) {
242 0 : AliWarning("No Track Reference information available");
243 0 : return;
244 : }
245 :
246 0 : if (i > -1 && i < fNparticles) {
247 0 : fTreeTR->GetEntry(fStack->TreeKEntry(i));
248 0 : } else {
249 0 : AliWarning("AliMCEvent::GetEntry: Index out of range");
250 : }
251 :
252 0 : Int_t nh = fTRBuffer->GetEntries();
253 :
254 :
255 0 : if (search) {
256 0 : while(nh <= search && i < fNparticles - 1) {
257 0 : i++;
258 0 : fTreeTR->GetEntry(fStack->TreeKEntry(i));
259 0 : nh = fTRBuffer->GetEntries();
260 : }
261 0 : printf("Found Hits at %5d\n", i);
262 0 : }
263 0 : TParticle* particle = fStack->Particle(i);
264 :
265 0 : TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
266 0 : Float_t x0 = particle->Vx();
267 0 : Float_t y0 = particle->Vy();
268 :
269 0 : Float_t x1 = particle->Vx() + particle->Px() * 50.;
270 0 : Float_t y1 = particle->Vy() + particle->Py() * 50.;
271 :
272 0 : TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
273 0 : h->Draw();
274 0 : a->SetLineColor(2);
275 :
276 0 : a->Draw();
277 :
278 0 : for (Int_t ih = 0; ih < nh; ih++) {
279 0 : AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
280 0 : TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
281 0 : m->Draw();
282 0 : m->SetMarkerSize(0.4);
283 :
284 : }
285 0 : }
286 :
287 :
288 : void AliMCEvent::ReorderAndExpandTreeTR()
289 : {
290 : //
291 : // Reorder and expand the track reference tree in order to match the kinematics tree.
292 : // Copy the information from different branches into one
293 : //
294 : // TreeTR
295 :
296 0 : fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
297 0 : fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
298 0 : if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
299 0 : fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 64000, 0);
300 :
301 :
302 : //
303 : // Activate the used branches only. Otherwisw we get a bad memory leak.
304 0 : if (fTreeTR) {
305 0 : fTreeTR->SetBranchStatus("*", 0);
306 0 : fTreeTR->SetBranchStatus("AliRun.*", 1);
307 0 : fTreeTR->SetBranchStatus("ITS.*", 1);
308 0 : fTreeTR->SetBranchStatus("TPC.*", 1);
309 0 : fTreeTR->SetBranchStatus("TRD.*", 1);
310 0 : fTreeTR->SetBranchStatus("TOF.*", 1);
311 0 : fTreeTR->SetBranchStatus("FRAME.*", 1);
312 0 : fTreeTR->SetBranchStatus("MUON.*", 1);
313 0 : }
314 :
315 : //
316 : // Connect the active branches
317 0 : TClonesArray* trefs[7];
318 0 : for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
319 0 : if (fTreeTR){
320 : // make branch for central track references
321 0 : if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
322 0 : if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
323 0 : if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
324 0 : if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
325 0 : if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
326 0 : if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
327 0 : if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
328 : }
329 :
330 0 : Int_t np = fStack->GetNprimary();
331 0 : Int_t nt = fTreeTR->GetEntries();
332 :
333 : //
334 : // Loop over tracks and find the secondaries with the help of the kine tree
335 : Int_t ifills = 0;
336 : Int_t it = 0;
337 : Int_t itlast = 0;
338 : TParticle* part;
339 :
340 0 : for (Int_t ip = np - 1; ip > -1; ip--) {
341 0 : part = fStack->Particle(ip);
342 : // printf("Particle %5d %5d %5d %5d %5d %5d \n",
343 : // ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
344 : // part->GetLastDaughter(), part->TestBit(kTransportBit));
345 :
346 : // Determine range of secondaries produced by this primary during transport
347 0 : Int_t dau1 = part->GetFirstDaughter();
348 0 : if (dau1 < np) continue; // This particle has no secondaries produced during transport
349 : Int_t dau2 = -1;
350 0 : if (dau1 > -1) {
351 0 : Int_t inext = ip - 1;
352 0 : while (dau2 < 0) {
353 0 : if (inext >= 0) {
354 0 : part = fStack->Particle(inext);
355 0 : dau2 = part->GetFirstDaughter();
356 0 : if (dau2 == -1 || dau2 < np) {
357 : dau2 = -1;
358 0 : } else {
359 0 : dau2--;
360 : }
361 : } else {
362 0 : dau2 = fStack->GetNtrack() - 1;
363 : }
364 0 : inext--;
365 : } // find upper bound
366 0 : } // dau2 < 0
367 :
368 :
369 : // printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
370 : //
371 : // Loop over reference hits and find secondary label
372 : // First the tricky part: find the entry in treeTR than contains the hits or
373 : // make sure that no hits exist.
374 : //
375 : Bool_t hasHits = kFALSE;
376 : Bool_t isOutside = kFALSE;
377 :
378 : it = itlast;
379 0 : while (!hasHits && !isOutside && it < nt) {
380 0 : fTreeTR->GetEntry(it++);
381 0 : for (Int_t ib = 0; ib < 7; ib++) {
382 0 : if (!trefs[ib]) continue;
383 0 : Int_t nh = trefs[ib]->GetEntries();
384 0 : for (Int_t ih = 0; ih < nh; ih++) {
385 0 : AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
386 0 : Int_t label = tr->Label();
387 0 : if (label >= dau1 && label <= dau2) {
388 : hasHits = kTRUE;
389 : itlast = it - 1;
390 0 : break;
391 : }
392 0 : if (label > dau2 || label < ip) {
393 : isOutside = kTRUE;
394 : itlast = it - 1;
395 0 : break;
396 : }
397 0 : } // hits
398 0 : if (hasHits || isOutside) break;
399 0 : } // branches
400 : } // entries
401 :
402 0 : if (!hasHits) {
403 : // Write empty entries
404 0 : for (Int_t id = dau1; (id <= dau2); id++) {
405 0 : fTmpTreeTR->Fill();
406 0 : ifills++;
407 : }
408 0 : } else {
409 : // Collect all hits
410 0 : fTreeTR->GetEntry(itlast);
411 0 : for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
412 0 : for (Int_t ib = 0; ib < 7; ib++) {
413 0 : if (!trefs[ib]) continue;
414 0 : Int_t nh = trefs[ib]->GetEntries();
415 0 : for (Int_t ih = 0; ih < nh; ih++) {
416 0 : AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
417 0 : Int_t label = tr->Label();
418 : // Skip primaries
419 0 : if (label == ip) continue;
420 0 : if (label > dau2 || label < dau1)
421 0 : printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
422 : itlast, label, dau1, dau2);
423 0 : if (label == id) {
424 : // secondary found
425 0 : tr->SetDetectorId(ib-1);
426 0 : Int_t nref = fTRBuffer->GetEntriesFast();
427 0 : TClonesArray &lref = *fTRBuffer;
428 0 : new(lref[nref]) AliTrackReference(*tr);
429 0 : }
430 0 : } // hits
431 0 : } // branches
432 0 : fTmpTreeTR->Fill();
433 0 : fTRBuffer->Delete();
434 0 : ifills++;
435 : } // daughters
436 : } // has hits
437 0 : } // tracks
438 :
439 : //
440 : // Now loop again and write the primaries
441 : //
442 0 : it = nt - 1;
443 0 : for (Int_t ip = 0; ip < np; ip++) {
444 : Int_t labmax = -1;
445 0 : while (labmax < ip && it > -1) {
446 0 : fTreeTR->GetEntry(it--);
447 0 : for (Int_t ib = 0; ib < 7; ib++) {
448 0 : if (!trefs[ib]) continue;
449 0 : Int_t nh = trefs[ib]->GetEntries();
450 : //
451 : // Loop over reference hits and find primary labels
452 0 : for (Int_t ih = 0; ih < nh; ih++) {
453 0 : AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
454 0 : Int_t label = tr->Label();
455 0 : if (label < np && label > labmax) {
456 : labmax = label;
457 0 : }
458 :
459 0 : if (label == ip) {
460 0 : tr->SetDetectorId(ib-1);
461 0 : Int_t nref = fTRBuffer->GetEntriesFast();
462 0 : TClonesArray &lref = *fTRBuffer;
463 0 : new(lref[nref]) AliTrackReference(*tr);
464 0 : }
465 : } // hits
466 0 : } // branches
467 : } // entries
468 0 : it++;
469 0 : fTmpTreeTR->Fill();
470 0 : fTRBuffer->Delete();
471 0 : ifills++;
472 : } // tracks
473 : // Check
474 :
475 :
476 : // Clean-up
477 0 : delete fTreeTR; fTreeTR = 0;
478 :
479 0 : for (Int_t ib = 0; ib < 7; ib++) {
480 0 : if (trefs[ib]) {
481 0 : trefs[ib]->Clear();
482 0 : delete trefs[ib];
483 0 : trefs[ib] = 0;
484 0 : }
485 : }
486 :
487 0 : if (ifills != fStack->GetNtrack())
488 0 : printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
489 0 : ifills, fStack->GetNtrack());
490 :
491 0 : fTmpTreeTR->Write();
492 0 : fTreeTR = fTmpTreeTR;
493 0 : }
494 :
495 : AliVParticle* AliMCEvent::GetTrack(Int_t i) const
496 : {
497 : // Get MC Particle i
498 : //
499 :
500 14454 : if (fExternal) {
501 0 : return ((AliVParticle*) (fMCParticles->At(i)));
502 : }
503 :
504 : //
505 : // Check first if this explicitely accesses the subsidiary event
506 :
507 7227 : if (i >= BgLabelOffset()) {
508 0 : if (fSubsidiaryEvents) {
509 0 : AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
510 0 : return (bgEvent->GetTrack(i - BgLabelOffset()));
511 : } else {
512 0 : return 0;
513 : }
514 : }
515 :
516 : //
517 : AliMCParticle *mcParticle = 0;
518 : TParticle *particle = 0;
519 : TClonesArray *trefs = 0;
520 : Int_t ntref = 0;
521 : TObjArray *rarray = 0;
522 :
523 :
524 :
525 : // Out of range check
526 14454 : if (i < 0 || i >= fNparticles) {
527 0 : AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
528 : mcParticle = 0;
529 0 : return (mcParticle);
530 : }
531 :
532 :
533 7227 : if (fSubsidiaryEvents) {
534 0 : AliMCEvent* mc;
535 0 : Int_t idx = FindIndexAndEvent(i, mc);
536 0 : return (mc->GetTrack(idx));
537 0 : }
538 :
539 : //
540 : // First check If the MC Particle has been already cached
541 7227 : if(!fMCParticleMap->At(i)) {
542 : // Get particle from the stack
543 2149 : particle = fStack->Particle(i);
544 : // Get track references from Tree TR
545 2149 : if (fTreeTR) {
546 2149 : fTreeTR->GetEntry(fStack->TreeKEntry(i));
547 2149 : trefs = fTRBuffer;
548 2149 : ntref = trefs->GetEntriesFast();
549 2149 : rarray = new TObjArray(ntref);
550 2149 : Int_t nen = fTrackReferences->GetEntriesFast();
551 13192 : for (Int_t j = 0; j < ntref; j++) {
552 : // Save the track references in a TClonesArray
553 13341 : AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
554 : // Save the pointer in a TRefArray
555 4447 : if (ref) {
556 4447 : new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
557 4447 : rarray->AddAt((*fTrackReferences)[nen], j);
558 4447 : nen++;
559 4447 : }
560 : } // loop over track references for entry i
561 2149 : } // if TreeTR available
562 2149 : Int_t nentries = fMCParticles->GetEntriesFast();
563 2149 : mcParticle = new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
564 2149 : fMCParticleMap->AddAt(mcParticle, i);
565 2149 : if (mcParticle) {
566 2149 : TParticle* part = mcParticle->Particle();
567 2149 : Int_t imo = part->GetFirstMother();
568 2149 : Int_t id1 = part->GetFirstDaughter();
569 2149 : Int_t id2 = part->GetLastDaughter();
570 4298 : if (fPrimaryOffset > 0 || fSecondaryOffset > 0) {
571 : // Remapping of the mother and daughter indices
572 0 : if (imo < fNprimaries) {
573 0 : mcParticle->SetMother(imo + fPrimaryOffset);
574 0 : } else {
575 0 : mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
576 : }
577 :
578 0 : if (id1 < fNprimaries) {
579 0 : mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
580 0 : mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
581 0 : } else {
582 0 : mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
583 0 : mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
584 : }
585 :
586 :
587 0 : if (i > fNprimaries) {
588 0 : mcParticle->SetLabel(i + fPrimaryOffset);
589 0 : } else {
590 0 : mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
591 : }
592 : } else {
593 2149 : mcParticle->SetFirstDaughter(id1);
594 2149 : mcParticle->SetLastDaughter (id2);
595 2149 : mcParticle->SetMother (imo);
596 : }
597 2149 : }
598 2149 : } else {
599 15234 : mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
600 : }
601 :
602 : //Printf("mcParticleGetMother %d",mcParticle->GetMother());
603 7227 : return mcParticle;
604 7227 : }
605 :
606 : AliGenEventHeader* AliMCEvent::GenEventHeader() const
607 : {
608 16 : if (!fExternal) {
609 : // ESD
610 8 : return (fHeader->GenEventHeader());
611 : } else {
612 : // AOD
613 0 : if (fAODMCHeader) {
614 0 : TList * lh = fAODMCHeader->GetCocktailHeaders();
615 0 : if (lh) {return ((AliGenEventHeader*) lh->At(0));}
616 0 : }
617 : }
618 0 : return 0;
619 8 : }
620 :
621 :
622 : void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event)
623 : {
624 : // Add a subsidiary event to the list; for example merged background event.
625 0 : if (!fSubsidiaryEvents) {
626 0 : TList* events = new TList();
627 0 : events->Add(new AliMCEvent(*this));
628 0 : fSubsidiaryEvents = events;
629 0 : }
630 :
631 0 : fSubsidiaryEvents->Add(event);
632 0 : }
633 :
634 : AliGenEventHeader *AliMCEvent::FindHeader(Int_t ipart) {
635 : //
636 : // Get Header belonging to this track;
637 : // only works for primaries (i.e. particles coming from the Generator)
638 : // Also sorts out the case of Cocktail event (get header of subevent in cocktail generetor header)
639 : //
640 :
641 0 : AliMCEvent *event = this;
642 :
643 0 : if (fSubsidiaryEvents) {
644 : // Get pointer to subevent if needed
645 0 : ipart = FindIndexAndEvent(ipart,event);
646 0 : }
647 :
648 0 : AliGenEventHeader* header = event->GenEventHeader();
649 0 : if (ipart >= header->NProduced()) {
650 0 : AliWarning(Form("Not a primary -- returning 0 (idx %d, nPrimary %d)",ipart,header->NProduced()));
651 0 : return 0;
652 : }
653 0 : AliGenCocktailEventHeader *coHeader = dynamic_cast<AliGenCocktailEventHeader*>(header);
654 0 : if (coHeader) { // Cocktail event
655 0 : TList* headerList = coHeader->GetHeaders();
656 0 : TIter headIt(headerList);
657 : Int_t nproduced = 0;
658 0 : do { // Go trhough all headers and look for the correct one
659 0 : header = (AliGenEventHeader*) headIt();
660 0 : if (header) nproduced += header->NProduced();
661 0 : } while (header && ipart >= nproduced);
662 0 : }
663 :
664 : return header;
665 0 : }
666 :
667 : Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
668 : {
669 : // Find the index and event in case of composed events like signal + background
670 0 : TIter next(fSubsidiaryEvents);
671 0 : next.Reset();
672 0 : if (oldidx < fNprimaries) {
673 0 : while((event = (AliMCEvent*)next())) {
674 0 : if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
675 : }
676 0 : if (event) {
677 0 : return (oldidx - event->GetPrimaryOffset());
678 : } else {
679 0 : return (-1);
680 : }
681 : } else {
682 0 : while((event = (AliMCEvent*)next())) {
683 0 : if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
684 : }
685 0 : if (event) {
686 0 : return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
687 : } else {
688 0 : return (-1);
689 : }
690 : }
691 0 : }
692 :
693 : Int_t AliMCEvent::BgLabelToIndex(Int_t label)
694 : {
695 : // Convert a background label to an absolute index
696 0 : if (fSubsidiaryEvents) {
697 0 : AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
698 0 : label -= BgLabelOffset();
699 0 : if (label < bgEvent->GetNumberOfPrimaries()) {
700 0 : label += bgEvent->GetPrimaryOffset();
701 0 : } else {
702 0 : label += (bgEvent->GetSecondaryOffset() - fNprimaries);
703 : }
704 0 : }
705 0 : return (label);
706 : }
707 :
708 :
709 : Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) const
710 : {
711 : //
712 : // Delegate to subevent if necesarry
713 :
714 :
715 438 : if (!fSubsidiaryEvents) {
716 219 : return fStack->IsPhysicalPrimary(i);
717 : } else {
718 0 : AliMCEvent* evt = 0;
719 0 : Int_t idx = FindIndexAndEvent(i, evt);
720 0 : return (evt->IsPhysicalPrimary(idx));
721 0 : }
722 219 : }
723 :
724 : Bool_t AliMCEvent::IsSecondaryFromWeakDecay(Int_t i)
725 : {
726 : //
727 : // Delegate to subevent if necesarry
728 438 : if (!fSubsidiaryEvents) {
729 219 : return fStack->IsSecondaryFromWeakDecay(i);
730 : } else {
731 0 : AliMCEvent* evt = 0;
732 0 : Int_t idx = FindIndexAndEvent(i, evt);
733 0 : return (evt->IsSecondaryFromWeakDecay(idx));
734 0 : }
735 219 : }
736 :
737 : Bool_t AliMCEvent::IsSecondaryFromMaterial(Int_t i)
738 : {
739 : //
740 : // Delegate to subevent if necesarry
741 438 : if (!fSubsidiaryEvents) {
742 219 : return fStack->IsSecondaryFromMaterial(i);
743 : } else {
744 0 : AliMCEvent* evt = 0;
745 0 : Int_t idx = FindIndexAndEvent(i, evt);
746 0 : return (evt->IsSecondaryFromMaterial(idx));
747 0 : }
748 219 : }
749 :
750 :
751 : void AliMCEvent::InitEvent()
752 : {
753 : //
754 : // Initialize the subsidiary event structure
755 0 : if (fSubsidiaryEvents) {
756 0 : TIter next(fSubsidiaryEvents);
757 : AliMCEvent* evt;
758 0 : fNprimaries = 0;
759 0 : fNparticles = 0;
760 :
761 0 : while((evt = (AliMCEvent*)next())) {
762 0 : fNprimaries += evt->GetNumberOfPrimaries();
763 0 : fNparticles += evt->GetNumberOfTracks();
764 : }
765 :
766 : Int_t ioffp = 0;
767 0 : Int_t ioffs = fNprimaries;
768 0 : next.Reset();
769 :
770 0 : while((evt = (AliMCEvent*)next())) {
771 0 : evt->SetPrimaryOffset(ioffp);
772 0 : evt->SetSecondaryOffset(ioffs);
773 0 : ioffp += evt->GetNumberOfPrimaries();
774 0 : ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
775 : }
776 0 : }
777 0 : }
778 :
779 : void AliMCEvent::PreReadAll()
780 : {
781 : // Preread the MC information
782 : Int_t i;
783 : // secondaries
784 4086 : for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
785 : {
786 2037 : GetTrack(i);
787 : }
788 : // primaries
789 232 : for (i = 0; i < fStack->GetNprimary(); i++)
790 : {
791 112 : GetTrack(i);
792 : }
793 4 : AssignGeneratorIndex();
794 4 : }
795 :
796 : const AliVVertex * AliMCEvent::GetPrimaryVertex() const
797 : {
798 : // Create a MCVertex object from the MCHeader information
799 0 : TArrayF v;
800 0 : GenEventHeader()->PrimaryVertex(v) ;
801 0 : if (!fVertex) {
802 0 : fVertex = new AliMCVertex(v[0], v[1], v[2]);
803 0 : } else {
804 0 : ((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);
805 : }
806 0 : return fVertex;
807 0 : }
808 :
809 : Bool_t AliMCEvent::IsFromBGEvent(Int_t index)
810 : {
811 : // Checks if a particle is from the background events
812 : // Works for HIJING inside Cocktail
813 0 : if (fNBG == -1) {
814 : AliGenCocktailEventHeader* coHeader =
815 0 : dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
816 0 : if (!coHeader) return (0);
817 0 : TList* list = coHeader->GetHeaders();
818 0 : AliGenHijingEventHeader* hijingH = dynamic_cast<AliGenHijingEventHeader*>(list->FindObject("Hijing"));
819 0 : if (!hijingH) return (0);
820 0 : fNBG = hijingH->NProduced();
821 0 : }
822 :
823 0 : return (index < fNBG);
824 0 : }
825 :
826 :
827 : TList* AliMCEvent::GetCocktailList()
828 : {
829 : //gives the CocktailHeaders when reading ESDs/AODs (corresponding to fExteral=kFALSE/kTRUE)
830 : //the AODMC header (and the aodmc array) is passed as an instance to MCEvent by the AliAODInputHandler
831 8 : if(fExternal==kFALSE) {
832 12 : AliGenCocktailEventHeader* coHeader =dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
833 4 : if(!coHeader) {
834 0 : return 0;
835 : } else {
836 4 : return (coHeader->GetHeaders());
837 : }
838 : } else {
839 0 : if(!fAODMCHeader) {
840 0 : return 0;
841 : } else {
842 0 : return (fAODMCHeader->GetCocktailHeaders());
843 : }
844 : }
845 4 : }
846 :
847 :
848 : TString AliMCEvent::GetGenerator(Int_t index)
849 : {
850 0 : Int_t nsumpart=fNprimaries;
851 0 : TList* lh = GetCocktailList();
852 0 : if(!lh){ TString noheader="nococktailheader";
853 0 : return noheader;}
854 0 : Int_t nh=lh->GetEntries();
855 0 : for (Int_t i = nh-1; i >= 0; i--){
856 0 : AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
857 0 : TString genname=gh->GetName();
858 0 : Int_t npart=gh->NProduced();
859 0 : if (i == 0) npart = nsumpart;
860 0 : if(index < nsumpart && index >= (nsumpart-npart)) return genname;
861 0 : nsumpart-=npart;
862 0 : }
863 0 : TString empty="";
864 0 : return empty;
865 0 : }
866 :
867 : void AliMCEvent::AssignGeneratorIndex() {
868 : //
869 : // Assign the generator index to each particle
870 : //
871 8 : TList* list = GetCocktailList();
872 4 : if (fNprimaries <= 0) {
873 0 : AliWarning(Form("AliMCEvent::AssignGeneratorIndex: no primaries %10d\n", fNprimaries));
874 0 : return;
875 : }
876 4 : if (!list) {
877 0 : return;
878 : } else {
879 4 : Int_t nh = list->GetEntries();
880 4 : Int_t nsumpart = fNprimaries;
881 40 : for(Int_t i = nh-1; i >= 0; i--){
882 16 : AliGenEventHeader* gh = (AliGenEventHeader*)list->At(i);
883 16 : Int_t npart = gh->NProduced();
884 16 : if (i==0) {
885 : if (npart != nsumpart) {
886 : // printf("Header inconsistent ! %5d %5d \n", npart, nsumpart);
887 : }
888 : npart = nsumpart;
889 4 : }
890 : //
891 : // Loop over primary particles for generator i
892 272 : for (Int_t j = nsumpart-1; j >= nsumpart-npart; j--) {
893 112 : AliVParticle* part = GetTrack(j);
894 112 : if (!part) {
895 0 : AliWarning(Form("AliMCEvent::AssignGeneratorIndex: 0-pointer to particle j %8d npart %8d nsumpart %8d Nprimaries %8d\n",
896 : j, npart, nsumpart, fNprimaries));
897 0 : break;
898 : }
899 112 : part->SetGeneratorIndex(i);
900 112 : Int_t dmin = part->GetFirstDaughter();
901 112 : Int_t dmax = part->GetLastDaughter();
902 135 : if (dmin == -1) continue;
903 89 : AssignGeneratorIndex(i, dmin, dmax);
904 89 : }
905 : nsumpart -= npart;
906 : }
907 : }
908 8 : }
909 : void AliMCEvent::AssignGeneratorIndex(Int_t index, Int_t dmin, Int_t dmax) {
910 6936 : for (Int_t k = dmin; k <= dmax; k++) {
911 2037 : AliVParticle* dpart = GetTrack(k);
912 2037 : dpart->SetGeneratorIndex(index);
913 2037 : Int_t d1 = dpart->GetFirstDaughter();
914 2037 : Int_t d2 = dpart->GetLastDaughter();
915 2037 : if (d1 > -1) {
916 865 : AssignGeneratorIndex(index, d1, d2);
917 865 : }
918 : }
919 954 : }
920 :
921 : Bool_t AliMCEvent::GetCocktailGenerator(Int_t index,TString &nameGen){
922 : //method that gives the generator for a given particle with label index (or that of the corresponding primary)
923 0 : AliVParticle* mcpart0 = (AliVParticle*) (GetTrack(index));
924 0 : if(!mcpart0){
925 0 : printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",index);
926 0 : return 0;
927 : }
928 : /*
929 : Int_t ig = mcpart0->GetGeneratorIndex();
930 : if (ig != -1) {
931 : nameGen = ((AliGenEventHeader*)GetCocktailList()->At(ig))->GetName();
932 : return 1;
933 : }
934 : */
935 0 : nameGen=GetGenerator(index);
936 0 : if(nameGen.Contains("nococktailheader") )return 0;
937 : Int_t lab=index;
938 :
939 0 : while(nameGen.IsWhitespace()){
940 :
941 :
942 0 : AliVParticle* mcpart = (AliVParticle*) (GetTrack(lab));
943 :
944 0 : if(!mcpart){
945 0 : printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",lab);
946 0 : break;}
947 : Int_t mother=0;
948 0 : mother = mcpart->GetMother();
949 :
950 0 : if(mother<0){
951 0 : printf("AliMCEvent - BREAK: Reached primary particle without valid mother\n");
952 0 : break;
953 : }
954 0 : AliVParticle* mcmom = (AliVParticle*) (GetTrack(mother));
955 0 : if(!mcmom){
956 0 : printf("AliMCEvent-BREAK: No valid AliMCParticle mother at label %i\n",mother);
957 0 : break;
958 : }
959 : lab=mother;
960 :
961 0 : nameGen=GetGenerator(mother);
962 0 : }
963 :
964 : return 1;
965 0 : }
966 :
967 : void AliMCEvent::SetParticleArray(TClonesArray* mcParticles)
968 : {
969 0 : fMCParticles = mcParticles;
970 0 : fNparticles = fMCParticles->GetEntries();
971 0 : fExternal = kTRUE;
972 0 : fNprimaries = 0;
973 : struct Local {
974 : static Int_t binaryfirst(TClonesArray* a, Int_t low, Int_t high)
975 : {
976 0 : Int_t mid = low + (high - low)/2;
977 0 : if (low > a->GetEntries()-1) return (a->GetEntries()-1);
978 0 : if (!((AliVParticle*) a->At(mid))->IsPrimary()) {
979 0 : if (mid > 1 && !((AliVParticle*) a->At(mid-1))->IsPrimary()) {
980 0 : return binaryfirst(a, low, mid-1);
981 : } else {
982 0 : return mid;
983 : }
984 : } else {
985 0 : return binaryfirst(a, mid+1, high);
986 : }
987 0 : }
988 : };
989 0 : fNprimaries = Local::binaryfirst(mcParticles, 0, mcParticles->GetEntries()-1);
990 0 : AssignGeneratorIndex();
991 0 : }
992 :
993 : AliVEvent::EDataLayoutType AliMCEvent::GetDataLayoutType() const
994 : {
995 0 : return AliVEvent::kMC;
996 : }
997 :
998 176 : ClassImp(AliMCEvent)
|
