STEER/AliMillePede2.cxx

   1 /**********************************************************************************************/
   2 /* General class for alignment with large number of degrees of freedom                        */
   3 /* Based on the original milliped2 by Volker Blobel                                           */
   4 /* http://www.desy.de/~blobel/mptalks.html                                                    */
   5 /*                                                                                            */
   6 /* Author: ruben.shahoyan@cern.ch                                                             */
   7 /*                                                                                            */
   8 /**********************************************************************************************/
   9
  10 #include "AliMillePede2.h"
  11 #include "AliLog.h"
  12 #include <TStopwatch.h>
  13 #include <TFile.h>
  14 #include <TMath.h>
  15 #include <TVectorD.h>
  16 #include "AliMatrixSq.h"
  17 #include "AliSymMatrix.h"
  18 #include "AliRectMatrix.h"
  19 #include "AliMatrixSparse.h"
  20
  21
  22 using namespace std;
  23
  24
  25 ClassImp(AliMillePede2)
  26
  27 Bool_t   AliMillePede2::fgInvChol        = kTRUE;     // Invert global matrix with Cholesky solver
  28 Bool_t   AliMillePede2::fgWeightSigma    = kTRUE;     // weight local constraint by module statistics
  29 Bool_t   AliMillePede2::fgIsMatGloSparse = kFALSE;    // use faster dense matrix by default
  30 Int_t    AliMillePede2::fgMinResCondType = 1;         // Jacoby preconditioner by default
  31 Double_t AliMillePede2::fgMinResTol      = 1.e-11;    // default tolerance
  32 Int_t    AliMillePede2::fgMinResMaxIter  = 10000;     // default max number of iterations
  33 Int_t    AliMillePede2::fgIterSol        = AliMinResSolve::kSolMinRes; // default iterative solver
  34 Int_t    AliMillePede2::fgNKrylovV       = 240;        // default number of Krylov vectors to keep
  35
  36 //_____________________________________________________________________________________________
  37 AliMillePede2::AliMillePede2()
  38 : fNLocPar(0),
  39   fNGloPar(0),
  40   fNGloSize(0),
  41 //
  42   fNLocEquations(0),
  43   fIter(0),
  44   fMaxIter(10),
  45   fNStdDev(3),
  46   fNGloConstraints(0),
  47   fNLagrangeConstraints(0),
  48   fNLocFits(0),
  49   fNLocFitsRejected(0),
  50   fNGloFix(0),
  51   fGloSolveStatus(gkFailed),
  52 //
  53   fChi2CutFactor(1.),
  54   fChi2CutRef(1.),
  55   fResCutInit(100.),
  56   fResCut(100.),
  57   fMinPntValid(1),
  58 //
  59   fNGroupsSet(0),
  60   fParamGrID(0),
  61   fProcPnt(0),
  62   fVecBLoc(0),
  63   fDiagCGlo(0),
  64   fVecBGlo(0),
  65   fInitPar(0),
  66   fDeltaPar(0),
  67   fSigmaPar(0),
  68   fIsLinear(0),
  69   fConstrUsed(0),
  70 //
  71   fGlo2CGlo(0),
  72   fCGlo2Glo(0),
  73 //
  74   fMatCLoc(0),
  75   fMatCGlo(0),
  76   fMatCGloLoc(0),
  77   //
  78   fFillIndex(0),
  79   fFillValue(0),
  80   //
  81   fDataRecFName("/tmp/mp2_data_records.root"),
  82   fRecord(0),
  83   fDataRecFile(0),
  84   fTreeData(0),
  85   //
  86   fConstrRecFName("/tmp/mp2_constraints_records.root"),
  87   fTreeConstr(0),
  88   fConsRecFile(0),
  89   fCurrRecDataID(0),
  90   fCurrRecConstrID(0),
  91   fLocFitAdd(kTRUE)
  92 {}
  93
  94 //_____________________________________________________________________________________________
  95 AliMillePede2::AliMillePede2(const AliMillePede2& src) :
  96   TObject(src),fNLocPar(0),fNGloPar(0),fNGloSize(0),fNLocEquations(0),fIter(0),
  97   fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLagrangeConstraints(0),
  98   fNLocFits(0),fNLocFitsRejected(0),
  99   fNGloFix(0),fGloSolveStatus(0),fChi2CutFactor(0),fChi2CutRef(0),fResCutInit(0),
 100   fResCut(0),fMinPntValid(1),fNGroupsSet(0),fParamGrID(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0),
 101   fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0),
 102   fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0),
 103   fDataRecFName(0),fRecord(0),fDataRecFile(0),
 104   fTreeData(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
 105   fCurrRecConstrID(0),fLocFitAdd(kTRUE)
 106 {printf("Dummy\n");}
 107
 108 //_____________________________________________________________________________________________
 109 AliMillePede2::~AliMillePede2()
 110 {
 111   CloseDataRecStorage();
 112   CloseConsRecStorage();
 113   //
 114   delete[] fParamGrID;
 115   delete[] fProcPnt;
 116   delete[] fVecBLoc;
 117   delete[] fDiagCGlo;
 118   delete[] fVecBGlo;
 119   delete[] fInitPar;
 120   delete[] fDeltaPar;
 121   delete[] fSigmaPar;
 122   delete[] fGlo2CGlo;
 123   delete[] fCGlo2Glo;
 124   delete[] fIsLinear;
 125   delete[] fConstrUsed;
 126   delete[] fFillIndex;
 127   delete[] fFillValue;
 128   //
 129   delete fRecord;
 130   delete fMatCLoc;
 131   delete fMatCGlo;
 132   delete fMatCGloLoc;
 133 }
 134
 135 //_____________________________________________________________________________________________
 136 Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit)
 137 {
 138   //
 139   if (nLoc>0)        fNLocPar = nLoc;
 140   if (nGlo>0)        fNGloPar = nGlo;
 141   if (lResCutInit>0) fResCutInit = lResCutInit;
 142   if (lResCut>0)     fResCut     = lResCut;
 143   if (lNStdDev>0)    fNStdDev    = lNStdDev;
 144   //
 145   fNGloSize = fNGloPar;
 146   //
 147   try {
 148     //
 149     if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);}
 150     else                   fMatCGlo = new AliSymMatrix(fNGloPar);
 151     //
 152     fFillIndex    = new Int_t[fNGloPar];
 153     fFillValue    = new Double_t[fNGloPar];
 154     //
 155     fMatCLoc      = new AliSymMatrix(fNLocPar);
 156     fMatCGloLoc   = new AliRectMatrix(fNGloPar,fNLocPar);
 157     //
 158     fParamGrID    = new Int_t[fNGloPar];
 159     fProcPnt      = new Int_t[fNGloPar];
 160     fVecBLoc      = new Double_t[fNLocPar];
 161     fDiagCGlo     = new Double_t[fNGloPar];
 162     //
 163     fInitPar      = new Double_t[fNGloPar];
 164     fDeltaPar     = new Double_t[fNGloPar];
 165     fSigmaPar     = new Double_t[fNGloPar];
 166     fIsLinear     = new Bool_t[fNGloPar];
 167     //
 168     fGlo2CGlo     = new Int_t[fNGloPar];
 169     fCGlo2Glo     = new Int_t[fNGloPar];
 170   }
 171   catch(bad_alloc&) {
 172     AliInfo(Form("Failed to allocate the memory for %d global and %d local parameters",fNGloPar,fNLocPar));
 173     return 0;
 174   }
 175   //
 176   memset(fVecBLoc   ,0,fNLocPar*sizeof(Double_t));
 177   memset(fDiagCGlo  ,0,fNGloPar*sizeof(Double_t));
 178   memset(fInitPar   ,0,fNGloPar*sizeof(Double_t));
 179   memset(fDeltaPar  ,0,fNGloPar*sizeof(Double_t));
 180   memset(fSigmaPar  ,0,fNGloPar*sizeof(Double_t));
 181   memset(fProcPnt   ,0,fNGloPar*sizeof(Int_t));
 182   //
 183   for (int i=fNGloPar;i--;) {
 184     fGlo2CGlo[i] = fCGlo2Glo[i] = -1;
 185     fIsLinear[i] = kTRUE;
 186     fParamGrID[i] = -1;
 187   }
 188   //
 189   return 1;
 190 }
 191
 192 //_____________________________________________________________________________________________
 193 Bool_t AliMillePede2::ImposeDataRecFile(const char* fname)
 194 {
 195   CloseDataRecStorage();
 196   SetDataRecFName(fname);
 197   return InitDataRecStorage(kTRUE); // open in read mode
 198 }
 199
 200 //_____________________________________________________________________________________________
 201 Bool_t AliMillePede2::ImposeConsRecFile(const char* fname)
 202 {
 203   CloseConsRecStorage();
 204   SetConsRecFName(fname);
 205   return InitConsRecStorage(kTRUE); // open in read mode
 206 }
 207
 208 //_____________________________________________________________________________________________
 209 Bool_t AliMillePede2::InitDataRecStorage(Bool_t read)
 210 {
 211   // initialize the buffer for processed measurements records
 212   //
 213   if (fDataRecFile) {AliInfo("Data Records File is already initialized"); return kFALSE;}
 214   //
 215   if (!fRecord) fRecord = new AliMillePedeRecord();
 216   //
 217   fDataRecFile = TFile::Open(GetDataRecFName(),read ? "":"recreate");
 218   if (!fDataRecFile) {AliInfo(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;}
 219   //
 220   AliInfo(Form("File %s used for derivatives records",GetDataRecFName()));
 221   if (read) {
 222     fTreeData = (TTree*)fDataRecFile->Get("AliMillePedeRecords_Data");
 223     if (!fTreeData) {AliInfo(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;}
 224     fTreeData->SetBranchAddress("Record_Data",&fRecord);
 225     AliInfo(Form("Found %d derivatives records",fTreeData->GetEntries()));
 226   }
 227   else {
 228     fTreeData = new TTree("AliMillePedeRecords_Data","Data Records for AliMillePede2");
 229     fTreeData->Branch("Record_Data","AliMillePedeRecord",&fRecord,32000,99);
 230   }
 231   fCurrRecDataID = -1;
 232   //
 233   return kTRUE;
 234 }
 235
 236 //_____________________________________________________________________________________________
 237 Bool_t AliMillePede2::InitConsRecStorage(Bool_t read)
 238 {
 239   // initialize the buffer for processed measurements records
 240   //
 241   if (fConsRecFile) {AliInfo("Constraints Records File is already initialized"); return kFALSE;}
 242   //
 243   if (!fRecord) fRecord = new AliMillePedeRecord();
 244   //
 245   fConsRecFile = TFile::Open(GetConsRecFName(),read ? "":"recreate");
 246   if (!fConsRecFile) {AliInfo(Form("Failed to initialize constraints records file %s",GetConsRecFName())); return kFALSE;}
 247   //
 248   AliInfo(Form("File %s used for constraints records",GetConsRecFName()));
 249   if (read) {
 250     fTreeConstr = (TTree*)fConsRecFile->Get("AliMillePedeRecords_Constraints");
 251     if (!fTreeConstr) {AliInfo(Form("Did not find constraints records tree in %s",GetConsRecFName())); return kFALSE;}
 252     fTreeConstr->SetBranchAddress("Record_Constraints",&fRecord);
 253     AliInfo(Form("Found %d constraints records",fTreeConstr->GetEntries()));
 254     //
 255   }
 256   else {
 257     //
 258     fTreeConstr = new TTree("AliMillePedeRecords_Constraints","Constraints Records for AliMillePede2");
 259     fTreeConstr->Branch("Record_Constraints","AliMillePedeRecord",&fRecord,32000,99);
 260   }
 261   fCurrRecConstrID = -1;
 262   //
 263   return kTRUE;
 264 }
 265
 266 //_____________________________________________________________________________________________
 267 void AliMillePede2::CloseDataRecStorage()
 268 {
 269   if (fTreeData) {
 270     if (fDataRecFile->IsWritable()) {
 271       fDataRecFile->cd();
 272       fTreeData->Write();
 273     }
 274     delete fTreeData;
 275     fTreeData = 0;
 276     fDataRecFile->Close();
 277     delete fDataRecFile;
 278     fDataRecFile = 0;
 279   }
 280   //
 281 }
 282
 283 //_____________________________________________________________________________________________
 284 void AliMillePede2::CloseConsRecStorage()
 285 {
 286   if (fTreeConstr) {
 287     if (fConsRecFile->IsWritable()) {
 288       fConsRecFile->cd();
 289       fTreeConstr->Write();
 290     }
 291     delete fTreeConstr;
 292     fTreeConstr = 0;
 293     fConsRecFile->Close();
 294     delete fConsRecFile;
 295     fConsRecFile = 0;
 296   }
 297   //
 298 }
 299
 300 //_____________________________________________________________________________________________
 301 Bool_t AliMillePede2::ReadNextRecordData()
 302 {
 303   // read next data record (if any)
 304   if (!fTreeData || ++fCurrRecDataID >= fTreeData->GetEntries()) { fCurrRecDataID--; return kFALSE;}
 305   fTreeData->GetEntry(fCurrRecDataID);
 306   return kTRUE;
 307 }
 308
 309 //_____________________________________________________________________________________________
 310 Bool_t AliMillePede2::ReadNextRecordConstraint()
 311 {
 312   // read next constraint record (if any)
 313   if (!fTreeConstr || ++fCurrRecConstrID >= fTreeConstr->GetEntries()) { fCurrRecConstrID--; return kFALSE;}
 314   fTreeConstr->GetEntry(fCurrRecConstrID);
 315   return kTRUE;
 316 }
 317
 318 //_____________________________________________________________________________________________
 319 void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma)
 320 {
 321   if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
 322   //
 323   // write data of single measurement
 324   if (lSigma<=0.0) { // If parameter is fixed, then no equation
 325     for (int i=fNLocPar; i--;) derlc[i] = 0.0;
 326     for (int i=fNGloPar; i--;) dergb[i] = 0.0;
 327     return;
 328   }
 329   //
 330   fRecord->AddResidual(lMeas);
 331   //
 332   // Retrieve local param interesting indices
 333   for (int i=0;i<fNLocPar;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(i,derlc[i]); derlc[i] = 0.0;}
 334   //
 335   fRecord->AddWeight( 1.0/lSigma/lSigma );
 336   //
 337   // Idem for global parameters
 338   for (int i=0;i<fNGloPar;i++) if (!IsZero(dergb[i])) {
 339     fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0;
 340     fRecord->MarkGroup(fParamGrID[i]);
 341   }
 342   //
 343 }
 344
 345 //_____________________________________________________________________________________________
 346 void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc,
 347                                      double *derlc,int nlc,double lMeas,double lSigma)
 348 {
 349   // write data of single measurement
 350   if (lSigma<=0.0) { // If parameter is fixed, then no equation
 351     for (int i=nlc;i--;)  derlc[i] = 0.0;
 352     for (int i=ngb;i--;)  dergb[i] = 0.0;
 353     return;
 354   }
 355   //
 356   if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
 357   //
 358   fRecord->AddResidual(lMeas);
 359   //
 360   // Retrieve local param interesting indices
 361   for (int i=0;i<nlc;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;}
 362   //
 363   fRecord->AddWeight( 1./lSigma/lSigma );
 364   //
 365   // Idem for global parameters
 366   for (int i=0;i<ngb;i++) if (!IsZero(dergb[i])) {fRecord->AddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;}
 367   //
 368 }
 369
 370
 371 //_____________________________________________________________________________________________
 372 void AliMillePede2::SetGlobalConstraint(double *dergb, double val, double sigma)
 373 {
 374   // Define a constraint equation.
 375   if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
 376   //
 377   fRecord->Reset();
 378   fRecord->AddResidual(val);
 379   fRecord->AddWeight(sigma);
 380   for (int i=0; i<fNGloPar; i++) if (!IsZero(dergb[i]))  fRecord->AddIndexValue(i,dergb[i]);
 381   fNGloConstraints++;
 382   if (IsZero(sigma)) fNLagrangeConstraints++;
 383   //  printf("NewConstraint:\n"); fRecord->Print(); //RRR
 384   SaveRecordConstraint();
 385   //
 386 }
 387
 388 //_____________________________________________________________________________________________
 389 void AliMillePede2::SetGlobalConstraint(const int *indgb, double *dergb, int ngb, double val,double sigma)
 390 {
 391   // Define a constraint equation.
 392   if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
 393   fRecord->Reset();
 394   fRecord->AddResidual(val);
 395   fRecord->AddWeight(sigma);   // dummy
 396   for (int i=0; i<ngb; i++) if (!IsZero(dergb[i]))  fRecord->AddIndexValue(indgb[i],dergb[i]);
 397   fNGloConstraints++;
 398   if (IsZero(sigma)) fNLagrangeConstraints++;
 399   SaveRecordConstraint();
 400   //
 401 }
 402
 403 //_____________________________________________________________________________________________
 404 Int_t AliMillePede2::LocalFit(double *localParams)
 405 {
 406   /*
 407     Perform local parameters fit once all the local equations have been set
 408     -----------------------------------------------------------
 409     localParams = (if !=0) will contain the fitted track parameters and
 410     related errors
 411   */
 412   static int     nrefSize = 0;
 413   //  static TArrayI refLoc,refGlo,nrefLoc,nrefGlo;
 414   static Int_t  *refLoc=0,*refGlo=0,*nrefLoc=0,*nrefGlo=0;
 415   int nPoints = 0;
 416   //
 417   AliSymMatrix    &matCLoc    = *fMatCLoc;
 418   AliMatrixSq     &matCGlo    = *fMatCGlo;
 419   AliRectMatrix   &matCGloLoc = *fMatCGloLoc;
 420   //
 421   memset(fVecBLoc,0,fNLocPar*sizeof(double));
 422   matCLoc.Reset();
 423   //
 424   int cnt=0;
 425   int recSz = fRecord->GetSize();
 426   //
 427   while(cnt<recSz) {  // Transfer the measurement records to matrices
 428     //
 429     // extract addresses of residual, weight and pointers on local and global derivatives for each point
 430     if (nrefSize<=nPoints) {
 431       int *tmpA = 0;
 432       nrefSize = 2*(nPoints+1);
 433       tmpA = refLoc;  refLoc  = new Int_t[nrefSize]; if (tmpA) memcpy(refLoc,tmpA,nPoints*sizeof(int));
 434       tmpA = refGlo;  refGlo  = new Int_t[nrefSize]; if (tmpA) memcpy(refGlo,tmpA,nPoints*sizeof(int));
 435       tmpA = nrefLoc; nrefLoc = new Int_t[nrefSize]; if (tmpA) memcpy(nrefLoc,tmpA,nPoints*sizeof(int));
 436       tmpA = nrefGlo; nrefGlo = new Int_t[nrefSize]; if (tmpA) memcpy(nrefGlo,tmpA,nPoints*sizeof(int));
 437     }
 438     //
 439     refLoc[nPoints]  = ++cnt;
 440     int nLoc = 0;
 441     while(!fRecord->IsWeight(cnt)) {nLoc++; cnt++;}
 442     nrefLoc[nPoints] = nLoc;
 443     //
 444     refGlo[nPoints]  = ++cnt;
 445     int nGlo = 0;
 446     while(!fRecord->IsResidual(cnt) && cnt<recSz) {nGlo++; cnt++;}
 447     nrefGlo[nPoints] = nGlo;
 448     //
 449     nPoints++;
 450   }
 451   double vl;
 452   //
 453   Int_t maxLocUsed = 0;
 454   //
 455   for (int ip=nPoints;ip--;) {  // Transfer the measurement records to matrices
 456     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
 457     double  weight = fRecord->GetValue( refGlo[ip]-1 );
 458     double *derLoc = fRecord->GetValue()+refLoc[ip];
 459     double *derGlo = fRecord->GetValue()+refGlo[ip];
 460     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
 461     int    *indGlo = fRecord->GetIndex()+refGlo[ip];
 462     //
 463     for (int i=nrefGlo[ip];i--;) {       // suppress the global part (only relevant with iterations)
 464       int iID = indGlo[i];              // Global param indice
 465       if (fSigmaPar[iID]<=0.) continue;                                    // fixed parameter RRRCheck
 466       if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);  // linear parameter
 467       else                resid -= derGlo[i]*fDeltaPar[iID];                  // nonlinear parameter
 468     }
 469     //
 470     // Symmetric matrix, don't bother j>i coeffs
 471     for (int i=nrefLoc[ip];i--;) {         // Fill local matrix and vector
 472       fVecBLoc[ indLoc[i] ] += weight*resid*derLoc[i];
 473       if (indLoc[i]>maxLocUsed) maxLocUsed = indLoc[i];
 474       for (int j=i+1;j--;) matCLoc(indLoc[i] ,indLoc[j]) += weight*derLoc[i]*derLoc[j];
 475     }
 476     //
 477   } // end of the transfer of the measurement record to matrices
 478   //
 479   matCLoc.SetSizeUsed(++maxLocUsed);   // data with B=0 may use less than declared nLocals
 480   //
 481   // first try to solve by faster Cholesky decomposition, then by Gaussian elimination
 482   if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) {
 483     AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
 484     if (!matCLoc.SolveSpmInv(fVecBLoc,kTRUE)) {
 485       AliInfo("Failed to solve locals by Gaussian Elimination, skip...");
 486       matCLoc.Print("d");
 487       return 0; // failed to solve
 488     }
 489   }
 490   //
 491   // If requested, store the track params and errors
 492   //  printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
 493   if (localParams) for (int i=maxLocUsed; i--;) {
 494       localParams[2*i]   = fVecBLoc[i];
 495       localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i)));
 496     }
 497   //
 498   float lChi2 = 0;
 499   int   nEq   = 0;
 500   //
 501   for (int ip=nPoints;ip--;) {  // Calculate residuals
 502     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
 503     double  weight = fRecord->GetValue( refGlo[ip]-1 );
 504     double *derLoc = fRecord->GetValue()+refLoc[ip];
 505     double *derGlo = fRecord->GetValue()+refGlo[ip];
 506     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
 507     int    *indGlo = fRecord->GetIndex()+refGlo[ip];
 508     //
 509     // Suppress local and global contribution in residuals;
 510     for (int i=nrefLoc[ip];i--;) resid -= derLoc[i]*fVecBLoc[ indLoc[i] ];     // local part
 511     //
 512     for (int i=nrefGlo[ip];i--;) {                                             // global part
 513       int iID = indGlo[i];
 514       if ( fSigmaPar[iID] <= 0.) continue;                                     // fixed parameter RRRCheck
 515       if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);   // linear parameter
 516       else                resid -= derGlo[i]*fDeltaPar[iID];                   // nonlinear parameter
 517     }
 518     //
 519     // reject the track if the residual is too large (outlier)
 520     double absres = TMath::Abs(resid);
 521     if ( (absres >= fResCutInit && fIter ==1 ) ||
 522          (absres >= fResCut     && fIter > 1)) {
 523       if (fLocFitAdd) fNLocFitsRejected++;
 524       //      printf("reject res %5ld %+e\n",fCurrRecDataID,resid);
 525       return 0;
 526     }
 527     //
 528     lChi2 += weight*resid*resid ; // total chi^2
 529     nEq++;                        // number of equations
 530   } // end of Calculate residuals
 531   //
 532   int nDoF = nEq-maxLocUsed;
 533   lChi2 = (nDoF>0) ? lChi2/nDoF : 0;  // Chi^2/dof
 534   //
 535   if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2
 536     if (fLocFitAdd) fNLocFitsRejected++;
 537     //    printf("reject chi2 %5ld: %+e\n",fCurrRecDataID, lChi2);
 538     return 0;
 539   }
 540   //
 541   if (fLocFitAdd) {
 542     fNLocFits++;
 543     fNLocEquations += nEq;
 544   }
 545   else {
 546     fNLocFits--;
 547     fNLocEquations -= nEq;
 548   }
 549   //
 550   //  local operations are finished, track is accepted
 551   //  We now update the global parameters (other matrices)
 552   //
 553   int nGloInFit = 0;
 554   //
 555   for (int ip=nPoints;ip--;) {  // Update matrices
 556     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
 557     double  weight = fRecord->GetValue( refGlo[ip]-1 );
 558     double *derLoc = fRecord->GetValue()+refLoc[ip];
 559     double *derGlo = fRecord->GetValue()+refGlo[ip];
 560     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
 561     int    *indGlo = fRecord->GetIndex()+refGlo[ip];
 562     //
 563     for (int i=nrefGlo[ip];i--;) {    // suppress the global part
 564       int iID = indGlo[i];           // Global param indice
 565       if ( fSigmaPar[iID] <= 0.) continue;                                         // fixed parameter RRRCheck
 566       if (fIsLinear[iID])  resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);      // linear parameter
 567       else                 resid -= derGlo[i]*fDeltaPar[iID];                      // nonlinear parameter
 568     }
 569     //
 570     for (int ig=nrefGlo[ip];ig--;) {
 571       int iIDg = indGlo[ig];   // Global param indice (the matrix line)
 572       if ( fSigmaPar[iIDg] <= 0.) continue;                                    // fixed parameter RRRCheck
 573       if (fLocFitAdd) fVecBGlo[ iIDg ] += weight*resid*derGlo[ig]; //!!!
 574       else            fVecBGlo[ iIDg ] -= weight*resid*derGlo[ig]; //!!!
 575       //
 576       // First of all, the global/global terms (exactly like local matrix)
 577       int nfill = 0;
 578       for (int jg=ig+1;jg--;) { // matCGlo is symmetric by construction
 579         int jIDg = indGlo[jg];
 580         if ( fSigmaPar[jIDg] <= 0.) continue;                                    // fixed parameter RRRCheck
 581         if ( !IsZero(vl = weight*derGlo[ig]*derGlo[jg]) ) {
 582           fFillIndex[nfill]   = jIDg;
 583           fFillValue[nfill++] = fLocFitAdd ? vl:-vl;
 584         }
 585       }
 586       if (nfill) matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill);
 587       //
 588       // Now we have also rectangular matrices containing global/local terms.
 589       int iCIDg = fGlo2CGlo[iIDg];  // compressed Index of index
 590       if (iCIDg == -1) {
 591         Double_t *rowGL = matCGloLoc(nGloInFit);
 592         for (int k=maxLocUsed;k--;) rowGL[k] = 0.0;  // reset the row
 593         iCIDg = fGlo2CGlo[iIDg] = nGloInFit;
 594         fCGlo2Glo[nGloInFit++] = iIDg;
 595       }
 596       //
 597       Double_t *rowGLIDg = matCGloLoc(iCIDg);
 598       for (int il=nrefLoc[ip];il--;) rowGLIDg[ indLoc[il] ] += weight*derGlo[ig]*derLoc[il];
 599       fProcPnt[iIDg] += fLocFitAdd ? 1:-1;       // update counter
 600       //
 601     }
 602   } // end of Update matrices
 603   //
 604   // calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T
 605   // and       fVecBGlo -= fMatCGloLoc * fVecBLoc
 606   //
 607   //-------------------------------------------------------------- >>>
 608   double vll;
 609   for (int iCIDg=0; iCIDg<nGloInFit; iCIDg++) {
 610     int iIDg = fCGlo2Glo[iCIDg];
 611     //
 612     vl = 0;
 613     Double_t *rowGLIDg =  matCGloLoc(iCIDg);
 614     for (int kl=0;kl<maxLocUsed;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
 615     if  (!IsZero(vl)) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl;
 616     //
 617     int nfill = 0;
 618     for (int jCIDg=0;jCIDg<=iCIDg; jCIDg++) {
 619       int jIDg = fCGlo2Glo[jCIDg];
 620       //
 621       vl = 0;
 622       Double_t *rowGLJDg =  matCGloLoc(jCIDg);
 623       for (int kl=0;kl<maxLocUsed;kl++) {
 624         // diag terms
 625         if ( (!IsZero(vll=rowGLIDg[kl]*rowGLJDg[kl])) ) vl += matCLoc.QueryDiag(kl)*vll;
 626         //
 627         // off-diag terms
 628         for (int ll=0;ll<kl;ll++) {
 629           if ( !IsZero(vll=rowGLIDg[kl]*rowGLJDg[ll]) ) vl += matCLoc(kl,ll)*vll;
 630           if ( !IsZero(vll=rowGLIDg[ll]*rowGLJDg[kl]) ) vl += matCLoc(kl,ll)*vll;
 631         }
 632       }
 633       if (!IsZero(vl)) {
 634         fFillIndex[nfill]   = jIDg;
 635         fFillValue[nfill++] = fLocFitAdd ? -vl : vl;
 636       }
 637     }
 638     if (nfill)  matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill);
 639   }
 640   //
 641   // reset compressed index array
 642   //
 643   for (int i=nGloInFit;i--;) {
 644     fGlo2CGlo[ fCGlo2Glo[i] ] = -1;
 645     fCGlo2Glo[i] = -1;
 646   }
 647   //
 648   //---------------------------------------------------- <<<
 649   return 1;
 650 }
 651
 652 //_____________________________________________________________________________________________
 653 Int_t AliMillePede2::GlobalFit(Double_t *par, Double_t *error, Double_t *pull)
 654 {
 655   // performs a requested number of global iterations
 656   fIter = 1;
 657   //
 658   TStopwatch sw; sw.Start();
 659   //
 660   int res = 0;
 661   AliInfo("Starting Global fit.");
 662   while (fIter<=fMaxIter) {
 663     //
 664     res = GlobalFitIteration();
 665     if (!res) break;
 666     //
 667     if (!IsZero(fChi2CutFactor-fChi2CutRef)) {
 668       fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
 669       if (fChi2CutFactor < 1.2*fChi2CutRef) {
 670         fChi2CutFactor = fChi2CutRef;
 671         fIter = fMaxIter - 1;     // Last iteration
 672       }
 673     }
 674     fIter++;
 675   }
 676   //
 677   sw.Stop();
 678   AliInfo(Form("Global fit %s, CPU time: %.1f",res ? "Converged":"Failed",sw.CpuTime()));
 679   if (!res) return 0;
 680   //
 681   if (par) for (int i=fNGloPar;i--;) par[i] = fInitPar[i]+fDeltaPar[i];
 682   //
 683   if (fGloSolveStatus==gkInvert) { // errors on params are available
 684     if (error) for (int i=fNGloPar;i--;) error[i] = fProcPnt[i]>0 ? TMath::Sqrt(TMath::Abs(fMatCGlo->QueryDiag(i))) : 0.;
 685     if (pull)  for (int i=fNGloPar;i--;) pull[i]  = fProcPnt[i]>0 && (fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i))>0. && fSigmaPar[i]>0
 686                                            ? fDeltaPar[i]/TMath::Sqrt(fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i)) : 0;
 687   }
 688   //
 689   return 1;
 690 }
 691
 692 //_____________________________________________________________________________________________
 693 Int_t AliMillePede2::GlobalFitIteration()
 694 {
 695   // perform global parameters fit once all the local equations have been fitted
 696   //
 697   AliInfo(Form("Global Fit Iteration#%2d (Local Fit Chi^2 cut factor: %.2f)",fIter,fChi2CutFactor));
 698   //
 699   if (!fNGloPar || !fTreeData) {
 700     AliInfo("No data was stored, aborting iteration");
 701     return 0;
 702   }
 703   TStopwatch sw,sws;
 704   sw.Start();
 705   sws.Stop();
 706   //
 707   if (!fConstrUsed) {
 708     fConstrUsed = new Bool_t[fNGloConstraints];
 709     memset(fConstrUsed,0,fNGloConstraints*sizeof(Bool_t));
 710   }
 711   // Reset all info specific for this step
 712   AliMatrixSq& matCGlo = *fMatCGlo;
 713   matCGlo.Reset();
 714   memset(fProcPnt,0,fNGloPar*sizeof(Int_t));
 715   //
 716   fNGloConstraints = fTreeConstr ? fTreeConstr->GetEntries() : 0;
 717   //
 718   // count number of Lagrange constraints: they need new row/cols to be added
 719   fNLagrangeConstraints = 0;
 720   for (int i=0; i<fNGloConstraints; i++) {
 721     ReadRecordConstraint(i);
 722     if ( IsZero(fRecord->GetValue(1)) ) fNLagrangeConstraints++; // exact constraint (no error) -> Lagrange multiplier
 723   }
 724   //
 725   // if needed, readjust the size of the global vector (for matrices this is done automatically)
 726   if (!fVecBGlo || fNGloSize!=fNGloPar+fNLagrangeConstraints) {
 727     delete[] fVecBGlo;   // in case some constraint was added between the two manual iterations
 728     fNGloSize = fNGloPar+fNLagrangeConstraints;
 729     fVecBGlo = new Double_t[fNGloSize];
 730   }
 731   memset(fVecBGlo,0,fNGloSize*sizeof(double));
 732   //
 733   fNLocFits         = 0;
 734   fNLocFitsRejected = 0;
 735   fNLocEquations    = 0;
 736   //
 737   //  Process data records and build the matrices
 738   Long_t ndr = fTreeData->GetEntries();
 739   AliInfo(Form("Building the Global matrix from %d data records",ndr));
 740   if (!ndr) return 0;
 741   //
 742   TStopwatch swt; swt.Start();
 743   fLocFitAdd = kTRUE;  // add contributions of matching tracks
 744   for (Long_t i=0;i<ndr;i++) {
 745     ReadRecordData(i);
 746     LocalFit();
 747     if ( (i%int(0.2*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
 748   }
 749   swt.Stop();
 750   printf("%ld local fits done: ", ndr);
 751   swt.Print();
 752   sw.Start(kFALSE);
 753   //
 754   //
 755   // ---------------------- Reject parameters with low statistics ------------>>
 756   fNGloFix = 0;
 757   if (fMinPntValid>1 && fNGroupsSet) {
 758     //
 759     printf("Checking parameters with statistics < %d\n",fMinPntValid);
 760     TStopwatch swsup;
 761     swsup.Start();
 762     // 1) build the list of parameters to fix
 763     Int_t fixArrSize = 10;
 764     Int_t nFixedGroups = 0;
 765     TArrayI fixGroups(fixArrSize);
 766     //
 767     for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
 768       int grID = fParamGrID[i];
 769       if (grID<0) continue; // not in the group
 770       if (fProcPnt[i]>=fMinPntValid) continue;
 771       fProcPnt[i] = 0;
 772       // check if the group is already accounted
 773       Bool_t fnd = kFALSE;
 774       for (int j=nFixedGroups;j--;) if (fixGroups[j]==grID) {fnd=kTRUE; break;}
 775       if (fnd) continue;
 776       if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);}
 777       fixGroups[nFixedGroups++] = grID; // add group to fix
 778       //
 779     }
 780     // 2) loop over records and add contributions of fixed groups with negative sign
 781     fLocFitAdd = kFALSE;
 782     //
 783     for (Long_t i=0;i<ndr;i++) {
 784       ReadRecordData(i);
 785       Bool_t suppr = kFALSE;
 786       for (int ifx=nFixedGroups;ifx--;)if (fRecord->IsGroupPresent(fixGroups[ifx])) suppr = kTRUE;
 787       if (suppr) LocalFit();
 788     }
 789     fLocFitAdd = kTRUE;
 790     //
 791     if (nFixedGroups) {
 792       printf("Suppressed contributions of groups with NPoints<%d :\n",fMinPntValid);
 793       for (int i=0;i<nFixedGroups;i++) printf("%d ",fixGroups[i]); printf("\n");
 794     }
 795     swsup.Stop();
 796     swsup.Print();
 797   }
 798   // ---------------------- Reject parameters with low statistics ------------<<
 799   //
 800   // add large number to diagonal of fixed params
 801   //
 802   for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
 803     //    printf("#%3d : Nproc : %5d   grp: %d\n",i,fProcPnt[i],fParamGrID[i]);
 804     if (fProcPnt[i]<1) {
 805       fNGloFix++;
 806       fVecBGlo[i] = 0.;
 807       matCGlo.DiagElem(i) = float(fNLocEquations*fNLocEquations);
 808     }
 809     else matCGlo.DiagElem(i) += (fgWeightSigma ? fProcPnt[i] : 1.)/(fSigmaPar[i]*fSigmaPar[i]);
 810   }
 811   //
 812   for (int i=fNGloPar;i--;) fDiagCGlo[i] = matCGlo.QueryDiag(i); // save the diagonal elements
 813   //
 814   // add constraint equations
 815   int nVar = fNGloPar;                    // Current size of global matrix
 816   for (int i=0; i<fNGloConstraints; i++) {
 817     ReadRecordConstraint(i);
 818     double val   = fRecord->GetValue(0);
 819     double sig   = fRecord->GetValue(1);
 820     int    *indV = fRecord->GetIndex()+2;
 821     double *der  = fRecord->GetValue()+2;
 822     int    csize = fRecord->GetSize()-2;
 823     //
 824     // check if after suppression of fixed variables there are non-0 derivatives
 825     // and determine the max statistics of involved params
 826     int nSuppressed = 0;
 827     int maxStat = 1;
 828     for (int j=csize;j--;) {
 829       if (fProcPnt[indV[j]]<1) nSuppressed++;
 830       else {
 831         maxStat = TMath::Max(maxStat,fProcPnt[indV[j]]);
 832       }
 833     }
 834     //
 835     if (nSuppressed==csize) {
 836       //      AliInfo(Form("Neglecting constraint %d of %d derivatives since no free parameters left",i,csize));
 837       //
 838       // was this constraint ever created ?
 839       if ( sig==0 && fConstrUsed[i] ) { // this is needed only for constraints with Lagrange multiplier
 840         // to avoid empty row impose dummy constraint on "Lagrange multiplier"
 841         matCGlo.DiagElem(nVar) = 1.;
 842         fVecBGlo[nVar++] = 0;
 843       }
 844       continue;
 845     }
 846     //
 847     // account for already accumulated corrections
 848     for (int j=csize; j--;) val -= der[j]*(fInitPar[ indV[j] ]+fDeltaPar[ indV[j] ]);
 849     //
 850     if (sig > 0) {  // this is a gaussian constriant: no Lagrange multipliers are added
 851       //
 852       double sig2i = (fgWeightSigma ? TMath::Sqrt(maxStat) : 1.)/sig/sig;
 853       for (int ir=0;ir<csize;ir++) {
 854         int iID = indV[ir];
 855         for (int ic=0;ic<=ir;ic++) { // matrix is symmetric
 856           int jID = indV[ic];
 857           double vl = der[ir]*der[ic]*sig2i;
 858           if (!IsZero(vl)) matCGlo(iID,jID) += vl;
 859         }
 860         fVecBGlo[iID] += val*der[ir]*sig2i;
 861       }
 862     }
 863     else {   // this is exact constriant:  Lagrange multipliers must be added
 864       for (int j=csize; j--;) {
 865         int jID = indV[j];
 866         if (fProcPnt[jID]<1) continue;                      // this parameter was fixed, don't put it into constraint
 867         matCGlo(nVar,jID) = float(fNLocEquations)*der[j];   // fMatCGlo is symmetric, only lower triangle is filled
 868       }
 869       //
 870       if (matCGlo.QueryDiag(nVar)) matCGlo.DiagElem(nVar) = 0.0;
 871       fVecBGlo[nVar++] = float(fNLocEquations)*val; //RS ? should we use here fNLocFits ?
 872       fConstrUsed[i] = kTRUE;
 873     }
 874   }
 875   //
 876   AliInfo(Form("Obtained %-7ld equations from %-7ld records (%-7ld rejected). Fixed %-4d globals",
 877                fNLocEquations,fNLocFits,fNLocFitsRejected,fNGloFix));
 878
 879   //
 880   sws.Start();
 881   fGloSolveStatus = SolveGlobalMatEq();                     // obtain solution for this step
 882   sws.Stop();
 883   printf("Solve %d |",fIter); sws.Print();
 884   //
 885   sw.Stop();
 886   AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==gkFailed ? "Failed":"Converged",sw.CpuTime()));
 887   if (fGloSolveStatus==gkFailed) return 0;
 888   //
 889   for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
 890   //
 891   //  PrintGlobalParameters();
 892   return 1;
 893 }
 894
 895 //_____________________________________________________________________________________________
 896 Int_t AliMillePede2::SolveGlobalMatEq()
 897 {
 898   //
 899   // solve global matrix equation MatCGlob*X=VecBGlo and store the result in the VecBGlo
 900   //
 901   /*
 902   printf("GlobalMatrix\n");
 903   fMatCGlo->Print();
 904   printf("RHS\n");
 905   for (int i=0;i<fNGloPar;i++) printf("%d %+e\n",i,fVecBGlo[i]);
 906   */
 907   //
 908   if (!fgIsMatGloSparse) {
 909     //
 910     if (fNLagrangeConstraints==0) { // pos-def systems are faster to solve by Cholesky
 911       if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, fgInvChol) ) return fgInvChol ? gkInvert:gkNoInversion;
 912       else AliInfo("Solution of Global Dense System by Cholesky failed, trying Gaussian Elimiation");
 913     }
 914     //
 915     if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return gkInvert;
 916     else AliInfo("Solution of Global Dense System by Gaussian Elimiation failed, trying iterative methods");
 917   }
 918   // try to solve by minres
 919   TVectorD sol(fNGloSize);
 920   //
 921   AliMinResSolve *slv = new AliMinResSolve(fMatCGlo,fVecBGlo);
 922   if (!slv) return gkFailed;
 923   //
 924   Bool_t res = kFALSE;
 925   if      (fgIterSol == AliMinResSolve::kSolMinRes)
 926     res =  slv->SolveMinRes(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol);
 927   else if (fgIterSol == AliMinResSolve::kSolFGMRes)
 928     res =  slv->SolveFGMRES(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV);
 929   else
 930     AliInfo(Form("Undefined Iteritive Solver ID=%d, only %d are defined",fgIterSol,AliMinResSolve::kNSolvers));
 931   //
 932   if (!res) return gkFailed;
 933   for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i];
 934   return gkNoInversion;
 935   //
 936 }
 937
 938 //_____________________________________________________________________________________________
 939 Float_t AliMillePede2::Chi2DoFLim(int nSig, int nDoF) const
 940 {
 941   /// return the limit in chi^2/nd for n sigmas stdev authorized
 942   // Only n=1, 2, and 3 are expected in input
 943   int lNSig;
 944   float sn[3]        =  {0.47523, 1.690140, 2.782170};
 945   float table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
 946                          1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
 947                          1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
 948                          1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
 949                         {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
 950                          1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
 951                          1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
 952                          1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
 953                         {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
 954                          2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
 955                          2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
 956                          1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
 957
 958   if (nDoF < 1) {
 959     return 0.0;
 960   }
 961   else {
 962     lNSig = TMath::Max(1,TMath::Min(nSig,3));
 963
 964     if (nDoF <= 30) {
 965       return table[lNSig-1][nDoF-1];
 966     }
 967     else { // approximation
 968       return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
 969               (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);
 970     }
 971   }
 972 }
 973
 974 //_____________________________________________________________________________________________
 975 Int_t AliMillePede2::SetIterations(double lChi2CutFac)
 976 {
 977   // Number of iterations is calculated from lChi2CutFac
 978   fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
 979   //
 980   AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
 981   fIter = 1; // Initializes the iteration process
 982   return 1;
 983 }
 984
 985 //_____________________________________________________________________________________________
 986 Double_t AliMillePede2::GetParError(int iPar) const
 987 {
 988   // return error for parameter iPar
 989   if (fGloSolveStatus==gkInvert) {
 990     double res = fMatCGlo->QueryDiag(iPar);
 991     if (res>=0) return TMath::Sqrt(res);
 992   }
 993   return -1.;
 994 }
 995
 996
 997 //_____________________________________________________________________________________________
 998 Int_t AliMillePede2::PrintGlobalParameters() const
 999 {
1000   ///  Print the final results into the logfile
1001   double lError = 0.;
1002   double lGlobalCor =0.;
1003
1004   AliInfo("");
1005   AliInfo("   Result of fit for global parameters");
1006   AliInfo("   ===================================");
1007   AliInfo("    I       initial       final       differ        lastcor        error       gcor       Npnt");
1008   AliInfo("----------------------------------------------------------------------------------------------");
1009   //
1010   for (int i=0; i<fNGloPar; i++) {
1011     lError = GetParError(i);
1012     lGlobalCor = 0.0;
1013     //
1014     double dg;
1015     if (fGloSolveStatus==gkInvert && TMath::Abs( (dg=fMatCGlo->QueryDiag(i)) *fDiagCGlo[i]) > 0) {
1016       lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(dg*fDiagCGlo[i])));
1017       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %6d",
1018                    i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i]));
1019     }
1020     else {
1021       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t OFF\t OFF\t %6d",i,fInitPar[i],fInitPar[i]+fDeltaPar[i],
1022                    fDeltaPar[i],fVecBGlo[i],fProcPnt[i]));
1023     }
1024   }
1025   return 1;
1026 }