Memory leak corrected.
[u/mrichter/AliRoot.git] / STEER / AliMCEventHandler.cxx
CommitLineData
5fe09262 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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// Class AliMCEventHandler
19// This class gives access to MC truth during the analysis.
20// Monte Carlo truth is containe in the kinematics tree (produced particles) and
21// the tree of reference hits.
22//
23// Origin: Andreas Morsch, CERN, andreas.morsch@cern.ch
24//---------------------------------------------------------------------------------
25
26
27
28#include "AliMCEventHandler.h"
29#include "AliTrackReference.h"
30#include "AliHeader.h"
31#include "AliStack.h"
2c36081e 32#include "AliLog.h"
5fe09262 33
34#include <TTree.h>
35#include <TFile.h>
36#include <TParticle.h>
0a05cd41 37#include <TString.h>
5fe09262 38#include <TClonesArray.h>
39#include <TDirectoryFile.h>
40#include <TArrow.h>
41#include <TMarker.h>
42#include <TH2F.h>
43
44
45ClassImp(AliMCEventHandler)
46
47AliMCEventHandler::AliMCEventHandler() :
d2f1d9ef 48 AliVEventHandler(),
5fe09262 49 fFileE(0),
50 fFileK(0),
51 fFileTR(0),
52 fTmpFileTR(0),
53 fTreeE(0),
54 fTreeK(0),
55 fTreeTR(0),
56 fTmpTreeTR(0),
57 fStack(0),
58 fHeader(0),
59 fTrackReferences(0),
60 fNEvent(-1),
61 fEvent(-1),
62 fNprimaries(-1),
63 fNparticles(-1),
0a05cd41 64 fPathName(new TString("./")),
9aea8469 65 fExtension(""),
66 fFileNumber(0),
67 fEventsPerFile(0)
5fe09262 68{
69 // Default constructor
70}
71
72AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
d2f1d9ef 73 AliVEventHandler(name, title),
5fe09262 74 fFileE(0),
75 fFileK(0),
76 fFileTR(0),
77 fTmpFileTR(0),
78 fTreeE(0),
79 fTreeK(0),
80 fTreeTR(0),
81 fTmpTreeTR(0),
82 fStack(0),
83 fHeader(new AliHeader()),
84 fTrackReferences(new TClonesArray("AliTrackReference", 200)),
85 fNEvent(-1),
86 fEvent(-1),
87 fNprimaries(-1),
88 fNparticles(-1),
0a05cd41 89 fPathName(new TString("./")),
9aea8469 90 fExtension(""),
91 fFileNumber(0),
92 fEventsPerFile(0)
5fe09262 93{
94 // Constructor
95}
96AliMCEventHandler::~AliMCEventHandler()
97{
98 // Destructor
99 delete fFileE;
100 delete fFileK;
101 delete fFileTR;
102}
103
104Bool_t AliMCEventHandler::InitIO(Option_t* /*opt*/)
105{
106 // Initialize input
107 //
0a05cd41 108 fFileE = TFile::Open(Form("%sgalice.root", fPathName->Data()));
109 if (!fFileE) AliFatal(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName->Data()));
5fe09262 110
111 fFileE->GetObject("TE", fTreeE);
112 fTreeE->SetBranchAddress("Header", &fHeader);
113 fNEvent = fTreeE->GetEntries();
114 //
115 // Tree K
0a05cd41 116 fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
2c36081e 117 if (!fFileK) AliFatal(Form("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName));
9aea8469 118 fEventsPerFile = fFileK->GetNkeys() - fFileK->GetNProcessIDs();
5fe09262 119 //
120 // Tree TR
0a05cd41 121 fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
122 if (!fFileTR) AliWarning(Form("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName->Data()));
5fe09262 123 //
124 // Reset the event number
2c36081e 125 fEvent = -1;
126 fFileNumber = 0;
9aea8469 127
2c36081e 128 AliInfo(Form("AliMCEventHandler:Number of events in this directory %5d \n", fNEvent));
5fe09262 129 return kTRUE;
130
131}
132
9aea8469 133Bool_t AliMCEventHandler::GetEvent(Int_t iev)
134{
135 // Load the event number iev
2c36081e 136 //
137 // Calculate the file number
9aea8469 138 Int_t inew = iev/fEventsPerFile;
139 if (inew != fFileNumber) {
140 fFileNumber = inew;
141 if (!OpenFile(fFileNumber)){
142 return kFALSE;
143 }
144 }
2c36081e 145 // Folder name
5fe09262 146 char folder[20];
9aea8469 147 sprintf(folder, "Event%d", iev);
2c36081e 148
5fe09262 149 // TreeE
9aea8469 150 fTreeE->GetEntry(iev);
5fe09262 151 fStack = fHeader->Stack();
152 // Tree K
153 TDirectoryFile* dirK = 0;
154 fFileK->GetObject(folder, dirK);
9aea8469 155 if (!dirK) {
2c36081e 156 AliWarning(Form("AliMCEventHandler: Event #%5d not found\n", iev));
9aea8469 157 return kFALSE;
158 }
5fe09262 159 dirK->GetObject("TreeK", fTreeK);
160 fStack->ConnectTree(fTreeK);
161 fStack->GetEvent();
2d8f26f6 162
5fe09262 163 //Tree TR
2d8f26f6 164
07bd4750 165 if (fFileTR) {
2c36081e 166 TDirectoryFile* dirTR = 0;
167 fFileTR->GetObject(folder, dirTR);
168 dirTR->GetObject("TreeTR", fTreeTR);
169 if (fTreeTR->GetBranch("AliRun")) {
170 // This is an old format with one branch per detector not in synch with TreeK
171 ReorderAndExpandTreeTR();
2d8f26f6 172 delete dirTR;
2c36081e 173 } else {
174 // New format
175 fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
176 }
5fe09262 177 }
2d8f26f6 178
5fe09262 179 //
180 fNparticles = fStack->GetNtrack();
181 fNprimaries = fStack->GetNprimary();
2d8f26f6 182 AliInfo(Form("AliMCEventHandler: Event#: %5d Number of particles: %5d (all) %5d (primaries)\n",
183 fEvent, fNparticles, fNprimaries));
5fe09262 184
185 return kTRUE;
9aea8469 186}
187
188Bool_t AliMCEventHandler::OpenFile(Int_t i)
189{
190 // Open file i
191 Bool_t ok = kTRUE;
192 if (i > 0) {
193 fExtension = Form("%d", i);
194 } else {
195 fExtension = "";
196 }
197
198
199 delete fFileK;
0a05cd41 200 fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
9aea8469 201 if (!fFileK) {
0a05cd41 202 AliFatal(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
9aea8469 203 ok = kFALSE;
204 }
205
206 delete fFileTR;
0a05cd41 207 fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
9aea8469 208 if (!fFileTR) {
0a05cd41 209 AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
9aea8469 210 ok = kFALSE;
211 }
212
213 return ok;
214}
215
216Bool_t AliMCEventHandler::BeginEvent()
217{
218 // Read the next event
219 fEvent++;
220 if (fEvent >= fNEvent) {
2c36081e 221 AliWarning(Form("AliMCEventHandler: Event number out of range %5d\n", fEvent));
9aea8469 222 return kFALSE;
223 }
224 return GetEvent(fEvent);
5fe09262 225}
226
227Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
228{
229 // Retrieve entry i
2d8f26f6 230 if (i < 0 || i >= fNparticles) {
2c36081e 231 AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
232 particle = 0;
233 trefs = 0;
234 return (-1);
235 }
236 particle = fStack->Particle(i);
237 if (fFileTR) {
5fe09262 238 fTreeTR->GetEntry(fStack->TreeKEntry(i));
2c36081e 239 trefs = fTrackReferences;
240 return trefs->GetEntries();
5fe09262 241 } else {
2c36081e 242 trefs = 0;
243 return -1;
5fe09262 244 }
5fe09262 245}
246
2d8f26f6 247void AliMCEventHandler::DrawCheck(Int_t i, Int_t search)
5fe09262 248{
249 // Retrieve entry i and draw momentum vector and hits
2c36081e 250 if (!fFileTR) {
251 AliWarning("No Track Reference information available");
252 return;
253 }
254
5fe09262 255 if (i > -1 && i < fNparticles) {
256 fTreeTR->GetEntry(fStack->TreeKEntry(i));
257 } else {
2c36081e 258 AliWarning("AliMCEventHandler::GetEntry: Index out of range");
5fe09262 259 }
2c36081e 260
36e4cfbb 261 Int_t nh = fTrackReferences->GetEntries();
262
5fe09262 263
36e4cfbb 264 if (search) {
2d8f26f6 265 while(nh <= search && i < fNparticles - 1) {
36e4cfbb 266 i++;
267 fTreeTR->GetEntry(fStack->TreeKEntry(i));
268 nh = fTrackReferences->GetEntries();
269 }
270 printf("Found Hits at %5d\n", i);
271 }
272 TParticle* particle = fStack->Particle(i);
273
5fe09262 274 TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
275 Float_t x0 = particle->Vx();
276 Float_t y0 = particle->Vy();
277
278 Float_t x1 = particle->Vx() + particle->Px() * 50.;
279 Float_t y1 = particle->Vy() + particle->Py() * 50.;
280
281 TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
282 h->Draw();
283 a->SetLineColor(2);
284
285 a->Draw();
286
287 for (Int_t ih = 0; ih < nh; ih++) {
288 AliTrackReference* ref = (AliTrackReference*) fTrackReferences->At(ih);
289 TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
290 m->Draw();
291 m->SetMarkerSize(0.4);
292
293 }
294}
295
890126ab 296Bool_t AliMCEventHandler::Notify(const char *path)
5fe09262 297{
298 // Notify about directory change
890126ab 299 // The directory is taken from the 'path' argument
36e4cfbb 300 // Reconnect trees
890126ab 301
302 printf("AliMCEventHandler::Notify() file: %s\n", path);
0a05cd41 303 delete fPathName;
304 fPathName = new TString(path);
5fe09262 305 ResetIO();
306 InitIO("");
307 return kTRUE;
308}
309
310void AliMCEventHandler::ResetIO()
311{
36e4cfbb 312 // Reset files
5fe09262 313 if (fFileE) delete fFileE;
314 if (fFileK) delete fFileK;
315 if (fFileTR) delete fFileTR;
316}
317
318
319Bool_t AliMCEventHandler::FinishEvent()
320{
2d8f26f6 321 // Reset the stack
2d8f26f6 322 Stack()->Reset();
323
5fe09262 324 return kTRUE;
325}
326
327Bool_t AliMCEventHandler::Terminate()
328{
329 // Dummy
330 return kTRUE;
331}
332
333Bool_t AliMCEventHandler::TerminateIO()
334{
335 // Dummy
336 return kTRUE;
337}
338
339void AliMCEventHandler::ReorderAndExpandTreeTR()
340{
341//
342// Reorder and expand the track reference tree in order to match the kinematics tree.
343// Copy the information from different branches into one
344//
345// TreeTR
2d8f26f6 346 if (fTmpTreeTR) {
347 fTmpTreeTR->Delete("all");
348 }
349
5fe09262 350 if (fTmpFileTR) {
351 fTmpFileTR->Close();
352 delete fTmpFileTR;
353 }
354
355 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
356 fTmpTreeTR = new TTree("TreeTR", "Track References");
357 if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
358 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 4000);
2d8f26f6 359
5fe09262 360//
67dcc1eb 361//
362// fTreeTR->SetBranchStatus("*", 0);
363 fTreeTR->SetBranchStatus("ABSO", 0);
364 fTreeTR->SetBranchStatus("BODY", 0);
365 fTreeTR->SetBranchStatus("DIPO", 0);
366 fTreeTR->SetBranchStatus("EMCAL", 0);
367 fTreeTR->SetBranchStatus("FMD", 0);
368 fTreeTR->SetBranchStatus("HALL", 0);
369 fTreeTR->SetBranchStatus("MAG", 0);
370 fTreeTR->SetBranchStatus("PHOS", 0);
371 fTreeTR->SetBranchStatus("PIPE", 0);
372 fTreeTR->SetBranchStatus("PMD", 0);
373 fTreeTR->SetBranchStatus("RICH", 0);
374 fTreeTR->SetBranchStatus("SHIL", 0);
375 fTreeTR->SetBranchStatus("START", 0);
376 fTreeTR->SetBranchStatus("VZERO", 0);
377 fTreeTR->SetBranchStatus("ZDC", 0);
378
808ba8bf 379 fTreeTR->SetBranchStatus("AliRun", 1);
380 fTreeTR->SetBranchStatus("ITS", 1);
381 fTreeTR->SetBranchStatus("TPC", 1);
382 fTreeTR->SetBranchStatus("TRD", 1);
383 fTreeTR->SetBranchStatus("TOF", 1);
384 fTreeTR->SetBranchStatus("FRAME", 1);
385 fTreeTR->SetBranchStatus("MUON", 1);
386
5fe09262 387 TClonesArray* trefs[7];
388 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
389 if (fTreeTR){
390 // make branch for central track references
391 if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
392 if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
393 if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
394 if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
395 if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
396 if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
397 if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
398 }
399
400 Int_t np = fStack->GetNprimary();
401 Int_t nt = fTreeTR->GetEntries();
402 //
403 // Loop over tracks and find the secondaries with the help of the kine tree
404 Int_t ifills = 0;
405 Int_t it = 0;
406 Int_t itlast = 0;
407
408 for (Int_t ip = np - 1; ip > -1; ip--) {
409 TParticle *part = fStack->Particle(ip);
410// printf("Particle %5d %5d %5d %5d %5d %5d \n",
411// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
412// part->GetLastDaughter(), part->TestBit(kTransportBit));
36e4cfbb 413
414 // Determine range of secondaries produced by this primary during transport
5fe09262 415 Int_t dau1 = part->GetFirstDaughter();
36e4cfbb 416 if (dau1 < np) continue; // This particle has no secondaries produced during transport
5fe09262 417 Int_t dau2 = -1;
5fe09262 418 if (dau1 > -1) {
419 Int_t inext = ip - 1;
420 while (dau2 < 0) {
421 if (inext >= 0) {
422 part = fStack->Particle(inext);
423 dau2 = part->GetFirstDaughter();
424 if (dau2 == -1 || dau2 < np) {
425 dau2 = -1;
426 } else {
427 dau2--;
428 }
429 } else {
430 dau2 = fStack->GetNtrack() - 1;
431 }
432 inext--;
433 } // find upper bound
434 } // dau2 < 0
36e4cfbb 435
5fe09262 436
437// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
36e4cfbb 438//
439// Loop over reference hits and find secondary label
440// First the tricky part: find the entry in treeTR than contains the hits or
441// make sure that no hits exist.
442//
5fe09262 443 Bool_t hasHits = kFALSE;
444 Bool_t isOutside = kFALSE;
36e4cfbb 445
5fe09262 446 it = itlast;
5fe09262 447 while (!hasHits && !isOutside && it < nt) {
448 fTreeTR->GetEntry(it++);
449 for (Int_t ib = 0; ib < 7; ib++) {
450 if (!trefs[ib]) continue;
451 Int_t nh = trefs[ib]->GetEntries();
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 >= dau1 && label <= dau2) {
456 hasHits = kTRUE;
457 itlast = it - 1;
36e4cfbb 458 break;
5fe09262 459 }
460 if (label > dau2 || label < ip) {
461 isOutside = kTRUE;
462 itlast = it - 1;
36e4cfbb 463 break;
5fe09262 464 }
465 } // hits
36e4cfbb 466 if (hasHits || isOutside) break;
5fe09262 467 } // branches
468 } // entries
469
470 if (!hasHits) {
36e4cfbb 471 // Write empty entries
5fe09262 472 for (Int_t id = dau1; (id <= dau2); id++) {
473 fTmpTreeTR->Fill();
474 ifills++;
475 }
476 } else {
36e4cfbb 477 // Collect all hits
5fe09262 478 fTreeTR->GetEntry(itlast);
479 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
480 for (Int_t ib = 0; ib < 7; ib++) {
481 if (!trefs[ib]) continue;
482 Int_t nh = trefs[ib]->GetEntries();
483 for (Int_t ih = 0; ih < nh; ih++) {
484 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
485 Int_t label = tr->Label();
486 // Skip primaries
487 if (label == ip) continue;
488 if (label > dau2 || label < dau1)
36e4cfbb 489 printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
490 itlast, label, dau1, dau2);
5fe09262 491 if (label == id) {
492 // secondary found
493 tr->SetDetectorId(ib-1);
494 Int_t nref = fTrackReferences->GetEntriesFast();
495 TClonesArray &lref = *fTrackReferences;
496 new(lref[nref]) AliTrackReference(*tr);
497 }
498 } // hits
499 } // branches
500 fTmpTreeTR->Fill();
501 fTrackReferences->Clear();
502 ifills++;
503 } // daughters
504 } // has hits
505 } // tracks
5fe09262 506 //
507 // Now loop again and write the primaries
36e4cfbb 508 //
5fe09262 509 it = nt - 1;
510 for (Int_t ip = 0; ip < np; ip++) {
511 Int_t labmax = -1;
512 while (labmax < ip && it > -1) {
513 fTreeTR->GetEntry(it--);
514 for (Int_t ib = 0; ib < 7; ib++) {
515 if (!trefs[ib]) continue;
516 Int_t nh = trefs[ib]->GetEntries();
517 //
518 // Loop over reference hits and find primary labels
519 for (Int_t ih = 0; ih < nh; ih++) {
520 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
521 Int_t label = tr->Label();
522 if (label < np && label > labmax) {
523 labmax = label;
524 }
525
526 if (label == ip) {
527 tr->SetDetectorId(ib-1);
528 Int_t nref = fTrackReferences->GetEntriesFast();
529 TClonesArray &lref = *fTrackReferences;
530 new(lref[nref]) AliTrackReference(*tr);
531 }
532 } // hits
533 } // branches
534 } // entries
535 it++;
536 fTmpTreeTR->Fill();
537 fTrackReferences->Clear();
538 ifills++;
539 } // tracks
540 // Check
808ba8bf 541 delete fTreeTR;
542 for (Int_t ib = 0; ib < 7; ib++) {
543 if (trefs[ib]) delete trefs[ib];
544 }
545
5fe09262 546 if (ifills != fStack->GetNtrack())
36e4cfbb 547 printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
548 ifills, fStack->GetNtrack());
5fe09262 549 fTreeTR = fTmpTreeTR;
5fe09262 550}
0a05cd41 551
552void AliMCEventHandler::SetInputPath(char* fname)
553{
554 // Set the input path name
555 delete fPathName;
556 fPathName = new TString(fname);
557}