3da2d88fd3bc29732d18d3c90baf5b92009a015f
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityCorrection.h
1 /* $Id$ */
2
3 #ifndef ALIMULTIPLICITYCORRECTION_H
4 #define ALIMULTIPLICITYCORRECTION_H
5
6 #include "TNamed.h"
7
8 //
9 // class that contains the correction matrix and the functions for
10 // correction the multiplicity spectrum
11 //
12
13 class TH1;
14 class TH2;
15 class TH1F;
16 class TH2F;
17 class TH3F;
18 class TF1;
19 class TCollection;
20
21 #include <TMatrixD.h>
22 #include <TVectorD.h>
23
24 class AliMultiplicityCorrection : public TNamed {
25   public:
26     enum EventType { kTrVtx = 0, kMB, kINEL };
27     enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature };
28     enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8, kQualityRegions = 3 };
29
30     AliMultiplicityCorrection();
31     AliMultiplicityCorrection(const Char_t* name, const Char_t* title);
32     virtual ~AliMultiplicityCorrection();
33
34     virtual Long64_t Merge(TCollection* list);
35
36     void FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20);
37     void FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll);
38
39     void FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20);
40
41     Bool_t LoadHistograms(const Char_t* dir = 0);
42     void SaveHistograms();
43     void DrawHistograms();
44     void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE);
45
46     Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0);
47     void SetRegularizationParameters(RegularizationType type, Float_t weight);
48
49     void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
50
51     void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* inputDist = 0, Bool_t determineError = kTRUE);
52     TH1* BayesianStatisticsEffect(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, Float_t regPar = 1, Int_t nIterations = 100, TH1* compareTo = 0);
53
54     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
55
56     void ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
57
58     TH2F* GetMultiplicityESD(Int_t i) { return fMultiplicityESD[i]; }
59     TH2F* GetMultiplicityVtx(Int_t i) { return fMultiplicityVtx[i]; }
60     TH2F* GetMultiplicityMB(Int_t i) { return fMultiplicityMB[i]; }
61     TH2F* GetMultiplicityINEL(Int_t i) { return fMultiplicityINEL[i]; }
62     TH2F* GetMultiplicityMC(Int_t i, EventType eventType);
63     TH3F* GetCorrelation(Int_t i) { return fCorrelation[i]; }
64     TH1F* GetMultiplicityESDCorrected(Int_t i) { return fMultiplicityESDCorrected[i]; }
65
66     void SetMultiplicityESD(Int_t i, TH2F* hist)  { fMultiplicityESD[i]  = hist; }
67     void SetMultiplicityVtx(Int_t i, TH2F* hist)  { fMultiplicityVtx[i]  = hist; }
68     void SetMultiplicityMB(Int_t i, TH2F* hist)   { fMultiplicityMB[i]   = hist; }
69     void SetMultiplicityINEL(Int_t i, TH2F* hist) { fMultiplicityINEL[i] = hist; }
70     void SetCorrelation(Int_t i, TH3F* hist) { fCorrelation[i] = hist; }
71     void SetMultiplicityESDCorrected(Int_t i, TH1F* hist) { fMultiplicityESDCorrected[i] = hist; }
72
73     void SetGenMeasFromFunc(TF1* inputMC, Int_t id);
74     TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap);
75
76     static void NormalizeToBinWidth(TH1* hist);
77     static void NormalizeToBinWidth(TH2* hist);
78
79     void GetComparisonResults(Float_t* mc = 0, Int_t* mcLimit = 0, Float_t* residuals = 0, Float_t* ratioAverage = 0);
80
81     TH1* GetEfficiency(Int_t inputRange, EventType eventType);
82
83     static void SetQualityRegions(Bool_t SPDStudy);
84     Float_t GetQuality(Int_t region) { return fQuality[region]; }
85
86     void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y);
87
88   protected:
89     static const Int_t fgMaxParams;  // bins in unfolded histogram = number of fit params
90     static const Int_t fgMaxInput;   // bins in measured histogram
91
92     static Double_t RegularizationPol0(TVectorD& params);
93     static Double_t RegularizationPol1(TVectorD& params);
94     static Double_t RegularizationTotalCurvature(TVectorD& params);
95     static Double_t RegularizationEntropy(TVectorD& params);
96     static Double_t RegularizationLog(TVectorD& params);
97
98     static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
99     static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
100
101     void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
102
103     Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
104     TH1* UnfoldWithBayesian(TH1* measured, Float_t regPar, Int_t nIterations, TH1* inputDist);
105
106     static TH1* fCurrentESD;         // static variable to be accessed by MINUIT
107     static TH1* fCurrentCorrelation; // static variable to be accessed by MINUIT
108     static TH1* fCurrentEfficiency;  // static variable to be accessed by MINUIT
109
110     static TMatrixD* fCorrelationMatrix;            // contains fCurrentCorrelation in matrix form
111     static TMatrixD* fCorrelationCovarianceMatrix;  // contains the errors of fCurrentESD
112     static TVectorD* fCurrentESDVector;             // contains fCurrentESD
113     static TVectorD* fEntropyAPriori;               // a-priori distribution for entropy regularization
114
115     static TF1* fNBD;   // negative binomial distribution
116
117     static RegularizationType fRegularizationType; // regularization that is used during Chi2 method
118     static Float_t fRegularizationWeight;          // factor for regularization term
119
120     TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2 (0..3)
121     TH2F* fMultiplicityVtx[kMCHists];  // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2, inf (0..4)
122     TH2F* fMultiplicityMB[kMCHists];   // multiplicity histogram of triggered events                        : vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2, inf (0..4)
123     TH2F* fMultiplicityINEL[kMCHists]; // multiplicity histogram of all (inelastic) events                  : vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2, inf (0..4)
124
125     TH3F* fCorrelation[kCorrHists];              // vtx vs. (gene multiplicity (trig+vtx)) vs. (meas multiplicity); array: |eta| < 0.5, 1, 1.5, 2 (0..3 and 4..7), the first corrects to the eta range itself, the second to full phase space
126     TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms
127
128     Float_t fLastChi2MC;        // last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
129     Int_t   fLastChi2MCLimit;   // bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
130     Float_t fLastChi2Residuals; // last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
131     Float_t fRatioAverage;      // last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
132
133     static Int_t   fgQualityRegionsB[kQualityRegions]; // begin, given in multiplicity units
134     static Int_t   fgQualityRegionsE[kQualityRegions]; // end
135     Float_t fQuality[kQualityRegions];        // stores the quality of the last comparison (calculated in DrawComparison). Contains 3 values that are averages of (MC - unfolded) / e(MC) in 3 regions, these are defined in fQualityRegionB,E
136
137  private:
138     AliMultiplicityCorrection(const AliMultiplicityCorrection&);
139     AliMultiplicityCorrection& operator=(const AliMultiplicityCorrection&);
140
141   ClassDef(AliMultiplicityCorrection, 1);
142 };
143
144 #endif
145