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 | |
44 | ClassImp(AliMCEventHandler) |
45 | |
46 | AliMCEventHandler::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 | |
68 | AliMCEventHandler::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 | } |
89 | AliMCEventHandler::~AliMCEventHandler() |
90 | { |
91 | // Destructor |
92 | delete fFileE; |
93 | delete fFileK; |
94 | delete fFileTR; |
95 | } |
96 | |
97 | Bool_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 | |
123 | Bool_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 | |
158 | Int_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 |
171 | void 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 | |
215 | Bool_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 | |
231 | void 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 | |
240 | Bool_t AliMCEventHandler::FinishEvent() |
241 | { |
242 | // Dummy |
243 | return kTRUE; |
244 | } |
245 | |
246 | Bool_t AliMCEventHandler::Terminate() |
247 | { |
248 | // Dummy |
249 | return kTRUE; |
250 | } |
251 | |
252 | Bool_t AliMCEventHandler::TerminateIO() |
253 | { |
254 | // Dummy |
255 | return kTRUE; |
256 | } |
257 | |
258 | void 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 | } |