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