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