Fixes in reconstruction:
[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 // 13 Nov 2007:
27 // Added a method to create a list of reconstructed AliMUONTrack objects from
28 // ESD data. This is necessary since the track objects that are actually created
29 // during offline reconstruction are no longer stored to disk.
30 //  - Artur Szostak <artursz@iafrica.com>
31 // 25 Jan 2008:
32 // Use the new ESDInterface to create MUON objects from ESD data
33 // - Philippe Pillot
34 // 
35
36 #include "AliMUONRecoCheck.h"
37 #include "AliMUONTrack.h"
38 #include "AliMUONVTrackStore.h"
39 #include "AliMUONVCluster.h"
40 #include "AliMUONVClusterStore.h"
41 #include "AliMUONConstants.h"
42 #include "AliMUONESDInterface.h"
43 #include "AliMUONTrackParam.h"
44 #include "AliMCEventHandler.h"
45 #include "AliMCEvent.h"
46 #include "AliStack.h"
47 #include "AliTrackReference.h"
48 #include "AliLog.h"
49 #include "AliESDEvent.h"
50 #include "AliESDMuonTrack.h"
51
52 #include <TFile.h>
53 #include <TTree.h>
54 #include <TParticle.h>
55 #include <TParticlePDG.h>
56 #include <Riostream.h>
57
58 /// \cond CLASSIMP
59 ClassImp(AliMUONRecoCheck)
60 /// \endcond
61
62 //_____________________________________________________________________________
63 AliMUONRecoCheck::AliMUONRecoCheck(const Char_t *esdFileName, const Char_t *pathSim)
64 : TObject(),
65 fMCEventHandler(new AliMCEventHandler()),
66 fESDEvent(new AliESDEvent()),
67 fESDTree (0x0),
68 fESDFile (0x0),
69 fCurrentEvent(0),
70 fTrackRefStore(0x0),
71 fRecoTrackRefStore(0x0),
72 fRecoTrackStore(0x0),
73 fESDEventOwner(kTRUE)
74 {
75   /// Normal ctor
76   
77   // TrackRefs and Particules
78   fMCEventHandler->SetInputPath(pathSim);
79   fMCEventHandler->InitIO("");
80   
81   // ESD MUON Tracks
82   fESDFile = TFile::Open(esdFileName); // open the file
83   if (!fESDFile || !fESDFile->IsOpen()) {
84     AliError(Form("opening ESD file %s failed", esdFileName));
85     fESDFile = 0x0;
86     return;
87   }
88   fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
89   if (!fESDTree) {
90     AliError("no ESD tree found");
91     fESDFile->Close();
92     fESDFile = 0x0;
93     return;
94   }
95   fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
96 }
97
98 //_____________________________________________________________________________
99 AliMUONRecoCheck::AliMUONRecoCheck(AliESDEvent *esdEvent, AliMCEventHandler *mcEventHandler)
100 : TObject(),
101 fMCEventHandler(0),
102 fESDEvent(0),
103 fESDTree (0x0),
104 fESDFile (0x0),
105 fCurrentEvent(0),
106 fTrackRefStore(0x0),
107 fRecoTrackRefStore(0x0),
108 fRecoTrackStore(0x0),
109 fESDEventOwner(kFALSE)
110 {
111   /// Normal ctor
112   
113   // TrackRefs and Particules
114   fMCEventHandler = mcEventHandler;
115   
116   // ESD MUON Tracks
117   fESDEvent = esdEvent;
118   
119 }
120
121 //_____________________________________________________________________________
122 AliMUONRecoCheck::~AliMUONRecoCheck()
123 {
124   /// Destructor
125   if (fESDEventOwner) {
126     delete fMCEventHandler;
127     delete fESDEvent;
128     if (fESDFile) fESDFile->Close();
129   }
130   ResetStores();
131 }
132
133 //_____________________________________________________________________________
134 void AliMUONRecoCheck::ResetStores()
135 {
136   /// Deletes all the store objects that have been created and resets the pointers to 0x0
137   delete fTrackRefStore;      fTrackRefStore = 0x0;
138   delete fRecoTrackRefStore;  fRecoTrackRefStore = 0x0;
139   delete fRecoTrackStore;     fRecoTrackStore = 0x0;
140 }
141
142 //_____________________________________________________________________________
143 Int_t AliMUONRecoCheck::GetRunNumber()
144 {
145   /// Return the run number of the current ESD event
146   
147   if (fESDEventOwner && fRecoTrackStore == 0x0) {
148     if (!fESDTree || fESDTree->GetEvent(fCurrentEvent) <= 0) {
149       AliError(Form("fails to read ESD object for event %d: cannot get the run number",fCurrentEvent));
150       return -1;
151     }
152   }
153   
154   return fESDEvent->GetRunNumber();
155 }
156
157 //_____________________________________________________________________________
158 Int_t AliMUONRecoCheck::NumberOfEvents() const
159 {
160   /// Return the number of events
161   if (fESDEventOwner && fESDTree) return fESDTree->GetEntries();
162   return 0;
163 }
164
165 //_____________________________________________________________________________
166 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event, Bool_t refit)
167 {
168   /// Return a track store containing the reconstructed tracks (converted into 
169   /// MUONTrack objects) for a given event.
170   /// Track parameters at each clusters are computed or not depending on the flag "refit".
171   /// If not, only the track parameters at first cluster are valid.
172   
173   if (!fESDEventOwner) {
174     if (fRecoTrackStore == 0x0) MakeReconstructedTracks(refit);
175     return fRecoTrackStore;
176   }
177
178   if (event != fCurrentEvent) {
179     ResetStores();
180     fCurrentEvent = event;
181   }
182   
183   if (fRecoTrackStore != 0x0) return fRecoTrackStore;
184   else {
185     if (!fESDTree) return 0x0;
186     if (fESDTree->GetEvent(event) <= 0) {
187       AliError(Form("fails to read ESD object for event %d", event));
188       return 0x0;
189     }
190     MakeReconstructedTracks(refit);
191     return fRecoTrackStore;
192   }
193 }
194
195 //_____________________________________________________________________________
196 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
197 {
198   /// Return a track store containing the track references (converted into 
199   /// MUONTrack objects) for a given event
200   
201   if (!fESDEventOwner) {
202     if (fTrackRefStore == 0x0) MakeTrackRefs();
203     return fTrackRefStore;
204   }
205
206   if (event != fCurrentEvent) {
207     ResetStores();
208     fCurrentEvent = event;
209   }
210   
211   if (fTrackRefStore != 0x0) return fTrackRefStore;
212   else {
213     if (!fMCEventHandler->GetEvent(event)) {
214       AliError(Form("fails to read MC objects for event %d", event));
215       return 0x0;
216     }
217     MakeTrackRefs();
218     return fTrackRefStore;
219   }
220 }
221
222 //_____________________________________________________________________________
223 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event, UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
224 {
225   /// Return a track store containing the reconstructible tracks for a given event,
226   /// according to the mask of requested stations and the minimum number of chambers hit in stations 4 & 5.
227
228   if (!fESDEventOwner) {
229     if (fRecoTrackRefStore == 0x0) {
230       if (TrackRefs(event) == 0x0) return 0x0;
231       MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
232     }
233     return fRecoTrackRefStore;
234   }
235
236   if (event != fCurrentEvent) {
237     ResetStores();
238     fCurrentEvent = event;
239   }
240   
241   if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
242   else {
243     if (TrackRefs(event) == 0x0) return 0x0;
244     MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
245     return fRecoTrackRefStore;
246   }
247 }
248
249 //_____________________________________________________________________________
250 void AliMUONRecoCheck::MakeReconstructedTracks(Bool_t refit)
251 {
252   /// Make reconstructed tracks
253   if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
254   
255   // loop over all reconstructed tracks and add them to the store (skip ghosts)
256   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
257   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
258     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
259     if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore, refit);
260   }
261   
262 }
263
264 //_____________________________________________________________________________
265 void AliMUONRecoCheck::MakeTrackRefs()
266 {
267   /// Make reconstructible tracks
268   AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
269   if (!tmpTrackRefStore) return;
270   
271   Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
272   TParticle* particle;
273   TClonesArray* trackRefs;
274   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
275   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
276   if (!cStore) return;
277   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
278   
279   // loop over simulated tracks
280   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
281     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
282     
283     // skip empty trackRefs
284     if (nHits < 1) continue;
285     
286     // get the particle charge for further calculation
287     TParticlePDG* ppdg = particle->GetPDG();
288     Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
289     
290     AliMUONTrack track;
291     
292     // loop over simulated track hits
293     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
294       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
295       
296       // skip trackRefs not in MUON
297       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
298     
299       // Get track parameters of current hit
300       x = trackReference->X();
301       y = trackReference->Y();
302       z = trackReference->Z();
303       pX = trackReference->Px();
304       pY = trackReference->Py();
305       pZ = trackReference->Pz();
306       
307       // check chamberId of current trackReference
308       Int_t detElemId = trackReference->UserId();
309       Int_t chamberId = detElemId / 100 - 1;
310       if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
311       
312       // set hit parameters
313       hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
314       hit->SetXYZ(x,y,z);
315       hit->SetErrXY(0.,0.);
316       
317       // compute track parameters at hit
318       bendingSlope = 0;
319       nonBendingSlope = 0;
320       inverseBendingMomentum = 0;
321       if (TMath::Abs(pZ) > 0) {
322         bendingSlope = pY/pZ;
323         nonBendingSlope = pX/pZ;
324       }
325       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
326       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
327       inverseBendingMomentum *= charge;
328       
329       // set track parameters at hit
330       AliMUONTrackParam trackParam;
331       trackParam.SetNonBendingCoor(x);
332       trackParam.SetBendingCoor(y);
333       trackParam.SetZ(z);
334       trackParam.SetBendingSlope(bendingSlope);
335       trackParam.SetNonBendingSlope(nonBendingSlope);
336       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
337       
338       // add track parameters at current hit to the track
339       track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
340     }
341     
342     // if none of the track hits was in MUON, goto the next track
343     if (track.GetNClusters() < 1) continue;
344     
345     // get track parameters at particle's vertex
346     x = particle->Vx();
347     y = particle->Vy();
348     z = particle->Vz();
349     pX = particle->Px();
350     pY = particle->Py();
351     pZ = particle->Pz();
352     
353     // compute rest of track parameters at particle's vertex
354     bendingSlope = 0;
355     nonBendingSlope = 0;
356     inverseBendingMomentum = 0;
357     if (TMath::Abs(pZ) > 0) {
358       bendingSlope = pY/pZ;
359       nonBendingSlope = pX/pZ;
360     }
361     Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
362     if (pYZ >0) inverseBendingMomentum = 1/pYZ;
363     inverseBendingMomentum *= charge;
364     
365     // set track parameters at particle's vertex
366     AliMUONTrackParam trackParamAtVertex;
367     trackParamAtVertex.SetNonBendingCoor(x);
368     trackParamAtVertex.SetBendingCoor(y);
369     trackParamAtVertex.SetZ(z);
370     trackParamAtVertex.SetBendingSlope(bendingSlope);
371     trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
372     trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
373     
374     // add track parameters at vertex
375     track.SetTrackParamAtVertex(&trackParamAtVertex);
376     
377     // store the track
378     track.SetUniqueID(iTrackRef);
379     tmpTrackRefStore->Add(track);
380   }
381   
382   CleanMuonTrackRef(tmpTrackRefStore);
383   
384   delete hit;
385   delete cStore;
386   delete tmpTrackRefStore;
387 }
388
389 //_____________________________________________________________________________
390 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
391 {
392   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
393   /// each chamber (one when particle is entering + one when particle is leaving 
394   /// the sensitive volume) 
395   if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
396   
397   Double_t maxGasGap = 1.; // cm 
398   Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
399   Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
400   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
401   if (!cStore) return;
402   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
403   
404   // create iterator
405   TIter next(tmpTrackRefStore->CreateIterator());
406   
407   // loop over tmpTrackRef
408   AliMUONTrack* track;
409   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
410     
411     AliMUONTrack newTrack;
412     
413     // loop over tmpTrackRef's hits
414     Int_t iHit1 = 0;
415     Int_t nTrackHits = track->GetNClusters();
416     while (iHit1 < nTrackHits) {
417       AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
418       
419       // get track parameters at hit1
420       x1  = trackParam1->GetNonBendingCoor();
421       y1  = trackParam1->GetBendingCoor();
422       z1  = trackParam1->GetZ();
423       pX1 = trackParam1->Px();
424       pY1 = trackParam1->Py();
425       pZ1 = trackParam1->Pz();
426       
427       // prepare new track parameters
428       x  = x1;
429       y  = y1;
430       z  = z1;
431       pX = pX1;
432       pY = pY1;
433       pZ = pZ1;
434       
435       // loop over next tmpTrackRef's hits
436       Int_t nCombinedHits = 1;
437       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
438         AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
439         
440         // get z position of hit2
441         z2 = trackParam2->GetZ();
442         
443         // complete new track parameters if hit2 is on the same detection element
444         if ( TMath::Abs(z2-z1) < maxGasGap ) {
445           x  += trackParam2->GetNonBendingCoor();
446           y  += trackParam2->GetBendingCoor();
447           z  += z2;
448           pX += trackParam2->Px();
449           pY += trackParam2->Py();
450           pZ += trackParam2->Pz();
451           nCombinedHits++;
452           iHit1 = iHit2;
453         }
454         
455       }
456       
457       // finalize new track parameters
458       x  /= (Double_t)nCombinedHits;
459       y  /= (Double_t)nCombinedHits;
460       z  /= (Double_t)nCombinedHits;
461       pX /= (Double_t)nCombinedHits;
462       pY /= (Double_t)nCombinedHits;
463       pZ /= (Double_t)nCombinedHits;
464       bendingSlope = 0;
465       nonBendingSlope = 0;
466       inverseBendingMomentum = 0;
467       if (TMath::Abs(pZ) > 0) {
468         bendingSlope = pY/pZ;
469         nonBendingSlope = pX/pZ;
470       }
471       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
472       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
473       inverseBendingMomentum *= trackParam1->GetCharge();
474       
475       // set hit parameters
476       hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
477       hit->SetXYZ(x,y,z);
478       hit->SetErrXY(0.,0.);
479       
480       // set new track parameters at new hit
481       AliMUONTrackParam trackParam;
482       trackParam.SetNonBendingCoor(x);
483       trackParam.SetBendingCoor(y);
484       trackParam.SetZ(z);
485       trackParam.SetBendingSlope(bendingSlope);
486       trackParam.SetNonBendingSlope(nonBendingSlope);
487       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
488       
489       // add track parameters at current hit to the track
490       newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
491       
492       iHit1++;
493     }
494     
495     newTrack.SetUniqueID(track->GetUniqueID());
496     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
497     fTrackRefStore->Add(newTrack);
498     
499   }
500   
501   delete hit;
502   delete cStore;
503 }
504
505 //_____________________________________________________________________________
506 void AliMUONRecoCheck::MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
507 {
508   /// Isolate the reconstructible tracks
509   if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
510   
511   // create iterator on trackRef
512   TIter next(fTrackRefStore->CreateIterator());
513   
514   // loop over trackRef and add reconstructible tracks to fRecoTrackRefStore
515   AliMUONTrack* track;
516   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
517     if (track->IsValid(requestedStationMask, request2ChInSameSt45)) fRecoTrackRefStore->Add(*track);
518   }
519
520 }
521
522 //_____________________________________________________________________________
523 AliMUONTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore,
524                                                     Int_t &nMatchClusters, Bool_t useLabel, Double_t sigmaCut)
525 {
526   /// Return the track from the store matched with the given track (or 0x0) and the number of matched clusters.
527   /// Matching is done by using the MC label of by comparing cluster/TrackRef positions according to the flag "useLabel".
528   /// WARNING: Who match who matters since the matching algorithm uses the *fraction* of matched clusters of the given track
529   
530   AliMUONTrack *matchedTrack = 0x0;
531   nMatchClusters = 0;
532   
533   if (useLabel) { // by using the MC label
534     
535     // get the corresponding simulated track if any
536     Int_t label = track.GetMCLabel();
537     matchedTrack = (AliMUONTrack*) trackStore.FindObject(label);
538     
539     // get the fraction of matched clusters
540     if (matchedTrack) {
541       Int_t nClusters = track.GetNClusters();
542       for (Int_t iCl = 0; iCl < nClusters; iCl++)
543         if (((AliMUONTrackParam*) track.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
544           nMatchClusters++;
545     }
546     
547   } else { // by comparing cluster/TrackRef positions
548     
549     // look for the corresponding simulated track if any
550     TIter next(trackStore.CreateIterator());
551     AliMUONTrack* track2;
552     while ( ( track2 = static_cast<AliMUONTrack*>(next()) ) ) {
553       
554       // check compatibility
555       Int_t n = 0;
556       if (track.Match(*track2, sigmaCut, n)) {
557         matchedTrack = track2;
558         nMatchClusters = n;
559         break;
560       }
561       
562     }
563     
564   }
565   
566   return matchedTrack;
567   
568 }
569