TagAndProbe/src/RooCBExGaussShape.cc

0001 #include "PhysicsTools/TagAndProbe/interface/RooCBExGaussShape.h"
0002
0003 ClassImp(RooCBExGaussShape);
0004
0005 RooCBExGaussShape::RooCBExGaussShape(const char* name,
0006                                      const char* title,
0007                                      RooAbsReal& _m,
0008                                      RooAbsReal& _m0,
0009                                      RooAbsReal& _sigma,
0010                                      RooAbsReal& _alpha,
0011                                      RooAbsReal& _n,
0012                                      RooAbsReal& _sigma_2,
0013                                      RooAbsReal& _frac)
0014     : RooAbsPdf(name, title),
0015       m("m", "m", this, _m),
0016       m0(" m0", " m0", this, _m0),
0017       sigma(" sigma", " sigma", this, _sigma),
0018       alpha(" alpha", " alpha", this, _alpha),
0019       n(" n", " n", this, _n),
0020       sigma_2(" sigma_2", " sigma_2", this, _sigma_2),
0021       frac(" frac", " frac", this, _frac) {}
0022
0023 RooCBExGaussShape::RooCBExGaussShape(const RooCBExGaussShape& other, const char* name)
0024     : RooAbsPdf(other, name),
0025       m("m", this, other.m),
0026       m0(" m0", this, other.m0),
0027       sigma(" sigma", this, other.sigma),
0028       alpha(" alpha", this, other.alpha),
0029       n(" n", this, other.n),
0030       sigma_2(" sigma_2", this, other.sigma_2),
0031       frac(" frac", this, other.frac) {}
0032
0033 Double_t RooCBExGaussShape::evaluate() const {
0034   Double_t rval = 0;
0035
0036   Double_t t = (m - m0) / sigma;
0037   Double_t t0 = (m - m0) / sigma_2;
0038   if (alpha < 0) {
0039     t = -t;
0040     t0 = -t0;
0041   }
0042
0043   Double_t absAlpha = fabs((Double_t)alpha);
0044
0045   if (t >= -absAlpha) {
0046     rval = frac * exp(-0.5 * t * t) + (1.0 - frac) * exp(-0.5 * t0 * t0);
0047   } else {
0048     Double_t a = TMath::Power(n / absAlpha, n) * exp(-0.5 * absAlpha * absAlpha);
0049     Double_t b = n / absAlpha - absAlpha;
0050     rval = a / TMath::Power(b - t, n);
0051   }
0052
0053   //std::cout<<"RooCBExGaussShape: m, evaluate= "<<m<<", "<<rval<<std::endl;
0054   return rval;
0055 }