Added some more scripts
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.h
CommitLineData
f8715167 1#ifndef ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTER_H
2#define ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTER_H
3#include <TNamed.h>
4#include <TH1D.h>
5#include <TAxis.h>
6#include <TList.h>
7#include <TObjArray.h>
8#include "AliForwardUtil.h"
9class AliESDFMD;
10class TFitResult;
11class TF1;
12class TArrayD;
13
14/**
15 * Class to fit the energy distribution.
16 *
17 * @par Input:
18 * - AliESDFMD object - from reconstruction
19 *
20 * @par Output:
21 * - Lists of histogram - one per ring. Each list has a number of
22 * histograms corresponding to the number of eta bins defined.
23 *
24 * @par Corrections used:
25 * - None
26 *
27 *
28 * @ingroup pwg2_forward_analysis
29 */
30class AliFMDEnergyFitter : public TNamed
31{
32public:
c389303e 33 enum {
34 kC = AliForwardUtil::ELossFitter::kC,
35 kDelta = AliForwardUtil::ELossFitter::kDelta,
36 kXi = AliForwardUtil::ELossFitter::kXi,
37 kSigma = AliForwardUtil::ELossFitter::kSigma,
38 kSigmaN = AliForwardUtil::ELossFitter::kSigmaN,
39 kN = AliForwardUtil::ELossFitter::kN,
40 kA = AliForwardUtil::ELossFitter::kA
41 };
42
f8715167 43 /**
44 * Destructor
45 */
46 virtual ~AliFMDEnergyFitter();
47 /**
48 * Default Constructor - do not use
49 */
50 AliFMDEnergyFitter();
51 /**
52 * Constructor
53 *
54 * @param title Title of object - not significant
55 */
56 AliFMDEnergyFitter(const char* title);
57 /**
58 * Copy constructor
59 *
60 * @param o Object to copy from
61 */
62 AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
63 /**
64 * Assignment operator
65 *
66 * @param o Object to assign from
67 *
68 * @return Reference to this
69 */
70 AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
71
4b9857f3 72 /**
73 * Initialise the task
74 *
75 * @param etaAxis The eta axis to use. Note, that if the eta axis
76 * has already been set (using SetEtaAxis), then this parameter will be
77 * ignored
78 */
f8715167 79 void Init(const TAxis& etaAxis);
4b9857f3 80 /**
81 * Set the eta axis to use. This will force the code to use this
82 * eta axis definition - irrespective of whatever axis is passed to
83 * the Init member function. Therefore, this member function can be
84 * used to force another eta axis than one found in the correction
85 * objects.
86 *
87 * @param nBins Number of bins
88 * @param etaMin Minimum of the eta axis
89 * @param etaMax Maximum of the eta axis
90 */
91 void SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax);
92 /**
93 * Set the eta axis to use. This will force the code to use this
94 * eta axis definition - irrespective of whatever axis is passed to
95 * the Init member function. Therefore, this member function can be
96 * used to force another eta axis than one found in the correction
97 * objects.
98 *
99 * @param etaAxis Eta axis to use
100 */
101 void SetEtaAxis(const TAxis& etaAxis);
f8715167 102 /**
103 * Set the low cut used for energy
104 *
105 * @param lowCut Low cut
106 */
107 void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
4b9857f3 108 /**
109 * Set the number of bins to subtract
110 *
111 * @param n
112 */
c389303e 113 void SetFitRangeBinWidth(UShort_t n=4) { fFitRangeBinWidth = n; }
f8715167 114 /**
115 * Whether or not to enable fitting of the final merged result.
116 * Note, fitting takes quite a while and one should be careful not to do
117 * this needlessly
118 *
119 * @param doFit Whether to do the fits or not
120 */
c389303e 121 void SetDoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
f8715167 122 /**
c389303e 123 * Set how many particles we will try to fit at most to the data
f8715167 124 *
c389303e 125 * @param n Max number of particle to try to fit
f8715167 126 */
c389303e 127 void SetNParticles(UShort_t n) { fNParticles = (n < 1 ? 1 : (n > 5 ? 5 : n)); }
f8715167 128 /**
c389303e 129 * Set the minimum number of entries each histogram must have
130 * before we try to fit our response function to it
f8715167 131 *
c389303e 132 * @param n Minimum number of entries
f8715167 133 */
134 void SetMinEntries(UShort_t n) { fMinEntries = (n < 1 ? 1 : n); }
4b9857f3 135 /**
136 * Set maximum energy loss to consider
137 *
138 * @param x Maximum energy loss to consider
139 */
140 void SetMaxE(Double_t x) { fMaxE = x; }
141 /**
142 * Set number of energy loss bins
143 *
144 * @param x Number of energy loss bins
145 */
146 void SetNEbins(Int_t x) { fNEbins = x; }
c389303e 147 void SetMaxRelativeParameterError(Double_t e) { fMaxRelParError = e; }
148 void SetMaxChi2PerNDF(Double_t c) { fMaxChi2PerNDF = c; }
4b9857f3 149 /**
150 * Set wheter to use increasing bin sizes
151 *
152 * @param x Wheter to use increasing bin sizes
153 */
154 void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
f8715167 155 /**
156 * Fitter the input AliESDFMD object
157 *
158 * @param input Input
159 * @param empty Whether the event is 'empty'
160 *
161 * @return True on success, false otherwise
162 */
163 Bool_t Accumulate(const AliESDFMD& input,
164 Bool_t empty);
165 /**
166 * Scale the histograms to the total number of events
167 *
c389303e 168 * @param dir Where the histograms are
f8715167 169 */
170 void Fit(TList* dir);
171
172 /**
173 * Define the output histograms. These are put in a sub list of the
174 * passed list. The histograms are merged before the parent task calls
175 * AliAnalysisTaskSE::Terminate
176 *
177 * @param dir Directory to add to
178 */
179 void DefineOutput(TList* dir);
180 /**
181 * Set the debug level. The higher the value the more output
182 *
183 * @param dbg Debug level
184 */
185 void SetDebug(Int_t dbg=1);
186protected:
187 /**
188 * Internal data structure to keep track of the histograms
189 */
190 struct RingHistos : public AliForwardUtil::RingHistos
191 {
192 /**
193 * Default CTOR
194 */
195 RingHistos();
196 /**
197 * Constructor
198 *
199 * @param d detector
200 * @param r ring
201 */
202 RingHistos(UShort_t d, Char_t r);
203 /**
204 * Copy constructor
205 *
206 * @param o Object to copy from
207 */
208 RingHistos(const RingHistos& o);
209 /**
210 * Assignment operator
211 *
212 * @param o Object to assign from
213 *
214 * @return Reference to this
215 */
216 RingHistos& operator=(const RingHistos& o);
217 /**
218 * Destructor
219 */
220 ~RingHistos();
221 /**
222 * Define outputs
223 *
224 * @param dir
225 */
226 void Output(TList* dir);
227 /**
228 * Initialise object
229 *
c389303e 230 * @param eAxis Eta axis
231 * @param maxDE Max energy loss to consider
232 * @param nDEbins Number of bins
233 * @param useIncrBin Whether to use an increasing bin size
f8715167 234 */
66b34429 235 void Init(const TAxis& eAxis,
4b9857f3 236 Double_t maxDE=10,
237 Int_t nDEbins=300,
238 Bool_t useIncrBin=true);
f8715167 239 /**
240 * Fill histogram
241 *
242 * @param empty True if event is empty
243 * @param ieta Eta bin
244 * @param mult Signal
245 */
246 void Fill(Bool_t empty, Int_t ieta, Double_t mult);
247 /**
c389303e 248 * Fit each histogram to up to @a nParticles particle responses.
f8715167 249 *
c389303e 250 * @param dir Output list
251 * @param eta Eta axis
252 * @param lowCut Lower cut
253 * @param nParticles Max number of convolved landaus to fit
254 * @param minEntries Minimum number of entries
255 * @param minusBins Number of bins from peak to subtract to
256 * get the fit range
257 * @param relErrorCut Cut applied to relative error of parameter.
258 * Note, for multi-particle weights, the cut
259 * is loosend by a factor of 2
260 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
261 * the reduced @f$\chi^2@f$
f8715167 262 */
4b9857f3 263 TObjArray* Fit(TList* dir,
264 const TAxis& eta,
265 Double_t lowCut,
c389303e 266 UShort_t nParticles,
4b9857f3 267 UShort_t minEntries,
c389303e 268 UShort_t minusBins,
269 Double_t relErrorCut,
270 Double_t chi2nuCut) const;
f8715167 271 /**
c389303e 272 * Fit a signal histogram. First, the bin @f% b_{min}@f$ with
273 * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
274 * found. Then the fit range is set to the bin range
275 * @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
276 * particle signal is fitted to that. The parameters of that fit
277 * is then used as seeds for a fit of the @f$ N@f$ particle response
278 * to the data in the range
279 * @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
f8715167 280 *
c389303e 281 * @param dist Histogram to fit
282 * @param lowCut Lower cut @f$ E_{min}@f$ on signal
283 * @param nParticles Max number @f$ N@f$ of convolved landaus to fit
284 * @param minusBins Number of bins @f$ \Delta b@f$ from peak to
285 * subtract to get the fit range
286 * @param relErrorCut Cut applied to relative error of parameter.
287 * Note, for multi-particle weights, the cut
288 * is loosend by a factor of 2
289 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
290 * the reduced @f$\chi^2@f$
f8715167 291 *
292 * @return The best fit function
293 */
4b9857f3 294 TF1* FitHist(TH1* dist,
295 Double_t lowCut,
c389303e 296 UShort_t nParticles,
297 UShort_t minusBins,
298 Double_t relErrorCut,
299 Double_t chi2nuCut) const;
f8715167 300 /**
4b9857f3 301 * Check the result of the fit. Returns true if
c389303e 302 * - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
303 * - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters. Note,
304 * for multi-particle fits, this requirement is relaxed by a
305 * factor of 2
306 * - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$
307 * particle response
f8715167 308 *
c389303e 309 * @param r Result to check
310 * @param relErrorCut Cut @f$ \delta_e@f$ applied to relative error
311 * of parameter.
312 * @param chi2nuCut Cut @f$ \max{\chi^2/\nu}@f$
f8715167 313 *
4b9857f3 314 * @return true if fit is good.
f8715167 315 */
c389303e 316 Bool_t CheckResult(TFitResult* r,
317 Double_t relErrorCut,
318 Double_t chi2nuCut) const;
f8715167 319 /**
4b9857f3 320 * Make an axis with increasing bins
f8715167 321 *
4b9857f3 322 * @param n Number of bins
323 * @param min Minimum
324 * @param max Maximum
f8715167 325 *
4b9857f3 326 * @return An axis with quadratically increasing bin size
f8715167 327 */
66b34429 328 TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
f8715167 329 /**
330 * Make E/E_mip histogram
331 *
332 * @param ieta Eta bin
333 * @param eMin Least signal
334 * @param eMax Largest signal
335 */
66b34429 336 void Make(Int_t ieta, Double_t eMin, Double_t eMax,
337 Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
f8715167 338 /**
339 * Make a parameter histogram
340 *
341 * @param name Name of histogram.
342 * @param title Title of histogram.
343 * @param eta Eta axis
344 *
345 * @return
346 */
347 TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
4b9857f3 348 /**
349 * Make a histogram that contains the results of the fit over the full ring
350 *
351 * @param name Name
352 * @param title Title
353 * @param eta Eta axis
354 * @param low Least bin
355 * @param high Largest bin
356 * @param val Value of parameter
357 * @param err Error on parameter
358 *
359 * @return The newly allocated histogram
360 */
f8715167 361 TH1D* MakeTotal(const char* name,
362 const char* title,
363 const TAxis& eta,
364 Int_t low,
365 Int_t high,
366 Double_t val,
367 Double_t err) const;
368 TH1D* fEDist; // Ring energy distribution
369 TH1D* fEmpty; // Ring energy distribution for empty events
370 TList fEtaEDists; // Energy distributions per eta bin.
371 TList* fList;
372 Int_t fDebug;
373 ClassDef(RingHistos,1);
374 };
375 /**
376 * Get the ring histogram container
377 *
378 * @param d Detector
379 * @param r Ring
380 *
381 * @return Ring histogram container
382 */
383 RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
384
385 TList fRingHistos; // List of histogram containers
386 Double_t fLowCut; // Low cut on energy
c389303e 387 UShort_t fNParticles; // Number of landaus to try to fit
f8715167 388 UShort_t fMinEntries; // Minimum number of entries
c389303e 389 UShort_t fFitRangeBinWidth;// Number of bins to subtract from found max
f8715167 390 Bool_t fDoFits; // Wheter to actually do the fits
391 TAxis fEtaAxis; // Eta axis
4b9857f3 392 Double_t fMaxE; // Maximum energy loss to consider
393 Int_t fNEbins; // Number of energy loss bins
394 Bool_t fUseIncreasingBins; // Wheter to use increasing bin sizes
c389303e 395 Double_t fMaxRelParError;// Relative error cut
396 Double_t fMaxChi2PerNDF; // chi^2/nu cit
f8715167 397 Int_t fDebug; // Debug level
398
399 ClassDef(AliFMDEnergyFitter,1); //
400};
401
402#endif
403// Local Variables:
404// mode: C++
405// End: