-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathRooPDF_DSCB.cxx
More file actions
157 lines (130 loc) · 6.81 KB
/
RooPDF_DSCB.cxx
File metadata and controls
157 lines (130 loc) · 6.81 KB
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
/*****************************************************************************
* Project: RooFit *
* *
* This code was autogenerated by RooClassFactory *
*****************************************************************************/
// Class for the Signal Double Sided Crystal Ball functions
#include "Riostream.h"
#include "RooPDF_DSCB.h"
#include "RooAbsReal.h"
#include "RooArgList.h"
#include "RooAbsCategory.h"
#include <math.h>
#include "TMath.h"
#include "RooFormulaVar.h"
ClassImp(RooPDF_DSCB);
RooPDF_DSCB::RooPDF_DSCB(const char *name, const char *title,
RooAbsReal& _x,
RooAbsReal& _realHiggsMass,
RooAbsReal& _branch_ratio_1,
RooAbsReal& _branch_ratio_2,
const std::vector<std::vector<double>>& _signal_params,
bool _multiplyBy2) :
RooAbsPdf(name,title),
x("x","x",this,_x),
realHiggsMass("realHiggsMass","realHiggsMass",this,_realHiggsMass),
branch_ratio_1("branch_ratio_1", "branch_ratio_1", this, _branch_ratio_1),
branch_ratio_2("branch_ratio_2", "branch_ratio_2", this, _branch_ratio_2),
signal_params(_signal_params),
multiplyBy2(_multiplyBy2)
{
//std:: cout << "Main constructor signal params size: " << signal_params.size() << '\n';
}
RooPDF_DSCB::RooPDF_DSCB(const RooPDF_DSCB& other, const char* name) :
RooAbsPdf(other,name),
x("x",this,other.x),
realHiggsMass("realHiggsMass",this,other.realHiggsMass),
branch_ratio_1("branch_ratio_1", this, other.branch_ratio_1),
branch_ratio_2("branch_ratio_2", this, other.branch_ratio_2),
signal_params(other.signal_params),
multiplyBy2(other.multiplyBy2)
{
//std:: cout << "Copy constructor signal params size: " << signal_params.size() << '\n';
}
RooFormulaVar RooPDF_DSCB::signal_norm(std::string channel_name)
{
std::string norm_string = std::to_string(signal_params[6][0]) + "* (@0 - " + std::to_string(signal_params[6][1]) + ")^" + std::to_string(signal_params[6][2]) + " + " + std::to_string(signal_params[6][3]);
RooFormulaVar norm((channel_name + "_norm").c_str(), (channel_name + "_norm").c_str(), norm_string.c_str(), RooArgList(*realHiggsMass.absArg()));
return norm;
}
// In the future, it would be nice if each of these calculated parameters was its own Roo Object that we pass in to handle systematics
Double_t RooPDF_DSCB::evaluate() const
{
//std::cout << "Signal params size: " << signal_params.size() << '\n';
double alpha_l = signal_params[0][0] * std::pow((realHiggsMass - signal_params[0][1]), signal_params[0][2]) + signal_params[0][3];
double alpha_h = signal_params[1][0] * std::pow((realHiggsMass - signal_params[1][1]), signal_params[1][2]) + signal_params[1][3];
double n_l = signal_params[2][0] * std::pow((realHiggsMass - signal_params[2][1]), signal_params[2][2]) + signal_params[2][3];
double n_h = signal_params[3][0] * std::pow((realHiggsMass - signal_params[3][1]), signal_params[3][2]) + signal_params[3][3];
double mean = signal_params[4][0] * std::pow((realHiggsMass - signal_params[4][1]), signal_params[4][2]) + signal_params[4][3];
double sigma = signal_params[5][0] * std::pow((realHiggsMass - signal_params[5][1]), signal_params[5][2]) + signal_params[5][3];
double norm = signal_params[6][0] * std::pow((realHiggsMass - signal_params[6][1]), signal_params[6][2]) + signal_params[6][3];
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);
if (multiplyBy2)
{
// return 2 * branch_ratio_1 * branch_ratio_2 * norm * functionNormalization * result;
return 2 * norm * functionNormalization * result;
}
else
{
// return branch_ratio_1 * branch_ratio_2 * norm * functionNormalization * result;
return norm * functionNormalization * result;
}
}
void RooPDF_DSCB::printParameters()
{
double alpha_l = signal_params[0][0] * std::pow((realHiggsMass - signal_params[0][1]), signal_params[0][2]) + signal_params[0][3];
double alpha_h = signal_params[1][0] * std::pow((realHiggsMass - signal_params[1][1]), signal_params[1][2]) + signal_params[1][3];
double n_l = signal_params[2][0] * std::pow((realHiggsMass - signal_params[2][1]), signal_params[2][2]) + signal_params[2][3];
double n_h = signal_params[3][0] * std::pow((realHiggsMass - signal_params[3][1]), signal_params[3][2]) + signal_params[3][3];
double mean = signal_params[4][0] * std::pow((realHiggsMass - signal_params[4][1]), signal_params[4][2]) + signal_params[4][3];
double sigma = signal_params[5][0] * std::pow((realHiggsMass - signal_params[5][1]), signal_params[5][2]) + signal_params[5][3];
double norm = signal_params[6][0] * std::pow((realHiggsMass - signal_params[6][1]), signal_params[6][2]) + signal_params[6][3];
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);
std::cout<< "FuncNorm: " << functionNormalization << "; norm: " << norm << "; result: " << result << std::endl;
}