PWGPP/TPC/macros/PIDCalib/checkShapeEtaTree.C

   1 #include "TTree.h"
   2 #include "TH2D.h"
   3 #include "TH1D.h"
   4 #include "TH3D.h"
   5 #include "TH3F.h"
   6 #include "THnSparse.h"
   7 #include "TProfile.h"
   8 #include "TProfile2D.h"
   9 #include "TF1.h"
  10 #include "TFitResultPtr.h"
  11 #include "TFitResult.h"
  12 #include "TCanvas.h"
  13 #include "TStyle.h"
  14 #include "TVirtualFitter.h"
  15 #include "TLinearFitter.h"
  16 #include "TList.h"
  17 #include "TString.h"
  18 #include "TLegend.h"
  19 #include "TLegendEntry.h"
  20 #include "TFile.h"
  21 #include "TGraphErrors.h"
  22 #include "TGraph.h"
  23 #include "TMath.h"
  24 #include "TDatime.h"
  25 #include "TSystem.h"
  26 #include "TSpline.h"
  27
  28 #include "AliPID.h"
  29
  30 #include <iostream>
  31 #include <iomanip>
  32
  33 #include "THnSparseDefinitions.h"
  34 #include "ProgressBar.h"
  35
  36 enum mapType {
  37   kNoNormalisation = 1,
  38   kSmallEtaNormalisation = 2,
  39   kLargeEtaNormalisation = 3
  40 };
  41 enum type {
  42   kTypeDelta = 1,
  43   kTypeDeltaPrime = 2,
  44   kTypeSigmaDeltaPrime = 3
  45 };
  46
  47 const Double_t binContentThreshold = 1e-4;
  48 const Double_t massProton = AliPID::ParticleMass(AliPID::kProton);
  49 const Double_t massPion = AliPID::ParticleMass(AliPID::kPion);
  50 const TString suffixMapType[kLargeEtaNormalisation + 1] = {"", "NoNormalisation", "SmallEtaNormalisation", "LargeEtaNormalisation"};
  51
  52 TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", 0.6, 1.6);
  53
  54
  55 //TODO NOW getBin... functions needed? Only, if correction (momentum) necessary
  56 //_________________________________________________________________________________________
  57 Int_t getBinX(TH2D* h, Double_t tanTheta)
  58 {
  59   Int_t binX = h->GetXaxis()->FindBin(tanTheta);
  60
  61   if (binX <= 0)
  62     binX = 1;
  63   if (binX > h->GetXaxis()->GetNbins())
  64     binX = h->GetXaxis()->GetNbins();
  65
  66     return binX;
  67 }
  68
  69
  70 //_________________________________________________________________________________________
  71 Int_t getBinY(TH2D* h, Double_t dEdxInv)
  72 {
  73   Int_t binY = h->GetYaxis()->FindBin(dEdxInv);
  74
  75   if (binY <= 0)
  76     binY = 1;
  77   if (binY > h->GetYaxis()->GetNbins())
  78     binY = h->GetYaxis()->GetNbins();
  79
  80   return binY;
  81 }
  82
  83
  84 //__________________________________________________________________________________________
  85 Bool_t FindFitRange(TH1* h, Double_t& lowThreshold, Double_t& highThreshold, Double_t fractionForRange = 0.1)
  86 {
  87     lowThreshold = 0.6;
  88     highThreshold = 1.6;
  89
  90     if (!h)
  91       return kFALSE;
  92
  93     // Average around maximum bin -> Might catch some outliers
  94     Int_t maxBin = h->GetMaximumBin();
  95     Double_t maxVal = h->GetBinContent(maxBin);
  96     UChar_t usedBins = 1;
  97     if (maxBin > 1) {
  98       maxVal += h->GetBinContent(maxBin - 1);
  99       usedBins++;
 100     }
 101     if (maxBin < h->GetNbinsX()) {
 102       maxVal += h->GetBinContent(maxBin + 1);
 103       usedBins++;
 104     }
 105     maxVal /= usedBins;
 106
 107     Double_t thresholdFraction = fractionForRange * maxVal;
 108     Int_t lowThrBin = TMath::Max(1, h->FindFirstBinAbove(thresholdFraction));
 109     Int_t highThrBin = TMath::Min(h->GetNbinsX(), h->FindLastBinAbove(thresholdFraction));
 110
 111     lowThreshold = h->GetBinCenter(lowThrBin);
 112     highThreshold = h->GetBinCenter(highThrBin);
 113
 114     return kTRUE;
 115 }
 116
 117
 118 //__________________________________________________________________________________________
 119 void SetHistAxis(TH2D* h, Int_t type)
 120 {
 121   h->GetXaxis()->SetTitle("tan(#Theta)");
 122   h->GetYaxis()->SetTitle("1/(dE/dx) (arb. unit)");
 123   if (type == kTypeDelta)
 124     h->GetZaxis()->SetTitle("#Delta (arb. unit)");
 125   else if (type == kTypeDeltaPrime)
 126     h->GetZaxis()->SetTitle("#Delta' (arb. unit)");
 127   else if (type == kTypeSigmaDeltaPrime)
 128     h->GetZaxis()->SetTitle("Par1(#sigma_{rel}(#Delta'))");
 129   else
 130     printf("SetHistAxis: Unknown type %d!\n", type);
 131   h->GetXaxis()->SetTitleSize(0.04);
 132   h->GetXaxis()->SetTitleOffset(2.2);
 133   h->GetXaxis()->SetLabelSize(0.03);
 134   h->GetYaxis()->SetTitleSize(0.04);
 135   h->GetYaxis()->SetTitleOffset(2.6);
 136   h->GetYaxis()->SetLabelSize(0.03);
 137   h->GetZaxis()->SetTitleSize(0.04);
 138   h->GetZaxis()->SetTitleOffset(1.3);
 139 }
 140
 141 //__________________________________________________________________________________________
 142 Double_t getMedianOfNonZeros(Double_t* input, const Int_t dim)
 143 {
 144   Double_t values[dim];
 145   for (Int_t i = 0; i < dim; i++)
 146     values[i] = 0.0;
 147
 148   Int_t numNonZero = 0;
 149
 150   for (Int_t i = 0; i < dim; i++) {
 151     if (input[i] > 0) {
 152       values[numNonZero] = input[i];
 153       numNonZero++;
 154     }
 155   }
 156
 157   return ((numNonZero > 0) ? TMath::Median(numNonZero, values) : 0);
 158 }
 159
 160
 161 //__________________________________________________________________________________________
 162 void normaliseHisto(TH2D* h, Double_t lowerLimit, Double_t upperLimit, Int_t type)
 163 {
 164   if (lowerLimit >= upperLimit) {
 165     printf("!!!Error normaliseHisto: lowerLimit >= upperLimit!\n");
 166     return;
 167   }
 168
 169   if (lowerLimit < 0 || upperLimit < 0) {
 170     printf("!!!Error normaliseHisto: lowerLimit, upperLimit >= 0 required!\n");
 171     return;
 172   }
 173
 174   Int_t binLow = h->GetXaxis()->FindBin(lowerLimit);
 175   if (binLow < 1)
 176     binLow = 1;
 177
 178   Int_t binHigh = h->GetXaxis()->FindBin(upperLimit);
 179   if (binHigh > h->GetXaxis()->GetNbins())
 180     binHigh = h->GetXaxis()->GetNbins();
 181
 182   Int_t binLow2 = h->GetXaxis()->FindBin(-upperLimit);
 183   if (binLow2 < 1)
 184     binLow2 = 1;
 185
 186   Int_t binHigh2 = h->GetXaxis()->FindBin(-lowerLimit);
 187   if (binHigh2 > h->GetXaxis()->GetNbins())
 188     binHigh2 = h->GetXaxis()->GetNbins();
 189
 190   // If 0 is included in range, it might happen that both ranges overlap -> Remove one of the double counted binsPforMap
 191     if (binHigh2 >= binLow) {
 192       binHigh2 = binLow - 1;
 193       if (binHigh2 < 1)   {
 194         printf("!!!Error: binHigh2 = binLow - 1 is < 1\n");
 195         return;
 196       }
 197     }
 198
 199     if (binLow2 > binHigh2)
 200       binLow2 = binHigh2;
 201
 202     // Normalise with respect to some given tan(theta)
 203       // To be robust against fluctuations: Take median of a few tan(theta) bins around the desired value for scaling
 204       const Int_t nThetaBins = (binHigh - binLow + 1) + (binHigh2 - binLow2 + 1);
 205
 206       if (nThetaBins <= 0)  {
 207         printf("Cannot renormalise histo due to bad limits for reference bins: %f -> %f\n", lowerLimit, upperLimit);
 208         return;
 209       }
 210
 211       Double_t* values;
 212       values = new Double_t[nThetaBins];
 213
 214       for (Int_t i = 1; i <= h->GetYaxis()->GetNbins(); i++) {
 215         // Reset values
 216         Int_t binIndex = 0;
 217         for (Int_t thetaBin = 1; thetaBin <= nThetaBins; thetaBin++)
 218           values[thetaBin - 1] = 0;
 219
 220         for (Int_t thetaBin = binLow; thetaBin <= binHigh; thetaBin++)  {
 221           Double_t temp = h->GetBinContent(thetaBin, i);
 222           values[binIndex] = (temp > 0 || type == kTypeDelta) ? temp : 0;
 223           binIndex++;
 224         }
 225         for (Int_t thetaBin = binLow2; thetaBin <= binHigh2; thetaBin++)  {
 226           Double_t temp = h->GetBinContent(thetaBin, i);
 227           values[binIndex] = (temp > 0 || type == kTypeDelta) ? temp : 0;
 228           binIndex++;
 229         }
 230
 231         Double_t temp = 0;
 232         Double_t scale = 0;
 233
 234         if (type == kTypeDeltaPrime)  {
 235           temp = getMedianOfNonZeros(values, nThetaBins);
 236           if (temp <= 0)
 237             continue;
 238           scale = 1. / temp;
 239         }
 240         else if (type == kTypeDelta)  {
 241           scale = TMath::Median(nThetaBins, values);
 242         }
 243         else  {
 244           printf("Error: Type %d not supported for normalisation!\n", type);
 245           return;
 246         }
 247
 248
 249         for (Int_t thetaBin = 1; thetaBin <= h->GetXaxis()->GetNbins(); thetaBin++)  {
 250           if (type == kTypeDeltaPrime)  {
 251             h->SetBinContent(thetaBin, i, h->GetBinContent(thetaBin, i) * scale);
 252             h->SetBinError(thetaBin, i, h->GetBinError(thetaBin, i) * scale);
 253           }
 254           else if (type == kTypeDelta)  {
 255             h->SetBinContent(thetaBin, i, h->GetBinContent(thetaBin, i) - scale);
 256             h->SetBinError(thetaBin, i, h->GetBinError(thetaBin, i) - scale);
 257           }
 258         }
 259       }
 260       delete values;
 261 }
 262
 263
 264 //__________________________________________________________________________________________
 265 void eliminateNonPositivePointsInHisto(TH2D* h)
 266 {
 267   const Int_t nNeighbours = 24;
 268   Double_t* values;
 269   values = new Double_t[nNeighbours];
 270
 271   // Search for bins with content <= 0 (bins with fit errors). Then take all surrounding points in +-2 rows and columns
 272   // and take their median (only those bins without fit errors!) as the new bin content of the considered bin.
 273
 274   for (Int_t binX = 1; binX <= h->GetXaxis()->GetNbins(); binX++) {
 275     Int_t firstBinLeft = TMath::Max(1, binX - 2);
 276     Int_t lastBinRight = TMath::Min(h->GetXaxis()->GetNbins(), binX + 2);
 277
 278     for (Int_t binY = 1; binY <= h->GetYaxis()->GetNbins(); binY++) {
 279       if (h->GetBinContent(binX, binY) <= 0) {
 280         Int_t firstBinBelow = TMath::Max(1, binY - 2);
 281         Int_t lastBinAbove = TMath::Min(h->GetYaxis()->GetNbins(), binY + 2);
 282
 283         // Reset values
 284         Int_t binIndex = 0;
 285         for (Int_t i = 0; i < nNeighbours; i++)
 286           values[i] = 0;
 287
 288         for (Int_t binX2 = firstBinLeft; binX2 <= lastBinRight; binX2++) {
 289           for (Int_t binY2 = firstBinBelow; binY2 <= lastBinAbove; binY2++) {
 290             if (binX2 == binX && binY2 == binY)
 291               continue; // skip bin that is currently under consideration
 292
 293             // Only take values from (hopefully) valid fits
 294             if (h->GetBinContent(binX2, binY2) > 0) {
 295               values[binIndex] = h->GetBinContent(binX2, binY2);
 296               binIndex++;
 297             }
 298           }
 299         }
 300
 301         Double_t temp = getMedianOfNonZeros(values, nNeighbours);
 302         if (temp <= 0) {
 303           // Only print error message, if there is at least one positive value
 304           if (binIndex > 0)
 305             printf("Error: Could not eliminate values <= 0 for bin at (%f, %f)!\n",
 306                   h->GetXaxis()->GetBinCenter(binX), h->GetYaxis()->GetBinCenter(binY));
 307           temp = -1;
 308         }
 309
 310         printf("Eliminated non-positive value at bin (%f, %f): %f -> %f!\n",
 311                h->GetXaxis()->GetBinCenter(binX), h->GetYaxis()->GetBinCenter(binY),
 312                h->GetBinContent(binX, binY), temp);
 313         h->SetBinContent(binX, binY, temp);
 314         h->SetBinError(binX, binY, 1000);
 315       }
 316     }
 317   }
 318
 319   delete values;
 320 }
 321
 322
 323 //__________________________________________________________________________________________
 324 void eliminateOutliers(TH2D* h, Double_t threshold)
 325 {
 326   const Int_t nNeighbours = 8;
 327   Double_t* values;
 328   values = new Double_t[nNeighbours];
 329
 330   // Search for outliers (most likely bad fits). Then take all surrounding points in +-1 rows and columns
 331   // and take their median (only those bins without (known) fit errors, i.e. bin content > 0).
 332   // If the current bin content deviates by more than "threshold" from the median, take the median as the new bin content
 333   // of the considered bin.
 334
 335   for (Int_t binY = 1; binY <= h->GetYaxis()->GetNbins(); binY++) {
 336     Int_t firstBinBelow = TMath::Max(1, binY - 1);
 337     Int_t lastBinAbove = TMath::Min(h->GetYaxis()->GetNbins(), binY + 1);
 338
 339     for (Int_t binX = 1; binX <= h->GetXaxis()->GetNbins(); binX++) {
 340       Int_t firstBinLeft = TMath::Max(1, binX - 1);
 341       Int_t lastBinRight = TMath::Min(h->GetXaxis()->GetNbins(), binX + 1);
 342
 343       Double_t binContent = h->GetBinContent(binX, binY);
 344
 345
 346
 347       // Reset values
 348       Int_t binIndex = 0;
 349       for (Int_t i = 0; i < nNeighbours; i++)
 350         values[i] = 0;
 351
 352       for (Int_t binX2 = firstBinLeft; binX2 <= lastBinRight; binX2++) {
 353         for (Int_t binY2 = firstBinBelow; binY2 <= lastBinAbove; binY2++) {
 354           if (binX2 == binX && binY2 == binY)
 355             continue; // skip bin that is currently under consideration
 356
 357           // Only take values from (hopefully) valid fits
 358           if (h->GetBinContent(binX2, binY2) > 0) {
 359             values[binIndex] = h->GetBinContent(binX2, binY2);
 360             binIndex++;
 361           }
 362         }
 363       }
 364
 365       Double_t temp = getMedianOfNonZeros(values, nNeighbours);
 366       if (temp <= 0) {
 367         // Only print error message, if there is at least one positive value
 368         if (binIndex > 0)
 369           printf("Error: Could not eliminate outlier for bin at (%f, %f)!\n",
 370                  h->GetXaxis()->GetBinCenter(binX), h->GetYaxis()->GetBinCenter(binY));
 371         temp = -1;
 372       }
 373       else {
 374         if (TMath::Abs(binContent - temp) > threshold) {
 375           printf("Eliminated outlier at bin (%f, %f): %f -> %f!\n",
 376                h->GetXaxis()->GetBinCenter(binX), h->GetYaxis()->GetBinCenter(binY),
 377                binContent, temp);
 378           h->SetBinContent(binX, binY, temp);
 379           h->SetBinError(binX, binY, 1000);
 380         }
 381       }
 382     }
 383   }
 384
 385   delete values;
 386 }
 387
 388
 389 //__________________________________________________________________________________________
 390 void addPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
 391 {
 392   if (h->GetBinContent(binX, binY) <= binContentThreshold)
 393     return; // Reject bins without content (within some numerical precision) or with strange content
 394
 395   Double_t coord[2] = {0, 0};
 396   coord[0] = h->GetXaxis()->GetBinCenter(binX);
 397   coord[1] = h->GetYaxis()->GetBinCenter(binY);
 398   Double_t binError = h->GetBinError(binX, binY);
 399   if (binError <= 0) {
 400     printf("Warning: Trying to add bin in addPointToHyperplane with error (%e) not set (%f, %f), but bin content %f. Maybe fit problem -> Setting large error\n",
 401            binError, coord[0], coord[1], h->GetBinContent(binX, binY));
 402     binError = 1000;
 403   }
 404   linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
 405 }
 406
 407 //__________________________________________________________________________________________
 408 TH2D* refineHistoViaLinearExtrapolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY, Int_t mapType,
 409                                         TFile* fSave, Bool_t sigmaMap = kFALSE, Bool_t skipBinsAtBoundaries = kFALSE)
 410 {
 411   if (!h)
 412     return 0x0;
 413
 414   if (h->GetEntries() <= 0)
 415     return 0x0;
 416
 417   Int_t nBinsX = h->GetXaxis()->GetNbins();
 418   Int_t nBinsY = h->GetYaxis()->GetNbins();
 419
 420   Int_t firstYbin = 1;
 421
 422   if (skipBinsAtBoundaries) {
 423     // Remove the two first and the last bin on y-axis
 424     nBinsY -= 3;
 425     firstYbin = 3;
 426   }
 427
 428   // Extrapolate map to larger 1/dEdx using the last 3 points (taking into account skipBinsAtBoundaries). Then: Refine
 429
 430   Int_t nBinsYextrapolated = nBinsY + TMath::Ceil((1. / 20. - h->GetYaxis()->GetBinUpEdge(nBinsY)) / h->GetYaxis()->GetBinWidth(nBinsY));
 431   // Make sure that the "old" bins stay the same and only some additional bins are added (i.e. bin width should NOT change)
 432   Double_t yUpBoundExtrapolated =  h->GetYaxis()->GetBinUpEdge(nBinsY) +  h->GetYaxis()->GetBinWidth(nBinsY) * (nBinsYextrapolated - nBinsY);
 433
 434   TH2D* hExtrapolated = new TH2D(Form("%s_%s_extrapolated", h->GetName(), suffixMapType[mapType].Data()), Form("%s (refined)",h->GetTitle()),
 435                                  nBinsX, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(nBinsX),
 436                                  nBinsYextrapolated, h->GetYaxis()->GetBinLowEdge(1), yUpBoundExtrapolated);
 437
 438   // Copy the old bins and take into account skipBinsAtBoundaries
 439   // Get rid of most non-positive bin contents
 440   eliminateNonPositivePointsInHisto(h);
 441
 442   for (Int_t binX = 1; binX <= nBinsX; binX++)  {
 443     for (Int_t binY = 1; binY <= nBinsY; binY++)    {
 444       if (h->GetBinContent(binX, binY) > 0 && h->GetBinError(binX, binY) > 0)   {
 445         hExtrapolated->SetBinContent(binX, binY, h->GetBinContent(binX, binY));
 446         hExtrapolated->SetBinError(binX, binY, h->GetBinError(binX, binY));
 447       }
 448     }
 449     // Default value for extrapolated points is 1
 450     for (Int_t binY = nBinsY + 1; binY <= nBinsYextrapolated; binY++)    {
 451       hExtrapolated->SetBinContent(binX, binY, 1);
 452       hExtrapolated->SetBinError(binX, binY, 1000);
 453     }
 454   }
 455
 456
 457   // Do the extrapolation
 458
 459   TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
 460
 461   for (Int_t binX = 1; binX <= nBinsX; binX++)  {
 462     linExtrapolation->ClearPoints();
 463
 464     // Last 3 points for fitting for eta map,
 465     // but only last 2 points for sigma map, since the resolution
 466     // usually levels off around MIP and 3 points used for fitting
 467     // would cause an underestimation of the resolution.
 468
 469     for (Int_t binY = (sigmaMap ? (nBinsY - 1) : (nBinsY - 2)); binY <= nBinsY; binY++)   {
 470       Double_t centerX = hExtrapolated->GetXaxis()->GetBinCenter(binX);
 471       Double_t centerY = hExtrapolated->GetYaxis()->GetBinCenter(binY);
 472
 473       Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
 474       if (oldBinX < 1)
 475         oldBinX = 1;
 476       if (oldBinX > nBinsX)
 477         oldBinX = nBinsX;
 478
 479       Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
 480       if (oldBinY < firstYbin)
 481         oldBinY = firstYbin;
 482       if (oldBinY > nBinsY)
 483         oldBinY = nBinsY;
 484
 485       // Neighbours left column
 486         if (oldBinX >= 2) {
 487           if (oldBinY >= 2) {
 488             addPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
 489           }
 490
 491           addPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
 492
 493           if (oldBinY < nBinsY) {
 494             addPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
 495           }
 496         }
 497
 498         // Neighbours (and point itself) same column
 499         if (oldBinY >= 2) {
 500           addPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
 501         }
 502
 503         addPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
 504
 505         if (oldBinY < nBinsY) {
 506           addPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
 507         }
 508
 509         // Neighbours right column
 510         if (oldBinX < nBinsX) {
 511           if (oldBinY >= 2) {
 512             addPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
 513           }
 514
 515           addPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
 516
 517           if (oldBinY < nBinsY) {
 518             addPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
 519           }
 520         }
 521     }
 522     // Fit 2D-hyperplane
 523     if (linExtrapolation->GetNpoints() <= 0)
 524       continue;
 525
 526     if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
 527         continue;
 528
 529     // ....extrapolation to the rest
 530     for (Int_t binY = nBinsY + 1; binY <= nBinsYextrapolated; binY++)    {
 531       Double_t centerX = hExtrapolated->GetXaxis()->GetBinCenter(binX);
 532       Double_t centerY = hExtrapolated->GetYaxis()->GetBinCenter(binY);
 533
 534       // Fill the bin of the refined histogram with the extrapolated value
 535       Double_t extrapolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
 536       + linExtrapolation->GetParameter(2) * centerY;
 537
 538       // Do not allow too small or even negative values
 539       extrapolatedValue = TMath::Max(binContentThreshold * 2, extrapolatedValue);
 540       hExtrapolated->SetBinContent(binX, binY, extrapolatedValue);
 541       hExtrapolated->SetBinError(binX, binY, 50); // High, but constant, value => should almost not influence not-extrapolated data
 542     }
 543   }
 544
 545   eliminateNonPositivePointsInHisto(hExtrapolated);
 546
 547   // Normalise map on demand
 548
 549   // Use kNoNormalisation for final QA
 550   if (mapType == kSmallEtaNormalisation) {
 551     normaliseHisto(hExtrapolated, 0., 0.11, kTypeDeltaPrime);
 552   }
 553   else if (mapType == kLargeEtaNormalisation) {
 554     normaliseHisto(hExtrapolated, 0.81, 0.99, kTypeDeltaPrime);
 555   }
 556
 557   // Interpolate to finer map
 558   Int_t nBinsXrefined = nBinsX * refineFactorX;
 559   Int_t nBinsYrefined = nBinsYextrapolated * refineFactorY;
 560
 561   TH2D* hRefined = new TH2D(Form("%s_%s_refined", h->GetName(), suffixMapType[mapType].Data()),  Form("%s (refined)",h->GetTitle()),
 562                             nBinsXrefined, hExtrapolated->GetXaxis()->GetBinLowEdge(1), hExtrapolated->GetXaxis()->GetBinUpEdge(nBinsX),
 563                             nBinsYrefined, hExtrapolated->GetYaxis()->GetBinLowEdge(firstYbin),
 564                             hExtrapolated->GetYaxis()->GetBinUpEdge(nBinsYextrapolated));
 565
 566   for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
 567     Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
 568
 569     for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
 570
 571       hRefined->SetBinContent(binX, binY, 1); // Default value is 1
 572
 573       Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
 574
 575       /*OLD
 576       linExtrapolation->ClearPoints();
 577
 578       // For interpolation: Just take the corresponding bin from the old histo.
 579       // For extrapolation: take the last available bin from the old histo.
 580       // If the boundaries are to be skipped, also skip the corresponding bins
 581       Int_t oldBinX = hExtrapolated->GetXaxis()->FindBin(centerX);
 582       if (oldBinX < 1)
 583         oldBinX = 1;
 584       if (oldBinX > nBinsX)
 585         oldBinX = nBinsX;
 586
 587       Int_t oldBinY = hExtrapolated->GetYaxis()->FindBin(centerY);
 588       if (oldBinY < firstYbin)
 589         oldBinY = firstYbin;
 590       if (oldBinY > nBinsYextrapolated)
 591         oldBinY = nBinsYextrapolated;
 592
 593       // Neighbours left column
 594       if (oldBinX >= 2) {
 595         if (oldBinY >= 2) {
 596           addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX - 1, oldBinY - 1);
 597         }
 598
 599         addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX - 1, oldBinY);
 600
 601         if (oldBinY < nBinsYextrapolated) {
 602           addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX - 1, oldBinY + 1);
 603         }
 604       }
 605
 606       // Neighbours (and point itself) same column
 607       if (oldBinY >= 2) {
 608         addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX, oldBinY - 1);
 609       }
 610
 611       addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX, oldBinY);
 612
 613       if (oldBinY < nBinsYextrapolated) {
 614         addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX, oldBinY + 1);
 615       }
 616
 617       // Neighbours right column
 618       if (oldBinX < nBinsX) {
 619         if (oldBinY >= 2) {
 620           addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX + 1, oldBinY - 1);
 621         }
 622
 623         addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX + 1, oldBinY);
 624
 625         if (oldBinY < nBinsYextrapolated) {
 626           addPointToHyperplane(hExtrapolated, linExtrapolation, oldBinX + 1, oldBinY + 1);
 627         }
 628       }
 629
 630
 631       // Fit 2D-hyperplane
 632       if (linExtrapolation->GetNpoints() <= 0)
 633         continue;
 634
 635       if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
 636         continue;
 637
 638       // Fill the bin of the refined histogram with the extrapolated value
 639       Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
 640                                  + linExtrapolation->GetParameter(2) * centerY;
 641       */
 642       Double_t interpolatedValue = hExtrapolated->Interpolate(centerX, centerY);
 643       hRefined->SetBinContent(binX, binY, interpolatedValue);
 644     }
 645   }
 646
 647
 648   // Problem: Interpolation does not work before/beyond center of first/last bin (as the name suggests).
 649   // Therefore, for each row in dEdx: Take last bin from old map and interpolate values from center and edge.
 650   // Assume line through these points and extropolate to last bin of refined map
 651   const Double_t firstOldXbinUpEdge = hExtrapolated->GetXaxis()->GetBinUpEdge(1);
 652   const Double_t firstOldXbinCenter = hExtrapolated->GetXaxis()->GetBinCenter(1);
 653
 654   const Double_t oldXbinHalfWidth = firstOldXbinUpEdge - firstOldXbinCenter;
 655
 656   const Double_t lastOldXbinLowEdge = hExtrapolated->GetXaxis()->GetBinLowEdge(hExtrapolated->GetNbinsX());
 657   const Double_t lastOldXbinCenter = hExtrapolated->GetXaxis()->GetBinCenter(hExtrapolated->GetNbinsX());
 658
 659   for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
 660     Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
 661
 662     const Double_t interpolatedCenterFirstXbin = hExtrapolated->Interpolate(firstOldXbinCenter, centerY);
 663     const Double_t interpolatedUpEdgeFirstXbin = hExtrapolated->Interpolate(firstOldXbinUpEdge, centerY);
 664
 665     const Double_t extrapolationSlopeFirstXbin = (interpolatedUpEdgeFirstXbin - interpolatedCenterFirstXbin) / oldXbinHalfWidth;
 666     const Double_t extrapolationOffsetFirstXbin = interpolatedCenterFirstXbin;
 667
 668
 669     const Double_t interpolatedCenterLastXbin = hExtrapolated->Interpolate(lastOldXbinCenter, centerY);
 670     const Double_t interpolatedLowEdgeLastXbin = hExtrapolated->Interpolate(lastOldXbinLowEdge, centerY);
 671
 672     const Double_t extrapolationSlopeLastXbin = (interpolatedCenterLastXbin - interpolatedLowEdgeLastXbin) / oldXbinHalfWidth;
 673     const Double_t extrapolationOffsetLastXbin = interpolatedCenterLastXbin;
 674
 675     for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
 676       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
 677
 678       if (centerX < firstOldXbinCenter) {
 679         Double_t extrapolatedValue = extrapolationOffsetFirstXbin + (centerX - firstOldXbinCenter) * extrapolationSlopeFirstXbin;
 680         hRefined->SetBinContent(binX, binY, extrapolatedValue);
 681       }
 682       else if (centerX <= lastOldXbinCenter) {
 683         continue;
 684       }
 685       else {
 686         Double_t extrapolatedValue = extrapolationOffsetLastXbin + (centerX - lastOldXbinCenter) * extrapolationSlopeLastXbin;
 687         hRefined->SetBinContent(binX, binY, extrapolatedValue);
 688       }
 689     }
 690   }
 691
 692   delete linExtrapolation;
 693
 694
 695   if (fSave)    {
 696     fSave->cd();
 697     hExtrapolated->Write();
 698   }
 699
 700   delete hExtrapolated;
 701
 702   return hRefined;
 703 }
 704
 705
 706 //__________________________________________________________________________________________
 707 TFitResult*  doubleGaussFit(TH1D* h, Double_t currentMeanMomentum, TSpline3* splPr, TSpline3* splPi, TString fitOption = "QNS")
 708 {
 709   Double_t lowThreshold = 0.6;
 710   Double_t highThreshold = 1.6;
 711   FindFitRange(h, lowThreshold, highThreshold);
 712   TFitResultPtr result = h->Fit("gaus", "QNS", "", lowThreshold, highThreshold);
 713
 714   Double_t contaminationPeakMean = splPi->Eval(currentMeanMomentum / massPion) / splPr->Eval(currentMeanMomentum / massProton);
 715
 716   if (contaminationPeakMean < h->GetXaxis()->GetBinLowEdge(1) ||
 717       contaminationPeakMean > h->GetXaxis()->GetBinUpEdge(h->GetNbinsX())) {
 718     return ((Int_t)result == 0) ? (TFitResult*)result.Get()->Clone() : 0x0;
 719   }
 720
 721   // Estimate parameters for doubleGauss fit
 722   Double_t estimatedMean = 0;
 723   Double_t estimatedSigma = 0;
 724   Double_t estimatedYield = 0;
 725
 726   Double_t chi2oneGauss = 0;
 727   Int_t NDFoneGauss = 0;
 728   Double_t reducedChi2oneGauss = 0;
 729
 730   if ((Int_t) result == 0) {
 731     estimatedMean = result->GetParams()[1];
 732     estimatedSigma = result->GetParams()[2];
 733     estimatedYield = result->GetParams()[0];
 734
 735     chi2oneGauss = result->Chi2();
 736     NDFoneGauss = result->Ndf();
 737     reducedChi2oneGauss = (NDFoneGauss > 0) ? chi2oneGauss / NDFoneGauss : -1;
 738   }
 739   else {
 740     estimatedMean = h->GetMean();
 741     estimatedSigma = h->GetRMS();
 742     estimatedYield = (estimatedSigma > 0) ? (h->Integral("width") / estimatedSigma) : h->GetEntries();
 743   }
 744
 745   TF1* doubleGaus = new TF1("doubleGaus", "[0]*TMath::Gaus(x,[1],[2],0)+[3]*TMath::Gaus(x,[4],[2],0)", 0.6, 1.6);
 746
 747   estimatedMean = TMath::Max(0.6, estimatedMean);
 748   estimatedYield = TMath::Max(1., estimatedYield);
 749
 750   Double_t newPars[5] = { estimatedYield, estimatedMean, estimatedSigma, estimatedYield / 10., contaminationPeakMean };
 751   doubleGaus->SetParameters(newPars);
 752   doubleGaus->SetParLimits(0, 0.01 * estimatedYield, 100. * estimatedYield);//TODO 0., 100. * estimatedYield);
 753   doubleGaus->SetParLimits(1, estimatedMean - 0.05, estimatedMean + 0.05);//TODO 0.6, 1.6);
 754   doubleGaus->SetParLimits(2, 0, 100. * estimatedSigma);
 755   doubleGaus->SetParLimits(3, 0, 0.8 * estimatedYield);
 756   doubleGaus->SetParLimits(4, contaminationPeakMean - 0.15, contaminationPeakMean + 0.15);//TODO doubleGaus->SetParLimits(4, 0.6, 1.6);
 757   TFitResultPtr result2 = h->Fit(doubleGaus, fitOption.Data());
 758
 759   printf("\n%f -> %f\n", currentMeanMomentum, contaminationPeakMean);//TODO NOW NOW NOW
 760   printf("%f, %f, %f, %f, %f\n%f, %f, %f\n", result2->GetParams()[0], result2->GetParams()[1], result2->GetParams()[2], result2->GetParams()[3], result2->GetParams()[4], result->GetParams()[0], result->GetParams()[1], result->GetParams()[2]); printedSomething = kTRUE;//TODO NOW NOW NOW
 761   if ((Int_t)result2 == 0) {
 762     Double_t chi2doubleGauss = 0;
 763     Int_t NDFdoubleGauss = 0;
 764     Double_t reducedChi2doubleGauss = 0;
 765
 766     chi2doubleGauss = result2->Chi2();
 767     NDFdoubleGauss = result2->Ndf();
 768     reducedChi2doubleGauss = (NDFdoubleGauss > 0) ? chi2doubleGauss / NDFdoubleGauss : -1;
 769
 770     printf("%f / %d = %f\n %f / %d = %f\n\n", chi2oneGauss, NDFoneGauss, reducedChi2oneGauss, chi2doubleGauss, NDFdoubleGauss, reducedChi2doubleGauss); printedSomething = kTRUE;//TODO NOW NOW NOW
 771     // Only take double gauss result, if (valid) reduced chi2 is better than that of the one gauss fit
 772     if (reducedChi2oneGauss < 0 || (reducedChi2doubleGauss > 0 && reducedChi2doubleGauss <= reducedChi2oneGauss))
 773       return (TFitResult*)result2.Get()->Clone();
 774   }
 775
 776   // If fit failed, return results of standard fit instead
 777   return ((Int_t)result == 0) ? (TFitResult*)result.Get()->Clone() : 0x0;
 778 }
 779
 780
 781 //__________________________________________________________________________________________
 782 TFitResult* extractClusterDependence(TH3D** hPreMapClusterResolved, Int_t numClusterBins, Int_t clusterLowBound, Int_t clusterUpBound,
 783                                      Int_t binX, Int_t binY, Int_t binYhigh, Double_t c0, TFile* fSave, TSpline3* splPr, TSpline3* splPi)
 784 {
 785   TH1D* hClusMean = new TH1D(Form("hClusMean_tanTheta_%f_pTPC_%f", hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
 786                                   hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY)),
 787                              Form("#Delta' vs. ncl (tan(#Theta) %f, p_{TPC} %f);ncl;#Delta'", hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
 788                                   hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY)),
 789                              numClusterBins, clusterLowBound, clusterUpBound);
 790
 791   TH1D* hClusSigmaRelGauss = new TH1D(Form("hClusSigmaRelGauss_tanTheta_%f_pTPC_%f", hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
 792                                            hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY)),
 793                                       Form("#sigma_{rel, Gauss} vs. ncl (tan(#Theta) %f, p_{TPC} %f);ncl;#sigma_{rel, Gauss}",
 794                                            hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
 795                                            hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY)),
 796                                       numClusterBins, clusterLowBound, clusterUpBound);
 797
 798   TH1D* hClusSigmaRel = new TH1D(Form("hClusSigmaRel_tanTheta_%f_pTPC_%f", hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
 799                                            hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY)),
 800                                       Form("#sigma_{rel} vs. ncl (tan(#Theta) %f, p_{TPC} %f);ncl;#sigma_{rel}",
 801                                            hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
 802                                            hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY)),
 803                                       numClusterBins, clusterLowBound, clusterUpBound);
 804
 805   /*//TODO ADDED 03.07.2013
 806   TH1D* hTest = hPreMapClusterResolved[0]->ProjectionZ("hTest", binX, binX, binY, binYhigh);
 807   Long64_t nEntries = hTest->GetEntries();
 808   Double_t clusterMean = (clusterLowBound + (clusterUpBound - clusterLowBound) / numClusterBins * 0) * nEntries;
 809
 810   for (Int_t clusterBin = 1; clusterBin < numClusterBins; clusterBin++) {
 811     TH1D* hTempProjectionZ = hPreMapClusterResolved[clusterBin]->ProjectionZ("hTempProjectionZ", binX, binX, binY, binYhigh);
 812
 813     if (hTempProjectionZ->GetEntries() < 10) {
 814       delete hTempProjectionZ;
 815       continue;
 816     }
 817
 818     hTest->Add(hTempProjectionZ);
 819     nEntries += hTempProjectionZ->GetEntries();
 820     clusterMean += (clusterLowBound + (clusterUpBound - clusterLowBound) / numClusterBins * clusterBin) * hTempProjectionZ->GetEntries();
 821     delete hTempProjectionZ;
 822   }
 823   if (nEntries > 0)
 824     clusterMean /= nEntries;
 825
 826   Double_t lowThreshold = 0.6;
 827   Double_t highThreshold = 1.6;
 828   FindFitRange(hTest, lowThreshold, highThreshold);
 829   TFitResultPtr fitResPtr3 = hTest->Fit("gaus", "QNS", "", lowThreshold, highThreshold);
 830   TFitResult* fitRes3 = ((Int_t)fitResPtr3 == 0) ? (TFitResult*)fitResPtr3.Get()->Clone() : 0x0;
 831
 832   if (!fitRes3 || fitRes3->GetParams()[1] <= 0) {
 833     printf("ERROR: Fit test failed\n");
 834     return 0x0;
 835   }
 836
 837   Double_t mean = fitRes3->GetParams()[1];
 838   Double_t sigma = fitRes3->GetParams()[2];
 839
 840   Double_t relSigma = sigma / mean;
 841
 842   TGraphErrors * grClusSigmaRel2 = new TGraphErrors(1);
 843   grClusSigmaRel2->SetPoint(0, clusterMean, relSigma);
 844
 845   TF1* fStats2 = new TF1("fStats2", "TMath::Sqrt([0] * [0] + [1] * [1] / x)", clusterLowBound, clusterUpBound);
 846   fStats2->SetLineColor(kBlue);
 847   fStats2->FixParameter(0, c0);
 848   //fStats2->SetParameter(1, 0.1);
 849   //fStats2->SetParLimits(1, 0, 2);
 850   fStats2->SetNpx(300);
 851
 852   TFitResultPtr fitRes14 = grClusSigmaRel2->Fit(fStats2, "QS", "", clusterLowBound, clusterUpBound);
 853
 854   delete fStats2;
 855   delete grClusSigmaRel2;
 856
 857   delete hTest;
 858
 859   if (((Int_t)fitRes14) != 0) {
 860     printf("\nError: Cluster fit: Dependence fit failed (%d) for (%f, %f - %f): skipped\n", (Int_t)fitRes14,
 861            hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX), hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY),
 862            hPreMapClusterResolved[0]->GetYaxis()->GetBinUpEdge(binYhigh));
 863     printedSomething = kTRUE;
 864     return 0x0;
 865   }
 866
 867   return fitRes14.Get();
 868   //TODO END*/
 869
 870   Double_t minLowThreshold = 0.6;
 871   Double_t maxHighThreshold = 1.6;
 872   fGauss.SetRange(minLowThreshold, maxHighThreshold);
 873
 874   for (Int_t clusterBin = 0; clusterBin < numClusterBins; clusterBin++) {
 875     TH1D* hTempProjectionZ = hPreMapClusterResolved[clusterBin]->ProjectionZ("hTempProjectionZ", binX, binX, binY, binYhigh);
 876
 877     if (hTempProjectionZ->GetEntries() < 10) {
 878       /*
 879       printf("\nWarning: Cluster fit: No entries for (%f, %f - %f, ncl %f): skipped\n",
 880              hPreMapClusterResolved[clusterBin]->GetXaxis()->GetBinCenter(binX), hPreMapClusterResolved[clusterBin]->GetYaxis()->GetBinLowEdge(binY),
 881              hPreMapClusterResolved[clusterBin]->GetYaxis()->GetBinUpEdge(binYhigh), hClusMean->GetXaxis()->GetBinCenter(clusterBin + 1));
 882       printedSomething = kTRUE;
 883
 884       */
 885       delete hTempProjectionZ;
 886       continue;
 887     }
 888
 889     TFitResult* fitRes = 0x0;
 890
 891     Double_t mean = -999.;
 892     Double_t errorMean = 999.;
 893
 894     Double_t sigma = 999.;
 895     Double_t errorSigma = 999.;
 896
 897     if (splPr && splPi) { // do DoubleGaussFit
 898       fitRes = doubleGaussFit(hTempProjectionZ,
 899                               0.5 * (hPreMapClusterResolved[0]->GetYaxis()->GetBinCenter(binY) +
 900                                      hPreMapClusterResolved[0]->GetYaxis()->GetBinCenter(binYhigh)),
 901                               splPr, splPi, "QNS");
 902
 903       if (fitRes) {
 904         mean = fitRes->GetParams()[1];
 905         errorMean = fitRes->GetErrors()[1];
 906
 907         sigma = fitRes->GetParams()[2];
 908         errorSigma = fitRes->GetErrors()[2];
 909       }
 910     }
 911     else {
 912       Double_t lowThreshold = minLowThreshold;
 913       Double_t highThreshold = maxHighThreshold;
 914
 915       // First, fit with the desired restriction via FindFitRange to fix the mean,
 916       // but to fit the width, fit the full range (otherwise one typically gets biased
 917       // to larger sigmas). However, fitting the full range would also result in a different
 918       // mean. Thus, to be consistent, these two steps are needed.
 919       // NOTE: If the restricted fit fails, just take the mean of the fit over the full range (if it doesn't fail),
 920       // since this is still better than no information
 921
 922       FindFitRange(hTempProjectionZ, lowThreshold, highThreshold);
 923       TFitResultPtr fitResPtrFirstStep = hTempProjectionZ->Fit("gaus", "QNS", "", lowThreshold, highThreshold);
 924
 925       if ((Int_t)fitResPtrFirstStep == 0) {
 926         fGauss.SetParameter(0, fitResPtrFirstStep->GetParams()[0]);
 927         fGauss.SetParError(0, fitResPtrFirstStep->GetErrors()[0]);
 928         fGauss.SetParameter(2, fitResPtrFirstStep->GetParams()[2]);
 929         fGauss.SetParError(2, fitResPtrFirstStep->GetErrors()[2]);
 930
 931         fGauss.FixParameter(1, fitResPtrFirstStep->GetParams()[1]);
 932         TFitResultPtr fitResPtrSecondStep = hTempProjectionZ->Fit(&fGauss, "QNS", "", minLowThreshold, maxHighThreshold);
 933
 934         fitRes = ((Int_t)fitResPtrSecondStep == 0) ? (TFitResult*)fitResPtrSecondStep.Get()->Clone() : 0x0;
 935
 936         if (fitRes) {
 937           // Mean and error from first step (error is zero in second step because mean is fixed there)
 938           mean = fitResPtrFirstStep->GetParams()[1];
 939           errorMean = fitResPtrFirstStep->GetErrors()[1];
 940
 941           sigma = fitRes->GetParams()[2];
 942           errorSigma = fitRes->GetErrors()[2];
 943         }
 944       }
 945       else {
 946         TFitResultPtr fitResPtrSecondStep = hTempProjectionZ->Fit("gaus", "QNS", "", minLowThreshold, maxHighThreshold);
 947
 948         fitRes = ((Int_t)fitResPtrSecondStep == 0) ? (TFitResult*)fitResPtrSecondStep.Get()->Clone() : 0x0;
 949
 950         if (fitRes) {
 951           mean = fitRes->GetParams()[1];
 952           errorMean = fitRes->GetErrors()[1];
 953
 954           sigma = fitRes->GetParams()[2];
 955           errorSigma = fitRes->GetErrors()[2];
 956         }
 957       }
 958     }
 959
 960     if (!fitRes) {
 961       printf("\nWarning: Cluster fit: Gauss fit failed for (%f, %f - %f, ncl %f): skipped\n",
 962              hPreMapClusterResolved[clusterBin]->GetXaxis()->GetBinCenter(binX), hPreMapClusterResolved[clusterBin]->GetYaxis()->GetBinLowEdge(binY),
 963              hPreMapClusterResolved[clusterBin]->GetYaxis()->GetBinUpEdge(binYhigh), hClusMean->GetXaxis()->GetBinCenter(clusterBin + 1));
 964       printedSomething = kTRUE;
 965
 966       delete hTempProjectionZ;
 967       continue;
 968     }
 969     else {
 970       hClusMean->SetBinContent(clusterBin + 1, mean);
 971       hClusMean->SetBinError(clusterBin + 1, errorMean);
 972
 973       hClusSigmaRelGauss->SetBinContent(clusterBin + 1, sigma);
 974       hClusSigmaRelGauss->SetBinError(clusterBin + 1, errorSigma);
 975
 976       if (mean > 1e-4)   {
 977         hClusSigmaRel->SetBinContent(clusterBin + 1, sigma / mean);
 978         hClusSigmaRel->SetBinError(clusterBin + 1,
 979                                    TMath::Sqrt(TMath::Power(errorSigma / mean, 2) + TMath::Power(errorMean * sigma / (mean * mean), 2)));
 980       }
 981       else {
 982         printf("\nWarning: Cluster fit: Gauss fit with strange mean value (%e) for (%f, %f - %f, ncl %f): skipped\n", mean,
 983              hPreMapClusterResolved[clusterBin]->GetXaxis()->GetBinCenter(binX), hPreMapClusterResolved[clusterBin]->GetYaxis()->GetBinLowEdge(binY),
 984              hPreMapClusterResolved[clusterBin]->GetYaxis()->GetBinUpEdge(binYhigh), hClusMean->GetXaxis()->GetBinCenter(clusterBin + 1));
 985         printedSomething = kTRUE;
 986       }
 987
 988       delete hTempProjectionZ;
 989     }
 990   }
 991
 992   TF1* fStats = new TF1("fStats", "TMath::Sqrt([0] * [0] + [1] * [1] / x)", clusterLowBound, clusterUpBound);
 993   fStats->SetLineColor(kBlue);
 994   fStats->FixParameter(0, c0);
 995   //fStats->SetParameter(1, 0.1);
 996   //fStats->SetParLimits(1, 0, 2);
 997   fStats->SetNpx(300);
 998
 999   TGraphErrors * grClusSigmaRel = new TGraphErrors(hClusSigmaRel);
1000   for(Int_t ip = 0; ip < grClusSigmaRel->GetN(); ip ++) {
1001     Bool_t removePoint = grClusSigmaRel->GetY()[ip] <= 0 || grClusSigmaRel->GetEY()[ip] / TMath::Max(1e-8, grClusSigmaRel->GetY()[ip]) > 0.20; //TODO NOW was > 0.10
1002                          //|| grClusSigmaRel->GetX()[ip] > 140; // TODO (was 160) For some reason it is better to skip the last clusters
1003     if (removePoint) {
1004       grClusSigmaRel->RemovePoint(ip);
1005       ip--;
1006       continue;
1007     }
1008   }
1009   grClusSigmaRel->SetMarkerStyle(20);
1010   grClusSigmaRel->SetMarkerColor(kMagenta);
1011
1012   if (grClusSigmaRel->GetN() <= 0) {
1013     printf("\nError: Cluster fit: No good data points for dependence fit for (%f, %f - %f): skipped\n",
1014            hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX), hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY),
1015            hPreMapClusterResolved[0]->GetYaxis()->GetBinUpEdge(binYhigh));
1016     printedSomething = kTRUE;
1017     return 0x0;
1018   }
1019
1020   TFitResultPtr fitRes = grClusSigmaRel->Fit(fStats, "QS", "", clusterLowBound, clusterUpBound);
1021   //printf("\nchi^2/NDF = %f / %d = %f\n\n", fitRes->Chi2(), fitRes->Ndf(), (fitRes->Ndf() > 0) ? fitRes->Chi2() / fitRes->Ndf() : -1);
1022
1023   TCanvas* canvCluster = new TCanvas(Form("canvCluster_tanTheta_%f_pTPC_%f_%f",
1024                                           hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
1025                                           hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY),
1026                                           hPreMapClusterResolved[0]->GetYaxis()->GetBinUpEdge(binYhigh)),
1027                                      Form("canvCluster_tanTheta_%f_pTPC_%f_%f",
1028                                           hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX),
1029                                           hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY),
1030                                           hPreMapClusterResolved[0]->GetYaxis()->GetBinUpEdge(binYhigh)),
1031                                      100,10,1380,800);
1032   canvCluster->Divide(3,1);
1033
1034   canvCluster->cd(1);
1035   hClusMean->DrawCopy("e");
1036
1037   canvCluster->cd(2);
1038   hClusSigmaRelGauss->DrawCopy("e");
1039
1040   canvCluster->cd(3);
1041   hClusSigmaRel->DrawCopy("e");
1042   grClusSigmaRel->Draw("same p");
1043
1044   if (fSave)  {
1045     fSave->cd("clusterQA");
1046     canvCluster->Write();
1047     //hClusMean->Write();
1048     //hClusSigmaRelGauss->Write();
1049     //hClusSigmaRel->Write();
1050   }
1051
1052   delete fStats;
1053   delete grClusSigmaRel;
1054
1055   delete hClusMean;
1056   delete hClusSigmaRel;
1057   delete hClusSigmaRelGauss;
1058
1059   delete canvCluster;
1060
1061   if (((Int_t)fitRes) != 0) {
1062     printf("\nError: Cluster fit: Dependence fit failed (%d) for (%f, %f - %f): skipped\n", (Int_t)fitRes,
1063            hPreMapClusterResolved[0]->GetXaxis()->GetBinCenter(binX), hPreMapClusterResolved[0]->GetYaxis()->GetBinLowEdge(binY),
1064            hPreMapClusterResolved[0]->GetYaxis()->GetBinUpEdge(binYhigh));
1065     printedSomething = kTRUE;
1066     return 0x0;
1067   }
1068
1069   return fitRes.Get();
1070 }
1071
1072 //__________________________________________________________________________________________
1073 // NOTE: Use smoothing only for creation of sigma map. For normal eta map this pulls down the edges too much <-> trade-off between
1074 // (normally very small) jumps in the map and good description at the edges.
1075 // For the sigma map, there are not that sharp edges, but the fluctuations are by far stronger. Here it makes sense to use the smoothing.
1076 TH2* checkShapeEtaTree(TTree* tree = 0x0, Double_t resolutionFactorX = 1/*2 or 2.1*/, Double_t resolutionFactorY = 1/*2 or 2.1*/,
1077                        Double_t pThresholdTOFcut = -1/*0.6*/, Int_t nMergeBinsAroundThreshold = 3, Double_t pThresholdV0cut = 1.2,
1078                        Double_t pThresholdV0plusTOFcut = 2.0, Double_t outlierThreshold = 0.04,
1079                        Bool_t doSmoothing = kFALSE, Double_t c0 = 0.0232,  Int_t nclCut = 60,
1080                        TString path = ".", TString suffixForFileName = "", Bool_t useDoubleGaussFit = kFALSE,
1081                        TString pathNameSplinesFile = "TPCPIDResponse.root", TString prSplinesName = "TSPLINE3_DATA_PROTON_LHC10H_PASS2_PBPB_MEAN",
1082                        Int_t returnMapType = kNoNormalisation, TString treeName = "fTree")
1083 {
1084   TString fileName = Form("bhess_PIDetaTree%s.root", suffixForFileName.Data());
1085
1086   if (resolutionFactorX <= 0 || resolutionFactorY <= 0)  {
1087           printf("Factor for lower resolution <= 0 is not allowed!\n");
1088           return 0x0;
1089   }
1090
1091   TFile* f = 0x0;
1092
1093   Bool_t saveResults = kFALSE;
1094
1095   if (!tree)  {
1096     f = TFile::Open(Form("%s/%s", path.Data(), fileName.Data()));
1097     if (!f)  {
1098       std::cout << "Failed to open file \"" << Form("%s/%s", path.Data(), fileName.Data()) << "\"!" << std::endl;
1099
1100
1101       return 0x0;
1102     }
1103
1104     // Extract the data Tree
1105     tree = dynamic_cast<TTree*>(f->Get(treeName.Data()));
1106     if (!tree) {
1107       std::cout << "Failed to load data tree!" << std::endl;
1108       return 0x0;
1109     }
1110
1111     saveResults = kTRUE;
1112   }
1113
1114
1115
1116   //TODO NOW start
1117   /*
1118   //TString pathNameMomentumMap = "pionTree_11b2/outputCheckShapeEtaTreePions_2012_10_16__11_38.root";
1119   //TString pathNameMomentumMap = "pionTree_10d1/outputCheckShapeEtaTreePions_2012_10_17__07_38.root";
1120   TString pathNameMomentumMap = "finalCuts/uncorrected/0.6GeVcut/pp/7TeV/LHC10d.pass2/correctedV0splines_newV0tree/outputCheckShapeEtaTreePions_2012_10_23__16_02_correctedWith_10_38.root";
1121   TFile* fPiMap = TFile::Open(pathNameMomentumMap.Data());
1122   if (!fPiMap)  {
1123     std::cout << "Failed to open pion momentum map file \"" << pathNameMomentumMap.Data() << "\"!" << std::endl;
1124     return 0x0;
1125   }
1126
1127   TH2D* hPiMap = 0x0;
1128
1129   if (fPiMap) {
1130     hPiMap = dynamic_cast<TH2D*>(fPiMap->Get("hRefinedSmallEtaNormalisation")); //TODO NOW NOW NOWForm(hRefined%s", suffixMapType[returnMapType].Data())));
1131     if (!hPiMap) {
1132       std::cout << "Failed to load pion momentum map!" << std::endl;
1133       return 0x0;
1134     }
1135   }*/
1136   //TODO NOW end
1137
1138
1139
1140
1141   TSpline3* splPr = 0x0;
1142   TSpline3* splPi = 0x0;
1143
1144   if (useDoubleGaussFit) {
1145     std::cout << "Using double gauss fit!" << std::endl << std::endl;
1146
1147     std::cout << "Loading splines from file \"" << pathNameSplinesFile.Data() << "\"" << std::endl;
1148
1149     TFile* fSpl = TFile::Open(pathNameSplinesFile.Data());
1150     if (!fSpl) {
1151       std::cout << "Failed to open spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
1152       return 0x0;
1153     }
1154
1155     TObjArray* TPCPIDResponse = (TObjArray*)fSpl->Get("TPCPIDResponse");
1156     if (!TPCPIDResponse) {
1157       std::cout << "Failed to load object array from spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
1158       return 0x0;
1159     }
1160
1161     splPr = (TSpline3*)TPCPIDResponse->FindObject(prSplinesName.Data());
1162     TString piSplinesName = prSplinesName.ReplaceAll("PROTON", "PION");
1163     splPi = (TSpline3*)TPCPIDResponse->FindObject(piSplinesName.Data());
1164
1165     if (!splPr || !splPi) {
1166       std::cout << "Failed to load splines from file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
1167       return 0x0;
1168     }
1169   }
1170   else
1171     std::cout << "NOT using double gauss fit!" << std::endl << std::endl;
1172
1173
1174   // Output file
1175
1176   TFile* fSave = 0x0;
1177   TString saveFileName = "";
1178   if (saveResults)  {
1179     TDatime daTime;
1180
1181     TString saveSuffix = "";
1182     if (suffixForFileName.Length() > 0)
1183       saveSuffix = Form("_%s", suffixForFileName.Data());
1184
1185     saveFileName = Form("outputCheckShapeEtaTree_%04d_%02d_%02d__%02d_%02d%s.root", daTime.GetYear(), daTime.GetMonth(),
1186                         daTime.GetDay(), daTime.GetHour(), daTime.GetMinute(), saveSuffix.Data());
1187
1188     fSave = TFile::Open(Form("%s/%s", path.Data(), saveFileName.Data()), "recreate");
1189     if (!fSave) {
1190       std::cout << "Failed to open save file \"" << Form("%s/%s", path.Data(), saveFileName.Data()) << "\"!" << std::endl;
1191       return 0x0;
1192     }
1193
1194     fSave->mkdir("clusterQA");
1195   }
1196
1197   printf("\nProcessing file %s\n", f->GetName());
1198
1199   // Only activate the branches of interest to save processing time
1200   tree->SetBranchStatus("*", 0); // Disable all branches
1201   tree->SetBranchStatus("pTPC", 1);
1202   //tree->SetBranchStatus("pT", 1);
1203   //tree->SetBranchStatus("phiPrime", 1);
1204   tree->SetBranchStatus("dEdx", 1);
1205   tree->SetBranchStatus("dEdxExpected", 1);
1206   tree->SetBranchStatus("tanTheta", 1);
1207   //if (isPbPb)
1208   //  tree->SetBranchStatus("fMultiplicity", 1);
1209   //tree->SetBranchStatus("sinAlpha", 1);
1210   tree->SetBranchStatus("tpcSignalN", 1);
1211   tree->SetBranchStatus("pidType", 1);
1212
1213   TString multiplicitySelectionString = "";
1214 // if (isPbPb)// Usually 20 bins from 0 to 4000 => bin width ~ 200
1215 //    multiplicitySelectionString = "&& fMultiplicity > 3000 && fMultiplicity <= 3200";
1216
1217   // Bins along tanTheta
1218   const Int_t numTanThetaBins = 30; //WARNING: Also change in AliPIDResponse, if changed here!
1219   Int_t binsX = numTanThetaBins;
1220
1221   if (resolutionFactorX != 1)
1222     binsX /= resolutionFactorX;
1223
1224   // Do not mix tan(theta) > 0 with tan(theta) < 0!
1225   // Since the eta range is symmetric around zero (i.e. |lowerTanTheta| = |upperTanTheta|),
1226   // this means that the number of bins for tan(theta) < 0 must equal that of tan(theta) > 0.
1227   // Thus, in order to have a bin edge at tan(theta) = 0, binsX must be even.
1228   if (binsX % 2 != 0) {
1229     binsX++;
1230     // Set resolutionFactorX correspondingly to the really used value
1231     resolutionFactorX = ((Double_t)numTanThetaBins) / ((Double_t)binsX);
1232     printf("\nNOTE: In order not to mix up tan(theta) < 0 with > 0, an even number of tanTheta-bins is required => ResolutionFactorX set to %f!\n\n",
1233            resolutionFactorX);
1234   }
1235
1236   // ---------------------------------------------------------------------------------
1237   // Use momentum bins and take expected dEdx to translate from momentum to dEdx in map
1238
1239   const Int_t clusterUpBound = 160;
1240   //TODO was 70 and 10; changed on 24.06.13
1241   const Int_t clusterLowBound = 60; // Other possibility: 80. Fits look better, but obviously the results for the PID-task are worse
1242   const Int_t clusterBinWidth = 10; // Other possibility: 20. Fits look better, but obviously the results for the PID-task are worse
1243   // NOTE: Too large cluster bins are bad because you fit a superposition of Gauss curves with different widths. This potentially
1244   // biases your fit to larger values. The effect for the resulting dependence 1/sqrt(ncl) is typically few permille for deltaNcl = 20
1245   // compared to deltaNcl = 10
1246   const Int_t numClusterBins = (clusterUpBound - clusterLowBound) / clusterBinWidth;
1247
1248   const Double_t pBoundLow = 0.15;
1249   const Double_t pBoundUp = 5.0;
1250   const Double_t lowerTanTheta = -1.0;
1251   const Double_t upperTanTheta = 1.0;
1252
1253
1254
1255
1256   const Int_t nPbinsForMap = 110;
1257   Double_t binsPforMap[nPbinsForMap + 1];
1258
1259   const Double_t factorBinning = TMath::Power(pBoundUp/pBoundLow, 1./nPbinsForMap);
1260
1261   // Log binning
1262   binsPforMap[0] = pBoundLow;
1263   for (Int_t i = 0 + 1; i <= nPbinsForMap; i++) {
1264     binsPforMap[i] = factorBinning * binsPforMap[i - 1];
1265   }
1266   binsPforMap[nPbinsForMap] = pBoundUp;
1267
1268   /*TODO changed on 23.06.13
1269   // 68 bins from 0.15 to 1.0 (width 0.0125)
1270   // 20 bins from 1.0 to 1.5 (width 0.025)
1271   // 10 bins from 1.5 to 2.0 (width 0.05)
1272   //  5 bins from 2.0 to 2.5 (width 0.1)
1273   //  4 bins from 2.5 to 3.5 (width 0.25)
1274   //  3 bins from 3.5 to 5.0 (width 0.5)
1275   const Int_t nPbinsForMap = 68 + 20 + 10 + 5 + 4 + 3;
1276   Double_t binsPforMap[nPbinsForMap + 1];
1277   for (Int_t i = 0; i < 68; i++)  {
1278     binsPforMap[i] = pBoundLow + i * 0.0125;
1279   }
1280   for (Int_t i = 68, j = 0; i < 68 + 20; i++, j++)  {
1281     binsPforMap[i] = 1.0 + j * 0.025;
1282   }
1283   for (Int_t i = 68 + 20, j = 0; i < 68 + 20 + 10; i++, j++)  {
1284     binsPforMap[i] = 1.5 + j * 0.05;
1285   }
1286   for (Int_t i = 68 + 20 + 10, j = 0; i < 68 + 20 + 10 + 5; i++, j++)  {
1287     binsPforMap[i] = 2.0 + j * 0.1;
1288   }
1289   for (Int_t i = 68 + 20 + 10 + 5, j = 0; i < 68 + 20 + 10 + 5 + 4; i++, j++)  {
1290     binsPforMap[i] = 2.5 + j * 0.25;
1291   }
1292   for (Int_t i = 68 + 20 + 10 + 5 + 4, j = 0; i < 68 + 20 + 10 + 5 + 4 + 3; i++, j++)  {
1293     binsPforMap[i] = 3.5 + j * 0.5;
1294   }
1295   binsPforMap[nPbinsForMap] = pBoundUp;
1296   */
1297
1298   progressbar(0.);
1299
1300   const Int_t nDim = 3;
1301
1302   Int_t thnBins[nDim] = {  binsX, nPbinsForMap, 200  };
1303   Double_t xmin[nDim] = {  lowerTanTheta, pBoundLow, 0.6  };
1304   Double_t xmax[nDim] = {  upperTanTheta, pBoundUp, 1.6  };
1305   THnSparseD* hDummy = new THnSparseD("hDummy", "", nDim, thnBins, xmin, xmax);
1306   hDummy->SetBinEdges(1, binsPforMap);
1307   TH3D* hPreMap = hDummy->Projection(0, 1, 2);
1308   hPreMap->SetName("hPreMap");
1309   hPreMap->SetTitle("#Delta' map for p_{TPC_Inner} vs. tan(#theta)");
1310   hPreMap->GetXaxis()->SetTitle("tan(#theta)");
1311   hPreMap->GetYaxis()->SetTitle("p_{TPC_inner} (GeV/c)");
1312   hPreMap->GetZaxis()->SetTitle("#Delta'");
1313
1314   delete hDummy;
1315
1316   // Use the expected dEdx to translate from pTPC to 1/dEdx -> Create a 2D histo with the correct momentum binning for this translation
1317   const Int_t nDim2 = 2;
1318
1319   Int_t thnBins2[nDim2] = {  nPbinsForMap, 200000  };
1320   Double_t xmin2[nDim2] = {  pBoundLow, 0.  };
1321   Double_t xmax2[nDim2] = {  pBoundUp, 2000.  };
1322
1323   hDummy = new THnSparseD("hDummy", "", nDim2, thnBins2, xmin2, xmax2);
1324   hDummy->SetBinEdges(0, binsPforMap);
1325   TH2D* hMomTranslation = hDummy->Projection(1, 0);
1326   hMomTranslation->SetName("hMomTranslation");
1327   hMomTranslation->SetTitle("Momentum to dEdxExpected translation map)");
1328   hMomTranslation->GetXaxis()->SetTitle("p_{TPC_inner} (GeV/c)");
1329   hMomTranslation->GetYaxis()->SetTitle("expected(dE/dx) (arb. unit)");
1330
1331   delete hDummy;
1332
1333   // Maps for different number of clusters - used to extract dependence of ncl and create corresponding map for paramatrisation
1334   TH3D* hPreMapClusterResolved[numClusterBins];
1335
1336   hDummy = new THnSparseD("hDummy", "", nDim, thnBins, xmin, xmax);
1337   hDummy->SetBinEdges(1, binsPforMap);
1338
1339   for (Int_t numClusters = clusterLowBound, clusterBin = 0; numClusters + clusterBinWidth <= clusterUpBound;
1340        numClusters += clusterBinWidth, clusterBin++)   {
1341     TH3D* hTemp = 0x0;
1342     hTemp = hDummy->Projection(0, 1, 2);
1343     hTemp->SetName("hTemp");
1344     hTemp->SetTitle(Form("#Delta' map for p_{TPC_Inner} vs. tan(#theta) for ncl around %d", numClusters));
1345     hTemp->GetXaxis()->SetTitle("tan(#theta)");
1346     hTemp->GetYaxis()->SetTitle("p_{TPC_inner} (GeV/c)");
1347     hTemp->GetZaxis()->SetTitle("#Delta'");
1348
1349     hTemp->SetName(Form("hPreMapClusterResolved_%d", clusterBin));
1350     hPreMapClusterResolved[clusterBin] = hTemp;
1351   }
1352
1353   delete hDummy;
1354
1355
1356
1357   Long64_t nTreeEntries = tree->GetEntriesFast();
1358
1359   Double_t dEdx = 0.; // Measured dE/dx
1360   Double_t dEdxExpected = 0.; // Expected dE/dx according to parametrisation
1361   Double_t tanTheta = 0.; // Tangens of (local) theta at TPC inner wall
1362   Double_t pTPC = 0.; // Momentum at TPC inner wall
1363   UShort_t tpcSignalN = 0; // Number of clusters used for dEdx
1364   //Int_t fMultiplicity = 0;
1365   UChar_t pidType = 0;
1366   //Double_t phiPrime = 0;
1367
1368
1369   tree->SetBranchAddress("dEdx", &dEdx);
1370   tree->SetBranchAddress("dEdxExpected", &dEdxExpected);
1371   tree->SetBranchAddress("tanTheta", &tanTheta);
1372   tree->SetBranchAddress("tpcSignalN", &tpcSignalN);
1373   tree->SetBranchAddress("pTPC", &pTPC);
1374   tree->SetBranchAddress("pidType", &pidType);
1375   //tree->SetBranchAddress("phiPrime", &phiPrime);
1376   //if (isPbPb)
1377   //  tree->SetBranchAddress("fMultiplicity", &fMultiplicity);
1378
1379
1380   /*TODO NOW
1381   TF1* fShapeSmallP = new TF1("fShapeSmallP", "pol5", -0.4, 0.4);
1382   fShapeSmallP->SetParameters(1.01712, -0.0202725, -0.260692, 0.261623, 0.671854, -1.14014);
1383   */
1384   for (Long64_t i = 0; i < nTreeEntries; i++) {
1385     tree->GetEntry(i);
1386
1387     if (dEdx <= 0 || dEdxExpected <= 0)
1388       continue;
1389
1390     if (nclCut >= 0 && tpcSignalN <= nclCut)
1391       continue;
1392
1393
1394     //if (pidType == kV0idNoTOF || pidType == kV0idPlusTOFaccepted || pidType == kV0idPlusTOFrejected) continue; //For testing to get the old result
1395
1396     // Use different PID techniques depending on the momentum - except for MC
1397     if (pidType != kMCid) {
1398       if ((pTPC < 0.6 && pidType != kTPCid)               // No V0s/TOF below 0.6 due to bad statistics
1399         //TODO NOW Maybe don't use this line since the tails get fewer V0's than the main peak, so the tails are overall reduced || (pTPC >= 0.6 && pTPC < pThresholdV0cut && pidType != kTPCandTOFid) || // Do not mix V0's and TOF in this range
1400           || (pTPC >= pThresholdV0cut && pTPC < pThresholdV0plusTOFcut &&
1401               pidType != kV0idNoTOF && pidType != kV0idPlusTOFaccepted) // Only V0's above pThresholdV0cut due to contamination
1402           || (pTPC >= pThresholdV0plusTOFcut && pidType != kV0idPlusTOFaccepted)) // Additionally require TOF
1403         continue;
1404     }
1405
1406     //Double_t pT = pTPC*TMath::Sin(-TMath::ATan(tanTheta)+TMath::Pi()/2.0);
1407     //if ((phiPrime > 0.072/pT+TMath::Pi()/18.0-0.035 && phiPrime < 0.07/pT/pT+0.1/pT+TMath::Pi()/18.0+0.035))
1408     //  continue;
1409
1410     /*TODO
1411     if (TMath::Abs(tanTheta) <= 0.4) {
1412       Double_t p0 = fShapeSmallP->Eval(tanTheta) - 1.0; // Strength of the correction
1413       Double_t p1 = -9.0; // How fast the correction is turned off
1414       Double_t p2 = -0.209; // Turn off correction around 0.2 GeV/c
1415       Double_t p3 = 1.0; // Delta' for large p should be 1
1416
1417       Double_t corrFactor = TMath::Erf((pTPC + p2) * p1) * p0 + p3 + p0; // Add p0 to have 1 for p3 = 1 and large pTPC
1418       dEdxExpected *= corrFactor;
1419     }*/
1420
1421     // TODO Old unsuccessful try
1422     //Double_t thetaGlobalTPC = -TMath::ATan(tanTheta) + TMath::Pi() / 2.;
1423     //Double_t pTtpc = pTPC * TMath::Sin(thetaGlobalTPC);
1424     //Double_t pTtpcInv = (pTtpc > 0) ? 1. / pTtpc : 0;
1425     //Double_t p0 = 1.0;
1426     //Double_t p1 = 1./ 0.5;//TODO 2.0;
1427     //Double_t p2 = -0.05;//TODO 0.1
1428     //Double_t pTcorrFactor = p0 + (pTtpcInv > p1) * p2 * (pTtpcInv - p1);
1429     //
1430     //dEdxExpected *= pTcorrFactor;
1431
1432
1433     // Correct for pure momentum dependence (azimuthal dependence comes into play)
1434     /*TODO NOW
1435     Double_t correctionFactor = hPiMap->GetBinContent(getBinX(hPiMap, tanTheta), getBinY(hPiMap, pTPC));
1436     correctionFactor -= 1.;
1437     correctionFactor *= 0.5*(TMath::Erf((0.5 - pTPC) / 0.05) + 1.); // Only correct up to 0.5 GeV/c and switch off within about 2*0.05 GeV/c
1438     //correctionFactor *=  0.5* splPr->Eval(pTPC / massProton) / splPi->Eval(pTPC / massPion);//TODO NOW
1439     correctionFactor += 1.;
1440     dEdx /= correctionFactor;*/
1441
1442
1443     // TODO Maybe needed(?) - maybe not for the maps, since this should be valid for all particles!: "acceptedByPhiPrimeCut == 1"
1444     //if (!(TMath::Abs(tpcSignalN - 150) <= 5 && (phiPrime < 0.072/pT+pi/18.0-0.035 || phiPrime > 0.07/pT/pT+0.1/pT+pi/18.0+0.035)))
1445     //  continue
1446
1447     hPreMap->Fill(tanTheta, pTPC, dEdx / dEdxExpected);
1448     hMomTranslation->Fill(pTPC, dEdxExpected);
1449
1450     for (Int_t numClusters = clusterLowBound, clusterBin = 0; numClusters + clusterBinWidth <= clusterUpBound;
1451          numClusters += clusterBinWidth, clusterBin++)   {
1452       if (tpcSignalN > numClusters && tpcSignalN <= numClusters + clusterBinWidth) {
1453         hPreMapClusterResolved[clusterBin]->Fill(tanTheta, pTPC, dEdx / dEdxExpected);
1454         break;
1455       }
1456     }
1457
1458     if (i % 1000 == 0)
1459       progressbar(100. * (((Double_t)i) / nTreeEntries));
1460   }
1461
1462   progressbar(100.);
1463
1464
1465   // Go to MIP region (around p = 3 GeV/c, tan(theta) small) and take the mean value as the end point of the map
1466   // Since the contamination with kaons already starts, it is safer to go only to up to p = 2.7 GeV/c. Or: Use V0s
1467   printf("\nDetermining upper and lower 1/dEdx bound of map...\n\n");
1468   ///*TODO NOW NOW
1469   tree->Project("hAdjustMIP(200000,0,100)", "dEdxExpected",
1470                 Form("pTPC > 1.5 && dEdx > 0 && dEdxExpected > 0 %s", multiplicitySelectionString.Data()));
1471
1472   TH1D* hAdjustMIP = (TH1D*)gDirectory->Get("hAdjustMIP");
1473
1474   Int_t binMIP = hAdjustMIP->FindFirstBinAbove(0.0);
1475
1476   Double_t tempMean = (binMIP > 0) ? hAdjustMIP->GetBinCenter(binMIP) : -1.;
1477   if (tempMean <= 0)  {
1478     printf("ERROR: Failed to determine expected dEdx in MIP region!\n");
1479     return 0x0;
1480   }
1481   Double_t upperMapBound = 1. / tempMean;
1482
1483   TProfile* hMomTranslationProfile = hMomTranslation->ProfileX();
1484         hMomTranslationProfile->GetXaxis()->SetRangeUser(0.7, 4.5);
1485   Int_t mipBin = hMomTranslationProfile->GetMinimumBin();
1486   Double_t mipMomentum = hMomTranslationProfile->GetXaxis()->GetBinUpEdge(mipBin);
1487
1488   printf("Upper 1/dEdx bound of map set to: 1./%f = %f (pTPC ~ %.2f GeV/c)\n", tempMean, upperMapBound, mipMomentum);
1489
1490   delete hAdjustMIP;
1491   delete hMomTranslationProfile;
1492
1493   // Take pre-map and loop through tanTheta-rows in momentum bins. Take the first row with sufficient statistics as
1494   // the lower bound and translate to the corresponding 1./dEdx
1495   TH2D* hPreMap_p_tanTheta = (TH2D*)hPreMap->Project3D("yx");
1496   Int_t lowMomentumBinWithSufficientStatistics = 1;
1497   for (Int_t binY = 1; binY <= hPreMap_p_tanTheta->GetNbinsY(); binY++) {
1498     Bool_t sufficientStatistics = kTRUE;
1499     for (Int_t binX = 1; binX <= hPreMap_p_tanTheta->GetNbinsX(); binX++) {
1500       // TODO Try to estimate the statistics with resolutionFactorY
1501       if (hPreMap_p_tanTheta->GetBinContent(binX, binY) < 100 / resolutionFactorY) {
1502         sufficientStatistics = kFALSE;
1503         break;
1504       }
1505     }
1506     if (sufficientStatistics) {
1507       lowMomentumBinWithSufficientStatistics = binY;
1508       break;
1509     }
1510   }
1511   TH1D* hDeDxExpectedLowPbound = hMomTranslation->ProjectionY("_pyLow", lowMomentumBinWithSufficientStatistics, lowMomentumBinWithSufficientStatistics);
1512
1513   Double_t meanDeDxExpectedLowPbound = hDeDxExpectedLowPbound->GetMean();
1514   delete hDeDxExpectedLowPbound;
1515
1516   if (meanDeDxExpectedLowPbound <= 0)  {
1517     printf("\nCould not determine mean dEdx_expected for lower p bound (pTPC = %f)!\n",
1518            hPreMap_p_tanTheta->GetYaxis()->GetBinCenter(lowMomentumBinWithSufficientStatistics));
1519     return 0x0;
1520   }
1521   Double_t lowerMapBound = 1. / meanDeDxExpectedLowPbound;
1522   printf("Lower 1/dEdx bound of map set to: 1./%f = %f (pTPC ~ %.2f GeV/c)\n", meanDeDxExpectedLowPbound, lowerMapBound,
1523          hPreMap_p_tanTheta->GetYaxis()->GetBinCenter(lowMomentumBinWithSufficientStatistics));
1524   //*/
1525
1526   /*
1527   printf("TODO: Fixed upper map bound to 0.02!!\n");//TODO NOW NOW
1528   Double_t upperMapBound = 0.02;//TODO NOW NOW
1529   printf("TODO: Fixed lower map bound to 0.002!!\n");//TODO NOW NOW
1530   Double_t lowerMapBound = 0.002;//TODO NOW NOW
1531   */
1532   //Double_t lowerMapBound = 0.0016;//TODO NOW NOW Should be used for PbPb due to bad splines
1533   //printf("\nTODO NOW: lowerMapBound fixed to 0.0016 - only needed for bad splines in PbPb!!!\n\n");//TODO NOW NOW
1534
1535
1536   // Binning was found to yield good results, if 40 bins are chosen for the range 0.0016 to 0.02. For the new variable range,
1537   // scale the number of bins correspondingly
1538
1539   Int_t binsY = TMath::Nint((upperMapBound - lowerMapBound) / (0.02 - 0.0016) * 40);//WARNING: Also change in AliPIDResponse, if changed here!
1540
1541   if (resolutionFactorY != 1)
1542     binsY /= resolutionFactorY;
1543
1544   printf("Used number of y bins: %d\n", binsY);
1545
1546   printf("\n\nExtracting histograms from tree. This may take some time...\n");
1547
1548   TH2D* hRes3DprimeFit = new TH2D("hRes3DprimeFit",
1549                                   "#Delta' map for 1./(dE/dx) vs. tan(#theta) via Gauss fit;tan(#theta);1./(dE/dx) (arb. unit);#Delta'",
1550                                   binsX, lowerTanTheta, upperTanTheta, binsY, lowerMapBound, upperMapBound);
1551   TH2D* hSigmaPar1 = (TH2D*)hRes3DprimeFit->Clone("hSigmaPar1");
1552   hSigmaPar1->SetTitle("Parameter 1 of #sigma_{rel}(#Delta') map for 1./(dE/dx) vs. tan(#theta) via Gauss fit");
1553   hSigmaPar1->GetZaxis()->SetTitle("Par1(#sigma_{rel}(#Delta'))");
1554
1555   TH2D* hRes3DprimeProfile = (TH2D*)hRes3DprimeFit->Clone("hRes3DprimeProfile");
1556   hRes3DprimeProfile->SetTitle("#Delta' map for 1./(dE/dx) vs. tan(#theta) via mean");
1557
1558
1559
1560   printf("\n\nStarted map creation....\n");
1561
1562
1563   const Int_t numTotalBins = hPreMap->GetXaxis()->GetNbins() * hPreMap->GetYaxis()->GetNbins();
1564   Bool_t specialHandlingApplied = kFALSE;
1565
1566   // If special handling for the threshold of the TOF cut is requested: Determine the p-Bin which contains the threshold
1567   Int_t binWithThreshold = -1;
1568   if (pThresholdTOFcut > 0) {
1569     for (Int_t binY = 1; binY <= hPreMap->GetYaxis()->GetNbins(); binY++)  {
1570       if (hPreMap->GetYaxis()->GetBinLowEdge(binY) <= pThresholdTOFcut &&
1571           hPreMap->GetYaxis()->GetBinUpEdge(binY) > pThresholdTOFcut)  {
1572         binWithThreshold = binY;
1573         break;
1574       }
1575     }
1576   }
1577
1578   for (Int_t binY = 1; binY <= hPreMap->GetYaxis()->GetNbins(); binY++) {
1579      // TODO added on 24.06.13
1580     // Do not go beyond MIP momentum
1581     if (hPreMap->GetYaxis()->GetBinLowEdge(binY) >= mipMomentum) {
1582       printf("\nReached momentum of MIPs. Stopping...\n");
1583       printedSomething = kTRUE;
1584       break;
1585     }
1586
1587     // Use the expected dEdx to translate from pTPC to 1/dEdx
1588     TH1D* hDeDxExpected = hMomTranslation->ProjectionY("_py", binY, binY);
1589
1590     Double_t meanDeDxExpected = hDeDxExpected->GetMean();
1591     delete hDeDxExpected;
1592
1593     if (meanDeDxExpected <= 0)  {
1594       printf("\nCould not determine mean dEdx_expected for (p = %f) - skipping!\n", hPreMap->GetYaxis()->GetBinCenter(binY));
1595       printedSomething = kTRUE;
1596       continue;
1597     }
1598
1599     Int_t dEdxBin = hRes3DprimeFit->GetYaxis()->FindBin(1. / meanDeDxExpected);
1600     if (dEdxBin < 1 || dEdxBin > hRes3DprimeFit->GetYaxis()->GetNbins())  {
1601       printf("\nSkipped bin (out of map range): p = %f GeV/c, 1./<dEdxExpected> = %f\n", hPreMap->GetYaxis()->GetBinCenter(binY),
1602              1. / meanDeDxExpected);
1603       printedSomething = kTRUE;
1604       continue;
1605     }
1606
1607     // If there are subsequent bins with same 1/dEdxExpected-bin, merge them to increase statistics.
1608     Int_t nMergeBins = 1;
1609     for (; ; nMergeBins++) {
1610       if (binY + nMergeBins > hPreMap->GetYaxis()->GetNbins())  {
1611         nMergeBins--;
1612         break;
1613       }
1614
1615       // TODO added on 24.06.13
1616       // Do not go beyond MIP momentum
1617       if (hPreMap->GetYaxis()->GetBinLowEdge(binY + nMergeBins) >= mipMomentum) {
1618         nMergeBins--;
1619         break;
1620       }
1621
1622       hDeDxExpected = hMomTranslation->ProjectionY("_py", binY + nMergeBins, binY + nMergeBins);
1623       Double_t meanDeDxExpectedTemp = hDeDxExpected->GetMean();
1624       delete hDeDxExpected;
1625
1626       if (meanDeDxExpectedTemp <= 0 || hRes3DprimeFit->GetYaxis()->FindBin(1. / meanDeDxExpectedTemp) != dEdxBin) {
1627         nMergeBins--;
1628         break;
1629       }
1630     }
1631
1632     Int_t binYhigh = binY + nMergeBins;
1633
1634     //if (binYhigh != binY) {
1635     //  printf("****Merged: %d for %f - %f\n", nMergeBins, hPreMap->GetYaxis()->GetBinLowEdge(binY),
1636     //         hPreMap->GetYaxis()->GetBinUpEdge(binYhigh));
1637     //}
1638
1639     // Average over some bins to get rid of threshold effects around the cut
1640     //Double_t referenceY = hPreMap->GetYaxis()->GetBinLowEdge(binYhigh);
1641     //Double_t binWidthY = hPreMap->GetYaxis()->GetBinWidth(binYhigh);
1642
1643     Int_t nMapRows = 1; // If no special handling: Only one row of the map filled
1644
1645     if (pThresholdTOFcut > 0) {
1646       // TOF cut causes step down in statistics behind the treshold towards larger p. Therefore: Merge the requested number of bins
1647       // above threshold - 1 (or only a part of them, if due to the same 1/dEdx-Bin some bins already got merged above)
1648       if (binYhigh >= binWithThreshold - 1 && binYhigh <= binWithThreshold - 1 + nMergeBinsAroundThreshold - 1) {
1649         // Check, if the next momentum bin after the merging, i.e. binYhigh + nMergeBinsAroundThreshold + 1, is apart in the row
1650         // of the map. If so, fill all rows before (usually at most 1 row) with the same value.
1651         // Take the bin AFTER the last merged bin because this bin should have better statistics
1652         hDeDxExpected = hMomTranslation->ProjectionY("_py", binWithThreshold + nMergeBinsAroundThreshold + 1,
1653                                                             binWithThreshold + nMergeBinsAroundThreshold + 1);
1654         Double_t meanDeDxExpectedTemp = hDeDxExpected->GetMean();
1655         delete hDeDxExpected;
1656
1657         if (meanDeDxExpectedTemp <= 0)  {
1658           printf("\nError in special handling: Failed to determine mean dEdxExpected of next row after \"the gap\"!\n");
1659           printedSomething = kTRUE;
1660         }
1661         else  {
1662           // nMapRows should be at least 1. If the row minus 1 of the bin after the map coincides with the "usual" row, then the second
1663           // argument in the following call should yield 1. For sanity, take Max(1,x) instead of x, although this should make no
1664           // difference, if everything is fine.
1665           nMapRows = TMath::Max(1, (hRes3DprimeFit->GetYaxis()->FindBin(1. / meanDeDxExpectedTemp) - 1) - dEdxBin + 1);
1666         }
1667
1668         Int_t nAdditionalMergeBins = binWithThreshold + nMergeBinsAroundThreshold - binYhigh;
1669         nMergeBins += nAdditionalMergeBins;
1670         binYhigh += nAdditionalMergeBins;
1671
1672         // Special handling is only really applied, if bins get merged due to it
1673         if (nAdditionalMergeBins > 0) {
1674           if (specialHandlingApplied) {
1675             printf("\nWarning in special handling: It is applied twice, which should not happen by design!\n");
1676             printedSomething = kTRUE;
1677           }
1678
1679           printf("\nAdditional merging due to special handling of threshold for TOF cut: %d p bins and %d 1/dEdx rows\n",
1680                  nAdditionalMergeBins, nMapRows);
1681           printedSomething = kTRUE;
1682
1683           specialHandlingApplied = kTRUE;
1684         }
1685       }
1686       /*OLD
1687       if (referenceY + nMergeBinsAroundThreshold / 4. * binWidthY >= pThresholdTOFcut && referenceY < pThresholdTOFcut) {
1688         // Check, if the next momentum bin after the merging, i.e. binYhigh + nMergeBinsAroundThreshold + 1, is appart in the row
1689         // of the map. If so, fill all rows before (usually at most 1 row) with the same value.
1690         // Take the bin AFTER the last merged bin because this bin should have better statistics
1691         hDeDxExpected = hMomTranslation->ProjectionY("_py", binYhigh + nMergeBinsAroundThreshold + 1,
1692                                                             binYhigh + nMergeBinsAroundThreshold + 1);
1693         Double_t meanDeDxExpectedTemp = hDeDxExpected->GetMean();
1694         delete hDeDxExpected;
1695
1696         if (meanDeDxExpectedTemp <= 0)  {
1697           printf("\nError in special handling: Failed to determine mean dEdxExpected of next row after \"the gap\"!\n");
1698           printedSomething = kTRUE;
1699         }
1700         else  {
1701           // nMapRows should be at least 1. If the row minus 1 of the bin after the map coincides with the "usual" row, then the second
1702           // argument in the following call should yield 1. For sanity, take Max(1,x) instead of x, although this should make no
1703           // difference, if everything is fine.
1704           nMapRows = TMath::Max(1, (hRes3DprimeFit->GetYaxis()->FindBin(1. / meanDeDxExpectedTemp) - 1) - dEdxBin + 1);
1705         }
1706
1707         nMergeBins += nMergeBinsAroundThreshold;
1708         binYhigh += nMergeBinsAroundThreshold;
1709         specialHandlingApplied = kTRUE;
1710       }*/
1711     }
1712
1713     for (Int_t binX = 1; binX <= hPreMap->GetXaxis()->GetNbins(); binX++) {
1714       TH1D* hTempProjectionZ = 0x0;
1715
1716       hTempProjectionZ = hPreMap->ProjectionZ("hTempProjectionZ", binX, binX, binY, binYhigh);
1717
1718       if (hTempProjectionZ->GetEntries() < 10) {
1719         printf("\nWarning: Too few entries for (%f, %f - %f (->%f)): skipped\n",
1720                hPreMap->GetXaxis()->GetBinCenter(binX), hPreMap->GetYaxis()->GetBinLowEdge(binY),
1721                hPreMap->GetYaxis()->GetBinUpEdge(binYhigh), hRes3DprimeFit->GetYaxis()->GetBinCenter(dEdxBin));
1722         printedSomething = kTRUE;
1723
1724         delete hTempProjectionZ;
1725         continue;
1726       }
1727
1728       Double_t meanWithoutFit = hTempProjectionZ->GetMean();
1729
1730       // Do not overwrite with data from lower statistics (larger momentum) -> Most of the cases should be caught by the merging!
1731       for (Int_t i = 0; i < nMapRows; i++)  {
1732         if (hRes3DprimeProfile->GetBinContent(binX, dEdxBin + i) <= 0) {
1733           hRes3DprimeProfile->SetBinContent(binX, dEdxBin + i, meanWithoutFit);
1734           hRes3DprimeProfile->SetBinError(binX, dEdxBin + i, (meanWithoutFit <= 0) ? 1000 : 1);
1735         }
1736         else  {
1737           printf("\nWarning: Data for mean (%f) already set (%f, %f (->%f) to %f - will NOT be overwritten\n",
1738                   meanWithoutFit, hPreMap->GetXaxis()->GetBinCenter(binX),
1739                   hPreMap->GetYaxis()->GetBinCenter(binY), hRes3DprimeFit->GetYaxis()->GetBinCenter(dEdxBin + i),
1740                   hRes3DprimeProfile->GetBinContent(binX, dEdxBin + i));
1741           printedSomething = kTRUE;
1742         }
1743       }
1744
1745
1746       TFitResult* fitRes = 0x0;
1747       if (useDoubleGaussFit) {
1748         fitRes = doubleGaussFit(hTempProjectionZ, 0.5 *(hPreMap->GetYaxis()->GetBinCenter(binY) + hPreMap->GetYaxis()->GetBinCenter(binYhigh)),
1749                                 splPr, splPi, "QNS") ;
1750       }
1751       else {
1752         Double_t lowThreshold = 0.6;
1753         Double_t highThreshold = 1.6;
1754         FindFitRange(hTempProjectionZ, lowThreshold, highThreshold);
1755         TFitResultPtr fitResPtr = hTempProjectionZ->Fit("gaus", "QNS", "", lowThreshold, highThreshold);
1756         fitRes = ((Int_t)fitResPtr == 0) ? (TFitResult*)fitResPtr.Get()->Clone() : 0x0;
1757       }
1758
1759       Double_t mean = 0;
1760       Double_t meanError = 0;
1761
1762       // If the fit failed, use the mean of the histogram as an approximation instead
1763       if (!fitRes) {
1764         mean = meanWithoutFit;
1765         meanError = 1000;
1766
1767         printf("\nWarning: Fit failed for (%f, %f (->%f)), entries: %.0f -> Using mean instead (%f +- %f)\n",
1768                hPreMap->GetXaxis()->GetBinCenter(binX), hPreMap->GetYaxis()->GetBinCenter(binY),
1769                hRes3DprimeFit->GetYaxis()->GetBinCenter(dEdxBin), hTempProjectionZ->GetEntries(), mean, meanError);
1770         printedSomething = kTRUE;
1771       }
1772       else {
1773         mean = fitRes->GetParams()[1];
1774         meanError = fitRes->GetErrors()[1];
1775       }
1776
1777       // Do not overwrite with data from lower statistics (larger momentum) -> Most of the cases should be caught by the merging!
1778       for (Int_t i = 0; i < nMapRows; i++)  {
1779         if (hRes3DprimeFit->GetBinContent(binX, dEdxBin + i) <= 0)  {
1780           hRes3DprimeFit->SetBinContent(binX, dEdxBin + i, mean);
1781           hRes3DprimeFit->SetBinError(binX, dEdxBin + i, meanError);
1782         }
1783         else  {
1784           printf("\nWarning: Data for fit (%f) already set (%f, %f (->%f) to %f - will NOT be overwritten\n",
1785                   mean, hPreMap->GetXaxis()->GetBinCenter(binX),
1786                   hPreMap->GetYaxis()->GetBinCenter(binY), hRes3DprimeFit->GetYaxis()->GetBinCenter(dEdxBin + i),
1787                   hRes3DprimeFit->GetBinContent(binX, dEdxBin + i));
1788           printedSomething = kTRUE;
1789         }
1790       }
1791
1792       delete fitRes;
1793       delete hTempProjectionZ;
1794
1795       //progressbar(100. * (((Double_t)(((binY - 1.) * hPreMap->GetXaxis()->GetNbins()) + binX)) / numTotalBins)); continue; //TODO NOW NOW
1796       // Extract the dependence on the number of clusters and fill the map with the resolution parameter
1797       TFitResult* fNclStats = extractClusterDependence(&hPreMapClusterResolved[0], numClusterBins, clusterLowBound, clusterUpBound,
1798                                                        binX, binY, binYhigh, c0, fSave, splPr, splPi);
1799
1800
1801       Double_t c1 = -1;
1802       Double_t c1Error = 1000;
1803
1804       if (fNclStats) {
1805         c1 = fNclStats->GetParams()[1];
1806         c1Error = fNclStats->GetErrors()[1];
1807       }
1808
1809       // Do not overwrite with data from lower statistics (larger momentum) -> Most of the cases should be caught by the merging!
1810       for (Int_t i = 0; i < nMapRows; i++)  {
1811         if (hSigmaPar1->GetBinContent(binX, dEdxBin + i) <= 0)  {
1812           hSigmaPar1->SetBinContent(binX, dEdxBin + i, c1);
1813           hSigmaPar1->SetBinError(binX, dEdxBin + i, c1Error);
1814         }
1815         else  {
1816           printf("\nWarning: Data (sigma par 1) for fit (%f) already set (%f, %f (->%f) to %f - will NOT be overwritten\n",
1817                  c1, hPreMap->GetXaxis()->GetBinCenter(binX), hPreMap->GetYaxis()->GetBinCenter(binY),
1818                  hSigmaPar1->GetYaxis()->GetBinCenter(dEdxBin + i),
1819                  hSigmaPar1->GetBinContent(binX, dEdxBin + i));
1820           printedSomething = kTRUE;
1821         }
1822       }
1823
1824       //delete fNclStats;
1825
1826
1827       progressbar(100. * (((Double_t)(((binY - 1.) * hPreMap->GetXaxis()->GetNbins()) + binX)) / numTotalBins));
1828     }
1829
1830     binY += nMergeBins;
1831   }
1832
1833   progressbar(100);
1834   printf("\n");
1835
1836   /*OLD
1837    // ---------------------------------------------------------------------------------
1838    // Use dEdx bins directly to create the map
1839
1840    //upperMapBound usually = 0.02
1841    tree->Project(Form("hRes3Dprime(%d,-1,1, %d,0.004,%f, 200,0.6,1.6)", binsX, binsY, upperMapBound),
1842    "dEdx/dEdxExpected:1./dEdx:tanTheta", "dEdx > 0 && dEdxExpected > 0 && pTPC < 3");
1843    //tree->Draw(Form("dEdx/dEdxExpected:1./dEdx:tanTheta>>hRes3Dprime(%d,-1,1, %d,0.004,%f, 200,0.6,1.6)", binsX, binsY, upperMapBound),
1844    //           "dEdx > 0 && dEdxExpected > 0 && pTPC < 3");
1845    //tree->Draw("dEdx/dEdxExpected:1./dEdx:tanTheta>>hRes3Dprime(40,-1,1, 40,0.002,0.0225, 200,0.6,1.6)");// "acceptedByPhiPrimeCut == 0");
1846    //tree->Draw("dEdx/dEdxExpected:dEdx:tanTheta>>hRes3Dprime(30,-1,1, 100,45.,600., 200,0.6,1.6)");
1847    TH3F* hRes3Dprime = (TH3F*)gDirectory->Get("hRes3Dprime");
1848    TProfile2D* hTemp = hRes3Dprime->Project3DProfile("yx");
1849    hTemp->SetName("hTemp");
1850    TH2D* hRes3DprimeProfile = hTemp->ProjectionXY();
1851    hRes3DprimeProfile->SetName("hRes3DprimeProfile");
1852    normaliseHisto(hRes3DprimeProfile, -0.1, 0.1, kTypeDeltaPrime);
1853    hRes3DprimeProfile->GetZaxis()->SetRangeUser(0.95, 1.15);
1854    SetHistAxis(hRes3DprimeProfile, kTypeDeltaPrime);
1855    hRes3DprimeProfile->DrawCopy("surf1");
1856    canv3Dprime->cd(2);
1857    //OLD
1858    //hRes3Dprime->FitSlicesZ();
1859    //TH2D *hRes3DprimeFit = (TH2D*)gDirectory->Get("hRes3Dprime_1");
1860
1861
1862    TH2D *hRes3DprimeFit = (TH2D*)hRes3DprimeProfile->Clone("hRes3DprimeFit");
1863    hRes3DprimeFit->Reset();
1864
1865    const Int_t numBinsSpecialHandling = 2;
1866    for (Int_t binY = 1; binY <= hRes3Dprime->GetYaxis()->GetNbins(); binY++) {
1867      Bool_t specialHandling = kFALSE;
1868      Double_t referenceY = hRes3Dprime->GetYaxis()->GetBinLowEdge(binY);
1869      Double_t binWidthY = hRes3Dprime->GetYaxis()->GetBinWidth(binY);
1870
1871      // Average over some bins to get rid of threshold effects around the cut
1872      if (TOFcutDeDx > 0) {
1873        if (referenceY + numBinsSpecialHandling / 2. * binWidthY >= 1. / TOFcutDeDx && referenceY < 1. / TOFcutDeDx)
1874          specialHandling = kTRUE;
1875      }
1876
1877      for (Int_t binX = 1; binX <= hRes3Dprime->GetXaxis()->GetNbins(); binX++) {
1878        TH1D* hTempProjectionZ = 0x0;
1879        if (specialHandling)
1880          hTempProjectionZ = hRes3Dprime->ProjectionZ("hTempProjectionZ", binX, binX, binY, binY + numBinsSpecialHandling, "e");
1881        else
1882          hTempProjectionZ = hRes3Dprime->ProjectionZ("hTempProjectionZ", binX, binX, binY, binY, "e");
1883
1884        Int_t binZ = 0;
1885        Int_t binZmax = 0;
1886
1887        // Start from the very right side and go to the left -> Remove some of the first bins with non-zero content from the fitting.
1888        // This is important, since these bin might be affected from binning effects and spoil the gaussian shape.
1889        for (binZ = hTempProjectionZ->GetXaxis()->GetNbins(); binZ >= 1; binZ--) {
1890          if (hTempProjectionZ->GetBinContent(binZ) > 0)  {
1891            binZmax = binZ;
1892
1893            break;
1894          }
1895        }
1896
1897        binZmax = TMath::Max(1, binZmax - 2);
1898
1899        TFitResultPtr fitRes = hTempProjectionZ->Fit("gaus", "QS", "", hTempProjectionZ->GetXaxis()->GetBinLowEdge(1),
1900        hTempProjectionZ->GetXaxis()->GetBinUpEdge(binZmax));
1901        Double_t mean = 0;
1902
1903        // If the fit failed, use the mean of the histogram as an approximation instead
1904        if (((Int_t)fitRes) != 0) {
1905          printf("Warning: Fit failed for (%f, %f), entries: %.0f -> Using mean instead (%f)\n",
1906          hRes3DprimeProfile->GetXaxis()->GetBinCenter(binX), hRes3DprimeProfile->GetYaxis()->GetBinCenter(binY),
1907          fitRes->GetParams()[1], hTempProjectionZ->GetEntries(), hTempProjectionZ->GetMean());
1908          mean = hTempProjectionZ->GetMean();
1909        }
1910        else
1911          mean = fitRes->GetParams()[1];
1912        hRes3DprimeFit->SetBinContent(binX, binY, mean);
1913
1914        if (specialHandling)  {
1915          for (Int_t i = 1; i <= numBinsSpecialHandling && i + binY <= hRes3DprimeFit->GetYaxis()->GetNbins(); i++)
1916            hRes3DprimeFit->SetBinContent(binX, binY + i, mean);
1917        }
1918
1919
1920        delete hTempProjectionZ;
1921    }
1922
1923    if (specialHandling)  {
1924      binY += numBinsSpecialHandling;
1925    }
1926    }
1927
1928    normaliseHisto(hRes3DprimeFit, -0.1, 0.1, kTypeDeltaPrime);
1929    hRes3DprimeFit->GetZaxis()->SetRangeUser(0.95, 1.15);
1930    SetHistAxis(hRes3DprimeFit, kTypeDeltaPrime);
1931    hRes3DprimeFit->DrawCopy("surf1");
1932
1933    delete hTemp;
1934    delete hRes3Dprime;
1935    */
1936
1937
1938
1939
1940   TH2D* hPureResults = (TH2D*)hRes3DprimeFit->Clone("hPureResults");
1941   TH2D* hPureResultsSigmaPar1 = (TH2D*)hSigmaPar1->Clone("hPureResultsSigmaPar1");
1942
1943   printf("\n\nEliminating outliers....\n\n");
1944   eliminateOutliers(hRes3DprimeFit, outlierThreshold);
1945
1946   printf("\n\nSame for the sigma map....\n\n");
1947   eliminateOutliers(hSigmaPar1, outlierThreshold);
1948
1949   TH2D* hRes3DprimeDiffSmooth = 0x0;
1950
1951   if (doSmoothing) {
1952     // --------------------------------------------------
1953     // Determine difference between smoothed and pure map (gauss fit)
1954     hRes3DprimeDiffSmooth = (TH2D*)hRes3DprimeFit->Clone("hRes3DprimeDiffSmooth");
1955     hRes3DprimeDiffSmooth->SetTitle("tan(#theta) and 1/(dE/dx) dependence of #Delta': Smooth(Gauss) - Gauss");
1956
1957     // Smooth the histos
1958
1959     // If the smoothing is to be applied, the bins with fit errors (content <= 0) must be eliminated in order not to bias the smoothing.
1960     // In case of no smoothing, the refinement will take care of such bins....
1961     eliminateNonPositivePointsInHisto(hRes3DprimeProfile);
1962     eliminateNonPositivePointsInHisto(hRes3DprimeFit);
1963     eliminateNonPositivePointsInHisto(hSigmaPar1);
1964
1965
1966     hRes3DprimeProfile->Smooth();
1967     hRes3DprimeFit->Smooth();
1968
1969     hSigmaPar1->Smooth();
1970
1971     hRes3DprimeDiffSmooth->Add(hRes3DprimeFit, -1);
1972     //hRes3DprimeDiffSmooth->GetZaxis()->SetRangeUser(-0.05, 0.05);
1973     SetHistAxis(hRes3DprimeDiffSmooth, kTypeDeltaPrime);
1974   }
1975
1976   /* Normalisation is now done after the extrapolation
1977    *   // Use kNoNormalisation for final QA
1978    *   if (mapType == kSmallEtaNormalisation) {
1979    *     normaliseHisto(hRes3DprimeProfile, 0., 0.11, kTypeDeltaPrime);
1980    *     normaliseHisto(hRes3DprimeFit, 0., 0.11, kTypeDeltaPrime);
1981    }
1982    else if (mapType == kLargeEtaNormalisation) {
1983     normaliseHisto(hRes3DprimeProfile, 0.81, 0.99, kTypeDeltaPrime);
1984   normaliseHisto(hRes3DprimeFit, 0.81, 0.99, kTypeDeltaPrime);
1985    }
1986    */
1987
1988
1989   //hRes3DprimeProfile->GetZaxis()->SetRangeUser(0.95, 1.15);
1990   SetHistAxis(hRes3DprimeProfile, kTypeDeltaPrime);
1991
1992   //hRes3DprimeFit->GetZaxis()->SetRangeUser(0.95, 1.15);
1993   SetHistAxis(hRes3DprimeFit, kTypeDeltaPrime);
1994
1995   SetHistAxis(hSigmaPar1, kTypeSigmaDeltaPrime);
1996
1997
1998   if (saveResults)  {
1999     fSave->cd();
2000     hPreMap->Write();
2001     hMomTranslation->Write();
2002   }
2003
2004   delete hPreMap;
2005   delete hMomTranslation;
2006
2007
2008   // --------------------------------------------------
2009   // Determine difference between gauss fit and pure mean
2010   TH2D* hRes3DprimeDiff = (TH2D*)hRes3DprimeFit->Clone("hRes3DprimeDiff");
2011   hRes3DprimeDiff->SetTitle("tan(#theta) and 1/(dE/dx) dependence of #Delta': Gauss - Profile");
2012   hRes3DprimeDiff->Add(hRes3DprimeProfile, -1);
2013   hRes3DprimeDiff->GetZaxis()->SetRangeUser(-0.05, 0.05);
2014   SetHistAxis(hRes3DprimeDiff, kTypeDeltaPrime);
2015
2016   // --------------------------------------------------
2017   // Refine the map
2018   printf("\n\nMap created. Started to refine....\n\n");
2019
2020   Double_t factorX = 6;
2021   Double_t factorY = 6;
2022
2023   if (resolutionFactorX != 1)
2024     factorX *= resolutionFactorX;
2025   if (resolutionFactorY != 1)
2026     factorY *= resolutionFactorY;
2027
2028
2029
2030   Bool_t skipBinsAtBoundaries = kFALSE; // Diese Bins vielleicht doch nicht rausnehmen. Der Anstieg in Sigma bei extrem kleinen Impulsen ist insofern richtig, da hier die Splines abweichen und somit eine große Streuung auftritt. Und bei großen Impulsen sollten im Moment die V0s Abhilfe von fehlerhaften Abweichungen schaffen. Beim Mean gilt dasselbe.
2031   TH2D* hRefinedNoNorm = refineHistoViaLinearExtrapolation(hRes3DprimeFit, factorX, factorY, kNoNormalisation, fSave, kFALSE,
2032                                                            skipBinsAtBoundaries);
2033   if (hRefinedNoNorm) {
2034     hRefinedNoNorm->SetName(Form("hRefined%s", suffixMapType[kNoNormalisation].Data()));
2035     SetHistAxis(hRefinedNoNorm, kTypeDeltaPrime);
2036   }
2037
2038   TH2D* hRefinedSmallEtaNorm = refineHistoViaLinearExtrapolation(hRes3DprimeFit, factorX, factorY, kSmallEtaNormalisation, fSave, kFALSE,
2039                                                                  skipBinsAtBoundaries);
2040   if (hRefinedSmallEtaNorm)  {
2041     hRefinedSmallEtaNorm->SetName(Form("hRefined%s", suffixMapType[kSmallEtaNormalisation].Data()));
2042     SetHistAxis(hRefinedSmallEtaNorm, kTypeDeltaPrime);
2043   }
2044
2045   TH2D* hRefinedLargeEtaNorm = refineHistoViaLinearExtrapolation(hRes3DprimeFit, factorX, factorY, kLargeEtaNormalisation, fSave, kFALSE,
2046                                                                  skipBinsAtBoundaries);
2047   if (hRefinedLargeEtaNorm) {
2048     hRefinedLargeEtaNorm->SetName(Form("hRefined%s", suffixMapType[kLargeEtaNormalisation].Data()));
2049     SetHistAxis(hRefinedLargeEtaNorm, kTypeDeltaPrime);
2050   }
2051
2052
2053
2054   printf("\nSame for the sigma map.....\n\n");
2055
2056   // Normalisation makes no sense for sigma map (will anyway be created for corrected data only)
2057   TH2D* hRefinedSigmaPar1 = refineHistoViaLinearExtrapolation(hSigmaPar1, factorX, factorY, kNoNormalisation, fSave, kTRUE,
2058                                                               skipBinsAtBoundaries);
2059   if (hRefinedSigmaPar1) {
2060     hRefinedSigmaPar1->SetName("hThetaMapSigmaPar1");
2061     SetHistAxis(hRefinedSigmaPar1, kTypeSigmaDeltaPrime);
2062   }
2063
2064
2065   printf("\nDone!\n\n");
2066
2067   if (resolutionFactorX != 1 || resolutionFactorY != 1) {
2068     printf("\n***WARNING***: Resolution factor != 1 used for map creation! Resolution scaled by factor (%f, %f)\n\n",
2069            resolutionFactorX, resolutionFactorY);
2070   }
2071
2072   if (pThresholdTOFcut > 0) {
2073     printf("\n***INFO***: Requested special handling for map around p = %f GeV/c to take care of potential TOF cut around this value!",
2074            pThresholdTOFcut);
2075     if (specialHandlingApplied)
2076       printf(" -> Special handling was used!\n\n");
2077     else
2078       printf(" -> Special handling was not used (seemed to be not necessary)!\n\n");
2079   }
2080   else  {
2081     printf("\n***WARNING***: No special care taken of the TOF cut. Fine for MC, but not for data. Please check, if it is correct!\n\n");
2082   }
2083
2084   printf("WARNING: Remember to switch on/off ncl cut w.r.t. geo cut!\n");
2085
2086   if (saveResults)  {
2087     TNamed* c0Info = new TNamed("c0", Form("%f", c0));
2088
2089     fSave->cd();
2090     hPureResults->Write();
2091     hRes3DprimeProfile->Write();
2092     hRes3DprimeFit->Write();
2093     hRes3DprimeDiff->Write();
2094     if (doSmoothing)
2095       hRes3DprimeDiffSmooth->Write();
2096
2097     hPureResultsSigmaPar1->Write();
2098     hSigmaPar1->Write();
2099
2100     if (hRefinedNoNorm)
2101       hRefinedNoNorm->Write();
2102     if (hRefinedSmallEtaNorm)
2103       hRefinedSmallEtaNorm->Write();
2104     if (hRefinedLargeEtaNorm)
2105       hRefinedLargeEtaNorm->Write();
2106
2107     if (hRefinedSigmaPar1)
2108       hRefinedSigmaPar1->Write();
2109
2110     c0Info->Write();
2111
2112     TNamed* settings = new TNamed(
2113       Form("Settings for map: resolutionFactor (%f, %f),  pThresholdTOFcut %f, nMergeBinsAroundThreshold %d, pThresholdV0cut %f, pThresholdV0plusTOFcut %f, outlierThreshold %f, doSmoothing %d, c0 %f, nclCut %d, useDoubleGaussFit %d\n",
2114            resolutionFactorX, resolutionFactorY, pThresholdTOFcut, nMergeBinsAroundThreshold, pThresholdV0cut, pThresholdV0plusTOFcut,
2115            outlierThreshold, doSmoothing, c0, nclCut, useDoubleGaussFit), "");
2116     settings->Write();
2117
2118     f->Close();
2119     fSave->Close();
2120
2121     delete c0Info;
2122
2123     gSystem->Exec(Form("echo \"%s\" >> %s/settingsForMapCreation.txt",
2124                        TString(settings->GetName()).ReplaceAll("Settings for map",
2125                                                                Form("Settings for map %s", saveFileName.Data())).Data(), path.Data()));
2126     return 0x0;
2127   }
2128
2129   if (returnMapType == kSmallEtaNormalisation)
2130     return hRefinedSmallEtaNorm;
2131   else if (returnMapType == kLargeEtaNormalisation)
2132     return hRefinedLargeEtaNorm;
2133   else if (returnMapType == kNoNormalisation)
2134     return hRefinedNoNorm;
2135
2136   printf("Map type %d is not supported. Returning 0x0...\n", returnMapType);
2137
2138   return 0x0;
2139 }