STEER/AliMillepede.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 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------------
  19 /// \class AliMillepede
  20 /// Detector independent alignment class
  21 ///
  22 /// This modified C++ version is based on a C++ translation of Millepede used
  23 /// for LHCb Vertex Detector alignment (lhcb-2005-101), available here
  24 /// http://isscvs.cern.ch/cgi-bin/cvsweb.cgi/Alignment/AlignmentTools/src/?cvsroot=lhcb
  25 /// The original millepede fortran package is available at:
  26 /// http://www.desy.de/~blobel/wwwmille.html
  27 ///
  28 /// \author Javier Castillo
  29 //-----------------------------------------------------------------------------
  30
  31 #include <TArrayI.h>
  32 #include <TArrayD.h>
  33 #include <TMath.h>
  34
  35 #include "AliLog.h"
  36
  37 #include "AliMillepede.h"
  38
  39 //=============================================================================
  40 AliMillepede::AliMillepede()
  41   : TObject(),
  42     fIndexLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
  43     fDerivLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
  44     fIndexAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
  45     fDerivAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
  46     fLocEqPlace(1000),
  47     fNIndexLocEq(0),
  48     fNDerivLocEq(0),
  49     fNIndexAllEqs(0),
  50     fNDerivAllEqs(0),
  51     fNLocEqPlace(0),
  52     fNLocalEquations(0),
  53     fResCutInit(0.),
  54     fResCut(0.),
  55     fChi2CutFactor(1.0),
  56     fChi2CutRef(1.0),
  57     fIter(0),
  58     fMaxIter(10),
  59     fNStdDev(3),
  60     fNGlobalConstraints(0),
  61     fNLocalFits(0),
  62     fNLocalFitsRejected(0),
  63     fNGlobalPar(0),
  64     fNLocalPar(0)
  65 {
  66   /// Standard constructor
  67
  68   AliInfo("                                           ");
  69   AliInfo("            * o o                   o      ");
  70   AliInfo("              o o                   o      ");
  71   AliInfo("   o ooooo  o o o  oo  ooo   oo   ooo  oo  ");
  72   AliInfo("    o  o  o o o o o  o o  o o  o o  o o  o ");
  73   AliInfo("    o  o  o o o o oooo o  o oooo o  o oooo ");
  74   AliInfo("    o  o  o o o o o    ooo  o    o  o o    ");
  75   AliInfo("    o  o  o o o o  oo  o     oo   ooo  oo  ++ starts");
  76   AliInfo("                                           ");
  77
  78 }
  79
  80 //=============================================================================
  81 AliMillepede::~AliMillepede() {
  82   /// Destructor
  83
  84 }
  85
  86 //=============================================================================
  87 Int_t AliMillepede::InitMille(int nGlo, int nLoc, int lNStdDev,
  88                               double lResCut, double lResCutInit)
  89 {
  90   /// Initialization of millepede
  91   AliDebug(1,"");
  92   AliDebug(1,"----------------------------------------------------");
  93   AliDebug(1,"");
  94   AliDebug(1,"    Entering InitMille");
  95   AliDebug(1,"");
  96   AliDebug(1,"-----------------------------------------------------");
  97   AliDebug(1,"");
  98
  99   fNGlobalConstraints = 0;
 100   fNLocalFits  = 0;                      // Total number of local fits
 101   fNLocalFitsRejected  = 0;              // Total number of local fits rejected
 102   fChi2CutRef  = 1.0;                    // Reference value for Chi^2/ndof cut
 103
 104   AliMillepede::SetNLocalEquations(0);     // Number of local fits (starts at 0)
 105
 106   fIter    = 0;  // By default iterations are turned off, turned on by use SetIterations
 107   fMaxIter = 10;
 108
 109   fResCutInit = lResCutInit;
 110   fResCut = lResCut;
 111
 112   fNGlobalPar = nGlo;     // Number of global derivatives
 113   fNLocalPar  = nLoc;     // Number of local derivatives
 114   fNStdDev    = lNStdDev; // Number of StDev for local fit chisquare cut
 115
 116   AliDebug(1,Form("Number of global parameters   : %d", fNGlobalPar));
 117   AliDebug(1,Form("Number of local parameters    : %d", fNLocalPar));
 118   AliDebug(1,Form("Number of standard deviations : %d", fNStdDev));
 119
 120   if (fNGlobalPar>fgkMaxGlobalPar || fNLocalPar>fgkMaxLocalPar) {
 121     AliDebug(1,"Two many parameters !!!!!");
 122     return 0;
 123   }
 124
 125   // Global parameters initializations
 126   for (int i=0; i<fNGlobalPar; i++) {
 127     fVecBGlo[i]=0.;
 128     fInitPar[i]=0.;
 129     fDeltaPar[i]=0.;
 130     fSigmaPar[i]=-1.;
 131     fIsNonLinear[i]=0;
 132     fGlo2CGLRow[i]=-1;
 133     fCGLRow2Glo[i]=-1;
 134
 135     for (int j=0; j<fNGlobalPar;j++) {
 136       fMatCGlo[i][j]=0.;
 137     }
 138   }
 139
 140   // Local parameters initializations
 141
 142   for (int i=0; i<fNLocalPar; i++) {
 143     fVecBLoc[i]=0.;
 144
 145     for (int j=0; j<fNLocalPar;j++) {
 146       fMatCLoc[i][j]=0.;
 147     }
 148   }
 149
 150   // Then we fix all parameters...
 151
 152   for (int j=0; j<fNGlobalPar; j++) {
 153     AliMillepede::SetParSigma(j,0.0);
 154   }
 155
 156   fDerivLocEq.Reset();  fNDerivLocEq=0;
 157   fIndexLocEq.Reset();  fNIndexLocEq=0;
 158
 159   fIndexAllEqs.Reset();  fNIndexAllEqs=0;
 160   fDerivAllEqs.Reset();  fNDerivAllEqs=0;
 161   fLocEqPlace.Reset();  fNLocEqPlace=0;
 162
 163   AliDebug(1,"");
 164   AliDebug(1,"----------------------------------------------------");
 165   AliDebug(1,"");
 166   AliDebug(1,"    InitMille has been successfully called!");
 167   AliDebug(1,"");
 168   AliDebug(1,"-----------------------------------------------------");
 169   AliDebug(1,"");
 170
 171   return 1;
 172 }
 173
 174 /*
 175 -----------------------------------------------------------
 176   PARGLO: initialization of global parameters
 177 -----------------------------------------------------------
 178
 179   param    = array of starting values
 180
 181 -----------------------------------------------------------
 182 */
 183 Int_t AliMillepede::SetGlobalParameters(double *param)
 184 {
 185   /// initialization of global parameters
 186   for(Int_t iPar=0; iPar<fNGlobalPar; iPar++){
 187     fInitPar[iPar] = param[iPar];
 188   }
 189
 190  return 1;
 191 }
 192
 193 /*
 194 -----------------------------------------------------------
 195   PARGLO: initialization of global parameters
 196 -----------------------------------------------------------
 197
 198   iPar    = the index of the global parameter in the
 199              result array (equivalent to fDeltaPar[]).
 200
 201   param    = the starting value
 202
 203 -----------------------------------------------------------
 204 */
 205 Int_t AliMillepede::SetGlobalParameter(int iPar, double param)
 206 {
 207   /// initialization of global parameter iPar
 208   if (iPar<0 || iPar>=fNGlobalPar) {
 209     return 0;
 210   }
 211   else {
 212     fInitPar[iPar] = param;
 213   }
 214   return 1;
 215 }
 216
 217
 218 /*
 219 -----------------------------------------------------------
 220   PARSIG: define a constraint for a single global param
 221           param is 'encouraged' to vary within [-sigma;sigma]
 222           range
 223 -----------------------------------------------------------
 224
 225   iPar    = the index of the global parameter in the
 226              result array (equivalent to fDeltaPar[]).
 227
 228   sigma    = value of the constraint (sigma <= 0. will
 229              mean that parameter is FIXED !!!)
 230
 231 -----------------------------------------------------------
 232 */
 233 Int_t AliMillepede::SetParSigma(int iPar, double sigma)
 234 {
 235   /// Define a range [-sigma;sigma] where iPar is encourage to vary
 236   if (iPar>=fNGlobalPar) {
 237     return 0;
 238   }
 239   else {
 240     fSigmaPar[iPar] = sigma;
 241   }
 242
 243   return 1;
 244 }
 245
 246
 247 /*
 248 -----------------------------------------------------------
 249   NONLIN: set nonlinear flag for a single global param
 250           update of param durin iterations will not
 251           consider initial starting value
 252 -----------------------------------------------------------
 253
 254   iPar    = the index of the global parameter in the
 255              result array (equivalent to fDeltaPar[]).
 256
 257 -----------------------------------------------------------
 258 */
 259 Int_t AliMillepede::SetNonLinear(int iPar)
 260 {
 261   /// Set nonlinear flag for iPar
 262   if (iPar<0 || iPar>=fNGlobalPar) {
 263     return 0;
 264   }
 265   else {
 266     fIsNonLinear[iPar] = 1;
 267   }
 268
 269   return 1;
 270 }
 271
 272
 273 /*
 274 -----------------------------------------------------------
 275   INITUN: unit for iteration
 276 -----------------------------------------------------------
 277
 278   lChi2CutFac is used by Fitloc to define the Chi^2/ndof cut value
 279
 280   A large cutfac value enables to take a wider range of tracks
 281   for first iterations, which might be useful if misalignments
 282   are large.
 283
 284   As soon as cutfac differs from 0 iteration are requested.
 285   cutfac is then reduced, from one iteration to the other,
 286   and iterations are stopped when it reaches the value 1.
 287
 288   At least one more iteration is often needed in order to remove
 289   tracks containing outliers.
 290
 291 -----------------------------------------------------------
 292 */
 293 Int_t AliMillepede::SetIterations(double lChi2CutFac)
 294 {
 295   /// Number of iterations is calculated from lChi2CutFac
 296   fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
 297
 298   AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
 299   fIter = 1; // Initializes the iteration process
 300   return 1;
 301 }
 302
 303 /*
 304 -----------------------------------------------------------
 305   CONSTF: define a constraint equation in AliMillepede
 306 -----------------------------------------------------------
 307
 308   dercs    = the row containing constraint equation
 309              derivatives (put into the final matrix)
 310
 311   lLagMult      = the lagrange multiplier value (sum of equation)
 312
 313 -----------------------------------------------------------
 314 */
 315 Int_t AliMillepede::SetGlobalConstraint(double dercs[], double lLagMult)
 316 {
 317   /// Define a constraint equation
 318   if (fNGlobalConstraints>=fgkMaxGloCsts) {
 319     AliInfo("Too many constraints !!!");
 320     return 0;
 321   }
 322
 323   for (int i=0; i<fNGlobalPar; i++) {
 324     fMatDerConstr[fNGlobalConstraints][i] = dercs[i];
 325   }
 326
 327   fLagMult[fNGlobalConstraints] = lLagMult;
 328   fNGlobalConstraints++ ;
 329   AliInfo(Form("Number of constraints increased to %d",fNGlobalConstraints));
 330   return 1;
 331 }
 332
 333
 334 /*
 335 -----------------------------------------------------------
 336   EQULOC: write ONE equation in the matrices
 337 -----------------------------------------------------------
 338
 339   dergb[1..fNGlobalPar] = global parameters derivatives
 340   derlc[1..fNLocalPar]  = local parameters derivatives
 341   rmeas                 = measured value
 342   sigma                 = error on measured value (nothing to do with SetParSigma!!!)
 343
 344 -----------------------------------------------------------
 345 */
 346 Int_t AliMillepede::SetLocalEquation(double dergb[], double derlc[], double lMeas, double lSigma)
 347 {
 348   /// Write one local equation
 349   if (lSigma<=0.0) { // If parameter is fixed, then no equation
 350     for (int i=0; i<fNLocalPar; i++) {
 351       derlc[i] = 0.0;
 352     }
 353     for (int i=0; i<fNGlobalPar; i++) {
 354       dergb[i] = 0.0;
 355     }
 356     return 1;
 357   }
 358
 359   // Serious equation, initialize parameters
 360   double lWeight =  1.0/(lSigma*lSigma);
 361   int iLocFirst  = -1;
 362   int iLocLast   = -1;
 363   int iGlobFirst = -1;
 364   int iGlobLast  = -1;
 365
 366   for (int i=0; i<fNLocalPar; i++) { // Retrieve local param interesting indices
 367     if (derlc[i]!=0.0) {
 368       if (iLocFirst == -1) {
 369         iLocFirst = i;  // first index
 370       }
 371       iLocLast = i;     // last index
 372     }
 373   }
 374   AliDebug(2,Form("%d / %d",iLocFirst, iLocLast));
 375
 376   for (int i=0; i<fNGlobalPar; i++) { // Idem for global parameters
 377     if (dergb[i]!=0.0) {
 378       if (iGlobFirst == -1) {
 379         iGlobFirst = i;  // first index
 380       }
 381       iGlobLast = i;     // last index
 382     }
 383   }
 384   AliDebug(2,Form("%d / %d",iGlobFirst,iGlobLast));
 385
 386   if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
 387   fIndexLocEq.AddAt(-1,fNIndexLocEq++);
 388   if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
 389   fDerivLocEq.AddAt(lMeas,fNDerivLocEq++);
 390
 391
 392   for (int i=iLocFirst; i<=iLocLast; i++) { // Store interesting local parameters
 393     if (derlc[i]!=0.0) {
 394       if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
 395       fIndexLocEq.AddAt(i,fNIndexLocEq++);
 396       if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
 397       fDerivLocEq.AddAt(derlc[i],fNDerivLocEq++);
 398       derlc[i] = 0.0;
 399     }
 400   }
 401
 402   if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
 403   fIndexLocEq.AddAt(-1,fNIndexLocEq++);
 404   if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
 405   fDerivLocEq.AddAt(lWeight,fNDerivLocEq++);
 406
 407   for (int i=iGlobFirst; i<=iGlobLast; i++) { // Store interesting global parameters
 408     if (dergb[i]!=0.0) {
 409       if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
 410       fIndexLocEq.AddAt(i,fNIndexLocEq++);
 411       if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
 412       fDerivLocEq.AddAt(dergb[i],fNDerivLocEq++);
 413       dergb[i] = 0.0;
 414     }
 415   }
 416
 417   AliDebug(2,Form("Out Equloc --  NST = %d",fNDerivLocEq));
 418   return 1;
 419 }
 420
 421 /*
 422 -----------------------------------------------------------
 423   FITLOC:  perform local params fit, once all the equations
 424            have been written by EquLoc
 425 -----------------------------------------------------------
 426
 427   iFit        = number of the fit, it is used to store
 428                 fit parameters and then retrieve them
 429                 for iterations (via FINDEXALLEQS and FDERIVALLEQS)
 430
 431   localParams = contains the fitted track parameters and
 432                 related errors
 433
 434   bSingleFit  = is an option, if it is set to 1, we don't
 435                 perform the last loop. It is used to update
 436                 the track parameters without modifying global
 437                 matrices
 438
 439 -----------------------------------------------------------
 440 */
 441 Int_t AliMillepede::LocalFit(int iFit, double localParams[], Bool_t bSingleFit)
 442 {
 443   /// Perform local parameters fit once all the local equations have been set
 444
 445   // Few initializations
 446   int iEqTerm = 0;
 447   int i, j, k;
 448   int iIdx, jIdx, kIdx;
 449   int iIdxIdx, jIdxIdx;
 450   int iMeas   = -1;
 451   int iWeight = 0;
 452   int iLocFirst = 0;
 453   int iLocLast  = 0;
 454   int iGloFirst = 0;
 455   int iGloLast  = 0;
 456   int nGloInFit = 0;
 457
 458   double lMeas    = 0.0;
 459   double lWeight  = 0.0;
 460
 461   double lChi2    = 0.0;
 462   double lRedChi2 = 0.0;
 463   double lChi2Cut = 0.0;
 464   int nEq  = 0;
 465   int nDoF = 0;
 466   int nEqTerms = fNDerivLocEq; // Number of terms (local + global derivatives,
 467                                // measurement and weight) involved in this local equation
 468
 469
 470   // Fill the track store at first pass
 471
 472   if (fIter < 2 && !bSingleFit) { // Do it only once
 473     AliDebug(1,Form("Store equation no: %d", iFit));
 474
 475     for (i=0; i<nEqTerms; i++) { // Store the track parameters
 476       if (fNIndexAllEqs==fIndexAllEqs.GetSize()) fIndexAllEqs.Set(2*fNIndexAllEqs);
 477       fIndexAllEqs.AddAt(fIndexLocEq[i],fNIndexAllEqs++);
 478       if (fNDerivAllEqs==fDerivAllEqs.GetSize()) fDerivAllEqs.Set(2*fNDerivAllEqs);
 479       fDerivAllEqs.AddAt(fDerivLocEq[i],fNDerivAllEqs++);
 480     }
 481     if (fNLocEqPlace==fLocEqPlace.GetSize()) fLocEqPlace.Set(2*fNLocEqPlace);
 482     fLocEqPlace.AddAt(fNIndexAllEqs,fNLocEqPlace++);
 483
 484
 485     AliDebug(2,Form("FLocEqPlace size = %d",fLocEqPlace[iFit]));
 486     AliDebug(2,Form("FIndexAllEqs size   = %d",fNIndexAllEqs));
 487   }
 488
 489
 490   for (i=0; i<fNLocalPar; i++) { // reset local params
 491     fVecBLoc[i] = 0.0;
 492
 493     for (j=0; j<fNLocalPar; j++) {
 494       fMatCLoc[i][j] = 0.0;
 495     }
 496   }
 497
 498   for (i=0; i<fNGlobalPar; i++) {
 499     fGlo2CGLRow[i] = -1;  // reset mixed params
 500   }
 501
 502
 503 /*
 504
 505   LOOPS : HOW DOES IT WORKS ?
 506
 507   Now start by reading the informations stored with EquLoc.
 508   Those informations are in vector FINDEXSTORE and FDERIVSTORE.
 509   Each -1 in FINDEXSTORE delimits the equation parameters:
 510
 511   First -1  ---> lMeas in FDERIVSTORE
 512   Then we have indices of local eq in FINDEXSTORE, and derivatives in FDERIVSTORE
 513   Second -1 ---> weight in FDERIVSTORE
 514   Then follows indices and derivatives of global eq.
 515   ....
 516
 517   We took them and store them into matrices.
 518
 519   As we want ONLY local params, we substract the part of the estimated value
 520   due to global params. Indeed we could have already an idea of these params,
 521   with previous alignment constants for example (set with PARGLO). Also if there
 522   are more than one iteration (FITLOC could be called by FITGLO)
 523
 524 */
 525
 526
 527 //
 528 // FIRST LOOP : local track fit
 529 //
 530
 531   iEqTerm = 0;
 532 //   fIndexLocEq.push_back(-1);
 533   if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
 534   fIndexLocEq.AddAt(-1,fNIndexLocEq++);
 535
 536   while (iEqTerm <= nEqTerms) {
 537     if (fIndexLocEq[iEqTerm] == -1) {
 538       if (iMeas == -1) {        // First  -1 : lMeas
 539         iMeas = iEqTerm;
 540         iLocFirst = iEqTerm+1;
 541       }
 542       else if (iWeight == 0) {  // Second -1 : weight
 543         iWeight = iEqTerm;
 544         iLocLast = iEqTerm-1;
 545         iGloFirst = iEqTerm+1;
 546       }
 547       else {                    // Third  -1 : end of equation; start of next
 548         iGloLast = iEqTerm-1;
 549         lMeas   = fDerivLocEq[iMeas];
 550         lWeight         = fDerivLocEq[iWeight];
 551 //      AliDebug(1,Form("lMeas = %f", lMeas));
 552 //      AliDebug(1,Form("lWeight = %f", lWeight));
 553
 554         // Now suppress the global part (only relevant with iterations)
 555         //
 556         for (i=iGloFirst; i<=iGloLast; i++) {
 557           iIdx = fIndexLocEq[i];              // Global param indice
 558 //        AliDebug(2,Form("fDeltaPar[%d] = %f", iIdx, fDeltaPar[iIdx]));
 559 //        AliDebug(2,Form("Starting misalignment = %f",fInitPar[iIdx]));
 560           if (fIsNonLinear[iIdx] == 0)
 561             lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]); // linear parameter
 562           else
 563             lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
 564         }
 565 //      AliDebug(2,Form("lMeas after global stuff removal = %f", lMeas));
 566
 567         for (i=iLocFirst; i<=iLocLast; i++) { // Finally fill local matrix and vector
 568           iIdx = fIndexLocEq[i];   // Local param indice (the matrix line)
 569           fVecBLoc[iIdx] += lWeight*lMeas*fDerivLocEq[i];
 570 //        AliDebug(2,Form("fVecBLoc[%d] = %f", iIdx, fVecBLoc[iIdx]));
 571
 572           for (j=iLocFirst; j<=i ; j++) { // Symmetric matrix, don't bother j>i coeffs
 573             jIdx = fIndexLocEq[j];
 574             fMatCLoc[iIdx][jIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[j];
 575 //          AliDebug(2,Form("fMatCLoc[%d][%d] = ", iIdx, jIdx, fMatCLoc[iIdx][jIdx]));
 576           }
 577         }
 578         iMeas   = -1;
 579         iWeight = 0;
 580         iEqTerm--;   // end of one equation is the beginning of next
 581       } // End of "end of equation" operations
 582     } // End of loop on equation
 583     iEqTerm++;
 584   } // End of loop on all equations used in the fit
 585
 586
 587 //
 588 // Local params matrix is completed, now invert to solve...
 589 //
 590
 591   Int_t nRank = AliMillepede::SpmInv(fMatCLoc, fVecBLoc, fNLocalPar);
 592   // nRank is the number of nonzero diagonal elements
 593
 594   AliDebug(1,"");
 595   AliDebug(1," __________________________________________________");
 596   AliDebug(1,Form(" Printout of local fit  (FITLOC)  with rank= %d", nRank));
 597   AliDebug(1," Result of local fit :      (index/parameter/error)");
 598
 599   for (i=0; i<fNLocalPar; i++) {
 600     AliDebug(1,Form("%d   /   %.6f   /   %.6f", i, fVecBLoc[i], TMath::Sqrt(fMatCLoc[i][i])));
 601   }
 602
 603
 604 // Store the track params and errors
 605
 606   for (i=0; i<fNLocalPar; i++) {
 607     localParams[2*i] = fVecBLoc[i];
 608     localParams[2*i+1] = TMath::Sqrt(TMath::Abs(fMatCLoc[i][i]));
 609   }
 610
 611
 612 //
 613 // SECOND LOOP : residual calculation
 614 //
 615
 616   iEqTerm = 0;
 617   iMeas = -1;
 618   iWeight = 0;
 619
 620   while (iEqTerm <= nEqTerms) {
 621     if (fIndexLocEq[iEqTerm] == -1) {
 622       if (iMeas == -1) {        // First  -1 : lMeas
 623         iMeas = iEqTerm;
 624         iLocFirst = iEqTerm+1;
 625       }
 626       else if (iWeight == 0) {  // Second -1 : weight
 627         iWeight = iEqTerm;
 628         iLocLast = iEqTerm-1;
 629         iGloFirst = iEqTerm+1;
 630       }
 631       else {                    // Third  -1 : end of equation; start of next
 632         iGloLast = iEqTerm-1;
 633         lMeas   = fDerivLocEq[iMeas];
 634         lWeight         = fDerivLocEq[iWeight];
 635
 636         // Print all (for debugging purposes)
 637
 638 //      int nDerLoc = iLocLast-iLocFirst+1;   // Number of local derivatives involved
 639 //      int nDerGlo = iGloLast-iGloFirst+1;   // Number of global derivatives involved
 640
 641 //      AliDebug(2,"");
 642 //      AliDebug(2,Form(". equation:  measured value %.6f +/- %.6f", lMeas, 1.0/TMath::Sqrt(lWeight)));
 643 //      AliDebug(2,Form("Number of derivatives (global, local): %d, %d",nDerGlo,nDerLoc));
 644 //      AliDebug(2,"Global derivatives are: (index/derivative/parvalue) ");
 645
 646 //      for (i=iGloFirst; i<=iGloLast; i++) {
 647 //        AliDebug(2,Form("%d / %.6f / %.6f",fIndexLocEq[i],fDerivLocEq[i],fInitPar[fIndexLocEq[i]]));
 648 //      }
 649
 650 //      AliDebug(2,"Local derivatives are: (index/derivative) ");
 651
 652 //      for (i=(ja+1); i<jb; i++) {AliDebug(2,Form("%d / %.6f",fIndexLocEq[i], fDerivLocEq[i]));}
 653
 654         // Now suppress local and global parts to LMEAS;
 655         //
 656         // First the local part
 657         for (i=iLocFirst; i<=iLocLast; i++) {
 658           iIdx = fIndexLocEq[i];
 659           lMeas -= fDerivLocEq[i]*fVecBLoc[iIdx];
 660         }
 661         // Then the global part
 662         for (i=iGloFirst; i<=iGloLast; i++) {
 663           iIdx = fIndexLocEq[i];
 664           if (fIsNonLinear[iIdx] == 0)
 665             lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]); // linear parameter
 666           else
 667             lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
 668         }
 669
 670         // lMeas contains now the residual value
 671         AliDebug(2,Form("Residual value : %.6f", lMeas));
 672
 673         // reject the track if lMeas is too important (outlier)
 674         if (TMath::Abs(lMeas) >= fResCutInit && fIter <= 1) {
 675           AliDebug(2,"Rejected track !!!!!");
 676           fNLocalFitsRejected++;
 677           fIndexLocEq.Reset();  fNIndexLocEq=0; // reset stores and go to the next track
 678           fDerivLocEq.Reset();  fNDerivLocEq=0;
 679           return 0;
 680         }
 681
 682         if (TMath::Abs(lMeas) >= fResCut && fIter > 1) {
 683           AliDebug(2,"Rejected track !!!!!");
 684           fNLocalFitsRejected++;
 685           fIndexLocEq.Reset();  fNIndexLocEq=0; // reset stores and go to the next track
 686           fDerivLocEq.Reset();  fNDerivLocEq=0;
 687           return 0;
 688         }
 689
 690         lChi2 += lWeight*lMeas*lMeas ; // total chi^2
 691         nEq++;                    // number of equations
 692         iMeas   = -1;
 693         iWeight = 0;
 694         iEqTerm--;
 695       } // End of "end of equation" operations
 696     }   // End of loop on equation
 697     iEqTerm++;
 698   } // End of loop on all equations used in the fit
 699
 700   nDoF = nEq-nRank;
 701   lRedChi2 = 0.0;
 702
 703   AliDebug(1,Form("Final chi square / degrees of freedom %.2f / %d",lChi2, nDoF));
 704
 705   if (nDoF > 0) lRedChi2 = lChi2/float(nDoF);  // Chi^2/dof
 706
 707   fNLocalFits++;
 708
 709   if (fNStdDev != 0 && nDoF > 0 && !bSingleFit) // Chisquare cut
 710   {
 711     lChi2Cut = AliMillepede::Chi2DoFLim(fNStdDev, nDoF)*fChi2CutFactor;
 712
 713     AliDebug(1,Form("Reject if Chisq/Ndf = %.4f > %.4f",lRedChi2,lChi2Cut));
 714
 715     if (lRedChi2 > lChi2Cut) // Reject the track if too much...
 716     {
 717       AliDebug(2,"Rejected track !!!!!");
 718       fNLocalFitsRejected++;
 719       fIndexLocEq.Reset();  fNIndexLocEq=0; // reset stores and go to the next track
 720       fDerivLocEq.Reset();  fNDerivLocEq=0;
 721       return 0;
 722     }
 723   }
 724
 725   if (bSingleFit) // Stop here if just updating the track parameters
 726   {
 727     fIndexLocEq.Reset();  fNIndexLocEq=0; // Reset store for the next track
 728     fDerivLocEq.Reset();  fNDerivLocEq=0;
 729     return 1;
 730   }
 731
 732 //
 733 // THIRD LOOP: local operations are finished, track is accepted
 734 // We now update the global parameters (other matrices)
 735 //
 736
 737   iEqTerm = 0;
 738   iMeas = -1;
 739   iWeight = 0;
 740
 741   while (iEqTerm <= nEqTerms)
 742   {
 743     if (fIndexLocEq[iEqTerm] == -1)
 744     {
 745       if (iMeas == -1) {        // First  -1 : lMeas
 746         iMeas = iEqTerm;
 747         iLocFirst = iEqTerm+1;
 748       }
 749       else if (iWeight == 0) {  // Second -1 : weight
 750         iWeight = iEqTerm;
 751         iLocLast = iEqTerm-1;
 752         iGloFirst = iEqTerm+1;
 753       }
 754       else {                    // Third  -1 : end of equation; start of next
 755         iGloLast = iEqTerm-1;
 756         lMeas   = fDerivLocEq[iMeas];
 757         lWeight         = fDerivLocEq[iWeight];
 758
 759         // Now suppress the global part
 760         for (i=iGloFirst; i<=iGloLast; i++) {
 761           iIdx = fIndexLocEq[i];   // Global param indice
 762           if (fIsNonLinear[iIdx] == 0)
 763             lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]); // linear parameter
 764           else
 765             lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
 766         }
 767
 768         for (i=iGloFirst; i<=iGloLast; i++) {
 769           iIdx = fIndexLocEq[i];   // Global param indice (the matrix line)
 770
 771           fVecBGlo[iIdx] += lWeight*lMeas*fDerivLocEq[i];
 772 //        AliDebug(2,Form("fVecBGlo[%d] = %.6f", j, fVecBGlo[j] ));
 773
 774           // First of all, the global/global terms (exactly like local matrix)
 775           //
 776           for (j=iGloFirst; j<=iGloLast; j++) {
 777             jIdx = fIndexLocEq[j];
 778             fMatCGlo[iIdx][jIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[j];
 779 //          AliDebug(2,Form("fMatCGlo[%d][%d] = %.6f",iIdx,jIdx,fMatCGlo[iIdx][jIdx]));
 780           }
 781
 782           // Now we have also rectangular matrices containing global/local terms.
 783           //
 784           iIdxIdx = fGlo2CGLRow[iIdx];  // Index of index
 785           if (iIdxIdx == -1) {    // New global variable
 786             for (k=0; k<fNLocalPar; k++) {
 787               fMatCGloLoc[nGloInFit][k] = 0.0;  // Initialize the row
 788             }
 789             fGlo2CGLRow[iIdx] = nGloInFit;
 790             fCGLRow2Glo[nGloInFit] = iIdx;
 791             iIdxIdx = nGloInFit;
 792             nGloInFit++;
 793           }
 794
 795           // Now fill the rectangular matrix
 796           for (k=iLocFirst; k<=iLocLast ; k++) {
 797             kIdx = fIndexLocEq[k];
 798             fMatCGloLoc[iIdxIdx][kIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[k];
 799 //          AliDebug(2,Form("fMatCGloLoc[%d][%d] = %.6f",iIdxIdx,kIdx,fMatCGloLoc[iIdxIdx][kIdx]));
 800           }
 801         }
 802         iMeas   = -1;
 803         iWeight =  0;
 804         iEqTerm--;
 805       } // End of "end of equation" operations
 806     }   // End of loop on equation
 807     iEqTerm++;
 808   } // End of loop on all equations used in the fit
 809
 810   // Third loop is finished, now we update the correction matrices
 811   AliMillepede::SpAVAt(fMatCLoc, fMatCGloLoc, fMatCGloCorr, fNLocalPar, nGloInFit);
 812   AliMillepede::SpAX(fMatCGloLoc, fVecBLoc, fVecBGloCorr, fNLocalPar, nGloInFit);
 813
 814   for (iIdxIdx=0; iIdxIdx<nGloInFit; iIdxIdx++) {
 815     iIdx = fCGLRow2Glo[iIdxIdx];
 816     fVecBGlo[iIdx] -= fVecBGloCorr[iIdxIdx];
 817
 818     for (jIdxIdx=0; jIdxIdx<=iIdxIdx; jIdxIdx++) {
 819       int jIdx = fCGLRow2Glo[jIdxIdx];
 820       fMatCGlo[iIdx][jIdx] -= fMatCGloCorr[iIdxIdx][jIdxIdx];
 821       fMatCGlo[jIdx][iIdx] = fMatCGlo[iIdx][jIdx];
 822     }
 823   }
 824
 825   fIndexLocEq.Reset();  fNIndexLocEq=0; // Reset store for the next track
 826   fDerivLocEq.Reset();  fNDerivLocEq=0;
 827
 828   return 1;
 829 }
 830
 831
 832 /*
 833 -----------------------------------------------------------
 834   MAKEGLOBALFIT:  perform global params fit, once all the 'tracks'
 835                   have been fitted by FitLoc
 836 -----------------------------------------------------------
 837
 838   par[]        = array containing the computed global
 839                  parameters (the misalignment constants)
 840
 841   error[]      = array containing the error on global
 842                  parameters (estimated by AliMillepede)
 843
 844   pull[]        = array containing the corresponding pulls
 845
 846 -----------------------------------------------------------
 847 */
 848 Int_t AliMillepede::GlobalFit(double par[], double error[], double pull[])
 849 {
 850   /// perform global parameters fit once all the local equations have been fitted
 851   int i, j;
 852   int nVar    = 0;
 853   int nGloFix = 0;
 854   double lConstraint;
 855
 856   double step[fgkMaxGlobalPar];
 857
 858   double localPars[2*fgkMaxLocalPar];
 859
 860   int nLocFitsGood = 0;
 861   int nLocFitsTot  = 0;
 862   int nLocFits     = 0;
 863
 864   AliInfo("..... Making global fit .....");
 865   AliInfo(Form("First Iteration with cut factor %.2f", fChi2CutFactor));
 866
 867   nLocFitsTot = AliMillepede::GetNLocalEquations();
 868
 869   while (fIter < fMaxIter)  // Iteration for the final loop
 870   {
 871
 872     nLocFits = AliMillepede::GetNLocalEquations();
 873     AliInfo(Form("...using %d local fits...",nLocFits));
 874
 875 // Start by saving the diagonal elements
 876
 877     for (i=0; i<fNGlobalPar; i++) {
 878       fDiagCGlo[i] = fMatCGlo[i][i];
 879     }
 880
 881 //  Then we retrieve the different constraints: fixed parameter or global equation
 882
 883     nGloFix = 0; // First look at the fixed global params
 884
 885     for (i=0; i<fNGlobalPar; i++) {
 886       if (fSigmaPar[i] <= 0.0) {  // fixed global param
 887         nGloFix++;
 888         for (j=0; j<fNGlobalPar; j++) {
 889           fMatCGlo[i][j] = 0.0;  // Reset row and column
 890           fMatCGlo[j][i] = 0.0;
 891         }
 892       }
 893       else {
 894         fMatCGlo[i][i] += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
 895       }
 896     }
 897
 898     nVar = fNGlobalPar;  // Current number of equations
 899     AliDebug(1,Form("Number of constraint equations : %d", fNGlobalConstraints));
 900
 901     for (i=0; i<fNGlobalConstraints; i++) { // Then the constraint equation
 902       lConstraint = fLagMult[i];
 903       for (j=0; j<fNGlobalPar; j++) {
 904         fMatCGlo[nVar][j] = float(nLocFits)*fMatDerConstr[i][j];
 905         fMatCGlo[j][nVar] = float(nLocFits)*fMatDerConstr[i][j];
 906         lConstraint -= fMatDerConstr[i][j]*(fInitPar[j]+fDeltaPar[j]);
 907       }
 908
 909       fMatCGlo[nVar][nVar] = 0.0;
 910       fVecBGlo[nVar] = float(nLocFits)*lConstraint;
 911       nVar++;
 912     }
 913
 914
 915     // Intended to compute the final global chisquare
 916
 917     double lFinalCor = 0.0;
 918
 919     if (fIter > 1) {
 920       for (i=0; i<fNGlobalPar; i++) {
 921         for (j=0; j<fNGlobalPar; j++) {
 922 //        printf("%d, %d, %.6f  %.6f  %.6f\n",i,j,step[i],fMatCGlo[i][j],step[j]);
 923           lFinalCor += step[i]*fMatCGlo[i][j]*step[j];
 924           if (i == j && fSigmaPar[i] != 0) {
 925             lFinalCor -= step[i]*step[i]/(fSigmaPar[i]*fSigmaPar[i]);
 926           }
 927         }
 928       }
 929     }
 930
 931     AliInfo(Form(" Final coeff is %.6f",lFinalCor));
 932     AliInfo(Form(" Final NDOFs = %d", fNGlobalPar));
 933
 934     //  The final matrix inversion
 935
 936     Int_t nRank = AliMillepede::SpmInv(fMatCGlo, fVecBGlo, nVar);
 937
 938     for (i=0; i<fNGlobalPar; i++) {
 939       fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
 940       AliDebug(1,Form("fDeltaPar[%d] = %.6f", i, fDeltaPar[i]));
 941       AliDebug(1,Form("fMatCGlo[%d][%d] = %.6f", i, i, fMatCGlo[i][i]));
 942       AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
 943
 944       step[i] = fVecBGlo[i];
 945
 946       if (fIter == 1) error[i] = fMatCGlo[i][i]; // Unfitted error
 947     }
 948     AliInfo(Form("The rank defect of the symmetric %d by %d matrix is %d (bad if non 0)",
 949             nVar, nVar, nVar-nGloFix-nRank));
 950
 951     AliInfo(Form("Total : %d local fits, %d rejected.", fNLocalFits, fNLocalFitsRejected));
 952     if (fIter == 0)  break;  // No iterations set
 953     fIter++;
 954
 955     if (fIter == fMaxIter)  break;  // End of story
 956
 957     // Reinitialize parameters for iteration
 958     //
 959     fNLocalFits = 0;
 960     fNLocalFitsRejected = 0;
 961
 962     if (fChi2CutFactor != fChi2CutRef) {
 963       fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
 964       if (fChi2CutFactor < 1.2*fChi2CutRef) {
 965         fChi2CutFactor = fChi2CutRef;
 966         fIter = fMaxIter - 1;     // Last iteration
 967       }
 968     }
 969     AliInfo("");
 970     AliInfo("----------------------------------------------");
 971     AliInfo(Form("Iteration %d with cut factor %.2f", fIter, fChi2CutFactor));
 972
 973     // Reset global variables
 974     //
 975     for (i=0; i<nVar; i++) {
 976       fVecBGlo[i] = 0.0;
 977       for (j=0; j<nVar; j++) {
 978         fMatCGlo[i][j] = 0.0;
 979       }
 980     }
 981
 982     //
 983     // We start a new iteration
 984     //
 985
 986     // First we read the stores for retrieving the local params
 987     //
 988     nLocFitsGood = 0;
 989
 990     for (i=0; i<nLocFitsTot; i++) {
 991       int iEqFirst = 0;
 992       int iEqLast = 0;
 993
 994       (i>0) ? iEqFirst = fLocEqPlace[i-1] : iEqFirst = 0;
 995       iEqLast = fLocEqPlace[i];
 996
 997       AliDebug(2,Form("Track %d : ",i));
 998       AliDebug(2,Form("Starts at %d", iEqFirst));
 999       AliDebug(2,Form("Ends at %d",iEqLast));
1000
1001       if (fIndexAllEqs[iEqFirst] != -999) { // Fit is still OK
1002         fIndexLocEq.Reset();  fNIndexLocEq=0;
1003         fDerivLocEq.Reset();  fNDerivLocEq=0;
1004
1005         for (j=iEqFirst; j<iEqLast; j++) {
1006           if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
1007           fIndexLocEq.AddAt(fIndexAllEqs[j],fNIndexLocEq++);
1008           if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
1009           fDerivLocEq.AddAt(fDerivAllEqs[j],fNDerivLocEq++);
1010         }
1011
1012         for (j=0; j<2*fNLocalPar; j++) {
1013           localPars[j] = 0.;
1014         }
1015
1016 //      Int_t sc = AliMillepede::LocalFit(i,localPars,0);
1017 //      (sc) ? nLocFitsGood++ : fIndexAllEqs[iEqFirst] = -999;
1018         AliMillepede::LocalFit(i,localPars,0);
1019         nLocFitsGood++;
1020       }
1021     } // End of loop on fits
1022
1023     AliMillepede::SetNLocalEquations(nLocFitsGood);
1024
1025   } // End of iteration loop
1026
1027   AliMillepede::PrintGlobalParameters(); // Print the final results
1028
1029   for (j=0; j<fNGlobalPar; j++) {
1030     par[j]   = fInitPar[j]+fDeltaPar[j];
1031     pull[j]  = (fSigmaPar[j] <= 0. || fSigmaPar[j]*fSigmaPar[j]-fMatCGlo[j][j] <=0.) ? 0. : fDeltaPar[j]/TMath::Sqrt(fSigmaPar[j]*fSigmaPar[j]-fMatCGlo[j][j]);
1032     error[j] = TMath::Sqrt(TMath::Abs(fMatCGlo[j][j]));
1033   }
1034
1035   AliInfo(" ");
1036   AliInfo("            * o o                   o      ");
1037   AliInfo("              o o                   o      ");
1038   AliInfo("   o ooooo  o o o  oo  ooo   oo   ooo  oo  ");
1039   AliInfo("    o  o  o o o o o  o o  o o  o o  o o  o ");
1040   AliInfo("    o  o  o o o o oooo o  o oooo o  o oooo ");
1041   AliInfo("    o  o  o o o o o    ooo  o    o  o o    ");
1042   AliInfo("    o  o  o o o o  oo  o     oo   ooo  oo ++ ends.");
1043   AliInfo("                       o                   ");
1044
1045   return 1;
1046 }
1047
1048 /*
1049 -----------------------------------------------------------
1050   ERRPAR: return error for parameter iPar
1051 -----------------------------------------------------------
1052
1053   iPar     = the index of the global parameter in the
1054              result array (equivalent to fDeltaPar[]).
1055
1056 -----------------------------------------------------------
1057 */
1058 Double_t AliMillepede::GetParError(int iPar) const
1059 {
1060   /// return error for parameter iPar
1061   Double_t lErr = -1.;
1062   if (iPar>=0 && iPar<fNGlobalPar) {
1063     lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
1064   }
1065   return lErr;
1066 }
1067
1068
1069 /*
1070 -----------------------------------------------------------
1071   SPMINV:  obtain solution of a system of linear equations with symmetric matrix
1072            and the inverse (using 'singular-value friendly' GAUSS pivot)
1073 -----------------------------------------------------------
1074
1075         Solve the equation :  V * X = B
1076
1077         V is replaced by inverse matrix and B by X, the solution vector
1078 -----------------------------------------------------------
1079 */
1080 int AliMillepede::SpmInv(double matV[][fgkMaxGloPC], double vecB[], int nGlo)
1081 {
1082   ///  Obtain solution of a system of linear equations with symmetric matrix
1083   ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
1084
1085   Int_t nRank = 0;
1086   int iPivot;
1087   double vPivot = 0.;
1088   double eps = 0.00000000000001;
1089
1090   bool *bUnUsed = new bool[nGlo];
1091   double *diagV = new double[nGlo];
1092   double *rowMax = new double[nGlo];
1093   double *colMax = new double[nGlo];
1094
1095   double *temp = new double[nGlo];
1096
1097   for (Int_t i=0; i<nGlo; i++) {
1098     rowMax[i] = 0.0;
1099     colMax[i] = 0.0;
1100     bUnUsed[i] = true;
1101
1102     for (Int_t j=0; j<i; j++) {
1103       if (matV[j][i] == 0) {
1104         matV[j][i] = matV[i][j];
1105       }
1106     }
1107   }
1108
1109   // Small loop for matrix equilibration (gives a better conditioning)
1110
1111   for (Int_t i=0; i<nGlo; i++) {
1112     for (Int_t j=0; j<nGlo; j++) {
1113       if (TMath::Abs(matV[i][j]) >= rowMax[i]) rowMax[i] = TMath::Abs(matV[i][j]); // Max elemt of row i
1114       if (TMath::Abs(matV[j][i]) >= colMax[i]) colMax[i] = TMath::Abs(matV[j][i]); // Max elemt of column i
1115     }
1116   }
1117
1118   for (Int_t i=0; i<nGlo; i++) {
1119     if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
1120     if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
1121   }
1122
1123   for (Int_t i=0; i<nGlo; i++) {
1124     for (Int_t j=0; j<nGlo; j++) {
1125       matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]); // Equilibrate the V matrix
1126     }
1127     diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
1128   }
1129
1130
1131   for (Int_t i=0; i<nGlo; i++) {
1132     vPivot = 0.0;
1133     iPivot = -1;
1134
1135     for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element
1136       if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1137         vPivot = matV[j][j];
1138         iPivot = j;
1139       }
1140     }
1141
1142     if (iPivot >= 0) {   // pivot found
1143       nRank++;
1144       bUnUsed[iPivot] = false; // This value is used
1145       vPivot = 1.0/vPivot;
1146       matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1147
1148       for (Int_t j=0; j<nGlo; j++) {
1149         for (Int_t jj=0; jj<nGlo; jj++) {
1150           if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
1151             matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1152           }
1153         }
1154       }
1155
1156       for (Int_t j=0; j<nGlo; j++) {
1157         if (j != iPivot) { // Pivot row or column elements
1158           matV[j][iPivot] = matV[j][iPivot]*vPivot;     // Column
1159           matV[iPivot][j] = matV[iPivot][j]*vPivot;     // Line
1160         }
1161       }
1162     }
1163     else {  // No more pivot value (clear those elements)
1164       for (Int_t j=0; j<nGlo; j++) {
1165         if (bUnUsed[j]) {
1166           vecB[j] = 0.0;
1167
1168           for (Int_t k=0; k<nGlo; k++) {
1169             matV[j][k] = 0.0;
1170             matV[k][j] = 0.0;
1171           }
1172         }
1173       }
1174       break;  // No more pivots anyway, stop here
1175     }
1176   }
1177
1178   for (Int_t i=0; i<nGlo; i++) {
1179     for (Int_t j=0; j<nGlo; j++) {
1180       matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]); // Correct matrix V
1181     }
1182   }
1183
1184   for (Int_t j=0; j<nGlo; j++) {
1185     temp[j] = 0.0;
1186
1187     for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
1188       matV[j][jj] = -matV[j][jj];
1189       temp[j] += matV[j][jj]*vecB[jj];
1190     }
1191   }
1192
1193   for (Int_t j=0; j<nGlo; j++) {
1194     vecB[j] = temp[j]; // The final result
1195   }
1196
1197   delete [] temp;
1198   delete [] bUnUsed;
1199   delete [] diagV;
1200   delete [] rowMax;
1201   delete [] colMax;
1202
1203   return nRank;
1204 }
1205
1206 //
1207 // Same method but for local fit, so heavily simplified
1208 //
1209 int AliMillepede::SpmInv(double matV[][fgkMaxLocalPar], double vecB[], int nLoc)
1210 {
1211   ///  Obtain solution of a system of linear equations with symmetric matrix
1212   ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
1213
1214   Int_t nRank = 0;
1215   Int_t iPivot = -1;
1216   double vPivot = 0.;
1217   double eps = 0.0000000000001;
1218
1219   bool *bUnUsed = new bool[nLoc];
1220   double *diagV  = new double[nLoc];
1221   double *temp  = new double[nLoc];
1222
1223   for (Int_t i=0; i<nLoc; i++) {
1224     bUnUsed[i] = true;
1225     diagV[i] = TMath::Abs(matV[i][i]);     // save diagonal elem absolute values
1226     for (Int_t j=0; j<i; j++) {
1227       matV[j][i] = matV[i][j] ;
1228     }
1229   }
1230
1231   for (Int_t i=0; i<nLoc; i++) {
1232     vPivot = 0.0;
1233     iPivot = -1;
1234
1235     for (Int_t j=0; j<nLoc; j++) { // First look for the pivot, ie max unused diagonal element
1236       if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1237         vPivot = matV[j][j];
1238         iPivot = j;
1239       }
1240     }
1241
1242     if (iPivot >= 0) {   // pivot found
1243       nRank++;
1244       bUnUsed[iPivot] = false;
1245       vPivot = 1.0/vPivot;
1246       matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1247
1248       for (Int_t j=0; j<nLoc; j++) {
1249         if (j != iPivot) {
1250           for (Int_t jj=0; jj<=j; jj++) {
1251             if (jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
1252               matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1253               matV[jj][j] = matV[j][jj];
1254             }
1255           }
1256         }
1257       }
1258
1259       for (Int_t j=0; j<nLoc; j++) {
1260         if (j != iPivot) {      // Pivot row or column elements
1261           matV[j][iPivot] = matV[j][iPivot]*vPivot; // Column
1262           matV[iPivot][j] = matV[j][iPivot];
1263         }
1264       }
1265     }
1266     else { // No more pivot value (clear those elements)
1267       for (Int_t j=0; j<nLoc; j++) {
1268         if (bUnUsed[j]) {
1269           vecB[j] = 0.0;
1270           matV[j][j] = 0.0;
1271           for (Int_t k=0; k<j; k++) {
1272             matV[j][k] = 0.0;
1273             matV[k][j] = 0.0;
1274           }
1275         }
1276       }
1277
1278       break;  // No more pivots anyway, stop here
1279     }
1280   }
1281
1282   for (Int_t j=0; j<nLoc; j++) {
1283     temp[j] = 0.0;
1284     for (Int_t jj=0; jj<nLoc; jj++) { // Reverse matrix elements
1285       matV[j][jj] = -matV[j][jj];
1286       temp[j] += matV[j][jj]*vecB[jj];
1287     }
1288   }
1289
1290   for (Int_t j=0; j<nLoc; j++) {
1291     vecB[j] = temp[j];
1292   }
1293
1294   delete [] bUnUsed;
1295   delete [] diagV;
1296   delete [] temp;
1297
1298   return nRank;
1299 }
1300
1301
1302 /*
1303 -----------------------------------------------------------
1304   SPAVAT
1305 -----------------------------------------------------------
1306
1307   multiply symmetric N-by-N matrix from the left with general M-by-N
1308   matrix and from the right with the transposed of the same  general
1309   matrix  to  form  symmetric  M-by-M   matrix.
1310
1311                                                        T
1312   CALL SPAVAT(V,A,W,N,M)      W   =   A   *   V   *   A
1313                            M*M     M*N     N*N     N*M
1314
1315   where V = symmetric N-by-N matrix
1316         A = general N-by-M matrix
1317         W = symmetric M-by-M matrix
1318 -----------------------------------------------------------
1319 */
1320 Int_t AliMillepede::SpAVAt(double matV[][fgkMaxLocalPar], double matA[][fgkMaxLocalPar], double matW[][fgkMaxGlobalPar], int nLoc, int nGlo)
1321 {
1322   ///  multiply symmetric N-by-N matrix from the left with general M-by-N
1323   ///  matrix and from the right with the transposed of the same general
1324   ///  matrix to form symmetric M-by-M matrix.
1325
1326   for (Int_t i=0; i<nGlo; i++) {
1327     for (Int_t j=0; j<=i; j++) {  // Matrix w is symmetric
1328       matW[i][j] = 0.0;      // Reset final matrix
1329       for (Int_t k=0; k<nLoc; k++) {
1330         matW[i][j] += matA[i][k]*matV[k][k]*matA[j][k];    // diagonal terms of v
1331         for (Int_t l=0; l<k; l++) {
1332           matW[i][j] += matA[i][k]*matV[k][l]*matA[j][l];  // Use symmetric properties
1333           matW[i][j] += matA[i][l]*matV[k][l]*matA[j][k];  // of matrix v
1334         }
1335       }
1336       if (i!=j){
1337         matW[j][i] = matW[i][j]; // Matrix w is symmetric
1338       }
1339     }
1340   }
1341
1342   return 1;
1343 }
1344
1345
1346 /*
1347 -----------------------------------------------------------
1348   SPAX
1349 -----------------------------------------------------------
1350
1351   multiply general M-by-N matrix A and N-vector X
1352
1353   CALL  SPAX(A,X,Y,M,N)          Y =  A * X
1354                                  M   M*N  N
1355
1356   where A = general M-by-N matrix (A11 A12 ... A1N  A21 A22 ...)
1357         X = N vector
1358         Y = M vector
1359 -----------------------------------------------------------
1360 */
1361 Int_t AliMillepede::SpAX(double matA[][fgkMaxLocalPar], double vecX[], double vecY[], int nCol, int nRow)
1362 {
1363   ///   multiply general M-by-N matrix A and N-vector X
1364   for (Int_t i=0; i<nRow; i++) {
1365     vecY[i] = 0.0;          // Reset final vector
1366     for (Int_t j=0; j<nCol; j++) {
1367       vecY[i] += matA[i][j]*vecX[j];  // fill the vector
1368     }
1369   }
1370
1371   return 1;
1372 }
1373
1374
1375 /*
1376 -----------------------------------------------------------
1377   PRTGLO
1378 -----------------------------------------------------------
1379
1380   Print the final results into the logfile
1381
1382 -----------------------------------------------------------
1383 */
1384 Int_t AliMillepede::PrintGlobalParameters() const
1385 {
1386   ///  Print the final results into the logfile
1387   double lError = 0.;
1388   double lGlobalCor =0.;
1389
1390   AliInfo("");
1391   AliInfo("   Result of fit for global parameters");
1392   AliInfo("   ===================================");
1393   AliInfo("    I       initial       final       differ        lastcor        error       gcor");
1394   AliInfo("-----------------------------------------------------------------------------------");
1395
1396   for (int i=0; i<fNGlobalPar; i++) {
1397     lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
1398     if (fMatCGlo[i][i] < 0.0) lError = -lError;
1399     lGlobalCor = 0.0;
1400
1401     if (TMath::Abs(fMatCGlo[i][i]*fDiagCGlo[i]) > 0) {
1402       lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(fMatCGlo[i][i]*fDiagCGlo[i])));
1403       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f",
1404                    i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor));
1405     }
1406     else {
1407       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t OFF\t OFF",i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i]));
1408     }
1409   }
1410   return 1;
1411 }
1412
1413
1414 /*
1415 ----------------------------------------------------------------
1416   CHI2DOFLIM:  return the limit in chi^2/nd for n sigmas stdev authorized
1417 ----------------------------------------------------------------
1418
1419   Only n=1, 2, and 3 are expected in input
1420 ----------------------------------------------------------------
1421 */
1422 double AliMillepede::Chi2DoFLim(int nSig, int nDoF)
1423 {
1424   /// return the limit in chi^2/nd for n sigmas stdev authorized
1425   int lNSig;
1426   double sn[3]        = {0.47523, 1.690140, 2.782170};
1427   double table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
1428                           1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
1429                           1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
1430                           1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
1431                          {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
1432                           1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
1433                           1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
1434                           1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
1435                          {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
1436                           2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
1437                           2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
1438                           1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
1439
1440   if (nDoF < 1) {
1441     return 0.0;
1442   }
1443   else {
1444     lNSig = TMath::Max(1,TMath::Min(nSig,3));
1445
1446     if (nDoF <= 30) {
1447       return table[lNSig-1][nDoF-1];
1448     }
1449     else { // approximation
1450       return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
1451               (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);
1452     }
1453   }
1454 }