merging post-CVS changes into SVN:
[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 // implements a several unfolding methods: e.g. chi2 minimization and bayesian unfolding
12 //
13
14 class TH1;
15 class TH2;
16 class TH1F;
17 class TH2F;
18 class TH3F;
19 class TF1;
20 class TCollection;
21
22 // defined here, because it does not seem possible to predeclare these (or i do not know how)
23 // -->
24 // $ROOTSYS/include/TVectorDfwd.h:21: conflicting types for `typedef struct TVectorT<Double_t> TVectorD'
25 // PWG0/dNdEta/AliMultiplicityCorrection.h:21: previous declaration as `struct TVectorD'
26
27 #include <TMatrixD.h>
28 #include <TVectorD.h>
29
30 class AliMultiplicityCorrection : public TNamed {
31   public:
32     enum EventType { kTrVtx = 0, kMB, kINEL };
33     enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature };
34     enum MethodType { kChi2Minimization = 0, kBayesian = 1 };
35     enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8, kQualityRegions = 3 };
36
37     AliMultiplicityCorrection();
38     AliMultiplicityCorrection(const Char_t* name, const Char_t* title);
39     virtual ~AliMultiplicityCorrection();
40
41     virtual Long64_t Merge(TCollection* list);
42
43     void FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20);
44     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);
45
46     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);
47
48     Bool_t LoadHistograms(const Char_t* dir = 0);
49     void SaveHistograms(const char* dir = 0);
50     void DrawHistograms();
51     void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE);
52
53     Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* initialConditions = 0);
54
55     static void SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams = -1);
56     static void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
57     static void SetCreateBigBin(Bool_t flag) { fgCreateBigBin = flag; }
58
59     void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
60
61     void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Bool_t determineError = kTRUE);
62
63     static TH1* CalculateStdDev(TH1** results, Int_t max);
64     TH1* StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0);
65
66     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
67
68     void ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
69
70     TH2F* GetMultiplicityESD(Int_t i) { return fMultiplicityESD[i]; }
71     TH2F* GetMultiplicityVtx(Int_t i) { return fMultiplicityVtx[i]; }
72     TH2F* GetMultiplicityMB(Int_t i) { return fMultiplicityMB[i]; }
73     TH2F* GetMultiplicityINEL(Int_t i) { return fMultiplicityINEL[i]; }
74     TH2F* GetMultiplicityMC(Int_t i, EventType eventType);
75     TH3F* GetCorrelation(Int_t i) { return fCorrelation[i]; }
76     TH1F* GetMultiplicityESDCorrected(Int_t i) { return fMultiplicityESDCorrected[i]; }
77
78     void SetMultiplicityESD(Int_t i, TH2F* hist)  { fMultiplicityESD[i]  = hist; }
79     void SetMultiplicityVtx(Int_t i, TH2F* hist)  { fMultiplicityVtx[i]  = hist; }
80     void SetMultiplicityMB(Int_t i, TH2F* hist)   { fMultiplicityMB[i]   = hist; }
81     void SetMultiplicityINEL(Int_t i, TH2F* hist) { fMultiplicityINEL[i] = hist; }
82     void SetCorrelation(Int_t i, TH3F* hist) { fCorrelation[i] = hist; }
83     void SetMultiplicityESDCorrected(Int_t i, TH1F* hist) { fMultiplicityESDCorrected[i] = hist; }
84
85     void SetGenMeasFromFunc(TF1* inputMC, Int_t id);
86     TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap);
87
88     void GetComparisonResults(Float_t* mc = 0, Int_t* mcLimit = 0, Float_t* residuals = 0, Float_t* ratioAverage = 0) const;
89
90     TH1* GetEfficiency(Int_t inputRange, EventType eventType);
91
92     static void SetQualityRegions(Bool_t SPDStudy);
93     Float_t GetQuality(Int_t region) const { return fQuality[region]; }
94
95     void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y);
96
97   protected:
98     static const Int_t fgkMaxParams;  //! bins in unfolded histogram = number of fit params
99     static const Int_t fgkMaxInput;   //! bins in measured histogram
100
101     static Double_t RegularizationPol0(TVectorD& params);
102     static Double_t RegularizationPol1(TVectorD& params);
103     static Double_t RegularizationTotalCurvature(TVectorD& params);
104     static Double_t RegularizationEntropy(TVectorD& params);
105     static Double_t RegularizationLog(TVectorD& params);
106
107     static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
108     static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
109
110     static void DrawGuess(Double_t *params);
111
112     void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
113
114     Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
115     static Int_t UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations);
116     static Int_t UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
117
118     TH1* fCurrentESD;         //! static variable to be accessed by MINUIT
119     TH1* fCurrentCorrelation; //! static variable to be accessed by MINUIT
120     TH1* fCurrentEfficiency;  //! static variable to be accessed by MINUIT
121
122     // static variable to be accessed by MINUIT
123     static TMatrixD* fgCorrelationMatrix;            //! contains fCurrentCorrelation in matrix form
124     static TMatrixD* fgCorrelationCovarianceMatrix;  //! contains the errors of fCurrentESD
125     static TVectorD* fgCurrentESDVector;             //! contains fCurrentESD
126     static TVectorD* fgEntropyAPriori;               //! a-priori distribution for entropy regularization
127
128     static TF1* fgNBD;   //! negative binomial distribution
129
130     // configuration params follow
131     static RegularizationType fgRegularizationType; //! regularization that is used during Chi2 method
132     static Float_t fgRegularizationWeight;          //! factor for regularization term
133     static Bool_t  fgCreateBigBin;                  //! to fix fluctuations at high multiplicities, all entries above a certain limit are summarized in one bin
134     static Int_t   fgNParamsMinuit;                 //! number of parameters minuit uses for unfolding (todo: to be merged w/ fgkMaxParams that has to be const. for the moment)
135
136     static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
137     static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
138     // end of configuration
139
140     TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2 (0..3)
141
142     TH2F* fMultiplicityVtx[kMCHists];  // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2, inf (0..4)
143     TH2F* fMultiplicityMB[kMCHists];   // multiplicity histogram of triggered events                        : vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2, inf (0..4)
144     TH2F* fMultiplicityINEL[kMCHists]; // multiplicity histogram of all (inelastic) events                  : vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2, inf (0..4)
145
146     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
147
148     TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms
149
150     Int_t fLastBinLimit;        //! last bin limit, determined in SetupCurrentHists()
151     Float_t fLastChi2MC;        //! last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
152     Int_t   fLastChi2MCLimit;   //! bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
153     Float_t fLastChi2Residuals; //! last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
154     Float_t fRatioAverage;      //! last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
155
156     static Int_t   fgQualityRegionsB[kQualityRegions]; //! begin, given in multiplicity units
157     static Int_t   fgQualityRegionsE[kQualityRegions]; //! end
158     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
159
160  private:
161     AliMultiplicityCorrection(const AliMultiplicityCorrection&);
162     AliMultiplicityCorrection& operator=(const AliMultiplicityCorrection&);
163
164   ClassDef(AliMultiplicityCorrection, 2);
165 };
166
167 #endif
168