Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.cxx
index d1f77f1..2a3e830 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.3  2000/06/25 13:23:28  hristov
-stdlib.h needed for non-Linux compilation
+/* $Id$ */
 
-Revision 1.2  2000/06/15 07:58:48  morsch
-Code from MUON-dev joined
-
-Revision 1.1.2.3  2000/06/12 10:11:34  morsch
-Dummy copy constructor and assignment operator added
-
-Revision 1.1.2.2  2000/06/09 12:58:05  gosset
-Removed comment beginnings in Log sections of .cxx files
-Suppressed most violations of coding rules
-
-Revision 1.1.2.1  2000/06/07 14:44:53  gosset
-Addition of files for track reconstruction in C++
-*/
-
-//__________________________________________________________________________
+///////////////////////////////////////////////////
 //
-// Reconstructed track in ALICE dimuon spectrometer
-//__________________________________________________________________________
-
-#include "AliMUONTrack.h"
+// Reconstructed track
+// in
+// ALICE
+// dimuon
+// spectrometer
+//
+///////////////////////////////////////////////////
 
-#include <iostream.h>
+#include <Riostream.h> // for cout
+#include <stdlib.h> // for exit()
 
 #include <TClonesArray.h>
-#include <TMinuit.h>
 #include <TMath.h>
-#include <TMatrix.h>
+#include <TMatrixD.h>
+#include <TObjArray.h>
+#include <TVirtualFitter.h>
 
 #include "AliMUONEventReconstructor.h" 
 #include "AliMUONHitForRec.h" 
 #include "AliMUONSegment.h" 
+#include "AliMUONTrack.h"
 #include "AliMUONTrackHit.h"
 
-#include <stdlib.h>
-
-// variables to be known from minimization functions
-static AliMUONTrack *trackBeingFitted;
-static AliMUONTrackParam *trackParamBeingFitted;
-
 // Functions to be minimized with Minuit
 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
 
+void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
+
 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
 
 ClassImp(AliMUONTrack) // Class implementation in ROOT context
 
+TVirtualFitter* AliMUONTrack::fgFitter = NULL; 
+
   //__________________________________________________________________________
 AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor)
 {
   // Constructor from two Segment's
   fEventReconstructor = EventReconstructor; // link back to EventReconstructor
-  // memory allocation for the TClonesArray of reconstructed TrackHit's
-  fTrackHitsPtr = new  TClonesArray("AliMUONTrackHit", 10);
+  // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
+  fTrackHitsPtr = new TObjArray(10);
   fNTrackHits = 0;
   AddSegment(BegSegment); // add hits from BegSegment
   AddSegment(EndSegment); // add hits from EndSegment
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
+  // set fit conditions...
   fFitMCS = 0;
+  fFitNParam = 3;
+  fFitStart = 1;
+  fFitFMin = -1.0;
   return;
 }
 
@@ -86,54 +77,145 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec,
 {
   // Constructor from one Segment and one HitForRec
   fEventReconstructor = EventReconstructor; // link back to EventReconstructor
-  // memory allocation for the TClonesArray of reconstructed TrackHit's
-  fTrackHitsPtr = new  TClonesArray("AliMUONTrackHit", 10);
+  // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
+  fTrackHitsPtr = new TObjArray(10);
   fNTrackHits = 0;
   AddSegment(Segment); // add hits from Segment
   AddHitForRec(HitForRec); // add HitForRec
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
+  // set fit conditions...
   fFitMCS = 0;
+  fFitNParam = 3;
+  fFitStart = 1;
+  fFitFMin = -1.0;
   return;
 }
 
-AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
+  //__________________________________________________________________________
+AliMUONTrack::~AliMUONTrack()
+{
+  // Destructor
+  if (fTrackHitsPtr) {
+    delete fTrackHitsPtr; // delete the TObjArray of pointers to TrackHit's
+    fTrackHitsPtr = NULL;
+  }
+}
+
+  //__________________________________________________________________________
+AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack):TObject(MUONTrack)
 {
-// Dummy copy constructor
+  fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor);
+  fTrackParamAtVertex = MUONTrack.fTrackParamAtVertex;
+  fTrackHitsPtr =  new TObjArray(*MUONTrack.fTrackHitsPtr);
+  fNTrackHits =  MUONTrack.fNTrackHits;
+  fFitMCS     =  MUONTrack.fFitMCS;
+  fFitNParam  =  MUONTrack.fFitNParam;
+  fFitFMin    =  MUONTrack.fFitFMin;
+  fFitStart   =  MUONTrack.fFitStart;
 }
 
+  //__________________________________________________________________________
 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
 {
-// Dummy assignment operator
+  if (this == &MUONTrack)
     return *this;
+
+  fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor);
+  fTrackParamAtVertex = MUONTrack.fTrackParamAtVertex;
+  fTrackHitsPtr =  new TObjArray(*MUONTrack.fTrackHitsPtr);
+  fNTrackHits =  MUONTrack.fNTrackHits;
+  fFitMCS     =  MUONTrack.fFitMCS;
+  fFitNParam  =  MUONTrack.fFitNParam;
+  fFitFMin    =  MUONTrack.fFitFMin;
+  fFitStart   =  MUONTrack.fFitStart;
+  return *this;
+}
+
+  //__________________________________________________________________________
+void AliMUONTrack::Remove()
+{
+  // Remove current track from array of tracks,
+  // and corresponding track hits from array of track hits.
+  // Compress the TClonesArray it belongs to.
+  AliMUONTrackHit *nextTrackHit;
+  AliMUONEventReconstructor *eventRec = this->fEventReconstructor;
+  TClonesArray *trackHitsPtr = eventRec->GetRecTrackHitsPtr();
+  // Loop over all track hits of track
+  AliMUONTrackHit *trackHit = (AliMUONTrackHit*) fTrackHitsPtr->First();
+  while (trackHit) {
+    nextTrackHit = (AliMUONTrackHit*) fTrackHitsPtr->After(trackHit);
+    // Remove TrackHit from event TClonesArray.
+    // Destructor is called,
+    // hence links between HitForRec's and TrackHit's are updated
+    trackHitsPtr->Remove(trackHit);
+    trackHit = nextTrackHit;
+  }
+  // Remove the track from event TClonesArray
+  // Destructor is called,
+  // hence space for TObjArray of pointers to TrackHit's is freed
+  eventRec->GetRecTracksPtr()->Remove(this);
+  // Number of tracks decreased by 1
+  eventRec->SetNRecTracks(eventRec->GetNRecTracks() - 1);
+  // Compress event TClonesArray of Track's:
+  // this is essential to retrieve the TClonesArray afterwards
+  eventRec->GetRecTracksPtr()->Compress();
+  // Compress event TClonesArray of TrackHit's:
+  // this is probably also essential to retrieve the TClonesArray afterwards
+  trackHitsPtr->Compress();
 }
 
   //__________________________________________________________________________
 void AliMUONTrack::SetFitMCS(Int_t FitMCS)
 {
-  // Set track fit option with or without multiple Coulomb scattering
+  // Set multiple Coulomb scattering option for track fit "fFitMCS"
   // from "FitMCS" argument: 0 without, 1 with
-  fFitMCS = FitMCS;
+  if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS;
   // better implementation with enum(with, without) ????
+  else {
+    cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl;
+    cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
+    exit(0);
+  }
+  return;
+}
+
+  //__________________________________________________________________________
+void AliMUONTrack::SetFitNParam(Int_t FitNParam)
+{
+  // Set number of parameters for track fit "fFitNParam" from "FitNParam":
+  // 3 for momentum, 5 for momentum and position
+  if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam;
+  else {
+    cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl;
+    cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl;
+    exit(0);
+  }
+  return;
+}
+
+  //__________________________________________________________________________
+void AliMUONTrack::SetFitStart(Int_t FitStart)
+{
+  // Set multiple Coulomb scattering option for track fit "fFitStart"
+  // from "FitStart" argument: 0 without, 1 with
+  if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart;
+  // better implementation with enum(vertex, firstHit) ????
+  else {
+    cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl;
+    cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
+    exit(0);
+  }
   return;
 }
 
-// Inline functions for Get and Set: inline removed because it does not work !!!!
-AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) {
-  // Get pointer to fTrackParamAtVertex
-  return &fTrackParamAtVertex;}
-AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
+  //__________________________________________________________________________
+AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) const {
   // Get pointer to TrackParamAtFirstHit
   return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
-TClonesArray* AliMUONTrack::GetTrackHitsPtr(void) {
-  // Get fTrackHitsPtr
-  return fTrackHitsPtr;}
-Int_t AliMUONTrack::GetNTrackHits(void) {
-  // Get fNTrackHits
-  return fNTrackHits;}
 
   //__________________________________________________________________________
-void AliMUONTrack::RecursiveDump(void)
+void AliMUONTrack::RecursiveDump(void) const
 {
   // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's
   AliMUONTrackHit *trackHit;
@@ -155,81 +237,112 @@ void AliMUONTrack::RecursiveDump(void)
 }
 
   //__________________________________________________________________________
-void AliMUONTrack::Fit(AliMUONTrackParam *TrackParam, Int_t NParam)
+Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track)
+{
+  // Returns the number of hits in common
+  // between the current track ("this")
+  // and the track pointed to by "Track".
+  Int_t hitsInCommon = 0;
+  AliMUONTrackHit *trackHit1, *trackHit2;
+  // Loop over hits of first track
+  trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->First();
+  while (trackHit1) {
+    // Loop over hits of second track
+    trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->First();
+    while (trackHit2) {
+      // Increment "hitsInCommon" if both TrackHits point to the same HitForRec
+      if ( (trackHit1->GetHitForRecPtr()) ==
+          (trackHit2->GetHitForRecPtr())    ) hitsInCommon++;
+      trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->After(trackHit2);
+    } // trackHit2
+    trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->After(trackHit1);
+  } // trackHit1
+  return hitsInCommon;
+}
+
+  //__________________________________________________________________________
+void AliMUONTrack::Fit()
 {
   // Fit the current track ("this"),
-  // starting with track parameters pointed to by "TrackParam",
-  // and with 3 or 5 parameters ("NParam"):
-  // 3 if one keeps X and Y fixed in "TrackParam",
-  // 5 if one lets them vary.
-  if ((NParam != 3) && (NParam != 5)) {
-    cout << "ERROR in AliMUONTrack::Fit, NParam = " << NParam;
-    cout << " , i.e. neither 3 nor 5 ====> EXIT" << endl;
-    exit(0); // right instruction for exit ????
-  }
-  Int_t error = 0;
+  // with or without multiple Coulomb scattering according to "fFitMCS",
+  // with the number of parameters given by "fFitNParam"
+  // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary),
+  // starting, according to "fFitStart",
+  // with track parameters at vertex or at the first TrackHit.
+  // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before
+  // by calling the corresponding Set methods.
   Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
-  TString parName;
-  TMinuit *minuit = new TMinuit(5);
-  trackBeingFitted = this; // for the track to be known from the function to minimize
-  trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ????
-  // try to use TVirtualFitter to get this feature automatically !!!!
-  minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ????
+  char parName[50];
+  AliMUONTrackParam *trackParam;
+  // Check if Minuit is initialized...
+  fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ???
+  fgFitter->Clear(); // necessary ???? probably yes
+  // how to reset the printout number at every fit ????
+  // is there any risk to leave it like that ????
   // how to go faster ???? choice of Minuit parameters like EDM ????
   // choice of function to be minimized according to fFitMCS
-  if (fFitMCS == 0) minuit->SetFCN(TrackChi2);
-  else minuit->SetFCN(TrackChi2MCS);
-  minuit->SetPrintLevel(1); // More printing !!!!
-  // set first 3 parameters
+  if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
+  else fgFitter->SetFCN(TrackChi2MCS);
+  arg[0] = -1;
+  fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
+  // Parameters according to "fFitStart"
+  // (should be a function to be used at every place where needed ????)
+  if (fFitStart == 0) trackParam = &fTrackParamAtVertex;
+  else trackParam = this->GetTrackParamAtFirstHit();
+  // set first 3 Minuit parameters
   // could be tried with no limits for the search (min=max=0) ????
-  minuit->mnparm(0, "InvBenP",
-                TrackParam->GetInverseBendingMomentum(),
-                0.003, -0.4, 0.4, error);
-  minuit->mnparm(1, "BenS",
-                TrackParam->GetBendingSlope(),
-                0.001, -0.5, 0.5, error);
-  minuit->mnparm(2, "NonBenS",
-                TrackParam->GetNonBendingSlope(),
-                0.001, -0.5, 0.5, error);
-  if (NParam == 5) {
-    // set last 2 parameters (no limits for the search: min=max=0)
-    minuit->mnparm(3, "X",
-                  TrackParam->GetNonBendingCoor(),
-                  0.03, 0.0, 0.0, error);
-    minuit->mnparm(4, "Y",
-                  TrackParam->GetBendingCoor(),
-                  0.10, 0.0, 0.0, error);
+  fgFitter->SetParameter(0, "InvBenP",
+                        trackParam->GetInverseBendingMomentum(),
+                        0.003, -0.4, 0.4);
+  fgFitter->SetParameter(1, "BenS",
+                        trackParam->GetBendingSlope(),
+                        0.001, -0.5, 0.5);
+  fgFitter->SetParameter(2, "NonBenS",
+                        trackParam->GetNonBendingSlope(),
+                        0.001, -0.5, 0.5);
+  if (fFitNParam == 5) {
+    // set last 2 Minuit parameters
+    // mandatory limits in Bending to avoid NaN values of parameters
+    fgFitter->SetParameter(3, "X",
+                          trackParam->GetNonBendingCoor(),
+                          0.03, -500.0, 500.0);
+    // mandatory limits in non Bending to avoid NaN values of parameters
+    fgFitter->SetParameter(4, "Y",
+                          trackParam->GetBendingCoor(),
+                          0.10, -500.0, 500.0);
   }
   // search without gradient calculation in the function
-  minuit->mnexcm("SET NOGRADIENT", arg, 0, error);
+  fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
   // minimization
-  minuit->mnexcm("MINIMIZE", arg, 0, error);
+  fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
   // exit from Minuit
-  minuit->mnexcm("EXIT", arg, 0, error); // necessary ????
-  // print results
-  minuit->mnpout(0, parName, invBenP, errorParam, lower, upper, error);
-  minuit->mnpout(1, parName, benC, errorParam, lower, upper, error);
-  minuit->mnpout(2, parName, nonBenC, errorParam, lower, upper, error);
-  if (NParam == 5) {
-    minuit->mnpout(3, parName, x, errorParam, lower, upper, error);
-    minuit->mnpout(4, parName, y, errorParam, lower, upper, error);
+  fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
+  // get results into "invBenP", "benC", "nonBenC" ("x", "y")
+  fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
+  fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
+  fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
+  if (fFitNParam == 5) {
+    fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
+    fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
   }
   // result of the fit into track parameters
-  TrackParam->SetInverseBendingMomentum(invBenP);
-  TrackParam->SetBendingSlope(benC);
-  TrackParam->SetNonBendingSlope(nonBenC);
-  if (NParam == 5) {
-    TrackParam->SetNonBendingCoor(x);
-    TrackParam->SetBendingCoor(y);
+  trackParam->SetInverseBendingMomentum(invBenP);
+  trackParam->SetBendingSlope(benC);
+  trackParam->SetNonBendingSlope(nonBenC);
+  if (fFitNParam == 5) {
+    trackParam->SetNonBendingCoor(x);
+    trackParam->SetBendingCoor(y);
   }
-  trackBeingFitted = NULL;
-  delete minuit;
+  // global result of the fit
+  Double_t fedm, errdef;
+  Int_t npari, nparx;
+  fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx);
 }
 
   //__________________________________________________________________________
 void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
 {
-  // Add Segment
+  // Add Segment to the track
   AddHitForRec(Segment->GetHitForRec1()); // 1st hit
   AddHitForRec(Segment->GetHitForRec2()); // 2nd hit
 }
@@ -237,17 +350,27 @@ void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
   //__________________________________________________________________________
 void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
 {
-  // Add HitForRec
-  new ((*fTrackHitsPtr)[fNTrackHits]) AliMUONTrackHit(HitForRec);
+  // Add HitForRec to the track:
+  // actual TrackHit into TClonesArray of TrackHit's for the event;
+  // pointer to actual TrackHit in TObjArray of pointers to TrackHit's for the track
+  TClonesArray *recTrackHitsPtr = this->fEventReconstructor->GetRecTrackHitsPtr();
+  Int_t eventTrackHits = this->fEventReconstructor->GetNRecTrackHits();
+  // event
+  AliMUONTrackHit* trackHit =
+    new ((*recTrackHitsPtr)[eventTrackHits]) AliMUONTrackHit(HitForRec);
+  this->fEventReconstructor->SetNRecTrackHits(eventTrackHits + 1);
+  // track
+  fTrackHitsPtr->Add(trackHit);
   fNTrackHits++;
 }
 
   //__________________________________________________________________________
-void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam)
+void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) const
 {
   // Set track parameters at TrackHit with index "indexHit"
   // from the track parameters pointed to by "TrackParam".
-  AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]);
+  //PH  AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]);
+  AliMUONTrackHit* trackHit = (AliMUONTrackHit*) (fTrackHitsPtr->At(indexHit));
   trackHit->SetTrackParam(TrackParam);
 }
 
@@ -295,7 +418,7 @@ void AliMUONTrack::SetTrackParamAtVertex()
 }
 
   //__________________________________________________________________________
-void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
+void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
 {
   // Return the "Chi2" to be minimized with Minuit for track fitting,
   // with "NParam" parameters
@@ -303,13 +426,17 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para
   // Assumes that the track hits are sorted according to increasing Z.
   // Track parameters at each TrackHit are updated accordingly.
   // Multiple Coulomb scattering is not taken into account
+  AliMUONTrack *trackBeingFitted;
   AliMUONTrackHit* hit;
   AliMUONTrackParam param1;
   Int_t hitNumber;
   Double_t zHit;
   Chi2 = 0.0; // initialize Chi2
   // copy of track parameters to be fitted
-  param1 = *trackParamBeingFitted;
+  trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
+  if (trackBeingFitted->GetFitStart() == 0)
+    param1 = *(trackBeingFitted->GetTrackParamAtVertex());
+  else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
   // Minuit parameters to be fitted into this copy
   param1.SetInverseBendingMomentum(Param[0]);
   param1.SetBendingSlope(Param[1]);
@@ -342,7 +469,7 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para
 }
 
   //__________________________________________________________________________
-void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
+void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
 {
   // Return the "Chi2" to be minimized with Minuit for track fitting,
   // with "NParam" parameters
@@ -350,10 +477,14 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
   // Assumes that the track hits are sorted according to increasing Z.
   // Track parameters at each TrackHit are updated accordingly.
   // Multiple Coulomb scattering is taken into account with covariance matrix.
+  AliMUONTrack *trackBeingFitted;
   AliMUONTrackParam param1;
   Chi2 = 0.0; // initialize Chi2
   // copy of track parameters to be fitted
-  param1 = *trackParamBeingFitted;
+  trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
+  if (trackBeingFitted->GetFitStart() == 0)
+    param1 = *(trackBeingFitted->GetTrackParamAtVertex());
+  else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
   // Minuit parameters to be fitted into this copy
   param1.SetInverseBendingMomentum(Param[0]);
   param1.SetBendingSlope(Param[1]);
@@ -363,104 +494,139 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
     param1.SetBendingCoor(Param[4]);
   }
 
-  AliMUONTrackHit* hit, hit1, hit2, hit3;
-  Bool_t GoodDeterminant;
-  Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3;
-  Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10];
-  Double_t hitBendingCoor[10], hitNonBendingCoor[10];
-  Double_t hitBendingReso2[10], hitNonBendingReso2[10];
-  Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
-  TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
-  TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
+  AliMUONTrackHit *hit;
+  Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
+  Double_t z, z1, z2, z3;
+  AliMUONTrackHit *hit1, *hit2, *hit3;
+  Double_t hbc1, hbc2, pbc1, pbc2;
+  Double_t hnbc1, hnbc2, pnbc1, pnbc2;
+  Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
+  TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
+  TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
+  Double_t *msa2 = new Double_t[numberOfHit];
 
   // Predicted coordinates and  multiple scattering angles are first calculated
   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
     hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
-    zHit[hitNumber] = hit->GetHitForRecPtr()->GetZ();
+    z = hit->GetHitForRecPtr()->GetZ();
     // do something special if 2 hits with same Z ????
     // security against infinite loop ????
-    (&param1)->ExtrapToZ(zHit[hitNumber]); // extrapolation
+    (&param1)->ExtrapToZ(z); // extrapolation
     hit->SetTrackParam(&param1);
-    paramBendingCoor[hitNumber]= (&param1)->GetBendingCoor();
-    paramNonBendingCoor[hitNumber]= (&param1)->GetNonBendingCoor();
-    hitBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetBendingCoor();
-    hitNonBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingCoor();
-    hitBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetBendingReso2();
-    hitNonBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingReso2();
-    ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2  
+    // square of multiple scattering angle at current hit, with one chamber
+    msa2[hitNumber] = MultipleScatteringAngle2(hit);
+    // correction for eventual missing hits or multiple hits in a chamber,
+    // according to the number of chambers
+    // between the current hit and the previous one
+    chCurrent = hit->GetHitForRecPtr()->GetChamberNumber();
+    if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
+    chPrev = chCurrent;
   }
 
   // Calculates the covariance matrix
-  for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {    
+  for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { 
+    hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
+    z1 = hit1->GetHitForRecPtr()->GetZ();
     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
-      (*covBending)(hitNumber1, hitNumber2) = 0;
-      (*covBending)(hitNumber2, hitNumber1) = 0;
-      if (hitNumber1 == hitNumber2){ // diagonal elements
-       (*covBending)(hitNumber2, hitNumber1) =
-         (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
-      }
-      // Multiple Scattering...  loop on upstream chambers ??
-      for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++){     
+      hit2 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
+      z2 = hit2->GetHitForRecPtr()->GetZ();
+      // initialization to 0 (diagonal plus upper triangular part)
+      (*covBending)(hitNumber2, hitNumber1) = 0.0;
+      // contribution from multiple scattering in bending plane:
+      // loop over upstream hits
+      for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {    
+       hit3 = (AliMUONTrackHit*)
+         (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber3];
+       z3 = hit3->GetHitForRecPtr()->GetZ();
        (*covBending)(hitNumber2, hitNumber1) =
          (*covBending)(hitNumber2, hitNumber1) +
-         ((zHit[hitNumber1] - zHit[hitNumber3]) *
-          (zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]); 
-      }  
-      (*covNonBending)(hitNumber1, hitNumber2) = 0;
+         ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); 
+      }
+      // equal contribution from multiple scattering in non bending plane
       (*covNonBending)(hitNumber2, hitNumber1) =
        (*covBending)(hitNumber2, hitNumber1);
-      if (hitNumber1 == hitNumber2) {  // diagonal elements
+      if (hitNumber1 == hitNumber2) {
+       // Diagonal elements: add contribution from position measurements
+       // in bending plane
+       (*covBending)(hitNumber2, hitNumber1) =
+         (*covBending)(hitNumber2, hitNumber1) +
+         hit1->GetHitForRecPtr()->GetBendingReso2();
+       // and in non bending plane
        (*covNonBending)(hitNumber2, hitNumber1) =
-         (*covNonBending)(hitNumber2, hitNumber1) -
-         hitBendingReso2[hitNumber1] + hitNonBendingReso2[hitNumber1] ;
-      }      
-    }
-  }
-
-  // Inverts covariance matrix 
-  GoodDeterminant = kTRUE;
-  if (covBending->Determinant() != 0) {
-    covBending->Invert();
-  } else {
-    GoodDeterminant = kFALSE;
-    cout << "Warning in ChiMCS  Determinant Bending=0: " << endl;  
-  }
-  if (covNonBending->Determinant() != 0){
-    covNonBending->Invert();
-  } else {
-    GoodDeterminant = kFALSE;
-    cout << "Warning in ChiMCS  Determinant non Bending=0: " << endl;  
-  }
+         (*covNonBending)(hitNumber2, hitNumber1) +
+         hit1->GetHitForRecPtr()->GetNonBendingReso2();
+      }
+      else {
+       // Non diagonal elements: symmetrization
+       // for bending plane
+       (*covBending)(hitNumber1, hitNumber2) =
+         (*covBending)(hitNumber2, hitNumber1);
+       // and non bending plane
+       (*covNonBending)(hitNumber1, hitNumber2) =
+         (*covNonBending)(hitNumber2, hitNumber1);
+      }
+    } // for (hitNumber2 = hitNumber1;...
+  } // for (hitNumber1 = 0;...
+    
+  // Inversion of covariance matrices
+  // with "mnvertLocal", local "mnvert" function of Minuit.
+  // One cannot use directly "mnvert" since "TVirtualFitter" does not know it.
+  // One will have to replace this local function by the right inversion function
+  // from a specialized Root package for symmetric positive definite matrices,
+  // when available!!!!
+  Int_t ifailBending;
+  mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit,
+             ifailBending);
+  Int_t ifailNonBending;
+  mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit,
+             ifailNonBending);
+
+  // It would be worth trying to calculate the inverse of the covariance matrix
+  // only once per fit, since it cannot change much in principle,
+  // and it would save a lot of computing time !!!!
   
   // Calculates Chi2
-  if (GoodDeterminant) { // with Multiple Scattering if inversion correct
-    for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){ 
-      for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){
-       Chi2 = Chi2 +
-         ((*covBending)(hitNumber2, hitNumber1) * 
-          (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
-          (hitBendingCoor[hitNumber2] - paramBendingCoor[hitNumber2]));
+  if ((ifailBending == 0) && (ifailNonBending == 0)) {
+    // with Multiple Scattering if inversion correct
+    for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
+      hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
+      hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
+      pbc1 = hit1->GetTrackParam()->GetBendingCoor();
+      hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
+      pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
+      for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
+       hit2 = (AliMUONTrackHit*)
+         (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
+       hbc2 = hit2->GetHitForRecPtr()->GetBendingCoor();
+       pbc2 = hit2->GetTrackParam()->GetBendingCoor();
+       hnbc2 = hit2->GetHitForRecPtr()->GetNonBendingCoor();
+       pnbc2 = hit2->GetTrackParam()->GetNonBendingCoor();
        Chi2 = Chi2 +
+         ((*covBending)(hitNumber2, hitNumber1) *
+          (hbc1 - pbc1) * (hbc2 - pbc2)) +
          ((*covNonBending)(hitNumber2, hitNumber1) *
-          (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
-          (hitNonBendingCoor[hitNumber2] - paramNonBendingCoor[hitNumber2]));
+          (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
       }
     }
-  } else {  // without Multiple Scattering if inversion impossible
-    for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++) { 
-      Chi2 = Chi2 +
-       ((hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
-        (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) /
-        hitBendingReso2[hitNumber1]);
-      Chi2 = Chi2 +
-       ((hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
-        (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) /
-        hitNonBendingReso2[hitNumber1]);      
+  } else {
+    // without Multiple Scattering if inversion impossible
+    for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
+      hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
+      hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
+      pbc1 = hit1->GetTrackParam()->GetBendingCoor();
+      hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
+      pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
+      Chi2 = Chi2 + 
+       ((hbc1 - pbc1) * (hbc1 - pbc1) /
+        hit1->GetHitForRecPtr()->GetBendingReso2()) +
+       ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) /
+        hit1->GetHitForRecPtr()->GetNonBendingReso2());
     }
   }
   
   delete covBending;
   delete covNonBending;
+  delete [] msa2;
 }
 
 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
@@ -469,7 +635,9 @@ Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
   // at TrackHit pointed to by "TrackHit"
   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
   Double_t varMultipleScatteringAngle;
+  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
   AliMUONTrackParam *param = TrackHit->GetTrackParam();
+  // Better implementation in AliMUONTrack class ????
   slopeBending = param->GetBendingSlope();
   slopeNonBending = param->GetNonBendingSlope();
   // thickness in radiation length for the current track,
@@ -481,10 +649,111 @@ Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
   inverseBendingMomentum2 = 
     param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
   inverseTotalMomentum2 =
-    inverseBendingMomentum2 * (1.0 + slopeBending*slopeBending) /
+    inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
     (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
+  // The velocity is assumed to be 1 !!!!
   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
     varMultipleScatteringAngle * varMultipleScatteringAngle;
   return varMultipleScatteringAngle;
 }
+
+//______________________________________________________________________________
+ void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
+{
+//*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
+//*-*                    ==========================
+//*-*        inverts a symmetric matrix.   matrix is first scaled to
+//*-*        have all ones on the diagonal (equivalent to change of units)
+//*-*        but no pivoting is done since matrix is positive-definite.
+//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
+
+  // taken from TMinuit package of Root (l>=n)
+  // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
+  //  Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
+  Double_t * localVERTs = new Double_t[n];
+  Double_t * localVERTq = new Double_t[n];
+  Double_t * localVERTpp = new Double_t[n];
+  // fMaxint changed to localMaxint
+  Int_t localMaxint = n;
+
+    /* System generated locals */
+    Int_t aOffset;
+
+    /* Local variables */
+    Double_t si;
+    Int_t i, j, k, kp1, km1;
+
+    /* Parameter adjustments */
+    aOffset = l + 1;
+    a -= aOffset;
+
+    /* Function Body */
+    ifail = 0;
+    if (n < 1) goto L100;
+    if (n > localMaxint) goto L100;
+//*-*-                  scale matrix by sqrt of diag elements
+    for (i = 1; i <= n; ++i) {
+        si = a[i + i*l];
+        if (si <= 0) goto L100;
+        localVERTs[i-1] = 1 / TMath::Sqrt(si);
+    }
+    for (i = 1; i <= n; ++i) {
+        for (j = 1; j <= n; ++j) {
+            a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
+        }
+    }
+//*-*-                                       . . . start main loop . . . .
+    for (i = 1; i <= n; ++i) {
+        k = i;
+//*-*-                  preparation for elimination step1
+        if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
+        else goto L100;
+        localVERTpp[k-1] = 1;
+        a[k + k*l] = 0;
+        kp1 = k + 1;
+        km1 = k - 1;
+        if (km1 < 0) goto L100;
+        else if (km1 == 0) goto L50;
+        else               goto L40;
+L40:
+        for (j = 1; j <= km1; ++j) {
+            localVERTpp[j-1] = a[j + k*l];
+            localVERTq[j-1]  = a[j + k*l]*localVERTq[k-1];
+            a[j + k*l]   = 0;
+        }
+L50:
+        if (k - n < 0) goto L51;
+        else if (k - n == 0) goto L60;
+        else                goto L100;
+L51:
+        for (j = kp1; j <= n; ++j) {
+            localVERTpp[j-1] = a[k + j*l];
+            localVERTq[j-1]  = -a[k + j*l]*localVERTq[k-1];
+            a[k + j*l]   = 0;
+        }
+//*-*-                  elimination proper
+L60:
+        for (j = 1; j <= n; ++j) {
+            for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
+        }
+    }
+//*-*-                  elements of left diagonal and unscaling
+    for (j = 1; j <= n; ++j) {
+        for (k = 1; k <= j; ++k) {
+            a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
+            a[j + k*l] = a[k + j*l];
+        }
+    }
+    delete localVERTs;
+    delete localVERTq;
+    delete localVERTpp;
+    return;
+//*-*-                  failure return
+L100:
+    delete localVERTs;
+    delete localVERTq;
+    delete localVERTpp;
+    ifail = 1;
+} /* mnvertLocal */
+