hpstr
The Heavy Photon Search Toolkit for Reconstruction (hpstr) provides an interface to physics data from the HPS experiment saved in the LCIO format and converts it into an ROOT based format.
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
SimpEquations.cxx
Go to the documentation of this file.
1#include "SimpEquations.h"
2#include <fstream>
3#include <iostream>
4#include <cmath>
5
7 year_ = year;
8
9 std::cout << "[SimpEquations]::Parameters " << std::endl;
10 std::cout << "mass ratio Ap:Pi_D = " << mass_ratio_Ap_to_Pid_ << std::endl;
11 std::cout << "mass ratio Ap:V_D = " << mass_ratio_Ap_to_Vd_ << std::endl;
12 std::cout << "ratio of dark pion mass to dark pion decay constant: " << ratio_mPi_to_fPi_ << std::endl;
13 std::cout << "Dark Vector decay lepton mass: " << m_l_ << std::endl;
14 std::cout << "Alpha Dark: " << alpha_dark_ << std::endl;
15
16}
17
18SimpEquations::SimpEquations(int year, const std::string paramsConfigFile){
19 year_ = year;
20 loadParametersConfig(paramsConfigFile);
21
22 std::cout << "[SimpEquations]::Parameters Loaded From Json Config " << std::endl;
23 std::cout << "mass ratio Ap:Pi_D = " << mass_ratio_Ap_to_Pid_ << std::endl;
24 std::cout << "mass ratio Ap:V_D = " << mass_ratio_Ap_to_Vd_ << std::endl;
25 std::cout << "ratio of dark pion mass to dark pion decay constant: " << ratio_mPi_to_fPi_ << std::endl;
26 std::cout << "Dark Vector decay lepton mass: " << m_l_ << std::endl;
27 std::cout << "Alpha Dark: " << alpha_dark_ << std::endl;
28
29}
30
31void SimpEquations::loadParametersConfig(const std::string paramsConfigFile){
32
33 //Read SIMP Model parameters from json
34 std::ifstream i_file(paramsConfigFile);
35 i_file >> params_config_;
36 for(auto& el : params_config_.items())
37 std::cout << el.key() << " : " << el.value() << "\n";
38 i_file.close();
39
40 //Set parameter values. All parameters must be specificed in json config
41 try
42 {
43 for(auto param : params_config_.items()){
44 mass_ratio_Ap_to_Vd_ = param.value().at("mass_ratio_Ap_to_Vd");
45 mass_ratio_Ap_to_Pid_ = param.value().at("mass_ratio_Ap_to_Pid");
46 ratio_mPi_to_fPi_ = param.value().at("ratio_mPi_to_fPi");
47 m_l_ = param.value().at("lepton_mass");
48 alpha_dark_ = param.value().at("alpha_dark");
49 }
50 }
51 catch(...)
52 {
53 std::cout << "[SimpEquations]::ERROR::REQUIRED PARAMETER CONFIGURATIONS NOT FOUND IN "
54 << paramsConfigFile << std::endl;
55 }
56}
57
58double SimpEquations::rate_2pi(double m_Ap, double m_pi, double m_V, double alpha_dark){
59 double coeff = (2.0*alpha_dark/3.0) * m_Ap;
60 double pow1 = std::pow((1-(4*m_pi*m_pi/(m_Ap*m_Ap))),3/2.);
61 double pow2 = std::pow(((m_V*m_V)/((m_Ap*m_Ap)-(m_V*m_V))),2);
62 return coeff*pow1*pow2;
63}
64
65double SimpEquations::rate_Vrho_pi(double m_Ap, double m_pi, double m_V,
66 double alpha_dark, double f_pi){
67 double x = m_pi / m_Ap;
68 double y = m_V / m_Ap;
69 double Tv = 3.0/4.0;
70 double coeff = alpha_dark*Tv/(192.*std::pow(M_PI,4));
71 return coeff * std::pow((m_Ap/m_pi),2) * std::pow(m_V/m_pi,2) * std::pow((m_pi/f_pi),4) * m_Ap*std::pow(Beta(x,y),3./2.);
72}
73
74double SimpEquations::rate_Vphi_pi(double m_Ap, double m_pi, double m_V,
75 double alpha_dark, double f_pi){
76 double x = m_pi / m_Ap;
77 double y = m_V / m_Ap;
78 double Tv = 3.0/2.0;
79 double coeff = alpha_dark*Tv/(192.*std::pow(M_PI,4));
80 return coeff * std::pow((m_Ap/m_pi),2) * std::pow(m_V/m_pi,2) * std::pow((m_pi/f_pi),4) * m_Ap*std::pow(Beta(x,y),3./2.);
81}
82double SimpEquations::rate_Vcharged_pi(double m_Ap, double m_pi, double m_V,
83 double alpha_dark, double f_pi){
84 double x = m_pi / m_Ap;
85 double y = m_V / m_Ap;
86 double Tv = 18.0-((3.0/2.0) +(3.0/4.0)) ;
87 double coeff = alpha_dark*Tv/(192.*std::pow(M_PI,4));
88 return coeff * std::pow((m_Ap/m_pi),2) * std::pow(m_V/m_pi,2) * std::pow((m_pi/f_pi),4) * m_Ap*std::pow(Beta(x,y),3./2.);
89}
90
91double SimpEquations::br_2pi(double m_Ap, double m_pi, double m_V, double alpha_dark,
92 double f_pi){
93 double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
94 rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
95 rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
96 rate_2pi(m_Ap, m_pi, m_V, alpha_dark);
97 if(m_Ap > 2.0*m_V)
98 total_rate + rate_2V(m_Ap, m_V, alpha_dark);
99 return rate_2pi(m_Ap, m_pi, m_V, alpha_dark)/total_rate;
100}
101
102double SimpEquations::br_Vrho_pi(double m_Ap, double m_pi, double m_V, double alpha_dark,
103 double f_pi){
104 double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
105 rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
106 rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
107 rate_2pi(m_Ap, m_pi, m_V, alpha_dark);
108 if(m_Ap > 2.0*m_V)
109 total_rate + rate_2V(m_Ap, m_V, alpha_dark);
110 return rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi)/total_rate;
111}
112
113double SimpEquations::br_Vphi_pi(double m_Ap, double m_pi, double m_V, double alpha_dark,
114 double f_pi){
115 double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
116 rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
117 rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
118 rate_2pi(m_Ap, m_pi, m_V, alpha_dark);
119 if(m_Ap > 2.0*m_V)
120 total_rate + rate_2V(m_Ap, m_V, alpha_dark);
121 return rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi)/total_rate;
122}
123
124double SimpEquations::br_Vcharged_pi(double m_Ap, double m_pi, double m_V,
125 double alpha_dark, double f_pi){
126 double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
127 rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
128 rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
129 rate_2pi(m_Ap, m_pi, m_V, alpha_dark);
130 if(m_Ap > 2.0*m_V)
131 total_rate + rate_2V(m_Ap, m_V, alpha_dark);
132 return rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark,f_pi)/total_rate;
133}
134
135double SimpEquations::br_2V(double m_Ap,double m_pi,double m_V,double alpha_dark,double f_pi,double rho,double phi){
136 double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
137 rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
138 rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) +
139 rate_2pi(m_Ap, m_pi, m_V, alpha_dark);
140 if(m_Ap > 2.0*m_V){
141 total_rate + rate_2V(m_Ap, m_V, alpha_dark);
142 return 0.0;
143 }
144 return rate_2V(m_Ap,m_V,alpha_dark)/total_rate;
145}
146
147double SimpEquations::Tv(bool rho,bool phi){
148 if(rho) return 3./4.;
149 else if(phi) return 3./2.;
150 else return 18.;
151}
152
153double SimpEquations::Beta(double x,double y){
154 return (1+std::pow(y,2)-std::pow(x,2)-2*y)*(1+std::pow(y,2)-std::pow(x,2)+2*y);
155}
156
157double SimpEquations::rate_2V(double m_Ap,double m_V,double alpha_dark){
158 double r = m_V/m_Ap;
159 return alpha_dark/6. * m_Ap * f(r);
160}
161
162double SimpEquations::f(double r){
163 double num = 1 + 16*std::pow(r,2) - 68*std::pow(r,4) - 48*std::pow(r,6);
164 double den = std::pow(1-std::pow(r,2),2);
165 return num/den * std::pow((1-4*std::pow(r,2)),0.5);
166}
167
168double SimpEquations::rate_2l(double m_Ap,double m_pi,double m_V,double eps,double alpha_dark,double f_pi,double m_l,bool rho){
169 double alpha = 1./137.0;
170 double coeff = 16*M_PI*alpha_dark*alpha*std::pow(eps,2)*std::pow(f_pi,2)/(3*std::pow(m_V,2));
171 double term1 = std::pow((std::pow(m_V,2)/(std::pow(m_Ap,2) - std::pow(m_V,2))),2);
172 double term2 = std::pow((1-(4*std::pow(m_l,2)/std::pow(m_V,2))),0.5);
173 double term3 = 1+(2*std::pow(m_l,2)/std::pow(m_V,2));
174 int constant = 1;
175 if(rho) constant = 2;
176 return coeff * term1 * term2 * term3 * m_V * constant;
177}
178
179double SimpEquations::getCtau(double m_Ap,double m_pi,double m_V,double eps,double alpha_dark,double f_pi,double m_l,bool rho){
180 double c = 3.00e11; //mm/s
181 double hbar = 6.58e-22; //MeV*sec
182 double rate = rate_2l(m_Ap,m_pi,m_V,eps,alpha_dark,f_pi,m_l,rho); //MeV
183 double tau = hbar/rate;
184 double ctau = c*tau;
185 return ctau;
186}
187
188double SimpEquations::gamma(double m_V,double E_V){
189 double gamma = E_V/m_V;
190 return gamma;
191}
192
193//Calculate expected signal by passing radFrac, radAcc, and dNdm
194double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho,
195 double E_V, TEfficiency* effCalc_h, double dNdm, double radFrac, double radAcc, double target_pos, double zcut){
196
197 //Signal mass dependent SIMP parameters
198 double m_Ap = m_V*(mass_ratio_Ap_to_Vd_);
199 double m_pi = m_Ap/mass_ratio_Ap_to_Pid_;
200 double f_pi = m_pi/ratio_mPi_to_fPi_;
201
202 //Mass in MeV
203 double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_dark_,f_pi,m_l_,rho);
204 double gcTau = ctau * gamma(m_V/1000.0, E_V); //E_V in GeV
205 std::cout << "gcTau: " << gcTau << std::endl;
206 std::cout << "radFrac: " << radFrac << std::endl;
207 std::cout << "radAcc: " << radAcc << std::endl;
208 std::cout << "dNdm:" << dNdm << std::endl;
209
210 //Calculate the Efficiency Vertex (Displaced VD Acceptance)
211 double effVtx = 0.0;
212 for(int zbin = 0; zbin < (effCalc_h->GetTotalHistogram()->GetNbinsX())+1; zbin++){
213 double zz = effCalc_h->GetTotalHistogram()->GetBinLowEdge(zbin);
214 if(zz < zcut) continue;
215 effVtx += (TMath::Exp((target_pos-zz)/gcTau)/gcTau)*
216 (effCalc_h->GetEfficiency(zbin) - effCalc_h->GetEfficiencyErrorLow(zbin))*
217 effCalc_h->GetTotalHistogram()->GetBinWidth(zbin);
218 }
219
220 //Total A' Production Rate
221 double apProduction = (3.*137/2.)*3.14159*(m_Ap*eps*eps*radFrac*dNdm)/radAcc;
222
223 //A' -> V+Pi Branching Ratio
224 double br_VPi = 0.0;
225 if(rho)
226 br_VPi = br_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi);
227 else
228 br_VPi = br_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi);
229
230 std::cout << "Branching ratio is " << br_VPi << std::endl;
231 //Vector to e+e- BR = 1
232 double br_V_ee = 1.0;
233
234 //Expected Signal
235 double expSignal = apProduction * effVtx * br_VPi * br_V_ee;
236
237 return expSignal;
238}
239
240double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho,
241 double E_V, TEfficiency* effCalc_h, double target_pos, double zcut){
242
243 //Signal mass dependent SIMP parameters
244 double m_Ap = m_V*(mass_ratio_Ap_to_Vd_);
245 double m_pi = m_Ap/mass_ratio_Ap_to_Pid_;
246 double f_pi = m_pi/ratio_mPi_to_fPi_;
247
248 //Mass in MeV
249 double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_dark_,f_pi,m_l_,rho);
250 double gcTau = ctau * gamma(m_V/1000.0, E_V); //E_V in GeV
251
252 //Calculate the Efficiency Vertex (Displaced VD Acceptance)
253 double effVtx = 0.0;
254 for(int zbin = 0; zbin < (effCalc_h->GetTotalHistogram()->GetNbinsX())+1; zbin++){
255 double zz = effCalc_h->GetTotalHistogram()->GetBinLowEdge(zbin);
256 if(zz < zcut) continue;
257 effVtx += (TMath::Exp((target_pos-zz)/gcTau)/gcTau)*
258 (effCalc_h->GetEfficiency(zbin) - effCalc_h->GetEfficiencyErrorLow(zbin))*
259 effCalc_h->GetTotalHistogram()->GetBinWidth(zbin);
260 }
261
262
263 //Total A' Production Rate
264 double apProduction = (3.*137/2.)*3.14159*(m_Ap*eps*eps*radiativeFraction(m_Ap)*controlRegionBackgroundRate(m_Ap))
265 /radiativeAcceptance(m_Ap);
266
267 //A' -> V+Pi Branching Ratio
268 double br_VPi;
269 if(rho)
270 br_VPi = br_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi);
271 else
272 br_VPi = br_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi);
273
274 //Vector to e+e- BR = 1
275 double br_V_ee = 1.0;
276
277 //Expected Signal
278 double expSignal = apProduction * effVtx * br_VPi * br_V_ee;
279
280 return expSignal;
281}
282
283double SimpEquations::expectedSignalCalculation(double m_Ap, double m_pi, double m_V, double eps, double alpha_dark,
284 double f_pi, double m_l, bool rho, double E_V, TEfficiency* effCalc_h, double target_pos, double zcut){
285
286 //Mass in MeV
287 double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_dark,f_pi,m_l,rho);
288 double gcTau = ctau * gamma(m_V/1000.0, E_V); //E_V in GeV
289
290 //Calculate the Efficiency Vertex (Displaced VD Acceptance)
291 double effVtx = 0.0;
292 for(int zbin = 0; zbin < (effCalc_h->GetTotalHistogram()->GetNbinsX())+1; zbin++){
293 double zz = effCalc_h->GetTotalHistogram()->GetBinLowEdge(zbin);
294 if(zz < zcut) continue;
295 effVtx += (TMath::Exp((target_pos-zz)/gcTau)/gcTau)*
296 (effCalc_h->GetEfficiency(zbin) - effCalc_h->GetEfficiencyErrorLow(zbin))*
297 effCalc_h->GetTotalHistogram()->GetBinWidth(zbin);
298 }
299
300
301 //Total A' Production Rate
302 double apProduction = (3.*137/2.)*3.14159*(m_Ap*eps*eps*radiativeFraction(m_Ap)*controlRegionBackgroundRate(m_Ap))
303 /radiativeAcceptance(m_Ap);
304
305 //A' -> V+Pi Branching Ratio
306 double br_VPi;
307 if(rho)
308 br_VPi = br_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi);
309 else
310 br_VPi = br_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi);
311
312 //Vector to e+e- BR = 1
313 double br_V_ee = 1.0;
314
315 //Expected Signal
316 double expSignal = apProduction * effVtx * br_VPi * br_V_ee;
317
318 return expSignal;
319}
320
322 double radFrac = 0.0;
323 if(year_ == 2016){
324 radFrac = ( -1.92497e-01 + 1.47144e-02*m_Ap + -2.91966e-04*std::pow(m_Ap,2) + 2.65603e-06*std::pow(m_Ap,3) + -1.12471e-8*std::pow(m_Ap,4) + 1.74765e-11*std::pow(m_Ap,5) + 2.235718e-15*std::pow(m_Ap,6)); //alic 2016 simp reach estimate
325 }
326
327 return radFrac;
328}
329
331 double radAcc = 0.0;
332 if(year_ == 2016){
333 radAcc = ( -7.93151e-01 + 1.04324e-01*m_Ap + -5.55225e-03*std::pow(m_Ap,2) + 1.55480e-04*std::pow(m_Ap,3) + -2.53281e-06*std::pow(m_Ap,4) + 2.54558e-08*std::pow(m_Ap,5) + -1.60877e-10*std::pow(m_Ap,6) + 6.24627e-13*std::pow(m_Ap,7) + -1.36375e-15*std::pow(m_Ap,8) + 1.28312e-18*std::pow(m_Ap,9) ); //alic 2016 simp reach est
334 }
335
336 return radAcc;
337}
338
340 double massRes = 0.0;
341 if(year_ == 2016){
342 massRes = 9.73217e-01 + 3.63659e-02*m_V + -7.32046e-05*m_V*m_V; //2016 simps alic
343 }
344
345 return massRes;
346}
347
349 double dNdm = 0.0;
350 m_Ap = m_Ap/1000.0;
351 /*
352 if(year_ == 2016){
353 std::cout << "[SimpEquations]::WARNING! BACKGROUND NOT DEFINED FOR MASS SPECTRUM!" << std::endl;
354 dNdm = 513800.0;
355 }*/
356
357 if(year_ == 2016){
358 std::cout << "[SimpEquations]::WARNING! USING TEMPORARY BACKGROUND MODEL" << std::endl;
359 //dNdm = 100000 + -1.02011e07*std::pow(m_Ap,1) + 4.02136e08*std::pow(m_Ap,2) + -8.02288e09*std::pow(m_Ap,3)
360 // + 9.16485e10*std::pow(m_Ap,4) + -6.2886e11*std::pow(m_Ap,5) + 2.57095e12*std::pow(m_Ap,6)
361 // + -5.78477e12*std::pow(m_Ap,7) + 5.52322e12*std::pow(m_Ap, 8);
362 dNdm = 204656 + -2.08037e7*std::pow(m_Ap,1) + 8.17706e8*std::pow(m_Ap,2) + -1.62741e10*std::pow(m_Ap,3)
363 + 1.85484e11*std::pow(m_Ap,4) + -1.26987e12*std::pow(m_Ap,5) + 5.17974e12*std::pow(m_Ap,6)
364 + -1.16277e13*std::pow(m_Ap,7) + 1.10758e13*std::pow(m_Ap, 8);
365 }
366
367 std::cout << "LOOK: Background Rate taken from Control Region InvMass is " << dNdm << std::endl;
368 return dNdm;
369}
370
371
double Tv(bool rho, bool phi)
double f(double r)
double rate_Vrho_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi)
double expectedSignalCalculation(double m_V, double eps, bool rho, double E_V, TEfficiency *effCalc_h, double target_pos, double zcut)
double getCtau(double m_Ap, double m_pi, double m_V, double eps, double alpha_dark, double f_pi, double m_l, bool rho)
double rate_Vphi_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi)
double rate_2V(double m_Ap, double m_V, double alpha_dark)
double controlRegionBackgroundRate(double m_Ap)
double rate_Vcharged_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi)
double radiativeFraction(double m_Ap)
void loadParametersConfig(const std::string paramsConfigFile)
double Beta(double x, double y)
double br_Vrho_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi)
int year_
year (used to specify polynomial choices)
double rate_2pi(double m_Ap, double m_pi, double m_V, double alpha_dark)
double radiativeAcceptance(double m_Ap)
SimpEquations(int year)
double br_Vphi_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi)
double m_l_
default lepton mass (ele/pos only)
double gamma(double m_V, double E_V)
double br_Vcharged_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi)
double br_2V(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi, double rho, double phi)
double mass_ratio_Ap_to_Pid_
default Ap to dark pion mass ratio
double alpha_dark_
default A' to dark meson coupling
double rate_2l(double m_Ap, double m_pi, double m_V, double eps, double alpha_dark, double f_pi, double m_l, bool rho)
double br_2pi(double m_Ap, double m_pi, double m_V, double alpha_dark, double f_pi)
double ratio_mPi_to_fPi_
defualt dark pion mass to decay constant ratio
json params_config_
read in simp params
double mass_ratio_Ap_to_Vd_
default Ap to dark vector mass ratio
double massResolution(double m_V)