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