StrumiaVissaniIBDPXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Corey Reed <cjreed \at nikhef.nl>
8  Nikhef - June 22, 2009
9 
10  For the class documentation see the corresponding header file.
11 */
12 //____________________________________________________________________________
13 
14 #include <TMath.h>
15 
23 
24 using namespace genie;
25 using namespace genie::constants;
26 
28 
29 //____________________________________________________________________________
32  fXSecIntegrator(0)
33 {
34 
35 }
36 //____________________________________________________________________________
38  XSecAlgorithmI("genie::StrumiaVissaniIBDPXSec", config),
39  fXSecIntegrator(0)
40 {
41 
42 }
43 //____________________________________________________________________________
45 {
46 
47 }
48 //____________________________________________________________________________
50  const Interaction * interaction, KinePhaseSpace_t kps) const
51 {
52  // compute the differential cross section ds/dt
53 
54  if(! this -> ValidProcess (interaction) ) return 0.;
55  if(! this -> ValidKinematics (interaction) ) return 0.;
56 
57  const InitialState & init_state = interaction -> InitState();
58  const double Ev = init_state.ProbeE(kRfHitNucRest);
59  const Target & target = init_state.Tgt();
60  const bool isProt = target.IsProton();
61  const Kinematics & kine = interaction->Kine();
62  const double q2 = kine.q2();
63  const double mp = (isProt) ? kProtonMass : kNeutronMass;
64  const double mp2 = (isProt) ? kProtonMass2 : kNeutronMass2;
65  const double mn2 = (isProt) ? kNeutronMass2 : kProtonMass2;
66 
67  // calculate s-u and s-m_nucleon^2 in nucleon rest frame
68  // need E_e
69  const double Ee = Ev + ( (q2 - mn2 + mp2) / (2.000*mp) );
70  const double sMinusU = ((2.000*mp*(Ev+Ee)) - kElectronMass2)
71  * (isProt ? 1.000 : -1.000);
72  const double sMinusMp2 = 2.000*mp*Ev;
73 
74 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
75  LOG("StrumiaVissani", pDEBUG) << "*** Ev = " << Ev
76  << ", q2 = " << q2
77  << ", Ee = " << Ee
78  << ", s-u = " << sMinusU
79  << ", s-Mp2 = " << sMinusMp2;
80 #endif
81 
82  double xsec = dSigDt(sMinusU, sMinusMp2, q2);
83 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
84  LOG("StrumiaVissani", pDEBUG) << "*** dSdt = " << xsec;
85 #endif
86 
87  // apply correction factors
88  const double rdcf = RadiativeCorr(Ee);
89  const double fscf = (isProt) ? 1.00 : FinalStateCorr(Ee);
90  xsec *= rdcf * fscf;
91 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
92  LOG("StrumiaVissani", pDEBUG) << "*** rad.corr. = " << rdcf
93  << ", fin.st.cor. = " << fscf
94  << ", xsec = " << xsec;
95 #endif
96 
97  //----- The algorithm computes dxsec/dt, t=q^2
98  // Check whether variable tranformation is needed
99  if(kps!=kPSq2fE && kps!=kPSQ2fE) {
100  const double J = utils::kinematics::Jacobian(interaction,kPSq2fE,kps);
101  xsec *= J;
102 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
103  LOG("StrumiaVissani", pDEBUG) << "*** Jacobian = " << J;
104 #endif
105  }
106 
107 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
108  LOG("StrumiaVissani", pDEBUG) << "*** xsec = " << xsec;
109 #endif
110  return xsec;
111 }
112 //____________________________________________________________________________
113 double StrumiaVissaniIBDPXSec::Integral(const Interaction * interaction) const
114 {
115  // Compute the total cross section for a free nucleon target
116 
117  assert(interaction!=0);
118  if(! this -> ValidProcess (interaction) ) return 0.;
119  if(! this -> ValidKinematics (interaction) ) return 0.;
120 
122  double xsec = fXSecIntegrator->Integrate(this, interaction);
123 
124  return xsec;
125 }
126 //____________________________________________________________________________
127 bool StrumiaVissaniIBDPXSec::ValidProcess(const Interaction * interaction) const
128 {
129  if(interaction->TestBit(kISkipProcessChk)) return true;
130 
131  // should be IBD and either nu_e + n or anu_e + p
132  if (interaction->ProcInfo().IsInverseBetaDecay()) {
133 
134  const InitialState & init_state = interaction -> InitState();
135  const Target & target = init_state.Tgt();
136  if ( (target.IsProton() && pdg::IsAntiNuE(init_state.ProbePdg())) || (target.IsNeutron() && pdg::IsNuE(init_state.ProbePdg()) ))
137  return true;
138  }
139  LOG("StrumiaVissani", pERROR) << "*** Should be IBD processes either nu_e + n or anu_e + p!";
140  return false;
141 }
142 //____________________________________________________________________________
144 {
145  if(interaction->TestBit(kISkipKinematicChk)) return true;
146 
147  const InitialState & init_state = interaction -> InitState();
148  const double Ev = init_state.ProbeE(kRfHitNucRest);
149  static const double Ecutoff = kNucleonMass / 2;
150  if (Ev > Ecutoff) {
151  LOG("StrumiaVissani", pERROR) << "*** Ev=" << Ev
152  << " is too large for VLE QEL!";
153  } else if (init_state.IsNuBarP() || init_state.IsNuN()) {
154  const KPhaseSpace & kps = interaction->PhaseSpace();
155  return kps.IsAboveThreshold();
156  }
157 
158  return false;
159 }
160 //____________________________________________________________________________
162 {
163  Algorithm::Configure(config);
164  this->LoadConfig();
165 }
166 //____________________________________________________________________________
168 {
169  Algorithm::Configure(config);
170  this->LoadConfig();
171 }
172 //____________________________________________________________________________
174 {
175 
176  // cabibbo angle
177  double cab ;
178  GetParam( "CabibboAngle", cab ) ;
179  const double cosCab = TMath::Cos(cab);
180  fCosCabibbo2 = cosCab*cosCab;
181 
182  // form factor params
183  GetParam( "QEL-FA0", fg1of0 ) ;
184  double ma, mv ;
185  GetParam( "QEL-Ma", ma ) ;
186  GetParam( "QEL-Mv", mv ) ;
187  fMa2 = ma*ma;
188  fMv2 = mv*mv;
189 
190  // magnetic moments
191 
192  double mup, mun ;
193  GetParam( "AnomMagnMoment-P", mup );
194  GetParam( "AnomMagnMoment-N", mun );
195  fNucleonMMDiff = (mup - 1.000) - mun; // *anamolous* mag. mom. diff.
196 
197  // numeric
198  int epmag ;
199  GetParam("EpsilonMag", epmag ) ;
200  fEpsilon = TMath::Power(10.000, -1.000 * static_cast<double>(epmag) );
201 
202  LOG("StrumiaVissani", pINFO) << "*** USING: cos2(Cabibbo)="
203  << fCosCabibbo2
204  << ", g1(0)=" << fg1of0
205  << ", Ma^2=" << fMa2
206  << ", Mv^2=" << fMv2
207  << ", xi=" << fNucleonMMDiff
208  << ", epsilon=" << fEpsilon;
209 
210  // load XSec Integrator
212  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
214 
215 }
216 //____________________________________________________________________________
217 double StrumiaVissaniIBDPXSec::dSigDt(const double sMinusU,
218  const double sMinusMnuc,
219  const double t) const
220 {
221  // return the differential cross section dS/dt. eqn 3 from reference
222  // t = q^2
223  //
224  // for anue + p -> e+ + n , sMinusU = s - u and sMinusMnuc = s - m_p^2
225  // for nue + n -> e- + p , sMinusU = u - s and sMinusMnuc = s - m_n^2
226  // where s = (p_nu + p_p)^2 and u = (p_nu - p_n)^2
227 
228  const double numer = kGF2 * fCosCabibbo2;
229  const double denom = (2.000 * kPi) * sMinusMnuc*sMinusMnuc;
230  assert(TMath::Abs(denom) > fEpsilon); // avoid divide by zero
231 
232  return ( (numer / denom) * MtxElm(sMinusU, t) );
233 }
234 //____________________________________________________________________________
235 double StrumiaVissaniIBDPXSec::MtxElm(const double sMinusU,
236  const double t) const
237 {
238  // return the square of the matrix element. eqn 5 from the reference paper
239  // |M^2| = A(t) - (s-u)B(t) + (s-u)^2 C(t)
240  //
241  // factor 16 comes from eqn 6
242  //
243  // for anue + p -> e+ + n , sMinusU = s - u
244  // for nue + n -> e- + p , sMinusU = u - s
245  // where s = (p_nu + p_p)^2 and u = (p_nu - p_n)^2
246 
247  // store multiply used variables to reduce number of calculations
248  const double t4m = t / k4NucMass2;
249  const double t2 = t * t;
250 
251  const double fdenomB = 1.000 - (t / fMv2);
252  const double fdenom = (1.000 - t4m)*fdenomB*fdenomB;
253  const double g1denom = 1.000 - (t / fMa2);
254  const double g2denom = kPionMass2 - t;
255  assert(TMath::Abs(fdenom) > fEpsilon); // avoid divide by zero
256  assert(TMath::Abs(g1denom) > fEpsilon); // avoid divide by zero
257  assert(TMath::Abs(g2denom) > fEpsilon); // avoid divide by zero
258 
259  const double f1 = (1.000 - ( (1.000+fNucleonMMDiff) * t4m)) / fdenom;
260  const double f2 = fNucleonMMDiff / fdenom;
261  const double g1 = fg1of0 / (g1denom*g1denom);
262  const double g2 = (2.000 * kNucleonMass2 * g1) / g2denom;
263 
264  const double f12 = f1 * f1;
265  const double f124 = 4.000 * f12;
266  const double f22 = f2 * f2;
267  const double g12 = g1 * g1;
268  const double g124 = 4.000 * g12;
269  const double g22 = g2 * g2;
270  const double g224meM2 = 4.000 * g22 * kElectronMass2 / kNucleonMass2;
271 
272  const double g1cFsumR = g1 * (f1+f2);
273  const double f1cf2 = f1 * f2;
274  const double g1cg2 = g1 * g2;
275  const double f1cf2R8 = f1cf2 * 8.000;
276  const double g1cg2R16me = g1cg2 * 16.000 * kElectronMass2;
277 
278 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
279  LOG("StrumiaVissani", pDEBUG) << "*** t = " << t
280  << ", g2 = " << g2
281  << ", g1 = " << g1
282  << ", f1 = " << f1
283  << ", f2 = " << f2;
284 #endif
285 
286  // ok - now calculate the terms of the matrix element
287  const double mat = MAterm(t, t2, f124, f22, g124,
288  g224meM2, f1cf2R8, g1cg2R16me, g1cFsumR);
289  const double mbt = MBterm(t, f1cf2, g1cg2, g1cFsumR, f22);
290  const double mct = MCterm(t, f124, f22, g124);
291 
292  const double M2 = mat - (sMinusU * (mbt - (sMinusU * mct)));
293 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
294  LOG("StrumiaVissani", pDEBUG) << "*** Matrix element = " << (M2 / 16.000)
295  << " [16A=" << mat << "] [16B=" << mbt
296  << "] [16C=" << mct << "]";
297 #endif
298 
299  return ( M2 / 16.000 );
300 }
301 //____________________________________________________________________________
302 double StrumiaVissaniIBDPXSec::MAterm(const double t,
303  const double t2,
304  const double f124,
305  const double f22,
306  const double g124,
307  const double g224meM2,
308  const double f1cf2R8,
309  const double g1cg2R16me,
310  const double g1cFsumR)
311 {
312  // return A(t) in the matrix element calculation
313  // from eqn 6 (without the factor 16)
314  //
315  // for speed purposes, no check that |t|^2 = t2 is performed
316 
317  // store multiply used terms
318  const double tmme2 = t - kElectronMass2;
319  const double tpme2 = t + kElectronMass2;
320 
321  double r1 = f124 * ( k4NucMass2 + tpme2);
322  r1 += g124 * (-k4NucMass2 + tpme2);
323  r1 += f22 * ( (t2/kNucleonMass2) + 4.000*tpme2 );
324  r1 += g224meM2 * t;
325  r1 += f1cf2R8 * (2.000*t + kElectronMass2);
326  r1 += g1cg2R16me;
327  r1 *= tmme2;
328 
329  double r2 = (f124 + (t * (f22 / kNucleonMass2))) * (k4NucMass2 + tmme2);
330  r2 += g124 * (k4NucMass2 - tmme2);
331  r2 += g224meM2 * tmme2;
332  r2 += f1cf2R8 * (2.000*t - kElectronMass2);
333  r2 += g1cg2R16me;
334  r2 *= kNucMassDiff2;
335 
336  const double r3 = 32.000 * kElectronMass2 * kNucleonMass * kNucMassDiff
337  * g1cFsumR;
338 
339  return ( r1 - r2 - r3 );
340 }
341 //____________________________________________________________________________
342 double StrumiaVissaniIBDPXSec::MBterm(const double t,
343  const double f1cf2,
344  const double g1cg2,
345  const double g1cFsumR,
346  const double f22)
347 {
348  // return C(t) in the matrix element calculation
349  // from eqn 6 (without the factor 16)
350 
351  double bterm = 16.000 * t * g1cFsumR;
352  bterm += ( k4EleMass2 * kNucMassDiff
353  * (f22 + (f1cf2 + 2.000*g1cg2)) ) / kNucleonMass2;
354  return bterm;
355 }
356 //____________________________________________________________________________
357 double StrumiaVissaniIBDPXSec::MCterm(const double t,
358  const double f124,
359  const double f22,
360  const double g124)
361 {
362  // return C(t) in the matrix element calculation
363  // from eqn 6 (without the factor 16)
364 
365  return ( f124 + g124 - (t * ( f22 / kNucleonMass2 )) );
366 }
367 //____________________________________________________________________________
368 double StrumiaVissaniIBDPXSec::RadiativeCorr(const double Ee) const
369 {
370  // radiative correction to the cross section. eqn 14 from the reference
371  // only valid for Ev<<m_p!
372 
373  assert(Ee > fEpsilon); // must be non-zero and positive
374  double rc = 6.000 + (1.500 * TMath::Log(kProtonMass / (2.000*Ee)));
375  rc += 1.200 * TMath::Power((kElectronMass / Ee), 1.500);
376  rc *= kAem / kPi;
377  rc += 1.000;
378  return rc;
379 
380 }
381 //____________________________________________________________________________
382 double StrumiaVissaniIBDPXSec::FinalStateCorr(const double Ee) const
383 {
384  // Sommerfeld factor; correction for final state interactions.
385  // eqn 15 of the reference
386 
387  assert(Ee > fEpsilon); // must be non-zero and positive
388  const double eta = 2.000*kPi*kAem
389  / TMath::Sqrt(1.000 - (kElectronMass2 / (Ee*Ee)));
390  const double expn = TMath::Exp(-1.000 * eta);
391  assert(expn < 1.000);
392  return ( eta / (1.000 - expn) );
393 }
394 //____________________________________________________________________________
395 
Cross Section Calculation Interface.
const double kPi
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
const XML_Char * target
Definition: expat.h:268
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
Cross Section Integrator Interface.
static const double kNucleonMass
Definition: Constants.h:78
bool IsNeutron(void) const
Definition: Target.cxx:284
static double MAterm(const double t, const double t2, const double f124, const double f22, const double g124, const double g224meM2, const double f1cf2R8, const double g1cg2R16me, const double g1cFsumR)
double RadiativeCorr(const double Ee) const
bool IsNuN(void) const
is neutrino + neutron?
static double MCterm(const double t, const double f124, const double f22, const double g124)
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
bool IsNuE(int pdgc)
Definition: PDGUtils.cxx:150
bool IsInverseBetaDecay(void) const
static double MBterm(const double t, const double f1cf2, const double g1cg2, const double g1cFsumR, const double f22)
ClassImp(StrumiaVissaniIBDPXSec) StrumiaVissaniIBDPXSec
static const double k4EleMass2
Definition: Constants.h:31
Definition: config.py:1
Float_t f2
TRandom3 * r3
enum genie::EKinePhaseSpace KinePhaseSpace_t
An implementation of the neutrino - (free) nucleon [inverse beta decay] cross section, valid from the threshold energy (1.806MeV) up to hundreds of MeV. Currently cut off at 1/2 nucleon mass. Based on the Strumia/Vissani paper Phys.Lett.B564:42-54,2003.
static const double kElectronMass
Definition: Constants.h:71
static const double kAem
Definition: Constants.h:57
Double_t q2[12][num]
Definition: f2_nu.C:137
Summary information for an interaction.
Definition: Interaction.h:56
double q2(bool selected=false) const
Definition: Kinematics.cxx:151
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double FinalStateCorr(const double Ee) const
static const double kElectronMass2
Definition: Constants.h:84
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
Kinematical phase space.
Definition: KPhaseSpace.h:34
static const double kNeutronMass
Definition: Constants.h:77
TString mp
Definition: loadincs.C:4
Float_t f1
double MtxElm(const double sMinusU, const double t) const
static const double kNucMassDiff
Definition: Constants.h:32
double Integral(const Interaction *i) const
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
#define pINFO
Definition: Messenger.h:63
double t2
static const double kNucMassDiff2
Definition: Constants.h:33
static const double kNucleonMass2
Definition: Constants.h:90
Double_t xsec[nknots]
Definition: testXsec.C:47
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
static const double k4NucMass2
Definition: Constants.h:30
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double dSigDt(const double sMinusU, const double sMinusMnuc, const double t) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool IsNuBarP(void) const
is anti-neutrino + proton?
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void Configure(const Registry &config)
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
static const double kProtonMass2
Definition: Constants.h:88
static const double kGF2
Definition: Constants.h:60
static const double kNeutronMass2
Definition: Constants.h:89
bool IsProton(void) const
Definition: Target.cxx:279
double ProbeE(RefFrame_t rf) const
bool IsAntiNuE(int pdgc)
Definition: PDGUtils.cxx:165
static const double kProtonMass
Definition: Constants.h:76
Eigen::MatrixXd mat
const XSecIntegratorI * fXSecIntegrator
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
static const double kPionMass2
Definition: Constants.h:87
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353