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