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