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