Corrected Doxygen warnings:
[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 /// 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   delete fMuonTrackRef;
121   if(fIsLoadConstructor){
122     fRunLoader->UnloadKinematics();
123     fRunLoader->UnloadTrackRefs();
124     fRunLoader->UnloadTracks();
125     delete fMUONData;
126   }
127   delete fRunLoader;
128 }
129
130 //_____________________________________________________________________________
131 void AliMUONRecoCheck::MakeTrackRef()
132 {
133 /// Make reconstructible tracks
134
135   AliTrackReference *trackReference;
136   AliMUONTrack *muonTrack;
137   AliMUONTrackParam *trackParam;
138   AliMUONHitForRec *hitForRec;
139   Float_t x, y, z, pX, pY, pZ, pYZ;
140   Int_t track, trackSave;
141   TParticle *particle;   
142   Float_t bendingSlope = 0;
143   Float_t  nonBendingSlope = 0;
144   Float_t inverseBendingMomentum = 0;
145   
146   TTree* treeTR = fRunLoader->TreeTR();
147   if (treeTR == NULL) return;
148
149   TBranch* branch = treeTR->GetBranch("MUON");
150   if (branch == NULL) return;
151
152   TClonesArray* trackRefs = 0;
153   branch->SetAddress(&trackRefs);
154   branch->SetAutoDelete(kTRUE);  
155
156   Int_t nTrackRef = (Int_t)branch->GetEntries();
157  
158   track = trackSave = -999;
159   Bool_t isNewTrack;
160   Int_t iHitMin, iChamber, detElemId;
161
162   trackParam = new AliMUONTrackParam();
163   hitForRec = new AliMUONHitForRec();
164
165   Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
166   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; iTrackRef++) {
167
168     branch->GetEntry(iTrackRef);
169     
170     iHitMin = 0;
171     isNewTrack = kTRUE;
172     
173     if (!trackRefs->GetEntries()) continue; 
174
175     while (isNewTrack) {
176
177       muonTrack = new AliMUONTrack();
178       
179       for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
180       
181         trackReference = (AliTrackReference*)trackRefs->At(iHit);
182         x = trackReference->X();
183         y = trackReference->Y();
184         z = trackReference->Z();
185         pX = trackReference->Px();
186         pY = trackReference->Py();
187         pZ = trackReference->Pz();
188
189         track = trackReference->GetTrack();
190
191         if(track >= max){
192           AliWarningStream()
193             << "Track ID " << track 
194             << " larger than max number of particles " << max << endl;
195           isNewTrack = kFALSE;
196           break;
197         }
198         if (track != trackSave && iHit != 0) {
199           iHitMin = iHit;
200           trackSave = track;
201           break;
202         }
203
204         // track parameters at hit
205         trackParam->SetBendingCoor(y);
206         trackParam->SetNonBendingCoor(x);
207         trackParam->SetZ(z);
208
209         if (TMath::Abs(pZ) > 0) {
210           bendingSlope = pY/pZ;
211           nonBendingSlope = pX/pZ;
212         }
213         pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
214         if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
215
216         trackParam->SetBendingSlope(bendingSlope);
217         trackParam->SetNonBendingSlope(nonBendingSlope);
218         trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
219
220         hitForRec->SetBendingCoor(y);
221         hitForRec->SetNonBendingCoor(x);
222         hitForRec->SetZ(z);
223         hitForRec->SetBendingReso2(0.0); 
224         hitForRec->SetNonBendingReso2(0.0);
225         detElemId = hitForRec->GetDetElemId();
226         if (detElemId) iChamber = detElemId / 100 - 1; 
227         else iChamber = AliMUONConstants::ChamberNumber(z);
228         hitForRec->SetChamberNumber(iChamber);
229
230         muonTrack->AddTrackParamAtHit(trackParam,0);
231         muonTrack->AddHitForRecAtHit(hitForRec);
232         muonTrack->SetTrackID(track);
233
234         trackSave = track;
235         if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
236       }
237
238       // track parameters at vertex 
239       particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
240
241       if (particle) {
242
243         x = particle->Vx();
244         y = particle->Vy();
245         z = particle->Vz();
246         pX = particle->Px();
247         pY = particle->Py();
248         pZ = particle->Pz();
249
250         trackParam->SetBendingCoor(y);
251         trackParam->SetNonBendingCoor(x);
252         trackParam->SetZ(z);
253         if (TMath::Abs(pZ) > 0) {
254           bendingSlope = pY/pZ;
255           nonBendingSlope = pX/pZ;
256         }
257         pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
258         if (pYZ >0) inverseBendingMomentum = 1/pYZ;       
259         trackParam->SetBendingSlope(bendingSlope);
260         trackParam->SetNonBendingSlope(nonBendingSlope);
261         trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
262       
263         muonTrack->SetTrackParamAtVertex(trackParam);
264       }
265
266       AddMuonTrackReference(muonTrack);
267       delete muonTrack;
268       muonTrack = NULL;
269       
270     } // end while isNewTrack
271   }
272   delete trackRefs;
273   CleanMuonTrackRef();
274   
275   ReconstructibleTracks();
276
277   delete trackParam;
278   delete hitForRec;
279
280 }
281
282 //____________________________________________________________________________
283 TClonesArray* AliMUONRecoCheck::GetTrackReco()
284 {
285 /// Return TClonesArray of reconstructed tracks
286
287   fMUONData->ResetRecTracks();
288   fMUONData->SetTreeAddress("RT");
289   fTrackReco = fMUONData->RecTracks(); 
290   fMUONData->GetRecTracks();
291   fRecoTracks = fTrackReco->GetEntriesFast();
292   return fTrackReco;
293 }
294 //_____________________________________________________________________________
295 void AliMUONRecoCheck::PrintEvent() const
296 {
297 /// Debug facility
298
299   AliMUONTrack *track;
300   AliMUONHitForRec *hitForRec;
301   TClonesArray *  hitForRecAtHit = 0;
302   Float_t xRec,yRec,zRec;
303
304   Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
305   
306   printf(" ******************************* \n");
307   printf(" nb of tracks %d \n",nTrackRef);
308
309   for (Int_t index = 0; index < nTrackRef; index++) {
310     track = (AliMUONTrack*)fMuonTrackRef->At(index);
311     hitForRecAtHit = track->GetHitForRecAtHit();
312     Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
313     printf(" track number %d \n",index);
314     for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
315       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
316       xRec  = hitForRec->GetNonBendingCoor();
317       yRec  = hitForRec->GetBendingCoor();
318       zRec  = hitForRec->GetZ();
319       printf("  x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
320     }
321   }
322 }
323
324 //_____________________________________________________________________________
325 void AliMUONRecoCheck::ResetTracks() const
326 {
327 /// Reset tracks
328
329   if (fMuonTrackRef) fMuonTrackRef->Delete();
330 }
331 //_____________________________________________________________________________
332 void AliMUONRecoCheck::CleanMuonTrackRef()
333 {
334 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
335 /// each chamber (one when particle is entering + one when particle is leaving 
336 /// the sensitive volume) 
337
338   Float_t maxGasGap = 1.; // cm 
339   AliMUONTrack *track, *trackNew;
340   AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
341   AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
342   TClonesArray *  hitForRecAtHit = 0;
343   TClonesArray *  trackParamAtHit = 0;
344   Float_t xRec,yRec,zRec;
345   Float_t xRec1,yRec1,zRec1;
346   Float_t xRec2,yRec2,zRec2;
347   Float_t bendingSlope,nonBendingSlope,bendingMomentum;
348   Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
349   Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
350   TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
351   Int_t iHit1;
352   Int_t iChamber = 0, detElemId = 0;
353   Int_t nRec = 0;
354   Int_t nTrackHits = 0;
355
356   hitForRec = new AliMUONHitForRec();
357   trackParam = new AliMUONTrackParam();
358
359   Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
360   for (Int_t index = 0; index < nTrackRef; index++) {
361     track = (AliMUONTrack*)fMuonTrackRef->At(index);
362     hitForRecAtHit = track->GetHitForRecAtHit();
363     trackParamAtHit = track->GetTrackParamAtHit();
364     trackParamAtVertex = track->GetTrackParamAtVertex();
365     nTrackHits = hitForRecAtHit->GetEntriesFast();
366     trackNew = new AliMUONTrack();
367     iHit1 = 0;
368     while (iHit1 < nTrackHits) {
369       hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1); 
370       trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1); 
371       xRec1  = hitForRec1->GetNonBendingCoor();
372       yRec1  = hitForRec1->GetBendingCoor();
373       zRec1  = hitForRec1->GetZ();      
374       xRec   = xRec1;
375       yRec   = yRec1;
376       zRec   = zRec1;
377       bendingSlope1 = trackParam1->GetBendingSlope();
378       nonBendingSlope1 = trackParam1->GetNonBendingSlope();
379       bendingMomentum1 = 0;
380       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
381         bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
382       bendingSlope = bendingSlope1;
383       nonBendingSlope = nonBendingSlope1;
384       bendingMomentum = bendingMomentum1;
385       nRec = 1;  
386       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
387         hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2); 
388         trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2); 
389         xRec2  = hitForRec2->GetNonBendingCoor();
390         yRec2  = hitForRec2->GetBendingCoor();
391         zRec2  = hitForRec2->GetZ();      
392         bendingSlope2 = trackParam2->GetBendingSlope();
393         nonBendingSlope2 = trackParam2->GetNonBendingSlope();
394         bendingMomentum2 = 0;
395         if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
396           bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
397         
398         if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
399
400           nRec++;
401           xRec += xRec2;
402           yRec += yRec2;
403           zRec += zRec2;
404           bendingSlope += bendingSlope2;
405           nonBendingSlope += nonBendingSlope2;
406           bendingMomentum += bendingMomentum2;
407           iHit1 = iHit2;
408         }
409         
410       } // end iHit2
411       xRec /= (Float_t)nRec;
412       yRec /= (Float_t)nRec;
413       zRec /= (Float_t)nRec;
414       bendingSlope /= (Float_t)nRec;
415       nonBendingSlope /= (Float_t)nRec;
416       bendingMomentum /= (Float_t)nRec;
417       
418       hitForRec->SetNonBendingCoor(xRec);
419       hitForRec->SetBendingCoor(yRec);
420       hitForRec->SetZ(zRec);
421       detElemId = hitForRec->GetDetElemId();
422       if (detElemId) iChamber = detElemId / 100 - 1;
423       else iChamber = AliMUONConstants::ChamberNumber(zRec);
424       hitForRec->SetChamberNumber(iChamber);
425       hitForRec->SetBendingReso2(0.0); 
426       hitForRec->SetNonBendingReso2(0.0); 
427       trackParam->SetNonBendingCoor(xRec);
428       trackParam->SetBendingCoor(yRec);
429       trackParam->SetZ(zRec);
430       trackParam->SetNonBendingSlope(nonBendingSlope);
431       trackParam->SetBendingSlope(bendingSlope);
432       if (TMath::Abs(bendingMomentum) > 0)
433         trackParam->SetInverseBendingMomentum(1./bendingMomentum);
434
435       trackNew->AddHitForRecAtHit(hitForRec);
436       trackNew->AddTrackParamAtHit(trackParam,0);
437       
438       iHit1++;
439     } // end iHit1
440
441     trackNew->SetTrackID(track->GetTrackID());
442     trackNew->SetTrackParamAtVertex(trackParamAtVertex);
443     {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}   
444     delete trackNew;
445     trackNew = NULL;
446     
447   } // end trackRef
448
449   fMuonTrackRef->Delete();
450   nTrackRef = newMuonTrackRef->GetEntriesFast();
451   for (Int_t index = 0; index < nTrackRef; index++) {
452     track = (AliMUONTrack*)newMuonTrackRef->At(index);
453     AddMuonTrackReference(track);
454   }
455   
456   delete hitForRec;
457   delete trackParam;
458   newMuonTrackRef->Delete();
459   delete newMuonTrackRef;
460   
461 }
462
463 //_____________________________________________________________________________
464 void AliMUONRecoCheck::ReconstructibleTracks()
465 {
466 /// Calculate the number of reconstructible tracks
467   
468   AliMUONTrack* track;
469   TClonesArray* hitForRecAtHit = NULL;
470   AliMUONHitForRec* hitForRec;
471   Float_t zRec;
472   Int_t nTrackRef, nTrackHits;
473   Int_t isChamberInTrack[10];
474   Int_t iChamber = 0;
475   Bool_t isTrackOK = kTRUE;
476
477   fReconstructibleTracks = 0;
478
479   nTrackRef = fMuonTrackRef->GetEntriesFast();
480   for (Int_t index = 0; index < nTrackRef; index++) {
481     track = (AliMUONTrack*)fMuonTrackRef->At(index);
482     hitForRecAtHit = track->GetHitForRecAtHit();
483     nTrackHits = hitForRecAtHit->GetEntriesFast();
484     for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
485  
486     for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
487       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
488       zRec  = hitForRec->GetZ();
489       iChamber = hitForRec->GetChamberNumber();
490       if (iChamber < 0 || iChamber > 10) continue;
491       isChamberInTrack[iChamber] = 1;
492     } 
493     // track is reconstructible if the particle is depositing a hit
494     // in the following chamber combinations:
495
496     isTrackOK = kTRUE;
497     if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
498     if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
499     if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
500     Int_t nHitsInLastStations=0;
501     for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
502       if (isChamberInTrack[ch]) nHitsInLastStations++; 
503     if(nHitsInLastStations < 3) isTrackOK = kFALSE;
504
505     if (isTrackOK) fReconstructibleTracks++;
506     if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
507   }
508   fMuonTrackRef->Compress();
509 }
510