- Added classes and macros for TPC PID calibration
[u/mrichter/AliRoot.git] / 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 }