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