Applying line fit to trigger record hits to get better idea of where the track traver...
[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{
53faeca4 302 // Notify about directory change
303 // The directory is taken from the 'path' argument
304 // Reconnect trees
305 TString fileName(path);
306 if(fileName.Contains("AliESDs.root")){
307 fileName.ReplaceAll("AliESDs.root", "");
308 }
309 else if(fileName.Contains("galice.root")){
310 // for running with galice and kinematics alone...
311 fileName.ReplaceAll("galice.root", "");
312 }
313
314 *fPathName = fileName;
315 printf("AliMCEventHandler::Notify() Path: %s\n", fPathName->Data());
316
5fe09262 317 ResetIO();
318 InitIO("");
319 return kTRUE;
320}
321
322void AliMCEventHandler::ResetIO()
323{
5efedd31 324// Clear header and stack
325
326 if (fHeader) {
327 delete fHeader;
328 fHeader = 0;
329 }
330
331 delete fStack;
332 delete fTreeE; fTreeE = 0;
333
334// Clear TR
335
336 if (fTrackReferences) {
337 fTrackReferences->Clear();
338 delete fTrackReferences;
339 fTrackReferences = 0;
340 }
341
342// Reset files
5fe09262 343 if (fFileE) delete fFileE;
344 if (fFileK) delete fFileK;
345 if (fFileTR) delete fFileTR;
346}
347
348
349Bool_t AliMCEventHandler::FinishEvent()
350{
5efedd31 351 // Clean-up after each event
352 delete fDirTR; fDirTR = 0;
353 delete fDirK; fDirK = 0;
354 Stack()->Reset(0);
5fe09262 355 return kTRUE;
356}
357
358Bool_t AliMCEventHandler::Terminate()
359{
360 // Dummy
361 return kTRUE;
362}
363
364Bool_t AliMCEventHandler::TerminateIO()
365{
366 // Dummy
367 return kTRUE;
368}
369
370void AliMCEventHandler::ReorderAndExpandTreeTR()
371{
372//
373// Reorder and expand the track reference tree in order to match the kinematics tree.
374// Copy the information from different branches into one
375//
376// TreeTR
5fe09262 377
378 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
5efedd31 379 fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
5fe09262 380 if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
5efedd31 381 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 32000, 0);
382
2d8f26f6 383
5fe09262 384//
5efedd31 385// Activate the used branches only. Otherwisw we get a bad memory leak.
386 fTreeTR->SetBranchStatus("*", 0);
387 fTreeTR->SetBranchStatus("AliRun.*", 1);
388 fTreeTR->SetBranchStatus("ITS.*", 1);
389 fTreeTR->SetBranchStatus("TPC.*", 1);
390 fTreeTR->SetBranchStatus("TRD.*", 1);
391 fTreeTR->SetBranchStatus("TOF.*", 1);
392 fTreeTR->SetBranchStatus("FRAME.*", 1);
393 fTreeTR->SetBranchStatus("MUON.*", 1);
67dcc1eb 394//
5efedd31 395// Connect the active branches
5fe09262 396 TClonesArray* trefs[7];
397 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
398 if (fTreeTR){
399 // make branch for central track references
400 if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
401 if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
402 if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
403 if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
404 if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
405 if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
406 if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
407 }
408
409 Int_t np = fStack->GetNprimary();
410 Int_t nt = fTreeTR->GetEntries();
5efedd31 411
5fe09262 412 //
413 // Loop over tracks and find the secondaries with the help of the kine tree
414 Int_t ifills = 0;
415 Int_t it = 0;
416 Int_t itlast = 0;
5efedd31 417 TParticle* part;
418
5fe09262 419 for (Int_t ip = np - 1; ip > -1; ip--) {
5efedd31 420 part = fStack->Particle(ip);
5fe09262 421// printf("Particle %5d %5d %5d %5d %5d %5d \n",
422// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
423// part->GetLastDaughter(), part->TestBit(kTransportBit));
36e4cfbb 424
425 // Determine range of secondaries produced by this primary during transport
5fe09262 426 Int_t dau1 = part->GetFirstDaughter();
36e4cfbb 427 if (dau1 < np) continue; // This particle has no secondaries produced during transport
5fe09262 428 Int_t dau2 = -1;
5fe09262 429 if (dau1 > -1) {
430 Int_t inext = ip - 1;
431 while (dau2 < 0) {
432 if (inext >= 0) {
433 part = fStack->Particle(inext);
434 dau2 = part->GetFirstDaughter();
435 if (dau2 == -1 || dau2 < np) {
436 dau2 = -1;
437 } else {
438 dau2--;
439 }
440 } else {
441 dau2 = fStack->GetNtrack() - 1;
442 }
443 inext--;
444 } // find upper bound
445 } // dau2 < 0
36e4cfbb 446
5fe09262 447
448// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
36e4cfbb 449//
450// Loop over reference hits and find secondary label
451// First the tricky part: find the entry in treeTR than contains the hits or
452// make sure that no hits exist.
453//
5fe09262 454 Bool_t hasHits = kFALSE;
455 Bool_t isOutside = kFALSE;
36e4cfbb 456
5fe09262 457 it = itlast;
5fe09262 458 while (!hasHits && !isOutside && it < nt) {
459 fTreeTR->GetEntry(it++);
460 for (Int_t ib = 0; ib < 7; ib++) {
461 if (!trefs[ib]) continue;
462 Int_t nh = trefs[ib]->GetEntries();
463 for (Int_t ih = 0; ih < nh; ih++) {
464 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
465 Int_t label = tr->Label();
466 if (label >= dau1 && label <= dau2) {
467 hasHits = kTRUE;
468 itlast = it - 1;
36e4cfbb 469 break;
5fe09262 470 }
471 if (label > dau2 || label < ip) {
472 isOutside = kTRUE;
473 itlast = it - 1;
36e4cfbb 474 break;
5fe09262 475 }
476 } // hits
36e4cfbb 477 if (hasHits || isOutside) break;
5fe09262 478 } // branches
479 } // entries
480
481 if (!hasHits) {
36e4cfbb 482 // Write empty entries
5fe09262 483 for (Int_t id = dau1; (id <= dau2); id++) {
484 fTmpTreeTR->Fill();
485 ifills++;
486 }
487 } else {
36e4cfbb 488 // Collect all hits
5fe09262 489 fTreeTR->GetEntry(itlast);
490 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
491 for (Int_t ib = 0; ib < 7; ib++) {
492 if (!trefs[ib]) continue;
493 Int_t nh = trefs[ib]->GetEntries();
494 for (Int_t ih = 0; ih < nh; ih++) {
495 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
496 Int_t label = tr->Label();
497 // Skip primaries
498 if (label == ip) continue;
499 if (label > dau2 || label < dau1)
36e4cfbb 500 printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
501 itlast, label, dau1, dau2);
5fe09262 502 if (label == id) {
503 // secondary found
504 tr->SetDetectorId(ib-1);
505 Int_t nref = fTrackReferences->GetEntriesFast();
506 TClonesArray &lref = *fTrackReferences;
507 new(lref[nref]) AliTrackReference(*tr);
508 }
509 } // hits
510 } // branches
511 fTmpTreeTR->Fill();
512 fTrackReferences->Clear();
513 ifills++;
514 } // daughters
515 } // has hits
516 } // tracks
5efedd31 517
5fe09262 518 //
519 // Now loop again and write the primaries
36e4cfbb 520 //
5fe09262 521 it = nt - 1;
522 for (Int_t ip = 0; ip < np; ip++) {
523 Int_t labmax = -1;
524 while (labmax < ip && it > -1) {
525 fTreeTR->GetEntry(it--);
526 for (Int_t ib = 0; ib < 7; ib++) {
527 if (!trefs[ib]) continue;
528 Int_t nh = trefs[ib]->GetEntries();
529 //
530 // Loop over reference hits and find primary labels
531 for (Int_t ih = 0; ih < nh; ih++) {
532 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
533 Int_t label = tr->Label();
534 if (label < np && label > labmax) {
535 labmax = label;
536 }
537
538 if (label == ip) {
539 tr->SetDetectorId(ib-1);
540 Int_t nref = fTrackReferences->GetEntriesFast();
541 TClonesArray &lref = *fTrackReferences;
542 new(lref[nref]) AliTrackReference(*tr);
543 }
544 } // hits
545 } // branches
546 } // entries
547 it++;
548 fTmpTreeTR->Fill();
549 fTrackReferences->Clear();
550 ifills++;
551 } // tracks
552 // Check
5efedd31 553
554
555 // Clean-up
556 delete fTreeTR; fTreeTR = 0;
557
808ba8bf 558 for (Int_t ib = 0; ib < 7; ib++) {
5efedd31 559 if (trefs[ib]) {
560 trefs[ib]->Clear();
561 delete trefs[ib];
562 trefs[ib] = 0;
563 }
808ba8bf 564 }
565
5fe09262 566 if (ifills != fStack->GetNtrack())
36e4cfbb 567 printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
568 ifills, fStack->GetNtrack());
5efedd31 569
570 fTmpTreeTR->Write();
5fe09262 571 fTreeTR = fTmpTreeTR;
5fe09262 572}
0a05cd41 573
574void AliMCEventHandler::SetInputPath(char* fname)
575{
576 // Set the input path name
577 delete fPathName;
578 fPathName = new TString(fname);
579}