PWGJE/macros/PID/AliTPCPIDmathFit.cxx

   1 #include "TH1.h"
   2 #include "TMath.h"
   3 #include "TMatrixDSym.h"
   4 #include "TMinuit.h"
   5 #include "TStopwatch.h"
   6 #include "TString.h"
   7
   8 #include "AliLog.h"
   9 #include <iostream>
  10 #include <limits>
  11
  12 #include "AliTPCPIDmathFit.h"
  13
  14 // Class for advanced fitting methods (template fits, (weighted) LL fits, regularised fits, simultaneous fits)
  15 //
  16 // Based on original code of:
  17 // Xianguo Lu
  18 // lu@physi.uni-heidelberg.de
  19 //
  20 // Modified by:
  21 // Benjamin Hess
  22 // Benjamin-Andreas.Hess@Uni-Tuebingen.de
  23
  24
  25 ClassImp(AliTPCPIDmathFit);
  26
  27 AliTPCPIDmathFit* AliTPCPIDmathFit::fgInstance = NULL;
  28
  29 //***Public functions***
  30 //______________________________________________________
  31 Int_t AliTPCPIDmathFit::AddRefHisto(TH1* refHisto)
  32 {
  33   // Add refHisto to the list of reference histos. The index of the added histo is returned.
  34
  35   TH1** newRefHistos = new TH1*[fNrefHistos + 1];
  36   for (Int_t i = 0; i < fNrefHistos; i++)   {
  37     newRefHistos[i] = (TH1*)fRefHistos[i]->Clone();
  38     newRefHistos[i]->SetName(Form("%s_%d", newRefHistos[i]->GetName(), i));
  39     newRefHistos[i]->SetDirectory(0);
  40     delete fRefHistos[i];
  41   }
  42
  43   newRefHistos[fNrefHistos] = (TH1*)refHisto->Clone();
  44
  45   delete [] fRefHistos;
  46   fRefHistos = newRefHistos;
  47
  48   return fNrefHistos++;
  49 }
  50
  51
  52 //______________________________________________________
  53 void AliTPCPIDmathFit::ClearRefHistos()
  54 {
  55   // Clear list of reference histos
  56
  57   for (Int_t i = 0; i < fNrefHistos; i++)
  58     delete fRefHistos[i];
  59
  60   delete [] fRefHistos;
  61   fRefHistos = 0x0;
  62   fNrefHistos = 0;
  63 }
  64
  65
  66 //______________________________________________________
  67 TH1* AliTPCPIDmathFit::GetRefHisto(Int_t index) const
  68 {
  69   if (index < 0 || index >= fNrefHistos || !fRefHistos)
  70     return 0x0;
  71
  72   return fRefHistos[index];
  73 }
  74
  75
  76 //______________________________________________________
  77 Int_t AliTPCPIDmathFit::GetIndexParametersToRegularise(Int_t index) const
  78 {
  79   if (index < 0 || index >= fNumParametersToRegularise || !fIndexParametersToRegularise)
  80     return -1;
  81
  82   return fIndexParametersToRegularise[index];
  83 }
  84
  85
  86 //______________________________________________________
  87 AliTPCPIDmathFit* AliTPCPIDmathFit::Instance(const Int_t numXbinsRegularisation, const Int_t numSimultaneousFits,
  88                                              const Int_t maxDataPoints)
  89 {
  90   // Get instance of AliTPCPIDmathFit. If numSimultaneousFits is > 0 and there is already an existing instance,
  91   // the instance will be deleted and a new one with the new number of simultaneous fits will be created.
  92   // The same applies applies to numXbinsRegularisation.
  93
  94   if(!fgInstance) {
  95     fgInstance = new AliTPCPIDmathFit(numXbinsRegularisation, numSimultaneousFits, maxDataPoints);
  96   }
  97
  98   if (numSimultaneousFits > 0 && numSimultaneousFits != fgInstance->GetNumSimultaneousFits()) {
  99     printf("Different number of simultaneous fits. Creating new instance with %d instead of %d simultaneous fits...\n",
 100            numSimultaneousFits, fgInstance->GetNumSimultaneousFits());
 101     delete fgInstance;
 102     fgInstance = new AliTPCPIDmathFit(numXbinsRegularisation, numSimultaneousFits, maxDataPoints);
 103   }
 104
 105   if (numXbinsRegularisation > 0 && numXbinsRegularisation != fgInstance->GetNumXbinsRegularisation()) {
 106     printf("Different number of x bins for regularisation. Creating new instance with %d instead of %d such bins...\n",
 107            numXbinsRegularisation, fgInstance->GetNumXbinsRegularisation());
 108     delete fgInstance;
 109     fgInstance = new AliTPCPIDmathFit(numXbinsRegularisation, numSimultaneousFits, maxDataPoints);
 110   }
 111
 112   return fgInstance;
 113 }
 114
 115
 116 //______________________________________________________
 117 void AliTPCPIDmathFit::InputData(const TH1 *hh, const Int_t indexXbinRegularisation, const Int_t indexSimultaneousFit,
 118                                  const Double_t leftBoundary, const Double_t rightBoundary, const Double_t threshold,
 119                                  const Bool_t kXerr)
 120 {
 121   // Take index of simultaneousFit and xBinRegularisation and fill the corresponding data points from the histo
 122
 123   if (indexSimultaneousFit < 0 || indexSimultaneousFit >= fNumSimultaneousFits) {
 124     printf("Error: Index %d of simultaneous fit does not exist. Must be 0 - %d!\n", indexSimultaneousFit, fNumSimultaneousFits - 1);
 125     return;
 126   }
 127
 128   if (indexXbinRegularisation < 0 || indexXbinRegularisation >= fNumXbinsRegularisation) {
 129     printf("Error: Index %d of x bin regularisation does not exist. Must be 0 - %d!\n", indexXbinRegularisation,
 130            fNumXbinsRegularisation - 1);
 131     return;
 132   }
 133
 134   fNdata[indexXbinRegularisation][indexSimultaneousFit] = 0;
 135
 136   for (Int_t i = 0; i < fkMaxNdata; i++) {
 137     fX[indexXbinRegularisation][indexSimultaneousFit][i] = -999;
 138     fY[indexXbinRegularisation][indexSimultaneousFit][i] = -999;
 139     fXerr[indexXbinRegularisation][indexSimultaneousFit][i] = -999;
 140     fYerr[indexXbinRegularisation][indexSimultaneousFit][i] = -999;
 141   }
 142
 143   //------
 144
 145   //const Double_t hmax =  hh->GetBinContent(hh->GetMaximumBin());
 146   const TAxis * ax = hh->GetXaxis();
 147
 148   // Only fill bins inside the desired range and only bins with content > threshold
 149   for (Int_t id = ax->GetFirst(); id <= ax->GetLast(); id++) {
 150     if (hh->GetBinCenter(id) < leftBoundary)  continue;
 151     if (hh->GetBinCenter(id) > rightBoundary)  break;
 152
 153     const Double_t yy = hh->GetBinContent(id);
 154     if (yy <= threshold)//"=" important to exclude 0-bins
 155         continue;
 156
 157     fX[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = hh->GetBinCenter(id);
 158     fY[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = yy;
 159
 160     fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = hh->GetBinError(id);
 161
 162     // Small yErr only matters for chi^2 fit. For likelihood fit, just set error to 1, if no error is set.
 163     if (fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]]
 164         < fkDoubleEpsilonLimit) {
 165       if (fUseLogLikelihood) {
 166         fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = 1.0;
 167       }
 168       else if (fY[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]]
 169                < fkDoubleEpsilonLimit) {
 170         Int_t firstNonZeroBin = hh->FindFirstBinAbove(0.0);
 171         Int_t lastNonZeroBin = hh->FindLastBinAbove(0.0);
 172
 173         // Bins with small content usually sit at the edges of the histo
 174         Double_t errorOfBinsWithSmallestContent = 0.0;
 175         if (firstNonZeroBin > 0 && lastNonZeroBin > 0) {
 176           errorOfBinsWithSmallestContent = (hh->GetBinError(firstNonZeroBin) + hh->GetBinError(lastNonZeroBin)) / 2.0;
 177
 178
 179           if (fDebugLevel > 0)
 180             printf("MathFit::InputData strange yerr %e for y %e -> Setting arbitrary error %e!\n",
 181                   fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]],
 182                   fY[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]],
 183                   errorOfBinsWithSmallestContent);
 184
 185           fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] =
 186             errorOfBinsWithSmallestContent;
 187         }
 188         else {
 189           // There is no information, exclude from fit
 190           fX[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = -999;
 191           fY[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = -999;
 192           fXerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = -999;
 193           fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = -999;
 194           continue;
 195         }
 196
 197         //// For empty bins assume bin content = 1, i.e. error = 1
 198         //fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = 1.0;
 199         /*
 200         const Double_t ww = 10.;
 201         fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = hmax / ww;
 202         //printf("MathFit::InputData yerr asigned arbitrary as %f = hamx %f / %f for 0-bin!!\n",
 203                  fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]], hmax, ww);
 204         */
 205       }
 206       else {
 207         printf("AliTPCPIDmathFit::InputData fYerr = %e for fY = %e and fNdata = %d\n",
 208                fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]],
 209                fY[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]], fNdata[indexXbinRegularisation][indexSimultaneousFit]);
 210         exit(1);
 211       }
 212     }
 213     /*
 214     if (abs(fYerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]])
 215         < std::numeric_limits<double>::min()){
 216       printf("AliTPCPIDmathFit::InputData fYerr = 0!! %d\n", id);
 217       exit(1);
 218     }*/
 219
 220     if (kXerr)
 221       fXerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] =
 222         hh->GetBinWidth(id) / 2.0;
 223     else
 224       fXerr[indexXbinRegularisation][indexSimultaneousFit][fNdata[indexXbinRegularisation][indexSimultaneousFit]] = 0.0;
 225
 226     fNdata[indexXbinRegularisation][indexSimultaneousFit]++;
 227     if(fNdata[indexXbinRegularisation][indexSimultaneousFit] >= fkMaxNdata) {
 228       printf("AliTPCPIDmathFit::InputData fNdata >= fkMaxNdata!! %d\n", fkMaxNdata);
 229       exit(1);
 230     }
 231   }
 232 }
 233
 234
 235 //______________________________________________________
 236 Int_t AliTPCPIDmathFit::MinuitFit(FitFunc_t *fn, FitFunc_t *fnError, const Int_t nPar, Double_t *par, Double_t *err,
 237                                   Double_t *covMatrix, Double_t &ChiSquareOrLogLikelihood, Int_t &NDF,
 238                                   const Double_t *stepSize, const Double_t *lowLim, const Double_t *hiLim)
 239 {
 240   // Setup minuit object with in put parameters and internally stored things (like data points) and do the fitting.
 241
 242   Bool_t hasData = kFALSE;
 243   for (Int_t indexXbinRegularisation = 0; indexXbinRegularisation < fNumXbinsRegularisation; indexXbinRegularisation++) {
 244     for (Int_t indexSimultaneousFit = 0; indexSimultaneousFit < fNumSimultaneousFits; indexSimultaneousFit++) {
 245       if (fNdata[indexXbinRegularisation][indexSimultaneousFit] > 0) {
 246         hasData = kTRUE;
 247         break;
 248       }
 249     }
 250   }
 251
 252   if (!hasData) {
 253     printf("AliTPCPIDmathFit::MinuitFit no data yet, call InputData() first!!\n");
 254     return -1;
 255   }
 256
 257   if (fRegularisation > 0 && fNumSimultaneousFits <= 1) {
 258     printf("Regularisation only implemented for simultaneous fit!\n");
 259     return -1;
 260   }
 261
 262   for (Int_t indexSimultaneousFit = 0; indexSimultaneousFit < fNumSimultaneousFits; indexSimultaneousFit++) {
 263     fFunc[indexSimultaneousFit] = fn ? fn[indexSimultaneousFit] : 0x0;
 264     fErrorFunc[indexSimultaneousFit] = fnError ? fnError[indexSimultaneousFit] : 0x0;
 265   }
 266
 267   fNpar = nPar;
 268   fPar = par;
 269   fErr = err;
 270
 271   const Double_t *inpar = fPar;
 272   Int_t errFlag = 0;
 273
 274   // Start the stopwatch. If reset is kTRUE reset the stopwatch before starting it (including the stopwatch counter).
 275   // Use kFALSE to continue timing after a Stop() without resetting the stopwatch.
 276   TStopwatch stwa;
 277   stwa.Start(kTRUE);
 278   errFlag = Fit(inpar, covMatrix, stepSize, lowLim, hiLim);
 279   stwa.Stop();
 280   stwa.Print("m");//"m": milli sec; "u": micro sec
 281
 282   //Store Chi^2/Loglikelihood and NDF
 283   if (fUseLogLikelihood)
 284     ChiSquareOrLogLikelihood = fNegLogLikelihood;
 285   else
 286     ChiSquareOrLogLikelihood = fChi2;
 287   NDF = fNDF;
 288
 289   // Reset internal variables after fit
 290   for (Int_t indexSimultaneousFit = 0; indexSimultaneousFit < fNumSimultaneousFits; indexSimultaneousFit++) {
 291     fFunc[indexSimultaneousFit] = 0x0;
 292     fErrorFunc[indexSimultaneousFit] = 0x0;
 293   }
 294
 295   fNpar = 0;
 296   fPar = 0x0;
 297   fErr = 0x0;
 298
 299   fNstep = 0;
 300
 301   fChi2 = 0;
 302   fNegLogLikelihood = 0;
 303   fNDF = 0;
 304
 305   return errFlag;
 306 }
 307
 308
 309 //______________________________________________________
 310 void AliTPCPIDmathFit::SetEpsilon(const Double_t eps)
 311 {
 312   // Set epsilon defining convergence of fit
 313
 314   if (eps > 0)
 315     fEpsilon = eps;
 316   else
 317     AliError("epsilon must be > 0");
 318 }
 319
 320
 321 //______________________________________________________
 322 void AliTPCPIDmathFit::SetMaxCalls(const Int_t maxCalls)
 323 {
 324   // Set maximum number of calls for fitting
 325
 326   if (maxCalls > 0)
 327     fMaxCalls = maxCalls;
 328   else
 329     AliError("maxCalls must be > 0");
 330 }
 331
 332
 333 //______________________________________________________
 334 Bool_t AliTPCPIDmathFit::SetParametersToRegularise(const Int_t numParams, const Int_t numParamsPerXbin,
 335                                                    const Int_t* indexParametersToRegularise, const Int_t* lastNotFixedIndexOfParameters,
 336                                                    const Double_t* xValuesForRegularisation, const Double_t* xStatisticalWeight,
 337                                                    const Double_t* xStatisticalWeightError)
 338 {
 339   // Sets the number of parameters to the new value and saves the indices of the parameters to regularise in an internal array.
 340   // Also, the number of parameters per xBin is set, the x values of the individual x bins and the last index of the x bin
 341   // in which the corresponding parameter is not fixed plus the statistical weights and errors for each xBin.
 342   // kFALSE is returned and the internal values are reset, if invalid values are provided.
 343
 344   fNumParametersPerXbin = 0;
 345   fNumParametersToRegularise = 0;
 346
 347   delete fIndexParametersToRegularise;
 348   fIndexParametersToRegularise = 0x0;
 349
 350   delete fLastNotFixedIndexOfParameters;
 351   fLastNotFixedIndexOfParameters = 0x0;
 352
 353   delete fXvaluesForRegularisation;
 354   fXvaluesForRegularisation = 0x0;
 355
 356   delete fXstatisticalWeight;
 357   fXstatisticalWeight = 0x0;
 358
 359   delete fXstatisticalWeightError;
 360   fXstatisticalWeightError = 0x0;
 361
 362   if (!indexParametersToRegularise || !xValuesForRegularisation || numParamsPerXbin <= 0 || numParams <= 0
 363       || !xStatisticalWeight || !xStatisticalWeightError || !lastNotFixedIndexOfParameters) {
 364     printf("Error: Cannot set parameters to regularise. Invalid input!\n");
 365
 366     return kFALSE;
 367   }
 368
 369   fIndexParametersToRegularise = new Int_t[numParams];
 370   for (Int_t i = 0; i < numParams; i++)
 371     fIndexParametersToRegularise[i] = indexParametersToRegularise[i];
 372
 373   fLastNotFixedIndexOfParameters = new Int_t[numParamsPerXbin];
 374   for (Int_t i = 0; i < numParamsPerXbin; i++)
 375     fLastNotFixedIndexOfParameters[i] = lastNotFixedIndexOfParameters[i];
 376
 377   fXvaluesForRegularisation = new Double_t[fNumXbinsRegularisation];
 378   fXstatisticalWeight = new Double_t[fNumXbinsRegularisation];
 379   fXstatisticalWeightError = new Double_t[fNumXbinsRegularisation];
 380
 381
 382   Bool_t success = kTRUE;
 383
 384   for (Int_t i = 0; i < fNumXbinsRegularisation; i++) {
 385     fXvaluesForRegularisation[i] = xValuesForRegularisation[i];
 386
 387     if (i > 0 && fXvaluesForRegularisation[i] <= fXvaluesForRegularisation[i - 1]) {
 388       printf("Error: x values for regularisation must be strictly increasing!\n");
 389       success = kFALSE;
 390       break;
 391     }
 392
 393     if (xStatisticalWeight[i] < 0 || xStatisticalWeightError[i] < 0) {
 394       printf("Error: Invalid statistical weight and/or error for bin %d: value %e, error %e\n",
 395              i, xStatisticalWeight[i], xStatisticalWeightError[i]);
 396       success = kFALSE;
 397       break;
 398     }
 399
 400     fXstatisticalWeight[i] = xStatisticalWeight[i];
 401     fXstatisticalWeightError[i] = xStatisticalWeightError[i];
 402   }
 403
 404   if (!success) {
 405     fNumParametersPerXbin = 0;
 406     fNumParametersToRegularise = 0;
 407
 408     delete fIndexParametersToRegularise;
 409     fIndexParametersToRegularise = 0x0;
 410
 411     delete fLastNotFixedIndexOfParameters;
 412     fLastNotFixedIndexOfParameters = 0x0;
 413
 414     delete fXvaluesForRegularisation;
 415     fXvaluesForRegularisation = 0x0;
 416
 417     delete fXstatisticalWeight;
 418     fXstatisticalWeight = 0x0;
 419
 420     delete fXstatisticalWeightError;
 421     fXstatisticalWeightError = 0x0;
 422
 423     return kFALSE;
 424   }
 425
 426   fNumParametersToRegularise = numParams;
 427   fNumParametersPerXbin = numParamsPerXbin;
 428
 429   return kTRUE;
 430 }
 431
 432
 433
 434 //***Private functions***
 435 //______________________________________________________
 436 AliTPCPIDmathFit::AliTPCPIDmathFit(const Int_t numXbinsRegularisation, const Int_t numSimultaneousFits, const Int_t maxDataPoints):
 437 fkDoubleEpsilonLimit(2. * std::numeric_limits<double>::min()),
 438 fFunc(0x0),
 439 fErrorFunc(0x0),
 440 fNpar(0),
 441 fPar(0x0),
 442 fErr(0x0),
 443 fChi2(0),
 444 fNegLogLikelihood(0),
 445 fNDF(0),
 446 fNstep(0),
 447 fNdata(0x0),
 448 fX(0x0),
 449 fY(0x0),
 450 fXerr(0x0),
 451 fYerr(0x0),
 452 fMaxCalls(2000),
 453 fEpsilon(1.e-4),
 454 fUseLogLikelihood(kFALSE),
 455 fUseWeightsForLoglikelihood(kFALSE),
 456 fUseWeightsForLoglikelihoodForCurrentFit(kFALSE),
 457 fScaleFactorError(1.0),
 458 fRegularisation(0),
 459 fXbinIndex(0),
 460 fNumParametersToRegularise(0),
 461 fNumParametersPerXbin(0),
 462 fIndexParametersToRegularise(0x0),
 463 fLastNotFixedIndexOfParameters(0x0),
 464 fXvaluesForRegularisation(0x0),
 465 fXstatisticalWeight(0x0),
 466 fXstatisticalWeightError(0x0),
 467 fRegularisationFactor(1),
 468 fMinimisationString("MINIMIZE"),
 469 fkMaxNdata(maxDataPoints),
 470 fNumSimultaneousFits(numSimultaneousFits),
 471 fNumXbinsRegularisation(numXbinsRegularisation),
 472 fMinuit(0x0),
 473 fPolIntX(0x0),
 474 fPolIntY(0x0),
 475 fPolIntC(0x0),
 476 fPolIntD(0x0),
 477 fDebugLevel(0),
 478 fRefHistos(0x0),
 479 fNrefHistos(0)
 480 {
 481   // Private ctor
 482
 483   if (fNumSimultaneousFits < 1) {
 484     printf("Error: numSimultaneousFits must be >= 1! Is now set to 1!\n\n");
 485     fNumSimultaneousFits = 1;
 486   }
 487
 488   fX = new Double_t**[fNumXbinsRegularisation];
 489   fY = new Double_t**[fNumXbinsRegularisation];
 490   fXerr = new Double_t**[fNumXbinsRegularisation];
 491   fYerr = new Double_t**[fNumXbinsRegularisation];
 492
 493   fNdata = new Int_t*[fNumXbinsRegularisation];
 494
 495   fFunc = new FitFunc_t[fNumSimultaneousFits];
 496   fErrorFunc = new FitFunc_t[fNumSimultaneousFits];
 497
 498   for (Int_t i = 0; i < fNumXbinsRegularisation; i++) {
 499     fX[i] = new Double_t*[fNumSimultaneousFits];
 500     fY[i] = new Double_t*[fNumSimultaneousFits];
 501     fXerr[i] = new Double_t*[fNumSimultaneousFits];
 502     fYerr[i] = new Double_t*[fNumSimultaneousFits];
 503
 504     fNdata[i] = new Int_t[fNumSimultaneousFits];
 505
 506     for (Int_t j = 0; j < fNumSimultaneousFits; j++) {
 507       fX[i][j] = new Double_t[fkMaxNdata];
 508       fY[i][j] = new Double_t[fkMaxNdata];
 509       fXerr[i][j] = new Double_t[fkMaxNdata];
 510       fYerr[i][j] = new Double_t[fkMaxNdata];
 511
 512       for (Int_t k = 0; k < fkMaxNdata; k++) {
 513         fX[i][j][k] = -999;
 514         fY[i][j][k] = -999;
 515         fXerr[i][j][k] = -999;
 516         fYerr[i][j][k] = -999;
 517       }
 518
 519       fNdata[i][j] = 0;
 520     }
 521   }
 522
 523   for (Int_t i = 0; i < fNumSimultaneousFits; i++) {
 524     fFunc[i] = 0x0;
 525     fErrorFunc[i] = 0x0;
 526   }
 527 }
 528
 529
 530 //______________________________________________________
 531 AliTPCPIDmathFit::~AliTPCPIDmathFit()
 532 {
 533   // dtor
 534
 535   ClearRefHistos();
 536
 537   for (Int_t i = 0; i < fNumXbinsRegularisation; i++) {
 538     for (Int_t j = 0; j < fNumSimultaneousFits; j++) {
 539       delete [] fX[i][j];
 540       delete [] fY[i][j];
 541       delete [] fXerr[i][j];
 542       delete [] fYerr[i][j];
 543
 544       fX[i][j]  = 0x0;
 545       fY[i][j]  = 0x0;
 546       fXerr[i][j]  = 0x0;
 547       fYerr[i][j]  = 0x0;
 548     }
 549
 550     delete [] fX[i];
 551     delete [] fY[i];
 552     delete [] fXerr[i];
 553     delete [] fYerr[i];
 554
 555     fX[i]  = 0x0;
 556     fY[i]  = 0x0;
 557     fXerr[i]  = 0x0;
 558     fYerr[i]  = 0x0;
 559   }
 560
 561
 562   delete [] fX;
 563   fX = 0x0;
 564
 565   delete [] fY;
 566   fY = 0x0;
 567
 568   delete [] fXerr;
 569   fXerr = 0x0;
 570
 571   delete [] fYerr;
 572   fYerr = 0x0;
 573
 574
 575   for (Int_t i = 0; i < fNumXbinsRegularisation; i++) {
 576     delete [] fNdata[i];
 577
 578     fNdata[i]  = 0x0;
 579   }
 580
 581   delete [] fNdata;
 582   fNdata = 0x0;
 583
 584
 585
 586   delete [] fFunc;
 587   fFunc = 0x0;
 588
 589   delete [] fErrorFunc;
 590   fErrorFunc = 0x0;
 591
 592   delete [] fIndexParametersToRegularise;
 593   fIndexParametersToRegularise = 0x0;
 594
 595   delete [] fLastNotFixedIndexOfParameters;
 596   fLastNotFixedIndexOfParameters = 0x0;
 597
 598   delete [] fXvaluesForRegularisation;
 599   fXvaluesForRegularisation = 0;
 600
 601   delete fMinuit;
 602   fMinuit = 0x0;
 603
 604   delete fXstatisticalWeight;
 605   fXstatisticalWeight = 0x0;
 606
 607   delete fXstatisticalWeightError;
 608   fXstatisticalWeightError = 0x0;
 609
 610   delete fPolIntX;
 611   fPolIntX = 0x0;
 612
 613   delete fPolIntY;
 614   fPolIntY = 0x0;
 615
 616   delete fPolIntC;
 617   fPolIntC = 0x0;
 618
 619   delete fPolIntD;
 620   fPolIntD = 0x0;
 621
 622   fNumXbinsRegularisation = 0;
 623   fNumSimultaneousFits = 0;
 624
 625   fgInstance = 0x0;
 626 }
 627
 628
 629 //______________________________________________________
 630 inline Double_t AliTPCPIDmathFit::EvalLog(Double_t x) const
 631 {
 632   // Safely evaluate logarithms -> Adapted root util.h
 633
 634   if (x < fkDoubleEpsilonLimit)
 635     return x / fkDoubleEpsilonLimit + TMath::Log(fkDoubleEpsilonLimit) - 1.;
 636   else
 637     return TMath::Log(x);
 638 }
 639
 640
 641 //______________________________________________________
 642 Int_t AliTPCPIDmathFit::Fit(const Double_t *inPar, Double_t *covMatrix, const Double_t *stepSize, const Double_t *lowLim,
 643                             const Double_t *hiLim)
 644 {
 645   // Fit the currently loaded data, intialise the parameters with inPar (stepSize, lowLim and hiLim).
 646   // Save covariance matrix to covMatrix.
 647
 648   Int_t fNdataTotal = 0;
 649
 650   for (Int_t iNumXbinRegularisation = 0; iNumXbinRegularisation < fNumXbinsRegularisation; iNumXbinRegularisation++) {
 651     for (Int_t iNumSimultaneousFits = 0; iNumSimultaneousFits < fNumSimultaneousFits; iNumSimultaneousFits++) {
 652       fNdataTotal += fNdata[iNumXbinRegularisation][iNumSimultaneousFits];
 653
 654       for (Int_t idata = 0; idata < fNdata[iNumXbinRegularisation][iNumSimultaneousFits]; idata++) {
 655         if (fDebugLevel > 1)
 656           printf("%3d %2d %7d %13.6e %13.6e %10.2e %10.2e\n", iNumXbinRegularisation, iNumSimultaneousFits,
 657                  idata, fX[iNumXbinRegularisation][iNumSimultaneousFits][idata],
 658                  fY[iNumXbinRegularisation][iNumSimultaneousFits][idata], fXerr[iNumXbinRegularisation][iNumSimultaneousFits][idata],
 659                  fYerr[iNumXbinRegularisation][iNumSimultaneousFits][idata]);
 660       }
 661     }
 662   }
 663
 664   AliTPCPIDmathFit* ptr = AliTPCPIDmathFit::Instance();
 665   printf("\n================= Starting MathFit useLogLikelihood %d ndata %d =================\n\n", ptr->GetUseLogLikelihood(),
 666          fNdataTotal);
 667
 668   fMinuit = new TMinuit(fNpar);
 669   fMinuit->SetFCN(FitFCN);
 670
 671   // Standard error is for chi^2 -> Need factor 1/2 for log likelihood
 672   fMinuit->SetErrorDef(ptr->GetUseLogLikelihood() ? 0.5 : 1);
 673
 674   // Define the start parameters
 675   for (Int_t i = 0; i < fNpar; i++)  {
 676     const Double_t start = inPar ? inPar[i] : 0.5;
 677     // Here is only the INITIAL step size
 678     const Double_t ss = stepSize ? stepSize[i] : 0.01;
 679     const Double_t low = lowLim ? lowLim[i] : 0;
 680     const Double_t hi = hiLim ? hiLim[i] : 0;
 681
 682     fMinuit->DefineParameter(i, Form("p%d", i), start, ss, low, hi);
 683   }
 684
 685   fMinuit->SetMaxIterations(10000);
 686
 687   Double_t arg[2] = {fMaxCalls, fEpsilon};
 688   Int_t errFlag = 0;
 689   Int_t hesseFlag = -999;
 690
 691
 692
 693   // If regularisation is to be used, initialise the function for interpolation.
 694   // Use a pol(2*n-1) for the interpolation, if +-n bins (set by fRegularisation) are taken into account,
 695   // i.e. 2 * fRegularisation points are used.
 696   if (GetUseRegularisation()) {
 697     if (!fXvaluesForRegularisation) {
 698       printf("Cannot calculate regularisation: x value array not initialised!\n");
 699       exit(-1);
 700     }
 701
 702     // NOTE: The interpolation routine is based on self-defined vectors, which have indices from 1 to n.
 703     // In order not to mess something up, just add one dummy entry at index zero
 704     fPolIntX = new Double_t[fRegularisation * 2 + 1];
 705     fPolIntY = new Double_t[fRegularisation * 2 + 1];
 706     fPolIntC = new Double_t[fRegularisation * 2 + 1];
 707     fPolIntD = new Double_t[fRegularisation * 2 + 1];
 708
 709     for (Int_t i = 0; i < fRegularisation * 2 + 1; i++) {
 710       fPolIntX[i] = 0.;
 711       fPolIntY[i] = 0.;
 712       fPolIntC[i] = 0.;
 713       fPolIntD[i] = 0.;
 714     }
 715
 716     printf("Regularisation enabled! Interpolating with pol%d between +-%d bins and regularisation factor %f.\n\n",
 717            2 * fRegularisation - 1, fRegularisation, fRegularisationFactor);
 718   }
 719   else
 720      printf("Regularisation disabled!\n\n");
 721
 722   // Using "MINIMIZE" means to try "MIGRAD". Only if MIGRAD fails,
 723   // SIMPLEX and then MIGRAD are applied. This yields to better convergence,
 724   // but in case of MIGRAD failing the first time, errors are not reliable and
 725   // also the minimum might not be correct. But this is most likely still better
 726   // than a failed fit.
 727
 728   if (fUseWeightsForLoglikelihood && fUseLogLikelihood) {
 729
 730     Double_t covMatrixTemp[fNpar][fNpar];
 731     Double_t hesseMatrixTemp[fNpar][fNpar];
 732     Double_t matrixTemp[fNpar][fNpar];
 733     Double_t covMatrixFinal[fNpar][fNpar];
 734     for (Int_t i = 0; i < fNpar; i++) {
 735       for (Int_t j = 0; j < fNpar; j++) {
 736         covMatrixTemp[i][j] = 0;
 737         hesseMatrixTemp[i][j] = 0;
 738         matrixTemp[i][j] = 0;
 739         covMatrixFinal[i][j] = 0;
 740       }
 741     }
 742
 743
 744     fUseWeightsForLoglikelihoodForCurrentFit = kFALSE;
 745
 746     // Firstly, minimise w/o weighting and obtain covariance matrix
 747     fMinuit->mnexcm(fMinimisationString.Data(), arg, 2, errFlag);
 748     GetCovarianceMatrix(&covMatrixTemp[0][0]);
 749
 750     // Save results (to be used if further steps fail)
 751     for(Int_t ipar = 0; ipar < fNpar; ipar++)
 752       fMinuit->GetParameter(ipar, fPar[ipar], fErr[ipar]);
 753
 754     if (covMatrix) {
 755       GetCovarianceMatrix(covMatrix);
 756
 757
 758       // If the parameter errors (NOT(!) the errors for the data points in fY!!!) are to be scaled, we need to do this here!
 759       for (Int_t i = 0; i < fNpar; ++i) {
 760         for (Int_t j = 0; j < fNpar; ++j) {
 761             covMatrix[i * fNpar + j] *= fScaleFactorError * fScaleFactorError;
 762         }
 763       }
 764     }
 765
 766     // If the parameter errors (NOT(!) the errors for the data points in fY!!!) are to be scaled, we need to do this here!
 767     for (Int_t i = 0; i < fNpar; ++i)
 768       fErr[i] *= fScaleFactorError;
 769
 770
 771     // Secondly, enable weithing and obtain HESSE matrix
 772     fUseWeightsForLoglikelihoodForCurrentFit = kTRUE;
 773
 774
 775     // TODO THIS CALL IS NOT NECESSARILY NEEDED fMinuit->mnexcm(fMinimisationString.Data(), arg, 2, errFlag2);
 776
 777     //improve error bars a lot!
 778     //only 1 argument: fMaxCalls
 779     fMinuit->mnexcm("HESSE", arg, 1, hesseFlag);
 780
 781     if (hesseFlag == 0) {
 782       // If HESSE call was successful, save results (to be used if further steps fail)
 783       for(Int_t ipar = 0; ipar < fNpar; ipar++)
 784         fMinuit->GetParameter(ipar, fPar[ipar], fErr[ipar]);
 785
 786       if (covMatrix) {
 787         GetCovarianceMatrix(covMatrix);
 788
 789         // If the parameter errors (NOT(!) the errors for the data points in fY!!!) are to be scaled, we need to do this here!
 790         for (Int_t i = 0; i < fNpar; ++i) {
 791           for (Int_t j = 0; j < fNpar; ++j) {
 792               covMatrix[i * fNpar + j] *= fScaleFactorError * fScaleFactorError;
 793           }
 794         }
 795       }
 796
 797       // If the parameter errors (NOT(!) the errors for the data points in fY!!!) are to be scaled, we need to do this here!
 798       for (Int_t i = 0; i < fNpar; ++i)
 799         fErr[i] *= fScaleFactorError;
 800     }
 801
 802
 803     // Hesse matrix is the inverse of the error matrix after the hesse step
 804
 805     Int_t nfree = fMinuit->GetNumFreePars();
 806     TMatrixDSym orig(nfree);
 807     fMinuit->mnemat(orig.GetMatrixArray(), nfree);
 808
 809     TMatrixDSym inv(orig);
 810     inv.Invert();
 811
 812     // If the inversion failed, the result is equal to the original
 813     Bool_t success = kFALSE;
 814     for (Int_t i = 0; i < nfree; ++i) {
 815       for (Int_t j = 0; j < nfree; ++j) {
 816         if (TMath::Abs(inv(i,j) - orig(i,j)) > 1e-9) {
 817           success = kTRUE;
 818           break;
 819         }
 820       }
 821
 822       if (success)
 823         break;
 824     }
 825
 826
 827     // Only re-weight errors, if the matrix inversion and the fits have been successful.
 828     if ((errFlag + hesseFlag) == 0 && success) {
 829       Int_t l = 0;
 830       for (Int_t i = 0; i < fNpar; ++i) {
 831         if (fMinuit->fNiofex[i] > 0) {  // not fixed ?
 832           Int_t m = 0;
 833           for (Int_t j = 0; j <= i; ++j) {
 834             if (fMinuit->fNiofex[j] > 0) {  //not fixed
 835                 hesseMatrixTemp[i][j] =  inv(l, m);
 836                 hesseMatrixTemp[j][i] = hesseMatrixTemp[i][j];
 837                 m++;
 838             }
 839           }
 840           l++;
 841         }
 842       }
 843
 844       // We want cov * hes * cov
 845       // -> Do hes * cov first
 846       for (Int_t i = 0; i < fNpar; ++i) {
 847         for (Int_t j = 0; j < fNpar; ++j) {
 848             for (Int_t k = 0; k < fNpar; ++k)
 849               matrixTemp[i][j] += hesseMatrixTemp[i][k] * covMatrixTemp[k][j];
 850         }
 851       }
 852
 853       // Now do cov * tmp
 854       for (Int_t i = 0; i < fNpar; ++i) {
 855         for (Int_t j = 0; j < fNpar; ++j) {
 856           for (Int_t k = 0; k < fNpar; ++k) {
 857             covMatrixFinal[i][j] += covMatrixTemp[i][k] * matrixTemp[k][j];
 858           }
 859
 860           // If the parameter errors (NOT(!) the errors for the data points in fY!!!) are to be scaled, we need to do this here!
 861           covMatrixFinal[i][j] *= fScaleFactorError * fScaleFactorError;
 862
 863           if (covMatrix)
 864             covMatrix[i * fNpar + j] = covMatrixFinal[i][j];
 865         }
 866       }
 867
 868       // Obtain the parameters (and old errors)
 869       for(Int_t ipar = 0; ipar < fNpar; ipar++)
 870         fMinuit->GetParameter(ipar, fPar[ipar], fErr[ipar]);
 871
 872       // Set the new parameter errors
 873       for(Int_t ipar = 0; ipar < fNpar; ipar++)
 874         fErr[ipar] = TMath::Sqrt(covMatrixFinal[ipar][ipar]);
 875     }
 876     else {
 877       // An error occurred -> Leave the errors as they are (still better than doing re-weighting with problematic matrices)
 878       printf("Error: Could not re-weight errors (MIGRAD %d, HESSE %d, matrix inversion %d) -> Errors and covariance matrix not reliable!\n",
 879              errFlag, hesseFlag, success);
 880     }
 881   }
 882   else {
 883     fMinuit->mnexcm(fMinimisationString.Data(), arg, 2, errFlag);
 884
 885     //improve error bars a lot!
 886     //only 1 argument: fMaxCalls
 887     fMinuit->mnexcm("HESSE", arg, 1, hesseFlag);
 888
 889     for(Int_t ipar = 0; ipar < fNpar; ipar++)
 890       fMinuit->GetParameter(ipar, fPar[ipar], fErr[ipar]);
 891
 892     if (covMatrix) {
 893       GetCovarianceMatrix(covMatrix);
 894
 895       // If the parameter errors (NOT(!) the errors for the data points in fY!!!) are to be scaled, we need to do this here!
 896       for (Int_t i = 0; i < fNpar; ++i) {
 897         for (Int_t j = 0; j < fNpar; ++j) {
 898             covMatrix[i * fNpar + j] *= fScaleFactorError * fScaleFactorError;
 899         }
 900       }
 901     }
 902
 903     // If the parameter errors (NOT(!) the errors for the data points in fY!!!) are to be scaled, we need to do this here!
 904     for (Int_t i = 0; i < fNpar; ++i)
 905       fErr[i] *= fScaleFactorError;
 906   }
 907
 908   fNDF = fNdataTotal - fNpar;
 909   Double_t fChi2orLogLikelihood = 0;
 910
 911   if (!fUseLogLikelihood) {
 912     GetChiSquare();
 913     fChi2orLogLikelihood = fChi2;
 914   }
 915   else {
 916     GetNegLogLikelihood();
 917     fChi2orLogLikelihood = fNegLogLikelihood;
 918   }
 919
 920
 921   if ((errFlag + hesseFlag) != 0) printf("*****Warning: Abnormal termination: MIGRAD flag %d, HESSE flag %d\n", errFlag, hesseFlag);
 922
 923
 924   printf("======> ChiSquare/NDF: %f / %d%s\nerrFlag %d (MIGRAD %d, HESSE %d)\n", fChi2orLogLikelihood, fNDF,
 925          fNDF > 0 ? Form(" = %f", fChi2orLogLikelihood / fNDF) : "", errFlag + hesseFlag, errFlag, hesseFlag);
 926   errFlag += hesseFlag;
 927
 928   delete fMinuit;
 929   fMinuit = 0x0;
 930
 931   delete fPolIntX;
 932   fPolIntX = 0x0;
 933
 934   delete fPolIntY;
 935   fPolIntY = 0x0;
 936
 937   delete fPolIntC;
 938   fPolIntC = 0x0;
 939
 940   delete fPolIntD;
 941   fPolIntD = 0x0;
 942
 943   fXbinIndex = 0;
 944
 945
 946   return errFlag;
 947 }
 948
 949
 950 //______________________________________________________
 951 void AliTPCPIDmathFit::FitFCN(Int_t &nPar, Double_t *, Double_t &chiSquareOrLogLikelihood, Double_t *par, Int_t)
 952 {
 953   // Perform fitting step
 954
 955   AliTPCPIDmathFit* ptr = AliTPCPIDmathFit::Instance();
 956   ptr->fNstep++;
 957
 958   for (Int_t ii = 0; ii < ptr->fNpar; ii++)
 959     ptr->fPar[ii] = par[ii];
 960
 961   if (ptr->GetUseLogLikelihood()) {
 962     chiSquareOrLogLikelihood = ptr->GetNegLogLikelihood();
 963   }
 964   else {
 965     chiSquareOrLogLikelihood = ptr->GetChiSquare();
 966   }
 967
 968   if (ptr->GetDebugLevel() >= 3) {
 969     printf("%d %f %d --- ", ptr->fNstep, chiSquareOrLogLikelihood, nPar);
 970
 971     for (Int_t i = 0; i < nPar; i++)
 972       printf("%f ", par[i]);
 973
 974     printf("-----------------\n");
 975
 976     if (ptr->fNstep == 2 && ptr->GetDebugLevel() >= 4)
 977        exit(1);
 978   }
 979 }
 980
 981
 982 //______________________________________________________
 983 Double_t AliTPCPIDmathFit::GetChiSquare()
 984 {
 985   // Calculates chi^2 and sets this value for the internal variable fChi2. In addition, the calculated value is returned.
 986
 987   fChi2 = 0;
 988
 989   for (Int_t xBin = 0; xBin < fNumXbinsRegularisation; xBin++) {
 990     // This parameter tells the following fFunc call which x bin (and which parameters) are to be used for the fit
 991     // of the current x bin
 992
 993     fXbinIndex = xBin;
 994
 995     for (Int_t iSimultaneousFit = 0; iSimultaneousFit < fNumSimultaneousFits; iSimultaneousFit++) {
 996       for(Int_t idata = 0; idata < fNdata[xBin][iSimultaneousFit]; idata++) {
 997         const Double_t evaly = fFunc[iSimultaneousFit](&(fX[xBin][iSimultaneousFit][idata]), fPar);
 998         const Double_t dy = fY[xBin][iSimultaneousFit][idata] - evaly;
 999
1000         Double_t dfdx = 0;
1001         //calculated df/dx if necessary
1002         if (fXerr[xBin][iSimultaneousFit][idata] != 0)  {
1003           const Double_t derix0 = fX[xBin][iSimultaneousFit][idata] - fXerr[xBin][iSimultaneousFit][idata];
1004           const Double_t derix1 = fX[xBin][iSimultaneousFit][idata] + fXerr[xBin][iSimultaneousFit][idata];
1005           const Double_t deriy0 = fFunc[iSimultaneousFit](&derix0, fPar);
1006           const Double_t deriy1 = fFunc[iSimultaneousFit](&derix1, fPar);
1007           if (abs(derix1 - derix0) < fkDoubleEpsilonLimit) {
1008             AliError("derix1-derix0 is 0 -> Devision by zero!");
1009             dfdx = 0;
1010           }
1011           else
1012             dfdx = (deriy1 - deriy0) / (derix1 - derix0);
1013         }
1014         Double_t den = TMath::Power(fYerr[xBin][iSimultaneousFit][idata], 2) +
1015                        TMath::Power(dfdx * fXerr[xBin][iSimultaneousFit][idata], 2);
1016
1017         if (fErrorFunc[iSimultaneousFit])  {
1018           den += fErrorFunc[iSimultaneousFit](&(fX[xBin][iSimultaneousFit][idata]), fPar);
1019         }
1020
1021         Double_t part = 0;
1022
1023         if (abs(den) < fkDoubleEpsilonLimit) {
1024           AliError("den is 0 -> Devision by zero!");
1025           part = 0;
1026         }
1027         else
1028           part = dy * dy / den;
1029
1030         if (fDebugLevel >= 4)  {
1031           for(Int_t ip = 0; ip < fNpar; ip++) {
1032             printf("par%d %e\n", ip, fPar[ip]);
1033           }
1034           printf("%d %d: %e %e %e -- %e %e -- %e %e %e -- %e\n", iSimultaneousFit, idata, fX[xBin][iSimultaneousFit][idata],
1035                  fY[xBin][iSimultaneousFit][idata], fYerr[xBin][iSimultaneousFit][idata], evaly, dy, dfdx, den, part, fChi2);
1036         }
1037
1038         fChi2 += part;
1039       }
1040     }
1041   }
1042
1043   AddPenaltyTermForRegularisation(kFALSE);
1044
1045   return fChi2;
1046 }
1047
1048
1049 //______________________________________________________
1050 Double_t AliTPCPIDmathFit::GetNegLogLikelihood()
1051 {
1052   // Implementation based on ROOT's FitUtil::EvaluatePoissonLogL in $ROOTSYS/math/mathcore/src/FitUtil.cxx.
1053   // Especially the weighted branch is based on the root version.
1054
1055   fNegLogLikelihood = 0;
1056
1057   if (fUseWeightsForLoglikelihoodForCurrentFit) {
1058     for (Int_t xBin = 0; xBin < fNumXbinsRegularisation; xBin++) {
1059       // This parameter tells the following fFunc call which x bin (and which parameters) are to be used for the fit
1060       // of the current x bin
1061
1062       fXbinIndex = xBin;
1063
1064       for (Int_t iSimultaneousFit = 0; iSimultaneousFit < fNumSimultaneousFits; iSimultaneousFit++) {
1065         for(Int_t idata = 0; idata < fNdata[xBin][iSimultaneousFit]; idata++) {
1066
1067           // Need to apply weight correction. Effective weight is error^2 / y
1068           // and expected events in bins is evalY/weight.
1069           // One can apply the correction only when y is not zero otherwise weight is undefined
1070           // (in case of weighted likelihood one doesn't care about the constant term due to
1071           // the saturated model).
1072
1073           if (TMath::Abs(fY[xBin][iSimultaneousFit][idata]) > fkDoubleEpsilonLimit) {
1074             Double_t evalY = fFunc[iSimultaneousFit](&(fX[xBin][iSimultaneousFit][idata]), fPar);
1075             evalY = TMath::Max(evalY, 0.0);  // Avoid negative or too small values
1076
1077             // Bin effective weight
1078             const Double_t weight = (fYerr[xBin][iSimultaneousFit][idata] * fYerr[xBin][iSimultaneousFit][idata]) /
1079                                     fY[xBin][iSimultaneousFit][idata];
1080             fNegLogLikelihood += weight * evalY - weight * fY[xBin][iSimultaneousFit][idata] * EvalLog(evalY);
1081           }
1082         }
1083       }
1084     }
1085   }
1086   else {
1087     // Binned Poissonian likelihood chi-square as in http://www.scribd.com/doc/65048561/1983-S-Baker-R-D-Cousins-NIM-221-437-442
1088     // -> Download http://www.sciencedirect.com/science/article/pii/0167508784900164#
1089     for (Int_t xBin = 0; xBin < fNumXbinsRegularisation; xBin++) {
1090       // This parameter tells the following fFunc call which x bin (and which parameters) are to be used for the fit
1091       // of the current x bin
1092       fXbinIndex = xBin;
1093
1094       for (Int_t iSimultaneousFit = 0; iSimultaneousFit < fNumSimultaneousFits; iSimultaneousFit++) {
1095         for(Int_t idata = 0; idata < fNdata[xBin][iSimultaneousFit]; idata++) {
1096           Double_t evalY = fFunc[iSimultaneousFit](&(fX[xBin][iSimultaneousFit][idata]), fPar);
1097           // Since the fit MINIMIZES, we need the negative value of the log likelihood here,
1098           // i.e. minus sign already taken into account in the following!
1099           evalY = TMath::Max(evalY, 0.0);  // Avoid negative or too small values
1100
1101           Double_t part = evalY - fY[xBin][iSimultaneousFit][idata];
1102           if (fY[xBin][iSimultaneousFit][idata] > 0.0) {
1103             part += fY[xBin][iSimultaneousFit][idata] * (EvalLog(fY[xBin][iSimultaneousFit][idata]) - EvalLog(evalY));
1104           }
1105
1106           if(fDebugLevel >= 4)
1107             printf("%d %d: %e -- %e %e %e\n", iSimultaneousFit, idata, fX[xBin][iSimultaneousFit][idata], evalY, part,
1108                    fNegLogLikelihood);
1109
1110           // NO factor 2 as in paper. This is because the error definition is set to 0.5, which exactly absorbs the parameter
1111           // 2 here (error means that changing the parameters by the error, the negLogLikelihood changes by the error definition)
1112           // -> So changing the negLogLikelihood as defined here by 1 sigma produces 0.5 of the change of the real negLoglikelihood
1113           // with the factor 2 in front.
1114           fNegLogLikelihood += part;
1115         }
1116       }
1117     }
1118   }
1119
1120   AddPenaltyTermForRegularisation(kTRUE);
1121
1122   return fNegLogLikelihood;
1123 }
1124
1125
1126 //______________________________________________________
1127 void AliTPCPIDmathFit::AddPenaltyTermForRegularisation(const Bool_t useLogLikelihood)
1128 {
1129   // In case of regularisation, add the corresponding penalty term to chi^2 or negLogLikelihood.
1130
1131   if (GetUseRegularisation()) {
1132     // In the following, the roles of x and y are different from the actual fit!
1133     // In particular, y is now the value of a certain parameter!
1134
1135     // Loop over all x bins to obtain the penalty term
1136     for (Int_t iPar = 0; iPar < fNumParametersToRegularise; iPar++) {
1137       const Int_t currParIndex = fIndexParametersToRegularise[iPar];
1138
1139       // the xBin of the current parameter can be obtained from the parameter index
1140       // and the number of parameters per xBin
1141       const Int_t xBin = (Int_t)(currParIndex / fNumParametersPerXbin);
1142
1143       if (fXstatisticalWeight[xBin] < fkDoubleEpsilonLimit)
1144           continue; // If there is no statistics in this bin, one can safely skip it (the term cannot be calculated anyway)
1145
1146       // Determine the number of points used for the regularisation (= order of interpolation polynomial)
1147       const Int_t nRegBins = 2 * fRegularisation; // Regularise with +-fRegularisation bins
1148
1149       // No regularisation at the edges, i.e. less than fRegularisation bins left or right of the current bin
1150       // => This is IMPORTANT. Consider for example fRegularisation = 1, i.e. +-1 bins. If the left bin is missing,
1151       // the "inter"polation polynomial would be a pol0, i.e. a constant. But for low pT there is some finite slope
1152       // at the edges. So, extrapolating to the edges would bias the results.
1153       // The edge points are constraint by their neighbours in the sense that they provide an interpolation value.
1154       // A bin is also an "edge bin", if the right neighbour bin is fixed
1155       if (xBin - fRegularisation < 0)
1156         continue; // Points left of the current bin are missing
1157       if (xBin + fRegularisation > fNumXbinsRegularisation - 1)
1158         continue; // Points right of the current bin are missing
1159       if (xBin + fRegularisation > fLastNotFixedIndexOfParameters[currParIndex % fNumParametersPerXbin])
1160         continue; // Points right of the current bin cannot be used, since they are fixed
1161
1162       /* OLD
1163       Int_t nRegBins = 2 * fRegularisation; // Regularise with +-fRegularisation bins
1164       if (xBin - fRegularisation < 0)
1165         nRegBins -= fRegularisation - xBin; // Points left of the current bin are missing
1166       if (xBin + fRegularisation > fNumXbinsRegularisation - 1)
1167         nRegBins -= (xBin + fRegularisation) - (fNumXbinsRegularisation - 1); // Points right of the current bin are missing
1168
1169       // To perform the interpolation, the interpolation function must be well determined.
1170       // This is only true for nRegBins >= DOF(pol(2 * fRegularisation - 1)) = 2 * fRegularisation.
1171       // Therefore, bins at the edges will not be regularised.
1172
1173       if (nRegBins < 2 * fRegularisation)
1174         continue;
1175
1176
1177
1178       Int_t xBinTemp = TMath::Max(0, xBin - fRegularisation);
1179       */
1180
1181       // Bin effective weight required for weighted log likelihood
1182       const Double_t weight = fUseWeightsForLoglikelihoodForCurrentFit
1183                                 ? (fXstatisticalWeightError[xBin] * fXstatisticalWeightError[xBin] / fXstatisticalWeight[xBin])
1184                                 : 1.0;
1185
1186       // Current values of this parameter
1187       Double_t currYerr = 0, currYvalue = 0;
1188       fMinuit->GetParameter(currParIndex, currYvalue, currYerr);
1189
1190
1191       Double_t y, yErr;
1192
1193       // Set the data points for the interpolation
1194
1195       // Reset first; i = 0 is not used!
1196       for (Int_t i = 1; i < fRegularisation * 2 + 1; i++) {
1197         fPolIntX[i] = 0.;
1198         fPolIntY[i] = 0.;
1199       }
1200
1201       Int_t xBinTemp = xBin - fRegularisation; // Above if statement already guarantees xBinTemp >= 0
1202       for (Int_t i = 0; i < nRegBins; xBinTemp++) {
1203         if (xBinTemp == xBin) // Do not add considered bin to interpolation
1204           continue;
1205
1206         y = 0;
1207         yErr = 0;
1208
1209         fMinuit->GetParameter(currParIndex % fNumParametersPerXbin + xBinTemp * fNumParametersPerXbin, y, yErr);
1210
1211         // NOTE: +1 to ignore dummy index 0!
1212         fPolIntX[i + 1] = fXvaluesForRegularisation[xBinTemp];
1213         fPolIntY[i + 1] = y;
1214
1215         i++;
1216       }
1217
1218       // Get the interpolated y value
1219       Double_t yInterpolated = 0., yInterPolatedErr = 0;
1220       if (!PolynomialInterpolation(fPolIntX, fPolIntY, nRegBins, fXvaluesForRegularisation[xBin], &yInterpolated, &yInterPolatedErr))
1221       {
1222         printf("Polynomial interpolation failed! Regularisation term for current x bin (%d) ignored\n", xBin);
1223         continue;
1224       }
1225
1226       // Calculate the penalty term...
1227       const Double_t yDiff = currYvalue - yInterpolated;
1228
1229       // For the error, just take the statistical error (assuming the statistical weight for this bin to have no error, since
1230       // this is only a number of measured counts, which is a fact). A weighting correction will be applied for the
1231       // penalty term.
1232
1233       const Double_t yError2 = (0.5 * (currYvalue + yInterpolated)) / fXstatisticalWeight[xBin];
1234
1235       const Double_t regPenaltyTerm = fRegularisationFactor * weight * (yDiff * yDiff) / yError2;
1236
1237       // ...and add it to the negLogLikelihood or the chi^2 depending on the fit method.
1238       // The factor 0.5 for the log likelihood is to take the proper definition of negLogLikelihood
1239       // (which is 0.5), i.e. changing by yError should change the penalty term by 0.5.
1240       if (useLogLikelihood)
1241         fNegLogLikelihood += 0.5 * regPenaltyTerm;
1242       else
1243         fChi2 += regPenaltyTerm;
1244     }
1245   }
1246 }
1247
1248
1249 //______________________________________________________
1250 void AliTPCPIDmathFit::GetCovarianceMatrix(Double_t* covMatrix) const
1251 {
1252   // Retrieve covariance matrix.
1253
1254   if (!fMinuit || !covMatrix)
1255     return;
1256
1257   // Take special care in case of fixed parameters, i.e. fill in zeros, compare
1258   // http://root.cern.ch/root/html/src/TMinuitMinimizer.cxx.html#Ilx.Z
1259
1260   const Int_t nfree = fMinuit->GetNumFreePars();
1261
1262   TMatrixDSym matfree(nfree);
1263   fMinuit->mnemat(matfree.GetMatrixArray(), nfree);
1264
1265   // r: row, c: column
1266   Int_t rfree = 0;
1267   for (Int_t rr = 0; rr < fNpar; rr++) {
1268     if (fMinuit->fNiofex[rr] <= 0) // fMinuit->fNiofex[rr] is zero in case of fix parameter
1269       continue;
1270
1271     Int_t cfree = 0;
1272     for (Int_t cc = 0; cc < fNpar; cc++) {
1273       if (fMinuit->fNiofex[cc] <= 0) // fMinuit->fNiofex[cc] is zero in case of fix parameter
1274         continue;
1275
1276       covMatrix[rr * fNpar + cc] = matfree[rfree][cfree];
1277       cfree++;
1278     }
1279     rfree++;
1280   }
1281 }
1282
1283
1284 //______________________________________________________
1285 Bool_t AliTPCPIDmathFit::PolynomialInterpolation(Double_t xa[], Double_t ya[], Int_t n, Double_t x, Double_t* y, Double_t* dy)
1286 {
1287   // Code from "Numerical Recipes in C" with slight modifications.
1288   // In the following, n must be <= 2*fRegularisation!
1289   // Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y and
1290   // an error estimate dy. If P(x) is the polynomial of degree N − 1 such that P(xai) = yai, i =
1291   // 1, . . . , n, then the returned value y = P(x).
1292
1293
1294   Int_t i, m, ns = 1;
1295   Double_t den, dif, dift, ho, hp, w;
1296
1297   dif = fabs(x - xa[1]);
1298
1299   for (Int_t j = 1; j <= n; j++) {
1300     fPolIntC[j] = 0.;
1301     fPolIntD[j] = 0.;
1302   }
1303
1304   for (i = 1; i <= n; i++) {
1305     // Find the index ns of the closest table entry
1306     if ( (dift = fabs(x - xa[i])) < dif) {
1307       ns = i;
1308       dif = dift;
1309     }
1310
1311     // Initialise tableau of c's and d's
1312     fPolIntC[i] = ya[i];
1313     fPolIntD[i] = ya[i];
1314   }
1315
1316   *y = ya[ns--]; // This is the initial approximation to y.
1317   for (m = 1; m < n; m++) { // For each column of the tableau,
1318     for (i = 1; i <= n - m; i++) { // loop over the current c’s and d’s and update them
1319       ho = xa[i] - x;
1320       hp = xa[i + m] - x;
1321       w = fPolIntC[i + 1] - fPolIntD[i];
1322
1323       if ( fabs((den = ho - hp)) < fkDoubleEpsilonLimit) {
1324         // This error can occur only if two input xa’s are (to within roundoff) identical.
1325         printf("Error in PolynomialInterpolation!\n");
1326
1327         return kFALSE;
1328       }
1329
1330       den = w / den;
1331
1332       // Update c's and d's
1333       fPolIntD[i] = hp * den;
1334       fPolIntC[i] = ho * den;
1335     }
1336
1337     // After each column in the tableau is completed, we decide which correction, c or d,
1338     // we want to add to our accumulating value of y, i.e., which path to take through the
1339     // tableau—forking up or down. We do this in such a way as to take the most “straight
1340     // line” route through the tableau to its apex, updating ns accordingly to keep track of
1341     // where we are. This route keeps the partial approximations centered (insofar as possible)
1342     // on the target x. The last dy added is thus the error indication.
1343     *y += (*dy = (2 * ns < (n-m) ? fPolIntC[ns + 1] : fPolIntD[ns--]));
1344
1345   }
1346
1347   return kTRUE;
1348 }