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