RooGausDExp.cxx

/*****************************************************************************
 * Project: RooFit                                                           *
 *                                                                           *
 * This code was autogenerated by RooClassFactory                            *
 *****************************************************************************/

#include "RooGausDExp.h"

#include <Riostream.h>
#include <TMath.h>

ClassImp(RooGausDExp)

RooGausDExp::RooGausDExp(const char *name, const char *title,
                       RooAbsReal& _x,
                       RooAbsReal& _mu,
                       RooAbsReal& _sig,
                       RooAbsReal& _tau0,
                       RooAbsReal& _tau1)
: RooAbsPdf(name,title)
, x("x","x",this,_x)
, mu("mu","mu",this,_mu)
, sig("sig","sig",this,_sig)
, tau0("tau0","tau",this,_tau0)
, tau1("tau1","tau",this,_tau1) {
}


RooGausDExp::RooGausDExp(const RooGausDExp& other, const char* name)
: RooAbsPdf(other,name)
, x("x",this,other.x)
, mu("mu",this,other.mu)
, sig("sig",this,other.sig)
, tau0("tau0",this,other.tau0)
, tau1("tau1",this,other.tau1) {
}



Double_t RooGausDExp::evaluate() const {
  double u = (x - mu) / sig;
  if (u < tau0)
    return TMath::Exp(-tau0 * (u - 0.5 * tau0));
  else if (u <= tau1)
    return TMath::Exp(-u * u * 0.5);
  else
    return TMath::Exp(-tau1 * (u - 0.5 * tau1));
}

#ifdef _AN_INT_
Int_t RooGausDExp::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char*) const {
  if (matchArgs(allVars,analVars,x)) return 1 ;
  return 0 ;
}

Double_t RooGausDExp::analyticalIntegral(Int_t code, const char* r) const
{
  double umin = (x.min(r) - mu) / sig;
  double umax = (x.max(r) - mu) / sig;
  R__ASSERT(code==1);
  double integral = 0.;
  
  if (umin < tau0) {
    integral -= IntExp(umin, tau0);
    if (umax <= tau0)
      integral += IntExp(umax, tau0);
    else {
      integral += IntExp(tau0, tau0) - IntGaus(tau0);
      if (umax <= tau1)
        integral += IntGaus(umax);
      else
        integral += IntGaus(tau1) - IntExp(tau1, tau1) + IntExp(umax, tau1);
    }
  } else if (umin < tau1) {
    integral -= IntGaus(umin);
    if (umax < tau1)
      integral += IntGaus(umax);
    else
      integral += IntGaus(tau1) - IntExp(tau1, tau1) + IntExp(umax, tau1);
  } else {
    integral = IntExp(umax, tau1) - IntExp(umin, tau1);
  }
  
  return sig * integral;
}
#endif

double RooGausDExp::IntExp(double x,double tau) const {
  return -1. * TMath::Exp(tau * (0.5 * tau - x)) / tau;
}

double RooGausDExp::IntGaus(double x) const {
  static const double rootPiBy2 = TMath::Sqrt(TMath::PiOver2());
  return rootPiBy2 * (TMath::Erf(x / TMath::Sqrt2()));
}