76eb487fd22128c8b89a44a79a5414f8aa3be145
[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 "AliMUONTrackParam.h" 
27 #include "AliMUONVCluster.h" 
28 #include "AliMUONObjectPair.h" 
29 #include "AliMUONConstants.h"
30 #include "AliMUONTrackExtrap.h" 
31
32 #include "AliLog.h"
33
34 #include <TMath.h>
35 #include <TMatrixD.h>
36
37 #include <Riostream.h>
38
39 /// \cond CLASSIMP
40 ClassImp(AliMUONTrack) // Class implementation in ROOT context
41 /// \endcond
42
43 //__________________________________________________________________________
44 AliMUONTrack::AliMUONTrack()
45   : TObject(),
46     fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
47     fFitWithVertex(kFALSE),
48     fVertexErrXY2(),
49     fFitWithMCS(kFALSE),
50     fClusterWeightsNonBending(0x0),
51     fClusterWeightsBending(0x0),
52     fGlobalChi2(-1.),
53     fImproved(kFALSE),
54     fMatchTrigger(-1),
55     floTrgNum(-1),
56     fChi2MatchTrigger(0.),
57     fTrackID(0),
58     fTrackParamAtVertex(0x0),
59     fHitsPatternInTrigCh(0),
60     fLocalTrigger(0)
61 {
62   /// Default constructor
63   fTrackParamAtCluster->SetOwner(kTRUE);
64   fVertexErrXY2[0] = 0.;
65   fVertexErrXY2[1] = 0.;
66 }
67
68   //__________________________________________________________________________
69 AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment)
70   : TObject(),
71     fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
72     fFitWithVertex(kFALSE),
73     fVertexErrXY2(),
74     fFitWithMCS(kFALSE),
75     fClusterWeightsNonBending(0x0),
76     fClusterWeightsBending(0x0),
77     fGlobalChi2(0.),
78     fImproved(kFALSE),
79     fMatchTrigger(-1),
80     floTrgNum(-1),    
81     fChi2MatchTrigger(0.),
82     fTrackID(0),
83     fTrackParamAtVertex(0x0),
84     fHitsPatternInTrigCh(0),
85     fLocalTrigger(0)
86 {
87   /// Constructor from two clusters
88   fTrackParamAtCluster->SetOwner(kTRUE);
89   
90   fVertexErrXY2[0] = 0.;
91   fVertexErrXY2[1] = 0.;
92   
93   // Pointers to clusters from the segment
94   AliMUONVCluster* cluster1 = (AliMUONVCluster*) segment->First();
95   AliMUONVCluster* cluster2 = (AliMUONVCluster*) segment->Second();
96   
97   // check sorting in -Z (spectro z<0)
98   if (cluster1->GetZ() < cluster2->GetZ()) {
99     cluster1 = cluster2;
100     cluster2 = (AliMUONVCluster*) segment->First();
101   }
102   
103   // order the clusters into the track according to the station the segment belong to
104   //(the cluster first attached is the one from which we will start the tracking procedure)
105   AliMUONVCluster *firstCluster, *lastCluster;
106   if (cluster1->GetChamberId() == 8) {
107     firstCluster = cluster1;
108     lastCluster = cluster2;
109   } else {
110     firstCluster = cluster2;
111     lastCluster = cluster1;
112   }
113   
114   // Compute track parameters
115   Double_t z1 = firstCluster->GetZ();
116   Double_t z2 = lastCluster->GetZ();
117   Double_t dZ = z1 - z2;
118   // Non bending plane
119   Double_t nonBendingCoor1 = firstCluster->GetX();
120   Double_t nonBendingCoor2 = lastCluster->GetX();
121   Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ;
122   // Bending plane
123   Double_t bendingCoor1 = firstCluster->GetY();
124   Double_t bendingCoor2 = lastCluster->GetY();
125   Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ;
126   // Inverse bending momentum
127   Double_t bendingImpact = bendingCoor1 - z1 * bendingSlope;
128   Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
129   
130   
131   // Set track parameters at first cluster
132   AliMUONTrackParam trackParamAtFirstCluster;
133   trackParamAtFirstCluster.SetZ(z1);
134   trackParamAtFirstCluster.SetNonBendingCoor(nonBendingCoor1);
135   trackParamAtFirstCluster.SetNonBendingSlope(nonBendingSlope);
136   trackParamAtFirstCluster.SetBendingCoor(bendingCoor1);
137   trackParamAtFirstCluster.SetBendingSlope(bendingSlope);
138   trackParamAtFirstCluster.SetInverseBendingMomentum(inverseBendingMomentum);
139   
140   
141   // Set track parameters at last cluster
142   AliMUONTrackParam trackParamAtLastCluster;
143   trackParamAtLastCluster.SetZ(z2);
144   trackParamAtLastCluster.SetNonBendingCoor(nonBendingCoor2);
145   trackParamAtLastCluster.SetNonBendingSlope(nonBendingSlope);
146   trackParamAtLastCluster.SetBendingCoor(bendingCoor2);
147   trackParamAtLastCluster.SetBendingSlope(bendingSlope);
148   trackParamAtLastCluster.SetInverseBendingMomentum(inverseBendingMomentum);
149   
150   
151   // Compute and set track parameters covariances at first cluster
152   TMatrixD paramCov1(5,5);
153   paramCov1.Zero();
154   // Non bending plane
155   paramCov1(0,0) = firstCluster->GetErrX2();
156   paramCov1(0,1) = firstCluster->GetErrX2() / dZ;
157   paramCov1(1,0) = paramCov1(0,1);
158   paramCov1(1,1) = ( firstCluster->GetErrX2() + lastCluster->GetErrX2() ) / dZ / dZ;
159   // Bending plane
160   paramCov1(2,2) = firstCluster->GetErrY2();
161   paramCov1(2,3) = firstCluster->GetErrY2() / dZ;
162   paramCov1(3,2) = paramCov1(2,3);
163   paramCov1(3,3) = ( firstCluster->GetErrY2() + lastCluster->GetErrY2() ) / dZ / dZ;
164   // Inverse bending momentum (50% error)
165   paramCov1(4,4) = 0.5*inverseBendingMomentum * 0.5*inverseBendingMomentum;
166   // Set covariances
167   trackParamAtFirstCluster.SetCovariances(paramCov1);
168   
169   
170   // Compute and set track parameters covariances at last cluster (as if the first cluster did not exist)
171   TMatrixD paramCov2(5,5);
172   paramCov2.Zero();
173   // Non bending plane
174   paramCov2(0,0) = paramCov1(0,0);
175   paramCov2(1,1) = 100.*paramCov1(1,1);
176   // Bending plane
177   paramCov2(2,2) = paramCov1(2,2);
178   paramCov2(3,3) = 100.*paramCov1(3,3);
179   // Inverse bending momentum
180   paramCov2(4,4) = paramCov1(4,4);
181   // Set covariances
182   trackParamAtLastCluster.SetCovariances(paramCov2);
183   
184   // Flag clusters as being removable
185   trackParamAtFirstCluster.SetRemovable(kTRUE);
186   trackParamAtLastCluster.SetRemovable(kTRUE);
187   
188   // Add track parameters at clusters
189   AddTrackParamAtCluster(trackParamAtFirstCluster,*firstCluster);
190   AddTrackParamAtCluster(trackParamAtLastCluster,*lastCluster);
191   
192 }
193
194   //__________________________________________________________________________
195 AliMUONTrack::AliMUONTrack(const AliMUONTrack& track)
196   : TObject(track),
197     fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)),
198     fFitWithVertex(track.fFitWithVertex),
199     fVertexErrXY2(),
200     fFitWithMCS(track.fFitWithMCS),
201     fClusterWeightsNonBending(0x0),
202     fClusterWeightsBending(0x0),
203     fGlobalChi2(track.fGlobalChi2),
204     fImproved(track.fImproved),
205     fMatchTrigger(track.fMatchTrigger),
206     floTrgNum(track.floTrgNum),    
207     fChi2MatchTrigger(track.fChi2MatchTrigger),
208     fTrackID(track.fTrackID),
209     fTrackParamAtVertex(0x0),
210     fHitsPatternInTrigCh(track.fHitsPatternInTrigCh),
211     fLocalTrigger(track.fLocalTrigger)
212 {
213   ///copy constructor
214   
215   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
216   AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First();
217   while (trackParamAtCluster) {
218     new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster);
219     trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster);
220   }
221   
222   // copy vertex resolution square used during the tracking procedure
223   fVertexErrXY2[0] = track.fVertexErrXY2[0];
224   fVertexErrXY2[1] = track.fVertexErrXY2[1];
225   
226   // copy cluster weights matrices if any
227   if (track.fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
228   if (track.fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
229   
230   // copy track parameters at vertex if any
231   if (track.fTrackParamAtVertex) fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
232   
233 }
234
235   //__________________________________________________________________________
236 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
237 {
238   /// Asignment operator
239   // check assignement to self
240   if (this == &track)
241     return *this;
242
243   // base class assignement
244   TObject::operator=(track);
245
246   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
247   fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10);
248   AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First();
249   while (trackParamAtCluster) {
250     new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster);
251     trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster);
252   }
253   
254   // copy cluster weights matrix if any
255   if (track.fClusterWeightsNonBending) {
256     if (fClusterWeightsNonBending) {
257       fClusterWeightsNonBending->ResizeTo(*(track.fClusterWeightsNonBending));
258       *fClusterWeightsNonBending = *(track.fClusterWeightsNonBending);
259     } else fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
260   } else if (fClusterWeightsNonBending) {
261     delete fClusterWeightsNonBending;
262     fClusterWeightsNonBending = 0x0;
263   }
264   
265   // copy cluster weights matrix if any
266   if (track.fClusterWeightsBending) {
267     if (fClusterWeightsBending) {
268       fClusterWeightsBending->ResizeTo(*(track.fClusterWeightsBending));
269       *fClusterWeightsBending = *(track.fClusterWeightsBending);
270     } else fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
271   } else if (fClusterWeightsBending) {
272     delete fClusterWeightsBending;
273     fClusterWeightsBending = 0x0;
274   }
275   
276   // copy track parameters at vertex if any
277   if (track.fTrackParamAtVertex) {
278     if (fTrackParamAtVertex) *fTrackParamAtVertex = *(track.fTrackParamAtVertex);
279     else fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
280   } else if (fTrackParamAtVertex) {
281     delete fTrackParamAtVertex;
282     fTrackParamAtVertex = 0x0;
283   }
284   
285   fFitWithVertex      =  track.fFitWithVertex;
286   fVertexErrXY2[0]    =  track.fVertexErrXY2[0];
287   fVertexErrXY2[1]    =  track.fVertexErrXY2[1];
288   fFitWithMCS         =  track.fFitWithMCS;
289   fGlobalChi2         =  track.fGlobalChi2;
290   fImproved           =  track.fImproved;
291   fMatchTrigger       =  track.fMatchTrigger;
292   floTrgNum           =  track.floTrgNum;
293   fChi2MatchTrigger   =  track.fChi2MatchTrigger;
294   fTrackID            =  track.fTrackID; 
295   fHitsPatternInTrigCh = track.fHitsPatternInTrigCh;
296   fLocalTrigger        = track.fLocalTrigger;
297
298   return *this;
299 }
300
301   //__________________________________________________________________________
302 AliMUONTrack::~AliMUONTrack()
303 {
304   /// Destructor
305   delete fTrackParamAtCluster;
306   delete fClusterWeightsNonBending;
307   delete fClusterWeightsBending;
308   delete fTrackParamAtVertex;
309 }
310
311   //__________________________________________________________________________
312 void AliMUONTrack::Clear(Option_t* opt)
313 {
314   /// Clear arrays
315   fTrackParamAtCluster->Clear(opt);
316   delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
317   delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
318   delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
319 }
320
321   //__________________________________________________________________________
322 void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy)
323 {
324   /// Copy given track parameters into a new TrackParamAtCluster
325   /// Link parameters with the associated cluster
326   /// If copy=kTRUE: the cluster is copied then passed the trackParam which become its owner 
327   ///     otherwise: make sure to do not delete the cluster until it is used by the track
328   
329   // check chamber ID of the associated cluster
330   if (cluster.GetChamberId() < 0 || cluster.GetChamberId() > AliMUONConstants::NTrackingCh()) {
331     AliError(Form("Chamber ID of the associated cluster is not valid (ChamberId=%d)",cluster.GetChamberId()));
332     return;
333   }
334   
335   // check whether track parameters are given at the correct cluster z position
336   if (cluster.GetZ() != trackParam.GetZ()) {
337     AliError("track parameters are given at a different z position than the one of the associated cluster");
338     return;
339   }
340   
341   // add parameters to the array of track parameters
342   AliMUONTrackParam* trackParamAtCluster = new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(trackParam);
343   
344   // link parameters with the associated cluster or its copy
345   if (copy) {
346     AliMUONVCluster *clusterCopy = cluster.CreateCopy();
347     trackParamAtCluster->SetClusterPtr(clusterCopy, kTRUE);
348   } else trackParamAtCluster->SetClusterPtr(&cluster);
349 }
350
351   //__________________________________________________________________________
352 void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam)
353 {
354   /// Remove trackParam from the array of TrackParamAtCluster
355   if (!fTrackParamAtCluster->Remove(trackParam)) {
356     AliWarning("object to remove does not exist in array fTrackParamAtCluster");
357     return;
358   }
359   
360   fTrackParamAtCluster->Compress();
361 }
362
363   //__________________________________________________________________________
364 void AliMUONTrack::UpdateTrackParamAtCluster()
365 {
366   /// Update track parameters at each attached cluster
367   
368   if (GetNClusters() == 0) {
369     AliWarning("no cluster attached to the track");
370     return;
371   }
372   
373   AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First();
374   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam);
375   while (trackParamAtCluster) {
376     
377     // reset track parameters and their covariances
378     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
379     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
380     
381     // extrapolation to the given z
382     AliMUONTrackExtrap::ExtrapToZ(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ());
383     
384     // prepare next step
385     startingTrackParam = trackParamAtCluster;
386     trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster));
387   }
388
389 }
390
391   //__________________________________________________________________________
392 void AliMUONTrack::UpdateCovTrackParamAtCluster()
393 {
394   /// Update track parameters and their covariances at each attached cluster
395   /// Include effects of multiple scattering in chambers
396   
397   if (GetNClusters() == 0) {
398     AliWarning("no cluster attached to the track");
399     return;
400   }
401   
402   AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First();
403   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam);
404   Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() + 1;
405   Int_t currentChamber;
406   while (trackParamAtCluster) {
407     
408     // reset track parameters and their covariances
409     trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
410     trackParamAtCluster->SetZ(startingTrackParam->GetZ());
411     trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
412     
413     // add MCS effect
414     AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
415     
416     // add MCS in missing chambers if any
417     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
418     while (currentChamber > expectedChamber) {
419       // extrapolation to the missing chamber
420       AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber));
421       // add MCS effect
422       AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
423       expectedChamber++;
424     }
425     
426     // extrapolation to the z of the current cluster
427     AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ());
428     
429     // prepare next step
430     expectedChamber = currentChamber + 1;
431     startingTrackParam = trackParamAtCluster;
432     trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster));
433   }
434   
435 }
436
437   //__________________________________________________________________________
438 Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
439 {
440   /// Compute each cluster contribution to the chi2 of the track
441   /// accounting for multiple scattering or not according to the flag
442   /// - Also recompute the weight matrices of the attached clusters if accountForMCS=kTRUE
443   /// - Assume that track parameters at each cluster are corrects
444   /// - Return kFALSE if computation failed
445   
446   if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
447       
448     // Compute MCS covariance matrix only once
449     Int_t nClusters = GetNClusters();
450     TMatrixD mcsCovariances(nClusters,nClusters);
451     ComputeMCSCovariances(mcsCovariances);
452     
453     // Make sure cluster weights are consistent with following calculations
454     if (!ComputeClusterWeights(&mcsCovariances)) {
455       AliWarning("cannot take into account the multiple scattering effects");
456       return ComputeLocalChi2(kFALSE);
457     }
458     
459     // Compute chi2 of the track
460     Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
461     if (globalChi2 < 0.) return kFALSE;
462     
463     // Loop over removable clusters and compute their local chi2
464     AliMUONTrackParam* trackParamAtCluster1;
465     AliMUONVCluster *cluster, *discardedCluster;
466     Int_t iCluster1, iCluster2, iCurrentCluster1, iCurrentCluster2;
467     TMatrixD ClusterWeightsNB(nClusters-1,nClusters-1);
468     TMatrixD ClusterWeightsB(nClusters-1,nClusters-1);
469     Double_t *dX = new Double_t[nClusters-1];
470     Double_t *dY = new Double_t[nClusters-1];
471     Double_t globalChi2b;
472     AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
473     while (trackParamAtCluster) {
474       
475       discardedCluster = trackParamAtCluster->GetClusterPtr();
476       
477       // Recompute cluster weights without the current cluster
478       if (!ComputeClusterWeights(ClusterWeightsNB, ClusterWeightsB, &mcsCovariances, discardedCluster)) {
479         AliWarning("cannot take into account the multiple scattering effects");
480         delete [] dX;
481         delete [] dY;
482         return ComputeLocalChi2(kFALSE);
483       }
484       
485       // Compute track chi2 without the current cluster
486       globalChi2b = 0.;
487       iCurrentCluster1 = 0;
488       for (iCluster1 = 0; iCluster1 < nClusters ; iCluster1++) { 
489         trackParamAtCluster1 = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
490         cluster = trackParamAtCluster1->GetClusterPtr();
491         
492         if (cluster == discardedCluster) continue;
493         
494         // Compute and save residuals
495         dX[iCurrentCluster1] = cluster->GetX() - trackParamAtCluster1->GetNonBendingCoor();
496         dY[iCurrentCluster1] = cluster->GetY() - trackParamAtCluster1->GetBendingCoor();
497         
498         iCurrentCluster2 = 0;
499         for (iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
500           cluster = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
501           
502           if (cluster == discardedCluster) continue;
503           
504           // Add contribution from covariances
505           globalChi2b += (ClusterWeightsNB(iCurrentCluster1, iCurrentCluster2) +
506                           ClusterWeightsNB(iCurrentCluster2, iCurrentCluster1)) * dX[iCurrentCluster1] * dX[iCurrentCluster2] +
507                          (ClusterWeightsB(iCurrentCluster1, iCurrentCluster2) +
508                           ClusterWeightsB(iCurrentCluster2, iCurrentCluster1)) * dY[iCurrentCluster1] * dY[iCurrentCluster2];
509           
510           iCurrentCluster2++;
511         }
512         
513         // Add contribution from variances
514         globalChi2b += ClusterWeightsNB(iCurrentCluster1, iCurrentCluster1) * dX[iCurrentCluster1] * dX[iCurrentCluster1] +
515                        ClusterWeightsB(iCurrentCluster1, iCurrentCluster1) * dY[iCurrentCluster1] * dY[iCurrentCluster1];
516         
517         iCurrentCluster1++;
518       }
519
520       // Set local chi2
521       trackParamAtCluster->SetLocalChi2(globalChi2 - globalChi2b);
522       
523       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
524     }
525     
526     delete [] dX;
527     delete [] dY;
528     
529   } else { // without multiple scattering effects
530     
531     AliMUONVCluster *discardedCluster;
532     Double_t dX, dY;
533     AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
534     while (trackParamAtCluster) {
535       
536       discardedCluster = trackParamAtCluster->GetClusterPtr();
537       
538       // Compute residuals
539       dX = discardedCluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
540       dY = discardedCluster->GetY() - trackParamAtCluster->GetBendingCoor();
541       
542       // Set local chi2
543       trackParamAtCluster->SetLocalChi2(dX * dX / discardedCluster->GetErrX2() + dY * dY / discardedCluster->GetErrY2());
544     
545       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
546     }
547   
548   }
549   
550   return kTRUE;
551   
552 }
553
554   //__________________________________________________________________________
555 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
556 {
557   /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
558   /// - Assume that track parameters at each cluster are corrects
559   /// - Assume the cluster weights matrices are corrects
560   /// - Return negative value if chi2 computation failed
561   
562   Double_t chi2 = 0.;
563   
564   if (accountForMCS) {
565     
566     // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
567     if (!fClusterWeightsNonBending || !fClusterWeightsBending) {
568       AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
569       return ComputeGlobalChi2(kFALSE);
570     }
571     Int_t nClusters = GetNClusters();
572     if (fClusterWeightsNonBending->GetNrows() != nClusters || fClusterWeightsBending->GetNcols() != nClusters) {
573       AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
574       return ComputeGlobalChi2(kFALSE);
575     }
576     
577     // Compute chi2
578     AliMUONVCluster *cluster;
579     Double_t *dX = new Double_t[nClusters];
580     Double_t *dY = new Double_t[nClusters];
581     AliMUONTrackParam* trackParamAtCluster;
582     for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
583       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
584       cluster = trackParamAtCluster->GetClusterPtr();
585       dX[iCluster1] = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
586       dY[iCluster1] = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
587       for (Int_t iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
588         chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster2) + (*fClusterWeightsNonBending)(iCluster2, iCluster1)) * dX[iCluster1] * dX[iCluster2] +
589                 ((*fClusterWeightsBending)(iCluster1, iCluster2) + (*fClusterWeightsBending)(iCluster2, iCluster1)) * dY[iCluster1] * dY[iCluster2];
590       }
591       chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster1) * dX[iCluster1] * dX[iCluster1]) +
592               ((*fClusterWeightsBending)(iCluster1, iCluster1) * dY[iCluster1] * dY[iCluster1]);
593     }
594     delete [] dX;
595     delete [] dY;
596     
597   } else {
598     
599     AliMUONVCluster *cluster;
600     Double_t dX, dY;
601     AliMUONTrackParam* trackParamAtCluster;
602     Int_t nClusters = GetNClusters();
603     for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) { 
604       trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
605       cluster = trackParamAtCluster->GetClusterPtr();
606       dX = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
607       dY = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
608       chi2 += dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
609     }
610     
611   }
612   
613   return chi2;
614   
615 }
616
617   //__________________________________________________________________________
618 Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances)
619 {
620   /// Compute the weight matrices of the attached clusters, in non bending and bending direction,
621   /// accounting for multiple scattering correlations and cluster resolution
622   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
623   /// - Assume that track parameters at each cluster are corrects
624   /// - Return kFALSE if computation failed
625   
626   // Alocate memory
627   Int_t nClusters = GetNClusters();
628   if (!fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(nClusters,nClusters);
629   if (!fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(nClusters,nClusters);
630   
631   // Compute weights matrices
632   if (!ComputeClusterWeights(*fClusterWeightsNonBending, *fClusterWeightsBending, mcsCovariances)) return kFALSE;
633   
634   return kTRUE;
635   
636 }
637
638   //__________________________________________________________________________
639 Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& ClusterWeightsNB, TMatrixD& ClusterWeightsB,
640                                            TMatrixD* mcsCovariances, AliMUONVCluster* discardedCluster) const
641 {
642   /// Compute the weight matrices, in non bending and bending direction,
643   /// of the other attached clusters assuming the discarded one does not exist
644   /// accounting for multiple scattering correlations and cluster resolution
645   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
646   /// - Return kFALSE if computation failed
647   
648   // Check MCS covariance matrix and recompute it if need
649   Int_t nClusters = GetNClusters();
650   Bool_t deleteMCSCov = kFALSE;
651   if (!mcsCovariances) {
652     mcsCovariances = new TMatrixD(nClusters,nClusters);
653     deleteMCSCov = kTRUE;
654     ComputeMCSCovariances(*mcsCovariances);
655   }
656   
657   // Resize the weights matrices; alocate memory
658   if (discardedCluster) {
659     ClusterWeightsNB.ResizeTo(nClusters-1,nClusters-1);
660     ClusterWeightsB.ResizeTo(nClusters-1,nClusters-1);
661   } else {
662     ClusterWeightsNB.ResizeTo(nClusters,nClusters);
663     ClusterWeightsB.ResizeTo(nClusters,nClusters);
664   }
665   
666   // Define variables
667   AliMUONVCluster *cluster1, *cluster2;
668   Int_t iCurrentCluster1, iCurrentCluster2;
669   
670   // Compute the covariance matrices
671   iCurrentCluster1 = 0;
672   for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
673     cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
674     
675     if (cluster1 == discardedCluster) continue;
676     
677     // Loop over next clusters
678     iCurrentCluster2 = iCurrentCluster1;
679     for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
680       cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
681       
682       if (cluster2 == discardedCluster) continue;
683       
684       // Fill with MCS covariances
685       ClusterWeightsNB(iCurrentCluster1, iCurrentCluster2) = (*mcsCovariances)(iCluster1,iCluster2);
686       
687       // Equal contribution from multiple scattering in non bending and bending directions
688       ClusterWeightsB(iCurrentCluster1, iCurrentCluster2) = ClusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
689       
690       // Add contribution from cluster resolution to diagonal element and symmetrize the matrix
691       if (iCurrentCluster1 == iCurrentCluster2) {
692         
693         // In non bending plane
694         ClusterWeightsNB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrX2();
695         // In bending plane
696         ClusterWeightsB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrY2();
697         
698       } else {
699         
700         // In non bending plane
701         ClusterWeightsNB(iCurrentCluster2, iCurrentCluster1) = ClusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
702         // In bending plane
703         ClusterWeightsB(iCurrentCluster2, iCurrentCluster1) = ClusterWeightsB(iCurrentCluster1, iCurrentCluster2);
704         
705       }
706       
707       iCurrentCluster2++;
708     }
709     
710     iCurrentCluster1++;
711   }
712     
713   // Inversion of covariance matrices to get the weights
714   if (ClusterWeightsNB.Determinant() != 0 && ClusterWeightsB.Determinant() != 0) {
715     ClusterWeightsNB.Invert();
716     ClusterWeightsB.Invert();
717   } else {
718     AliWarning(" Determinant = 0");
719     ClusterWeightsNB.ResizeTo(0,0);
720     ClusterWeightsB.ResizeTo(0,0);
721     if(deleteMCSCov) delete mcsCovariances;
722     return kFALSE;
723   }
724   
725   if(deleteMCSCov) delete mcsCovariances;
726   
727   return kTRUE;
728   
729 }
730
731   //__________________________________________________________________________
732 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
733 {
734   /// Compute the multiple scattering covariance matrix
735   /// (assume that track parameters at each cluster are corrects)
736   
737   // Reset the size of the covariance matrix if needed
738   Int_t nClusters = GetNClusters();
739   if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters);
740   
741   // Define variables
742   Int_t nChambers = AliMUONConstants::NTrackingCh();
743   AliMUONTrackParam* trackParamAtCluster;
744   AliMUONTrackParam extrapTrackParam;
745   Int_t currentChamber = 0, expectedChamber = 0, size = 0;
746   Double_t *mcsAngle2 = new Double_t[2*nChambers];
747   Double_t *zMCS = new Double_t[2*nChambers];
748   Int_t *indices = new Int_t[2*nClusters];
749   
750   // Compute multiple scattering dispersion angle at each chamber
751   // and save the z position where it is calculated
752   for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
753     trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
754     
755     // look for missing chambers if any
756     currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
757     while (currentChamber > expectedChamber) {
758       
759       // Save the z position where MCS dispersion is calculated
760       zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
761       
762       // Do not take into account MCS in chambers prior the first cluster
763       if (iCluster > 0) {
764         
765         // Get track parameters at missing chamber z
766         extrapTrackParam = *trackParamAtCluster;
767         AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]);
768         
769         // Save multiple scattering dispersion angle in missing chamber
770         mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
771         
772       } else mcsAngle2[size] = 0.;
773       
774       expectedChamber++;
775       size++;
776     }
777     
778     // Save z position where MCS dispersion is calculated
779     zMCS[size] = trackParamAtCluster->GetZ();
780     
781     // Save multiple scattering dispersion angle in current chamber
782     mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
783     
784     // Save indice in zMCS array corresponding to the current cluster
785     indices[iCluster] = size;
786     
787     expectedChamber = currentChamber + 1;
788     size++;
789   }
790   
791   // complete array of z if last cluster is on the last but one chamber
792   if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
793   
794   // Compute the covariance matrix
795   for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { 
796     
797     for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
798       
799       // Initialization to 0 (diagonal plus upper triangular part)
800       mcsCovariances(iCluster1,iCluster2) = 0.;
801       
802       // Compute contribution from multiple scattering in upstream chambers
803       for (Int_t k = 0; k < indices[iCluster1]; k++) {  
804         mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k];
805       }
806       
807       // Symetrize the matrix
808       mcsCovariances(iCluster2,iCluster1) = mcsCovariances(iCluster1,iCluster2);
809     }
810     
811   }
812     
813   delete [] mcsAngle2;
814   delete [] zMCS;
815   delete [] indices;
816   
817 }
818
819   //__________________________________________________________________________
820 Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track) const
821 {
822   /// Returns the number of clusters in common between the current track ("this")
823   /// and the track pointed to by "track".
824   Int_t clustersInCommon = 0;
825   AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
826   // Loop over clusters of first track
827   trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First();
828   while (trackParamAtCluster1) {
829     // Loop over clusters of second track
830     trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First();
831     while (trackParamAtCluster2) {
832       // Increment "clustersInCommon" if both trackParamAtCluster1 & 2 point to the same cluster
833       if ((trackParamAtCluster1->GetClusterPtr()) == (trackParamAtCluster2->GetClusterPtr())) {
834         clustersInCommon++;
835         break;
836       }
837       trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2);
838     } // trackParamAtCluster2
839     trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1);
840   } // trackParamAtCluster1
841   return clustersInCommon;
842 }
843
844   //__________________________________________________________________________
845 Double_t AliMUONTrack::GetNormalizedChi2() const
846 {
847   /// return the chi2 value divided by the number of degrees of freedom (or 1.e10 if ndf < 0)
848   
849   Double_t numberOfDegFree = (2. * GetNClusters() - 5.);
850   if (numberOfDegFree > 0.) return fGlobalChi2 / numberOfDegFree;
851   else return 1.e10;
852 }
853
854   //__________________________________________________________________________
855 Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack *track, Double_t sigma2Cut) const
856 {
857   /// for each chamber: return kTRUE (kFALSE) if clusters are compatible (not compatible)
858   AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
859   AliMUONVCluster *cluster1, *cluster2;
860   Double_t chi2, dX, dY, dZ;
861   Double_t chi2Max = sigma2Cut * sigma2Cut;
862   Double_t dZMax = 1.; // 1 cm
863   
864   Bool_t *compatibleCluster = new Bool_t[AliMUONConstants::NTrackingCh()]; 
865   for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE;
866
867   // Loop over clusters of first track
868   trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First();
869   while (trackParamAtCluster1) {
870     
871     cluster1 = trackParamAtCluster1->GetClusterPtr();
872     
873     // Loop over clusters of second track
874     trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First();
875     while (trackParamAtCluster2) {
876       
877       cluster2 = trackParamAtCluster2->GetClusterPtr();
878       
879       //prepare next step
880       trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2);
881       
882       // z direction
883       dZ = cluster1->GetZ() - cluster2->GetZ();
884       if (dZ > dZMax) continue;
885       
886       // non bending direction
887       dX = cluster1->GetX() - cluster2->GetX();
888       chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2());
889       if (chi2 > chi2Max) continue;
890       
891       // bending direction
892       dY = cluster1->GetY() - cluster2->GetY();
893       chi2 = dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2());
894       if (chi2 > chi2Max) continue;
895       
896       compatibleCluster[cluster1->GetChamberId()] = kTRUE;
897       break;
898     }
899     
900     trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1);
901   }
902   
903   return compatibleCluster;
904 }
905
906 //__________________________________________________________________________
907 AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex()
908 {
909   /// return reference to track parameters at vertex (create it before if needed)
910   if (!fTrackParamAtVertex) fTrackParamAtVertex = new AliMUONTrackParam();
911   return fTrackParamAtVertex;
912 }
913
914 //__________________________________________________________________________
915 void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam)
916 {
917   /// set track parameters at vertex
918   if (trackParam == 0x0) return;
919   if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam;
920   else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam);
921 }
922
923 //__________________________________________________________________________
924 void AliMUONTrack::RecursiveDump() const
925 {
926   /// Recursive dump of AliMUONTrack, i.e. with dump of trackParamAtCluster and attached clusters
927   AliMUONTrackParam *trackParamAtCluster;
928   AliMUONVCluster *cluster;
929   cout << "Recursive dump of Track: " << this << endl;
930   // Track
931   this->Dump();
932   for (Int_t iCluster = 0; iCluster < GetNClusters(); iCluster++) {
933     trackParamAtCluster = (AliMUONTrackParam*) ((*fTrackParamAtCluster)[iCluster]);
934     // trackParamAtCluster
935     cout << "trackParamAtCluster: " << trackParamAtCluster << " (index: " << iCluster << ")" << endl;
936     trackParamAtCluster->Dump();
937     cluster = trackParamAtCluster->GetClusterPtr();
938     // cluster
939     cout << "cluster: " << cluster << endl;
940     cluster->Print();
941   }
942   return;
943 }
944   
945 //_____________________________________________-
946 void AliMUONTrack::Print(Option_t*) const
947 {
948   /// Printing Track information 
949
950   cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNClusters() << 
951       ", Match2Trig=" << setw(1) << GetMatchTrigger()  << 
952       ", LoTrgNum=" << setw(3) << GetLoTrgNum()  << 
953     ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) <<  GetChi2MatchTrigger();
954   cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl;
955   fTrackParamAtCluster->First()->Print("FULL");
956 }
957
958 //__________________________________________________________________________
959 void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt)
960 {
961   /// pack the local trigger information and store
962
963   if (loCirc < 0) return;
964
965   fLocalTrigger = 0;
966   fLocalTrigger += loCirc;
967   fLocalTrigger += loStripX << 8;
968   fLocalTrigger += loStripY << 13;
969   fLocalTrigger += loDev    << 17;
970   fLocalTrigger += loLpt    << 22;
971   fLocalTrigger += loHpt    << 24;
972
973 }
974