PWGHF/vertexingHF/AliHFPtSpectrum.h

   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   }
  81   void SetIsEventPlaneAnalysis(Bool_t flag){ fIsEventPlane = flag; }
  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
 297   Bool_t fIsEventPlane;             // flag : when the analysis is done for In/Out of plane, divide the B-prediction by two
 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
 305   ClassDef(AliHFPtSpectrum,4) // Class for Heavy Flavor spectra corrections
 306 };
 307
 308 #endif