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