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