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"
24 #include "AliLandauGaus.h"
33 * Class to fit the energy distribution.
36 * - AliESDFMD object - from reconstruction
39 * - Lists of histogram - one per ring. Each list has a number of
40 * histograms corresponding to the number of eta bins defined.
42 * @par Corrections used:
45 * @image html alice-int-2012-040-eloss_fits.png "Summary of fits"
47 * @ingroup pwglf_forward_algo
48 * @ingroup pwglf_forward_eloss
50 class AliFMDEnergyFitter : public TNamed
54 * Enumeration of parameters
57 /** Index of pre-constant @f$ C@f$ */
58 kC = AliLandauGaus::kC,
59 /** Index of most probable value @f$ \Delta_p@f$ */
60 kDelta = AliLandauGaus::kDelta,
61 /** Index of Landau width @f$ \xi@f$ */
62 kXi = AliLandauGaus::kXi,
63 /** Index of Gaussian width @f$ \sigma@f$ */
64 kSigma = AliLandauGaus::kSigma,
65 /** Index of Gaussian additional width @f$ \sigma_n@f$ */
66 kSigmaN = AliLandauGaus::kSigmaN,
67 /** Index of Number of particles @f$ N@f$ */
68 kN = AliLandauGaus::kN,
69 /** Base index of particle strengths @f$ a_i@f$ for
71 kA = AliLandauGaus::kA
74 * Enumeration of residual methods
76 enum EResidualMethod {
77 /** Do not calculate residuals */
79 /** The residuals stored are the difference, and the errors are
80 stored in the error bars of the histogram. */
82 /** The residuals stored are the differences scaled to the error
84 kResidualScaledDifference,
85 /** The residuals stored are the square difference scale to the
86 square error on the data. */
87 kResidualSquareDifference
91 * FMD ring bits for skipping
103 kFMD2 =kFMD2I|kFMD2O,
114 virtual ~AliFMDEnergyFitter();
116 * Default Constructor - do not use
118 AliFMDEnergyFitter();
122 * @param title Title of object - not significant
124 AliFMDEnergyFitter(const char* title);
126 // -----------------------------------------------------------------
129 * @name Setters of options and parameters
132 * Set the eta axis to use. This will force the code to use this
133 * eta axis definition - irrespective of whatever axis is passed to
134 * the Init member function. Therefore, this member function can be
135 * used to force another eta axis than one found in the correction
138 * @param nBins Number of bins
139 * @param etaMin Minimum of the eta axis
140 * @param etaMax Maximum of the eta axis
142 void SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax);
144 * Set the eta axis to use. This will force the code to use this
145 * eta axis definition - irrespective of whatever axis is passed to
146 * the Init member function. Therefore, this member function can be
147 * used to force another eta axis than one found in the correction
150 * @param etaAxis Eta axis to use
152 void SetEtaAxis(const TAxis& etaAxis);
154 * Set the centrality bins. E.g.,
157 * Double_t bins[] = { 0., 5., 10., 15., 20., 30.,
158 * 40., 50., 60., 70., 80., 100. };
159 * task->GetFitter().SetCentralityBins(n, bins);
162 * @param nBins Size of @a bins
163 * @param bins Bin limits.
165 void SetCentralityAxis(UShort_t nBins, Double_t* bins);
167 * Set the low cut used for energy
169 * @param lowCut Low cut
171 void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
173 * Set the number of bins to subtract
177 void SetFitRangeBinWidth(UShort_t n=4) { fFitRangeBinWidth = n; }
179 * Whether or not to enable fitting of the final merged result.
180 * Note, fitting takes quite a while and one should be careful not to do
183 * @param doFit Whether to do the fits or not
185 void SetDoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
187 * Set whether to make the corrections object on the output. Note,
188 * fits should be enable for this to have any effect.
190 * @param doMake If true (false is default), do make the corrections object.
192 void SetDoMakeObject(Bool_t doMake=kTRUE) { fDoMakeObject = doMake; }
194 * Set how many particles we will try to fit at most to the data
196 * @param n Max number of particle to try to fit
198 void SetNParticles(UShort_t n) { fNParticles = (n<1 ? 1 : (n>5 ? 5 : n)); }
200 * Set the minimum number of entries each histogram must have
201 * before we try to fit our response function to it
203 * @param n Minimum number of entries
205 void SetMinEntries(UShort_t n) { fMinEntries = (n < 1 ? 1 : n); }
207 * Set maximum energy loss to consider
209 * @param x Maximum energy loss to consider
211 void SetMaxE(Double_t x) { fMaxE = x; }
213 * Set number of energy loss bins
215 * @param x Number of energy loss bins
217 void SetNEbins(Int_t x) { fNEbins = x; }
219 * Set the maximum relative error
221 * @param e Maximum relative error
223 void SetMaxRelativeParameterError(Double_t e=0.2) { fMaxRelParError = e; }
225 * Set the maximum @f$ \chi^2/\nu@f$
227 * @param c Maximum @f$ \chi^2/\nu@f$
229 void SetMaxChi2PerNDF(Double_t c=10) { fMaxChi2PerNDF = c; }
231 * Set the least weight
233 * @param c Least weight
235 void SetMinWeight(Double_t c=1e-7) { fMinWeight = c; }
237 * Set wheter to use increasing bin sizes
239 * @param x Wheter to use increasing bin sizes
241 void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
243 * Set whether to make residuals, and in that case how.
245 * - Square difference: @f$chi_i^2=(h_i - f(x_i))^2/\delta_i^2@f$
246 * - Scaled difference: @f$(h_i - f(x_i))/\delta_i@f$
247 * - Difference: @f$(h_i - f(x_i)) \pm\delta_i@f$
249 * where @f$h_i, x_i, \delta_i@f$ is the bin content, bin center,
250 * and bin error for bin @f$i@f$ respectively, and @f$ f@f$ is the
253 * @param x Residual method
255 void SetStoreResiduals(EResidualMethod x=kResidualDifference)
260 * Set the regularization cut @f$c_{R}@f$. If a @f$\Delta@f$
261 * distribution has more entries @f$ N_{dist}@f$ than @f$c_{R}@f$,
262 * then we modify the errors of the the distribution with the factor
265 * \sqrt{N_{dist}/c_{R}}
268 * to keep the @f$\chi^2/\nu@f$ within resonable limits.
270 * The large residuals @f$chi_i^2=(h_i - f(x_i))^2/\delta_i^2@f$
271 * (see also SetStoreResiduals) comes about on the boundary between
272 * the @f$N@f$ and @f$N+1@f$ particle contributions, and seems to
273 * fall off for larger @f$N@f$. This may indicate that there's a
274 * component in the distributions that the function
277 * f(\Delta;\Delta_p,\xi,\sigma,\mathbf{a}) = \sum_i=1^{n} a_i\int
278 * d\Delta' L(\Delta;\Delta',\xi) G(\Delta';\Delta_p,\sigma)
285 void SetRegularizationCut(Double_t cut=3e6)
287 fRegularizationCut = cut;
289 void SetSkips(UShort_t skip) { fSkips = skip; }
291 * Set the debug level. The higher the value the more output
293 * @param dbg Debug level
295 void SetDebug(Int_t dbg=1);
297 * Whether to enable the extra shift in the MPV from @f$ \sigma/\xi@f$
299 * @param use If true, enable extra shift @f$\delta\Delta_p(\sigma/\xi)@f$
301 void SetEnableDeltaShift(Bool_t use=true);
304 // -----------------------------------------------------------------
311 * Define the output histograms. These are put in a sub list of the
312 * passed list. The histograms are merged before the parent task calls
313 * AliAnalysisTaskSE::Terminate
315 * @param dir Directory to add to
317 virtual void CreateOutputObjects(TList* dir);
319 * Initialise the task
321 * @param etaAxis The eta axis to use. Note, that if the eta axis
322 * has already been set (using SetEtaAxis), then this parameter will be
325 virtual void SetupForData(const TAxis& etaAxis);
327 * Fitter the input AliESDFMD object
330 * @param cent Event centrality (or < 0 if not valid)
331 * @param empty Whether the event is 'empty'
333 * @return True on success, false otherwise
335 virtual Bool_t Accumulate(const AliESDFMD& input,
339 * Scale the histograms to the total number of events
341 * @param dir Where the histograms are
343 virtual void Fit(const TList* dir);
345 * Generate the corrections object
347 * @param dir List to analyse
349 void MakeCorrectionsObject(TList* dir);
354 * @param option Not used
356 void Print(Option_t* option="") const;
358 * Read the parameters from a list - used when re-running the code
360 * @param list Input list
362 * @return true if the parameter where read
364 Bool_t ReadParameters(const TCollection* list);
369 * @param o Object to copy from
371 AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
373 * Assignment operator
375 * @param o Object to assign from
377 * @return Reference to this
379 AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
382 * Internal data structure to keep track of the histograms
384 struct RingHistos : public AliForwardUtil::RingHistos
386 typedef AliFMDCorrELossFit::ELossFit ELossFit_t;
397 RingHistos(UShort_t d, Char_t r);
399 * Copy constructor - not defined
401 * @param o Object to copy from
403 RingHistos(const RingHistos& o);
405 * Assignment operator - not defined
407 * @param o Object to assign from
409 * @return Reference to this
411 RingHistos& operator=(const RingHistos& o);
417 * Make an axis with increasing bins
419 * @param n Number of bins
423 * @return An axis with quadratically increasing bin size
425 virtual TArrayD MakeIncreasingAxis(Int_t n,
429 * Make E/E_mip histogram
431 * @param name Name of histogram
432 * @param title Title of histogram
433 * @param eAxis @f$\eta@f$ axis
434 * @param deMax Maximum energy loss
435 * @param nDeBins Number energy loss bins
436 * @param incr Whether to make bins of increasing size
438 TH2* Make(const char* name,
449 virtual void CreateOutputObjects(TList* dir);
453 * @param eAxis Eta axis
454 * @param cAxis Centrality axis
455 * @param maxDE Max energy loss to consider
456 * @param nDEbins Number of bins
457 * @param useIncrBin Whether to use an increasing bin size
459 virtual void SetupForData(const TAxis& eAxis,
463 Bool_t useIncrBin=true);
467 * @param empty True if event is empty
468 * @param eta @f$ Eta@f$
469 * @param icent Centrality bin (1 based)
472 virtual void Fill(Bool_t empty, Double_t eta, Int_t icent, Double_t mult);
474 * Get the the 2D histogram eloss name from our sub-list of @a dir
475 * and call the Fit function described below (with &fBest) as last
478 * @param dir Output list
479 * @param lowCut Lower cut
480 * @param nParticles Max number of convolved landaus to fit
481 * @param minEntries Minimum number of entries
482 * @param minusBins Number of bins from peak to subtract to
484 * @param relErrorCut Cut applied to relative error of parameter.
485 * Note, for multi-particle weights, the cut
486 * is loosend by a factor of 2
487 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
488 * the reduced @f$\chi^2@f$
489 * @param minWeight Least weight ot consider
490 * @param regCut Regularization cut-off
491 * @param residuals Mode for residual plots
493 * @return List of fit parameters
495 virtual TObjArray* Fit(TList* dir,
500 Double_t relErrorCut,
504 EResidualMethod residuals) const;
506 * Get the the 2D histogram @a name from our sub-list of @a
507 * dir. Then for each eta slice, try to fit the energu loss
508 * distribution up to @a nParticles particle responses.
510 * The fitted distributions (along with the functions fitted) are
511 * stored in a newly created sublist (<i>name</i>Dists).
513 * The fit parameters are also recorded in the newly created sub-list
514 * <i>name</i>Results.
516 * If @a residuals is not equal to kNoResiduals, then the
517 * residuals of the fits will be stored in the newly created
518 * sub-list <i>name</i>Residuals.
520 * A histogram named <i>name</i>Status is also generated and
521 * stored in the output list.
523 * @param dir Output list
524 * @param name Name of 2D base histogram in list
525 * @param lowCut Lower cut
526 * @param nParticles Max number of convolved landaus to fit
527 * @param minEntries Minimum number of entries
528 * @param minusBins Number of bins from peak to subtract to
530 * @param relErrorCut Cut applied to relative error of parameter.
531 * Note, for multi-particle weights, the cut
532 * is loosend by a factor of 2
533 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
534 * the reduced @f$\chi^2@f$
535 * @param minWeight Least weight ot consider
536 * @param regCut Regularization cut-off
537 * @param residuals Mode for residual plots
538 * @param scaleToPeak If true, scale distribution to peak value
539 * @param best Optional array to store fits in
541 * @return List of fit parameters
543 virtual TObjArray* FitSlices(TList* dir,
549 Double_t relErrorCut,
553 EResidualMethod residuals,
554 Bool_t scaleToPeak=true,
555 TObjArray* best=0) const;
558 * Fit a signal histogram. First, the bin @f$ b_{min}@f$ with
559 * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
560 * found. Then the fit range is set to the bin range
561 * @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
562 * particle signal is fitted to that. The parameters of that fit
563 * is then used as seeds for a fit of the @f$ N@f$ particle response
564 * to the data in the range
565 * @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
567 * @param dist Histogram to fit
568 * @param lowCut Lower cut @f$ E_{min}@f$ on signal
569 * @param minEntries Least number of entries required
570 * @param nParticles Max number @f$ N@f$ of convolved landaus to fit
571 * @param minusBins Number of bins @f$ \Delta b@f$ from peak to
572 * subtract to get the fit range
573 * @param relErrorCut Cut applied to relative error of parameter.
574 * Note, for multi-particle weights, the cut
575 * is loosend by a factor of 2
576 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
577 * the reduced @f$\chi^2@f$
578 * @param minWeight Least weight ot consider
579 * @param regCut Regularization cut-off
580 * @param scaleToPeak If true, scale distribution to peak value
581 * @param status On return, contain the status code (0: OK, 1:
582 * empty, 2: low statistics, 3: fit failed)
584 * @return The best fit function
586 virtual ELossFit_t* FitHist(TH1* dist,
591 Double_t relErrorCut,
596 UShort_t& status) const;
600 * @param dist Histogram
601 * @param relErrorCut Cut applied to relative error of parameter.
602 * Note, for multi-particle weights, the cut
603 * is loosend by a factor of 2
604 * @param chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
605 * the reduced @f$\chi^2@f$
606 * @param minWeightCut Least valid @f$ a_i@f$
610 virtual ELossFit_t* FindBestFit(const TH1* dist,
611 Double_t relErrorCut,
613 Double_t minWeightCut) const;
615 * Calculate residuals of the fit
617 * @param mode How to calculate
618 * @param lowCut Lower cut
619 * @param dist Distribution
620 * @param fit Function fitted to distribution
621 * @param out Output list to store residual histogram in
623 virtual void CalculateResiduals(EResidualMethod mode,
627 TCollection* out) const;
629 * Find the best fits. This assumes that the array fBest has been
630 * filled with the best possible fits for each eta bin, and that
631 * the fits are placed according to the bin number of the eta bin.
633 * This is called by the parent class when generating the corretion
636 * @param d Parent list
637 * @param obj Object to add fits to
638 * @param eta Eta axis
640 virtual void FindBestFits(const TList* d,
641 AliFMDCorrELossFit& obj,
644 * Make a parameter histogram
646 * @param name Name of histogram.
647 * @param title Title of histogram.
648 * @param eta Eta axis
652 TH1* MakePar(const char* name, const char* title, const TAxis& eta) const;
654 * Make a histogram that contains the results of the fit over the
659 * @param eta Eta axis
660 * @param low Least bin
661 * @param high Largest bin
662 * @param val Value of parameter
663 * @param err Error on parameter
665 * @return The newly allocated histogram
667 TH1* MakeTotal(const char* name,
674 TH1* fEDist; // Ring energy distribution
675 TH1* fEmpty; // Ring energy dist for empty events
676 TH2* fHist; // Two dimension Delta distribution
677 // TList* fEtaEDists; // Energy distributions per eta bin.
679 mutable TObjArray fBest;
680 mutable TClonesArray fFits;
682 ClassDef(RingHistos,4);
684 virtual RingHistos* CreateRingHistos(UShort_t d, Char_t r) const;
686 * Get the ring histogram container
691 * @return Ring histogram container
693 RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
695 * Check if the detector @a d, ring @a r is listed <i>in</i> the @a
696 * skips bit mask. If the detector/ring is in the mask, return true.
698 * That is, use case is
700 * for (UShort_t d=1. d<=3, d++) {
701 * UShort_t nr = (d == 1 ? 1 : 2);
702 * for (UShort_t q = 0; q < nr; q++) {
703 * Char_t r = (q == 0 ? 'I' : 'O');
704 * if (CheckSkips(d, r, skips)) continue;
705 * // Process detector/ring
712 * @param skips Mask of detector/rings to skip
714 * @return True if detector @a d, ring @a r is in the mask @a skips
716 static Bool_t CheckSkip(UShort_t d, Char_t r, UShort_t skips);
718 TList fRingHistos; // List of histogram containers
719 Double_t fLowCut; // Low cut on energy
720 UShort_t fNParticles; // Number of landaus to try to fit
721 UShort_t fMinEntries; // Minimum number of entries
722 UShort_t fFitRangeBinWidth; // N-bins to subtract from found max
723 Bool_t fDoFits; // Whether to actually do the fits
724 Bool_t fDoMakeObject; // Whether to make corrections object
725 TAxis fEtaAxis; // Eta axis
726 TAxis fCentralityAxis; // Centrality axis
727 Double_t fMaxE; // Maximum energy loss to consider
728 Int_t fNEbins; // Number of energy loss bins
729 Bool_t fUseIncreasingBins; // Wheter to use increasing bin sizes
730 Double_t fMaxRelParError; // Relative error cut
731 Double_t fMaxChi2PerNDF; // chi^2/nu cit
732 Double_t fMinWeight; // Minimum weight value
733 Int_t fDebug; // Debug level
734 EResidualMethod fResidualMethod; // Whether to store residuals (debugging)
735 UShort_t fSkips; // Rings to skip when fitting
736 Double_t fRegularizationCut; // When to regularize the chi^2
738 ClassDef(AliFMDEnergyFitter,8); //