Fix coverity defect
[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,
343                                                             Bool_t request2ChInSameSt45, Bool_t hitInFrontOfPad)
344 {
345   /// Return a track store containing the reconstructible tracks for a given event,
346   /// according to the mask of requested stations and the minimum number of chambers hit in stations 4 & 5.
347
348   if (!fESDEventOwner) {
349     if (fRecoTrackRefStore == 0x0) {
350       if (TrackRefs(event) == 0x0) return 0x0;
351       MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45, hitInFrontOfPad);
352     }
353     return fRecoTrackRefStore;
354   }
355
356   if (event != fCurrentEvent) {
357     ResetStores();
358     fCurrentEvent = event;
359   }
360   
361   if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
362   else {
363     if (TrackRefs(event) == 0x0) return 0x0;
364     MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45, hitInFrontOfPad);
365     return fRecoTrackRefStore;
366   }
367 }
368
369
370 //_____________________________________________________________________________
371 void AliMUONRecoCheck::MakeReconstructedTracks(Bool_t refit)
372 {
373   /// Make reconstructed tracks
374   if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
375   
376   // loop over all reconstructed tracks and add them to the store (skip ghosts)
377   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
378   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
379     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
380     if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore, refit);
381   }
382   
383 }
384
385
386 //_____________________________________________________________________________
387 void AliMUONRecoCheck::MakeTriggeredTracks()
388 {
389   /// Make reconstructed trigger tracks
390   if (!(fRecoTriggerTrackStore = AliMUONESDInterface::NewTriggerTrackStore())) return;
391         
392   AliMUONVTriggerStore* tmpTriggerStore = AliMUONESDInterface::NewTriggerStore();
393   if ( ! tmpTriggerStore ) return;
394   
395   // loop over all reconstructed tracks and add them to the store (include ghosts)
396   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
397   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
398     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
399     if (esdTrack->ContainTriggerData()) AliMUONESDInterface::Add(*esdTrack, *tmpTriggerStore);
400   }
401         
402   if ( InitCircuit() ) {
403   
404     AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
405     tracker->EventReconstructTrigger(*fTriggerCircuit, *tmpTriggerStore, *fRecoTriggerTrackStore);
406   }
407   
408   delete tmpTriggerStore;
409 }
410
411 //_____________________________________________________________________________
412 void AliMUONRecoCheck::TriggerToTrack(const AliMUONLocalTrigger& locTrg, AliMUONTriggerTrack& triggerTrack)
413 {
414   /// Make trigger track from local trigger info
415   if ( ! InitCircuit() ) return;
416   AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
417   tracker->TriggerToTrack(*fTriggerCircuit, locTrg, triggerTrack);
418 }
419
420
421 //_____________________________________________________________________________
422 void AliMUONRecoCheck::MakeTrackRefs()
423 {
424   /// Make reconstructible tracks
425   AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
426   if (!tmpTrackRefStore) return;
427   
428   Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
429   TParticle* particle;
430   TClonesArray* trackRefs;
431   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
432   AliMUONVCluster* hit = AliMUONESDInterface::NewCluster();
433   
434   // loop over simulated tracks
435   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
436     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
437     
438     // skip empty trackRefs
439     if (nHits < 1) continue;
440     
441     // get the particle charge for further calculation
442     TParticlePDG* ppdg = particle->GetPDG();
443     Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
444     
445     AliMUONTrack track;
446     
447     // loop over simulated track hits
448     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
449       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
450       
451       // skip trackRefs not in MUON
452       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
453     
454       // Get track parameters of current hit
455       x = trackReference->X();
456       y = trackReference->Y();
457       z = trackReference->Z();
458       pX = trackReference->Px();
459       pY = trackReference->Py();
460       pZ = trackReference->Pz();
461       
462       // check chamberId of current trackReference
463       Int_t detElemId = trackReference->UserId();
464       Int_t chamberId = detElemId / 100 - 1;
465       if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
466       
467       // set hit parameters
468       hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
469       hit->SetXYZ(x,y,z);
470       hit->SetErrXY(0.,0.);
471       
472       // compute track parameters at hit
473       bendingSlope = 0;
474       nonBendingSlope = 0;
475       inverseBendingMomentum = 0;
476       if (TMath::Abs(pZ) > 0) {
477         bendingSlope = pY/pZ;
478         nonBendingSlope = pX/pZ;
479       }
480       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
481       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
482       inverseBendingMomentum *= charge;
483       
484       // set track parameters at hit
485       AliMUONTrackParam trackParam;
486       trackParam.SetNonBendingCoor(x);
487       trackParam.SetBendingCoor(y);
488       trackParam.SetZ(z);
489       trackParam.SetBendingSlope(bendingSlope);
490       trackParam.SetNonBendingSlope(nonBendingSlope);
491       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
492       
493       // add track parameters at current hit to the track
494       track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
495     }
496     
497     // if none of the track hits was in MUON, goto the next track
498     if (track.GetNClusters() < 1) continue;
499     
500     // get track parameters at particle's vertex
501     x = particle->Vx();
502     y = particle->Vy();
503     z = particle->Vz();
504     pX = particle->Px();
505     pY = particle->Py();
506     pZ = particle->Pz();
507     
508     // compute rest of track parameters at particle's vertex
509     bendingSlope = 0;
510     nonBendingSlope = 0;
511     inverseBendingMomentum = 0;
512     if (TMath::Abs(pZ) > 0) {
513       bendingSlope = pY/pZ;
514       nonBendingSlope = pX/pZ;
515     }
516     Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
517     if (pYZ >0) inverseBendingMomentum = 1/pYZ;
518     inverseBendingMomentum *= charge;
519     
520     // set track parameters at particle's vertex
521     AliMUONTrackParam trackParamAtVertex;
522     trackParamAtVertex.SetNonBendingCoor(x);
523     trackParamAtVertex.SetBendingCoor(y);
524     trackParamAtVertex.SetZ(z);
525     trackParamAtVertex.SetBendingSlope(bendingSlope);
526     trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
527     trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
528     
529     // add track parameters at vertex
530     track.SetTrackParamAtVertex(&trackParamAtVertex);
531     
532     // store the track
533     track.SetUniqueID(iTrackRef);
534     tmpTrackRefStore->Add(track);
535   }
536   
537   CleanMuonTrackRef(tmpTrackRefStore);
538   
539   delete hit;
540   delete tmpTrackRefStore;
541 }
542
543 //_____________________________________________________________________________
544 void AliMUONRecoCheck::MakeTriggerableTracks()
545 {
546   /// Make triggerable tracks
547  if (!(fRecoTriggerRefStore = AliMUONESDInterface::NewTriggerTrackStore()))
548    return;
549
550   Double_t x, y, z, slopeX, slopeY, pZ, xLoc, yLoc, zLoc;
551   TParticle* particle;
552   TClonesArray* trackRefs;
553   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
554         
555   // loop over simulated tracks
556   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
557     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
558                 
559     // skip empty trackRefs
560     if (nHits < 1) continue;
561                                 
562     AliMUONTriggerTrack track;
563     AliMUONDigitStoreV2S digitStore;
564     Int_t hitsOnTrigger = 0;
565     Int_t currCh = -1;
566                 
567     // loop over simulated track hits
568     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
569       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
570                         
571       // skip trackRefs not in MUON
572       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
573
574       // check chamberId of current trackReference
575       Int_t detElemId = trackReference->UserId();
576       Int_t chamberId = detElemId / 100 - 1;
577       if (chamberId < AliMUONConstants::NTrackingCh() || chamberId >= AliMUONConstants::NCh() ) continue;
578
579       // Get track parameters of current hit
580       x = trackReference->X();
581       y = trackReference->Y();
582       z = trackReference->Z();
583       
584       if ( InitTriggerResponse() ) {
585         fGeometryTransformer->Global2Local(detElemId, x, y, z, xLoc, yLoc, zLoc);
586       
587         Int_t nboard = 0;
588         for ( Int_t cath = AliMp::kCath0; cath <= AliMp::kCath1; ++cath )
589         {
590           const AliMpVSegmentation* seg 
591           = AliMpSegmentation::Instance()
592           ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cath));
593         
594           AliMpPad pad = seg->PadByPosition(xLoc,yLoc,kFALSE);
595           Int_t ix = pad.GetIx();
596           Int_t iy = pad.GetIy();
597         
598           if ( !pad.IsValid() ) continue;
599         
600           if ( cath == AliMp::kCath0 ) nboard = pad.GetLocalBoardId(0);
601         
602           AliMUONDigit* digit = new AliMUONDigit(detElemId,nboard,
603                                                  pad.GetLocalBoardChannel(0),cath);
604           digit->SetPadXY(ix,iy);
605           digit->SetCharge(1.);
606           digitStore.Add(*digit,AliMUONVDigitStore::kDeny);
607         }
608       }
609     
610       if ( hitsOnTrigger == 0 ) {
611         pZ = trackReference->Pz();
612         slopeX = ( pZ == 0. ) ? 99999. : trackReference->Px() / pZ;
613         slopeY = ( pZ == 0. ) ? 99999. : trackReference->Py() / pZ;
614
615         track.SetX11(x);
616         track.SetY11(y);
617         track.SetZ11(z);
618         track.SetSlopeX(slopeX);
619         track.SetSlopeY(slopeY);
620       }
621
622       if ( currCh != chamberId ) {
623         hitsOnTrigger++;
624         currCh = chamberId;
625       }
626
627     } // loop on hits
628     
629     if ( hitsOnTrigger < 3 ) continue;
630
631     // Check if the track passes the trigger algorithm
632     if ( InitTriggerResponse() ) {
633       AliMUONTriggerStoreV1 triggerStore;
634       fTriggerElectronics->Digits2Trigger(digitStore,triggerStore);
635     
636       TIter next(triggerStore.CreateIterator());
637       AliMUONLocalTrigger* locTrg(0x0);
638     
639       Int_t ptCutLevel = 0;
640     
641       while ( ( locTrg = static_cast<AliMUONLocalTrigger*>(next()) ) )
642       {
643         if ( locTrg->IsTrigX() && locTrg->IsTrigY() ) 
644         {
645           ptCutLevel = TMath::Max(ptCutLevel, 1);
646           if ( locTrg->LoHpt() ) ptCutLevel = TMath::Max(ptCutLevel, 3);
647           else if ( locTrg->LoLpt() ) ptCutLevel = TMath::Max(ptCutLevel, 2);
648         } // board is fired 
649       } // end of loop on Local Trigger
650       track.SetPtCutLevel(ptCutLevel);
651     }
652     track.SetUniqueID(iTrackRef);
653     fRecoTriggerRefStore->Add(track);
654   }
655 }
656
657
658 //_____________________________________________________________________________
659 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
660 {
661   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
662   /// each chamber (one when particle is entering + one when particle is leaving 
663   /// the sensitive volume) 
664   if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
665   
666   Double_t maxGasGap = 1.; // cm 
667   Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
668   Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
669   AliMUONVCluster* hit = AliMUONESDInterface::NewCluster();
670   
671   // create iterator
672   TIter next(tmpTrackRefStore->CreateIterator());
673   
674   // loop over tmpTrackRef
675   AliMUONTrack* track;
676   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
677     
678     AliMUONTrack newTrack;
679     
680     // loop over tmpTrackRef's hits
681     Int_t iHit1 = 0;
682     Int_t nTrackHits = track->GetNClusters();
683     while (iHit1 < nTrackHits) {
684       AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
685       
686       // get track parameters at hit1
687       x1  = trackParam1->GetNonBendingCoor();
688       y1  = trackParam1->GetBendingCoor();
689       z1  = trackParam1->GetZ();
690       pX1 = trackParam1->Px();
691       pY1 = trackParam1->Py();
692       pZ1 = trackParam1->Pz();
693       
694       // prepare new track parameters
695       x  = x1;
696       y  = y1;
697       z  = z1;
698       pX = pX1;
699       pY = pY1;
700       pZ = pZ1;
701       
702       // loop over next tmpTrackRef's hits
703       Int_t nCombinedHits = 1;
704       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
705         AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
706         
707         // get z position of hit2
708         z2 = trackParam2->GetZ();
709         
710         // complete new track parameters if hit2 is on the same detection element
711         if ( TMath::Abs(z2-z1) < maxGasGap ) {
712           x  += trackParam2->GetNonBendingCoor();
713           y  += trackParam2->GetBendingCoor();
714           z  += z2;
715           pX += trackParam2->Px();
716           pY += trackParam2->Py();
717           pZ += trackParam2->Pz();
718           nCombinedHits++;
719           iHit1 = iHit2;
720         }
721         
722       }
723       
724       // finalize new track parameters
725       x  /= (Double_t)nCombinedHits;
726       y  /= (Double_t)nCombinedHits;
727       z  /= (Double_t)nCombinedHits;
728       pX /= (Double_t)nCombinedHits;
729       pY /= (Double_t)nCombinedHits;
730       pZ /= (Double_t)nCombinedHits;
731       bendingSlope = 0;
732       nonBendingSlope = 0;
733       inverseBendingMomentum = 0;
734       if (TMath::Abs(pZ) > 0) {
735         bendingSlope = pY/pZ;
736         nonBendingSlope = pX/pZ;
737       }
738       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
739       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
740       inverseBendingMomentum *= trackParam1->GetCharge();
741       
742       // set hit parameters
743       hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
744       hit->SetXYZ(x,y,z);
745       hit->SetErrXY(0.,0.);
746       
747       // set new track parameters at new hit
748       AliMUONTrackParam trackParam;
749       trackParam.SetNonBendingCoor(x);
750       trackParam.SetBendingCoor(y);
751       trackParam.SetZ(z);
752       trackParam.SetBendingSlope(bendingSlope);
753       trackParam.SetNonBendingSlope(nonBendingSlope);
754       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
755       
756       // add track parameters at current hit to the track
757       newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
758       
759       iHit1++;
760     }
761     
762     newTrack.SetUniqueID(track->GetUniqueID());
763     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
764     fTrackRefStore->Add(newTrack);
765     
766   }
767   
768   delete hit;
769 }
770
771 //_____________________________________________________________________________
772 void AliMUONRecoCheck::MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45,
773                                                  Bool_t hitInFrontOfPad)
774 {
775   /// Isolate the reconstructible tracks
776   if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
777   
778   // need geometry transformer and segmentation to check if trackrefs are located in front of pad(s)
779   if (hitInFrontOfPad) {
780     if (!InitGeometryTransformer()) return;
781     if (!AliMpSegmentation::Instance(kFALSE) && !AliMUONCDB::LoadMapping(kTRUE)) return;
782   }
783   
784   // create iterator on trackRef
785   TIter next(fTrackRefStore->CreateIterator());
786   
787   // loop over trackRef and add reconstructible tracks to fRecoTrackRefStore
788   AliMUONTrack* track;
789   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
790     
791     // check if the track is valid as it is
792     if (track->IsValid(requestedStationMask, request2ChInSameSt45)) {
793       
794       if (hitInFrontOfPad) {
795         
796         AliMUONTrack trackTmp(*track);
797         
798         // loop over clusters and remove those which are not in front of pad(s)
799         AliMUONTrackParam *param = static_cast<AliMUONTrackParam*>(trackTmp.GetTrackParamAtCluster()->First());
800         while (param) {
801           AliMUONTrackParam *nextParam = static_cast<AliMUONTrackParam*>(trackTmp.GetTrackParamAtCluster()->After(param));
802           if (!IsHitInFrontOfPad(param)) trackTmp.RemoveTrackParamAtCluster(param);
803           param = nextParam;
804         }
805         
806         // add the track if it is still valid
807         if (trackTmp.IsValid(requestedStationMask, request2ChInSameSt45)) fRecoTrackRefStore->Add(trackTmp);
808         
809       } else {
810         
811         // add the track
812         fRecoTrackRefStore->Add(*track);
813         
814       }
815       
816     }
817       
818   }
819
820 }
821
822 //_____________________________________________________________________________
823 AliMUONTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore,
824                                                     Int_t &nMatchClusters, Bool_t useLabel, Double_t sigmaCut)
825 {
826   /// Return the track from the store matched with the given track (or 0x0) and the number of matched clusters.
827   /// Matching is done by using the MC label of by comparing cluster/TrackRef positions according to the flag "useLabel".
828   /// WARNING: Who match who matters since the matching algorithm uses the *fraction* of matched clusters of the given track
829   
830   AliMUONTrack *matchedTrack = 0x0;
831   nMatchClusters = 0;
832   
833   if (useLabel) { // by using the MC label
834     
835     // get the corresponding simulated track if any
836     Int_t label = track.GetMCLabel();
837     matchedTrack = (AliMUONTrack*) trackStore.FindObject(label);
838     
839     // get the fraction of matched clusters
840     if (matchedTrack) {
841       Int_t nClusters = track.GetNClusters();
842       for (Int_t iCl = 0; iCl < nClusters; iCl++)
843         if (((AliMUONTrackParam*) track.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
844           nMatchClusters++;
845     }
846     
847   } else { // by comparing cluster/TrackRef positions
848     
849     // look for the corresponding simulated track if any
850     TIter next(trackStore.CreateIterator());
851     AliMUONTrack* track2;
852     while ( ( track2 = static_cast<AliMUONTrack*>(next()) ) ) {
853       
854       // check compatibility
855       Int_t n = 0;
856       if (track.Match(*track2, sigmaCut, n)) {
857         matchedTrack = track2;
858         nMatchClusters = n;
859         break;
860       }
861       
862     }
863     
864   }
865   
866   return matchedTrack;
867   
868 }
869
870
871 //_____________________________________________________________________________
872 AliMUONTriggerTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTriggerTrack &track, const AliMUONVTriggerTrackStore &triggerTrackStore,
873                                                            Double_t sigmaCut)
874 {
875   /// Return the trigger track from the store matched with the given track (or 0x0).
876   /// Matching is done by comparing cluster/TrackRef positions.
877   
878   AliMUONTriggerTrack *matchedTrack = 0x0;
879   
880   // look for the corresponding simulated track if any
881   TIter next(triggerTrackStore.CreateIterator());
882   AliMUONTriggerTrack* track2;
883   while ( ( track2 = static_cast<AliMUONTriggerTrack*>(next()) ) ) {
884     // check compatibility
885     if (track.Match(*track2, sigmaCut)) {
886       matchedTrack = track2;
887       break;
888     }
889   }
890   
891   return matchedTrack;
892   
893 }
894
895
896 //____________________________________________________________________________ 
897 Bool_t AliMUONRecoCheck::InitTriggerResponse()
898 {
899   /// Initialize trigger electronics
900   /// for building of triggerable tracks from MC
901   
902   if ( fTriggerElectronics ) return kTRUE;
903   
904   if ( ! InitGeometryTransformer() ) return kFALSE;
905   
906   if ( ! AliMpDDLStore::Instance(false) ) AliMpCDB::LoadDDLStore();
907   
908   if ( ! InitCalibrationData() ) return kFALSE;
909   
910   fTriggerElectronics = new AliMUONTriggerElectronics(fCalibrationData);
911   
912   return kTRUE;
913 }
914
915
916 //____________________________________________________________________________ 
917 Bool_t AliMUONRecoCheck::InitCalibrationData()
918 {
919   /// Initialize calibration data
920   if ( ! fCalibrationData ) {
921     if ( !AliMUONCDB::CheckOCDB() ) return kFALSE;
922     fCalibrationData = new AliMUONCalibrationData(AliCDBManager::Instance()->GetRun());
923   }
924   return kTRUE;
925 }
926
927
928 //____________________________________________________________________________ 
929 Bool_t AliMUONRecoCheck::InitGeometryTransformer()
930 {
931   /// Return geometry transformer
932   /// (create it if necessary)
933   if ( ! fGeometryTransformer ) {
934     
935     if ( !AliMUONCDB::CheckOCDB() ) return kFALSE;
936     
937     if (!AliGeomManager::GetGeometry()) {
938       AliGeomManager::LoadGeometry();
939       if (!AliGeomManager::GetGeometry()) return kFALSE;  
940       if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return kFALSE;
941     }
942     
943     fGeometryTransformer = new AliMUONGeometryTransformer();
944     fGeometryTransformer->LoadGeometryData();
945   }
946   
947   return kTRUE;
948 }
949
950 //________________________________________________________________________
951 Bool_t AliMUONRecoCheck::IsHitInFrontOfPad(AliMUONTrackParam *param) const
952 {
953   /// Return kTRUE if the hit is located in front of at least 1 pad
954   
955   Int_t deId = param->GetClusterPtr()->GetDetElemId();
956   
957   Double_t xL, yL, zL;
958   fGeometryTransformer->Global2Local(deId, param->GetNonBendingCoor(), param->GetBendingCoor(), param->GetZ(), xL, yL, zL);
959   
960   const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, AliMp::kCath0);
961   const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, AliMp::kCath1);
962   if (!seg1 || !seg2) return kFALSE;
963   
964   AliMpPad pad1 = seg1->PadByPosition(xL, yL, kFALSE);
965   AliMpPad pad2 = seg2->PadByPosition(xL, yL, kFALSE);
966   
967   return (pad1.IsValid() || pad2.IsValid());
968 }
969