STEER/AliMillePede2.cxx

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