Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / final_spectra / convertToMyBins.C
CommitLineData
4ebdd20e 1#include <TFile.h>
2#include <TH1.h>
3#include <TMath.h>
4#include <TGraphErrors.h>
5
6#include "my_tools.C"
7
8#include <iostream>
9
10using namespace std;
11
12/*
13 Info:
14
15 Method to convert charged pT histograms to have the same histogram binning
16 as our histograms have. They in general go to 100 GeV/c in pT while we only
17 go to 50. In that range the binning is EXACTLY the same (and the code checks
18 it).
19
20 In charged spectra there is NO NORMALIZATION uncertainty. We add it by hand
21 to the syst error.
22
23 * Need at some point to get the final spectra
24 * Need at some point to get the final normalization uncertainty
25
26
27 To run:
28
29 gSystem->AddIncludePath("-I../macros")
30 gSystem->AddIncludePath("-I../lib")
31 gROOT->SetMacroPath(".:../macros")
32 .L my_tools.C+
33 .L convertToMyBins.C+
34
35
36 ConvertPPfile("charged_spectra/dNdPt_pp_7TeV.root", "NOPID_pp_7TeV.root", 0.034);
37 ConvertPPfile("charged_spectra/dNdPt_pp_2.76TeV.root", "NOPID_pp_2.76TeV.root", 0.034);
38 ConvertPPfile("charged_spectra/dNdPt_pp_900GeV.root", "NOPID_pp_900GeV.root", 0.034);
39
40
41
42
43 Older methods that can be useful for AA:
44 The method convertToMyBinsAA is for pp and AA spectra.
45 The method convertToMyBinsRaa is for RAA.
46 The method convertToMyBinsOld is for AA data when the centralities had to be merged.
47
48 */
49
50
51
52TH1D* Convert(TH1D* histIn, Double_t addSyst=0);
53TH1D* ConvertGraph(TGraphErrors* graphIn);
54void convertToMyBinsAA();
55void convertToMyBinsRaa(const Char_t* fileName = "PbPb_RAA_2760GeV_200511.root");
56void convertToMyBinsOld(const Char_t* fileName = "dndpt_all_2011-05-15.root");
57
58
59void ConvertPPfile(const Char_t* inFileName, const Char_t* outFileName, Double_t normSyst=0.0)
60{
61 TFile* fileDataCharged = FindFileFresh(inFileName);
62
63 if(!fileDataCharged)
64 return;
65
66 TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat");
67 TH1D* hPtCharged = Convert(hPtHelp);
68 hPtCharged->SetName("hDnDptCharged_PP");
69 hPtCharged->SetMarkerStyle(29);
70
71 hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst");
72 TH1D* hPtChargedSyst = Convert(hPtHelp, normSyst);
73 hPtChargedSyst->SetName("hDnDptChargedSyst_PP");
74 hPtChargedSyst->SetMarkerStyle(29);
75
76 TFile* fileOut = new TFile(outFileName, "RECREATE");
77 hPtCharged->Write();
78 hPtChargedSyst->Write();
79 fileOut->Close();
80}
81
82//________________________________________________________________________________________
83TH1D* Convert(TH1D* histIn, Double_t addSyst)
84{
85 const Int_t nPtBins = 68;
86 Double_t xBins[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
87 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
88 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 ,
89 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 ,
90 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0,
91 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
92 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
93
94
95 TH1D* hist = new TH1D("hist", "", nPtBins, xBins);
96 hist->SetTitle(histIn->GetTitle());
97 hist->GetXaxis()->SetTitle(histIn->GetXaxis()->GetTitle());
98 hist->GetYaxis()->SetTitle(histIn->GetYaxis()->GetTitle());
99
100 const Double_t deltapt = 0.0001;
101
102 for(Int_t bin = 1; bin <= nPtBins; bin++) {
103
104 // check bin size
105 if(TMath::Abs(hist->GetXaxis()->GetBinLowEdge(bin) -
106 histIn->GetXaxis()->GetBinLowEdge(bin)) > deltapt) {
107 cout << "pt edge low does not agree!!!!!!!!!!!!!!!" << endl;
108 return 0;
109 }
110 if(TMath::Abs(hist->GetXaxis()->GetBinUpEdge(bin) -
111 histIn->GetXaxis()->GetBinUpEdge(bin)) > deltapt) {
112 cout << "pt edge high does not agree!!!!!!!!!!!!!!!" << endl;
113 return 0;
114 }
115 if(TMath::Abs(hist->GetXaxis()->GetBinCenter(bin) -
116 histIn->GetXaxis()->GetBinCenter(bin)) > deltapt) {
117 cout << "pt center does not agree!!!!!!!!!!!!!!!" << endl;
118 return 0;
119 }
120
121
122 hist->SetBinContent(bin, histIn->GetBinContent(bin));
123
124 Double_t error = histIn->GetBinError(bin);
125 if(addSyst) {
126
127 Double_t extra_error = addSyst*histIn->GetBinContent(bin);
128 error = TMath::Sqrt(error*error + extra_error*extra_error);
129 }
130 hist->SetBinError(bin, error);
131 }
132
133 return hist;
134}
135
136//________________________________________________________________________________________
137TH1D* ConvertGraph(TGraphErrors* graphIn)
138{
139 const Int_t nPtBins = 68;
140 Double_t xBins[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
141 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
142 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 ,
143 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 ,
144 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0,
145 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
146 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
147
148
149 TH1D* hist = new TH1D("hist", "", nPtBins, xBins);
150 // hist->SetTitle(graphIn->GetTitle());
151 // hist->GetXaxis()->SetTitle(graphIn->GetXaxis()->GetTitle());
152 // hist->GetYaxis()->SetTitle(graphIn->GetYaxis()->GetTitle());
153
154 const Double_t deltapt = 0.0001;
155
156 // Jacek does not fill the first 3 bins of his graph?
157 for(Int_t bin = 4; bin <= nPtBins; bin++) {
158
159 // check bin size
160 if(TMath::Abs(hist->GetXaxis()->GetBinCenter(bin) -
161 graphIn->GetX()[bin-1]) > deltapt) {
162 cout << "pt center does not agree!!!!!!!!!!!!!!!" << endl;
163 return 0;
164 }
165
166
167 hist->SetBinContent(bin, graphIn->GetY()[bin-1]);
168 hist->SetBinError(bin, graphIn->GetEY()[bin-1]);
169 }
170
171 return hist;
172}
173
174//________________________________________________________________________________________
175void convertToMyBinsAA()
176{
177 const Int_t nCent = 7;
178 const Int_t cent[nCent+1] = {-5, 0, 5, 10, 20, 40, 60, 80};
179 const Int_t colors[nCent] = {kBlack, kRed+1, kOrange+7, kOrange, kGreen+2, kCyan+2, kBlue+1};
180
181 TH1D* hPtCharged[nCent] = {0, 0, 0, 0, 0, 0, 0};
182 TH1D* hPtChargedSyst[nCent] = {0, 0, 0, 0, 0, 0, 0};
183
184 for (Int_t i = 0; i < nCent; i++) {
185
186 if(cent[i]==-5) {
187
188 TFile* fileDataCharged = FindFileFresh("dNdPt_pp_2.76TeV.root");
189
190 if(!fileDataCharged)
191 return;
192
193 TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat");
194 hPtCharged[i] = Convert(hPtHelp);
195 hPtCharged[i]->SetName("hDnDptCharged_PP");
196 hPtCharged[i]->SetMarkerColor(colors[i]);
197 hPtCharged[i]->SetMarkerStyle(29);
198 hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst");
199 hPtChargedSyst[i] = Convert(hPtHelp, 0.034); // 3.4 % norm error
200 /*
201 Dear Peter,
202
203 for 2.76TeV the systematic error for MB->INEL is 3.40%,, we decided to
204 drop the uncertainty from VdM scan since we don't really use it...
205
206 Michael
207 */
208 hPtChargedSyst[i]->SetName("hDnDptChargedSyst_PP");
209 hPtChargedSyst[i]->SetMarkerColor(colors[i]);
210 hPtChargedSyst[i]->SetMarkerStyle(29);
211 } else {
212
213 TFile* fileDataCharged = FindFileFresh(Form("dNdPt_PbPb_2.76ATeV_centrality_%d-%d.root", cent[i], cent[i+1]));
214
215 if(!fileDataCharged)
216 return;
217
218 TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat");
219 hPtCharged[i] = Convert(hPtHelp);
220 hPtCharged[i]->SetName(Form("hDnDptCharged_%d_%d", cent[i], cent[i+1]));
221 hPtCharged[i]->SetMarkerColor(colors[i]);
222 hPtCharged[i]->SetMarkerStyle(29);
223 hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst");
224 hPtChargedSyst[i] = Convert(hPtHelp);
225 hPtChargedSyst[i]->SetName(Form("hDnDptChargedSyst_%d_%d", cent[i], cent[i+1]));
226 hPtChargedSyst[i]->SetMarkerColor(colors[i]);
227 hPtChargedSyst[i]->SetMarkerStyle(29);
228 }
229 }
230
231 TFile* fileOut = new TFile("dNdPt_Charged.root", "RECREATE");
232 for (Int_t i = 0; i < nCent; i++) {
233 hPtCharged[i]->Write();
234 hPtChargedSyst[i]->Write();
235 }
236 fileOut->Close();
237}
238
239//________________________________________________________________________________________
240void convertToMyBinsRaa(const Char_t* fileName)
241{
242 TFile* fileData = FindFileFresh(fileName);
243
244 if(!fileData)
245 return;
246
247 TGraphErrors* hRaa_0_5 = (TGraphErrors*)fileData->Get("raa_c0_5_stat");
248 TH1D* hRaaCharged_0_5 = ConvertGraph(hRaa_0_5);
249 hRaaCharged_0_5->SetName("hRaaCharged_0_5");
250
251 TGraphErrors* hRaa_5_10 = (TGraphErrors*)fileData->Get("raa_c5_10_stat");
252 TH1D* hRaaCharged_5_10 = ConvertGraph(hRaa_5_10);
253 hRaaCharged_5_10->SetName("hRaaCharged_5_10");
254
255 TGraphErrors* hRaa_10_20 = (TGraphErrors*)fileData->Get("raa_c10_20_stat");
256 TH1D* hRaaCharged_10_20 = ConvertGraph(hRaa_10_20);
257 hRaaCharged_10_20->SetName("hRaaCharged_10_20");
258
259 TGraphErrors* hRaa_20_40 = (TGraphErrors*)fileData->Get("raa_c20_40_stat");
260 TH1D* hRaaCharged_20_40 = ConvertGraph(hRaa_20_40);
261 hRaaCharged_20_40->SetName("hRaaCharged_20_40");
262
263 TGraphErrors* hRaa_40_60 = (TGraphErrors*)fileData->Get("raa_c40_60_stat");
264 TH1D* hRaaCharged_40_60 = ConvertGraph(hRaa_40_60);
265 hRaaCharged_40_60->SetName("hRaaCharged_40_60");
266
267 TGraphErrors* hRaa_60_80 = (TGraphErrors*)fileData->Get("raa_c60_80_stat");
268 TH1D* hRaaCharged_60_80 = ConvertGraph(hRaa_60_80);
269 hRaaCharged_60_80->SetName("hRaaCharged_60_80");
270
271 TGraphErrors* hRaa_0_20 = (TGraphErrors*)fileData->Get("raa_c0_20_stat");
272 TH1D* hRaaCharged_0_20 = ConvertGraph(hRaa_0_20);
273 hRaaCharged_0_20->SetName("hRaaCharged_0_20");
274
275 TGraphErrors* hRaa_40_80 = (TGraphErrors*)fileData->Get("raa_c40_80_stat");
276 TH1D* hRaaCharged_40_80 = ConvertGraph(hRaa_40_80);
277 hRaaCharged_40_80->SetName("hRaaCharged_40_80");
278
279
280 TFile* fileOut = new TFile("Raa_Charged.root", "RECREATE");
281 hRaaCharged_0_5->Write();
282 hRaaCharged_5_10->Write();
283 hRaaCharged_10_20->Write();
284 hRaaCharged_20_40->Write();
285 hRaaCharged_40_60->Write();
286 hRaaCharged_60_80->Write();
287 hRaaCharged_0_20->Write();
288 hRaaCharged_40_80->Write();
289 fileOut->Close();
290}
291
292//________________________________________________________________________________________
293
294void convertToMyBinsOld(const Char_t* fileName)
295{
296 TFile* fileData = FindFileFresh(fileName);
297
298 if(!fileData)
299 return;
300
301 TH1D* hDnDPt_0_5 = (TH1D*)fileData->Get("dNdPt_PbPb_c0_5");
302 TH1D* hDnDptCharged_0_5 = Convert(hDnDPt_0_5);
303 hDnDptCharged_0_5->SetName("hDnDptCharged_0_5");
304
305 TH1D* hDnDPt_5_10 = (TH1D*)fileData->Get("dNdPt_PbPb_c5_10");
306 TH1D* hDnDptCharged_5_10 = Convert(hDnDPt_5_10);
307 hDnDptCharged_5_10->SetName("hDnDptCharged_5_10");
308
309 TH1D* hDnDPt_10_20 = (TH1D*)fileData->Get("dNdPt_PbPb_c10_20");
310 TH1D* hDnDptCharged_10_20 = Convert(hDnDPt_10_20);
311 hDnDptCharged_10_20->SetName("hDnDptCharged_10_20");
312
313 TH1D* hDnDPt_20_40 = (TH1D*)fileData->Get("dNdPt_PbPb_c20_30");
314 TH1D* hHelp = (TH1D*)fileData->Get("dNdPt_PbPb_c30_40");
315 hDnDPt_20_40->Add(hHelp);
316 hDnDPt_20_40->Scale(1.0/2.0);
317 TH1D* hDnDptCharged_20_40 = Convert(hDnDPt_20_40);
318 hDnDptCharged_20_40->SetName("hDnDptCharged_20_40");
319
320 TH1D* hDnDPt_40_60 = (TH1D*)fileData->Get("dNdPt_PbPb_c40_50");
321 hHelp = (TH1D*)fileData->Get("dNdPt_PbPb_c50_60");
322 hDnDPt_40_60->Add(hHelp);
323 hDnDPt_40_60->Scale(1.0/2.0);
324 TH1D* hDnDptCharged_40_60 = Convert(hDnDPt_40_60);
325 hDnDptCharged_40_60->SetName("hDnDptCharged_40_60");
326
327 TH1D* hDnDPt_60_80 = (TH1D*)fileData->Get("dNdPt_PbPb_c60_70");
328 hHelp = (TH1D*)fileData->Get("dNdPt_PbPb_c70_80");
329 hDnDPt_60_80->Add(hHelp);
330 hDnDPt_60_80->Scale(1.0/2.0);
331 TH1D* hDnDptCharged_60_80 = Convert(hDnDPt_60_80);
332 hDnDptCharged_60_80->SetName("hDnDptCharged_60_80");
333
334
335 TFile* fileOut = new TFile("dNdPt_Charged.root", "RECREATE");
336 hDnDptCharged_0_5->Write();
337 hDnDptCharged_5_10->Write();
338 hDnDptCharged_10_20->Write();
339 hDnDptCharged_20_40->Write();
340 hDnDptCharged_40_60->Write();
341 hDnDptCharged_60_80->Write();
342 fileOut->Close();
343}