]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/QATasks/post/PostProcessQAHighPtDeDx.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] / PWGLF / QATasks / post / PostProcessQAHighPtDeDx.C
CommitLineData
21cfe1dd 1/*
2
3 Use AliRoot because of AliXRDPROOFtoolkit:
4
5
6
7 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/Base")
8 gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/macros")
9 gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/grid")
10 gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/lib")
11 gROOT->SetMacroPath(".:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/macros:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/grid:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/lib/")
12 .L my_functions.C+
13 .L my_tools.C+
14 .L PostProcessQAHighPtDeDx.C+
15
16 PlotQA("FileRoot/AnalysisResults.root")
17 MakeFitsExternalData("FileRoot/AnalysisResults.root", "HistosForBB")
18
19 MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",0)
20 MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",1)
21 MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",2)
22 MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",3)
23
24
25 PlotParametrizations("HistosForBB")
26
27
28
29
30 FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE, 0, 2, 0,1, 0, "HistosForBB/hres_0_5_02.root",27);
31 FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE, 2, 4, 0,1, 0, "HistosForBB/hres_0_5_24.root",27);
32 FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE, 4, 6, 0,1, 0, "HistosForBB/hres_0_5_46.root",27);
33 FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE, 6, 8, 0,1, 0, "HistosForBB/hres_0_5_68.root",27);
34
35
36 MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/02_dataPbPb.root",2,50, 0, "02_dataPbPb.root");
37 MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/24_dataPbPb.root",2,50, 1, "24_dataPbPb.root");
38 MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/46_dataPbPb.root",2,50, 2, "46_dataPbPb.root");
39 MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/68_dataPbPb.root",2,50, 3, "68_dataPbPb.root");
40
41
42 PlotNSigma("nsigma_results")
43
44
45
46 //add particle fractions vs p
47
48
49 ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,0, "results/eta02","fractions");
50 ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,1, "results/eta24","fractions");
51 ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,2, "results/eta46","fractions");
52 ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,3, "results/eta68","fractions");
53
54
55*/
56#include <TLegend.h>
57#include <TFile.h>
58#include <TList.h>
59#include <TTree.h>
60#include <TH2.h>
61#include <TF1.h>
62#include <TF2.h>
63#include <TH1.h>
64#include <TCanvas.h>
65#include <TProfile.h>
66#include <TClonesArray.h>
67#include <TObjString.h>
68#include <TString.h>
69#include <TROOT.h>
70#include <TStyle.h>
71#include <TLatex.h>
72#include <TGraphErrors.h>
73#include <TSpectrum.h>
74#include <TFitResultPtr.h>
75#include <TFitResult.h>
76
77#include "my_tools.C"
78#include "my_functions.C"
79
80#include <iostream>
81#include <fstream>
82#include <string>
83
84using namespace std;
85
86
87
88
89
90//const Char_t* ending[4] = {"02", "24", "46", "68"};
91
92const Char_t * legEtaCut[4]=
93 {"|#eta_{lab}|<0.2", "0.2<|#eta_{lab}|<0.4", "0.4<|#eta_{lab}|<0.6", "0.6<|#eta_{lab}|<0.8"};
94
95
96const Char_t* endingCent[1] =
97 {"MB"};
98const Char_t* CentLatex[1] =
99 {"MB"};
100
101 const Char_t *Pid[4]={"Pion", "Kaon", "Proton", "Electron"};
102//const Char_t *Pid[3]={"Pion","Kaon","Proton"};
103const Char_t *PidLatex[3]={"#pi^{+} + #pi^{-}","K^{+} + K^{-}","p + #bar{p}"};
104
105
106const Color_t PidColor[3]={2,kGreen,4};
107
108TF1* piFunc = 0;
109TF1* kFunc = 0;
110TF1* pFunc = 0;
111TF1* eFunc = 0;
112TF1* sigmaFunc = 0;
113
114const Int_t nPtBins = 68;
115
116Double_t xBins[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
117 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
118 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 ,
119 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 ,
120 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0,
121 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
122 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
123
124
125
126
127const Int_t nPtBinsV0s = 25;
128Double_t ptBinsV0s[nPtBinsV0s+1] = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 ,
129 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
130 9.0 , 12.0, 15.0, 20.0 };
131
132
133TH2D * AddTwoSameBinningTH2D(TH2D *hPos, TH2D *hNeg, const Char_t *nameHist);
134
135
136void PlotNSigma(const Char_t *path){
137
138 gStyle->SetCanvasColor(10);
139 gStyle->SetFrameFillColor(10);
140 gStyle->SetOptStat(0);
141 gStyle->SetOptTitle(0);
142 gStyle->SetOptFit(0);
143
144
145
146 TFile *fin[2]={0,0};
147
148 fin[0] = FindFileFresh(Form("%s/02_dataPbPb.root",path));
149 if(!fin[0])
150 return;
151
152 fin[1] = FindFileFresh(Form("%s/68_dataPbPb.root",path));
153 if(!fin[1])
154 return;
155
156 TH1D *hNPiP[2]={0,0};
157 TH1D *hNPiK[2]={0,0};
158 TH1D *hNPK[2]={0,0};
159 Int_t colorestmp[2]={2,4};
160 for(Int_t i=0; i<2; ++i){
161
162 hNPiP[i]=(TH1D *)fin[i]->Get("hPionProtonSeparation");
163 hNPiP[i]->SetLineColor(colorestmp[i]);
164 hNPiP[i]->SetLineWidth(2);
165 hNPiK[i]=(TH1D *)fin[i]->Get("hPionKaonSeparation");
166 hNPiK[i]->SetLineColor(colorestmp[i]);
167 hNPiK[i]->SetLineWidth(2);
168 hNPK[i]=(TH1D *)fin[i]->Get("hKaonProtonSeparation");
169 hNPK[i]->SetLineColor(colorestmp[i]);
170 hNPK[i]->SetLineWidth(2);
171
172 }
173
174
175 TH1D *hframe[3]={0,0,0};
176 for(Int_t i=0; i<3; ++i){
177
178 hframe[i] = new TH1D(Form("hframe%d",i),";#it{p} (GeV/#it{c});",10,3,15);
179 hframe[i]->GetYaxis()->SetRangeUser(-0.01,7.01);
180 hframe[i]->GetXaxis()->SetTitleSize(0.05);
181 hframe[i]->GetYaxis()->SetTitleSize(0.05);
182 }
183
184
185 TCanvas* vC1 = new TCanvas("vC1","vC1",200,10,1050,400);
186 TPad *pad[3]={0,0,0};
187 for(Int_t i=0; i<3; ++i){
188
189 pad[i]=new TPad(Form("pad%d",i),"pad1",0.01+0.33*i,0.01,0.33+0.33*i,0.99,0);
190 pad[i]->SetLeftMargin(0.12);
191 pad[i]->SetRightMargin(0.01);
192 pad[i]->SetTopMargin(0.01);
193 pad[i]->SetBottomMargin(0.1);
194 pad[i]->SetBorderSize(0);
195 pad[i]->SetBorderMode(0);
196
197 }
198
199
200
201
202 for(Int_t i =0 ; i<3; ++i){
203
204
205
206
207 vC1->cd();
208 pad[i]->Draw();
209
210
211 pad[i]->cd();
212
213
214 pad[i]->SetTickx(1);
215 pad[i]->SetTicky(1);
216
217 if(i==0)
218 hframe[i]->GetYaxis()->SetTitle("Separation #pi/p, N_{#sigma}");
219 if(i==1)
220 hframe[i]->GetYaxis()->SetTitle("Separation #pi/K, N_{#sigma}");
221 if(i==2)
222 hframe[i]->GetYaxis()->SetTitle("Separation p/K, N_{#sigma}");
223
224
225 hframe[i]->Draw();
226
227 if(i==0){
228 for(Int_t j =0 ; j<2; ++j){
229
230 hNPiP[j]->Draw("samel");
231
232 }
233 TF1 *fpip1=new TF1("fpip1","5.0+pol0",5,15);
234 fpip1->SetLineColor(1);
235 fpip1->SetLineWidth(2);
236 fpip1->SetLineStyle(2);
237 fpip1->Draw("same");
238
239 TF1 *fpip2=new TF1("fpip2","3.5+pol0",5,15);
240 fpip2->SetLineColor(1);
241 fpip2->SetLineWidth(2);
242 fpip2->SetLineStyle(4);
243 fpip2->Draw("same");
244
245 TLatex* latex0 = new TLatex();
246 latex0->SetNDC();
247 latex0->SetTextAlign(22);
248 latex0->SetTextFont(42);
249 latex0->SetTextSize(0.05);
250 latex0->DrawLatex(0.5,0.77,"pp @ 2.76 TeV, 0.6<|#eta|<0.8 (final)");
251
252 latex0->DrawLatex(0.5,0.57,"Pb-Pb, 0-5%, 0.6<|#eta|<0.8 (final)");
253
254 }
255
256 if(i==1){
257 for(Int_t j =0 ; j<2; ++j){
258 hNPiK[j]->Draw("samel");
259
260 }
261
262
263 TF1 *fpik1=new TF1("fpik1","3.5+pol0",5,15);
264 fpik1->SetLineColor(1);
265 fpik1->SetLineWidth(2);
266 fpik1->SetLineStyle(2);
267 fpik1->Draw("same");
268
269 TF1 *fpik2=new TF1("fpik2","2.2+pol0",5,15);
270 fpik2->SetLineColor(1);
271 fpik2->SetLineWidth(2);
272 fpik2->SetLineStyle(4);
273 fpik2->Draw("same");
274
275 TLegend* legend = new TLegend(0.51, 0.65, 0.85, 0.85);
276 legend->SetBorderSize(0);
277 legend->SetFillColor(0);
278 legend->SetTextFont(42);
279 legend->SetTextSize(0.05);
280 legend->SetLineColor(kWhite);
281 legend->SetLineStyle(3);
282 legend->SetShadowColor(kWhite);
283 legend->SetFillColor(kWhite);
284 legend->SetFillStyle(0);
285
286 legend->AddEntry(hNPiK[0], "|#eta|<0.2", "L");
287 legend->AddEntry(hNPiK[1], "0.6<|#eta|<0.8", "L");
288 legend->Draw();
289
290
291 }
292
293 if(i==2){
294 for(Int_t j =0 ; j<2; ++j){
295 hNPK[j]->Draw("samel");
296
297 }
298
299 TF1 *fpk1=new TF1("fpk1","1.8+pol0",5,15);
300 fpk1->SetLineColor(1);
301 fpk1->SetLineWidth(2);
302 fpk1->SetLineStyle(2);
303 fpk1->Draw("same");
304
305 TF1 *fpk2=new TF1("fpk2","1.2+pol0",5,15);
306 fpk2->SetLineColor(1);
307 fpk2->SetLineWidth(2);
308 fpk2->SetLineStyle(4);
309 fpk2->Draw("same");
310
311 }
312
313 }
314
315
316
317 vC1->SaveAs("NSigma.png");
318 vC1->SaveAs("NSigma.eps");
319
320
321
322
323
324
325
326}
327
328
329
330//____________________________________________________________________________
331void MakeNSigmaPlot(const Char_t* inFile, const Char_t* fitFileName,
332 Double_t ptStart, Double_t ptStop, Int_t i_eta,
333 const Char_t* outFileName=0)
334{
335
336 const Char_t* ending[4] = {"02", "24", "46", "68"};
337 gStyle->SetOptStat(0);
338
339 if(fitFileName) {
340
341 TFile* fitFile = FindFileFresh(fitFileName);
342 if(!fitFile)
343 return;
344 DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo");
345 fitPar->Print();
346
347 fixMIP = fitPar->MIP;
348 fixPlateau = fitPar->plateau;
349
350 Double_t dedxPar[6] = {0, 0, 0, 0, 0, 0};
351 Double_t sigmaPar[8] = {0, 0, 0, 0, 0, 0, 0, 0};
352
353 dedxPar[0] = fitPar->optionDeDx;
354 for(Int_t i = 0; i < fitPar->nDeDxPar; i++) {
355 dedxPar[i+1] = fitPar->parDeDx[i];
356 }
357
358 sigmaPar[0] = fitPar->optionSigma;
359 for(Int_t i = 0; i < fitPar->nSigmaPar; i++) {
360 sigmaPar[i+1] = fitPar->parSigma[i];
361 }
362
363 piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
364 piFunc->SetParameters(dedxPar);
365
366
367 kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
368 kFunc->SetParameters(dedxPar);
369 kFunc->SetParameter(0, kFunc->GetParameter(0)+10);
370
371 pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
372 pFunc->SetParameters(dedxPar);
373 pFunc->SetParameter(0, pFunc->GetParameter(0)+20);
374
375 eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
376 eFunc->SetParameters(dedxPar);
377 eFunc->SetParameter(0, eFunc->GetParameter(0)+30);
378
379 sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1);
380 sigmaFunc->SetParameters(sigmaPar);
381 }
382
383 kFunc->SetLineColor(kGreen+2);
384 pFunc->SetLineColor(4);
385 eFunc->SetLineColor(kMagenta);
386
387
388 /*
389 TString baseName(gSystem->BaseName(calibFileName));
390 baseName.ReplaceAll(".root", "");
391
392 TFile* calibFile = FindFileFresh(calibFileName);
393 if(!calibFile)
394 return;
395 AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, filter, usePhiCut, run, etaAbs, etaLow, etaHigh);
396 if(calib) {
397 fixMIP = calib->GetHistDeDx(kTRUE, 0)->GetMean();
398 fixPlateau = calib->GetHistDeDx(kFALSE, 0)->GetMean();
399 cout << "Plateau: " << fixPlateau << endl;
400 } else {
401 calib = (AliHighPtDeDxCalib*)calibFile->Get("v0phicut");
402 fixMIP = 50.0;
403 fixPlateau = 83.461;
404 }
405 calib->Print();
406
407 hDeDxVsP = calib->GetHistDeDxVsP(0);
408 */
409
410 TFile* dedxFile = FindFileFresh(inFile);
411 if(!dedxFile)
412 return;
413
414 TList * list = (TList *)dedxFile->Get("outputdedx");
415 TH2D *hDeDxVsPPlus = 0;
416 hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
417 hDeDxVsPPlus->Sumw2();
418
419 TH2D *hDeDxVsPMinus = 0;
420 hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
421 hDeDxVsPMinus->Sumw2();
422
423 gSystem->Exec("rm *.root");
424 TFile *fplus = new TFile("hist2Dplus.root","recreate");
425 fplus->cd();
426 hDeDxVsPPlus->SetName(Form("histAllCh%s", ending[i_eta]));
427 hDeDxVsPPlus->Write();
428 fplus->Close();
429
430 TFile *fminus = new TFile("hist2Dminus.root","recreate");
431 fminus->cd();
432 hDeDxVsPMinus->SetName(Form("histAllCh%s", ending[i_eta]));
433 hDeDxVsPMinus->Write();
434 fminus->Close();
435
436 gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
437
438 TFile *ffinal = TFile::Open("hist2Ddedx.root");
439 hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%s", ending[i_eta]));
440
441
442
443
444 //in this case, sigma is not the relative resolution
445
446 // Root is a bit stupid with finidng bins so we have to add and subtract a
447 // little to be sure we get the right bin as we typically put edges as
448 // limits
449 const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(ptStart+0.01);
450 ptStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
451 const Int_t binStop = hDeDxVsP->GetXaxis()->FindBin(ptStop-0.01);
452 ptStop = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
453 // const Int_t nBins = binStop - binStart + 1;
454
455 cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
456 << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
457
458
459 TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
460 hPionRatio->Reset();
461
462 TH1D* hPionKaonSeparation = (TH1D*)hPionRatio->Clone("hPionKaonSeparation");
463 hPionKaonSeparation->SetTitle(Form("#pi/K separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
464 TH1D* hPionProtonSeparation = (TH1D*)hPionRatio->Clone("hPionProtonSeparation");
465 hPionProtonSeparation->SetTitle(Form("#pi/p separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
466 TH1D* hKaonProtonSeparation = (TH1D*)hPionRatio->Clone("hKaonProtonSeparation");
467 hKaonProtonSeparation->SetTitle(Form("K/p separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
468
469
470
471
472 for(Int_t bin = binStart; bin <= binStop; bin++){
473
474 Double_t pt=hPionRatio->GetBinCenter(bin);
475
476 Double_t dEdx_pi = piFunc->Eval(pt);
477 Double_t Sigma_pi = sigmaFunc->Eval(piFunc->Eval(pt));
478
479 Double_t dEdx_k = kFunc->Eval(pt);
480 Double_t Sigma_k = sigmaFunc->Eval(kFunc->Eval(pt));
481
482 Double_t dEdx_p = pFunc->Eval(pt);
483 Double_t Sigma_p = sigmaFunc->Eval(pFunc->Eval(pt));
484
485 Double_t N1 = (dEdx_pi - dEdx_k)/( (Sigma_pi + Sigma_k)/2.0 );
486 Double_t N2 = (dEdx_pi - dEdx_p)/( (Sigma_pi + Sigma_p)/2.0 );
487 Double_t N3 = (dEdx_k - dEdx_p)/( (Sigma_k + Sigma_p)/2.0 );
488
489 hPionKaonSeparation->SetBinContent(bin,N1);
490 hPionProtonSeparation->SetBinContent(bin,N2);
491 hKaonProtonSeparation->SetBinContent(bin,N3);
492
493 }
494
495
496
497 if(outFileName) {
498
499 CreateDir("nsigma_results");
500 TFile* fileOut = new TFile(Form("nsigma_results/%s", outFileName), "RECREATE");
501
502 hPionKaonSeparation->Write();
503 hPionProtonSeparation->Write();
504 hKaonProtonSeparation->Write();
505
506
507
508 fileOut->Close();
509 }
510
511
512
513}
514
515
516
517//________________________________________________________________________
518void PlotParametrizations(const Char_t *path){
519
520
521 gStyle->SetCanvasColor(10);
522 gStyle->SetFrameFillColor(10);
523 gStyle->SetOptStat(0);
524 gStyle->SetOptTitle(0);
525 gStyle->SetOptFit(0);
526
527 TFile *fin[2]={0,0};
528
529
530 fin[0] = FindFileFresh(Form("%s/hres_0_5_02.root",path));
531 if(!fin[0])
532 return;
533
534 fin[1] = FindFileFresh(Form("%s/hres_0_5_68.root",path));
535 if(!fin[1])
536 return;
537
538
539 TH1D *hframeRes = new TH1D("hframeRes","; #LT dE/dx #GT; #sigma/#LT dE/dx #GT",50,40,79);
540 hframeRes->GetYaxis()->SetRangeUser(0,0.105);
541 hframeRes->GetYaxis()->SetTitleSize(0.05);
542 hframeRes->GetXaxis()->SetTitleSize(0.05);
543
544 TH1D *hframeBB = new TH1D("hframeBB","; #beta#gamma;#LT dE/dx #GT",50,2,59);
545 hframeBB->GetYaxis()->SetRangeUser(40,81);
546 hframeBB->GetYaxis()->SetTitleSize(0.05);
547 hframeBB->GetXaxis()->SetTitleSize(0.05);
548
549 TGraphErrors *gRes[2] = {0,0};
550
551 gRes[0] = (TGraphErrors *)fin[0]->Get("res_allpions");
552 gRes[1] = (TGraphErrors *)fin[1]->Get("res_allpions");
553
554
555 TGraphErrors *gBB[2] = {0,0};
556
557 gBB[0] = (TGraphErrors *)fin[0]->Get("gBB");
558 gBB[1] = (TGraphErrors *)fin[1]->Get("gBB");
559
560 Int_t colorestmp[2]={2,4};
561
562
563
564 TCanvas* vC1 = new TCanvas("vC1","vC1",200,10,700,400);
565 TPad * pad1=new TPad("pad1","pad1",0.01,0.01,0.5,0.99,0);
566 TPad * pad2=new TPad("pad2","pad2",0.5,0.01,0.99,0.99,0);
567
568 pad1->SetLeftMargin(0.12);
569 pad1->SetRightMargin(0.01);
570 pad1->SetTopMargin(0.01);
571 pad1->SetBottomMargin(0.1);
572 pad1->SetBorderSize(0);
573 pad1->SetBorderMode(0);
574
575 pad2->SetLeftMargin(0.12);
576 pad2->SetRightMargin(0.01);
577 pad2->SetBottomMargin(0.1);
578 pad2->SetTopMargin(0);
579 pad2->SetBorderSize(0);
580 vC1->cd();
581 pad1->Draw();
582
583
584 pad1->cd();
585
586
587 pad1->SetTickx(1);
588 pad1->SetTicky(1);
589
590
591
592
593
594
595 hframeRes->Draw();
596
597 for(Int_t i =0 ; i<2; ++i){
598 gRes[i]->SetMarkerStyle(20);
599 gRes[i]->SetLineWidth(2);
600 gRes[i]->SetLineColor(colorestmp[i]);
601 gRes[i]->SetMarkerColor(colorestmp[i]);
602 gRes[i]->Draw("samep");
603 }
604
605
606 TF1 *f0005 = new TF1("f0005","pol2",50,80);
607 // f0005->SetParameter(0,1.20000000000000000e+01);
608 f0005->SetParameter(0,6.08803411659195118e-02);
609 f0005->SetParameter(1,2.34760597152924448e-04);
610 f0005->SetParameter(2,-2.81841683363273653e-06);
611
612 f0005->SetLineColor(1);
613 f0005->SetLineWidth(2);
614 f0005->SetLineStyle(2);
615 f0005->Draw("samel");
616
617
618 TF1 *fpp = new TF1("fpp","pol2",50,80);
619 fpp->SetParameter(0,1.06773399595441187e-01);
620 fpp->SetParameter(1,-1.55087701882609280e-03);
621 fpp->SetParameter(2,1.01579672586848698e-05);
622
623 fpp->SetLineColor(1);
624 fpp->SetLineWidth(2);
625 fpp->SetLineStyle(3);
626 fpp->Draw("samel");
627
628
629 TLegend* legendRes = new TLegend(0.21, 0.15, 0.65, 0.35);
630 legendRes->SetBorderSize(0);
631 legendRes->SetFillColor(0);
632 legendRes->SetTextFont(42);
633 legendRes->SetTextSize(0.05);
634 legendRes->SetLineColor(kWhite);
635 legendRes->SetLineStyle(3);
636 legendRes->SetShadowColor(kWhite);
637 legendRes->SetFillColor(kWhite);
638 legendRes->SetFillStyle(0);
639 legendRes->SetHeader("Final, 0.6<|#eta|<0.8");
640 legendRes->AddEntry(f0005, "Pb-Pb 0-5%", "L");
641 legendRes->AddEntry(fpp, "pp @ 2.76 TeV", "L");
642 legendRes->Draw();
643
644 vC1->cd();
645 pad2->Draw();
646
647
648 pad2->cd();
649
650
651 pad2->SetTickx(1);
652 pad2->SetTicky(1);
653
654
655
656
657
658
659 hframeBB->Draw();
660
661 for(Int_t i =0 ; i<2; ++i){
662 gBB[i]->SetMarkerStyle(20);
663 gBB[i]->SetLineWidth(2);
664 gBB[i]->SetLineColor(colorestmp[i]);
665 gBB[i]->SetMarkerColor(colorestmp[i]);
666 gBB[i]->Draw("samep");
667 }
668
669 TLegend* legend = new TLegend(0.51, 0.8, 0.85, 0.95);
670 legend->SetBorderSize(0);
671 legend->SetFillColor(0);
672 legend->SetTextFont(42);
673 legend->SetTextSize(0.05);
674 legend->SetLineColor(kWhite);
675 legend->SetLineStyle(3);
676 legend->SetShadowColor(kWhite);
677 legend->SetFillColor(kWhite);
678 legend->SetFillStyle(0);
679
680 legend->AddEntry(gBB[0], "|#eta|<0.2", "P");
681 legend->AddEntry(gBB[1], "0.6<|#eta|<0.8", "P");
682 legend->Draw();
683
684
685
686
687 vC1->SaveAs("Parametrizations.png");
688 vC1->SaveAs("Parametrizations.eps");
689
690
691
692}
693//____________________________________________________________________________
694void ExtractUncPartFractvsP(const Char_t * inFile, Double_t ptStart, Double_t ptStop,
695 Int_t i_cent=1,
696 Int_t i_eta=1,
697 const Char_t* dirName="debugfits", const Char_t* outFileName=0)
698
699
700{
701
702
703 const Char_t* ending[4] = {"02", "24", "46", "68"};
704 const Double_t etaHigh[4]={2,4,6,8};
705 gStyle->SetOptStat(0);
706
707
708
709 TFile* fitFile = FindFileFresh(Form("fitparameters/MB/%s_dataPbPb.root",ending[i_eta]));
710 if(!fitFile)
711 return;
712
713
714 DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo");
715 fitPar->Print();
716
717 fixMIP = fitPar->MIP;
718
719
720
721 //return;
722
723 Double_t dedxPar[6] = {0, 0, 0, 0, 0, 0};
724 Double_t sigmaPar[8] = {0, 0, 0, 0, 0, 0, 0, 0};
725
726 dedxPar[0] = fitPar->optionDeDx;
727 for(Int_t i = 0; i < fitPar->nDeDxPar; i++) {
728 dedxPar[i+1] = fitPar->parDeDx[i];
729 cout<<"idedx_par"<<i+1<<"="<<dedxPar[i+1]<<endl;
730
731 }
732 fixPlateau = dedxPar[5];
733 cout<<"fixPlateau="<<fixPlateau<<" fixMIP="<<fixMIP<<endl;
734 //return;
735
736 sigmaPar[0] = fitPar->optionSigma;
737 for(Int_t i = 0; i < fitPar->nSigmaPar; i++) {
738 sigmaPar[i+1] = fitPar->parSigma[i];
739 }
740
741 piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
742 piFunc->SetParameters(dedxPar);
743
744
745 kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
746 kFunc->SetParameters(dedxPar);
747 kFunc->SetParameter(0, kFunc->GetParameter(0)+10);
748
749 pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
750 pFunc->SetParameters(dedxPar);
751 pFunc->SetParameter(0, pFunc->GetParameter(0)+20);
752
753 eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
754 eFunc->SetParameters(dedxPar);
755 eFunc->SetParameter(0, eFunc->GetParameter(0)+30);
756
757 sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1);
758 sigmaFunc->SetParameters(sigmaPar);
759
760
761 kFunc->SetLineColor(kGreen+2);
762 pFunc->SetLineColor(4);
763 eFunc->SetLineColor(kMagenta);
764
765 /*
766 TFile* calibFile = FindFileFresh(inFile);
767 if(!calibFile)
768 return;
769 hDeDxVsP = (TH2D *)calibFile->Get(Form("hDeDxVsP%s_%s", ending[i_eta], endingCent[i_cent]));
770 hDeDxVsP->Sumw2();
771 */
772
773 TFile* dedxFile = FindFileFresh(inFile);
774 if(!dedxFile)
775 return;
776
777 TList * list = (TList *)dedxFile->Get("outputdedx");
778 TH2D *hDeDxVsPPlus = 0;
779 hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
780 hDeDxVsPPlus->Sumw2();
781
782 TH2D *hDeDxVsPMinus = 0;
783 hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
784 hDeDxVsPMinus->Sumw2();
785
786 gSystem->Exec("rm *.root");
787 TFile *fplus = new TFile("hist2Dplus.root","recreate");
788 fplus->cd();
789 hDeDxVsPPlus->SetName(Form("histAllCh%s", ending[i_eta]));
790 hDeDxVsPPlus->Write();
791 fplus->Close();
792
793 TFile *fminus = new TFile("hist2Dminus.root","recreate");
794 fminus->cd();
795 hDeDxVsPMinus->SetName(Form("histAllCh%s", ending[i_eta]));
796 hDeDxVsPMinus->Write();
797 fminus->Close();
798
799 gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
800
801 TFile *ffinal = TFile::Open("hist2Ddedx.root");
802 hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%s", ending[i_eta]));
803
804
805
806
807
808
809
810 CreateDir(Form("fit_results/%s_%s",dirName,endingCent[i_cent]));
811
812 TCanvas* cDeDxVsPLogX = new TCanvas("cDeDxVsPLogX", "; #it{p} (GeV/#it{c}) ;d#it{E}/d#it{x}", 500, 400);
813 cDeDxVsPLogX->Clear();
814 cDeDxVsPLogX->cd();
815 cDeDxVsPLogX->SetLogz();
816
817 TH2D *hDeDxVsP2=(TH2D *)hDeDxVsP->Clone("h2DP");
818 hDeDxVsP2->SetTitle("; #it{p} (GeV/#it{c}) ;d#it{E}/d#it{x}");
819 hDeDxVsP2->GetXaxis()->SetRangeUser(2,19.5);
820 hDeDxVsP2->Draw("COLZ");
821
822 piFunc->SetLineColor(1);
823 piFunc->Draw("samel");
824
825 eFunc->SetLineColor(1);
826 eFunc->Draw("samel");
827
828 kFunc->SetLineColor(1);
829 kFunc->Draw("samel");
830
831 pFunc->SetLineColor(1);
832 pFunc->Draw("samel");
833
834 TLatex* latex0 = new TLatex();
835 latex0->SetNDC();
836 latex0->SetTextAlign(22);
837 latex0->SetTextSize(0.05);
838 latex0->DrawLatex(0.7,0.94+0.035,Form("p-Pb, %s",legEtaCut[i_eta]));
839 latex0->DrawLatex(0.71,0.89+0.035,Form("#sqrt{s_{NN}}=5.02 TeV, %s",CentLatex[i_cent]));
840
841
842
843 cDeDxVsPLogX->SaveAs(Form("fit_results/%s_%s/%s.png",dirName,endingCent[i_cent],cDeDxVsPLogX->GetName()));
844 cDeDxVsPLogX->SaveAs(Form("fit_results/%s_%s/%s.eps",dirName,endingCent[i_cent],cDeDxVsPLogX->GetName()));
845 //in this case, sigma is not the relative resolution
846
847 // Root is a bit stupid with finidng bins so we have to add and subtract a
848 // little to be sure we get the right bin as we typically put edges as
849 // limits
850 const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(ptStart+0.01);
851 ptStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
852 const Int_t binStop = hDeDxVsP->GetXaxis()->FindBin(ptStop-0.01);
853 ptStop = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
854 // const Int_t nBins = binStop - binStart + 1;
855
856 cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
857 << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
858
859
860 //cross check
861 TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
862 cFits->Clear();
863 cFits->Divide(7, 4);
864
865
866 TF1* pion = new TF1("pion", "gausn", 40, 100);
867
868 TF1* kaon = new TF1("kaon", "gausn", 40, 100);
869
870 TF1* proton = new TF1("proton", "gausn", 40, 100);
871 //proton->SetLineWidth(2);
872 //proton->SetLineColor(kBlue);
873 TF1* electron = new TF1("electron", "gausn", 40, 100);
874 //electron->SetLineWidth(2);
875 //electron->SetLineColor(kMagenta);
876 // TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)+gausn(9)", -30, 30);
877 TF1* total = 0;
878 total = new TF1("total", "([0]+[12])*gausn(1)+[4]*gausn(5)+[8]*gausn(9)+[12]*gausn(13)", 40, 100);
879
880 total->SetLineColor(kBlack);
881 total->SetLineWidth(2);
882 total->SetLineStyle(2);
883
884 TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);
885 legend->SetBorderSize(0);
886 legend->SetFillColor(0);
887 legend->AddEntry(total, "4-Gauss fit", "L");
888 legend->AddEntry(pion, "#pi", "L");
889 legend->AddEntry(kaon, "K", "L");
890 legend->AddEntry(proton, "p", "L");
891 legend->AddEntry(electron, "e", "L");
892
893 TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
894 cSingleFit->SetLeftMargin(0.11);
895 cSingleFit->SetRightMargin(0.03);
896 cSingleFit->SetTopMargin(0.01);
897 cSingleFit->SetBottomMargin(0.1);
898 cSingleFit->SetBorderSize(0);
899 cSingleFit->SetBorderMode(0);
900
901
902 TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
903 hPionRatio->SetTitle("particle fractions; p [GeV/c]; particle fraction");
904 hPionRatio->Reset();
905 TH1D* hKaonRatio = (TH1D*)hPionRatio->Clone("hKaonRatio");
906 TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
907 TH1D* hElectronRatio = (TH1D*)hPionRatio->Clone("hElectronRatio");
908
909
910
911 TH1D* hPionYield =(TH1D*)hDeDxVsP->ProjectionX("hPionYield", 1, 1);
912 hPionYield->SetTitle("particle fractions; p [GeV/c]; particle fraction");
913 hPionYield->Reset();
914 TH1D* hKaonYield = (TH1D*)hPionYield->Clone("hKaonYield");
915 TH1D* hProtonYield = (TH1D*)hPionYield->Clone("hProtonYield");
916 TH1D* hElectronYield = (TH1D*)hPionYield->Clone("hElectronYield");
917
918
919 TF1 *fmip=new TF1("fmip","pol0",0,1);
920 fmip->SetParameter(0,fixMIP);
921
922
923
924 TF1* fElectronFraction = 0;
925 if(etaHigh[i_eta]==8||etaHigh[i_eta]==-6)
926 fElectronFraction = new TF1("fElectronFraction", "[0]+(x<9.0)*[1]*(x-9.0)", 0, ptStop);
927 if(etaHigh[i_eta]==6||etaHigh[i_eta]==-4)
928 fElectronFraction = new TF1("fElectronFraction", "[0]+(x<8.0)*[1]*(x-8.0)", 0, ptStop);
929 if(etaHigh[i_eta]==4||etaHigh[i_eta]==-2)
930 fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.5)*[1]*(x-7.5)", 0, ptStop);
931 if(etaHigh[i_eta]==2||etaHigh[i_eta]==0)
932 fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.0)*[1]*(x-7.0)", 0, ptStop);
933
934 //TF1* fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.0)*[1]*(x-7.0)", 0, ptStop);
935 fElectronFraction->SetParameters(0.01, 0.0);
936
937
938 TH1D *dEdxpoints[binStop-binStart];
939 TF1 *fPion[binStop-binStart];
940 TF1 *fKaon[binStop-binStart];
941 TF1 *fProton[binStop-binStart];
942 TF1 *fElectron[binStop-binStart];
943 TF1 *fTotal[binStop-binStart];
944
945
946 for(Int_t bin = binStart; bin <= binStop; bin++){
947
948 Double_t pt=hPionRatio->GetBinCenter(bin);
949
950 cout << "Making projection for bin: " << bin << endl;
951
952 const Int_t j = bin-binStart;
953
954 TH1D* hDeltaPiVsPtProj =(TH1D*)hDeDxVsP->ProjectionY(Form("hDeltaPiVsPtProj%d", bin), bin, bin);
955 // hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(-25, 20);
956 hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(40, 90);
957 hDeltaPiVsPtProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
958 hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
959 hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
960
961 dEdxpoints[j]=0;
962 dEdxpoints[j]=hDeltaPiVsPtProj;
963 fPion[j]=0;
964 fPion[j]=new TF1(Form("fPiongauss_%d",bin), "gausn", 40, 100);
965
966 fKaon[j]=0;
967 fKaon[j]=new TF1(Form("fKaongauss_%d",bin), "gausn", 40, 100);
968
969 fProton[j]=0;
970 fProton[j]=new TF1(Form("fProtongauss_%d",bin), "gausn", 40, 100);
971
972 fElectron[j]=0;
973 fElectron[j]=new TF1(Form("fElectrongauss_%d",bin), "gausn", 40, 100);
974
975 fTotal[j]=0;
976 fTotal[j]=new TF1(Form("fTotalgauss_%d",bin), "([0]+[12])*gausn(1)+[4]*gausn(5)+[8]*gausn(9)+[12]*gausn(13)", 40, 100);
977
978
979 Double_t all = hDeltaPiVsPtProj->GetEntries();
980
981
982 const Int_t nPar = 16;
983 Double_t gausParams[nPar] = {
984 0.59,
985 all,
986 piFunc->Eval(pt),
987 //fpi->Eval(pt),
988 sigmaFunc->Eval(piFunc->Eval(pt)),
989
990 0.3,
991 all,
992 kFunc->Eval(pt),
993 //fp->Eval(pt),
994 sigmaFunc->Eval(kFunc->Eval(pt)),
995
996 0.1,
997 all,
998 pFunc->Eval(pt),
999 //fp->Eval(pt),
1000 sigmaFunc->Eval(pFunc->Eval(pt)),
1001
1002 0.01,
1003 all,
1004 eFunc->Eval(pt),
1005 //fpi->Eval(pt),
1006 sigmaFunc->Eval(eFunc->Eval(pt)),
1007 };
1008
1009
1010 for(Int_t ipar=0;ipar<nPar;++ipar)
1011 cout<<gausParams[ipar]<<endl;
1012
1013
1014
1015 cFits->cd();
1016 cFits->cd(j + 1);
1017
1018 total->SetParameters(gausParams);
1019
1020
1021 for(Int_t i = 0; i < nPar; i++) {
1022 if(i!=0 && i!=4 && i!=8 && i!=12)
1023 total->FixParameter(i, gausParams[i]);
1024 else
1025 continue;
1026 }
1027
1028 total->SetParLimits(7, gausParams[7]-0.05*gausParams[7],gausParams[7]+0.05*gausParams[7]);
1029 total->SetParLimits(8, 0.005,0.4);
1030 total->SetParLimits(4, 0.005,0.6);
1031
1032 if(bin==48) {
1033 if(etaHigh[i_eta]==8||etaHigh[i_eta]==-6)
1034 hElectronRatio->Fit(fElectronFraction, "N", "", 4.0, 10.0);
1035 if(etaHigh[i_eta]==6||etaHigh[i_eta]==-4)
1036 hElectronRatio->Fit(fElectronFraction, "N", "", 3.66, 10.0);
1037 if(etaHigh[i_eta]==4||etaHigh[i_eta]==-2)
1038 hElectronRatio->Fit(fElectronFraction, "N", "", 3.33, 10.0);
1039 if(etaHigh[i_eta]==2||etaHigh[i_eta]==0)
1040 hElectronRatio->Fit(fElectronFraction, "N", "", 3.0, 10.0);
1041 fElectronFraction->SetRange(0, ptStop);
1042 }
1043
1044 if(bin>=48) {
1045 total->FixParameter(12, fElectronFraction->Eval(hDeDxVsP->GetXaxis()->GetBinCenter(bin)));
1046 }
1047
1048
1049 hDeltaPiVsPtProj->Fit(total, "0L");
1050
1051
1052 hDeltaPiVsPtProj->Draw();
1053 total->DrawCopy("same");
1054
1055 Double_t parametersOut[nPar];
1056 total->GetParameters(parametersOut);
1057 const Double_t* parameterErrorsOut = total->GetParErrors();
1058
1059 for(Int_t i = 0; i < nPar; i++)
1060 cout << parametersOut[i] << ", ";
1061 cout << endl;
1062 fTotal[j]->SetParameters(parametersOut);
1063
1064
1065
1066 pion->SetParameters(&parametersOut[1]);
1067 pion->SetParameter(0,all*(parametersOut[0]+parametersOut[12]));
1068
1069 fPion[j]->SetParameters(&parametersOut[1]);;
1070 fPion[j]->SetParameter(0,all*(parametersOut[0]+parametersOut[12]));
1071
1072
1073 hPionRatio->SetBinContent(bin, parametersOut[0]);
1074 hPionRatio->SetBinError(bin, parameterErrorsOut[0]);
1075 hPionYield->SetBinContent(bin, parametersOut[0]*all);
1076 hPionYield->SetBinError(bin, parameterErrorsOut[0]*all);
1077
1078 kaon->SetParameters(&parametersOut[5]);
1079 kaon->SetParameter(0,all*parametersOut[4]);
1080
1081 fKaon[j]->SetParameters(&parametersOut[5]);
1082 fKaon[j]->SetParameter(0,all*parametersOut[4]);
1083
1084
1085 hKaonRatio->SetBinContent(bin, parametersOut[4]);
1086 hKaonRatio->SetBinError(bin, parameterErrorsOut[4]);
1087 hKaonYield->SetBinContent(bin, parametersOut[4]*all);
1088 hKaonYield->SetBinError(bin, parameterErrorsOut[4]*all);
1089
1090 proton->SetParameters(&parametersOut[9]);
1091 proton->SetParameter(0,all*parametersOut[8]);
1092 //proton->DrawCopy("same");
1093
1094 fProton[j]->SetParameters(&parametersOut[9]);
1095 fProton[j]->SetParameter(0,all*parametersOut[8]);
1096
1097
1098 hProtonRatio->SetBinContent(bin, parametersOut[8]);
1099 hProtonRatio->SetBinError(bin, parameterErrorsOut[8]);
1100 hProtonYield->SetBinContent(bin, parametersOut[8]*all);
1101 hProtonYield->SetBinError(bin, parameterErrorsOut[8]*all);
1102
1103 electron->SetParameters(&parametersOut[13]);
1104 electron->SetParameter(0,all*parametersOut[12]);
1105
1106 fElectron[j]->SetParameters(&parametersOut[13]);
1107 fElectron[j]->SetParameter(0,all*parametersOut[12]);
1108
1109
1110 //electron->DrawCopy("same");
1111 hElectronRatio->SetBinContent(bin, parametersOut[12]);
1112 hElectronRatio->SetBinError(bin, parameterErrorsOut[12]);
1113 hElectronYield->SetBinContent(bin, parametersOut[12]*all);
1114 hElectronYield->SetBinError(bin, parameterErrorsOut[12]*all);
1115
1116 //DrawALICELogo(kFALSE, 0.71, 0.76, 0.82, 0.98);
1117 cSingleFit->cd();
1118 cSingleFit->Clear();
1119 cSingleFit->SetTickx(1);
1120 cSingleFit->SetTicky(1);
1121 cSingleFit->SetLogy(0);
1122
1123 // cSingleFit->SetLogy();
1124 hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(40,90);
1125 hDeltaPiVsPtProj->SetMarkerStyle(25);
1126 hDeltaPiVsPtProj->SetMarkerColor(1);
1127 hDeltaPiVsPtProj->SetLineColor(1);
1128 hDeltaPiVsPtProj->SetLineWidth(2);
1129 hDeltaPiVsPtProj->GetYaxis()->SetLabelFont(42);
1130 hDeltaPiVsPtProj->GetXaxis()->SetLabelFont(42);
1131 hDeltaPiVsPtProj->GetZaxis()->SetLabelFont(42);
1132 hDeltaPiVsPtProj->GetYaxis()->SetTitleFont(42);
1133 hDeltaPiVsPtProj->GetXaxis()->SetTitleFont(42);
1134 hDeltaPiVsPtProj->GetZaxis()->SetTitleFont(42);
1135 hDeltaPiVsPtProj->GetYaxis()->SetTitle("Entries");
1136 hDeltaPiVsPtProj->GetXaxis()->SetTitle("d#it{E}/d#it{x} in TPC (arb. units)");
1137 hDeltaPiVsPtProj->GetXaxis()->SetTitleSize(0.05);
1138 hDeltaPiVsPtProj->GetXaxis()->SetTitleOffset(0.9);
1139 hDeltaPiVsPtProj->GetYaxis()->SetTitleSize(0.05);
1140 hDeltaPiVsPtProj->GetYaxis()->SetTitleOffset(0.9);
1141 hDeltaPiVsPtProj->Draw();
1142
1143 total->SetLineColor(kGray+1);
1144 total->SetLineWidth(3);
1145 total->SetLineStyle(1);
1146 total->DrawCopy("same");
1147
1148 pion->GetXaxis()->SetRangeUser(40,90);
1149 pion->SetLineColor(2);
1150 pion->SetLineWidth(2);
1151 //pion->SetLineStyle(2);
1152 pion->DrawCopy("same");
1153
1154 kaon->SetLineColor(kGreen+2);
1155 kaon->SetLineWidth(2);
1156 //kaon->SetLineStyle(2);
1157 kaon->DrawCopy("same");
1158
1159
1160 proton->SetLineColor(4);
1161 proton->SetLineWidth(2);
1162 //proton->SetLineStyle(2);
1163 proton->DrawCopy("same");
1164
1165
1166 TLatex* latex = new TLatex();
1167 // latex->SetTextFont(132);
1168 latex->SetNDC();
1169 latex->SetTextAlign(22);
1170
1171 latex->SetTextSize(0.05);
1172 latex->SetTextFont(42);
1173 latex->SetTextColor(kRed+2);
1174 //latex->DrawLatex(0.8,0.3,"pp");
1175 //latex->DrawLatex(0.81,0.25,"#sqrt{s}=2.76 TeV");
1176 latex->DrawLatex(0.8,0.25+0.035,Form("p-Pb, %s",CentLatex[i_cent]));
1177 latex->DrawLatex(0.81,0.2+0.035,"#sqrt{s_{NN}}=5.02 TeV");
1178
1179 latex->SetTextColor(1);
1180 latex->SetTextSize(0.04);
1181 //latex->DrawLatex(0.81,0.2,"25/07/2012");
1182 //latex->DrawLatex(0.80,0.3+0.036,"31/10/2012");
1183
1184 latex->SetTextSize(0.05);
1185 latex->DrawLatex(0.8,0.71,legEtaCut[i_eta]);
1186 latex->DrawLatex(0.8,0.64,Form("%.1f<#it{p}<%.1f GeV/#it{c}", hDeDxVsP->GetXaxis()->GetBinLowEdge(bin), hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
1187
1188
1189
1190 //cSingleFit->SaveAs(Form("%s/ptspectrum_bin%d.gif", dirName, bin));
1191 cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.png", dirName, endingCent[i_cent],bin));
1192 cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.pdf", dirName, endingCent[i_cent],bin));
1193 cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.eps", dirName, endingCent[i_cent],bin));
1194 // legend->Draw();
1195
1196
1197
1198
1199 }
1200
1201
1202 TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p; pt", 600, 400);
1203
1204 cRatio->Clear();
1205 cRatio->SetGridy();
1206 hElectronRatio->GetYaxis()->SetRangeUser(-0.05, 0.1);
1207 hElectronRatio->DrawCopy("P E");
1208 fElectronFraction->DrawCopy("SAME");
1209
1210 gROOT->ProcessLine(".x drawText.C");
1211 cRatio->SaveAs(Form("fit_results/%s_%s/electron_ratio.gif", dirName,endingCent[i_cent]));
1212 cRatio->SaveAs(Form("fit_results/%s_%s/electron_ratio.pdf", dirName,endingCent[i_cent]));
1213
1214 cRatio->Clear();
1215 cRatio->SetGridy(kFALSE);
1216 hPionRatio->SetMarkerStyle(20);
1217 hPionRatio->SetMarkerColor(2);
1218 hPionRatio->SetLineColor(2);
1219 hPionRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1220 hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
1221 hPionRatio->DrawCopy("P E");
1222 hPionYield->SetMarkerStyle(20);
1223 hPionYield->SetMarkerColor(2);
1224 hPionYield->SetLineColor(2);
1225 hPionYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1226
1227 hKaonRatio->SetMarkerStyle(20);
1228 hKaonRatio->SetMarkerColor(kGreen);
1229 hKaonRatio->SetLineColor(kGreen);
1230 hKaonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1231 hKaonYield->SetMarkerStyle(20);
1232 hKaonYield->SetMarkerColor(kGreen);
1233 hKaonYield->SetLineColor(kGreen);
1234 hKaonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1235 hKaonRatio->DrawCopy("SAME P E");
1236
1237 hProtonRatio->SetMarkerStyle(20);
1238 hProtonRatio->SetMarkerColor(4);
1239 hProtonRatio->SetLineColor(4);
1240 hProtonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1241 hProtonYield->SetMarkerStyle(20);
1242 hProtonYield->SetMarkerColor(4);
1243 hProtonYield->SetLineColor(4);
1244 hProtonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1245 hProtonRatio->DrawCopy("SAME P E");
1246
1247 hElectronRatio->SetMarkerStyle(20);
1248 hElectronRatio->SetMarkerColor(6);
1249 hElectronRatio->SetLineColor(6);
1250 hElectronRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1251 hElectronRatio->DrawCopy("SAME P E");
1252 hElectronYield->SetMarkerStyle(20);
1253 hElectronYield->SetMarkerColor(6);
1254 hElectronYield->SetLineColor(6);
1255 hElectronYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1256
1257 TH1D *hallRec=0;
1258
1259 hallRec=(TH1D *)hPionRatio->Clone("hRecAll");
1260 hallRec->Add(hKaonRatio);
1261 hallRec->Add(hProtonRatio);
1262 hallRec->SetLineColor(kYellow-3);
1263 hallRec->SetMarkerColor(kYellow-3);
1264 hallRec->DrawCopy("SAME P E");
1265
1266
1267 gROOT->ProcessLine(".x drawText.C");
1268 cRatio->SaveAs(Form("fit_results/%s_%s/particle_ratios.gif", dirName,endingCent[i_cent]));
1269 cRatio->SaveAs(Form("fit_results/%s_%s/particle_ratios.pdf", dirName,endingCent[i_cent]));
1270
1271
1272
1273 if(outFileName) {
1274
1275 CreateDir(Form("fit_results/fit_yields_results_%s",endingCent[i_cent]));
1276 TFile* fileOut = new TFile(Form("fit_results/fit_yields_results_%s/%s%s.root", endingCent[i_cent], outFileName, ending[i_eta]), "RECREATE");
1277
1278
1279 hPionRatio->Write();
1280 hKaonRatio->Write();
1281 hProtonRatio->Write();
1282 hElectronRatio->Write();
1283
1284 hPionYield->Write();
1285 hKaonYield->Write();
1286 hProtonYield->Write();
1287 hElectronYield->Write();
1288
1289
1290
1291 fElectronFraction->Write();
1292
1293 for(Int_t bin = binStart; bin <= binStop; bin++){
1294
1295 Int_t j = bin-binStart;
1296 dEdxpoints[j]->Write();
1297 fPion[j]->Write();
1298 fKaon[j]->Write();
1299 fProton[j]->Write();
1300 fElectron[j]->Write();
1301 fTotal[j]->Write();
1302
1303
1304 }
1305 fmip->Write();
1306 fileOut->Close();
1307 }
1308
1309
1310
1311
1312}
1313
1314
1315
1316
1317
1318
1319//____________________________________________________________________________
1320void FitDeDxVsP(const Char_t * dedxFileName,
1321 Double_t pStart, Double_t pStop, Int_t i_cent,
1322 Int_t optionDeDx, Int_t optionSigma,
1323 Bool_t performFit = kFALSE,
1324 Int_t etaLow=0, Int_t etaHigh=8,
1325 Bool_t fixAllParBB=kFALSE,
1326 Bool_t fixAllParSigma=kFALSE,
1327 Bool_t fixKtoZero=kFALSE,
1328 const Char_t* initialFitFileName=0,
1329 Int_t fixParametersDedx=-1,
1330 Int_t fixParametersSigma=-1)
1331{
1332 //gStyle->SetOptStat(0);
1333
1334
1335
1336 TFile* dedxFile = FindFileFresh(dedxFileName);
1337 if(!dedxFile)
1338 return;
1339
1340 TList * list = (TList *)dedxFile->Get("outputdedx");
1341 TH2D *hDeDxVsPPlus = 0;
1342 hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%d%d", etaLow, etaHigh));
1343 hDeDxVsPPlus->Sumw2();
1344
1345 TH2D *hDeDxVsPMinus = 0;
1346 hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%d%d", etaHigh, etaLow));
1347 hDeDxVsPMinus->Sumw2();
1348
1349 gSystem->Exec("rm *.root");
1350 TFile *fplus = new TFile("hist2Dplus.root","recreate");
1351 fplus->cd();
1352 hDeDxVsPPlus->SetName(Form("histAllCh%d%d", etaLow, etaHigh));
1353 hDeDxVsPPlus->Write();
1354 fplus->Close();
1355
1356 TFile *fminus = new TFile("hist2Dminus.root","recreate");
1357 fminus->cd();
1358 hDeDxVsPMinus->SetName(Form("histAllCh%d%d", etaLow, etaHigh));
1359 hDeDxVsPMinus->Write();
1360 fminus->Close();
1361
1362 gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
1363
1364 TFile *ffinal = TFile::Open("hist2Ddedx.root");
1365 hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%d%d", etaLow, etaHigh));
1366
1367
1368
1369 CreateDir("old");
1370 gSystem->Exec(Form("mv results/calibdedx_%d%d/* old/",etaLow, etaHigh));
1371
1372
1373 TCanvas* cDeDxVsP = new TCanvas("cDeDxVsP", "dE/dx vs p", 400, 300);
1374 cDeDxVsP->Clear();
1375 cDeDxVsP->cd();
1376 cDeDxVsP->SetLogz();
1377 hDeDxVsP->SetTitle(0);
1378 hDeDxVsP->Draw("COLZ");
1379
1380 TCanvas* cDeDxVsPLogX = new TCanvas("cDeDxVsPLogX", "dE/dx vs p", 400, 300);
1381 cDeDxVsPLogX->Clear();
1382 cDeDxVsPLogX->cd();
1383 cDeDxVsPLogX->SetLogz();
1384 cDeDxVsPLogX->SetLogx();
1385 hDeDxVsP->Draw("COLZ");
1386
1387 // Root is a bit stupid with finidng bins so we have to add and subtract a
1388 // little to be sure we get the right bin as we typically put edges as
1389 // limits
1390 const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(pStart+0.01);
1391 pStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
1392 const Int_t binStop = hDeDxVsP->GetXaxis()->FindBin(pStop-0.01);
1393 pStop = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
1394 const Int_t nBins = binStop - binStart + 1;
1395
1396 cout << "Doing 2d fit from pTlow = " << pStart << " (bin: " << binStart
1397 << ") to pThigh = " << pStop << " (bin: " << binStop << ")" << endl;
1398
1399 // Double_t binSize = (histdEdxvsP->GetXaxis()->GetXmax() - histdEdxvsP->GetXaxis()->GetXmin())/ histdEdxvsP->GetXaxis()->GetNbins();
1400
1401 Double_t parDeDx[6] = {0, 0, 0, 0, 0, 0};
1402 Double_t parSigma[8] = {0, 0, 0, 0, 0, 0, 0, 0};
1403
1404 const Int_t nLocalParams = 3; // pi, K, p yields
1405 Int_t nDeDxPar = 0;
1406 Int_t nSigmaPar = 0;
1407
1408 switch(optionDeDx) {
1409
1410 case 1:
1411 nDeDxPar = 2;
1412 parDeDx[0] = 39.7;
1413 parDeDx[1] = 6.3;
1414 break;
1415 case 2:
1416 nDeDxPar = 1;
1417 parDeDx[0] = 7.3;
1418 break;
1419 case 3:
1420 nDeDxPar = 2;
1421 parDeDx[0] = 6.85097;
1422 parDeDx[1] = 29.5933;
1423 break;
1424 case 4:
1425 nDeDxPar = 3;
1426 parDeDx[0] = 31.2435;
1427 parDeDx[1] = 9.73037;
1428 parDeDx[2] = 1.95832;
1429 break;
1430 case 5:
1431 nDeDxPar = 4;
1432 parDeDx[0] = 31.384;
1433 parDeDx[1] = 9.5253;
1434 parDeDx[2] = 7.3542;
1435 parDeDx[3] = 1.4;
1436 break;
1437 case 6:
1438 nDeDxPar = 5;
1439 parDeDx[0] = 31.384;
1440 parDeDx[1] = 9.5253;
1441 parDeDx[2] = 7.3542;
1442 parDeDx[3] = 1.5;
1443 parDeDx[4] = 81;
1444 break;
1445 case 7:
1446 nDeDxPar = 3;
1447 parDeDx[0] = 31.3;
1448 parDeDx[1] = 9.5;
1449 parDeDx[2] = 1.4;
1450 break;
1451 default:
1452
1453 cout << "optionDeDx does not support option: " << optionDeDx << endl;
1454 return;
1455 break;
1456 }
1457
1458 switch(optionSigma) {
1459
1460 case 1:
1461 nSigmaPar = 1;
1462 parSigma[0] = 3.0;
1463 break;
1464 case 2:
1465 nSigmaPar = 1;
1466 parSigma[0] = 0.0655688;
1467 break;
1468 case 3:
1469 nSigmaPar = 2;
1470 parSigma[0] = 0.06;
1471 parSigma[1] = -0.001;
1472 case 4:
1473 nSigmaPar = 2;
1474 parSigma[0] = 0.06;
1475 parSigma[1] = 1.0;
1476 break;
1477 case 5:
1478 nSigmaPar = 2;
1479 parSigma[0] = 0.06;
1480 parSigma[1] = 1.0;
1481 break;
1482 case 6:
1483 nSigmaPar = 2;
1484 parSigma[0] = 0.06;
1485 parSigma[1] = 0.0;
1486 break;
1487 case 7:
1488 nSigmaPar = 3;
1489 parSigma[0] = 0.06;
1490 parSigma[1] = 0.0;
1491 parSigma[2] = 2.0;
1492 break;
1493 case 8:
1494 nSigmaPar = 3;
1495 parSigma[0] = 0.06;
1496 parSigma[1] = 0.0;
1497 parSigma[2] = 2.0;
1498 break;
1499 case 9://works for 0-5 %
1500 nSigmaPar = 3;
1501 parSigma[0] = 1.88409e-01;
1502 parSigma[1] = -3.84073e-03;
1503 parSigma[2] = 3.03110e-05;
1504 break;
1505 case 10://works for 0-5 %
1506 nSigmaPar = 6;
1507 parSigma[0]=7.96380e-03;
1508 parSigma[1]=8.09869e-05;
1509 parSigma[2]=6.29707e-06;
1510 parSigma[3]=-3.27295e-08;
1511 parSigma[4]=-1.20200e+03;
1512 parSigma[5]=-3.97089e+01;
1513 break;
1514 case 11://works for 0-5 %
1515 nSigmaPar = 6;
1516 parSigma[0]=7.96380e-03;
1517 parSigma[1]=8.09869e-05;
1518 parSigma[2]=6.29707e-06;
1519 parSigma[3]=-3.27295e-08;
1520 parSigma[4]=-1.20200e+03;
1521 parSigma[5]=-3.97089e+01;
1522 break;
1523 case 12://works for 0-5 %
1524 nSigmaPar = 3;
1525 parSigma[0]=7.96380e-03;
1526 parSigma[1]=8.09869e-05;
1527 parSigma[2]=6.29707e-06;
1528 break;
1529 case 13://works for 0-5 %
1530 nSigmaPar = 3;
1531 parSigma[0]=7.96380e-03;
1532 parSigma[1]=8.09869e-05;
1533 parSigma[2]=6.29707e-06;
1534 break;
1535 case 14:
1536 nSigmaPar = 2;
1537 parSigma[0]=7.96380e-03;
1538 parSigma[1]=8.09869e-05;
1539 break;
1540 case 15:
1541 nSigmaPar = 3;
1542 parSigma[0]=7.96380e-03;
1543 parSigma[1]=8.09869e-05;
1544 break;
1545
1546
1547
1548
1549 default:
1550
1551 cout << "optionSigma does not support option: " << optionSigma << endl;
1552 return;
1553 break;
1554 }
1555
1556 if(initialFitFileName) {
1557
1558
1559
1560 TFile* fresBB = FindFileFresh(initialFitFileName);
1561 if(!fresBB)
1562 return;
1563 TF1 *fres=(TF1 *)fresBB->Get("sigmaRes");
1564 TF1 *fBB=(TF1 *)fresBB->Get("dedxFunc");
1565
1566 Int_t nparres=fres->GetNpar();
1567 Int_t nparBB=fBB->GetNpar();
1568 for(Int_t ires=0;ires<nparres;++ires){
1569 parSigma[ires] = fres->GetParameter(ires+1);
1570 //cout<<fres->GetParameter(ires+1)<<endl;
1571 }
1572 for(Int_t iBB=0;iBB<nparBB;++iBB){
1573 parDeDx[iBB] = fBB->GetParameter(iBB+1);
1574 cout<<parDeDx[iBB]<<endl;
1575
1576 }
1577 //Esto lo puse apenas
1578 fixPlateau = parDeDx[4];
1579 fixMIP = fBB->Eval(3.5);
1580 cout<<"Hello!!"<<" fixPlateau="<<fixPlateau<<" fixMIP="<<fixMIP<<endl;
1581
1582 //fixMIP = fitPar->MIP;
1583 //return;
1584 }
1585
1586
1587
1588 const Int_t nGlobalParams = 2 // binStart, binStop,
1589 + 2 + nDeDxPar // optionDeDx, nDeDxPar, dedxpar0 ....
1590 + 2 + nSigmaPar; // nSigmaPar, sigmapar0 .....
1591
1592 const Int_t nParams = nBins*nLocalParams + nGlobalParams;
1593
1594
1595 TF2* fitAll = new TF2("fitAll", fitf3G, pStart, pStop, 30, 90, nParams);
1596 Double_t parametersIn[nParams];
1597
1598 parametersIn[0] = binStart;
1599 fitAll->SetParName(0,"binStart");
1600 fitAll->FixParameter(0, parametersIn[0]);
1601
1602 parametersIn[1] = binStop;
1603 fitAll->SetParName(1,"binStop");
1604 fitAll->FixParameter(1, parametersIn[1]);
1605
1606 // dE/dx parameters
1607 parametersIn[2] = nDeDxPar;
1608 fitAll->SetParName(2,"nDeDxPar");
1609 fitAll->FixParameter(2, parametersIn[2]);
1610
1611 parametersIn[3] = optionDeDx;
1612 fitAll->SetParName(3,"optionDeDx");
1613 fitAll->FixParameter(3, parametersIn[3]);
1614
1615 for(Int_t i = 0; i < nDeDxPar; i++) {
1616
1617 parametersIn[4+i] = parDeDx[i];
1618 fitAll->SetParName(4+i,Form("dEdXpar%d", i));
1619 Int_t bit = TMath::Nint(TMath::Power(2, i));
1620 if(fixParametersDedx>=0)
1621 if (fixParametersDedx==0 || fixParametersDedx&bit) {
1622
1623 cout << "Fixing dE/dx parameter " << i << endl;
1624 fitAll->FixParameter(4+i, parametersIn[4+i]);
1625 }
1626
1627 if(fixAllParBB) {
1628
1629 fitAll->FixParameter(4+i, parametersIn[4+i]);
1630 }
1631 }
1632
1633 // sigma parameters
1634 parametersIn[4+nDeDxPar] = nSigmaPar;
1635 fitAll->SetParName(4+nDeDxPar,"nSigmaPar");
1636 fitAll->FixParameter(4+nDeDxPar, parametersIn[4+nDeDxPar]);
1637
1638 parametersIn[5+nDeDxPar] = optionSigma;
1639 fitAll->SetParName(5+nDeDxPar,"optionSigma");
1640 fitAll->FixParameter(5+nDeDxPar, parametersIn[5+nDeDxPar]);
1641
1642 for(Int_t i = 0; i < nSigmaPar; i++) {
1643
1644 parametersIn[6+nDeDxPar+i] = parSigma[i];
1645 fitAll->SetParName(6+nDeDxPar+i,Form("sigmapar%d", i));
1646 Int_t bit = TMath::Nint(TMath::Power(2, i));
1647 if(fixParametersSigma>=0)
1648 if (fixParametersSigma==0 || fixParametersSigma&bit) {
1649
1650 cout << "Fixing sigma parameter " << i << endl;
1651 fitAll->FixParameter(6+nDeDxPar+i, parametersIn[6+nDeDxPar+i]);
1652 }
1653
1654 if(fixAllParSigma) {
1655
1656 fitAll->FixParameter(6+nDeDxPar+i, parametersIn[6+nDeDxPar+i]);
1657 }
1658 }
1659
1660 // Initial yields
1661
1662 for(Int_t bin = binStart; bin <= binStop; bin++) {
1663
1664 TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY("hDeDxVsPProj", bin, bin);
1665
1666 const Int_t offset = nGlobalParams + (bin - binStart)*nLocalParams;
1667 const Double_t all = hDeDxVsPProj->Integral();
1668
1669 fitAll->SetParName(offset + 0, Form("piYield%d", bin));
1670 fitAll->SetParName(offset + 1, Form("kYield%d", bin));
1671 fitAll->SetParName(offset + 2, Form("pYield%d", bin));
1672 fitAll->SetParLimits(offset + 0, 0, 10*all);
1673 fitAll->SetParLimits(offset + 2, 0, 10*all);
1674
1675 if(fixKtoZero) {
1676 parametersIn[offset + 0] = (all)*0.5;
1677 parametersIn[offset + 1] = (all)*0.0;
1678 fitAll->FixParameter(offset + 1, 0.0);
1679 parametersIn[offset + 2] = (all)*0.5;
1680 } else {
1681 parametersIn[offset + 0] = (all)*0.6;
1682 parametersIn[offset + 1] = (all)*0.3;
1683 fitAll->SetParLimits(offset + 1, 0, 10*all);
1684 parametersIn[offset + 2] = (all)*0.1;
1685 }
1686 // fitAll->SetParLimits(offset + 0, 0., histdEdxvsPpy->GetEntries());
1687 // fitAll->SetParLimits(offset + 1, 0., histdEdxvsPpy->GetEntries());
1688 // fitAll->SetParLimits(offset + 2, 0., histdEdxvsPpy->GetEntries());
1689 }
1690
1691 fitAll->SetParameters(parametersIn);
1692 fitAll->Print();
1693
1694 Bool_t converged = kFALSE;
1695
1696 if(performFit) {
1697 for(Int_t i = 0; i < 4; i++) {
1698 TFitResultPtr fitResultPtr = hDeDxVsP->Fit(fitAll, "0S", "", pStart+0.01, pStop-0.01);
1699 if(!fitResultPtr->Status()) {
1700 // if(!fitResultPtr->Status() && !fitResultPtr->CovMatrixStatus()) {
1701
1702 converged = kTRUE;
1703 break;
1704 }
1705 }
1706 }
1707 // else we just draw how the results looks with the input parameters
1708
1709 Double_t parametersOut[nParams];
1710 fitAll->GetParameters(parametersOut);
1711 const Double_t* parameterErrorsOut = fitAll->GetParErrors();
1712
1713 // overlay the fitfunction
1714
1715
1716 TF1* fDeDxPi = new TF1("fDeDxPi", FitFunc, 0, 50, nDeDxPar+1); // +1 because of option!
1717 fDeDxPi->SetParameters(&parametersOut[3]);
1718 fDeDxPi->SetLineWidth(2);
1719 cDeDxVsP->cd();
1720 fDeDxPi->Draw("SAME");
1721 cDeDxVsPLogX->cd();
1722 fDeDxPi->Draw("SAME");
1723
1724 TF1* fDeDxK = new TF1("fDeDxK", FitFunc, 0, 50, nDeDxPar+1);
1725 fDeDxK->SetParameters(&parametersOut[3]);
1726 fDeDxK->SetParameter(0, fDeDxK->GetParameter(0)+10);
1727 fDeDxK->SetLineWidth(2);
1728 cDeDxVsP->cd();
1729 fDeDxK->Draw("SAME");
1730 cDeDxVsPLogX->cd();
1731 fDeDxK->Draw("SAME");
1732
1733 TF1* fDeDxP = new TF1("fDeDxP", FitFunc, 0, 50, nDeDxPar+1);
1734 fDeDxP->SetParameters(&parametersOut[3]);
1735 fDeDxP->SetParameter(0, fDeDxP->GetParameter(0)+20);
1736 fDeDxP->SetLineWidth(2);
1737 cDeDxVsP->cd();
1738 fDeDxP->Draw("SAME");
1739 cDeDxVsPLogX->cd();
1740 fDeDxP->Draw("SAME");
1741
1742 TF1* fDeDxE = new TF1("fDeDxE", FitFunc, 0, 50, nDeDxPar+1);
1743 fDeDxE->SetParameters(&parametersOut[3]);
1744 fDeDxE->SetParameter(0, fDeDxE->GetParameter(0)+30);
1745 fDeDxE->SetLineWidth(2);
1746 cDeDxVsP->cd();
1747 fDeDxE->Draw("SAME");
1748 cDeDxVsPLogX->cd();
1749 fDeDxE->Draw("SAME");
1750
1751 TF1* fSigma = new TF1("fSigma", SigmaFunc, 0, 50, nSigmaPar+1);
1752 fSigma->SetParameters(&parametersOut[5 + nDeDxPar]);
1753
1754 // fitAll->Draw("same cont3");
1755
1756 CreateDir(Form("results/%s/calibdedx_%d%d",endingCent[i_cent],etaLow,etaHigh));
1757
1758 //CreateDir(Form("results/calibdedx_%d%d",etaLow,etaHigh));
1759
1760 cDeDxVsP->cd();
1761 cDeDxVsP->Modified();
1762 cDeDxVsP->Update();
1763 gROOT->ProcessLine(".x drawText.C");
1764 cDeDxVsP->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p.gif",endingCent[i_cent],etaLow,etaHigh));
1765 cDeDxVsP->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p.pdf",endingCent[i_cent],etaLow,etaHigh));
1766
1767 cDeDxVsPLogX->cd();
1768 cDeDxVsPLogX->Modified();
1769 cDeDxVsPLogX->Update();
1770 gROOT->ProcessLine(".x drawText.C");
1771 cDeDxVsPLogX->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p_logx.gif",endingCent[i_cent],etaLow,etaHigh));
1772 cDeDxVsPLogX->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p_logx.pdf",endingCent[i_cent],etaLow,etaHigh));
1773
1774 //cross check
1775 TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
1776
1777 cFits->Clear();
1778 cFits->Divide(7, 5);
1779
1780 TF1* pion = new TF1("pion", "gausn", 30, 90);
1781 pion->SetLineWidth(2);
1782 pion->SetLineColor(kRed);
1783 TF1* kaon = new TF1("kaon", "gausn", 30, 90);
1784 kaon->SetLineWidth(2);
1785 kaon->SetLineColor(kGreen);
1786 TF1* proton = new TF1("proton", "gausn", 30, 90);
1787 proton->SetLineWidth(2);
1788 proton->SetLineColor(kBlue);
1789 TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)", 30, 90);
1790 total->SetLineColor(kBlack);
1791 total->SetLineWidth(2);
1792 total->SetLineStyle(2);
1793
1794 TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);
1795 legend->SetBorderSize(0);
1796 legend->SetFillColor(0);
1797 legend->AddEntry(total, "3-Gauss fit", "L");
1798 legend->AddEntry(pion, "#pi", "L");
1799 legend->AddEntry(kaon, "K", "L");
1800 legend->AddEntry(proton, "p", "L");
1801
1802
1803 TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
1804 cSingleFit->SetLogy(1);
1805
1806 TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
1807 hPionRatio->Reset();
1808 hPionRatio->GetXaxis()->SetRangeUser(pStart+0.001, pStop-0.001);
1809 hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
1810 hPionRatio->SetTitle("particle fractions; p [GeV/c]; particle fraction");
1811 TH1D* hKaonRatio = (TH1D*)hPionRatio->Clone("hKaonRatio");
1812 TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
1813
1814 for(Int_t bin = binStart; bin <= binStop; bin++){
1815
1816 cout << "Making projection for bin: " << bin << endl;
1817
1818 const Int_t j = bin-binStart;
1819 cFits->cd();
1820 cFits->cd(j + 1);
1821
1822 TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY(Form("hDeDxVsPProj%d", bin), bin, bin);
1823 // hDeDxVsPProj->GetXaxis()->SetRangeUser(40, 90);
1824 hDeDxVsPProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
1825 hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
1826 hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
1827 hDeDxVsPProj->Draw();
1828
1829 const Int_t offset = nGlobalParams + j*nLocalParams;
1830 const Double_t p = hDeDxVsP->GetXaxis()->GetBinCenter(bin);
1831 const Double_t pKeff = p*piMass/kMass; // corresponding p of a pion with same dE/dx
1832 const Double_t pPeff = p*piMass/pMass; // corresponding p of a pion with same dE/dx
1833 const Double_t meanDeDxPi = fDeDxPi->Eval(p);
1834 const Double_t meanDeDxK = fDeDxPi->Eval(pKeff);
1835 const Double_t meanDeDxP = fDeDxPi->Eval(pPeff);
1836 Double_t gausParams[9] = {
1837 parametersOut[offset + 0],
1838 meanDeDxPi,
1839 fSigma->Eval(meanDeDxPi) ,
1840 parametersOut[offset + 1],
1841 meanDeDxK,
1842 fSigma->Eval(meanDeDxK) ,
1843 parametersOut[offset + 2],
1844 meanDeDxP,
1845 fSigma->Eval(meanDeDxP) ,
1846 };
1847
1848 for(Int_t i = 0; i < 9; i++)
1849 cout << gausParams[i] << ", ";
1850
1851 cout << endl;
1852
1853 pion->SetParameters(&gausParams[0]);
1854 pion->DrawCopy("same");
1855 Double_t all = hDeDxVsPProj->Integral();
1856 hPionRatio->SetBinContent(bin, parametersOut[offset + 0]/all);
1857 hPionRatio->SetBinError(bin, parameterErrorsOut[offset + 0]/all);
1858
1859 kaon->SetParameters(&gausParams[3]);
1860 kaon->DrawCopy("same");
1861 hKaonRatio->SetBinContent(bin, parametersOut[offset + 1]/all);
1862 hKaonRatio->SetBinError(bin, parameterErrorsOut[offset + 1]/all);
1863
1864 proton->SetParameters(&gausParams[6]);
1865 proton->DrawCopy("same");
1866 hProtonRatio->SetBinContent(bin, parametersOut[offset + 2]/all);
1867 hProtonRatio->SetBinError(bin, parameterErrorsOut[offset + 2]/all);
1868
1869 total->SetParameters(gausParams);
1870 total->DrawCopy("same");
1871
1872 cSingleFit->cd();
1873 cSingleFit->Clear();
1874 // cSingleFit->SetLogy();
1875 hDeDxVsPProj->Draw();
1876 pion->DrawCopy("same");
1877 kaon->DrawCopy("same");
1878 proton->DrawCopy("same");
1879 total->DrawCopy("same");
1880
1881 gROOT->ProcessLine(".x drawText.C(2)");
1882
1883
1884 cSingleFit->SaveAs(Form("results/%s/calibdedx_%d%d/ptspectrum_bin%d.gif",endingCent[i_cent], etaLow,etaHigh, bin));
1885 cSingleFit->SaveAs(Form("results/%s/calibdedx_%d%d/ptspectrum_bin%d.pdf",endingCent[i_cent],etaLow, etaHigh, bin));
1886 // legend->Draw();
1887
1888
1889 }
1890
1891 TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p", 600, 400);
1892 cRatio->Clear();
1893 hPionRatio->SetMaximum(0.8);
1894 hPionRatio->SetMarkerStyle(20);
1895 hPionRatio->SetMarkerColor(2);
1896 hPionRatio->Draw("P E");
1897
1898 hKaonRatio->SetMarkerStyle(20);
1899 hKaonRatio->SetMarkerColor(3);
1900 hKaonRatio->Draw("SAME P E");
1901
1902 hProtonRatio->SetMarkerStyle(20);
1903 hProtonRatio->SetMarkerColor(4);
1904 hProtonRatio->Draw("SAME P E");
1905 gROOT->ProcessLine(".x drawText.C(2)");
1906 cRatio->SaveAs(Form("results/%s/calibdedx_%d%d/particle_ratios.gif",endingCent[i_cent],etaLow,etaHigh));
1907 cRatio->SaveAs(Form("results/%s/calibdedx_%d%d/particle_ratios.pdf",endingCent[i_cent],etaLow,etaHigh));
1908
1909
1910
1911 //
1912 // Store the <dE/dx> parameters in a file that we can get them back to use for the Delta-pi!
1913 //
1914 DeDxFitInfo* fitInfo = new DeDxFitInfo();
1915 fitInfo->MIP = fixMIP;
1916 //fitInfo->plateau = 80;
1917 fitInfo->optionDeDx = optionDeDx;
1918 fitInfo->nDeDxPar = nDeDxPar;
1919 for(Int_t i = 0; i < nDeDxPar; i++) {
1920 fitInfo->parDeDx[i] = fDeDxPi->GetParameter(i+1); // 1st parameter is option
1921 }
1922 //fitInfo->plateau = fDeDxPi->GetParameter(5);//lo puse 16/09/13
1923 fitInfo->optionSigma = optionSigma;
1924 fitInfo->nSigmaPar = nSigmaPar;
1925 for(Int_t i = 0; i < nSigmaPar; i++) {
1926 cout<<"my sugma param="<<fSigma->GetParameter(i+1)<<endl;
1927 fitInfo->parSigma[i] = fSigma->GetParameter(i+1); // 1st parameter is option
1928 }
1929
1930
1931 fitInfo->Print();
1932
1933
1934 if(converged) {
1935
1936 cout << "Fit converged and error matrix was ok" << endl;
1937 } else {
1938
1939 cout << "WARNING!!!!!!!!!!!!!!! Fit did not converge, or error matrix was not ok!" << endl;
1940 }
1941
1942
1943 CreateDir(Form("fitparameters/%s",endingCent[i_cent]));
1944 //cout<<"flag 0"<<" gSystem->BaseName(calibFileName))="<< gSystem->BaseName(calibFileName)<<endl;
1945 TFile* outFile = new TFile(Form("fitparameters/%s/%d%d_dataPbPb.root", endingCent[i_cent], etaLow, etaHigh), "RECREATE");
1946 outFile->cd();
1947 fitInfo->Write("fitInfo");
1948 cout<<"flag 1"<<endl;
1949
1950 outFile->Close();
1951
1952 if(converged) {
1953
1954 cout << "Fit converged and error matrix was ok" << endl;
1955 } else {
1956
1957 cout << "WARNING!!!!!!!!!!!!!!! Fit did not converge, or error matrix was not ok!" << endl;
1958 }
1959}
1960
1961
1962
1963
1964
1965void MakeFitsV0s(const Char_t * fileInputName, const Char_t * fileE, const Char_t *outDir, Int_t index_eta){
1966 /*
1967 gStyle->SetPalette(1);
1968 gStyle->SetCanvasColor(10);
1969 gStyle->SetFrameFillColor(10);
1970 gStyle->SetOptStat(0);
1971 gStyle->SetOptTitle(0);
1972 gStyle->SetOptFit(0);
1973 */
1974 TFile *f12=TFile::Open(fileE);
1975 TF1 *fPlateau = (TF1 *)f12->Get(Form("fGauss%d",index_eta));
1976 /*
1977 TH1D * hEV0s=(TH1D *)f12->Get("hdedx_vs_p_e");
1978 TGraphErrors *gResE=(TGraphErrors *)f12->Get("gRes_Electrons");
1979 */
1980
1981 const Char_t * endNamePos[4]={"02","24","46","68"};
1982 const Char_t * endNameNeg[4]={"20","42","64","86"};
1983
1984
1985
1986 TFile *f1=TFile::Open(fileInputName);
1987 TList * list = (TList *)f1->Get("outputdedx");
1988
1989
1990 //positive eta
1991 TH2D *hPi2DPos=(TH2D *)list->FindObject(Form("histPiTof%s",endNamePos[index_eta]));
1992 TH2D *hPi2DV0sPos=(TH2D *)list->FindObject(Form("histPiV0%s",endNamePos[index_eta]));
1993 TH2D *hP2DPos=(TH2D *)list->FindObject(Form("histPV0%s",endNamePos[index_eta]));
1994 TH1D *hPi2DpPos=(TH1D *)list->FindObject(Form("histPiTof1D%s",endNamePos[index_eta]));
1995 TH1D *hPi2DV0spPos=(TH1D *)list->FindObject(Form("histPiV01D%s",endNamePos[index_eta]));
1996 TH1D *hP2DpPos=(TH1D *)list->FindObject(Form("histPV01D%s",endNamePos[index_eta]));
1997
1998
1999 gSystem->Exec("rm *.root");
2000
2001 TFile *fpos = new TFile("histtmppos.root","recreate");
2002 fpos->cd();
2003 hPi2DPos->SetName(Form("histPiTof%s",endNamePos[index_eta]));
2004 hPi2DV0sPos->SetName(Form("histPiV0%s",endNamePos[index_eta]));
2005 hP2DPos->SetName(Form("histPV0%s",endNamePos[index_eta]));
2006 hPi2DpPos->SetName(Form("histPiTof1D%s",endNamePos[index_eta]));
2007 hPi2DV0spPos->SetName(Form("histPiV01D%s",endNamePos[index_eta]));
2008 hP2DpPos->SetName(Form("histPV01D%s",endNamePos[index_eta]));
2009 hPi2DPos->Write();
2010 hPi2DV0sPos->Write();
2011 hP2DPos->Write();
2012 hPi2DpPos->Write();
2013 hPi2DV0spPos->Write();
2014 hP2DpPos->Write();
2015 fpos->Close();
2016
2017
2018 //negative eta
2019 TH2D *hPi2DNeg=(TH2D *)list->FindObject(Form("histPiTof%s",endNameNeg[index_eta]));
2020 TH2D *hPi2DV0sNeg=(TH2D *)list->FindObject(Form("histPiV0%s",endNameNeg[index_eta]));
2021 TH2D *hP2DNeg=(TH2D *)list->FindObject(Form("histPV0%s",endNameNeg[index_eta]));
2022 TH1D *hPi2DpNeg=(TH1D *)list->FindObject(Form("histPiTof1D%s",endNameNeg[index_eta]));
2023 TH1D *hPi2DV0spNeg=(TH1D *)list->FindObject(Form("histPiV01D%s",endNameNeg[index_eta]));
2024 TH1D *hP2DpNeg=(TH1D *)list->FindObject(Form("histPV01D%s",endNameNeg[index_eta]));
2025
2026
2027
2028 TFile *fneg = new TFile("histtmpneg.root","recreate");
2029 fneg->cd();
2030 hPi2DNeg->SetName(Form("histPiTof%s",endNamePos[index_eta]));
2031 hPi2DV0sNeg->SetName(Form("histPiV0%s",endNamePos[index_eta]));
2032 hP2DNeg->SetName(Form("histPV0%s",endNamePos[index_eta]));
2033 hPi2DpNeg->SetName(Form("histPiTof1D%s",endNamePos[index_eta]));
2034 hPi2DV0spNeg->SetName(Form("histPiV01D%s",endNamePos[index_eta]));
2035 hP2DpNeg->SetName(Form("histPV01D%s",endNamePos[index_eta]));
2036 hPi2DNeg->Write();
2037 hPi2DV0sNeg->Write();
2038 hP2DNeg->Write();
2039 hPi2DpNeg->Write();
2040 hPi2DV0spNeg->Write();
2041 hP2DpNeg->Write();
2042 fneg->Close();
2043
2044 //pos+neg eta
2045 TH2D *hPi2D=0;
2046 TH2D *hPi2DV0s=0;
2047 TH2D *hP2D=0;
2048 TH1D *hPi2Dp=0;
2049 TH1D *hPi2DV0sp=0;
2050 TH1D *hP2Dp=0;
2051
2052
2053 gSystem->Exec(Form("hadd tmp%d.root histtmppos.root histtmpneg.root",index_eta));
2054
2055 TFile *ffinalv0 = TFile::Open(Form("tmp%d.root",index_eta));
2056 hPi2D=(TH2D *)ffinalv0->Get(Form("histPiTof%s",endNamePos[index_eta]));
2057 hPi2DV0s=(TH2D *)ffinalv0->Get(Form("histPiV0%s",endNamePos[index_eta]));
2058 hP2D=(TH2D *)ffinalv0->Get(Form("histPV0%s",endNamePos[index_eta]));
2059 hPi2Dp=(TH1D *)ffinalv0->Get(Form("histPiTof1D%s",endNamePos[index_eta]));
2060 hPi2DV0sp=(TH1D *)ffinalv0->Get(Form("histPiV01D%s",endNamePos[index_eta]));
2061 hP2Dp=(TH1D *)ffinalv0->Get(Form("histPV01D%s",endNamePos[index_eta]));
2062
2063
2064
2065 TList *l1=new TList();
2066 l1->SetOwner();
2067
2068
2069 TGraphErrors* graphRes = new TGraphErrors();
2070 Int_t nERes = 0;
2071
2072 TGraphErrors* graphBB = new TGraphErrors();
2073 Int_t nBB = 0;
2074
2075 TGraphErrors* graphBBpitof = new TGraphErrors();
2076 Int_t nBBpitof = 0;
2077 TGraphErrors* graphBBpiv0s = new TGraphErrors();
2078 Int_t nBBpiv0s = 0;
2079 TGraphErrors* graphBBpv0s = new TGraphErrors();
2080 Int_t nBBpv0s = 0;
2081
2082
2083
2084
2085 TF1* fFit1D = new TF1("fFit1D", "gausn(0)", 40, 90);
2086
2087
2088
2089 //TF1* pionsv0s = new TF1("pionsv0s", "gausn", 56, 80);
2090 TF1* pionsv0s = new TF1("pionsv0s", "gausn", 40, 80);
2091 pionsv0s->SetLineColor(kRed);
2092 pionsv0s->SetLineWidth(2);
2093 pionsv0s->SetLineStyle(2);
2094
2095 //TF1* protonsv0s = new TF1("protonsv0s", "gausn", 40, 55);
2096 TF1* protonsv0s = new TF1("protonsv0s", "gausn", 40, 80);
2097 protonsv0s->SetLineColor(kBlue);
2098 protonsv0s->SetLineWidth(2);
2099 protonsv0s->SetLineStyle(2);
2100
2101
2102 TF1* total = new TF1("total", "gausn(0)+gausn(3)", 40, 80);
2103 total->SetLineColor(1);
2104 total->SetLineWidth(2);
2105 total->SetLineStyle(2);
2106
2107
2108
2109 Int_t nbins=hPi2D->GetXaxis()->GetNbins();
2110
2111
2112
2113 for(Int_t bin=1;bin<=nbins;++bin){
2114
2115
2116
2117 TH1D *hproy_pion=0;
2118 TH1D *hproy_proton=0;
2119
2120 Double_t p=hPi2D->GetXaxis()->GetBinCenter(bin);
2121 Double_t lowp = hPi2D->GetXaxis()->GetBinLowEdge(bin);
2122 Double_t widthp = hPi2D->GetXaxis()->GetBinWidth(bin);
2123 hPi2Dp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
2124 Double_t meanp= hPi2Dp->GetMean();
2125 Double_t meanpE= hPi2Dp->GetMeanError();
2126
2127
2128 cout<<"p="<<p<<endl;
2129
2130 if(p>2.0&&p<9.0){
2131
2132 cout<<bin<<endl;
2133
2134 hproy_pion=(TH1D *)hPi2D->ProjectionY(Form("hproy_pionTof_%d",bin),bin,bin);
2135
2136 cout<<bin<<endl;
2137
2138 hproy_pion->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2139 hPi2D->GetXaxis()->GetBinLowEdge(bin),
2140 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2141 cout<<bin<<endl;
2142 hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
2143
2144 cout<<"nbinsproj="<<hproy_proton->GetNbinsX()<<endl;
2145 //hproy_pion->Draw();
2146
2147
2148 TSpectrum *s = new TSpectrum(2);
2149 Int_t nfound = s->Search(hproy_pion,2,"",0.15);
2150 cout<<"!!!!!!!!!!!!!! nfound="<<nfound<<endl;
2151 Int_t npeaks = 0;
2152 Float_t *xpeaks = s->GetPositionX();
2153
2154
2155 for (Int_t p1=0;p1<nfound;p1++) {
2156 //Float_t xp = xpeaks[p1];
2157 npeaks++;
2158 }
2159
2160 if(npeaks==0)
2161 continue;
2162 if(npeaks==1){
2163
2164 fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[0], hproy_pion->GetRMS());
2165 hproy_pion->Fit(fFit1D, "Q");
2166 hproy_pion->Fit(fFit1D, "Q", "", xpeaks[0]-hproy_pion->GetRMS(),
2167 xpeaks[0]+hproy_pion->GetRMS());
2168
2169 hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2170 fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2171
2172
2173
2174
2175 graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2176 graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2177 nBB++;
2178
2179 graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2180 graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2181 nBBpitof++;
2182
2183 if(p>3&&p<9){
2184 graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
2185 graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
2186 nERes++;
2187 }
2188
2189
2190
2191
2192 }
2193 else if(npeaks==2){
2194 Int_t bin1 = hproy_pion->GetXaxis()->FindBin(xpeaks[0]);
2195 Int_t bin2 = hproy_pion->GetXaxis()->FindBin(xpeaks[1]);
2196
2197 Int_t bin1C = hproy_pion->GetBinContent(bin1);
2198 Int_t bin2C = hproy_pion->GetBinContent(bin2);
2199
2200 cout<<"bin1="<<bin1<<" content="<<bin1C<<endl;
2201 cout<<"bin2="<<bin2<<" content="<<bin2C<<endl;
2202
2203
2204 if(bin==24){
2205
2206 hproy_pion->Fit(fFit1D, "R", "", 65, 85);
2207 graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2208 graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2209 nBB++;
2210
2211 graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2212 graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2213 nBBpitof++;
2214
2215
2216
2217 }
2218
2219
2220 else if(bin1C>bin2C){//pion signal larger
2221
2222
2223 fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[0], hproy_pion->GetRMS());
2224 hproy_pion->Fit(fFit1D, "Q");
2225 hproy_pion->Fit(fFit1D, "Q", "", xpeaks[0]-hproy_pion->GetRMS(),
2226 xpeaks[0]+hproy_pion->GetRMS());
2227
2228
2229
2230
2231 hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2232 fFit1D->GetParameter(1)+2.5*fFit1D->GetParameter(2));
2233
2234
2235
2236
2237 cout<<"hello"<<endl;
2238 //protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[1], hproy_proton->GetRMS());
2239
2240
2241 hproy_proton->Fit(protonsv0s, "Q");
2242
2243
2244 hproy_proton->Fit(protonsv0s, "Q", "", xpeaks[1]-hproy_proton->GetRMS(),
2245 xpeaks[1]+hproy_proton->GetRMS());
2246
2247
2248 hproy_proton->Fit(protonsv0s, "R", "", xpeaks[1]-2*protonsv0s->GetParameter(2),
2249 xpeaks[1]+2*protonsv0s->GetParameter(2));
2250
2251 protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[1], protonsv0s->GetParameter(2));
2252 hproy_proton->Fit(protonsv0s, "R", "", xpeaks[1]-2.0*protonsv0s->GetParameter(2),
2253 xpeaks[1]+2.0*protonsv0s->GetParameter(2));
2254
2255 total->SetParameter(1,protonsv0s->GetParameter(1));
2256 total->SetParameter(2,protonsv0s->GetParameter(2));
2257 total->SetParameter(4,fFit1D->GetParameter(1));
2258 total->SetParameter(5,fFit1D->GetParameter(2));
2259 hproy_pion->Fit(total,"R");
2260
2261 Float_t reducedCh2=total->GetChisquare()/total->GetNDF();
2262 if(reducedCh2>1.0){
2263 hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-2*fFit1D->GetParameter(2),
2264 fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2265 //hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-0.5*fFit1D->GetParameter(2),
2266 hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2267 fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2268
2269
2270 graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2271 graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2272 nBB++;
2273
2274
2275 graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2276 graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2277 nBBpitof++;
2278
2279 if(p>3.0&&p<9){
2280 graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
2281 graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
2282 nERes++;
2283 }
2284
2285
2286
2287 }else{
2288
2289 graphBB->SetPoint(nBB,meanp/0.14, total->GetParameter(4));
2290 graphBB->SetPointError(nBB,meanpE, total->GetParError(4));
2291 nBB++;
2292
2293 graphBBpitof->SetPoint(nBBpitof,meanp/0.14, total->GetParameter(4));
2294 graphBBpitof->SetPointError(nBBpitof,meanpE, total->GetParError(4));
2295 nBBpitof++;
2296
2297 if(p>3.0&&p<9){
2298 graphRes->SetPoint(nERes,total->GetParameter(4),total->GetParameter(5)/total->GetParameter(4));
2299
2300 graphRes->SetPointError(nERes,total->GetParError(4),total->GetParError(5)/total->GetParameter(4));
2301 nERes++;
2302 }
2303
2304
2305
2306
2307 }
2308
2309
2310
2311
2312 }else{
2313
2314 fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[1], hproy_pion->GetRMS());
2315 hproy_pion->Fit(fFit1D, "Q");
2316 hproy_pion->Fit(fFit1D, "Q", "", xpeaks[1]-hproy_pion->GetRMS(),
2317 xpeaks[0]+hproy_pion->GetRMS());
2318
2319 hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-2*fFit1D->GetParameter(2),
2320 fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2321
2322
2323 protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[0], hproy_proton->GetRMS());
2324 hproy_proton->Fit(protonsv0s, "Q");
2325 hproy_proton->Fit(protonsv0s, "Q", "", xpeaks[0]-hproy_proton->GetRMS(),
2326 xpeaks[0]+hproy_proton->GetRMS());
2327
2328 hproy_proton->Fit(protonsv0s, "R", "", xpeaks[0]-2*protonsv0s->GetParameter(2),
2329 xpeaks[0]+2*protonsv0s->GetParameter(2));
2330
2331 protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[0], protonsv0s->GetParameter(2));
2332 hproy_proton->Fit(protonsv0s, "R", "", xpeaks[0]-2*protonsv0s->GetParameter(2),
2333 xpeaks[0]+2*protonsv0s->GetParameter(2));
2334
2335
2336 //graphBB->SetPoint(nBB,p/0.14, xpeaks[1]);
2337 //graphBB->SetPointError(nBB,0, fFit1D->GetParError(1));
2338 //nBB++;
2339 total->SetParameter(1,protonsv0s->GetParameter(1));
2340 total->SetParameter(2,protonsv0s->GetParameter(2));
2341 total->SetParameter(4,fFit1D->GetParameter(1));
2342 total->SetParameter(5,fFit1D->GetParameter(2));
2343 hproy_pion->Fit(total,"R");
2344
2345 //graphBB->SetPoint(nBB,p/0.14, xpeaks[1]);
2346 //graphBB->SetPointError(nBB,0, total->GetParError(4));
2347 graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2348 graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2349 nBB++;
2350
2351 graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2352 graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2353 nBBpitof++;
2354
2355
2356 if(p>3.0&&p<9){
2357 graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
2358 graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
2359 nERes++;
2360 }
2361
2362
2363
2364 }
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376 }
2377
2378 l1->Add(hproy_pion);
2379
2380
2381 }
2382 else
2383 continue;
2384
2385 }
2386
2387
2388
2389 graphRes->SetTitle(";#LT dE/dx #GT;#sigma/#LT dE/dx #GT");
2390 graphRes->SetName("gRes_PionsTof");
2391
2392 graphBBpitof->SetTitle(";#beta*#gamma; #LT dE/dx #GT");
2393 graphBBpitof->SetName("BB_piTOF");
2394 graphBBpitof->SetMarkerColor(kGreen+2);
2395 graphBBpitof->SetLineColor(kGreen+2);
2396 graphBBpitof->SetMarkerStyle(21);
2397
2398 l1->Add(graphRes);
2399 l1->Add(graphBBpitof);
2400
2401
2402
2403 //pi by V0s
2404
2405 TGraphErrors* graphResV0s = new TGraphErrors();
2406 nERes = 0;
2407
2408 nbins=hPi2DV0s->GetXaxis()->GetNbins();
2409
2410
2411
2412 for(Int_t bin=1;bin<=nbins;++bin){
2413
2414
2415
2416 TH1D *hproy_pionV0s=0;
2417 TH1D *hproy_proton=0;
2418
2419 Double_t p=hPi2DV0s->GetXaxis()->GetBinCenter(bin);
2420
2421 Double_t lowp = hPi2DV0s->GetXaxis()->GetBinLowEdge(bin);
2422 Double_t widthp = hPi2DV0s->GetXaxis()->GetBinWidth(bin);
2423 hPi2DV0sp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
2424 Double_t meanp= hPi2DV0sp->GetMean();
2425 Double_t meanpE= hPi2DV0sp->GetMeanError();
2426
2427
2428
2429 if(p>2.0&&p<9){
2430
2431 hproy_pionV0s=(TH1D *)hPi2DV0s->ProjectionY(Form("hproy_pionV0s_%d",bin),bin,bin);
2432 hproy_pionV0s->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2433 hPi2DV0s->GetXaxis()->GetBinLowEdge(bin),
2434 hPi2DV0s->GetXaxis()->GetBinUpEdge(bin)));
2435
2436 fFit1D->SetParameters(hproy_pionV0s->GetEntries(), hproy_pionV0s->GetMean(), hproy_pionV0s->GetRMS());
2437 hproy_pionV0s->Fit(fFit1D, "Q");
2438 hproy_pionV0s->Fit(fFit1D, "Q", "", hproy_pionV0s->GetMean()-hproy_pionV0s->GetRMS(),
2439 hproy_pionV0s->GetMean()+hproy_pionV0s->GetRMS());
2440
2441
2442 // hproy_pionV0s->Fit(fFit1D, "Q", "", fFit1D->GetParameter(1)-1.5*fFit1D->GetParameter(2),
2443 // fFit1D->GetParameter(1)+3.0*fFit1D->GetParameter(2));
2444
2445 hproy_pionV0s->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2446 fFit1D->GetParameter(1)+2.0*fFit1D->GetParameter(2));
2447
2448
2449 hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
2450 hproy_proton->SetTitle(Form("protonsv0s: %.1f < p < %.1f GeV/c; dE/dx; Entries",
2451 hPi2D->GetXaxis()->GetBinLowEdge(bin),
2452 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2453
2454 //protonsv0s->SetParameters(hproy_proton->GetEntries(), hproy_proton->GetMean(), hproy_proton->GetRMS());
2455 hproy_proton->Fit(protonsv0s,"R");
2456 //hproy_proton->Fit(protonsv0s,"R", "",protonsv0s->GetParameter(1)-3.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+1.0*protonsv0s->GetParameter(2));
2457 hproy_proton->Fit(protonsv0s,"R","",protonsv0s->GetParameter(1)-2.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+0.8*protonsv0s->GetParameter(2));
2458
2459 total->SetParameter(0,0.1*hproy_pionV0s->GetEntries());
2460 total->SetParameter(1,protonsv0s->GetParameter(1));
2461 total->SetParameter(2,protonsv0s->GetParameter(2));
2462 total->SetParameter(3,0.9*hproy_pionV0s->GetEntries());
2463 total->SetParameter(4,fFit1D->GetParameter(1));
2464 total->SetParameter(5,fFit1D->GetParameter(2));
2465
2466 hproy_pionV0s->Fit(total,"0R");
2467
2468 cout<<"p="<<p<<" yield0="<<total->GetParameter(0)<<" yield1="<<total->GetParameter(3)<<endl;
2469 cout<<"p="<<p<<" mean0="<<total->GetParameter(1)<<" mean1="<<total->GetParameter(4)<<endl;
2470 cout<<"p="<<p<<" sigma0="<<total->GetParameter(1)<<" sigma1="<<total->GetParameter(5)<<endl;
2471
2472 total->SetParameter(1,total->GetParameter(1));
2473 total->SetParameter(2,total->GetParameter(2));
2474 total->SetParameter(4,total->GetParameter(4));
2475 total->SetParameter(5,total->GetParameter(5));
2476
2477 hproy_pionV0s->Fit(total,"R");
2478
2479
2480
2481 if(p>4.0){
2482
2483
2484 //return;
2485
2486 graphResV0s->SetPoint(nERes,total->GetParameter(4),total->GetParameter(5)/total->GetParameter(4));
2487 graphResV0s->SetPointError(nERes,total->GetParError(4),total->GetParError(5)/total->GetParameter(4));
2488 nERes++;
2489
2490 //graphBB->SetPoint(nBB,p/0.14, total->GetParameter(4));
2491 //graphBB->SetPointError(nBB,0, total->GetParError(4));
2492 //nBB++;
2493
2494
2495 graphBBpiv0s->SetPoint(nBBpiv0s,meanp/0.14, total->GetParameter(4));
2496 graphBBpiv0s->SetPointError(nBBpiv0s,meanpE,total->GetParError(4));
2497 nBBpiv0s++;
2498
2499 l1->Add(hproy_pionV0s);
2500
2501 }
2502
2503
2504
2505 }
2506
2507 }
2508
2509
2510 graphResV0s->SetTitle(";#LT dE/dx #GT;#sigma/#LT dE/dx #GT");
2511 graphResV0s->SetName("gRes_PionsV0sArmenteros");
2512 l1->Add(graphResV0s);
2513
2514 graphBBpiv0s->SetTitle(";#beta*#gamma; #LT dE/dx #GT");
2515 graphBBpiv0s->SetName("BB_piV0s");
2516 graphBBpiv0s->SetMarkerColor(2);
2517 graphBBpiv0s->SetLineColor(2);
2518 graphBBpiv0s->SetMarkerStyle(24);
2519 l1->Add(graphBBpiv0s);
2520
2521 //protons
2522
2523
2524
2525
2526 TGraphErrors* graphPRes = new TGraphErrors();
2527 nERes = 0;
2528 TGraphErrors* graphdEdxvsP = new TGraphErrors();
2529 TGraphErrors* graphSigmavsP = new TGraphErrors();
2530 graphdEdxvsP->SetName("graphdEdxvsP_p");
2531 graphSigmavsP->SetName("graphSigmavsP_p");
2532
2533 nbins=hPi2D->GetXaxis()->GetNbins();
2534 for(Int_t bin=1;bin<=nbins;++bin){
2535
2536 TH1D *hproy_pion=0;
2537 TH1D *hproy_proton=0;
2538
2539
2540 Double_t p=hPi2D->GetXaxis()->GetBinCenter(bin);
2541 Double_t lowp = hP2D->GetXaxis()->GetBinLowEdge(bin);
2542 Double_t widthp = hP2D->GetXaxis()->GetBinWidth(bin);
2543 hP2Dp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
2544 Double_t meanp= hP2Dp->GetMean();
2545 Double_t meanpE= hP2Dp->GetMeanError();
2546
2547
2548 if(p<2.0||p>5.0)
2549 continue;
2550
2551
2552 if(p>2){
2553
2554 cout<<"p="<<hPi2D->GetXaxis()->GetBinCenter(bin)<<endl;
2555
2556 hproy_pion=(TH1D *)hPi2D->ProjectionY(Form("hproy_pion_%d",bin),bin,bin);
2557 hproy_pion->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2558 hPi2D->GetXaxis()->GetBinLowEdge(bin),
2559 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2560
2561 hproy_pion->Fit(pionsv0s,"R");
2562
2563 hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
2564 hproy_proton->SetTitle(Form("protonsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2565 hPi2D->GetXaxis()->GetBinLowEdge(bin),
2566 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2567
2568
2569
2570 hproy_proton->Fit(protonsv0s,"R");
2571
2572 hproy_proton->Fit(protonsv0s,"R","",protonsv0s->GetParameter(1)-2.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+0.8*protonsv0s->GetParameter(2));
2573
2574
2575 Float_t reducedCh2=protonsv0s->GetChisquare()/protonsv0s->GetNDF();
2576 //if(reducedCh2<3)
2577 if(reducedCh2>3){
2578 total->SetParameter(1,protonsv0s->GetParameter(1));
2579 total->SetParameter(2,protonsv0s->GetParameter(2));
2580 total->SetParameter(4,pionsv0s->GetParameter(1));
2581 total->SetParameter(5,pionsv0s->GetParameter(2));
2582
2583 hproy_proton->Fit(total,"R");
2584 }
2585 //if(total->GetParameter(1)<=0)continue;
2586
2587 if(p>2&&p<7){
2588 if(reducedCh2<5){
2589
2590
2591
2592 if(p>3.0&&p<5){
2593 graphPRes->SetPoint(nERes,protonsv0s->GetParameter(1),protonsv0s->GetParameter(2)/protonsv0s->GetParameter(1));
2594 graphPRes->SetPointError(nERes,protonsv0s->GetParError(1),protonsv0s->GetParError(2)/protonsv0s->GetParameter(1));
2595 nERes++;
2596 }
2597 graphBB->SetPoint(nBB,meanp/0.943, protonsv0s->GetParameter(1));
2598 graphBB->SetPointError(nBB,meanpE,protonsv0s->GetParError(1));
2599
2600 graphBBpv0s->SetPoint(nBBpv0s,meanp/0.943, protonsv0s->GetParameter(1));
2601 graphBBpv0s->SetPointError(nBBpv0s,meanpE,protonsv0s->GetParError(1));
2602 nBBpv0s++;
2603
2604 graphdEdxvsP->SetPoint(nBB,p, protonsv0s->GetParameter(1));
2605 graphdEdxvsP->SetPointError(nBB,0.001, protonsv0s->GetParError(1));
2606 graphSigmavsP->SetPoint(nBB,p, protonsv0s->GetParameter(2));
2607 graphSigmavsP->SetPointError(nBB,0.0001, protonsv0s->GetParError(2));
2608 nBB++;
2609
2610 }
2611 else{
2612
2613 if(p>3.0&&p<5){
2614 graphPRes->SetPoint(nERes,total->GetParameter(1),total->GetParameter(2)/total->GetParameter(1));
2615 graphPRes->SetPointError(nERes,total->GetParError(1),total->GetParError(2)/total->GetParameter(1));
2616 nERes++;
2617 }
2618
2619 //graphBB->SetPoint(nBB,p/0.943, total->GetParameter(1));
2620 //graphBB->SetPointError(nBB,0,total->GetParError(1));
2621
2622 graphBBpv0s->SetPoint(nBBpv0s,meanp/0.943, protonsv0s->GetParameter(1));
2623 graphBBpv0s->SetPointError(nBBpv0s,meanpE,protonsv0s->GetParError(1));
2624 nBBpv0s++;
2625
2626
2627 graphdEdxvsP->SetPoint(nBB,p, total->GetParameter(1));
2628 graphdEdxvsP->SetPointError(nBB,0.001, total->GetParError(1));
2629 graphSigmavsP->SetPoint(nBB,p, total->GetParameter(2));
2630 graphSigmavsP->SetPointError(nBB,0.0001, total->GetParError(2));
2631
2632 graphBB->SetPoint(nBB,meanp/0.943, protonsv0s->GetParameter(1));
2633 graphBB->SetPointError(nBB,meanpE,protonsv0s->GetParError(1));
2634
2635 nBB++;
2636 }
2637
2638
2639 }
2640 /*else{
2641
2642 graphPRes->SetPoint(nERes,protonsv0s->GetParameter(1),protonsv0s->GetParameter(2)/protonsv0s->GetParameter(1));
2643 graphPRes->SetPointError(nERes,protonsv0s->GetParError(1),protonsv0s->GetParError(2)/protonsv0s->GetParameter(1));
2644
2645
2646 nERes++;
2647
2648 }*/
2649
2650 l1->Add(hproy_proton);
2651
2652
2653 }
2654
2655
2656
2657 }
2658
2659 graphPRes->SetTitle(";#LT dE/dx #GT;#sigma/#LT dE/dx #GT");
2660 graphPRes->SetName("gRes_Protons");
2661
2662 l1->Add(graphPRes);
2663
2664 graphBBpv0s->SetTitle(";#beta*#gamma; #LT dE/dx #GT");
2665 graphBBpv0s->SetName("BB_pv0s");
2666 graphBBpv0s->SetMarkerColor(4);
2667 graphBBpv0s->SetLineColor(4);
2668 graphBBpv0s->SetMarkerStyle(24);
2669 l1->Add(graphBBpv0s);
2670
2671
2672
2673 graphBB->SetTitle(";#beta*#gamma; #LT dE/dx #GT");
2674 graphBB->SetName("gBB");
2675
2676
2677 Double_t dedxPar[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
2678
2679 dedxPar[0] = 0;
2680 dedxPar[1] = 31.384;
2681 dedxPar[2] = 9.5253;
2682 dedxPar[3] = 7.3542;
2683 dedxPar[4] = 1.5;
2684
2685 TF1 *dedxFunc = new TF1("dedxFunc", FitFunc, 2.5, 60, 6);
2686 dedxFunc->SetParameters(dedxPar);
2687 // Fit betagamma (option+40)
2688 dedxFunc->SetParameter(0, 6+40);
2689 dedxFunc->FixParameter(0, 6+40);
2690 //parameters for MB
2691 dedxFunc->SetParameter(1, 33.7538);
2692 dedxFunc->SetParameter(2, 8.05926);
2693 dedxFunc->SetParameter(3, 1.69954);
2694 dedxFunc->SetParameter(4, 1.37793);
2695 dedxFunc->SetParameter(5, 78.0465);
2696
2697
2698 cout << "!!!!!!!!!!!!!!!!! Setting new plateau=" << fPlateau->GetParameter(1) << endl;
2699 dedxFunc->SetParameter(5, fPlateau->GetParameter(1));
2700 dedxFunc->FixParameter(5, fPlateau->GetParameter(1));
2701
2702
2703 graphBB->Fit(dedxFunc, "0R", "SAME");
2704 graphBB->Fit(dedxFunc, "0R", "SAME");
2705 graphBB->Fit(dedxFunc, "0R", "SAME");
2706
2707 graphBB->Draw();
2708
2709 for(Int_t i=0;i<dedxFunc->GetNpar();++i)
2710 dedxFunc->SetParameter(i,dedxFunc->GetParameter(i));
2711
2712 l1->Add(dedxFunc);
2713 l1->Add(graphBB);
2714
2715
2716
2717 //Fit to get resolution
2718 TGraphErrors *g1=0;
2719 g1=new TGraphErrors();
2720 g1->SetName("res_allpions");
2721
2722 Int_t npionts=0;
2723
2724
2725 for(Int_t ipoint = 0; ipoint < graphPRes->GetN(); ipoint++){
2726 g1->SetPoint(npionts,graphPRes->GetX()[ipoint],graphPRes->GetY()[ipoint]);
2727 g1->SetPointError(npionts,graphPRes->GetEX()[ipoint],graphPRes->GetEY()[ipoint]);
2728 npionts++;
2729 }
2730
2731 // pions TOF
2732
2733 for(Int_t ipoint = 0; ipoint < graphRes->GetN(); ipoint++){
2734
2735 g1->SetPoint(npionts,graphRes->GetX()[ipoint],graphRes->GetY()[ipoint]);
2736 g1->SetPointError(npionts,graphRes->GetEX()[ipoint],graphRes->GetEY()[ipoint]);
2737 npionts++;
2738 }
2739 //electrons
2740 g1->SetPoint(npionts,fPlateau->GetParameter(1),fPlateau->GetParameter(2)/fPlateau->GetParameter(1));
2741 g1->SetPointError(npionts,fPlateau->GetParError(1),fPlateau->GetParError(2)/fPlateau->GetParameter(1));
2742
2743
2744 Float_t sigmaPar[7];
2745 sigmaPar[0]=12;
2746 sigmaPar[1]=7.96380e-03;
2747 sigmaPar[2]=8.09869e-05;
2748 sigmaPar[3]=6.29707e-06;
2749 sigmaPar[4]=-3.27295e-08;
2750 sigmaPar[5]=-1.20200e+03;
2751 sigmaPar[6]=-3.97089e+01;
2752
2753
2754
2755 TF1 *sigmarrrrr=new TF1("sigmaRes",SigmaFunc,47.5, 85,4);
2756 sigmarrrrr->FixParameter(0,sigmaPar[0]);
2757 g1->Fit(sigmarrrrr, "0R", "SAME");
2758 g1->Fit(sigmarrrrr, "0R", "SAME");
2759 g1->Fit(sigmarrrrr, "0R", "SAME");
2760
2761 TF1 *fdedxvspP=new TF1("fdedxvspP","pol0",3.0,5.5);
2762 graphdEdxvsP->Fit(fdedxvspP, "R", "SAME");
2763 TF1 *fsigmavspP=new TF1("fsigmavspP","pol0",3.0,5.5);
2764 graphSigmavsP->Fit(fsigmavspP, "R", "SAME");
2765
2766
2767
2768 for(Int_t i=0;i<6;++i)
2769 sigmarrrrr->SetParameter(i,sigmarrrrr->GetParameter(i));
2770
2771 l1->Add(graphdEdxvsP);
2772 l1->Add(graphSigmavsP);
2773
2774 l1->Add(fdedxvspP);
2775 l1->Add(fsigmavspP);
2776
2777 l1->Add(sigmarrrrr);
2778 l1->Add(g1);
2779
2780 TFile *fout=new TFile(Form("%s/hres_0_5_%s.root",outDir,endNamePos[index_eta]),"recreate");
2781 fout->cd();
2782 l1->Write();
2783 fout->Close();
2784
2785
2786
2787
2788}
2789
2790
2791
2792void MakeFitsExternalData(const Char_t * inFile, const Char_t * outDir){
2793
2794
2795 TFile *fIn = TFile::Open(inFile);
2796 TList * list = (TList *)fIn->Get("outputdedx");
2797 //plateau, electrons 0.4<p<0.6 GeV/c
2798 TH2D *hEVsEta = (TH2D *)list->FindObject("hPlateauVsEta");
2799 hEVsEta->Sumw2();
2800 //four th1, projections
2801 TH1D *hdEdxE[4]={0,0,0,0};// |eta|
2802 TH1D *hdEdxEpos[4]={0,0,0,0};// eta>0
2803 TF1 *fGaussE[4]={0,0,0,0};
2804 Int_t index_max = 1;
2805 for(Int_t i=0; i<4; ++i){
2806
2807 Int_t index_min_negeta = i+index_max;
2808 Int_t index_max_negeta = i+index_max+1;
2809 Int_t index_min_poseta = 17-(i+index_max+1);
2810 Int_t index_max_poseta = 17-(i+index_max);
2811
2812 hdEdxE[i] = (TH1D *)hEVsEta->ProjectionY(Form("hdEdxNegEta%d",i),index_min_negeta,index_max_negeta);
2813 hdEdxEpos[i] = (TH1D *)hEVsEta->ProjectionY(Form("hdEdxPosEta%d",i),index_min_poseta,index_max_poseta);
2814 hdEdxE[i]->Add(hdEdxEpos[i]);
2815
2816 //Fitting the signal
2817 fGaussE[i]=new TF1(Form("fGauss%d",i),"gaus");
2818 hdEdxE[i]->Fit(fGaussE[i],"0","");
2819 hdEdxE[i]->Fit(fGaussE[i],"","", fGaussE[i]->GetParameter(1)-1.5*fGaussE[i]->GetParameter(2), fGaussE[i]->GetParameter(1)+1.5*fGaussE[i]->GetParameter(2));
2820
2821 index_max++;
2822
2823 }
2824
2825 CreateDir(outDir);
2826
2827 TFile *fout = 0;
2828
2829 if(outDir){
2830
2831 fout = new TFile(Form("%s/PrimaryElectrons.root",outDir),"recreate");
2832 fout->cd();
2833 for(Int_t i=0; i<4; ++i){
2834 hdEdxE[i]->Write();
2835 fGaussE[i]->Write();
2836 }
2837 fout->Close();
2838
2839
2840 }
2841
2842
2843
2844}
2845//__________________________________________
2846void PlotQA(const Char_t * inFile){
2847
2848 //gStyle->SetPalette(1);
2849 gStyle->SetCanvasColor(10);
2850 gStyle->SetFrameFillColor(10);
2851 gStyle->SetOptStat(0);
2852 gStyle->SetOptTitle(0);
2853 gStyle->SetOptFit(0);
2854
2855 const Char_t* ending[8] = { "86", "64", "42", "20", "02", "24", "46", "68" };
2856 const Char_t* LatexEta[8] = { "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0",
2857 "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8" };
2858 TFile *fIn = TFile::Open(inFile);
2859 if(!fIn)
2860 return;
2861
2862 TList * list = (TList *)fIn->Get("outputdedx");
2863
2864
2865 TH2D *hMIP[8];
2866 TH2D *hPlateau[8];
2867 TProfile *pMIP[8];
2868 TProfile *pPlateau[8];
2869
2870 TH2D *hMIPVsNch[8];
2871 TProfile *pMIPVsNch[8];
2872
2873
2874 //Secondary pions
2875 TH2D *hPiV0s2D[8];
2876 TH1D *hPion3[8];
2877 TH1D *hPion2GeV = new TH1D("hPion2GeV","hPion2GeV",8,-0.8,0.8);
2878
2879
2880 for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
2881 hMIP[index_eta] = (TH2D *)list->FindObject(Form("hMIPVsPhi%s",ending[index_eta]));
2882 hPlateau[index_eta] = (TH2D *)list->FindObject(Form("hPlateauVsPhi%s",ending[index_eta]));
2883 pMIP[index_eta] = (TProfile *)list->FindObject(Form("pMIPVsPhi%s",ending[index_eta]));
2884 pPlateau[index_eta] = (TProfile *)list->FindObject(Form("pPlateauVsPhi%s",ending[index_eta]));
2885 hPiV0s2D[index_eta] = (TH2D *)list->FindObject(Form("histPiV0%s",ending[index_eta]));
2886 hPion3[index_eta] = (TH1D *)hPiV0s2D[index_eta]->ProjectionY(Form("Pion3%s",ending[index_eta]),16,16);
2887
2888 hMIPVsNch[index_eta] = (TH2D *)list->FindObject(Form("hMIPVsNch%s",ending[index_eta]));
2889 pMIPVsNch[index_eta] = (TProfile *)list->FindObject(Form("pMIPVsNch%s",ending[index_eta]));
2890
2891
2892
2893 TF1 *fgaus = 0;
2894 fgaus = new TF1("fgaus","gausn",40,80);
2895 hPion3[index_eta]->Fit(fgaus,"0","",40,80);
2896
2897 hPion2GeV->SetBinContent(index_eta+1,fgaus->GetParameter(1));
2898 hPion2GeV->SetBinError(index_eta+1,fgaus->GetParError(1));
2899 }
2900
2901
2902
2903 TH2D *hMIPVsEta = (TH2D *)list->FindObject("hMIPVsEta");
2904 TProfile *pMIPVsEta = (TProfile *)list->FindObject("pMIPVsEta");
2905 TProfile *pMIPVsEtaV0s = (TProfile *)list->FindObject("pMIPVsEtaV0s");
2906 TProfile *pPlateauVsEta = (TProfile *)list->FindObject("pPlateauVsEta");
2907 TH2D *hPlateauVsEta = (TH2D *)list->FindObject("hPlateauVsEta");
2908
2909
2910
2911 const Int_t nPadX = 4;
2912 const Int_t nPadY = 2;
2913 Float_t noMargin = 0.000;
2914 Float_t topMargin = 0.01;
2915 Float_t botMargin = 0.1;
2916 Float_t leftMargin = 0.05;
2917 Float_t rightMargin = 0.01;
2918 Float_t width = (1-rightMargin-leftMargin)/nPadX;
2919 Float_t height = (1-botMargin-topMargin)/nPadY;
2920 Float_t shift = 0.05;
2921
2922
2923 TCanvas* cMIP = new TCanvas("cMIP2", "Raa Pions", 900, 500);
2924 TPad* padMIP[nPadX*nPadY];
2925 //Float_t shift = 0.1;
2926 //Float_t shift = 0.05;
2927 const char* yTitleMIP = "d#it{E}/d#it{x}_{MIP}";
2928 const char* xTitleMIP = "#phi (rad)";
2929 cMIP->cd();
2930
2931
2932 TLatex* latex = new TLatex();
2933 latex->SetNDC();
2934 latex->SetTextAlign(22);
2935 latex->SetTextAngle(0);
2936 latex->SetTextFont(62);
2937
2938 latex->DrawLatex(0.53,0.04,xTitleMIP);
2939 latex->SetTextAngle(90);
2940 latex->SetTextSize(0.055);
2941 latex->DrawLatex(0.025,0.5,yTitleMIP);
2942 latex->SetTextSize(0.05);
2943
2944
2945 for (Int_t i = 1; i < 9; i++) {
2946
2947 Int_t iX = (i-1)%nPadX;
2948 Int_t iY = (i-1)/nPadX;
2949 Float_t x1 = leftMargin + iX*width;
2950 if(iX==2)
2951 x1-=0.01;
2952 else if(iX==1)
2953 x1+=0.01;
2954 Float_t x2 = leftMargin + (iX + 1)*width;
2955 if(iX==0)
2956 x2+=0.01;
2957 else if(iX==1)
2958 x2-=0.01;
2959 Float_t y1 = 1 - topMargin - (iY +1)*height;
2960 Float_t y2 = 1 - topMargin - iY*height;
2961 padMIP[i-1] = new TPad(Form("padRaa%d",i),"", x1, y1, x2, y2, 0, 0, 0);
2962 Float_t mBot = noMargin;
2963 Float_t mTop = noMargin;
2964 Float_t mLeft = noMargin;
2965 Float_t mRight = noMargin;
2966 if(iY==0) mTop = shift;
2967 if(iY==nPadY-1) mBot = shift;
2968 if(iX==0) mLeft = 0.08;
2969 if(iX==nPadX-1) mRight = 0.08;
2970 padMIP[i-1]->SetBottomMargin(mBot);
2971 padMIP[i-1]->SetTopMargin(mTop);
2972 padMIP[i-1]->SetLeftMargin(mLeft);
2973 padMIP[i-1]->SetRightMargin(mRight);
2974
2975 cMIP->cd();
2976
2977 padMIP[i-1]->Draw();
2978
2979 }
2980
2981 latex->SetTextAngle(0);
2982 latex->SetTextSize(0.085);
2983
2984 Int_t index_dec = 7;
2985
2986 for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
2987 padMIP[index_eta]->cd();
2988 padMIP[index_eta]->SetTickx(1);
2989 padMIP[index_eta]->SetTicky(1);
2990 padMIP[index_eta]->SetLogy(0);
2991 padMIP[index_eta]->SetLogx(0);
2992 hMIP[index_eta]->GetXaxis()->SetLabelSize(0.06);
2993 hMIP[index_eta]->GetYaxis()->SetLabelSize(0.06);
2994 hMIP[index_eta]->GetXaxis()->SetTitleSize(0);
2995 hMIP[index_eta]->GetYaxis()->SetTitleSize(0);
2996
2997 if(index_eta<4){
2998 hMIP[index_eta]->Draw("colz");
2999 pMIP[index_eta]->Draw("samep");
3000 latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
3001 }
3002 else{
3003 hMIP[index_dec]->Draw("colz");
3004 pMIP[index_dec]->Draw("samep");
3005 latex->DrawLatex(0.5,0.1,LatexEta[index_dec]);
3006 index_dec--;
3007 }
3008
3009 }
3010
3011
3012
3013 cMIP->SaveAs("MIPvsPhi.png");
3014
3015 //Plateau vs phi, different eta intervals
3016 TCanvas* cPlateau = new TCanvas("cPlateau2", "Raa Pions", 900, 500);
3017 TPad* padPlateau[nPadX*nPadY];
3018
3019 const char* yTitlePlateau = "d#it{E}/d#it{x}_{Plateau}";
3020 const char* xTitlePlateau = "#phi (rad)";
3021 cPlateau->cd();
3022
3023
3024 //TLatex* latex = new TLatex();
3025 latex->SetNDC();
3026 latex->SetTextAlign(22);
3027 latex->SetTextAngle(0);
3028 latex->SetTextFont(62);
3029
3030 latex->DrawLatex(0.53,0.04,xTitlePlateau);
3031 latex->SetTextAngle(90);
3032 latex->SetTextSize(0.055);
3033 latex->DrawLatex(0.025,0.5,yTitlePlateau);
3034 latex->SetTextSize(0.05);
3035
3036
3037 for (Int_t i = 1; i < 9; i++) {
3038
3039 Int_t iX = (i-1)%nPadX;
3040 Int_t iY = (i-1)/nPadX;
3041 Float_t x1 = leftMargin + iX*width;
3042 if(iX==2)
3043 x1-=0.01;
3044 else if(iX==1)
3045 x1+=0.01;
3046 Float_t x2 = leftMargin + (iX + 1)*width;
3047 if(iX==0)
3048 x2+=0.01;
3049 else if(iX==1)
3050 x2-=0.01;
3051 Float_t y1 = 1 - topMargin - (iY +1)*height;
3052 Float_t y2 = 1 - topMargin - iY*height;
3053 padPlateau[i-1] = new TPad(Form("padPlateau%d",i),"", x1, y1, x2, y2, 0, 0, 0);
3054 Float_t mBot = noMargin;
3055 Float_t mTop = noMargin;
3056 Float_t mLeft = noMargin;
3057 Float_t mRight = noMargin;
3058 if(iY==0) mTop = shift;
3059 if(iY==nPadY-1) mBot = shift;
3060 if(iX==0) mLeft = 0.08;
3061 if(iX==nPadX-1) mRight = 0.08;
3062 padPlateau[i-1]->SetBottomMargin(mBot);
3063 padPlateau[i-1]->SetTopMargin(mTop);
3064 padPlateau[i-1]->SetLeftMargin(mLeft);
3065 padPlateau[i-1]->SetRightMargin(mRight);
3066
3067 cPlateau->cd();
3068
3069 padPlateau[i-1]->Draw();
3070
3071 }
3072
3073 latex->SetTextAngle(0);
3074 latex->SetTextSize(0.085);
3075
3076 index_dec = 7;
3077
3078 for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
3079 padPlateau[index_eta]->cd();
3080 padPlateau[index_eta]->SetTickx(1);
3081 padPlateau[index_eta]->SetTicky(1);
3082 padPlateau[index_eta]->SetLogy(0);
3083 padPlateau[index_eta]->SetLogx(0);
3084 hPlateau[index_eta]->GetXaxis()->SetLabelSize(0.06);
3085
3086 hPlateau[index_eta]->GetYaxis()->SetLabelSize(0.06);
3087 hPlateau[index_eta]->GetXaxis()->SetTitleSize(0);
3088 hPlateau[index_eta]->GetYaxis()->SetTitleSize(0);
3089
3090 if(index_eta<4){
3091 hPlateau[index_eta]->Draw("colz");
3092 pPlateau[index_eta]->Draw("samep");
3093 latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
3094 }
3095 else{
3096 hPlateau[index_dec]->Draw("colz");
3097 pPlateau[index_dec]->Draw("samep");
3098 latex->DrawLatex(0.5,0.1,LatexEta[index_dec]);
3099 index_dec--;
3100 }
3101
3102 }
3103
3104
3105
3106 cPlateau->SaveAs("PlateauVsPhi.png");
3107
3108
3109 //MIP vs Nch
3110 TCanvas* cMIPVsNch = new TCanvas("cMIPVsNch2", "Raa Pions", 900, 500);
3111 TPad* padMIPVsNch[nPadX*nPadY];
3112
3113 const char* yTitleMIPVsNch = "d#it{E}/d#it{x}_{MIP}";
3114 const char* xTitleMIPVsNch = "TPC-track multiplicity";
3115 cMIPVsNch->cd();
3116
3117
3118 latex->SetNDC();
3119 latex->SetTextAlign(22);
3120 latex->SetTextAngle(0);
3121 latex->SetTextFont(62);
3122
3123 latex->DrawLatex(0.53,0.04,xTitleMIPVsNch);
3124 latex->SetTextAngle(90);
3125 latex->SetTextSize(0.055);
3126 latex->DrawLatex(0.025,0.5,yTitleMIPVsNch);
3127 latex->SetTextSize(0.05);
3128
3129
3130 for (Int_t i = 1; i < 9; i++) {
3131
3132 Int_t iX = (i-1)%nPadX;
3133 Int_t iY = (i-1)/nPadX;
3134 Float_t x1 = leftMargin + iX*width;
3135 if(iX==2)
3136 x1-=0.01;
3137 else if(iX==1)
3138 x1+=0.01;
3139 Float_t x2 = leftMargin + (iX + 1)*width;
3140 if(iX==0)
3141 x2+=0.01;
3142 else if(iX==1)
3143 x2-=0.01;
3144 Float_t y1 = 1 - topMargin - (iY +1)*height;
3145 Float_t y2 = 1 - topMargin - iY*height;
3146 padMIPVsNch[i-1] = new TPad(Form("padMIPVsNch%d",i),"", x1, y1, x2, y2, 0, 0, 0);
3147 Float_t mBot = noMargin;
3148 Float_t mTop = noMargin;
3149 Float_t mLeft = noMargin;
3150 Float_t mRight = noMargin;
3151 if(iY==0) mTop = shift;
3152 if(iY==nPadY-1) mBot = shift;
3153 if(iX==0) mLeft = 0.08;
3154 if(iX==nPadX-1) mRight = 0.08;
3155 padMIPVsNch[i-1]->SetBottomMargin(mBot);
3156 padMIPVsNch[i-1]->SetTopMargin(mTop);
3157 padMIPVsNch[i-1]->SetLeftMargin(mLeft);
3158 padMIPVsNch[i-1]->SetRightMargin(mRight);
3159
3160 if(i>=5&&i<9)
3161 padMIPVsNch[i-1]->SetBottomMargin(mBot+0.05);
3162
3163 cMIPVsNch->cd();
3164
3165 padMIPVsNch[i-1]->Draw();
3166
3167 }
3168
3169 latex->SetTextAngle(0);
3170 latex->SetTextSize(0.085);
3171
3172 index_dec = 7;
3173
3174 for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
3175 padMIPVsNch[index_eta]->cd();
3176
3177 padMIPVsNch[index_eta]->SetTickx(1);
3178 padMIPVsNch[index_eta]->SetTicky(1);
3179 padMIPVsNch[index_eta]->SetLogy(0);
3180 padMIPVsNch[index_eta]->SetLogx(1);
3181 hMIPVsNch[index_eta]->GetXaxis()->SetLabelSize(0.08);
3182
3183 hMIPVsNch[index_eta]->GetXaxis()->SetRangeUser(10,2000);
3184 hMIPVsNch[index_eta]->GetYaxis()->SetRangeUser(43.5,54.5);
3185 hMIPVsNch[index_eta]->GetYaxis()->SetLabelSize(0.06);
3186 hMIPVsNch[index_eta]->GetXaxis()->SetTitleSize(0);
3187 hMIPVsNch[index_eta]->GetYaxis()->SetTitleSize(0);
3188
3189 if(index_eta<4){
3190 hMIPVsNch[index_eta]->Draw("colz");
3191 pMIPVsNch[index_eta]->Draw("samep");
3192 latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
3193 }
3194 else{
3195 hMIPVsNch[index_dec]->Draw("colz");
3196 pMIPVsNch[index_dec]->Draw("samep");
3197 latex->DrawLatex(0.5,0.2,LatexEta[index_dec]);
3198 index_dec--;
3199 }
3200
3201 }
3202
3203
3204
3205
3206
3207
3208
3209
3210 cMIPVsNch->SaveAs("MIPvsNch.png");
3211
3212
3213
3214
3215
3216
3217
3218
3219
3220 //eta dependence, scaling
3221 const Int_t nPadXS = 3;
3222 const Int_t nPadYS = 1;
3223 Float_t noMarginS = 0.0;
3224 Float_t topMarginS = 0.01;
3225 Float_t botMarginS = 0.1;
3226 Float_t leftMarginS = 0.05;
3227 Float_t rightMarginS = 0.01;
3228 Float_t widthS = (1-rightMarginS-leftMarginS)/nPadXS;
3229 Float_t heightS = (1-botMarginS-topMarginS)/nPadYS;
3230 Float_t shiftS = 0.05;
3231
3232
3233 TCanvas* cMIPEta = new TCanvas("cMIPEta", "MIP eta", 900, 500);
3234 TPad* padMIPEta[nPadXS*nPadYS];
3235 //Float_t shift = 0.1;
3236 //Float_t shift = 0.05;
3237 const char* yTitleMIPS = "d#it{E}/d#it{x}";
3238 const char* xTitleMIPS = "#eta";
3239 cMIPEta->cd();
3240
3241
3242 //TLatexlatex = new TLatex();
3243 latex->SetNDC();
3244 latex->SetTextAlign(22);
3245 latex->SetTextAngle(0);
3246 latex->SetTextFont(42);
3247
3248 latex->DrawLatex(0.53,0.04,xTitleMIPS);
3249 latex->SetTextAngle(90);
3250 latex->SetTextSize(0.055);
3251 latex->DrawLatex(0.025,0.5,yTitleMIPS);
3252 latex->SetTextSize(0.05);
3253
3254
3255 for (Int_t i = 1; i < 4; i++) {
3256
3257 Int_t iX = (i-1)%nPadXS;
3258 Int_t iY = (i-1)/nPadXS;
3259 Float_t x1 = leftMarginS + iX*widthS;
3260 if(iX==2)
3261 x1-=0.01;
3262 else if(iX==1)
3263 x1+=0.01;
3264 Float_t x2 = leftMarginS + (iX + 1)*widthS;
3265 if(iX==0)
3266 x2+=0.01;
3267 else if(iX==1)
3268 x2-=0.01;
3269 Float_t y1 = 1 - topMarginS - (iY +1)*heightS;
3270 Float_t y2 = 1 - topMarginS - iY*heightS;
3271 padMIPEta[i-1] = new TPad(Form("padMIPEta%d",i),"", x1, y1, x2, y2, 0, 0, 0);
3272 Float_t mBot = noMarginS;
3273 Float_t mTop = noMarginS;
3274 Float_t mLeft = noMarginS;
3275 Float_t mRight = noMarginS;
3276 if(iY==0) mTop = shiftS;
3277 if(iY==nPadYS-1) mBot = shiftS;
3278 if(iX==0) mLeft = 0.08;
3279 if(iX==nPadXS-1) mRight = 0.08;
3280 padMIPEta[i-1]->SetBottomMargin(mBot);
3281 padMIPEta[i-1]->SetTopMargin(mTop);
3282 padMIPEta[i-1]->SetLeftMargin(mLeft);
3283 padMIPEta[i-1]->SetRightMargin(mRight);
3284
3285 cMIPEta->cd();
3286
3287 padMIPEta[i-1]->Draw();
3288
3289 }
3290
3291 latex->SetTextAngle(0);
3292 latex->SetTextSize(0.085);
3293
3294
3295 TF1 *fetapos = new TF1("fetapos","pol3",0.0,0.8);
3296 TF1 *fetaneg = new TF1("fetaneg","pol3",-0.8,0.0);
3297
3298 pMIPVsEta->Fit(fetaneg,"0","",-0.8,0.0);
3299 pMIPVsEta->Fit(fetapos,"0","",0.0,0.8);
3300 TLegend *leg1=new TLegend(0.15,0.5,0.4,0.7);
3301
3302 fetapos->SetLineColor(1);
3303 fetapos->SetLineStyle(1);
3304 fetapos->SetLineWidth(2);
3305
3306 fetaneg->SetLineColor(1);
3307 fetaneg->SetLineStyle(1);
3308 fetaneg->SetLineWidth(2);
3309
3310 TH1D *hframe20 = new TH1D("hframe20","hframe20",8,-0.8,0.8);
3311
3312 TH1D *hMIPPiPrim = new TH1D("hMIPPiPrim","",16,-0.8,0.8);
3313 for(Int_t bin =1; bin <= hMIPPiPrim->GetNbinsX(); ++bin){
3314 Double_t eta = pMIPVsEta->GetBinCenter(bin);
3315 Double_t dedx = pMIPVsEta->GetBinContent(bin);
3316 Double_t e_dedx = pMIPVsEta->GetBinError(bin);
3317 if(eta<0){
3318 dedx *= 50/fetaneg->Eval(eta);
3319 e_dedx *= 50/fetaneg->Eval(eta);
3320 }else{
3321 dedx *= 50/fetapos->Eval(eta);
3322 e_dedx *= 50/fetapos->Eval(eta);
3323 }
3324 hMIPPiPrim->SetBinContent(bin,dedx);
3325 hMIPPiPrim->SetBinError(bin,e_dedx);
3326
3327 }
3328
3329 TH1D *hPlateauPiPrim = new TH1D("hPlateauPiPrim","",16,-0.8,0.8);
3330 for(Int_t bin =1; bin <= hPlateauPiPrim->GetNbinsX(); ++bin){
3331 Double_t eta = pPlateauVsEta->GetBinCenter(bin);
3332 Double_t dedx = pPlateauVsEta->GetBinContent(bin);
3333 Double_t e_dedx = pPlateauVsEta->GetBinError(bin);
3334 if(eta<0){
3335 dedx *= 50/fetaneg->Eval(eta);
3336 e_dedx *= 50/fetaneg->Eval(eta);
3337 }else{
3338 dedx *= 50/fetapos->Eval(eta);
3339 e_dedx *= 50/fetapos->Eval(eta);
3340 }
3341 hPlateauPiPrim->SetBinContent(bin,dedx);
3342 hPlateauPiPrim->SetBinError(bin,e_dedx);
3343
3344 }
3345
3346
3347 TH1D *h2GeVPiSec = new TH1D("h2GeVPiSec","",8,-0.8,0.8);
3348 for(Int_t bin =1; bin <= h2GeVPiSec->GetNbinsX(); ++bin){
3349 Double_t eta = hPion2GeV->GetBinCenter(bin);
3350 Double_t dedx = hPion2GeV->GetBinContent(bin);
3351 Double_t e_dedx = hPion2GeV->GetBinError(bin);
3352 if(eta<0){
3353 dedx *= 50/fetaneg->Eval(eta);
3354 e_dedx *= 50/fetaneg->Eval(eta);
3355 }else{
3356 dedx *= 50/fetapos->Eval(eta);
3357 e_dedx *= 50/fetapos->Eval(eta);
3358 }
3359 h2GeVPiSec->SetBinContent(bin,dedx);
3360 h2GeVPiSec->SetBinError(bin,e_dedx);
3361
3362 }
3363
3364
3365
3366
3367
3368
3369
3370
3371
3372
3373 for(Int_t index_eta = 0; index_eta < 3; ++index_eta){
3374 padMIPEta[index_eta]->cd();
3375 padMIPEta[index_eta]->SetTickx(1);
3376 padMIPEta[index_eta]->SetTicky(1);
3377 padMIPEta[index_eta]->SetLogy(0);
3378 padMIPEta[index_eta]->SetLogx(0);
3379
3380 if(index_eta==0){
3381 hframe20->GetYaxis()->SetRangeUser(30,95);
3382 hframe20->GetXaxis()->SetLabelSize(0.06);
3383 hframe20->GetXaxis()->SetLabelOffset(0.004);
3384 hframe20->GetYaxis()->SetLabelSize(0.06);
3385 hframe20->GetXaxis()->SetTitleSize(0);
3386 hframe20->GetYaxis()->SetTitleSize(0);
3387 hframe20->Draw();
3388 hMIPVsEta->Draw("colz same");
3389 hPlateauVsEta->Draw("colz same");
3390
3391 latex->SetTextSize(0.05);
3392 latex->DrawLatex(0.55,0.45,"MIP Pions, 0.4<#it{p}<0.6 GeV/c");
3393 latex->DrawLatex(0.55,0.85,"Electrons, 0.4<#it{p}<0.6 GeV/c");
3394
3395 }
3396
3397
3398 if(index_eta==1){
3399 pMIPVsEtaV0s->GetXaxis()->SetLabelSize(0.06);
3400 pMIPVsEtaV0s->GetYaxis()->SetLabelSize(0.06);
3401 pMIPVsEtaV0s->GetXaxis()->SetLabelOffset(0.004);
3402 pMIPVsEtaV0s->GetXaxis()->SetTitleSize(0);
3403 pMIPVsEtaV0s->GetYaxis()->SetTitleSize(0);
3404 pMIPVsEtaV0s->GetYaxis()->SetRangeUser(30,95);
3405 pMIPVsEtaV0s->SetMarkerStyle(29);
3406 pMIPVsEtaV0s->SetMarkerColor(4);
3407 pMIPVsEtaV0s->SetLineColor(4);
3408 pMIPVsEtaV0s->Draw();
3409
3410 pMIPVsEta->SetMarkerStyle(25);
3411 pMIPVsEta->SetMarkerColor(1);
3412 pMIPVsEta->SetLineColor(1);
3413 pMIPVsEta->Draw("samep");
3414 pPlateauVsEta->SetMarkerStyle(21);
3415 pPlateauVsEta->SetMarkerColor(kBlue);
3416 pPlateauVsEta->SetLineColor(4);
3417 pPlateauVsEta->Draw("samep");
3418
3419 //fetapos->Draw("same");
3420 //fetaneg->Draw("same");
3421
3422
3423 pMIPVsEta->Draw("samep");
3424 hPion2GeV->SetMarkerStyle(20);
3425 hPion2GeV->SetMarkerColor(4);
3426 hPion2GeV->SetLineColor(4);
3427 hPion2GeV->Draw("samep");
3428 pPlateauVsEta->Draw("samep");
3429
3430
3431 leg1->SetTextFont(42);
3432 leg1->SetTextSize(0.045);
3433 leg1->SetLineColor(kWhite);
3434 leg1->SetLineStyle(3);
3435 leg1->SetShadowColor(kWhite);
3436 leg1->SetFillColor(kWhite);
3437 leg1->SetFillStyle(0);
3438 leg1->AddEntry(pMIPVsEta,"Primary pions at MIP","P");
3439 leg1->AddEntry(pMIPVsEtaV0s,"Secondary pions at MIP","P");
3440 leg1->AddEntry(hPion2GeV,"Secondary pions at #beta#gamma ~ 14","P");
3441 leg1->AddEntry(pPlateauVsEta,"Electrons, at #beta#gamma ~ 800","P");
3442 leg1->Draw();
3443
3444
3445
3446 }
3447
3448
3449
3450 if(index_eta==2){
3451 hMIPPiPrim->GetXaxis()->SetLabelSize(0.06);
3452 hMIPPiPrim->GetYaxis()->SetLabelSize(0.06);
3453 hMIPPiPrim->GetXaxis()->SetLabelOffset(0.004);
3454 hMIPPiPrim->GetXaxis()->SetTitleSize(0);
3455 hMIPPiPrim->GetYaxis()->SetTitleSize(0);
3456 hMIPPiPrim->GetYaxis()->SetRangeUser(30,95);
3457 hMIPPiPrim->SetMarkerStyle(29);
3458 hMIPPiPrim->SetMarkerColor(4);
3459 hMIPPiPrim->SetLineColor(4);
3460 hMIPPiPrim->Draw();
3461
3462
3463 hPlateauPiPrim->SetMarkerStyle(21);
3464 hPlateauPiPrim->SetLineColor(4);
3465 hPlateauPiPrim->SetMarkerColor(4);
3466 hPlateauPiPrim->Draw("samep");
3467 h2GeVPiSec->SetMarkerStyle(20);
3468 h2GeVPiSec->SetLineColor(4);
3469 h2GeVPiSec->SetMarkerColor(4);
3470 h2GeVPiSec->Draw("samep");
3471
3472 latex->DrawLatex(0.45,0.85,"#LT d#it{E}/d#it{x} #GT #times 50 / #LT d#it{E}/d#it{x} #GT_{primary #pi at MIP}");
3473
3474 }
3475
3476
3477
3478
3479
3480
3481
3482
3483
3484
3485
3486 }
3487
3488
3489 cMIPEta->SaveAs("MIPvsEta.png");
3490
3491
3492
3493
3494}
3495TH2D * AddTwoSameBinningTH2D(TH2D *hPos, TH2D *hNeg, const Char_t *nameHist){
3496
3497 TH2D *hout = (TH2D *)hPos->Clone(nameHist);
3498 hout->Reset();
3499
3500
3501 for(Int_t binx = 1; binx <= hPos->GetNbinsX(); ++binx){
3502
3503 for(Int_t biny = 1; biny <= hPos->GetNbinsY(); ++biny){
3504
3505 Double_t yield1 = hPos->GetBinContent(binx,biny);
3506 Double_t yield2 = hNeg->GetBinContent(binx,biny);
3507 Double_t e_yield1 = hPos->GetBinError(binx,biny);
3508 Double_t e_yield2 = hNeg->GetBinError(binx,biny);
3509
3510 Double_t yield = yield1 + yield2;
3511 Double_t e_yield = TMath::Sqrt( e_yield1*e_yield1 + e_yield2*e_yield2 );
3512
3513 hout->SetBinContent(binx,biny,yield);
3514 hout->SetBinError(binx,biny,e_yield);
3515
3516 }
3517
3518 }
3519
3520 hout->SetEntries( hPos->GetEntries()+hNeg->GetEntries() );
3521
3522 return hout;
3523
3524}