PWG2/FORWARD/analysis2/AliFMDCorrELossFit.h

   1 // Object holding the Energy loss fit 'correction'
   2 //
   3 // These are generated from Monte-Carlo or real ESDs.
   4 #ifndef ALIFMDCORRELOSSFIT_H
   5 #define ALIFMDCORRELOSSFIT_H
   6 /**
   7  * @file   AliFMDCorrELossFit.h
   8  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
   9  * @date   Wed Mar 23 14:01:15 2011
  10  *
  11  * @brief
  12  *
  13  * @ingroup pwg2_forward_eloss
  14  *
  15  */
  16 #include <TObject.h>
  17 #include <TAxis.h>
  18 #include <TObjArray.h>
  19 class TF1;
  20 class TBrowser;
  21
  22 /**
  23  * @defgroup pwg2_forward_corr Corrections
  24  *
  25  * @ingroup pwg2_forward
  26  */
  27 /**
  28  * Object holding the Energy loss fit 'correction'
  29  *
  30  * These are generated from Monte-Carlo or real ESDs.
  31  *
  32  * @ingroup pwg2_forward_corr
  33  * @ingroup pwg2_forward_eloss
  34  */
  35 class AliFMDCorrELossFit : public TObject
  36 {
  37 public:
  38   /**
  39    * POD structure to hold data from fits
  40    *
  41    * @ingroup pwg2_forward_corr
  42    */
  43   struct ELossFit : public TObject
  44   {
  45     Int_t     fN;      // Number of peaks fitted
  46     UShort_t  fNu;     // Number degrees of freedom
  47     Double_t  fChi2;   // Chi square from fit
  48     Double_t  fC;      // Scaling constant
  49     Double_t  fDelta;  // Most probable value
  50     Double_t  fXi;     // Width parameter of Landau
  51     Double_t  fSigma;  // Sigma on folded gaussian
  52     Double_t  fSigmaN; // Sigma of detector noise
  53     Double_t* fA;      // [fN] Weights
  54     Double_t  fEC;     // Error on C
  55     Double_t  fEDelta; // Error on Delta
  56     Double_t  fEXi;    // Error on Xi
  57     Double_t  fESigma; // Error on sigma
  58     Double_t  fESigmaN;// Error on sigma (noise)
  59     Double_t* fEA;     // [fN] Error on weights
  60     Int_t     fQuality;// Assigned quality
  61     UShort_t  fDet;    // Detector
  62     Char_t    fRing;   // Ring
  63     UShort_t  fBin;    // Eta bin
  64
  65     static Double_t fgMaxRelError;  // Global default max relative error
  66     static Double_t fgLeastWeight;  // Global default least weight
  67     static Double_t fgMaxChi2nu;    // Global default maximum reduced chi^2
  68     /**
  69      * Default constructor
  70      *
  71      */
  72     ELossFit();
  73     /**
  74      * Construct from a function
  75      *
  76      * @param quality Quality flag
  77      * @param f       Function
  78      */
  79     ELossFit(Int_t quality,const TF1& f);
  80     /**
  81      * Constructor with full parameter set
  82      *
  83      * @param quality   Quality flag
  84      * @param n         @f$ N@f$ - Number of fitted peaks
  85      * @param chi2      @f$ \chi^2 @f$
  86      * @param nu        @f$ \nu @f$ - number degrees of freedom
  87      * @param c         @f$ C@f$ - scale constant
  88      * @param ec        @f$ \delta C@f$ - error on @f$ C@f$
  89      * @param delta     @f$ \Delta@f$ - Most probable value
  90      * @param edelta    @f$ \delta\Delta@f$ - error on @f$\Delta@f$
  91      * @param xi        @f$ \xi@f$ - width
  92      * @param exi       @f$ \delta\xi@f$ - error on @f$\xi@f$
  93      * @param sigma     @f$ \sigma@f$ - Width of Gaussian
  94      * @param esigma    @f$ \delta\sigma@f$ - error on @f$\sigma@f$
  95      * @param sigman    @f$ \sigma_n@f$ - Noise width
  96      * @param esigman   @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
  97      * @param a         Array of @f$ N-1@f$ weights @f$ a_i@f$ for
  98      *                  @f$ i=2,\ldots@f$
  99      * @param ea        Array of @f$ N-1@f$ error on the weights @f$ a_i@f$ for
 100      *                  @f$ i=2,\ldots@f$
 101      */
 102     ELossFit(Int_t     quality,UShort_t  n,
 103              Double_t  chi2,   UShort_t  nu,
 104              Double_t  c,      Double_t  ec,
 105              Double_t  delta,  Double_t  edelta,
 106              Double_t  xi,     Double_t  exi,
 107              Double_t  sigma,  Double_t  esigma,
 108              Double_t  sigman, Double_t  esigman,
 109              const Double_t* a,const Double_t* ea);
 110     /**
 111      * Copy constructor
 112      *
 113      * @param o Object to copy from
 114      */
 115     ELossFit(const ELossFit& o);
 116     /**
 117      * Assignment operator
 118      *
 119      * @param o Object to assign from
 120      *
 121      * @return Reference to this object
 122      */
 123     ELossFit& operator=(const ELossFit& o);
 124     /**
 125      * Destructor
 126      */
 127     ~ELossFit();
 128     /**
 129      * @{
 130      * @name Access to parameters
 131      */
 132     /**
 133      * @return Number of peaks fitted
 134      */
 135     Int_t GetN() const { return fN; }
 136     /**
 137      * @return Number degrees of freedom
 138      */
 139     UShort_t GetNu() const { return fNu; }
 140     /**
 141      * @return Chi square from fit
 142      */
 143     Double_t GetChi2() const { return fChi2; }
 144     /**
 145      * @return Scaling constant
 146      */
 147     Double_t GetC() const { return fC; }
 148     /**
 149      * @return Most probable value
 150      */
 151     Double_t GetDelta() const { return fDelta; }
 152     /**
 153      * @return Width parameter of Landau
 154      */
 155     Double_t GetXi() const { return fXi; }
 156     /**
 157      * @return Sigma on folded gaussian
 158      */
 159     Double_t GetSigma() const { return fSigma; }
 160     /**
 161      * @return Sigma of detector noise
 162      */
 163     Double_t GetSigmaN() const { return fSigmaN; }
 164     /**
 165      * @return Weights
 166      */
 167     Double_t* GetAs() const { return fA; }
 168     /**
 169      * @return Weights
 170      */
 171     Double_t GetA(UShort_t i) const;
 172     /**
 173      * @return Error on C
 174      */
 175     Double_t GetEC() const { return fEC; }
 176     /**
 177      * @return Error on Delta
 178      */
 179     Double_t GetEDelta() const { return fEDelta; }
 180     /**
 181      * @return Error on Xi
 182      */
 183     Double_t GetEXi() const { return fEXi; }
 184     /**
 185      * @return Error on sigma
 186      */
 187     Double_t GetESigma() const { return fESigma; }
 188     /**
 189      * @return Error on sigma (noise)
 190      */
 191     Double_t GetESigmaN() const { return fESigmaN; }
 192     /**
 193      * @return Error on weights
 194      */
 195     Double_t* GetEAs() const { return fEA; }
 196     /**
 197      * @return Error on weights
 198      */
 199     Double_t GetEA(UShort_t i) const;
 200     /**
 201      * @return Assigned quality
 202      */
 203     Int_t GetQuality() const { return fQuality; }
 204     /**
 205      * @return Detector
 206      */
 207     UShort_t GetDet() const { return fDet; }
 208     /**
 209      * @return Ring
 210      */
 211     Char_t GetRing() const { return fRing; }
 212     /**
 213      * @return Eta bin
 214      */
 215     UShort_t GetBin() const { return fBin; }
 216     /* @} */
 217
 218     /**
 219      * @{
 220      * @name Evaluation
 221      */
 222     /**
 223      * Evaluate
 224      * @f[
 225      *  f_N(x;\Delta,\xi,\sigma') =
 226      *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
 227      * @f]
 228      *
 229      * (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
 230      * that fulfills the requirements
 231      *
 232      * @param x           Where to evaluate
 233      * @param maxN        @f$ \max{N}@f$
 234      *
 235      * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
 236      */
 237     Double_t Evaluate(Double_t x,
 238                       UShort_t maxN=999) const;
 239     /**
 240      * Evaluate
 241      * @f[
 242      *   f_W(x;\Delta,\xi,\sigma') =
 243      *   \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma')}{
 244      *     f_N(x;\Delta,\xi,\sigma')} =
 245      *   \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i')}{
 246      *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')}
 247      * @f]
 248      * where @f$ n@f$ fulfills the requirements (see FindMaxWeight).
 249      *
 250      * If the denominator is zero, then 1 is returned.
 251      *
 252      * See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
 253      * for more information on the evaluated functions.
 254      *
 255      * @param x           Where to evaluate
 256      * @param maxN        @f$ \max{N}@f$
 257      *
 258      * @return @f$ f_W(x;\Delta,\xi,\sigma')@f$.
 259      */
 260     Double_t EvaluateWeighted(Double_t x,
 261                               UShort_t maxN=9999) const;
 262     /**
 263      * Find the maximum weight to use.  The maximum weight is the
 264      * largest i for which
 265      *
 266      * - @f$ i \leq \max{N}@f$
 267      * - @f$ a_i > \min{a}@f$
 268      * - @f$ \delta a_i/a_i > \delta_{max}@f$
 269      *
 270      * @param maxRelError @f$ \min{a}@f$
 271      * @param leastWeight @f$ \delta_{max}@f$
 272      * @param maxN        @f$ \max{N}@f$
 273      *
 274      * @return The largest index @f$ i@f$ for which the above
 275      * conditions hold.  Will never return less than 1.
 276      */
 277     Int_t FindMaxWeight(Double_t maxRelError=2*fgMaxRelError,
 278                         Double_t leastWeight=fgLeastWeight,
 279                         UShort_t maxN=999) const;
 280     /* @} */
 281     /**
 282      * @{
 283      * @name TObject Sortable interface
 284      */
 285     /**
 286      * Declare this object as sortable
 287      *
 288      * @return Always true
 289      */
 290     Bool_t IsSortable() const { return kTRUE; }
 291     /**
 292      * Compare to another ELossFit object.
 293      *
 294      * - +1, if this quality is better (larger) than other objects quality
 295      * - -1, if this quality is worse (smaller) than other objects quality
 296      * - +1, if this @f$|\chi^2/\nu-1|@f$ is smaller than the same for other
 297      * - -1, if this @f$|\chi^2/\nu-1|@f$ is larger than the same for other
 298      * - 0 otherwise
 299      *
 300      * @param o Other object to compare to
 301      */
 302     Int_t Compare(const TObject* o) const;
 303     /* @} */
 304     /**
 305      * @{
 306      * name Auxiliary member functions
 307      */
 308     /**
 309      * Information to standard output
 310      *
 311      * @param option Not used
 312      */
 313     void Print(Option_t* option) const; // *MENU*
 314     /**
 315      * Draw this fit
 316      *
 317      * @param option Options
 318      *  - COMP  Draw components too
 319      */
 320     void Draw(Option_t* option="comp"); // *MENU*
 321     /**
 322      * Browse this object
 323      *
 324      * @param b Browser
 325      */
 326     void Browse(TBrowser* b);
 327     /**
 328      * Get the name of this object
 329      *
 330      *
 331      * @return
 332      */
 333     const Char_t* GetName() const;
 334     /**
 335      * Calculate the lower bound
 336      *
 337      * @param f             Width factor
 338      * @param includeSigma  Whether to include sigma
 339      *
 340      * @return @f$ \Delta - f (\xi + \sigma)@f$
 341      */
 342     Double_t GetLowerBound(Double_t f, Bool_t includeSigma) const;
 343     /**
 344      * Calculate the quality
 345      */
 346     void CalculateQuality(Double_t maxChi2nu=fgMaxChi2nu,
 347                           Double_t maxRelError=fgMaxRelError,
 348                           Double_t leastWeight=fgLeastWeight);
 349     /* @} */
 350     ClassDef(ELossFit,1); // Result of fit
 351   };
 352
 353   /**
 354    * Default constructor
 355    */
 356   AliFMDCorrELossFit();
 357   /**
 358    * Copy constructor
 359    *
 360    * @param o Object to copy from
 361    */
 362   AliFMDCorrELossFit(const AliFMDCorrELossFit& o);
 363   /**
 364    * Destructor
 365    */
 366   virtual ~AliFMDCorrELossFit();
 367   /**
 368    * Assignment operator
 369    *
 370    * @param o Object to assign from
 371    *
 372    * @return Reference to this object
 373    */
 374   AliFMDCorrELossFit& operator=(const AliFMDCorrELossFit& o);
 375
 376   /**
 377    * @{
 378    * @name Set fits
 379    */
 380   /**
 381    * Set the fit parameters from a function
 382    *
 383    * @param d        Detector
 384    * @param r        Ring
 385    * @param eta      Eta
 386    * @param quality  Quality flag
 387    * @param f        Function from fit
 388    */
 389   Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, Int_t quality,
 390                 const TF1& f);
 391   /**
 392    * Set the fit parameters from a function
 393    *
 394    * @param d    Detector
 395    * @param r    Ring
 396    * @param eta  Eta
 397    * @param f    ELoss fit result - note, the object will take ownership
 398    */
 399   Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* f);
 400   /**
 401    * Set the fit parameters from a function
 402    *
 403    * @param d       Detector
 404    * @param r       Ring
 405    * @param etaBin  Eta (bin number, 1->nBins)
 406    * @param f       ELoss fit result - note, the object will take ownership
 407    */
 408   Bool_t SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* f);
 409   /**
 410    * Set the fit parameters from a function
 411    *
 412    * @param d         Detector number
 413    * @param r         Ring identifier
 414    * @param eta       Eta value
 415    * @param quality   Quality flag
 416    * @param n         @f$ N@f$ - Number of fitted peaks
 417    * @param chi2      @f$ \chi^2 @f$
 418    * @param nu        @f$ \nu @f$ - number degrees of freedom
 419    * @param c         @f$ C@f$ - scale constant
 420    * @param ec        @f$ \delta C@f$ - error on @f$ C@f$
 421    * @param delta     @f$ \Delta@f$ - most probable value
 422    * @param edelta    @f$ \delta\Delta@f$ - error on @f$\Delta@f$
 423    * @param xi        @f$ \xi@f$ - Landau width
 424    * @param exi       @f$ \delta\xi@f$ - error on @f$\xi@f$
 425    * @param sigma     @f$ \sigma@f$ - Gaussian width
 426    * @param esigma    @f$ \delta\sigma@f$ - error on @f$\sigma@f$
 427    * @param sigman    @f$ \sigma_n@f$ - Noise width
 428    * @param esigman   @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
 429    * @param a         Array of @f$ N-1@f$ weights @f$ a_i@f$ for
 430    *                  @f$ i=2,\ldots@f$
 431    * @param ea        Array of @f$ N-1@f$ errors on weights @f$ a_i@f$ for
 432    *                  @f$ i=2,\ldots@f$
 433    */
 434   Bool_t SetFit(UShort_t  d,      Char_t    r, Double_t eta,
 435                 Int_t     quality,UShort_t  n,
 436                 Double_t  chi2,   UShort_t  nu,
 437                 Double_t  c,      Double_t  ec,
 438                 Double_t  delta,  Double_t  edelta,
 439                 Double_t  xi,     Double_t  exi,
 440                 Double_t  sigma,  Double_t  esigma,
 441                 Double_t  sigman, Double_t  esigman,
 442                 Double_t* a,      Double_t* ea);
 443   /* @} */
 444
 445   /**
 446    * @{
 447    * @name Set and get eta axis
 448    */
 449   /**
 450    * Set the eta axis to use
 451    *
 452    * @param axis Eta axis
 453    */
 454   void SetEtaAxis(const TAxis& axis);
 455   /**
 456    * Set the eta axis to use
 457    *
 458    * @param nBins Number of bins
 459    * @param min   Minimum @f$ \eta@f$
 460    * @param max   maximum @f$ \eta@f$
 461    */
 462   void SetEtaAxis(Int_t nBins, Double_t min, Double_t max);
 463   /**
 464    * Get the eta axis used
 465    *
 466    * @return
 467    */
 468   const TAxis& GetEtaAxis() const { return fEtaAxis; }
 469   /**
 470    * Set the low cut used when fitting
 471    *
 472    * @param cut Cut value
 473    */
 474   void SetLowCut(Double_t cut) { fLowCut = cut; }
 475   /**
 476    * Get the low cut used when fitting
 477    *
 478    * @return Low cut used for fitting
 479    */
 480   Double_t GetLowCut() const { return fLowCut; }
 481   /**
 482    * Find the eta bin corresponding to the given eta
 483    *
 484    * @param eta  Eta value
 485    *
 486    * @return Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
 487    * eta, or 0 if out of range.
 488    */
 489   Int_t FindEtaBin(Double_t eta) const;
 490   /* @} */
 491
 492   /**
 493    * @{
 494    * @name Find fits
 495    */
 496   /**
 497    * Find the fit corresponding to the specified parameters
 498    *
 499    * @param d   Detector
 500    * @param r   Ring
 501    * @param eta Eta value
 502    *
 503    * @return Fit parameters or null in case of problems
 504    */
 505   ELossFit* FindFit(UShort_t d, Char_t r, Double_t eta) const;
 506   /**
 507    * Find the fit corresponding to the specified parameters
 508    *
 509    * @param d      Detector
 510    * @param r      Ring
 511    * @param etabin Eta bin (1 based)
 512    *
 513    * @return Fit parameters or null in case of problems
 514    */
 515   ELossFit* FindFit(UShort_t d, Char_t r, Int_t etabin) const;
 516   /**
 517    * Find the fit corresponding to the specified parameters
 518    *
 519    * @param d   Detector
 520    * @param r   Ring
 521    * @param eta Eta value
 522    *
 523    * @return Fit parameters or null in case of problems
 524    */
 525   ELossFit* GetFit(UShort_t d, Char_t r, Double_t eta) const;
 526   /**
 527    * Find the fit corresponding to the specified parameters
 528    *
 529    * @param d      Detector
 530    * @param r      Ring
 531    * @param etabin Eta bin (1 based)
 532    *
 533    * @return Fit parameters or null in case of problems
 534    */
 535   ELossFit* GetFit(UShort_t d, Char_t r, Int_t etabin) const;
 536   /* @} */
 537
 538   /**
 539    * @{
 540    * @name Lower bounds on fits
 541    */
 542   /**
 543    * Get the lower validity bound of the fit
 544    *
 545    * @param d            Detector
 546    * @param r            Ring
 547    * @param etaBin       Eta bin (1-based)
 548    * @param f            Factor on xi (and sigma)
 549    * @param showErrors   Show errors
 550    * @param includeSigma Whether to include sigma
 551    *
 552    * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$
 553    */
 554   Double_t GetLowerBound(UShort_t d, Char_t r, Int_t etaBin,
 555                          Double_t f, Bool_t showErrors=true,
 556                          Bool_t includeSigma=true) const;
 557   /**
 558    * Get the lower validity bound of the fit
 559    *
 560    * @param d            Detector
 561    * @param r            Ring
 562    * @param eta          Eta value
 563    * @param f            Factor on xi (and sigma)
 564    * @param showErrors   Show errors
 565    * @param includeSigma Whether to include sigma
 566    *
 567    * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$
 568    */
 569   Double_t GetLowerBound(UShort_t d, Char_t r, Double_t eta,
 570                          Double_t f, Bool_t showErrors=true,
 571                          Bool_t includeSigma=true) const;
 572
 573
 574   /**
 575    * @{
 576    * @name Miscellaneous
 577    */
 578   /**
 579    * Get the ring array corresponding to the specified ring
 580    *
 581    * @param d Detector
 582    * @param r Ring
 583    *
 584    * @return Pointer to ring array, or null in case of problems
 585    */
 586   TObjArray* GetRingArray(UShort_t d, Char_t r) const;
 587   /**
 588    * Signal that this is a folder
 589    *
 590    * @return Always true
 591    */
 592   Bool_t IsFolder() const { return true; }
 593   /**
 594    * Browse this object
 595    *
 596    * @param b
 597    */
 598   void Browse(TBrowser* b);
 599   /**
 600    * Draw this object
 601    *
 602    * @param option Options.  Possible values are
 603    *  - error Plot error bars
 604    *  - relative Plot relative errors
 605    */
 606   void Draw(Option_t* option=""); //*MENU*
 607   /**
 608    * Print this object.
 609    *
 610    * @param option Options
 611    *   - R   Print recursive
 612    *
 613    */
 614   void Print(Option_t* option="R") const; //*MENU*
 615   /* @} */
 616 protected:
 617   /**
 618    * Get the ring array corresponding to the specified ring
 619    *
 620    * @param d Detector
 621    * @param r Ring
 622    *
 623    * @return Pointer to ring array, or newly created container
 624    */
 625   TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r);
 626
 627   TObjArray  fRings;    // Array of rings
 628   TAxis      fEtaAxis;  // Eta axis used
 629   Double_t   fLowCut;   // Low cut used when fitting
 630
 631   ClassDef(AliFMDCorrELossFit,2);
 632 };
 633
 634 //____________________________________________________________________
 635 inline void
 636 AliFMDCorrELossFit::SetEtaAxis(Int_t nBins, Double_t min, Double_t max)
 637 {
 638   fEtaAxis.Set(nBins, min, max);
 639 }
 640 //____________________________________________________________________
 641 inline void
 642 AliFMDCorrELossFit::SetEtaAxis(const TAxis& e)
 643 {
 644   fEtaAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
 645 }
 646 //____________________________________________________________________
 647 inline Double_t
 648 AliFMDCorrELossFit::ELossFit::GetA(UShort_t i) const
 649 {
 650   if (i <  1)   return 0;
 651   if (i >  fN)  return 0;
 652   if (i == 1)   return 1;
 653   return fA[i-2];
 654 }
 655 //____________________________________________________________________
 656 inline Double_t
 657 AliFMDCorrELossFit::ELossFit::GetEA(UShort_t i) const
 658 {
 659   if (i <  1)   return 0;
 660   if (i >  fN)  return 0;
 661   if (i == 1)   return 1;
 662   return fEA[i-2];
 663 }
 664
 665
 666 #endif
 667 // Local Variables:
 668 //   mode: C++
 669 // End:
 670