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