0daff813097ecc89e1704f8b7f89e14f9ee97747
[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)
224 {
225   /// Return a track store containing the reconstructible tracks for a given event
226
227   if (!fESDEventOwner) {
228     if (fRecoTrackRefStore == 0x0) {
229       if (TrackRefs(event) == 0x0) return 0x0;
230       MakeReconstructibleTracks();
231     }
232     return fRecoTrackRefStore;
233   }
234
235   if (event != fCurrentEvent) {
236     ResetStores();
237     fCurrentEvent = event;
238   }
239   
240   if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
241   else {
242     if (TrackRefs(event) == 0x0) return 0x0;
243     MakeReconstructibleTracks();
244     return fRecoTrackRefStore;
245   }
246 }
247
248 //_____________________________________________________________________________
249 void AliMUONRecoCheck::MakeReconstructedTracks(Bool_t refit)
250 {
251   /// Make reconstructed tracks
252   if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
253   
254   // loop over all reconstructed tracks and add them to the store (skip ghosts)
255   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
256   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
257     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
258     if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore, refit);
259   }
260   
261 }
262
263 //_____________________________________________________________________________
264 void AliMUONRecoCheck::MakeTrackRefs()
265 {
266   /// Make reconstructible tracks
267   AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
268   if (!tmpTrackRefStore) return;
269   
270   Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
271   TParticle* particle;
272   TClonesArray* trackRefs;
273   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
274   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
275   if (!cStore) return;
276   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
277   
278   // loop over simulated tracks
279   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
280     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
281     
282     // skip empty trackRefs
283     if (nHits < 1) continue;
284     
285     // get the particle charge for further calculation
286     TParticlePDG* ppdg = particle->GetPDG();
287     Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
288     
289     AliMUONTrack track;
290     
291     // loop over simulated track hits
292     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
293       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
294       
295       // skip trackRefs not in MUON
296       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
297     
298       // Get track parameters of current hit
299       x = trackReference->X();
300       y = trackReference->Y();
301       z = trackReference->Z();
302       pX = trackReference->Px();
303       pY = trackReference->Py();
304       pZ = trackReference->Pz();
305       
306       // check chamberId of current trackReference
307       Int_t detElemId = trackReference->UserId();
308       Int_t chamberId = detElemId / 100 - 1;
309       if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
310       
311       // set hit parameters
312       hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
313       hit->SetXYZ(x,y,z);
314       hit->SetErrXY(0.,0.);
315       
316       // compute track parameters at hit
317       bendingSlope = 0;
318       nonBendingSlope = 0;
319       inverseBendingMomentum = 0;
320       if (TMath::Abs(pZ) > 0) {
321         bendingSlope = pY/pZ;
322         nonBendingSlope = pX/pZ;
323       }
324       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
325       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
326       inverseBendingMomentum *= charge;
327       
328       // set track parameters at hit
329       AliMUONTrackParam trackParam;
330       trackParam.SetNonBendingCoor(x);
331       trackParam.SetBendingCoor(y);
332       trackParam.SetZ(z);
333       trackParam.SetBendingSlope(bendingSlope);
334       trackParam.SetNonBendingSlope(nonBendingSlope);
335       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
336       
337       // add track parameters at current hit to the track
338       track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
339     }
340     
341     // if none of the track hits was in MUON, goto the next track
342     if (track.GetNClusters() < 1) continue;
343     
344     // get track parameters at particle's vertex
345     x = particle->Vx();
346     y = particle->Vy();
347     z = particle->Vz();
348     pX = particle->Px();
349     pY = particle->Py();
350     pZ = particle->Pz();
351     
352     // compute rest of track parameters at particle's vertex
353     bendingSlope = 0;
354     nonBendingSlope = 0;
355     inverseBendingMomentum = 0;
356     if (TMath::Abs(pZ) > 0) {
357       bendingSlope = pY/pZ;
358       nonBendingSlope = pX/pZ;
359     }
360     Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
361     if (pYZ >0) inverseBendingMomentum = 1/pYZ;
362     inverseBendingMomentum *= charge;
363     
364     // set track parameters at particle's vertex
365     AliMUONTrackParam trackParamAtVertex;
366     trackParamAtVertex.SetNonBendingCoor(x);
367     trackParamAtVertex.SetBendingCoor(y);
368     trackParamAtVertex.SetZ(z);
369     trackParamAtVertex.SetBendingSlope(bendingSlope);
370     trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
371     trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
372     
373     // add track parameters at vertex
374     track.SetTrackParamAtVertex(&trackParamAtVertex);
375     
376     // store the track
377     track.SetUniqueID(iTrackRef);
378     tmpTrackRefStore->Add(track);
379   }
380   
381   CleanMuonTrackRef(tmpTrackRefStore);
382   
383   delete hit;
384   delete cStore;
385   delete tmpTrackRefStore;
386 }
387
388 //_____________________________________________________________________________
389 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
390 {
391   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
392   /// each chamber (one when particle is entering + one when particle is leaving 
393   /// the sensitive volume) 
394   if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
395   
396   Double_t maxGasGap = 1.; // cm 
397   Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
398   Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
399   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
400   if (!cStore) return;
401   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
402   
403   // create iterator
404   TIter next(tmpTrackRefStore->CreateIterator());
405   
406   // loop over tmpTrackRef
407   AliMUONTrack* track;
408   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
409     
410     AliMUONTrack newTrack;
411     
412     // loop over tmpTrackRef's hits
413     Int_t iHit1 = 0;
414     Int_t nTrackHits = track->GetNClusters();
415     while (iHit1 < nTrackHits) {
416       AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
417       
418       // get track parameters at hit1
419       x1  = trackParam1->GetNonBendingCoor();
420       y1  = trackParam1->GetBendingCoor();
421       z1  = trackParam1->GetZ();
422       pX1 = trackParam1->Px();
423       pY1 = trackParam1->Py();
424       pZ1 = trackParam1->Pz();
425       
426       // prepare new track parameters
427       x  = x1;
428       y  = y1;
429       z  = z1;
430       pX = pX1;
431       pY = pY1;
432       pZ = pZ1;
433       
434       // loop over next tmpTrackRef's hits
435       Int_t nCombinedHits = 1;
436       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
437         AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
438         
439         // get z position of hit2
440         z2 = trackParam2->GetZ();
441         
442         // complete new track parameters if hit2 is on the same detection element
443         if ( TMath::Abs(z2-z1) < maxGasGap ) {
444           x  += trackParam2->GetNonBendingCoor();
445           y  += trackParam2->GetBendingCoor();
446           z  += z2;
447           pX += trackParam2->Px();
448           pY += trackParam2->Py();
449           pZ += trackParam2->Pz();
450           nCombinedHits++;
451           iHit1 = iHit2;
452         }
453         
454       }
455       
456       // finalize new track parameters
457       x  /= (Double_t)nCombinedHits;
458       y  /= (Double_t)nCombinedHits;
459       z  /= (Double_t)nCombinedHits;
460       pX /= (Double_t)nCombinedHits;
461       pY /= (Double_t)nCombinedHits;
462       pZ /= (Double_t)nCombinedHits;
463       bendingSlope = 0;
464       nonBendingSlope = 0;
465       inverseBendingMomentum = 0;
466       if (TMath::Abs(pZ) > 0) {
467         bendingSlope = pY/pZ;
468         nonBendingSlope = pX/pZ;
469       }
470       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
471       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
472       inverseBendingMomentum *= trackParam1->GetCharge();
473       
474       // set hit parameters
475       hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
476       hit->SetXYZ(x,y,z);
477       hit->SetErrXY(0.,0.);
478       
479       // set new track parameters at new hit
480       AliMUONTrackParam trackParam;
481       trackParam.SetNonBendingCoor(x);
482       trackParam.SetBendingCoor(y);
483       trackParam.SetZ(z);
484       trackParam.SetBendingSlope(bendingSlope);
485       trackParam.SetNonBendingSlope(nonBendingSlope);
486       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
487       
488       // add track parameters at current hit to the track
489       newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
490       
491       iHit1++;
492     }
493     
494     newTrack.SetUniqueID(track->GetUniqueID());
495     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
496     fTrackRefStore->Add(newTrack);
497     
498   }
499   
500   delete hit;
501   delete cStore;
502 }
503
504 //_____________________________________________________________________________
505 void AliMUONRecoCheck::MakeReconstructibleTracks()
506 {
507   /// Isolate the reconstructible tracks
508   if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
509   
510   // create iterator on trackRef
511   TIter next(fTrackRefStore->CreateIterator());
512   
513   // loop over trackRef
514   AliMUONTrack* track;
515   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
516     
517     Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
518     for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
519     
520     // loop over trackRef's hits to get hit chambers
521     Int_t nTrackHits = track->GetNClusters();
522     for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
523       AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr(); 
524       chamberInTrack[hit->GetChamberId()] = kTRUE;
525     } 
526     
527     // track is reconstructible if the particle is depositing a hit
528     // in the following chamber combinations:
529     Bool_t trackOK = kTRUE;
530     if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
531     if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
532     if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
533     Int_t nHitsInLastStations = 0;
534     for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
535       if (chamberInTrack[iCh]) nHitsInLastStations++; 
536     if(nHitsInLastStations < 3) trackOK = kFALSE;
537     
538     // Add reconstructible tracks to fRecoTrackRefStore
539     if (trackOK) fRecoTrackRefStore->Add(*track);
540     
541     delete [] chamberInTrack;
542   }
543
544 }
545