Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.cxx
index a3d4211..2a3e830 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.8  2001/01/08 11:01:02  gosset
-Modifications used for addendum to Dimuon TDR (JP Cussonneau):
-*. MaxBendingMomentum to make both a segment and a track (default 500)
-*. MaxChi2 per degree of freedom to make a track (default 100)
-*. MinBendingMomentum used also to make a track
-   and not only a segment (default 3)
-*. wider roads for track search in stations 1 to 3
-*. extrapolation to actual Z instead of Z(chamber) in FollowTracks
-*. in track fit:
-   - limits on parameters X and Y (+/-500)
-   - covariance matrices in double precision
-   - normalization of covariance matrices before inversion
-   - suppression of Minuit printouts
-*. correction against memory leak (delete extrapHit) in FollowTracks
-*. RMax to 10 degrees with Z(chamber) instead of fixed values;
-   RMin and Rmax cuts suppressed in NewHitForRecFromGEANT,
-   because useless with realistic geometry
-
-Revision 1.7  2000/09/19 15:50:46  gosset
-TrackChi2MCS function: covariance matrix better calculated,
-taking into account missing planes...
-
-Revision 1.6  2000/07/20 12:45:27  gosset
-New "EventReconstructor..." structure,
-       hopefully more adapted to tree/streamer.
-"AliMUONEventReconstructor::RemoveDoubleTracks"
-       to keep only one track among similar ones.
-
-Revision 1.5  2000/07/18 16:04:06  gosset
-AliMUONEventReconstructor package:
-* a few minor modifications and more comments
-* a few corrections
-  * right sign for Z of raw clusters
-  * right loop over chambers inside station
-  * symmetrized covariance matrix for measurements (TrackChi2MCS)
-  * right sign of charge in extrapolation (ExtrapToZ)
-  * right zEndAbsorber for Branson correction below 3 degrees
-* use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit
-* no parameter for AliMUONTrack::Fit() but more fit parameters in Track object
-
-Revision 1.4  2000/06/30 10:15:48  gosset
-Changes to EventReconstructor...:
-precision fit with multiple Coulomb scattering;
-extrapolation to vertex with Branson correction in absorber (JPC)
-
-Revision 1.3  2000/06/25 13:23:28  hristov
-stdlib.h needed for non-Linux compilation
-
-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
-//__________________________________________________________________________
+/* $Id$ */
 
-#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 <TMath.h>
@@ -95,14 +37,15 @@ Addition of files for track reconstruction in C++
 #include "AliMUONEventReconstructor.h" 
 #include "AliMUONHitForRec.h" 
 #include "AliMUONSegment.h" 
+#include "AliMUONTrack.h"
 #include "AliMUONTrackHit.h"
 
-#include <stdlib.h>
-
 // 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
@@ -160,16 +103,33 @@ AliMUONTrack::~AliMUONTrack()
 }
 
   //__________________________________________________________________________
-AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
+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;
 }
 
   //__________________________________________________________________________
@@ -250,12 +210,12 @@ void AliMUONTrack::SetFitStart(Int_t FitStart)
 }
 
   //__________________________________________________________________________
-AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
+AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) const {
   // Get pointer to TrackParamAtFirstHit
   return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
 
   //__________________________________________________________________________
-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;
@@ -405,11 +365,12 @@ void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
 }
 
   //__________________________________________________________________________
-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);
 }
 
@@ -457,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
@@ -508,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
@@ -534,7 +495,6 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
   }
 
   AliMUONTrackHit *hit;
-  Bool_t goodDeterminant;
   Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
   Double_t z, z1, z2, z3;
   AliMUONTrackHit *hit1, *hit2, *hit3;
@@ -607,50 +567,27 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
       }
     } // for (hitNumber2 = hitNumber1;...
   } // for (hitNumber1 = 0;...
-  // Normalization of covariance matrices
-  Double_t normCovBending2 = covBending->E2Norm();
-  Double_t normCovNonBending2 = covNonBending->E2Norm();
-  (*covBending) *= 1/normCovBending2;
-  (*covNonBending) *= 1/normCovNonBending2;
-//   if (covBending->Determinant() < 1.e-33) {
-//     printf(" *** covBending *** \n");
-//     covBending->Print();
-//     printf(" *** covNonBending *** \n");
-//     covNonBending->Print();
-//     cout << " number of hits " <<  numberOfHit << endl;
-//     cout << "Momentum = " << 1/Param[0] <<endl;
-//     cout << "normCovBending = " << normCovBending2 << endl; 
-//     cout << "normCovNonBending = " << normCovNonBending2 << endl; 
-//     exit(0);
     
-//   }
-  // Inverts covariance matrix 
-  goodDeterminant = kTRUE;
-  // check whether the Invert method returns flag if matrix cannot be inverted,
-  // and do not calculate the Determinant in that case !!!!
-  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;  
-  }
+  // 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) {
+  if ((ifailBending == 0) && (ifailNonBending == 0)) {
     // with Multiple Scattering if inversion correct
-    // Inverse  matrices without normalization
-    (*covBending) *= 1/normCovBending2;
-    (*covNonBending) *= 1/normCovNonBending2;
     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
       hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
       hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
@@ -720,3 +657,103 @@ Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
     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 */
+