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