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