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