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