]> 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
add cluster occupancy to the ESD friend
[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{
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;
69fcabe8 499 TObjArray *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();
69fcabe8 527 rarray = new TObjArray(ntref);
93df0e9b 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
95121282 597AliGenEventHeader *AliMCEvent::FindHeader(Int_t ipart) {
598 //
599 // Get Header belonging to this track;
600 // only works for primaries (i.e. particles coming from the Generator)
601 // Also sorts out the case of Cocktail event (get header of subevent in cocktail generetor header)
602 //
603
604 AliMCEvent *event = this;
605
606 if (fSubsidiaryEvents) {
607 // Get pointer to subevent if needed
608 ipart = FindIndexAndEvent(ipart,event);
609 }
610
611 AliGenEventHeader* header = event->GenEventHeader();
612 if (ipart >= header->NProduced()) {
613 AliWarning(Form("Not a primary -- returning 0 (idx %d, nPrimary %d)",ipart,header->NProduced()));
614 return 0;
615 }
616 AliGenCocktailEventHeader *coHeader = dynamic_cast<AliGenCocktailEventHeader*>(header);
617 if (coHeader) { // Cocktail event
618 TList* headerList = coHeader->GetHeaders();
619 TIter headIt(headerList);
620 Int_t nproduced = 0;
621 do { // Go trhough all headers and look for the correct one
622 header = (AliGenEventHeader*) headIt();
4636bbd8 623 if (header) nproduced += header->NProduced();
95121282 624 } while (header && ipart >= nproduced);
625 }
626
627 return header;
628}
629
93836e1b 630Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
631{
632 // Find the index and event in case of composed events like signal + background
633 TIter next(fSubsidiaryEvents);
634 next.Reset();
635 if (oldidx < fNprimaries) {
636 while((event = (AliMCEvent*)next())) {
637 if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
638 }
58c338d2 639 if (event) {
640 return (oldidx - event->GetPrimaryOffset());
641 } else {
642 return (-1);
643 }
93836e1b 644 } else {
645 while((event = (AliMCEvent*)next())) {
646 if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
647 }
58c338d2 648 if (event) {
649 return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
650 } else {
651 return (-1);
652 }
93836e1b 653 }
654}
655
656Int_t AliMCEvent::BgLabelToIndex(Int_t label)
657{
658 // Convert a background label to an absolute index
659 if (fSubsidiaryEvents) {
660 AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
661 label -= BgLabelOffset();
662 if (label < bgEvent->GetNumberOfPrimaries()) {
663 label += bgEvent->GetPrimaryOffset();
664 } else {
665 label += (bgEvent->GetSecondaryOffset() - fNprimaries);
666 }
667 }
668 return (label);
669}
670
671
69fcabe8 672Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) const
93836e1b 673{
674//
675// Delegate to subevent if necesarry
676
677
678 if (!fSubsidiaryEvents) {
679 return fStack->IsPhysicalPrimary(i);
680 } else {
681 AliMCEvent* evt = 0;
682 Int_t idx = FindIndexAndEvent(i, evt);
683 return (evt->IsPhysicalPrimary(idx));
684 }
685}
686
1081342d 687Bool_t AliMCEvent::IsSecondaryFromWeakDecay(Int_t i)
688{
689//
690// Delegate to subevent if necesarry
691 if (!fSubsidiaryEvents) {
692 return fStack->IsSecondaryFromWeakDecay(i);
693 } else {
694 AliMCEvent* evt = 0;
695 Int_t idx = FindIndexAndEvent(i, evt);
696 return (evt->IsSecondaryFromWeakDecay(idx));
697 }
698}
699
700Bool_t AliMCEvent::IsSecondaryFromMaterial(Int_t i)
701{
702//
703// Delegate to subevent if necesarry
704 if (!fSubsidiaryEvents) {
705 return fStack->IsSecondaryFromMaterial(i);
706 } else {
707 AliMCEvent* evt = 0;
708 Int_t idx = FindIndexAndEvent(i, evt);
709 return (evt->IsSecondaryFromMaterial(idx));
710 }
711}
712
713
93836e1b 714void AliMCEvent::InitEvent()
715{
35152e4b 716//
717// Initialize the subsidiary event structure
93836e1b 718 if (fSubsidiaryEvents) {
719 TIter next(fSubsidiaryEvents);
720 AliMCEvent* evt;
721 fNprimaries = 0;
722 fNparticles = 0;
723
724 while((evt = (AliMCEvent*)next())) {
725 fNprimaries += evt->GetNumberOfPrimaries();
726 fNparticles += evt->GetNumberOfTracks();
727 }
728
729 Int_t ioffp = 0;
730 Int_t ioffs = fNprimaries;
731 next.Reset();
732
733 while((evt = (AliMCEvent*)next())) {
734 evt->SetPrimaryOffset(ioffp);
735 evt->SetSecondaryOffset(ioffs);
736 ioffp += evt->GetNumberOfPrimaries();
737 ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
738 }
739 }
740}
091fa189 741
35152e4b 742void AliMCEvent::PreReadAll()
743{
744 // Preread the MC information
745 Int_t i;
746 // secondaries
747 for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
748 {
749 GetTrack(i);
750 }
751 // primaries
752 for (i = 0; i < fStack->GetNprimary(); i++)
753 {
754 GetTrack(i);
755 }
35152e4b 756}
93836e1b 757
5df3496b 758const AliVVertex * AliMCEvent::GetPrimaryVertex() const
759{
760 // Create a MCVertex object from the MCHeader information
570308e4 761 TArrayF v;
762 GenEventHeader()->PrimaryVertex(v) ;
5df3496b 763 if (!fVertex) {
5df3496b 764 fVertex = new AliMCVertex(v[0], v[1], v[2]);
570308e4 765 } else {
766 ((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);
5df3496b 767 }
5df3496b 768 return fVertex;
769}
770
2e7ca177 771Bool_t AliMCEvent::IsFromBGEvent(Int_t index)
772{
773 // Checks if a particle is from the background events
774 // Works for HIJING inside Cocktail
775 if (fNBG == -1) {
776 AliGenCocktailEventHeader* coHeader =
777 dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
778 if (!coHeader) return (0);
779 TList* list = coHeader->GetHeaders();
780 AliGenHijingEventHeader* hijingH = dynamic_cast<AliGenHijingEventHeader*>(list->FindObject("Hijing"));
781 if (!hijingH) return (0);
782 fNBG = hijingH->NProduced();
783 }
784
785 return (index < fNBG);
786}
787
5df3496b 788
415d9f5c 789ClassImp(AliMCEvent)