LimitSetting/RooPDF_DSCB_2.cxx at main · IMSA-CMS/LimitSetting · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
/*****************************************************************************
 * Project: RooFit                                                           *
 *                                                                           *
 * This code was autogenerated by RooClassFactory                            *
 *****************************************************************************/

// Your description goes here...

#include "Riostream.h"

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

ClassImp(RooPDF_DSCB_2);

 RooPDF_DSCB_2::RooPDF_DSCB_2(const char *name, const char *title,
                        RooAbsReal& _x,
                        RooAbsReal& _realHiggsMass,
                        RooAbsReal& _alpha_l,
                        RooAbsReal& _alpha_h,
                        RooAbsReal& _n_l,
                        RooAbsReal& _n_h,
                        RooAbsReal& _mean,
                        RooAbsReal& _sigma,
                        RooAbsReal& _norm) :
   RooAbsPdf(name,title),
   x("x","x",this,_x),
   realHiggsMass("realHiggsMass","realHiggsMass",this,_realHiggsMass),
   alpha_l("alpha_l","alpha_l",this,_alpha_l),
   alpha_h("alpha_h","alpha_h",this,_alpha_h),
   n_l("n_l","n_l",this,_n_l),
   n_h("n_h","n_h",this,_n_h),
   mean("mean","mean",this,_mean),
   sigma("sigma","sigma",this,_sigma),
   norm("norm","norm",this,_norm)
 {
 }


 RooPDF_DSCB_2::RooPDF_DSCB_2(const RooPDF_DSCB_2& other, const char* name) :
   RooAbsPdf(other,name),
   x("x",this,other.x),
   realHiggsMass("realHiggsMass",this,other.realHiggsMass),
   alpha_l("alpha_l",this,other.alpha_l),
   alpha_h("alpha_h",this,other.alpha_h),
   n_l("n_l",this,other.n_l),
   n_h("n_h",this,other.n_h),
   mean("mean",this,other.mean),
   sigma("sigma",this,other.sigma),
   norm("norm",this,other.norm)
 {
 }


 Double_t RooPDF_DSCB_2::evaluate() const
 {
  double t = (x - mean) / sigma;
	double result;
	double fact1TLessMinosAlphaL = alpha_l/n_l;
	double fact2TLessMinosAlphaL = (n_l/alpha_l) - alpha_l -t;
	double fact1THhigerAlphaH = alpha_h/n_h;
	double fact2THigherAlphaH = (n_h/alpha_h) - alpha_h +t;

	double root2 = std::pow(2,0.5);
	if (-alpha_l <= t && alpha_h >= t)
	{
		result = exp(-0.5*t*t);
	}
	else if (t < -alpha_l)
	{
		result = exp(-0.5*alpha_l*alpha_l)*pow(fact1TLessMinosAlphaL*fact2TLessMinosAlphaL, -n_l);
	}
	else
	{
		result = exp(-0.5*alpha_h*alpha_h)*pow(fact1THhigerAlphaH*fact2THigherAlphaH, -n_h);
	}

	double lowTailNorm = (n_l/std::abs(alpha_l)) * 1/(n_l - 1) * std::exp(-0.5 * alpha_l * alpha_l);
	double highTailNorm = (n_h/std::abs(alpha_h)) * 1/(n_h - 1) * std::exp(-0.5 * alpha_h * alpha_h);
	double gaussianNormA = erf(std::abs(alpha_l/root2)) + erf(std::abs(alpha_h/root2));
	double gaussianNormB = std::pow(M_PI/2, 0.5) * gaussianNormA;
	double functionNormalization = std::pow(sigma * (gaussianNormB + lowTailNorm + highTailNorm ), -1);
  return  norm * functionNormalization * result;
 }