PWGLF/FORWARD/analysis2/AliLandauGaus.h

   1 #ifndef ALILANDAUGAUS_H
   2 #define ALILANDAUGAUS_H
   3 /**
   4  * @file   AliLandauGaus.h
   5  * @author Christian Holm Christensen <cholm@nbi.dk>
   6  * @date   Tue Mar 11 08:53:47 2014
   7  *
   8  * @brief  Declaration and implementation of Landau-Gauss distributions.
   9  *
  10  * @ingroup pwglf_forward
  11  */
  12
  13 #include <TObject.h>
  14 #include <TF1.h>
  15 #include <TMath.h>
  16
  17 /**
  18  * This class contains static member functions to calculate the energy
  19  * loss stragling - most notably the N-particle energy loss as a sum
  20  * of convolutions of a Landau and Gauss distribution.
  21  *
  22  * That is, for a single particle we have the function @f$ f(x)@f$:
  23  *
  24  * @f[
  25  * f(x;\Delta_p,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
  26  *    \int_{-\infty}^{+\infty} dx' f'_{L}(x',\Delta_p,\xi)
  27  *    \exp{-\frac{(x-x')^2}{2\sigma'^2}}}
  28  * @f]
  29  *
  30  * where @f$ f'_{L}@f$ is the Landau distribution, @f$\Delta_p@f$ the
  31  * most probable energy loss, @f$ \xi@f$ the width of the Landau, and
  32  * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
  33  * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter
  34  * modelling noise in the detector.
  35  *
  36  * for @f$ i@f$ particles this is modified to
  37  *
  38  * @f[
  39  *    f_i(x;\Delta_{p},\xi,\sigma')=f(x;\Delta_{p,i},\xi_i,\sigma_i')
  40  * @f]
  41  *
  42  * corresponding to @f$ i@f$ particles i.e., with the substitutions
  43  * @f{eqnarray*}{
  44  *    \Delta_p  \rightarrow \Delta_{p,i}&=& i(\Delta_p + \xi\log(i))\\
  45  *    \xi       \rightarrow \xi_i       &=& i \xi\\
  46  *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
  47  *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
  48  * @f}
  49  *
  50  * Because of the convolution with a Gaussian, the most-probable-value
  51  * @f$\Delta_p'@f$ of the resulting distribution is not really at the
  52  * Landau most-probable-value @f$\Delta_p@f$.  In fact we find that
  53  * @f$\Delta_p' > \Delta_p@f$.
  54  *
  55  * Ideally, one would find an analytic expression for this shift by
  56  * solving
  57  *
  58  * @f{eqnarray}
  59  *   0 &=& \frac{d f_i(x;\Delta_p,\xi,\sigma)}{d\Delta}\\
  60  *     &=& \frac{d}{d\Delta}\frac{1}{\sigma' \sqrt{2 \pi}}
  61  *      \int_{-\infty}^{+\infty} dx' f'_{L}(x',\Delta_p,\xi)
  62  *      \exp{-\frac{(x-x')^2}{2\sigma'^2}}}
  63  * @f{eqnarray}
  64  *
  65  * for @f$x@f$ as a function of @f$\Delta_p,\xi,\sigma,i@f$. However,
  66  * do to the complex nature of the Landau distribution this is not
  67  * really feasible.
  68  *
  69  * Instead, the shift was studied numerically. Landau-Gauss
  70  * distributions for @f$i=1,\ldots@f$ where generated with varying
  71  * @f$\xi@f$ and @f$\sigma@f$.  The distributions was then numerically
  72  * differentiated and the root @f$\Delta_p'@f$ of that derivative
  73  * found numerically.  The difference
  74  * @f$\delta\Delta_p=\Delta_p'-\Delta_p@f$ was then studied as a
  75  * function of the @f$\sigma,\xi@f$ parameters and an approximate
  76  * expression was found
  77  *
  78  * @f[
  79  *  \delta\Delta_p \approx \frac{c \sigma u}{(1+1/i)^{p u^{3/2}}}
  80  * @f]
  81  *
  82  * where @f$ u=\sigma/\xi@f$.  The parameters @f$c@f$ and @f$p@f$ is
  83  * found to depend on @f$ u@f$ only weakly, and for practical
  84  * applications where @f$u\approx1@f$, we set @f$ c=p=1/2@f$.
  85  *
  86  * For the evaluating the full energy loss distribution from f@$
  87  * 1+2+\ldots,n@f$ particles, we evaluate
  88  *
  89  * @f[
  90  *   f_N(x;\Delta_p,\xi,\sigma')=\sum_{i=1}^N a_i f_i(x;\Delta_p,\xi,\sigma',a)
  91  * @f]
  92  *
  93  * where @f$ f(x;\Delta_p,\xi,\sigma')@f$ is the convolution of a
  94  * Landau with a Gaussian (see LandauGaus), and @f$ a@f$ is a vector of
  95  * weights for each @f$ f_i@f$. Note that @f$ a_1 = 1@f$.
  96  *
  97  * Everything is defined in this header file to make it easy to move
  98  * this code around. Nothing here's meant to be persistent, so we
  99  * can easily do that.
 100  *
 101  * References:
 102  *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">
 103  *        Nucl.Instrum.Meth.B1:16</a>
 104  *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 105  *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">
 106  *        ROOT implementation</a>
 107  *
 108  * @ingroup pwglf_forward
 109  */
 110 class AliLandauGaus
 111 {
 112 public:
 113   /** Enumeration of parameters */
 114   enum {
 115     kC     = 0,
 116     kDelta,
 117     kXi,
 118     kSigma,
 119     kSigmaN,
 120     kN,
 121     kA
 122   };
 123   /** Enumeration of colors */
 124
 125   //__________________________________________________________________
 126   /**
 127    * @{
 128    * @name Constants
 129    */
 130   //------------------------------------------------------------------
 131   /**
 132    * The shift of the most probable value for the ROOT function
 133    * TMath::Landau
 134    *
 135    * @return  Shift of TMath::Landau
 136    */
 137   static Double_t MPShift() { return -0.22278298; }
 138   /**
 139    * Constant of @f$\sigma@f$ shift
 140    *
 141    * @return Constant of @f$\sigma@f$ shift
 142    */
 143   static Double_t SigmaShiftC() { return 0.5446; }
 144   /**
 145    * Power factor of @f$\sigma@f$ shift
 146    *
 147    * @return Power factor of @f$\sigma@f$ shift
 148    */
 149   static Double_t SigmaShiftP() { return 0.5746; }
 150   /**
 151    * Normalization constant
 152    *
 153    * @return The landau-gauss normalization constant
 154    */
 155   static Double_t InvSq2Pi() { return 1. / TMath::Sqrt(TMath::TwoPi()); }
 156   /**
 157    * How many sigma's of the Gaussian in the Landau, Gaussian
 158    * convolution to integrate over
 159    */
 160   static Double_t NSigma() { return 5; }
 161   /**
 162    * Number of steps to do in the Landau, Gaussiam convolution
 163    */
 164   static Int_t NSteps() { return 100; }
 165   /* @} */
 166
 167   //__________________________________________________________________
 168   /**
 169    * @{
 170    * @name Function calculations
 171    */
 172   //------------------------------------------------------------------
 173   /**
 174    * Calculate the shifted Landau
 175    * @f[
 176    *    f'_{L}(x;\Delta_p,\xi) = f_L(x;\Delta_p+0.22278298\xi)
 177    * @f]
 178    *
 179    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 180    * distribution (known to have @f$\Delta_{p}=-0.22278298@f$ for
 181    * @f$\Delta_p=0,\xi=1@f$.
 182    *
 183    * @param x      Where to evaluate @f$ f'_{L}@f$
 184    * @param delta  Most probable value
 185    * @param xi     The 'width' of the distribution
 186    *
 187    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
 188    */
 189   static Double_t Fl(Double_t x, Double_t delta, Double_t xi);
 190   //------------------------------------------------------------------
 191   /**
 192    * Calculate the value of a Landau convolved with a Gaussian
 193    *
 194    * @f[
 195    * f(x;\Delta_p,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 196      *    \int_{-\infty}^{+\infty} dx' f'_{L}(x';\Delta_p,\xi)
 197    *    \exp{-\frac{(x-x')^2}{2\sigma'^2}}
 198    * @f]
 199    *
 200    * Note that this function uses the constants NSteps() and
 201    * NSigma()
 202    *
 203    * @param x         where to evaluate @f$ f@f$
 204    * @param delta     @f$ \Delta_p@f$ of @f$ f(x;\Delta_p,\xi,\sigma')@f$
 205    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta_p,\xi,\sigma')@f$
 206    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 207    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 208    *
 209    * @return @f$ f@f$ evaluated at @f$ x@f$.
 210    */
 211   static Double_t F(Double_t x, Double_t delta, Double_t xi,
 212                     Double_t sigma, Double_t sigma_n);
 213   //------------------------------------------------------------------
 214   /**
 215    * Evaluate
 216    * @f[
 217    *    f_i(x;\Delta_p,\xi,\sigma') = f(x;\Delta_{p,i},\xi_i,\sigma_i')
 218    * @f]
 219    * corresponding to @f$ i@f$ particles.
 220    *
 221    * @param x        Where to evaluate
 222    * @param delta    @f$ \Delta@f$
 223    * @param xi       @f$ \xi@f$
 224    * @param sigma    @f$ \sigma@f$
 225    * @param sigma_n  @f$ \sigma_n@f$
 226    * @param i        @f$ i @f$
 227    *
 228    * @return @f$ f_i @f$ evaluated
 229    */
 230   static Double_t Fi(Double_t x, Double_t delta, Double_t xi,
 231                      Double_t sigma, Double_t sigma_n, Int_t i);
 232
 233   //------------------------------------------------------------------
 234   /**
 235    * Numerically evaluate
 236    * @f[
 237    *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 238    * @f]
 239    * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 240    * of the parameters is given by
 241    *
 242    * - 0: @f$\Delta@f$
 243    * - 1: @f$\xi@f$
 244    * - 2: @f$\sigma@f$
 245    * - 3: @f$\sigma_n@f$
 246    *
 247    * @param x        Where to evaluate
 248    * @param ipar     Parameter number
 249    * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 250    * @param delta    @f$ \Delta_p@f$
 251    * @param xi       @f$ \xi@f$
 252    * @param sigma    @f$ \sigma@f$
 253    * @param sigma_n  @f$ \sigma_n@f$
 254    * @param i        @f$ i@f$
 255    *
 256    * @return @f$ f_i@f$ evaluated
 257    */
 258   static Double_t DFidPar(Double_t x, UShort_t ipar, Double_t dp,
 259                           Double_t delta, Double_t xi,
 260                           Double_t sigma, Double_t sigma_n, Int_t i);
 261
 262   //------------------------------------------------------------------
 263   /**
 264    * Evaluate
 265    * @f[
 266    * f_N(x;\Delta_p,\xi,\sigma')=\sum_{i=1}^N a_i f_i(x;\Delta_p,\xi,\sigma',a)
 267    * @f]
 268    *
 269    * where @f$ f(x;\Delta_p,\xi,\sigma')@f$ is the convolution of a
 270      * Landau with a Gaussian (see LandauGaus).
 271      *
 272      * @param x        Where to evaluate @f$ f_N@f$
 273      * @param delta    @f$ \Delta_1@f$
 274      * @param xi       @f$ \xi_1@f$
 275      * @param sigma    @f$ \sigma_1@f$
 276      * @param sigma_n  @f$ \sigma_n@f$
 277      * @param n        @f$ N@f$ in the sum above.
 278      * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 279      *                 @f$ i > 1@f$
 280      *
 281      * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
 282    */
 283   static Double_t Fn(Double_t x, Double_t delta, Double_t xi,
 284                      Double_t sigma, Double_t sigma_n, Int_t n,
 285                      const Double_t* a);
 286   /**
 287    * Get parameters for the @f$ i@f$ particle response.
 288    *
 289    * @f{eqnarray*}{
 290    *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\
 291    *    \xi       \rightarrow \xi_i       &=& i \xi\\
 292    *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
 293    *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 294    * @f}
 295    *
 296    *
 297    * @param i     Number of particles
 298    * @param delta Input: single particle, output @f$ i@f$ particle
 299    * @param xi    Input: single particle, output @f$ i@f$ particle
 300    * @param sigma Input: single particle, output @f$ i@f$ particle
 301    */
 302   static void IPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma);
 303   /**
 304    * Set and check if sigma shift is enabled
 305    *
 306    * @param val if <0, then only check.  Otherwise set enabled (>0) or not (=0)
 307    *
 308    * @return whether the sigma shift is enabled or not
 309    */
 310   static Bool_t EnableSigmaShift(Short_t val=-1);
 311   /**
 312    * Get the shift of the MPV due to convolution with a Gaussian.
 313    *
 314    * @f[
 315    *  \delta\Delta_p \approx \frac{c \sigma u}{(1+1/i)^{p u^{3/2}}}
 316    * @f]
 317    *
 318    * where @f$ u=\sigma/\xi@f$.
 319    *
 320    * @param i       Number of particles
 321    * @param xi      Landau width @f$\xi@f$
 322    * @param sigma   Gaussian variance @f$\sigma@f$
 323    *
 324    * @return The shift
 325    */
 326   static Double_t SigmaShift(Int_t i, Double_t xi, Double_t sigma);
 327   /* @} */
 328
 329
 330   //__________________________________________________________________
 331   /**
 332    * @{
 333    * @name Utilities for defining TF1 objects
 334    */
 335   //------------------------------------------------------------------
 336   /**
 337    * Generate a TF1 object of @f$ f_1@f$
 338    *
 339    * @param c        Constant
 340    * @param delta    @f$ \Delta_1@f$
 341    * @param xi       @f$ \xi_1@f$
 342    * @param sigma    @f$ \sigma_1@f$
 343    * @param sigma_n  @f$ \sigma_n@f$
 344    * @param xmin     Least value of range
 345    * @param xmax     Largest value of range
 346    *
 347    * @return Newly allocated TF1 object
 348    */
 349   static TF1* MakeF1(Double_t c,
 350                      Double_t delta, Double_t xi,
 351                      Double_t sigma, Double_t sigma_n,
 352                      Double_t xmin,  Double_t  xmax);
 353   //------------------------------------------------------------------
 354   /**
 355    * Generate a TF1 object of @f$ f_I@f$
 356    *
 357    * @param c        Constant
 358    * @param delta    @f$ \Delta_1@f$
 359    * @param xi       @f$ \xi_1@f$
 360    * @param sigma    @f$ \sigma_1@f$
 361    * @param sigma_n  @f$ \sigma_n@f$
 362    * @param i        @f$ i@f$ - the number of particles
 363    * @param xmin     Least value of range
 364    * @param xmax     Largest value of range
 365    *
 366    * @return Newly allocated TF1 object
 367    */
 368   static TF1* MakeFi(Double_t c,
 369                      Double_t delta, Double_t xi,
 370                      Double_t sigma, Double_t sigma_n,
 371                      Int_t    i,
 372                      Double_t xmin,  Double_t  xmax);
 373   //------------------------------------------------------------------
 374   /**
 375    * Generate a TF1 object of @f$ f_N@f$
 376    *
 377    * @param c         Constant
 378    * @param delta     @f$ \Delta_1@f$
 379    * @param xi        @f$ \xi_1@f$
 380    * @param sigma     @f$ \sigma_1@f$
 381    * @param sigma_n   @f$ \sigma_n@f$
 382    * @param n         @f$ N@f$ - how many particles to sum to
 383    * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 384    *                  @f$ i > 1@f$
 385    * @param xmin      Least value of range
 386    * @param xmax      Largest value of range
 387    *
 388    * @return Newly allocated TF1 object
 389    */
 390   static TF1* MakeFn(Double_t c,
 391                      Double_t delta, Double_t  xi,
 392                      Double_t sigma, Double_t  sigma_n,
 393                      Int_t    n,     const Double_t* a,
 394                      Double_t xmin,  Double_t  xmax);
 395   /**
 396    * Make a TF1 object that describes the convoluted energy loss of
 397    * primary and secondary particles.
 398    *
 399    * @param c1     Primary weight
 400    * @param delta  Most-probable value
 401    * @param xi1    Primary width
 402    * @param sigma  Gaussian variance
 403    * @param c2     Secondary weight
 404    * @param xi2    Secondary with
 405    *
 406    * @return Pointer to newly allocated TF1 object
 407    */
 408   static TF1* MakeComposite(Double_t c1,
 409                             Double_t delta,
 410                             Double_t xi1,
 411                             Double_t sigma,
 412                             Double_t c2,
 413                             Double_t xi2,
 414                             Double_t xmin,
 415                             Double_t xmax);
 416
 417   //------------------------------------------------------------------
 418   /**
 419    * Get the color of the @f$ i@f$ particle response
 420    *
 421    * @param i Particle number
 422    *
 423    * @return Color
 424    */
 425   static Color_t GetIColor(Int_t i);
 426   //------------------------------------------------------------------
 427   /**
 428    * Utility function for TF1 definition
 429    *
 430    * @param pp Array of parameters
 431    * @param xp Pointer to independent variable
 432    *
 433    * @see AliLandauGaus::F
 434    *
 435    * @return Landau convolved with a Gauss
 436    */
 437   static Double_t F1Func(Double_t* xp, Double_t* pp);
 438   //------------------------------------------------------------------
 439   /**
 440    * Utility function for TF1 definition
 441    *
 442    * @param pp Array of parameters
 443    * @param xp Pointer to independent variable
 444    *
 445    * @see AliLandauGaus::Fn
 446    *
 447    * @return Landau convolved with a Gauss
 448    */
 449   static Double_t FnFunc(Double_t* xp, Double_t* pp);
 450   //------------------------------------------------------------------
 451   /**
 452    * Utility function for TF1 definition
 453    *
 454    * @param pp Array of parameters
 455    * @param xp Pointer to independent variable
 456    *
 457    * @see AliLandauGaus::Fn
 458    *
 459    * @return Landau convolved with a Gauss
 460    */
 461   static Double_t FiFunc(Double_t* xp, Double_t* pp);
 462   //------------------------------------------------------------------
 463   /**
 464    * Utility function for TF1 definition
 465    *
 466    * @param pp Array of parameters
 467    * @param xp Pointer to independent variable
 468    *
 469    * @see AliLandauGaus::F
 470    *
 471    * @return Landau convolved with a Gauss
 472    */
 473   static Double_t CompFunc(Double_t* xp, Double_t* pp);
 474   /* @} */
 475 };
 476 //____________________________________________________________________
 477 inline Bool_t
 478 AliLandauGaus::EnableSigmaShift(Short_t val)
 479 {
 480   static Bool_t enabled = true;
 481   if (val >= 0) enabled = val == 1;
 482   return enabled;
 483 }
 484 //____________________________________________________________________
 485 inline void
 486 AliLandauGaus::IPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma)
 487 {
 488 #ifndef NO_SIGMA_SHIFT
 489   Double_t dDelta = EnableSigmaShift() ? SigmaShift(i, xi, sigma) : 0;
 490 #else
 491   // This is for testing
 492   Double_t dDelta = 0;
 493 #endif
 494   if (i == 1) {
 495     delta -= dDelta;
 496     return;
 497   }
 498
 499   delta = i * (delta + xi * TMath::Log(i)) - dDelta;
 500   xi    = i * xi;
 501   sigma = TMath::Sqrt(Double_t(i)) * sigma;
 502 }
 503 //____________________________________________________________________
 504 inline Double_t
 505 AliLandauGaus::SigmaShift(Int_t i, Double_t xi, Double_t sigma)
 506 {
 507   if (xi <= 0) return 0;
 508   if (sigma <= 0) return 0;
 509   const Double_t c = SigmaShiftC();
 510   const Double_t p = SigmaShiftP();
 511   const Double_t u = sigma / xi;
 512   const Double_t q = p*u*TMath::Sqrt(u);
 513   if (q > 100) return 0;
 514   return c * sigma / TMath::Power(1+1./i, q);
 515 }
 516 //____________________________________________________________________
 517 inline Double_t
 518 AliLandauGaus::Fl(Double_t x, Double_t delta, Double_t xi)
 519 {
 520   Double_t deltaP = delta - xi * MPShift();
 521   return TMath::Landau(x, deltaP, xi, true);
 522 }
 523 //____________________________________________________________________
 524 inline Double_t
 525 AliLandauGaus::F(Double_t x, Double_t delta, Double_t xi,
 526                  Double_t sigma, Double_t sigmaN)
 527 {
 528   if (xi <= 0) return 0;
 529
 530   const Int_t    nSteps = NSteps();
 531   const Double_t nSigma = NSigma();
 532   const Double_t deltaP = delta; // - sigma * sigmaShift; // + sigma * mpshift;
 533   const Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
 534   const Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
 535   const Double_t xlow   = x - nSigma * sigma1;
 536   const Double_t xhigh  = x + nSigma * sigma1;
 537   const Double_t step   = (xhigh - xlow) / nSteps;
 538   Double_t       sum    = 0;
 539
 540   for (Int_t i = 0; i <= nSteps/2; i++) {
 541     const Double_t x1 = xlow  + (i - .5) * step;
 542     const Double_t x2 = xhigh - (i - .5) * step;
 543     sum += Fl(x1, deltaP, xi) * TMath::Gaus(x, x1, sigma1);
 544     sum += Fl(x2, deltaP, xi) * TMath::Gaus(x, x2, sigma1);
 545   }
 546   return step * sum * InvSq2Pi() / sigma1;
 547 }
 548
 549 //____________________________________________________________________
 550 inline Double_t
 551 AliLandauGaus::Fi(Double_t x, Double_t delta, Double_t xi,
 552                   Double_t sigma, Double_t sigmaN, Int_t i)
 553 {
 554   Double_t deltaI = delta;
 555   Double_t xiI    = xi;
 556   Double_t sigmaI = sigma;
 557   IPars(i, deltaI, xiI, sigmaI);
 558   if (sigmaI < 1e-10)
 559     // Fall back to landau
 560     return Fl(x, deltaI, xiI);
 561
 562   return F(x, deltaI, xiI, sigmaI, sigmaN);
 563 }
 564 //____________________________________________________________________
 565 inline Double_t
 566 AliLandauGaus::Fn(Double_t x, Double_t delta, Double_t xi,
 567                   Double_t sigma, Double_t sigmaN, Int_t n,
 568                   const Double_t* a)
 569 {
 570   Double_t result = Fi(x, delta, xi, sigma, sigmaN, 1);
 571   for (Int_t i = 2; i <= n; i++)
 572     result += a[i-2] * Fi(x,delta,xi,sigma,sigmaN,i);
 573   return result;
 574 }
 575
 576 //____________________________________________________________________
 577 inline Double_t
 578 AliLandauGaus::DFidPar(Double_t x,
 579                        UShort_t par,   Double_t dPar,
 580                        Double_t delta, Double_t xi,
 581                        Double_t sigma, Double_t sigmaN,
 582                        Int_t    i)
 583 {
 584   if (dPar == 0) return 0;
 585   Double_t dp      = dPar;
 586   Double_t d2      = dPar / 2;
 587   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
 588   Double_t xiI     =  i * xi;
 589   Double_t si      =  TMath::Sqrt(Double_t(i));
 590   Double_t sigmaI  =  si*sigma;
 591   Double_t y1      = 0;
 592   Double_t y2      = 0;
 593   Double_t y3      = 0;
 594   Double_t y4      = 0;
 595   switch (par) {
 596   case 0:
 597     y1 = Fi(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
 598     y2 = Fi(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
 599     y3 = Fi(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
 600     y4 = Fi(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
 601     break;
 602   case 1:
 603     y1 = Fi(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
 604     y2 = Fi(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
 605     y3 = Fi(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
 606     y4 = Fi(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
 607     break;
 608   case 2:
 609     y1 = Fi(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
 610     y2 = Fi(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
 611     y3 = Fi(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
 612     y4 = Fi(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
 613     break;
 614   case 3:
 615     y1 = Fi(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
 616     y2 = Fi(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
 617     y3 = Fi(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
 618     y4 = Fi(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
 619     break;
 620   default:
 621     return 0;
 622   }
 623
 624   Double_t d0  = y1 - y4;
 625   Double_t d1  = 2 * (y2 - y3);
 626
 627   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
 628
 629   return g;
 630 }
 631
 632
 633 //____________________________________________________________________
 634 inline Color_t
 635 AliLandauGaus::GetIColor(Int_t i)
 636 {
 637   const Int_t kColors[] = { kRed+1,
 638                             kPink+3,
 639                             kMagenta+2,
 640                             kViolet+2,
 641                             kBlue+1,
 642                             kAzure+3,
 643                             kCyan+1,
 644                             kTeal+2,
 645                             kGreen+2,
 646                             kSpring+3,
 647                             kYellow+2,
 648                             kOrange+2 };
 649   return kColors[((i-1) % 12)];
 650 }
 651 //____________________________________________________________________
 652 inline Double_t
 653 AliLandauGaus::F1Func(Double_t* xp, Double_t* pp)
 654 {
 655   Double_t x        = xp[0];
 656   Double_t constant = pp[kC];
 657   Double_t delta    = pp[kDelta];
 658   Double_t xi       = pp[kXi];
 659   Double_t sigma    = pp[kSigma];
 660   Double_t sigmaN   = pp[kSigmaN];
 661
 662   return constant * F(x, delta, xi, sigma, sigmaN);
 663 }
 664 //____________________________________________________________________
 665 inline Double_t
 666 AliLandauGaus::FiFunc(Double_t* xp, Double_t* pp)
 667 {
 668   Double_t x        = xp[0];
 669   Double_t constant = pp[kC];
 670   Double_t delta    = pp[kDelta];
 671   Double_t xi       = pp[kXi];
 672   Double_t sigma    = pp[kSigma];
 673   Double_t sigmaN   = pp[kSigmaN];
 674   Int_t    i        = Int_t(pp[kN]);
 675
 676   return constant * Fi(x, delta, xi, sigma, sigmaN, i);
 677 }
 678 //____________________________________________________________________
 679 inline Double_t
 680 AliLandauGaus::FnFunc(Double_t* xp, Double_t* pp)
 681 {
 682   Double_t  x        = xp[0];
 683   Double_t constant  = pp[kC];
 684   Double_t delta     = pp[kDelta];
 685   Double_t xi        = pp[kXi];
 686   Double_t sigma     = pp[kSigma];
 687   Double_t sigmaN    = pp[kSigmaN];
 688   Int_t     n        = Int_t(pp[kN]);
 689   Double_t* a        = &(pp[kA]);
 690
 691   return constant * Fn(x, delta, xi, sigma, sigmaN, n, a);
 692 }
 693 //____________________________________________________________________
 694 inline Double_t
 695 AliLandauGaus::CompFunc(Double_t* xp, Double_t* pp)
 696 {
 697   Double_t x           = xp[0];
 698   Double_t cP          = pp[kC];
 699   Double_t deltaP      = pp[kDelta];
 700   Double_t xiP         = pp[kXi];
 701   Double_t sigmaP      = pp[kSigma];
 702   Double_t cS          = pp[kSigma+1];
 703   Double_t deltaS      = deltaP; // pp[kSigma+2];
 704   Double_t xiS         = pp[kSigma+2/*3*/];
 705   Double_t sigmaS      = sigmaP; // pp[kSigma+4];
 706
 707   return (cP * F(x,deltaP,xiP,sigmaP,0) +
 708           cS * F(x,deltaS,xiS,sigmaS,0));
 709
 710 }
 711 //____________________________________________________________________
 712 inline TF1*
 713 AliLandauGaus::MakeF1(Double_t  c,
 714                       Double_t  delta, Double_t xi,
 715                       Double_t  sigma, Double_t sigmaN,
 716                       Double_t  xmin, Double_t xmax)
 717 {
 718   // Define the function to fit
 719   TF1* f = new TF1("landau1", &F1Func, xmin,xmax,kSigmaN+1);
 720
 721   // Set initial guesses, parameter names, and limits
 722   f->SetParameters(c,delta,xi,sigma,sigmaN);
 723   f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
 724   f->SetNpx(500);
 725   f->SetLineColor(GetIColor(1));
 726
 727   return f;
 728 }
 729
 730 //____________________________________________________________________
 731 inline TF1*
 732 AliLandauGaus::MakeFn(Double_t  c,
 733                       Double_t  delta, Double_t xi,
 734                       Double_t  sigma, Double_t sigmaN, Int_t n,
 735                       const Double_t* a,
 736                       Double_t  xmin, Double_t xmax)
 737 {
 738   Int_t npar = kN+n;
 739   TF1*  f    = new TF1(Form("nlandau%d", n), &FnFunc,xmin,xmax,npar);
 740   f->SetLineColor(GetIColor(n));
 741   f->SetLineWidth(2);
 742   f->SetNpx(500);
 743   f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
 744   f->SetParameter(kC,      c);
 745   f->SetParameter(kDelta,  delta);
 746   f->SetParameter(kXi,     xi);
 747   f->SetParameter(kSigma,  sigma);
 748   f->SetParameter(kSigmaN, sigmaN);
 749   f->FixParameter(kN,      n);
 750   for (UShort_t i = 2; i <= n; i++) {
 751     f->SetParameter(kA+i-2, a[i-2]);
 752     f->SetParName(kA+i-2, Form("a_{%d}", i));
 753   }
 754   return f;
 755 }
 756 //____________________________________________________________________
 757 inline TF1*
 758 AliLandauGaus::MakeFi(Double_t  c,
 759                       Double_t  delta, Double_t xi,
 760                       Double_t  sigma, Double_t sigmaN, Int_t i,
 761                       Double_t  xmin, Double_t xmax)
 762 {
 763   Int_t npar = kN+1;
 764   TF1*  f    = new TF1(Form("ilandau%d", i), &FiFunc,xmin,xmax,npar);
 765   f->SetLineColor(GetIColor(i));
 766   f->SetLineWidth(1);
 767   f->SetNpx(500);
 768   f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
 769   f->SetParameter(kC,      c);
 770   f->SetParameter(kDelta,  delta);
 771   f->SetParameter(kXi,     xi);
 772   f->SetParameter(kSigma,  sigma);
 773   f->SetParameter(kSigmaN, sigmaN);
 774   f->FixParameter(kN,      i);
 775
 776   return f;
 777 }
 778 //____________________________________________________________________
 779 inline TF1*
 780 AliLandauGaus::MakeComposite(Double_t c1,
 781                              Double_t delta,
 782                              Double_t xi1,
 783                              Double_t sigma,
 784                              Double_t c2,
 785                              Double_t xi2,
 786                              Double_t xmin,
 787                              Double_t xmax)
 788 {
 789   TF1* comp = new TF1("composite", &CompFunc, xmin, xmax, kSigma+1+2);
 790   comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
 791                     "C#prime", "#xi#prime");
 792   comp->SetParameters(c1,     // 0 Primary weight
 793                       delta,  // 1 Primary Delta
 794                       xi1,    // 2 primary Xi
 795                       sigma,  // 3 primary sigma
 796                       c2,     // 5 Secondary weight
 797                       xi2);   // 6 secondary Xi
 798   return comp;
 799 }
 800 #endif
 801 // Local Variables:
 802 //  mode: C++
 803 // End:
 804
 805
 806