-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathRooPDF_DSCB_test.cxx
More file actions
176 lines (141 loc) · 8.23 KB
/
RooPDF_DSCB_test.cxx
File metadata and controls
176 lines (141 loc) · 8.23 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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
/*****************************************************************************
* Project: RooFit *
* *
* This code was autogenerated by RooClassFactory *
*****************************************************************************/
// Your description goes here...
#include "Riostream.h"
#include "RooPDF_DSCB_test.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include <math.h>
#include "TMath.h"
#include "RooFormulaVar.h"
// Indices in signal_params for each parameter
const int ALPHA_L = 1;
const int ALPHA_H = 0;
const int N_L = 5;
const int N_H = 4;
const int MEAN = 2;
const int SIGMA = 3;
const int NORM = 6;
ClassImp(RooPDF_DSCB_test);
RooPDF_DSCB_test::RooPDF_DSCB_test(const char *name, const char *title,
RooAbsReal& _x,
RooAbsReal& _realHiggsMass,
RooAbsReal& _branch_ratio_1,
RooAbsReal& _branch_ratio_2,
RooAbsReal& _norm_Systematic,
RooAbsReal& _shape_Systematic,
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),
norm_Systematic("norm_Systematic","norm_Systematic",this,_norm_Systematic),
shape_Systematic("shape_Systematic","shape_Systematic",this,_shape_Systematic),
signal_params(_signal_params),
multiplyBy2(_multiplyBy2)
{
}
RooPDF_DSCB_test::RooPDF_DSCB_test(const RooPDF_DSCB_test& 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),
norm_Systematic("norm_Systematic",this,other.norm_Systematic),
shape_Systematic("shape_Systematic",this,other.shape_Systematic),
signal_params(other.signal_params),
multiplyBy2(other.multiplyBy2)
{
}
RooFormulaVar RooPDF_DSCB_test::signal_norm(std::string channel_name)
{
// std::cout << "signal_params[NORM][0] = " << signal_params[NORM][0] << "\n";
// std::cout << "signal_params[NORM][1] = " << signal_params[NORM][1] << "\n";
// std::cout << "signal_params[NORM][2] = " << signal_params[NORM][2] << "\n";
// std::cout << "signal_params[NORM][3] = " << signal_params[NORM][3] << "\n";
//std::string norm_string = "@1 * " + std::to_string(signal_params[NORM][0]) + "* (@0 - " + std::to_string(signal_params[NORM][1]) + ")^" + std::to_string(signal_params[NORM][2]) + " + " + std::to_string(signal_params[NORM][3]);
std::string norm_string = "@1 * " + std::to_string(signal_params[NORM][0]) + "* (@0 - " + std::to_string(signal_params[NORM][1]) + ")^" + std::to_string(signal_params[NORM][2]);
RooFormulaVar norm((channel_name + "_norm").c_str(), (channel_name + "_norm").c_str(), norm_string.c_str(), RooArgList(*realHiggsMass.absArg(), *norm_Systematic.absArg()));
return norm;
}
double RooPDF_DSCB_test::evaluatePowerLaw(double p0, double p1, double p2) const
{
return p0 * std::pow((realHiggsMass - p1), p2);
}
Double_t RooPDF_DSCB_test::evaluate() const
{
// double alpha_l = signal_params[ALPHA_L][0] * std::pow((realHiggsMass - signal_params[ALPHA_L][1]), signal_params[ALPHA_L][2]) + signal_params[ALPHA_L][3];
// double alpha_h = signal_params[ALPHA_H][0] * std::pow((realHiggsMass - signal_params[ALPHA_H][1]), signal_params[ALPHA_H][2]) + signal_params[ALPHA_H][3];
// double n_l = signal_params[N_L][0] * std::pow((realHiggsMass - signal_params[N_L][1]), signal_params[N_L][2]) + signal_params[N_L][3];
// double n_h = signal_params[N_H][0] * std::pow((realHiggsMass - signal_params[N_H][1]), signal_params[N_H][2]) + signal_params[N_H][3];
// double mean = shape_Systematic * signal_params[MEAN][0] * std::pow((realHiggsMass - signal_params[MEAN][1]), signal_params[MEAN][2]) + signal_params[MEAN][3];
// double sigma = signal_params[SIGMA][0] * std::pow((realHiggsMass - signal_params[SIGMA][1]), signal_params[SIGMA][2]) + signal_params[SIGMA][3];
// double norm = signal_params[NORM][0] * std::pow((realHiggsMass - signal_params[NORM][1]), signal_params[NORM][2]) + signal_params[NORM][3];
// double t = (x - mean) / sigma;
double alpha_l = evaluatePowerLaw(signal_params[ALPHA_L][0], signal_params[ALPHA_L][1], signal_params[ALPHA_L][2]);
double alpha_h = evaluatePowerLaw(signal_params[ALPHA_H][0], signal_params[ALPHA_H][1], signal_params[ALPHA_H][2]);
double n_l = evaluatePowerLaw(signal_params[N_L][0], signal_params[N_L][1], signal_params[N_L][2]);
double n_h = evaluatePowerLaw(signal_params[N_H][0], signal_params[N_H][1], signal_params[N_H][2]);
double mean = evaluatePowerLaw(signal_params[MEAN][0], signal_params[MEAN][1], signal_params[MEAN][2]);
double sigma = evaluatePowerLaw(signal_params[SIGMA][0], signal_params[SIGMA][1], signal_params[SIGMA][2]);
double norm = evaluatePowerLaw(signal_params[NORM][0], signal_params[NORM][1], signal_params[NORM][2]);
double t = (x - mean) / sigma;
//std::cout << GetName() << " parameters: " << "\n";
//std::cout << "x: " << x << "\n";
// std::cout << "Mean: " << mean << "\n";
// std::cout << "Norm: " << norm << "\n";
// std::cout << "n_l: " << n_l << "\n";
// std::cout << "n_h: " << n_h << "\n";
// std::cout << "alpha_l: " << alpha_l << "\n";
// std::cout << "alpha_h: " << alpha_h << "\n";
// std::cout << "sigma: " << sigma << "\n";
// std::cout << "Mean parameters: " << signal_params[MEAN][0] << ", " << signal_params[MEAN][1] << ", " << signal_params[MEAN][2] << "\n";
// std::cout << "Norm parameters: " << signal_params[NORM][0] << ", " << signal_params[NORM][1] << ", " << signal_params[NORM][2] << "\n";
// std::cout << "n_l parameters: " << signal_params[N_L][0] << ", " << signal_params[N_L][1] << ", " << signal_params[N_L][2] << "\n";
// std::cout << "n_h parameters: " << signal_params[N_H][0] << ", " << signal_params[N_H][1] << ", " << signal_params[N_H][2] << "\n";
// std::cout << "alpha_l parameters: " << signal_params[ALPHA_L][0] << ", " << signal_params[ALPHA_L][1] << ", " << signal_params[ALPHA_L][2] << "\n";
// std::cout << "alpha_h parameters: " << signal_params[ALPHA_H][0] << ", " << signal_params[ALPHA_H][1] << ", " << signal_params[ALPHA_H][2] << "\n";
// std::cout << "sigma parameters: " << signal_params[SIGMA][0] << ", " << signal_params[SIGMA][1] << ", " << signal_params[SIGMA][2] << "\n";
// std::cout << "realHiggsMass: " << realHiggsMass << "\n \n \n";
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 << "calculated norm: " << norm << "\n";
// std::cout << "functionNormalization: " << functionNormalization << "\n";
// std::cout << "result: " << result << "\n \n \n";
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;
}
}