PWGPP/TPC/macros/PIDCalib/extractMultiplicityDependence.C

   1 #include "TTree.h"
   2 #include "TCanvas.h"
   3 #include "TF1.h"
   4 #include "TFitResult.h"
   5 #include "TH1D.h"
   6 #include "TH2D.h"
   7 #include "TGraphErrors.h"
   8 #include "TList.h"
   9 #include "TString.h"
  10 #include "TFile.h"
  11 #include "TMath.h"
  12 #include "TDatime.h"
  13
  14 #include "TSpline.h"
  15 #include "AliPID.h"
  16
  17 #include <iostream>
  18 #include "ProgressBar.h"
  19 #include "THnSparseDefinitions.h"
  20
  21 //________________________________________________________________________
  22 void FitSlicesY(TH2 *hist, Double_t heightFractionForRange, Int_t cutThreshold, TString fitOption, TObjArray *arr)
  23 {
  24   //heightPercentageForRange
  25   // custom slices fit
  26   //
  27
  28   if (!hist) return;
  29   if (!arr) return;
  30
  31   // If old style is to be used
  32   /*
  33   hist->FitSlicesY(0, 0, -1, cutThreshold, fitOption.Data(), arr);
  34   return;
  35   */
  36
  37
  38   arr->Clear();
  39   arr->SetOwner();
  40   arr->Expand(4);
  41
  42   TAxis *axis=hist->GetXaxis();
  43   const TArrayD *bins = axis->GetXbins();
  44   TH1D** hList = new TH1D*[4];
  45
  46   for (Int_t i = 0; i < 4; i++) {
  47     delete gDirectory->FindObject(Form("%s_%d", hist->GetName(), i));
  48
  49     if (bins->fN == 0) {
  50       hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF",
  51                           hist->GetNbinsX(), axis->GetXmin(), axis->GetXmax());
  52     } else {
  53       hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF", hist->GetNbinsX(), bins->fArray);
  54     }
  55
  56     (*arr)[i] = hList[i];
  57   }
  58
  59   for (Int_t ibin=axis->GetFirst(); ibin<=axis->GetLast(); ++ibin){
  60     TH1 *h=hist->ProjectionY("_temp",ibin,ibin);
  61     if (!h)
  62       continue;
  63
  64     if (h->GetEntries() < cutThreshold) {
  65       delete h;
  66       continue;
  67     }
  68
  69     // Average around maximum bin -> Might catch some outliers
  70     Int_t maxBin = h->GetMaximumBin();
  71     Double_t maxVal = h->GetBinContent(maxBin);
  72
  73     if (maxVal < 2) { // It could happen that all entries are in overflow/underflow; don't fit in this case
  74       delete h;
  75       continue;
  76     }
  77
  78     UChar_t usedBins = 1;
  79     if (maxBin > 1) {
  80       maxVal += h->GetBinContent(maxBin - 1);
  81       usedBins++;
  82     }
  83     if (maxBin < h->GetNbinsX()) {
  84       maxVal += h->GetBinContent(maxBin + 1);
  85       usedBins++;
  86     }
  87     maxVal /= usedBins;
  88
  89     Double_t thresholdFraction = heightFractionForRange * maxVal;
  90     Int_t lowThrBin = TMath::Max(1, h->FindFirstBinAbove(thresholdFraction));
  91     Int_t highThrBin = TMath::Min(h->GetNbinsX(), h->FindLastBinAbove(thresholdFraction));
  92
  93     Double_t lowThreshold = h->GetBinCenter(lowThrBin);
  94     Double_t highThreshold = h->GetBinCenter(highThrBin);
  95
  96     TFitResultPtr res = h->Fit("gaus", Form("%sS", fitOption.Data()), "", lowThreshold, highThreshold);
  97
  98     if ((Int_t)res == 0) {
  99       Int_t resBin = ibin;
 100       hList[0]->SetBinContent(resBin,res->GetParams()[0]);
 101       hList[0]->SetBinError(resBin,res->GetErrors()[0]);
 102       hList[1]->SetBinContent(resBin,res->GetParams()[1]);
 103       hList[1]->SetBinError(resBin,res->GetErrors()[1]);
 104       hList[2]->SetBinContent(resBin,res->GetParams()[2]);
 105       hList[2]->SetBinError(resBin,res->GetErrors()[2]);
 106       hList[3]->SetBinContent(resBin,res->Ndf()>0?res->Chi2()/res->Ndf():0);
 107     }
 108
 109     delete h;
 110   }
 111
 112   delete [] hList;
 113 }
 114
 115 //__________________________________________________________________________________________
 116 Int_t extractMultiplicityDependence(TString pathTree, Double_t widthFactor /*=1*/, Int_t multStepSize /*= 400*/, Bool_t isMC,
 117                                     Bool_t correct = kFALSE, TString fileNameTree = "bhess_PIDetaTree.root", TString treeName = "fTree",
 118                                     TString pathNameSplines = "splines_10h.pass2.root", TString splineName = "TSPLINE3_DATA_PROTON_LHC10H_PASS2_PBPB_MEAN")
 119 {
 120   const Double_t massProton = AliPID::ParticleMass(AliPID::kProton);
 121   const Double_t massPion = AliPID::ParticleMass(AliPID::kPion);
 122
 123   Double_t mass = massProton;
 124
 125
 126   // Extract the data tree
 127   TFile* fTree = TFile::Open(Form("%s/%s", pathTree.Data(), fileNameTree.Data()));
 128   if (!fTree)  {
 129     std::cout << "Failed to open file \"" << Form("%s/%s", pathTree.Data(), fileNameTree.Data()) << "\"!" << std::endl;
 130     return -1;
 131   }
 132
 133   TTree* tree = dynamic_cast<TTree*>(fTree->Get(treeName.Data()));
 134   if (!tree) {
 135     std::cout << "Failed to load data tree!" << std::endl;
 136     return -1;
 137   }
 138
 139   // Output file
 140   TDatime daTime;
 141   TString savefileName = Form("%s_multiplicityDependence_%04d_%02d_%02d__%02d_%02d.root", fileNameTree.ReplaceAll(".root", "").Data(),
 142                               daTime.GetYear(), daTime.GetMonth(), daTime.GetDay(), daTime.GetHour(), daTime.GetMinute());
 143
 144   TFile* fSave = TFile::Open(Form("%s/%s", pathTree.Data(), savefileName.Data()), "recreate");
 145   if (!fSave) {
 146     std::cout << "Failed to open save file \"" << Form("%s/%s", pathTree.Data(), savefileName.Data()) << "\"!" << std::endl;
 147     return -1;
 148   }
 149
 150   Long64_t nTreeEntries = tree->GetEntriesFast();
 151
 152   Double_t pTPC = 0.;
 153   Double_t dEdx = 0.;
 154   //Double_t dEdxExpected = 0.;
 155   Double_t tanTheta = 0.;
 156   UShort_t tpcSignalN = 0.;
 157   Int_t fMultiplicity = 0.;
 158   UChar_t pidType = 0;
 159
 160   tree->SetBranchStatus("*", 0); // Disable all branches
 161   tree->SetBranchStatus("pTPC", 1);
 162   tree->SetBranchStatus("dEdx", 1);
 163   //tree->SetBranchStatus("dEdxExpected", 1);
 164   tree->SetBranchStatus("tanTheta", 1);
 165   tree->SetBranchStatus("fMultiplicity", 1);
 166   tree->SetBranchStatus("tpcSignalN", 1);
 167   tree->SetBranchStatus("pidType", 1);
 168
 169
 170   tree->SetBranchAddress("pTPC", &pTPC);
 171   tree->SetBranchAddress("dEdx", &dEdx);
 172   //tree->SetBranchAddress("dEdxExpecttanThetaCentered", &dEdxExpected);
 173   tree->SetBranchAddress("tanTheta", &tanTheta);
 174   tree->SetBranchAddress("fMultiplicity", &fMultiplicity);
 175   tree->SetBranchAddress("tpcSignalN", &tpcSignalN);
 176   tree->SetBranchAddress("pidType", &pidType);
 177
 178   const Int_t numMultBins = TMath::Ceil((20000. - 0.) / multStepSize);
 179
 180   const Int_t numPbins = 24;
 181   Double_t pBinEdges[numPbins * 2] = {
 182     0.3, 0.3 + 0.005*widthFactor,
 183     0.35, 0.35 + 0.005*widthFactor,
 184     0.4, 0.4 + 0.005*widthFactor,
 185     0.45, 0.45 + 0.005*widthFactor,
 186     0.5, 0.5 + 0.005*widthFactor,
 187     0.55, 0.55 + 0.005*widthFactor,
 188     0.6, 0.6 + 0.005*widthFactor,
 189     0.65, 0.65 + 0.005*widthFactor,
 190     0.7, 0.7 + 0.005*widthFactor,
 191     0.8, 0.8 + 0.01*widthFactor,
 192     0.9, 0.9 + 0.01*widthFactor,
 193     1.0, 1.0 + 0.01*widthFactor,
 194     1.1, 1.1 + 0.01*widthFactor,
 195     1.2, 1.2 + 0.01*widthFactor,
 196     1.3, 1.3 + 0.01*widthFactor,
 197     1.4, 1.4 + 0.01*widthFactor,
 198     1.5, 1.5 + 0.01*widthFactor,
 199     1.6, 1.6 + 0.02*widthFactor,
 200     1.7, 1.7 + 0.02*widthFactor,
 201     1.8, 1.8 + 0.02*widthFactor,
 202     1.9, 1.9 + 0.02*widthFactor,
 203     2.0, 2.0 + 0.1*widthFactor,
 204     2.5, 2.5 + 0.1*widthFactor,
 205     3.0, 3.0 + 0.15*widthFactor };
 206
 207   const Int_t numTanThetaBins = 10;
 208
 209   Double_t tanThetaBinEdges[numTanThetaBins * 2] = {
 210     -1.0, -0.8,
 211     -0.8, -0.6,
 212     -0.6, -0.4,
 213     -0.4, -0.2,
 214     -0.2, 0.0,
 215     0.0, 0.2,
 216     0.2, 0.4,
 217     0.4, 0.6,
 218     0.6, 0.8,
 219     0.8, 1.0 };
 220     /*TODO tanThetaAbs
 221   Double_t tanThetaBinEdges[numTanThetaBins * 2] = {
 222     0.0, 0.1,
 223     0.1, 0.2,
 224     0.2, 0.3,
 225     0.3, 0.4,
 226     0.4, 0.5,
 227     0.5, 0.6,
 228     0.6, 0.7,
 229     0.7, 0.8,
 230     0.8, 0.9,
 231     0.9, 1.0 };*/
 232
 233   TH2D* h2[numPbins][numTanThetaBins];
 234   TH2D* h2AllEta[numPbins];
 235
 236   for (Int_t i = 0; i < numPbins; i++) {
 237     h2AllEta[i] = new TH2D(Form("h2AllEta_%d", i), Form("p %.2f - %.2f GeV/c, all #eta; Multiplicity; dEdx", pBinEdges[2*i], pBinEdges[2*i+1]),
 238                            numMultBins, 0, 20000, 590, 10.0, 600.0);
 239
 240     for (Int_t j = 0; j < numTanThetaBins; j++) {
 241       h2[i][j] = new TH2D(Form("h2_%d_%d", i, j), Form("p %.2f - %.2f GeV/c & %.2f < tanTheta #leq %.2f; Multiplicity; dEdx",
 242       //TODO tanThetaAbs h2[i][j] = new TH2D(Form("h2_%d_%d", i, j), Form("p %.2f - %.2f GeV/c & %.2f < |tanTheta| #leq %.2f; Multiplicity; dEdx",
 243                                                        pBinEdges[2*i], pBinEdges[2*i+1], tanThetaBinEdges[2*j], tanThetaBinEdges[2*j+1]),
 244                           numMultBins, 0, 20000, 590, 10.0, 600.0);
 245     }
 246   }
 247
 248   TSpline3* splProton = 0x0;
 249
 250   if (correct) {
 251     printf("CORRECTING data for multiplicity dependence...TODO OBSOLETE - DO NOT USE!!!!\n");
 252
 253     TFile* f = 0x0;
 254     f = TFile::Open(pathNameSplines.Data());
 255     if (!f)  {
 256       std::cout << "Failed to open spline file \"" << pathNameSplines.Data() << "\"!" << std::endl;
 257       return -1;
 258     }
 259
 260     splProton = (TSpline3*)f->Get(splineName.Data());
 261     if (!splProton) {
 262       std::cout << "Failed to load splines: " << splineName.Data() << "!" << std::endl;
 263       return -1;
 264     }
 265
 266     if (splineName.Contains("PION")){
 267       printf("Pion mass assumed: %f...\n", massPion);
 268       mass = massPion;
 269     }
 270     else {
 271       printf("Proton mass assumed: %f...\n", massProton);
 272       mass = massProton;
 273     }
 274   }
 275   else
 276     printf("NOT CORRECTING data for multiplicity dependence...\n");
 277
 278   progressbar(0.);
 279
 280   TF1 cutFunc("cutFunc","50./TMath::Power(x,1.3)", 0.05, 6);
 281
 282   TF1 corrFunc("corrFunc", "pol2", 0, 0.01);
 283   corrFunc.SetParameter(0, 5.5e-6);
 284   corrFunc.SetParameter(1, -0.00436);
 285   corrFunc.SetParameter(2, 0.103);
 286
 287   for (Long64_t i = 0; i < nTreeEntries; i++) {
 288     tree->GetEntry(i);
 289
 290     if (i % 100000 == 0)
 291       progressbar(100. * (((Double_t)i) / nTreeEntries));
 292
 293     // Filter tracks according to shape recognition
 294     //if (tpcSignalN < 60 || (!isMC && dEdx < cutFunc.Eval(pTPC))) continue;
 295     if (dEdx <= 0 || tpcSignalN < 60) {
 296       continue;
 297     }
 298
 299     if (!isMC) {
 300       if ((pTPC < 0.6 && pidType != kTPCid) ||               // No V0s/TOF below 0.6 due to bad statistics
 301       //if ((pTPC < 0.6 && pidType == kTPCandTOFid) ||         // No TOF below 0.6 GeV/c due to bad statistics
 302           //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 < pThreholdV0cut && pidType != kTPCandTOFid) || // Do not mix V0's and TOF in this range
 303           (pTPC >= 0.6 && pTPC < 0.9 && pidType != kTPCandTOFid) ||
 304           (pTPC >= 0.9 && pidType != kV0idPlusTOFaccepted && pidType != kV0idNoTOF))      // Only V0's WITH TOF above pThreholdV0cut due to contamination
 305         continue;
 306     }
 307
 308     // CORRECT multiplicity dependence -> WARNING: Better take into account the eta dependence for dEdxSplines for the determination of corrFactor!!!
 309     if (correct) {
 310       Double_t dEdxSplines = 50. * splProton->Eval(pTPC / mass);
 311       Double_t dEdxSplinesInv = 1. / dEdxSplines;
 312       Double_t relSlope = corrFunc.Eval(dEdxSplinesInv);
 313       Double_t corrFactor = relSlope * fMultiplicity;
 314       dEdx /= 1. + corrFactor;
 315     }
 316
 317
 318
 319     //TODO tanThetaAbs Double_t absTanTheta = TMath::Abs(tanTheta);
 320
 321     for (Int_t j = 0; j < numPbins; j++) {
 322       if (pTPC >= pBinEdges[2*j] && pTPC < pBinEdges[2*j+1]) {
 323         h2AllEta[j]->Fill(fMultiplicity, dEdx);
 324
 325         for (Int_t k = 0; k < numTanThetaBins; k++) {
 326           if (tanTheta >= tanThetaBinEdges[2*k] && tanTheta < tanThetaBinEdges[2*k+1]) {
 327           //TODO tanThetaAbs if (absTanTheta >= tanThetaBinEdges[2*k] && absTanTheta < tanThetaBinEdges[2*k+1]) {
 328             h2[j][k]->Fill(fMultiplicity, dEdx);
 329             break;
 330           }
 331         }
 332         break;
 333       }
 334     }
 335   }
 336
 337   progressbar(100.);
 338
 339   printf("\n\n");
 340
 341
 342   const Double_t MCfitThreshold = 0.006;
 343   const Double_t heightFractionForRange = 0.1;
 344   TObjArray arr;
 345
 346   { // Fitting for all eta
 347     TGraphErrors* gSlvsOffset = new TGraphErrors(numPbins);
 348     TGraphErrors* gSlvsInvOffset = new TGraphErrors(numPbins);
 349
 350     TGraphErrors* gSigmaSlopevsdEdx = new TGraphErrors(numPbins);
 351     TGraphErrors* gSigmaSlopevsInvdEdx = new TGraphErrors(numPbins);
 352
 353     fSave->mkdir("allTanTheta");
 354
 355     for (Int_t i = 0; i < numPbins; i++) {
 356       gSlvsOffset->SetPoint(i, -10, 0);
 357       gSlvsOffset->SetPointError(i, 0, 0);
 358
 359       gSlvsInvOffset->SetPoint(i, -10, 0);
 360       gSlvsInvOffset->SetPointError(i, 0, 0);
 361
 362       gSigmaSlopevsdEdx->SetPoint(i, -10, 0);
 363       gSigmaSlopevsdEdx->SetPointError(i, 0, 0);
 364
 365       gSigmaSlopevsInvdEdx->SetPoint(i, -10, 0);
 366       gSigmaSlopevsInvdEdx->SetPointError(i, 0, 0);
 367
 368       if (h2AllEta[i]->GetEntries() < 10)
 369         continue;
 370
 371       printf("Momentum %.2f - %.2f GeV/c, all eta:\n", pBinEdges[2*i], pBinEdges[2*i+1]);
 372       TCanvas* c = new TCanvas(Form("canv_pBin%d_allEta", i), Form("canv_pBin%d_allEta", i), 100,10,1380,800);
 373       FitSlicesY(h2AllEta[i], heightFractionForRange, 10, "QN", &arr);
 374       h2AllEta[i]->GetYaxis()->SetRange(h2AllEta[i]->FindFirstBinAbove(1, 2), h2AllEta[i]->FindLastBinAbove(1,2));
 375       h2AllEta[i]->Draw("colz");
 376
 377       TH1D* hMean = (TH1D*)arr.At(1);//new TH1D((TH1D&)*arr.At(1));
 378       TH1D* hSigma = (TH1D*)arr.At(2);//new TH1D((TH1D&)*arr.At(2));
 379       hMean->SetName(Form("hMean_allEta_pBin%d", i));
 380       hSigma->SetName(Form("hSigma_allEta_pBin%d", i));
 381       hSigma->SetLineColor(kMagenta);
 382
 383       TH1D* hSigmaRel = new TH1D(*hSigma);
 384       hSigmaRel->SetName(Form("hSigmaRel_allEta_pBin%d", i));
 385       hSigmaRel->Divide(hMean);
 386       hSigmaRel->SetLineColor(kMagenta + 2);
 387
 388       TFitResultPtr fitRes = hMean->Fit("pol1", "s", "same", 0, 12000);
 389       TFitResultPtr fitResSigma = hSigmaRel->Fit("pol1", "s", "same", 0, 12000);
 390
 391       hMean->DrawClone("same");
 392
 393       if (((Int_t)fitRes) != 0) {
 394         printf("Fit failed!\n");
 395       }
 396       else {
 397         Double_t p0 = fitRes->GetParams()[0];
 398         Double_t p1 = fitRes->GetParams()[1];
 399         Double_t err0 = fitRes->GetErrors()[0];
 400         Double_t err1 = fitRes->GetErrors()[1];
 401
 402         Double_t totError = 99999;
 403         if (p0 != 0 && p1 != 0) {
 404           totError = TMath::Abs(p1 / p0 * TMath::Sqrt(TMath::Power(err0 / p0, 2) + TMath::Power(err1 / p1, 2)));
 405         }
 406         printf("Result: p1 / p0 = %f / %f = %e\n\n", p1, p0, (p0 != 0 ? p1 / p0 : -99999));
 407         if (p0 > 0) {
 408           gSlvsOffset->SetPoint(i, p0, p1/p0);
 409           gSlvsOffset->SetPointError(i, err0, totError);
 410
 411           gSlvsInvOffset->SetPoint(i, 1./p0, p1/p0);
 412           gSlvsInvOffset->SetPointError(i, TMath::Abs(err0/p0/p0), totError);
 413
 414           // Only makes sense, if valid p0 available
 415           if (((Int_t)fitResSigma) != 0) {
 416             printf("Sigma fit failed!\n");
 417           }
 418           else {
 419             Double_t pSigma0 = fitResSigma->GetParams()[0];
 420             Double_t errSigma0 = fitResSigma->GetErrors()[0];
 421
 422             Double_t pSigma1 = fitResSigma->GetParams()[1];
 423             Double_t errSigma1 = fitResSigma->GetErrors()[1];
 424
 425             printf("Result sigma: pSigma1 / pSigma0 = %f / %f = %e\n\n", pSigma1, pSigma0, (pSigma0 != 0 ? pSigma1 / pSigma0 : -99999));
 426
 427             Double_t totErrorSigma = 99999;
 428             if (pSigma0 != 0 && pSigma1 != 0) {
 429               totErrorSigma = TMath::Abs(pSigma1 / pSigma0 * TMath::Sqrt(TMath::Power(errSigma0 / pSigma0, 2) + TMath::Power(errSigma1 / pSigma1, 2)));
 430             }
 431
 432             if (pSigma0 > 0) {
 433               gSigmaSlopevsdEdx->SetPoint(i, p0, pSigma1/pSigma0);
 434               gSigmaSlopevsdEdx->SetPointError(i, err0, totErrorSigma);
 435
 436               gSigmaSlopevsInvdEdx->SetPoint(i, 1./p0, pSigma1/pSigma0);
 437               gSigmaSlopevsInvdEdx->SetPointError(i, TMath::Abs(err0/p0/p0), totErrorSigma);
 438
 439                // Set artificially large errors for points with positive slope vs. multiplicity
 440               // (negative slope makes physically no sense and is most likely due to uncertainties)
 441               if (pSigma1 < 0) {
 442                 gSigmaSlopevsdEdx->SetPointError(i, err0, 1000);
 443                 gSigmaSlopevsInvdEdx->SetPointError(i, TMath::Abs(err0/p0/p0), 1000);
 444               }
 445             }
 446           }
 447         }
 448       }
 449
 450       fSave->cd("allTanTheta");
 451       c->Write();
 452       hSigma->Write();
 453       hSigmaRel->Write();
 454       arr.Clear();
 455     }
 456
 457     gSlvsOffset->SetName("slopes_AllEta");
 458     gSlvsOffset->SetTitle("Rel. slope vs. offset (all #eta)");
 459     gSlvsOffset->SetMarkerColor(kBlack);
 460     gSlvsOffset->SetLineColor(kBlack);
 461     gSlvsOffset->SetMarkerStyle(20);
 462     gSlvsOffset->Draw("Ap");
 463     gSlvsOffset->GetHistogram()->GetXaxis()->SetTitle("dEdx (offset)");
 464     gSlvsOffset->GetHistogram()->GetYaxis()->SetTitle("Multiplicity slope/offset");
 465     gSlvsOffset->Fit("pol2", "W EX0", "same", 50, 300);
 466     TF1* fitFuncResult = dynamic_cast<TF1*>(gSlvsOffset->GetListOfFunctions()->At(0));
 467     if (fitFuncResult)
 468       fitFuncResult->SetLineColor(kBlack);
 469
 470     gSlvsInvOffset->SetName("slopes2_AllEta");
 471     gSlvsInvOffset->SetTitle("Rel. slope vs. 1./offset (all #eta)");
 472     gSlvsInvOffset->SetMarkerColor(kBlack);
 473     gSlvsInvOffset->SetLineColor(kBlack);
 474     gSlvsInvOffset->SetMarkerStyle(20);
 475     gSlvsInvOffset->Draw("Ap");
 476     gSlvsInvOffset->GetHistogram()->GetXaxis()->SetTitle("1./dEdx (1./offset)");
 477     gSlvsInvOffset->GetHistogram()->GetYaxis()->SetTitle("Multiplicity slope/offset");
 478
 479     TFitResultPtr fitRes = gSlvsInvOffset->Fit("pol2", "S 0 W EX0", "", 0, 0.02);//isMC ? MCfitThreshold : 0., 0.018); // Estimate pol2 params for subsequent fit
 480     //TODO NEW Introduced parameter [4]
 481     TF1 fitFunc("fitFunc", "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)", 0., 0.2);
 482     //TF1 fitFunc("fitFunc", "[0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2)", 0., 0.1);
 483
 484
 485     //TF1* fitFuncResult2 = dynamic_cast<TF1*>(gSlvsInvOffset->GetListOfFunctions()->At(0));
 486     if ((Int_t)fitRes == 0)
 487       fitFunc.SetParameters(fitRes->GetParams());
 488     fitFunc.SetParameter(3, isMC ? 0.1 : 0.018);
 489     fitFunc.SetParLimits(3, isMC ? 0.014 : 0.01, 0.2);//TODO NEW
 490     // No lower threshold for data
 491     if (isMC) {
 492       fitFunc.SetParameter(4, MCfitThreshold);
 493       fitFunc.SetParLimits(4, 0., 0.01);
 494     }
 495     else
 496       fitFunc.FixParameter(4, 0.);
 497
 498     fitFunc.SetLineColor(kBlack);
 499
 500     gSlvsInvOffset->Fit(&fitFunc, "EX0", "same", /*TODO NEW isMC ? MCfitThreshold : */0., 0.2);
 501
 502
 503
 504     gSigmaSlopevsdEdx->SetName("sigmaSlopes_AllEta");
 505     gSigmaSlopevsdEdx->SetTitle("Sigma rel slope vs. dEdx offset (all #eta)");
 506     gSigmaSlopevsdEdx->SetMarkerColor(kBlack);
 507     gSigmaSlopevsdEdx->SetLineColor(kBlack);
 508     gSigmaSlopevsdEdx->SetMarkerStyle(20);
 509     gSigmaSlopevsdEdx->Draw("Ap");
 510     gSigmaSlopevsdEdx->GetHistogram()->GetXaxis()->SetTitle("dEdx (offset)");
 511     gSigmaSlopevsdEdx->GetHistogram()->GetYaxis()->SetTitle("Multiplicity sigma rel. slope");
 512
 513     gSigmaSlopevsInvdEdx->SetName("sigmaSlopes2_AllEta");
 514     gSigmaSlopevsInvdEdx->SetTitle("Sigma rel slope vs. 1./dEdx offset (all #eta)");
 515     gSigmaSlopevsInvdEdx->SetMarkerColor(kBlack);
 516     gSigmaSlopevsInvdEdx->SetLineColor(kBlack);
 517     gSigmaSlopevsInvdEdx->SetMarkerStyle(20);
 518     gSigmaSlopevsInvdEdx->Draw("Ap");
 519     gSigmaSlopevsInvdEdx->GetHistogram()->GetXaxis()->SetTitle("1./dEdx (1./offset)");
 520     gSigmaSlopevsInvdEdx->GetHistogram()->GetYaxis()->SetTitle("Multiplicity sigma rel. slope");
 521
 522     TFitResultPtr fitRes2 = gSigmaSlopevsInvdEdx->Fit("pol2", "S 0 W EX0", "same", 0., 0.012); // Estimate pol2 params for subsequent fit
 523     TF1 fitFuncSigma("fitFuncSigma_AllEta", "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2);
 524     if ((Int_t)fitRes2 == 0)
 525       fitFuncSigma.SetParameters(fitRes2->GetParams());
 526     fitFuncSigma.SetParameter(3, 0.012);
 527     fitFuncSigma.SetParLimits(3, 0.01, 0.2);
 528
 529     fitFuncSigma.SetLineColor(kBlack);
 530
 531     gSigmaSlopevsInvdEdx->Fit(&fitFuncSigma, "EX0", "same", 0., 0.2);
 532
 533
 534     fSave->cd("allTanTheta");
 535     gSlvsOffset->Write();
 536     gSlvsInvOffset->Write();
 537
 538     gSigmaSlopevsdEdx->Write();
 539     gSigmaSlopevsInvdEdx->Write();
 540   }
 541
 542
 543   TGraphErrors* gSlvsTanTheta[numPbins];
 544   for (Int_t j = 0; j < numPbins; j++) {
 545     gSlvsTanTheta[j] = new TGraphErrors(numTanThetaBins);
 546     gSlvsTanTheta[j]->SetName(Form("tanThetaSlopes_pBin%d", j));
 547     gSlvsTanTheta[j]->SetTitle(Form("Rel. slope vs. tanTheta (%.4f #leq pTPC/(GeV/c) < %.4f)", pBinEdges[2*j], pBinEdges[2*j+1]));
 548     gSlvsTanTheta[j]->SetMarkerColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 549     gSlvsTanTheta[j]->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 550     gSlvsTanTheta[j]->SetMarkerStyle(20);
 551   }
 552
 553   TGraphErrors* gSigmaSlopevsTanTheta[numPbins];
 554   for (Int_t j = 0; j < numPbins; j++) {
 555     gSigmaSlopevsTanTheta[j] = new TGraphErrors(numTanThetaBins);
 556     gSigmaSlopevsTanTheta[j]->SetName(Form("tanThetaSigmaSlopes_pBin%d", j));
 557     gSigmaSlopevsTanTheta[j]->SetTitle(Form("Rel. sigma slope vs. tanTheta (%.4f #leq pTPC/(GeV/c) < %.4f)", pBinEdges[2*j], pBinEdges[2*j+1]));
 558     gSigmaSlopevsTanTheta[j]->SetMarkerColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 559     gSigmaSlopevsTanTheta[j]->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 560     gSigmaSlopevsTanTheta[j]->SetMarkerStyle(20);
 561   }
 562
 563   printf("\n\n\n");
 564
 565   // Fitting for eta bins
 566   for (Int_t j = 0; j < numTanThetaBins; j++) {
 567     printf("%.2f < tanTheta <= %.2f:\n", tanThetaBinEdges[2*j], tanThetaBinEdges[2*j+1]); //TODO absTanTheta
 568
 569     TGraphErrors* gSlvsOffset = new TGraphErrors(numPbins);
 570     TGraphErrors* gSlvsInvOffset = new TGraphErrors(numPbins);
 571
 572     TGraphErrors* gSigmaSlopevsdEdx = new TGraphErrors(numPbins);
 573     TGraphErrors* gSigmaSlopevsInvdEdx = new TGraphErrors(numPbins);
 574
 575     fSave->mkdir(Form("tanThetaBin%d", j));
 576
 577     for (Int_t i = 0; i < numPbins; i++) {
 578       gSlvsOffset->SetPoint(i, -10, 0);
 579       gSlvsOffset->SetPointError(i, 0, 0);
 580
 581       gSlvsInvOffset->SetPoint(i, -10, 0);
 582       gSlvsInvOffset->SetPointError(i, 0, 0);
 583
 584       gSlvsTanTheta[i]->SetPoint(j, -10, 0);
 585       gSlvsTanTheta[i]->SetPointError(j, 0, 0);
 586
 587       gSigmaSlopevsTanTheta[i]->SetPoint(j, -10, 0);
 588       gSigmaSlopevsTanTheta[i]->SetPointError(j, 0, 0);
 589
 590       gSigmaSlopevsdEdx->SetPoint(i, -10, 0);
 591       gSigmaSlopevsdEdx->SetPointError(i, 0, 0);
 592
 593       gSigmaSlopevsInvdEdx->SetPoint(i, -10, 0);
 594       gSigmaSlopevsInvdEdx->SetPointError(i, 0, 0);
 595
 596       if (h2[i][j]->GetEntries() < 10)
 597         continue;
 598
 599       printf("Momentum %.2f - %.2f GeV/c:\n", pBinEdges[2*i], pBinEdges[2*i+1]);
 600       TCanvas* c = new TCanvas(Form("canv_pBin%d_tanThetaBin%d", i, j), Form("canv_pBin%d_tanThetaBin_%d", i, j), 100,10,1380,800);
 601       FitSlicesY(h2[i][j], heightFractionForRange, 10, "QN", &arr);
 602       h2[i][j]->GetYaxis()->SetRange(h2[i][j]->FindFirstBinAbove(1, 2), h2[i][j]->FindLastBinAbove(1,2));
 603       h2[i][j]->Draw("colz");
 604
 605       TH1D* hMean = (TH1D*)arr.At(1);//new TH1D((TH1D&)*arr.At(1));
 606       TH1D* hSigma = (TH1D*)arr.At(2);//new TH1D((TH1D&)*arr.At(2));
 607       hMean->SetName(Form("hMean_pBin%d_tanThetaBin_%d", i, j));
 608       hSigma->SetName(Form("hSigma_pBin%d_tanThetaBin_%d", i, j));
 609       hSigma->SetLineColor(kMagenta);
 610
 611       TH1D* hSigmaRel = new TH1D(*hSigma);
 612       hSigmaRel->SetName(Form("hSigmaRel_pBin%d_tanThetaBin_%d", i, j));
 613       hSigmaRel->Divide(hMean);
 614       hSigmaRel->SetLineColor(kMagenta + 2);
 615       TFitResultPtr fitRes = hMean->Fit("pol1", "s", "same", 0, 12000);
 616       TFitResultPtr fitResSigma = hSigmaRel->Fit("pol1", "s", "same", 0, 12000);
 617
 618       hMean->DrawClone("same");
 619
 620       if (((Int_t)fitRes) != 0) {
 621         printf("Fit failed!\n");
 622       }
 623       else {
 624         Double_t p0 = fitRes->GetParams()[0];
 625         Double_t p1 = fitRes->GetParams()[1];
 626         Double_t err0 = fitRes->GetErrors()[0];
 627         Double_t err1 = fitRes->GetErrors()[1];
 628
 629         Double_t totError = 99999;
 630         if (p0 != 0 && p1 != 0) {
 631           totError = TMath::Abs(p1 / p0 * TMath::Sqrt(TMath::Power(err0 / p0, 2) + TMath::Power(err1 / p1, 2)));
 632         }
 633         printf("Result: p1 / p0 = %f / %f = %e\n\n", p1, p0, (p0 != 0 ? p1 / p0 : -99999));
 634         if (p0 > 0) {
 635           gSlvsOffset->SetPoint(i, p0, p1/pSigma0);
 636           gSlvsOffset->SetPointError(i, err0, totError);
 637
 638           gSlvsInvOffset->SetPoint(i, 1./p0, p1/pSigma0);
 639           gSlvsInvOffset->SetPointError(i, TMath::Abs(err0/p0/p0), totError);
 640
 641
 642           gSlvsTanTheta[i]->SetPoint(j, (tanThetaBinEdges[2*j] + tanThetaBinEdges[2*j+1]) / 2., p1/pSigma0);
 643           gSlvsTanTheta[i]->SetPointError(j, (tanThetaBinEdges[2*j+1] - tanThetaBinEdges[2*j]) / 2., totError);
 644
 645           // Only makes sense, if valid p0 available
 646           if (((Int_t)fitResSigma) != 0) {
 647             printf("Sigma fit failed!\n");
 648           }
 649           else {
 650             Double_t pSigma0 = fitResSigma->GetParams()[0];
 651             Double_t errSigma0 = fitResSigma->GetErrors()[0];
 652
 653             Double_t pSigma1 = fitResSigma->GetParams()[1];
 654             Double_t errSigma1 = fitResSigma->GetErrors()[1];
 655
 656             printf("Result sigma: pSigma1 / pSigma0 = %f / %f = %e\n\n", pSigma1, pSigma0, (pSigma0 != 0 ? pSigma1 / pSigma0 : -99999));
 657
 658             Double_t totErrorSigma = 99999;
 659             if (pSigma0 != 0 && pSigma1 != 0) {
 660               totErrorSigma = TMath::Abs(pSigma1 / pSigma0 * TMath::Sqrt(TMath::Power(errSigma0 / pSigma0, 2) + TMath::Power(errSigma1 / pSigma1, 2)));
 661             }
 662
 663             if (pSigma0 > 0) {
 664               gSigmaSlopevsdEdx->SetPoint(i, p0, pSigma1/pSigma0);
 665               gSigmaSlopevsdEdx->SetPointError(i, err0, totErrorSigma);
 666
 667               gSigmaSlopevsInvdEdx->SetPoint(i, 1./p0, pSigma1/pSigma0);
 668               gSigmaSlopevsInvdEdx->SetPointError(i, TMath::Abs(err0/p0/p0), totErrorSigma);
 669
 670               gSigmaSlopevsTanTheta[i]->SetPoint(j, (tanThetaBinEdges[2*j] + tanThetaBinEdges[2*j+1]) / 2., pSigma1/pSigma0);
 671               gSigmaSlopevsTanTheta[i]->SetPointError(j, (tanThetaBinEdges[2*j+1] - tanThetaBinEdges[2*j]) / 2., totErrorSigma);
 672
 673               // Set artificially large errors for points with positive slope vs. multiplicity
 674               // (negative slope makes physically no sense and is most likely due to uncertainties)
 675               if (pSigma1 < 0) {
 676                 gSigmaSlopevsdEdx->SetPointError(i, err0, 1000);
 677                 gSigmaSlopevsInvdEdx->SetPointError(i, TMath::Abs(err0/p0/p0), 1000);
 678                 gSigmaSlopevsTanTheta[i]->SetPointError(j, (tanThetaBinEdges[2*j+1] - tanThetaBinEdges[2*j]) / 2., 1000);
 679               }
 680             }
 681           }
 682         }
 683       }
 684
 685       fSave->cd(Form("tanThetaBin%d", j));
 686       c->Write();
 687       hSigma->Write();
 688       hSigmaRel->Write();
 689       arr.Clear();
 690     }
 691
 692
 693     gSlvsOffset->SetName(Form("slopes_tanTheta%d", j));
 694     gSlvsOffset->SetTitle(Form("Rel. slope vs. offset (%.2f < tanTheta #leq %.2f)", tanThetaBinEdges[2*j], tanThetaBinEdges[2*j+1]));//TODO tanThetaAbs
 695     gSlvsOffset->SetMarkerColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));//   (2 + ((j > 2) ? j + 1 : j));
 696     gSlvsOffset->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));//(2 + ((j > 2) ? j + 1 : j));
 697     gSlvsOffset->SetMarkerStyle(20);
 698     gSlvsOffset->Draw("Ap");
 699     gSlvsOffset->GetHistogram()->GetXaxis()->SetTitle("dEdx (offset)");
 700     gSlvsOffset->GetHistogram()->GetYaxis()->SetTitle("Multiplicity slope/offset");
 701     gSlvsOffset->Fit("pol2", "W Ex0", "same", 10, 2000);
 702     TF1* fitFuncResult = dynamic_cast<TF1*>(gSlvsOffset->GetListOfFunctions()->At(0));
 703     if (fitFuncResult)
 704       fitFuncResult->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));//(2 + ((j > 2) ? j + 1 : j));
 705
 706     gSlvsInvOffset->SetName(Form("slopes2_tanTheta%d", j));
 707     gSlvsInvOffset->SetTitle(Form("Rel. slope vs. 1./offset (%.2f < tanTheta #leq %.2f)", tanThetaBinEdges[2*j], tanThetaBinEdges[2*j+1]));//TODO tanThetaAbs
 708     gSlvsInvOffset->SetMarkerColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));//(2 + ((j > 2) ? j + 1 : j));
 709     gSlvsInvOffset->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));//(2 + ((j > 2) ? j + 1 : j));
 710     gSlvsInvOffset->SetMarkerStyle(20);
 711     gSlvsInvOffset->Draw("Ap");
 712     gSlvsInvOffset->GetHistogram()->GetXaxis()->SetTitle("1./dEdx (1./offset)");
 713     gSlvsInvOffset->GetHistogram()->GetYaxis()->SetTitle("Multiplicity slope/offset");
 714
 715     TFitResultPtr fitRes = gSlvsInvOffset->Fit("pol2", "S 0 W EX0", "same", isMC ? MCfitThreshold : 0., 0.018); // Estimate pol2 params for subsequent fit
 716     //TODO NEW Introduced parameter [4]
 717     TF1 fitFunc(Form("fitFunc_tanTheta%d", j), "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)", 0., 0.2);
 718     //TF1 fitFunc(Form("fitFunc_tanTheta%d", j), "[0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2)", 0., 0.1);
 719     if ((Int_t)fitRes == 0)
 720       fitFunc.SetParameters(fitRes->GetParams());
 721     fitFunc.SetParameter(3, isMC ? 0.1 : 0.018);
 722     fitFunc.SetParLimits(3, isMC ? 0.014 : 0.01, 0.2);//TODO NEW
 723     // No lower threshold for data
 724     if (isMC) {
 725       fitFunc.SetParameter(4, MCfitThreshold);
 726       fitFunc.SetParLimits(4, 0., 0.01);
 727     }
 728     else
 729       fitFunc.FixParameter(4, 0.);
 730
 731     fitFunc.SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));//(2 + ((j > 2) ? j + 1 : j));
 732
 733     gSlvsInvOffset->Fit(&fitFunc, "EX0", "same", /*TODO NEW isMC ? MCfitThreshold : */0., 0.2);
 734
 735
 736
 737     gSigmaSlopevsdEdx->SetName(Form("sigmaSlopes_tanTheta%d", j));
 738     gSigmaSlopevsdEdx->SetTitle(Form("Sigma rel slope vs. dEdx offset (%.2f < tanTheta #leq %.2f)", tanThetaBinEdges[2*j],
 739                                      tanThetaBinEdges[2*j+1]));//TODO tanThetaAbs
 740     gSigmaSlopevsdEdx->SetMarkerColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 741     gSigmaSlopevsdEdx->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 742     gSigmaSlopevsdEdx->SetMarkerStyle(20);
 743     gSigmaSlopevsdEdx->Draw("Ap");
 744     gSigmaSlopevsdEdx->GetHistogram()->GetXaxis()->SetTitle("dEdx (offset)");
 745     gSigmaSlopevsdEdx->GetHistogram()->GetYaxis()->SetTitle("Multiplicity sigma rel. slope");
 746
 747     gSigmaSlopevsInvdEdx->SetName(Form("sigmaSlopes2_tanTheta%d", j));
 748     gSigmaSlopevsInvdEdx->SetTitle(Form("Sigma rel slope vs. 1./dEdx offset (%.2f < tanTheta #leq %.2f)", tanThetaBinEdges[2*j],
 749                                         tanThetaBinEdges[2*j+1]));//TODO tanThetaAbs
 750     gSigmaSlopevsInvdEdx->SetMarkerColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 751     gSigmaSlopevsInvdEdx->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 752     gSigmaSlopevsInvdEdx->SetMarkerStyle(20);
 753     gSigmaSlopevsInvdEdx->Draw("Ap");
 754     gSigmaSlopevsInvdEdx->GetHistogram()->GetXaxis()->SetTitle("1./dEdx (1./offset)");
 755     gSigmaSlopevsInvdEdx->GetHistogram()->GetYaxis()->SetTitle("Multiplicity sigma rel. slope");
 756
 757
 758     TFitResultPtr fitRes2 = gSigmaSlopevsInvdEdx->Fit("pol2", "S 0 W EX0", "same", 0., 0.012); // Estimate pol2 params for subsequent fit
 759     TF1 fitFuncSigma(Form("fitFuncSigma_tanTheta%d", j), "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))",
 760                      0., 0.2);
 761     if ((Int_t)fitRes2 == 0)
 762       fitFuncSigma.SetParameters(fitRes2->GetParams());
 763     fitFuncSigma.SetParameter(3, 0.012);
 764     fitFuncSigma.SetParLimits(3, 0.006, 0.2); //TODO was 3, 0.01, 0.2
 765
 766     fitFuncSigma.SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));//(2 + ((j > 2) ? j + 1 : j));
 767
 768     gSigmaSlopevsInvdEdx->Fit(&fitFuncSigma, "EX0", "same", 0., 0.2);
 769
 770
 771     fSave->cd(Form("tanThetaBin%d", j));
 772     gSlvsOffset->Write();
 773     gSlvsInvOffset->Write();
 774
 775     gSigmaSlopevsdEdx->Write();
 776     gSigmaSlopevsInvdEdx->Write();
 777
 778     printf("\n\n\n");
 779   }
 780
 781
 782
 783
 784   fSave->mkdir("tanThetaFits");
 785   fSave->cd("tanThetaFits");
 786
 787   for (Int_t j = 0; j < numPbins; j++) {
 788     /*
 789     // Scale graphs, just that first bin sits at +-1
 790     Double_t scale = 1.0;
 791     for (Int_t k = 0; k < gSlvsTanTheta[j]->GetN(); k++) {
 792       if (k == 0) {
 793         scale = 1. / TMath::Max(1e-10, TMath::Abs(gSlvsTanTheta[j]->GetY()[k]));
 794       }
 795
 796       gSlvsTanTheta[j]->SetPoint(k, gSlvsTanTheta[j]->GetX()[k], gSlvsTanTheta[j]->GetY()[k] * scale);
 797       gSlvsTanTheta[j]->SetPointError(k, gSlvsTanTheta[j]->GetEX()[k], gSlvsTanTheta[j]->GetEY()[k] * scale);
 798     }
 799     */
 800
 801     // Shift graphs, just that reference bin sits at 0
 802     Double_t shift = 0.0;
 803     const Int_t refBin = 7;
 804     shift = gSlvsTanTheta[j]->GetY()[refBin];
 805
 806     for (Int_t k = 0; k < gSlvsTanTheta[j]->GetN(); k++) {
 807       gSlvsTanTheta[j]->SetPoint(k, gSlvsTanTheta[j]->GetX()[k], gSlvsTanTheta[j]->GetY()[k] - shift);
 808     }
 809
 810     gSlvsTanTheta[j]->Draw("Ap");
 811     gSlvsTanTheta[j]->GetHistogram()->GetXaxis()->SetTitle("tan(#Theta)");//TODO tanThetaAbs
 812     gSlvsTanTheta[j]->GetHistogram()->GetYaxis()->SetTitle("Multiplicity slope/offset");
 813
 814     gSlvsTanTheta[j]->Fit("pol2", "EX0", "same", tanThetaBinEdges[0], tanThetaBinEdges[numTanThetaBins * 2 - 1]);
 815
 816     TF1* fitFuncResult = dynamic_cast<TF1*>(gSlvsTanTheta[j]->GetListOfFunctions()->At(0));
 817     if (fitFuncResult)
 818       fitFuncResult->SetLineColor((j == 7) ? kBlack : (2 + ((j > 2) ? j + 1 : j)));
 819     gSlvsTanTheta[j]->Write();
 820
 821     gSigmaSlopevsTanTheta[j]->Draw("Ap");
 822     gSigmaSlopevsTanTheta[j]->GetHistogram()->GetXaxis()->SetTitle("tanTheta");//TODO tanThetaAbs
 823     gSigmaSlopevsTanTheta[j]->GetHistogram()->GetYaxis()->SetTitle("Multiplicity rel. sigma slope");
 824     gSigmaSlopevsTanTheta[j]->Write();
 825   }
 826
 827   fSave->cd();
 828   if (correct)
 829     corrFunc.Write();
 830
 831
 832   TNamed* settings = new TNamed(Form("Settings: widthFactor %f, multStepSize %d, isMC %d, correct %d, pathNameSplines %s, splineName %s\n",
 833                                      widthFactor, multStepSize, isMC, correct, pathNameSplines.Data(), splineName.Data()), "");
 834   settings->Write();
 835   fSave->Close();
 836
 837
 838   return 0;
 839 }