- Added classes and macros for TPC PID calibration
[u/mrichter/AliRoot.git] / PWGPP / TPC / macros / PIDCalib / checkPullTree.C
1 #include "TTree.h"
2 #include "TH2D.h"
3 #include "TH3D.h"
4 #include "TH3F.h"
5 #include "TCanvas.h"
6 #include "TStyle.h"
7 #include "TList.h"
8 #include "TString.h"
9 #include "TFile.h"
10 #include "TMath.h"
11 #include "TDatime.h"
12 #include "TRandom3.h"
13 #include "TSpline.h"
14
15 #include "TF1.h" //TODO NOW
16
17 #include <iostream>
18 #include <iomanip>
19
20 #include "THnSparseDefinitions.h"
21
22 //__________________________________________________________________________________________
23 void normaliseHisto(TH2D* h)
24 {
25   h->Sumw2();
26   //TODO NOW: Disabled normalisation (doesn't make sense here)
27   /*
28   for (Int_t x = 1; x <= h->GetXaxis()->GetNbins(); x++) {
29     Double_t binWidth = h->GetXaxis()->GetBinWidth(x);
30     for (Int_t y = 1; y <= h->GetYaxis()->GetNbins(); y++) {
31       Double_t binError = h->GetBinError(x,y);
32       h->SetBinContent(x, y, h->GetBinContent(x, y) / binWidth);
33       h->SetBinError(x, y, binError / binWidth);
34     }
35   }*/
36 }
37
38
39 //_________________________________________________________________________________________
40 Int_t getBinX(TH2D* h, Double_t tanTheta)
41 {
42   Int_t binX = h->GetXaxis()->FindBin(tanTheta);
43   
44   if (binX <= 0)
45     binX = 1;
46   if (binX > h->GetXaxis()->GetNbins())
47     binX = h->GetXaxis()->GetNbins();
48
49     return binX;
50 }
51
52
53 //_________________________________________________________________________________________
54 Int_t getBinY(TH2D* h, Double_t dEdxInv)
55 {
56   Int_t binY = h->GetYaxis()->FindBin(dEdxInv);
57   
58   if (binY <= 0)
59     binY = 1;
60   if (binY > h->GetYaxis()->GetNbins())
61     binY = h->GetYaxis()->GetNbins();
62
63   return binY;
64 }
65
66
67 //_________________________________________________________________________________________
68 Int_t checkPullTree(TString pathTree,  TString pathNameThetaMap, TString pathNameSigmaMap,
69                     TString mapSuffix, const Int_t collType /*0: pp, 1: pPb, 2: PbPb*/, const Bool_t plotPull = kTRUE,
70                     const Double_t downScaleFactor = 1,
71                     TString pathNameSplinesFile = "", TString prSplinesName = "",
72                     TString fileNameTree = "bhess_PIDetaTree.root", TString treeName = "fTree")
73 {
74   const Bool_t isNonPP = collType != 0;
75   const Double_t massProton = AliPID::ParticleMass(AliPID::kProton);
76   
77   Bool_t recalculateExpecteddEdx = pathNameSplinesFile != "";
78   
79   TFile* f = 0x0;
80         
81   f = TFile::Open(Form("%s/%s", pathTree.Data(), fileNameTree.Data()));
82   if (!f)  {
83     std::cout << "Failed to open tree file \"" << Form("%s/%s", pathTree.Data(), fileNameTree.Data()) << "\"!" << std::endl;
84     return -1;
85   }
86       
87   // Extract the data Tree
88   TTree* tree = dynamic_cast<TTree*>(f->Get(treeName.Data()));
89   if (!tree) {
90     std::cout << "Failed to load data tree!" << std::endl;
91     return -1;
92   }
93   
94   // Extract the splines, if desired
95   TSpline3* splPr = 0x0;
96   if (recalculateExpecteddEdx) {
97     std::cout << "Loading splines to recalculate expected dEdx!" << std::endl << std::endl;
98     
99     TFile* fSpl = TFile::Open(pathNameSplinesFile.Data());
100     if (!fSpl) {
101       std::cout << "Failed to open spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
102       return 0x0;
103     }
104     
105     TObjArray* TPCPIDResponse = (TObjArray*)fSpl->Get("TPCPIDResponse");
106     if (!TPCPIDResponse) {
107       splPr = (TSpline3*)fSpl->Get(prSplinesName.Data());
108       
109       // If splines are in file directly, without TPCPIDResponse object, try to load them
110       if (!splPr) {
111         std::cout << "Failed to load object array from spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
112         return 0x0;
113       }
114     }
115     else {
116       splPr = (TSpline3*)TPCPIDResponse->FindObject(prSplinesName.Data());
117       
118       if (!splPr) {
119         std::cout << "Failed to load splines from file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
120         return 0x0;
121       }
122     }
123   }
124   else
125     std::cout << "Taking dEdxExpected from Tree..." << std::endl << std::endl;
126
127   // Extract the correction maps
128   TFile* fMap = TFile::Open(pathNameThetaMap.Data());
129   if (!fMap)  {
130     std::cout << "Failed to open thetaMap file \"" << pathNameThetaMap.Data() << "\"! Will not additionally correct data...." << std::endl;
131   }
132
133   TH2D* hMap = 0x0;
134   
135   if (fMap) {
136     hMap = dynamic_cast<TH2D*>(fMap->Get(Form("hRefined%s", mapSuffix.Data())));
137     if (!hMap) {
138       std::cout << "Failed to load theta map!" << std::endl;
139       return -1;
140     }
141   }
142
143   TFile* fSigmaMap = TFile::Open(pathNameSigmaMap.Data());
144   if (!fSigmaMap)  {
145     std::cout << "Failed to open simgaMap file \"" << pathNameSigmaMap.Data() << "\"!" << std::endl;
146     return -1;
147   }
148
149   TH2D* hThetaMapSigmaPar1 = dynamic_cast<TH2D*>(fSigmaMap->Get("hThetaMapSigmaPar1"));
150   if (!hThetaMapSigmaPar1) {
151     std::cout << "Failed to load sigma map for par 1!" << std::endl;
152     return -1;
153   }
154
155   Double_t c0 = -1;
156   TNamed* c0Info = dynamic_cast<TNamed*>(fSigmaMap->Get("c0"));
157   if (!c0Info) {
158     std::cout << "Failed to extract c0 from file with sigma map!" << std::endl;
159     return -1;
160   }
161
162   TString c0String = c0Info->GetTitle();
163   c0 = c0String.Atof();
164   printf("Loaded parameter 0 for sigma: %f\n\n", c0);
165
166   if (plotPull)
167     std::cout << "Plotting pull..." << std::endl << std::endl;
168   else
169     std::cout << "Plotting delta'..." << std::endl << std::endl;
170
171   Long64_t nTreeEntries = tree->GetEntriesFast();
172
173   Double_t dEdx = 0.; // Measured dE/dx
174   Double_t dEdxExpected = 0.; // Expected dE/dx according to parametrisation
175   Double_t tanTheta = 0.; // Tangens of (local) theta at TPC inner wall
176   Double_t pTPC = 0.; // Momentum at TPC inner wall
177   UShort_t tpcSignalN = 0; // Number of clusters used for dEdx
178   UChar_t  pidType = 0;
179   Int_t    fMultiplicity = 0;
180   //Double_t phiPrime = 0;
181
182   // Only activate the branches of interest to save processing time
183   tree->SetBranchStatus("*", 0); // Disable all branches
184   tree->SetBranchStatus("pTPC", 1);
185   tree->SetBranchStatus("dEdx", 1);
186   tree->SetBranchStatus("dEdxExpected", 1);
187   tree->SetBranchStatus("tanTheta", 1);
188   tree->SetBranchStatus("tpcSignalN", 1);
189   tree->SetBranchStatus("pidType", 1);
190   //tree->SetBranchStatus("phiPrime", 1);
191   if (isNonPP)
192     tree->SetBranchStatus("fMultiplicity", 1);
193
194   
195   tree->SetBranchAddress("dEdx", &dEdx);
196   tree->SetBranchAddress("dEdxExpected", &dEdxExpected);
197   tree->SetBranchAddress("tanTheta", &tanTheta);
198   tree->SetBranchAddress("tpcSignalN", &tpcSignalN);
199   tree->SetBranchAddress("pTPC", &pTPC);
200   tree->SetBranchAddress("pidType", &pidType);
201   //tree->SetBranchAddress("phiPrime", &phiPrime);
202   if (isNonPP)
203     tree->SetBranchAddress("fMultiplicity", &fMultiplicity);
204
205   
206   // Output file
207   TDatime daTime;
208   TString savefileName = Form("%s%s_checkPullSigma_%04d_%02d_%02d__%02d_%02d.root", fileNameTree.ReplaceAll(".root", "").Data(),
209                               recalculateExpecteddEdx ? "_recalcdEdx" : "",
210                               daTime.GetYear(), daTime.GetMonth(), daTime.GetDay(), daTime.GetHour(), daTime.GetMinute());
211
212   TFile* fSave = TFile::Open(Form("%s/%s", pathTree.Data(), savefileName.Data()), "recreate");
213   if (!fSave) {
214     std::cout << "Failed to open save file \"" << Form("%s/%s", pathTree.Data(), savefileName.Data()) << "\"!" << std::endl;
215     return -1;
216   }
217   
218   const Double_t pBoundLow = 0.1;
219   const Double_t pBoundUp = 5;
220
221   const Int_t nBins1 = TMath::Ceil(180 / downScaleFactor);
222   const Int_t nBins2 = TMath::Ceil(100 / downScaleFactor);
223   const Int_t nBins3 = TMath::Ceil(60 / downScaleFactor);
224   
225   const Int_t nPbinsForMap = nBins1 + nBins2 + nBins3;
226   Double_t binsPforMap[nPbinsForMap + 1];
227   
228   Double_t binWidth1 = (1.0 - pBoundLow) / nBins1;
229   Double_t binWidth2 = (2.0 - 1.0 ) / nBins2;
230   Double_t binWidth3 = (pBoundUp - 2.0) / nBins3;
231   
232   for (Int_t i = 0; i < nBins1; i++)  {
233     binsPforMap[i] = pBoundLow + i * binWidth1;
234   }
235   for (Int_t i = nBins1, j = 0; i < nBins1 + nBins2; i++, j++)  {
236     binsPforMap[i] = 1.0 + j * binWidth2;
237   }
238   for (Int_t i = nBins1 + nBins2, j = 0; i < nBins1 + nBins2 + nBins3; i++, j++)  {
239     binsPforMap[i] = 2.0 + j * binWidth3;
240   }
241   binsPforMap[nPbinsForMap] = pBoundUp;
242
243   TH2D* hPull = new TH2D("hPull", "Pull vs. p_{TPC} integrated over tan(#Theta);p_{TPC} (GeV/c);Pull", nPbinsForMap, binsPforMap, 
244                          plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);
245   TH2D* hPullAdditionalCorr = (TH2D*)hPull->Clone("hPullAdditionalCorr");
246   hPullAdditionalCorr->SetTitle("Pull vs. p_{TPC} integrated over tan(#Theta) with additional dEdx correction w.r.t. tan(#Theta)");
247   /*
248   const Int_t nThetaHistos = 3;
249   TH2D* hPullTheta[nThetaHistos];
250   TH2D* hPullAdditionalCorrTheta[nThetaHistos];
251   Double_t tThetaLow[nThetaHistos] = { 0.0, 0.4, 0.9 };
252   Double_t tThetaHigh[nThetaHistos] = { 0.1, 0.5, 1.0 };
253   */
254   const Int_t nThetaHistos = 10;
255   TH2D* hPullTheta[nThetaHistos];
256   TH2D* hPullAdditionalCorrTheta[nThetaHistos];
257   Double_t tThetaLow[nThetaHistos] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
258   Double_t tThetaHigh[nThetaHistos] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
259   
260
261   for (Int_t i = 0; i < nThetaHistos; i++)    {
262     hPullTheta[i] = new TH2D(Form("hPullTheta_%d", i),
263                              Form("Pull vs. p_{TPC} for %.2f <= |tan(#Theta)| < %.2f;p_{TPC} (GeV/c);Pull", tThetaLow[i], tThetaHigh[i]),
264                              nPbinsForMap, binsPforMap, plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);
265
266     hPullAdditionalCorrTheta[i] =
267       new TH2D(Form("hPullAdditionalCorrTheta_%d", i),
268                Form("Pull vs. p_{TPC} for %.2f <= |tan(#Theta)| < %.2f with additional dEdx correction w.r.t. tan(#Theta);p_{TPC} (GeV/c);Pull",
269                     tThetaLow[i], tThetaHigh[i]),
270                nPbinsForMap, binsPforMap, plotPull ? 120 : 240, plotPull ? -6 : -0.6, plotPull ? 6 : 0.6);
271   }
272
273   
274   
275   
276   
277   
278   TF1 corrFuncMult("corrFuncMult", "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
279                    0., 0.2);
280   TF1 corrFuncMultTanTheta("corrFuncMultTanTheta", "[0] * (x -[2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
281   TF1 corrFuncSigmaMult("corrFuncSigmaMul", "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2);
282   
283   
284   // LHC13b.pass2
285   if (isNonPP)
286     printf("Using corr Parameters for 13b.pass2\n!");
287   
288   corrFuncMult.SetParameter(0, -5.906e-06);
289   corrFuncMult.SetParameter(1, -5.064e-04);
290   corrFuncMult.SetParameter(2, -3.521e-02);
291   corrFuncMult.SetParameter(3,  2.469e-02);
292   corrFuncMult.SetParameter(4, 0);
293   
294   corrFuncMultTanTheta.SetParameter(0, -5.32e-06);
295   corrFuncMultTanTheta.SetParameter(1,  1.177e-05);
296   corrFuncMultTanTheta.SetParameter(2, -0.5);
297   
298   corrFuncSigmaMult.SetParameter(0, 0.);
299   corrFuncSigmaMult.SetParameter(1, 0.);
300   corrFuncSigmaMult.SetParameter(2, 0.);
301   corrFuncSigmaMult.SetParameter(3, 0.);
302   
303   
304   /* OK, but PID task was not very satisfying
305   corrFuncMult.SetParameter(0, -6.27187e-06);
306   corrFuncMult.SetParameter(1, -4.60649e-04);
307   corrFuncMult.SetParameter(2, -4.26450e-02);
308   corrFuncMult.SetParameter(3, 2.40590e-02);
309   corrFuncMult.SetParameter(4, 0);
310   
311   corrFuncMultTanTheta.SetParameter(0, -5.338e-06);
312   corrFuncMultTanTheta.SetParameter(1,  1.220e-05);
313   corrFuncMultTanTheta.SetParameter(2, -0.5);
314   
315   corrFuncSigmaMult.SetParameter(0, 7.89237e-05);
316   corrFuncSigmaMult.SetParameter(1, -1.30662e-02);
317   corrFuncSigmaMult.SetParameter(2, 8.91548e-01);
318   corrFuncSigmaMult.SetParameter(3, 1.47931e-02);
319   */
320   
321   
322   /*
323   // LHC11a10a
324   if (isNonPP)
325     printf("Using corr Parameters for 11a10a\n!");
326   
327   corrFuncMult.SetParameter(0, 6.90133e-06);
328   corrFuncMult.SetParameter(1, -1.22123e-03);
329   corrFuncMult.SetParameter(2, 1.80220e-02);
330   corrFuncMult.SetParameter(3, 0.1);
331   corrFuncMult.SetParameter(4, 6.45306e-03);
332   
333   corrFuncMultTanTheta.SetParameter(0, -2.85505e-07);
334   corrFuncMultTanTheta.SetParameter(1, -1.31911e-06);
335   corrFuncMultTanTheta.SetParameter(2, -0.5);
336   
337   corrFuncSigmaMult.SetParameter(0, -4.29665e-05);
338   corrFuncSigmaMult.SetParameter(1, 1.37023e-02);
339   corrFuncSigmaMult.SetParameter(2, -6.36337e-01);
340   corrFuncSigmaMult.SetParameter(3, 1.13479e-02);
341   */
342   
343   /* OLD without saturation and large error for negative slopes
344   corrFuncSigmaMult.SetParameter(0, -4.79684e-05);
345   corrFuncSigmaMult.SetParameter(1, 1.49938e-02);
346   corrFuncSigmaMult.SetParameter(2, -7.15269e-01);
347   corrFuncSigmaMult.SetParameter(3, 1.06855e-02);
348   */
349   
350   /* OLD very good try, but with fewer pBins for the fitting
351   corrFuncMult.SetParameter(0, 6.88365e-06);
352   corrFuncMult.SetParameter(1, -1.22324e-03);
353   corrFuncMult.SetParameter(2, 1.81625e-02);
354   corrFuncMult.SetParameter(3, 0.1);
355   corrFuncMult.SetParameter(4, 6.36890e-03);
356   
357   corrFuncMultTanTheta.SetParameter(0, -2.85505e-07);
358   corrFuncMultTanTheta.SetParameter(1, -1.31911e-06);
359   corrFuncMultTanTheta.SetParameter(2, -0.5);
360   
361   corrFuncSigmaMult.SetParameter(0, -4.28401e-05);
362   corrFuncSigmaMult.SetParameter(1, 1.24812e-02);
363   corrFuncSigmaMult.SetParameter(2, -5.28531e-01);
364   corrFuncSigmaMult.SetParameter(3, 1.25147e-02);
365   */
366   /*OLD good try
367   corrFuncMult.SetParameter(0, 7.50321e-06);
368   corrFuncMult.SetParameter(1, -1.25250e-03);
369   corrFuncMult.SetParameter(2, 1.85437e-02);
370   corrFuncMult.SetParameter(3, 0.1);
371   corrFuncMult.SetParameter(4, 6.21192e-03);
372   
373   corrFuncMultTanTheta.SetParameter(0, -1.43112e-07);
374   corrFuncMultTanTheta.SetParameter(1, -1.53e-06);
375   corrFuncMultTanTheta.SetParameter(2, 0.3);
376   
377   corrFuncSigmaMult.SetParameter(0, -2.54019e-05);
378   corrFuncSigmaMult.SetParameter(1, 8.68883e-03);
379   corrFuncSigmaMult.SetParameter(2, -3.36176e-01);
380   corrFuncSigmaMult.SetParameter(3, 1.29230e-02);
381   */
382   
383   /*
384   // LHC10h.pass2
385   if (isNonPP)
386     printf("Using corr Parameters for 10h.pass2\n!");
387   
388   corrFuncMult.SetParameter(0, 3.21636e-07);
389   corrFuncMult.SetParameter(1, -6.65876e-04);
390   corrFuncMult.SetParameter(2, 1.28786e-03);
391   corrFuncMult.SetParameter(3, 1.47677e-02);
392   corrFuncMult.SetParameter(4, 0.);
393   
394   corrFuncMultTanTheta.SetParameter(0, 7.23591e-08);
395   corrFuncMultTanTheta.SetParameter(1, 2.7469e-06);
396   corrFuncMultTanTheta.SetParameter(2, -0.5);
397   
398   corrFuncSigmaMult.SetParameter(0, -1.22590e-05);
399   corrFuncSigmaMult.SetParameter(1, 6.88888e-03);
400   corrFuncSigmaMult.SetParameter(2, -3.20788e-01);
401   corrFuncSigmaMult.SetParameter(3, 1.07345e-02);
402   */
403   
404   /*OLD bad try
405   corrFuncMult.SetParameter(0, 2.71514e-07);
406   corrFuncMult.SetParameter(1, -6.92031e-04);
407   corrFuncMult.SetParameter(2, 3.56042e-03);
408   corrFuncMult.SetParameter(3, 1.47497e-02);
409   corrFuncMult.SetParameter(4, 0.);
410   
411   corrFuncMultTanTheta.SetParameter(0, 8.53204e-08);
412   corrFuncMultTanTheta.SetParameter(1, 2.85591e-06);
413   corrFuncMultTanTheta.SetParameter(2, -0.5);
414   
415   corrFuncSigmaMult.SetParameter(0, -6.82477e-06);
416   corrFuncSigmaMult.SetParameter(1, 4.97051e-03);
417   corrFuncSigmaMult.SetParameter(2, -1.64954e-01);
418   corrFuncSigmaMult.SetParameter(3, 9.21061e-03);
419   */
420   
421
422   //TODO NOW
423   TF1* fShapeSmallP = new TF1("fShapeSmallP", "pol5", -0.4, 0.4);
424   fShapeSmallP->SetParameters(1.01712, -0.0202725, -0.260692, 0.261623, 0.671854, -1.14014);
425     
426   for (Long64_t i = 0; i < nTreeEntries; i++) {
427     tree->GetEntry(i);
428
429     if (dEdx <= 0 || dEdxExpected <= 0 || tpcSignalN <= 10)
430       continue;
431     /*
432     Double_t pT = pTPC*TMath::Sin(-TMath::ATan(tanTheta)+TMath::Pi()/2.0);
433     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)) 
434       continue;
435     */
436       
437     if (pidType != kMCid) {
438       if (pidType == kTPCid && pTPC > 0.6)
439         continue;
440       if (pidType == kTPCandTOFid && (pTPC < 0.6 || pTPC > 2.0))
441         continue;
442       if ((collType == 2) && pidType == kTPCandTOFid && pTPC > 1.0)
443         continue;// Only V0's in case of PbPb above 1.0 GeV/c
444       if (pidType == kV0idPlusTOFrejected) //TODO NOW NEW
445         continue;
446     }
447     
448     if (recalculateExpecteddEdx) {
449       dEdxExpected = 50. * splPr->Eval(pTPC / massProton); //WARNING: What, if MIP is different from 50.? Seems not to be used (tested for pp, MC_pp, PbPb and MC_PbPb), but can in principle happen
450     }
451       
452     //TODO NOW
453     /*
454     if (TMath::Abs(tanTheta) <= 0.4) {
455       Double_t p0 = fShapeSmallP->Eval(tanTheta) - 1.0; // Strength of the correction
456       Double_t p1 = -9.0; // How fast the correction is turned off
457       Double_t p2 = -0.209; // Turn off correction around 0.2 GeV/c
458       Double_t p3 = 1.0; // Delta' for large p should be 1
459
460       Double_t corrFactor = TMath::Erf((pTPC + p2) * p1) * p0 + p3 + p0; // Add p0 to have 1 for p3 = 1 and large pTPC
461       dEdxExpected *= corrFactor;
462     }*/
463     
464      /*TODO old unsuccessful try 
465     Double_t thetaGlobalTPC = -TMath::ATan(tanTheta) + TMath::Pi() / 2.;
466     Double_t pTtpc = pTPC * TMath::Sin(thetaGlobalTPC);
467     Double_t pTtpcInv = (pTtpc > 0) ? 1. / pTtpc : 0;
468     Double_t p0 = 1.0;
469     Double_t p1 = 1./ 0.5;//TODO 2.0;
470     Double_t p2 = -0.2;//TODO 0.1
471     Double_t pTcorrFactor = p0 + (pTtpcInv > p1) * p2 * (pTtpcInv - p1);
472     
473     dEdxExpected *= pTcorrFactor;
474     */
475     
476       
477     // From the momentum (via dEdxExpected) and the tanTheta of the track, the expected dEdx can be calculated (correctedDeDxExpected).
478     // If the splines are correct, this should give in average the same value as dEdx. 
479     // Now valid: Maps created from corrected data with splines adopted to corrected data, so lookup should be for dEdxExpected=dEdxSplines (no further
480     // eta correction) or the corrected dEdx from the track (which should ideally be = dEdxSplines)
481     
482     // Tested with corrected data for LHC10d.pass2: using dEdx for the lookup (which is the corrected value and should ideally be = dEdxSplines):
483     // Results almost the same. Maybe slightly better for dEdxExpected.
484     
485     // No longer valid: Note that the maps take always the uncorrected dEdx w.r.t.
486     // tanTheta, so that correctedDeDxExpected is needed here normally. However, the information for the correction will be lost at some point.
487     // Therefore, dEdxExpected can be used instead and should provide a good approximation.
488     Double_t c1FromSigmaMap = hThetaMapSigmaPar1->GetBinContent(getBinX(hThetaMapSigmaPar1, tanTheta), getBinY(hThetaMapSigmaPar1, 1./dEdxExpected));
489     
490     Double_t expectedSigma = dEdxExpected * TMath::Sqrt( c0 * c0 + (c1FromSigmaMap * c1FromSigmaMap) / tpcSignalN);
491     Double_t pull = (dEdx - dEdxExpected) / (plotPull ? expectedSigma: dEdxExpected);
492
493     // Fill pull histo
494     hPull->Fill(pTPC, pull);
495
496     Double_t tanThetaAbs = TMath::Abs(tanTheta);
497     
498     for (Int_t j = 0; j < nThetaHistos; j++)    {
499       if (tanThetaAbs  >= tThetaLow[j] && tanThetaAbs < tThetaHigh[j])  {
500         hPullTheta[j]->Fill(pTPC, pull);
501       }
502     }
503
504     if (!hMap)
505       continue;
506
507     Double_t correctionFactor = 1.;
508     
509     if (isNonPP) {
510       // 1. Correct eta dependence
511       correctionFactor = hMap->GetBinContent(getBinX(hMap, tanTheta), getBinY(hMap, 1./dEdxExpected));
512       
513       // 2. Correct for multiplicity dependence:
514       Double_t multCorrectionFactor = 1.;
515       
516       if (fMultiplicity > 0) {
517         Double_t relSlope = corrFuncMult.Eval(1. / (dEdxExpected * correctionFactor));
518         relSlope += corrFuncMultTanTheta.Eval(tanTheta);
519
520         multCorrectionFactor = 1. + relSlope * fMultiplicity;
521       }
522
523       c1FromSigmaMap = hThetaMapSigmaPar1->GetBinContent(getBinX(hThetaMapSigmaPar1, tanTheta), getBinY(hThetaMapSigmaPar1, 1./dEdxExpected));
524       
525       // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity, i.e. the eta (only) corrected dEdxExpected value has to be used
526       // since all maps etc. have been created for ~zero multiplicity
527       Double_t relSigmaSlope = corrFuncSigmaMult.Eval(1. / (dEdxExpected * correctionFactor));
528       Double_t multSigmaCorrectionFactor = 1. + relSigmaSlope * fMultiplicity;
529       
530       dEdxExpected *= correctionFactor * multCorrectionFactor; 
531       
532       expectedSigma = dEdxExpected * TMath::Sqrt( c0 * c0 + (c1FromSigmaMap * c1FromSigmaMap) / tpcSignalN);
533       expectedSigma *= multSigmaCorrectionFactor;
534       
535       pull = (dEdx - dEdxExpected) / (plotPull ? expectedSigma: dEdxExpected);
536     }
537     else {
538       correctionFactor = hMap->GetBinContent(getBinX(hMap, tanTheta), getBinY(hMap, 1./dEdxExpected));
539       
540       c1FromSigmaMap = hThetaMapSigmaPar1->GetBinContent(getBinX(hThetaMapSigmaPar1, tanTheta), getBinY(hThetaMapSigmaPar1, 1./dEdxExpected));
541    
542       dEdxExpected *= correctionFactor; // If data is not corrected, but the sigma map is for corrected data, re-do analysis with corrected dEdx
543       
544       expectedSigma = dEdxExpected * TMath::Sqrt( c0 * c0 + (c1FromSigmaMap * c1FromSigmaMap) / tpcSignalN);
545       pull = (dEdx - dEdxExpected) / (plotPull ? expectedSigma: dEdxExpected);
546     }
547
548     pull = (dEdx - dEdxExpected) / (plotPull ? expectedSigma: dEdxExpected);
549
550     hPullAdditionalCorr->Fill(pTPC, pull);
551
552     for (Int_t j = 0; j < nThetaHistos; j++)    {
553       if (tanThetaAbs  >= tThetaLow[j] && tanThetaAbs < tThetaHigh[j])  {
554         hPullAdditionalCorrTheta[j]->Fill(pTPC, pull);
555       }
556     }
557   }
558 /*
559   // Mean, Sigma, chi^2/NDF of pull of different theta bins and all in one plot
560   TCanvas* canvPullMean = new TCanvas("canvPullMean", "canvPullMean", 100,10,1380,800);
561   canvPullMean->SetLogx(kTRUE);
562   canvPullMean->SetGridx(kTRUE);
563   canvPullMean->SetGridy(kTRUE);
564   TCanvas* canvPullSigma = new TCanvas("canvPullSigma", "canvPullSigma", 100,10,1380,800);
565   canvPullSigma->SetLogx(kTRUE);
566   canvPullSigma->SetGridx(kTRUE);
567   canvPullSigma->SetGridy(kTRUE);
568   TCanvas* canvPullChi2 = new TCanvas("canvPullChi2", "canvPullChi2", 100,10,1380,800);
569   canvPullChi2->SetLogx(kTRUE);
570   canvPullChi2->SetGridx(kTRUE);
571   canvPullChi2->SetGridy(kTRUE);
572   
573
574   TCanvas* canvPull[nThetaHistos + 1];
575   for (Int_t i = 0, j = nThetaHistos; i < nThetaHistos + 1; i++, j--)  {
576     canvPull[i] = new TCanvas(Form("canvPull_%d", i), "canvPull", 100,10,1380,800);
577     canvPull[i]->cd();
578     canvPull[i]->SetLogx(kTRUE);
579     canvPull[i]->SetLogz(kTRUE);
580     canvPull[i]->SetGrid(kTRUE, kTRUE);
581
582     TH2D* hTemp = 0x0;
583     TString thetaString = "";
584     if (i == nThetaHistos)  {
585       hTemp = hPull;
586       thetaString = "tan(#Theta) integrated";
587     }
588     else {
589       hTemp = hPullTheta[i];
590       thetaString = Form("%.2f #leq |tan(#Theta)| < %.2f", tThetaLow[i], tThetaHigh[i]);
591     }
592     
593     normaliseHisto(hTemp);
594     hTemp->FitSlicesY();
595     hTemp->GetYaxis()->SetNdivisions(12);
596     hTemp->GetXaxis()->SetMoreLogLabels(kTRUE);
597     TH1D* hTempMean = (TH1D*)gDirectory->Get(Form("%s_1", hTemp->GetName()));
598     hTempMean->SetTitle(Form("mean(pull), %s", thetaString.Data()));
599     hTempMean->GetXaxis()->SetMoreLogLabels(kTRUE);
600     hTempMean->SetLineWidth(2);
601     hTempMean->SetMarkerStyle(20);
602     TH1D* hTempSigma = (TH1D*)gDirectory->Get(Form("%s_2", hTemp->GetName()));
603     hTempSigma->SetTitle(Form("#sigma(pull), %s", thetaString.Data()));
604     hTempSigma->GetXaxis()->SetMoreLogLabels(kTRUE);
605     hTempSigma->SetLineColor(kMagenta);
606     hTempSigma->SetMarkerStyle(20);
607     hTempSigma->SetMarkerColor(kMagenta);
608     hTempSigma->SetLineWidth(2);
609     TH1D* hTempChi2 = (TH1D*)gDirectory->Get(Form("%s_chi2", hTemp->GetName()));
610     hTempChi2->SetTitle(Form("#chi^{2} / NDF (pull), %s", thetaString.Data()));
611     hTempChi2->GetXaxis()->SetMoreLogLabels(kTRUE);
612     hTempChi2->SetLineColor(kMagenta + 2);
613     hTempChi2->SetMarkerStyle(20);
614     hTempChi2->SetMarkerColor(kMagenta + 2);
615     hTempChi2->SetLineWidth(2);
616
617     hTemp->DrawCopy("colz");
618     hTempMean->DrawCopy("same");
619     hTempSigma->DrawCopy("same");
620     hTempChi2->Scale(-1./10.);
621     hTempChi2->DrawCopy("same");
622     hTempChi2->Scale(-10.);
623
624     canvPullMean->cd();
625     hTempMean->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
626     hTempMean->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
627     hTempMean->DrawCopy((i == 0 ? "" : "same"));
628
629     canvPullSigma->cd();
630     hTempSigma->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
631     hTempSigma->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
632     hTempSigma->DrawCopy((i == 0 ? "" : "same"));
633
634     canvPullChi2->cd();
635     hTempChi2->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
636     hTempChi2->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
637     hTempChi2->DrawCopy((i == 0 ? "" : "same"));
638   }
639
640   canvPullMean->BuildLegend();
641   canvPullSigma->BuildLegend();
642   canvPullChi2->BuildLegend();
643 */
644   // Histograms with additional correction
645   TCanvas* canvPullMeanCorr = 0x0;
646   TCanvas* canvPullSigmaCorr = 0x0;
647   TCanvas* canvPullChi2Corr = 0x0;
648   TCanvas* canvPullCorr[nThetaHistos + 1];
649   for (Int_t i = 0; i < nThetaHistos + 1; i++) 
650     canvPullCorr[i] = 0x0;
651   
652   if (hMap) {
653     // Mean, Sigma, chi^2/NDF of pull of different theta bins and all in one plot
654     canvPullMeanCorr = new TCanvas("canvPullMeanCorr", "canvPullMeanCorr", 100,10,1380,800);
655     canvPullMeanCorr->SetLogx(kTRUE);
656     canvPullMeanCorr->SetGridx(kTRUE);
657     canvPullMeanCorr->SetGridy(kTRUE);
658     canvPullSigmaCorr = new TCanvas("canvPullSigmaCorr", "canvPullSigmaCorr", 100,10,1380,800);
659     canvPullSigmaCorr->SetLogx(kTRUE);
660     canvPullSigmaCorr->SetGridx(kTRUE);
661     canvPullSigmaCorr->SetGridy(kTRUE);
662     canvPullChi2Corr = new TCanvas("canvPullChi2Corr", "canvPullChi2Corr", 100,10,1380,800);
663     canvPullChi2Corr->SetLogx(kTRUE);
664     canvPullChi2Corr->SetGridx(kTRUE);
665     canvPullChi2Corr->SetGridy(kTRUE);
666     
667     for (Int_t i = 0, j = nThetaHistos; i < nThetaHistos + 1; i++, j--)  {
668       canvPullCorr[i] = new TCanvas(Form("canvPullCorr_%d", i), "canvPullCorr", 100,10,1380,800);
669       canvPullCorr[i]->cd();
670       canvPullCorr[i]->SetLogx(kTRUE);
671       canvPullCorr[i]->SetLogz(kTRUE);
672       canvPullCorr[i]->SetGrid(kTRUE, kTRUE);
673
674       TH2D* hTemp = 0x0;
675       TString thetaString = "";
676       
677       if (i == nThetaHistos)  {
678         hTemp = hPullAdditionalCorr;
679         thetaString = "tan(#Theta) integrated";
680       }
681       else    {
682         hTemp = hPullAdditionalCorrTheta[i];
683         thetaString = Form("%.2f #leq |tan(#Theta)| < %.2f", tThetaLow[i], tThetaHigh[i]);
684       }
685
686       normaliseHisto(hTemp);
687       hTemp->FitSlicesY();
688       hTemp->GetYaxis()->SetNdivisions(12);
689       hTemp->GetXaxis()->SetMoreLogLabels(kTRUE);
690       TH1D* hTempMean = (TH1D*)gDirectory->Get(Form("%s_1", hTemp->GetName()));
691       hTempMean->SetTitle(Form("mean(pull), %s", thetaString.Data()));
692       hTempMean->GetXaxis()->SetMoreLogLabels(kTRUE);
693       hTempMean->SetLineWidth(2);
694       hTempMean->SetMarkerStyle(20);
695       TH1D* hTempSigma = (TH1D*)gDirectory->Get(Form("%s_2", hTemp->GetName()));
696       hTempSigma->SetTitle(Form("#sigma(pull), %s", thetaString.Data()));
697       hTempSigma->GetXaxis()->SetMoreLogLabels(kTRUE);
698       hTempSigma->SetLineColor(kMagenta);
699       hTempSigma->SetMarkerStyle(20);
700       hTempSigma->SetMarkerColor(kMagenta);
701       hTempSigma->SetLineWidth(2);
702       TH1D* hTempChi2 = (TH1D*)gDirectory->Get(Form("%s_chi2", hTemp->GetName()));
703       hTempChi2->SetTitle(Form("#chi^{2} / NDF (pull), %s", thetaString.Data()));
704       hTempChi2->GetXaxis()->SetMoreLogLabels(kTRUE);
705       hTempChi2->SetLineColor(kMagenta + 2);
706       hTempChi2->SetMarkerStyle(20);
707       hTempChi2->SetMarkerColor(kMagenta + 2);
708       hTempChi2->SetLineWidth(2);
709
710       hTemp->DrawCopy("colz");
711       hTempMean->DrawCopy("same");
712       hTempSigma->DrawCopy("same");
713       hTempChi2->Scale(-1./10.);
714       hTempChi2->DrawCopy("same");
715       hTempChi2->Scale(-10.);
716   
717       canvPullMeanCorr->cd();
718       hTempMean->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
719       hTempMean->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
720       hTempMean->DrawCopy((i == 0 ? "" : "same"));
721       
722       canvPullSigmaCorr->cd();
723       hTempSigma->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
724       hTempSigma->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
725       hTempSigma->DrawCopy((i == 0 ? "" : "same"));
726       
727       canvPullChi2Corr->cd();
728       hTempChi2->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
729       hTempChi2->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
730       hTempChi2->DrawCopy((i == 0 ? "" : "same"));
731     }
732     
733     canvPullMeanCorr->BuildLegend();
734     canvPullSigmaCorr->BuildLegend();
735     canvPullChi2Corr->BuildLegend();
736   }
737   
738   
739   
740   
741   
742   fSave->cd();
743   /*canvPullMean->Write();
744   canvPullSigma->Write();
745   canvPullChi2->Write();
746   
747   for (Int_t  i = 0; i < nThetaHistos + 1; i++) {
748     canvPull[i]->Write();
749   }*/
750   
751   canvPullMeanCorr->Write();
752   canvPullSigmaCorr->Write();
753   canvPullChi2Corr->Write();
754   
755   for (Int_t  i = 0; i < nThetaHistos + 1; i++) {
756     canvPullCorr[i]->Write();
757   }
758
759   TNamed* info = new TNamed(Form("Theta map: %s\n\nSigma map: %s\n\nSplines file: %s\n\nSplines name: %s", pathNameThetaMap.Data(), 
760                                  pathNameSigmaMap.Data(), pathNameSplinesFile.Data(), prSplinesName.Data()),
761                             "info");
762   info->Write();
763   fSave->Close();
764   
765   return 0;
766 }