AliHFPtSpectrum.h

Go to the documentation of this file.

#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$ */ 

//***********************************************************************
//***********************************************************************

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

#include "AliLog.h"

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


class AliHFPtSpectrum: public TNamed
{

 public:
  
  AliHFPtSpectrum(const char* name="AliHFPtSpectrum", const char* title="HF feed down correction class", Int_t option=1);
  AliHFPtSpectrum(const AliHFPtSpectrum &rhs);
  AliHFPtSpectrum& operator=(const AliHFPtSpectrum &source);
  virtual ~AliHFPtSpectrum();

  void SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown);
  void SetFeedDownMCptSpectra(TH1D *hFeedDown);
  void SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin);
  void SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin);
  void SetDirectAccEffCorrection(TH1D *hDirectEff);
  void SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff);
  void SetReconstructedSpectrum(TH1D *hRec);
  void SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec); 
  void SetFeedDownCalculationOption(Int_t option){ fFeedDownOption = option; }
  void SetComputeAsymmetricUncertainties(Bool_t flag){ fAsymUncertainties = flag; }
  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; }
  void SetComputeElossHypothesis(Bool_t flag){ fPbPbElossHypothesis = flag; }
  void SetLuminosity(Double_t luminosity, Double_t unc){
    fLuminosity[0]=luminosity;  fLuminosity[1]=unc;
  }
  void SetTriggerEfficiency(Double_t efficiency, Double_t unc){
    fTrigEfficiency[0]=efficiency; fTrigEfficiency[1]=unc;
  }
  void SetAccEffPercentageUncertainty(Double_t globalEffUnc, Double_t globalBCEffRatioUnc){
    fGlobalEfficiencyUncertainties[0] = globalEffUnc;
    fGlobalEfficiencyUncertainties[1] = globalBCEffRatioUnc;
  }
  void SetUsePtDependentEffUncertainty(Bool_t flag) {
    fGlobalEfficiencyPtDependent = flag;
  }
  void SetSystematicUncertainty(AliHFSystErr *syst){
    if(fSystematics) { delete fSystematics; fSystematics=NULL; }
    fSystematics = syst;
  }
  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;
  }
  //
  void SetTabParameter(Double_t tabvalue, Double_t uncertainty){
    fTab[0] = tabvalue;
    fTab[1] = uncertainty;
  }
  void SetCollisionType(Int_t ct){
    fCollisionType=ct;
  }

  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); }
  TH1D * GetDirectAccEffCorrection() const { return (fhDirectEffpt ? (TH1D*)fhDirectEffpt : NULL); }
  TH1D * GetFeedDownAccEffCorrection() const { return (fhFeedDownEffpt ? (TH1D*)fhFeedDownEffpt : NULL); }
  Bool_t IsElossHypothesisCalculated(){ return fPbPbElossHypothesis; }
  TGraphAsymmErrors * GetFeedDownCorrectionFcExtreme() const { return (fgFcExtreme ?  fgFcExtreme : NULL); }
  TGraphAsymmErrors * GetFeedDownCorrectionFcConservative() const { return (fgFcConservative ?  fgFcConservative : NULL); }
  TH1D * GetHistoFeedDownCorrectionFc() const { return (fhFc ?  (TH1D*)fhFc : NULL); }
  TH1D * GetHistoUpperLimitFeedDownCorrectionFc() const { return (fhFcMax ? (TH1D*)fhFcMax : NULL); }
  TH1D * GetHistoLowerLimitFeedDownCorrectionFc() const { return (fhFcMin ? (TH1D*)fhFcMin : NULL); }
  TH2D * GetHistoFeedDownCorrectionFcVsEloss() const { return (fhFcRcb ?  (TH2D*)fhFcRcb : NULL); }
  TGraphAsymmErrors * GetFeedDownCorrectedSpectrum() const { return (fgYieldCorr ? fgYieldCorr : NULL); }
  TGraphAsymmErrors * GetFeedDownCorrectedSpectrumExtreme() const { return (fgYieldCorrExtreme ? fgYieldCorrExtreme : NULL); }
  TGraphAsymmErrors * GetFeedDownCorrectedSpectrumConservative() const { return (fgYieldCorrConservative ? fgYieldCorrConservative : NULL); }
  TH1D * GetHistoFeedDownCorrectedSpectrum() const { return (fhYieldCorr ? (TH1D*)fhYieldCorr : NULL); }
  TH1D * GetHistoUpperLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMax ? (TH1D*)fhYieldCorrMax : NULL); }
  TH1D * GetHistoLowerLimitFeedDownCorrectedSpectrum() const { return (fhYieldCorrMin ? (TH1D*)fhYieldCorrMin : NULL); }
  TH2D * GetHistoFeedDownCorrectedSpectrumVsEloss() const { return (fhYieldCorrRcb ? (TH2D*)fhYieldCorrRcb : NULL); }
  TGraphAsymmErrors * GetCrossSectionFromYieldSpectrum() const { return (fgSigmaCorr ? fgSigmaCorr : NULL); }
  TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumExtreme() const { return (fgSigmaCorrExtreme ? fgSigmaCorrExtreme : NULL); }
  TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumConservative() const { return (fgSigmaCorrConservative ? fgSigmaCorrConservative : NULL); }
  TH1D * GetHistoCrossSectionFromYieldSpectrum() const { return (fhSigmaCorr ? (TH1D*)fhSigmaCorr : NULL); }
  TH1D * GetHistoUpperLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMax ? (TH1D*)fhSigmaCorrMax : NULL); }
  TH1D * GetHistoLowerLimitCrossSectionFromYieldSpectrum() const { return (fhSigmaCorrMin ? (TH1D*)fhSigmaCorrMin : NULL); }
  TH1D * GetHistoCrossSectionDataSystematics() const { return (fhSigmaCorrDataSyst ? (TH1D*)fhSigmaCorrDataSyst : NULL); }
  TH2D * GetHistoCrossSectionFromYieldSpectrumVsEloss() const { return (fhSigmaCorrRcb ? (TH2D*)fhSigmaCorrRcb : NULL); }
  TNtuple * GetNtupleCrossSectionVsEloss() { return (fnSigma ? (TNtuple*)fnSigma : NULL); }
  TH1D * GetDirectStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffcSigma; }
  TH1D * GetFeedDownStatEffUncOnSigma() const { return (TH1D*)fhStatUncEffbSigma; }
  TH1D * GetDirectStatEffUncOnFc() const { return (TH1D*)fhStatUncEffcFD; }
  TH1D * GetFeedDownStatEffUncOnFc() const { return (TH1D*)fhStatUncEffbFD; }


  void ComputeHFPtSpectrum(Double_t deltaY=1.0, Double_t branchingRatioC=1.0, Double_t branchingRatioBintoFinalDecay=1.0);

//    void ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown);
    void ComputeSystUncertainties(Bool_t combineFeedDown);
  void DrawSpectrum(TGraphAsymmErrors *gPrediction);

  void EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco);
  void EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco);
  
  TH1D * ReweightHisto(TH1D *hToReweight, TH1D *hReference);
  TH1D * ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference);
  Int_t FindTH2YBin(TH2D *histo, Float_t yvalue);


 protected:

  Bool_t Initialize();
  
  void CalculateCorrectedSpectrumNoFeedDown();
  void CalculateFeedDownCorrectionFc();
  void CalculateFeedDownCorrectedSpectrumFc();
  void CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay);
  Double_t CalculateEfficiencyPtDepedentUncertainty(Double_t pt, Bool_t useOnlyCutVar);

  Bool_t CheckHistosConsistency(TH1D *h1, TH1D *h2);
  TH1D * RebinTheoreticalSpectra(TH1D *hTheory, const char *name);
  TH1D * EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name);
  void ResetStatUncEff();


  TH1D *fhDirectMCpt;            
  TH1D *fhFeedDownMCpt;          
  TH1D *fhDirectMCptMax;         
  TH1D *fhDirectMCptMin;         
  TH1D *fhFeedDownMCptMax;       
  TH1D *fhFeedDownMCptMin;       
  TH1D *fhDirectEffpt;           
  TH1D *fhFeedDownEffpt;         
  TH1D *fhRECpt;                 
  //
  TGraphAsymmErrors *fgRECSystematics; 
  //
  Int_t fNevts;                      
  Double_t fLuminosity[2];           
  Double_t fTrigEfficiency[2];       
  Double_t fGlobalEfficiencyUncertainties[2]; 
  Bool_t fGlobalEfficiencyPtDependent;
  Double_t fTab[2];                   
  //
  AliHFSystErr *fSystematics;          

  //
  //
  TH1D *fhFc;                            
  TH1D *fhFcMax;                         
  TH1D *fhFcMin;                         
  TH2D *fhFcRcb;                         
  TGraphAsymmErrors * fgFcExtreme;       
  TGraphAsymmErrors * fgFcConservative;  
  TH1D *fhYieldCorr;                     
  TH1D *fhYieldCorrMax;                  
  TH1D *fhYieldCorrMin;                  
  TH2D *fhYieldCorrRcb;                  
  TGraphAsymmErrors * fgYieldCorr;              
  TGraphAsymmErrors * fgYieldCorrExtreme;       
  TGraphAsymmErrors * fgYieldCorrConservative;  
  TH1D *fhSigmaCorr;                     
  TH1D *fhSigmaCorrMax;                  
  TH1D *fhSigmaCorrMin;                  
  TH1D *fhSigmaCorrDataSyst;             
  TH2D *fhSigmaCorrRcb;                  
  TGraphAsymmErrors * fgSigmaCorr;              
  TGraphAsymmErrors * fgSigmaCorrExtreme;       
  TGraphAsymmErrors * fgSigmaCorrConservative;  
  //
  TNtuple *fnSigma;          
  TNtuple *fnHypothesis;     
  Int_t fCollisionType;      
  //
  Int_t fFeedDownOption;            
  Bool_t fAsymUncertainties;        
  Bool_t fPbPbElossHypothesis;      
  Bool_t fIsStatUncEff;             
  Int_t fParticleAntiParticle;      
  Bool_t fIsEventPlane;             

  Int_t     fnPtBins;      
  Double_t *fPtBinLimits;  //[fnPtBins+1] limits of the pt intervals
  Double_t *fPtBinWidths;  //[fnPtBins] width of the pt intervals

  //
  TH1D *fhStatUncEffcSigma;       
  TH1D *fhStatUncEffbSigma;       
  TH1D *fhStatUncEffcFD;          
  TH1D *fhStatUncEffbFD;          

    ClassDef(AliHFPtSpectrum,7); 
};

#endif