GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MECUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6 
7  For documentation see the corresponding header file.
8 
9 */
10 //____________________________________________________________________________
11 
12 #include <TMath.h>
13 #include <TLorentzVector.h>
14 
28 
29 using namespace genie;
30 using namespace genie::constants;
31 
32 //______________________________________________________________________
34  double dq0, double dq3, double Enu, double lmass,
35  double &tmu, double &cost, double &area)
36 {
37  tmu = Enu - dq0 - lmass;
38  if(tmu < 0.0){
39  cost = -999;
40  tmu = -999;
41  area = 0.0;
42  return -999;
43  }
44 
45  double thisE = tmu + lmass;
46  double thisp = sqrt( thisE * thisE - lmass * lmass);
47  double numerator = Enu * Enu + thisp * thisp - dq3 * dq3;
48  double denominator = 2.0 * thisp * Enu;
49  if(denominator <= 0.0 ){
50  cost = 0.0;
51  if(denominator < 0.0){
52  return -999;
53  }
54  }
55  else cost = numerator / denominator;
56 
57  if(TMath::Abs(numerator) > TMath::Abs(denominator)){
58  cost = -999;
59  tmu = -999;
60  area = 0.0;
61  return -999;
62  }
63 
64  // xCrossSect is not yet in the right units for this particular case.
65  // need areaElement to go dsigma/dTmudcost to dsigma/dq0dq3
66  // Recompute the area element jacobian
67 
68  // dT/dq0 x dC/dq3 - dT/dq3 x dC/dq0
69  double areaElement = 0.0;
70  //double veryCloseToZero = 0.000001; // in GeV, this should be safe.
71  double insqrt = 0.0;
72  numerator = 0.0;
73  insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - lmass * lmass;
74  numerator = dq3 / Enu;
75  if(insqrt < 0.0)areaElement=0.0;
76  else areaElement = numerator / TMath::Sqrt(insqrt);
77  area = areaElement;
78 
79  return 0;
80 }
81 //______________________________________________________________________
83  double q0, double q3, double Enu, double ml, double& Tl, double& costl)
84 {
85  Tl = Enu - q0 - ml;
86  if(Tl < 0.) {
87  costl = -999;
88  Tl = -999;
89  return false;
90  }
91 
92  double El = Tl + ml;
93  double plsq = El * El - ml * ml;
94  if(plsq < 0.) {
95  costl = -999;
96  Tl = -999;
97  return false;
98  }
99  double pl = TMath::Sqrt(plsq);
100  double numerator = Enu * Enu + pl * pl - q3 * q3;
101  double denominator = 2.0 * pl * Enu;
102  if(denominator <= 0.0) {
103  costl = 0.0;
104  if(denominator < 0.0) {
105  return false;
106  }
107  }
108  else {
109  costl = numerator / denominator;
110  }
111 
112  if(TMath::Abs(numerator) > TMath::Abs(denominator)){
113  costl = -999;
114  Tl = -999;
115  return false;
116  }
117 
118  return true;
119 }
120 //______________________________________________________________________
122  double Tl, double costl, double Enu, double ml, double& q0, double& q3)
123 {
124  q0 = -999;
125  q3 = -999;
126 
127  double El = Tl + ml;
128  q0 = Enu - El;
129  if(q0 < 0) {
130  q0 = -999;
131  q3 = -999;
132  return false;
133  }
134 
135  double pl = TMath::Sqrt(El * El - ml * ml);
136  double q3sq = pl * pl + Enu * Enu - 2.0 * pl * Enu * costl;
137  if(q3sq < 0) {
138  q0 = -999;
139  q3 = -999;
140  return false;
141  }
142  q3 = TMath::Sqrt(q3sq);
143 
144  return true;
145 }
146 //______________________________________________________________________
148  double dq0, double dq3, double Enu, double ml)
149 {
150  // dT/dq0 x dC/dq3 - dT/dq3 x dC/dq0
151  double area = 0.0;
152  double insqrt = 0.0;
153  insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - ml * ml;
154  double numerator = dq3 / Enu;
155  if(insqrt < 0.0) {
156  area=0.0;
157  }
158  else {
159  area = numerator / TMath::Sqrt(insqrt);
160  }
161  return area;
162 }
163 //______________________________________________________________________
164 double genie::utils::mec::Qvalue(int targetpdg, int nupdg)
165 {
166 // The had tensor calculation requires parameters that describe the
167 // nuclear density. For the Valencia model they are taken by
168 // C. Garcia-Recio, Nieves, Oset Nucl.Phys.A 547 (1992) 473-487 Table I
169 // and simplifies that the p and n density are not necessarily the same?
170 // Standard tables such as deVries et al are similar ~5%
171 // This is the only nuclear input on the hadron side of the Hadron Tensor.
172 // and is what makes PseudoPb and PseudoFe different than Ni56 and Rf208
173 //
174 
175  int nearpdg = targetpdg;
176 
177  if(nupdg > 0){
178  // get Z+1 for CC interaction e.g. 1000210400 from 1000200400
179  nearpdg = targetpdg + 10000;
180  } else {
181  // get Z-1 for CC interaction e.g. 1000210400 from 1000200400
182  nearpdg = targetpdg - 10000;
183  }
184 
185  TParticlePDG *partf = PDGLibrary::Instance()->Find(nearpdg);
186  TParticlePDG *parti = PDGLibrary::Instance()->Find(targetpdg);
187 
188  if(NULL == partf || NULL == parti){
189  // maybe not every valid nucleus in the table has Z+1 or Z-1
190  // for example, Ca40 did not.
191  LOG("MECUtils", pFATAL)
192  << "Can't get qvalue, nucleus " << targetpdg << " or " << nearpdg
193  << " is not in the table of nuclei in /data/evgen/pdg ";
194  exit(-1);
195  }
196 
197  double massf = partf->Mass();
198  double massi = parti->Mass();
199 
200  // keep this here as a reminder of what not to do
201  // +/- emass; // subtract electron from Z+1, add it to Z-1
202  // the lookup table gives the nucleus mass, not the atomic
203  // mass, so don't need this.
204 
205  return massf - massi; // in GeV.
206 }
207 //______________________________________________________________________
209  int nupdg, int targetpdg,
210  double Enu, double Ml, double Tl, double costhl,
211  int tensorpdg,
212  HadronTensorType_t tensor_type,
213  char* tensor_model)
214 {
215  TLorentzVector v4lep;
216  TLorentzVector v4Nu(0,0,Enu,Enu); // assuming traveling along z:
217  TLorentzVector v4q;
218  double q0nucleus;
219  double facconv = 0.0389391289e15; // std::pow(0.19733,2)*1e15;
220 
221  double myQvalue = 0.0;
222 
223  // Angles
224  double sinthl = 1. - costhl * costhl;
225  if(sinthl < 0.0) sinthl = 0.0;
226  else sinthl = TMath::Sqrt(sinthl);
227 
228  double Cosh = TMath::Cos(TMath::ACos(costhl)/2.);
229  double Sinh = TMath::Sin(TMath::ACos(costhl)/2.);
230 
231  // Lepton
232  v4lep.SetE( Tl + Ml );
233  // energy transfer from the lepton
234  double q0 = v4Nu.E() - v4lep.E();
235  // energy transfer that actually gets to the nucleons
236  q0nucleus = q0 - myQvalue;
237 
238  //Get q3
239  double pl = TMath::Sqrt(v4lep.E() * v4lep.E() - Ml * Ml);
240  double q3sq = pl * pl + Enu * Enu - 2.0 * pl * Enu * costhl;
241  double q3 = sqrt(q3sq);
242 
243  // Define some calculation placeholders
244  double part1, part2;
245  double modkprime ;
246  double w1, w2, w3, w4, w5;
247  double wtotd[5];
248 
249  double xsec = 0;
250  if (q0nucleus <= 0){
251  //nothing transfered to nucleus thus no 2p2h
252  xsec = 0.;
253  } else {
254  // energy was transfered to the nucleus
255  modkprime = std::sqrt(TMath::Power(v4lep.E(),2)-TMath::Power(Ml,2));
256  v4lep.SetX(modkprime*sinthl);
257  v4lep.SetY(0);
258  v4lep.SetZ(modkprime*costhl);
259 
260  //q: v4q = v4Nu - v4lep;
261  v4q.SetE(q0nucleus);
262  v4q.SetX(v4Nu.X() - v4lep.X());
263  v4q.SetY(v4Nu.Y() - v4lep.Y());
264  v4q.SetZ(v4Nu.Z() - v4lep.Z());
265 
266 
267  // Get the appropriate hadron tensor model object
268  const genie::HadronTensorModelI* ht_model
269  = dynamic_cast<const genie::HadronTensorModelI*>(
270  genie::AlgFactory::Instance()->GetAlgorithm( tensor_model, "Default" ));
271 
272  const LabFrameHadronTensorI* tensor_table
273  = dynamic_cast<const LabFrameHadronTensorI*>( ht_model->GetTensor(targetpdg,
274  tensor_type) );
275 
276  double W00 = (tensor_table->tt(q0nucleus,q3)).real();
277  double W0Z = (tensor_table->tz(q0nucleus,q3)).real();
278  double WXX = (tensor_table->xx(q0nucleus,q3)).real();
279  double WXY = -1.0*(tensor_table->xy(q0nucleus,q3)).imag();
280  double WZZ = (tensor_table->zz(q0nucleus,q3)).real();
281 
282  w1=WXX/2.;
283  w2=(W00+WXX+(q0*q0/(v4q.Vect().Mag()*v4q.Vect().Mag())
284  *(WZZ-WXX))
285  -(2.*q0/v4q.Vect().Mag()*W0Z))/2.;
286  w3=WXY/v4q.Vect().Mag();
287  w4=(WZZ-WXX)/(2.*v4q.Vect().Mag()*v4q.Vect().Mag());
288  w5=(W0Z-(q0/v4q.Vect().Mag()*(WZZ-WXX)))/v4q.Vect().Mag();
289  //w6 we have no need for w6, noted at the end of IIA.
290 
291  // adjust for anti neutrinos
292 
293  if (nupdg < 0) w3 = -1. * w3;
294 
295  // calculate cross section, in parts
296  double xw1 = w1*costhl;
297  double xw2 = w2/2.*costhl;
298  double xw3 = w3/2.*(v4lep.E()+modkprime-(v4lep.E()+v4Nu.E())*costhl);
299  double xw4 = w4/2.*(Ml*Ml*costhl+2.*v4lep.E()*(v4lep.E()+modkprime)*Sinh*Sinh);
300  double xw5 = w5*(v4lep.E()+modkprime)/2.;
301  part1 = xw1 - xw2 + xw3 + xw4 - xw5;
302 
303  double yw1 = 2.*w1*Sinh*Sinh;
304  double yw2 = w2*Cosh*Cosh;
305  double yw3 = w3*(v4lep.E()+v4Nu.E())*Sinh*Sinh;
306  double yw4 = Ml*Ml*part1/(v4lep.E()*(v4lep.E()+modkprime));
307  part2 = yw1 + yw2 - yw3 + yw4;
308 
309  xsec = modkprime*v4lep.E()*kGF2*2./kPi*part2*facconv;
310 
311  if( ! (xsec >= 0.0) ){
312  // Traps negative numbers and nan's.
313  // Sometimes costhl is just over 1.0 due to fluke or numerical
314  // precision, so sinthl would be undefined.
315  xsec = 0;
316  }
317  }
318  if((Enu>1.15 && Enu<1.25) && (v4q.Vect().Mag()>0.615 && v4q.Vect().Mag()<0.625) && (q0<0.500 && q0>0.490))
319  LOG("MECUtils", pFATAL)//SB
320  << "Got xsec = " << xsec
321  << "\n " << part1 << " " << part2
322  << "\n w[] : " << w1 << ", " << w2 << ", " << w3 << ", " << w4 << ", " << w5
323  << "\n wtotd[] : " << wtotd[0] << ", " << wtotd[1] << ", " << wtotd[2] << ", " << wtotd[3] << ", " << wtotd[4]
324  << "\n vec " << v4q.Vect().Mag() << ", " << q0 << ", " << v4q.Px() << ", " << v4q.Py() << ", " << v4q.Pz()
325  << "\n v4qX " << v4Nu.X() << ", " << v4lep.X() << ", " << costhl << ", " << sinthl << ", " << modkprime
326  << "\n input " << Enu << ", " << Ml << ", " << Tl << ", " <<costhl << ", " << v4q.Vect().Mag() << ", " << q0;
327 
328  //LOG("MECUtils", pFATAL) << "xsec(Enu = " << Enu << " GeV, Ml = " << Ml << " GeV; " << "Tl = " << Tl << " GeV, costhl = " << costhl << ") = " << xsec << " x 1E-41 cm^2";
329 
330  //return (xsec * (1.0E-39 * units::cm2));
331  return (xsec);
332 }
333 //___________________________________________________________________________
335  const Interaction& inter, const double tolerance, const double safety_factor,
336  const int max_n_layers )
337 {
338  // Clone the input interaction so that we can modify the Tl and ctl values
339  // without messing anything up
340  Interaction* interaction = new Interaction( inter );
341 
342  // Choose the appropriate minimum Q^2 value based on the interaction
343  // mode (this is important for EM interactions since the differential
344  // cross section blows up as Q^2 --> 0)
345  double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit
346  if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics
347  ::electromagnetic::kMinQ2Limit; // EM limit
348 
349  const double Enu = interaction->InitState().ProbeE( kRfLab );
350  const double ProbeMass = interaction->InitState().Probe()->Mass();
351  const double LepMass = interaction->FSPrimLepton()->Mass();
352 
353  // Determine the bounds for the current layer being scanned.
354 
355  // Maximum energy transfer to consider
356  const double Q0Max = std::min( utils::mec::Q0LimitMaxXSec, Enu );
357 
358  // Maximum momentum transfer to consider
359  const double QMagMax = std::min( utils::mec::QMagLimitMaxXSec, 2.*Enu );
360 
361  // Translate these bounds into limits for the lepton
362  // kinetic energy and Q^2
363  const double TMax = std::max( 0., Enu - LepMass );
364  double CurTMin = std::max( 0., TMax - Q0Max );
365  double CurTMax = TMax;
366  double CurQ2Min = Q2min;
367  // Overshoots the true maximum Q^2, but not so severely that it's expected to
368  // be a problem.
369  double CurQ2Max = QMagMax * QMagMax;
370 
371  // This is based on the technique used to compute the maximum differential
372  // cross section for rejection sampling in QELEventGenerator. A brute-force
373  // scan over the two-dimensional kPSTlctl phase space is used to find the
374  // maximum cross section. Multiple layers are used to "zoom in" on the
375  // vicinity of the maximum. Tighter and tighter layers are made until
376  // the maximum number of iterations is reached or the xsec stabilizes
377  // to within a given tolerance.
378  const int N_T = 10; // number of kinetic energy steps per layer
379  const int N_Q2 = 10; // number of scattering cosine steps per layer
380  double T_at_xsec_max = CurTMax;
381  double Q2_at_xsec_max = CurQ2Max;
382  double cur_max_xsec = -1.;
383 
384  for ( int ilayer = 0; ilayer < max_n_layers; ++ilayer ) {
385 
386  double last_layer_xsec_max = cur_max_xsec;
387 
388  double T_increment = ( CurTMax - CurTMin ) / ( N_T - 1 );
389  double Q2_increment = ( CurQ2Max - CurQ2Min ) / ( N_Q2 - 1 );
390 
391  // Scan through the 2D phase space using the bounds of the current layer
392  for ( int iT = 0; iT < N_T; ++iT ) {
393  double T = CurTMin + iT * T_increment;
394  for ( int iQ2 = 0; iQ2 < N_Q2; ++iQ2 ) {
395  double Q2 = CurQ2Min + iQ2 * Q2_increment;
396 
397  // Compute the scattering cosine at fixed Tl given a value of Q^2.
398  double Elep = T + LepMass;
399  double pnu = std::sqrt( std::max(0., Enu*Enu - ProbeMass*ProbeMass) );
400  double plep = std::sqrt( std::max(0., std::pow(T + LepMass, 2)
401  - LepMass*LepMass) );
402 
403  double Costh = plep==0?0:( 2.*Enu*Elep - ProbeMass*ProbeMass - LepMass*LepMass
404  - Q2 ) / ( 2. * pnu * plep );
405  // Respect the bounds of the cosine function.
406  Costh = std::min( std::max(-1., Costh), 1. );
407 
408  // Update the values of Tl and ctl in the interaction, then
409  // compute the differential cross section
410  interaction->KinePtr()->SetKV( kKVTl, T );
411  interaction->KinePtr()->SetKV( kKVctl, Costh );
412  double xsec = xsec_model.XSec( interaction, kPSTlctl );
413 
414  // If the current differential cross section exceeds the previous
415  // maximum value, update the stored value and information about where it
416  // lies in the kPSTlctl phase space.
417  if ( xsec > cur_max_xsec ) {
418  T_at_xsec_max = T;
419  Q2_at_xsec_max = Q2;
420  cur_max_xsec = xsec;
421  }
422 
423  } // Done with cos(theta) scan
424  } // Done with lepton kinetic energy scan
425 
426  LOG("MEC", pDEBUG) << "Layer " << ilayer << " has max xsec = "
427  << cur_max_xsec << " for T = " << T_at_xsec_max << ", Q2 = "
428  << Q2_at_xsec_max;
429 
430  // ** Calculate the range for the next layer **
431 
432  // Don't let the minimum kinetic energy fall below zero
433  CurTMin = std::max( 0., T_at_xsec_max - T_increment );
434 
435  // We know a priori that the maximum cross section should occur near the
436  // maximum lepton kinetic energy, so keep the old maximum value when
437  // recalculating the new range. It was originally set to something near the
438  // maximum allowed value, so this should zero in on the appropriate region
439  // in future scans. The updated maximum from the old code is shown below,
440  // but due to problems with round-off in the kinematic limits, it can
441  // sometimes miss the edge of the phase space if the T_increment value is
442  // coarse enough. - S. Gardiner, 1 July 2020
443  CurTMax = TMax; // T_at_xsec_max + T_increment;
444 
445  // We use the same idea here. A priori, we know that the maximum cross section
446  // should lie near the minimum Q^2 value. Round-off gives us the same problems
447  // in finding the low Q^2 edge of phase space, so just keep the minimum equal
448  // to the cut for the next scan. The old assignment (which mostly worked
449  // but had problems for a coarse Q2_increment value) is commented out below.
450  // - S. Gardiner, 1 July 2020
451  CurQ2Min = Q2min; // std::max( Q2min, Q2_at_xsec_max - Q2_increment );
452  CurQ2Max = Q2_at_xsec_max + Q2_increment;
453 
454  // If the xsec has stabilized within the requested tolerance, then
455  // stop adding more layers
456  double improvement_factor = cur_max_xsec / last_layer_xsec_max;
457  if ( ilayer && std::abs(improvement_factor - 1.) < tolerance ) break;
458 
459  } // Done with iterations over layers
460 
461 
462  // We're done. Delete the cloned interaction object before returning the
463  // estimate of the maximum differential cross section.
464  delete interaction;
465 
466  double XSecMax = cur_max_xsec * safety_factor;
467  return XSecMax;
468 }
469 //___________________________________________________________________________
471  const double Enu, const double LepMass, const double Factor ) :
472 ROOT::Math::IBaseFunctionMultiDim(),
473 fModel(m),
474 fInteraction(i),
475 fEnu(Enu),
476 fLepMass(LepMass),
477 fFactor(Factor)
478 {
479 
480 }
481 //____________________________________________________________________________
483 {
484 
485 }
486 //____________________________________________________________________________
488 {
489  return 2;
490 }
491 //____________________________________________________________________________
492 double genie::utils::mec::gsl::d2Xsec_dTCosth::DoEval(const double * xin) const
493 {
494 // inputs:
495 // T [GeV]
496 // cos(theta)
497 // outputs:
498 // differential cross section (hbar=c=1 units)
499 //
500 
501  double T = xin[0];
502  double costh = xin[1];
503 
504  Kinematics * kinematics = fInteraction.KinePtr();
505  kinematics->SetKV(kKVTl, T);
506  kinematics->SetKV(kKVctl, costh);
507 
508  double Q0 = 0 ;
509  double Q3 = 0 ;
510  genie::utils::mec::Getq0q3FromTlCostl(T, costh, fEnu, fLepMass, Q0, Q3);
511 
512  kinematics ->SetKV(kKVQ0, Q0) ;
513  kinematics ->SetKV(kKVQ3, Q3) ;
514 
515  double xsec = fModel->XSec( &fInteraction, kPSTlctl);
516  return fFactor * xsec;
517 
518 }
519 //____________________________________________________________________________
520 ROOT::Math::IBaseFunctionMultiDim *
522 {
523  return
524  new genie::utils::mec::gsl::d2Xsec_dTCosth(fModel,fInteraction, fEnu, fLepMass, fFactor );
525 }
526 //____________________________________________________________________________
527 
Cross Section Calculation Interface.
virtual std::complex< double > zz(double q0, double q_mag) const =0
The tensor element .
bool GetTlCostlFromq0q3(double q0, double q3, double Enu, double ml, double &Tl, double &costl)
Definition: MECUtils.cxx:82
virtual const HadronTensorI * GetTensor(int tensor_pdg, HadronTensorType_t type) const =0
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
#define pFATAL
Definition: Messenger.h:56
double GetMaxXSecTlctl(const XSecAlgorithmI &xsec_model, const Interaction &inter, const double tolerance=0.01, const double safety_factor=1.2, const int max_n_layers=100)
Definition: MECUtils.cxx:334
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
TParticlePDG * Probe(void) const
Abstract interface for an object that computes the elements ( , , etc.) and structure functions ( ...
const double Q0LimitMaxXSec
Definition: MECUtils.h:83
enum genie::HadronTensorType HadronTensorType_t
static const double kMinQ2Limit
Definition: Controls.h:41
double Qvalue(int targetpdg, int nupdg)
Definition: MECUtils.cxx:164
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual std::complex< double > xy(double q0, double q_mag) const =0
The tensor element .
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
bool Getq0q3FromTlCostl(double Tl, double costl, double Enu, double ml, double &q0, double &q3)
Definition: MECUtils.cxx:121
unsigned int NDim(void) const
Definition: MECUtils.cxx:487
double DoEval(const double *xin) const
Definition: MECUtils.cxx:492
virtual std::complex< double > tz(double q0, double q_mag) const =0
The tensor element .
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
d2Xsec_dTCosth(const XSecAlgorithmI *m, const Interaction &i, const double Enu, const double LepMass, const double Factor=1.)
Definition: MECUtils.cxx:470
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
Definition: MECUtils.cxx:521
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
Creates hadron tensor objects for use in cross section calculations.
const double QMagLimitMaxXSec
Definition: MECUtils.h:86
double GetTmuCostFromq0q3(double dq0, double dq3, double Enu, double lmass, double &tmu, double &cost, double &area)
Definition: MECUtils.cxx:33
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
double OldTensorContraction(int nupdg, int targetpdg, double Enu, double Ml, double Tl, double costhl, int tensorpdg, genie::HadronTensorType_t tensor_type, char *tensor_model)
Definition: MECUtils.cxx:208
double ProbeE(RefFrame_t rf) const
virtual std::complex< double > tt(double q0, double q_mag) const =0
The tensor element .
static constexpr double m
Definition: Units.h:71
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
virtual std::complex< double > xx(double q0, double q_mag) const =0
The tensor element .
#define pDEBUG
Definition: Messenger.h:63