- Apply trigger algorithm to MC tracks when defining a triggerable track (for perform...
[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 "AliMUONDigit.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   fTriggerCircuit = new AliMUONTriggerCircuit(fGeometryTransformer);
194
195   // reset tracker for local trigger to trigger track conversion
196   if ( ! AliMUONESDInterface::GetTracker() )
197     AliMUONESDInterface::ResetTracker();
198   
199   return kTRUE;
200 }
201
202
203 //_____________________________________________________________________________
204 Int_t AliMUONRecoCheck::GetRunNumber()
205 {
206   /// Return the run number of the current ESD event
207   
208   if (fESDEventOwner && fRecoTrackStore == 0x0 && fRecoTriggerTrackStore == 0x0) {
209     if (!fESDTree || fESDTree->GetEvent(fCurrentEvent) <= 0) {
210       AliError(Form("fails to read ESD object for event %d: cannot get the run number",fCurrentEvent));
211       return -1;
212     }
213   }
214   
215   return fESDEvent->GetRunNumber();
216 }
217
218 //_____________________________________________________________________________
219 Int_t AliMUONRecoCheck::NumberOfEvents() const
220 {
221   /// Return the number of events
222   if (fESDEventOwner && fESDTree) return fESDTree->GetEntries();
223   return 0;
224 }
225
226 //_____________________________________________________________________________
227 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event, Bool_t refit)
228 {
229   /// Return a track store containing the reconstructed tracks (converted into 
230   /// MUONTrack objects) for a given event.
231   /// Track parameters at each clusters are computed or not depending on the flag "refit".
232   /// If not, only the track parameters at first cluster are valid.
233   
234   if (!fESDEventOwner) {
235     if (fRecoTrackStore == 0x0) MakeReconstructedTracks(refit);
236     return fRecoTrackStore;
237   }
238
239   if (event != fCurrentEvent) {
240     ResetStores();
241     fCurrentEvent = event;
242   }
243   
244   if (fRecoTrackStore != 0x0) return fRecoTrackStore;
245   else {
246     if (!fESDTree) return 0x0;
247     if (fESDTree->GetEvent(event) <= 0) {
248       AliError(Form("fails to read ESD object for event %d", event));
249       return 0x0;
250     }
251     MakeReconstructedTracks(refit);
252     return fRecoTrackStore;
253   }
254 }
255
256
257 //_____________________________________________________________________________
258 AliMUONVTriggerTrackStore* AliMUONRecoCheck::TriggeredTracks(Int_t event)
259 {
260   /// Return a track store containing the reconstructed trigger tracks (converted into 
261   /// MUONTriggerTrack objects) for a given event.
262         
263   if (!fESDEventOwner) {
264     if (fRecoTriggerTrackStore == 0x0) MakeTriggeredTracks();
265     return fRecoTriggerTrackStore;
266   }
267         
268   if (event != fCurrentEvent) {
269     ResetStores();
270     fCurrentEvent = event;
271   }
272         
273   if (fRecoTriggerTrackStore != 0x0) return fRecoTriggerTrackStore;
274   else {
275     if (!fESDTree) return 0x0;
276     if (fESDTree->GetEvent(event) <= 0) {
277       AliError(Form("fails to read ESD object for event %d", event));
278       return 0x0;
279     }
280     MakeTriggeredTracks();
281     return fRecoTriggerTrackStore;
282   }
283 }
284
285
286 //_____________________________________________________________________________
287 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
288 {
289   /// Return a track store containing the track references (converted into 
290   /// MUONTrack objects) for a given event
291   
292   if (!fESDEventOwner) {
293     if (fTrackRefStore == 0x0) MakeTrackRefs();
294     return fTrackRefStore;
295   }
296
297   if (event != fCurrentEvent) {
298     ResetStores();
299     fCurrentEvent = event;
300   }
301   
302   if (fTrackRefStore != 0x0) return fTrackRefStore;
303   else {
304     if (!fMCEventHandler->LoadEvent(event)) {
305       AliError(Form("fails to read MC objects for event %d", event));
306       return 0x0;
307     }
308     MakeTrackRefs();
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->LoadEvent(event)) {
332       AliError(Form("fails to read MC objects for event %d", event));
333       return 0x0;
334     }
335     MakeTriggerableTracks();
336     return fRecoTriggerRefStore;
337   }
338 }
339
340
341 //_____________________________________________________________________________
342 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event, UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
343 {
344   /// Return a track store containing the reconstructible tracks for a given event,
345   /// according to the mask of requested stations and the minimum number of chambers hit in stations 4 & 5.
346
347   if (!fESDEventOwner) {
348     if (fRecoTrackRefStore == 0x0) {
349       if (TrackRefs(event) == 0x0) return 0x0;
350       MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
351     }
352     return fRecoTrackRefStore;
353   }
354
355   if (event != fCurrentEvent) {
356     ResetStores();
357     fCurrentEvent = event;
358   }
359   
360   if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
361   else {
362     if (TrackRefs(event) == 0x0) return 0x0;
363     MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
364     return fRecoTrackRefStore;
365   }
366 }
367
368
369 //_____________________________________________________________________________
370 void AliMUONRecoCheck::MakeReconstructedTracks(Bool_t refit)
371 {
372   /// Make reconstructed tracks
373   if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
374   
375   // loop over all reconstructed tracks and add them to the store (skip ghosts)
376   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
377   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
378     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
379     if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore, refit);
380   }
381   
382 }
383
384
385 //_____________________________________________________________________________
386 void AliMUONRecoCheck::MakeTriggeredTracks()
387 {
388   /// Make reconstructed trigger tracks
389   if (!(fRecoTriggerTrackStore = AliMUONESDInterface::NewTriggerTrackStore())) return;
390         
391   AliMUONVTriggerStore* tmpTriggerStore = AliMUONESDInterface::NewTriggerStore();
392   if ( ! tmpTriggerStore ) return;
393   
394   // loop over all reconstructed tracks and add them to the store (include ghosts)
395   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
396   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
397     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
398     if (esdTrack->ContainTriggerData()) AliMUONESDInterface::Add(*esdTrack, *tmpTriggerStore);
399   }
400         
401   if ( InitCircuit() ) {
402   
403     AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
404     tracker->EventReconstructTrigger(*fTriggerCircuit, *tmpTriggerStore, *fRecoTriggerTrackStore);
405   }
406   
407   delete tmpTriggerStore;
408 }
409
410 //_____________________________________________________________________________
411 void AliMUONRecoCheck::TriggerToTrack(const AliMUONLocalTrigger& locTrg, AliMUONTriggerTrack& triggerTrack)
412 {
413   /// Make trigger track from local trigger info
414   if ( ! InitCircuit() ) return;
415   AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
416   tracker->TriggerToTrack(*fTriggerCircuit, locTrg, triggerTrack);
417 }
418
419
420 //_____________________________________________________________________________
421 void AliMUONRecoCheck::MakeTrackRefs()
422 {
423   /// Make reconstructible tracks
424   AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
425   if (!tmpTrackRefStore) return;
426   
427   Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
428   TParticle* particle;
429   TClonesArray* trackRefs;
430   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
431   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
432   if (!cStore) return;
433   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
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 cStore;
542   delete tmpTrackRefStore;
543 }
544
545 //_____________________________________________________________________________
546 void AliMUONRecoCheck::MakeTriggerableTracks()
547 {
548   /// Make triggerable tracks
549  if (!(fRecoTriggerRefStore = AliMUONESDInterface::NewTriggerTrackStore()))
550    return;
551
552   Double_t x, y, z, slopeX, slopeY, pZ, xLoc, yLoc, zLoc;
553   TParticle* particle;
554   TClonesArray* trackRefs;
555   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
556         
557   // loop over simulated tracks
558   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
559     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
560                 
561     // skip empty trackRefs
562     if (nHits < 1) continue;
563                                 
564     AliMUONTriggerTrack track;
565     AliMUONDigitStoreV2S digitStore;
566     Int_t hitsOnTrigger = 0;
567     Int_t currCh = -1;
568                 
569     // loop over simulated track hits
570     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
571       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
572                         
573       // skip trackRefs not in MUON
574       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
575
576       // check chamberId of current trackReference
577       Int_t detElemId = trackReference->UserId();
578       Int_t chamberId = detElemId / 100 - 1;
579       if (chamberId < AliMUONConstants::NTrackingCh() || chamberId >= AliMUONConstants::NCh() ) continue;
580
581       // Get track parameters of current hit
582       x = trackReference->X();
583       y = trackReference->Y();
584       z = trackReference->Z();
585       
586       if ( InitTriggerResponse() ) {
587         fGeometryTransformer->Global2Local(detElemId, x, y, z, xLoc, yLoc, zLoc);
588       
589         Int_t nboard = 0;
590         for ( Int_t cath = AliMp::kCath0; cath <= AliMp::kCath1; ++cath )
591         {
592           const AliMpVSegmentation* seg 
593           = AliMpSegmentation::Instance()
594           ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cath));
595         
596           AliMpPad pad = seg->PadByPosition(xLoc,yLoc,kFALSE);
597           Int_t ix = pad.GetIx();
598           Int_t iy = pad.GetIy();
599         
600           if ( !pad.IsValid() ) continue;
601         
602           if ( cath == AliMp::kCath0 ) nboard = pad.GetLocalBoardId(0);
603         
604           AliMUONDigit* digit = new AliMUONDigit(detElemId,nboard,
605                                                  pad.GetLocalBoardChannel(0),cath);
606           digit->SetPadXY(ix,iy);
607           digit->SetCharge(1.);
608           digitStore.Add(*digit,AliMUONVDigitStore::kDeny);
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   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
672   if (!cStore) return;
673   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
674   
675   // create iterator
676   TIter next(tmpTrackRefStore->CreateIterator());
677   
678   // loop over tmpTrackRef
679   AliMUONTrack* track;
680   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
681     
682     AliMUONTrack newTrack;
683     
684     // loop over tmpTrackRef's hits
685     Int_t iHit1 = 0;
686     Int_t nTrackHits = track->GetNClusters();
687     while (iHit1 < nTrackHits) {
688       AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
689       
690       // get track parameters at hit1
691       x1  = trackParam1->GetNonBendingCoor();
692       y1  = trackParam1->GetBendingCoor();
693       z1  = trackParam1->GetZ();
694       pX1 = trackParam1->Px();
695       pY1 = trackParam1->Py();
696       pZ1 = trackParam1->Pz();
697       
698       // prepare new track parameters
699       x  = x1;
700       y  = y1;
701       z  = z1;
702       pX = pX1;
703       pY = pY1;
704       pZ = pZ1;
705       
706       // loop over next tmpTrackRef's hits
707       Int_t nCombinedHits = 1;
708       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
709         AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
710         
711         // get z position of hit2
712         z2 = trackParam2->GetZ();
713         
714         // complete new track parameters if hit2 is on the same detection element
715         if ( TMath::Abs(z2-z1) < maxGasGap ) {
716           x  += trackParam2->GetNonBendingCoor();
717           y  += trackParam2->GetBendingCoor();
718           z  += z2;
719           pX += trackParam2->Px();
720           pY += trackParam2->Py();
721           pZ += trackParam2->Pz();
722           nCombinedHits++;
723           iHit1 = iHit2;
724         }
725         
726       }
727       
728       // finalize new track parameters
729       x  /= (Double_t)nCombinedHits;
730       y  /= (Double_t)nCombinedHits;
731       z  /= (Double_t)nCombinedHits;
732       pX /= (Double_t)nCombinedHits;
733       pY /= (Double_t)nCombinedHits;
734       pZ /= (Double_t)nCombinedHits;
735       bendingSlope = 0;
736       nonBendingSlope = 0;
737       inverseBendingMomentum = 0;
738       if (TMath::Abs(pZ) > 0) {
739         bendingSlope = pY/pZ;
740         nonBendingSlope = pX/pZ;
741       }
742       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
743       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
744       inverseBendingMomentum *= trackParam1->GetCharge();
745       
746       // set hit parameters
747       hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
748       hit->SetXYZ(x,y,z);
749       hit->SetErrXY(0.,0.);
750       
751       // set new track parameters at new hit
752       AliMUONTrackParam trackParam;
753       trackParam.SetNonBendingCoor(x);
754       trackParam.SetBendingCoor(y);
755       trackParam.SetZ(z);
756       trackParam.SetBendingSlope(bendingSlope);
757       trackParam.SetNonBendingSlope(nonBendingSlope);
758       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
759       
760       // add track parameters at current hit to the track
761       newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
762       
763       iHit1++;
764     }
765     
766     newTrack.SetUniqueID(track->GetUniqueID());
767     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
768     fTrackRefStore->Add(newTrack);
769     
770   }
771   
772   delete hit;
773   delete cStore;
774 }
775
776 //_____________________________________________________________________________
777 void AliMUONRecoCheck::MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
778 {
779   /// Isolate the reconstructible tracks
780   if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
781   
782   // create iterator on trackRef
783   TIter next(fTrackRefStore->CreateIterator());
784   
785   // loop over trackRef and add reconstructible tracks to fRecoTrackRefStore
786   AliMUONTrack* track;
787   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
788     if (track->IsValid(requestedStationMask, request2ChInSameSt45)) fRecoTrackRefStore->Add(*track);
789   }
790
791 }
792
793 //_____________________________________________________________________________
794 AliMUONTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore,
795                                                     Int_t &nMatchClusters, Bool_t useLabel, Double_t sigmaCut)
796 {
797   /// Return the track from the store matched with the given track (or 0x0) and the number of matched clusters.
798   /// Matching is done by using the MC label of by comparing cluster/TrackRef positions according to the flag "useLabel".
799   /// WARNING: Who match who matters since the matching algorithm uses the *fraction* of matched clusters of the given track
800   
801   AliMUONTrack *matchedTrack = 0x0;
802   nMatchClusters = 0;
803   
804   if (useLabel) { // by using the MC label
805     
806     // get the corresponding simulated track if any
807     Int_t label = track.GetMCLabel();
808     matchedTrack = (AliMUONTrack*) trackStore.FindObject(label);
809     
810     // get the fraction of matched clusters
811     if (matchedTrack) {
812       Int_t nClusters = track.GetNClusters();
813       for (Int_t iCl = 0; iCl < nClusters; iCl++)
814         if (((AliMUONTrackParam*) track.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
815           nMatchClusters++;
816     }
817     
818   } else { // by comparing cluster/TrackRef positions
819     
820     // look for the corresponding simulated track if any
821     TIter next(trackStore.CreateIterator());
822     AliMUONTrack* track2;
823     while ( ( track2 = static_cast<AliMUONTrack*>(next()) ) ) {
824       
825       // check compatibility
826       Int_t n = 0;
827       if (track.Match(*track2, sigmaCut, n)) {
828         matchedTrack = track2;
829         nMatchClusters = n;
830         break;
831       }
832       
833     }
834     
835   }
836   
837   return matchedTrack;
838   
839 }
840
841
842 //_____________________________________________________________________________
843 AliMUONTriggerTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTriggerTrack &track, const AliMUONVTriggerTrackStore &triggerTrackStore,
844                                                            Double_t sigmaCut)
845 {
846   /// Return the trigger track from the store matched with the given track (or 0x0).
847   /// Matching is done by comparing cluster/TrackRef positions.
848   
849   AliMUONTriggerTrack *matchedTrack = 0x0;
850   
851   // look for the corresponding simulated track if any
852   TIter next(triggerTrackStore.CreateIterator());
853   AliMUONTriggerTrack* track2;
854   while ( ( track2 = static_cast<AliMUONTriggerTrack*>(next()) ) ) {
855     // check compatibility
856     if (track.Match(*track2, sigmaCut)) {
857       matchedTrack = track2;
858       break;
859     }
860   }
861   
862   return matchedTrack;
863   
864 }
865
866
867 //____________________________________________________________________________ 
868 Bool_t AliMUONRecoCheck::InitTriggerResponse()
869 {
870   /// Initialize trigger electronics
871   /// for building of triggerable tracks from MC
872   
873   if ( fTriggerElectronics ) return kTRUE;
874   
875   if ( ! InitGeometryTransformer() ) return kFALSE;
876   
877   if ( ! InitCalibrationData() ) return kFALSE;
878   
879   fTriggerElectronics = new AliMUONTriggerElectronics(fCalibrationData);
880   
881   return kTRUE;
882 }
883
884
885 //____________________________________________________________________________ 
886 Bool_t AliMUONRecoCheck::InitCalibrationData()
887 {
888   /// Initialize calibration data
889   if ( ! fCalibrationData ) {
890     if ( !AliMUONCDB::CheckOCDB() ) return kFALSE;
891     fCalibrationData = new AliMUONCalibrationData(AliCDBManager::Instance()->GetRun());
892   }
893   return kTRUE;
894 }
895
896
897 //____________________________________________________________________________ 
898 Bool_t AliMUONRecoCheck::InitGeometryTransformer()
899 {
900   /// Return calibration data
901   /// (create it if necessary)
902   if ( ! fGeometryTransformer ) {
903     
904     if ( !AliMUONCDB::CheckOCDB() ) return kFALSE;
905     
906     if ( !AliGeomManager::GetGeometry() )
907       AliGeomManager::LoadGeometry();
908     
909     if ( !AliMpDDLStore::Instance(false) )
910       AliMpCDB::LoadDDLStore();
911     
912     fGeometryTransformer = new AliMUONGeometryTransformer();
913     fGeometryTransformer->LoadGeometryData();
914   }
915   
916   return kTRUE;
917 }
918