[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / CombineFeedDownMCSubtractionMethodsUncertainties.C

#include "TFile.h"
#include "TH1.h"
#include "TH1D.h"
#include "TH2.h"
#include "TH2F.h"
#include "TGraphAsymmErrors.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TMath.h"
#include "TROOT.h"
#include "TStyle.h"

#include "AliHFSystErr.h"

#include <Riostream.h>

//_________________________________________________________________________________________
//
//  Macro to combine the the MonteCarlo B feed-down subtraction uncertainties
//
//   Take as input the output files from the HFPtSpectrum class 
//    from both fc & Nb subtraction methods and combine the uncertainties. 
//   The final central value is set as the one from the Nb-method. 
//   The final uncertainties are defined as the envelope of both fc & Nb
//      uncertainties with respect to the new central-value.
//   The final global uncertainties are also defined and a preliminary drawing done. 
//
//
//   Usage parameters:
//      1. HFPtSpectrum fc subtraction file 
//      2. HFPtSpectrum Nb subtraction file 
//      3. Output file name
//      4. FONLL theoretical predictions file to draw on top
//      5. Decay channel as defined in the AliHFSystErr class
//
//_________________________________________________________________________________________


void CombineFeedDownMCSubtractionMethodsUncertainties(const char *fcfilename="HFPtSpectrum_D0Kpi_method1_221110_newnorm.root",
						      const char *nbfilename="HFPtSpectrum_D0Kpi_method2_221110_newnorm.root",
						      const char *outfilename="HFPtSpectrum_D0Kpi_combinedFD.root",
						      const char *thfilename="D0DplusDstarPredictions_y05.root",
						      Int_t decay=1)
{
  
  // 
  // Get fc file inputs
  TFile * fcfile = new TFile(fcfilename,"read");
  TH1D * histoSigmaCorrFc = (TH1D*)fcfile->Get("histoSigmaCorr");
  histoSigmaCorrFc->SetNameTitle("histoSigmaCorrFc","histoSigmaCorrFc");
  TGraphAsymmErrors * gSigmaCorrFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorr");
  gSigmaCorrFc->SetNameTitle("gSigmaCorrFc","gSigmaCorrFc");
  TGraphAsymmErrors * gSigmaCorrConservativeFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorrConservative");
  gSigmaCorrConservativeFc->SetNameTitle("gSigmaCorrConservativeFc","gSigmaCorrConservativeFc");

  // 
  // Get Nb file inputs
  TFile * nbfile = new TFile(nbfilename,"read");
  TH1D * histoSigmaCorrNb = (TH1D*)nbfile->Get("histoSigmaCorr");
  histoSigmaCorrNb->SetNameTitle("histoSigmaCorrNb","histoSigmaCorrNb");
  TGraphAsymmErrors * gSigmaCorrNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorr");
  gSigmaCorrNb->SetNameTitle("gSigmaCorrNb","gSigmaCorrNb");
  TGraphAsymmErrors * gSigmaCorrConservativeNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorrConservative");
  gSigmaCorrConservativeNb->SetNameTitle("gSigmaCorrConservativeNb","gSigmaCorrConservativeNb");

  //
  // Get the predictions input
  TFile *thfile = new TFile(thfilename,"read");
  TGraphAsymmErrors * thD0KpifromBprediction = (TGraphAsymmErrors*)thfile->Get("D0Kpiprediction");
  TGraphAsymmErrors * thDpluskpipiprediction = (TGraphAsymmErrors*)thfile->Get("Dpluskpipiprediction");
  TGraphAsymmErrors * thDstarD0piprediction = (TGraphAsymmErrors*)thfile->Get("DstarD0piprediction");
  thD0KpifromBprediction->SetLineColor(4);
  thD0KpifromBprediction->SetFillColor(kAzure+9);
  thDpluskpipiprediction->SetLineColor(4);
  thDpluskpipiprediction->SetFillColor(kAzure+9);
  thDstarD0piprediction->SetLineColor(4);
  thDstarD0piprediction->SetFillColor(kAzure+9);

  //
  // Get the spectra bins & limits
  Int_t nbins = histoSigmaCorrFc->GetNbinsX();
  Double_t *limits = new Double_t[nbins+1];
  Double_t xlow=0., binwidth=0.;
  for (Int_t i=1; i<=nbins; i++) {
    binwidth = histoSigmaCorrFc->GetBinWidth(i);
    xlow = histoSigmaCorrFc->GetBinLowEdge(i);
    limits[i-1] = xlow;
  }
  limits[nbins] = xlow + binwidth;


  //
  // Define a new histogram with the real-data reconstructed spectrum binning 
  //   they will be filled with central value equal to the Nb result
  //   and uncertainties taken from the envelope of the result uncertainties
  //   The systematical unc. (but FD) will also be re-calculated
  //
  TH1D * histoSigmaCorr = new TH1D("histoSigmaCorr","corrected cross-section (combined fc and Nb MC feed-down subtraction)",nbins,limits);
  TGraphAsymmErrors * gSigmaCorr = new TGraphAsymmErrors(nbins+1);
  gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (combined fc and Nb MC FD)");
  TGraphAsymmErrors * gSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
  gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (combined fc and Nb MC FD)");
  TGraphAsymmErrors * gSigmaCorrConservativePC = new TGraphAsymmErrors(nbins+1);
  gSigmaCorrConservativePC->SetNameTitle("gSigmaCorrConservativePC","Conservative gSigmaCorr (combined fc and Nb MC FD) in percentages [for drawing with AliHFSystErr]");

  //
  // Call the systematics uncertainty class for a given decay
  //   will help to compute the systematical unc. (but FD) 
  AliHFSystErr systematics(decay);

  // 
  // Loop on all the bins to do the calculations
  //
  Double_t pt=0., average = 0., averageStatUnc=0., avErrx=0., avErryl=0., avErryh=0., avErryfdl=0., avErryfdh=0.;
  Double_t avErrylPC=0., avErryhPC=0., avErryfdlPC=0., avErryfdhPC=0.;
  Double_t valFc = 0., valFcErrstat=0., valFcErrx=0., valFcErryl=0., valFcErryh=0., valFcErryfdl=0., valFcErryfdh=0.;
  Double_t valNb = 0., valNbErrstat=0., valNbErrx=0., valNbErryl=0., valNbErryh=0., valNbErryfdl=0., valNbErryfdh=0.;
  //
  for(Int_t ibin=1; ibin<=nbins; ibin++){

    // Get input values from fc method
    valFc = histoSigmaCorrFc->GetBinContent(ibin);
    pt = histoSigmaCorrFc->GetBinCenter(ibin);
    valFcErrstat = histoSigmaCorrFc->GetBinError(ibin);
    Double_t value =0., ptt=0.;
    gSigmaCorrConservativeFc->GetPoint(ibin,ptt,value);
    if ( TMath::Abs(valFc-value)>0.1 || TMath::Abs(pt-ptt)>0.1 ) 
      cout << "Hey you ! There might be a problem with the fc input file, please, have a look !" << endl;
    valFcErrx = gSigmaCorrFc->GetErrorXlow(ibin);
    valFcErryl = gSigmaCorrFc->GetErrorYlow(ibin);
    valFcErryh = gSigmaCorrFc->GetErrorYhigh(ibin);
    valFcErryfdl = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYlow(ibin) );
    valFcErryfdh = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYhigh(ibin) );

    // Get input values from Nb method
    valNb = histoSigmaCorrNb->GetBinContent(ibin);
    pt = histoSigmaCorrNb->GetBinCenter(ibin);
    valNbErrstat = histoSigmaCorrNb->GetBinError(ibin);
    gSigmaCorrConservativeNb->GetPoint(ibin,ptt,value);
    if ( TMath::Abs(valNb-value)>0.1 || TMath::Abs(pt-ptt)>0.1 ) 
      cout << "Hey you ! There might be a problem with the Nb input file, please, have a look !" << endl;
    valNbErrx = gSigmaCorrNb->GetErrorXlow(ibin);
    valNbErryl = gSigmaCorrNb->GetErrorYlow(ibin);
    valNbErryh = gSigmaCorrNb->GetErrorYhigh(ibin);
    valNbErryfdl = gSigmaCorrConservativeNb->GetErrorYlow(ibin);
    valNbErryfdh = gSigmaCorrConservativeNb->GetErrorYhigh(ibin);
    

    // Compute the FD combined value
    //    average = valNb
    average = valNb ;
    avErrx = valFcErrx;
    if ( TMath::Abs( valFcErrx - valNbErrx ) > 0.1 ) 
      cout << "Hey you ! There might be consistency problem with the fc & Nb input files, please, have a look !" << endl;
    averageStatUnc = valNbErrstat ;
//     cout << " pt=" << pt << ", average="<<average<<endl;
//     cout << "   stat unc (pc)=" << averageStatUnc/average << ", stat-fc (pc)="<<(valFcErrstat/valFc) << ", stat-Nb (pc)="<<(valNbErrstat/valNb)<<endl;
    
    // now estimate the new feed-down combined uncertainties
    Double_t minimum[2] = { (valFc - valFcErryfdl), (valNb - valNbErryfdl) };
    Double_t maximum[2] = { (valFc + valFcErryfdh), (valNb + valNbErryfdh) };
    avErryfdl = average - TMath::MinElement(2,minimum);
    avErryfdh = TMath::MaxElement(2,maximum) - average;
    avErryfdlPC = avErryfdl / average ; // in percentage
    avErryfdhPC = avErryfdh / average ; // in percentage
//     cout << " fc : val " << valFc << " + " << valFcErryfdh <<" - " << valFcErryfdl <<endl;
//     cout << " Nb : val " << valNb << " + " << valNbErryfdh <<" - " << valNbErryfdl <<endl;
//     cout << " fc  & Nb: val " << average << " + " << avErryfdh <<" - " << avErryfdl <<endl;


    // compute the global systematics
    avErrylPC = systematics.GetTotalSystErr(pt,avErryfdlPC); // in percentage
    avErryhPC = systematics.GetTotalSystErr(pt,avErryfdhPC); // in percentage
    avErryl = avErrylPC * average ;
    avErryh = avErryhPC * average ;
//     cout << "   syst av-l="<<avErryl<<", av-h="<<avErryh<<endl;
//     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;

    // fill in the histos and TGraphs
    //   fill them only when for non empty bins
    if ( average > 0.1 ) {
      histoSigmaCorr->SetBinContent(ibin,average);
      histoSigmaCorr->SetBinError(ibin,averageStatUnc);
      gSigmaCorr->SetPoint(ibin,pt,average);
      gSigmaCorr->SetPointError(ibin,valFcErrx,valFcErrx,avErryl,avErryh);
      gSigmaCorrConservative->SetPoint(ibin,pt,average);
      gSigmaCorrConservative->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdl,avErryfdh);
      gSigmaCorrConservativePC->SetPoint(ibin,pt,0.);
      gSigmaCorrConservativePC->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdlPC,avErryfdhPC);
    }

  }


  gROOT->SetStyle("Plain");
  gStyle->SetOptTitle(0);

  //
  // Plot the results
  TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,20.,100,1e4,1e8);
  histo2Draw->SetStats(0);
  histo2Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
  histo2Draw->GetXaxis()->SetTitleSize(0.05);
  histo2Draw->GetXaxis()->SetTitleOffset(0.95);
  histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}");
  histo2Draw->GetYaxis()->SetTitleSize(0.05);
  //
  TCanvas *combinefdunc = new TCanvas("combinefdunc","show the FD results combination");
  //
  histo2Draw->Draw();
  //
  histoSigmaCorrFc->SetMarkerStyle(20);
  histoSigmaCorrFc->SetMarkerColor(kGreen+2);
  histoSigmaCorrFc->SetLineColor(kGreen+2);
  histoSigmaCorrFc->Draw("esame");
  gSigmaCorrConservativeFc->SetMarkerStyle(20);
  gSigmaCorrConservativeFc->SetMarkerColor(kGreen+2);
  gSigmaCorrConservativeFc->SetLineColor(kGreen+2);
  gSigmaCorrConservativeFc->SetFillStyle(3002);
  gSigmaCorrConservativeFc->SetFillColor(kGreen);
  gSigmaCorrConservativeFc->Draw("2[]same");
  //
  histoSigmaCorrNb->SetMarkerStyle(25);
  histoSigmaCorrNb->SetMarkerColor(kViolet+5);
  histoSigmaCorrNb->SetLineColor(kViolet+5);
  histoSigmaCorrNb->Draw("esame");
  gSigmaCorrConservativeNb->SetMarkerStyle(25);
  gSigmaCorrConservativeNb->SetMarkerColor(kViolet+5);
  gSigmaCorrConservativeNb->SetLineColor(kViolet+5);
  gSigmaCorrConservativeNb->SetFillStyle(3002);
  gSigmaCorrConservativeNb->SetFillColor(kMagenta);
  gSigmaCorrConservativeNb->Draw("2[]same");
  //
  gSigmaCorrConservative->SetLineColor(kRed);
  gSigmaCorrConservative->SetLineWidth(1);
  gSigmaCorrConservative->SetFillColor(kRed);
  gSigmaCorrConservative->SetFillStyle(0);
  gSigmaCorrConservative->Draw("2");
  histoSigmaCorr->SetMarkerColor(kRed);
  histoSigmaCorr->Draw("esame");
  //
  combinefdunc->SetLogy();
  combinefdunc->Update();

  //
  // Plot the results
  TCanvas *finalresults = new TCanvas("finalresults","show all combined results");
  //
  if ( decay==1 ) {
    thD0KpifromBprediction->SetLineColor(kGreen+2);
    thD0KpifromBprediction->SetLineWidth(3);
    thD0KpifromBprediction->SetFillColor(kGreen-6);
    thD0KpifromBprediction->Draw("3CA");
    thD0KpifromBprediction->Draw("CX");
  }
  else if ( decay==2 ) {
    thDpluskpipiprediction->SetLineColor(kGreen+2);
    thDpluskpipiprediction->SetLineWidth(3);
    thDpluskpipiprediction->SetFillColor(kGreen-6);
    thDpluskpipiprediction->Draw("3CA");
    thDpluskpipiprediction->Draw("CX");
  }
  else if ( decay==3 ) {
    thDstarD0piprediction->SetLineColor(kGreen+2);
    thDstarD0piprediction->SetLineWidth(3);
    thDstarD0piprediction->SetFillColor(kGreen-6);
    thDstarD0piprediction->Draw("3CA");
    thDstarD0piprediction->Draw("CX");
  }
  //
  gSigmaCorr->SetLineColor(kRed);
  gSigmaCorr->SetLineWidth(1);
  gSigmaCorr->SetFillColor(kRed);
  gSigmaCorr->SetFillStyle(0);
  gSigmaCorr->Draw("2");
  histoSigmaCorr->SetMarkerStyle(21);
  histoSigmaCorr->SetMarkerColor(kRed);
  histoSigmaCorr->Draw("esame");
  //
  TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
  leg->SetBorderSize(0);
  leg->SetLineColor(0);
  leg->SetFillColor(0);
  leg->SetTextFont(42);
  if ( decay==1 ) leg->AddEntry(thD0KpifromBprediction,"FONLL ","fl");
  else if ( decay==2 ) leg->AddEntry(thDpluskpipiprediction,"FONLL ","fl");
  else if ( decay==3 ) leg->AddEntry(thDstarD0piprediction,"FONLL ","fl");
  leg->AddEntry(histoSigmaCorr,"data stat. unc.","pl");
  leg->AddEntry(gSigmaCorr,"data syst. unc.","f");
  leg->Draw();
  //
  finalresults->SetLogy();
  finalresults->Update();


  //
  // Draw all the systematics independently
  systematics.DrawErrors(gSigmaCorrConservativePC);


  // Write the output to a file
  TFile * out = new TFile(outfilename,"recreate");
  histoSigmaCorr->Write();
  gSigmaCorr->Write();
  gSigmaCorrConservative->Write();
  gSigmaCorrConservativePC->Write();
  out->Write();

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