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