2 // Class to fit the energy distribution.
4 #ifndef ALIFMDENERGYFITTER_H
5 #define ALIFMDENERGYFITTER_H
7 * @file AliFMDEnergyFitter.h
8 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
9 * @date Wed Mar 23 14:02:23 2011
14 * @ingroup pwglf_forward_eloss
20 #include <TObjArray.h>
21 #include <TClonesArray.h>
22 #include "AliFMDCorrELossFit.h"
23 #include "AliForwardUtil.h"
32 * Class to fit the energy distribution.
35 * - AliESDFMD object - from reconstruction
38 * - Lists of histogram - one per ring. Each list has a number of
39 * histograms corresponding to the number of eta bins defined.
41 * @par Corrections used:
44 * @image html alice-int-2012-040-eloss_fits.png "Summary of fits"
46 * @ingroup pwglf_forward_algo
47 * @ingroup pwglf_forward_eloss
49 class AliFMDEnergyFitter : public TNamed
53 * Enumeration of parameters
56 /** Index of pre-constant @f$ C@f$ */
57 kC = AliForwardUtil::ELossFitter::kC,
58 /** Index of most probable value @f$ \Delta_p@f$ */
59 kDelta = AliForwardUtil::ELossFitter::kDelta,
60 /** Index of Landau width @f$ \xi@f$ */
61 kXi = AliForwardUtil::ELossFitter::kXi,
62 /** Index of Gaussian width @f$ \sigma@f$ */
63 kSigma = AliForwardUtil::ELossFitter::kSigma,
64 /** Index of Gaussian additional width @f$ \sigma_n@f$ */
65 kSigmaN = AliForwardUtil::ELossFitter::kSigmaN,
66 /** Index of Number of particles @f$ N@f$ */
67 kN = AliForwardUtil::ELossFitter::kN,
68 /** Base index of particle strengths @f$ a_i@f$ for
70 kA = AliForwardUtil::ELossFitter::kA
73 * Enumeration of residual methods
75 enum EResidualMethod {
76 /** Do not calculate residuals */
78 /** The residuals stored are the difference, and the errors are
79 stored in the error bars of the histogram. */
81 /** The residuals stored are the differences scaled to the error
83 kResidualScaledDifference,
84 /** The residuals stored are the square difference scale to the
85 square error on the data. */
86 kResidualSquareDifference
90 * FMD ring bits for skipping
102 kFMD2 =kFMD2I|kFMD2O,
113 virtual ~AliFMDEnergyFitter();
115 * Default Constructor - do not use
117 AliFMDEnergyFitter();
121 * @param title Title of object - not significant
123 AliFMDEnergyFitter(const char* title);
125 // -----------------------------------------------------------------
128 * @name Setters of options and parameters
131 * Set the eta axis to use. This will force the code to use this
132 * eta axis definition - irrespective of whatever axis is passed to
133 * the Init member function. Therefore, this member function can be
134 * used to force another eta axis than one found in the correction
137 * @param nBins Number of bins
138 * @param etaMin Minimum of the eta axis
139 * @param etaMax Maximum of the eta axis
141 void SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax);
143 * Set the eta axis to use. This will force the code to use this
144 * eta axis definition - irrespective of whatever axis is passed to
145 * the Init member function. Therefore, this member function can be
146 * used to force another eta axis than one found in the correction
149 * @param etaAxis Eta axis to use
151 void SetEtaAxis(const TAxis& etaAxis);
153 * Set the centrality bins. E.g.,
156 * Double_t bins[] = { 0., 5., 10., 15., 20., 30.,
157 * 40., 50., 60., 70., 80., 100. };
158 * task->GetFitter().SetCentralityBins(n, bins);
161 * @param nBins Size of @a bins
162 * @param bins Bin limits.
164 void SetCentralityAxis(UShort_t nBins, Double_t* bins);
166 * Set the low cut used for energy
168 * @param lowCut Low cut
170 void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
172 * Set the number of bins to subtract
176 void SetFitRangeBinWidth(UShort_t n=4) { fFitRangeBinWidth = n; }
178 * Whether or not to enable fitting of the final merged result.
179 * Note, fitting takes quite a while and one should be careful not to do
182 * @param doFit Whether to do the fits or not
184 void SetDoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
186 * Set whether to make the corrections object on the output. Note,
187 * fits should be enable for this to have any effect.
189 * @param doMake If true (false is default), do make the corrections object.
191 void SetDoMakeObject(Bool_t doMake=kTRUE) { fDoMakeObject = doMake; }
193 * Set how many particles we will try to fit at most to the data
195 * @param n Max number of particle to try to fit
197 void SetNParticles(UShort_t n) { fNParticles = (n<1 ? 1 : (n>5 ? 5 : n)); }
199 * Set the minimum number of entries each histogram must have
200 * before we try to fit our response function to it
202 * @param n Minimum number of entries
204 void SetMinEntries(UShort_t n) { fMinEntries = (n < 1 ? 1 : n); }
206 * Set maximum energy loss to consider
208 * @param x Maximum energy loss to consider
210 void SetMaxE(Double_t x) { fMaxE = x; }
212 * Set number of energy loss bins
214 * @param x Number of energy loss bins
216 void SetNEbins(Int_t x) { fNEbins = x; }
218 * Set the maximum relative error
220 * @param e Maximum relative error
222 void SetMaxRelativeParameterError(Double_t e=0.2) { fMaxRelParError = e; }
224 * Set the maximum @f$ \chi^2/\nu@f$
226 * @param c Maximum @f$ \chi^2/\nu@f$
228 void SetMaxChi2PerNDF(Double_t c=10) { fMaxChi2PerNDF = c; }
230 * Set the least weight
232 * @param c Least weight
234 void SetMinWeight(Double_t c=1e-7) { fMinWeight = c; }
236 * Set wheter to use increasing bin sizes
238 * @param x Wheter to use increasing bin sizes
240 void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
242 * Set whether to make residuals, and in that case how.
244 * - Square difference: @f$chi_i^2=(h_i - f(x_i))^2/\delta_i^2@f$
245 * - Scaled difference: @f$(h_i - f(x_i))/\delta_i@f$
246 * - Difference: @f$(h_i - f(x_i)) \pm\delta_i@f$
248 * where @f$h_i, x_i, \delta_i@f$ is the bin content, bin center,
249 * and bin error for bin @f$i@f$ respectively, and @f$ f@f$ is the
252 * @param x Residual method
254 void SetStoreResiduals(EResidualMethod x=kResidualDifference)
259 * Set the regularization cut @f$c_{R}@f$. If a @f$\Delta@f$
260 * distribution has more entries @f$ N_{dist}@f$ than @f$c_{R}@f$,
261 * then we modify the errors of the the distribution with the factor
264 * \sqrt{N_{dist}/c_{R}}
267 * to keep the @f$\chi^2/\nu@f$ within resonable limits.
269 * The large residuals @f$chi_i^2=(h_i - f(x_i))^2/\delta_i^2@f$
270 * (see also SetStoreResiduals) comes about on the boundary between
271 * the @f$N@f$ and @f$N+1@f$ particle contributions, and seems to
272 * fall off for larger @f$N@f$. This may indicate that there's a
273 * component in the distributions that the function
276 * f(\Delta;\Delta_p,\xi,\sigma,\mathbf{a}) = \sum_i=1^{n} a_i\int
277 * d\Delta' L(\Delta;\Delta',\xi) G(\Delta';\Delta_p,\sigma)
284 void SetRegularizationCut(Double_t cut=3e6)
286 fRegularizationCut = cut;
288 void SetSkips(UShort_t skip) { fSkips = skip; }
290 * Set the debug level. The higher the value the more output
292 * @param dbg Debug level
294 void SetDebug(Int_t dbg=1);
296 // -----------------------------------------------------------------
303 * Define the output histograms. These are put in a sub list of the
304 * passed list. The histograms are merged before the parent task calls
305 * AliAnalysisTaskSE::Terminate
307 * @param dir Directory to add to
309 virtual void CreateOutputObjects(TList* dir);
311 * Initialise the task
313 * @param etaAxis The eta axis to use. Note, that if the eta axis
314 * has already been set (using SetEtaAxis), then this parameter will be
317 virtual void SetupForData(const TAxis& etaAxis);
319 * Fitter the input AliESDFMD object
322 * @param cent Event centrality (or < 0 if not valid)
323 * @param empty Whether the event is 'empty'
325 * @return True on success, false otherwise
327 virtual Bool_t Accumulate(const AliESDFMD& input,
331 * Scale the histograms to the total number of events
333 * @param dir Where the histograms are
335 virtual void Fit(const TList* dir);
337 * Generate the corrections object
339 * @param dir List to analyse
341 void MakeCorrectionsObject(TList* dir);
346 * @param option Not used
348 void Print(Option_t* option="") const;
350 * Read the parameters from a list - used when re-running the code
352 * @param list Input list
354 * @return true if the parameter where read
356 Bool_t ReadParameters(const TCollection* list);
361 * @param o Object to copy from
363 AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
365 * Assignment operator
367 * @param o Object to assign from
369 * @return Reference to this
371 AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
374 * Internal data structure to keep track of the histograms
376 struct RingHistos : public AliForwardUtil::RingHistos
378 typedef AliFMDCorrELossFit::ELossFit ELossFit_t;
389 RingHistos(UShort_t d, Char_t r);
391 * Copy constructor - not defined
393 * @param o Object to copy from
395 RingHistos(const RingHistos& o);
397 * Assignment operator - not defined
399 * @param o Object to assign from
401 * @return Reference to this
403 RingHistos& operator=(const RingHistos& o);
409 * Make an axis with increasing bins
411 * @param n Number of bins
415 * @return An axis with quadratically increasing bin size
417 virtual TArrayD MakeIncreasingAxis(Int_t n,
421 * Make E/E_mip histogram
423 * @param name Name of histogram
424 * @param title Title of histogram
425 * @param eAxis @f$\eta@f$ axis
426 * @param deMax Maximum energy loss
427 * @param nDeBins Number energy loss bins
428 * @param incr Whether to make bins of increasing size
430 TH2* Make(const char* name,
441 virtual void CreateOutputObjects(TList* dir);
445 * @param eAxis Eta axis
446 * @param cAxis Centrality axis
447 * @param maxDE Max energy loss to consider
448 * @param nDEbins Number of bins
449 * @param useIncrBin Whether to use an increasing bin size
451 virtual void SetupForData(const TAxis& eAxis,
455 Bool_t useIncrBin=true);
459 * @param empty True if event is empty
460 * @param eta @f$ Eta@f$
461 * @param icent Centrality bin (1 based)
464 virtual void Fill(Bool_t empty, Double_t eta, Int_t icent, Double_t mult);
466 * Get the the 2D histogram eloss name from our sub-list of @a dir
467 * and call the Fit function described below (with &fBest) as last
470 * @param dir Output list
471 * @param lowCut Lower cut
472 * @param nParticles Max number of convolved landaus to fit
473 * @param minEntries Minimum number of entries
474 * @param minusBins Number of bins from peak to subtract to
476 * @param relErrorCut Cut applied to relative error of parameter.
477 * Note, for multi-particle weights, the cut
478 * is loosend by a factor of 2
479 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
480 * the reduced @f$\chi^2@f$
481 * @param minWeight Least weight ot consider
482 * @param regCut Regularization cut-off
483 * @param residuals Mode for residual plots
485 * @return List of fit parameters
487 virtual TObjArray* Fit(TList* dir,
492 Double_t relErrorCut,
496 EResidualMethod residuals) const;
498 * Get the the 2D histogram @a name from our sub-list of @a
499 * dir. Then for each eta slice, try to fit the energu loss
500 * distribution up to @a nParticles particle responses.
502 * The fitted distributions (along with the functions fitted) are
503 * stored in a newly created sublist (<i>name</i>Dists).
505 * The fit parameters are also recorded in the newly created sub-list
506 * <i>name</i>Results.
508 * If @a residuals is not equal to kNoResiduals, then the
509 * residuals of the fits will be stored in the newly created
510 * sub-list <i>name</i>Residuals.
512 * A histogram named <i>name</i>Status is also generated and
513 * stored in the output list.
515 * @param dir Output list
516 * @param name Name of 2D base histogram in list
517 * @param lowCut Lower cut
518 * @param nParticles Max number of convolved landaus to fit
519 * @param minEntries Minimum number of entries
520 * @param minusBins Number of bins from peak to subtract to
522 * @param relErrorCut Cut applied to relative error of parameter.
523 * Note, for multi-particle weights, the cut
524 * is loosend by a factor of 2
525 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
526 * the reduced @f$\chi^2@f$
527 * @param minWeight Least weight ot consider
528 * @param regCut Regularization cut-off
529 * @param residuals Mode for residual plots
530 * @param scaleToPeak If true, scale distribution to peak value
531 * @param best Optional array to store fits in
533 * @return List of fit parameters
535 virtual TObjArray* FitSlices(TList* dir,
541 Double_t relErrorCut,
545 EResidualMethod residuals,
546 Bool_t scaleToPeak=true,
547 TObjArray* best=0) const;
550 * Fit a signal histogram. First, the bin @f$ b_{min}@f$ with
551 * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
552 * found. Then the fit range is set to the bin range
553 * @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
554 * particle signal is fitted to that. The parameters of that fit
555 * is then used as seeds for a fit of the @f$ N@f$ particle response
556 * to the data in the range
557 * @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
559 * @param dist Histogram to fit
560 * @param lowCut Lower cut @f$ E_{min}@f$ on signal
561 * @param minEntries Least number of entries required
562 * @param nParticles Max number @f$ N@f$ of convolved landaus to fit
563 * @param minusBins Number of bins @f$ \Delta b@f$ from peak to
564 * subtract to get the fit range
565 * @param relErrorCut Cut applied to relative error of parameter.
566 * Note, for multi-particle weights, the cut
567 * is loosend by a factor of 2
568 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
569 * the reduced @f$\chi^2@f$
570 * @param minWeight Least weight ot consider
571 * @param regCut Regularization cut-off
572 * @param scaleToPeak If true, scale distribution to peak value
573 * @param status On return, contain the status code (0: OK, 1:
574 * empty, 2: low statistics, 3: fit failed)
576 * @return The best fit function
578 virtual ELossFit_t* FitHist(TH1* dist,
583 Double_t relErrorCut,
588 UShort_t& status) const;
592 * @param dist Histogram
593 * @param relErrorCut Cut applied to relative error of parameter.
594 * Note, for multi-particle weights, the cut
595 * is loosend by a factor of 2
596 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
597 * the reduced @f$\chi^2@f$
598 * @param minWeightCut Least valid @f$ a_i@f$
602 virtual ELossFit_t* FindBestFit(const TH1* dist,
603 Double_t relErrorCut,
605 Double_t minWeightCut) const;
607 * Calculate residuals of the fit
609 * @param mode How to calculate
610 * @param lowCut Lower cut
611 * @param dist Distribution
612 * @param fit Function fitted to distribution
613 * @param out Output list to store residual histogram in
615 virtual void CalculateResiduals(EResidualMethod mode,
619 TCollection* out) const;
621 * Find the best fits. This assumes that the array fBest has been
622 * filled with the best possible fits for each eta bin, and that
623 * the fits are placed according to the bin number of the eta bin.
625 * This is called by the parent class when generating the corretion
628 * @param d Parent list
629 * @param obj Object to add fits to
630 * @param eta Eta axis
632 virtual void FindBestFits(const TList* d,
633 AliFMDCorrELossFit& obj,
636 * Make a parameter histogram
638 * @param name Name of histogram.
639 * @param title Title of histogram.
640 * @param eta Eta axis
644 TH1* MakePar(const char* name, const char* title, const TAxis& eta) const;
646 * Make a histogram that contains the results of the fit over the
651 * @param eta Eta axis
652 * @param low Least bin
653 * @param high Largest bin
654 * @param val Value of parameter
655 * @param err Error on parameter
657 * @return The newly allocated histogram
659 TH1* MakeTotal(const char* name,
666 TH1* fEDist; // Ring energy distribution
667 TH1* fEmpty; // Ring energy dist for empty events
668 TH2* fHist; // Two dimension Delta distribution
669 // TList* fEtaEDists; // Energy distributions per eta bin.
671 mutable TObjArray fBest;
672 mutable TClonesArray fFits;
674 ClassDef(RingHistos,4);
676 virtual RingHistos* CreateRingHistos(UShort_t d, Char_t r) const;
678 * Get the ring histogram container
683 * @return Ring histogram container
685 RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
687 * Check if the detector @a d, ring @a r is listed <i>in</i> the @a
688 * skips bit mask. If the detector/ring is in the mask, return true.
690 * That is, use case is
692 * for (UShort_t d=1. d<=3, d++) {
693 * UShort_t nr = (d == 1 ? 1 : 2);
694 * for (UShort_t q = 0; q < nr; q++) {
695 * Char_t r = (q == 0 ? 'I' : 'O');
696 * if (CheckSkips(d, r, skips)) continue;
697 * // Process detector/ring
704 * @param skips Mask of detector/rings to skip
706 * @return True if detector @a d, ring @a r is in the mask @a skips
708 static Bool_t CheckSkip(UShort_t d, Char_t r, UShort_t skips);
710 TList fRingHistos; // List of histogram containers
711 Double_t fLowCut; // Low cut on energy
712 UShort_t fNParticles; // Number of landaus to try to fit
713 UShort_t fMinEntries; // Minimum number of entries
714 UShort_t fFitRangeBinWidth; // N-bins to subtract from found max
715 Bool_t fDoFits; // Whether to actually do the fits
716 Bool_t fDoMakeObject; // Whether to make corrections object
717 TAxis fEtaAxis; // Eta axis
718 TAxis fCentralityAxis; // Centrality axis
719 Double_t fMaxE; // Maximum energy loss to consider
720 Int_t fNEbins; // Number of energy loss bins
721 Bool_t fUseIncreasingBins; // Wheter to use increasing bin sizes
722 Double_t fMaxRelParError; // Relative error cut
723 Double_t fMaxChi2PerNDF; // chi^2/nu cit
724 Double_t fMinWeight; // Minimum weight value
725 Int_t fDebug; // Debug level
726 EResidualMethod fResidualMethod; // Whether to store residuals (debugging)
727 UShort_t fSkips; // Rings to skip when fitting
728 Double_t fRegularizationCut; // When to regularize the chi^2
730 ClassDef(AliFMDEnergyFitter,6); //