New book-keeping class
[u/mrichter/AliRoot.git] / 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 class TH3F;
17 class TF1;
18 class TF2;
19 class TCollection;
20
21 #include "TNamed.h"
22 #include <TH2F.h>
23 #include <THnSparse.h>
24
25 class AliJetSpectrumUnfolding : public TNamed {
26   public:
27     AliJetSpectrumUnfolding();
28     AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title);
29     virtual ~AliJetSpectrumUnfolding();
30
31     virtual Long64_t Merge(TCollection* list);
32
33     Bool_t LoadHistograms(const Char_t* dir = 0);
34     void SaveHistograms();
35     void DrawHistograms();
36     void DrawComparison(const char* name, TH2* const genSpec);
37
38     void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
39     void ApplyBayesianMethod(Float_t regPar = 1, Int_t nIterations = 100, TH2* initialConditions = 0, Bool_t determineError = kTRUE);
40
41     TH2F* GetRecSpectrum() const      { return fRecSpectrum; }
42     TH2F* GetGenSpectrum() const     { return fGenSpectrum; }
43     TH2F* GetUnfSpectrum() const     { return fUnfSpectrum; }
44     THnSparseF* GetCorrelation() const { return fCorrelation; }
45
46     void SetRecSpectrum(TH2F* const hist)   { if(fRecSpectrum) delete fRecSpectrum;
47       fRecSpectrum  = hist; }
48     void SetGenSpectrum(TH2F* const hist)      { if(fGenSpectrum)delete fGenSpectrum;
49       fGenSpectrum  = hist; }
50     void SetUnfSpectrum(TH2F*  const hist)      { if(fUnfSpectrum)delete fUnfSpectrum;
51       fUnfSpectrum  = hist; }
52     void SetCorrelation(THnSparseF* const hist){ if(fCorrelation)delete fCorrelation;
53       fCorrelation  = hist; }
54
55     void SetGenRecFromFunc(TF2* inputGen);
56     TH2F* CalculateRecSpectrum(TH2* inputGen);
57
58     static void NormalizeToBinWidth(TH2* const hist);
59
60     void GetComparisonResults(Float_t* const gen = 0, Int_t* const genLimit = 0, Float_t* const residuals = 0, Float_t* const ratioAverage = 0) const;
61
62   protected:
63     void SetupCurrentHists(Bool_t createBigBin);
64
65     static Double_t BayesCov(THnSparseF* const M, THnSparseF* const correlation, Int_t* const binTM, Int_t* const binTM1);
66     static Double_t BayesUncertaintyTerms(THnSparseF* const M, THnSparseF *const C, Int_t* const binTM, Int_t* const binTM1, Double_t nt);
67     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);
68
69     static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
70     static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
71
72     TH2F* fCurrentRec;
73     THnSparseF* fCurrentCorrelation;
74
75     TH2F* fRecSpectrum;
76     TH2F* fGenSpectrum;
77     TH2F* fUnfSpectrum;
78     THnSparseF* fCorrelation;
79
80     Float_t fLastChi2MC;        //! last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
81     Int_t   fLastChi2MCLimit;   //! bin where the last chi2 breached a certain threshold
82     Float_t fLastChi2Residuals; //! last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
83     Float_t fRatioAverage;      //! last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
84
85  private:
86
87     static const Int_t fgkNBINSE;
88     static const Int_t fgkNBINSZ;
89     static const Int_t fgkNEVENTS;
90     static const Double_t fgkaxisLowerLimitE;
91     static const Double_t fgkaxisLowerLimitZ;
92     static const Double_t fgkaxisUpperLimitE;
93     static const Double_t fgkaxisUpperLimitZ;
94
95
96
97     AliJetSpectrumUnfolding(const AliJetSpectrumUnfolding&);
98     AliJetSpectrumUnfolding& operator=(const AliJetSpectrumUnfolding&);
99
100   ClassDef(AliJetSpectrumUnfolding, 2);
101 };
102
103 #endif
104