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