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