Bug fix. Removed delete statement
[u/mrichter/AliRoot.git] / STEER / AliMCEvent.cxx
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 <TRefArray.h>
30
31 #include "AliLog.h"
32 #include "AliMCEvent.h"
33 #include "AliStack.h"
34 #include "AliTrackReference.h"
35 #include "AliHeader.h"
36 #include "AliGenEventHeader.h"
37
38
39 AliMCEvent::AliMCEvent():
40     AliVEvent(),
41     fStack(0),
42     fMCParticles(new TClonesArray("AliMCParticle",1000)),
43     fMCParticleMap(0),
44     fHeader(new AliHeader()),
45     fTRBuffer(0),
46     fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
47     fTreeTR(0),
48     fTmpTreeTR(0),
49     fTmpFileTR(0),
50     fNprimaries(-1),
51     fNparticles(-1)
52 {
53     // Default constructor
54 }
55
56 AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
57     AliVEvent(mcEvnt),
58     fStack(0),
59     fMCParticles(0),
60     fMCParticleMap(0),
61     fHeader(new AliHeader()),
62     fTRBuffer(0),
63     fTrackReferences(0),
64     fTreeTR(0),
65     fTmpTreeTR(0),
66     fTmpFileTR(0),
67     fNprimaries(-1),
68     fNparticles(-1)
69
70 // Copy constructor
71 }
72
73
74 AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
75 {
76     // assignment operator
77     if (this!=&mcEvnt) { 
78         AliVEvent::operator=(mcEvnt); 
79     }
80   
81     return *this; 
82 }
83
84 void AliMCEvent::ConnectTreeE (TTree* tree)
85 {
86     // Connect the event header tree
87     tree->SetBranchAddress("Header", &fHeader);
88 }
89
90 void 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();
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));
105  
106     // This is a cache for the TParticles converted to MCParticles on user request
107     if (fMCParticleMap) {
108         fMCParticleMap->Clear();
109         fMCParticles->Delete();
110         if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
111     }
112     else
113         fMCParticleMap = new TRefArray(fNparticles);
114 }
115
116 void 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 
130         fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
131     }
132 }
133
134 Int_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));
146         trefs    = fTRBuffer;
147         return trefs->GetEntries();
148     } else {
149         trefs = 0;
150         return -1;
151     }
152 }
153
154
155 void AliMCEvent::Clean()
156 {
157     // Clean-up before new trees are connected
158 //    if (fHeader) {
159 //      delete fHeader;
160 //      fHeader = 0;
161 //    }
162
163     delete fStack; fStack = 0;
164
165     // Clear TR
166     if (fTRBuffer) {
167         fTRBuffer->Delete();
168         delete fTRBuffer;
169         fTRBuffer = 0;
170     }
171 }
172
173 void AliMCEvent::FinishEvent()
174 {
175   // Clean-up after event
176   //    
177     fStack->Reset(0);
178     fMCParticles->Delete();
179     fMCParticleMap->Clear();
180     if (fTRBuffer)
181       fTRBuffer->Delete();
182     fTrackReferences->Delete();
183     fNparticles = -1;
184     fNprimaries = -1;    
185     fStack      =  0;
186     
187 }
188
189
190
191 void 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     
206     Int_t nh = fTRBuffer->GetEntries();
207     
208     
209     if (search) {
210         while(nh <= search && i < fNparticles - 1) {
211             i++;
212             fTreeTR->GetEntry(fStack->TreeKEntry(i));
213             nh =  fTRBuffer->GetEntries();
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++) {
233         AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
234         TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
235         m->Draw();
236         m->SetMarkerSize(0.4);
237         
238     }
239 }
240
241
242 void 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");
252     if (!fTRBuffer)  fTRBuffer = new TClonesArray("AliTrackReference", 100);
253     fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
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);
377                             Int_t nref =  fTRBuffer->GetEntriesFast();
378                             TClonesArray &lref = *fTRBuffer;
379                             new(lref[nref]) AliTrackReference(*tr);
380                         }
381                     } // hits
382                 } // branches
383                 fTmpTreeTR->Fill();
384                 fTRBuffer->Delete();
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);
412                         Int_t nref = fTRBuffer->GetEntriesFast();
413                         TClonesArray &lref = *fTRBuffer;
414                         new(lref[nref]) AliTrackReference(*tr);
415                     }
416                 } // hits
417             } // branches
418         } // entries
419         it++;
420         fTmpTreeTR->Fill();
421         fTRBuffer->Delete();
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
446 AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
447 {
448     // Get MC Particle i
449     AliMCParticle *mcParticle = 0;
450     TParticle     *particle   = 0;
451     TClonesArray  *trefs      = 0;
452     Int_t          ntref      = 0;
453     TRefArray     *rarray     = 0;
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
462     //
463     // First check of the MC Particle has been already cached
464     if(!fMCParticleMap->At(i)) {
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
483         Int_t nentries = fMCParticles->GetEntriesFast();
484         new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
485         mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
486         fMCParticleMap->AddAt(mcParticle, i);
487     } else {
488         mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
489     }
490     return mcParticle;
491 }
492
493  AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
494
495 ClassImp(AliMCEvent)