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