Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / final_spectra / convertToMyBins.C
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
10 using 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
52 TH1D* Convert(TH1D* histIn, Double_t addSyst=0);
53 TH1D* ConvertGraph(TGraphErrors* graphIn);
54 void convertToMyBinsAA();
55 void convertToMyBinsRaa(const Char_t* fileName = "PbPb_RAA_2760GeV_200511.root");
56 void convertToMyBinsOld(const Char_t* fileName = "dndpt_all_2011-05-15.root");
57
58
59 void  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 //________________________________________________________________________________________
83 TH1D* 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 //________________________________________________________________________________________
137 TH1D* 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 //________________________________________________________________________________________
175 void 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 //________________________________________________________________________________________
240 void 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
294 void 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 }