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

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include "TH1D.h"
#include "TH1.h"
#include "TH2F.h"
#include "TNtuple.h"
#include "TFile.h"
#include "TGraphAsymmErrors.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TMath.h"
#endif

/******************************************************
 *
 *  Macro to compute fprompt for pPb and PbPb data
 *    using the outputs of the HFPtSpectrumRaa maco
 *  Questions and complains to Z.Conesa del Valle
 *
 *******************************************************/

Bool_t elossFDQuadSum = true;

TGraphAsymmErrors * DrawFpromptVsRaaElossHypo(const char *infile="", //const char *outfile="",
					      Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t CentralHypo = 1.0, Bool_t isYieldUnc=true);


//_______________________________________________________________________
void DrawFpromptVsRaaElossHypoCombined(const char *infileNb="", const char *infilefc="", const char *outfile="", Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t CentralHypo = 1.0, Bool_t isYieldUnc=true)
{

  TH2F *hdraw = new TH2F("hdraw","fprompt vs pt; p_{T} (GeV/c); f_{prompt}",36,0.,36.,10.,0.,1);
  hdraw->Draw();

  cout<<endl<<endl<<" ********* Checking Nb method ********* "<<endl<<endl;
  TGraphAsymmErrors *gfPromptNb = DrawFpromptVsRaaElossHypo(infileNb,MinHypo,MaxHypo,CentralHypo,isYieldUnc);
  gfPromptNb->SetNameTitle("gfPromptNb","fprompt Nb method, Eloss range considered");
  gfPromptNb->SetLineColor(kBlack);
  gfPromptNb->SetMarkerColor(kGreen+2);
  gfPromptNb->SetMarkerStyle(25);
  gfPromptNb->Draw("P");
					     
  cout<<endl<<endl<<" ********* Checking fc method ********* "<<endl<<endl;
  TGraphAsymmErrors *gfPromptfc = DrawFpromptVsRaaElossHypo(infilefc,MinHypo,MaxHypo,CentralHypo,isYieldUnc);
  gfPromptfc->SetNameTitle("gfPromptfc","fprompt fc method, Eloss range considered");
  gfPromptfc->SetLineColor(kBlack);
  gfPromptfc->SetMarkerColor(kRed);
  gfPromptfc->SetMarkerStyle(21);
  gfPromptfc->Draw("P");

  cout<<endl<<endl<<" ********* Computing the envelope ******** "<<endl<<endl;
  TGraphAsymmErrors *gFPromptCombined = new TGraphAsymmErrors(0);
  gFPromptCombined->SetNameTitle("gFPromptCombined","fprompt Nb + fc, Eloss range considered");

  Int_t nbins=gfPromptNb->GetN();
  for(Int_t i=0; i<nbins; i++) {
    Double_t pt=0.,value=0.,ptfc=0.,valuefc=0.,exlow=0.,eylow=0.,eyNblow=0.,eyfclow=0.,eyhigh=0.,eyNbhigh=0.,eyfchigh=0.;
    gfPromptNb->GetPoint(i,pt,value);
    gfPromptfc->GetPoint(i,ptfc,valuefc);
    exlow = gfPromptNb->GetErrorXlow(i);
    eyNblow = gfPromptNb->GetErrorYlow(i);
    eyNbhigh = gfPromptNb->GetErrorYhigh(i);
    eyfclow = gfPromptfc->GetErrorYlow(i);
    eyfchigh = gfPromptfc->GetErrorYhigh(i);
    Double_t uncertainties[5]={ value, value-eyNblow, value+eyNbhigh, valuefc-eyfclow, valuefc+eyfchigh };
    eylow = value - TMath::MinElement(5,uncertainties);
    eyhigh = TMath::MaxElement(5,uncertainties) - value;
    gFPromptCombined->SetPoint(i,pt,value);
    gFPromptCombined->SetPointError(i,exlow,exlow,eylow,eyhigh);
  }
  
  gFPromptCombined->SetFillStyle(0);
  gFPromptCombined->SetLineColor(kBlack);
  gFPromptCombined->SetMarkerColor(kBlack);
  gFPromptCombined->SetMarkerStyle(20);
  gFPromptCombined->Draw("2P");

  cout<<endl<<endl;

  TFile* fout= new TFile(outfile,"recreate");
  gfPromptNb->Write();
  gfPromptfc->Write();
  gFPromptCombined->Write();
  fout->Write();
}

//_______________________________________________________________________
TGraphAsymmErrors * DrawFpromptVsRaaElossHypo(const char *infile, 
					      //const char *outfile="",
					      Double_t MinHypo, Double_t MaxHypo, Double_t CentralHypo, Bool_t isYieldUnc)
{

  Float_t pt=0.,TAB=0.,sigmaPP=0.,invyieldAB=0.,RaaCharm=0.,RaaBeauty=0., fc=0.;
  Float_t yieldFdHigh=0., yieldFdLow=0.,RaaCharmFdHigh=0.,RaaCharmFdLow=0.;

  TFile *fRaa = new TFile(infile,"read");
  TH1D *hRABvsPt = (TH1D*)fRaa->Get("hRABvsPt");
  TNtuple *ntRaa = (TNtuple*)fRaa->Get("ntupleRAB");
  ntRaa->SetBranchAddress("pt",&pt);
  ntRaa->SetBranchAddress("TAB",&TAB);
  ntRaa->SetBranchAddress("sigmaPP",&sigmaPP);
  ntRaa->SetBranchAddress("invyieldAB",&invyieldAB);
  ntRaa->SetBranchAddress("RABCharm",&RaaCharm);
  ntRaa->SetBranchAddress("RABBeauty",&RaaBeauty);
  ntRaa->SetBranchAddress("RABCharmFDHigh",&RaaCharmFdHigh);
  ntRaa->SetBranchAddress("RABCharmFDLow",&RaaCharmFdLow);
  ntRaa->SetBranchAddress("invyieldABFDHigh",&yieldFdHigh);
  ntRaa->SetBranchAddress("invyieldABFDLow",&yieldFdLow);
  ntRaa->SetBranchAddress("fc",&fc);

  // define the binning
  Int_t nbins = hRABvsPt->GetNbinsX();

  //
  //
  // Search the central value of the energy loss hypothesis Rb = Rc (bin)
  //
  Double_t ElossCentral[nbins+1];
  Double_t fcV[nbins+1], fcVmax[nbins+1], fcVmin[nbins+1];
  Double_t fdMax[nbins+1], fdMin[nbins+1];
  for(Int_t i=0; i<=nbins; i++) { 
    ElossCentral[i]=0.; 
    fcV[i]=0.; fcVmax[i]=0.; fcVmin[i]=6.;
    fdMax[i]=0.; fdMin[i]=6.;
  }
  //
  Int_t netries = ntRaa->GetEntries();
  //
  for(Int_t ientry=0; ientry<=netries; ientry++){
    ntRaa->GetEntry(ientry);
    Double_t ElossHypo =  1. / (RaaCharm / RaaBeauty) ;
    //
    // Find the bin for the central Eloss hypo
    //
    if( TMath::Abs( ElossHypo - CentralHypo ) < 0.075 ){
      Int_t hABbin = hRABvsPt->FindBin( pt );
      Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - CentralHypo );
      Double_t DeltaV = TMath::Abs( ElossHypo - CentralHypo );
      //      cout << " pt " << ptTot << " ECentral Tot " << ElossCentralTot[ hABbin ] << " Ehypo "<< ElossHypo ;
      if ( DeltaV < DeltaIni ) ElossCentral[ hABbin ] = ElossHypo;
      //      cout << " final ECentral " << ElossCentralTot[ hABbin ] << endl;
    }
  }

  //
  for(Int_t ientry=0; ientry<=netries; ientry++){

    ntRaa->GetEntry(ientry);
    Double_t ElossHypo =  1. / (RaaCharm / RaaBeauty) ;
    Int_t hABbin = hRABvsPt->FindBin( pt );

    // Get the central value of Raa
    if ( ElossHypo == ElossCentral[ hABbin ] ) {

      fcV[hABbin] = fc;

      // Now look at FONLL band for the central value
      Double_t elems[3]= {0., 0., 0.};
      // 	fdMax[hABbin]= (RaaCharmFdHigh-RaaCharm)/RaaCharm;
      // 	fdMin[hABbin]= (RaaCharm-RaaCharmFdLow)/RaaCharm;
      if(isYieldUnc) {
	elems[0]=invyieldAB;
	elems[1]=yieldFdHigh;
	elems[2]=yieldFdLow;
      } else {
	elems[0]=RaaCharm;
	elems[1]=RaaCharmFdHigh;
	elems[2]=RaaCharmFdLow;
      }
      fdMax[hABbin]= (TMath::MaxElement(3,elems)-elems[0])/elems[0];
      fdMin[hABbin]= (elems[0]-TMath::MinElement(3,elems))/elems[0];
      //      cout<<" Check FONLL band giving for pt="<< pt <<" +"<< fdMax[hABbin]<< " -"<< fdMin[hABbin]<<" (relative variation)"<<endl;
    }

    // Here check the Eloss band
    if( RaaCharm>0 && ElossHypo >= MinHypo && ElossHypo <=MaxHypo ) {
      if( fc < fcVmin[hABbin] ) fcVmin[hABbin] = fc;
      if( fc > fcVmax[hABbin] ) fcVmax[hABbin] = fc;
    }
  }

  TGraphAsymmErrors *gFPrompt = new TGraphAsymmErrors(0);
  Int_t j=0;
  for(Int_t i=1; i<=nbins; i++) {
    Double_t xpt = hRABvsPt->GetBinCenter(i);
    Double_t width  = hRABvsPt->GetBinWidth(i)/2;
    Double_t elossLow = (fcV[i]-fcVmin[i])/fcV[i];
    Double_t elossHigh = (fcVmax[i]-fcV[i])/fcV[i];
    Double_t fcLow = fcV[i]-fcVmin[i];
    Double_t fcHigh = fcVmax[i]-fcV[i];
    //    cout<<" pt="<<xpt<<" , fprompt="<<fcV[i]<<" ["<< fcVmin[i] <<","<<fcVmax[i] <<"] ( Eloss) "<<endl;

    if(elossFDQuadSum){ // add quadratically
      fcLow = TMath::Sqrt( fdMin[i]*fdMin[i] + elossLow*elossLow )*fcV[i];
      fcHigh = TMath::Sqrt( fdMax[i]*fdMax[i] + elossHigh*elossHigh )*fcV[i];
    } else { // add linearly
      fcLow = ( fdMin[i] + elossLow )*fcV[i];
      fcHigh = ( fdMax[i] + elossHigh )*fcV[i];
    }

    gFPrompt->SetPoint(j,xpt,fcV[i]);
    //    gFPrompt->SetPointError(j,width,width,fcV[i]-fcVmin[i],fcVmax[i]-fcV[i]);
    gFPrompt->SetPointError(j,width,width,fcLow,fcHigh);
    j++;

    cout<<" pt="<<xpt<<" , fprompt="<<fcV[i]<<" +"<<fcHigh<<" -"<<fcLow<<" ( splitted it is +"<<elossLow*fcV[i]<<" -"<<elossHigh*fcV[i]<<" Eloss, +"<<fdMax[i]*fcV[i]<<" -"<<fdMin[i]*fcV[i]<<" FONLL-band )"<<endl;
  }
  //  gFPrompt->Draw("AP");

  // if(strcmp(outfile,"")!=0) {
  //   TFile *out = new TFile(outfile,"recreate");
  //   gFPrompt->Write();
  // }

  fRaa->Close();

  return gFPrompt;
  
}
Commit	Line	Data
4988649c	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <Riostream.h>
	3	#include "TH1D.h"
	4	#include "TH1.h"
	5	#include "TH2F.h"
	6	#include "TNtuple.h"
	7	#include "TFile.h"
	8	#include "TGraphAsymmErrors.h"
	9	#include "TCanvas.h"
	10	#include "TROOT.h"
	11	#include "TStyle.h"
	12	#include "TLegend.h"
	13	#include "TMath.h"
	14	#endif
	15
	16	/******************************************************
	17	*
	18	* Macro to compute fprompt for pPb and PbPb data
	19	* using the outputs of the HFPtSpectrumRaa maco
	20	* Questions and complains to Z.Conesa del Valle
	21	*
	22	*******************************************************/
	23
	24	Bool_t elossFDQuadSum = true;
	25
	26	TGraphAsymmErrors * DrawFpromptVsRaaElossHypo(const char infile="", //const char outfile="",
	27	Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t CentralHypo = 1.0, Bool_t isYieldUnc=true);
	28
	29
	30	//_______________________________________________________________________
	31	void DrawFpromptVsRaaElossHypoCombined(const char infileNb="", const char infilefc="", const char *outfile="", Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t CentralHypo = 1.0, Bool_t isYieldUnc=true)
	32	{
	33
	34	TH2F *hdraw = new TH2F("hdraw","fprompt vs pt; p_{T} (GeV/c); f_{prompt}",36,0.,36.,10.,0.,1);
	35	hdraw->Draw();
	36
	37	cout<<endl<<endl<<" ******* Checking Nb method ******* "<<endl<<endl;
	38	TGraphAsymmErrors *gfPromptNb = DrawFpromptVsRaaElossHypo(infileNb,MinHypo,MaxHypo,CentralHypo,isYieldUnc);
	39	gfPromptNb->SetNameTitle("gfPromptNb","fprompt Nb method, Eloss range considered");
	40	gfPromptNb->SetLineColor(kBlack);
	41	gfPromptNb->SetMarkerColor(kGreen+2);
	42	gfPromptNb->SetMarkerStyle(25);
	43	gfPromptNb->Draw("P");
	44
	45	cout<<endl<<endl<<" ******* Checking fc method ******* "<<endl<<endl;
	46	TGraphAsymmErrors *gfPromptfc = DrawFpromptVsRaaElossHypo(infilefc,MinHypo,MaxHypo,CentralHypo,isYieldUnc);
	47	gfPromptfc->SetNameTitle("gfPromptfc","fprompt fc method, Eloss range considered");
	48	gfPromptfc->SetLineColor(kBlack);
	49	gfPromptfc->SetMarkerColor(kRed);
	50	gfPromptfc->SetMarkerStyle(21);
	51	gfPromptfc->Draw("P");
	52
	53	cout<<endl<<endl<<" ******* Computing the envelope ****** "<<endl<<endl;
	54	TGraphAsymmErrors *gFPromptCombined = new TGraphAsymmErrors(0);
	55	gFPromptCombined->SetNameTitle("gFPromptCombined","fprompt Nb + fc, Eloss range considered");
	56
	57	Int_t nbins=gfPromptNb->GetN();
	58	for(Int_t i=0; i<nbins; i++) {
	59	Double_t pt=0.,value=0.,ptfc=0.,valuefc=0.,exlow=0.,eylow=0.,eyNblow=0.,eyfclow=0.,eyhigh=0.,eyNbhigh=0.,eyfchigh=0.;
	60	gfPromptNb->GetPoint(i,pt,value);
	61	gfPromptfc->GetPoint(i,ptfc,valuefc);
	62	exlow = gfPromptNb->GetErrorXlow(i);
	63	eyNblow = gfPromptNb->GetErrorYlow(i);
	64	eyNbhigh = gfPromptNb->GetErrorYhigh(i);
65	eyfclow = gfPromptfc->GetErrorYlow(i);
66	eyfchigh = gfPromptfc->GetErrorYhigh(i);
67	Double_t uncertainties[5]={ value, value-eyNblow, value+eyNbhigh, valuefc-eyfclow, valuefc+eyfchigh };
68	eylow = value - TMath::MinElement(5,uncertainties);
69	eyhigh = TMath::MaxElement(5,uncertainties) - value;
70	gFPromptCombined->SetPoint(i,pt,value);
71	gFPromptCombined->SetPointError(i,exlow,exlow,eylow,eyhigh);
72	}
73
74	gFPromptCombined->SetFillStyle(0);
75	gFPromptCombined->SetLineColor(kBlack);
76	gFPromptCombined->SetMarkerColor(kBlack);
77	gFPromptCombined->SetMarkerStyle(20);
78	gFPromptCombined->Draw("2P");
79
80	cout<<endl<<endl;
81
82	TFile* fout= new TFile(outfile,"recreate");
83	gfPromptNb->Write();
84	gfPromptfc->Write();
85	gFPromptCombined->Write();
86	fout->Write();
87	}
88
89	//_______________________________________________________________________
90	TGraphAsymmErrors * DrawFpromptVsRaaElossHypo(const char *infile,
91	//const char *outfile="",
92	Double_t MinHypo, Double_t MaxHypo, Double_t CentralHypo, Bool_t isYieldUnc)
93	{
94
95	Float_t pt=0.,TAB=0.,sigmaPP=0.,invyieldAB=0.,RaaCharm=0.,RaaBeauty=0., fc=0.;
96	Float_t yieldFdHigh=0., yieldFdLow=0.,RaaCharmFdHigh=0.,RaaCharmFdLow=0.;
97
98	TFile *fRaa = new TFile(infile,"read");
99	TH1D hRABvsPt = (TH1D)fRaa->Get("hRABvsPt");
100	TNtuple ntRaa = (TNtuple)fRaa->Get("ntupleRAB");
101	ntRaa->SetBranchAddress("pt",&pt);
102	ntRaa->SetBranchAddress("TAB",&TAB);
103	ntRaa->SetBranchAddress("sigmaPP",&sigmaPP);
104	ntRaa->SetBranchAddress("invyieldAB",&invyieldAB);
105	ntRaa->SetBranchAddress("RABCharm",&RaaCharm);
106	ntRaa->SetBranchAddress("RABBeauty",&RaaBeauty);
107	ntRaa->SetBranchAddress("RABCharmFDHigh",&RaaCharmFdHigh);
108	ntRaa->SetBranchAddress("RABCharmFDLow",&RaaCharmFdLow);
109	ntRaa->SetBranchAddress("invyieldABFDHigh",&yieldFdHigh);
110	ntRaa->SetBranchAddress("invyieldABFDLow",&yieldFdLow);
111	ntRaa->SetBranchAddress("fc",&fc);
112
113	// define the binning
114	Int_t nbins = hRABvsPt->GetNbinsX();
115
116	//
117	//
118	// Search the central value of the energy loss hypothesis Rb = Rc (bin)
119	//
120	Double_t ElossCentral[nbins+1];
121	Double_t fcV[nbins+1], fcVmax[nbins+1], fcVmin[nbins+1];
122	Double_t fdMax[nbins+1], fdMin[nbins+1];
123	for(Int_t i=0; i<=nbins; i++) {
124	ElossCentral[i]=0.;
125	fcV[i]=0.; fcVmax[i]=0.; fcVmin[i]=6.;
126	fdMax[i]=0.; fdMin[i]=6.;
127	}
128	//
129	Int_t netries = ntRaa->GetEntries();
130	//
131	for(Int_t ientry=0; ientry<=netries; ientry++){
132	ntRaa->GetEntry(ientry);
133	Double_t ElossHypo = 1. / (RaaCharm / RaaBeauty) ;
134	//
135	// Find the bin for the central Eloss hypo
136	//
137	if( TMath::Abs( ElossHypo - CentralHypo ) < 0.075 ){
138	Int_t hABbin = hRABvsPt->FindBin( pt );
139	Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - CentralHypo );
140	Double_t DeltaV = TMath::Abs( ElossHypo - CentralHypo );
141	// cout << " pt " << ptTot << " ECentral Tot " << ElossCentralTot[ hABbin ] << " Ehypo "<< ElossHypo ;
142	if ( DeltaV < DeltaIni ) ElossCentral[ hABbin ] = ElossHypo;
143	// cout << " final ECentral " << ElossCentralTot[ hABbin ] << endl;
144	}
145	}
146
147	//
148	for(Int_t ientry=0; ientry<=netries; ientry++){
149
150	ntRaa->GetEntry(ientry);
151	Double_t ElossHypo = 1. / (RaaCharm / RaaBeauty) ;
152	Int_t hABbin = hRABvsPt->FindBin( pt );
153
154	// Get the central value of Raa
155	if ( ElossHypo == ElossCentral[ hABbin ] ) {
156
157	fcV[hABbin] = fc;
158
159	// Now look at FONLL band for the central value
160	Double_t elems[3]= {0., 0., 0.};
161	// fdMax[hABbin]= (RaaCharmFdHigh-RaaCharm)/RaaCharm;
162	// fdMin[hABbin]= (RaaCharm-RaaCharmFdLow)/RaaCharm;
163	if(isYieldUnc) {
164	elems[0]=invyieldAB;
165	elems[1]=yieldFdHigh;
166	elems[2]=yieldFdLow;
167	} else {
168	elems[0]=RaaCharm;
169	elems[1]=RaaCharmFdHigh;
170	elems[2]=RaaCharmFdLow;
171	}
172	fdMax[hABbin]= (TMath::MaxElement(3,elems)-elems[0])/elems[0];
173	fdMin[hABbin]= (elems[0]-TMath::MinElement(3,elems))/elems[0];
174	// cout<<" Check FONLL band giving for pt="<< pt <<" +"<< fdMax[hABbin]<< " -"<< fdMin[hABbin]<<" (relative variation)"<<endl;
175	}
176
177	// Here check the Eloss band
178	if( RaaCharm>0 && ElossHypo >= MinHypo && ElossHypo <=MaxHypo ) {
179	if( fc < fcVmin[hABbin] ) fcVmin[hABbin] = fc;
180	if( fc > fcVmax[hABbin] ) fcVmax[hABbin] = fc;
181	}
182	}
183
184	TGraphAsymmErrors *gFPrompt = new TGraphAsymmErrors(0);
185	Int_t j=0;
186	for(Int_t i=1; i<=nbins; i++) {
187	Double_t xpt = hRABvsPt->GetBinCenter(i);
188	Double_t width = hRABvsPt->GetBinWidth(i)/2;
189	Double_t elossLow = (fcV[i]-fcVmin[i])/fcV[i];
190	Double_t elossHigh = (fcVmax[i]-fcV[i])/fcV[i];
191	Double_t fcLow = fcV[i]-fcVmin[i];
192	Double_t fcHigh = fcVmax[i]-fcV[i];
193	// cout<<" pt="<<xpt<<" , fprompt="<<fcV[i]<<" ["<< fcVmin[i] <<","<<fcVmax[i] <<"] ( Eloss) "<<endl;
194
195	if(elossFDQuadSum){ // add quadratically
196	fcLow = TMath::Sqrt( fdMin[i]fdMin[i] + elossLowelossLow )*fcV[i];
197	fcHigh = TMath::Sqrt( fdMax[i]fdMax[i] + elossHighelossHigh )*fcV[i];
198	} else { // add linearly
199	fcLow = ( fdMin[i] + elossLow )*fcV[i];
200	fcHigh = ( fdMax[i] + elossHigh )*fcV[i];
201	}
202
203	gFPrompt->SetPoint(j,xpt,fcV[i]);
204	// gFPrompt->SetPointError(j,width,width,fcV[i]-fcVmin[i],fcVmax[i]-fcV[i]);
205	gFPrompt->SetPointError(j,width,width,fcLow,fcHigh);
206	j++;
207
208	cout<<" pt="<<xpt<<" , fprompt="<<fcV[i]<<" +"<<fcHigh<<" -"<<fcLow<<" ( splitted it is +"<<elossLowfcV[i]<<" -"<<elossHighfcV[i]<<" Eloss, +"<<fdMax[i]fcV[i]<<" -"<<fdMin[i]fcV[i]<<" FONLL-band )"<<endl;
209	}
210	// gFPrompt->Draw("AP");
211
212	// if(strcmp(outfile,"")!=0) {
213	// TFile *out = new TFile(outfile,"recreate");
214	// gFPrompt->Write();
215	// }
216
217	fRaa->Close();
218
219	return gFPrompt;
220
221	}