1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-------------------------------------------------------------------------
19 // Class for Kinematic Events
20 // Author: Andreas Morsch, CERN
21 //-------------------------------------------------------------------------
27 #include <TParticle.h>
28 #include <TClonesArray.h>
29 #include <TRefArray.h>
34 #include "AliMCEvent.h"
35 #include "AliMCVertex.h"
37 #include "AliTrackReference.h"
38 #include "AliHeader.h"
39 #include "AliGenEventHeader.h"
42 Int_t AliMCEvent::fgkBgLabelOffset(10000000);
45 AliMCEvent::AliMCEvent():
50 fHeader(new AliHeader()),
52 fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
64 // Default constructor
67 AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
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),
84 fVertex(mcEvnt.fVertex)
90 AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
92 // assignment operator
94 AliVEvent::operator=(mcEvnt);
100 void AliMCEvent::ConnectTreeE (TTree* tree)
102 // Connect the event header tree
103 tree->SetBranchAddress("Header", &fHeader);
106 void AliMCEvent::ConnectTreeK (TTree* tree)
108 // Connect the kinematics tree to the stack
109 if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
111 fStack = fHeader->Stack();
112 fStack->ConnectTree(tree);
116 fNparticles = fStack->GetNtrack();
117 fNprimaries = fStack->GetNprimary();
119 Int_t iev = fHeader->GetEvent();
120 Int_t ievr = fHeader->GetEventNrInRun();
121 AliDebug(1, Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
122 iev, ievr, fNparticles, fNprimaries));
124 // This is a cache for the TParticles converted to MCParticles on user request
125 if (fMCParticleMap) {
126 fMCParticleMap->Clear();
127 fMCParticles->Delete();
128 if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
131 fMCParticleMap = new TRefArray(fNparticles);
134 void AliMCEvent::ConnectTreeTR (TTree* tree)
136 // Connect the track reference tree
139 if (fTreeTR->GetBranch("AliRun")) {
144 // This is an old format with one branch per detector not in synch with TreeK
145 ReorderAndExpandTreeTR();
148 fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
152 Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
155 if (i < 0 || i >= fNparticles) {
156 AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
161 particle = fStack->Particle(i);
163 fTreeTR->GetEntry(fStack->TreeKEntry(i));
165 return trefs->GetEntries();
173 void AliMCEvent::Clean()
175 // Clean-up before new trees are connected
176 delete fStack; fStack = 0;
188 void AliMCEvent::FinishEvent()
190 // Clean-up after event
192 if (fStack) fStack->Reset(0);
193 fMCParticles->Delete();
196 fMCParticleMap->Clear();
200 // fTrackReferences->Delete();
201 fTrackReferences->Clear();
205 // fSubsidiaryEvents->Clear();
206 fSubsidiaryEvents = 0;
211 void AliMCEvent::DrawCheck(Int_t i, Int_t search)
214 // Simple event display for debugging
216 AliWarning("No Track Reference information available");
220 if (i > -1 && i < fNparticles) {
221 fTreeTR->GetEntry(fStack->TreeKEntry(i));
223 AliWarning("AliMCEvent::GetEntry: Index out of range");
226 Int_t nh = fTRBuffer->GetEntries();
230 while(nh <= search && i < fNparticles - 1) {
232 fTreeTR->GetEntry(fStack->TreeKEntry(i));
233 nh = fTRBuffer->GetEntries();
235 printf("Found Hits at %5d\n", i);
237 TParticle* particle = fStack->Particle(i);
239 TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
240 Float_t x0 = particle->Vx();
241 Float_t y0 = particle->Vy();
243 Float_t x1 = particle->Vx() + particle->Px() * 50.;
244 Float_t y1 = particle->Vy() + particle->Py() * 50.;
246 TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
252 for (Int_t ih = 0; ih < nh; ih++) {
253 AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
254 TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
256 m->SetMarkerSize(0.4);
262 void AliMCEvent::ReorderAndExpandTreeTR()
265 // Reorder and expand the track reference tree in order to match the kinematics tree.
266 // Copy the information from different branches into one
270 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
271 fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
272 if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
273 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 64000, 0);
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);
287 // Connect the active branches
288 TClonesArray* trefs[7];
289 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
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]);
301 Int_t np = fStack->GetNprimary();
302 Int_t nt = fTreeTR->GetEntries();
305 // Loop over tracks and find the secondaries with the help of the kine tree
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));
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
322 Int_t inext = ip - 1;
325 part = fStack->Particle(inext);
326 dau2 = part->GetFirstDaughter();
327 if (dau2 == -1 || dau2 < np) {
333 dau2 = fStack->GetNtrack() - 1;
336 } // find upper bound
340 // printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
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.
346 Bool_t hasHits = kFALSE;
347 Bool_t isOutside = kFALSE;
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) {
363 if (label > dau2 || label < ip) {
369 if (hasHits || isOutside) break;
374 // Write empty entries
375 for (Int_t id = dau1; (id <= dau2); id++) {
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();
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);
396 tr->SetDetectorId(ib-1);
397 Int_t nref = fTRBuffer->GetEntriesFast();
398 TClonesArray &lref = *fTRBuffer;
399 new(lref[nref]) AliTrackReference(*tr);
411 // Now loop again and write the primaries
414 for (Int_t ip = 0; ip < np; ip++) {
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();
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) {
431 tr->SetDetectorId(ib-1);
432 Int_t nref = fTRBuffer->GetEntriesFast();
433 TClonesArray &lref = *fTRBuffer;
434 new(lref[nref]) AliTrackReference(*tr);
448 delete fTreeTR; fTreeTR = 0;
450 for (Int_t ib = 0; ib < 7; ib++) {
458 if (ifills != fStack->GetNtrack())
459 printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
460 ifills, fStack->GetNtrack());
463 fTreeTR = fTmpTreeTR;
466 AliVParticle* AliMCEvent::GetTrack(Int_t i) const
472 return ((AliVParticle*) (fMCParticles->At(i)));
476 // Check first if this explicitely accesses the subsidiary event
478 if (i >= BgLabelOffset()) {
479 if (fSubsidiaryEvents) {
480 AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
481 return (bgEvent->GetTrack(i - BgLabelOffset()));
488 AliMCParticle *mcParticle = 0;
489 TParticle *particle = 0;
490 TClonesArray *trefs = 0;
492 TRefArray *rarray = 0;
496 // Out of range check
497 if (i < 0 || i >= fNparticles) {
498 AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
504 if (fSubsidiaryEvents) {
506 Int_t idx = FindIndexAndEvent(i, mc);
507 return (mc->GetTrack(idx));
511 // First check If the MC Particle has been already cached
512 if(!fMCParticleMap->At(i)) {
513 // Get particle from the stack
514 particle = fStack->Particle(i);
515 // Get track references from Tree TR
517 fTreeTR->GetEntry(fStack->TreeKEntry(i));
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);
529 } // loop over track references for entry i
530 } // if TreeTR available
531 Int_t nentries = fMCParticles->GetEntriesFast();
532 new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
533 mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
534 fMCParticleMap->AddAt(mcParticle, i);
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);
545 mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
548 if (id1 < fNprimaries) {
549 mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
550 mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
552 mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
553 mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
557 if (i > fNprimaries) {
558 mcParticle->SetLabel(i + fPrimaryOffset);
560 mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
565 mcParticle->SetFirstDaughter(id1);
566 mcParticle->SetLastDaughter (id2);
567 mcParticle->SetMother (imo);
571 mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
578 AliGenEventHeader* AliMCEvent::GenEventHeader() const {return (fHeader->GenEventHeader());}
581 void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event)
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;
590 fSubsidiaryEvents->Add(event);
593 Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
595 // Find the index and event in case of composed events like signal + background
596 TIter next(fSubsidiaryEvents);
598 if (oldidx < fNprimaries) {
599 while((event = (AliMCEvent*)next())) {
600 if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
602 return (oldidx - event->GetPrimaryOffset());
604 while((event = (AliMCEvent*)next())) {
605 if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
607 return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
611 Int_t AliMCEvent::BgLabelToIndex(Int_t label)
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();
620 label += (bgEvent->GetSecondaryOffset() - fNprimaries);
627 Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i)
630 // Delegate to subevent if necesarry
633 if (!fSubsidiaryEvents) {
634 return fStack->IsPhysicalPrimary(i);
637 Int_t idx = FindIndexAndEvent(i, evt);
638 return (evt->IsPhysicalPrimary(idx));
642 void AliMCEvent::InitEvent()
645 // Initialize the subsidiary event structure
646 if (fSubsidiaryEvents) {
647 TIter next(fSubsidiaryEvents);
652 while((evt = (AliMCEvent*)next())) {
653 fNprimaries += evt->GetNumberOfPrimaries();
654 fNparticles += evt->GetNumberOfTracks();
658 Int_t ioffs = fNprimaries;
661 while((evt = (AliMCEvent*)next())) {
662 evt->SetPrimaryOffset(ioffp);
663 evt->SetSecondaryOffset(ioffs);
664 ioffp += evt->GetNumberOfPrimaries();
665 ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
669 void AliMCEvent::PreReadAll()
671 // Preread the MC information
674 for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
679 for (i = 0; i < fStack->GetNprimary(); i++)
688 const AliVVertex * AliMCEvent::GetPrimaryVertex() const
690 // Create a MCVertex object from the MCHeader information
692 GenEventHeader()->PrimaryVertex(v) ;
694 fVertex = new AliMCVertex(v[0], v[1], v[2]);
696 ((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);