Changes for #93172: Improve trigger chamber efficiency calculation
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.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 AliMUONTrack
20 //-------------------
21 // Reconstructed track in ALICE dimuon spectrometer
22 //-----------------------------------------------------------------------------
23
24 #include "AliMUONTrack.h"
25
26 #include "AliMUONReconstructor.h"
27 #include "AliMUONVCluster.h"
28 #include "AliMUONVClusterStore.h"
29 #include "AliMUONObjectPair.h"
30 #include "AliMUONTrackExtrap.h"
31 #include "AliMUONConstants.h"
32 #include "AliMUONTrackParam.h"
33
34 #include "AliLog.h"
35
36 #include <TMath.h>
37
38 #include <Riostream.h>
39
40 /// \cond CLASSIMP
41 ClassImp(AliMUONTrack) // Class implementation in ROOT context
42 /// \endcond
43
44
45 const Double_t AliMUONTrack::fgkMaxChi2 = 1.e10; ///< maximum chi2 above which the track can be considered as abnormal
46
47
48 //__________________________________________________________________________
49 AliMUONTrack::AliMUONTrack()
50   : TObject(),
51     fTrackParamAtCluster(0x0),
52     fFitWithVertex(kFALSE),
53     fVertexErrXY2(),
54     fFitWithMCS(kFALSE),
55     fClusterWeightsNonBending(0x0),
56     fClusterWeightsBending(0x0),
57     fGlobalChi2(-1.),
58     fImproved(kFALSE),
59     fMatchTrigger(-1),
60     fChi2MatchTrigger(0.),
61     fTrackID(-1),
62     fTrackParamAtVertex(0x0),
63     fHitsPatternInTrigCh(0),
64     fHitsPatternInTrigChTrk(0),
65     fLocalTrigger(0),
66     fConnected(kFALSE)
67 {
68   /// Default constructor
69   fVertexErrXY2[0] = 0.;
70   fVertexErrXY2[1] = 0.;
71 }
72
73   //__________________________________________________________________________
74 AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment, Double_t bendingVertexDispersion)
75   : TObject(),
76     fTrackParamAtCluster(new TObjArray(20)),
77     fFitWithVertex(kFALSE),
78     fVertexErrXY2(),
79     fFitWithMCS(kFALSE),
80     fClusterWeightsNonBending(0x0),
81     fClusterWeightsBending(0x0),
82     fGlobalChi2(0.),
83     fImproved(kFALSE),
84     fMatchTrigger(-1),
85     fChi2MatchTrigger(0.),
86     fTrackID(-1),
87     fTrackParamAtVertex(0x0),
88     fHitsPatternInTrigCh(0),
89     fHitsPatternInTrigChTrk(0),
90     fLocalTrigger(0),
91     fConnected(kFALSE)
92 {
93   /// Constructor from two clusters
94   
95   fTrackParamAtCluster->SetOwner(kTRUE);
96   
97   fVertexErrXY2[0] = 0.;
98   fVertexErrXY2[1] = 0.;
99   
100   // Pointers to clusters from the segment
101   AliMUONVCluster* firstCluster = (AliMUONVCluster*) segment->First();
102   AliMUONVCluster* lastCluster = (AliMUONVCluster*) segment->Second();
103   
104   // Compute track parameters
105   Double_t z1 = firstCluster->GetZ();
106   Double_t z2 = lastCluster->GetZ();
107   Double_t dZ = z1 - z2;
108   // Non bending plane
109   Double_t nonBendingCoor1 = firstCluster->GetX();
110   Double_t nonBendingCoor2 = lastCluster->GetX();
111   Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ;
112   // Bending plane
113   Double_t bendingCoor1 = firstCluster->GetY();
114   Double_t bendingCoor2 = lastCluster->GetY();
115   Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ;
116   // Inverse bending momentum
117   Double_t bendingImpact = bendingCoor1 - z1 * bendingSlope;
118   Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
119   
120   // Set track parameters at first cluster
121   AliMUONTrackParam trackParamAtFirstCluster;
122   trackParamAtFirstCluster.SetZ(z1);
123   trackParamAtFirstCluster.SetNonBendingCoor(nonBendingCoor1);
124   trackParamAtFirstCluster.SetNonBendingSlope(nonBendingSlope);
125   trackParamAtFirstCluster.SetBendingCoor(bendingCoor1);
126   trackParamAtFirstCluster.SetBendingSlope(bendingSlope);
127   trackParamAtFirstCluster.SetInverseBendingMomentum(inverseBendingMomentum);
128   
129   // Set track parameters at last cluster
130   AliMUONTrackParam trackParamAtLastCluster;
131   trackParamAtLastCluster.SetZ(z2);
132   trackParamAtLastCluster.SetNonBendingCoor(nonBendingCoor2);
133   trackParamAtLastCluster.SetNonBendingSlope(nonBendingSlope);
134   trackParamAtLastCluster.SetBendingCoor(bendingCoor2);
135   trackParamAtLastCluster.SetBendingSlope(bendingSlope);
136   trackParamAtLastCluster.SetInverseBendingMomentum(inverseBendingMomentum);
137   
138   // Compute and set track parameters covariances at first cluster
139   TMatrixD paramCov(5,5);
140   paramCov.Zero();
141   // Non bending plane
142   paramCov(0,0) = firstCluster->GetErrX2();
143   paramCov(0,1) = firstCluster->GetErrX2() / dZ;
144   paramCov(1,0) = paramCov(0,1);
145   paramCov(1,1) = ( firstCluster->GetErrX2() + lastCluster->GetErrX2() ) / dZ / dZ;
146   // Bending plane
147   paramCov(2,2) = firstCluster->GetErrY2();
148   paramCov(2,3) = firstCluster->GetErrY2() / dZ;
149   paramCov(3,2) = paramCov(2,3);
150   paramCov(3,3) = ( firstCluster->GetErrY2() + lastCluster->GetErrY2() ) / dZ / dZ;
151   // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
152   if (AliMUONTrackExtrap::IsFieldON()) {
153     paramCov(4,4) = ( ( bendingVertexDispersion*bendingVertexDispersion +
154                        (z1 * z1 * lastCluster->GetErrY2() + z2 * z2 * firstCluster->GetErrY2()) / dZ / dZ) /
155                      bendingImpact / bendingImpact + 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum ;
156     paramCov(2,4) = - z2 * firstCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
157     paramCov(4,2) = paramCov(2,4);
158     paramCov(3,4) = - (z1 * lastCluster->GetErrY2() + z2 * firstCluster->GetErrY2()) * inverseBendingMomentum / bendingImpact / dZ / dZ;
159     paramCov(4,3) = paramCov(3,4);
160   } else paramCov(4,4) = inverseBendingMomentum*inverseBendingMomentum;
161   trackParamAtFirstCluster.SetCovariances(paramCov);
162   
163   // Compute and set track parameters covariances at last cluster
164   // Non bending plane
165   paramCov(0,0) = lastCluster->GetErrX2();
166   paramCov(0,1) = - lastCluster->GetErrX2() / dZ;
167   paramCov(1,0) = paramCov(0,1);
168   // Bending plane
169   paramCov(2,2) = lastCluster->GetErrY2();
170   paramCov(2,3) = - lastCluster->GetErrY2() / dZ;
171   paramCov(3,2) = paramCov(2,3);
172   // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
173   if (AliMUONTrackExtrap::IsFieldON()) {
174     paramCov(2,4) = z1 * lastCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
175     paramCov(4,2) = paramCov(2,4);
176   }
177   trackParamAtLastCluster.SetCovariances(paramCov);
178   
179   // Add track parameters at clusters
180   AddTrackParamAtCluster(trackParamAtFirstCluster,*firstCluster);
181   AddTrackParamAtCluster(trackParamAtLastCluster,*lastCluster);
182   
183 }
184
185 //__________________________________________________________________________
186 AliMUONTrack::AliMUONTrack(const AliMUONTrack& track)
187   : TObject(track),
188     fTrackParamAtCluster(0x0),
189     fFitWithVertex(track.fFitWithVertex),
190     fVertexErrXY2(),
191     fFitWithMCS(track.fFitWithMCS),
192     fClusterWeightsNonBending(0x0),
193     fClusterWeightsBending(0x0),
194     fGlobalChi2(track.fGlobalChi2),
195     fImproved(track.fImproved),
196     fMatchTrigger(track.fMatchTrigger),
197     fChi2MatchTrigger(track.fChi2MatchTrigger),
198     fTrackID(track.fTrackID),
199     fTrackParamAtVertex(0x0),
200     fHitsPatternInTrigCh(track.fHitsPatternInTrigCh),
201     fHitsPatternInTrigChTrk(track.fHitsPatternInTrigChTrk),
202     fLocalTrigger(track.fLocalTrigger),
203     fConnected(track.fConnected)
204 {
205   ///copy constructor
206   
207   // necessary to make a copy of the objects and not only the pointers in TObjArray.
208   if (track.fTrackParamAtCluster) {
209     fTrackParamAtCluster = new TObjArray(track.fTrackParamAtCluster->GetSize());
210     fTrackParamAtCluster->SetOwner(kTRUE);
211     for (Int_t i = 0; i < track.GetNClusters(); i++)
212       fTrackParamAtCluster->AddLast(new AliMUONTrackParam(*static_cast<AliMUONTrackParam*>(track.fTrackParamAtCluster->UncheckedAt(i))));
213   }
214   
215   // copy vertex resolution square used during the tracking procedure
216   fVertexErrXY2[0] = track.fVertexErrXY2[0];
217   fVertexErrXY2[1] = track.fVertexErrXY2[1];
218   
219   // copy cluster weights matrices if any
220   if (track.fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
221   if (track.fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
222   
223   // copy track parameters at vertex if any
224   if (track.fTrackParamAtVertex) fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
225   
226 }
227
228   //__________________________________________________________________________
229 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
230 {
231   /// Asignment operator
232   // check assignement to self
233   if (this == &track)
234     return *this;
235
236   // base class assignement
237   TObject::operator=(track);
238   
239   // clear memory
240   Clear();
241   
242   // necessary to make a copy of the objects and not only the pointers in TObjArray
243   if (track.fTrackParamAtCluster) {
244     fTrackParamAtCluster = new TObjArray(track.fTrackParamAtCluster->GetSize());
245     fTrackParamAtCluster->SetOwner(kTRUE);
246     for (Int_t i = 0; i < track.GetNClusters(); i++)
247       fTrackParamAtCluster->AddLast(new AliMUONTrackParam(*static_cast<AliMUONTrackParam*>(track.fTrackParamAtCluster->UncheckedAt(i))));
248   }
249   
250   // copy cluster weights matrix if any
251   if (track.fClusterWeightsNonBending) {
252     if (fClusterWeightsNonBending) {
253       fClusterWeightsNonBending->ResizeTo(*(track.fClusterWeightsNonBending));
254       *fClusterWeightsNonBending = *(track.fClusterWeightsNonBending);
255     } else fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
256   }
257   
258   // copy cluster weights matrix if any
259   if (track.fClusterWeightsBending) {
260     if (fClusterWeightsBending) {
261       fClusterWeightsBending->ResizeTo(*(track.fClusterWeightsBending));
262       *fClusterWeightsBending = *(track.fClusterWeightsBending);
263     } else fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
264   }
265   
266   // copy track parameters at vertex if any
267   if (track.fTrackParamAtVertex) {
268     if (fTrackParamAtVertex) *fTrackParamAtVertex = *(track.fTrackParamAtVertex);
269     else fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
270   }
271   
272   fFitWithVertex      =  track.fFitWithVertex;
273   fVertexErrXY2[0]    =  track.fVertexErrXY2[0];
274   fVertexErrXY2[1]    =  track.fVertexErrXY2[1];
275   fFitWithMCS         =  track.fFitWithMCS;
276   fGlobalChi2         =  track.fGlobalChi2;
277   fImproved           =  track.fImproved;
278   fMatchTrigger       =  track.fMatchTrigger;
279   fChi2MatchTrigger   =  track.fChi2MatchTrigger;
280   fTrackID            =  track.fTrackID; 
281   fHitsPatternInTrigCh = track.fHitsPatternInTrigCh;
282   fHitsPatternInTrigChTrk = track.fHitsPatternInTrigChTrk;
283   fLocalTrigger        = track.fLocalTrigger;
284   fConnected          =  track.fConnected;
285
286   return *this;
287 }
288
289   //__________________________________________________________________________
290 AliMUONTrack::~AliMUONTrack()
291 {
292   /// Destructor
293   delete fTrackParamAtCluster;
294   delete fClusterWeightsNonBending;
295   delete fClusterWeightsBending;
296   delete fTrackParamAtVertex;
297 }
298
299   //__________________________________________________________________________
300 void AliMUONTrack::Clear(Option_t* /*opt*/)
301 {
302   /// Clear arrays
303   delete fTrackParamAtCluster; fTrackParamAtCluster = 0x0;
304   delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
305   delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
306   delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
307 }
308
309   //__________________________________________________________________________
310 void AliMUONTrack::Reset()
311 {
312   /// Reset to default values
313   SetUniqueID(0);
314   fFitWithVertex = kFALSE;
315   fVertexErrXY2[0] = 0.;
316   fVertexErrXY2[1] = 0.;
317   fFitWithMCS = kFALSE;
318   fGlobalChi2 = -1.;
319   fImproved = kFALSE;
320   fMatchTrigger = -1;
321   fChi2MatchTrigger = 0.;
322   fTrackID = -1;
323   fHitsPatternInTrigCh = 0;
324   fHitsPatternInTrigChTrk = 0;
325   fLocalTrigger = 0;
326   fConnected = kFALSE;
327   delete fTrackParamAtCluster; fTrackParamAtCluster = 0x0;
328   delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
329   delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
330   delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
331 }
332
333   //__________________________________________________________________________
334 TObjArray* AliMUONTrack::GetTrackParamAtCluster() const
335 {
336   /// return array of track parameters at cluster (create it if needed)
337   if (!fTrackParamAtCluster) {
338     fTrackParamAtCluster = new TObjArray(20);
339     fTrackParamAtCluster->SetOwner(kTRUE);
340   }
341   return fTrackParamAtCluster;
342 }
343
344   //__________________________________________________________________________
345 void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy)
346 {
347   /// Copy given track parameters into a new TrackParamAtCluster
348   /// Link parameters with the associated cluster
349   /// If copy=kTRUE: the cluster is copied then passed the trackParam which become its owner 
350   ///     otherwise: make sure to do not delete the cluster until it is used by the track
351   
352   // check chamber ID of the associated cluster
353   if (cluster.GetChamberId() < 0 || cluster.GetChamberId() > AliMUONConstants::NTrackingCh()) {
354     AliError(Form("Chamber ID of the associated cluster is not valid (ChamberId=%d)",cluster.GetChamberId()));
355     return;
356   }
357   
358   // check whether track parameters are given at the correct cluster z position
359   if (TMath::Abs(cluster.GetZ() - trackParam.GetZ())>1.e-5) {   // AU
360     AliError("track parameters are given at a different z position than the one of the associated cluster");
361     return;
362   }
363   
364   // add parameters to the array of track parameters
365   if (!fTrackParamAtCluster) {
366     fTrackParamAtCluster = new TObjArray(20);
367     fTrackParamAtCluster->SetOwner(kTRUE);
368   }
369   AliMUONTrackParam* trackParamAtCluster = new AliMUONTrackParam(trackParam);
370   fTrackParamAtCluster->AddLast(trackParamAtCluster);
371   
372   // link parameters with the associated cluster or its copy
373   if (copy) {
374     AliMUONVCluster *clusterCopy = static_cast<AliMUONVCluster*>(cluster.Clone());
375     trackParamAtCluster->SetClusterPtr(clusterCopy, kTRUE);
376   } else trackParamAtCluster->SetClusterPtr(&cluster);
377   
378   // sort the array of track parameters
379   fTrackParamAtCluster->Sort();
380 }
381
382   //__________________________________________________________________________
383 void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam)
384 {
385   /// Remove trackParam from the array of TrackParamAtCluster and delete it since the array is owner
386   
387   if (fTrackParamAtCluster) {
388     
389     AliMUONTrackParam* trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->Remove(trackParam));
390     
391     if (trackParamAtCluster) {
392       
393       // clean memory
394       delete trackParamAtCluster;
395       
396       // remove hole
397       fTrackParamAtCluster->Compress();
398       
399     } else AliWarning("object to remove does not exist in array fTrackParamAtCluster");
400     
401   } else AliWarning("array fTrackParamAtCluster does not exist");
402   
403 }
404
405   //__________________________________________________________________________
406 Bool_t AliMUONTrack::UpdateTrackParamAtCluster()
407 {
408   /// Update track parameters at each attached cluster
409   /// Return kFALSE in case of failure (i.e. extrapolation problem)
410   
411   Int_t nClusters = GetNClusters();
412   if (nClusters == 0) {
413     AliWarning("no cluster attached to the track");
414     return kFALSE;
415   }
416   
417   Bool_t extrapStatus = kTRUE;
418   AliMUONTrackParam* startingTrackParam = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(0));
419   
420   for (Int_t i = 1; i < nClusters; i++) {
421     AliMUONTrackParam* trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(i));
422     
423     // reset track parameters and their covariances
424     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
425     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
426     
427     // extrapolation to the given z
428     if (!AliMUONTrackExtrap::ExtrapToZ(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ())) extrapStatus = kFALSE;
429     
430     // prepare next step
431     startingTrackParam = trackParamAtCluster;
432   }
433
434   // set global chi2 to max value in case of problem during track extrapolation
435   if (!extrapStatus) SetGlobalChi2(2.*MaxChi2());
436   return extrapStatus;
437   
438 }
439
440   //__________________________________________________________________________
441 Bool_t AliMUONTrack::UpdateCovTrackParamAtCluster()
442 {
443   /// Update track parameters and their covariances at each attached cluster
444   /// Include effects of multiple scattering in chambers
445   /// Return kFALSE in case of failure (i.e. extrapolation problem)
446   
447   Int_t nClusters = GetNClusters();
448   if (nClusters == 0) {
449     AliWarning("no cluster attached to the track");
450     return kFALSE;
451   }
452   
453   Bool_t extrapStatus = kTRUE;
454   AliMUONTrackParam* startingTrackParam = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(0));
455   Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() + 1;
456   Int_t currentChamber;
457   
458   for (Int_t i = 1; i < nClusters; i++) {
459     AliMUONTrackParam* trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(i));
460     
461     // reset track parameters and their covariances
462     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
463     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
464     trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
465     
466     // add MCS effect
467     AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber-1),-1.);
468     
469     // add MCS in missing chambers if any
470     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
471     while (currentChamber > expectedChamber) {
472       // extrapolation to the missing chamber
473       if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber))) extrapStatus = kFALSE;
474       // add MCS effect
475       AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber),-1.);
476       expectedChamber++;
477     }
478     
479     // extrapolation to the z of the current cluster
480     if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ())) extrapStatus = kFALSE;
481     
482     // prepare next step
483     expectedChamber = currentChamber + 1;
484     startingTrackParam = trackParamAtCluster;
485   }
486   
487   // set global chi2 to max value in case of problem during track extrapolation
488   if (!extrapStatus) SetGlobalChi2(2.*MaxChi2());
489   return extrapStatus;
490   
491 }
492
493   //__________________________________________________________________________
494 Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
495 {
496   /// check the validity of the current track:
497   /// at least one cluster per requested station
498   /// and at least 2 chambers in stations 4 & 5 that contain cluster(s)
499   /// + if request2ChInSameSt45 = kTRUE: 2 chambers hit in the same station (4 or 5)
500   
501   Int_t nClusters = GetNClusters();
502   AliMUONTrackParam *trackParam;
503   Int_t currentCh, currentSt, previousCh = -1, nChHitInSt4 = 0, nChHitInSt5 = 0;
504   UInt_t presentStationMask(0);
505   
506   // first loop over clusters
507   for (Int_t i = 0; i < nClusters; i++) {
508     trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
509     
510     currentCh = trackParam->GetClusterPtr()->GetChamberId();
511     currentSt = currentCh/2;
512     
513     // build present station mask
514     presentStationMask |= ( 1 << currentSt );
515     
516     // count the number of chambers hit in station 4 that contain cluster(s)
517     if (currentSt == 3 && currentCh != previousCh) {
518       nChHitInSt4++;
519       previousCh = currentCh;
520     }
521     
522     // count the number of chambers hit in station 5 that contain cluster(s)
523     if (currentSt == 4 && currentCh != previousCh) {
524       nChHitInSt5++;
525       previousCh = currentCh;
526     }
527     
528   }
529   
530   // at least one cluster per requested station
531   if ((requestedStationMask & presentStationMask) != requestedStationMask) return kFALSE;
532   
533   // 2 chambers hit in the same station (4 or 5)
534   if (request2ChInSameSt45) return (nChHitInSt4 == 2 || nChHitInSt5 == 2);
535   // or 2 chambers hit in station 4 & 5 together
536   else return (nChHitInSt4+nChHitInSt5 >= 2);
537   
538 }
539
540   //__________________________________________________________________________
541 void AliMUONTrack::TagRemovableClusters(UInt_t requestedStationMask) {
542   /// Identify clusters that can be removed from the track,
543   /// with the only requirements to have at least 1 cluster per requested station
544   /// and at least 2 chambers over 4 in stations 4 & 5 that contain cluster(s)
545   
546   Int_t nClusters = GetNClusters();
547   AliMUONTrackParam *trackParam, *nextTrackParam;
548   Int_t currentCh, nextCh, currentSt, nextSt, previousCh = -1, nChHitInSt45 = 0;
549   
550   // first loop over clusters
551   for (Int_t i = 0; i < nClusters; i++) {
552     trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
553     
554     currentCh = trackParam->GetClusterPtr()->GetChamberId();
555     currentSt = currentCh/2;
556     
557     // reset flags to kFALSE for all clusters in required station
558     if ((1 << currentSt) & requestedStationMask) trackParam->SetRemovable(kFALSE);
559     else trackParam->SetRemovable(kTRUE);
560     
561     // count the number of chambers in station 4 & 5 that contain cluster(s)
562     if (currentCh > 5 && currentCh != previousCh) {
563       nChHitInSt45++;
564       previousCh = currentCh;
565     }
566     
567   }
568   
569   // second loop over clusters
570   for (Int_t i = 0; i < nClusters; i++) {
571     trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
572     
573     currentCh = trackParam->GetClusterPtr()->GetChamberId();
574     currentSt = currentCh/2;
575     
576     // make sure they are more than 2 clusters in 2 different chambers of stations 4 & 5
577     // but 2 clusters in he same chamber will still be flagged as removable
578     if (nChHitInSt45 < 3 && currentSt > 2) {
579       
580       if (i == nClusters-1) {
581         
582         trackParam->SetRemovable(kFALSE);
583       
584       } else {
585         
586         nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i+1);
587         nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
588         
589         // set clusters in the same chamber as being removable
590         if (nextCh == currentCh) {
591           trackParam->SetRemovable(kTRUE);
592           nextTrackParam->SetRemovable(kTRUE);
593           i++; // skip cluster already checked
594         } else {
595           trackParam->SetRemovable(kFALSE);
596         }
597         
598       }
599       
600     } else {
601       
602       // skip clusters already flag as removable
603       if (trackParam->IsRemovable()) continue;
604       
605       // loop over next track parameters
606       for (Int_t j = i+1; j < nClusters; j++) {
607         nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(j);
608         
609         nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
610         nextSt = nextCh/2;
611         
612         // set clusters in the same station as being removable
613         if (nextSt == currentSt) {
614           trackParam->SetRemovable(kTRUE);
615           nextTrackParam->SetRemovable(kTRUE);
616           i++; // skip cluster already checked
617         }
618         
619       }
620       
621     }
622       
623   }
624     
625 }
626
627   //__________________________________________________________________________
628 Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
629 {
630   /// Compute each cluster contribution to the chi2 of the track
631   /// accounting for multiple scattering or not according to the flag
632   /// - Also recompute the weight matrices of the attached clusters if accountForMCS=kTRUE
633   /// - Assume that track parameters at each cluster are corrects
634   /// - Return kFALSE if computation failed
635   AliDebug(1,"Enter ComputeLocalChi2");
636   
637   if (!fTrackParamAtCluster) {
638     AliWarning("no cluster attached to this track");
639     return kFALSE;
640   }
641   
642   if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
643       
644     // Compute MCS covariance matrix only once
645     Int_t nClusters = GetNClusters();
646     TMatrixD mcsCovariances(nClusters,nClusters);
647     ComputeMCSCovariances(mcsCovariances);
648     
649     // Make sure cluster weights are consistent with following calculations
650     if (!ComputeClusterWeights(&mcsCovariances)) {
651       AliWarning("cannot take into account the multiple scattering effects");
652       return ComputeLocalChi2(kFALSE);
653     }
654     
655     // Compute chi2 of the track
656     Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
657     if (globalChi2 < 0.) return kFALSE;
658     
659     // Loop over removable clusters and compute their local chi2
660     AliMUONTrackParam* trackParamAtCluster;
661     AliMUONTrackParam* trackParamAtCluster1;
662     AliMUONVCluster *cluster, *discardedCluster;
663     Int_t iCluster1, iCluster2, iCurrentCluster1, iCurrentCluster2;
664     TMatrixD clusterWeightsNB(nClusters-1,nClusters-1);
665     TMatrixD clusterWeightsB(nClusters-1,nClusters-1);
666     Double_t *dX = new Double_t[nClusters-1];
667     Double_t *dY = new Double_t[nClusters-1];
668     Double_t globalChi2b;
669     for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) { 
670       trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(iCluster));
671       
672       discardedCluster = trackParamAtCluster->GetClusterPtr();
673       
674       // Recompute cluster weights without the current cluster
675       if (!ComputeClusterWeights(clusterWeightsNB, clusterWeightsB, &mcsCovariances, discardedCluster)) {
676         AliWarning("cannot take into account the multiple scattering effects");
677         delete [] dX;
678         delete [] dY;
679         return ComputeLocalChi2(kFALSE);
680       }
681       
682       // Compute track chi2 without the current cluster
683       globalChi2b = 0.;
684       iCurrentCluster1 = 0;
685       for (iCluster1 = 0; iCluster1 < nClusters ; iCluster1++) { 
686         trackParamAtCluster1 = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
687         cluster = trackParamAtCluster1->GetClusterPtr();
688         
689         if (cluster == discardedCluster) continue;
690         
691         // Compute and save residuals
692         dX[iCurrentCluster1] = cluster->GetX() - trackParamAtCluster1->GetNonBendingCoor();
693         dY[iCurrentCluster1] = cluster->GetY() - trackParamAtCluster1->GetBendingCoor();
694         
695         iCurrentCluster2 = 0;
696         for (iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
697           cluster = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
698           
699           if (cluster == discardedCluster) continue;
700           
701           // Add contribution from covariances
702           globalChi2b += (clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) +
703                           clusterWeightsNB(iCurrentCluster2, iCurrentCluster1)) * dX[iCurrentCluster1] * dX[iCurrentCluster2] +
704                          (clusterWeightsB(iCurrentCluster1, iCurrentCluster2) +
705                           clusterWeightsB(iCurrentCluster2, iCurrentCluster1)) * dY[iCurrentCluster1] * dY[iCurrentCluster2];
706           
707           iCurrentCluster2++;
708         }
709         
710         // Add contribution from variances
711         globalChi2b += clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) * dX[iCurrentCluster1] * dX[iCurrentCluster1] +
712                        clusterWeightsB(iCurrentCluster1, iCurrentCluster1) * dY[iCurrentCluster1] * dY[iCurrentCluster1];
713         
714         iCurrentCluster1++;
715       }
716
717       // Set local chi2
718       trackParamAtCluster->SetLocalChi2(globalChi2 - globalChi2b);
719     }
720     
721     delete [] dX;
722     delete [] dY;
723     
724   } else { // without multiple scattering effects
725     
726     Int_t nClusters = GetNClusters();
727     AliMUONTrackParam* trackParamAtCluster;
728     AliMUONVCluster *discardedCluster;
729     Double_t dX, dY;
730     for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) { 
731       trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(iCluster));
732       
733       discardedCluster = trackParamAtCluster->GetClusterPtr();
734       
735       // Compute residuals
736       dX = discardedCluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
737       dY = discardedCluster->GetY() - trackParamAtCluster->GetBendingCoor();
738       
739       // Set local chi2
740       trackParamAtCluster->SetLocalChi2(dX * dX / discardedCluster->GetErrX2() + dY * dY / discardedCluster->GetErrY2());
741     }
742   
743   }
744   
745   return kTRUE;
746   
747 }
748
749   //__________________________________________________________________________
750 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
751 {
752   /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
753   /// - Assume that track parameters at each cluster are corrects
754   /// - Assume the cluster weights matrices are corrects
755   /// - Return a value of chi2 higher than the maximum allowed if computation failed
756   AliDebug(1,"Enter ComputeGlobalChi2");
757   
758   if (!fTrackParamAtCluster) {
759     AliWarning("no cluster attached to this track");
760     return 2.*MaxChi2();
761   }
762   
763   Double_t chi2 = 0.;
764   
765   if (accountForMCS) {
766     
767     // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
768     if (!fClusterWeightsNonBending || !fClusterWeightsBending) {
769       AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
770       return ComputeGlobalChi2(kFALSE);
771     }
772     Int_t nClusters = GetNClusters();
773     if (fClusterWeightsNonBending->GetNrows() != nClusters || fClusterWeightsBending->GetNcols() != nClusters) {
774       AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
775       return ComputeGlobalChi2(kFALSE);
776     }
777     
778     // Compute chi2
779     AliMUONVCluster *cluster;
780     Double_t *dX = new Double_t[nClusters];
781     Double_t *dY = new Double_t[nClusters];
782     AliMUONTrackParam* trackParamAtCluster;
783     for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
784       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
785       cluster = trackParamAtCluster->GetClusterPtr();
786       dX[iCluster1] = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
787       dY[iCluster1] = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
788       for (Int_t iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
789         chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster2) + (*fClusterWeightsNonBending)(iCluster2, iCluster1)) * dX[iCluster1] * dX[iCluster2] +
790                 ((*fClusterWeightsBending)(iCluster1, iCluster2) + (*fClusterWeightsBending)(iCluster2, iCluster1)) * dY[iCluster1] * dY[iCluster2];
791       }
792       chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster1) * dX[iCluster1] * dX[iCluster1]) +
793               ((*fClusterWeightsBending)(iCluster1, iCluster1) * dY[iCluster1] * dY[iCluster1]);
794     }
795     delete [] dX;
796     delete [] dY;
797     
798   } else {
799     
800     AliMUONVCluster *cluster;
801     Double_t dX, dY;
802     AliMUONTrackParam* trackParamAtCluster;
803     Int_t nClusters = GetNClusters();
804     for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) { 
805       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
806       cluster = trackParamAtCluster->GetClusterPtr();
807       dX = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
808       dY = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
809       chi2 += dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
810     }
811     
812   }
813   
814   return chi2;
815   
816 }
817
818   //__________________________________________________________________________
819 Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances)
820 {
821   /// Compute the weight matrices of the attached clusters, in non bending and bending direction,
822   /// accounting for multiple scattering correlations and cluster resolution
823   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
824   /// - Assume that track parameters at each cluster are corrects
825   /// - Return kFALSE if computation failed
826   AliDebug(1,"Enter ComputeClusterWeights1");
827   
828   if (!fTrackParamAtCluster) {
829     AliWarning("no cluster attached to this track");
830     return kFALSE;
831   }
832   
833   // Alocate memory
834   Int_t nClusters = GetNClusters();
835   if (!fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(nClusters,nClusters);
836   if (!fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(nClusters,nClusters);
837   
838   // Compute weights matrices
839   if (!ComputeClusterWeights(*fClusterWeightsNonBending, *fClusterWeightsBending, mcsCovariances)) return kFALSE;
840   
841   return kTRUE;
842   
843 }
844
845   //__________________________________________________________________________
846 Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& clusterWeightsNB, TMatrixD& clusterWeightsB,
847                                            TMatrixD* mcsCovariances, const AliMUONVCluster* discardedCluster) const
848 {
849   /// Compute the weight matrices, in non bending and bending direction,
850   /// of the other attached clusters assuming the discarded one does not exist
851   /// accounting for multiple scattering correlations and cluster resolution
852   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
853   /// - Return kFALSE if computation failed
854   AliDebug(1,"Enter ComputeClusterWeights2");
855   
856   // Check MCS covariance matrix and recompute it if need
857   Int_t nClusters = GetNClusters();
858   Bool_t deleteMCSCov = kFALSE;
859   if (!mcsCovariances) {
860     mcsCovariances = new TMatrixD(nClusters,nClusters);
861     deleteMCSCov = kTRUE;
862     ComputeMCSCovariances(*mcsCovariances);
863   }
864   
865   // Resize the weights matrices; alocate memory
866   if (discardedCluster) {
867     clusterWeightsNB.ResizeTo(nClusters-1,nClusters-1);
868     clusterWeightsB.ResizeTo(nClusters-1,nClusters-1);
869   } else {
870     clusterWeightsNB.ResizeTo(nClusters,nClusters);
871     clusterWeightsB.ResizeTo(nClusters,nClusters);
872   }
873   
874   // Define variables
875   AliMUONVCluster *cluster1, *cluster2;
876   Int_t iCurrentCluster1, iCurrentCluster2;
877   
878   // Compute the covariance matrices
879   iCurrentCluster1 = 0;
880   for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
881     cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
882     
883     if (cluster1 == discardedCluster) continue;
884     
885     // Loop over next clusters
886     iCurrentCluster2 = iCurrentCluster1;
887     for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
888       cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
889       
890       if (cluster2 == discardedCluster) continue;
891       
892       // Fill with MCS covariances
893       clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) = (*mcsCovariances)(iCluster1,iCluster2);
894       
895       // Equal contribution from multiple scattering in non bending and bending directions
896       clusterWeightsB(iCurrentCluster1, iCurrentCluster2) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
897       
898       // Add contribution from cluster resolution to diagonal element and symmetrize the matrix
899       if (iCurrentCluster1 == iCurrentCluster2) {
900         
901         // In non bending plane
902         clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrX2();
903         // In bending plane
904         clusterWeightsB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrY2();
905         
906       } else {
907         
908         // In non bending plane
909         clusterWeightsNB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
910         // In bending plane
911         clusterWeightsB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsB(iCurrentCluster1, iCurrentCluster2);
912         
913       }
914       
915       iCurrentCluster2++;
916     }
917     
918     iCurrentCluster1++;
919   }
920     
921   // Inversion of covariance matrices to get the weights
922   if (clusterWeightsNB.Determinant() != 0 && clusterWeightsB.Determinant() != 0) {
923     clusterWeightsNB.Invert();
924     clusterWeightsB.Invert();
925   } else {
926     AliWarning(" Determinant = 0");
927     clusterWeightsNB.ResizeTo(0,0);
928     clusterWeightsB.ResizeTo(0,0);
929     if(deleteMCSCov) delete mcsCovariances;
930     return kFALSE;
931   }
932   
933   if(deleteMCSCov) delete mcsCovariances;
934   
935   return kTRUE;
936   
937 }
938
939   //__________________________________________________________________________
940 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
941 {
942   /// Compute the multiple scattering covariance matrix
943   /// (assume that track parameters at each cluster are corrects)
944   AliDebug(1,"Enter ComputeMCSCovariances");
945   
946   // Reset the size of the covariance matrix if needed
947   Int_t nClusters = GetNClusters();
948   if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters);
949   
950   // Define variables
951   Int_t nChambers = AliMUONConstants::NTrackingCh();
952   AliMUONTrackParam* trackParamAtCluster;
953   AliMUONTrackParam extrapTrackParam;
954   Int_t currentChamber = 0, expectedChamber = 0, size = 0;
955   Double_t *mcsAngle2 = new Double_t[2*nChambers];
956   Double_t *zMCS = new Double_t[2*nChambers];
957   Int_t *indices = new Int_t[2*nClusters];
958   
959   // Compute multiple scattering dispersion angle at each chamber
960   // and save the z position where it is calculated
961   for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
962     trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
963     
964     // look for missing chambers if any
965     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
966     while (currentChamber > expectedChamber) {
967       
968       // Save the z position where MCS dispersion is calculated
969       zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
970       
971       // Do not take into account MCS in chambers prior the first cluster
972       if (iCluster > 0) {
973         
974         // Get track parameters at missing chamber z
975         extrapTrackParam = *trackParamAtCluster;
976         AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]);
977         
978         // Save multiple scattering dispersion angle in missing chamber
979         mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(expectedChamber),1.);
980         
981       } else mcsAngle2[size] = 0.;
982       
983       expectedChamber++;
984       size++;
985     }
986     
987     // Save z position where MCS dispersion is calculated
988     zMCS[size] = trackParamAtCluster->GetZ();
989     
990     // Save multiple scattering dispersion angle in current chamber
991     mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(currentChamber),1.);
992     
993     // Save indice in zMCS array corresponding to the current cluster
994     indices[iCluster] = size;
995     
996     expectedChamber = currentChamber + 1;
997     size++;
998   }
999   
1000   // complete array of z if last cluster is on the last but one chamber
1001   if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
1002   
1003   // Compute the covariance matrix
1004   for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
1005     
1006     for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
1007       
1008       // Initialization to 0 (diagonal plus upper triangular part)
1009       mcsCovariances(iCluster1,iCluster2) = 0.;
1010       
1011       // Compute contribution from multiple scattering in upstream chambers
1012       for (Int_t k = 0; k < indices[iCluster1]; k++) {  
1013         mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k];
1014       }
1015       
1016       // Symetrize the matrix
1017       mcsCovariances(iCluster2,iCluster1) = mcsCovariances(iCluster1,iCluster2);
1018     }
1019     
1020   }
1021     
1022   delete [] mcsAngle2;
1023   delete [] zMCS;
1024   delete [] indices;
1025   
1026 }
1027
1028   //__________________________________________________________________________
1029 Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track, Int_t stMin, Int_t stMax) const
1030 {
1031   /// Returns the number of clusters in common in stations [stMin, stMax]
1032   /// between the current track ("this") and the track pointed to by "track".
1033   
1034   if (!track || !track->fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0;
1035   
1036   Int_t chMin = 2 * stMin;
1037   Int_t chMax = 2 * stMax + 1;
1038   Int_t clustersInCommon = 0;
1039   
1040   // Loop over clusters of first track
1041   Int_t nCl1 = this->GetNClusters();
1042   for(Int_t iCl1 = 0; iCl1 < nCl1; iCl1++) {
1043     AliMUONVCluster* cl1 = ((AliMUONTrackParam*) this->fTrackParamAtCluster->UncheckedAt(iCl1))->GetClusterPtr();
1044     
1045     Int_t chCl1 = cl1->GetChamberId();
1046     if (chCl1 < chMin || chCl1 > chMax) continue;
1047     
1048     // Loop over clusters of second track
1049     Int_t nCl2 = track->GetNClusters();
1050     for(Int_t iCl2 = 0; iCl2 < nCl2; iCl2++) {
1051       AliMUONVCluster* cl2 = ((AliMUONTrackParam*) track->fTrackParamAtCluster->UncheckedAt(iCl2))->GetClusterPtr();
1052       
1053       Int_t chCl2 = cl2->GetChamberId();
1054       if (chCl2 < chMin || chCl2 > chMax) continue;
1055       
1056       // Increment "clustersInCommon" if both clusters have the same ID
1057       if (cl1->GetUniqueID() == cl2->GetUniqueID()) {
1058         clustersInCommon++;
1059         break;
1060       }
1061       
1062     }
1063     
1064   }
1065   
1066   return clustersInCommon;
1067 }
1068
1069   //__________________________________________________________________________
1070 Int_t AliMUONTrack::GetNDF() const
1071 {
1072   /// return the number of degrees of freedom
1073   
1074   Int_t ndf = 2 * GetNClusters() - 5;
1075   return (ndf > 0) ? ndf : 0;
1076 }
1077
1078   //__________________________________________________________________________
1079 Double_t AliMUONTrack::GetNormalizedChi2() const
1080 {
1081   /// return the chi2 value divided by the number of degrees of freedom (or FLT_MAX if ndf <= 0)
1082   
1083   Double_t ndf = (Double_t) GetNDF();
1084   return (ndf > 0.) ? fGlobalChi2 / ndf : 2.*MaxChi2();
1085 }
1086
1087   //__________________________________________________________________________
1088 Int_t AliMUONTrack::FindCompatibleClusters(const AliMUONTrack &track, Double_t sigmaCut, Bool_t compatibleCluster[10]) const
1089 {
1090   /// Try to match clusters from this track with clusters from the given track within the provided sigma cut:
1091   /// - Fill the array compatibleCluster[iCh] with kTRUE if a compatible cluster has been found in chamber iCh.
1092   /// - Return the number of clusters of "this" track matched with one cluster of the given track.
1093   AliMUONVCluster *cluster1, *cluster2;
1094   Double_t chi2, dX, dY;
1095   Double_t chi2Max = sigmaCut * sigmaCut;
1096   
1097   // initialization
1098   Int_t nMatchClusters = 0;
1099   for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE;
1100
1101   if (!track.fTrackParamAtCluster || !this->fTrackParamAtCluster) return nMatchClusters;
1102   
1103   // Loop over clusters of first track
1104   Int_t nCl1 = this->GetNClusters();
1105   for(Int_t iCl1 = 0; iCl1 < nCl1; iCl1++) {
1106     cluster1 = static_cast<AliMUONTrackParam*>(this->fTrackParamAtCluster->UncheckedAt(iCl1))->GetClusterPtr();
1107     
1108     // Loop over clusters of second track
1109     Int_t nCl2 = track.GetNClusters();
1110     for(Int_t iCl2 = 0; iCl2 < nCl2; iCl2++) {
1111       cluster2 = static_cast<AliMUONTrackParam*>(track.fTrackParamAtCluster->UncheckedAt(iCl2))->GetClusterPtr();
1112       
1113       // check DE Id
1114       if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue;
1115       
1116       // check local chi2
1117       dX = cluster1->GetX() - cluster2->GetX();
1118       dY = cluster1->GetY() - cluster2->GetY();
1119       chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()) + dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2());
1120       if (chi2 > 2. * chi2Max) continue; // 2 because 2 quantities in chi2
1121       
1122       compatibleCluster[cluster1->GetChamberId()] = kTRUE;
1123       nMatchClusters++;
1124       break;
1125     }
1126     
1127   }
1128   
1129   return nMatchClusters;
1130 }
1131
1132 //__________________________________________________________________________
1133 Bool_t AliMUONTrack::Match(AliMUONTrack &track, Double_t sigmaCut, Int_t &nMatchClusters) const
1134 {
1135   /// Try to match this track with the given track. Matching conditions:
1136   /// - more than 50% of clusters from this track matched with clusters from the given track
1137   /// - at least 1 cluster matched before and 1 cluster matched after the dipole
1138   
1139   Bool_t compTrack[10];
1140   nMatchClusters = FindCompatibleClusters(track, sigmaCut, compTrack);
1141   
1142   if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
1143       (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
1144       2 * nMatchClusters > GetNClusters()) return kTRUE;                // more than 50% of clusters matched
1145   else return kFALSE;
1146   
1147 }
1148
1149 //__________________________________________________________________________
1150 void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam)
1151 {
1152   /// set track parameters at vertex
1153   if (trackParam == 0x0) return;
1154   if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam;
1155   else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam);
1156 }
1157
1158 //__________________________________________________________________________
1159 void AliMUONTrack::RecursiveDump() const
1160 {
1161   /// Recursive dump of AliMUONTrack, i.e. with dump of trackParamAtCluster and attached clusters
1162   AliMUONTrackParam *trackParamAtCluster;
1163   AliMUONVCluster *cluster;
1164   cout << "Recursive dump of Track: " << this << endl;
1165   // Track
1166   this->Dump();
1167   for (Int_t iCluster = 0; iCluster < GetNClusters(); iCluster++) {
1168     trackParamAtCluster = (AliMUONTrackParam*) ((*fTrackParamAtCluster)[iCluster]);
1169     // trackParamAtCluster
1170     cout << "trackParamAtCluster: " << trackParamAtCluster << " (index: " << iCluster << ")" << endl;
1171     trackParamAtCluster->Dump();
1172     cluster = trackParamAtCluster->GetClusterPtr();
1173     // cluster
1174     cout << "cluster: " << cluster << endl;
1175     cluster->Print();
1176   }
1177   return;
1178 }
1179   
1180 //_____________________________________________-
1181 void AliMUONTrack::Print(Option_t*) const
1182 {
1183   /// Printing Track information 
1184
1185   streamsize curW = cout.width();
1186   streamsize curPrecision = cout.precision();
1187   cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNClusters() << 
1188       ", Match2Trig=" << setw(1) << GetMatchTrigger()  << 
1189       ", LoTrgNum=" << setw(3) << LoCircuit()  << 
1190     ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) <<  GetChi2MatchTrigger();
1191   cout << Form(" HitTriggerPattern trig %x  track %x",fHitsPatternInTrigCh, fHitsPatternInTrigChTrk);
1192   cout << Form(" MClabel=%d",fTrackID) << endl;
1193   if (fTrackParamAtCluster) fTrackParamAtCluster->First()->Print("FULL");
1194   cout.width(curW);
1195   cout.precision(curPrecision);
1196 }
1197
1198 //__________________________________________________________________________
1199 void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt, UChar_t respWithoutChamber)
1200 {
1201   /// pack the local trigger information and store
1202
1203   if (loCirc < 0) return;
1204
1205   fLocalTrigger = 0;
1206   fLocalTrigger += loCirc;
1207   fLocalTrigger += loStripX << 8;
1208   fLocalTrigger += loStripY << 13;
1209   fLocalTrigger += loDev    << 17;
1210   fLocalTrigger += loLpt    << 22;
1211   fLocalTrigger += loHpt    << 24;
1212   fLocalTrigger += respWithoutChamber << 26;
1213
1214 }
1215
1216 //__________________________________________________________________________
1217 void AliMUONTrack::FindMCLabel()
1218 {
1219   /// Determine the MC label from the label of the attached clusters and fill fMCLabel data member:
1220   /// More than 50% of clusters, including 1 before and 1 after the dipole, must share the same label
1221   
1222   Int_t nClusters = GetNClusters();
1223   Int_t halfCluster = nClusters/2;
1224   
1225   // reset MC label
1226   fTrackID = -1;
1227   
1228   // loop over first clusters (if nClusters left < (nClusters-halfCluster) the conditions cannot be fulfilled)
1229   for (Int_t iCluster1 = 0; iCluster1 < nClusters-halfCluster; iCluster1++) {
1230     AliMUONVCluster* cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
1231     
1232     // if the first cluster is not on station 1 or 2 the conditions cannot be fulfilled
1233     if (cluster1->GetChamberId() > 3) return;
1234     
1235     Int_t label1 = cluster1->GetMCLabel();
1236     if (label1 < 0) continue;
1237     
1238     Int_t nIdenticalLabel = 1;
1239     
1240     // Loop over next clusters
1241     for (Int_t iCluster2 = iCluster1+1; iCluster2 < nClusters; iCluster2++) {
1242       AliMUONVCluster* cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
1243       
1244       if (cluster2->GetMCLabel() != label1) continue;
1245       
1246       nIdenticalLabel++;
1247       
1248       // stop as soon as conditions are fulfilled
1249       if (nIdenticalLabel > halfCluster && cluster2->GetChamberId() > 5) {
1250         fTrackID = label1;
1251         return;
1252       }
1253       
1254     }
1255     
1256   }
1257   
1258 }
1259