Arrays with variable size created by new and deleted at the end
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.cxx
index 2f0474cc3ce0b2f47c4e879f394db32f9ddfbe4c..919c11fd5f8c43f407f3d58820985962117e5ab5 100644 (file)
 
 /*
 $Log$
+Revision 1.10  2001/04/09 12:25:09  gosset
+Inversion of covariance matrices with local copy of TMinuit::mnvert,
+for symmetric positive definite matrices, instead of TMatrixD::Invert
+
+Revision 1.9  2001/01/17 20:59:24  hristov
+chPrev initialised
+
+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
@@ -60,7 +95,7 @@ Addition of files for track reconstruction in C++
 
 #include <TClonesArray.h>
 #include <TMath.h>
-#include <TMatrix.h>
+#include <TMatrixD.h>
 #include <TObjArray.h>
 #include <TVirtualFitter.h>
 
@@ -75,6 +110,8 @@ Addition of files for track reconstruction in C++
 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
@@ -295,7 +332,7 @@ void AliMUONTrack::Fit()
   // choice of function to be minimized according to fFitMCS
   if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
   else fgFitter->SetFCN(TrackChi2MCS);
-  arg[0] = 1;
+  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 ????)
@@ -313,13 +350,15 @@ void AliMUONTrack::Fit()
                         trackParam->GetNonBendingSlope(),
                         0.001, -0.5, 0.5);
   if (fFitNParam == 5) {
-    // set last 2 Minuit parameters (no limits for the search: min=max=0)
+    // set last 2 Minuit parameters
+    // mandatory limits in Bending to avoid NaN values of parameters
     fgFitter->SetParameter(3, "X",
                           trackParam->GetNonBendingCoor(),
-                          0.03, 0.0, 0.0);
+                          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, 0.0, 0.0);
+                          0.10, -500.0, 500.0);
   }
   // search without gradient calculation in the function
   fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
@@ -504,48 +543,52 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
   }
 
   AliMUONTrackHit *hit;
-  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];
-  // dimension 10 in parameter ??? related to AliMUONConstants::NTrackingCh() !!!!
-  Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
-  TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
-  TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
+  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
-  // One chamber is taken into account between successive hits.
-  // "ap" should be changed for taking into account the eventual missing hits
-  // by defining an "equivalent" chamber thickness !!!!
-  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++) {
+      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]); 
+         ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); 
       }
       // equal contribution from multiple scattering in non bending plane
       (*covNonBending)(hitNumber2, hitNumber1) =
@@ -554,10 +597,12 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
        // Diagonal elements: add contribution from position measurements
        // in bending plane
        (*covBending)(hitNumber2, hitNumber1) =
-         (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
+         (*covBending)(hitNumber2, hitNumber1) +
+         hit1->GetHitForRecPtr()->GetBendingReso2();
        // and in non bending plane
        (*covNonBending)(hitNumber2, hitNumber1) =
-         (*covNonBending)(hitNumber2, hitNumber1) + hitNonBendingReso2[hitNumber1];
+         (*covNonBending)(hitNumber2, hitNumber1) +
+         hit1->GetHitForRecPtr()->GetNonBendingReso2();
       }
       else {
        // Non diagonal elements: symmetrization
@@ -570,57 +615,66 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
       }
     } // for (hitNumber2 = hitNumber1;...
   } // for (hitNumber1 = 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) { // 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)
@@ -651,3 +705,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 a_offset;
+
+    /* Local variables */
+    Double_t si;
+    Int_t i, j, k, kp1, km1;
+
+    /* Parameter adjustments */
+    a_offset = l + 1;
+    a -= a_offset;
+
+    /* 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 */
+