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>
33 #include "AliMCEvent.h"
34 #include "AliMCVertex.h"
36 #include "AliTrackReference.h"
37 #include "AliHeader.h"
38 #include "AliGenEventHeader.h"
39 #include "AliGenHijingEventHeader.h"
40 #include "AliGenCocktailEventHeader.h"
43 Int_t AliMCEvent::fgkBgLabelOffset(10000000);
46 AliMCEvent::AliMCEvent():
51 fHeader(new AliHeader()),
53 fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
66 // Default constructor
69 AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
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),
86 fVertex(mcEvnt.fVertex),
93 AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
95 // assignment operator
97 AliVEvent::operator=(mcEvnt);
103 void AliMCEvent::ConnectTreeE (TTree* tree)
105 // Connect the event header tree
106 tree->SetBranchAddress("Header", &fHeader);
109 void AliMCEvent::ConnectTreeK (TTree* tree)
111 // Connect the kinematics tree to the stack
112 if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
114 fStack = fHeader->Stack();
115 fStack->ConnectTree(tree);
119 fNparticles = fStack->GetNtrack();
120 fNprimaries = fStack->GetNprimary();
122 Int_t iev = fHeader->GetEvent();
123 Int_t ievr = fHeader->GetEventNrInRun();
124 AliDebug(1, Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
125 iev, ievr, fNparticles, fNprimaries));
127 // This is a cache for the TParticles converted to MCParticles on user request
128 if (fMCParticleMap) {
129 fMCParticleMap->Clear();
130 fMCParticles->Delete();
131 if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
134 fMCParticleMap = new TObjArray(fNparticles);
137 void AliMCEvent::ConnectTreeTR (TTree* tree)
139 // Connect the track reference tree
142 if (fTreeTR->GetBranch("AliRun")) {
147 // This is an old format with one branch per detector not in synch with TreeK
148 ReorderAndExpandTreeTR();
151 fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
155 Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
158 if (i < 0 || i >= fNparticles) {
159 AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
164 particle = fStack->Particle(i);
166 fTreeTR->GetEntry(fStack->TreeKEntry(i));
168 return trefs->GetEntries();
176 void AliMCEvent::Clean()
178 // Clean-up before new trees are connected
179 delete fStack; fStack = 0;
191 void AliMCEvent::FinishEvent()
193 // Clean-up after event
195 if (fStack) fStack->Reset(0);
196 fMCParticles->Delete();
199 fMCParticleMap->Clear();
203 // fTrackReferences->Delete();
204 fTrackReferences->Clear();
208 // fSubsidiaryEvents->Clear();
209 fSubsidiaryEvents = 0;
215 void AliMCEvent::DrawCheck(Int_t i, Int_t search)
218 // Simple event display for debugging
220 AliWarning("No Track Reference information available");
224 if (i > -1 && i < fNparticles) {
225 fTreeTR->GetEntry(fStack->TreeKEntry(i));
227 AliWarning("AliMCEvent::GetEntry: Index out of range");
230 Int_t nh = fTRBuffer->GetEntries();
234 while(nh <= search && i < fNparticles - 1) {
236 fTreeTR->GetEntry(fStack->TreeKEntry(i));
237 nh = fTRBuffer->GetEntries();
239 printf("Found Hits at %5d\n", i);
241 TParticle* particle = fStack->Particle(i);
243 TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
244 Float_t x0 = particle->Vx();
245 Float_t y0 = particle->Vy();
247 Float_t x1 = particle->Vx() + particle->Px() * 50.;
248 Float_t y1 = particle->Vy() + particle->Py() * 50.;
250 TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
256 for (Int_t ih = 0; ih < nh; ih++) {
257 AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
258 TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
260 m->SetMarkerSize(0.4);
266 void AliMCEvent::ReorderAndExpandTreeTR()
269 // Reorder and expand the track reference tree in order to match the kinematics tree.
270 // Copy the information from different branches into one
274 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
275 fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
276 if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
277 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 64000, 0);
281 // Activate the used branches only. Otherwisw we get a bad memory leak.
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);
294 // Connect the active branches
295 TClonesArray* trefs[7];
296 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
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]);
308 Int_t np = fStack->GetNprimary();
309 Int_t nt = fTreeTR->GetEntries();
312 // Loop over tracks and find the secondaries with the help of the kine tree
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));
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
329 Int_t inext = ip - 1;
332 part = fStack->Particle(inext);
333 dau2 = part->GetFirstDaughter();
334 if (dau2 == -1 || dau2 < np) {
340 dau2 = fStack->GetNtrack() - 1;
343 } // find upper bound
347 // printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
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.
353 Bool_t hasHits = kFALSE;
354 Bool_t isOutside = kFALSE;
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) {
370 if (label > dau2 || label < ip) {
376 if (hasHits || isOutside) break;
381 // Write empty entries
382 for (Int_t id = dau1; (id <= dau2); id++) {
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();
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);
403 tr->SetDetectorId(ib-1);
404 Int_t nref = fTRBuffer->GetEntriesFast();
405 TClonesArray &lref = *fTRBuffer;
406 new(lref[nref]) AliTrackReference(*tr);
418 // Now loop again and write the primaries
421 for (Int_t ip = 0; ip < np; ip++) {
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();
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) {
438 tr->SetDetectorId(ib-1);
439 Int_t nref = fTRBuffer->GetEntriesFast();
440 TClonesArray &lref = *fTRBuffer;
441 new(lref[nref]) AliTrackReference(*tr);
455 delete fTreeTR; fTreeTR = 0;
457 for (Int_t ib = 0; ib < 7; ib++) {
465 if (ifills != fStack->GetNtrack())
466 printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
467 ifills, fStack->GetNtrack());
470 fTreeTR = fTmpTreeTR;
473 AliVParticle* AliMCEvent::GetTrack(Int_t i) const
479 return ((AliVParticle*) (fMCParticles->At(i)));
483 // Check first if this explicitely accesses the subsidiary event
485 if (i >= BgLabelOffset()) {
486 if (fSubsidiaryEvents) {
487 AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
488 return (bgEvent->GetTrack(i - BgLabelOffset()));
495 AliMCParticle *mcParticle = 0;
496 TParticle *particle = 0;
497 TClonesArray *trefs = 0;
499 TObjArray *rarray = 0;
503 // Out of range check
504 if (i < 0 || i >= fNparticles) {
505 AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
511 if (fSubsidiaryEvents) {
513 Int_t idx = FindIndexAndEvent(i, mc);
514 return (mc->GetTrack(idx));
518 // First check If the MC Particle has been already cached
519 if(!fMCParticleMap->At(i)) {
520 // Get particle from the stack
521 particle = fStack->Particle(i);
522 // Get track references from Tree TR
524 fTreeTR->GetEntry(fStack->TreeKEntry(i));
526 ntref = trefs->GetEntriesFast();
527 rarray = new TObjArray(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
534 new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
535 rarray->AddAt((*fTrackReferences)[nen], j);
538 } // loop over track references for entry i
539 } // if TreeTR available
540 Int_t nentries = fMCParticles->GetEntriesFast();
541 mcParticle = new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
542 fMCParticleMap->AddAt(mcParticle, i);
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);
553 mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
556 if (id1 < fNprimaries) {
557 mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
558 mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
560 mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
561 mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
565 if (i > fNprimaries) {
566 mcParticle->SetLabel(i + fPrimaryOffset);
568 mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
571 mcParticle->SetFirstDaughter(id1);
572 mcParticle->SetLastDaughter (id2);
573 mcParticle->SetMother (imo);
577 mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
582 AliGenEventHeader* AliMCEvent::GenEventHeader() const {return (fHeader->GenEventHeader());}
585 void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event)
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;
594 fSubsidiaryEvents->Add(event);
597 Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
599 // Find the index and event in case of composed events like signal + background
600 TIter next(fSubsidiaryEvents);
602 if (oldidx < fNprimaries) {
603 while((event = (AliMCEvent*)next())) {
604 if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
607 return (oldidx - event->GetPrimaryOffset());
612 while((event = (AliMCEvent*)next())) {
613 if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
616 return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
623 Int_t AliMCEvent::BgLabelToIndex(Int_t label)
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();
632 label += (bgEvent->GetSecondaryOffset() - fNprimaries);
639 Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) const
642 // Delegate to subevent if necesarry
645 if (!fSubsidiaryEvents) {
646 return fStack->IsPhysicalPrimary(i);
649 Int_t idx = FindIndexAndEvent(i, evt);
650 return (evt->IsPhysicalPrimary(idx));
654 void AliMCEvent::InitEvent()
657 // Initialize the subsidiary event structure
658 if (fSubsidiaryEvents) {
659 TIter next(fSubsidiaryEvents);
664 while((evt = (AliMCEvent*)next())) {
665 fNprimaries += evt->GetNumberOfPrimaries();
666 fNparticles += evt->GetNumberOfTracks();
670 Int_t ioffs = fNprimaries;
673 while((evt = (AliMCEvent*)next())) {
674 evt->SetPrimaryOffset(ioffp);
675 evt->SetSecondaryOffset(ioffs);
676 ioffp += evt->GetNumberOfPrimaries();
677 ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
682 void AliMCEvent::PreReadAll()
684 // Preread the MC information
687 for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
692 for (i = 0; i < fStack->GetNprimary(); i++)
698 const AliVVertex * AliMCEvent::GetPrimaryVertex() const
700 // Create a MCVertex object from the MCHeader information
702 GenEventHeader()->PrimaryVertex(v) ;
704 fVertex = new AliMCVertex(v[0], v[1], v[2]);
706 ((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);
711 Bool_t AliMCEvent::IsFromBGEvent(Int_t index)
713 // Checks if a particle is from the background events
714 // Works for HIJING inside Cocktail
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();
725 return (index < fNBG);