//
// (the corrected yields per bin are divided by the bin-width)
//
+//
+// In HIC you can also evaluate how the feed-down correction is influenced by an energy loss hypothesis:
+// Raa(c-->D) / Raa(b-->D) defined here as Rcb for the "fc" method
+// Raa(b-->D) defined here as Rb for the "Nb" method
+//
// Author: Z.Conesa, zconesa@in2p3.fr
//***********************************************************************
#include "TMath.h"
#include "TH1.h"
#include "TH1D.h"
+#include "TH2.h"
+#include "TH2D.h"
+#include "TNtuple.h"
#include "TGraphAsymmErrors.h"
#include "TNamed.h"
#include "TCanvas.h"
fhFeedDownEffpt(NULL),
fhRECpt(NULL),
fgRECSystematics(NULL),
+ fNevts(1),
fLuminosity(),
fTrigEfficiency(),
fGlobalEfficiencyUncertainties(),
+ fTab(),
fhFc(NULL),
fhFcMax(NULL),
fhFcMin(NULL),
+ fhFcRcb(NULL),
fgFcExtreme(NULL),
fgFcConservative(NULL),
fhYieldCorr(NULL),
fhYieldCorrMax(NULL),
fhYieldCorrMin(NULL),
+ fhYieldCorrRcb(NULL),
fgYieldCorr(NULL),
fgYieldCorrExtreme(NULL),
fgYieldCorrConservative(NULL),
fhSigmaCorr(NULL),
fhSigmaCorrMax(NULL),
fhSigmaCorrMin(NULL),
+ fhSigmaCorrDataSyst(NULL),
+ fhSigmaCorrRcb(NULL),
fgSigmaCorr(NULL),
fgSigmaCorrExtreme(NULL),
fgSigmaCorrConservative(NULL),
+ fnSigma(NULL),
fFeedDownOption(option),
- fAsymUncertainties(kTRUE)
+ fAsymUncertainties(kTRUE),
+ fPbPbElossHypothesis(kFALSE),
+ fhStatUncEffcSigma(NULL),
+ fhStatUncEffbSigma(NULL),
+ fhStatUncEffcFD(NULL),
+ fhStatUncEffbFD(NULL)
{
//
// Default constructor
fLuminosity[0]=1.; fLuminosity[1]=0.;
fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
+ fTab[0]=1.; fTab[1]=0.;
}
fhFeedDownEffpt(rhs.fhFeedDownEffpt),
fhRECpt(rhs.fhRECpt),
fgRECSystematics(rhs.fgRECSystematics),
+ fNevts(rhs.fNevts),
fLuminosity(),
fTrigEfficiency(),
fGlobalEfficiencyUncertainties(),
+ fTab(),
fhFc(rhs.fhFc),
fhFcMax(rhs.fhFcMax),
fhFcMin(rhs.fhFcMin),
+ fhFcRcb(rhs.fhFcRcb),
fgFcExtreme(rhs.fgFcExtreme),
fgFcConservative(rhs.fgFcConservative),
fhYieldCorr(rhs.fhYieldCorr),
fhYieldCorrMax(rhs.fhYieldCorrMax),
fhYieldCorrMin(rhs.fhYieldCorrMin),
+ fhYieldCorrRcb(rhs.fhYieldCorrRcb),
fgYieldCorr(rhs.fgYieldCorr),
fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
fgYieldCorrConservative(rhs.fgYieldCorrConservative),
fhSigmaCorr(rhs.fhSigmaCorr),
fhSigmaCorrMax(rhs.fhSigmaCorrMax),
fhSigmaCorrMin(rhs.fhSigmaCorrMin),
+ fhSigmaCorrDataSyst(rhs.fhSigmaCorrDataSyst),
+ fhSigmaCorrRcb(rhs.fhSigmaCorrRcb),
fgSigmaCorr(rhs.fgSigmaCorr),
fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
+ fnSigma(rhs.fnSigma),
fFeedDownOption(rhs.fFeedDownOption),
- fAsymUncertainties(rhs.fAsymUncertainties)
+ fAsymUncertainties(rhs.fAsymUncertainties),
+ fPbPbElossHypothesis(rhs.fPbPbElossHypothesis),
+ fhStatUncEffcSigma(NULL),
+ fhStatUncEffbSigma(NULL),
+ fhStatUncEffcFD(NULL),
+ fhStatUncEffbFD(NULL)
{
//
// Copy constructor
fLuminosity[i] = rhs.fLuminosity[i];
fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
+ fTab[i] = rhs.fTab[i];
}
}
fhFeedDownEffpt = source.fhFeedDownEffpt;
fhRECpt = source.fhRECpt;
fgRECSystematics = source.fgRECSystematics;
+ fNevts = source.fNevts;
fhFc = source.fhFc;
fhFcMax = source.fhFcMax;
fhFcMin = source.fhFcMin;
+ fhFcRcb = source.fhFcRcb;
fgFcExtreme = source.fgFcExtreme;
fgFcConservative = source.fgFcConservative;
fhYieldCorr = source.fhYieldCorr;
fhYieldCorrMax = source.fhYieldCorrMax;
fhYieldCorrMin = source.fhYieldCorrMin;
+ fhYieldCorrRcb = source.fhYieldCorrRcb;
fgYieldCorr = source.fgYieldCorr;
fgYieldCorrExtreme = source.fgYieldCorrExtreme;
fgYieldCorrConservative = source.fgYieldCorrConservative;
fhSigmaCorr = source.fhSigmaCorr;
fhSigmaCorrMax = source.fhSigmaCorrMax;
fhSigmaCorrMin = source.fhSigmaCorrMin;
+ fhSigmaCorrDataSyst = source.fhSigmaCorrDataSyst;
+ fhSigmaCorrRcb = source.fhSigmaCorrRcb;
fgSigmaCorr = source.fgSigmaCorr;
fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
fgSigmaCorrConservative = source.fgSigmaCorrConservative;
+ fnSigma = source.fnSigma;
fFeedDownOption = source.fFeedDownOption;
fAsymUncertainties = source.fAsymUncertainties;
+ fPbPbElossHypothesis = source.fPbPbElossHypothesis;
for(Int_t i=0; i<2; i++){
fLuminosity[i] = source.fLuminosity[i];
fTrigEfficiency[i] = source.fTrigEfficiency[i];
fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
+ fTab[i] = source.fTab[i];
}
return *this;
if (fhFc) delete fhFc;
if (fhFcMax) delete fhFcMax;
if (fhFcMin) delete fhFcMin;
+ if (fhFcRcb) delete fhFcRcb;
if (fgFcExtreme) delete fgFcExtreme;
if (fgFcConservative) delete fgFcConservative;
if (fhYieldCorr) delete fhYieldCorr;
if (fhYieldCorrMax) delete fhYieldCorrMax;
- if (fhYieldCorrMin) delete fhYieldCorrMin;
+ if (fhYieldCorrMin) delete fhYieldCorrMin;
+ if (fhYieldCorrRcb) delete fhYieldCorrRcb;
if (fgYieldCorr) delete fgYieldCorr;
if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
if (fgYieldCorrConservative) delete fgYieldCorrConservative;
if (fhSigmaCorr) delete fhSigmaCorr;
if (fhSigmaCorrMax) delete fhSigmaCorrMax;
- if (fhSigmaCorrMin) delete fhSigmaCorrMin;
+ if (fhSigmaCorrMin) delete fhSigmaCorrMin;
+ if (fhSigmaCorrDataSyst) delete fhSigmaCorrDataSyst;
if (fgSigmaCorr) delete fgSigmaCorr;
if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
+ if (fnSigma) delete fnSigma;
}
// Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
// (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 )
// (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
+ //
+ // In HIC the feed-down correction varies with an energy loss hypothesis:
+ // Raa(c-->D) / Raa(b-->D) for the "fc" method, Raa(b-->D) for the "Nb" method (see exact formulas in the functions)
+ //
//
// First: Initialization
}
// Print out information
- // printf("\n\n Correcting the spectra with : \n luminosity = %2.2e +- %2.2e, trigger efficiency = %2.2e +- %2.2e, \n delta_y = %2.2f, BR_c = %2.2e, BR_b_decay = %2.2e \n\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
printf("\n\n Correcting the spectra with : \n luminosity = %2.2e +- %2.2e, trigger efficiency = %2.2e +- %2.2e, \n delta_y = %2.2f, BR_c = %2.2e, BR_b_decay = %2.2e \n %2.2f percent uncertainty on the efficiencies, and %2.2f percent uncertainty on the b/c efficiencies ratio \n\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay,fGlobalEfficiencyUncertainties[0],fGlobalEfficiencyUncertainties[1]);
+ if (fPbPbElossHypothesis) printf("\n\n The considered Tab is %4.2e +- %2.2e \n\n",fTab[0],fTab[1]);
//
// Finally: Correct from yields to cross-section
fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
+ fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",nbins,limits);
+ if (fPbPbElossHypothesis && fFeedDownOption==1) {
+ fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
+ fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma");
+ }
+ if (fPbPbElossHypothesis && fFeedDownOption==2) {
+ fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
+ fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma");
+ }
// and the output TGraphAsymmErrors
if (fAsymUncertainties){
fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
}
+ fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",nbins,limits);
+ fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",nbins,limits);
+
// protect against null denominator
if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
+ Double_t errvalueStatUncEffc=0., errvalueStatUncEffb=0.;
for(Int_t ibin=1; ibin<=nbins; ibin++){
// Sigma calculation
errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
+
+ // stat unc of the efficiencies, separately
+ errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
+ errvalueStatUncEffb = 0.;
+
}
else {
// protect against null denominator
errvalueMin = errvalueMax;
}
+ //
// Fill the histograms
+ //
fhSigmaCorr->SetBinContent(ibin,value);
fhSigmaCorr->SetBinError(ibin,errValue);
+ //
+ // Fill the histos and ntuple vs the Eloss hypothesis
+ //
+ if (fPbPbElossHypothesis) {
+ // Loop over the Eloss hypothesis
+ for (Float_t rval=0.0025; rval<4.0; rval+=0.005) {
+ Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
+ Double_t yieldRcbvalue = fhYieldCorrRcb->GetBinContent(ibin,rbin);
+ // Sigma calculation
+ // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
+ Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
+ ( yieldRcbvalue / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
+ : 0. ;
+ fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , rval, sigmaRcbvalue );
+ // if(ibin==3)
+ // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
+ fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
+ rval, fhFcRcb->GetBinContent(ibin,rbin),
+ yieldRcbvalue, sigmaRcbvalue );
+ }
+ }
+ //
// Fill the TGraphAsymmErrors
if (fAsymUncertainties) {
Double_t x = fhYieldCorr->GetBinCenter(ibin);
fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
+ fhStatUncEffcSigma->SetBinContent(ibin,0.);
+ if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
+ fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
+ // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
}
}
// (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
// systematic uncertainty on the acceptance x efficiency b/c ratio are included
//
+ // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
+ // fc (Rcb) = ( 1. / ( 1 + (eff_b/eff_c)*(N_b/N_c)* (1/Rcb) ) );
+ //
// define the variables
Int_t nbins = fhRECpt->GetNbinsX();
fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
+ if(fPbPbElossHypothesis) fhFcRcb = new TH2D("fhFcRcb","fc correction factor vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; fc correction",nbins,limits,800,0.,4.);
// two local control histograms
TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
}
+ fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
+ fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
+ Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
+ Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
//
// Compute fc
// fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
- // delta(eff_b/eff_c) is a percentage = fGlobalEfficiencyUncertainties[1]*effRatio
- // correctionUnc = correction*correction * theoryRatio * fGlobalEfficiencyUncertainties[1]*effRatio;
- // correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * fGlobalEfficiencyUncertainties[1]*effRatio;
- // correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * fGlobalEfficiencyUncertainties[1]*effRatio;
- // correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA * fGlobalEfficiencyUncertainties[1]*effRatio;
- // correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB * fGlobalEfficiencyUncertainties[1]*effRatio;
+ // delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
correctionUnc = correction*correction * theoryRatio * effRatio *
TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
);
+ //
+ correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
+ (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
+ correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
+ (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
+
correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
);
+ correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
+ (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
+ correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
+ (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
// Fill in the histograms
hTheoryRatio->SetBinContent(ibin,theoryRatio);
hEffRatio->SetBinContent(ibin,effRatio);
fhFc->SetBinContent(ibin,correction);
+ //
+ // Estimate how the result varies vs charm/beauty Eloss hypothesis
+ //
+ if ( TMath::Abs(correction-1.0)>0.01 && fPbPbElossHypothesis){
+ // Loop over the Eloss hypothesis
+ // Int_t rbin=0;
+ for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
+ Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
+ fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
+ // if(ibin==3){
+ // cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
+ // rbin++;
+ // }
+ }
+ }
+ //
+ // Fill the rest of (asymmetric) histograms
+ //
if (fAsymUncertainties) {
Double_t x = fhDirectMCpt->GetBinCenter(ibin);
Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
}
+
+ Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
+ correctionConservativeBUncStatEffc/correctionConservativeB };
+ Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
+ correctionConservativeBUncStatEffb/correctionConservativeB };
+ Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
+ Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
+ fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
+ fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
+ // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
+
}
}
// (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
//
// ( Calculation done bin by bin )
+ //
+ // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
if (!fhFc || !fhRECpt) {
AliError(" Reconstructed or fc distributions are not defined");
// declare the output histograms
fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
- fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
+ fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
+ if(fPbPbElossHypothesis) fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by fc) vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; corrected yield",nbins,limits,800,0.,4.);
// and the output TGraphAsymmErrors
if (fAsymUncertainties){
fgYieldCorr = new TGraphAsymmErrors(nbins+1);
}
else { errvalueMax = 0.; errvalueMin = 0.; }
- // fill in the histograms
+ //
+ // Fill in the histograms
+ //
fhYieldCorr->SetBinContent(ibin,value);
fhYieldCorr->SetBinError(ibin,errvalue);
+ //
+ // Fill the histos and ntuple vs the Eloss hypothesis
+ //
+ if (fPbPbElossHypothesis) {
+ // Loop over the Eloss hypothesis
+ for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
+ Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
+ Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
+ // physics = reco * fcRcb / bin-width
+ Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
+ fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
+ Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
+ fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
+ // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
+ }
+ }
if (fAsymUncertainties) {
Double_t center = fhYieldCorr->GetBinCenter(ibin);
fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
// + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
// where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
//
+ // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(b-->D) = Rb
+ // physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
+ //
Int_t nbins = fhRECpt->GetNbinsX();
Double_t binwidth = fhRECpt->GetBinWidth(1);
fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
+ if(fPbPbElossHypothesis) {
+ fhFcRcb = new TH2D("fhFcRcb","fc correction factor (Nb method) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; fc correction",nbins,limits,800,0.,4.);
+ fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by Nb) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; corrected yield",nbins,limits,800,0.,4.);
+ }
// and the output TGraphAsymmErrors
if (fAsymUncertainties){
fgYieldCorr = new TGraphAsymmErrors(nbins+1);
// variables to define fc-conservative
double correction=0, correctionMax=0., correctionMin=0.;
+ fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
+ fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
+ Double_t correctionUncStatEffc=0.;
+ Double_t correctionUncStatEffb=0.;
+
+
//
// Do the calculation
//
// Calculate the value
// physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
+ // In HIC : physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
+ //
+ //
+ Double_t frac = 1.0, errfrac =0.;
+ if(fPbPbElossHypothesis) {
+ frac = fTab[0]*fNevts;
+ errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
+ } else {
+ frac = fLuminosity[0];
+ errfrac = fLuminosity[1];
+ }
+
value = ( fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0. &&
fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
- fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) )
+ fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) )
: 0. ;
value /= fhRECpt->GetBinWidth(ibin);
+ if (value<0.) value =0.;
// Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
// Correction (fc) : Estimate of the relative amount feed-down subtracted
// correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
- correction = 1 - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ) / fhRECpt->GetBinContent(ibin) ;
+ // in HIC: correction = [ 1 - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
+ correction = 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ) / fhRECpt->GetBinContent(ibin) ;
+ if (correction<0.) correction = 0.;
// Systematic uncertainties
// (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
// (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
// + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
// where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
- kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
-
+ kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
//
if (fAsymUncertainties) {
Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
// Feed-down systematics
// min value with the maximum Nb
- errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
+ errvalueExtremeMin = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) +
( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
) / fhRECpt->GetBinWidth(ibin);
// max value with the minimum Nb
- errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
+ errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) +
( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
// Correction systematics (fc)
// min value with the maximum Nb
- correctionMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
+ correctionMin = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) +
( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
) / fhRECpt->GetBinContent(ibin) ;
// max value with the minimum Nb
- correctionMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
+ correctionMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) +
( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
) / fhRECpt->GetBinContent(ibin) ;
+ //
+ correctionUncStatEffb = TMath::Sqrt( ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) )
+ ) / fhRECpt->GetBinContent(ibin) ;
+ correctionUncStatEffc = 0.;
}
else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
- errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
+ errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
// fill in histograms
fhYieldCorr->SetBinContent(ibin,value);
- fhYieldCorr->SetBinError(ibin,errvalue);
+ fhYieldCorr->SetBinError(ibin,errvalue);
+ //
+ // Estimate how the result varies vs charm/beauty Eloss hypothesis
+ //
+ if ( correction>0.0001 && fPbPbElossHypothesis){
+ // Loop over the Eloss hypothesis
+ // Int_t rbin=0;
+ for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
+ // correction = [ 1 - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)* (rval) /reco ]
+ Double_t fcRcbvalue = 1 - (fTab[0]*fNevts*deltaY*branchingRatioBintoFinalDecay*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
+ if(fcRcbvalue<0.) fcRcbvalue=0.;
+ fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
+ // physics = reco * fcRcb / bin-width
+ Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
+ fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
+ Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
+ fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
+ // if(ibin==3){
+ // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
+ // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
+ // rbin++;
+ // }
+ }
+ }
+ //
+ // Fill the rest of (asymmetric) histograms
+ //
if (fAsymUncertainties) {
Double_t x = fhYieldCorr->GetBinCenter(ibin);
fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
if(correction>0.){
fgFcConservative->SetPoint(ibin,x,correction);
fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
+
+ fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
+ fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
+ // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
}
else{
fgFcConservative->SetPoint(ibin,x,0.);
// (in quadrature) with the feed-down subtraction uncertainties
//
-
// Estimate the feed-down uncertainty in percentage
Int_t nentries = fgSigmaCorrConservative->GetN();
TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
erryhcomb = systematics->GetTotalSystErr(x) * y ;
}
fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
+ //
+ fhSigmaCorrDataSyst->SetBinContent(i,y);
+ erryl = systematics->GetTotalSystErr(x) * y ;
+ fhSigmaCorrDataSyst->SetBinError(i,erryl);
}
}
return (TH1D*)hRecReweighted;
}
+
+
+//_________________________________________________________________________________________________________
+Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
+ //
+ // Function to find the y-axis bin of a TH2 for a given y-value
+ //
+
+ Int_t nbins = histo->GetNbinsY();
+ Int_t ybin=0;
+ for (int j=0; j<=nbins; j++) {
+ Float_t value = histo->GetYaxis()->GetBinCenter(j);
+ Float_t width = histo->GetYaxis()->GetBinWidth(j);
+ // if( TMath::Abs(yvalue-value)<= width/2. ) {
+ if( TMath::Abs(yvalue-value)<= width ) {
+ ybin =j;
+ // cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
+ break;
+ }
+ }
+
+ return ybin;
+}
//
// (the corrected yields per bin are divided by the bin-width)
//
+//
+// In HIC you can also evaluate how the feed-down correction is influenced by an energy loss hypothesis:
+// Raa(c-->D) / Raa(b-->D) defined here as Rcb for the "fc" method
+// Raa(b-->D) defined here as Rb for the "Nb" method
+//
// Author: Z.Conesa, zconesa@in2p3.fr
//***********************************************************************
#include "TMath.h"
class TH1;
+class TH2;
+class TNtuple;
class TGraphAsymmErrors;
-class AliHFSystErr;
class AliHFPtSpectrum: public TNamed
void SetFeedDownCalculationOption(Int_t option){ fFeedDownOption = option; }
// Set if the calculation has to consider asymmetric uncertaInt_ties or not
void SetComputeAsymmetricUncertainties(Bool_t flag){ fAsymUncertainties = flag; }
+ // Set if the calculation has to consider Ratio(c/b eloss) hypothesis
+ void SetComputeElossHypothesis(Bool_t flag){ fPbPbElossHypothesis = flag; }
// Set the luminosity and its uncertainty
void SetLuminosity(Double_t luminosity, Double_t unc){
fLuminosity[0]=luminosity; fLuminosity[1]=unc;
}
// Set the normalization factors
void SetNormalization(Double_t normalization){
- fLuminosity[0]=normalization; fTrigEfficiency[0]=1.0;
+ fLuminosity[0]=normalization;
}
- void SetNormalization(Double_t nevents, Double_t sigma){
- fLuminosity[0]=nevents/sigma; fTrigEfficiency[0]=1.0;
+ void SetNormalization(Int_t nevents, Double_t sigma){
+ fLuminosity[0]=nevents/sigma;
+ fNevts = nevents;
}
- void SetNormalization(Double_t nevents, Double_t sigma, Double_t sigmaunc){
+ void SetNormalization(Int_t nevents, Double_t sigma, Double_t sigmaunc){
fLuminosity[0] = nevents/sigma;
- fTrigEfficiency[0] = 1.0;
fLuminosity[1] = fLuminosity[0] * TMath::Sqrt( (1/nevents) + (sigmaunc/sigma)*(sigmaunc/sigma) );
+ fNevts = nevents;
+ }
+ //
+ // Set the Tab parameter and its uncertainty
+ void SetTabParameter(Double_t tabvalue, Double_t uncertainty){
+ fTab[0] = tabvalue;
+ fTab[1] = uncertainty;
}
+
//
// Getters
//
// Return the acceptance and efficiency corrections (rebinned if needed)
TH1D * GetDirectAccEffCorrection() const { return (fhDirectEffpt ? (TH1D*)fhDirectEffpt : NULL); }
TH1D * GetFeedDownAccEffCorrection() const { return (fhFeedDownEffpt ? (TH1D*)fhFeedDownEffpt : NULL); }
+ // Return whether the Ratio(c/b eloss) hypothesis has been considered
+ Bool_t IsElossHypothesisCalculated(){ return fPbPbElossHypothesis; }
// Return the TGraphAsymmErrors of the feed-down correction (extreme systematics)
TGraphAsymmErrors * GetFeedDownCorrectionFcExtreme() const { return (fgFcExtreme ? fgFcExtreme : NULL); }
// Return the TGraphAsymmErrors of the feed-down correction (conservative systematics)
// Return the histograms of the feed-down correction bounds
TH1D * GetHistoUpperLimitFeedDownCorrectionFc() const { return (fhFcMax ? (TH1D*)fhFcMax : NULL); }
TH1D * GetHistoLowerLimitFeedDownCorrectionFc() const { return (fhFcMin ? (TH1D*)fhFcMin : NULL); }
+ // Return the histogram of the feed-down correction times the Ratio(c/b eloss)
+ TH2D * GetHistoFeedDownCorrectionFcVsEloss() const { return (fhFcRcb ? (TH2D*)fhFcRcb : NULL); }
// Return the TGraphAsymmErrors of the yield after feed-down correction (systematics but feed-down)
TGraphAsymmErrors * GetFeedDownCorrectedSpectrum() const { return (fgYieldCorr ? fgYieldCorr : NULL); }
// Return the TGraphAsymmErrors of the yield after feed-down correction (feed-down extreme systematics)
// Return the histogram of the yield after feed-down correction bounds
TH1D * GetHistoUpperLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMax ? (TH1D*)fhYieldCorrMax : NULL); }
TH1D * GetHistoLowerLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMin ? (TH1D*)fhYieldCorrMin : NULL); }
+ // Return the histogram of the yield after feed-down correction vs the Ratio(c/b eloss)
+ TH2D * GetHistoFeedDownCorrectedSpectrumVsEloss() const { return (fhYieldCorrRcb ? (TH2D*)fhYieldCorrRcb : NULL); }
// Return the equivalent invariant cross-section TGraphAsymmErrors (systematics but feed-down)
TGraphAsymmErrors * GetCrossSectionFromYieldSpectrum() const { return (fgSigmaCorr ? fgSigmaCorr : NULL); }
// Return the equivalent invariant cross-section TGraphAsymmErrors (feed-down extreme systematics)
// Return the equivalent invariant cross-section histogram bounds
TH1D * GetHistoUpperLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMax ? (TH1D*)fhSigmaCorrMax : NULL); }
TH1D * GetHistoLowerLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMin ? (TH1D*)fhSigmaCorrMin : NULL); }
+ // Return the cross section systematics from data systematics
+ TH1D * GetHistoCrossSectionDataSystematics() const { return (fhSigmaCorrDataSyst ? (TH1D*)fhSigmaCorrDataSyst : NULL); }
+ //
+ // PbPb special calculations
+ // Return the equivalent invariant cross-section histogram vs the Ratio(c/b eloss)
+ TH2D * GetHistoCrossSectionFromYieldSpectrumVsEloss() const { return (fhSigmaCorrRcb ? (TH2D*)fhSigmaCorrRcb : NULL); }
+ // Return the ntuple of the calculation vs the Ratio(c/b eloss)
+ TNtuple * GetNtupleCrossSectionVsEloss() { return (fnSigma ? (TNtuple*)fnSigma : NULL); }
+ //
+ //
+ // Histograms to keep track of the influence of the efficiencies statistical uncertainty on the cross-section
+ TH1D * GetDirectStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffcSigma; }
+ TH1D * GetFeedDownStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffbSigma; }
+ // Histograms to keep track of the influence of the efficiencies statistical uncertainty on the feed-down correction factor
+ TH1D * GetDirectStatEffUncOnFc() const { return (TH1D*)fhStatUncEffcFD; }
+ TH1D * GetFeedDownStatEffUncOnFc() const { return (TH1D*)fhStatUncEffbFD; }
+
//
// Main function:
TH1D * ReweightHisto(TH1D *hToReweight, TH1D *hReference);
// to reweight the reco-histos: hRecToReweight is reweighted as hReference/hMCToReweight
TH1D * ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference);
+ // Functionality to find the y-axis bin of a TH2 for a given y-value
+ Int_t FindTH2YBin(TH2D *histo, Float_t yvalue);
protected:
TGraphAsymmErrors *fgRECSystematics; // all reconstructed D Systematic uncertainties
//
// Normalization factors
+ Int_t fNevts; // nb of analyzed events
Double_t fLuminosity[2]; // analyzed luminosity & uncertainty
Double_t fTrigEfficiency[2]; // trigger efficiency & uncertainty
Double_t fGlobalEfficiencyUncertainties[2]; // uncertainties on the efficiency [0]=c, b, [1]=b/c
+ Double_t fTab[2]; // Tab parameter and its uncertainty
//
// Output spectra
TH1D *fhFc; // Correction histo fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
TH1D *fhFcMax; // Maximum fc histo
TH1D *fhFcMin; // Minimum fc histo
+ TH2D *fhFcRcb; // Correction histo fc vs the Ratio(c/b eloss)
TGraphAsymmErrors * fgFcExtreme; // Extreme correction as TGraphAsymmErrors
TGraphAsymmErrors * fgFcConservative; // Extreme correction as TGraphAsymmErrors
TH1D *fhYieldCorr; // Corrected yield (stat unc. only)
TH1D *fhYieldCorrMax; // Maximum corrected yield
TH1D *fhYieldCorrMin; // Minimum corrected yield
+ TH2D *fhYieldCorrRcb; // Corrected yield (stat unc. only) vs the Ratio(c/b eloss)
TGraphAsymmErrors * fgYieldCorr; // Corrected yield as TGraphAsymmErrors (syst but feed-down)
TGraphAsymmErrors * fgYieldCorrExtreme; // Extreme corrected yield as TGraphAsymmErrors (syst from feed-down)
TGraphAsymmErrors * fgYieldCorrConservative; // Conservative corrected yield as TGraphAsymmErrors (syst from feed-down)
TH1D *fhSigmaCorr; // Corrected cross-section (stat unc. only)
TH1D *fhSigmaCorrMax; // Maximum corrected cross-section
TH1D *fhSigmaCorrMin; // Minimum corrected cross-section
+ TH1D *fhSigmaCorrDataSyst; // Corrected cross-section (syst. unc. from data only)
+ TH2D *fhSigmaCorrRcb; // Corrected cross-section (stat unc. only) vs the Ratio(c/b eloss)
TGraphAsymmErrors * fgSigmaCorr; // Corrected cross-section as TGraphAsymmErrors (syst but feed-down)
TGraphAsymmErrors * fgSigmaCorrExtreme; // Extreme corrected cross-section as TGraphAsymmErrors (syst from feed-down)
TGraphAsymmErrors * fgSigmaCorrConservative; // Conservative corrected cross-section as TGraphAsymmErrors (syst from feed-down)
+ //
+ TNtuple *fnSigma; // Ntuple of the calculation vs the Ratio(c/b eloss)
//
Int_t fFeedDownOption; // feed-down correction flag: 0=none, 1=fc, 2=Nb
Bool_t fAsymUncertainties; // flag: asymmetric uncertainties are (1) or not (0) considered
+ Bool_t fPbPbElossHypothesis; // flag: whether to do estimates vs Ratio(c/b eloss) hypothesis
+ //
+ TH1D *fhStatUncEffcSigma; // Uncertainty on the cross-section due to the prompt efficiency statistical uncertainty
+ TH1D *fhStatUncEffbSigma; // Uncertainty on the cross-section due to the feed-down efficiency statistical uncertainty
+ TH1D *fhStatUncEffcFD; // Uncertainty on the feed-down correction due to the prompt efficiency statistical uncertainty
+ TH1D *fhStatUncEffbFD; // Uncertainty on the feed-down correction due to the feed-down efficiency statistical uncertainty
- ClassDef(AliHFPtSpectrum,1) // Class for Heavy Flavor spectra corrections
+ ClassDef(AliHFPtSpectrum,2) // Class for Heavy Flavor spectra corrections
};
#endif
#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 "AliHFSystErr.h"
+
#include "AliHFPtSpectrum.h"
//
// 2) reconstructed spectra file name
// 3) output file name
// 4) Set the feed-down calculation option flag: 0=none, 1=fc only, 2=Nb only
-// 5) Set the luminosity
-// 6) Set the trigger efficiency
-// 7-14) If the efficiency histos do not have the right bin width, set the files & histo-names, they'll be computed, if the efficiencies are in file (6), don't set this parameters
+// 5-6) Set the luminosity: the number of events analyzed, and the cross-section of the sample
+// 7) Set the trigger efficiency
+// 8) Set the centrality class
+// 9) Flag to decide if there is need to evaluate the dependence on the energy loss
+// 10-17) If the efficiency histos do not have the right bin width, set the files & histo-names, they'll be computed, if the efficiencies are in file (6), don't set these parameters
//
+
+enum centrality{ kpp7, kpp276, k010, k020, k2040, k4060, k6080, k4080, k80100 };
+
void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
const char *efffilename="Efficiencies.root",
const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
const char *outfilename="HFPtSpectrum.root",
- Int_t option=1, Double_t lumi=1.0, Double_t effTrig=1.0,
+ // Int_t option=1, Double_t lumi=1.0, Double_t effTrig=1.0, Int_t cc=kpp, Bool_t PbPbEloss=false,
+ Int_t option=1, Double_t nevents=1.0, Double_t sigma=1.0, // sigma[nb]
+ Double_t effTrig=1.0, Int_t cc=kpp7, Bool_t PbPbEloss=false,
const char *directsimufilename="", const char *directsimuhistoname="CFHFccontainer0_New_3Prong_SelStep0_proj-pt",
const char *directrecofilename="", const char *directrecohistoname="CFHFccontainer0_New_3Prong_SelStep8_proj-pt",
const char *feeddownsimufilename="", const char *feeddownsimuhistoname="CFHFccontainer0allD_New_3Prong_SelStep0_proj-pt",
const char *feeddownrecofilename="", const char *feeddownrecohistoname="CFHFccontainer0allD_New_3Prong_SelStep8_proj-pt") {
+
+ gROOT->Macro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
+
// Set if calculation considers asymmetric uncertainties or not
- bool asym = true;
+ Bool_t asym = true;
// Set the meson and decay
// (only D0 -> K pi, D+--> K pi pi & D* --> D0 pi implemented here)
- bool isD0Kpi = true;
- bool isDplusKpipi = false;
- bool isDstarD0pi = false;
+ Bool_t isD0Kpi = true;
+ Bool_t isDplusKpipi = false;
+ Bool_t isDstarD0pi = false;
if (isD0Kpi && isDplusKpipi && isDstarD0pi) {
cout << "Sorry, can not deal with more than one correction at the same time"<<endl;
return;
}
if (option==0) asym = false;
+
+ //
+ // Defining the Tab values for the given centrality class
+ //
+ Double_t tab = 1., tabUnc = 0.;
+ if ( cc == k010 ) {
+ tab = 23.48; tabUnc = 0.97;
+ } else if ( cc == k020 ) {
+ tab = 18.93; tabUnc = 0.74;
+ } else if ( cc == k2040 ) {
+ tab = 6.86; tabUnc = 0.28;
+ } else if ( cc == k4060 ) {
+ tab = 2.00; tabUnc= 0.11;
+ } else if ( cc == k4080 ) {
+ tab = 1.20451; tabUnc = 0.071843;
+ } else if ( cc == k6080 ) {
+ tab = 0.419; tabUnc = 0.033;
+ } else if ( cc == k80100 ){
+ tab = 0.0690; tabUnc = 0.0062;
+ }
+ tab *= 1e-9; // to pass from mb^{-1} to pb^{-1}
+ tabUnc *= 1e-9;
+
+
+
//
// Get the histograms from the files
//
TH1D *hDirectMCptMin=0; // Input MC minimum c-->D spectra
TH1D *hFeedDownMCptMax=0; // Input MC maximum b-->D spectra
TH1D *hFeedDownMCptMin=0; // Input MC minimum b-->D spectra
- TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
+ // TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
TH1D *hDirectEffpt=0; // c-->D Acceptance and efficiency correction
TH1D *hFeedDownEffpt=0; // b-->D Acceptance and efficiency correction
TH1D *hRECpt=0; // all reconstructed D
hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min");
hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max");
hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min");
- gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
+ // gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
}
else if (isDplusKpipi){
decay = 2;
hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min");
hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max");
hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min");
- gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
+ // gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
}
else if(isDstarD0pi){
decay = 3;
hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min");
hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max");
hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min");
- gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
+ // gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
}
//
hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
}
else {
TFile * efffile = new TFile(efffilename,"read");
- hDirectEffpt = (TH1D*)efffile->Get("hEffD");
+ hDirectEffpt = (TH1D*)efffile->Get("hEffD_rebin");//hDirectEffpt");//hRecoPIDGenLimAcc");//hDirectEffpt");//hEffD");
hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
- hFeedDownEffpt = (TH1D*)efffile->Get("hEffB");
+ hFeedDownEffpt = (TH1D*)efffile->Get("hEffB_rebin");//hFeedDownEffpt");//hRecoPIDGenLimAccFromB");//hFeedDownEffpt");//hEffB");
hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
}
//
TH1D *histoSigmaCorrMax=0;
TH1D *histoSigmaCorrMin=0;
//
- int nbins = hRECpt->GetNbinsX();
+ TH2D *histofcRcb=0;
+ TH1D *histofcRcb_px=0;
+ TH2D *histoYieldCorrRcb=0;
+ TH2D *histoSigmaCorrRcb=0;
+ //
+ Int_t nbins = hRECpt->GetNbinsX();
TGraphAsymmErrors * gYieldCorr = 0;
TGraphAsymmErrors * gSigmaCorr = 0;
TGraphAsymmErrors * gFcExtreme = 0;
TGraphAsymmErrors * gSigmaCorrExtreme = 0;
TGraphAsymmErrors * gYieldCorrConservative = 0;
TGraphAsymmErrors * gSigmaCorrConservative = 0;
+ //
+ TNtuple * nSigma = 0;
//
AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
spectra->SetFeedDownCalculationOption(option);
spectra->SetComputeAsymmetricUncertainties(asym);
+ // Set flag on whether to additional PbPb Eloss hypothesis have to be computed
+ spectra->SetComputeElossHypothesis(PbPbEloss);
// Feed the input histograms
// reconstructed spectra
cout << " and the normalization" <<endl;
// Set normalization factors (uncertainties set to 0. as example)
- double lumiUnc = 0.10*lumi; // 10% uncertainty on the luminosity
+ spectra->SetNormalization(nevents,sigma);
+ Double_t lumi = nevents / sigma ;
+ Double_t lumiUnc = 0.07*lumi; // 10% uncertainty on the luminosity
spectra->SetLuminosity(lumi,lumiUnc);
spectra->SetTriggerEfficiency(effTrig,0.);
// Set the global uncertainties on the efficiencies (in percent)
- double globalEffUnc = 0.15;
- double globalBCEffRatioUnc = 0.15;
- // double globalEffUnc = 0.;
- // double globalBCEffRatioUnc = 0.;
+ Double_t globalEffUnc = 0.15;
+ Double_t globalBCEffRatioUnc = 0.15;
spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc);
+ // Set the Tab parameter and uncertainties
+ if ( (cc != kpp7) && (cc != kpp276) ) {
+ spectra->SetTabParameter(tab,tabUnc);
+ }
+
// Do the calculations
cout << " Doing the calculation... "<< endl;
- double deltaY = 1.0;
- double branchingRatioC = 1.0;
- double branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
+ Double_t deltaY = 1.0;
+ Double_t branchingRatioC = 1.0;
+ Double_t branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
cout << " ended the calculation, getting the histograms back " << endl;
// Set the systematics externally
+ Bool_t combineFeedDown = true;
AliHFSystErr *systematics = new AliHFSystErr();
+ if( cc==kpp276 ) {
+ systematics->SetIsLowEnergy(true);
+ } else if( cc!=kpp7 ) {
+ systematics->SetCollisionType(1);
+ if ( cc == k020 ) {
+ systematics->SetCentrality("020");
+ }
+ else if ( cc == k4080 ) {
+ systematics->SetCentrality("4080");
+ }
+ else {
+ cout << " Systematics not yet implemented " << endl;
+ return;
+ }
+ } else { systematics->SetCollisionType(0); }
systematics->Init(decay);
- bool combineFeedDown = true;
spectra->ComputeSystUncertainties(systematics,combineFeedDown);
-
//
// Get the output histograms
//
// the efficiencies
if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
+ // Get the PbPb Eloss hypothesis histograms
+ if(PbPbEloss){
+ histofcRcb = spectra->GetHistoFeedDownCorrectionFcVsEloss();
+ histoYieldCorrRcb = spectra->GetHistoFeedDownCorrectedSpectrumVsEloss();
+ histoSigmaCorrRcb = spectra->GetHistoCrossSectionFromYieldSpectrumVsEloss();
+ histofcRcb->SetName("histofcRcb");
+ histoYieldCorrRcb->SetName("histoYieldCorrRcb");
+ histoSigmaCorrRcb->SetName("histoSigmaCorrRcb");
+ }
// Get & Rename the TGraphs
if (asym) {
gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
}
+ if(PbPbEloss){
+ nSigma = spectra->GetNtupleCrossSectionVsEloss();
+ }
+
//
// Now, plot the results ! :)
//
+ gROOT->SetStyle("Plain");
+
cout << " Drawing the results ! " << endl;
// control plots
}
+ // Draw the PbPb Eloss hypothesis histograms
+ if(PbPbEloss){
+ AliHFPtSpectrum *CalcBins;
+ gStyle->SetPalette(1);
+ TCanvas *canvasfcRcb = new TCanvas("canvasfcRcb","fc vs pt vs Rcb");
+ // histofcRcb->Draw("cont4z");
+ histofcRcb->Draw("colz");
+ canvasfcRcb->Update();
+ canvasfcRcb->cd(2);
+ TCanvas *canvasfcRcb1 = new TCanvas("canvasfcRcb1","fc vs pt vs Rcb=1");
+ histofcRcb_px = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px",40,40);
+ histofcRcb_px->SetLineColor(2);
+ if (option==1) {
+ histofc->Draw();
+ histofcRcb_px->Draw("same");
+ } else histofcRcb_px->Draw("");
+ canvasfcRcb1->Update();
+ TCanvas *canvasfcRcb2 = new TCanvas("canvasfcRcb2","fc vs pt vs Rcb fixed Rcb");
+ Int_t bin0 = CalcBins->FindTH2YBin(histofcRcb,0.25);
+ Int_t bin1 = CalcBins->FindTH2YBin(histofcRcb,0.5);
+ Int_t bin2 = CalcBins->FindTH2YBin(histofcRcb,1.0);
+ Int_t bin3 = CalcBins->FindTH2YBin(histofcRcb,1.5);
+ Int_t bin4 = CalcBins->FindTH2YBin(histofcRcb,2.0);
+ Int_t bin5 = CalcBins->FindTH2YBin(histofcRcb,3.0);
+ Int_t bin6 = CalcBins->FindTH2YBin(histofcRcb,4.0);
+ TH1D * histofcRcb_px0a = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0a",bin0,bin0);
+ TH1D * histofcRcb_px0 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0",bin1,bin1);
+ TH1D * histofcRcb_px1 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px1",bin2,bin2);
+ TH1D * histofcRcb_px2 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px2",bin3,bin3);
+ TH1D * histofcRcb_px3 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px3",bin4,bin4);
+ TH1D * histofcRcb_px4 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px4",bin5,bin5);
+ TH1D * histofcRcb_px5 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px5",bin6,bin6);
+ if (option==1) {
+ histofc->Draw();
+ // histofcRcb_px->Draw("same");
+ } else {
+ // histofcRcb_px->Draw("");
+ histofcRcb_px0a->SetLineColor(2);
+ histofcRcb_px0a->Draw("");
+ }
+ histofcRcb_px0a->SetLineColor(2);
+ histofcRcb_px0a->Draw("same");
+ histofcRcb_px0->SetLineColor(4);
+ histofcRcb_px0->Draw("same");
+ histofcRcb_px1->SetLineColor(3);
+ histofcRcb_px1->Draw("same");
+ histofcRcb_px2->SetLineColor(kCyan);
+ histofcRcb_px2->Draw("same");
+ histofcRcb_px3->SetLineColor(kMagenta+1);
+ histofcRcb_px3->Draw("same");
+ histofcRcb_px4->SetLineColor(kOrange+7);
+ histofcRcb_px4->Draw("same");
+ histofcRcb_px5->SetLineColor(kGreen+3);
+ histofcRcb_px5->Draw("same");
+ TLegend *legrcc = new TLegend(0.8,0.8,0.95,0.9);
+ legrcc->SetFillColor(0);
+ if (option==1) {
+ legrcc->AddEntry(histofcRcb_px0a,"Rc/b=0.25","l");
+ legrcc->AddEntry(histofcRcb_px0,"Rc/b=0.5","l");
+ legrcc->AddEntry(histofcRcb_px1,"Rc/b=1.0","l");
+ legrcc->AddEntry(histofcRcb_px2,"Rc/b=1.5","l");
+ legrcc->AddEntry(histofcRcb_px3,"Rc/b=2.0","l");
+ legrcc->AddEntry(histofcRcb_px4,"Rc/b=3.0","l");
+ legrcc->AddEntry(histofcRcb_px5,"Rc/b=4.0","l");
+ }else{
+ legrcc->AddEntry(histofcRcb_px0a,"Rb=0.25","l");
+ legrcc->AddEntry(histofcRcb_px0,"Rb=0.5","l");
+ legrcc->AddEntry(histofcRcb_px1,"Rb=1.0","l");
+ legrcc->AddEntry(histofcRcb_px2,"Rb=1.5","l");
+ legrcc->AddEntry(histofcRcb_px3,"Rb=2.0","l");
+ legrcc->AddEntry(histofcRcb_px4,"Rb=3.0","l");
+ legrcc->AddEntry(histofcRcb_px5,"Rb=4.0","l");
+ }
+ legrcc->Draw();
+ canvasfcRcb2->Update();
+ TCanvas *canvasYRcb = new TCanvas("canvasYRcb","corrected yield vs pt vs Rcb");
+ histoYieldCorrRcb->Draw("cont4z");
+ canvasYRcb->Update();
+ TCanvas *canvasSRcb = new TCanvas("canvasSRcb","sigma vs pt vs Rcb");
+ histoSigmaCorrRcb->Draw("cont4z");
+ canvasSRcb->Update();
+ TCanvas *canvasSRcb1 = new TCanvas("canvasSRcb1","sigma vs pt vs Rcb fixed Rcb");
+ TH1D * histoSigmaCorrRcb_px0a = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0a",bin0,bin0);
+ TH1D * histoSigmaCorrRcb_px0 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0",bin1,bin1);
+ TH1D * histoSigmaCorrRcb_px1 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px1",bin2,bin2);
+ TH1D * histoSigmaCorrRcb_px2 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px2",bin3,bin3);
+ TH1D * histoSigmaCorrRcb_px3 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px3",bin4,bin4);
+ TH1D * histoSigmaCorrRcb_px4 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px4",bin5,bin5);
+ TH1D * histoSigmaCorrRcb_px5 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px5",bin6,bin6);
+ histoSigmaCorr->Draw();
+ histoSigmaCorrRcb_px0a->SetLineColor(2);
+ histoSigmaCorrRcb_px0a->Draw("hsame");
+ histoSigmaCorrRcb_px0->SetLineColor(4);
+ histoSigmaCorrRcb_px0->Draw("hsame");
+ histoSigmaCorrRcb_px1->SetLineColor(3);
+ histoSigmaCorrRcb_px1->Draw("hsame");
+ histoSigmaCorrRcb_px2->SetLineColor(kCyan);
+ histoSigmaCorrRcb_px2->Draw("hsame");
+ histoSigmaCorrRcb_px3->SetLineColor(kMagenta+1);
+ histoSigmaCorrRcb_px3->Draw("hsame");
+ histoSigmaCorrRcb_px4->SetLineColor(kOrange+7);
+ histoSigmaCorrRcb_px4->Draw("same");
+ histoSigmaCorrRcb_px5->SetLineColor(kGreen+3);
+ histoSigmaCorrRcb_px5->Draw("same");
+ TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
+ legrcb->SetFillColor(0);
+ if (option==1) {
+ legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rc/b=0.25","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rc/b=0.5","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rc/b=1.0","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rc/b=1.5","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rc/b=2.0","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rc/b=3.0","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rc/b=4.0","l");
+ }else{
+ legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rb=0.25","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rb=0.5","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rb=1.0","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rb=1.5","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rb=2.0","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rb=3.0","l");
+ legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rb=4.0","l");
+ }
+ legrcb->Draw();
+ canvasSRcb1->Update();
+ }
//
histoSigmaCorr->Write();
histoSigmaCorrMax->Write(); histoSigmaCorrMin->Write();
+ if(PbPbEloss){
+ histofcRcb->Write(); histofcRcb_px->Write();
+ histoYieldCorrRcb->Write();
+ histoSigmaCorrRcb->Write();
+ nSigma->Write();
+ }
+
if(asym){
gYieldCorr->Write();
gSigmaCorr->Write();
}
+ TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma();
+ TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma();
+ TH1D * hStatUncEffcFD = spectra->GetDirectStatEffUncOnFc();
+ TH1D * hStatUncEffbFD = spectra->GetFeedDownStatEffUncOnFc();
+ hStatUncEffcSigma->Write();
+ hStatUncEffbSigma->Write();
+ hStatUncEffcFD->Write();
+ hStatUncEffbFD->Write();
+
// Draw the cross-section
// spectra->DrawSpectrum(gPrediction);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff