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