PWG4/JetTasks/AliJetSpectrumUnfolding.h

   1 #ifndef ALIJETSPECTRUMUNFOLDING_H
   2 #define ALIJETSPECTRUMUNFOLDING_H
   3 //
   4 // class that contains the correction matrix and the functions for
   5 // correction the jet spectrum
   6 // implements 2-dim bayesian method
   7 //
   8
   9
  10
  11
  12 class TH1;
  13 class TH2;
  14 class TH3;
  15 class TH1F;
  16
  17 class TH3F;
  18 class TF1;
  19 class TF2;
  20 class TCollection;
  21
  22 #include "TNamed.h"
  23 #include <THnSparse.h>
  24 #include <TH2F.h> // need to included for delete
  25
  26
  27 class AliJetSpectrumUnfolding : public TNamed {
  28   public:
  29     AliJetSpectrumUnfolding();
  30     AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title);
  31     virtual ~AliJetSpectrumUnfolding();
  32
  33     virtual Long64_t Merge(TCollection* list);
  34
  35     Bool_t LoadHistograms(const Char_t* dir = 0);
  36     void SaveHistograms();
  37     void DrawHistograms();
  38     void DrawComparison(const char* name, TH2* const genSpec);
  39
  40     void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
  41     void ApplyBayesianMethod(Float_t regPar = 1, Int_t nIterations = 100, TH2* initialConditions = 0, Bool_t determineError = kTRUE);
  42
  43     TH2F* GetRecSpectrum() const      { return fRecSpectrum; }
  44     TH2F* GetGenSpectrum() const     { return fGenSpectrum; }
  45     TH2F* GetUnfSpectrum() const     { return fUnfSpectrum; }
  46     THnSparseF* GetCorrelation() const { return fCorrelation; }
  47
  48     void SetRecSpectrum(TH2F* const hist)   { if(fRecSpectrum) delete fRecSpectrum;
  49       fRecSpectrum  = hist; }
  50     void SetGenSpectrum(TH2F* const hist)      { if(fGenSpectrum)delete fGenSpectrum;
  51       fGenSpectrum  = hist; }
  52     void SetUnfSpectrum(TH2F*  const hist)      { if(fUnfSpectrum)delete fUnfSpectrum;
  53       fUnfSpectrum  = hist; }
  54     void SetCorrelation(THnSparseF* const hist){ if(fCorrelation)delete fCorrelation;
  55       fCorrelation  = hist; }
  56
  57     void SetGenRecFromFunc(const TF2* const inputGen);
  58     TH2F* CalculateRecSpectrum(TH2* inputGen);
  59
  60     static void NormalizeToBinWidth(TH2* const hist);
  61
  62
  63   protected:
  64     void SetupCurrentHists(Bool_t createBigBin);
  65
  66     static Double_t BayesCov(THnSparseF* const M, THnSparseF* const correlation,const  Int_t* const binTM,const Int_t* const binTM1);
  67     static Double_t BayesUncertaintyTerms(THnSparseF* const M, THnSparseF *const C,const Int_t* const binTM,const Int_t* const binTM1, Double_t nt);
  68     static Int_t UnfoldWithBayesian(THnSparseF* const correlation, TH2* const measured, TH2* const initialConditions, TH2* const aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors = kFALSE);
  69
  70     static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
  71     static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
  72
  73     TH2F* fCurrentRec;   // current rec
  74     THnSparseF* fCurrentCorrelation; // current correlat
  75
  76     TH2F* fRecSpectrum; // rec spectrum
  77     TH2F* fGenSpectrum; // gen spectrum
  78     TH2F* fUnfSpectrum; // Unfolded spectrum
  79     THnSparseF* fCorrelation; // corealtion matrix
  80
  81     Float_t fLastChi2MC;        //! last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
  82     Int_t   fLastChi2MCLimit;   //! bin where the last chi2 breached a certain threshold
  83     Float_t fLastChi2Residuals; //! last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
  84     Float_t fRatioAverage;      //! last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
  85
  86  private:
  87
  88     static const Int_t fgkNBINSE; // bins energy
  89     static const Int_t fgkNBINSZ; // bins Z
  90     static const Int_t fgkNEVENTS; // bins events
  91     static const Double_t fgkaxisLowerLimitE; // lower limit e
  92     static const Double_t fgkaxisLowerLimitZ; // lower limit Z
  93     static const Double_t fgkaxisUpperLimitE; // upper limit E
  94     static const Double_t fgkaxisUpperLimitZ; // upper limit Z
  95
  96
  97
  98     AliJetSpectrumUnfolding(const AliJetSpectrumUnfolding&);
  99     AliJetSpectrumUnfolding& operator=(const AliJetSpectrumUnfolding&);
 100
 101   ClassDef(AliJetSpectrumUnfolding, 2);
 102 };
 103
 104 #endif
 105