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 }