Fixed phi calculation to ALICE convention
[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//
808ba8bf 361 fTreeTR->SetBranchStatus("*", 0);
362 fTreeTR->SetBranchStatus("AliRun", 1);
363 fTreeTR->SetBranchStatus("ITS", 1);
364 fTreeTR->SetBranchStatus("TPC", 1);
365 fTreeTR->SetBranchStatus("TRD", 1);
366 fTreeTR->SetBranchStatus("TOF", 1);
367 fTreeTR->SetBranchStatus("FRAME", 1);
368 fTreeTR->SetBranchStatus("MUON", 1);
369
5fe09262 370 TClonesArray* trefs[7];
371 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
372 if (fTreeTR){
373 // make branch for central track references
374 if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
375 if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
376 if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
377 if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
378 if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
379 if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
380 if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
381 }
382
383 Int_t np = fStack->GetNprimary();
384 Int_t nt = fTreeTR->GetEntries();
385 //
386 // Loop over tracks and find the secondaries with the help of the kine tree
387 Int_t ifills = 0;
388 Int_t it = 0;
389 Int_t itlast = 0;
390
391 for (Int_t ip = np - 1; ip > -1; ip--) {
392 TParticle *part = fStack->Particle(ip);
393// printf("Particle %5d %5d %5d %5d %5d %5d \n",
394// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
395// part->GetLastDaughter(), part->TestBit(kTransportBit));
36e4cfbb 396
397 // Determine range of secondaries produced by this primary during transport
5fe09262 398 Int_t dau1 = part->GetFirstDaughter();
36e4cfbb 399 if (dau1 < np) continue; // This particle has no secondaries produced during transport
5fe09262 400 Int_t dau2 = -1;
5fe09262 401 if (dau1 > -1) {
402 Int_t inext = ip - 1;
403 while (dau2 < 0) {
404 if (inext >= 0) {
405 part = fStack->Particle(inext);
406 dau2 = part->GetFirstDaughter();
407 if (dau2 == -1 || dau2 < np) {
408 dau2 = -1;
409 } else {
410 dau2--;
411 }
412 } else {
413 dau2 = fStack->GetNtrack() - 1;
414 }
415 inext--;
416 } // find upper bound
417 } // dau2 < 0
36e4cfbb 418
5fe09262 419
420// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
36e4cfbb 421//
422// Loop over reference hits and find secondary label
423// First the tricky part: find the entry in treeTR than contains the hits or
424// make sure that no hits exist.
425//
5fe09262 426 Bool_t hasHits = kFALSE;
427 Bool_t isOutside = kFALSE;
36e4cfbb 428
5fe09262 429 it = itlast;
5fe09262 430 while (!hasHits && !isOutside && it < nt) {
431 fTreeTR->GetEntry(it++);
432 for (Int_t ib = 0; ib < 7; ib++) {
433 if (!trefs[ib]) continue;
434 Int_t nh = trefs[ib]->GetEntries();
435 for (Int_t ih = 0; ih < nh; ih++) {
436 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
437 Int_t label = tr->Label();
438 if (label >= dau1 && label <= dau2) {
439 hasHits = kTRUE;
440 itlast = it - 1;
36e4cfbb 441 break;
5fe09262 442 }
443 if (label > dau2 || label < ip) {
444 isOutside = kTRUE;
445 itlast = it - 1;
36e4cfbb 446 break;
5fe09262 447 }
448 } // hits
36e4cfbb 449 if (hasHits || isOutside) break;
5fe09262 450 } // branches
451 } // entries
452
453 if (!hasHits) {
36e4cfbb 454 // Write empty entries
5fe09262 455 for (Int_t id = dau1; (id <= dau2); id++) {
456 fTmpTreeTR->Fill();
457 ifills++;
458 }
459 } else {
36e4cfbb 460 // Collect all hits
5fe09262 461 fTreeTR->GetEntry(itlast);
462 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
463 for (Int_t ib = 0; ib < 7; ib++) {
464 if (!trefs[ib]) continue;
465 Int_t nh = trefs[ib]->GetEntries();
466 for (Int_t ih = 0; ih < nh; ih++) {
467 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
468 Int_t label = tr->Label();
469 // Skip primaries
470 if (label == ip) continue;
471 if (label > dau2 || label < dau1)
36e4cfbb 472 printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
473 itlast, label, dau1, dau2);
5fe09262 474 if (label == id) {
475 // secondary found
476 tr->SetDetectorId(ib-1);
477 Int_t nref = fTrackReferences->GetEntriesFast();
478 TClonesArray &lref = *fTrackReferences;
479 new(lref[nref]) AliTrackReference(*tr);
480 }
481 } // hits
482 } // branches
483 fTmpTreeTR->Fill();
484 fTrackReferences->Clear();
485 ifills++;
486 } // daughters
487 } // has hits
488 } // tracks
5fe09262 489 //
490 // Now loop again and write the primaries
36e4cfbb 491 //
5fe09262 492 it = nt - 1;
493 for (Int_t ip = 0; ip < np; ip++) {
494 Int_t labmax = -1;
495 while (labmax < ip && it > -1) {
496 fTreeTR->GetEntry(it--);
497 for (Int_t ib = 0; ib < 7; ib++) {
498 if (!trefs[ib]) continue;
499 Int_t nh = trefs[ib]->GetEntries();
500 //
501 // Loop over reference hits and find primary labels
502 for (Int_t ih = 0; ih < nh; ih++) {
503 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
504 Int_t label = tr->Label();
505 if (label < np && label > labmax) {
506 labmax = label;
507 }
508
509 if (label == ip) {
510 tr->SetDetectorId(ib-1);
511 Int_t nref = fTrackReferences->GetEntriesFast();
512 TClonesArray &lref = *fTrackReferences;
513 new(lref[nref]) AliTrackReference(*tr);
514 }
515 } // hits
516 } // branches
517 } // entries
518 it++;
519 fTmpTreeTR->Fill();
520 fTrackReferences->Clear();
521 ifills++;
522 } // tracks
523 // Check
808ba8bf 524 delete fTreeTR;
525 for (Int_t ib = 0; ib < 7; ib++) {
526 if (trefs[ib]) delete trefs[ib];
527 }
528
5fe09262 529 if (ifills != fStack->GetNtrack())
36e4cfbb 530 printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
531 ifills, fStack->GetNtrack());
5fe09262 532 fTreeTR = fTmpTreeTR;
5fe09262 533}
0a05cd41 534
535void AliMCEventHandler::SetInputPath(char* fname)
536{
537 // Set the input path name
538 delete fPathName;
539 fPathName = new TString(fname);
540}