AliMUONEventReconstructor package:
[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 /*
17 $Log$
18 Revision 1.4  2000/06/30 10:15:48  gosset
19 Changes to EventReconstructor...:
20 precision fit with multiple Coulomb scattering;
21 extrapolation to vertex with Branson correction in absorber (JPC)
22
23 Revision 1.3  2000/06/25 13:23:28  hristov
24 stdlib.h needed for non-Linux compilation
25
26 Revision 1.2  2000/06/15 07:58:48  morsch
27 Code from MUON-dev joined
28
29 Revision 1.1.2.3  2000/06/12 10:11:34  morsch
30 Dummy copy constructor and assignment operator added
31
32 Revision 1.1.2.2  2000/06/09 12:58:05  gosset
33 Removed comment beginnings in Log sections of .cxx files
34 Suppressed most violations of coding rules
35
36 Revision 1.1.2.1  2000/06/07 14:44:53  gosset
37 Addition of files for track reconstruction in C++
38 */
39
40 //__________________________________________________________________________
41 //
42 // Reconstructed track in ALICE dimuon spectrometer
43 //__________________________________________________________________________
44
45 #include "AliMUONTrack.h"
46
47 #include <iostream.h>
48
49 #include <TClonesArray.h>
50 #include <TMath.h>
51 #include <TMatrix.h>
52 #include <TVirtualFitter.h>
53
54 #include "AliMUONEventReconstructor.h" 
55 #include "AliMUONHitForRec.h" 
56 #include "AliMUONSegment.h" 
57 #include "AliMUONTrackHit.h"
58
59 #include <stdlib.h>
60
61 // Functions to be minimized with Minuit
62 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
63 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
64
65 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
66
67 ClassImp(AliMUONTrack) // Class implementation in ROOT context
68
69 TVirtualFitter* AliMUONTrack::fgFitter = NULL; 
70
71   //__________________________________________________________________________
72 AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor)
73 {
74   // Constructor from two Segment's
75   fEventReconstructor = EventReconstructor; // link back to EventReconstructor
76   // memory allocation for the TClonesArray of reconstructed TrackHit's
77   fTrackHitsPtr = new  TClonesArray("AliMUONTrackHit", 10);
78   fNTrackHits = 0;
79   AddSegment(BegSegment); // add hits from BegSegment
80   AddSegment(EndSegment); // add hits from EndSegment
81   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
82   SetTrackParamAtVertex(); // set track parameters at vertex
83   // set fit conditions
84   fFitMCS = 0;
85   fFitNParam = 3;
86   fFitStart = 1;
87   return;
88 }
89
90   //__________________________________________________________________________
91 AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor)
92 {
93   // Constructor from one Segment and one HitForRec
94   fEventReconstructor = EventReconstructor; // link back to EventReconstructor
95   // memory allocation for the TClonesArray of reconstructed TrackHit's
96   fTrackHitsPtr = new  TClonesArray("AliMUONTrackHit", 10);
97   fNTrackHits = 0;
98   AddSegment(Segment); // add hits from Segment
99   AddHitForRec(HitForRec); // add HitForRec
100   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
101   SetTrackParamAtVertex(); // set track parameters at vertex
102   // set fit conditions
103   fFitMCS = 0;
104   fFitNParam = 3;
105   fFitStart = 1;
106   return;
107 }
108
109   //__________________________________________________________________________
110 AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
111 {
112 // Dummy copy constructor
113 }
114
115   //__________________________________________________________________________
116 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
117 {
118 // Dummy assignment operator
119     return *this;
120 }
121
122   //__________________________________________________________________________
123 void AliMUONTrack::SetFitMCS(Int_t FitMCS)
124 {
125   // Set multiple Coulomb scattering option for track fit "fFitMCS"
126   // from "FitMCS" argument: 0 without, 1 with
127   if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS;
128   // better implementation with enum(with, without) ????
129   else {
130     cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl;
131     cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
132     exit(0);
133   }
134   return;
135 }
136
137   //__________________________________________________________________________
138 void AliMUONTrack::SetFitNParam(Int_t FitNParam)
139 {
140   // Set number of parameters for track fit "fFitNParam" from "FitNParam":
141   // 3 for momentum, 5 for momentum and position
142   if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam;
143   else {
144     cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl;
145     cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl;
146     exit(0);
147   }
148   return;
149 }
150
151   //__________________________________________________________________________
152 void AliMUONTrack::SetFitStart(Int_t FitStart)
153 {
154   // Set multiple Coulomb scattering option for track fit "fFitStart"
155   // from "FitStart" argument: 0 without, 1 with
156   if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart;
157   // better implementation with enum(vertex, firstHit) ????
158   else {
159     cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl;
160     cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
161     exit(0);
162   }
163   return;
164 }
165
166   //__________________________________________________________________________
167 AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
168   // Get pointer to TrackParamAtFirstHit
169   return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
170
171   //__________________________________________________________________________
172 void AliMUONTrack::RecursiveDump(void)
173 {
174   // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's
175   AliMUONTrackHit *trackHit;
176   AliMUONHitForRec *hitForRec;
177   cout << "Recursive dump of Track: " << this << endl;
178   // Track
179   this->Dump();
180   for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
181     trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]);
182     // TrackHit
183     cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl;
184     trackHit->Dump();
185     hitForRec = trackHit->GetHitForRecPtr();
186     // HitForRec
187     cout << "HitForRec: " << hitForRec << endl;
188     hitForRec->Dump();
189   }
190   return;
191 }
192
193   //__________________________________________________________________________
194 void AliMUONTrack::Fit()
195 {
196   // Fit the current track ("this"),
197   // with or without multiple Coulomb scattering according to "fFitMCS",
198   // with the number of parameters given by "fFitNParam"
199   // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary),
200   // starting, according to "fFitStart",
201   // with track parameters at vertex or at the first TrackHit.
202   // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before
203   // by calling the corresponding Set methods.
204   Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
205   char parName[50];
206   AliMUONTrackParam *trackParam;
207   // Check if Minuit is initialized...
208   fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ???
209   fgFitter->Clear(); // necessary ???? probably yes
210   // how to reset the printout number at every fit ????
211   // is there any risk to leave it like that ????
212   // how to go faster ???? choice of Minuit parameters like EDM ????
213   // choice of function to be minimized according to fFitMCS
214   if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
215   else fgFitter->SetFCN(TrackChi2MCS);
216   arg[0] = 1;
217   fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
218   // Parameters according to "fFitStart"
219   // (should be a function to be used at every place where needed ????)
220   if (fFitStart == 0) trackParam = &fTrackParamAtVertex;
221   else trackParam = this->GetTrackParamAtFirstHit();
222   // set first 3 Minuit parameters
223   // could be tried with no limits for the search (min=max=0) ????
224   fgFitter->SetParameter(0, "InvBenP",
225                          trackParam->GetInverseBendingMomentum(),
226                          0.003, -0.4, 0.4);
227   fgFitter->SetParameter(1, "BenS",
228                          trackParam->GetBendingSlope(),
229                          0.001, -0.5, 0.5);
230   fgFitter->SetParameter(2, "NonBenS",
231                          trackParam->GetNonBendingSlope(),
232                          0.001, -0.5, 0.5);
233   if (fFitNParam == 5) {
234     // set last 2 Minuit parameters (no limits for the search: min=max=0)
235     fgFitter->SetParameter(3, "X",
236                            trackParam->GetNonBendingCoor(),
237                            0.03, 0.0, 0.0);
238     fgFitter->SetParameter(4, "Y",
239                            trackParam->GetBendingCoor(),
240                            0.10, 0.0, 0.0);
241   }
242   // search without gradient calculation in the function
243   fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
244   // minimization
245   fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
246   // exit from Minuit
247   fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
248   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
249   fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
250   fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
251   fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
252   if (fFitNParam == 5) {
253     fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
254     fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
255   }
256   // result of the fit into track parameters
257   trackParam->SetInverseBendingMomentum(invBenP);
258   trackParam->SetBendingSlope(benC);
259   trackParam->SetNonBendingSlope(nonBenC);
260   if (fFitNParam == 5) {
261     trackParam->SetNonBendingCoor(x);
262     trackParam->SetBendingCoor(y);
263   }
264 }
265
266   //__________________________________________________________________________
267 void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
268 {
269   // Add Segment
270   AddHitForRec(Segment->GetHitForRec1()); // 1st hit
271   AddHitForRec(Segment->GetHitForRec2()); // 2nd hit
272 }
273
274   //__________________________________________________________________________
275 void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
276 {
277   // Add HitForRec
278   new ((*fTrackHitsPtr)[fNTrackHits]) AliMUONTrackHit(HitForRec);
279   fNTrackHits++;
280 }
281
282   //__________________________________________________________________________
283 void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam)
284 {
285   // Set track parameters at TrackHit with index "indexHit"
286   // from the track parameters pointed to by "TrackParam".
287   AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]);
288   trackHit->SetTrackParam(TrackParam);
289 }
290
291   //__________________________________________________________________________
292 void AliMUONTrack::SetTrackParamAtVertex()
293 {
294   // Set track parameters at vertex.
295   // TrackHit's are assumed to be only in stations(1..) 4 and 5,
296   // and sorted according to increasing Z..
297   // Parameters are calculated from information in HitForRec's
298   // of first and last TrackHit's.
299   AliMUONTrackParam *trackParam =
300     &fTrackParamAtVertex; // pointer to track parameters
301   // Pointer to HitForRec of first TrackHit
302   AliMUONHitForRec *firstHit =
303     ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr();
304   // Pointer to HitForRec of last TrackHit
305   AliMUONHitForRec *lastHit =
306     ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr();
307   // Z difference between first and last hits
308   Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
309   // bending slope in stations(1..) 4 and 5
310   Double_t bendingSlope =
311     (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
312   trackParam->SetBendingSlope(bendingSlope);
313   // impact parameter
314   Double_t impactParam =
315     firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and  lastHit ????
316   // signed bending momentum
317   Double_t signedBendingMomentum =
318     fEventReconstructor->GetBendingMomentumFromImpactParam(impactParam);
319   trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum);
320   // bending slope at vertex
321   trackParam->
322     SetBendingSlope(bendingSlope +
323                     impactParam / fEventReconstructor->GetSimpleBPosition());
324   // non bending slope
325   Double_t nonBendingSlope =
326     (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ;
327   trackParam->SetNonBendingSlope(nonBendingSlope);
328   // vertex coordinates at (0,0,0)
329   trackParam->SetZ(0.0);
330   trackParam->SetBendingCoor(0.0);
331   trackParam->SetNonBendingCoor(0.0);
332 }
333
334   //__________________________________________________________________________
335 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
336 {
337   // Return the "Chi2" to be minimized with Minuit for track fitting,
338   // with "NParam" parameters
339   // and their current values in array pointed to by "Param".
340   // Assumes that the track hits are sorted according to increasing Z.
341   // Track parameters at each TrackHit are updated accordingly.
342   // Multiple Coulomb scattering is not taken into account
343   AliMUONTrack *trackBeingFitted;
344   AliMUONTrackHit* hit;
345   AliMUONTrackParam param1;
346   Int_t hitNumber;
347   Double_t zHit;
348   Chi2 = 0.0; // initialize Chi2
349   // copy of track parameters to be fitted
350   trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
351   if (trackBeingFitted->GetFitStart() == 0)
352     param1 = *(trackBeingFitted->GetTrackParamAtVertex());
353   else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
354   // Minuit parameters to be fitted into this copy
355   param1.SetInverseBendingMomentum(Param[0]);
356   param1.SetBendingSlope(Param[1]);
357   param1.SetNonBendingSlope(Param[2]);
358   if (NParam == 5) {
359     param1.SetNonBendingCoor(Param[3]);
360     param1.SetBendingCoor(Param[4]);
361   }
362   // Follow track through all planes of track hits
363   for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) {
364     hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
365     zHit = hit->GetHitForRecPtr()->GetZ();
366     // do something special if 2 hits with same Z ????
367     // security against infinite loop ????
368     (&param1)->ExtrapToZ(zHit); // extrapolation
369     hit->SetTrackParam(&param1);
370     // Increment Chi2
371     // done hit per hit, with hit resolution,
372     // and not with point and angle like in "reco_muon.F" !!!!
373     // Needs to add multiple scattering contribution ????
374     Double_t dX =
375       hit->GetHitForRecPtr()->GetNonBendingCoor() - (&param1)->GetNonBendingCoor();
376     Double_t dY =
377       hit->GetHitForRecPtr()->GetBendingCoor() - (&param1)->GetBendingCoor();
378     Chi2 =
379       Chi2 +
380       dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() +
381       dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
382   }
383 }
384
385   //__________________________________________________________________________
386 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
387 {
388   // Return the "Chi2" to be minimized with Minuit for track fitting,
389   // with "NParam" parameters
390   // and their current values in array pointed to by "Param".
391   // Assumes that the track hits are sorted according to increasing Z.
392   // Track parameters at each TrackHit are updated accordingly.
393   // Multiple Coulomb scattering is taken into account with covariance matrix.
394   AliMUONTrack *trackBeingFitted;
395   AliMUONTrackParam param1;
396   Chi2 = 0.0; // initialize Chi2
397   // copy of track parameters to be fitted
398   trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
399   if (trackBeingFitted->GetFitStart() == 0)
400     param1 = *(trackBeingFitted->GetTrackParamAtVertex());
401   else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
402   // Minuit parameters to be fitted into this copy
403   param1.SetInverseBendingMomentum(Param[0]);
404   param1.SetBendingSlope(Param[1]);
405   param1.SetNonBendingSlope(Param[2]);
406   if (NParam == 5) {
407     param1.SetNonBendingCoor(Param[3]);
408     param1.SetBendingCoor(Param[4]);
409   }
410
411   AliMUONTrackHit *hit;
412   Bool_t goodDeterminant;
413   Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3;
414   Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10];
415   Double_t hitBendingCoor[10], hitNonBendingCoor[10];
416   Double_t hitBendingReso2[10], hitNonBendingReso2[10];
417   // dimension 10 in parameter ??? related to AliMUONConstants::NTrackingCh() !!!!
418   Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
419   TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
420   TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
421
422   // Predicted coordinates and  multiple scattering angles are first calculated
423   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
424     hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
425     zHit[hitNumber] = hit->GetHitForRecPtr()->GetZ();
426     // do something special if 2 hits with same Z ????
427     // security against infinite loop ????
428     (&param1)->ExtrapToZ(zHit[hitNumber]); // extrapolation
429     hit->SetTrackParam(&param1);
430     paramBendingCoor[hitNumber] = (&param1)->GetBendingCoor();
431     paramNonBendingCoor[hitNumber] = (&param1)->GetNonBendingCoor();
432     hitBendingCoor[hitNumber] = hit->GetHitForRecPtr()->GetBendingCoor();
433     hitNonBendingCoor[hitNumber] = hit->GetHitForRecPtr()->GetNonBendingCoor();
434     hitBendingReso2[hitNumber] = hit->GetHitForRecPtr()->GetBendingReso2();
435     hitNonBendingReso2[hitNumber] = hit->GetHitForRecPtr()->GetNonBendingReso2();
436     ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2  
437   }
438
439   // Calculates the covariance matrix
440   // One chamber is taken into account between successive hits.
441   // "ap" should be changed for taking into account the eventual missing hits
442   // by defining an "equivalent" chamber thickness !!!!
443   for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {    
444     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
445       // initialization to 0 (diagonal plus upper triangular part)
446       (*covBending)(hitNumber2, hitNumber1) = 0.0;
447       // contribution from multiple scattering in bending plane:
448       // loop over upstream hits
449       for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {     
450         (*covBending)(hitNumber2, hitNumber1) =
451           (*covBending)(hitNumber2, hitNumber1) +
452           ((zHit[hitNumber1] - zHit[hitNumber3]) *
453            (zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]); 
454       }
455       // equal contribution from multiple scattering in non bending plane
456       (*covNonBending)(hitNumber2, hitNumber1) =
457         (*covBending)(hitNumber2, hitNumber1);
458       if (hitNumber1 == hitNumber2) {
459         // Diagonal elements: add contribution from position measurements
460         // in bending plane
461         (*covBending)(hitNumber2, hitNumber1) =
462           (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
463         // and in non bending plane
464         (*covNonBending)(hitNumber2, hitNumber1) =
465           (*covNonBending)(hitNumber2, hitNumber1) + hitNonBendingReso2[hitNumber1];
466       }
467       else {
468         // Non diagonal elements: symmetrization
469         // for bending plane
470         (*covBending)(hitNumber1, hitNumber2) =
471           (*covBending)(hitNumber2, hitNumber1);
472         // and non bending plane
473         (*covNonBending)(hitNumber1, hitNumber2) =
474           (*covNonBending)(hitNumber2, hitNumber1);
475       }
476     } // for (hitNumber2 = hitNumber1;...
477   } // for (hitNumber1 = 0;...
478
479   // Inverts covariance matrix 
480   goodDeterminant = kTRUE;
481   // check whether the Invert method returns flag if matrix cannot be inverted,
482   // and do not calculate the Determinant in that case !!!!
483   if (covBending->Determinant() != 0) {
484     covBending->Invert();
485   } else {
486     goodDeterminant = kFALSE;
487     cout << "Warning in ChiMCS  Determinant Bending=0: " << endl;  
488   }
489   if (covNonBending->Determinant() != 0) {
490     covNonBending->Invert();
491   } else {
492     goodDeterminant = kFALSE;
493     cout << "Warning in ChiMCS  Determinant non Bending=0: " << endl;  
494   }
495
496   // It would be worth trying to calculate the inverse of the covariance matrix
497   // only once per fit, since it cannot change much in principle,
498   // and it would save a lot of computing time !!!!
499   
500   // Calculates Chi2
501   if (goodDeterminant) { // with Multiple Scattering if inversion correct
502     for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){ 
503       for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){
504         Chi2 = Chi2 +
505           ((*covBending)(hitNumber2, hitNumber1) * 
506            (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
507            (hitBendingCoor[hitNumber2] - paramBendingCoor[hitNumber2]));
508         Chi2 = Chi2 +
509           ((*covNonBending)(hitNumber2, hitNumber1) *
510            (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
511            (hitNonBendingCoor[hitNumber2] - paramNonBendingCoor[hitNumber2]));
512       }
513     }
514   } else {  // without Multiple Scattering if inversion impossible
515     for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++) { 
516       Chi2 = Chi2 +
517         ((hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
518          (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) /
519          hitBendingReso2[hitNumber1]);
520       Chi2 = Chi2 +
521         ((hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
522          (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) /
523          hitNonBendingReso2[hitNumber1]);      
524     }
525   }
526   
527   delete covBending;
528   delete covNonBending;
529 }
530
531 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
532 {
533   // Returns square of multiple Coulomb scattering angle
534   // at TrackHit pointed to by "TrackHit"
535   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
536   Double_t varMultipleScatteringAngle;
537   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
538   AliMUONTrackParam *param = TrackHit->GetTrackParam();
539   // Better implementation in AliMUONTrack class ????
540   slopeBending = param->GetBendingSlope();
541   slopeNonBending = param->GetNonBendingSlope();
542   // thickness in radiation length for the current track,
543   // taking local angle into account
544   radiationLength =
545     trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() *
546     TMath::Sqrt(1.0 +
547                 slopeBending * slopeBending + slopeNonBending * slopeNonBending);
548   inverseBendingMomentum2 = 
549     param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
550   inverseTotalMomentum2 =
551     inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
552     (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
553   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
554   // The velocity is assumed to be 1 !!!!
555   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
556     varMultipleScatteringAngle * varMultipleScatteringAngle;
557   return varMultipleScatteringAngle;
558 }