Reset object count after each event.
[u/mrichter/AliRoot.git] / STEER / 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>
93df0e9b 29#include <TRefArray.h>
415d9f5c 30
31#include "AliLog.h"
32#include "AliMCEvent.h"
33#include "AliStack.h"
34#include "AliTrackReference.h"
35#include "AliHeader.h"
83626b6b 36#include "AliGenEventHeader.h"
415d9f5c 37
38
39AliMCEvent::AliMCEvent():
40 AliVEvent(),
41 fStack(0),
42 fMCParticles(new TClonesArray("AliMCParticle",1000)),
43 fMCParticleMap(0),
44 fHeader(0),
93df0e9b 45 fTRBuffer(0),
46 fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
415d9f5c 47 fTreeTR(0),
48 fTmpTreeTR(0),
49 fTmpFileTR(0),
50 fNprimaries(-1),
51 fNparticles(-1)
52{
53 // Default constructor
54}
55
56AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
57 AliVEvent(mcEvnt),
58 fStack(0),
59 fMCParticles(0),
60 fMCParticleMap(0),
61 fHeader(0),
93df0e9b 62 fTRBuffer(0),
415d9f5c 63 fTrackReferences(0),
64 fTreeTR(0),
65 fTmpTreeTR(0),
66 fTmpFileTR(0),
67 fNprimaries(-1),
68 fNparticles(-1)
69{
70// Copy constructor
71}
72
73
74AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
75{
76 // assignment operator
77 if (this!=&mcEvnt) {
78 AliVEvent::operator=(mcEvnt);
79 }
80
81 return *this;
82}
83
84void AliMCEvent::ConnectTreeE (TTree* tree)
85{
86 // Connect the event header tree
87 tree->SetBranchAddress("Header", &fHeader);
88}
89
90void AliMCEvent::ConnectTreeK (TTree* tree)
91{
92 // Connect the kinematics tree to the stack
93 fStack = fHeader->Stack();
94 fStack->ConnectTree(tree);
95 //
96 // Load the event
97 fStack->GetEvent();
98 fNparticles = fStack->GetNtrack();
99 fNprimaries = fStack->GetNprimary();
ed97dc98 100 Int_t iev = fHeader->GetEvent();
101 Int_t ievr = fHeader->GetEventNrInRun();
102
103 AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
104 iev, ievr, fNparticles, fNprimaries));
415d9f5c 105
106 // This is a cache for the TParticles converted to MCParticles on user request
107 if (fMCParticleMap) {
ee2774b7 108 fMCParticleMap->Delete();
109 fMCParticles->Delete();
415d9f5c 110 if (fNparticles>0) fMCParticleMap->Expand(fNparticles);}
111 else
93df0e9b 112 fMCParticleMap = new TRefArray(fNparticles);
113
415d9f5c 114}
115
116void AliMCEvent::ConnectTreeTR (TTree* tree)
117{
118 // Connect the track reference tree
119 fTreeTR = tree;
120
121 if (fTreeTR->GetBranch("AliRun")) {
122 if (fTmpFileTR) {
123 fTmpFileTR->Close();
124 delete fTmpFileTR;
125 }
126 // This is an old format with one branch per detector not in synch with TreeK
127 ReorderAndExpandTreeTR();
128 } else {
129 // New format
93df0e9b 130 fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
415d9f5c 131 }
132}
133
134Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
135{
136 // Retrieve entry i
137 if (i < 0 || i >= fNparticles) {
138 AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
139 particle = 0;
140 trefs = 0;
141 return (-1);
142 }
143 particle = fStack->Particle(i);
144 if (fTreeTR) {
145 fTreeTR->GetEntry(fStack->TreeKEntry(i));
93df0e9b 146 trefs = fTRBuffer;
415d9f5c 147 return trefs->GetEntries();
148 } else {
149 trefs = 0;
150 return -1;
151 }
152}
153
154
155void AliMCEvent::Clean()
156{
157 // Clean-up before new trees are connected
415d9f5c 158 if (fHeader) {
159 delete fHeader;
160 fHeader = 0;
161 }
c6360f33 162
163 delete fStack; fStack = 0;
415d9f5c 164
165 // Clear TR
93df0e9b 166 if (fTRBuffer) {
ee2774b7 167 fTRBuffer->Delete();
93df0e9b 168 delete fTRBuffer;
169 fTRBuffer = 0;
415d9f5c 170 }
171}
172
173void AliMCEvent::FinishEvent()
174{
f1ac8a97 175 // Clean-up after event
c6360f33 176 //
177 fStack->Reset(0);
178 fMCParticles->Delete();
179 fMCParticleMap->Delete();
ee2774b7 180 fTRBuffer->Delete();
c6360f33 181 fTrackReferences->Delete();
415d9f5c 182}
183
184
93df0e9b 185
415d9f5c 186void AliMCEvent::DrawCheck(Int_t i, Int_t search)
187{
188 //
189 // Simple event display for debugging
190 if (!fTreeTR) {
191 AliWarning("No Track Reference information available");
192 return;
193 }
194
195 if (i > -1 && i < fNparticles) {
196 fTreeTR->GetEntry(fStack->TreeKEntry(i));
197 } else {
198 AliWarning("AliMCEvent::GetEntry: Index out of range");
199 }
200
93df0e9b 201 Int_t nh = fTRBuffer->GetEntries();
415d9f5c 202
203
204 if (search) {
205 while(nh <= search && i < fNparticles - 1) {
206 i++;
207 fTreeTR->GetEntry(fStack->TreeKEntry(i));
93df0e9b 208 nh = fTRBuffer->GetEntries();
415d9f5c 209 }
210 printf("Found Hits at %5d\n", i);
211 }
212 TParticle* particle = fStack->Particle(i);
213
214 TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
215 Float_t x0 = particle->Vx();
216 Float_t y0 = particle->Vy();
217
218 Float_t x1 = particle->Vx() + particle->Px() * 50.;
219 Float_t y1 = particle->Vy() + particle->Py() * 50.;
220
221 TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
222 h->Draw();
223 a->SetLineColor(2);
224
225 a->Draw();
226
227 for (Int_t ih = 0; ih < nh; ih++) {
93df0e9b 228 AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
415d9f5c 229 TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
230 m->Draw();
231 m->SetMarkerSize(0.4);
232
233 }
234}
235
236
237void AliMCEvent::ReorderAndExpandTreeTR()
238{
239//
240// Reorder and expand the track reference tree in order to match the kinematics tree.
241// Copy the information from different branches into one
242//
243// TreeTR
244
245 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
246 fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
93df0e9b 247 if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
248 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
415d9f5c 249
250
251//
252// Activate the used branches only. Otherwisw we get a bad memory leak.
253 fTreeTR->SetBranchStatus("*", 0);
254 fTreeTR->SetBranchStatus("AliRun.*", 1);
255 fTreeTR->SetBranchStatus("ITS.*", 1);
256 fTreeTR->SetBranchStatus("TPC.*", 1);
257 fTreeTR->SetBranchStatus("TRD.*", 1);
258 fTreeTR->SetBranchStatus("TOF.*", 1);
259 fTreeTR->SetBranchStatus("FRAME.*", 1);
260 fTreeTR->SetBranchStatus("MUON.*", 1);
261//
262// Connect the active branches
263 TClonesArray* trefs[7];
264 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
265 if (fTreeTR){
266 // make branch for central track references
267 if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
268 if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
269 if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
270 if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
271 if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
272 if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
273 if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
274 }
275
276 Int_t np = fStack->GetNprimary();
277 Int_t nt = fTreeTR->GetEntries();
278
279 //
280 // Loop over tracks and find the secondaries with the help of the kine tree
281 Int_t ifills = 0;
282 Int_t it = 0;
283 Int_t itlast = 0;
284 TParticle* part;
285
286 for (Int_t ip = np - 1; ip > -1; ip--) {
287 part = fStack->Particle(ip);
288// printf("Particle %5d %5d %5d %5d %5d %5d \n",
289// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
290// part->GetLastDaughter(), part->TestBit(kTransportBit));
291
292 // Determine range of secondaries produced by this primary during transport
293 Int_t dau1 = part->GetFirstDaughter();
294 if (dau1 < np) continue; // This particle has no secondaries produced during transport
295 Int_t dau2 = -1;
296 if (dau1 > -1) {
297 Int_t inext = ip - 1;
298 while (dau2 < 0) {
299 if (inext >= 0) {
300 part = fStack->Particle(inext);
301 dau2 = part->GetFirstDaughter();
302 if (dau2 == -1 || dau2 < np) {
303 dau2 = -1;
304 } else {
305 dau2--;
306 }
307 } else {
308 dau2 = fStack->GetNtrack() - 1;
309 }
310 inext--;
311 } // find upper bound
312 } // dau2 < 0
313
314
315// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
316//
317// Loop over reference hits and find secondary label
318// First the tricky part: find the entry in treeTR than contains the hits or
319// make sure that no hits exist.
320//
321 Bool_t hasHits = kFALSE;
322 Bool_t isOutside = kFALSE;
323
324 it = itlast;
325 while (!hasHits && !isOutside && it < nt) {
326 fTreeTR->GetEntry(it++);
327 for (Int_t ib = 0; ib < 7; ib++) {
328 if (!trefs[ib]) continue;
329 Int_t nh = trefs[ib]->GetEntries();
330 for (Int_t ih = 0; ih < nh; ih++) {
331 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
332 Int_t label = tr->Label();
333 if (label >= dau1 && label <= dau2) {
334 hasHits = kTRUE;
335 itlast = it - 1;
336 break;
337 }
338 if (label > dau2 || label < ip) {
339 isOutside = kTRUE;
340 itlast = it - 1;
341 break;
342 }
343 } // hits
344 if (hasHits || isOutside) break;
345 } // branches
346 } // entries
347
348 if (!hasHits) {
349 // Write empty entries
350 for (Int_t id = dau1; (id <= dau2); id++) {
351 fTmpTreeTR->Fill();
352 ifills++;
353 }
354 } else {
355 // Collect all hits
356 fTreeTR->GetEntry(itlast);
357 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
358 for (Int_t ib = 0; ib < 7; ib++) {
359 if (!trefs[ib]) continue;
360 Int_t nh = trefs[ib]->GetEntries();
361 for (Int_t ih = 0; ih < nh; ih++) {
362 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
363 Int_t label = tr->Label();
364 // Skip primaries
365 if (label == ip) continue;
366 if (label > dau2 || label < dau1)
367 printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
368 itlast, label, dau1, dau2);
369 if (label == id) {
370 // secondary found
371 tr->SetDetectorId(ib-1);
93df0e9b 372 Int_t nref = fTRBuffer->GetEntriesFast();
373 TClonesArray &lref = *fTRBuffer;
415d9f5c 374 new(lref[nref]) AliTrackReference(*tr);
375 }
376 } // hits
377 } // branches
378 fTmpTreeTR->Fill();
ee2774b7 379 fTRBuffer->Delete();
415d9f5c 380 ifills++;
381 } // daughters
382 } // has hits
383 } // tracks
384
385 //
386 // Now loop again and write the primaries
387 //
388 it = nt - 1;
389 for (Int_t ip = 0; ip < np; ip++) {
390 Int_t labmax = -1;
391 while (labmax < ip && it > -1) {
392 fTreeTR->GetEntry(it--);
393 for (Int_t ib = 0; ib < 7; ib++) {
394 if (!trefs[ib]) continue;
395 Int_t nh = trefs[ib]->GetEntries();
396 //
397 // Loop over reference hits and find primary labels
398 for (Int_t ih = 0; ih < nh; ih++) {
399 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
400 Int_t label = tr->Label();
401 if (label < np && label > labmax) {
402 labmax = label;
403 }
404
405 if (label == ip) {
406 tr->SetDetectorId(ib-1);
93df0e9b 407 Int_t nref = fTRBuffer->GetEntriesFast();
408 TClonesArray &lref = *fTRBuffer;
415d9f5c 409 new(lref[nref]) AliTrackReference(*tr);
410 }
411 } // hits
412 } // branches
413 } // entries
414 it++;
415 fTmpTreeTR->Fill();
ee2774b7 416 fTRBuffer->Delete();
415d9f5c 417 ifills++;
418 } // tracks
419 // Check
420
421
422 // Clean-up
423 delete fTreeTR; fTreeTR = 0;
424
425 for (Int_t ib = 0; ib < 7; ib++) {
426 if (trefs[ib]) {
427 trefs[ib]->Clear();
428 delete trefs[ib];
429 trefs[ib] = 0;
430 }
431 }
432
433 if (ifills != fStack->GetNtrack())
434 printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
435 ifills, fStack->GetNtrack());
436
437 fTmpTreeTR->Write();
438 fTreeTR = fTmpTreeTR;
439}
440
441AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
442{
443 // Get MC Particle i
93df0e9b 444 AliMCParticle *mcParticle = 0;
445 TParticle *particle = 0;
446 TClonesArray *trefs = 0;
447 Int_t ntref = 0;
448 TRefArray *rarray = 0;
415d9f5c 449 // Out of range check
450 if (i < 0 || i >= fNparticles) {
451 AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
452 mcParticle = 0;
453 return (mcParticle);
454 }
455
456
415d9f5c 457 //
458 // First check of the MC Particle has been already cached
459 if(!fMCParticleMap->At(i)) {
93df0e9b 460 // Get particle from the stack
461 particle = fStack->Particle(i);
462 // Get track references from Tree TR
463 if (fTreeTR) {
464 fTreeTR->GetEntry(fStack->TreeKEntry(i));
465 trefs = fTRBuffer;
466 ntref = trefs->GetEntriesFast();
467 rarray = new TRefArray(ntref);
468 Int_t nen = fTrackReferences->GetEntriesFast();
469 for (Int_t j = 0; j < ntref; j++) {
470 // Save the track references in a TClonesArray
471 AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
472 // Save the pointer in a TRefArray
473 new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
474 rarray->AddAt((*fTrackReferences)[nen], j);
475 nen++;
476 } // loop over track references for entry i
477 } // if TreeTR available
415d9f5c 478 Int_t nentries = fMCParticles->GetEntriesFast();
eee13e8d 479 new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
415d9f5c 480 mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
481 fMCParticleMap->AddAt(mcParticle, i);
482 } else {
483 mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
484 }
415d9f5c 485 return mcParticle;
486}
487
83626b6b 488 AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
415d9f5c 489
490ClassImp(AliMCEvent)