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