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>
32 #include "AliMCEvent.h"
34 #include "AliTrackReference.h"
35 #include "AliHeader.h"
36 #include "AliGenEventHeader.h"
39 AliMCEvent::AliMCEvent():
42 fMCParticles(new TClonesArray("AliMCParticle",1000)),
44 fHeader(new AliHeader()),
46 fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
53 // Default constructor
56 AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
61 fHeader(new AliHeader()),
74 AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
76 // assignment operator
78 AliVEvent::operator=(mcEvnt);
84 void AliMCEvent::ConnectTreeE (TTree* tree)
86 // Connect the event header tree
87 tree->SetBranchAddress("Header", &fHeader);
90 void AliMCEvent::ConnectTreeK (TTree* tree)
92 // Connect the kinematics tree to the stack
93 fStack = fHeader->Stack();
94 fStack->ConnectTree(tree);
98 fNparticles = fStack->GetNtrack();
99 fNprimaries = fStack->GetNprimary();
100 Int_t iev = fHeader->GetEvent();
101 Int_t ievr = fHeader->GetEventNrInRun();
103 AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
104 iev, ievr, fNparticles, fNprimaries));
106 // This is a cache for the TParticles converted to MCParticles on user request
107 if (fMCParticleMap) {
108 fMCParticleMap->Clear();
109 fMCParticles->Delete();
110 if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
113 fMCParticleMap = new TRefArray(fNparticles);
116 void AliMCEvent::ConnectTreeTR (TTree* tree)
118 // Connect the track reference tree
121 if (fTreeTR->GetBranch("AliRun")) {
126 // This is an old format with one branch per detector not in synch with TreeK
127 ReorderAndExpandTreeTR();
130 fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
134 Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
137 if (i < 0 || i >= fNparticles) {
138 AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
143 particle = fStack->Particle(i);
145 fTreeTR->GetEntry(fStack->TreeKEntry(i));
147 return trefs->GetEntries();
155 void AliMCEvent::Clean()
157 // Clean-up before new trees are connected
163 delete fStack; fStack = 0;
173 void AliMCEvent::FinishEvent()
175 // Clean-up after event
178 fMCParticles->Delete();
179 fMCParticleMap->Clear();
182 fTrackReferences->Delete();
191 void AliMCEvent::DrawCheck(Int_t i, Int_t search)
194 // Simple event display for debugging
196 AliWarning("No Track Reference information available");
200 if (i > -1 && i < fNparticles) {
201 fTreeTR->GetEntry(fStack->TreeKEntry(i));
203 AliWarning("AliMCEvent::GetEntry: Index out of range");
206 Int_t nh = fTRBuffer->GetEntries();
210 while(nh <= search && i < fNparticles - 1) {
212 fTreeTR->GetEntry(fStack->TreeKEntry(i));
213 nh = fTRBuffer->GetEntries();
215 printf("Found Hits at %5d\n", i);
217 TParticle* particle = fStack->Particle(i);
219 TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
220 Float_t x0 = particle->Vx();
221 Float_t y0 = particle->Vy();
223 Float_t x1 = particle->Vx() + particle->Px() * 50.;
224 Float_t y1 = particle->Vy() + particle->Py() * 50.;
226 TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
232 for (Int_t ih = 0; ih < nh; ih++) {
233 AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
234 TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
236 m->SetMarkerSize(0.4);
242 void AliMCEvent::ReorderAndExpandTreeTR()
245 // Reorder and expand the track reference tree in order to match the kinematics tree.
246 // Copy the information from different branches into one
250 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
251 fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
252 if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
253 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
257 // Activate the used branches only. Otherwisw we get a bad memory leak.
258 fTreeTR->SetBranchStatus("*", 0);
259 fTreeTR->SetBranchStatus("AliRun.*", 1);
260 fTreeTR->SetBranchStatus("ITS.*", 1);
261 fTreeTR->SetBranchStatus("TPC.*", 1);
262 fTreeTR->SetBranchStatus("TRD.*", 1);
263 fTreeTR->SetBranchStatus("TOF.*", 1);
264 fTreeTR->SetBranchStatus("FRAME.*", 1);
265 fTreeTR->SetBranchStatus("MUON.*", 1);
267 // Connect the active branches
268 TClonesArray* trefs[7];
269 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
271 // make branch for central track references
272 if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
273 if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
274 if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
275 if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
276 if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
277 if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
278 if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
281 Int_t np = fStack->GetNprimary();
282 Int_t nt = fTreeTR->GetEntries();
285 // Loop over tracks and find the secondaries with the help of the kine tree
291 for (Int_t ip = np - 1; ip > -1; ip--) {
292 part = fStack->Particle(ip);
293 // printf("Particle %5d %5d %5d %5d %5d %5d \n",
294 // ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
295 // part->GetLastDaughter(), part->TestBit(kTransportBit));
297 // Determine range of secondaries produced by this primary during transport
298 Int_t dau1 = part->GetFirstDaughter();
299 if (dau1 < np) continue; // This particle has no secondaries produced during transport
302 Int_t inext = ip - 1;
305 part = fStack->Particle(inext);
306 dau2 = part->GetFirstDaughter();
307 if (dau2 == -1 || dau2 < np) {
313 dau2 = fStack->GetNtrack() - 1;
316 } // find upper bound
320 // printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
322 // Loop over reference hits and find secondary label
323 // First the tricky part: find the entry in treeTR than contains the hits or
324 // make sure that no hits exist.
326 Bool_t hasHits = kFALSE;
327 Bool_t isOutside = kFALSE;
330 while (!hasHits && !isOutside && it < nt) {
331 fTreeTR->GetEntry(it++);
332 for (Int_t ib = 0; ib < 7; ib++) {
333 if (!trefs[ib]) continue;
334 Int_t nh = trefs[ib]->GetEntries();
335 for (Int_t ih = 0; ih < nh; ih++) {
336 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
337 Int_t label = tr->Label();
338 if (label >= dau1 && label <= dau2) {
343 if (label > dau2 || label < ip) {
349 if (hasHits || isOutside) break;
354 // Write empty entries
355 for (Int_t id = dau1; (id <= dau2); id++) {
361 fTreeTR->GetEntry(itlast);
362 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
363 for (Int_t ib = 0; ib < 7; ib++) {
364 if (!trefs[ib]) continue;
365 Int_t nh = trefs[ib]->GetEntries();
366 for (Int_t ih = 0; ih < nh; ih++) {
367 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
368 Int_t label = tr->Label();
370 if (label == ip) continue;
371 if (label > dau2 || label < dau1)
372 printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
373 itlast, label, dau1, dau2);
376 tr->SetDetectorId(ib-1);
377 Int_t nref = fTRBuffer->GetEntriesFast();
378 TClonesArray &lref = *fTRBuffer;
379 new(lref[nref]) AliTrackReference(*tr);
391 // Now loop again and write the primaries
394 for (Int_t ip = 0; ip < np; ip++) {
396 while (labmax < ip && it > -1) {
397 fTreeTR->GetEntry(it--);
398 for (Int_t ib = 0; ib < 7; ib++) {
399 if (!trefs[ib]) continue;
400 Int_t nh = trefs[ib]->GetEntries();
402 // Loop over reference hits and find primary labels
403 for (Int_t ih = 0; ih < nh; ih++) {
404 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
405 Int_t label = tr->Label();
406 if (label < np && label > labmax) {
411 tr->SetDetectorId(ib-1);
412 Int_t nref = fTRBuffer->GetEntriesFast();
413 TClonesArray &lref = *fTRBuffer;
414 new(lref[nref]) AliTrackReference(*tr);
428 delete fTreeTR; fTreeTR = 0;
430 for (Int_t ib = 0; ib < 7; ib++) {
438 if (ifills != fStack->GetNtrack())
439 printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
440 ifills, fStack->GetNtrack());
443 fTreeTR = fTmpTreeTR;
446 AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
449 AliMCParticle *mcParticle = 0;
450 TParticle *particle = 0;
451 TClonesArray *trefs = 0;
453 TRefArray *rarray = 0;
454 // Out of range check
455 if (i < 0 || i >= fNparticles) {
456 AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
463 // First check of the MC Particle has been already cached
464 if(!fMCParticleMap->At(i)) {
465 // Get particle from the stack
466 particle = fStack->Particle(i);
467 // Get track references from Tree TR
469 fTreeTR->GetEntry(fStack->TreeKEntry(i));
471 ntref = trefs->GetEntriesFast();
472 rarray = new TRefArray(ntref);
473 Int_t nen = fTrackReferences->GetEntriesFast();
474 for (Int_t j = 0; j < ntref; j++) {
475 // Save the track references in a TClonesArray
476 AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
477 // Save the pointer in a TRefArray
478 new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
479 rarray->AddAt((*fTrackReferences)[nen], j);
481 } // loop over track references for entry i
482 } // if TreeTR available
483 Int_t nentries = fMCParticles->GetEntriesFast();
484 new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
485 mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
486 fMCParticleMap->AddAt(mcParticle, i);
488 mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
493 AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}