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