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

// Your description goes here... 

#include "Riostream.h"
#include "AsymBifurcatedCB.h" 
#include "RooAbsReal.h" 
#include "RooAbsCategory.h" 
#include <math.h> 

#include "TMath.h" 

//ClassImp(AsymBifurcatedCB) 

AsymBifurcatedCB::AsymBifurcatedCB(const char *name, const char *title, 
                                    RooAbsReal& _x,
                                   RooAbsReal& _mean,
                                   RooAbsReal& _sigma1, // left
                                   RooAbsReal& _sigma2, // right
                                   RooAbsReal& _a1,
                                   RooAbsReal& _n1,
                                   RooAbsReal& _a2,
                                   RooAbsReal& _n2) :
   RooAbsPdf(name,title), 
   x("x","x",this,_x),
   mean("mean","mean",this,_mean),
   sigma1("sigma1","sigmaL",this,_sigma1),
   sigma2("sigma2","sigmaR",this,_sigma2),
   a1("a1","a1",this,_a1),
   n1("n1","n1",this,_n1),
   a2("a2","a2",this,_a2),
   n2("n2","n2",this,_n2)
 { 
 } 


 AsymBifurcatedCB::AsymBifurcatedCB(const AsymBifurcatedCB& other, const char* name) :  
   RooAbsPdf(other,name), 
   x("x",this,other.x),
   mean("mean",this,other.mean),
   sigma1("sigma1",this,other.sigma1),
   sigma2("sigma2",this,other.sigma2),
   a1("a1",this,other.a1),
   n1("n1",this,other.n1),
   a2("a2",this,other.a2),
   n2("n2",this,other.n2)
 { 
 } 



Double_t AsymBifurcatedCB::evaluate() const {

  double sig = (x  < mean ) ? sigma1 : sigma2;

  double t = (x-mean)/sig;
  double func = 0;

  if (a1 < 0) return 0.;
  if (a2 > 0) return 0.;

  if (n1 < 0) return 0.;
  if (n2 < 0) return 0.;


  double absAlpha = fabs((double)a1);
  double absAlpha2 = fabs((double)a2);

  if (t >= -absAlpha && t <= +absAlpha2) {
    func= exp(-0.5*t*t);
  }
  else if (t < -absAlpha){
    double a = n1/absAlpha;
    double b= n1/absAlpha - absAlpha - t;
    func=exp(-0.5*absAlpha*absAlpha) *pow(a/b,n1);
  }
  else if( t > absAlpha2) {
    double a =  pow(n2/absAlpha2,n2)*exp(-0.5*absAlpha2*absAlpha2);
    double b= n2/absAlpha2 - absAlpha2;
    func= (a/pow(b + t, n2));
  }
  return func;
}