From 65e55bbd03e8512ac8683d5d5ea4854910f18d0b Mon Sep 17 00:00:00 2001 From: dainese Date: Fri, 17 Sep 2010 23:04:14 +0000 Subject: [PATCH] New class to produce corrected dsigma/dpt for charm hadrons, with efficiency and B-feed-down correction (Zaida) --- PWG3/PWG3vertexingHFLinkDef.h | 1 + PWG3/libPWG3vertexingHF.pkg | 1 + PWG3/vertexingHF/AliHFPtSpectrum.cxx | 839 +++++++++++++++++++++++++ PWG3/vertexingHF/AliHFPtSpectrum.h | 164 +++++ PWG3/vertexingHF/macros/HFPtSpectrum.C | 327 ++++++++++ 5 files changed, 1332 insertions(+) create mode 100644 PWG3/vertexingHF/AliHFPtSpectrum.cxx create mode 100644 PWG3/vertexingHF/AliHFPtSpectrum.h create mode 100644 PWG3/vertexingHF/macros/HFPtSpectrum.C diff --git a/PWG3/PWG3vertexingHFLinkDef.h b/PWG3/PWG3vertexingHFLinkDef.h index a9743a4937a..9fa7d16935b 100644 --- a/PWG3/PWG3vertexingHFLinkDef.h +++ b/PWG3/PWG3vertexingHFLinkDef.h @@ -41,6 +41,7 @@ #pragma link C++ class AliMultiDimVector+; #pragma link C++ class AliSignificanceCalculator+; #pragma link C++ class AliHFMassFitter+; +#pragma link C++ class AliHFPtSpectrum+; #pragma link C++ class AliAnalysisTaskSEBkgLikeSignJPSI+; #pragma link C++ class AliAnalysisTaskSEBkgLikeSignD0+; #pragma link C++ class AliAnalysisTaskSEJPSItoEle+; diff --git a/PWG3/libPWG3vertexingHF.pkg b/PWG3/libPWG3vertexingHF.pkg index 0ef320fb762..577e52a1213 100644 --- a/PWG3/libPWG3vertexingHF.pkg +++ b/PWG3/libPWG3vertexingHF.pkg @@ -34,6 +34,7 @@ SRCS:= vertexingHF/AliAODRecoDecayHF.cxx \ vertexingHF/AliAnalysisTaskSESignificance.cxx \ vertexingHF/AliMultiDimVector.cxx vertexingHF/AliSignificanceCalculator.cxx \ vertexingHF/AliHFMassFitter.cxx \ + vertexingHF/AliHFPtSpectrum.cxx \ vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.cxx \ vertexingHF/AliAnalysisTaskSEBkgLikeSignD0.cxx \ vertexingHF/AliAnalysisTaskSEJPSItoEle.cxx \ diff --git a/PWG3/vertexingHF/AliHFPtSpectrum.cxx b/PWG3/vertexingHF/AliHFPtSpectrum.cxx new file mode 100644 index 00000000000..7cb3b2d7259 --- /dev/null +++ b/PWG3/vertexingHF/AliHFPtSpectrum.cxx @@ -0,0 +1,839 @@ +/************************************************************************** + * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//*********************************************************************** +// Class AliHFPtSpectrum +// Base class for feed-down corrections on heavy-flavour decays +// computes the cross-section via one of the three implemented methods: +// 0) Consider no feed-down prediction +// 1) Subtract the feed-down with the "fc" method +// Yield = Reco * fc; where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ; +// 2) Subtract the feed-down with the "Nb" method +// Yield = Reco - Feed-down (exact formula on the function implementation) +// +// Author: Z.Conesa, zconesa@in2p3.fr +//*********************************************************************** + +#include + +#include "TMath.h" +#include "TH1.h" +#include "TH1D.h" +#include "TGraphAsymmErrors.h" + +#include "AliLog.h" +#include "AliHFPtSpectrum.h" + +ClassImp(AliHFPtSpectrum) + +//_________________________________________________________________________________________________________ +AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option): + TNamed(name,title), + fhDirectMCpt(), + fhFeedDownMCpt(), + fhDirectMCpt_max(), + fhDirectMCpt_min(), + fhFeedDownMCpt_max(), + fhFeedDownMCpt_min(), + fhDirectEffpt(), + fhFeedDownEffpt(), + fhRECpt(), + fLuminosity(), + fTrigEfficiency(), + fhFc(), + fhFc_max(), + fhFc_min(), + fgFc(), + fhYieldCorr(), + fhYieldCorr_max(), + fhYieldCorr_min(), + fgYieldCorr(), + fhSigmaCorr(), + fhSigmaCorr_max(), + fhSigmaCorr_min(), + fgSigmaCorr(), + fFeedDownOption(option), + fAsymUncertainties(kTRUE) +{ + // + // Default constructor + // + + fLuminosity[0]=1.; fLuminosity[1]=0.; + fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.; + +} + +//_________________________________________________________________________________________________________ +AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs): + TNamed(rhs), + fhDirectMCpt(rhs.fhDirectMCpt), + fhFeedDownMCpt(rhs.fhFeedDownMCpt), + fhDirectMCpt_max(rhs.fhDirectMCpt_max), + fhDirectMCpt_min(rhs.fhDirectMCpt_min), + fhFeedDownMCpt_max(rhs.fhFeedDownMCpt_max), + fhFeedDownMCpt_min(rhs.fhFeedDownMCpt_min), + fhDirectEffpt(rhs.fhDirectEffpt), + fhFeedDownEffpt(rhs.fhFeedDownEffpt), + fhRECpt(rhs.fhRECpt), + fLuminosity(), + fTrigEfficiency(), + fhFc(rhs.fhFc), + fhFc_max(rhs.fhFc_max), + fhFc_min(rhs.fhFc_min), + fgFc(rhs.fgFc), + fhYieldCorr(rhs.fhYieldCorr), + fhYieldCorr_max(rhs.fhYieldCorr_max), + fhYieldCorr_min(rhs.fhYieldCorr_min), + fgYieldCorr(rhs.fgYieldCorr), + fhSigmaCorr(rhs.fhSigmaCorr), + fhSigmaCorr_max(rhs.fhSigmaCorr_max), + fhSigmaCorr_min(rhs.fhSigmaCorr_min), + fgSigmaCorr(rhs.fgSigmaCorr), + fFeedDownOption(rhs.fFeedDownOption), + fAsymUncertainties(rhs.fAsymUncertainties) +{ + // + // Copy constructor + // + + for(Int_t i=0; i<2; i++){ + fLuminosity[i] = rhs.fLuminosity[i]; + fTrigEfficiency[i] = rhs.fTrigEfficiency[i]; + } + +} + +//_________________________________________________________________________________________________________ +AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){ + // + // Assignment operator + // + + if (&source == this) return *this; + + fhDirectMCpt = source.fhDirectMCpt; + fhFeedDownMCpt = source.fhFeedDownMCpt; + fhDirectMCpt_max = source.fhDirectMCpt_max; + fhDirectMCpt_min = source.fhDirectMCpt_min; + fhFeedDownMCpt_max = source.fhFeedDownMCpt_max; + fhFeedDownMCpt_min = source.fhFeedDownMCpt_min; + fhDirectEffpt = source.fhDirectEffpt; + fhFeedDownEffpt = source.fhFeedDownEffpt; + fhRECpt = source.fhRECpt; + fhFc = source.fhFc; + fhFc_max = source.fhFc_max; + fhFc_min = source.fhFc_min; + fgFc = source.fgFc; + fhYieldCorr = source.fhYieldCorr; + fhYieldCorr_max = source.fhYieldCorr_max; + fhYieldCorr_min = source.fhYieldCorr_min; + fgYieldCorr = source.fgYieldCorr; + fhSigmaCorr = source.fhSigmaCorr; + fhSigmaCorr_max = source.fhSigmaCorr_max; + fhSigmaCorr_min = source.fhSigmaCorr_min; + fgSigmaCorr = source.fgSigmaCorr; + fFeedDownOption = source.fFeedDownOption; + fAsymUncertainties = source.fAsymUncertainties; + + for(Int_t i=0; i<2; i++){ + fLuminosity[i] = source.fLuminosity[i]; + fTrigEfficiency[i] = source.fTrigEfficiency[i]; + } + + return *this; +} + +//_________________________________________________________________________________________________________ +AliHFPtSpectrum::~AliHFPtSpectrum(){ + // + // Destructor + // + ; +} + + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown){ + // + // Set the MonteCarlo or Theoretical spectra + // both for direct and feed-down contributions + // + + if (!hDirect || !hFeedDown) { + AliError("One or both (direct, feed-down) 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 = hDirect; + fhFeedDownMCpt = hFeedDown; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1 *hFeedDown){ + // + // Set the MonteCarlo or Theoretical spectra + // for feed-down contribution + // + + if (!hFeedDown) { + AliError("Feed-down spectra don't exist"); + return; + } + fhFeedDownMCpt = hFeedDown; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1 *hDirectMax, TH1 *hDirectMin, TH1 *hFeedDownMax, TH1 *hFeedDownMin){ + // + // Set the maximum and minimum MonteCarlo or Theoretical spectra + // both for direct and feed-down contributions + // used in case uncertainties are asymmetric and ca not be on the "basic histograms" + // + + if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin) { + AliError("One or all of the max/min direct/feed-down 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; + } + + fhDirectMCpt_max = hDirectMax; + fhDirectMCpt_min = hDirectMin; + fhFeedDownMCpt_max = hFeedDownMax; + fhFeedDownMCpt_min = hFeedDownMin; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1 *hFeedDownMax, TH1 *hFeedDownMin){ + // + // Set the maximum and minimum MonteCarlo or Theoretical spectra + // for feed-down contributions + // used in case uncertainties are asymmetric and can not be on the "basic histogram" + // + + if (!hFeedDownMax || !hFeedDownMin) { + 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; + } + + fhFeedDownMCpt_max = hFeedDownMax; + fhFeedDownMCpt_min = hFeedDownMin; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1 *hDirectEff){ + // + // Set the Acceptance and Efficiency corrections + // for the direct contribution + // + + if (!hDirectEff) { + AliError("The direct acceptance and efficiency corrections doesn't exist"); + return; + } + + fhDirectEffpt = hDirectEff; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::SetAccEffCorrection(TH1 *hDirectEff, TH1 *hFeedDownEff){ + // + // Set the Acceptance and Efficiency corrections + // both for direct and feed-down contributions + // + + 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 = hDirectEff; + fhFeedDownEffpt = hFeedDownEff; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::SetReconstructedSpectrum(TH1 *hRec) { + // + // Set the reconstructed spectrum + // + + if (!hRec) { + AliError("The reconstructed spectrum doesn't exist"); + return; + } + + fhRECpt = hRec; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t delta_y, Double_t BR_c, Double_t BR_b){ + // + // Main function to compute the corrected cross-section: + // + // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down) + // + // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 ) + + // + // First: Initialization + // + Bool_t areHistosOk = Initialize(); + if (!areHistosOk) { + AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?"); + return; + } + + // + // Second: Correct for feed-down + // + if (fFeedDownOption==1) { + // Compute the feed-down correction via fc-method + CalculateFeedDownCorrection_fc(); + // Correct the yield for feed-down correction via fc-method + CalculateFeedDownCorrectedSpectrum_fc(); + } + else if (fFeedDownOption==2) { + // Correct the yield for feed-down correction via Nb-method + CalculateFeedDownCorrectedSpectrum_Nb(delta_y,BR_b); + } + else if (fFeedDownOption==0) { + // If there is no need for feed-down correction, + // the "corrected" yield is equal to the raw yield + fhYieldCorr = fhRECpt; + fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield"); + fhYieldCorr_max = fhRECpt; + fhYieldCorr_min = fhRECpt; + fhYieldCorr_max->SetNameTitle("fhYieldCorr_max","un-corrected yield"); + fhYieldCorr_min->SetNameTitle("fhYieldCorr_min","un-corrected yield"); + fAsymUncertainties=kFALSE; + } + else { + AliInfo(" Are you sure the feed-down correction option is right ?"); + } + + // 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],delta_y,BR_c,BR_b); + + // + // Finally: Correct from yields to cross-section + // + Int_t nbins = fhRECpt->GetNbinsX(); + Double_t binwidth = fhRECpt->GetBinWidth(1); + Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ; + Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ; + + // declare the output histograms + TH1D *hSigmaCorr = new TH1D("hSigmaCorr","corrected sigma",nbins,xmin,xmax); + TH1D *hSigmaCorr_max = new TH1D("hSigmaCorr_max","max corrected sigma",nbins,xmin,xmax); + TH1D *hSigmaCorr_min = new TH1D("hSigmaCorr_min","min corrected sigma",nbins,xmin,xmax); + // and the output TGraphAsymmErrors + if (fAsymUncertainties & !fgSigmaCorr) fgSigmaCorr = new TGraphAsymmErrors(nbins); + + // protect against null denominator + if (delta_y==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || BR_c==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, value_max=0., value_min=0.; + for(Int_t ibin=0; ibin<=nbins; ibin++){ + + // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down) + value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? + ( fhYieldCorr->GetBinContent(ibin) / ( delta_y * BR_c * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) ) + : 0. ; + + // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 ) + if (fAsymUncertainties) { + value_max = 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]) ); + value_min = 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]) ); + } + else { + // protect against null denominator + value_max = (value!=0.) ? + value * TMath::Sqrt( (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))* (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) + + (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) + + (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) ) + : 0. ; + value_min = value_max; + } + + // Fill the histograms + hSigmaCorr->SetBinContent(ibin,value); + hSigmaCorr_max->SetBinContent(ibin,value_max); + hSigmaCorr_min->SetBinContent(ibin,value_min); + // Fill the TGraphAsymmErrors + if (fAsymUncertainties) { + Double_t x = fhYieldCorr->GetBinCenter(ibin); + fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y + fgSigmaCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),value_min,value_max); // i,xl,xh,yl,yh + } + + } + + fhSigmaCorr = hSigmaCorr ; + fhSigmaCorr_max = hSigmaCorr_max ; + fhSigmaCorr_min = hSigmaCorr_min ; +} + +//_________________________________________________________________________________________________________ +Bool_t AliHFPtSpectrum::Initialize(){ + // + // Initialization of the variables (histograms) + // + + 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; } + + // Start checking the input histograms consistency + Bool_t areconsistent=kTRUE; + + // General checks + 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; + + // + // Common checks for options 1 (fc) & 2(Nb) + 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 (!fhFeedDownMCpt_max || !fhFeedDownMCpt_min) { + AliError(" Max/Min theoretical Nb distributions are not defined"); + return kFALSE; + } + areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCpt_max); + } + if (!areconsistent) { + AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)"); + return kFALSE; + } + if (fFeedDownOption>1) return kTRUE; + + // + // Now checks for option 1 (fc correction) + if (!fhDirectMCpt) { + AliError("Theoretical Nc distributions is not defined"); + return kFALSE; + } + areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt); + areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt); + if (fAsymUncertainties) { + if (!fhDirectMCpt_max || !fhDirectMCpt_min) { + AliError(" Max/Min theoretical Nc distributions are not defined"); + return kFALSE; + } + areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCpt_max); + } + if (!areconsistent) { + AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)"); + return kFALSE; + } + + return kTRUE; +} + +//_________________________________________________________________________________________________________ +Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1 *h1, TH1 *h2){ + // + // Check the histograms consistency (bins, limits) + // + + 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 max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ; + Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ; +// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ; + + if (binwidth1!=binwidth2) { + AliInfo(" histograms with different bin width"); + return kFALSE; + } + if (min1!=min2) { + AliInfo(" histograms with different minimum"); + return kFALSE; + } +// if (max1!=max2) { +// AliInfo(" histograms with different maximum"); +// return kFALSE; +// } + + return kTRUE; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::CalculateFeedDownCorrection_fc(){ + // + // Compute fc factor and its uncertainties bin by bin + // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) + // + + // define the variables + Int_t nbins = fhRECpt->GetNbinsX(); + Double_t binwidth = fhRECpt->GetBinWidth(1); + Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ; + Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ; + Double_t correction=1.; + Double_t correction_max=1., correction_min=1.; + Double_t theory_ratio=1.; + Double_t eff_ratio=1.; + + // declare the output histograms + TH1D *hfc = new TH1D("hfc","fc correction factor",nbins,xmin,xmax); + TH1D *hfc_max = new TH1D("hfc_max","max fc correction factor",nbins,xmin,xmax); + TH1D *hfc_min = new TH1D("hfc_min","min fc correction factor",nbins,xmin,xmax); + // two local control histograms + TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax); + TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax); + // and the output TGraphAsymmErrors + if (fAsymUncertainties & !fgFc) fgFc = new TGraphAsymmErrors(nbins); + + // + // Compute fc + // + for (Int_t ibin=0; ibin<=nbins; ibin++) { + + // theory_ratio = (N_b/N_c) + theory_ratio = (fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0.) ? fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ; + // eff_ratio = (eff_b/eff_c) + eff_ratio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ; + // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) + correction = (eff_ratio && theory_ratio) ? ( 1. / ( 1 + ( eff_ratio * theory_ratio ) ) ) : 1.0 ; + + // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ] + // delta_fc = fc^2 * (Nb/Nc) * sqrt ( (delta_Nb/Nb)^2 + (delta_Nc/Nc)^2 ) + Double_t delta_Nb_max = fhFeedDownMCpt_max->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin) ; + Double_t delta_Nb_min = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCpt_min->GetBinContent(ibin) ; + Double_t delta_Nc_max = fhDirectMCpt_max->GetBinContent(ibin) - fhDirectMCpt->GetBinContent(ibin) ; + Double_t delta_Nc_min = fhDirectMCpt->GetBinContent(ibin) - fhDirectMCpt_min->GetBinContent(ibin) ; + + // Protect against null denominator. If so, define uncertainty as null + if (fhFeedDownMCpt->GetBinContent(ibin) && fhFeedDownMCpt->GetBinContent(ibin)!=0. && + fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0. ) { + correction_max = correction*correction*theory_ratio * + TMath::Sqrt( + (delta_Nb_max/fhFeedDownMCpt->GetBinContent(ibin))*(delta_Nb_max/fhFeedDownMCpt->GetBinContent(ibin)) + + (delta_Nc_max/fhDirectMCpt->GetBinContent(ibin))*(delta_Nc_max/fhDirectMCpt->GetBinContent(ibin)) + ); + correction_min = correction*correction*theory_ratio * + TMath::Sqrt( + (delta_Nb_min/fhFeedDownMCpt->GetBinContent(ibin))*(delta_Nb_min/fhFeedDownMCpt->GetBinContent(ibin)) + + (delta_Nc_min/fhDirectMCpt->GetBinContent(ibin))*(delta_Nc_min/fhDirectMCpt->GetBinContent(ibin)) + ); + } + else { correction_max = 0.; correction_min = 0.; } + + + // Fill in the histograms + hTheoryRatio->SetBinContent(ibin,theory_ratio); + hEffRatio->SetBinContent(ibin,eff_ratio); + hfc->SetBinContent(ibin,correction); + hfc_max->SetBinContent(ibin,correction+correction_max); + hfc_min->SetBinContent(ibin,correction-correction_min); + if (fAsymUncertainties) { + Double_t x = fhDirectMCpt->GetBinCenter(ibin); + fgFc->SetPoint(ibin,x,correction); // i,x,y + fgFc->SetPointError(ibin,(binwidth/2.),(binwidth/2.),correction_min,correction_max); // i,xl,xh,yl,yh + } + + } + + fhFc = hfc; + fhFc_max = hfc_max; + fhFc_min = hfc_min; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrum_fc(){ + // + // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin) + // physics = reco * fc + // + // uncertainty: delta_physics = physics * sqrt ( (delta_reco/reco)^2 + (delta_fc/fc)^2 ) + // + // ( Calculation done bin by bin ) + + if (!fhFc || !fhRECpt) { + AliError(" Reconstructed or fc distributions are not defined"); + return; + } + + Int_t nbins = fhRECpt->GetNbinsX(); + Double_t value = 0., value_dmax= 0., value_dmin= 0.; + Double_t binwidth = fhRECpt->GetBinWidth(1); + Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ; + Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ; + + // declare the output histograms + TH1D *hYield = new TH1D("hYield","corrected yield (by fc)",nbins,xmin,xmax); + TH1D *hYield_max = new TH1D("hYield_max","max corrected yield (by fc)",nbins,xmin,xmax); + TH1D *hYield_min = new TH1D("hYield_min","min corrected yield (by fc)",nbins,xmin,xmax); + // and the output TGraphAsymmErrors + if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins); + + // + // Do the calculation + // + for (Int_t ibin=0; ibin<=nbins; ibin++) { + + // calculate the value + value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ; + + // calculate the value uncertainty + // Protect against null denominator. If so, define uncertainty as null + if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) { + + if (fAsymUncertainties) { + + if (fhFc->GetBinContent(ibin) && fhFc->GetBinContent(ibin)!=0.) { + value_dmax = value * TMath::Sqrt( ( (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin))*(fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ) + ( (fgFc->GetErrorYhigh(ibin)/fhFc->GetBinContent(ibin))*(fgFc->GetErrorYhigh(ibin)/fhFc->GetBinContent(ibin)) ) ); + value_dmin = value * TMath::Sqrt( ( (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin))*(fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ) + ( (fgFc->GetErrorYlow(ibin)/fhFc->GetBinContent(ibin))*(fgFc->GetErrorYlow(ibin)/fhFc->GetBinContent(ibin)) ) ); + } + else { value_dmax = 0.; value_dmin = 0.; } + + } + else { // Don't consider fc uncertainty in this case [ to be tested!!! ] + value_dmax = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ; + value_dmin = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ; + } + + } + else { value_dmax = 0.; value_dmin = 0.; } + + // fill in the histograms + hYield->SetBinContent(ibin,value); + hYield_max->SetBinContent(ibin,value+value_dmax); + hYield_min->SetBinContent(ibin,value-value_dmin); + if (fAsymUncertainties) { + Double_t center = hYield->GetBinCenter(ibin); + fgYieldCorr->SetPoint(ibin,center,value); // i,x,y + fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),value_dmin,value_dmax); // i,xl,xh,yl,yh + } + + } + + fhYieldCorr = hYield; + fhYieldCorr_max = hYield_max; + fhYieldCorr_min = hYield_min; +} + +//_________________________________________________________________________________________________________ +void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrum_Nb(Float_t delta_y, Double_t BR_b){ + // + // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin) + // physics = reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) + // + // uncertainty: delta_physics = sqrt ( (delta_reco)^2 + (k*delta_lumi/lumi)^2 + + // (k*delta_eff_trig/eff_trig)^2 + (k*delta_Nb/Nb)^2 ) + // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th + // + + Int_t nbins = fhRECpt->GetNbinsX(); + Double_t binwidth = fhRECpt->GetBinWidth(1); + Double_t value = 0., value_dmax= 0., value_dmin= 0., kfactor=0.; + Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ; + Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ; + + // declare the output histograms + TH1D *hYield = new TH1D("hYield","corrected yield (by Nb)",nbins,xmin,xmax); + TH1D *hYield_max = new TH1D("hYield_max","max corrected yield (by Nb)",nbins,xmin,xmax); + TH1D *hYield_min = new TH1D("hYield_min","min corrected yield (by Nb)",nbins,xmin,xmax); + // and the output TGraphAsymmErrors + if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins); + + // + // Do the calculation + // + for (Int_t ibin=0; ibin<=nbins; ibin++) { + + // calculate the value + value = fhRECpt->GetBinContent(ibin) - (delta_y*BR_b*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ); + + kfactor = delta_y*BR_b*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ; + + // calculate the value uncertainty + if (fAsymUncertainties) { + Double_t Nb = fhFeedDownMCpt->GetBinContent(ibin); + Double_t Nb_dmax = fhFeedDownMCpt_max->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin); + Double_t Nb_dmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCpt_min->GetBinContent(ibin); + value_dmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) + + ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) + + ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) + + ( (kfactor*Nb_dmax/Nb)*(kfactor*Nb_dmax/Nb) ) ); + value_dmin = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) + + ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) + + ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) + + ( (kfactor*Nb_dmin/Nb)*(kfactor*Nb_dmin/Nb) ) ); + } + else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ] + value_dmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) + + ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) + + ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) ); + value_dmin = value_dmax ; + } + + // fill in histograms + hYield->SetBinContent(ibin,value); + hYield_max->SetBinContent(ibin,value+value_dmax); + hYield_min->SetBinContent(ibin,value-value_dmin); + if (fAsymUncertainties) { + Double_t x = hYield->GetBinCenter(ibin); + fgYieldCorr->SetPoint(ibin,x,value); // i,x,y + fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),value_dmin,value_dmax); // i,xl,xh,yl,yh + } + + } + + fhYieldCorr = hYield; + fhYieldCorr_max = hYield_max; + fhYieldCorr_min = hYield_min; +} + + +//_________________________________________________________________________________________________________ +TH1 * AliHFPtSpectrum::ReweightHisto(TH1 *hToReweight, TH1 *hReference){ + // + // Function to reweight histograms for testing purposes: + // This function takes the histo hToReweight and reweights + // it (its pt shape) with respect to hReference + // + + // check histograms consistency + Bool_t areconsistent=kTRUE; + areconsistent &= CheckHistosConsistency(hToReweight,hReference); + if (!areconsistent) { + AliInfo("the histograms to reweight are not consistent (bin width, bounds)"); + return NULL; + } + + // define a new empty histogram + TH1 *hReweighted = (TH1*)hToReweight->Clone("hReweighted"); + hReweighted->Reset(); + Double_t weight=1.0; + Double_t yvalue=1.0; + Double_t integral_ref = hReference->Integral(); + Double_t integral_h = hToReweight->Integral(); + + // now reweight the spectra + // + // the weight is the relative probability of the given pt bin in the reference histo + // divided by its relative probability (to normalize it) on the histo to re-weight + for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) { + weight = (hReference->GetBinContent(i)/integral_ref) / (hToReweight->GetBinContent(i)/integral_h) ; + yvalue = hToReweight->GetBinContent(i); + hReweighted->SetBinContent(i,yvalue*weight); + } + + return (TH1*)hReweighted; +} + +//_________________________________________________________________________________________________________ +TH1 * AliHFPtSpectrum::ReweightRecHisto(TH1 *hRecToReweight, TH1 *hMCToReweight, TH1 *hMCReference){ + // + // Function to reweight histograms for testing purposes: + // This function takes the histo hToReweight and reweights + // it (its pt shape) with respect to hReference /hMCToReweight + // + + // check histograms consistency + 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; + } + + // define a new empty histogram + TH1 *hReweighted = (TH1*)hMCToReweight->Clone("hReweighted"); + hReweighted->Reset(); + TH1 *hRecReweighted = (TH1*)hRecToReweight->Clone("hRecReweighted"); + hRecReweighted->Reset(); + Double_t weight=1.0; + Double_t yvalue=1.0, yrecvalue=1.0; + Double_t integral_ref = hMCReference->Integral(); + Double_t integral_h = hMCToReweight->Integral(); + + // now reweight the spectra + // + // the weight is the relative probability of the given pt bin + // that should be applied in the MC histo to get the reference histo shape + // Probabilities are properly normalized. + for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) { + weight = (hMCReference->GetBinContent(i)/integral_ref) / (hMCToReweight->GetBinContent(i)/integral_h) ; + yvalue = hMCToReweight->GetBinContent(i); + hReweighted->SetBinContent(i,yvalue*weight); + yrecvalue = hRecToReweight->GetBinContent(i); + hRecReweighted->SetBinContent(i,yrecvalue*weight); + } + + return (TH1*)hRecReweighted; +} + diff --git a/PWG3/vertexingHF/AliHFPtSpectrum.h b/PWG3/vertexingHF/AliHFPtSpectrum.h new file mode 100644 index 00000000000..e8c0f0c6df1 --- /dev/null +++ b/PWG3/vertexingHF/AliHFPtSpectrum.h @@ -0,0 +1,164 @@ +#ifndef ALIHFPTSPECTRUM_H +#define ALIHFPTSPECTRUM_H + +/* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//*********************************************************************** +// Class AliHFPtSpectrum +// Base class for feed-down corrections on heavy-flavour decays +// computes the cross-section via one of the three implemented methods: +// 0) Consider no feed-down prediction +// 1) Subtract the feed-down with the "fc" method +// Yield = Reco * fc; where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ; +// 2) Subtract the feed-down with the "Nb" method +// Yield = Reco - Feed-down (exact formula on the function implementation) +// +// Author: Z.Conesa, zconesa@in2p3.fr +//*********************************************************************** + +#include "TNamed.h" +#include "TH1.h" +#include "TGraphAsymmErrors.h" + +class AliHFPtSpectrum: public TNamed +{ + + public: + + // Constructor + AliHFPtSpectrum(const char* name="AliHFPtSpectrum", const char* title="HF feed down correction class", Int_t option=1); + // Copy constructor + AliHFPtSpectrum(const AliHFPtSpectrum &rhs); + // Assignment operator + AliHFPtSpectrum& operator=(const AliHFPtSpectrum &source); + // Destructor + virtual ~AliHFPtSpectrum(); + + // + // Setters + // + // Set the theoretical direct & feeddown pt spectrum + void SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown); + // Set the theoretical feeddown pt spectrum + void SetFeedDownMCptSpectra(TH1 *hFeedDown); + // Set the theoretical direct & feeddown pt spectrum upper and lower bounds + void SetMCptDistributionsBounds(TH1 *hDirectMax, TH1 *hDirectMin, TH1 *hFeedDownMax, TH1 *hFeedDownMin); + // Set the theoretical feeddown pt spectrum upper and lower bounds + void SetFeedDownMCptDistributionsBounds(TH1 *hFeedDownMax, TH1 *hFeedDownMin); + // Set the acceptance and efficiency corrections for direct + void SetDirectAccEffCorrection(TH1 *hDirectEff); + // Set the acceptance and efficiency corrections for direct & feeddown + void SetAccEffCorrection(TH1 *hDirectEff, TH1 *hFeedDownEff); + // Set the reconstructed spectrum + void SetReconstructedSpectrum(TH1 *hRec); + // Set the calculation option flag for feed-down correction: 0=none, 1=fc , 2=Nb + void SetFeedDownCalculationOption(Int_t option){ fFeedDownOption = option; } + // Set if the calculation has to consider asymmetric uncertaInt_ties or not + void SetComputeAsymmetricUncertainties(bool flag){ fAsymUncertainties = flag; } + // Set the luminosity and its uncertainty + void SetLuminosity(Double_t luminosity, Double_t unc){ + fLuminosity[0]=luminosity; fLuminosity[1]=unc; + } + // Set the trigger efficiency and its uncertainty + void SetTriggerEfficiency(Double_t efficiency, Double_t unc){ + fTrigEfficiency[0]=efficiency; fTrigEfficiency[1]=unc; + } + + // + // Getters + // + // Return the TGraphAsymmErrors of the feed-down correction + TGraphAsymmErrors * GetFeedDownCorrection_fc() { return (fgFc ? fgFc : NULL); } + // Return the histogram of the feed-down correction + TH1 * GetHistoFeedDownCorrection_fc() { return (fhFc ? fhFc : NULL); } + // Return the histograms of the feed-down correction bounds + TH1 * GetHistoUpperLimitFeedDownCorrection_fc() { return (fhFc_max ? fhFc_max : NULL); } + TH1 * GetHistoLowerLimitFeedDownCorrection_fc() { return (fhFc_min ? fhFc_min : NULL); } + // Return the TGraphAsymmErrors of the yield after feed-down correction + TGraphAsymmErrors * GetFeedDownCorrectedSpectrum() { return (fgYieldCorr ? fgYieldCorr : NULL); } + // Return the histogram of the yield after feed-down correction + TH1 * GetHistoFeedDownCorrectedSpectrum() { return (fhYieldCorr ? fhYieldCorr : NULL); } + // Return the histogram of the yield after feed-down correction bounds + TH1 * GetHistoUpperLimitFeedDownCorrectedSpectrum() { return (fhYieldCorr_max ? fhYieldCorr_max : NULL); } + TH1 * GetHistoLowerLimitFeedDownCorrectedSpectrum() { return (fhYieldCorr_min ? fhYieldCorr_min : NULL); } + // Return the equivalent invariant cross-section TGraphAsymmErrors + TGraphAsymmErrors * GetCrossSectionFromYieldSpectrum() { return (fgSigmaCorr ? fgSigmaCorr : NULL); } + // Return the equivalent invariant cross-section histogram + TH1 * GetHistoCrossSectionFromYieldSpectrum() { return (fhSigmaCorr ? fhSigmaCorr : NULL); } + // Return the equivalent invariant cross-section histogram bounds + TH1 * GetHistoUpperLimitCrossSectionFromYieldSpectrum() { return (fhSigmaCorr_max ? fhSigmaCorr_max : NULL); } + TH1 * GetHistoLowerLimitCrossSectionFromYieldSpectrum() { return (fhSigmaCorr_min ? fhSigmaCorr_min : NULL); } + + // + // Main function: + // Compute the invariant cross-section from the yield (correct it) + void ComputeHFPtSpectrum(Double_t delta_y=1.0, Double_t BR_c=1.0, Double_t BR_b_decay=1.0); + + // + // Functions to reweight histograms for testing purposes: + // to reweight the simulation: hToReweight is reweighted as hReference/hToReweight + TH1 * ReweightHisto(TH1 *hToReweight, TH1 *hReference); + // to reweight the reco-histos: hRecToReweight is reweighted as hReference/hMCToReweight + TH1 * ReweightRecHisto(TH1 *hRecToReweight, TH1 *hMCToReweight, TH1 *hMCReference); + + + protected: + + // Initialization + Bool_t Initialize(); + + // Basic functions + // + // Compute the feed-down correction via fc-method + void CalculateFeedDownCorrection_fc(); + // Correct the yield for feed-down correction via fc-method + void CalculateFeedDownCorrectedSpectrum_fc(); + // Correct the yield for feed-down correction via Nb-method + void CalculateFeedDownCorrectedSpectrum_Nb(Float_t delta_y, Double_t BR_b_decay); + + // Check histograms consistency function + Bool_t CheckHistosConsistency(TH1 *h1, TH1 *h2); + + // + // Input spectra + // + TH1 *fhDirectMCpt; // Input MC c-->D spectra + TH1 *fhFeedDownMCpt; // Input MC b-->D spectra + TH1 *fhDirectMCpt_max; // Input MC maximum c-->D spectra + TH1 *fhDirectMCpt_min; // Input MC minimum c-->D spectra + TH1 *fhFeedDownMCpt_max; // Input MC maximum b-->D spectra + TH1 *fhFeedDownMCpt_min; // Input MC minimum b-->D spectra + TH1 *fhDirectEffpt; // c-->D Acceptance and efficiency correction + TH1 *fhFeedDownEffpt; // b-->D Acceptance and efficiency correction + TH1 *fhRECpt; // all reconstructed D + // + // Normalization factors + Double_t fLuminosity[2]; // analyzed luminosity & uncertainty + Double_t fTrigEfficiency[2]; // trigger efficiency & uncertainty + + // + // Output spectra + // + TH1 *fhFc; // Correction histo fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) + TH1 *fhFc_max; // Maximum fc histo + TH1 *fhFc_min; // Minimum fc histo + TGraphAsymmErrors * fgFc; // Correction as TGraphAsymmErrors + TH1 *fhYieldCorr; // Corrected yield + TH1 *fhYieldCorr_max; // Maximum corrected yield + TH1 *fhYieldCorr_min; // Minimum corrected yield + TGraphAsymmErrors * fgYieldCorr; // Corrected yield as TGraphAsymmErrors + TH1 *fhSigmaCorr; // Corrected cross-section + TH1 *fhSigmaCorr_max; // Maximum corrected cross-section + TH1 *fhSigmaCorr_min; // Minimum corrected cross-section + TGraphAsymmErrors * fgSigmaCorr; // Corrected cross-section as TGraphAsymmErrors + + // + 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 + + + ClassDef(AliHFPtSpectrum,1) // Class for Heavy Flavor spectra corrections +}; + +#endif diff --git a/PWG3/vertexingHF/macros/HFPtSpectrum.C b/PWG3/vertexingHF/macros/HFPtSpectrum.C new file mode 100644 index 00000000000..b8f309f7869 --- /dev/null +++ b/PWG3/vertexingHF/macros/HFPtSpectrum.C @@ -0,0 +1,327 @@ +// +// Macro to use the AliHFPtSpectrum class +// to compute the feed-down corrections for heavy-flavor +// +// Z.Conesa, September 2010 (zconesa@in2p3.fr) +// + +#include +#include "TH1D.h" +#include "TH1.h" +#include "TH2F.h" +#include "TFile.h" +#include "TGraphAsymmErrors.h" +#include "TCanvas.h" + +#include "AliHFPtSpectrum.h" + +// +// Macro execution parameters: +// 0) filename with the theoretical predictions (direct & feed-down) +// 1) acceptance and reconstruction efficiencies file name (direct & feed-down) +// 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 +// +void HFPtSpectrum(const char *mcfilename="FeedDownCorrectionMC.root", + const char *efffilename="Efficiencies.root", + const char *recofilename="Reconstructed.root", + const char *outfilename="HFPtSpectrum.root", + int option=1, double lumi=1.0, double eff_trig=1.0){ + + // Set if calculation considers asymmetric uncertainties or not + bool asym = true; + + // Set the meson and decay + // (only D0 -> K pi & D+-->K pi pi implemented here) + bool isD0Kpi = true; + bool isDplusKpipi = false; + if (isD0Kpi && isDplusKpipi) { + cout << "Sorry, can not deal with the two corrections at the same time"<2) { + cout<< "Bad calculation option, should be <=2"<D spectra + TH1D *hFeedDownMCpt; // Input MC b-->D spectra + TH1D *hDirectMCpt_max; // Input MC maximum c-->D spectra + TH1D *hDirectMCpt_min; // Input MC minimum c-->D spectra + TH1D *hFeedDownMCpt_max; // Input MC maximum b-->D spectra + TH1D *hFeedDownMCpt_min; // Input MC minimum b-->D spectra + TH1D *hDirectEffpt; // c-->D Acceptance and efficiency correction + TH1D *hFeedDownEffpt; // b-->D Acceptance and efficiency correction + TH1D *hRECpt; // all reconstructed D + + // + // Define/Get the input histograms + // + TFile * mcfile = new TFile(mcfilename,"read"); + if (isD0Kpi){ + hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central"); + hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central"); + hDirectMCpt_max = (TH1D*)mcfile->Get("hD0Kpipred_max"); + hDirectMCpt_min = (TH1D*)mcfile->Get("hD0Kpipred_min"); + hFeedDownMCpt_max = (TH1D*)mcfile->Get("hD0KpifromBpred_max"); + hFeedDownMCpt_min = (TH1D*)mcfile->Get("hD0KpifromBpred_min"); + } + else if (isDplusKpipi){ + hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central"); + hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central"); + hDirectMCpt_max = (TH1D*)mcfile->Get("hDpluskpipipred_max"); + hDirectMCpt_min = (TH1D*)mcfile->Get("hDpluskpipipred_min"); + hFeedDownMCpt_max = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max"); + hFeedDownMCpt_min = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min"); + } + // + hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra"); + hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra"); + hDirectMCpt_max->SetNameTitle("hDirectMCpt_max","max direct MC spectra"); + hDirectMCpt_min->SetNameTitle("hDirectMCpt_min","min direct MC spectra"); + hFeedDownMCpt_max->SetNameTitle("hFeedDownMCpt_max","max feed-down MC spectra"); + hFeedDownMCpt_min->SetNameTitle("hFeedDownMCpt_min","min feed-down MC spectra"); + // + TFile * efffile = new TFile(efffilename,"read"); + hDirectEffpt = (TH1D*)efffile->Get("hDirectEffpt"); + hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff"); + hFeedDownEffpt = (TH1D*)efffile->Get("hFeedDownEffpt"); + hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff"); + // + TFile * recofile = new TFile(recofilename,"read"); + hRECpt = (TH1D*)recofile->Get("hRecoAll"); + hRECpt->SetNameTitle("hRECpt","Reconstructed spectra"); + + // + // Define the output histograms + // + TFile *out = new TFile(outfilename,"recreate"); + // + // + TH1D *histofc = (TH1D*)hRECpt->Clone(); + histofc->SetNameTitle("histofc","fc correction factor"); + histofc->Reset(); + TH1D *histofc_max = (TH1D*)histofc->Clone(); + TH1D *histofc_min = (TH1D*)histofc->Clone(); + TH1D *histoYieldCorr = (TH1D*)hRECpt->Clone(); + histoYieldCorr->SetNameTitle("histoYieldCorrFc","corrected yield"); + histoYieldCorr->Reset(); + TH1D *histoYieldCorr_max = (TH1D*)histoYieldCorr->Clone(); + histoYieldCorr_max->SetNameTitle("histoYieldCorr_max","max corrected yield"); + TH1D *histoYieldCorr_min = (TH1D*)histoYieldCorr->Clone(); + histoYieldCorr_min->SetNameTitle("histoYieldCorr_min","min corrected yield"); + TH1D *histoSigmaCorr = (TH1D*)hRECpt->Clone(); + histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section"); + histoSigmaCorr->Reset(); + // + int nbins = hRECpt->GetNbinsX(); + TGraphAsymmErrors * gFc = new TGraphAsymmErrors(nbins); + TGraphAsymmErrors * gYieldCorr = new TGraphAsymmErrors(nbins); + TGraphAsymmErrors * gSigmaCorr = new TGraphAsymmErrors(nbins); + + + // + // Main functionalities for the calculation + // + + // Define and set the basic option flags + AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option); + spectra->SetFeedDownCalculationOption(option); + spectra->SetComputeAsymmetricUncertainties(asym); + + // Feed the input histograms + spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt); + spectra->SetReconstructedSpectrum(hRECpt); + // option specific histos + if(option==1){ + spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt); + if(asym) spectra->SetMCptDistributionsBounds(hDirectMCpt_max,hDirectMCpt_min,hFeedDownMCpt_max,hFeedDownMCpt_min); + } + else if(option==2){ + spectra->SetFeedDownMCptSpectra(hFeedDownMCpt); + if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCpt_max,hFeedDownMCpt_min); + } + + // Set normalization factors (uncertainties set to 0. as example) + spectra->SetLuminosity(lumi,0.); + spectra->SetTriggerEfficiency(eff_trig,0.); + + // Do the calculations + cout << " Doing the calculation... "<< endl; + double delta_y = 1.0; + double BR_c=1.0, BR_b_decay=1.0; + spectra->ComputeHFPtSpectrum(delta_y,BR_c,BR_b_decay); + cout << " ended the calculation, getting the histograms back " << endl; + + // + // Get the output histograms + // + // the corrected yield and cross-section + histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum(); + histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum(); + histoYieldCorr_max = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum(); + histoYieldCorr_min = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum(); + if (asym) { + gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum(); + gYieldCorr = spectra->GetFeedDownCorrectedSpectrum(); + } + + if (option==0 && asym){ + gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)"); + gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)"); + } + if (option==1){ + // fc histos + histofc = (TH1D*)spectra->GetHistoFeedDownCorrection_fc(); + histofc_max = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrection_fc(); + histofc_min = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrection_fc(); + histofc_max->SetNameTitle("hfc_max","max fc correction factor"); + histofc_min->SetNameTitle("histofc_min","min fc correction factor"); + if (asym) { + gFc = spectra->GetFeedDownCorrection_fc(); + gFc->SetNameTitle("gFc","gFc"); + gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)"); + gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)"); + } + } + if (option==2 && asym) { + gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)"); + gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)"); + } + + // + // Now, plot the results ! :) + // + + cout << " Drawing the results ! " << endl; + + if (option==1) { + + TCanvas * cfc = new TCanvas("cfc","Fc"); + histofc_max->Draw("c"); + histofc->Draw("csame"); + histofc_min->Draw("csame"); + cfc->Update(); + + if (asym) { + TH2F *histofc_draw= new TH2F("histofc_draw","histofc (for drawing)",100,0,33.25,100,0.01,1.25); + histofc_draw->SetStats(0); + histofc_draw->GetXaxis()->SetTitle("p_{T} [GeV]"); + histofc_draw ->GetXaxis()->SetTitleSize(0.05); + histofc_draw->GetXaxis()->SetTitleOffset(0.95); + histofc_draw->GetYaxis()->SetTitle(" fc "); + histofc_draw->GetYaxis()->SetTitleSize(0.05); + TCanvas *cfc_asym = new TCanvas("cfc_asym","Asymmetric fc (TGraphAsymmErr)"); + gFc->SetFillStyle(3001); + gFc->SetLineWidth(3); + gFc->SetMarkerStyle(20); + gFc->SetFillColor(3); + histofc_draw->Draw(); + gFc->Draw("3Csame"); + gFc->Draw("Xsame"); + cfc_asym->Update(); + } + + } + + // + // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra) + // + TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma"); + hDirectMCpt->SetMarkerStyle(20); + hDirectMCpt->SetMarkerColor(4); + hDirectMCpt->Draw("p"); + histoSigmaCorr->SetMarkerStyle(21); + histoSigmaCorr->SetMarkerColor(2); + histoSigmaCorr->Draw("psame"); + histoYieldCorr->SetMarkerStyle(22); + histoYieldCorr->SetMarkerColor(6); + histoYieldCorr->Draw("psame"); + hRECpt->SetMarkerStyle(23); + hRECpt->SetMarkerColor(3); + hRECpt->Draw("psame"); + cresult->SetLogy(); + cresult->Update(); + + if (asym) { + + TH2F *histo_draw = new TH2F("histo_draw","histo (for drawing)",100,0,33.25,100,50.,1e7); + float max = 1.1*gYieldCorr->GetMaximum(); + histo_draw->SetAxisRange(0.1,max,"Y"); + histo_draw->SetStats(0); + histo_draw->GetXaxis()->SetTitle("p_{T} [GeV]"); + histo_draw->GetXaxis()->SetTitleSize(0.05); + histo_draw->GetXaxis()->SetTitleOffset(0.95); + histo_draw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]"); + histo_draw->GetYaxis()->SetTitleSize(0.05); + TCanvas * cyield_asym = new TCanvas("cyield_asym","Asymmetric corrected yield (TGraphAsymmErr)"); + gYieldCorr->SetFillStyle(3001); + gYieldCorr->SetLineWidth(3); + gYieldCorr->SetMarkerStyle(20); + gYieldCorr->SetFillColor(3); + histo_draw->Draw(); + gYieldCorr->Draw("3Csame"); + gYieldCorr->Draw("Xsame"); + cyield_asym->SetLogy(); + cyield_asym->Update(); + + TH2F *histo2_draw = new TH2F("histo2_draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9); + max = 1.1*gSigmaCorr->GetMaximum(); + histo2_draw->SetAxisRange(0.1,max,"Y"); + histo2_draw->SetStats(0); + histo2_draw->GetXaxis()->SetTitle("p_{T} [GeV]"); + histo2_draw->GetXaxis()->SetTitleSize(0.05); + histo2_draw->GetXaxis()->SetTitleOffset(0.95); + histo2_draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}"); + histo2_draw->GetYaxis()->SetTitleSize(0.05); + TCanvas * csigma_asym = new TCanvas("csigma_asym","Asymmetric corrected sigma (TGraphAsymmErr)"); + gSigmaCorr->SetFillStyle(3001); + gSigmaCorr->SetLineWidth(3); + gSigmaCorr->SetMarkerStyle(21); + gSigmaCorr->SetFillColor(4); + histo2_draw->Draw(); + gSigmaCorr->Draw("3Csame"); + gSigmaCorr->Draw("Xsame"); + csigma_asym->SetLogy(); + csigma_asym->Update(); + } + + + + // + // Write the histograms to the output file + // + out->cd(); + // + hDirectMCpt->Write(); hFeedDownMCpt->Write(); + hDirectMCpt_max->Write(); hDirectMCpt_min->Write(); + hFeedDownMCpt_max->Write(); hFeedDownMCpt_min->Write(); + hDirectEffpt->Write(); hFeedDownEffpt->Write(); + hRECpt->Write(); + // + histoYieldCorr->Write(); + histoYieldCorr_max->Write(); histoYieldCorr_min->Write(); + histoSigmaCorr->Write(); + + if(asym){ + gYieldCorr->Write(); + gSigmaCorr->Write(); + } + + if(option==1){ + histofc->Write(); + histofc_max->Write(); histofc_min->Write(); + if(asym) gFc->Write(); + } + + // out->Close(); + +} -- 2.39.3