]> 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
Temporary splines for LHC12b-e anchored MC
[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) {
4e0b0eb8 501
502
7aad0c47 503 return ((AliVParticle*) (fMCParticles->At(i)));
4e0b0eb8 504
7aad0c47 505 }
506
93836e1b 507 //
508 // Check first if this explicitely accesses the subsidiary event
7aad0c47 509
93836e1b 510 if (i >= BgLabelOffset()) {
511 if (fSubsidiaryEvents) {
512 AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
513 return (bgEvent->GetTrack(i - BgLabelOffset()));
514 } else {
515 return 0;
516 }
517 }
518
519 //
93df0e9b 520 AliMCParticle *mcParticle = 0;
521 TParticle *particle = 0;
522 TClonesArray *trefs = 0;
523 Int_t ntref = 0;
69fcabe8 524 TObjArray *rarray = 0;
93836e1b 525
526
527
415d9f5c 528 // Out of range check
529 if (i < 0 || i >= fNparticles) {
530 AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
531 mcParticle = 0;
532 return (mcParticle);
533 }
93836e1b 534
415d9f5c 535
93836e1b 536 if (fSubsidiaryEvents) {
537 AliMCEvent* mc;
538 Int_t idx = FindIndexAndEvent(i, mc);
539 return (mc->GetTrack(idx));
540 }
415d9f5c 541
415d9f5c 542 //
93836e1b 543 // First check If the MC Particle has been already cached
415d9f5c 544 if(!fMCParticleMap->At(i)) {
93df0e9b 545 // Get particle from the stack
546 particle = fStack->Particle(i);
547 // Get track references from Tree TR
548 if (fTreeTR) {
549 fTreeTR->GetEntry(fStack->TreeKEntry(i));
550 trefs = fTRBuffer;
551 ntref = trefs->GetEntriesFast();
69fcabe8 552 rarray = new TObjArray(ntref);
93df0e9b 553 Int_t nen = fTrackReferences->GetEntriesFast();
554 for (Int_t j = 0; j < ntref; j++) {
555 // Save the track references in a TClonesArray
556 AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
557 // Save the pointer in a TRefArray
58c338d2 558 if (ref) {
559 new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
560 rarray->AddAt((*fTrackReferences)[nen], j);
561 nen++;
562 }
93df0e9b 563 } // loop over track references for entry i
564 } // if TreeTR available
415d9f5c 565 Int_t nentries = fMCParticles->GetEntriesFast();
091fa189 566 mcParticle = new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
415d9f5c 567 fMCParticleMap->AddAt(mcParticle, i);
58c338d2 568 if (mcParticle) {
569 TParticle* part = mcParticle->Particle();
570 Int_t imo = part->GetFirstMother();
571 Int_t id1 = part->GetFirstDaughter();
572 Int_t id2 = part->GetLastDaughter();
573 if (fPrimaryOffset > 0 || fSecondaryOffset > 0) {
574 // Remapping of the mother and daughter indices
575 if (imo < fNprimaries) {
576 mcParticle->SetMother(imo + fPrimaryOffset);
577 } else {
578 mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
579 }
580
581 if (id1 < fNprimaries) {
582 mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
583 mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
584 } else {
585 mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
586 mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
587 }
588
589
590 if (i > fNprimaries) {
591 mcParticle->SetLabel(i + fPrimaryOffset);
592 } else {
593 mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
594 }
93836e1b 595 } else {
58c338d2 596 mcParticle->SetFirstDaughter(id1);
597 mcParticle->SetLastDaughter (id2);
598 mcParticle->SetMother (imo);
93836e1b 599 }
93836e1b 600 }
415d9f5c 601 } else {
602 mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
603 }
4e0b0eb8 604
605 //Printf("mcParticleGetMother %d",mcParticle->GetMother());
415d9f5c 606 return mcParticle;
607}
608
5df3496b 609AliGenEventHeader* AliMCEvent::GenEventHeader() const {return (fHeader->GenEventHeader());}
93836e1b 610
611
612void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event)
613{
614 // Add a subsidiary event to the list; for example merged background event.
615 if (!fSubsidiaryEvents) {
616 TList* events = new TList();
617 events->Add(new AliMCEvent(*this));
618 fSubsidiaryEvents = events;
619 }
620
621 fSubsidiaryEvents->Add(event);
622}
623
95121282 624AliGenEventHeader *AliMCEvent::FindHeader(Int_t ipart) {
625 //
626 // Get Header belonging to this track;
627 // only works for primaries (i.e. particles coming from the Generator)
628 // Also sorts out the case of Cocktail event (get header of subevent in cocktail generetor header)
629 //
630
631 AliMCEvent *event = this;
632
633 if (fSubsidiaryEvents) {
634 // Get pointer to subevent if needed
635 ipart = FindIndexAndEvent(ipart,event);
636 }
637
638 AliGenEventHeader* header = event->GenEventHeader();
639 if (ipart >= header->NProduced()) {
640 AliWarning(Form("Not a primary -- returning 0 (idx %d, nPrimary %d)",ipart,header->NProduced()));
641 return 0;
642 }
643 AliGenCocktailEventHeader *coHeader = dynamic_cast<AliGenCocktailEventHeader*>(header);
644 if (coHeader) { // Cocktail event
645 TList* headerList = coHeader->GetHeaders();
646 TIter headIt(headerList);
647 Int_t nproduced = 0;
648 do { // Go trhough all headers and look for the correct one
649 header = (AliGenEventHeader*) headIt();
4636bbd8 650 if (header) nproduced += header->NProduced();
95121282 651 } while (header && ipart >= nproduced);
652 }
653
654 return header;
655}
656
93836e1b 657Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
658{
659 // Find the index and event in case of composed events like signal + background
660 TIter next(fSubsidiaryEvents);
661 next.Reset();
662 if (oldidx < fNprimaries) {
663 while((event = (AliMCEvent*)next())) {
664 if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
665 }
58c338d2 666 if (event) {
667 return (oldidx - event->GetPrimaryOffset());
668 } else {
669 return (-1);
670 }
93836e1b 671 } else {
672 while((event = (AliMCEvent*)next())) {
673 if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
674 }
58c338d2 675 if (event) {
676 return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
677 } else {
678 return (-1);
679 }
93836e1b 680 }
681}
682
683Int_t AliMCEvent::BgLabelToIndex(Int_t label)
684{
685 // Convert a background label to an absolute index
686 if (fSubsidiaryEvents) {
687 AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
688 label -= BgLabelOffset();
689 if (label < bgEvent->GetNumberOfPrimaries()) {
690 label += bgEvent->GetPrimaryOffset();
691 } else {
692 label += (bgEvent->GetSecondaryOffset() - fNprimaries);
693 }
694 }
695 return (label);
696}
697
698
69fcabe8 699Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) const
93836e1b 700{
701//
702// Delegate to subevent if necesarry
703
704
705 if (!fSubsidiaryEvents) {
706 return fStack->IsPhysicalPrimary(i);
707 } else {
708 AliMCEvent* evt = 0;
709 Int_t idx = FindIndexAndEvent(i, evt);
710 return (evt->IsPhysicalPrimary(idx));
711 }
712}
713
1081342d 714Bool_t AliMCEvent::IsSecondaryFromWeakDecay(Int_t i)
715{
716//
717// Delegate to subevent if necesarry
718 if (!fSubsidiaryEvents) {
719 return fStack->IsSecondaryFromWeakDecay(i);
720 } else {
721 AliMCEvent* evt = 0;
722 Int_t idx = FindIndexAndEvent(i, evt);
723 return (evt->IsSecondaryFromWeakDecay(idx));
724 }
725}
726
727Bool_t AliMCEvent::IsSecondaryFromMaterial(Int_t i)
728{
729//
730// Delegate to subevent if necesarry
731 if (!fSubsidiaryEvents) {
732 return fStack->IsSecondaryFromMaterial(i);
733 } else {
734 AliMCEvent* evt = 0;
735 Int_t idx = FindIndexAndEvent(i, evt);
736 return (evt->IsSecondaryFromMaterial(idx));
737 }
738}
739
740
93836e1b 741void AliMCEvent::InitEvent()
742{
35152e4b 743//
744// Initialize the subsidiary event structure
93836e1b 745 if (fSubsidiaryEvents) {
746 TIter next(fSubsidiaryEvents);
747 AliMCEvent* evt;
748 fNprimaries = 0;
749 fNparticles = 0;
750
751 while((evt = (AliMCEvent*)next())) {
752 fNprimaries += evt->GetNumberOfPrimaries();
753 fNparticles += evt->GetNumberOfTracks();
754 }
755
756 Int_t ioffp = 0;
757 Int_t ioffs = fNprimaries;
758 next.Reset();
759
760 while((evt = (AliMCEvent*)next())) {
761 evt->SetPrimaryOffset(ioffp);
762 evt->SetSecondaryOffset(ioffs);
763 ioffp += evt->GetNumberOfPrimaries();
764 ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
765 }
766 }
767}
091fa189 768
35152e4b 769void AliMCEvent::PreReadAll()
770{
771 // Preread the MC information
772 Int_t i;
773 // secondaries
774 for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
775 {
776 GetTrack(i);
777 }
778 // primaries
779 for (i = 0; i < fStack->GetNprimary(); i++)
780 {
781 GetTrack(i);
782 }
35152e4b 783}
93836e1b 784
5df3496b 785const AliVVertex * AliMCEvent::GetPrimaryVertex() const
786{
787 // Create a MCVertex object from the MCHeader information
570308e4 788 TArrayF v;
789 GenEventHeader()->PrimaryVertex(v) ;
5df3496b 790 if (!fVertex) {
5df3496b 791 fVertex = new AliMCVertex(v[0], v[1], v[2]);
570308e4 792 } else {
793 ((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);
5df3496b 794 }
5df3496b 795 return fVertex;
796}
797
2e7ca177 798Bool_t AliMCEvent::IsFromBGEvent(Int_t index)
799{
800 // Checks if a particle is from the background events
801 // Works for HIJING inside Cocktail
802 if (fNBG == -1) {
803 AliGenCocktailEventHeader* coHeader =
804 dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
805 if (!coHeader) return (0);
806 TList* list = coHeader->GetHeaders();
807 AliGenHijingEventHeader* hijingH = dynamic_cast<AliGenHijingEventHeader*>(list->FindObject("Hijing"));
808 if (!hijingH) return (0);
809 fNBG = hijingH->NProduced();
810 }
811
812 return (index < fNBG);
813}
814
5df3496b 815
4e0b0eb8 816 Int_t AliMCEvent::GetCocktailList(TList*& lh){
817 //gives the CocktailHeaders when reading ESDs/AODs (corresponding to fExteral=kFALSE/kTRUE)
818 //the AODMC header (and the aodmc array) is passed as an instance to MCEvent by the AliAODInputHandler
819 if(fExternal==kFALSE) {
820 AliGenCocktailEventHeader* coHeader =dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
821 if(!coHeader) return 0;
822 lh=coHeader->GetHeaders();}
823 if(fExternal==kTRUE){
824 if(!fAODMCHeader) return 0;
825 lh=fAODMCHeader->GetCocktailHeaders();}
826 return 1; }
16632e6b 827
828
829
4e0b0eb8 830
831 TString AliMCEvent::GetGenerator(Int_t index){
832 Int_t nsumpart=0;
833
834 TList* lh;
835 Int_t nt= GetCocktailList(lh);
836 if(nt==0){ TString noheader="nococktailheader";
837 return noheader;}
838 Int_t nh=lh->GetEntries();
839 for(Int_t i=0;i<nh;i++){
840 AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
841 TString genname=gh->GetName();
842 Int_t npart=gh->NProduced();
843 if(index>=nsumpart && index<(nsumpart+npart)) return genname;
844 nsumpart+=npart;}
16632e6b 845 TString empty="";
4e0b0eb8 846 return empty;}
847
11762dd0 848void AliMCEvent::AssignGeneratorIndex() {
849 //
850 // Assign the generator index to each particle
851 //
852 TList* list;
853 Int_t nt = GetCocktailList(list);
854 if (nt == 0) {
855 } else {
856 Int_t nh = list->GetEntries();
857 Int_t nsumpart = 0;
858 for(Int_t i = 0; i < nh; i++){
859 AliGenEventHeader* gh = (AliGenEventHeader*)list->At(i);
860 Int_t npart = gh->NProduced();
861 for (Int_t j = nsumpart; j < npart; j++) {
862 AliVParticle* part = GetTrack(j);
863 part->SetGeneratorIndex(i);
864 Int_t dmin = part->GetFirstDaughter();
865 Int_t dmax = part->GetLastDaughter();
866 for (Int_t k = dmin; k <= dmax; k++) {
867 AliVParticle* dpart = GetTrack(k);
868 dpart->SetGeneratorIndex(i);
869 }
870 }
871 nsumpart += npart;
872 }
873 }
874}
4e0b0eb8 875
16632e6b 876
4e0b0eb8 877 Bool_t AliMCEvent::GetCocktailGenerator(Int_t index,TString &nameGen){
16632e6b 878 //method that gives the generator for a given particle with label index (or that of the corresponding primary)
4e0b0eb8 879
880 nameGen=GetGenerator(index);
881 if(nameGen.Contains("nococktailheader") )return 0;
882 Int_t lab=index;
883
884 while(nameGen.IsWhitespace()){
885
886
887 AliVParticle* mcpart = (AliVParticle*) (GetTrack(lab));
888
889 if(!mcpart){
890 printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",lab);
891 break;}
892 Int_t mother=0;
893 mother = mcpart->GetMother();
d597de9e 894
16632e6b 895 if(mother<0){
896 printf("AliMCEvent - BREAK: Reached primary particle without valid mother\n");
897 break;
898 }
4e0b0eb8 899 AliVParticle* mcmom = (AliVParticle*) (GetTrack(mother));
900 if(!mcmom){
901 printf("AliMCEvent-BREAK: No valid AliMCParticle mother at label %i\n",mother);
902 break;
903 }
904 lab=mother;
905
16632e6b 906 nameGen=GetGenerator(mother);
907 }
908
909 return 1;
4e0b0eb8 910 }
911
16632e6b 912
913
914
915
415d9f5c 916ClassImp(AliMCEvent)