AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliLandauGaus.h
Go to the documentation of this file.
1 #ifndef ALILANDAUGAUS_H
2 #define ALILANDAUGAUS_H
3 
13 #include <TObject.h>
14 #include <TF1.h>
15 #include <TMath.h>
16 
111 {
112 public:
114  enum {
115  kC = 0,
117  kXi,
120  kN,
122  };
125  //__________________________________________________________________
130  //------------------------------------------------------------------
137  static Double_t MPShift() { return -0.22278298; }
143  static Double_t SigmaShiftC() { return 0.5446; }
149  static Double_t SigmaShiftP() { return 0.5746; }
155  static Double_t InvSq2Pi() { return 1. / TMath::Sqrt(TMath::TwoPi()); }
160  static Double_t NSigma() { return 5; }
164  static Int_t NSteps() { return 100; }
165  /* @} */
166 
167  //__________________________________________________________________
172  //------------------------------------------------------------------
189  static Double_t Fl(Double_t x, Double_t delta, Double_t xi);
190  //------------------------------------------------------------------
211  static Double_t F(Double_t x, Double_t delta, Double_t xi,
212  Double_t sigma, Double_t sigma_n);
213  //------------------------------------------------------------------
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  //------------------------------------------------------------------
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  //------------------------------------------------------------------
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);
302  static void IPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma);
310  static Bool_t EnableSigmaShift(Short_t val=-1);
327  /* @} */
328 
329 
330  //__________________________________________________________________
335  //------------------------------------------------------------------
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  //------------------------------------------------------------------
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  //------------------------------------------------------------------
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);
410  static TF1* MakeComposite(Double_t c1,
411  Double_t delta,
412  Double_t xi1,
413  Double_t sigma,
414  Double_t c2,
415  Double_t xi2,
416  Double_t xmin,
417  Double_t xmax);
418 
419  //------------------------------------------------------------------
427  static Color_t GetIColor(Int_t i);
428  //------------------------------------------------------------------
439  static Double_t F1Func(Double_t* xp, Double_t* pp);
440  //------------------------------------------------------------------
451  static Double_t FnFunc(Double_t* xp, Double_t* pp);
452  //------------------------------------------------------------------
463  static Double_t FiFunc(Double_t* xp, Double_t* pp);
464  //------------------------------------------------------------------
475  static Double_t CompFunc(Double_t* xp, Double_t* pp);
476  /* @} */
477 };
478 //____________________________________________________________________
479 inline Bool_t
481 {
482  static Bool_t enabled = true;
483  if (val >= 0) enabled = val == 1;
484  return enabled;
485 }
486 //____________________________________________________________________
487 inline void
489 {
490 #ifndef NO_SIGMA_SHIFT
491  Double_t dDelta = EnableSigmaShift() ? SigmaShift(i, xi, sigma) : 0;
492 #else
493  // This is for testing
494  Double_t dDelta = 0;
495 #endif
496  if (i == 1) {
497  delta -= dDelta;
498  return;
499  }
500 
501  delta = i * (delta + xi * TMath::Log(i)) - dDelta;
502  xi = i * xi;
503  sigma = TMath::Sqrt(Double_t(i)) * sigma;
504 }
505 //____________________________________________________________________
506 inline Double_t
508 {
509  if (xi <= 0) return 0;
510  if (sigma <= 0) return 0;
511  const Double_t c = SigmaShiftC();
512  const Double_t p = SigmaShiftP();
513  const Double_t u = sigma / xi;
514  const Double_t q = p*u*TMath::Sqrt(u);
515  if (q > 100) return 0;
516  return c * sigma / TMath::Power(1+1./i, q);
517 }
518 //____________________________________________________________________
519 inline Double_t
521 {
522  Double_t deltaP = delta - xi * MPShift();
523  return TMath::Landau(x, deltaP, xi, true);
524 }
525 //____________________________________________________________________
526 inline Double_t
528  Double_t sigma, Double_t sigmaN)
529 {
530  if (xi <= 0) return 0;
531 
532  const Int_t nSteps = NSteps();
533  const Double_t nSigma = NSigma();
534  const Double_t deltaP = delta; // - sigma * sigmaShift; // + sigma * mpshift;
535  const Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
536  const Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
537  const Double_t xlow = x - nSigma * sigma1;
538  const Double_t xhigh = x + nSigma * sigma1;
539  const Double_t step = (xhigh - xlow) / nSteps;
540  Double_t sum = 0;
541 
542  for (Int_t i = 0; i <= nSteps/2; i++) {
543  const Double_t x1 = xlow + (i - .5) * step;
544  const Double_t x2 = xhigh - (i - .5) * step;
545  sum += Fl(x1, deltaP, xi) * TMath::Gaus(x, x1, sigma1);
546  sum += Fl(x2, deltaP, xi) * TMath::Gaus(x, x2, sigma1);
547  }
548  return step * sum * InvSq2Pi() / sigma1;
549 }
550 
551 //____________________________________________________________________
552 inline Double_t
554  Double_t sigma, Double_t sigmaN, Int_t i)
555 {
556  Double_t deltaI = delta;
557  Double_t xiI = xi;
558  Double_t sigmaI = sigma;
559  IPars(i, deltaI, xiI, sigmaI);
560  if (sigmaI < 1e-10)
561  // Fall back to landau
562  return Fl(x, deltaI, xiI);
563 
564  return F(x, deltaI, xiI, sigmaI, sigmaN);
565 }
566 //____________________________________________________________________
567 inline Double_t
569  Double_t sigma, Double_t sigmaN, Int_t n,
570  const Double_t* a)
571 {
572  Double_t result = Fi(x, delta, xi, sigma, sigmaN, 1);
573  for (Int_t i = 2; i <= n; i++)
574  result += a[i-2] * Fi(x,delta,xi,sigma,sigmaN,i);
575  return result;
576 }
577 
578 //____________________________________________________________________
579 inline Double_t
581  UShort_t par, Double_t dPar,
582  Double_t delta, Double_t xi,
583  Double_t sigma, Double_t sigmaN,
584  Int_t i)
585 {
586  if (dPar == 0) return 0;
587  Double_t dp = dPar;
588  Double_t d2 = dPar / 2;
589  Double_t deltaI = i * (delta + xi * TMath::Log(i));
590  Double_t xiI = i * xi;
591  Double_t si = TMath::Sqrt(Double_t(i));
592  Double_t sigmaI = si*sigma;
593  Double_t y1 = 0;
594  Double_t y2 = 0;
595  Double_t y3 = 0;
596  Double_t y4 = 0;
597  switch (par) {
598  case 0:
599  y1 = Fi(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
600  y2 = Fi(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
601  y3 = Fi(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
602  y4 = Fi(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
603  break;
604  case 1:
605  y1 = Fi(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
606  y2 = Fi(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
607  y3 = Fi(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
608  y4 = Fi(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
609  break;
610  case 2:
611  y1 = Fi(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
612  y2 = Fi(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
613  y3 = Fi(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
614  y4 = Fi(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
615  break;
616  case 3:
617  y1 = Fi(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
618  y2 = Fi(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
619  y3 = Fi(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
620  y4 = Fi(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
621  break;
622  default:
623  return 0;
624  }
625 
626  Double_t d0 = y1 - y4;
627  Double_t d1 = 2 * (y2 - y3);
628 
629  Double_t g = 1/(2*dp) * (4*d1 - d0) / 3;
630 
631  return g;
632 }
633 
634 
635 //____________________________________________________________________
636 inline Color_t
638 {
639  const Int_t kColors[] = { kRed+1,
640  kPink+3,
641  kMagenta+2,
642  kViolet+2,
643  kBlue+1,
644  kAzure+3,
645  kCyan+1,
646  kTeal+2,
647  kGreen+2,
648  kSpring+3,
649  kYellow+2,
650  kOrange+2 };
651  return kColors[((i-1) % 12)];
652 }
653 //____________________________________________________________________
654 inline Double_t
656 {
657  Double_t x = xp[0];
658  Double_t constant = pp[kC];
659  Double_t delta = pp[kDelta];
660  Double_t xi = pp[kXi];
661  Double_t sigma = pp[kSigma];
662  Double_t sigmaN = pp[kSigmaN];
663 
664  return constant * F(x, delta, xi, sigma, sigmaN);
665 }
666 //____________________________________________________________________
667 inline Double_t
669 {
670  Double_t x = xp[0];
671  Double_t constant = pp[kC];
672  Double_t delta = pp[kDelta];
673  Double_t xi = pp[kXi];
674  Double_t sigma = pp[kSigma];
675  Double_t sigmaN = pp[kSigmaN];
676  Int_t i = Int_t(pp[kN]);
677 
678  return constant * Fi(x, delta, xi, sigma, sigmaN, i);
679 }
680 //____________________________________________________________________
681 inline Double_t
683 {
684  Double_t x = xp[0];
685  Double_t constant = pp[kC];
686  Double_t delta = pp[kDelta];
687  Double_t xi = pp[kXi];
688  Double_t sigma = pp[kSigma];
689  Double_t sigmaN = pp[kSigmaN];
690  Int_t n = Int_t(pp[kN]);
691  Double_t* a = &(pp[kA]);
692 
693  return constant * Fn(x, delta, xi, sigma, sigmaN, n, a);
694 }
695 //____________________________________________________________________
696 inline Double_t
698 {
699  Double_t x = xp[0];
700  Double_t cP = pp[kC];
701  Double_t deltaP = pp[kDelta];
702  Double_t xiP = pp[kXi];
703  Double_t sigmaP = pp[kSigma];
704  Double_t cS = pp[kSigma+1];
705  Double_t deltaS = deltaP; // pp[kSigma+2];
706  Double_t xiS = pp[kSigma+2/*3*/];
707  Double_t sigmaS = sigmaP; // pp[kSigma+4];
708 
709  return (cP * F(x,deltaP,xiP,sigmaP,0) +
710  cS * F(x,deltaS,xiS,sigmaS,0));
711 
712 }
713 //____________________________________________________________________
714 inline TF1*
716  Double_t delta, Double_t xi,
717  Double_t sigma, Double_t sigmaN,
718  Double_t xmin, Double_t xmax)
719 {
720  // Define the function to fit
721  TF1* f = new TF1("landau1", &F1Func, xmin,xmax,kSigmaN+1);
722 
723  // Set initial guesses, parameter names, and limits
724  f->SetParameters(c,delta,xi,sigma,sigmaN);
725  f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
726  f->SetNpx(500);
727  f->SetLineColor(GetIColor(1));
728 
729  return f;
730 }
731 
732 //____________________________________________________________________
733 inline TF1*
735  Double_t delta, Double_t xi,
736  Double_t sigma, Double_t sigmaN, Int_t n,
737  const Double_t* a,
738  Double_t xmin, Double_t xmax)
739 {
740  Int_t npar = kN+n;
741  TF1* f = new TF1(Form("nlandau%d", n), &FnFunc,xmin,xmax,npar);
742  f->SetLineColor(GetIColor(n));
743  f->SetLineWidth(2);
744  f->SetNpx(500);
745  f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
746  f->SetParameter(kC, c);
747  f->SetParameter(kDelta, delta);
748  f->SetParameter(kXi, xi);
749  f->SetParameter(kSigma, sigma);
750  f->SetParameter(kSigmaN, sigmaN);
751  f->FixParameter(kN, n);
752  for (UShort_t i = 2; i <= n; i++) {
753  f->SetParameter(kA+i-2, a[i-2]);
754  f->SetParName(kA+i-2, Form("a_{%d}", i));
755  }
756  return f;
757 }
758 //____________________________________________________________________
759 inline TF1*
761  Double_t delta, Double_t xi,
762  Double_t sigma, Double_t sigmaN, Int_t i,
763  Double_t xmin, Double_t xmax)
764 {
765  Int_t npar = kN+1;
766  TF1* f = new TF1(Form("ilandau%d", i), &FiFunc,xmin,xmax,npar);
767  f->SetLineColor(GetIColor(i));
768  f->SetLineWidth(1);
769  f->SetNpx(500);
770  f->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
771  f->SetParameter(kC, c);
772  f->SetParameter(kDelta, delta);
773  f->SetParameter(kXi, xi);
774  f->SetParameter(kSigma, sigma);
775  f->SetParameter(kSigmaN, sigmaN);
776  f->FixParameter(kN, i);
777 
778  return f;
779 }
780 //____________________________________________________________________
781 inline TF1*
783  Double_t delta,
784  Double_t xi1,
785  Double_t sigma,
786  Double_t c2,
787  Double_t xi2,
788  Double_t xmin,
789  Double_t xmax)
790 {
791  TF1* comp = new TF1("composite", &CompFunc, xmin, xmax, kSigma+1+2);
792  comp->SetParNames("C", "#Delta_{p}", "#xi", "#sigma",
793  "C#prime", "#xi#prime");
794  comp->SetParameters(c1, // 0 Primary weight
795  delta, // 1 Primary Delta
796  xi1, // 2 primary Xi
797  sigma, // 3 primary sigma
798  c2, // 5 Secondary weight
799  xi2); // 6 secondary Xi
800  return comp;
801 }
802 #endif
803 // Local Variables:
804 // mode: C++
805 // End:
806 
807 
808 
static Double_t SigmaShift(Int_t i, Double_t xi, Double_t sigma)
static Double_t SigmaShiftP()
double Double_t
Definition: External.C:58
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 Double_t SigmaShiftC()
static Double_t NSigma()
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)
TCanvas * c
Definition: TestFitELoss.C:172
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 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 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)
Double_t * sigma
static Double_t InvSq2Pi()
int Int_t
Definition: External.C:63
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 Double_t F1Func(Double_t *xp, Double_t *pp)
static Double_t FnFunc(Double_t *xp, Double_t *pp)
static Color_t GetIColor(Int_t i)
static Double_t FiFunc(Double_t *xp, Double_t *pp)
short Short_t
Definition: External.C:23
static Double_t MPShift()
static void IPars(Int_t i, Double_t &delta, Double_t &xi, Double_t &sigma)
static Double_t F(Double_t x, Double_t delta, Double_t xi, Double_t sigma, Double_t sigma_n)
static Double_t Fl(Double_t x, Double_t delta, Double_t xi)
static Int_t NSteps()
static Bool_t EnableSigmaShift(Short_t val=-1)
static Double_t Fi(Double_t x, Double_t delta, Double_t xi, Double_t sigma, Double_t sigma_n, Int_t i)
unsigned short UShort_t
Definition: External.C:28
bool Bool_t
Definition: External.C:53
static Double_t CompFunc(Double_t *xp, Double_t *pp)