GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Namespaces | Functions | Variables
genie::utils::mec Namespace Reference

MEC utilities. More...

Namespaces

 gsl
 

Functions

double GetTmuCostFromq0q3 (double dq0, double dq3, double Enu, double lmass, double &tmu, double &cost, double &area)
 
bool GetTlCostlFromq0q3 (double q0, double q3, double Enu, double ml, double &Tl, double &costl)
 
bool Getq0q3FromTlCostl (double Tl, double costl, double Enu, double ml, double &q0, double &q3)
 
double J (double q0, double q3, double Enu, double ml)
 
double Qvalue (int targetpdg, int nupdg)
 
double OldTensorContraction (int nupdg, int targetpdg, double Enu, double Ml, double Tl, double costhl, int tensorpdg, genie::HadronTensorType_t tensor_type, char *tensor_model)
 
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)
 

Variables

const double Q0LimitMaxXSec = 2.
 
const double QMagLimitMaxXSec = 2.
 

Detailed Description

MEC utilities.

Author
Created:
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Function Documentation

double genie::utils::mec::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 at line 334 of file MECUtils.cxx.

References genie::Interaction::FSPrimLepton(), genie::Interaction::InitState(), genie::ProcessInfo::IsEM(), genie::Interaction::KinePtr(), genie::kKVctl, genie::kKVTl, genie::controls::kMinQ2Limit, genie::kPSTlctl, genie::kRfLab, LOG, pDEBUG, genie::InitialState::Probe(), genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), Q0LimitMaxXSec, genie::utils::kinematics::Q2(), QMagLimitMaxXSec, genie::Kinematics::SetKV(), and genie::XSecAlgorithmI::XSec().

Referenced by genie::QELEventGeneratorSuSA::ComputeMaxXSec(), and genie::MECGenerator::SelectSuSALeptonKinematics().

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 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
TParticlePDG * Probe(void) const
const double Q0LimitMaxXSec
Definition: MECUtils.h:83
static const double kMinQ2Limit
Definition: Controls.h:41
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 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
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
const double QMagLimitMaxXSec
Definition: MECUtils.h:86
double ProbeE(RefFrame_t rf) const
#define pDEBUG
Definition: Messenger.h:63
bool genie::utils::mec::Getq0q3FromTlCostl ( double  Tl,
double  costl,
double  Enu,
double  ml,
double &  q0,
double &  q3 
)

Definition at line 121 of file MECUtils.cxx.

Referenced by genie::utils::mec::gsl::d2Xsec_dTCosth::DoEval(), genie::SuSAv2MECPXSec::PairRatio(), genie::QELEventGeneratorSuSA::SelectLeptonKinematics(), genie::MECGenerator::SelectNSVLeptonKinematics(), genie::SuSAv2QELPXSec::XSec(), and genie::SuSAv2MECPXSec::XSec().

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 }
bool genie::utils::mec::GetTlCostlFromq0q3 ( double  q0,
double  q3,
double  Enu,
double  ml,
double &  Tl,
double &  costl 
)

Definition at line 82 of file MECUtils.cxx.

Referenced by genie::QELEventGeneratorSuSA::SelectLeptonKinematics().

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 }
double genie::utils::mec::GetTmuCostFromq0q3 ( double  dq0,
double  dq3,
double  Enu,
double  lmass,
double &  tmu,
double &  cost,
double &  area 
)

Definition at line 33 of file MECUtils.cxx.

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 }
double genie::utils::mec::J ( double  q0,
double  q3,
double  Enu,
double  ml 
)

Definition at line 147 of file MECUtils.cxx.

Referenced by genie::SKKinematicsGenerator::CalculateKin_AtharSingleKaon(), genie::utils::gsl::d3Xsec_dTldTkdCosThetal::DoEval(), genie::utils::gsl::d2Xsec_dQ2dv::DoEval(), genie::utils::kinematics::Jacobian(), genie::DFRKinematicsGenerator::ProcessEventRecord(), genie::QELKinematicsGenerator::ProcessEventRecord(), genie::DMELKinematicsGenerator::ProcessEventRecord(), genie::DISKinematicsGenerator::ProcessEventRecord(), genie::DMDISKinematicsGenerator::ProcessEventRecord(), genie::MECGenerator::SelectEmpiricalKinematics(), genie::SlowRsclCharmDISPXSecLO::XSec(), genie::PattonCEvNSPXSec::XSec(), genie::EmpiricalMECPXSec2015::XSec(), genie::RosenbluthPXSec::XSec(), genie::AhrensNCELPXSec::XSec(), genie::StrumiaVissaniIBDPXSec::XSec(), genie::ReinDFRPXSec::XSec(), genie::IMDAnnihilationPXSec::XSec(), genie::AivazisCharmPXSecLO::XSec(), genie::QPMDISPXSec::XSec(), genie::AhrensDMELPXSec::XSec(), genie::P33PaschosLalakulichPXSec::XSec(), genie::BergerSehgalCOHPiPXSec2015::XSec(), genie::QPMDMDISPXSec::XSec(), genie::LwlynSmithQELCCPXSec::XSec(), genie::NuElectronPXSec::XSec(), genie::PaisQELLambdaPXSec::XSec(), genie::DMElectronPXSec::XSec(), genie::BardinIMDRadCorPXSec::XSec(), genie::ReinSehgalCOHPiPXSec::XSec(), genie::NievesSimoVacasMECPXSec2016::XSec(), genie::BostedChristyEMPXSec::XSec(), genie::ReinSehgalRESPXSec::XSec(), genie::KovalenkoQELCharmPXSec::XSec(), genie::NievesQELCCPXSec::XSec(), genie::SmithMonizQELCCPXSec::XSec(), genie::BSKLNBaseRESPXSec2014::XSec(), and genie::MKSPPPXSec2020::XSec().

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 }
double genie::utils::mec::OldTensorContraction ( int  nupdg,
int  targetpdg,
double  Enu,
double  Ml,
double  Tl,
double  costhl,
int  tensorpdg,
genie::HadronTensorType_t  tensor_type,
char *  tensor_model 
)

Definition at line 208 of file MECUtils.cxx.

References genie::AlgFactory::GetAlgorithm(), genie::HadronTensorModelI::GetTensor(), genie::AlgFactory::Instance(), genie::constants::kGF2, genie::constants::kPi, LOG, pFATAL, genie::HadronTensorI::tt(), genie::HadronTensorI::tz(), genie::HadronTensorI::xx(), genie::HadronTensorI::xy(), and genie::HadronTensorI::zz().

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 }
virtual std::complex< double > zz(double q0, double q_mag) const =0
The tensor element .
virtual const HadronTensorI * GetTensor(int tensor_pdg, HadronTensorType_t type) const =0
#define pFATAL
Definition: Messenger.h:56
Abstract interface for an object that computes the elements ( , , etc.) and structure functions ( ...
#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
virtual std::complex< double > tz(double q0, double q_mag) const =0
The tensor element .
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
Creates hadron tensor objects for use in cross section calculations.
virtual std::complex< double > tt(double q0, double q_mag) const =0
The tensor element .
virtual std::complex< double > xx(double q0, double q_mag) const =0
The tensor element .
double genie::utils::mec::Qvalue ( int  targetpdg,
int  nupdg 
)

Definition at line 164 of file MECUtils.cxx.

References genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), LOG, and pFATAL.

Referenced by genie::utils::BindHitNucleon(), genie::SuSAv2MECPXSec::Qvalue(), genie::SuSAv2QELPXSec::XSec(), and genie::NievesSimoVacasMECPXSec2016::XSec().

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 }
#define pFATAL
Definition: Messenger.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96

Variable Documentation

const double genie::utils::mec::Q0LimitMaxXSec = 2.

Definition at line 83 of file MECUtils.h.

Referenced by GetMaxXSecTlctl().

const double genie::utils::mec::QMagLimitMaxXSec = 2.

Definition at line 86 of file MECUtils.h.

Referenced by GetMaxXSecTlctl().