Removing small memory leaks (Gines)
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
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 AliMUONRecoCheck
19 /// \brief Utility class to check reconstruction
20 /// Reconstructed tracks are compared to reference tracks. 
21 /// The reference tracks are built from AliTrackReference for the
22 /// hit in chamber (0..9) and from kinematics for the vertex parameters.     
23
24 #include "AliMUON.h"
25 #include "AliMUONRecoCheck.h"
26 #include "AliMUONTrack.h"
27 #include "AliMUONData.h"
28 #include "AliMUONConstants.h"
29
30 #include "AliLoader.h" 
31 #include "AliRunLoader.h" 
32 #include "AliHeader.h"
33 #include "AliStack.h"
34 #include "AliTrackReference.h"
35 #include "AliLog.h" 
36
37 #include <TParticle.h>
38
39 /// \cond CLASSIMP
40 ClassImp(AliMUONRecoCheck)
41 /// \endcond
42
43 //_____________________________________________________________________________
44   AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
45   : TObject(),
46   fRunLoader(0x0),
47   fMUONData(0x0),
48   fMuonTrackRef(0x0),
49   fTrackReco(0x0),
50   fReconstructibleTracks(0),
51   fRecoTracks(0),
52   fIsLoadConstructor(kTRUE)
53 {
54 /// Constructor using "galice.root",
55 /// takes care of loading/unloading Kinematics, TrackRefs and Tracks
56
57   fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
58
59   // run loader
60   fRunLoader = AliRunLoader::Open(chLoader);
61   if (!fRunLoader) {
62     AliError(Form("no run loader found " ));
63     return;
64   }
65
66  // initialize loader    
67   AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
68
69   // container
70   fMUONData  = new AliMUONData(loader,"MUON","MUON");
71   if (!fMUONData) {
72     AliError(Form("no MUONData found " ));
73     return;
74   }
75
76   fRunLoader->LoadKinematics("READ");    
77   fRunLoader->LoadTrackRefs("READ");     
78   loader->LoadTracks("READ");
79
80 }
81
82 //_____________________________________________________________________________
83   AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader *runloader, AliMUONData *muondata)
84   : TObject(),
85   fRunLoader(0x0),
86   fMUONData(0x0),
87   fMuonTrackRef(0x0),
88   fTrackReco(0x0),
89   fReconstructibleTracks(0),
90   fRecoTracks(0),
91   fIsLoadConstructor(kFALSE)
92 {
93 /// Constructor from AliRunLoader and AliMUONData,
94 /// does not load/unload Kinematics, TrackRefs and Tracks internally
95 /// it should be done in the execution program (e.g. see MUONRecoCheck.C)
96
97   fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
98
99   // run loader
100   fRunLoader = runloader;
101   if (!fRunLoader) {
102     AliError(Form("no run loader found " ));
103     return;
104   }
105
106   // container
107   fMUONData  = muondata;
108   if (!fMUONData) {
109     AliError(Form("no MUONData found " ));
110     return;
111   }
112
113 }
114
115
116 //_____________________________________________________________________________
117 AliMUONRecoCheck::~AliMUONRecoCheck()
118 {
119 /// Destructor
120
121   delete fMuonTrackRef;
122   if(fIsLoadConstructor){
123     fRunLoader->UnloadKinematics();
124     fRunLoader->UnloadTrackRefs();
125     fRunLoader->UnloadTracks();
126     delete fMUONData;
127   }
128   delete fRunLoader;
129 }
130
131 //_____________________________________________________________________________
132 void AliMUONRecoCheck::MakeTrackRef()
133 {
134 /// Make reconstructible tracks
135
136   AliTrackReference *trackReference;
137   AliMUONTrack *muonTrack;
138   AliMUONTrackParam *trackParam;
139   AliMUONHitForRec *hitForRec;
140   Float_t x, y, z, pX, pY, pZ, pYZ;
141   Int_t track, trackSave;
142   TParticle *particle;   
143   Float_t bendingSlope = 0;
144   Float_t  nonBendingSlope = 0;
145   Float_t inverseBendingMomentum = 0;
146   
147   TTree* treeTR = fRunLoader->TreeTR();
148   if (treeTR == NULL) return;
149
150   TBranch* branch = treeTR->GetBranch("MUON");
151   if (branch == NULL) return;
152
153   TClonesArray* trackRefs = 0;
154   branch->SetAddress(&trackRefs);
155   branch->SetAutoDelete(kTRUE);  
156
157   Int_t nTrackRef = (Int_t)branch->GetEntries();
158  
159   track = trackSave = -999;
160   Bool_t isNewTrack;
161   Int_t iHitMin, iChamber, detElemId;
162
163   trackParam = new AliMUONTrackParam();
164   hitForRec = new AliMUONHitForRec();
165
166   Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
167   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; iTrackRef++) {
168
169     branch->GetEntry(iTrackRef);
170     
171     iHitMin = 0;
172     isNewTrack = kTRUE;
173     
174     if (!trackRefs->GetEntries()) continue; 
175
176     while (isNewTrack) {
177
178       muonTrack = new AliMUONTrack();
179       
180       for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
181       
182         trackReference = (AliTrackReference*)trackRefs->At(iHit);
183         x = trackReference->X();
184         y = trackReference->Y();
185         z = trackReference->Z();
186         pX = trackReference->Px();
187         pY = trackReference->Py();
188         pZ = trackReference->Pz();
189
190         track = trackReference->GetTrack();
191
192         if(track >= max){
193           AliWarningStream()
194             << "Track ID " << track 
195             << " larger than max number of particles " << max << endl;
196           isNewTrack = kFALSE;
197           break;
198         }
199         if (track != trackSave && iHit != 0) {
200           iHitMin = iHit;
201           trackSave = track;
202           break;
203         }
204
205         // track parameters at hit
206         trackParam->SetBendingCoor(y);
207         trackParam->SetNonBendingCoor(x);
208         trackParam->SetZ(z);
209
210         if (TMath::Abs(pZ) > 0) {
211           bendingSlope = pY/pZ;
212           nonBendingSlope = pX/pZ;
213         }
214         pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
215         if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
216
217         trackParam->SetBendingSlope(bendingSlope);
218         trackParam->SetNonBendingSlope(nonBendingSlope);
219         trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
220
221         hitForRec->SetBendingCoor(y);
222         hitForRec->SetNonBendingCoor(x);
223         hitForRec->SetZ(z);
224         hitForRec->SetBendingReso2(0.0); 
225         hitForRec->SetNonBendingReso2(0.0);
226         detElemId = hitForRec->GetDetElemId();
227         if (detElemId) iChamber = detElemId / 100 - 1; 
228         else iChamber = AliMUONConstants::ChamberNumber(z);
229         hitForRec->SetChamberNumber(iChamber);
230
231         muonTrack->AddTrackParamAtHit(trackParam,0);
232         muonTrack->AddHitForRecAtHit(hitForRec);
233         muonTrack->SetTrackID(track);
234
235         trackSave = track;
236         if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
237       }
238
239       // track parameters at vertex 
240       particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
241
242       if (particle) {
243
244         x = particle->Vx();
245         y = particle->Vy();
246         z = particle->Vz();
247         pX = particle->Px();
248         pY = particle->Py();
249         pZ = particle->Pz();
250
251         trackParam->SetBendingCoor(y);
252         trackParam->SetNonBendingCoor(x);
253         trackParam->SetZ(z);
254         if (TMath::Abs(pZ) > 0) {
255           bendingSlope = pY/pZ;
256           nonBendingSlope = pX/pZ;
257         }
258         pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
259         if (pYZ >0) inverseBendingMomentum = 1/pYZ;       
260         trackParam->SetBendingSlope(bendingSlope);
261         trackParam->SetNonBendingSlope(nonBendingSlope);
262         trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
263       
264         muonTrack->SetTrackParamAtVertex(trackParam);
265       }
266
267       AddMuonTrackReference(muonTrack);
268       delete muonTrack;
269       muonTrack = NULL;
270       
271     } // end while isNewTrack
272   }
273   
274   CleanMuonTrackRef();
275   
276   ReconstructibleTracks();
277
278   delete trackParam;
279   delete hitForRec;
280
281 }
282
283 //____________________________________________________________________________
284 TClonesArray* AliMUONRecoCheck::GetTrackReco()
285 {
286 /// Return TClonesArray of reconstructed tracks
287
288   fMUONData->ResetRecTracks();
289   fMUONData->SetTreeAddress("RT");
290   fTrackReco = fMUONData->RecTracks(); 
291   fMUONData->GetRecTracks();
292   fRecoTracks = fTrackReco->GetEntriesFast();
293   return fTrackReco;
294 }
295 //_____________________________________________________________________________
296 void AliMUONRecoCheck::PrintEvent() const
297 {
298 /// Debug facility
299
300   AliMUONTrack *track;
301   AliMUONHitForRec *hitForRec;
302   TClonesArray *  hitForRecAtHit = 0;
303   Float_t xRec,yRec,zRec;
304
305   Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
306   
307   printf(" ******************************* \n");
308   printf(" nb of tracks %d \n",nTrackRef);
309
310   for (Int_t index = 0; index < nTrackRef; index++) {
311     track = (AliMUONTrack*)fMuonTrackRef->At(index);
312     hitForRecAtHit = track->GetHitForRecAtHit();
313     Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
314     printf(" track number %d \n",index);
315     for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
316       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
317       xRec  = hitForRec->GetNonBendingCoor();
318       yRec  = hitForRec->GetBendingCoor();
319       zRec  = hitForRec->GetZ();
320       printf("  x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
321     }
322   }
323 }
324
325 //_____________________________________________________________________________
326 void AliMUONRecoCheck::ResetTracks() const
327 {
328 /// Reset tracks
329
330   if (fMuonTrackRef) fMuonTrackRef->Clear();
331 }
332 //_____________________________________________________________________________
333 void AliMUONRecoCheck::CleanMuonTrackRef()
334 {
335 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
336 /// each chamber (one when particle is entering + one when particle is leaving 
337 /// the sensitive volume) 
338
339   Float_t maxGasGap = 1.; // cm 
340   AliMUONTrack *track, *trackNew;
341   AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
342   AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
343   TClonesArray *  hitForRecAtHit = 0;
344   TClonesArray *  trackParamAtHit = 0;
345   Float_t xRec,yRec,zRec;
346   Float_t xRec1,yRec1,zRec1;
347   Float_t xRec2,yRec2,zRec2;
348   Float_t bendingSlope,nonBendingSlope,bendingMomentum;
349   Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
350   Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
351   TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
352   Int_t iHit1;
353   Int_t iChamber = 0, detElemId = 0;
354   Int_t nRec = 0;
355   Int_t nTrackHits = 0;
356
357   hitForRec = new AliMUONHitForRec();
358   trackParam = new AliMUONTrackParam();
359
360   Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
361   for (Int_t index = 0; index < nTrackRef; index++) {
362     track = (AliMUONTrack*)fMuonTrackRef->At(index);
363     hitForRecAtHit = track->GetHitForRecAtHit();
364     trackParamAtHit = track->GetTrackParamAtHit();
365     trackParamAtVertex = track->GetTrackParamAtVertex();
366     nTrackHits = hitForRecAtHit->GetEntriesFast();
367     trackNew = new AliMUONTrack();
368     iHit1 = 0;
369     while (iHit1 < nTrackHits) {
370       hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1); 
371       trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1); 
372       xRec1  = hitForRec1->GetNonBendingCoor();
373       yRec1  = hitForRec1->GetBendingCoor();
374       zRec1  = hitForRec1->GetZ();      
375       xRec   = xRec1;
376       yRec   = yRec1;
377       zRec   = zRec1;
378       bendingSlope1 = trackParam1->GetBendingSlope();
379       nonBendingSlope1 = trackParam1->GetNonBendingSlope();
380       bendingMomentum1 = 0;
381       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
382         bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
383       bendingSlope = bendingSlope1;
384       nonBendingSlope = nonBendingSlope1;
385       bendingMomentum = bendingMomentum1;
386       nRec = 1;  
387       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
388         hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2); 
389         trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2); 
390         xRec2  = hitForRec2->GetNonBendingCoor();
391         yRec2  = hitForRec2->GetBendingCoor();
392         zRec2  = hitForRec2->GetZ();      
393         bendingSlope2 = trackParam2->GetBendingSlope();
394         nonBendingSlope2 = trackParam2->GetNonBendingSlope();
395         bendingMomentum2 = 0;
396         if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
397           bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
398         
399         if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
400
401           nRec++;
402           xRec += xRec2;
403           yRec += yRec2;
404           zRec += zRec2;
405           bendingSlope += bendingSlope2;
406           nonBendingSlope += nonBendingSlope2;
407           bendingMomentum += bendingMomentum2;
408           iHit1 = iHit2;
409         }
410         
411       } // end iHit2
412       xRec /= (Float_t)nRec;
413       yRec /= (Float_t)nRec;
414       zRec /= (Float_t)nRec;
415       bendingSlope /= (Float_t)nRec;
416       nonBendingSlope /= (Float_t)nRec;
417       bendingMomentum /= (Float_t)nRec;
418       
419       hitForRec->SetNonBendingCoor(xRec);
420       hitForRec->SetBendingCoor(yRec);
421       hitForRec->SetZ(zRec);
422       detElemId = hitForRec->GetDetElemId();
423       if (detElemId) iChamber = detElemId / 100 - 1;
424       else iChamber = AliMUONConstants::ChamberNumber(zRec);
425       hitForRec->SetChamberNumber(iChamber);
426       hitForRec->SetBendingReso2(0.0); 
427       hitForRec->SetNonBendingReso2(0.0); 
428       trackParam->SetNonBendingCoor(xRec);
429       trackParam->SetBendingCoor(yRec);
430       trackParam->SetZ(zRec);
431       trackParam->SetNonBendingSlope(nonBendingSlope);
432       trackParam->SetBendingSlope(bendingSlope);
433       if (TMath::Abs(bendingMomentum) > 0)
434         trackParam->SetInverseBendingMomentum(1./bendingMomentum);
435
436       trackNew->AddHitForRecAtHit(hitForRec);
437       trackNew->AddTrackParamAtHit(trackParam,0);
438       
439       iHit1++;
440     } // end iHit1
441
442     trackNew->SetTrackID(track->GetTrackID());
443     trackNew->SetTrackParamAtVertex(trackParamAtVertex);
444     {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}   
445     delete trackNew;
446     trackNew = NULL;
447     
448   } // end trackRef
449
450   fMuonTrackRef->Clear();
451   nTrackRef = newMuonTrackRef->GetEntriesFast();
452   for (Int_t index = 0; index < nTrackRef; index++) {
453     track = (AliMUONTrack*)newMuonTrackRef->At(index);
454     AddMuonTrackReference(track);
455   }
456   
457   delete hitForRec;
458   delete trackParam;
459   newMuonTrackRef->Delete();
460   delete newMuonTrackRef;
461   
462 }
463
464 //_____________________________________________________________________________
465 void AliMUONRecoCheck::ReconstructibleTracks()
466 {
467 /// Calculate the number of reconstructible tracks
468   
469   AliMUONTrack* track;
470   TClonesArray* hitForRecAtHit = NULL;
471   AliMUONHitForRec* hitForRec;
472   Float_t zRec;
473   Int_t nTrackRef, nTrackHits;
474   Int_t isChamberInTrack[10];
475   Int_t iChamber = 0;
476   Bool_t isTrackOK = kTRUE;
477
478   fReconstructibleTracks = 0;
479
480   nTrackRef = fMuonTrackRef->GetEntriesFast();
481   for (Int_t index = 0; index < nTrackRef; index++) {
482     track = (AliMUONTrack*)fMuonTrackRef->At(index);
483     hitForRecAtHit = track->GetHitForRecAtHit();
484     nTrackHits = hitForRecAtHit->GetEntriesFast();
485     for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
486  
487     for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
488       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
489       zRec  = hitForRec->GetZ();
490       iChamber = hitForRec->GetChamberNumber();
491       if (iChamber < 0 || iChamber > 10) continue;
492       isChamberInTrack[iChamber] = 1;
493     } 
494     // track is reconstructible if the particle is depositing a hit
495     // in the following chamber combinations:
496
497     isTrackOK = kTRUE;
498     if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
499     if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
500     if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
501     Int_t nHitsInLastStations=0;
502     for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
503       if (isChamberInTrack[ch]) nHitsInLastStations++; 
504     if(nHitsInLastStations < 3) isTrackOK = kFALSE;
505
506     if (isTrackOK) fReconstructibleTracks++;
507     if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
508   }
509   fMuonTrackRef->Compress();
510 }
511