Adding comment lines to class description needed for Root documentation,
[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(234)
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(234)
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(AliMUONConstants::NTrackingCh(),AliMUONConstants::NTrackingCh());
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   
604   return kTRUE;
605   
606 }
607
608   //__________________________________________________________________________
609 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
610 {
611   /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
612   /// - Assume that track parameters at each hit are corrects
613   /// - Assume the hits weights matrices are corrects
614   /// - Return negative value if chi2 computation failed
615   
616   // Check hits (if the first one exist, assume that the other ones exit too!)
617   AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
618   if (!trackParamAtHit->GetHitForRecPtr()) {
619     AliWarning("hit is missing");
620     return -1.;
621   }
622   
623   Double_t chi2 = 0.;
624   
625   if (accountForMCS) {
626     
627     // Check the weight matrices
628     Bool_t weightsAvailable = kTRUE;
629     if (!fHitWeightsNonBending || !fHitWeightsBending) weightsAvailable = kFALSE;
630     else if (fHitWeightsNonBending->GetNrows() != fNTrackHits || fHitWeightsNonBending->GetNcols() != fNTrackHits ||
631              fHitWeightsBending->GetNrows()    != fNTrackHits || fHitWeightsBending->GetNcols()    != fNTrackHits) weightsAvailable = kFALSE;
632     
633     // if weight matrices are not available compute chi2 without MCS
634     if (!weightsAvailable) {
635       AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
636       return ComputeGlobalChi2(kFALSE);
637     }
638     
639     // Compute chi2
640     AliMUONHitForRec *hitForRec;
641     Double_t *dX = new Double_t[fNTrackHits];
642     Double_t *dY = new Double_t[fNTrackHits];
643     Int_t hitNumber1, hitNumber2;
644     for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) { 
645       trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
646       hitForRec = trackParamAtHit->GetHitForRecPtr();
647       dX[hitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
648       dY[hitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
649       for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
650         chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber2) + (*fHitWeightsNonBending)(hitNumber2, hitNumber1)) * dX[hitNumber1] * dX[hitNumber2] +
651                 ((*fHitWeightsBending)(hitNumber1, hitNumber2) + (*fHitWeightsBending)(hitNumber2, hitNumber1)) * dY[hitNumber1] * dY[hitNumber2];
652       }
653       chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber1) * dX[hitNumber1] * dX[hitNumber1]) +
654               ((*fHitWeightsBending)(hitNumber1, hitNumber1) * dY[hitNumber1] * dY[hitNumber1]);
655     }
656     delete [] dX;
657     delete [] dY;
658     
659   } else {
660     
661     AliMUONHitForRec *hitForRec;
662     Double_t dX, dY;
663     for (Int_t hitNumber = 0; hitNumber < fNTrackHits ; hitNumber++) { 
664       trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
665       hitForRec = trackParamAtHit->GetHitForRecPtr();
666       dX = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
667       dY = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
668       chi2 += dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
669     }
670     
671   }
672   
673   return chi2;
674   
675 }
676
677   //__________________________________________________________________________
678 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD* mcsCovariances)
679 {
680   /// Compute the weight matrices of the attached hits, in non bending and bending direction,
681   /// accounting for multiple scattering correlations and hits resolution
682   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
683   /// - Assume that track parameters at each hit are corrects
684   /// - Return kFALSE if computation failed
685   
686   // Alocate memory
687   if (!fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(fNTrackHits,fNTrackHits);
688   if (!fHitWeightsBending) fHitWeightsBending = new TMatrixD(fNTrackHits,fNTrackHits);
689   
690   // Check hits (if the first one exist, assume that the other ones exit too!)
691   if (!((AliMUONTrackParam*) fTrackParamAtHit->First())->GetHitForRecPtr()) {
692     AliWarning("hit is missing");
693     fHitWeightsNonBending->ResizeTo(0,0);
694     fHitWeightsBending->ResizeTo(0,0);
695     return kFALSE;
696   }
697   
698   // Compute weights matrices
699   if (!ComputeHitWeights(*fHitWeightsNonBending, *fHitWeightsBending, mcsCovariances)) return kFALSE;
700   
701   return kTRUE;
702   
703 }
704
705   //__________________________________________________________________________
706 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD& hitWeightsNB, TMatrixD& hitWeightsB, TMatrixD* mcsCovariances, AliMUONHitForRec* discardedHit) const
707 {
708   /// Compute the weight matrices, in non bending and bending direction,
709   /// of the other attached hits assuming the discarded one does not exist
710   /// accounting for multiple scattering correlations and hits resolution
711   /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
712   /// - Return kFALSE if computation failed
713   
714   // Check MCS covariance matrix and recompute it if need
715   Bool_t deleteMCSCov = kFALSE;
716   if (!mcsCovariances) {
717     
718     // build MCS covariance matrix
719     mcsCovariances = new TMatrixD(AliMUONConstants::NTrackingCh(),AliMUONConstants::NTrackingCh());
720     deleteMCSCov = kTRUE;
721     ComputeMCSCovariances(*mcsCovariances);
722     
723   } else {
724     
725     // check MCS covariance matrix size
726     if (mcsCovariances->GetNrows() != AliMUONConstants::NTrackingCh() || mcsCovariances->GetNcols() != AliMUONConstants::NTrackingCh()) {
727       ComputeMCSCovariances(*mcsCovariances);
728     }
729     
730   }
731   
732   // Resize the weights matrices; alocate memory
733   if (discardedHit) {
734     hitWeightsNB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
735     hitWeightsB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
736   } else {
737     hitWeightsNB.ResizeTo(fNTrackHits,fNTrackHits);
738     hitWeightsB.ResizeTo(fNTrackHits,fNTrackHits);
739   }
740   
741   // Define variables
742   AliMUONHitForRec *hitForRec1, *hitForRec2;
743   Int_t chamber1, chamber2, currentHitNumber1, currentHitNumber2;
744   
745   // Compute the covariance matrices
746   currentHitNumber1 = 0;
747   for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) { 
748     hitForRec1 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1))->GetHitForRecPtr();
749     
750     if (hitForRec1 == discardedHit) continue;
751     
752     chamber1 = hitForRec1->GetChamberNumber();
753     
754     // Loop over next hits
755     currentHitNumber2 = currentHitNumber1;
756     for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
757       hitForRec2 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
758       
759       if (hitForRec2 == discardedHit) continue;
760       
761       chamber2 = hitForRec2->GetChamberNumber();
762     
763       // Fill with MCS covariances
764       hitWeightsNB(currentHitNumber1, currentHitNumber2) = (*mcsCovariances)(chamber1,chamber2);
765       
766       // Equal contribution from multiple scattering in non bending and bending directions
767       hitWeightsB(currentHitNumber1, currentHitNumber2) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
768       
769       // Add contribution from hit resolution to diagonal element and symmetrize the matrix
770       if (currentHitNumber1 == currentHitNumber2) {
771         
772         // In non bending plane
773         hitWeightsNB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetNonBendingReso2();
774         // In bending plane
775         hitWeightsB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetBendingReso2();
776         
777       } else {
778         
779         // In non bending plane
780         hitWeightsNB(currentHitNumber2, currentHitNumber1) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
781         // In bending plane
782         hitWeightsB(currentHitNumber2, currentHitNumber1) = hitWeightsB(currentHitNumber1, currentHitNumber2);
783         
784       }
785       
786       currentHitNumber2++;
787     }
788     
789     currentHitNumber1++;
790   }
791     
792   // Inversion of covariance matrices to get the weights
793   if (hitWeightsNB.Determinant() != 0 && hitWeightsB.Determinant() != 0) {
794     hitWeightsNB.Invert();
795     hitWeightsB.Invert();
796   } else {
797     AliWarning(" Determinant = 0");
798     hitWeightsNB.ResizeTo(0,0);
799     hitWeightsB.ResizeTo(0,0);
800     if(deleteMCSCov) delete mcsCovariances;
801     return kFALSE;
802   }
803   
804   if(deleteMCSCov) delete mcsCovariances;
805   
806   return kTRUE;
807   
808 }
809
810   //__________________________________________________________________________
811 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
812 {
813   /// Compute the multiple scattering covariance matrix
814   /// - Assume that track parameters at each hit are corrects
815   /// - Return kFALSE if computation failed
816   
817   // Make sure the size of the covariance matrix is correct
818   Int_t nChambers = AliMUONConstants::NTrackingCh();
819   mcsCovariances.ResizeTo(nChambers,nChambers);
820   
821   // check for too many track hits
822   if (fNTrackHits > nChambers) {
823     AliWarning("more than 1 hit per chamber!!");
824     mcsCovariances.Zero();
825     return;
826   }
827   
828   // Define variables
829   AliMUONTrackParam* trackParamAtHit;
830   AliMUONHitForRec *hitForRec;
831   AliMUONTrackParam extrapTrackParam;
832   Int_t currentChamber, expectedChamber;
833   Double_t *mcsAngle2 = new Double_t[nChambers];
834   Double_t *zMCS = new Double_t[nChambers];
835   
836   // Compute multiple scattering dispersion angle at each chamber
837   // and save the z position where it is calculated
838   currentChamber = 0;
839   expectedChamber = 0;
840   for (Int_t hitNumber = 0; hitNumber < fNTrackHits; hitNumber++) {
841     trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
842     hitForRec = trackParamAtHit->GetHitForRecPtr();
843     
844     // look for missing chambers if any
845     currentChamber = hitForRec->GetChamberNumber();
846     while (currentChamber > expectedChamber) {
847       
848       // Save the z position where MCS dispersion is calculated
849       zMCS[expectedChamber] = AliMUONConstants::DefaultChamberZ(expectedChamber);
850       
851       // Do not take into account MCS in chambers prior the first hit
852       if (hitNumber > 0) {
853         
854         // Get track parameters at missing chamber z
855         extrapTrackParam = *trackParamAtHit;
856         AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[expectedChamber]);
857         
858         // Save multiple scattering dispersion angle in missing chamber
859         mcsAngle2[expectedChamber] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
860         
861       } else mcsAngle2[expectedChamber] = 0.;
862       
863       expectedChamber++;
864     }
865     
866     // Save z position where MCS dispersion is calculated
867     zMCS[currentChamber] = trackParamAtHit->GetZ();
868     
869     // Save multiple scattering dispersion angle in current chamber
870     mcsAngle2[currentChamber] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.);
871     
872     expectedChamber++;
873   }
874   
875   // complete array of z if last hit is on the last but one chamber
876   if (currentChamber != nChambers-1) zMCS[nChambers-1] = AliMUONConstants::DefaultChamberZ(nChambers-1);
877   
878   
879   // Compute the covariance matrix
880   for (Int_t chamber1 = 0; chamber1 < nChambers; chamber1++) { 
881     
882     for (Int_t chamber2 = chamber1; chamber2 < nChambers; chamber2++) {
883       
884       // Initialization to 0 (diagonal plus upper triangular part)
885       mcsCovariances(chamber1, chamber2) = 0.;
886       
887       // Compute contribution from multiple scattering in upstream chambers
888       for (currentChamber = 0; currentChamber < chamber1; currentChamber++) {   
889         mcsCovariances(chamber1, chamber2) += (zMCS[chamber1] - zMCS[currentChamber]) * (zMCS[chamber2] - zMCS[currentChamber]) * mcsAngle2[currentChamber];
890       }
891       
892       // Symetrize the matrix
893       mcsCovariances(chamber2, chamber1) = mcsCovariances(chamber1, chamber2);
894     }
895     
896   }
897     
898   delete [] mcsAngle2;
899   delete [] zMCS;
900   
901 }
902
903   //__________________________________________________________________________
904 Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* track) const
905 {
906   /// Returns the number of hits in common between the current track ("this")
907   /// and the track pointed to by "track".
908   Int_t hitsInCommon = 0;
909   AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2;
910   // Loop over hits of first track
911   trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->First();
912   while (trackParamAtHit1) {
913     // Loop over hits of second track
914     trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->First();
915     while (trackParamAtHit2) {
916       // Increment "hitsInCommon" if both TrackParamAtHits point to the same HitForRec
917       if ((trackParamAtHit1->GetHitForRecPtr()) == (trackParamAtHit2->GetHitForRecPtr())) {
918         hitsInCommon++;
919         break;
920       }
921       trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->After(trackParamAtHit2);
922     } // trackParamAtHit2
923     trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->After(trackParamAtHit1);
924   } // trackParamAtHit1
925   return hitsInCommon;
926 }
927
928   //__________________________________________________________________________
929 Double_t AliMUONTrack::GetNormalizedChi2() const
930 {
931   /// return the chi2 value divided by the number of degrees of freedom (or 1.e10 if ndf < 0)
932   
933   Double_t numberOfDegFree = (2. * fNTrackHits - 5.);
934   if (numberOfDegFree > 0.) return fGlobalChi2 / numberOfDegFree;
935   else return 1.e10;
936 }
937
938   //__________________________________________________________________________
939 Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * track, Double_t sigma2Cut) const
940 {
941   /// Return kTRUE/kFALSE for each chamber if hit is compatible or not 
942   TClonesArray *hitArray, *thisHitArray;
943   AliMUONHitForRec *hit, *thisHit;
944   Int_t chamberNumber;
945   Float_t deltaZ;
946   Float_t deltaZMax = 1.; // 1 cm
947   Float_t chi2 = 0;
948   Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()]; 
949
950   for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
951     nCompHit[ch] = kFALSE;
952   }
953
954   thisHitArray = this->GetHitForRecAtHit();
955
956   hitArray =  track->GetHitForRecAtHit();
957
958   for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
959     thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
960     chamberNumber = thisHit->GetChamberNumber();
961     if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue; 
962     nCompHit[chamberNumber] = kFALSE;
963     for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
964       hit = (AliMUONHitForRec*) hitArray->At(iH);
965       deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
966       chi2 = thisHit->NormalizedChi2WithHitForRec(hit,sigma2Cut); // set cut to 4 sigmas
967       if (chi2 < 3. && deltaZ < deltaZMax) {
968         nCompHit[chamberNumber] = kTRUE;
969         break;
970       }
971     }  
972   }
973   
974   return nCompHit;
975 }
976
977   //__________________________________________________________________________
978 void AliMUONTrack::RecursiveDump(void) const
979 {
980   /// Recursive dump of AliMUONTrack, i.e. with dump of TrackParamAtHit's and attached HitForRec's
981   AliMUONTrackParam *trackParamAtHit;
982   AliMUONHitForRec *hitForRec;
983   cout << "Recursive dump of Track: " << this << endl;
984   // Track
985   this->Dump();
986   for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
987     trackParamAtHit = (AliMUONTrackParam*) ((*fTrackParamAtHit)[trackHitIndex]);
988     // TrackHit
989     cout << "TrackParamAtHit: " << trackParamAtHit << " (index: " << trackHitIndex << ")" << endl;
990     trackParamAtHit->Dump();
991     hitForRec = trackParamAtHit->GetHitForRecPtr();
992     // HitForRec
993     cout << "HitForRec: " << hitForRec << endl;
994     hitForRec->Dump();
995   }
996   return;
997 }
998   
999 //_____________________________________________-
1000 void AliMUONTrack::Print(Option_t*) const
1001 {
1002   /// Printing Track information 
1003
1004   cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNTrackHits() << 
1005       ", Match2Trig=" << setw(1) << GetMatchTrigger()  << 
1006       ", LoTrgNum=" << setw(3) << GetLoTrgNum()  << 
1007     ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) <<  GetChi2MatchTrigger();
1008   cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl;
1009   GetTrackParamAtHit()->First()->Print("FULL");
1010 }
1011
1012 //__________________________________________________________________________
1013 void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt)
1014 {
1015   /// pack the local trigger information and store
1016
1017   if (loCirc < 0 || loCirc > 233) return;
1018
1019   fLocalTrigger = 0;
1020   fLocalTrigger += loCirc;
1021   fLocalTrigger += loStripX << 8;
1022   fLocalTrigger += loStripY << 13;
1023   fLocalTrigger += loDev    << 17;
1024   fLocalTrigger += loLpt    << 22;
1025   fLocalTrigger += loHpt    << 24;
1026
1027 }
1028