[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliHFPtSpectrum.h

#ifndef ALIHFPTSPECTRUM_H
#define ALIHFPTSPECTRUM_H

/* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
 * See cxx source for full Copyright notice                               */

/* $Id$ */ 

//***********************************************************************
// 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)
//
//  (the corrected yields per bin are divided by the bin-width)
//
//
//  In HIC you can also evaluate how the feed-down correction is influenced by an energy loss hypothesis: 
//      Raa(c-->D) / Raa(b-->D) defined here as Rcb for the "fc" method
//      Raa(b-->D) defined here as Rb for the "Nb" method
//
// Author: Z.Conesa, zconesa@in2p3.fr
//***********************************************************************

#include "TNamed.h"
#include "TMath.h"

#include "AliLog.h"

class TH1;
class TH2;
class TNtuple;
class TGraphAsymmErrors;


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(TH1D *hDirect, TH1D *hFeedDown);
  // Set the theoretical feeddown pt spectrum
  void SetFeedDownMCptSpectra(TH1D *hFeedDown);
  // Set the theoretical direct & feeddown pt spectrum upper and lower bounds
  void SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin);
  // Set the theoretical feeddown pt spectrum upper and lower bounds
  void SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin);
  // Set the acceptance and efficiency corrections for direct  
  void SetDirectAccEffCorrection(TH1D *hDirectEff);
  // Set the acceptance and efficiency corrections for direct & feeddown 
  void SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff);
  // Set the reconstructed spectrum
  void SetReconstructedSpectrum(TH1D *hRec);
  void SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec); 
  // 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_t flag){ fAsymUncertainties = flag; }
  // Set if the yield is for particle plus anti-particle or not
  void SetIsParticlePlusAntiParticleYield(Bool_t flag){
    if (flag) { fParticleAntiParticle = 2; AliInfo(" Setting for particle + anti-particle yields"); }
    else { fParticleAntiParticle = 1; AliInfo(" Setting for only (anti)particle yields, not the sum of both"); }
  }
  void SetIsEventPlaneAnalysis(Bool_t flag){ fIsEventPlane = flag; }
  //
  void SetfIsStatUncEff(Bool_t flag){ fIsStatUncEff = flag; }
  // Set if the calculation has to consider Ratio(c/b eloss) hypothesis 
  void SetComputeElossHypothesis(Bool_t flag){ fPbPbElossHypothesis = flag; }
  // Set the luminosity and its uncertainty
  void SetLuminosity(Double_t luminosity, Double_t unc){
    fLuminosity[0]=luminosity;  fLuminosity[1]=unc;
  }
  // Set the trigger efficiency and its uncertainty
  void SetTriggerEfficiency(Double_t efficiency, Double_t unc){
    fTrigEfficiency[0]=efficiency; fTrigEfficiency[1]=unc;
  }
  // Set global acceptance x efficiency correction uncertainty (in percentages)
  void SetAccEffPercentageUncertainty(Double_t globalEffUnc, Double_t globalBCEffRatioUnc){
    fGlobalEfficiencyUncertainties[0] = globalEffUnc;
    fGlobalEfficiencyUncertainties[1] = globalBCEffRatioUnc;
  }
  // Set the normalization factors
  void SetNormalization(Double_t normalization){
    fLuminosity[0]=normalization;
  }
  void SetNormalization(Int_t nevents, Double_t sigma){
    fLuminosity[0]=nevents/sigma;
    fNevts = nevents;
  }
  void SetNormalization(Int_t nevents, Double_t sigma, Double_t sigmaunc){
    fLuminosity[0] = nevents/sigma; 
    fLuminosity[1] = fLuminosity[0] * TMath::Sqrt( (1/nevents) + (sigmaunc/sigma)*(sigmaunc/sigma) );
    fNevts = nevents;
  }
  //
  // Set the Tab parameter and its uncertainty
  void SetTabParameter(Double_t tabvalue, Double_t uncertainty){
    fTab[0] = tabvalue;
    fTab[1] = uncertainty;
  }


  //
  // Getters
  //
  // Return the theoretical predictions used for the calculation (rebinned if needed)
  TH1D * GetDirectTheoreticalSpectrum() const { return (fhDirectMCpt ? (TH1D*)fhDirectMCpt : NULL); }
  TH1D * GetDirectTheoreticalUpperLimitSpectrum() const { return (fhDirectMCptMax ? (TH1D*)fhDirectMCptMax : NULL); }
  TH1D * GetDirectTheoreticalLowerLimitSpectrum() const { return (fhDirectMCptMin ? (TH1D*)fhDirectMCptMin : NULL); }
  TH1D * GetFeedDownTheoreticalSpectrum() const { return (fhFeedDownMCpt ? (TH1D*)fhFeedDownMCpt : NULL); }
  TH1D * GetFeedDownTheoreticalUpperLimitSpectrum() const { return (fhFeedDownMCptMax ? (TH1D*)fhFeedDownMCptMax : NULL); }
  TH1D * GetFeedDownTheoreticalLowerLimitSpectrum() const { return (fhFeedDownMCptMin ? (TH1D*)fhFeedDownMCptMin : NULL); }
  // Return the acceptance and efficiency corrections (rebinned if needed)
  TH1D * GetDirectAccEffCorrection() const { return (fhDirectEffpt ? (TH1D*)fhDirectEffpt : NULL); }
  TH1D * GetFeedDownAccEffCorrection() const { return (fhFeedDownEffpt ? (TH1D*)fhFeedDownEffpt : NULL); }
  // Return whether the Ratio(c/b eloss) hypothesis has been considered
  Bool_t IsElossHypothesisCalculated(){ return fPbPbElossHypothesis; }
  // Return the TGraphAsymmErrors of the feed-down correction (extreme systematics)
  TGraphAsymmErrors * GetFeedDownCorrectionFcExtreme() const { return (fgFcExtreme ?  fgFcExtreme : NULL); }
  // Return the TGraphAsymmErrors of the feed-down correction (conservative systematics)
  TGraphAsymmErrors * GetFeedDownCorrectionFcConservative() const { return (fgFcConservative ?  fgFcConservative : NULL); }
  // Return the histogram of the feed-down correction
  TH1D * GetHistoFeedDownCorrectionFc() const { return (fhFc ?  (TH1D*)fhFc : NULL); }
  // Return the histograms of the feed-down correction bounds
  TH1D * GetHistoUpperLimitFeedDownCorrectionFc() const { return (fhFcMax ? (TH1D*)fhFcMax : NULL); }
  TH1D * GetHistoLowerLimitFeedDownCorrectionFc() const { return (fhFcMin ? (TH1D*)fhFcMin : NULL); }
  // Return the histogram of the feed-down correction times the Ratio(c/b eloss)
  TH2D * GetHistoFeedDownCorrectionFcVsEloss() const { return (fhFcRcb ?  (TH2D*)fhFcRcb : NULL); }
  // Return the TGraphAsymmErrors of the yield after feed-down correction (systematics but feed-down) 
  TGraphAsymmErrors * GetFeedDownCorrectedSpectrum() const { return (fgYieldCorr ? fgYieldCorr : NULL); }
  // Return the TGraphAsymmErrors of the yield after feed-down correction (feed-down extreme systematics)
  TGraphAsymmErrors * GetFeedDownCorrectedSpectrumExtreme() const { return (fgYieldCorrExtreme ? fgYieldCorrExtreme : NULL); }
  // Return the TGraphAsymmErrors of the yield after feed-down correction (feed-down conservative systematics)
  TGraphAsymmErrors * GetFeedDownCorrectedSpectrumConservative() const { return (fgYieldCorrConservative ? fgYieldCorrConservative : NULL); }
  // Return the histogram of the yield after feed-down correction 
  TH1D * GetHistoFeedDownCorrectedSpectrum() const { return (fhYieldCorr ? (TH1D*)fhYieldCorr : NULL); }
  // Return the histogram of the yield after feed-down correction bounds
  TH1D * GetHistoUpperLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMax ? (TH1D*)fhYieldCorrMax : NULL); }
  TH1D * GetHistoLowerLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMin ? (TH1D*)fhYieldCorrMin : NULL); }
  // Return the histogram of the yield after feed-down correction vs the Ratio(c/b eloss)
  TH2D * GetHistoFeedDownCorrectedSpectrumVsEloss() const { return (fhYieldCorrRcb ? (TH2D*)fhYieldCorrRcb : NULL); }
  // Return the equivalent invariant cross-section TGraphAsymmErrors (systematics but feed-down) 
  TGraphAsymmErrors * GetCrossSectionFromYieldSpectrum() const { return (fgSigmaCorr ? fgSigmaCorr : NULL); }
  // Return the equivalent invariant cross-section TGraphAsymmErrors (feed-down extreme systematics)
  TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumExtreme() const { return (fgSigmaCorrExtreme ? fgSigmaCorrExtreme : NULL); }
  // Return the equivalent invariant cross-section TGraphAsymmErrors (feed-down conservative systematics)
  TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumConservative() const { return (fgSigmaCorrConservative ? fgSigmaCorrConservative : NULL); }
  // Return the equivalent invariant cross-section histogram
  TH1D * GetHistoCrossSectionFromYieldSpectrum() const { return (fhSigmaCorr ? (TH1D*)fhSigmaCorr : NULL); }
  // Return the equivalent invariant cross-section histogram bounds
  TH1D * GetHistoUpperLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMax ? (TH1D*)fhSigmaCorrMax : NULL); }
  TH1D * GetHistoLowerLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMin ? (TH1D*)fhSigmaCorrMin : NULL); }
  // Return the cross section systematics from data systematics
  TH1D * GetHistoCrossSectionDataSystematics() const { return (fhSigmaCorrDataSyst ? (TH1D*)fhSigmaCorrDataSyst : NULL); }
  //
  // PbPb special calculations 
  // Return the equivalent invariant cross-section histogram vs the Ratio(c/b eloss)
  TH2D * GetHistoCrossSectionFromYieldSpectrumVsEloss() const { return (fhSigmaCorrRcb ? (TH2D*)fhSigmaCorrRcb : NULL); }
  // Return the ntuple of the calculation vs the Ratio(c/b eloss)
  TNtuple * GetNtupleCrossSectionVsEloss() { return (fnSigma ? (TNtuple*)fnSigma : NULL); }
  //
  //
  // Histograms to keep track of the influence of the efficiencies statistical uncertainty on the cross-section
  TH1D * GetDirectStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffcSigma; }
  TH1D * GetFeedDownStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffbSigma; }
  // Histograms to keep track of the influence of the efficiencies statistical uncertainty on the feed-down correction factor
  TH1D * GetDirectStatEffUncOnFc() const { return (TH1D*)fhStatUncEffcFD; }
  TH1D * GetFeedDownStatEffUncOnFc() const { return (TH1D*)fhStatUncEffbFD; }


  //
  // Main function:
  //    Compute the invariant cross-section from the yield (correct it)
  // variables : analysed delta_y, BR for the final correction, BR b --> decay (relative to the input theoretical prediction)
  void ComputeHFPtSpectrum(Double_t deltaY=1.0, Double_t branchingRatioC=1.0, Double_t branchingRatioBintoFinalDecay=1.0);

  // Compute the systematic uncertainties
  //   taking as input the AliHFSystErr uncertainties
  void ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown);
  //
  // Drawing the corrected spectrum comparing to theoretical prediction
  void DrawSpectrum(TGraphAsymmErrors *gPrediction);

  //
  // Basic functions
  // 
  void EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco);
  void EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco);
  
  //
  // Functions to  reweight histograms for testing purposes: 
  //   to reweight the simulation: hToReweight is reweighted as hReference/hToReweight
  TH1D * ReweightHisto(TH1D *hToReweight, TH1D *hReference);
  //   to reweight the reco-histos: hRecToReweight is reweighted as hReference/hMCToReweight
  TH1D * ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference);
  // Functionality to find the y-axis bin of a TH2 for a given y-value
  Int_t FindTH2YBin(TH2D *histo, Float_t yvalue);


 protected:

  // Initialization 
  Bool_t Initialize();
  
  // Basic functions
  //
  // Compute the feed-down correction via fc-method
  void CalculateFeedDownCorrectionFc(); 
  // Correct the yield for feed-down correction via fc-method
  void CalculateFeedDownCorrectedSpectrumFc(); 
  // Correct the yield for feed-down correction via Nb-method
  void CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay); 

  // Check histograms consistency function
  Bool_t CheckHistosConsistency(TH1D *h1, TH1D *h2);
  // Function to rebin the theoretical spectra in the data-reconstructed spectra binning
  TH1D * RebinTheoreticalSpectra(TH1D *hTheory, const char *name);
  // Function to estimate the efficiency in the data-reconstructed spectra binning
  TH1D * EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name);
  // Reset stat unc on the efficiencies
  void ResetStatUncEff();


  //
  // Input spectra
  //
  TH1D *fhDirectMCpt;            // Input MC c-->D spectra
  TH1D *fhFeedDownMCpt;          // Input MC b-->D spectra
  TH1D *fhDirectMCptMax;         // Input MC maximum c-->D spectra
  TH1D *fhDirectMCptMin;         // Input MC minimum c-->D spectra
  TH1D *fhFeedDownMCptMax;       // Input MC maximum b-->D spectra
  TH1D *fhFeedDownMCptMin;       // Input MC minimum b-->D spectra
  TH1D *fhDirectEffpt;           // c-->D Acceptance and efficiency correction
  TH1D *fhFeedDownEffpt;         // b-->D Acceptance and efficiency correction
  TH1D *fhRECpt;                 // all reconstructed D
  //
  TGraphAsymmErrors *fgRECSystematics; // all reconstructed D Systematic uncertainties
  //
  // Normalization factors
  Int_t fNevts;                      // nb of analyzed events
  Double_t fLuminosity[2];           // analyzed luminosity & uncertainty
  Double_t fTrigEfficiency[2];       // trigger efficiency & uncertainty
  Double_t fGlobalEfficiencyUncertainties[2]; // uncertainties on the efficiency [0]=c, b, [1]=b/c
  Double_t fTab[2];                   // Tab parameter and its uncertainty

  //
  // Output spectra
  //
  TH1D *fhFc;                            // Correction histo fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) 
  TH1D *fhFcMax;                         // Maximum fc histo
  TH1D *fhFcMin;                         // Minimum fc histo
  TH2D *fhFcRcb;                         // Correction histo fc vs the Ratio(c/b eloss)
  TGraphAsymmErrors * fgFcExtreme;       // Extreme correction as TGraphAsymmErrors
  TGraphAsymmErrors * fgFcConservative;  // Extreme correction as TGraphAsymmErrors
  TH1D *fhYieldCorr;                     // Corrected yield (stat unc. only)
  TH1D *fhYieldCorrMax;                  // Maximum corrected yield  
  TH1D *fhYieldCorrMin;                  // Minimum corrected yield  
  TH2D *fhYieldCorrRcb;                  // Corrected yield (stat unc. only) vs the Ratio(c/b eloss)
  TGraphAsymmErrors * fgYieldCorr;              // Corrected yield as TGraphAsymmErrors  (syst but feed-down)
  TGraphAsymmErrors * fgYieldCorrExtreme;       // Extreme corrected yield as TGraphAsymmErrors  (syst from feed-down)
  TGraphAsymmErrors * fgYieldCorrConservative;  // Conservative corrected yield as TGraphAsymmErrors  (syst from feed-down) 
  TH1D *fhSigmaCorr;                     // Corrected cross-section (stat unc. only)
  TH1D *fhSigmaCorrMax;                  // Maximum corrected cross-section  
  TH1D *fhSigmaCorrMin;                  // Minimum corrected cross-section
  TH1D *fhSigmaCorrDataSyst;             // Corrected cross-section (syst. unc. from data only)
  TH2D *fhSigmaCorrRcb;                  // Corrected cross-section (stat unc. only) vs the Ratio(c/b eloss)
  TGraphAsymmErrors * fgSigmaCorr;              // Corrected cross-section as TGraphAsymmErrors (syst but feed-down)
  TGraphAsymmErrors * fgSigmaCorrExtreme;       // Extreme corrected cross-section as TGraphAsymmErrors (syst from feed-down)
  TGraphAsymmErrors * fgSigmaCorrConservative;  // Conservative corrected cross-section as TGraphAsymmErrors  (syst from feed-down)
  //
  TNtuple *fnSigma;          // Ntuple of the calculation vs the Ratio(c/b eloss)
  TNtuple *fnHypothesis;     // Ntuple of the calculation vs the Ratio(c/b eloss)

  //
  Int_t fFeedDownOption;            // feed-down correction flag: 0=none, 1=fc, 2=Nb 
  Bool_t fAsymUncertainties;        // flag: asymmetric uncertainties are (1) or not (0) considered
  Bool_t fPbPbElossHypothesis;      // flag: whether to do estimates vs Ratio(c/b eloss) hypothesis
  Bool_t fIsStatUncEff;             // flag : consider (1) or not (0) the stat unc on the efficiencies
  Int_t fParticleAntiParticle;      // 1: only one sign, 2: yield is for particle+anti-particle
  Bool_t fIsEventPlane;             // flag : when the analysis is done for In/Out of plane, divide the B-prediction by two

  //
  TH1D *fhStatUncEffcSigma;       // Uncertainty on the cross-section due to the prompt efficiency statistical uncertainty
  TH1D *fhStatUncEffbSigma;       // Uncertainty on the cross-section due to the feed-down efficiency statistical uncertainty
  TH1D *fhStatUncEffcFD;          // Uncertainty on the feed-down correction due to the prompt efficiency statistical uncertainty
  TH1D *fhStatUncEffbFD;          // Uncertainty on the feed-down correction due to the feed-down efficiency statistical uncertainty

  ClassDef(AliHFPtSpectrum,4) // Class for Heavy Flavor spectra corrections
};

#endif
Commit	Line	Data
b188dc47	1	#ifndef ALIHFPTSPECTRUM_H
	2	#define ALIHFPTSPECTRUM_H
	3
	4	/* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
	5	* See cxx source for full Copyright notice */
	6
	7	/* $Id$ */
	8
	9	//***********************************************************************
	10	// Class AliHFPtSpectrum
	11	// Base class for feed-down corrections on heavy-flavour decays
	12	// computes the cross-section via one of the three implemented methods:
	13	// 0) Consider no feed-down prediction
	14	// 1) Subtract the feed-down with the "fc" method
	15	// Yield = Reco * fc; where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ;
	16	// 2) Subtract the feed-down with the "Nb" method
	17	// Yield = Reco - Feed-down (exact formula on the function implementation)
	18	//
	19	// (the corrected yields per bin are divided by the bin-width)
	20	//
	21	//
	22	// In HIC you can also evaluate how the feed-down correction is influenced by an energy loss hypothesis:
	23	// Raa(c-->D) / Raa(b-->D) defined here as Rcb for the "fc" method
	24	// Raa(b-->D) defined here as Rb for the "Nb" method
	25	//
	26	// Author: Z.Conesa, zconesa@in2p3.fr
	27	//***********************************************************************
	28
	29	#include "TNamed.h"
	30	#include "TMath.h"
	31
	32	#include "AliLog.h"
	33
	34	class TH1;
	35	class TH2;
	36	class TNtuple;
	37	class TGraphAsymmErrors;
	38
	39
	40	class AliHFPtSpectrum: public TNamed
	41	{
	42
	43	public:
	44
	45	// Constructor
	46	AliHFPtSpectrum(const char* name="AliHFPtSpectrum", const char* title="HF feed down correction class", Int_t option=1);
	47	// Copy constructor
	48	AliHFPtSpectrum(const AliHFPtSpectrum &rhs);
	49	// Assignment operator
	50	AliHFPtSpectrum& operator=(const AliHFPtSpectrum &source);
	51	// Destructor
	52	virtual ~AliHFPtSpectrum();
	53
	54	//
	55	// Setters
	56	//
	57	// Set the theoretical direct & feeddown pt spectrum
	58	void SetMCptSpectra(TH1D hDirect, TH1D hFeedDown);
	59	// Set the theoretical feeddown pt spectrum
	60	void SetFeedDownMCptSpectra(TH1D *hFeedDown);
	61	// Set the theoretical direct & feeddown pt spectrum upper and lower bounds
	62	void SetMCptDistributionsBounds(TH1D hDirectMax, TH1D hDirectMin, TH1D hFeedDownMax, TH1D hFeedDownMin);
	63	// Set the theoretical feeddown pt spectrum upper and lower bounds
	64	void SetFeedDownMCptDistributionsBounds(TH1D hFeedDownMax, TH1D hFeedDownMin);
65	// Set the acceptance and efficiency corrections for direct
66	void SetDirectAccEffCorrection(TH1D *hDirectEff);
67	// Set the acceptance and efficiency corrections for direct & feeddown
68	void SetAccEffCorrection(TH1D hDirectEff, TH1D hFeedDownEff);
69	// Set the reconstructed spectrum
70	void SetReconstructedSpectrum(TH1D *hRec);
71	void SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec);
72	// Set the calculation option flag for feed-down correction: 0=none, 1=fc , 2=Nb
73	void SetFeedDownCalculationOption(Int_t option){ fFeedDownOption = option; }
74	// Set if the calculation has to consider asymmetric uncertaInt_ties or not
75	void SetComputeAsymmetricUncertainties(Bool_t flag){ fAsymUncertainties = flag; }
76	// Set if the yield is for particle plus anti-particle or not
77	void SetIsParticlePlusAntiParticleYield(Bool_t flag){
78	if (flag) { fParticleAntiParticle = 2; AliInfo(" Setting for particle + anti-particle yields"); }
79	else { fParticleAntiParticle = 1; AliInfo(" Setting for only (anti)particle yields, not the sum of both"); }
80	}
ddd86f95	81	void SetIsEventPlaneAnalysis(Bool_t flag){ fIsEventPlane = flag; }
b188dc47	82	//
	83	void SetfIsStatUncEff(Bool_t flag){ fIsStatUncEff = flag; }
	84	// Set if the calculation has to consider Ratio(c/b eloss) hypothesis
	85	void SetComputeElossHypothesis(Bool_t flag){ fPbPbElossHypothesis = flag; }
	86	// Set the luminosity and its uncertainty
	87	void SetLuminosity(Double_t luminosity, Double_t unc){
	88	fLuminosity[0]=luminosity; fLuminosity[1]=unc;
	89	}
	90	// Set the trigger efficiency and its uncertainty
	91	void SetTriggerEfficiency(Double_t efficiency, Double_t unc){
	92	fTrigEfficiency[0]=efficiency; fTrigEfficiency[1]=unc;
	93	}
	94	// Set global acceptance x efficiency correction uncertainty (in percentages)
	95	void SetAccEffPercentageUncertainty(Double_t globalEffUnc, Double_t globalBCEffRatioUnc){
	96	fGlobalEfficiencyUncertainties[0] = globalEffUnc;
	97	fGlobalEfficiencyUncertainties[1] = globalBCEffRatioUnc;
	98	}
	99	// Set the normalization factors
	100	void SetNormalization(Double_t normalization){
	101	fLuminosity[0]=normalization;
	102	}
	103	void SetNormalization(Int_t nevents, Double_t sigma){
	104	fLuminosity[0]=nevents/sigma;
	105	fNevts = nevents;
	106	}
	107	void SetNormalization(Int_t nevents, Double_t sigma, Double_t sigmaunc){
	108	fLuminosity[0] = nevents/sigma;
	109	fLuminosity[1] = fLuminosity[0] * TMath::Sqrt( (1/nevents) + (sigmaunc/sigma)*(sigmaunc/sigma) );
	110	fNevts = nevents;
	111	}
	112	//
	113	// Set the Tab parameter and its uncertainty
	114	void SetTabParameter(Double_t tabvalue, Double_t uncertainty){
	115	fTab[0] = tabvalue;
	116	fTab[1] = uncertainty;
	117	}
	118
	119
	120	//
	121	// Getters
	122	//
	123	// Return the theoretical predictions used for the calculation (rebinned if needed)
	124	TH1D * GetDirectTheoreticalSpectrum() const { return (fhDirectMCpt ? (TH1D*)fhDirectMCpt : NULL); }
	125	TH1D * GetDirectTheoreticalUpperLimitSpectrum() const { return (fhDirectMCptMax ? (TH1D*)fhDirectMCptMax : NULL); }
	126	TH1D * GetDirectTheoreticalLowerLimitSpectrum() const { return (fhDirectMCptMin ? (TH1D*)fhDirectMCptMin : NULL); }
	127	TH1D * GetFeedDownTheoreticalSpectrum() const { return (fhFeedDownMCpt ? (TH1D*)fhFeedDownMCpt : NULL); }
	128	TH1D * GetFeedDownTheoreticalUpperLimitSpectrum() const { return (fhFeedDownMCptMax ? (TH1D*)fhFeedDownMCptMax : NULL); }
	129	TH1D * GetFeedDownTheoreticalLowerLimitSpectrum() const { return (fhFeedDownMCptMin ? (TH1D*)fhFeedDownMCptMin : NULL); }
	130	// Return the acceptance and efficiency corrections (rebinned if needed)
	131	TH1D * GetDirectAccEffCorrection() const { return (fhDirectEffpt ? (TH1D*)fhDirectEffpt : NULL); }
	132	TH1D * GetFeedDownAccEffCorrection() const { return (fhFeedDownEffpt ? (TH1D*)fhFeedDownEffpt : NULL); }
	133	// Return whether the Ratio(c/b eloss) hypothesis has been considered
	134	Bool_t IsElossHypothesisCalculated(){ return fPbPbElossHypothesis; }
	135	// Return the TGraphAsymmErrors of the feed-down correction (extreme systematics)
	136	TGraphAsymmErrors * GetFeedDownCorrectionFcExtreme() const { return (fgFcExtreme ? fgFcExtreme : NULL); }
	137	// Return the TGraphAsymmErrors of the feed-down correction (conservative systematics)
	138	TGraphAsymmErrors * GetFeedDownCorrectionFcConservative() const { return (fgFcConservative ? fgFcConservative : NULL); }
	139	// Return the histogram of the feed-down correction
	140	TH1D * GetHistoFeedDownCorrectionFc() const { return (fhFc ? (TH1D*)fhFc : NULL); }
	141	// Return the histograms of the feed-down correction bounds
	142	TH1D * GetHistoUpperLimitFeedDownCorrectionFc() const { return (fhFcMax ? (TH1D*)fhFcMax : NULL); }
	143	TH1D * GetHistoLowerLimitFeedDownCorrectionFc() const { return (fhFcMin ? (TH1D*)fhFcMin : NULL); }
	144	// Return the histogram of the feed-down correction times the Ratio(c/b eloss)
	145	TH2D * GetHistoFeedDownCorrectionFcVsEloss() const { return (fhFcRcb ? (TH2D*)fhFcRcb : NULL); }
146	// Return the TGraphAsymmErrors of the yield after feed-down correction (systematics but feed-down)
147	TGraphAsymmErrors * GetFeedDownCorrectedSpectrum() const { return (fgYieldCorr ? fgYieldCorr : NULL); }
148	// Return the TGraphAsymmErrors of the yield after feed-down correction (feed-down extreme systematics)
149	TGraphAsymmErrors * GetFeedDownCorrectedSpectrumExtreme() const { return (fgYieldCorrExtreme ? fgYieldCorrExtreme : NULL); }
150	// Return the TGraphAsymmErrors of the yield after feed-down correction (feed-down conservative systematics)
151	TGraphAsymmErrors * GetFeedDownCorrectedSpectrumConservative() const { return (fgYieldCorrConservative ? fgYieldCorrConservative : NULL); }
152	// Return the histogram of the yield after feed-down correction
153	TH1D * GetHistoFeedDownCorrectedSpectrum() const { return (fhYieldCorr ? (TH1D*)fhYieldCorr : NULL); }
154	// Return the histogram of the yield after feed-down correction bounds
155	TH1D * GetHistoUpperLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMax ? (TH1D*)fhYieldCorrMax : NULL); }
156	TH1D * GetHistoLowerLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMin ? (TH1D*)fhYieldCorrMin : NULL); }
157	// Return the histogram of the yield after feed-down correction vs the Ratio(c/b eloss)
158	TH2D * GetHistoFeedDownCorrectedSpectrumVsEloss() const { return (fhYieldCorrRcb ? (TH2D*)fhYieldCorrRcb : NULL); }
159	// Return the equivalent invariant cross-section TGraphAsymmErrors (systematics but feed-down)
160	TGraphAsymmErrors * GetCrossSectionFromYieldSpectrum() const { return (fgSigmaCorr ? fgSigmaCorr : NULL); }
161	// Return the equivalent invariant cross-section TGraphAsymmErrors (feed-down extreme systematics)
162	TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumExtreme() const { return (fgSigmaCorrExtreme ? fgSigmaCorrExtreme : NULL); }
163	// Return the equivalent invariant cross-section TGraphAsymmErrors (feed-down conservative systematics)
164	TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumConservative() const { return (fgSigmaCorrConservative ? fgSigmaCorrConservative : NULL); }
165	// Return the equivalent invariant cross-section histogram
166	TH1D * GetHistoCrossSectionFromYieldSpectrum() const { return (fhSigmaCorr ? (TH1D*)fhSigmaCorr : NULL); }
167	// Return the equivalent invariant cross-section histogram bounds
168	TH1D * GetHistoUpperLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMax ? (TH1D*)fhSigmaCorrMax : NULL); }
169	TH1D * GetHistoLowerLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMin ? (TH1D*)fhSigmaCorrMin : NULL); }
170	// Return the cross section systematics from data systematics
171	TH1D * GetHistoCrossSectionDataSystematics() const { return (fhSigmaCorrDataSyst ? (TH1D*)fhSigmaCorrDataSyst : NULL); }
172	//
173	// PbPb special calculations
174	// Return the equivalent invariant cross-section histogram vs the Ratio(c/b eloss)
175	TH2D * GetHistoCrossSectionFromYieldSpectrumVsEloss() const { return (fhSigmaCorrRcb ? (TH2D*)fhSigmaCorrRcb : NULL); }
176	// Return the ntuple of the calculation vs the Ratio(c/b eloss)
177	TNtuple * GetNtupleCrossSectionVsEloss() { return (fnSigma ? (TNtuple*)fnSigma : NULL); }
178	//
179	//
180	// Histograms to keep track of the influence of the efficiencies statistical uncertainty on the cross-section
181	TH1D * GetDirectStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffcSigma; }
182	TH1D * GetFeedDownStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffbSigma; }
183	// Histograms to keep track of the influence of the efficiencies statistical uncertainty on the feed-down correction factor
184	TH1D * GetDirectStatEffUncOnFc() const { return (TH1D*)fhStatUncEffcFD; }
185	TH1D * GetFeedDownStatEffUncOnFc() const { return (TH1D*)fhStatUncEffbFD; }
186
187
188	//
189	// Main function:
190	// Compute the invariant cross-section from the yield (correct it)
191	// variables : analysed delta_y, BR for the final correction, BR b --> decay (relative to the input theoretical prediction)
192	void ComputeHFPtSpectrum(Double_t deltaY=1.0, Double_t branchingRatioC=1.0, Double_t branchingRatioBintoFinalDecay=1.0);
193
194	// Compute the systematic uncertainties
195	// taking as input the AliHFSystErr uncertainties
196	void ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown);
197	//
198	// Drawing the corrected spectrum comparing to theoretical prediction
199	void DrawSpectrum(TGraphAsymmErrors *gPrediction);
200
201	//
202	// Basic functions
203	//
204	void EstimateAndSetDirectEfficiencyRecoBin(TH1D hSimu, TH1D hReco);
205	void EstimateAndSetFeedDownEfficiencyRecoBin(TH1D hSimu, TH1D hReco);
206
207	//
208	// Functions to reweight histograms for testing purposes:
209	// to reweight the simulation: hToReweight is reweighted as hReference/hToReweight
210	TH1D * ReweightHisto(TH1D hToReweight, TH1D hReference);
211	// to reweight the reco-histos: hRecToReweight is reweighted as hReference/hMCToReweight
212	TH1D * ReweightRecHisto(TH1D hRecToReweight, TH1D hMCToReweight, TH1D *hMCReference);
213	// Functionality to find the y-axis bin of a TH2 for a given y-value
214	Int_t FindTH2YBin(TH2D *histo, Float_t yvalue);
215
216
217	protected:
218
219	// Initialization
220	Bool_t Initialize();
221
222	// Basic functions
223	//
224	// Compute the feed-down correction via fc-method
225	void CalculateFeedDownCorrectionFc();
226	// Correct the yield for feed-down correction via fc-method
227	void CalculateFeedDownCorrectedSpectrumFc();
228	// Correct the yield for feed-down correction via Nb-method
229	void CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay);
230
231	// Check histograms consistency function
232	Bool_t CheckHistosConsistency(TH1D h1, TH1D h2);
233	// Function to rebin the theoretical spectra in the data-reconstructed spectra binning
234	TH1D * RebinTheoreticalSpectra(TH1D hTheory, const char name);
235	// Function to estimate the efficiency in the data-reconstructed spectra binning
236	TH1D * EstimateEfficiencyRecoBin(TH1D hSimu, TH1D hReco, const char *name);
237	// Reset stat unc on the efficiencies
238	void ResetStatUncEff();
239
240
241	//
242	// Input spectra
243	//
244	TH1D *fhDirectMCpt; // Input MC c-->D spectra
245	TH1D *fhFeedDownMCpt; // Input MC b-->D spectra
246	TH1D *fhDirectMCptMax; // Input MC maximum c-->D spectra
247	TH1D *fhDirectMCptMin; // Input MC minimum c-->D spectra
248	TH1D *fhFeedDownMCptMax; // Input MC maximum b-->D spectra
249	TH1D *fhFeedDownMCptMin; // Input MC minimum b-->D spectra
250	TH1D *fhDirectEffpt; // c-->D Acceptance and efficiency correction
251	TH1D *fhFeedDownEffpt; // b-->D Acceptance and efficiency correction
252	TH1D *fhRECpt; // all reconstructed D
253	//
254	TGraphAsymmErrors *fgRECSystematics; // all reconstructed D Systematic uncertainties
255	//
256	// Normalization factors
257	Int_t fNevts; // nb of analyzed events
258	Double_t fLuminosity[2]; // analyzed luminosity & uncertainty
259	Double_t fTrigEfficiency[2]; // trigger efficiency & uncertainty
260	Double_t fGlobalEfficiencyUncertainties[2]; // uncertainties on the efficiency [0]=c, b, [1]=b/c
261	Double_t fTab[2]; // Tab parameter and its uncertainty
262
263	//
264	// Output spectra
265	//
266	TH1D fhFc; // Correction histo fc = 1 / ( 1 + (eff_b/eff_c)(N_b/N_c) )
267	TH1D *fhFcMax; // Maximum fc histo
268	TH1D *fhFcMin; // Minimum fc histo
269	TH2D *fhFcRcb; // Correction histo fc vs the Ratio(c/b eloss)
270	TGraphAsymmErrors * fgFcExtreme; // Extreme correction as TGraphAsymmErrors
271	TGraphAsymmErrors * fgFcConservative; // Extreme correction as TGraphAsymmErrors
272	TH1D *fhYieldCorr; // Corrected yield (stat unc. only)
273	TH1D *fhYieldCorrMax; // Maximum corrected yield
274	TH1D *fhYieldCorrMin; // Minimum corrected yield
275	TH2D *fhYieldCorrRcb; // Corrected yield (stat unc. only) vs the Ratio(c/b eloss)
276	TGraphAsymmErrors * fgYieldCorr; // Corrected yield as TGraphAsymmErrors (syst but feed-down)
277	TGraphAsymmErrors * fgYieldCorrExtreme; // Extreme corrected yield as TGraphAsymmErrors (syst from feed-down)
278	TGraphAsymmErrors * fgYieldCorrConservative; // Conservative corrected yield as TGraphAsymmErrors (syst from feed-down)
279	TH1D *fhSigmaCorr; // Corrected cross-section (stat unc. only)
280	TH1D *fhSigmaCorrMax; // Maximum corrected cross-section
281	TH1D *fhSigmaCorrMin; // Minimum corrected cross-section
282	TH1D *fhSigmaCorrDataSyst; // Corrected cross-section (syst. unc. from data only)
283	TH2D *fhSigmaCorrRcb; // Corrected cross-section (stat unc. only) vs the Ratio(c/b eloss)
284	TGraphAsymmErrors * fgSigmaCorr; // Corrected cross-section as TGraphAsymmErrors (syst but feed-down)
285	TGraphAsymmErrors * fgSigmaCorrExtreme; // Extreme corrected cross-section as TGraphAsymmErrors (syst from feed-down)
286	TGraphAsymmErrors * fgSigmaCorrConservative; // Conservative corrected cross-section as TGraphAsymmErrors (syst from feed-down)
287	//
288	TNtuple *fnSigma; // Ntuple of the calculation vs the Ratio(c/b eloss)
289	TNtuple *fnHypothesis; // Ntuple of the calculation vs the Ratio(c/b eloss)
290
291	//
292	Int_t fFeedDownOption; // feed-down correction flag: 0=none, 1=fc, 2=Nb
293	Bool_t fAsymUncertainties; // flag: asymmetric uncertainties are (1) or not (0) considered
294	Bool_t fPbPbElossHypothesis; // flag: whether to do estimates vs Ratio(c/b eloss) hypothesis
295	Bool_t fIsStatUncEff; // flag : consider (1) or not (0) the stat unc on the efficiencies
296	Int_t fParticleAntiParticle; // 1: only one sign, 2: yield is for particle+anti-particle
ddd86f95	297	Bool_t fIsEventPlane; // flag : when the analysis is done for In/Out of plane, divide the B-prediction by two
b188dc47	298
	299	//
	300	TH1D *fhStatUncEffcSigma; // Uncertainty on the cross-section due to the prompt efficiency statistical uncertainty
	301	TH1D *fhStatUncEffbSigma; // Uncertainty on the cross-section due to the feed-down efficiency statistical uncertainty
	302	TH1D *fhStatUncEffcFD; // Uncertainty on the feed-down correction due to the prompt efficiency statistical uncertainty
	303	TH1D *fhStatUncEffbFD; // Uncertainty on the feed-down correction due to the feed-down efficiency statistical uncertainty
	304
ddd86f95	305	ClassDef(AliHFPtSpectrum,4) // Class for Heavy Flavor spectra corrections
b188dc47	306	};
	307
	308	#endif