]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TPC/macros/PIDCalib/checkPullTree.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGPP / TPC / macros / PIDCalib / checkPullTree.C
CommitLineData
88b71b9f 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//__________________________________________________________________________________________
23void 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//_________________________________________________________________________________________
40Int_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//_________________________________________________________________________________________
54Int_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//_________________________________________________________________________________________
68Int_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}