Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / corrected_fraction / corrected_fraction.C
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
41 void 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");
44 TH1D* GetSystErrorHist(TH1D* hRatio, Int_t centBin, TF1* electronFraction);
45
46
47 //________________________________________________________________________________________
48 void 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 //__________________________________________________________________________________
149 TH1D* 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