Update (Zaida)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / CombineFeedDownMCSubtractionMethodsUncertainties.C
1 #include "TFile.h"
2 #include "TH1.h"
3 #include "TH1D.h"
4 #include "TH2.h"
5 #include "TH2F.h"
6 #include "TGraphAsymmErrors.h"
7 #include "TCanvas.h"
8 #include "TLegend.h"
9 #include "TMath.h"
10 #include "TROOT.h"
11 #include "TStyle.h"
12
13 #include "AliHFSystErr.h"
14
15 #include <Riostream.h>
16
17 //_________________________________________________________________________________________
18 //
19 //  Macro to combine the the MonteCarlo B feed-down subtraction uncertainties
20 //
21 //   Take as input the output files from the HFPtSpectrum class 
22 //    from both fc & Nb subtraction methods and combine the uncertainties. 
23 //   The final central value is set as the one from the Nb-method. 
24 //   The final uncertainties are defined as the envelope of both fc & Nb
25 //      uncertainties with respect to the new central-value.
26 //   The final global uncertainties are also defined and a preliminary drawing done. 
27 //
28 //
29 //   Usage parameters:
30 //      1. HFPtSpectrum fc subtraction file 
31 //      2. HFPtSpectrum Nb subtraction file 
32 //      3. Output file name
33 //      4. FONLL theoretical predictions file to draw on top
34 //      5. Decay channel as defined in the AliHFSystErr class
35 //
36 //_________________________________________________________________________________________
37
38 enum centrality{ kpp, k010, k020, k2040, k4060, k6080, k4080, k80100 };
39
40 void CombineFeedDownMCSubtractionMethodsUncertainties(const char *fcfilename="HFPtSpectrum_D0Kpi_method1_221110_newnorm.root",
41                                                       const char *nbfilename="HFPtSpectrum_D0Kpi_method2_221110_newnorm.root",
42                                                       const char *outfilename="HFPtSpectrum_D0Kpi_combinedFD.root",
43                                                       const char *thfilename="D0DplusDstarPredictions_y05.root",
44                                                       Int_t decay=1, Int_t centrality=kpp)
45 {
46   
47   // 
48   // Get fc file inputs
49   TFile * fcfile = new TFile(fcfilename,"read");
50   TH1D * histoSigmaCorrFc = (TH1D*)fcfile->Get("histoSigmaCorr");
51   histoSigmaCorrFc->SetNameTitle("histoSigmaCorrFc","histoSigmaCorrFc");
52   TGraphAsymmErrors * gSigmaCorrFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorr");
53   gSigmaCorrFc->SetNameTitle("gSigmaCorrFc","gSigmaCorrFc");
54   TGraphAsymmErrors * gSigmaCorrConservativeFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorrConservative");
55   gSigmaCorrConservativeFc->SetNameTitle("gSigmaCorrConservativeFc","gSigmaCorrConservativeFc");
56
57   // 
58   // Get Nb file inputs
59   TFile * nbfile = new TFile(nbfilename,"read");
60   TH1D * histoSigmaCorrNb = (TH1D*)nbfile->Get("histoSigmaCorr");
61   histoSigmaCorrNb->SetNameTitle("histoSigmaCorrNb","histoSigmaCorrNb");
62   TGraphAsymmErrors * gSigmaCorrNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorr");
63   gSigmaCorrNb->SetNameTitle("gSigmaCorrNb","gSigmaCorrNb");
64   TGraphAsymmErrors * gSigmaCorrConservativeNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorrConservative");
65   gSigmaCorrConservativeNb->SetNameTitle("gSigmaCorrConservativeNb","gSigmaCorrConservativeNb");
66
67   //
68   // Get the predictions input
69   TFile *thfile = new TFile(thfilename,"read");
70   TGraphAsymmErrors * thD0KpifromBprediction = (TGraphAsymmErrors*)thfile->Get("D0Kpiprediction");
71   TGraphAsymmErrors * thDpluskpipiprediction = (TGraphAsymmErrors*)thfile->Get("Dpluskpipiprediction");
72   TGraphAsymmErrors * thDstarD0piprediction = (TGraphAsymmErrors*)thfile->Get("DstarD0piprediction");
73   thD0KpifromBprediction->SetLineColor(4);
74   thD0KpifromBprediction->SetFillColor(kAzure+9);
75   thDpluskpipiprediction->SetLineColor(4);
76   thDpluskpipiprediction->SetFillColor(kAzure+9);
77   thDstarD0piprediction->SetLineColor(4);
78   thDstarD0piprediction->SetFillColor(kAzure+9);
79
80   //
81   // Get the spectra bins & limits
82   Int_t nbins = histoSigmaCorrFc->GetNbinsX();
83   Double_t *limits = new Double_t[nbins+1];
84   Double_t xlow=0., binwidth=0.;
85   for (Int_t i=1; i<=nbins; i++) {
86     binwidth = histoSigmaCorrFc->GetBinWidth(i);
87     xlow = histoSigmaCorrFc->GetBinLowEdge(i);
88     limits[i-1] = xlow;
89   }
90   limits[nbins] = xlow + binwidth;
91
92
93   //
94   // Define a new histogram with the real-data reconstructed spectrum binning 
95   //   they will be filled with central value equal to the Nb result
96   //   and uncertainties taken from the envelope of the result uncertainties
97   //   The systematical unc. (but FD) will also be re-calculated
98   //
99   TH1D * histoSigmaCorr = new TH1D("histoSigmaCorr","corrected cross-section (combined fc and Nb MC feed-down subtraction)",nbins,limits);
100   TGraphAsymmErrors * gSigmaCorr = new TGraphAsymmErrors(nbins+1);
101   gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (combined fc and Nb MC FD)");
102   TGraphAsymmErrors * gSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
103   gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (combined fc and Nb MC FD)");
104   TGraphAsymmErrors * gSigmaCorrConservativePC = new TGraphAsymmErrors(nbins+1);
105   gSigmaCorrConservativePC->SetNameTitle("gSigmaCorrConservativePC","Conservative gSigmaCorr (combined fc and Nb MC FD) in percentages [for drawing with AliHFSystErr]");
106
107   //
108   // Call the systematics uncertainty class for a given decay
109   //   will help to compute the systematical unc. (but FD) 
110   AliHFSystErr systematics;
111   if( centrality!=kpp)  {
112     systematics.SetCollisionType(1);
113     if ( centrality == k020 ) {
114       systematics.SetCentrality(020);
115     }
116     else if ( centrality == k4080 ) {
117       systematics.SetCentrality(4080);
118     }
119     else { 
120       cout << " Systematics not yet implemented " << endl;
121       return;
122     }
123   }
124   else { systematics.SetCollisionType(0); }
125   systematics.Init(decay);
126
127   // 
128   // Loop on all the bins to do the calculations
129   //
130   Double_t pt=0., average = 0., averageStatUnc=0., avErrx=0., avErryl=0., avErryh=0., avErryfdl=0., avErryfdh=0.;
131   Double_t avErrylPC=0., avErryhPC=0., avErryfdlPC=0., avErryfdhPC=0.;
132   Double_t valFc = 0., valFcErrstat=0., valFcErrx=0., valFcErryl=0., valFcErryh=0., valFcErryfdl=0., valFcErryfdh=0.;
133   Double_t valNb = 0., valNbErrstat=0., valNbErrx=0., valNbErryl=0., valNbErryh=0., valNbErryfdl=0., valNbErryfdh=0.;
134   //
135   for(Int_t ibin=1; ibin<=nbins; ibin++){
136
137     // Get input values from fc method
138     valFc = histoSigmaCorrFc->GetBinContent(ibin);
139     pt = histoSigmaCorrFc->GetBinCenter(ibin);
140     valFcErrstat = histoSigmaCorrFc->GetBinError(ibin);
141     Double_t value =0., ptt=0.;
142     gSigmaCorrConservativeFc->GetPoint(ibin,ptt,value);
143     if ( TMath::Abs(valFc-value)>0.1 || TMath::Abs(pt-ptt)>0.1 ) 
144       cout << "Hey you ! There might be a problem with the fc input file, please, have a look !" << endl;
145     valFcErrx = gSigmaCorrFc->GetErrorXlow(ibin);
146     valFcErryl = gSigmaCorrFc->GetErrorYlow(ibin);
147     valFcErryh = gSigmaCorrFc->GetErrorYhigh(ibin);
148     valFcErryfdl = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYlow(ibin) );
149     valFcErryfdh = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYhigh(ibin) );
150
151     // Get input values from Nb method
152     valNb = histoSigmaCorrNb->GetBinContent(ibin);
153     pt = histoSigmaCorrNb->GetBinCenter(ibin);
154     valNbErrstat = histoSigmaCorrNb->GetBinError(ibin);
155     gSigmaCorrConservativeNb->GetPoint(ibin,ptt,value);
156     if ( TMath::Abs(valNb-value)>0.1 || TMath::Abs(pt-ptt)>0.1 ) 
157       cout << "Hey you ! There might be a problem with the Nb input file, please, have a look !" << endl;
158     valNbErrx = gSigmaCorrNb->GetErrorXlow(ibin);
159     valNbErryl = gSigmaCorrNb->GetErrorYlow(ibin);
160     valNbErryh = gSigmaCorrNb->GetErrorYhigh(ibin);
161     valNbErryfdl = gSigmaCorrConservativeNb->GetErrorYlow(ibin);
162     valNbErryfdh = gSigmaCorrConservativeNb->GetErrorYhigh(ibin);
163     
164
165     // Compute the FD combined value
166     //    average = valNb
167     average = valNb ;
168     avErrx = valFcErrx;
169     if ( TMath::Abs( valFcErrx - valNbErrx ) > 0.1 ) 
170       cout << "Hey you ! There might be consistency problem with the fc & Nb input files, please, have a look !" << endl;
171     averageStatUnc = valNbErrstat ;
172 //     cout << " pt=" << pt << ", average="<<average<<endl;
173 //     cout << "   stat unc (pc)=" << averageStatUnc/average << ", stat-fc (pc)="<<(valFcErrstat/valFc) << ", stat-Nb (pc)="<<(valNbErrstat/valNb)<<endl;
174     
175     // now estimate the new feed-down combined uncertainties
176     Double_t minimum[2] = { (valFc - valFcErryfdl), (valNb - valNbErryfdl) };
177     Double_t maximum[2] = { (valFc + valFcErryfdh), (valNb + valNbErryfdh) };
178     avErryfdl = average - TMath::MinElement(2,minimum);
179     avErryfdh = TMath::MaxElement(2,maximum) - average;
180     avErryfdlPC = avErryfdl / average ; // in percentage
181     avErryfdhPC = avErryfdh / average ; // in percentage
182 //     cout << " fc : val " << valFc << " + " << valFcErryfdh <<" - " << valFcErryfdl <<endl;
183 //     cout << " Nb : val " << valNb << " + " << valNbErryfdh <<" - " << valNbErryfdl <<endl;
184 //     cout << " fc  & Nb: val " << average << " + " << avErryfdh <<" - " << avErryfdl <<endl;
185
186
187     // compute the global systematics
188     avErrylPC = systematics.GetTotalSystErr(pt,avErryfdlPC); // in percentage
189     avErryhPC = systematics.GetTotalSystErr(pt,avErryfdhPC); // in percentage
190     avErryl = avErrylPC * average ;
191     avErryh = avErryhPC * average ;
192 //     cout << "   syst av-l="<<avErryl<<", av-h="<<avErryh<<endl;
193 //     cout << "     fd-l-pc="<<avErryfdlPC<<", fd-h-pc="<<avErryfdhPC<<", syst err(no fd)-pc="<<systematics.GetTotalSystErr(pt)<<", av-l-pc="<<avErrylPC<<", av-h-pc="<<avErryhPC<<endl;
194
195     // fill in the histos and TGraphs
196     //   fill them only when for non empty bins
197     if ( average > 0.1 ) {
198       histoSigmaCorr->SetBinContent(ibin,average);
199       histoSigmaCorr->SetBinError(ibin,averageStatUnc);
200       gSigmaCorr->SetPoint(ibin,pt,average);
201       gSigmaCorr->SetPointError(ibin,valFcErrx,valFcErrx,avErryl,avErryh);
202       gSigmaCorrConservative->SetPoint(ibin,pt,average);
203       gSigmaCorrConservative->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdl,avErryfdh);
204       gSigmaCorrConservativePC->SetPoint(ibin,pt,0.);
205       gSigmaCorrConservativePC->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdlPC,avErryfdhPC);
206     }
207
208   }
209
210
211   gROOT->SetStyle("Plain");
212   gStyle->SetOptTitle(0);
213
214   //
215   // Plot the results
216   TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,20.,100,1e4,1e10);
217   histo2Draw->SetStats(0);
218   histo2Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
219   histo2Draw->GetXaxis()->SetTitleSize(0.05);
220   histo2Draw->GetXaxis()->SetTitleOffset(0.95);
221   histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}");
222   histo2Draw->GetYaxis()->SetTitleSize(0.05);
223   //
224   TCanvas *combinefdunc = new TCanvas("combinefdunc","show the FD results combination");
225   //
226   histo2Draw->Draw();
227   //
228   histoSigmaCorrFc->SetMarkerStyle(20);
229   histoSigmaCorrFc->SetMarkerColor(kGreen+2);
230   histoSigmaCorrFc->SetLineColor(kGreen+2);
231   histoSigmaCorrFc->Draw("esame");
232   gSigmaCorrConservativeFc->SetMarkerStyle(20);
233   gSigmaCorrConservativeFc->SetMarkerColor(kGreen+2);
234   gSigmaCorrConservativeFc->SetLineColor(kGreen+2);
235   gSigmaCorrConservativeFc->SetFillStyle(3002);
236   gSigmaCorrConservativeFc->SetFillColor(kGreen);
237   gSigmaCorrConservativeFc->Draw("2[]same");
238   //
239   histoSigmaCorrNb->SetMarkerStyle(25);
240   histoSigmaCorrNb->SetMarkerColor(kViolet+5);
241   histoSigmaCorrNb->SetLineColor(kViolet+5);
242   histoSigmaCorrNb->Draw("esame");
243   gSigmaCorrConservativeNb->SetMarkerStyle(25);
244   gSigmaCorrConservativeNb->SetMarkerColor(kViolet+5);
245   gSigmaCorrConservativeNb->SetLineColor(kViolet+5);
246   gSigmaCorrConservativeNb->SetFillStyle(3002);
247   gSigmaCorrConservativeNb->SetFillColor(kMagenta);
248   gSigmaCorrConservativeNb->Draw("2[]same");
249   //
250   gSigmaCorrConservative->SetLineColor(kRed);
251   gSigmaCorrConservative->SetLineWidth(1);
252   gSigmaCorrConservative->SetFillColor(kRed);
253   gSigmaCorrConservative->SetFillStyle(0);
254   gSigmaCorrConservative->Draw("2");
255   histoSigmaCorr->SetMarkerColor(kRed);
256   histoSigmaCorr->Draw("esame");
257   //
258   combinefdunc->SetLogy();
259   combinefdunc->Update();
260
261   //
262   // Plot the results
263   TCanvas *finalresults = new TCanvas("finalresults","show all combined results");
264   //
265   if ( decay==1 ) {
266     thD0KpifromBprediction->SetLineColor(kGreen+2);
267     thD0KpifromBprediction->SetLineWidth(3);
268     thD0KpifromBprediction->SetFillColor(kGreen-6);
269     thD0KpifromBprediction->Draw("3CA");
270     thD0KpifromBprediction->Draw("CX");
271   }
272   else if ( decay==2 ) {
273     thDpluskpipiprediction->SetLineColor(kGreen+2);
274     thDpluskpipiprediction->SetLineWidth(3);
275     thDpluskpipiprediction->SetFillColor(kGreen-6);
276     thDpluskpipiprediction->Draw("3CA");
277     thDpluskpipiprediction->Draw("CX");
278   }
279   else if ( decay==3 ) {
280     thDstarD0piprediction->SetLineColor(kGreen+2);
281     thDstarD0piprediction->SetLineWidth(3);
282     thDstarD0piprediction->SetFillColor(kGreen-6);
283     thDstarD0piprediction->Draw("3CA");
284     thDstarD0piprediction->Draw("CX");
285   }
286   //
287   gSigmaCorr->SetLineColor(kRed);
288   gSigmaCorr->SetLineWidth(1);
289   gSigmaCorr->SetFillColor(kRed);
290   gSigmaCorr->SetFillStyle(0);
291   gSigmaCorr->Draw("2");
292   histoSigmaCorr->SetMarkerStyle(21);
293   histoSigmaCorr->SetMarkerColor(kRed);
294   histoSigmaCorr->Draw("esame");
295   //
296   TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
297   leg->SetBorderSize(0);
298   leg->SetLineColor(0);
299   leg->SetFillColor(0);
300   leg->SetTextFont(42);
301   if ( decay==1 ) leg->AddEntry(thD0KpifromBprediction,"FONLL ","fl");
302   else if ( decay==2 ) leg->AddEntry(thDpluskpipiprediction,"FONLL ","fl");
303   else if ( decay==3 ) leg->AddEntry(thDstarD0piprediction,"FONLL ","fl");
304   leg->AddEntry(histoSigmaCorr,"data stat. unc.","pl");
305   leg->AddEntry(gSigmaCorr,"data syst. unc.","f");
306   leg->Draw();
307   //
308   finalresults->SetLogy();
309   finalresults->Update();
310
311
312   //
313   // Draw all the systematics independently
314   systematics.DrawErrors(gSigmaCorrConservativePC);
315
316
317   // Write the output to a file
318   TFile * out = new TFile(outfilename,"recreate");
319   histoSigmaCorr->Write();
320   gSigmaCorr->Write();
321   gSigmaCorrConservative->Write();
322   gSigmaCorrConservativePC->Write();
323   out->Write();
324
325 }