GENIEGenerator
Loading...
Searching...
No Matches
HNLBRCalculator.cxx
Go to the documentation of this file.
1
2//----------------------------------------------------------------------------
3/*!
4 Copyright (c) 2003-2025, The GENIE Collaboration
5 For the full text of the license visit http://copyright.genie-mc.org
6
7 Author: John Plows <komninos-john.plows \at physics.ox.ac.uk>
8 University of Oxford
9 */
10//----------------------------------------------------------------------------
11
12#include "Physics/BeamHNL/HNLBRCalculator.h"
13
14using namespace genie;
15using namespace genie::hnl;
16
17//----------------------------------------------------------------------------
18BRCalculator::BRCalculator() :
19 ChannelCalculatorI("genie::hnl::BRCalculator")
20{
21
22}
23//----------------------------------------------------------------------------
24BRCalculator::BRCalculator(string name) :
25 ChannelCalculatorI(name)
26{
27
28}
29//----------------------------------------------------------------------------
30BRCalculator::BRCalculator(string name, string config) :
31 ChannelCalculatorI(name, config)
32{
33
34}
35//----------------------------------------------------------------------------
40//----------------------------------------------------------------------------
41void BRCalculator::Configure(const Registry & config)
42{
43 Algorithm::Configure(config);
44 this->LoadConfig();
45}
46//----------------------------------------------------------------------------
47void BRCalculator::Configure(string config)
48{
49 Algorithm::Configure(config);
50 this->LoadConfig();
51}
52//----------------------------------------------------------------------------
53void BRCalculator::LoadConfig(void)
54{
55 if( fIsConfigLoaded ) return;
56
57 this->GetParam( "HNL-Mass", fMass );
58 this->GetParamVect( "HNL-LeptonMixing", fCouplings );
59 fUe42 = fCouplings.at(0);
60 fUm42 = fCouplings.at(1);
61 fUt42 = fCouplings.at(2);
62 this->GetParam( "HNL-Majorana", fMajorana );
63
64 mPi0 = PDGLibrary::Instance()->Find(genie::kPdgPi0)->Mass();
65 mPi = PDGLibrary::Instance()->Find(genie::kPdgPiP)->Mass();
66 mMu = PDGLibrary::Instance()->Find(genie::kPdgMuon)->Mass();
67 mK = PDGLibrary::Instance()->Find(genie::kPdgKP)->Mass();
68 mK0 = PDGLibrary::Instance()->Find(genie::kPdgK0)->Mass();
69 mE = PDGLibrary::Instance()->Find(genie::kPdgElectron)->Mass();
70
71 this->GetParam( "WeinbergAngle", wAng );
72 s2w = std::pow( std::sin( wAng ), 2.0 );
73
74 this->GetParam( "CKM-Vud", Vud );
75 Vud2 = Vud * Vud;
76
77 this->GetParam( "Pion-FFactor", fpi );
78 fpi2 = fpi * fpi;
79
80 BR_C1 = 1./4. * ( 1. - 4. * s2w + 8. * s2w * s2w );
81 BR_C2 = 1./2. * ( -s2w + 2. * s2w * s2w );
82
83 this->GetParam( "PMNS-Ue1", Ue1 );
84 this->GetParam( "PMNS-Ue2", Ue2 );
85 this->GetParam( "PMNS-Ue3", Ue3 );
86 this->GetParam( "PMNS-Um1", Um1 );
87 this->GetParam( "PMNS-Um2", Um2 );
88 this->GetParam( "PMNS-Um3", Um3 );
89 this->GetParam( "PMNS-Ut1", Ut1 );
90 this->GetParam( "PMNS-Ut2", Ut2 );
91 this->GetParam( "PMNS-Ut3", Ut3 );
92
93 kscale_K3e = {
94 { 0.0, 1.0 }, { 0.01, (2.0 + 0.968309)/3.0 },
95 { 0.019970, 0.968309 }, { 0.029963, 0.952842 }, { 0.040037, 0.922646 }, { 0.049839, 0.907908 },
96 { 0.060075, 0.879136 }, { 0.070117, 0.824297 }, { 0.080272, 0.772880 }, { 0.090139, 0.736432 },
97 { 0.099924, 0.679466 }, { 0.110059, 0.607039 }, { 0.120445, 0.560082 }, { 0.130123, 0.508503 },
98 { 0.140579, 0.461674 }, { 0.150900, 0.405874 }, { 0.159906, 0.356819 }, { 0.170544, 0.303751 },
99 { 0.180722, 0.267038 }, { 0.190278, 0.227323 }, { 0.200340, 0.193514 }, { 0.209579, 0.159513 },
100 { 0.219244, 0.129386 }, { 0.230837, 0.101623 }, { 0.239932, 0.081113 }, { 0.249386, 0.060704 },
101 { 0.260887, 0.043990 }, { 0.269425, 0.031878 }, { 0.280041, 0.021660 }, { 0.289206, 0.014251 },
102 { 0.300601, 0.007729 }, { 0.310440, 0.004059 }, { 0.320600, 0.001935 }, { 0.330028, 0.000713 },
103 { 0.339733, 0.000184 }, { 0.350852, 0.00000953 }
104 };
105
106 kscale_K3mu = {
107 { 0.0, 1.0 }, { 0.01, (2.0 + 0.968309)/3.0 },
108 { 0.019970, 0.968309 }, { 0.029963, 0.937622 }, { 0.040037, 0.907908 }, { 0.049839, 0.879136 },
109 { 0.060075, 0.824297 }, { 0.070117, 0.772880 }, { 0.079757, 0.713094 }, { 0.090139, 0.647424 },
110 { 0.100570, 0.578413 }, { 0.110059, 0.525146 }, { 0.120045, 0.454300 }, { 0.130964, 0.393012 },
111 { 0.140127, 0.334561 }, { 0.149932, 0.275778 }, { 0.159906, 0.227323 }, { 0.170544, 0.181443 },
112 { 0.180722, 0.137994 }, { 0.190278, 0.101623 }, { 0.200340, 0.071309 }, { 0.210933, 0.046917 },
113 { 0.220661, 0.025445 }, { 0.229355, 0.013362 }, { 0.239932, 0.003930 }, { 0.250997, 0.000037 }
114 };
115
116 kscale_mu3e = {
117 { 0.0, 1.0 }, { 0.01, (2.0 + 0.772880)/3.0 },
118 { 0.020099, 0.772880 }, { 0.029963, 0.560082 }, { 0.040037, 0.356819 }, { 0.050161, 0.193514 },
119 { 0.060075, 0.089341 }, { 0.069667, 0.032922 }, { 0.080272, 0.007247 }, { 0.090139, 0.000713 },
120 { 0.100570, 0.00000363 }
121 };
122
123 fIsConfigLoaded = true;
124}
125//----------------------------------------------------------------------------
126double BRCalculator::KinematicScaling( HNLProd_t hnlprod ) const
127{
128 return this->KScale_Global( hnlprod, fMass );
129}
130//----------------------------------------------------------------------------
131double BRCalculator::DecayWidth( HNLDecayMode_t hnldm ) const
132{
133 return this->DWidth_Global( hnldm, fMass );
134}
135//----------------------------------------------------------------------------
136double BRCalculator::GetFormfactorF1( double x ) const {
137 if( x < 0. || x > 0.5 ) { LOG( "HNL", pERROR ) << "BRCalculator::GetFormfactorF1:: Illegal x = " << x; exit( 3 ); }
138 if( x == 0.5 ) return 0.;
139 int i = x/selector::PARTWIDTH;
140 if( x - i*selector::PARTWIDTH ==0 ) return selector::FormfactorF1[i];
141 return 1./2. * ( selector::FormfactorF1[i] + selector::FormfactorF1[i+1] );
142}
143//----------------------------------------------------------------------------
144double BRCalculator::GetFormfactorF2( double x ) const {
145 if( x < 0. || x > 0.5 ) { LOG( "HNL", pERROR ) << "BRCalculator::GetFormfactorF2:: Illegal x = " << x; exit( 3 ); }
146 if( x == 0.5 ) return 0.;
147 int i = x/selector::PARTWIDTH;
148 if( x - i*selector::PARTWIDTH==0 ) return selector::FormfactorF2[i];
149 return 1./2. * ( selector::FormfactorF2[i] + selector::FormfactorF2[i+1] );
150}
151//----------------------------------------------------------------------------
152// interface to scale factors
153double BRCalculator::KScale_Global( HNLProd_t hnldm, const double M ) const {
154 if( !utils::hnl::IsProdKinematicallyAllowed( hnldm ) ){
155 return 0.0;
156 }
157
158 switch( hnldm ){
159 case kHNLProdPion2Muon: return KScale_PseudoscalarToLepton( mPi, M, mMu );
160 case kHNLProdPion2Electron: return KScale_PseudoscalarToLepton( mPi, M, mE );
161 case kHNLProdKaon2Muon: return KScale_PseudoscalarToLepton( mK, M, mMu );
162 case kHNLProdKaon2Electron: return KScale_PseudoscalarToLepton( mK, M, mE );
163 case kHNLProdKaon3Muon: return KScale_PseudoscalarToPiLepton( mK, M, mMu );
164 case kHNLProdKaon3Electron: return KScale_PseudoscalarToPiLepton( mK, M, mE );
165 case kHNLProdNeuk3Muon: return KScale_PseudoscalarToPiLepton( mK0, M, mMu );
166 case kHNLProdNeuk3Electron: return KScale_PseudoscalarToPiLepton( mK0, M, mE );
167 case kHNLProdMuon3Numu:
168 case kHNLProdMuon3Nue:
169 case kHNLProdMuon3Nutau:
170 return KScale_MuonToNuAndElectron( M );
171 default: return 0.0;
172 }
173
174 return 0.0;
175}
176//----------------------------------------------------------------------------
177// HNL production widths
178double BRCalculator::KScale_PseudoscalarToLepton( const double mP, const double M, const double ma ) const {
179 double da = std::pow( utils::hnl::MassX( ma, mP ) , 2.0 );
180 double di = std::pow( utils::hnl::MassX( M, mP ) , 2.0 );
181 double num = utils::hnl::RhoFunc( da, di );
182 double den = da * std::pow( (1.0 - da), 2.0 );
183 return num/den;
184}
185//----------------------------------------------------------------------------
186double BRCalculator::DWidth_PseudoscalarToLepton( const double mP, const double M, const double Ua42, const double ma ) const {
187 assert( M + ma <= mP );
188
189 double KScale = KScale_PseudoscalarToLepton( mP, M, ma );
190 return Ua42 * KScale;
191}
192//----------------------------------------------------------------------------
193double BRCalculator::KScale_PseudoscalarToPiLepton( const double mP, const double M, const double ma ) const {
194 assert( mP == mK || mP == mK0 ); // RETHERE remove this when/if heavier pseudoscalars are considered
195 assert( ma == mE || ma == mMu );
196
197 std::map< double, double > scaleMap = ( ma == mE ) ? kscale_K3e : kscale_K3mu;
198
199 std::map< double, double >::iterator scmit = scaleMap.begin();
200 // iterate until we know between which two map points M is
201 // if we're very lucky, M will coincide with a map point
202 while( (*scmit).first <= M && scmit != scaleMap.end() ){ ++scmit; }
203 std::map< double, double >::iterator scpit = std::prev( scmit, 1 );
204 //LOG( "HNL", pDEBUG )
205 // << "Requested map for M = " << M << ": iter at ( " << (*scpit).first << ", " << (*scmit).first << " ]";
206 assert( scmit != scaleMap.end() );
207 // if coincide then return scale there
208 if( scaleMap.find( M ) != scaleMap.end() ) return (*scmit).second;
209 // otherwise transform scmit-1 and scmit second to log, do a linear extrapolation and return
210 double l1 = TMath::Log( (*scpit).second );
211 double l2 = TMath::Log( (*scmit).second );
212 double t = ( M - (*scpit).first ) / ( (*scmit).first - (*scpit).first );
213 return TMath::Exp( l1 + ( l2 - l1 ) * t );
214}
215//----------------------------------------------------------------------------
216double BRCalculator::DWidth_PseudoscalarToPiLepton( const double mP, const double M, const double Ua42, const double ma ) const {
217 assert( M + ma + mPi0 <= mP );
218
219 double KScale = KScale_PseudoscalarToPiLepton( mP, M, ma );
220 return Ua42 * KScale;
221}
222//----------------------------------------------------------------------------
223double BRCalculator::KScale_MuonToNuAndElectron( const double M ) const {
224 std::map< double, double > scaleMap = kscale_mu3e;
225 std::map< double, double >::iterator scmit = scaleMap.begin();
226 while( (*scmit).first <= M && scmit != scaleMap.end() ){ ++scmit; }
227 std::map< double, double >::iterator scpit = std::prev( scmit, 1 );
228 //LOG( "HNL", pDEBUG )
229 // << "Requested map for M = " << M << ": iter at ( " << (*scpit).first << ", " << (*scmit).first << " ]";
230 assert( scmit != scaleMap.end() );
231
232 if( scaleMap.find( M ) != scaleMap.end() ) return (*scmit).second;
233
234 double l1 = TMath::Log( (*scpit).second );
235 double l2 = TMath::Log( (*scmit).second );
236 double t = ( M - (*scpit).first ) / ( (*scmit).first - (*scpit).first );
237 return TMath::Exp( l1 + ( l2 - l1 ) * t );
238}
239//----------------------------------------------------------------------------
240double BRCalculator::DWidth_MuonToNuAndElectron( const double M, const double Ue42, const double Umu42, const double Ut42 ) const {
241 assert( M + mE <= mMu );
242
243 double KScale = KScale_MuonToNuAndElectron( M );
244 return ( Ue42 + Umu42 + Ut42 ) * KScale;
245}
246//----------------------------------------------------------------------------
247// interface to decay widths
248double BRCalculator::DWidth_Global( HNLDecayMode_t hnldm, const double M ) const {
249 if( !utils::hnl::IsKinematicallyAllowed( hnldm, M ) ){
250 return 0.0;
251 }
252
253 switch( hnldm ){
254 case kHNLDcyNuNuNu : return this->DWidth_Invisible( M, fUe42, fUm42, fUt42 );
255 case kHNLDcyNuEE : return this->DWidth_SameLepton( M, fUe42, fUm42, fUt42, mE, false );
256 case kHNLDcyNuMuE : return (this->DWidth_DiffLepton( M, fUe42, fUm42, fMajorana ) + this->DWidth_DiffLepton( M, fUm42, fUe42, fMajorana ));
257 case kHNLDcyPi0Nu : return this->DWidth_PiZeroAndNu( M, fUe42, fUm42, fUt42 );
258 case kHNLDcyPiE : return this->DWidth_PiAndLepton( M, fUe42, mE );
259 case kHNLDcyNuMuMu : return this->DWidth_SameLepton( M, fUe42, fUm42, fUt42, mMu, true );
260 case kHNLDcyPiMu : return this->DWidth_PiAndLepton( M, fUm42, mMu );
261 case kHNLDcyPi0Pi0Nu : return this->DWidth_Pi0Pi0Nu( M, fUe42, fUm42, fUt42 );
262 case kHNLDcyPiPi0E : return this->DWidth_PiPi0Ell( M, mE, fUe42, fUm42, fUt42, true );
263 case kHNLDcyPiPi0Mu : return this->DWidth_PiPi0Ell( M, mMu, fUe42, fUm42, fUt42, false );
264 default: return 0.0;
265 }
266
267 return 0.0;
268}
269//----------------------------------------------------------------------------
270// total decay widths, various channels
271double BRCalculator::DWidth_PiZeroAndNu( const double M, const double Ue42, const double Umu42, const double Ut42 ) const {
272 const double x = genie::utils::hnl::MassX( mPi0, M );
273 const double preFac = GF2 * M*M*M / ( 32. * pi );
274 const double kinPart = ( 1. - x*x ) * ( 1. - x*x );
275 return preFac * ( Ue42 + Umu42 + Ut42 ) * fpi2 * kinPart;
276}
277//----------------------------------------------------------------------------
278double BRCalculator::DWidth_PiAndLepton( const double M, const double Ua42, const double ma ) const {
279 const double xPi = genie::utils::hnl::MassX( mPi, M );
280 const double xLep = genie::utils::hnl::MassX( ma, M );
281 const double preFac = GF2 * M*M*M / ( 16. * pi );
282 const double kalPart = TMath::Sqrt( genie::utils::hnl::Kallen( 1, xPi*xPi, xLep*xLep ) );
283 const double othPart = 1. - xPi*xPi - xLep*xLep * ( 2. + xPi*xPi - xLep*xLep );
284 return preFac * fpi2 * Ua42 * Vud2 * kalPart * othPart;
285}
286//----------------------------------------------------------------------------
287double BRCalculator::DWidth_Invisible( const double M, const double Ue42, const double Umu42, const double Ut42 ) const {
288 const double preFac = GF2 * TMath::Power( M, 5. ) / ( 192. * pi*pi*pi );
289 return preFac * ( Ue42 + Umu42 + Ut42 );
290}
291//----------------------------------------------------------------------------
292double BRCalculator::DWidth_SameLepton( const double M, const double Ue42, const double Umu42, const double Ut42, const double mb, bool bIsMu ) const {
293 const double preFac = GF2 * TMath::Power( M, 5. ) / ( 192. * pi*pi*pi );
294 const double x = genie::utils::hnl::MassX( mb, M );
295 const double f1 = GetFormfactorF1( x );
296 const double f2 = GetFormfactorF2( x );
297 const double C1Part = ( Ue42 + Umu42 + Ut42 ) * f1 * BR_C1;
298 const double C2Part = ( Ue42 + Umu42 + Ut42 ) * f2 * BR_C2;
299 const double D1Part = bIsMu ? 2. * s2w * Umu42 * f1 : 2. * s2w * Ue42 * f1;
300 const double D2Part = bIsMu ? s2w * Umu42 * f2 : s2w * Ue42 * f2;
301 return preFac * ( C1Part + C2Part + D1Part + D2Part );
302}
303//----------------------------------------------------------------------------
304double BRCalculator::DWidth_DiffLepton( const double M, const double Ua42, const double Ub42, const int IsMajorana ) const {
305 const double preFac = GF2 * TMath::Power( M, 5. ) / ( 192. * pi*pi*pi );
306 const double x = genie::utils::hnl::MassX( mMu, M );
307 const double kinPol = 1. - 8. * x*x + 8. * TMath::Power( x, 6. ) - TMath::Power( x, 8. );
308 const double kinLn = -12. * TMath::Power( x, 4. ) * TMath::Log( x*x );
309 const double kinPart = kinPol + kinLn;
310 const double coupPart = (!IsMajorana) ? Ua42 : Ua42 + Ub42; // 2nd diagram in Majorana case!
311 return preFac * kinPart * coupPart;
312}
313//----------------------------------------------------------------------------
314// note that these BR are very very tiny.
315double BRCalculator::DWidth_PiPi0Ell( const double M, const double ml,
316 const double Ue42, const double Umu42, const double Ut42,
317 const bool isElectron) const
318{
319 // because the actual decay width is very hard to integrate onto a full DWidth,
320 // build 2Differential and then integrate numerically
321 // using Simpson's method for 2D.
322
323 const double preFac = fpi2 * fpi2 * GF2 * GF2 * Vud2 * M / ( 32.0 * pi*pi*pi );
324 const double Ua1 = isElectron ? Ue1 : Um1;
325 const double Ua2 = isElectron ? Ue2 : Um2;
326 const double Ua3 = isElectron ? Ue3 : Um3;
327 __attribute__((unused)) const double Ua4 = isElectron ? std::sqrt( Ue42 ) : std::sqrt( Umu42 );
328
329 const double Ue4 = std::sqrt( Ue42 );
330 const double Um4 = std::sqrt( Umu42 );
331 const double Ut4 = std::sqrt( Ut42 );
332 // assume all these to be real
333 const double bigMats =
334 Ua1 * ( Ue4 * Ue1 + Um4 * Um1 + Ut4 * Ut1 ) +
335 Ua2 * ( Ue4 * Ue2 + Um4 * Um2 + Ut4 * Ut2 ) +
336 Ua3 * ( Ue4 * Ue3 + Um4 * Um3 + Ut4 * Ut3 );
337
338 // now limits
339 const double maxMu =
340 ( ( M - mPi0 ) * ( M - mPi0 ) - mPi*mPi + ml*ml ) / ( 2.0 * ( M - mPi0 ) );
341 const double maxPi =
342 ( ( M - mPi0 ) * ( M - mPi0 ) + mPi*mPi - ml*ml ) / ( 2.0 * ( M - mPi0 ) );
343
344 // gotta put in the formula
345 TF2 * f = new TF2( "fPiPi0Ell", PiPi0EllForm, mPi, maxPi, ml, maxMu, 4 );
346 f->SetParameter( 0, M );
347 f->SetParameter( 1, ml );
348 f->SetParameter( 2, mPi );
349 f->SetParameter( 3, mPi0 );
350
351 // now we can use composite Simpson, iterating on both axes simultaneously
352 // This is like using Fubini over and over again for sampled Emu ==> integrate
353 // out Epi ==> Simpson again for Emu. Can see more at
354 // https://math.stackexchange.com/questions/1319892/simpsons-rule-for-double-integrals.
355
356 const int nSteps = 10000 + 1;
357 const double hEMu = ( maxMu - ml ) / ( nSteps - 1 );
358 const double hEPi = ( maxPi - mPi ) / ( nSteps - 1 );
359 const double preSimp = hEMu * hEPi / ( 9.0 * ( nSteps - 1 ) * ( nSteps - 1 ) );
360
361 double intNow = 0.0;
362 for( int i = 0; i < nSteps; i++ ){
363 for( int j = 0; j < nSteps; j++ ){
364 double midW = 0.0;
365 //determine midpoint coefficient for this step
366 if( i % (nSteps - 1) == 0 ){ // edge case i
367 if( j % (nSteps - 1) == 0 ){ midW = 1.0; } // edge case j
368 else if( j % 2 == 0 ){ midW = 2.0; } // even j
369 else{ midW = 4.0; } // odd j
370 }
371 else if( i % 2 == 0 ){ // even i
372 if( j % (nSteps - 1) == 0 ){ midW = 2.0; } // edge case j
373 else if( j % 2 == 0 ){ midW = 4.0; } // even j
374 else{ midW = 8.0; } // odd j
375 }
376 else{ // odd i
377 if( j % (nSteps - 1) == 0 ){ midW = 4.0; } // edge case j
378 else if( j % 2 == 0 ){ midW = 8.0; } // even j
379 else{ midW = 16.0; } // odd j
380 }
381 // finally, evaluate f at this point
382 const double xev = mPi + i * hEPi;
383 const double yev = ml + j * hEMu;
384 const double fev = f->Eval( xev, yev );
385
386 // and add to integral
387 intNow += std::abs( preSimp * midW * fev );
388 }
389 }
390
391 delete f;
392
393 intNow *= preFac * bigMats;
394
395 return intNow;
396
397}
398//----------------------------------------------------------------------------
399// *especially* this channel, there's N4 in the propagator so it emits *both* the pi-zeros!!!
400// It is subleading in |U_\ell 4|^2, therefore not important to get this exactly right
401double BRCalculator::DWidth_Pi0Pi0Nu( const double M,
402 const double Ue42, const double Umu42, const double Ut42 ) const
403{
404 const double preFac = fpi2 * fpi2 * GF2 * GF2 * std::pow( M, 5.0 ) / ( 64.0 * pi*pi*pi );
405
406 const double Ue4 = std::sqrt( Ue42 );
407 const double Um4 = std::sqrt( Umu42 );
408 const double Ut4 = std::sqrt( Ut42 );
409
410 // once again, assume all PMNS matrix elements real
411 const double bigMats = std::pow( Ue4 * ( Ue1 + Ue2 + Ue3 ) +
412 Um4 * ( Um1 + Um2 + Um3 ) +
413 Ut4 * ( Ut1 + Ut2 + Ut3 ), 2.0 );
414 const double smallMats = std::pow( Ue42 + Umu42 + Ut42 , 2.0 );
415
416 // let's make the limits
417 const double maxNu =
418 ( ( M - mPi0 ) * ( M - mPi0 ) - mPi0*mPi0 ) / ( 2.0 * ( M - mPi0 ) );
419 const double maxPi =
420 ( ( M - mPi0 ) * ( M - mPi0 ) + mPi0*mPi0 ) / ( 2.0 * ( M - mPi0 ) );
421
422 // gotta put in the formula
423 TF2 * f = new TF2( "fPi0Pi0Nu", Pi0Pi0NuForm, mPi0, maxPi, 0.0, maxNu, 2 );
424 f->SetParameter( 0, M );
425 f->SetParameter( 1, mPi0 );
426
427 // using composite Simpson to evaluate
428
429 const int nSteps = 10000 + 1;
430 const double hENu = ( maxNu - 0.0 ) / ( nSteps - 1 );
431 const double hEPi = ( maxPi - mPi0 ) / ( nSteps - 1 );
432 const double preSimp = hENu * hEPi / ( 9.0 * ( nSteps - 1 ) * ( nSteps - 1 ) );
433
434 double intNow = 0.0;
435 for( int i = 0; i < nSteps; i++ ){
436 for( int j = 0; j < nSteps; j++ ){
437 double midW = 0.0;
438 //determine midpoint coefficient for this step
439 if( i % (nSteps - 1) == 0 ){ // edge case i
440 if( j % (nSteps - 1) == 0 ){ midW = 1.0; } // edge case j
441 else if( j % 2 == 0 ){ midW = 2.0; } // even j
442 else{ midW = 4.0; } // odd j
443 }
444 else if( i % 2 == 0 ){ // even i
445 if( j % (nSteps - 1) == 0 ){ midW = 2.0; } // edge case j
446 else if( j % 2 == 0 ){ midW = 4.0; } // even j
447 else{ midW = 8.0; } // odd j
448 }
449 else{ // odd i
450 if( j % (nSteps - 1) == 0 ){ midW = 4.0; } // edge case j
451 else if( j % 2 == 0 ){ midW = 8.0; } // even j
452 else{ midW = 16.0; } // odd j
453 }
454 // finally, evaluate f at this point
455 const double xev = mPi0 + i * hEPi;
456 const double yev = 0.0 + j * hENu;
457 const double fev = f->Eval( xev, yev );
458
459 // and add to integral
460 intNow += std::abs( preSimp * midW * fev );
461 }
462 }
463
464 delete f;
465
466 intNow *= preFac * bigMats * smallMats;
467
468 return intNow;
469}
470//----------------------------------------------------------------------------
471// formula for N --> pi pi0 ell decay rate
472double BRCalculator::PiPi0EllForm( double *x, double *par ){
473 double MN = par[0];
474 double MMu = par[1];
475 double MPi = par[2];
476 double MPi0 = par[3];
477
478 double Epi = x[0];
479 double Emu = x[1];
480
481 double pi0Term = ( MN - Emu - Epi > MPi0 ) ?
482 std::sqrt( std::pow( ( MN - Emu - Epi ), 2.0 ) - MPi0 * MPi0 ) : 0.0;
483
484 double ETerm =
485 std::sqrt( Emu*Emu - MMu*MMu ) *
486 std::sqrt( Epi*Epi - MPi*MPi ) *
487 pi0Term / ( MN - Emu - Epi );
488
489 double FracNum1 = MN*MN - 2.0*( MN-Emu-Epi )*MN + MPi0*MPi0;
490 double FracNum2 = MN*MN - 2.0*Emu*MN + 2.0*MMu*MMu;
491 double FracNum3 = MN*MN - MPi0*MPi0;
492 double FracNum4 = MN*MN - 2.0*( MN-Emu-Epi )*MN + MPi0*MPi0 + MMu*MMu - MPi*MPi;
493 double FracNum = FracNum1*FracNum2 - FracNum3*FracNum4;
494 double FracDen = std::pow( MN*MN - 2.0*( MN - Emu - Epi ) * MN + MPi0*MPi0 , 2.0 );
495
496 return ETerm * FracNum / FracDen;
497}
498//----------------------------------------------------------------------------
499// formula for N --> pi0 pi0 nu decay rate
500double BRCalculator::Pi0Pi0NuForm( double *x, double *par ){
501 double MN = par[0];
502 double MPi0 = par[1];
503
504 double Epi = x[0]; // leading pi-zero energy
505 double Enu = x[1];
506
507 double ETerm =
508 std::sqrt( Epi*Epi - MPi0*MPi0 ) *
509 (Enu + MN) * Enu * Enu *
510 (MN - Enu - Epi);
511
512 double Frac1 = 1.0 / ( Enu * ( MN - Enu - Epi ) + MPi0 * MPi0 - MN * MN );
513 double Frac2 = 1.0 / ( Enu * Epi + MPi0 * MPi0 - MN * MN );
514
515 return ETerm * std::pow( ( Frac1 + Frac2 ), 2.0 );
516}
pERROR
#define pERROR
Definition Messenger.h:59
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::Algorithm::GetParamVect
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition PDGLibrary.cxx:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition PDGLibrary.cxx:86
genie::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition Registry.h:65
genie::hnl::BRCalculator::DWidth_SameLepton
double DWidth_SameLepton(const double M, const double Ue42, const double Umu42, const double Ut42, const double mb, bool bIsMu) const
Definition HNLBRCalculator.cxx:292
genie::hnl::BRCalculator::DWidth_PseudoscalarToPiLepton
double DWidth_PseudoscalarToPiLepton(const double mP, const double M, const double Ua42, const double ma) const
Definition HNLBRCalculator.cxx:216
genie::hnl::BRCalculator::DWidth_PiAndLepton
double DWidth_PiAndLepton(const double M, const double Ua42, const double ma) const
Definition HNLBRCalculator.cxx:278
genie::hnl::BRCalculator::GetFormfactorF1
double GetFormfactorF1(double x) const
Definition HNLBRCalculator.cxx:136
genie::hnl::BRCalculator::fCouplings
std::vector< double > fCouplings
Definition HNLBRCalculator.h:147
genie::hnl::BRCalculator::DWidth_Global
double DWidth_Global(genie::hnl::HNLDecayMode_t hnldm, const double M) const
Definition HNLBRCalculator.cxx:248
genie::hnl::BRCalculator::DecayWidth
double DecayWidth(genie::hnl::HNLDecayMode_t hnldm) const
Definition HNLBRCalculator.cxx:131
genie::hnl::BRCalculator::DWidth_PseudoscalarToLepton
double DWidth_PseudoscalarToLepton(const double mP, const double M, const double Ua42, const double ma) const
Definition HNLBRCalculator.cxx:186
genie::hnl::BRCalculator::KScale_PseudoscalarToLepton
double KScale_PseudoscalarToLepton(const double mP, const double M, const double ma) const
Definition HNLBRCalculator.cxx:178
genie::hnl::BRCalculator::kscale_K3mu
std::map< double, double > kscale_K3mu
Definition HNLBRCalculator.h:158
genie::hnl::BRCalculator::Configure
void Configure(const Registry &config)
Definition HNLBRCalculator.cxx:41
genie::hnl::BRCalculator::Pi0Pi0NuForm
static double Pi0Pi0NuForm(double *x, double *par)
Definition HNLBRCalculator.cxx:500
genie::hnl::BRCalculator::KScale_Global
double KScale_Global(genie::hnl::HNLProd_t hnldm, const double M) const
Definition HNLBRCalculator.cxx:153
genie::hnl::BRCalculator::DWidth_MuonToNuAndElectron
double DWidth_MuonToNuAndElectron(const double M, const double Ue42, const double Umu42, const double Ut42) const
Definition HNLBRCalculator.cxx:240
genie::hnl::BRCalculator::DWidth_PiZeroAndNu
double DWidth_PiZeroAndNu(const double M, const double Ue42, const double Umu42, const double Ut42) const
Definition HNLBRCalculator.cxx:271
genie::hnl::BRCalculator::GetFormfactorF2
double GetFormfactorF2(double x) const
Definition HNLBRCalculator.cxx:144
genie::hnl::BRCalculator::DWidth_DiffLepton
double DWidth_DiffLepton(const double M, const double Ua42, const double Ub42, const int IsMajorana) const
Definition HNLBRCalculator.cxx:304
genie::hnl::BRCalculator::DWidth_Invisible
double DWidth_Invisible(const double M, const double Ue42, const double Umu42, const double Ut42) const
Definition HNLBRCalculator.cxx:287
genie::hnl::BRCalculator::kscale_K3e
std::map< double, double > kscale_K3e
Definition HNLBRCalculator.h:158
genie::hnl::BRCalculator::PiPi0EllForm
static double PiPi0EllForm(double *x, double *par)
Definition HNLBRCalculator.cxx:472
genie::hnl::BRCalculator::KScale_PseudoscalarToPiLepton
double KScale_PseudoscalarToPiLepton(const double mP, const double M, const double ma) const
Definition HNLBRCalculator.cxx:193
genie::hnl::BRCalculator::kscale_mu3e
std::map< double, double > kscale_mu3e
Definition HNLBRCalculator.h:158
genie::hnl::BRCalculator::KinematicScaling
double KinematicScaling(genie::hnl::HNLProd_t hnlprod) const
Definition HNLBRCalculator.cxx:126
genie::hnl::BRCalculator::KScale_MuonToNuAndElectron
double KScale_MuonToNuAndElectron(const double M) const
Definition HNLBRCalculator.cxx:223
genie::hnl::BRCalculator::DWidth_PiPi0Ell
double DWidth_PiPi0Ell(const double M, const double ml, const double Ue42, const double Umu42, const double Ut42, const bool isElectron=false) const
Definition HNLBRCalculator.cxx:315
genie::hnl::BRCalculator::DWidth_Pi0Pi0Nu
double DWidth_Pi0Pi0Nu(const double M, const double Ue42, const double Umu42, const double Ut42) const
Definition HNLBRCalculator.cxx:401
genie::hnl::selector::PARTWIDTH
static const double PARTWIDTH
Definition FormFactorTables.h:30
genie::hnl::selector::FormfactorF1
static const double FormfactorF1[]
Definition FormFactorTables.h:58
genie::hnl::selector::FormfactorF2
static const double FormfactorF2[]
Definition FormFactorTables.h:109
genie::hnl::HNLProd_t
enum genie::hnl::t_HNLProd HNLProd_t
genie::hnl::HNLDecayMode_t
enum genie::hnl::EHNLDecayMode HNLDecayMode_t
genie::utils::hnl::MassX
double MassX(double m1, double m2)
Definition HNLKinUtils.h:32
genie::utils::hnl::IsKinematicallyAllowed
bool IsKinematicallyAllowed(genie::hnl::HNLDecayMode_t hnldm, double Mhnl)
genie::utils::hnl::IsProdKinematicallyAllowed
bool IsProdKinematicallyAllowed(genie::hnl::HNLProd_t hnlprod)
genie::utils::hnl::Kallen
double Kallen(double x, double y, double z)
Definition HNLKinUtils.h:37
genie::utils::hnl::RhoFunc
double RhoFunc(double x, double y)
Definition HNLKinUtils.h:45
genie
THE MAIN GENIE PROJECT NAMESPACE
Definition AlgCmp.h:25
genie::kPdgPi0
const int kPdgPi0
Definition PDGCodes.h:160
genie::kPdgKP
const int kPdgKP
Definition PDGCodes.h:172
genie::kPdgPiP
const int kPdgPiP
Definition PDGCodes.h:158
genie::kPdgMuon
const int kPdgMuon
Definition PDGCodes.h:37
genie::kPdgElectron
const int kPdgElectron
Definition PDGCodes.h:35
genie::kPdgK0
const int kPdgK0
Definition PDGCodes.h:174
Generated on for GENIEGenerator by doxygen 1.14.0