Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / corrected_fraction / corrected_fraction.C
CommitLineData
4ebdd20e 1#include <TFile.h>
2#include <TH1.h>
3#include <TF1.h>
4#include <TCanvas.h>
5#include <TList.h>
6#include <TLegend.h>
7#include <TStyle.h>
8#include <TMath.h>
9#include <TSystem.h>
10#include <TLatex.h>
11
12#include "my_tools.C"
13
14
15/*
16 Info:
17 * Efficiencies are currently set by hand!!!!!
18 * Systematic errors needs an update also.
19 * Need to accomodate k+p at soem point!
20
21 To run:
22
23 gSystem->AddIncludePath("-I../macros")
24 gSystem->AddIncludePath("-I../lib")
25 gROOT->SetMacroPath(".:../macros")
26 .L my_tools.C++
27 .L corrected_fraction.C+
28
29
30 CorrectedFraction("../ratios_7tevb/fit_yields_results/7tev_b.root", "fraction_7tev_b_pi.root", 1, "PP", 0, 10.0);
31 CorrectedFraction("../ratios_7tevb/fit_yields_results/7tev_b.root", "fraction_7tev_b_k.root", 2, "PP", 0, 10.0);
32 CorrectedFraction("../ratios_7tevb/fit_yields_results/7tev_b.root", "fraction_7tev_b_p.root", 3, "PP", 0, 10.0);
33
34
35 And in the shell one can add files like:
36
37 hadd fraction_7tev_b_all.root fraction_7tev_b_pi.root fraction_7tev_b_k.root fraction_7tev_b_p.root
38
39 */
40
41void CorrectedFraction(const Char_t* inFileName, const Char_t* outFileName, Int_t pid,
42 const Char_t* endName="PP", Int_t centBin=0,
43 Double_t ptMax=20.0, const Char_t* outDir="plots");
44TH1D* GetSystErrorHist(TH1D* hRatio, Int_t centBin, TF1* electronFraction);
45
46
47//________________________________________________________________________________________
48void CorrectedFraction(const Char_t* inFileName, const Char_t* outFileName, Int_t pid,
49 const Char_t* endName, Int_t centBin,
50 Double_t ptMax, const Char_t* outDir)
51{
52 gStyle->SetOptStat(0);
53
54
55 TFile* fileData = FindFileFresh(inFileName);
56
57 if(!fileData)
58 return;
59
60 TF1* fElectronFraction = (TF1*)fileData->Get("fElectronFraction");
61 cout << "Electron fraction found!" << endl;
62 TH1D* hPionRatio = 0;
63 if(pid==1) {
64 hPionRatio = (TH1D*)fileData->Get("hPionRatio");
65 CutHistogram(hPionRatio, 2.0, ptMax);
66 hPionRatio->SetMarkerStyle(24);
67 } else if(pid==2) {
68 hPionRatio = (TH1D*)fileData->Get("hKaonRatio");
69 CutHistogram(hPionRatio, 3.0, ptMax);
70 hPionRatio->SetMarkerStyle(25);
71 } else if(pid==3) {
72 hPionRatio = (TH1D*)fileData->Get("hProtonRatio");
73 CutHistogram(hPionRatio, 3.0, ptMax);
74 hPionRatio->SetMarkerStyle(25);
75 }
76
77 // Global variable
78 TH1D* hSystFraction = 0;
79 if(pid==1)
80 hSystFraction = GetSystErrorHist(hPionRatio, centBin, fElectronFraction);
81 else
82 hSystFraction = GetSystErrorHist(hPionRatio, centBin, 0);
83 hSystFraction->SetLineColor(1);
84 hSystFraction->SetMarkerStyle(1);
85 hSystFraction->SetFillStyle(1001);
86 hSystFraction->SetFillColor(kGray);
87
88 TH1D* hPionFraction = (TH1D*)hPionRatio->Clone("hPionFraction");
89 if(pid==1 && !fElectronFraction) {
90 TF1 f1("f1", "1.0", 0.0, 50.0);
91 cout << "NO ELECTRON FRACTION!!!" << endl;
92 hPionFraction->Add(&f1, -0.01); // correct for muons and electrons
93 CutHistogram(hPionFraction, 3.0, ptMax);
94 hSystFraction->Add(&f1, -0.01); // correct for muons and electrons
95 CutHistogram(hSystFraction, 3.0, ptMax);
96 }
97
98 if(pid==1 || pid ==3) {
99
100 hPionFraction->Scale(0.94); // correct for efficiency ratio
101 hSystFraction->Scale(0.94); // correct for efficiency ratio
102 } else {
103
104 TF1* effRatioK = new TF1("effRatioK", "exp(-[1]*x)+[0]", 0.0, 50.0);
105 effRatioK->SetParameters(9.82065e-01, 1.28157e+00);
106 hPionFraction->Multiply(effRatioK); // correct for efficiency ratio
107 hSystFraction->Multiply(effRatioK); // correct for efficiency ratio
108 }
109
110 if(pid==1)
111 hPionFraction->SetMarkerStyle(20);
112 else
113 hPionFraction->SetMarkerStyle(20);
114 TLatex latex;
115 latex.SetNDC();
116 latex.SetTextSize(0.05);
117
118 TCanvas* cRatios = new TCanvas("cRatios", "Particle fractions");
119 cRatios->Clear();
120 hSystFraction->GetXaxis()->SetRangeUser(0, ptMax);
121 hSystFraction->GetYaxis()->SetRangeUser(0, 1);
122 hSystFraction->Draw("E5");
123 hPionRatio->Draw("SAME");
124 hPionFraction->Draw("SAME");
125 latex.DrawLatex(0.6, 0.95, Form("%s", endName));
126
127 CreateDir(outDir);
128 if(pid==1) {
129 hPionFraction->SetName(Form("hPionFraction_%s", endName));
130 hSystFraction->SetName(Form("hPionFractionSyst_%s", endName));
131 cRatios->SaveAs(Form("%s/fractions_%s_pions.gif", outDir, endName));
132 } else if (pid==2) {
133 hPionFraction->SetName(Form("hKaonFraction_%s", endName));
134 hSystFraction->SetName(Form("hKaonFractionSyst_%s", endName));
135 cRatios->SaveAs(Form("%s/fractions_%s_kaons.gif", outDir, endName));
136 } else if (pid==3) {
137 hPionFraction->SetName(Form("hProtonFraction_%s", endName));
138 hSystFraction->SetName(Form("hProtonFractionSyst_%s", endName));
139 cRatios->SaveAs(Form("%s/fractions_%s_protons.gif", outDir, endName));
140 }
141
142 TFile* outFile = new TFile(outFileName, "RECREATE");
143 hPionFraction->Write();
144 hSystFraction->Write();
145 outFile->Close();
146}
147
148//__________________________________________________________________________________
149TH1D* GetSystErrorHist(TH1D* hRatio, Int_t centBin, TF1* electronFraction)
150{
151 TFile* fileSyst = FindFile("syst_vs_pt.root");
152
153 TF1* fSyst = (TF1*)fileSyst->Get(Form("systHigh%d", centBin));
154 fSyst->SetRange(2.0, 20.0);
155 fSyst->Print();
156
157 Double_t syst_error_mc = 0.03; // pp
158 if (centBin>0)
159 syst_error_mc = 0.05; // Pb+Pb
160
161 Double_t syst_error_correction = 0.01; // pp
162 if (centBin>0)
163 syst_error_correction = 0.03; // Pb+Pb
164
165 TH1D* hSystError = (TH1D*)hRatio->Clone("hPionFractionSyst");
166 hSystError->Multiply(fSyst);
167
168 const Int_t nBins = hSystError->GetNbinsX();
169 for(Int_t bin = 1; bin <= nBins; bin++) {
170
171 Double_t value = hRatio->GetBinContent(bin);
172 Double_t stat_error = hSystError->GetBinContent(bin);
173
174 if(value==0)
175 continue;
176
177 Double_t syst_error = stat_error*stat_error;
178 syst_error += value*value*syst_error_mc*syst_error_mc;
179 syst_error += value*value*syst_error_correction*syst_error_correction;
180
181
182 if (electronFraction) {
183
184 Double_t systEandMu = electronFraction->Eval(hRatio->GetXaxis()->GetBinCenter(bin));
185 syst_error += systEandMu*systEandMu;
186 }
187
188 syst_error = TMath::Sqrt(syst_error);
189 hSystError->SetBinContent(bin, value);
190 hSystError->SetBinError(bin, syst_error);
191 }
192
193 return hSystError;
194}
195