]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TPC/macros/PIDCalib/extractMultiplicityDependence.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 / extractMultiplicityDependence.C
CommitLineData
88b71b9f 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//________________________________________________________________________
22void 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//__________________________________________________________________________________________
116Int_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}