- Adding check and flagging for HG present
[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 "AliMUONTriggerTrack.h"
45 #include "AliMUONVTriggerTrackStore.h"
46 #include "AliMCEventHandler.h"
47 #include "AliMCEvent.h"
48 #include "AliStack.h"
49 #include "AliTrackReference.h"
50 #include "AliLog.h"
51 #include "AliESDEvent.h"
52 #include "AliESDMuonTrack.h"
53
54 #include "AliGeomManager.h"
55 #include "AliMpCDB.h"
56 #include "AliMpDDLStore.h"
57 #include "AliMUONCDB.h"
58 #include "AliMUONGeometryTransformer.h"
59 #include "AliMUONTriggerCircuit.h"
60 #include "AliMUONVTrackReconstructor.h"
61 #include "AliMUONVTriggerStore.h"
62
63 #include "TGeoManager.h"
64
65 #include <TFile.h>
66 #include <TTree.h>
67 #include <TParticle.h>
68 #include <TParticlePDG.h>
69 #include <Riostream.h>
70
71 #include "AliMUONRecoCheck.h"
72
73 /// \cond CLASSIMP
74 ClassImp(AliMUONRecoCheck)
75 /// \endcond
76
77 //_____________________________________________________________________________
78 AliMUONRecoCheck::AliMUONRecoCheck(const Char_t *esdFileName, const Char_t *pathSim)
79 : TObject(),
80 fMCEventHandler(new AliMCEventHandler()),
81 fESDEvent(new AliESDEvent()),
82 fESDTree (0x0),
83 fESDFile (0x0),
84 fCurrentEvent(0),
85 fTrackRefStore(0x0),
86 fRecoTrackRefStore(0x0),
87 fRecoTriggerRefStore(0x0),
88 fRecoTrackStore(0x0),
89 fRecoTriggerTrackStore(0x0),
90 fGeometryTransformer(0x0),
91 fTriggerCircuit(0x0),
92 fESDEventOwner(kTRUE)
93 {
94   /// Normal ctor
95   
96   // TrackRefs and Particules
97   fMCEventHandler->SetInputPath(pathSim);
98   fMCEventHandler->InitIO("");
99   
100   // ESD MUON Tracks
101   fESDFile = TFile::Open(esdFileName); // open the file
102   if (!fESDFile || !fESDFile->IsOpen()) {
103     AliError(Form("opening ESD file %s failed", esdFileName));
104     fESDFile = 0x0;
105     return;
106   }
107   fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
108   if (!fESDTree) {
109     AliError("no ESD tree found");
110     fESDFile->Close();
111     fESDFile = 0x0;
112     return;
113   }
114   fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
115 }
116
117 //_____________________________________________________________________________
118 AliMUONRecoCheck::AliMUONRecoCheck(AliESDEvent *esdEvent, AliMCEventHandler *mcEventHandler)
119 : TObject(),
120 fMCEventHandler(0),
121 fESDEvent(0),
122 fESDTree (0x0),
123 fESDFile (0x0),
124 fCurrentEvent(0),
125 fTrackRefStore(0x0),
126 fRecoTrackRefStore(0x0),
127 fRecoTriggerRefStore(0x0),
128 fRecoTrackStore(0x0),
129 fRecoTriggerTrackStore(0x0),
130 fGeometryTransformer(0x0),
131 fTriggerCircuit(0x0),
132 fESDEventOwner(kFALSE)
133 {
134   /// Normal ctor
135   
136   // TrackRefs and Particules
137   fMCEventHandler = mcEventHandler;
138   
139   // ESD MUON Tracks
140   fESDEvent = esdEvent;
141 }
142
143 //_____________________________________________________________________________
144 AliMUONRecoCheck::~AliMUONRecoCheck()
145 {
146   /// Destructor
147   if (fESDEventOwner) {
148     delete fMCEventHandler;
149     delete fESDEvent;
150     if (fESDFile) fESDFile->Close();
151   }
152   ResetStores();
153   delete fGeometryTransformer;
154   delete fTriggerCircuit;
155 }
156
157 //_____________________________________________________________________________
158 void AliMUONRecoCheck::ResetStores()
159 {
160   /// Deletes all the store objects that have been created and resets the pointers to 0x0
161   delete fTrackRefStore;      fTrackRefStore = 0x0;
162   delete fRecoTrackRefStore;  fRecoTrackRefStore = 0x0;
163   delete fRecoTriggerRefStore;  fRecoTriggerRefStore = 0x0;
164   delete fRecoTrackStore;     fRecoTrackStore = 0x0;
165   delete fRecoTriggerTrackStore; fRecoTriggerTrackStore = 0x0;
166 }
167
168 //_____________________________________________________________________________
169 Bool_t AliMUONRecoCheck::InitCircuit()
170 {
171
172   if ( fTriggerCircuit ) return kTRUE;
173
174   if ( !AliMUONCDB::CheckOCDB() ) return kFALSE;
175
176   if ( !AliGeomManager::GetGeometry() )
177     AliGeomManager::LoadGeometry();
178
179   if ( !AliMpDDLStore::Instance(false) )
180     AliMpCDB::LoadDDLStore();
181         
182   fGeometryTransformer = new AliMUONGeometryTransformer();
183   fGeometryTransformer->LoadGeometryData();
184   
185   fTriggerCircuit = new AliMUONTriggerCircuit(fGeometryTransformer);
186
187   // reset tracker for local trigger to trigger track conversion
188   if ( ! AliMUONESDInterface::GetTracker() )
189     AliMUONESDInterface::ResetTracker();
190   
191   return kTRUE;
192 }
193
194
195 //_____________________________________________________________________________
196 Int_t AliMUONRecoCheck::GetRunNumber()
197 {
198   /// Return the run number of the current ESD event
199   
200   if (fESDEventOwner && fRecoTrackStore == 0x0 && fRecoTriggerTrackStore == 0x0) {
201     if (!fESDTree || fESDTree->GetEvent(fCurrentEvent) <= 0) {
202       AliError(Form("fails to read ESD object for event %d: cannot get the run number",fCurrentEvent));
203       return -1;
204     }
205   }
206   
207   return fESDEvent->GetRunNumber();
208 }
209
210 //_____________________________________________________________________________
211 Int_t AliMUONRecoCheck::NumberOfEvents() const
212 {
213   /// Return the number of events
214   if (fESDEventOwner && fESDTree) return fESDTree->GetEntries();
215   return 0;
216 }
217
218 //_____________________________________________________________________________
219 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event, Bool_t refit)
220 {
221   /// Return a track store containing the reconstructed tracks (converted into 
222   /// MUONTrack objects) for a given event.
223   /// Track parameters at each clusters are computed or not depending on the flag "refit".
224   /// If not, only the track parameters at first cluster are valid.
225   
226   if (!fESDEventOwner) {
227     if (fRecoTrackStore == 0x0) MakeReconstructedTracks(refit);
228     return fRecoTrackStore;
229   }
230
231   if (event != fCurrentEvent) {
232     ResetStores();
233     fCurrentEvent = event;
234   }
235   
236   if (fRecoTrackStore != 0x0) return fRecoTrackStore;
237   else {
238     if (!fESDTree) return 0x0;
239     if (fESDTree->GetEvent(event) <= 0) {
240       AliError(Form("fails to read ESD object for event %d", event));
241       return 0x0;
242     }
243     MakeReconstructedTracks(refit);
244     return fRecoTrackStore;
245   }
246 }
247
248
249 //_____________________________________________________________________________
250 AliMUONVTriggerTrackStore* AliMUONRecoCheck::TriggeredTracks(Int_t event)
251 {
252   /// Return a track store containing the reconstructed trigger tracks (converted into 
253   /// MUONTriggerTrack objects) for a given event.
254         
255   if (!fESDEventOwner) {
256     if (fRecoTriggerTrackStore == 0x0) MakeTriggeredTracks();
257     return fRecoTriggerTrackStore;
258   }
259         
260   if (event != fCurrentEvent) {
261     ResetStores();
262     fCurrentEvent = event;
263   }
264         
265   if (fRecoTriggerTrackStore != 0x0) return fRecoTriggerTrackStore;
266   else {
267     if (!fESDTree) return 0x0;
268     if (fESDTree->GetEvent(event) <= 0) {
269       AliError(Form("fails to read ESD object for event %d", event));
270       return 0x0;
271     }
272     MakeTriggeredTracks();
273     return fRecoTriggerTrackStore;
274   }
275 }
276
277
278 //_____________________________________________________________________________
279 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
280 {
281   /// Return a track store containing the track references (converted into 
282   /// MUONTrack objects) for a given event
283   
284   if (!fESDEventOwner) {
285     if (fTrackRefStore == 0x0) MakeTrackRefs();
286     return fTrackRefStore;
287   }
288
289   if (event != fCurrentEvent) {
290     ResetStores();
291     fCurrentEvent = event;
292   }
293   
294   if (fTrackRefStore != 0x0) return fTrackRefStore;
295   else {
296     if (!fMCEventHandler->GetEvent(event)) {
297       AliError(Form("fails to read MC objects for event %d", event));
298       return 0x0;
299     }
300     MakeTrackRefs();
301     return fTrackRefStore;
302   }
303 }
304
305 //_____________________________________________________________________________
306 AliMUONVTriggerTrackStore* AliMUONRecoCheck::TriggerableTracks(Int_t event)
307 {
308   /// Return a trigger track store containing the triggerable track references (converted into 
309   /// AliMUONTriggerTrack objects) for a given event
310         
311   if (!fESDEventOwner) {
312     if (fRecoTriggerRefStore == 0x0) MakeTriggerableTracks();
313     return fRecoTriggerRefStore;
314   }
315         
316   if (event != fCurrentEvent) {
317     ResetStores();
318     fCurrentEvent = event;
319   }
320         
321   if (fRecoTriggerRefStore != 0x0) return fRecoTriggerRefStore;
322   else {
323     if (!fMCEventHandler->GetEvent(event)) {
324       AliError(Form("fails to read MC objects for event %d", event));
325       return 0x0;
326     }
327     MakeTriggerableTracks();
328     return fRecoTriggerRefStore;
329   }
330 }
331
332
333 //_____________________________________________________________________________
334 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event, UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
335 {
336   /// Return a track store containing the reconstructible tracks for a given event,
337   /// according to the mask of requested stations and the minimum number of chambers hit in stations 4 & 5.
338
339   if (!fESDEventOwner) {
340     if (fRecoTrackRefStore == 0x0) {
341       if (TrackRefs(event) == 0x0) return 0x0;
342       MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
343     }
344     return fRecoTrackRefStore;
345   }
346
347   if (event != fCurrentEvent) {
348     ResetStores();
349     fCurrentEvent = event;
350   }
351   
352   if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
353   else {
354     if (TrackRefs(event) == 0x0) return 0x0;
355     MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
356     return fRecoTrackRefStore;
357   }
358 }
359
360
361 //_____________________________________________________________________________
362 void AliMUONRecoCheck::MakeReconstructedTracks(Bool_t refit)
363 {
364   /// Make reconstructed tracks
365   if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
366   
367   // loop over all reconstructed tracks and add them to the store (skip ghosts)
368   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
369   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
370     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
371     if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore, refit);
372   }
373   
374 }
375
376
377 //_____________________________________________________________________________
378 void AliMUONRecoCheck::MakeTriggeredTracks()
379 {
380   /// Make reconstructed trigger tracks
381   if (!(fRecoTriggerTrackStore = AliMUONESDInterface::NewTriggerTrackStore())) return;
382         
383   AliMUONVTriggerStore* tmpTriggerStore = AliMUONESDInterface::NewTriggerStore();
384   if ( ! tmpTriggerStore ) return;
385   
386   // loop over all reconstructed tracks and add them to the store (include ghosts)
387   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
388   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
389     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
390     if (esdTrack->ContainTriggerData()) AliMUONESDInterface::Add(*esdTrack, *tmpTriggerStore);
391   }
392         
393   if ( ! InitCircuit() ) return;
394   
395   AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
396   tracker->EventReconstructTrigger(*fTriggerCircuit, *tmpTriggerStore, *fRecoTriggerTrackStore);
397   
398   delete tmpTriggerStore;
399 }
400
401 //_____________________________________________________________________________
402 void AliMUONRecoCheck::TriggerToTrack(const AliMUONLocalTrigger& locTrg, AliMUONTriggerTrack& triggerTrack)
403 {
404   /// Make trigger track from local trigger info
405   if ( ! InitCircuit() ) return;
406   AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
407   tracker->TriggerToTrack(*fTriggerCircuit, locTrg, triggerTrack);
408 }
409
410
411 //_____________________________________________________________________________
412 void AliMUONRecoCheck::MakeTrackRefs()
413 {
414   /// Make reconstructible tracks
415   AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
416   if (!tmpTrackRefStore) return;
417   
418   Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
419   TParticle* particle;
420   TClonesArray* trackRefs;
421   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
422   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
423   if (!cStore) return;
424   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
425   
426   // loop over simulated tracks
427   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
428     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
429     
430     // skip empty trackRefs
431     if (nHits < 1) continue;
432     
433     // get the particle charge for further calculation
434     TParticlePDG* ppdg = particle->GetPDG();
435     Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
436     
437     AliMUONTrack track;
438     
439     // loop over simulated track hits
440     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
441       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
442       
443       // skip trackRefs not in MUON
444       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
445     
446       // Get track parameters of current hit
447       x = trackReference->X();
448       y = trackReference->Y();
449       z = trackReference->Z();
450       pX = trackReference->Px();
451       pY = trackReference->Py();
452       pZ = trackReference->Pz();
453       
454       // check chamberId of current trackReference
455       Int_t detElemId = trackReference->UserId();
456       Int_t chamberId = detElemId / 100 - 1;
457       if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
458       
459       // set hit parameters
460       hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
461       hit->SetXYZ(x,y,z);
462       hit->SetErrXY(0.,0.);
463       
464       // compute track parameters at hit
465       bendingSlope = 0;
466       nonBendingSlope = 0;
467       inverseBendingMomentum = 0;
468       if (TMath::Abs(pZ) > 0) {
469         bendingSlope = pY/pZ;
470         nonBendingSlope = pX/pZ;
471       }
472       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
473       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
474       inverseBendingMomentum *= charge;
475       
476       // set track parameters at hit
477       AliMUONTrackParam trackParam;
478       trackParam.SetNonBendingCoor(x);
479       trackParam.SetBendingCoor(y);
480       trackParam.SetZ(z);
481       trackParam.SetBendingSlope(bendingSlope);
482       trackParam.SetNonBendingSlope(nonBendingSlope);
483       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
484       
485       // add track parameters at current hit to the track
486       track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
487     }
488     
489     // if none of the track hits was in MUON, goto the next track
490     if (track.GetNClusters() < 1) continue;
491     
492     // get track parameters at particle's vertex
493     x = particle->Vx();
494     y = particle->Vy();
495     z = particle->Vz();
496     pX = particle->Px();
497     pY = particle->Py();
498     pZ = particle->Pz();
499     
500     // compute rest of track parameters at particle's vertex
501     bendingSlope = 0;
502     nonBendingSlope = 0;
503     inverseBendingMomentum = 0;
504     if (TMath::Abs(pZ) > 0) {
505       bendingSlope = pY/pZ;
506       nonBendingSlope = pX/pZ;
507     }
508     Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
509     if (pYZ >0) inverseBendingMomentum = 1/pYZ;
510     inverseBendingMomentum *= charge;
511     
512     // set track parameters at particle's vertex
513     AliMUONTrackParam trackParamAtVertex;
514     trackParamAtVertex.SetNonBendingCoor(x);
515     trackParamAtVertex.SetBendingCoor(y);
516     trackParamAtVertex.SetZ(z);
517     trackParamAtVertex.SetBendingSlope(bendingSlope);
518     trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
519     trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
520     
521     // add track parameters at vertex
522     track.SetTrackParamAtVertex(&trackParamAtVertex);
523     
524     // store the track
525     track.SetUniqueID(iTrackRef);
526     tmpTrackRefStore->Add(track);
527   }
528   
529   CleanMuonTrackRef(tmpTrackRefStore);
530   
531   delete hit;
532   delete cStore;
533   delete tmpTrackRefStore;
534 }
535
536 //_____________________________________________________________________________
537 void AliMUONRecoCheck::MakeTriggerableTracks()
538 {
539   /// Make triggerable tracks
540  if (!(fRecoTriggerRefStore = AliMUONESDInterface::NewTriggerTrackStore()))
541    return;
542
543   Double_t x, y, z, slopeX, slopeY, pZ;
544   TParticle* particle;
545   TClonesArray* trackRefs;
546   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
547         
548   // loop over simulated tracks
549   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
550     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
551                 
552     // skip empty trackRefs
553     if (nHits < 1) continue;
554                                 
555     AliMUONTriggerTrack track;
556     Int_t hitsOnTrigger = 0;
557     Int_t currCh = -1;
558                 
559     // loop over simulated track hits
560     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
561       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
562                         
563       // skip trackRefs not in MUON
564       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
565
566       // check chamberId of current trackReference
567       Int_t detElemId = trackReference->UserId();
568       Int_t chamberId = detElemId / 100 - 1;
569       if (chamberId < AliMUONConstants::NTrackingCh() || chamberId >= AliMUONConstants::NCh() ) continue;
570                         
571                         
572       if ( hitsOnTrigger == 0 ) {
573         // Get track parameters of current hit
574         x = trackReference->X();
575         y = trackReference->Y();
576         z = trackReference->Z();
577         pZ = trackReference->Pz();
578         slopeX = ( pZ == 0. ) ? 99999. : trackReference->Px() / pZ;
579         slopeY = ( pZ == 0. ) ? 99999. : trackReference->Py() / pZ;
580
581         track.SetX11(x);
582         track.SetY11(y);
583         track.SetZ11(z);
584         track.SetSlopeX(slopeX);
585         track.SetSlopeY(slopeY);
586       }
587
588       if ( currCh != chamberId ) {
589         hitsOnTrigger++;
590         currCh = chamberId;
591       }
592
593     } // loop on hits
594                 
595     if ( hitsOnTrigger >= 3 ){
596       // store the track
597       track.SetUniqueID(iTrackRef);
598       fRecoTriggerRefStore->Add(track);
599     }
600   }
601 }
602
603
604 //_____________________________________________________________________________
605 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
606 {
607   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
608   /// each chamber (one when particle is entering + one when particle is leaving 
609   /// the sensitive volume) 
610   if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
611   
612   Double_t maxGasGap = 1.; // cm 
613   Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
614   Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
615   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
616   if (!cStore) return;
617   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
618   
619   // create iterator
620   TIter next(tmpTrackRefStore->CreateIterator());
621   
622   // loop over tmpTrackRef
623   AliMUONTrack* track;
624   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
625     
626     AliMUONTrack newTrack;
627     
628     // loop over tmpTrackRef's hits
629     Int_t iHit1 = 0;
630     Int_t nTrackHits = track->GetNClusters();
631     while (iHit1 < nTrackHits) {
632       AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
633       
634       // get track parameters at hit1
635       x1  = trackParam1->GetNonBendingCoor();
636       y1  = trackParam1->GetBendingCoor();
637       z1  = trackParam1->GetZ();
638       pX1 = trackParam1->Px();
639       pY1 = trackParam1->Py();
640       pZ1 = trackParam1->Pz();
641       
642       // prepare new track parameters
643       x  = x1;
644       y  = y1;
645       z  = z1;
646       pX = pX1;
647       pY = pY1;
648       pZ = pZ1;
649       
650       // loop over next tmpTrackRef's hits
651       Int_t nCombinedHits = 1;
652       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
653         AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
654         
655         // get z position of hit2
656         z2 = trackParam2->GetZ();
657         
658         // complete new track parameters if hit2 is on the same detection element
659         if ( TMath::Abs(z2-z1) < maxGasGap ) {
660           x  += trackParam2->GetNonBendingCoor();
661           y  += trackParam2->GetBendingCoor();
662           z  += z2;
663           pX += trackParam2->Px();
664           pY += trackParam2->Py();
665           pZ += trackParam2->Pz();
666           nCombinedHits++;
667           iHit1 = iHit2;
668         }
669         
670       }
671       
672       // finalize new track parameters
673       x  /= (Double_t)nCombinedHits;
674       y  /= (Double_t)nCombinedHits;
675       z  /= (Double_t)nCombinedHits;
676       pX /= (Double_t)nCombinedHits;
677       pY /= (Double_t)nCombinedHits;
678       pZ /= (Double_t)nCombinedHits;
679       bendingSlope = 0;
680       nonBendingSlope = 0;
681       inverseBendingMomentum = 0;
682       if (TMath::Abs(pZ) > 0) {
683         bendingSlope = pY/pZ;
684         nonBendingSlope = pX/pZ;
685       }
686       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
687       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
688       inverseBendingMomentum *= trackParam1->GetCharge();
689       
690       // set hit parameters
691       hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
692       hit->SetXYZ(x,y,z);
693       hit->SetErrXY(0.,0.);
694       
695       // set new track parameters at new hit
696       AliMUONTrackParam trackParam;
697       trackParam.SetNonBendingCoor(x);
698       trackParam.SetBendingCoor(y);
699       trackParam.SetZ(z);
700       trackParam.SetBendingSlope(bendingSlope);
701       trackParam.SetNonBendingSlope(nonBendingSlope);
702       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
703       
704       // add track parameters at current hit to the track
705       newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
706       
707       iHit1++;
708     }
709     
710     newTrack.SetUniqueID(track->GetUniqueID());
711     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
712     fTrackRefStore->Add(newTrack);
713     
714   }
715   
716   delete hit;
717   delete cStore;
718 }
719
720 //_____________________________________________________________________________
721 void AliMUONRecoCheck::MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
722 {
723   /// Isolate the reconstructible tracks
724   if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
725   
726   // create iterator on trackRef
727   TIter next(fTrackRefStore->CreateIterator());
728   
729   // loop over trackRef and add reconstructible tracks to fRecoTrackRefStore
730   AliMUONTrack* track;
731   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
732     if (track->IsValid(requestedStationMask, request2ChInSameSt45)) fRecoTrackRefStore->Add(*track);
733   }
734
735 }
736
737 //_____________________________________________________________________________
738 AliMUONTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore,
739                                                     Int_t &nMatchClusters, Bool_t useLabel, Double_t sigmaCut)
740 {
741   /// Return the track from the store matched with the given track (or 0x0) and the number of matched clusters.
742   /// Matching is done by using the MC label of by comparing cluster/TrackRef positions according to the flag "useLabel".
743   /// WARNING: Who match who matters since the matching algorithm uses the *fraction* of matched clusters of the given track
744   
745   AliMUONTrack *matchedTrack = 0x0;
746   nMatchClusters = 0;
747   
748   if (useLabel) { // by using the MC label
749     
750     // get the corresponding simulated track if any
751     Int_t label = track.GetMCLabel();
752     matchedTrack = (AliMUONTrack*) trackStore.FindObject(label);
753     
754     // get the fraction of matched clusters
755     if (matchedTrack) {
756       Int_t nClusters = track.GetNClusters();
757       for (Int_t iCl = 0; iCl < nClusters; iCl++)
758         if (((AliMUONTrackParam*) track.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
759           nMatchClusters++;
760     }
761     
762   } else { // by comparing cluster/TrackRef positions
763     
764     // look for the corresponding simulated track if any
765     TIter next(trackStore.CreateIterator());
766     AliMUONTrack* track2;
767     while ( ( track2 = static_cast<AliMUONTrack*>(next()) ) ) {
768       
769       // check compatibility
770       Int_t n = 0;
771       if (track.Match(*track2, sigmaCut, n)) {
772         matchedTrack = track2;
773         nMatchClusters = n;
774         break;
775       }
776       
777     }
778     
779   }
780   
781   return matchedTrack;
782   
783 }
784
785
786 //_____________________________________________________________________________
787 AliMUONTriggerTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTriggerTrack &track, const AliMUONVTriggerTrackStore &triggerTrackStore,
788                                                            Double_t sigmaCut)
789 {
790   /// Return the trigger track from the store matched with the given track (or 0x0).
791   /// Matching is done by comparing cluster/TrackRef positions.
792   
793   AliMUONTriggerTrack *matchedTrack = 0x0;
794   
795   // look for the corresponding simulated track if any
796   TIter next(triggerTrackStore.CreateIterator());
797   AliMUONTriggerTrack* track2;
798   while ( ( track2 = static_cast<AliMUONTriggerTrack*>(next()) ) ) {
799       
800     // check compatibility
801     if (track.Match(*track2, sigmaCut)) {
802       matchedTrack = track2;
803       break;
804     }
805   }
806   
807   return matchedTrack;
808   
809 }
810