CrystalBallEfficiency.cxx

/*****************************************************************************
 * Project: RooFit                                                           *
 *                                                                           *
 * This code was autogenerated by RooClassFactory                            *
 *****************************************************************************/
#include "Riostream.h"

#include "CrystalBallEfficiency.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include <math.h>
#include "TMath.h"

ClassImp(CrystalBallEfficiency)

 CrystalBallEfficiency::CrystalBallEfficiency(const char *name, const char *title,
                        RooAbsReal& _m,
                        RooAbsReal& _m0,
                        RooAbsReal& _sigma,
                        RooAbsReal& _alpha,
                        RooAbsReal& _n,
                        RooAbsReal& _norm) :
   RooAbsReal(name,title),
   m("m","m",this,_m),
   m0("m0","m0",this,_m0),
   sigma("sigma","sigma",this,_sigma),
   alpha("alpha","alpha",this,_alpha),
   n("n","n",this,_n),
   norm("norm","norm",this,_norm)
 {
 }


 CrystalBallEfficiency::CrystalBallEfficiency(const CrystalBallEfficiency& other, const char* name) :
   RooAbsReal(other,name),
   m("m",this,other.m),
   m0("m0",this,other.m0),
   sigma("sigma",this,other.sigma),
   alpha("alpha",this,other.alpha),
   n("n",this,other.n),
   norm("norm",this,other.norm)
 {
 }



 Double_t CrystalBallEfficiency::evaluate() const
 {
   double sqrtPiOver2 = std::sqrt(TMath::PiOver2());
   double sqrt2       = std::sqrt(2.);
   double sig         = std::abs(sigma);
   double t           = (m - m0)/sig * alpha / std::abs(alpha);
   double absAlpha    = std::abs(alpha/sig);
   double a           = TMath::Power(n/absAlpha, n) * TMath::Exp(-0.5 * absAlpha * absAlpha);
   double b           = absAlpha - n/absAlpha;
   double arg         = absAlpha / sqrt2;

   double ApproxErf = 0.;
   if (arg >  5.) {
      ApproxErf =  1.;
   } else if (arg < -5.) {
      ApproxErf = -1.;
   } else {
      ApproxErf = TMath::Erf(arg);
   }

   double leftArea    = (1. + ApproxErf) * sqrtPiOver2;
   double rightArea   = ( a * 1. / TMath::Power(absAlpha-b, n-1) ) / (n - 1.);
   double area        = leftArea + rightArea;

   if (t <= absAlpha) {
      arg = t / sqrt2;
      if (arg >  5.) {
        ApproxErf =  1.;
      } else if (arg < -5.) {
        ApproxErf = -1.;
      } else {
        ApproxErf = TMath::Erf(arg);
      }
      return norm * (1. + ApproxErf) * sqrtPiOver2 / area;
   } else {
      return norm * (leftArea + a * (1./TMath::Power(t-b,n-1) - \
                                     1./TMath::Power(absAlpha - b,n-1)) / (1 - n)) / area;
   }
 }