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

#if !defined(__CINT__) || defined(__MAKECINT__)
#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>
#endif

/* $Id$ */ 


//_________________________________________________________________________________________
//
//  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
//
//_________________________________________________________________________________________

enum centrality{ kpp7, kpp276, k010, k020, k2040, k4060, k6080, k4080, k80100 };

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, Int_t centrality=kpp7)
{
  
  // 
  // 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","Cross section (fc prompt fraction)");
  TGraphAsymmErrors * gFcConservativeFc = (TGraphAsymmErrors*)fcfile->Get("gFcConservative");
  gFcConservativeFc->SetNameTitle("gFcConservativeFc","fc prompt fraction");

  // 
  // 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","Cross section (Nb prompt fraction)");
  TGraphAsymmErrors * gFcConservativeNb = (TGraphAsymmErrors*)nbfile->Get("gFcConservative");
  gFcConservativeNb->SetNameTitle("gFcConservativeNb","Nb prompt fraction");

  //
  // 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");
  TGraphAsymmErrors * thDsKKpiprediction = (TGraphAsymmErrors*)thfile->Get("DsKkpiprediction");

  thD0KpifromBprediction->SetLineColor(4);
  thD0KpifromBprediction->SetFillColor(kAzure+9);
  thDpluskpipiprediction->SetLineColor(4);
  thDpluskpipiprediction->SetFillColor(kAzure+9);
  thDstarD0piprediction->SetLineColor(4);
  thDstarD0piprediction->SetFillColor(kAzure+9);
  thDsKKpiprediction->SetLineColor(4);
  thDsKKpiprediction->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 * gFcCorrConservative = new TGraphAsymmErrors(nbins+1);
  gFcCorrConservative->SetNameTitle("gFcCorrConservative","Combined prompt fraction");
  TGraphAsymmErrors * gSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
  gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Cross section (combined prompt fraction)");
  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;
  if( centrality==kpp276 ) {
    systematics.SetIsLowEnergy(true);
  } else if( centrality!=kpp7 )  {
    systematics.SetCollisionType(1);
    if ( centrality == k020 ) {
      systematics.SetCentrality("020");
    }
    else if ( centrality == k4080 ) {
      systematics.SetCentrality("4080");
    }
    else { 
      cout << " Systematics not yet implemented " << endl;
      return;
    }
  }
  else { systematics.SetCollisionType(0); }
  systematics.Init(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.;
  Double_t corrfd = 0., corrfdl=0., corrfdh=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 (value<=0.) continue;
    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) );
    Double_t valfdFc = 0., x=0.;
    gFcConservativeFc->GetPoint(ibin,x,valfdFc);
    Double_t valfdFch = gFcConservativeFc->GetErrorYhigh(ibin);
    Double_t valfdFcl = gFcConservativeFc->GetErrorYlow(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);
    Double_t valfdNb = 0.;
    gFcConservativeNb->GetPoint(ibin,x,valfdNb);
    Double_t valfdNbh = gFcConservativeNb->GetErrorYhigh(ibin);
    Double_t valfdNbl = gFcConservativeNb->GetErrorYlow(ibin);
    

    // Compute the FD combined value
    //    average = valNb
    average = valNb ;
    corrfd = valfdNb;
    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;
    Double_t minimumfc[2] = { (valfdNb - valfdNbl), (valfdFc - valfdFcl) };
    Double_t maximumfc[2] = { (valfdNb + valfdNbh), (valfdFc + valfdFch) };
    corrfdl = corrfd - TMath::MinElement(2,minimumfc);
    corrfdh = TMath::MaxElement(2,maximumfc) - corrfd;


    // 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);
      gFcCorrConservative->SetPoint(ibin,pt,corrfd);
      gFcCorrConservative->SetPointError(ibin,valFcErrx,valFcErrx,corrfdl,corrfdh);
    }

  }


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

  //
  // Plot the results
  TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,20.,100,1e3,5e7);
  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(3004);//2);
  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(kOrange+7);//kViolet+5);
  gSigmaCorrConservativeNb->SetLineColor(kOrange+7);//kOrange+7);//kViolet+5);
  gSigmaCorrConservativeNb->SetFillStyle(3018);//02);
  gSigmaCorrConservativeNb->SetFillColor(kMagenta);
  gSigmaCorrConservativeNb->Draw("2[]same");
  //
  gSigmaCorrConservative->SetLineColor(kRed);
  gSigmaCorrConservative->SetLineWidth(2);
  gSigmaCorrConservative->SetFillColor(kRed);
  gSigmaCorrConservative->SetFillStyle(0);
  gSigmaCorrConservative->Draw("2");
  histoSigmaCorr->SetMarkerColor(kRed);
  histoSigmaCorr->Draw("esame");
  //
  //
  TLegend* leg=combinefdunc->BuildLegend();
  leg->SetFillStyle(0);
  combinefdunc->SetLogy();
  combinefdunc->Update();

  TCanvas *combinefcunc = new TCanvas("combinefcunc","show the fc FD results combination");
  //
  TH2F *histo3Draw = new TH2F("histo3Draw","histo3 (for drawing)",100,0,20.,10,0.,1.);
  histo3Draw->SetStats(0);
  histo3Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
  histo3Draw->GetXaxis()->SetTitleSize(0.05);
  histo3Draw->GetXaxis()->SetTitleOffset(0.95);
  histo3Draw->GetYaxis()->SetTitle("Prompt fraction of the raw yields");
  histo3Draw->GetYaxis()->SetTitleSize(0.05);
  histo3Draw->Draw();
  //
  gFcConservativeFc->SetMarkerStyle(20);
  gFcConservativeFc->SetMarkerColor(kGreen+2);
  gFcConservativeFc->SetLineColor(kGreen+2);
  gFcConservativeFc->SetFillStyle(3004);
  gFcConservativeFc->SetFillColor(kGreen);
  gFcConservativeFc->Draw("2P");
  //
  gFcConservativeNb ->SetMarkerStyle(25);
  gFcConservativeNb ->SetMarkerSize(1.3);
  gFcConservativeNb->SetMarkerColor(kOrange+7);//kViolet+5);
  gFcConservativeNb->SetLineColor(kOrange+7);//kViolet+5);
  gFcConservativeNb->SetFillStyle(3018);
  gFcConservativeNb->SetFillColor(kMagenta);
  gFcConservativeNb->Draw("2P");
  //
  gFcCorrConservative->SetMarkerStyle(21);
  gFcCorrConservative->SetLineColor(kRed);
  gFcCorrConservative->SetLineWidth(2);
  gFcCorrConservative->SetFillColor(kRed);
  gFcCorrConservative->SetFillStyle(0);
  gFcCorrConservative->Draw("2P");
  //
  leg=combinefcunc->BuildLegend();
  leg->SetFillStyle(0);
  //
  combinefcunc->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");
  }
  else if ( decay==4 ) {
    thDsKKpiprediction->SetLineColor(kGreen+2);
    thDsKKpiprediction->SetLineWidth(3);
    thDsKKpiprediction->SetFillColor(kGreen-6);
    thDsKKpiprediction->Draw("3CA");
    thDsKKpiprediction->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");
  //
  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");
  else if ( decay==4 ) leg->AddEntry(thDsKKpiprediction,"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();
  gFcCorrConservative->Write();
  out->Write();

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