1 // Object holding the Energy loss fit 'correction'
3 // These are generated from Monte-Carlo or real ESDs.
4 #ifndef ALIFMDCORRELOSSFIT_H
5 #define ALIFMDCORRELOSSFIT_H
7 * @file AliFMDCorrELossFit.h
8 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
9 * @date Wed Mar 23 14:01:15 2011
13 * @ingroup pwglf_forward_eloss
18 #include <TObjArray.h>
25 * @defgroup pwglf_forward_corr Corrections
27 * Correction objects used
29 * @ingroup pwglf_forward
32 * Object holding the Energy loss fit 'correction'
34 * These are generated from Monte-Carlo or real ESDs.
36 * @ingroup pwglf_forward_corr
37 * @ingroup pwglf_forward_eloss
39 class AliFMDCorrELossFit : public TObject
46 * POD structure to hold data from fits
48 * @ingroup pwglf_forward_corr
50 struct ELossFit : public TObject
52 Int_t fN; // Number of peaks fitted
53 UShort_t fNu; // Number degrees of freedom
54 Double_t fChi2; // Chi square from fit
55 Double_t fC; // Scaling constant
56 Double_t fDelta; // Most probable value
57 Double_t fXi; // Width parameter of Landau
58 Double_t fSigma; // Sigma on folded gaussian
59 Double_t fSigmaN; // Sigma of detector noise
60 Double_t* fA; // [fN] Weights
61 Double_t fEC; // Error on C
62 Double_t fEDelta; // Error on Delta
63 Double_t fEXi; // Error on Xi
64 Double_t fESigma; // Error on sigma
65 Double_t fESigmaN;// Error on sigma (noise)
66 Double_t* fEA; // [fN] Error on weights
67 Int_t fQuality;// Assigned quality
68 UShort_t fDet; // Detector
70 UShort_t fBin; // Eta bin
72 mutable UShort_t fMaxWeight; //!Cached maximum weight
74 static Double_t fgMaxRelError; // Global default max relative error
75 static Double_t fgLeastWeight; // Global default least weight
76 static Double_t fgMaxChi2nu; // Global default maximum reduced chi^2
83 * Construct from a function
85 * @param quality Quality flag
88 ELossFit(Int_t quality,const TF1& f);
90 * Constructor with full parameter set
92 * @param quality Quality flag
93 * @param n @f$ N@f$ - Number of fitted peaks
94 * @param chi2 @f$ \chi^2 @f$
95 * @param nu @f$ \nu @f$ - number degrees of freedom
96 * @param c @f$ C@f$ - scale constant
97 * @param ec @f$ \delta C@f$ - error on @f$ C@f$
98 * @param delta @f$ \Delta@f$ - Most probable value
99 * @param edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
100 * @param xi @f$ \xi@f$ - width
101 * @param exi @f$ \delta\xi@f$ - error on @f$\xi@f$
102 * @param sigma @f$ \sigma@f$ - Width of Gaussian
103 * @param esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
104 * @param sigman @f$ \sigma_n@f$ - Noise width
105 * @param esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
106 * @param a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
108 * @param ea Array of @f$ N-1@f$ error on the weights @f$ a_i@f$ for
111 ELossFit(Int_t quality,UShort_t n,
112 Double_t chi2, UShort_t nu,
113 Double_t c, Double_t ec,
114 Double_t delta, Double_t edelta,
115 Double_t xi, Double_t exi,
116 Double_t sigma, Double_t esigma,
117 Double_t sigman, Double_t esigman,
118 const Double_t* a,const Double_t* ea);
122 * @param o Object to copy from
124 ELossFit(const ELossFit& o);
126 * Assignment operator
128 * @param o Object to assign from
130 * @return Reference to this object
132 ELossFit& operator=(const ELossFit& o);
139 * @name Access to parameters
142 * @return Number of peaks fitted
144 Int_t GetN() const { return fN; }
146 * @return Number degrees of freedom
148 UShort_t GetNu() const { return fNu; }
150 * @return Chi square from fit
152 Double_t GetChi2() const { return fChi2; }
154 * @return Scaling constant
156 Double_t GetC() const { return fC; }
158 * @return Most probable value
160 Double_t GetDelta() const { return fDelta; }
162 * @return Width parameter of Landau
164 Double_t GetXi() const { return fXi; }
166 * @return Sigma on folded gaussian
168 Double_t GetSigma() const { return fSigma; }
170 * @return Sigma of detector noise
172 Double_t GetSigmaN() const { return fSigmaN; }
176 Double_t* GetAs() const { return fA; }
178 * @param i Which weight to get
182 Double_t GetA(UShort_t i) const;
186 Double_t GetEC() const { return fEC; }
188 * @return Error on Delta
190 Double_t GetEDelta() const { return fEDelta; }
192 * @return Error on Xi
194 Double_t GetEXi() const { return fEXi; }
196 * @return Error on sigma
198 Double_t GetESigma() const { return fESigma; }
200 * @return Error on sigma (noise)
202 Double_t GetESigmaN() const { return fESigmaN; }
204 * @return Error on weights
206 Double_t* GetEAs() const { return fEA; }
208 * @param i Which weight to get
210 * @return Error on weights
212 Double_t GetEA(UShort_t i) const;
214 * @return Assigned quality
216 Int_t GetQuality() const { return fQuality; }
220 UShort_t GetDet() const { return fDet; }
224 Char_t GetRing() const { return fRing; }
228 UShort_t GetBin() const { return fBin; }
238 * f_N(x;\Delta,\xi,\sigma') =
239 * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
242 * (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
243 * that fulfills the requirements
245 * @param x Where to evaluate
246 * @param maxN @f$ \max{N}@f$
248 * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
250 Double_t Evaluate(Double_t x,
251 UShort_t maxN=999) const;
255 * f_W(x;\Delta,\xi,\sigma') =
256 * \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma')}{
257 * f_N(x;\Delta,\xi,\sigma')} =
258 * \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i')}{
259 * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')}
261 * where @f$ n@f$ fulfills the requirements (see FindMaxWeight).
263 * If the denominator is zero, then 1 is returned.
265 * See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
266 * for more information on the evaluated functions.
268 * @param x Where to evaluate
269 * @param maxN @f$ \max{N}@f$
271 * @return @f$ f_W(x;\Delta,\xi,\sigma')@f$.
273 Double_t EvaluateWeighted(Double_t x,
274 UShort_t maxN=9999) const;
276 * Find the maximum weight to use. The maximum weight is the
277 * largest i for which
279 * - @f$ i \leq \max{N}@f$
280 * - @f$ a_i > \min{a}@f$
281 * - @f$ \delta a_i/a_i > \delta_{max}@f$
283 * @param maxRelError @f$ \min{a}@f$
284 * @param leastWeight @f$ \delta_{max}@f$
285 * @param maxN @f$ \max{N}@f$
287 * @return The largest index @f$ i@f$ for which the above
288 * conditions hold. Will never return less than 1.
290 Int_t FindMaxWeight(Double_t maxRelError=2*fgMaxRelError,
291 Double_t leastWeight=fgLeastWeight,
292 UShort_t maxN=999) const;
294 * Get a function that expresses this fit.
297 * f_N(x;\Delta,\xi,\sigma') =
298 * \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
300 * (see AliForwardUtil::NLandauGaus) or, if @a i is 1 or larger
302 * f_i(x;\Delta,\xi,\sigma') = a_i f(x;\Delta_i,\xi_i,\sigma_i')
304 * (see AliForwardUtil::ILandauGaus).
306 * @param i Component to get. If @a i is 0 or less, then the full
307 * function is returned, otherwise the specified component (if
309 * @param max Upper bound on function
311 * @return Pointer to newly allocated function. The caller owns
312 * this object, and must clean it up.
314 TF1* GetF1(Int_t i=0, Double_t max=20) const;
316 * Find the x value that corresponds to a (normalized) probability
317 * of @a low or less. That is, we can use this to say: "Give me
318 * the x value under which it is unlikely that a particle gave a
321 * @param low Threshold (between 0 and 1)
323 * @return Cut value, or 1000 in case of problems
325 Double_t FindProbabilityCut(Double_t low) const;
329 * @name TObject Sortable interface
332 * Declare this object as sortable
334 * @return Always true
336 Bool_t IsSortable() const { return kTRUE; }
338 * Compare to another ELossFit object.
340 * - +1, if this quality is better (larger) than other objects quality
341 * - -1, if this quality is worse (smaller) than other objects quality
342 * - +1, if this @f$|\chi^2/\nu-1|@f$ is smaller than the same for other
343 * - -1, if this @f$|\chi^2/\nu-1|@f$ is larger than the same for other
346 * @param o Other object to compare to
350 Int_t Compare(const TObject* o) const;
354 * name Auxiliary member functions
357 * Information to standard output
359 * @param option Not used
361 void Print(Option_t* option) const; // *MENU*
365 * @param option Options
366 * - COMP Draw components too
368 void Draw(Option_t* option="comp"); // *MENU*
374 void Browse(TBrowser* b);
376 * Get the name of this object
381 const Char_t* GetName() const;
383 * Calculate the lower bound
385 * @param f Width factor
386 * @param includeSigma Whether to include sigma
388 * @return @f$ \Delta - f (\xi + \sigma)@f$
390 Double_t GetLowerBound(Double_t f, Bool_t includeSigma) const;
392 * Calculate the lower bound
394 * @param f fraction of @f$\Delta@f$
396 * @return @f$ f\Delta@f$
398 Double_t GetLowerBound(Double_t f) const;
400 * Calculate the quality
402 * @param maxChi2nu Maximum reduced @f$\chi^2@f$
403 * @param maxRelError Maximum relative error
404 * @param leastWeight Least weight to use
406 void CalculateQuality(Double_t maxChi2nu=fgMaxChi2nu,
407 Double_t maxRelError=fgMaxRelError,
408 Double_t leastWeight=fgLeastWeight);
410 ClassDef(ELossFit,2); // Result of fit
414 * Default constructor
416 AliFMDCorrELossFit();
420 * @param o Object to copy from
422 AliFMDCorrELossFit(const AliFMDCorrELossFit& o);
426 virtual ~AliFMDCorrELossFit();
428 * Assignment operator
430 * @param o Object to assign from
432 * @return Reference to this object
434 AliFMDCorrELossFit& operator=(const AliFMDCorrELossFit& o);
441 * Set the fit parameters from a function
446 * @param quality Quality flag
447 * @param f Function from fit
449 * @return true on success
451 Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, Int_t quality,
454 * Set the fit parameters from a function
459 * @param f ELoss fit result - note, the object will take ownership
461 * @return true on success
463 Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* f);
465 * Set the fit parameters from a function
469 * @param etaBin Eta (bin number, 1->nBins)
470 * @param f ELoss fit result - note, the object will take ownership
472 * @return true on success
474 Bool_t SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* f);
476 * Set the fit parameters from a function
478 * @param d Detector number
479 * @param r Ring identifier
480 * @param eta Eta value
481 * @param quality Quality flag
482 * @param n @f$ N@f$ - Number of fitted peaks
483 * @param chi2 @f$ \chi^2 @f$
484 * @param nu @f$ \nu @f$ - number degrees of freedom
485 * @param c @f$ C@f$ - scale constant
486 * @param ec @f$ \delta C@f$ - error on @f$ C@f$
487 * @param delta @f$ \Delta@f$ - most probable value
488 * @param edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
489 * @param xi @f$ \xi@f$ - Landau width
490 * @param exi @f$ \delta\xi@f$ - error on @f$\xi@f$
491 * @param sigma @f$ \sigma@f$ - Gaussian width
492 * @param esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
493 * @param sigman @f$ \sigma_n@f$ - Noise width
494 * @param esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
495 * @param a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
497 * @param ea Array of @f$ N-1@f$ errors on weights @f$ a_i@f$ for
500 * @return true on success
502 Bool_t SetFit(UShort_t d, Char_t r, Double_t eta,
503 Int_t quality,UShort_t n,
504 Double_t chi2, UShort_t nu,
505 Double_t c, Double_t ec,
506 Double_t delta, Double_t edelta,
507 Double_t xi, Double_t exi,
508 Double_t sigma, Double_t esigma,
509 Double_t sigman, Double_t esigman,
510 Double_t* a, Double_t* ea);
515 * @name Set and get eta axis
518 * Set the eta axis to use
520 * @param axis Eta axis
522 void SetEtaAxis(const TAxis& axis);
524 * Set the eta axis to use
526 * @param nBins Number of bins
527 * @param min Minimum @f$ \eta@f$
528 * @param max maximum @f$ \eta@f$
530 void SetEtaAxis(Int_t nBins, Double_t min, Double_t max);
532 * Get the eta axis used
536 const TAxis& GetEtaAxis() const { return fEtaAxis; }
538 * Set the low cut used when fitting
540 * @param cut Cut value
542 void SetLowCut(Double_t cut) { fLowCut = cut; }
544 * Get the low cut used when fitting
546 * @return Low cut used for fitting
548 Double_t GetLowCut() const { return fLowCut; }
550 * Find the eta bin corresponding to the given eta
552 * @param eta Eta value
554 * @return Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
555 * eta, or 0 if out of range.
557 Int_t FindEtaBin(Double_t eta) const;
565 * Find the fit corresponding to the specified parameters. This
566 * uses the cache map of good fits for the look-up. For un-cached
567 * look-up see GetFit.
571 * @param eta Eta value
572 * @param minQ Minimum quality
574 * @return Fit parameters or null in case of problems
576 ELossFit* FindFit(UShort_t d, Char_t r, Double_t eta,
577 UShort_t minQ) const;
579 * Find the fit corresponding to the specified parameters. This
580 * uses the cache map of good fits for the look-up. For un-cached
581 * look-up see GetFit.
585 * @param etabin Eta bin (1 based)
586 * @param minQ Minimum quality
588 * @return Fit parameters or null in case of problems
590 ELossFit* FindFit(UShort_t d, Char_t r, Int_t etabin,
591 UShort_t minQ) const;
593 * Find the fit corresponding to the specified parameters. Note,
594 * the a cache-map of good fits isn't used for this look-up. To use
595 * the cache, use FindFit.
599 * @param eta Eta value
601 * @return Fit parameters or null in case of problems
603 ELossFit* GetFit(UShort_t d, Char_t r, Double_t eta) const;
605 * Find the fit corresponding to the specified parameters. Note,
606 * the a cache-map of good fits isn't used for this look-up. To use
607 * the cache, use FindFit.
611 * @param etabin Eta bin (1 based)
613 * @return Fit parameters or null in case of problems
615 ELossFit* GetFit(UShort_t d, Char_t r, Int_t etabin) const;
620 * @name Lower bounds on fits
623 * Get the lower validity bound of the fit.
627 * @param etaBin Eta bin (1-based)
628 * @param f Fraction of @f$\Delta_{mp}@f$
630 * @return @f$ f\Delta_{mp}@f$
632 Double_t GetLowerBound(UShort_t d, Char_t r, Int_t etaBin,
635 * Get the lower validity bound of the fit.
639 * @param eta Eta value
640 * @param f Fraction of @f$\Delta_{mp}@f$
642 * @return @f$ f\Delta_{mp}@f$
644 Double_t GetLowerBound(UShort_t d, Char_t r, Double_t eta,
647 * Get the lower validity bound of the fit.
651 * @param etaBin Eta bin (1-based)
652 * @param p Probability cut
653 * @param dummy Not used
655 * @return @f$ x@f$ for which @f$ P(x>p)@f$
657 Double_t GetLowerBound(UShort_t d, Char_t r, Int_t etaBin,
658 Double_t p, Bool_t dummy) const;
660 * Get the lower validity bound of the fit.
664 * @param eta Eta value
665 * @param p Probability cut
666 * @param dummy Not used
668 * @return @f$ x@f$ for which @f$ P(x>p)@f$
670 Double_t GetLowerBound(UShort_t d, Char_t r, Double_t eta,
671 Double_t p, Bool_t dummy) const;
673 * Get the lower validity bound of the fit
677 * @param etaBin Eta bin (1-based)
678 * @param f Factor on xi (and sigma)
679 * @param showErrors Show errors
680 * @param includeSigma Whether to include sigma
682 * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$
684 Double_t GetLowerBound(UShort_t d, Char_t r, Int_t etaBin,
685 Double_t f, Bool_t showErrors,
686 Bool_t includeSigma) const;
688 * Get the lower validity bound of the fit
692 * @param eta Eta value
693 * @param f Factor on xi (and sigma)
694 * @param showErrors Show errors
695 * @param includeSigma Whether to include sigma
697 * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$
699 Double_t GetLowerBound(UShort_t d, Char_t r, Double_t eta,
700 Double_t f, Bool_t showErrors,
701 Bool_t includeSigma) const;
706 * @name Miscellaneous
708 void CacheBins(UShort_t minQuality) const;
710 * Get the ring array corresponding to the specified ring
715 * @return Pointer to ring array, or null in case of problems
717 TObjArray* GetRingArray(UShort_t d, Char_t r) const;
719 * Signal that this is a folder
721 * @return Always true
723 Bool_t IsFolder() const { return true; }
729 void Browse(TBrowser* b);
733 * @param option Options. Possible values are
734 * - error Plot error bars
735 * - relative Plot relative errors
737 void Draw(Option_t* option=""); //*MENU*
741 * @param option Options
742 * - R Print recursive
745 void Print(Option_t* option="R") const; //*MENU*
747 * Get a list of THStack - one for each parameter.
749 * If @a err is true, then error bars are set too. If @a rel is
750 * true, then the relative error (rather than the absolute value) is
751 * filled into the histograms. If @a good is true, then we use the
752 * cache-map of good fits rather than all fits.
754 * @param err Show errors
755 * @param rel Show relative errors
756 * @param good Only show good fits
757 * @param maxN Maximum weight to use
759 * @return List of THStack
761 TList* GetStacks(Bool_t err, Bool_t rel, Bool_t good, UShort_t maxN=5) const;
765 * Get the ring array corresponding to the specified ring
770 * @return Pointer to ring array, or newly created container
772 TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r);
774 * Fill fit values into stack histograms
777 * @param rel If true, fill in relative errors
778 * @param used The bin begin used for this fit
779 * @param hChi @f$\chi^2/\nu@f$ histogram
780 * @param hN @f$ N_{a}@f$ - number of components - histogram
781 * @param hC @f$ C@f$ - prefactor - histogram
782 * @param hDelta @f$ \Delta_p@f$ - most-probably value - histogram
783 * @param hXi @f$ \xi@f$ - Landau 'width' - histogram
784 * @param hSigma @f$ \sigma@f$ - Gaussian smear - histogram
785 * @param maxN @f$ N_{a,max}@f$ Largest possible @f$ N@f$
786 * @param hA @f$ a_{i}, i=\{2,..,N_{a,max}\}@f$ - histogram
788 void UpdateStackHist(ELossFit* f, Bool_t rel,
791 TH1* hC, TH1* hDelta,
792 TH1* hXi, TH1* hSigma,
793 Int_t maxN, TH1** hA) const;
795 TObjArray fRings; // Array of rings
796 TAxis fEtaAxis; // Eta axis used
797 Double_t fLowCut; // Low cut used when fitting
798 mutable TArrayI fCache;
800 ClassDef(AliFMDCorrELossFit,3);
803 //____________________________________________________________________
805 AliFMDCorrELossFit::SetEtaAxis(Int_t nBins, Double_t min, Double_t max)
807 fEtaAxis.Set(nBins, min, max);
809 //____________________________________________________________________
811 AliFMDCorrELossFit::SetEtaAxis(const TAxis& e)
813 fEtaAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
815 //____________________________________________________________________
817 AliFMDCorrELossFit::ELossFit::GetA(UShort_t i) const
820 if (i > fN) return 0;
821 if (i == 1) return 1;
824 //____________________________________________________________________
826 AliFMDCorrELossFit::ELossFit::GetEA(UShort_t i) const
829 if (i > fN) return 0;
830 if (i == 1) return 1;