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