0d7d8ab40dacbb9bc506bc4ed81b435c4a222be2
[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 //-----------------------------------------------------------------------------
19 /// \class AliMUONRecoCheck
20 /// Utility class to check reconstruction
21 /// Reconstructed tracks are compared to reference tracks. 
22 /// The reference tracks are built from AliTrackReference for the
23 /// hit in chamber (0..9) and from kinematics for the vertex parameters.     
24 //-----------------------------------------------------------------------------
25
26 #include "AliMUONRecoCheck.h"
27 #include "AliMUONHitForRec.h"
28 #include "AliMUONTrack.h"
29 #include "AliMUONConstants.h"
30 #include "AliMUONMCDataInterface.h"
31 #include "AliMUONDataInterface.h"
32 #include "AliStack.h"
33 #include "AliTrackReference.h"
34 #include "AliLog.h" 
35 #include "AliMUONTrackStoreV1.h"
36 #include <Riostream.h>
37 #include <TParticle.h>
38 #include <TParticlePDG.h>
39
40 /// \cond CLASSIMP
41 ClassImp(AliMUONRecoCheck)
42 /// \endcond
43
44 //_____________________________________________________________________________
45 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
46 : TObject(),
47 fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
48 fDataInterface(new AliMUONDataInterface(chLoader))
49 {
50   /// Normal ctor
51 }
52
53 //_____________________________________________________________________________
54 AliMUONRecoCheck::~AliMUONRecoCheck()
55 {
56   /// Destructor
57   delete fMCDataInterface;
58   delete fDataInterface;
59 }
60
61 //_____________________________________________________________________________
62 AliMUONVTrackStore* 
63 AliMUONRecoCheck::ReconstructedTracks(Int_t event)
64 {
65   /// Return the reconstructed track store for a given event
66   return fDataInterface->TrackStore(event);
67 }
68
69 //_____________________________________________________________________________
70 AliMUONVTrackStore* 
71 AliMUONRecoCheck::TrackRefs(Int_t event)
72 {
73   /// Return a track store containing the track references (converted into 
74   /// MUONTrack objects) for a given event
75   return MakeTrackRefs(event);
76 }
77
78 //_____________________________________________________________________________
79 AliMUONVTrackStore* 
80 AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
81 {
82   /// Return a track store containing the reconstructible tracks for a given event
83   AliMUONVTrackStore* tmp = MakeTrackRefs(event);
84   AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
85   delete tmp;
86   return reconstructible;
87 }
88
89 //_____________________________________________________________________________
90 AliMUONVTrackStore*
91 AliMUONRecoCheck::MakeTrackRefs(Int_t event)
92 {
93   /// Make reconstructible tracks
94     
95   AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
96   
97   Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
98   
99   Int_t trackSave(-999);
100   Int_t iHitMin(0);
101   
102   Bool_t isNewTrack;
103     
104   AliStack* stack = fMCDataInterface->Stack(event);
105   Int_t max = stack->GetNtrack();
106   
107   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) 
108   {
109     TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
110     
111     iHitMin = 0;
112     isNewTrack = kTRUE;
113     
114     if (!trackRefs->GetEntries()) continue; 
115     
116     while (isNewTrack) {
117       
118       AliMUONTrack muonTrack;
119       
120       for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit) 
121       {        
122         AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
123         
124         Float_t x = trackReference->X();
125         Float_t y = trackReference->Y();
126         Float_t z = trackReference->Z();
127         Float_t pX = trackReference->Px();
128         Float_t pY = trackReference->Py();
129         Float_t pZ = trackReference->Pz();
130         
131         Int_t track = trackReference->GetTrack();
132         
133         if ( track >= max )
134         {
135           AliWarningStream()
136           << "Track ID " << track 
137           << " larger than max number of particles " << max << endl;
138           isNewTrack = kFALSE;
139           break;
140         }
141         if (track != trackSave && iHit != 0) {
142           iHitMin = iHit;
143           trackSave = track;
144           break;
145         }
146         
147         Float_t bendingSlope = 0;
148         Float_t nonBendingSlope = 0;
149         Float_t inverseBendingMomentum = 0;
150         
151         AliMUONTrackParam trackParam;
152         
153         // track parameters at hit
154         trackParam.SetBendingCoor(y);
155         trackParam.SetNonBendingCoor(x);
156         trackParam.SetZ(z);
157         
158         if (TMath::Abs(pZ) > 0) 
159         {
160           bendingSlope = pY/pZ;
161           nonBendingSlope = pX/pZ;
162         }
163         Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
164         if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
165         
166         trackParam.SetBendingSlope(bendingSlope);
167         trackParam.SetNonBendingSlope(nonBendingSlope);
168         trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
169         
170         AliMUONHitForRec hitForRec;
171         
172         hitForRec.SetBendingCoor(y);
173         hitForRec.SetNonBendingCoor(x);
174         hitForRec.SetZ(z);
175         hitForRec.SetBendingReso2(0.0); 
176         hitForRec.SetNonBendingReso2(0.0);
177         Int_t detElemId = hitForRec.GetDetElemId();
178         Int_t iChamber;
179         if (detElemId) 
180         {
181           iChamber = detElemId / 100 - 1; 
182         }
183         else 
184         {
185           iChamber = AliMUONConstants::ChamberNumber(z);
186         }
187         hitForRec.SetChamberNumber(iChamber);
188         
189         muonTrack.AddTrackParamAtHit(&trackParam,0);
190         muonTrack.AddHitForRecAtHit(&hitForRec);
191         muonTrack.SetTrackID(track);
192         
193         trackSave = track;
194         if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
195       }
196       
197       // track parameters at vertex 
198       TParticle* particle = stack->Particle(muonTrack.GetTrackID());
199       
200       if (particle) 
201       {        
202         Float_t x = particle->Vx();
203         Float_t y = particle->Vy();
204         Float_t z = particle->Vz();
205         Float_t pX = particle->Px();
206         Float_t pY = particle->Py();
207         Float_t pZ = particle->Pz();
208         
209         AliMUONTrackParam trackParam;
210         
211         trackParam.SetBendingCoor(y);
212         trackParam.SetNonBendingCoor(x);
213         trackParam.SetZ(z);
214         
215         Float_t bendingSlope = 0;
216         Float_t nonBendingSlope = 0;
217         Float_t inverseBendingMomentum = 0;
218                 
219         if (TMath::Abs(pZ) > 0) 
220         {
221           bendingSlope = pY/pZ;
222           nonBendingSlope = pX/pZ;
223         }
224         
225         Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
226         if (pYZ >0) inverseBendingMomentum = 1/pYZ;       
227         
228         TParticlePDG* ppdg = particle->GetPDG(1);
229         Int_t charge = (Int_t)(ppdg->Charge()/3.0);
230         inverseBendingMomentum *= charge;
231         
232         trackParam.SetBendingSlope(bendingSlope);
233         trackParam.SetNonBendingSlope(nonBendingSlope);
234         trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
235         
236         muonTrack.SetTrackParamAtVertex(&trackParam);
237       }
238
239       trackRefStore->Add(muonTrack);
240     } // end while isNewTrack
241   }
242   
243   AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
244   
245   delete trackRefStore;
246
247   return rv;
248 }
249
250 //_____________________________________________________________________________
251 Int_t
252 AliMUONRecoCheck::NumberOfEvents() const
253 {
254   /// Return the number of events
255   if ( fDataInterface ) 
256   {
257     return fDataInterface->NumberOfEvents();
258   }
259   return 0;
260 }
261
262 //_____________________________________________________________________________
263 AliMUONVTrackStore*
264 AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
265 {
266   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
267   /// each chamber (one when particle is entering + one when particle is leaving 
268   /// the sensitive volume) 
269   
270   Float_t maxGasGap = 1.; // cm 
271   AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
272   AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
273   TClonesArray *  hitForRecAtHit = 0;
274   TClonesArray *  trackParamAtHit = 0;
275   Float_t xRec,yRec,zRec;
276   Float_t xRec1,yRec1,zRec1;
277   Float_t xRec2,yRec2,zRec2;
278   Float_t bendingSlope,nonBendingSlope,bendingMomentum;
279   Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
280   Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
281   
282   AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
283   Int_t iHit1;
284   Int_t iChamber = 0, detElemId = 0;
285   Int_t nRec = 0;
286   Int_t nTrackHits = 0;
287   
288   hitForRec = new AliMUONHitForRec();
289   trackParam = new AliMUONTrackParam();
290   
291   TIter next(trackRefs.CreateIterator());
292   AliMUONTrack* track;
293   
294   while ( ( track = static_cast<AliMUONTrack*>(next())) ) 
295   {
296     hitForRecAtHit = track->GetHitForRecAtHit();
297     trackParamAtHit = track->GetTrackParamAtHit();
298     trackParamAtVertex = track->GetTrackParamAtVertex();
299     nTrackHits = hitForRecAtHit->GetEntriesFast();
300     AliMUONTrack trackNew;
301     iHit1 = 0;
302     while (iHit1 < nTrackHits) 
303     {
304       hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1); 
305       trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1); 
306       xRec1  = hitForRec1->GetNonBendingCoor();
307       yRec1  = hitForRec1->GetBendingCoor();
308       zRec1  = hitForRec1->GetZ();      
309       xRec   = xRec1;
310       yRec   = yRec1;
311       zRec   = zRec1;
312       bendingSlope1 = trackParam1->GetBendingSlope();
313       nonBendingSlope1 = trackParam1->GetNonBendingSlope();
314       bendingMomentum1 = 0;
315       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
316         bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
317       bendingSlope = bendingSlope1;
318       nonBendingSlope = nonBendingSlope1;
319       bendingMomentum = bendingMomentum1;
320       nRec = 1;  
321       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
322       {
323         hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2); 
324         trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2); 
325         xRec2  = hitForRec2->GetNonBendingCoor();
326         yRec2  = hitForRec2->GetBendingCoor();
327         zRec2  = hitForRec2->GetZ();      
328         bendingSlope2 = trackParam2->GetBendingSlope();
329         nonBendingSlope2 = trackParam2->GetNonBendingSlope();
330         bendingMomentum2 = 0;
331         if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
332           bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
333         
334         if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
335           
336           nRec++;
337           xRec += xRec2;
338           yRec += yRec2;
339           zRec += zRec2;
340           bendingSlope += bendingSlope2;
341           nonBendingSlope += nonBendingSlope2;
342           bendingMomentum += bendingMomentum2;
343           iHit1 = iHit2;
344         }
345         
346       } // end iHit2
347       xRec /= (Float_t)nRec;
348       yRec /= (Float_t)nRec;
349       zRec /= (Float_t)nRec;
350       bendingSlope /= (Float_t)nRec;
351       nonBendingSlope /= (Float_t)nRec;
352       bendingMomentum /= (Float_t)nRec;
353       
354       hitForRec->SetNonBendingCoor(xRec);
355       hitForRec->SetBendingCoor(yRec);
356       hitForRec->SetZ(zRec);
357       detElemId = hitForRec->GetDetElemId();
358       if (detElemId) iChamber = detElemId / 100 - 1;
359       else iChamber = AliMUONConstants::ChamberNumber(zRec);
360       hitForRec->SetChamberNumber(iChamber);
361       hitForRec->SetBendingReso2(0.0); 
362       hitForRec->SetNonBendingReso2(0.0); 
363       trackParam->SetNonBendingCoor(xRec);
364       trackParam->SetBendingCoor(yRec);
365       trackParam->SetZ(zRec);
366       trackParam->SetNonBendingSlope(nonBendingSlope);
367       trackParam->SetBendingSlope(bendingSlope);
368       if (TMath::Abs(bendingMomentum) > 0)
369         trackParam->SetInverseBendingMomentum(1./bendingMomentum);
370       
371       trackNew.AddHitForRecAtHit(hitForRec);
372       trackNew.AddTrackParamAtHit(trackParam,0);
373       
374       iHit1++;
375     } // end iHit1
376     
377     trackNew.SetTrackID(track->GetTrackID());
378     trackNew.SetTrackParamAtVertex(trackParamAtVertex);
379     newMuonTrackRef->Add(trackNew);
380     
381   } // end trackRef
382   
383   delete hitForRec;
384   delete trackParam;
385   return newMuonTrackRef;  
386 }
387
388 //_____________________________________________________________________________
389 AliMUONVTrackStore*
390 AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
391 {
392   /// Calculate the number of reconstructible tracks
393   
394   AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
395   
396   TClonesArray* hitForRecAtHit = NULL;
397   AliMUONHitForRec* hitForRec;
398   Float_t zRec;
399   Int_t nTrackHits;
400   Int_t isChamberInTrack[10];
401   Int_t iChamber = 0;
402   Bool_t isTrackOK = kTRUE;
403     
404   TIter next(trackRefs.CreateIterator());
405   AliMUONTrack* track;
406
407   while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
408   {
409     hitForRecAtHit = track->GetHitForRecAtHit();
410     nTrackHits = hitForRecAtHit->GetEntriesFast();
411     for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
412     
413     for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
414       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
415       zRec  = hitForRec->GetZ();
416       iChamber = hitForRec->GetChamberNumber(); 
417       if (iChamber < 0 || iChamber > 10) continue;
418       isChamberInTrack[iChamber] = 1;
419     } 
420     // track is reconstructible if the particle is depositing a hit
421     // in the following chamber combinations:
422     
423     isTrackOK = kTRUE;
424     if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
425     if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
426     if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
427     Int_t nHitsInLastStations=0;
428     for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
429       if (isChamberInTrack[ch]) nHitsInLastStations++; 
430     if(nHitsInLastStations < 3) isTrackOK = kFALSE;
431     
432     if (isTrackOK) 
433     {
434       reconstructibleStore->Add(*track);
435     }
436   }
437
438   return reconstructibleStore;
439 }
440