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