AliLandauGaus.h

Go to the documentation of this file.

#ifndef ALILANDAUGAUS_H
#define ALILANDAUGAUS_H

#include <TObject.h>
#include <TF1.h>
#include <TMath.h>

class AliLandauGaus
{
public:
  enum { 
    kC     = 0,
    kDelta, 
    kXi, 
    kSigma, 
    kSigmaN, 
    kN, 
    kA
  };
  //__________________________________________________________________
  //------------------------------------------------------------------
  static Double_t MPShift() { return -0.22278298; }
  static Double_t SigmaShiftC() { return 0.5446; }
  static Double_t SigmaShiftP() { return 0.5746; }
  static Double_t InvSq2Pi() { return 1. / TMath::Sqrt(TMath::TwoPi()); }
  static Double_t NSigma() { return 5; }
  static Int_t NSteps() { return 100; }
  /* @} */

  //__________________________________________________________________
  //------------------------------------------------------------------
  static Double_t Fl(Double_t x, Double_t delta, Double_t xi);
  //------------------------------------------------------------------
  static Double_t F(Double_t x, Double_t delta, Double_t xi, 
        Double_t sigma, Double_t sigma_n);
  //------------------------------------------------------------------
  static Double_t Fi(Double_t x, Double_t delta, Double_t xi, 
         Double_t sigma, Double_t sigma_n, Int_t i);

  //------------------------------------------------------------------
  static Double_t DFidPar(Double_t x, UShort_t ipar, Double_t dp,
        Double_t delta, Double_t xi, 
        Double_t sigma, Double_t sigma_n, Int_t i);

  //------------------------------------------------------------------
  static Double_t Fn(Double_t x, Double_t delta, Double_t xi, 
         Double_t sigma, Double_t sigma_n, Int_t n, 
         const Double_t* a);
  static void IPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma);
  static Bool_t EnableSigmaShift(Short_t val=-1);
  static Double_t SigmaShift(Int_t i, Double_t xi, Double_t sigma);
  /* @} */

  
  //__________________________________________________________________
  //------------------------------------------------------------------
  static TF1* MakeF1(Double_t c, 
         Double_t delta, Double_t xi, 
         Double_t sigma, Double_t sigma_n,
         Double_t xmin,  Double_t  xmax);
  //------------------------------------------------------------------
  static TF1* MakeFi(Double_t c, 
         Double_t delta, Double_t xi, 
         Double_t sigma, Double_t sigma_n,
         Int_t    i, 
         Double_t xmin,  Double_t  xmax);
  //------------------------------------------------------------------
  static TF1* MakeFn(Double_t c, 
         Double_t delta, Double_t  xi, 
         Double_t sigma, Double_t  sigma_n,
         Int_t    n,     const Double_t* a, 
         Double_t xmin,  Double_t  xmax);
  static TF1* MakeComposite(Double_t c1, 
          Double_t delta, 
          Double_t xi1, 
          Double_t sigma, 
          Double_t c2,
          Double_t xi2, 
          Double_t xmin, 
          Double_t xmax);
  
  //------------------------------------------------------------------
  static Color_t GetIColor(Int_t i);
  //------------------------------------------------------------------
  static Double_t F1Func(Double_t* xp, Double_t* pp);
  //------------------------------------------------------------------
  static Double_t FnFunc(Double_t* xp, Double_t* pp);
  //------------------------------------------------------------------
  static Double_t FiFunc(Double_t* xp, Double_t* pp);
  //------------------------------------------------------------------
  static Double_t CompFunc(Double_t* xp, Double_t* pp);
  /* @} */
};
//____________________________________________________________________
inline Bool_t
AliLandauGaus::EnableSigmaShift(Short_t val)
{
  static Bool_t enabled = true;
  if (val >= 0) enabled = val == 1;
  return enabled;
}
//____________________________________________________________________
inline void
AliLandauGaus::IPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma)
{
#ifndef NO_SIGMA_SHIFT
  Double_t dDelta = EnableSigmaShift() ? SigmaShift(i, xi, sigma) : 0;
#else 
  // This is for testing 
  Double_t dDelta = 0;
#endif
  if (i == 1) {
    delta -= dDelta;
    return;
  }
  
  delta = i * (delta + xi * TMath::Log(i)) - dDelta;
  xi    = i * xi;
  sigma = TMath::Sqrt(Double_t(i)) * sigma;
}
//____________________________________________________________________
inline Double_t
AliLandauGaus::SigmaShift(Int_t i, Double_t xi, Double_t sigma)
{
  if (xi <= 0) return 0;
  if (sigma <= 0) return 0;
  const Double_t c = SigmaShiftC();
  const Double_t p = SigmaShiftP();
  const Double_t u = sigma / xi;
  const Double_t q = p*u*TMath::Sqrt(u);
  if (q > 100) return 0;
  return c * sigma / TMath::Power(1+1./i, q);
}
//____________________________________________________________________
inline Double_t 
AliLandauGaus::Fl(Double_t x, Double_t delta, Double_t xi)
{
  Double_t deltaP = delta - xi * MPShift();
  return TMath::Landau(x, deltaP, xi, true);
}
//____________________________________________________________________
inline Double_t 
AliLandauGaus::F(Double_t x, Double_t delta, Double_t xi,
     Double_t sigma, Double_t sigmaN)
{
  if (xi <= 0) return 0;

  const Int_t    nSteps = NSteps();
  const Double_t nSigma = NSigma();
  const Double_t deltaP = delta; // - sigma * sigmaShift; // + sigma * mpshift;
  const Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
  const Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
  const Double_t xlow   = x - nSigma * sigma1;
  const Double_t xhigh  = x + nSigma * sigma1;
  const Double_t step   = (xhigh - xlow) / nSteps;
  Double_t       sum    = 0;
  
  for (Int_t i = 0; i <= nSteps/2; i++) { 
    const Double_t x1 = xlow  + (i - .5) * step;
    const Double_t x2 = xhigh - (i - .5) * step;
    sum += Fl(x1, deltaP, xi) * TMath::Gaus(x, x1, sigma1);
    sum += Fl(x2, deltaP, xi) * TMath::Gaus(x, x2, sigma1);
  }
  return step * sum * InvSq2Pi() / sigma1;
}

//____________________________________________________________________
inline Double_t 
AliLandauGaus::Fi(Double_t x, Double_t delta, Double_t xi, 
      Double_t sigma, Double_t sigmaN, Int_t i)
{
  Double_t deltaI = delta;
  Double_t xiI    = xi;
  Double_t sigmaI = sigma;
  IPars(i, deltaI, xiI, sigmaI);
  if (sigmaI < 1e-10) 
    // Fall back to landau 
    return Fl(x, deltaI, xiI);
  
  return F(x, deltaI, xiI, sigmaI, sigmaN);
}
//____________________________________________________________________
inline Double_t 
AliLandauGaus::Fn(Double_t x, Double_t delta, Double_t xi, 
      Double_t sigma, Double_t sigmaN, Int_t n, 
      const Double_t* a)
{
  Double_t result = Fi(x, delta, xi, sigma, sigmaN, 1);
  for (Int_t i = 2; i <= n; i++) 
    result += a[i-2] * Fi(x,delta,xi,sigma,sigmaN,i);
  return result;
}

//____________________________________________________________________
inline Double_t 
AliLandauGaus::DFidPar(Double_t x, 
           UShort_t par,   Double_t dPar, 
           Double_t delta, Double_t xi, 
           Double_t sigma, Double_t sigmaN, 
           Int_t    i)
{
  if (dPar == 0) return 0;
  Double_t dp      = dPar;
  Double_t d2      = dPar / 2;
  Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
  Double_t xiI     =  i * xi;
  Double_t si      =  TMath::Sqrt(Double_t(i));
  Double_t sigmaI  =  si*sigma;
  Double_t y1      = 0;
  Double_t y2      = 0;
  Double_t y3      = 0;
  Double_t y4      = 0;
  switch (par) {
  case 0: 
    y1 = Fi(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
    y2 = Fi(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
    y3 = Fi(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
    y4 = Fi(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
    break;
  case 1: 
    y1 = Fi(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
    y2 = Fi(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
    y3 = Fi(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
    y4 = Fi(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
    break;
  case 2: 
    y1 = Fi(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
    y2 = Fi(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
    y3 = Fi(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
    y4 = Fi(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
    break;
  case 3: 
    y1 = Fi(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
    y2 = Fi(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
    y3 = Fi(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
    y4 = Fi(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
    break;
  default:
    return 0;
  } 
  
  Double_t d0  = y1 - y4;
  Double_t d1  = 2 * (y2 - y3);
  
  Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
   
  return g;
}


//____________________________________________________________________
inline Color_t
AliLandauGaus::GetIColor(Int_t i) 
{
  const Int_t kColors[] = { kRed+1, 
          kPink+3, 
          kMagenta+2, 
          kViolet+2, 
          kBlue+1, 
          kAzure+3, 
          kCyan+1, 
          kTeal+2, 
          kGreen+2, 
          kSpring+3, 
          kYellow+2, 
          kOrange+2 };
  return kColors[((i-1) % 12)];
}
//____________________________________________________________________
inline Double_t 
AliLandauGaus::F1Func(Double_t* xp, Double_t* pp) 
{ 
  Double_t x        = xp[0];
  Double_t constant = pp[kC];
  Double_t delta    = pp[kDelta];
  Double_t xi       = pp[kXi];
  Double_t sigma    = pp[kSigma];
  Double_t sigmaN   = pp[kSigmaN];
  
  return constant * F(x, delta, xi, sigma, sigmaN);
}
//____________________________________________________________________
inline Double_t 
AliLandauGaus::FiFunc(Double_t* xp, Double_t* pp) 
{
  Double_t x        = xp[0];
  Double_t constant = pp[kC];
  Double_t delta    = pp[kDelta];
  Double_t xi       = pp[kXi];
  Double_t sigma    = pp[kSigma];
  Double_t sigmaN   = pp[kSigmaN];
  Int_t    i        = Int_t(pp[kN]);

  return constant * Fi(x, delta, xi, sigma, sigmaN, i);
}
//____________________________________________________________________
inline Double_t 
AliLandauGaus::FnFunc(Double_t* xp, Double_t* pp) 
{ 
  Double_t  x        = xp[0];
  Double_t constant  = pp[kC];
  Double_t delta     = pp[kDelta];
  Double_t xi        = pp[kXi];
  Double_t sigma     = pp[kSigma];
  Double_t sigmaN    = pp[kSigmaN];
  Int_t     n        = Int_t(pp[kN]);
  Double_t* a        = &(pp[kA]);
  
  return constant * Fn(x, delta, xi, sigma, sigmaN, n, a);
}
//____________________________________________________________________
inline Double_t 
AliLandauGaus::CompFunc(Double_t* xp, Double_t* pp) 
{
  Double_t x           = xp[0];
  Double_t cP          = pp[kC];
  Double_t deltaP      = pp[kDelta];
  Double_t xiP         = pp[kXi];
  Double_t sigmaP      = pp[kSigma];
  Double_t cS          = pp[kSigma+1];
  Double_t deltaS      = deltaP; // pp[kSigma+2];
  Double_t xiS         = pp[kSigma+2/*3*/];
  Double_t sigmaS      = sigmaP; // pp[kSigma+4];
  
  return (cP * F(x,deltaP,xiP,sigmaP,0) + 
    cS * F(x,deltaS,xiS,sigmaS,0));

}
//____________________________________________________________________
inline TF1*
AliLandauGaus::MakeF1(Double_t  c, 
          Double_t  delta, Double_t xi, 
          Double_t  sigma, Double_t sigmaN,
          Double_t  xmin, Double_t xmax)
{
  // Define the function to fit 
  TF1* f = new TF1("landau1", &F1Func, xmin,xmax,kSigmaN+1);

  // Set initial guesses, parameter names, and limits  
  f->SetParameters(c,delta,xi,sigma,sigmaN);
  f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
  f->SetNpx(500);
  f->SetLineColor(GetIColor(1));
  
  return f;
}
  
//____________________________________________________________________
inline TF1*
AliLandauGaus::MakeFn(Double_t  c, 
          Double_t  delta, Double_t xi, 
          Double_t  sigma, Double_t sigmaN, Int_t n, 
          const Double_t* a, 
          Double_t  xmin, Double_t xmax)
{
  Int_t npar = kN+n;
  TF1*  f    = new TF1(Form("nlandau%d", n), &FnFunc,xmin,xmax,npar);
  f->SetLineColor(GetIColor(n)); 
  f->SetLineWidth(2);
  f->SetNpx(500);
  f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
  f->SetParameter(kC,      c);       
  f->SetParameter(kDelta,  delta);   
  f->SetParameter(kXi,     xi);      
  f->SetParameter(kSigma,  sigma);   
  f->SetParameter(kSigmaN, sigmaN); 
  f->FixParameter(kN,      n);       
  for (UShort_t i = 2; i <= n; i++) {
    f->SetParameter(kA+i-2, a[i-2]);
    f->SetParName(kA+i-2, Form("a_{%d}", i));
  }
  return f;
}
//____________________________________________________________________
inline TF1*
AliLandauGaus::MakeFi(Double_t  c, 
          Double_t  delta, Double_t xi, 
          Double_t  sigma, Double_t sigmaN, Int_t i, 
          Double_t  xmin, Double_t xmax)
{
  Int_t npar = kN+1;
  TF1*  f    = new TF1(Form("ilandau%d", i), &FiFunc,xmin,xmax,npar);
  f->SetLineColor(GetIColor(i));
  f->SetLineWidth(1);
  f->SetNpx(500);
  f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
  f->SetParameter(kC,      c);       
  f->SetParameter(kDelta,  delta);   
  f->SetParameter(kXi,     xi);      
  f->SetParameter(kSigma,  sigma);   
  f->SetParameter(kSigmaN, sigmaN); 
  f->FixParameter(kN,      i);       

  return f;
}
//____________________________________________________________________
inline TF1*
AliLandauGaus::MakeComposite(Double_t c1, 
           Double_t delta, 
           Double_t xi1, 
           Double_t sigma, 
           Double_t c2, 
           Double_t xi2, 
           Double_t xmin, 
           Double_t xmax)
{
  TF1* comp = new TF1("composite", &CompFunc, xmin, xmax, kSigma+1+2);
  comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
        "C#prime", "#xi#prime");
  comp->SetParameters(c1,     // 0 Primary weight 
          delta,  // 1 Primary Delta
          xi1,    // 2 primary Xi
          sigma,  // 3 primary sigma
          c2,     // 5 Secondary weight
          xi2);   // 6 secondary Xi
  return comp;
}
#endif
// Local Variables:
//  mode: C++ 
// End: