GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HNLBRCalculator.cxx
Go to the documentation of this file.
1 
2 //----------------------------------------------------------------------------
3 /*!
4  Copyright (c) 2003-2024, 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 
14 using namespace genie;
15 using namespace genie::hnl;
16 
17 //----------------------------------------------------------------------------
18 BRCalculator::BRCalculator() :
19  ChannelCalculatorI("genie::hnl::BRCalculator")
20 {
21 
22 }
23 //----------------------------------------------------------------------------
24 BRCalculator::BRCalculator(string name) :
25  ChannelCalculatorI(name)
26 {
27 
28 }
29 //----------------------------------------------------------------------------
30 BRCalculator::BRCalculator(string name, string config) :
31  ChannelCalculatorI(name, config)
32 {
33 
34 }
35 //----------------------------------------------------------------------------
36 BRCalculator::~BRCalculator()
37 {
38 
39 }
40 //----------------------------------------------------------------------------
41 void BRCalculator::Configure(const Registry & config)
42 {
43  Algorithm::Configure(config);
44  this->LoadConfig();
45 }
46 //----------------------------------------------------------------------------
47 void BRCalculator::Configure(string config)
48 {
49  Algorithm::Configure(config);
50  this->LoadConfig();
51 }
52 //----------------------------------------------------------------------------
53 void 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 //----------------------------------------------------------------------------
126 double BRCalculator::KinematicScaling( HNLProd_t hnlprod ) const
127 {
128  return this->KScale_Global( hnlprod, fMass );
129 }
130 //----------------------------------------------------------------------------
131 double BRCalculator::DecayWidth( HNLDecayMode_t hnldm ) const
132 {
133  return this->DWidth_Global( hnldm, fMass );
134 }
135 //----------------------------------------------------------------------------
136 double 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 //----------------------------------------------------------------------------
144 double 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
153 double 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
178 double 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 //----------------------------------------------------------------------------
186 double 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 //----------------------------------------------------------------------------
193 double 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 //----------------------------------------------------------------------------
216 double 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 //----------------------------------------------------------------------------
223 double 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 //----------------------------------------------------------------------------
240 double 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
248 double 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
271 double 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 //----------------------------------------------------------------------------
278 double 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 //----------------------------------------------------------------------------
287 double 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 //----------------------------------------------------------------------------
292 double 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 //----------------------------------------------------------------------------
304 double 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.
315 double 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
401 double 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
472 double 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
500 double 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 }
genie::hnl::HNLProd_t
enum genie::hnl::t_HNLProd HNLProd_t
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_PseudoscalarToLepton
double DWidth_PseudoscalarToLepton(const double mP, const double M, const double Ua42, const double ma) const
Definition: HNLBRCalculator.cxx:186
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
pERROR
#define pERROR
Definition: Messenger.h:59
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::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::kscale_K3e
std::map< double, double > kscale_K3e
Definition: HNLBRCalculator.h:158
genie::hnl::selector::PARTWIDTH
static const double PARTWIDTH
Definition: FormFactorTables.h:30
genie::kPdgElectron
const int kPdgElectron
Definition: PDGCodes.h:35
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::utils::hnl::Kallen
double Kallen(double x, double y, double z)
Definition: HNLKinUtils.h:37
genie::kPdgK0
const int kPdgK0
Definition: PDGCodes.h:174
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_PiAndLepton
double DWidth_PiAndLepton(const double M, const double Ua42, const double ma) const
Definition: HNLBRCalculator.cxx:278
genie::hnl::selector::FormfactorF2
static const double FormfactorF2[]
Definition: FormFactorTables.h:109
genie::hnl::BRCalculator::KinematicScaling
double KinematicScaling(genie::hnl::HNLProd_t hnlprod) const
Definition: HNLBRCalculator.cxx:126
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::hnl::BRCalculator::KScale_PseudoscalarToPiLepton
double KScale_PseudoscalarToPiLepton(const double mP, const double M, const double ma) const
Definition: HNLBRCalculator.cxx:193
genie::units::mb
static constexpr double mb
Definition: Units.h:79
genie::utils::hnl::RhoFunc
double RhoFunc(double x, double y)
Definition: HNLKinUtils.h:45
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::hnl::HNLDecayMode_t
enum genie::hnl::EHNLDecayMode HNLDecayMode_t
genie::hnl::selector::FormfactorF1
static const double FormfactorF1[]
Definition: FormFactorTables.h:58
genie::hnl::BRCalculator::DWidth_Global
double DWidth_Global(genie::hnl::HNLDecayMode_t hnldm, const double M) const
Definition: HNLBRCalculator.cxx:248
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::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::KScale_MuonToNuAndElectron
double KScale_MuonToNuAndElectron(const double M) const
Definition: HNLBRCalculator.cxx:223
genie::hnl::BRCalculator::PiPi0EllForm
static double PiPi0EllForm(double *x, double *par)
Definition: HNLBRCalculator.cxx:472
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::utils::hnl::IsProdKinematicallyAllowed
bool IsProdKinematicallyAllowed(genie::hnl::HNLProd_t hnlprod)
genie::hnl::selector::__attribute__
static __attribute__((unused)) double fDecayGammas[]
genie::hnl::BRCalculator::DecayWidth
double DecayWidth(genie::hnl::HNLDecayMode_t hnldm) const
Definition: HNLBRCalculator.cxx:131
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::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
genie::hnl::BRCalculator::kscale_K3mu
std::map< double, double > kscale_K3mu
Definition: HNLBRCalculator.h:158
genie::hnl::BRCalculator::kscale_mu3e
std::map< double, double > kscale_mu3e
Definition: HNLBRCalculator.h:158
genie::hnl::BRCalculator::GetFormfactorF2
double GetFormfactorF2(double x) const
Definition: HNLBRCalculator.cxx:144
genie::hnl::ChannelCalculatorI
Pure abstract base class. Defines the ChannelCalculatorI interface to be implemented by BRCalculator ...
Definition: HNLChannelCalculatorI.h:39
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::hnl::BRCalculator::Configure
void Configure(const Registry &config)
Definition: HNLBRCalculator.cxx:41
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::hnl::BRCalculator::fCouplings
std::vector< double > fCouplings
Definition: HNLBRCalculator.h:147
genie::hnl::BRCalculator::GetFormfactorF1
double GetFormfactorF1(double x) const
Definition: HNLBRCalculator.cxx:136
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::KScale_PseudoscalarToLepton
double KScale_PseudoscalarToLepton(const double mP, const double M, const double ma) const
Definition: HNLBRCalculator.cxx:178
genie::kPdgMuon
const int kPdgMuon
Definition: PDGCodes.h:37
genie::utils::hnl::MassX
double MassX(double m1, double m2)
Definition: HNLKinUtils.h:32
genie::hnl::BRCalculator::Pi0Pi0NuForm
static double Pi0Pi0NuForm(double *x, double *par)
Definition: HNLBRCalculator.cxx:500
genie::utils::hnl::IsKinematicallyAllowed
bool IsKinematicallyAllowed(genie::hnl::HNLDecayMode_t hnldm, double Mhnl)
Generated on Wed Jul 17 2024 14:46:42 for GENIEGenerator by   doxygen 1.8.5