#include <Riostream.h>
#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"
#include "TLegend.h"
#include "AliHFSystErr.h"
#include "AliHFPtSpectrum.h"
ClassImp(AliHFPtSpectrum)
AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
TNamed(name,title),
fhDirectMCpt(NULL),
fhFeedDownMCpt(NULL),
fhDirectMCptMax(NULL),
fhDirectMCptMin(NULL),
fhFeedDownMCptMax(NULL),
fhFeedDownMCptMin(NULL),
fhDirectEffpt(NULL),
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),
fnHypothesis(NULL),
fFeedDownOption(option),
fAsymUncertainties(kTRUE),
fPbPbElossHypothesis(kFALSE),
fIsStatUncEff(kTRUE),
fParticleAntiParticle(2),
fIsEventPlane(kFALSE),
fnPtBins(0),
fPtBinLimits(NULL),
fPtBinWidths(NULL),
fhStatUncEffcSigma(NULL),
fhStatUncEffbSigma(NULL),
fhStatUncEffcFD(NULL),
fhStatUncEffbFD(NULL)
{
fLuminosity[0]=1.; fLuminosity[1]=0.;
fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
fTab[0]=1.; fTab[1]=0.;
}
AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
TNamed(rhs),
fhDirectMCpt(rhs.fhDirectMCpt),
fhFeedDownMCpt(rhs.fhFeedDownMCpt),
fhDirectMCptMax(rhs.fhDirectMCptMax),
fhDirectMCptMin(rhs.fhDirectMCptMin),
fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
fhDirectEffpt(rhs.fhDirectEffpt),
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),
fnHypothesis(rhs.fnHypothesis),
fFeedDownOption(rhs.fFeedDownOption),
fAsymUncertainties(rhs.fAsymUncertainties),
fPbPbElossHypothesis(rhs.fPbPbElossHypothesis),
fIsStatUncEff(rhs.fIsStatUncEff),
fParticleAntiParticle(rhs.fParticleAntiParticle),
fIsEventPlane(rhs.fIsEventPlane),
fnPtBins(rhs.fnPtBins),
fPtBinLimits(),
fPtBinWidths(),
fhStatUncEffcSigma(NULL),
fhStatUncEffbSigma(NULL),
fhStatUncEffcFD(NULL),
fhStatUncEffbFD(NULL)
{
for(Int_t i=0; i<2; i++){
fLuminosity[i] = rhs.fLuminosity[i];
fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
fTab[i] = rhs.fTab[i];
}
for(Int_t i=0; i<fnPtBins; i++){
fPtBinLimits[i] = rhs.fPtBinLimits[i];
fPtBinWidths[i] = rhs.fPtBinWidths[i];
}
fPtBinLimits[fnPtBins] = rhs.fPtBinLimits[fnPtBins];
}
AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
if (&source == this) return *this;
fhDirectMCpt = source.fhDirectMCpt;
fhFeedDownMCpt = source.fhFeedDownMCpt;
fhDirectMCptMax = source.fhDirectMCptMax;
fhDirectMCptMin = source.fhDirectMCptMin;
fhFeedDownMCptMax = source.fhFeedDownMCptMax;
fhFeedDownMCptMin = source.fhFeedDownMCptMin;
fhDirectEffpt = source.fhDirectEffpt;
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;
fnHypothesis = source.fnHypothesis;
fFeedDownOption = source.fFeedDownOption;
fAsymUncertainties = source.fAsymUncertainties;
fPbPbElossHypothesis = source.fPbPbElossHypothesis;
fIsStatUncEff = source.fIsStatUncEff;
fParticleAntiParticle = source.fParticleAntiParticle;
fIsEventPlane = source.fIsEventPlane;
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];
}
fnPtBins = source.fnPtBins;
for(Int_t i=0; i<fnPtBins; i++){
fPtBinLimits[i] = source.fPtBinLimits[i];
fPtBinWidths[i] = source.fPtBinWidths[i];
}
fPtBinLimits[fnPtBins] = source.fPtBinLimits[fnPtBins];
return *this;
}
AliHFPtSpectrum::~AliHFPtSpectrum(){
if (fhDirectMCpt) delete fhDirectMCpt;
if (fhFeedDownMCpt) delete fhFeedDownMCpt;
if (fhDirectMCptMax) delete fhDirectMCptMax;
if (fhDirectMCptMin) delete fhDirectMCptMin;
if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
if (fhDirectEffpt) delete fhDirectEffpt;
if (fhFeedDownEffpt) delete fhFeedDownEffpt;
if (fhRECpt) delete fhRECpt;
if (fgRECSystematics) delete fgRECSystematics;
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 (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 (fhSigmaCorrDataSyst) delete fhSigmaCorrDataSyst;
if (fgSigmaCorr) delete fgSigmaCorr;
if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
if (fnSigma) delete fnSigma;
if (fnHypothesis) delete fnHypothesis;
if (fPtBinLimits) delete [] fPtBinLimits;
if (fPtBinWidths) delete [] fPtBinWidths;
}
TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
if (!hTheory || !fhRECpt) {
AliError("Feed-down or reconstructed spectra don't exist");
return NULL;
}
Int_t nbinsMC = hTheory->GetNbinsX();
Double_t thbinwidth = hTheory->GetBinWidth(1);
for (Int_t i=1; i<=fnPtBins; i++) {
if ( thbinwidth > fPtBinWidths[i-1] ) {
AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
}
}
TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",fnPtBins,fPtBinLimits);
Double_t sum[fnPtBins], items[fnPtBins];
for (Int_t ibin=0; ibin<fnPtBins; ibin++) {
sum[ibin]=0.; items[ibin]=0.;
}
for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++){
if (hTheory->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
hTheory->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1]){
sum[ibinrec]+=hTheory->GetBinContent(ibin);
items[ibinrec]+=1.;
}
}
}
for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++) {
hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
}
return (TH1D*)hTheoryRebin;
}
void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
if (!hDirect || !hFeedDown || !fhRECpt) {
AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
return;
}
Bool_t areconsistent = kTRUE;
areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
if (!areconsistent) {
AliInfo("Histograms are not consistent (bin width, bounds)");
return;
}
fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
}
void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
if (!hFeedDown || !fhRECpt) {
AliError("Feed-down or reconstructed spectra don't exist");
return;
}
fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
}
void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
return;
}
Bool_t areconsistent = kTRUE;
areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
if (!areconsistent) {
AliInfo("Histograms are not consistent (bin width, bounds)");
return;
}
fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
}
void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
AliError("One or all of the max/min direct/feed-down spectra don't exist");
return;
}
Bool_t areconsistent = kTRUE;
areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
if (!areconsistent) {
AliInfo("Histograms are not consistent (bin width, bounds)");
return;
}
fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
}
void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
if (!hDirectEff) {
AliError("The direct acceptance and efficiency corrections doesn't exist");
return;
}
fhDirectEffpt = (TH1D*)hDirectEff->Clone();
fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
}
void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
if (!hDirectEff || !hFeedDownEff) {
AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
return;
}
Bool_t areconsistent=kTRUE;
areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
if (!areconsistent) {
AliInfo("Histograms are not consistent (bin width, bounds)");
return;
}
fhDirectEffpt = (TH1D*)hDirectEff->Clone();
fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
}
void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
if (!hRec) {
AliError("The reconstructed spectrum doesn't exist");
return;
}
fhRECpt = (TH1D*)hRec->Clone();
fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
fnPtBins = fhRECpt->GetNbinsX();
fPtBinLimits = new Double_t[fnPtBins+1];
fPtBinWidths = new Double_t[fnPtBins];
Double_t xlow=0., binwidth=0.;
for (Int_t i=1; i<=fnPtBins; i++) {
binwidth = fhRECpt->GetBinWidth(i);
xlow = fhRECpt->GetBinLowEdge(i);
fPtBinLimits[i-1] = xlow;
fPtBinWidths[i-1] = binwidth;
}
fPtBinLimits[fnPtBins] = xlow + binwidth;
}
void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
Double_t hbinwidth = fhRECpt->GetBinWidth(1);
Double_t gxbincenter=0., gybincenter=0.;
gRec->GetPoint(1,gxbincenter,gybincenter);
Double_t hbincenter = fhRECpt->GetBinCenter(1);
if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
AliError(" The reconstructed spectrum and its systematics don't seem compatible");
return;
}
fgRECSystematics = gRec;
}
void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
Bool_t areHistosOk = Initialize();
if (!areHistosOk) {
AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
return;
}
if(!fIsStatUncEff) ResetStatUncEff();
if (fFeedDownOption==1) {
CalculateFeedDownCorrectionFc();
CalculateFeedDownCorrectedSpectrumFc();
}
else if (fFeedDownOption==2) {
CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
}
else if (fFeedDownOption==0) {
CalculateCorrectedSpectrumNoFeedDown();
}
else {
AliInfo(" Are you sure the feed-down correction option is right ?");
}
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]);
fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",fnPtBins,fPtBinLimits);
fgSigmaCorr = new TGraphAsymmErrors(fnPtBins+1);
fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",fnPtBins,fPtBinLimits);
fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",fnPtBins,fPtBinLimits);
fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",fnPtBins,fPtBinLimits);
if (fPbPbElossHypothesis && fFeedDownOption==1) {
fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",fnPtBins,fPtBinLimits,800,0.,4.);
fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
}
if (fPbPbElossHypothesis && fFeedDownOption==2) {
fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",fnPtBins,fPtBinLimits,800,0.,4.);
fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
}
if (fAsymUncertainties){
fgSigmaCorrExtreme = new TGraphAsymmErrors(fnPtBins+1);
fgSigmaCorrConservative = new TGraphAsymmErrors(fnPtBins+1);
}
fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",fnPtBins,fPtBinLimits);
fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",fnPtBins,fPtBinLimits);
if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
return ;
}
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.;
for(Int_t ibin=1; ibin<=fnPtBins; ibin++){
value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
errvalueExtremeMax=0.; errvalueExtremeMin=0.;
errvalueConservativeMax=0.; errvalueConservativeMin=0.;
errvalueStatUncEffc=0.;
value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
: 0. ;
errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
if (fAsymUncertainties && value>0.) {
errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
(fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
(fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
(fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
(fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
}
else {
errvalueMax = (value!=0.) ?
value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
(fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
: 0. ;
errvalueMin = errvalueMax;
}
fhSigmaCorr->SetBinContent(ibin,value);
fhSigmaCorr->SetBinError(ibin,errValue);
Double_t x = fhYieldCorr->GetBinCenter(ibin);
fgSigmaCorr->SetPoint(ibin,x,value);
fgSigmaCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax);
fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
if (fPbPbElossHypothesis) {
if(!fnHypothesis) {
AliError("Error reading the fc hypothesis no ntuple, please check !!");
return ;
}
Int_t entriesHypo = fnHypothesis->GetEntries();
Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
fnHypothesis->SetBranchAddress("pt",&pt);
if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
fnHypothesis->SetBranchAddress("fc",&fc);
fnHypothesis->SetBranchAddress("fcMax",&fcMax);
fnHypothesis->SetBranchAddress("fcMin",&fcMin);
for (Int_t item=0; item<entriesHypo; item++) {
fnHypothesis->GetEntry(item);
if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
: 0. ;
Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
: 0. ;
Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
: 0. ;
Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) ) : 0. ;
fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
}
}
if (fAsymUncertainties) {
fgSigmaCorrExtreme->SetPoint(ibin,x,value);
fgSigmaCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax);
fgSigmaCorrConservative->SetPoint(ibin,x,value);
fgSigmaCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax);
fhStatUncEffcSigma->SetBinContent(ibin,0.);
if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
}
}
}
TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
if(!fhRECpt){
AliInfo("Hey, the reconstructed histogram was not set yet !");
return NULL;
}
TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",fnPtBins,fPtBinLimits);
Double_t *sumSimu=new Double_t[fnPtBins];
Double_t *sumReco=new Double_t[fnPtBins];
for (Int_t ibin=0; ibin<fnPtBins; ibin++){
sumSimu[ibin]=0.; sumReco[ibin]=0.;
}
for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++){
if ( hSimu->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
hSimu->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1] ) {
sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
}
if ( hReco->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
hReco->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1] ) {
sumReco[ibinrec]+=hReco->GetBinContent(ibin);
}
}
}
Double_t eff=0., erreff=0.;
for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++) {
if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
}
else { eff=0.0; erreff=0.; }
hEfficiency->SetBinContent(ibinrec+1,eff);
hEfficiency->SetBinError(ibinrec+1,erreff);
}
delete [] sumSimu;
delete [] sumReco;
return (TH1D*)hEfficiency;
}
void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
if(!fhRECpt || !hSimu || !hReco){
AliError("Hey, the reconstructed histogram was not set yet !");
return;
}
fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
}
void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
if(!fhRECpt || !hSimu || !hReco){
AliError("Hey, the reconstructed histogram was not set yet !");
return;
}
fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
}
Bool_t AliHFPtSpectrum::Initialize(){
if (fFeedDownOption==0) {
AliInfo("Getting ready for the corrections without feed-down consideration");
} else if (fFeedDownOption==1) {
AliInfo("Getting ready for the fc feed-down correction calculation");
} else if (fFeedDownOption==2) {
AliInfo("Getting ready for the Nb feed-down correction calculation");
} else { AliError("The calculation option must be <=2"); return kFALSE; }
Bool_t areconsistent=kTRUE;
if (!fhDirectEffpt || !fhRECpt) {
AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
return kFALSE;
}
areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
if (!areconsistent) {
AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
return kFALSE;
}
if (fFeedDownOption==0) return kTRUE;
if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
return kFALSE;
}
areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
if (fAsymUncertainties) {
if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
AliError(" Max/Min theoretical Nb distributions are not defined");
return kFALSE;
}
areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
}
if (!areconsistent) {
AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
return kFALSE;
}
if (fFeedDownOption>1) return kTRUE;
if (!fhDirectMCpt) {
AliError("Theoretical Nc distributions is not defined");
return kFALSE;
}
areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
if (fAsymUncertainties) {
if (!fhDirectMCptMax || !fhDirectMCptMin) {
AliError(" Max/Min theoretical Nc distributions are not defined");
return kFALSE;
}
areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
}
if (!areconsistent) {
AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
return kFALSE;
}
return kTRUE;
}
Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
if (!h1 || !h2) {
AliError("One or both histograms don't exist");
return kFALSE;
}
Double_t binwidth1 = h1->GetBinWidth(1);
Double_t binwidth2 = h2->GetBinWidth(1);
Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
if (binwidth1!=binwidth2) {
AliInfo(" histograms with different bin width");
return kFALSE;
}
if (min1!=min2) {
AliInfo(" histograms with different minimum");
return kFALSE;
}
return kTRUE;
}
void AliHFPtSpectrum::CalculateCorrectedSpectrumNoFeedDown(){
fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
fgYieldCorr = new TGraphAsymmErrors(fnPtBins+1);
for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
Double_t value = fhRECpt->GetBinContent(ibin) /fhRECpt->GetBinWidth(ibin);
Double_t errvalue = fhRECpt->GetBinError(ibin) /fhRECpt->GetBinWidth(ibin);
fhYieldCorr->SetBinContent(ibin,value);
fhYieldCorr->SetBinError(ibin,errvalue);
fhYieldCorrMax->SetBinContent(ibin,value);
fhYieldCorrMax->SetBinError(ibin,errvalue);
fhYieldCorrMin->SetBinContent(ibin,value);
fhYieldCorrMin->SetBinError(ibin,errvalue);
}
fAsymUncertainties=kFALSE;
}
void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
AliInfo("Calculating the feed-down correction factor (fc method)");
Double_t correction=1.;
Double_t theoryRatio=1.;
Double_t effRatio=1.;
Double_t correctionExtremeA=1., correctionExtremeB=1.;
Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
Double_t correctionConservativeA=1., correctionConservativeB=1.;
Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
fhFc = new TH1D("fhFc","fc correction factor",fnPtBins,fPtBinLimits);
fhFcMax = new TH1D("fhFcMax","max fc correction factor",fnPtBins,fPtBinLimits);
fhFcMin = new TH1D("fhFcMin","min fc correction factor",fnPtBins,fPtBinLimits);
if(fPbPbElossHypothesis) {
fhFcRcb = new TH2D("fhFcRcb","fc correction factor vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; fc correction",fnPtBins,fPtBinLimits,800,0.,4.);
fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
}
TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",fnPtBins,fPtBinLimits);
TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",fnPtBins,fPtBinLimits);
if (fAsymUncertainties) {
fgFcExtreme = new TGraphAsymmErrors(fnPtBins+1);
fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
fgFcConservative = new TGraphAsymmErrors(fnPtBins+1);
fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
}
fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
correction=1.; theoryRatio=1.; effRatio=1.;
correctionExtremeA=1.; correctionExtremeB=1.;
theoryRatioExtremeA=1.; theoryRatioExtremeB=1.;
correctionConservativeA=1.; correctionConservativeB=1.;
theoryRatioConservativeA=1.; theoryRatioConservativeB=1.;
correctionExtremeAUnc=0.; correctionExtremeBUnc=0.;
correctionConservativeAUnc=0.; correctionConservativeBUnc=0.;
correctionConservativeAUncStatEffc=0.; correctionConservativeBUncStatEffc=0.;
correctionConservativeAUncStatEffb=0.; correctionConservativeBUncStatEffb=0.;
theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
if( TMath::Abs(effRatio - 1.0)<0.0001 || TMath::Abs(theoryRatio - 1.0)<0.0001 ) {
correction = 1.0;
correctionExtremeA = 1.0;
correctionExtremeB = 1.0;
correctionConservativeA = 1.0;
correctionConservativeB = 1.0;
}
else {
correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
}
Double_t relEffUnc = 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))
);
correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio * relEffUnc;
correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio * relEffUnc;
correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio * relEffUnc;
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 * relEffUnc;
correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
hTheoryRatio->SetBinContent(ibin,theoryRatio);
hEffRatio->SetBinContent(ibin,effRatio);
fhFc->SetBinContent(ibin,correction);
if ( correction>1.0e-16 && fPbPbElossHypothesis){
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 );
Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
Double_t uncConservativeRcbMax = TMath::MaxElement(4,consvalRcb) - correctionRcb;
fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
}
}
if (fAsymUncertainties) {
Double_t x = fhDirectMCpt->GetBinCenter(ibin);
Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
fgFcExtreme->SetPoint(ibin,x,correction);
fgFcExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax);
fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
fgFcConservative->SetPoint(ibin,x,correction);
fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax);
if( !(correction>0.) ){
fgFcExtreme->SetPoint(ibin,x,0.);
fgFcExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.);
fgFcConservative->SetPoint(ibin,x,0.);
fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.);
}
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.);
}
}
}
void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
if (!fhFc || !fhRECpt) {
AliError(" Reconstructed or fc distributions are not defined");
return;
}
Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",fnPtBins,fPtBinLimits);
fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",fnPtBins,fPtBinLimits);
fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",fnPtBins,fPtBinLimits);
if(fPbPbElossHypothesis) fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by fc) vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; corrected yield",fnPtBins,fPtBinLimits,800,0.,4.);
fgYieldCorr = new TGraphAsymmErrors(fnPtBins+1);
if (fAsymUncertainties){
fgYieldCorrExtreme = new TGraphAsymmErrors(fnPtBins+1);
fgYieldCorrConservative = new TGraphAsymmErrors(fnPtBins+1);
}
for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
value = 0.; errvalue = 0.; errvalueMax= 0.; errvalueMin= 0.;
valueExtremeMax= 0.; valueExtremeMin= 0.;
valueConservativeMax= 0.; valueConservativeMin= 0.;
value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
value /= fhRECpt->GetBinWidth(ibin) ;
errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
if (fAsymUncertainties) {
if (fgRECSystematics) {
errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
}
else { errvalueMax = 0.; errvalueMin = 0.; }
valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
}
}
else { errvalueMax = 0.; errvalueMin = 0.; }
fhYieldCorr->SetBinContent(ibin,value);
fhYieldCorr->SetBinError(ibin,errvalue);
if (fPbPbElossHypothesis) {
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);
Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
}
}
if (fAsymUncertainties) {
Double_t center = fhYieldCorr->GetBinCenter(ibin);
fgYieldCorr->SetPoint(ibin,center,value);
fgYieldCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax);
fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
fgYieldCorrExtreme->SetPoint(ibin,center,value);
fgYieldCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
fgYieldCorrConservative->SetPoint(ibin,center,value);
fgYieldCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
}
}
}
void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",fnPtBins,fPtBinLimits);
fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",fnPtBins,fPtBinLimits);
fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",fnPtBins,fPtBinLimits);
if(fPbPbElossHypothesis) {
fhFcRcb = new TH2D("fhFcRcb","fc correction factor (Nb method) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; fc correction",fnPtBins,fPtBinLimits,800,0.,4.);
fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by Nb) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; corrected yield",fnPtBins,fPtBinLimits,800,0.,4.);
fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
}
fgYieldCorr = new TGraphAsymmErrors(fnPtBins+1);
if (fAsymUncertainties){
fgYieldCorrExtreme = new TGraphAsymmErrors(fnPtBins+1);
fgYieldCorrConservative = new TGraphAsymmErrors(fnPtBins+1);
fgFcConservative = new TGraphAsymmErrors(fnPtBins+1);
AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
}
double correction=0, correctionMax=0., correctionMin=0.;
fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
Double_t correctionUncStatEffc=0.;
Double_t correctionUncStatEffb=0.;
for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
Double_t frac = 1.0, errfrac =0.;
value = 0.; errvalue = 0.; errvalueMax = 0.; errvalueMin = 0.; kfactor = 0.;
errvalueExtremeMax = 0.; errvalueExtremeMin = 0.;
correction=0; correctionMax=0.; correctionMin=0.;
correctionUncStatEffc=0.; correctionUncStatEffb=0.;
if(fPbPbElossHypothesis) {
frac = fTab[0]*fNevts;
if(fIsEventPlane) frac = frac/2.0;
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)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
: 0. ;
value /= fhRECpt->GetBinWidth(ibin);
if (value<0.) value =0.;
errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
fhRECpt->GetBinError(ibin) : 0.;
errvalue /= fhRECpt->GetBinWidth(ibin);
correction = (value>0.) ?
1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) : 0. ;
if (correction<0.) correction = 0.;
kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
if (fAsymUncertainties && value>0.) {
Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
if (fgRECSystematics){
errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
}
else { errvalueMax = 0.; errvalueMin = 0.; }
Double_t errCom = ( (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]) ) ;
errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
errvalueExtremeMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
correctionMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / 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{
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]) )
) / fhRECpt->GetBinWidth(ibin);
errvalueExtremeMin = errvalueExtremeMax ;
}
fhYieldCorr->SetBinContent(ibin,value);
fhYieldCorr->SetBinError(ibin,errvalue);
if ( correction>1.0e-16 && fPbPbElossHypothesis){
for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
Double_t fcRcbvalue = 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
if(fcRcbvalue<1.e-16) fcRcbvalue=0.;
fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
Double_t correctionMaxRcb = correctionMax*rval;
Double_t correctionMinRcb = correctionMin*rval;
fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
}
}
if (fAsymUncertainties) {
Double_t x = fhYieldCorr->GetBinCenter(ibin);
fgYieldCorr->SetPoint(ibin,x,value);
fgYieldCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax);
fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
fgYieldCorrExtreme->SetPoint(ibin,x,value);
fgYieldCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax);
fgYieldCorrConservative->SetPoint(ibin,x,value);
fgYieldCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax);
if(correction>0.){
fgFcConservative->SetPoint(ibin,x,correction);
fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[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.);
}
else{
fgFcConservative->SetPoint(ibin,x,0.);
fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.);
}
}
}
}
void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
Int_t nentries = 0;
TGraphAsymmErrors *grErrFeeddown = 0;
Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
if(fFeedDownOption!=0) {
nentries = fgSigmaCorrConservative->GetN();
grErrFeeddown = new TGraphAsymmErrors(nentries);
for(Int_t i=0; i<nentries; i++) {
x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
fgSigmaCorrConservative->GetPoint(i,x,y);
if(y>0.){
errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
}
grErrFeeddown->SetPoint(i,x,0.);
grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh);
}
}
systematics->DrawErrors(grErrFeeddown);
Double_t errylcomb=0., erryhcomb=0;
for(Int_t i=1; i<nentries; i++) {
fgSigmaCorr->GetPoint(i,x,y);
errx = grErrFeeddown->GetErrorXlow(i) ;
erryl = grErrFeeddown->GetErrorYlow(i);
erryh = grErrFeeddown->GetErrorYhigh(i);
if (combineFeedDown) {
errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
} else {
errylcomb = systematics->GetTotalSystErr(x) * y ;
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);
}
}
void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
csigma->SetFillColor(0);
gPrediction->GetXaxis()->SetTitleSize(0.05);
gPrediction->GetXaxis()->SetTitleOffset(0.95);
gPrediction->GetYaxis()->SetTitleSize(0.05);
gPrediction->GetYaxis()->SetTitleOffset(0.95);
gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
gPrediction->SetLineColor(kGreen+2);
gPrediction->SetLineWidth(3);
gPrediction->SetFillColor(kGreen+1);
gPrediction->Draw("3CA");
fgSigmaCorr->SetLineColor(kRed);
fgSigmaCorr->SetLineWidth(1);
fgSigmaCorr->SetFillColor(kRed);
fgSigmaCorr->SetFillStyle(0);
fgSigmaCorr->Draw("2");
fhSigmaCorr->SetMarkerColor(kRed);
fhSigmaCorr->Draw("esame");
csigma->SetLogy();
TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
leg->SetBorderSize(0);
leg->SetLineColor(0);
leg->SetFillColor(0);
leg->SetTextFont(42);
leg->AddEntry(gPrediction,"FONLL ","fl");
leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
leg->Draw();
csigma->Draw();
}
TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
Bool_t areconsistent=kTRUE;
areconsistent &= CheckHistosConsistency(hToReweight,hReference);
if (!areconsistent) {
AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
return NULL;
}
TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
hReweighted->Reset();
Double_t weight=1.0;
Double_t yvalue=1.0;
Double_t integralRef = hReference->Integral();
Double_t integralH = hToReweight->Integral();
for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
yvalue = hToReweight->GetBinContent(i);
hReweighted->SetBinContent(i,yvalue*weight);
}
return (TH1D*)hReweighted;
}
TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
Bool_t areconsistent=kTRUE;
areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
if (!areconsistent) {
AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
return NULL;
}
TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
hReweighted->Reset();
TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
hRecReweighted->Reset();
Double_t weight=1.0;
Double_t yvalue=1.0, yrecvalue=1.0;
Double_t integralRef = hMCReference->Integral();
Double_t integralH = hMCToReweight->Integral();
for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
yvalue = hMCToReweight->GetBinContent(i);
hReweighted->SetBinContent(i,yvalue*weight);
yrecvalue = hRecToReweight->GetBinContent(i);
hRecReweighted->SetBinContent(i,yrecvalue*weight);
}
return (TH1D*)hRecReweighted;
}
Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
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 ) {
ybin =j;
break;
}
}
return ybin;
}
void AliHFPtSpectrum::ResetStatUncEff(){
Int_t entries = fhDirectEffpt->GetNbinsX();
for(Int_t i=0; i<=entries; i++){
fhDirectEffpt->SetBinError(i,0.);
fhFeedDownEffpt->SetBinError(i,0.);
}
}