The MC Event part has been separated from the AliMCEventHandler.
[u/mrichter/AliRoot.git] / STEER / AliMCEvent.cxx
CommitLineData
415d9f5c 1/**************************************************************************
2 * Copyright(c) 1998-2007, 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//-------------------------------------------------------------------------
19// Class for Kinematic Events
20// Author: Andreas Morsch, CERN
21//-------------------------------------------------------------------------
22#include <TArrow.h>
23#include <TMarker.h>
24#include <TH2F.h>
25#include <TTree.h>
26#include <TFile.h>
27#include <TParticle.h>
28#include <TClonesArray.h>
29#include <TObjArray.h>
30
31#include "AliLog.h"
32#include "AliMCEvent.h"
33#include "AliStack.h"
34#include "AliTrackReference.h"
35#include "AliHeader.h"
36
37
38AliMCEvent::AliMCEvent():
39 AliVEvent(),
40 fStack(0),
41 fMCParticles(new TClonesArray("AliMCParticle",1000)),
42 fMCParticleMap(0),
43 fHeader(0),
44 fTrackReferences(0),
45 fTreeTR(0),
46 fTmpTreeTR(0),
47 fTmpFileTR(0),
48 fNprimaries(-1),
49 fNparticles(-1)
50{
51 // Default constructor
52}
53
54AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
55 AliVEvent(mcEvnt),
56 fStack(0),
57 fMCParticles(0),
58 fMCParticleMap(0),
59 fHeader(0),
60 fTrackReferences(0),
61 fTreeTR(0),
62 fTmpTreeTR(0),
63 fTmpFileTR(0),
64 fNprimaries(-1),
65 fNparticles(-1)
66{
67// Copy constructor
68}
69
70
71AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
72{
73 // assignment operator
74 if (this!=&mcEvnt) {
75 AliVEvent::operator=(mcEvnt);
76 }
77
78 return *this;
79}
80
81void AliMCEvent::ConnectTreeE (TTree* tree)
82{
83 // Connect the event header tree
84 tree->SetBranchAddress("Header", &fHeader);
85}
86
87void AliMCEvent::ConnectTreeK (TTree* tree)
88{
89 // Connect the kinematics tree to the stack
90 fStack = fHeader->Stack();
91 fStack->ConnectTree(tree);
92 //
93 // Load the event
94 fStack->GetEvent();
95 fNparticles = fStack->GetNtrack();
96 fNprimaries = fStack->GetNprimary();
97 AliInfo(Form("AliMCEvent: Number of particles: %5d (all) %5d (primaries)\n",
98 fNparticles, fNprimaries));
99
100 // This is a cache for the TParticles converted to MCParticles on user request
101 if (fMCParticleMap) {
102 fMCParticleMap->Clear();
103 if (fNparticles>0) fMCParticleMap->Expand(fNparticles);}
104 else
105 fMCParticleMap = new TObjArray(fNparticles);
106 fMCParticles->Clear();
107}
108
109void AliMCEvent::ConnectTreeTR (TTree* tree)
110{
111 // Connect the track reference tree
112 fTreeTR = tree;
113
114 if (fTreeTR->GetBranch("AliRun")) {
115 if (fTmpFileTR) {
116 fTmpFileTR->Close();
117 delete fTmpFileTR;
118 }
119 // This is an old format with one branch per detector not in synch with TreeK
120 ReorderAndExpandTreeTR();
121 } else {
122 // New format
123 fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
124 }
125}
126
127Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
128{
129 // Retrieve entry i
130 if (i < 0 || i >= fNparticles) {
131 AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
132 particle = 0;
133 trefs = 0;
134 return (-1);
135 }
136 particle = fStack->Particle(i);
137 if (fTreeTR) {
138 fTreeTR->GetEntry(fStack->TreeKEntry(i));
139 trefs = fTrackReferences;
140 return trefs->GetEntries();
141 } else {
142 trefs = 0;
143 return -1;
144 }
145}
146
147
148void AliMCEvent::Clean()
149{
150 // Clean-up before new trees are connected
151
152 if (fHeader) {
153 delete fHeader;
154 fHeader = 0;
155 }
156
157 delete fStack;
158
159 // Clear TR
160 if (fTrackReferences) {
161 fTrackReferences->Clear();
162 delete fTrackReferences;
163 fTrackReferences = 0;
164 }
165}
166
167void AliMCEvent::FinishEvent()
168{
169 // Clean-up after event
170 fStack->Reset(0);
171}
172
173
174void AliMCEvent::DrawCheck(Int_t i, Int_t search)
175{
176 //
177 // Simple event display for debugging
178 if (!fTreeTR) {
179 AliWarning("No Track Reference information available");
180 return;
181 }
182
183 if (i > -1 && i < fNparticles) {
184 fTreeTR->GetEntry(fStack->TreeKEntry(i));
185 } else {
186 AliWarning("AliMCEvent::GetEntry: Index out of range");
187 }
188
189 Int_t nh = fTrackReferences->GetEntries();
190
191
192 if (search) {
193 while(nh <= search && i < fNparticles - 1) {
194 i++;
195 fTreeTR->GetEntry(fStack->TreeKEntry(i));
196 nh = fTrackReferences->GetEntries();
197 }
198 printf("Found Hits at %5d\n", i);
199 }
200 TParticle* particle = fStack->Particle(i);
201
202 TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
203 Float_t x0 = particle->Vx();
204 Float_t y0 = particle->Vy();
205
206 Float_t x1 = particle->Vx() + particle->Px() * 50.;
207 Float_t y1 = particle->Vy() + particle->Py() * 50.;
208
209 TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
210 h->Draw();
211 a->SetLineColor(2);
212
213 a->Draw();
214
215 for (Int_t ih = 0; ih < nh; ih++) {
216 AliTrackReference* ref = (AliTrackReference*) fTrackReferences->At(ih);
217 TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
218 m->Draw();
219 m->SetMarkerSize(0.4);
220
221 }
222}
223
224
225void AliMCEvent::ReorderAndExpandTreeTR()
226{
227//
228// Reorder and expand the track reference tree in order to match the kinematics tree.
229// Copy the information from different branches into one
230//
231// TreeTR
232
233 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
234 fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
235 if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
236 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 32000, 0);
237
238
239//
240// Activate the used branches only. Otherwisw we get a bad memory leak.
241 fTreeTR->SetBranchStatus("*", 0);
242 fTreeTR->SetBranchStatus("AliRun.*", 1);
243 fTreeTR->SetBranchStatus("ITS.*", 1);
244 fTreeTR->SetBranchStatus("TPC.*", 1);
245 fTreeTR->SetBranchStatus("TRD.*", 1);
246 fTreeTR->SetBranchStatus("TOF.*", 1);
247 fTreeTR->SetBranchStatus("FRAME.*", 1);
248 fTreeTR->SetBranchStatus("MUON.*", 1);
249//
250// Connect the active branches
251 TClonesArray* trefs[7];
252 for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
253 if (fTreeTR){
254 // make branch for central track references
255 if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
256 if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
257 if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
258 if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
259 if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
260 if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
261 if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
262 }
263
264 Int_t np = fStack->GetNprimary();
265 Int_t nt = fTreeTR->GetEntries();
266
267 //
268 // Loop over tracks and find the secondaries with the help of the kine tree
269 Int_t ifills = 0;
270 Int_t it = 0;
271 Int_t itlast = 0;
272 TParticle* part;
273
274 for (Int_t ip = np - 1; ip > -1; ip--) {
275 part = fStack->Particle(ip);
276// printf("Particle %5d %5d %5d %5d %5d %5d \n",
277// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
278// part->GetLastDaughter(), part->TestBit(kTransportBit));
279
280 // Determine range of secondaries produced by this primary during transport
281 Int_t dau1 = part->GetFirstDaughter();
282 if (dau1 < np) continue; // This particle has no secondaries produced during transport
283 Int_t dau2 = -1;
284 if (dau1 > -1) {
285 Int_t inext = ip - 1;
286 while (dau2 < 0) {
287 if (inext >= 0) {
288 part = fStack->Particle(inext);
289 dau2 = part->GetFirstDaughter();
290 if (dau2 == -1 || dau2 < np) {
291 dau2 = -1;
292 } else {
293 dau2--;
294 }
295 } else {
296 dau2 = fStack->GetNtrack() - 1;
297 }
298 inext--;
299 } // find upper bound
300 } // dau2 < 0
301
302
303// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
304//
305// Loop over reference hits and find secondary label
306// First the tricky part: find the entry in treeTR than contains the hits or
307// make sure that no hits exist.
308//
309 Bool_t hasHits = kFALSE;
310 Bool_t isOutside = kFALSE;
311
312 it = itlast;
313 while (!hasHits && !isOutside && it < nt) {
314 fTreeTR->GetEntry(it++);
315 for (Int_t ib = 0; ib < 7; ib++) {
316 if (!trefs[ib]) continue;
317 Int_t nh = trefs[ib]->GetEntries();
318 for (Int_t ih = 0; ih < nh; ih++) {
319 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
320 Int_t label = tr->Label();
321 if (label >= dau1 && label <= dau2) {
322 hasHits = kTRUE;
323 itlast = it - 1;
324 break;
325 }
326 if (label > dau2 || label < ip) {
327 isOutside = kTRUE;
328 itlast = it - 1;
329 break;
330 }
331 } // hits
332 if (hasHits || isOutside) break;
333 } // branches
334 } // entries
335
336 if (!hasHits) {
337 // Write empty entries
338 for (Int_t id = dau1; (id <= dau2); id++) {
339 fTmpTreeTR->Fill();
340 ifills++;
341 }
342 } else {
343 // Collect all hits
344 fTreeTR->GetEntry(itlast);
345 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
346 for (Int_t ib = 0; ib < 7; ib++) {
347 if (!trefs[ib]) continue;
348 Int_t nh = trefs[ib]->GetEntries();
349 for (Int_t ih = 0; ih < nh; ih++) {
350 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
351 Int_t label = tr->Label();
352 // Skip primaries
353 if (label == ip) continue;
354 if (label > dau2 || label < dau1)
355 printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
356 itlast, label, dau1, dau2);
357 if (label == id) {
358 // secondary found
359 tr->SetDetectorId(ib-1);
360 Int_t nref = fTrackReferences->GetEntriesFast();
361 TClonesArray &lref = *fTrackReferences;
362 new(lref[nref]) AliTrackReference(*tr);
363 }
364 } // hits
365 } // branches
366 fTmpTreeTR->Fill();
367 fTrackReferences->Clear();
368 ifills++;
369 } // daughters
370 } // has hits
371 } // tracks
372
373 //
374 // Now loop again and write the primaries
375 //
376 it = nt - 1;
377 for (Int_t ip = 0; ip < np; ip++) {
378 Int_t labmax = -1;
379 while (labmax < ip && it > -1) {
380 fTreeTR->GetEntry(it--);
381 for (Int_t ib = 0; ib < 7; ib++) {
382 if (!trefs[ib]) continue;
383 Int_t nh = trefs[ib]->GetEntries();
384 //
385 // Loop over reference hits and find primary labels
386 for (Int_t ih = 0; ih < nh; ih++) {
387 AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
388 Int_t label = tr->Label();
389 if (label < np && label > labmax) {
390 labmax = label;
391 }
392
393 if (label == ip) {
394 tr->SetDetectorId(ib-1);
395 Int_t nref = fTrackReferences->GetEntriesFast();
396 TClonesArray &lref = *fTrackReferences;
397 new(lref[nref]) AliTrackReference(*tr);
398 }
399 } // hits
400 } // branches
401 } // entries
402 it++;
403 fTmpTreeTR->Fill();
404 fTrackReferences->Clear();
405 ifills++;
406 } // tracks
407 // Check
408
409
410 // Clean-up
411 delete fTreeTR; fTreeTR = 0;
412
413 for (Int_t ib = 0; ib < 7; ib++) {
414 if (trefs[ib]) {
415 trefs[ib]->Clear();
416 delete trefs[ib];
417 trefs[ib] = 0;
418 }
419 }
420
421 if (ifills != fStack->GetNtrack())
422 printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
423 ifills, fStack->GetNtrack());
424
425 fTmpTreeTR->Write();
426 fTreeTR = fTmpTreeTR;
427}
428
429AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
430{
431 // Get MC Particle i
432 AliMCParticle *mcParticle;
433 TParticle *particle;
434
435 // Out of range check
436 if (i < 0 || i >= fNparticles) {
437 AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
438 mcParticle = 0;
439 return (mcParticle);
440 }
441
442
443 particle = fStack->Particle(i);
444 //
445 // First check of the MC Particle has been already cached
446 if(!fMCParticleMap->At(i)) {
447 Int_t nentries = fMCParticles->GetEntriesFast();
448 new ((*fMCParticles)[nentries]) AliMCParticle(particle);
449 mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
450 fMCParticleMap->AddAt(mcParticle, i);
451 } else {
452 mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
453 }
454
455 return mcParticle;
456}
457
458inline AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
459
460ClassImp(AliMCEvent)