New class to produce corrected dsigma/dpt for charm hadrons, with efficiency and...
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 17 Sep 2010 23:04:14 +0000 (23:04 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 17 Sep 2010 23:04:14 +0000 (23:04 +0000)
PWG3/PWG3vertexingHFLinkDef.h
PWG3/libPWG3vertexingHF.pkg
PWG3/vertexingHF/AliHFPtSpectrum.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliHFPtSpectrum.h [new file with mode: 0644]
PWG3/vertexingHF/macros/HFPtSpectrum.C [new file with mode: 0644]

index a9743a4..9fa7d16 100644 (file)
@@ -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+;
index 0ef320f..577e52a 100644 (file)
@@ -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 (file)
index 0000000..7cb3b2d
--- /dev/null
@@ -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 <Riostream.h>
+
+#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 (file)
index 0000000..e8c0f0c
--- /dev/null
@@ -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 (file)
index 0000000..b8f309f
--- /dev/null
@@ -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 <Riostream.h>
+#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"<<endl;
+    return;
+  }
+
+  if (option>2) { 
+    cout<< "Bad calculation option, should be <=2"<<endl;
+    return;
+  }
+  if (option==0) asym = false;
+
+  //
+  // Get the histograms from the files
+  //
+  TH1D *hDirectMCpt;            // Input MC c-->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();
+
+}