Functions
genie::utils::mec Namespace Reference

MEC utilities. More...

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 TensorContraction (const Interaction *interaction, int tensor_pdg, MECHadronTensor::MECHadronTensorType_t tensor_type)
 
double TensorContraction (int nu_pdg, int target_pdg, double Enu, double M_l, double T_l, double costh_l, int tensor_pdg, MECHadronTensor::MECHadronTensorType_t tensor_type)
 

Detailed Description

MEC utilities.

Author

Copyright (c) 2003-2019, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Function Documentation

bool genie::utils::mec::Getq0q3FromTlCostl ( double  Tl,
double  costl,
double  Enu,
double  ml,
double &  q0,
double &  q3 
)

Definition at line 115 of file MECUtils.cxx.

References std_candles::pl, and ana::Sqrt().

Referenced by genie::NievesSimoVacasMECPXSec2016::XSec().

117 {
118  q0 = -999;
119  q3 = -999;
120 
121  double El = Tl + ml;
122  q0 = Enu - El;
123  if(q0 < 0) {
124  q0 = -999;
125  q3 = -999;
126  return false;
127  }
128 
129  double pl = TMath::Sqrt(El * El - ml * ml);
130  double q3sq = pl * pl + Enu * Enu - 2.0 * pl * Enu * costl;
131  if(q3sq < 0) {
132  q0 = -999;
133  q3 = -999;
134  return false;
135  }
136  q3 = TMath::Sqrt(q3sq);
137 
138  return true;
139 }
ana::Sqrt
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
std_candles.pl
pl
Definition: std_candles.py:675
bool genie::utils::mec::GetTlCostlFromq0q3 ( double  q0,
double  q3,
double  Enu,
double  ml,
double &  Tl,
double &  costl 
)

Definition at line 76 of file MECUtils.cxx.

References std_candles::pl, and ana::Sqrt().

78 {
79  Tl = Enu - q0 - ml;
80  if(Tl < 0.) {
81  costl = -999;
82  Tl = -999;
83  return false;
84  }
85 
86  double El = Tl + ml;
87  double plsq = El * El - ml * ml;
88  if(plsq < 0.) {
89  costl = -999;
90  Tl = -999;
91  return false;
92  }
93  double pl = TMath::Sqrt(plsq);
94  double numerator = Enu * Enu + pl * pl - q3 * q3;
95  double denominator = 2.0 * pl * Enu;
96  if(denominator <= 0.0) {
97  costl = 0.0;
98  if(denominator < 0.0) {
99  return false;
100  }
101  }
102  else {
103  costl = numerator / denominator;
104  }
105 
106  if(TMath::Abs(numerator) > TMath::Abs(denominator)){
107  costl = -999;
108  Tl = -999;
109  return false;
110  }
111 
112  return true;
113 }
ana::Sqrt
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
std_candles.pl
pl
Definition: std_candles.py:675
double genie::utils::mec::GetTmuCostFromq0q3 ( double  dq0,
double  dq3,
double  Enu,
double  lmass,
double &  tmu,
double &  cost,
double &  area 
)

Definition at line 27 of file MECUtils.cxx.

References std::sqrt(), and ana::Sqrt().

30 {
31  tmu = Enu - dq0 - lmass;
32  if(tmu < 0.0){
33  cost = -999;
34  tmu = -999;
35  area = 0.0;
36  return -999;
37  }
38 
39  double thisE = tmu + lmass;
40  double thisp = sqrt( thisE * thisE - lmass * lmass);
41  double numerator = Enu * Enu + thisp * thisp - dq3 * dq3;
42  double denominator = 2.0 * thisp * Enu;
43  if(denominator <= 0.0 ){
44  cost = 0.0;
45  if(denominator < 0.0){
46  return -999;
47  }
48  }
49  else cost = numerator / denominator;
50 
51  if(TMath::Abs(numerator) > TMath::Abs(denominator)){
52  cost = -999;
53  tmu = -999;
54  area = 0.0;
55  return -999;
56  }
57 
58  // xCrossSect is not yet in the right units for this particular case.
59  // need areaElement to go dsigma/dTmudcost to dsigma/dq0dq3
60  // Recompute the area element jacobian
61 
62  // dT/dq0 x dC/dq3 - dT/dq3 x dC/dq0
63  double areaElement = 0.0;
64  //double veryCloseToZero = 0.000001; // in GeV, this should be safe.
65  double insqrt = 0.0;
66  numerator = 0.0;
67  insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - lmass * lmass;
68  numerator = dq3 / Enu;
69  if(insqrt < 0.0)areaElement=0.0;
70  else areaElement = numerator / TMath::Sqrt(insqrt);
71  area = areaElement;
72 
73  return 0;
74 }
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
ana::Sqrt
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double genie::utils::mec::J ( double  q0,
double  q3,
double  Enu,
double  ml 
)

Definition at line 141 of file MECUtils.cxx.

References ana::Sqrt().

Referenced by bpfit::Lutz::CalcAinv(), genie::SKKinematicsGenerator::CalculateKin_AtharSingleKaon(), genie::utils::gsl::d3Xsec_dTldTkdCosThetal::DoEval(), genie::utils::gsl::d2Xsec_dQ2dv::DoEval(), rb::Cluster::Exclude(), lem::FindMatchesAlg::FindMatchesHeads(), stan::math::jacobian(), genie::utils::kinematics::Jacobian(), stan::math::multiply_lower_tri_self_transpose(), genie::DFRKinematicsGenerator::ProcessEventRecord(), genie::QELKinematicsGenerator::ProcessEventRecord(), genie::RESKinematicsGenerator::ProcessEventRecord(), genie::DMELKinematicsGenerator::ProcessEventRecord(), genie::DISKinematicsGenerator::ProcessEventRecord(), genie::DMDISKinematicsGenerator::ProcessEventRecord(), genie::QELEventGeneratorSM::ProcessEventRecord(), genie::MECGenerator::SelectEmpiricalKinematics(), stan::math::student_t_cdf(), stan::math::student_t_lccdf(), stan::math::student_t_lcdf(), stan::math::unit_vector_constrain(), genie::COHElasticPXSec::XSec(), genie::EmpiricalMECPXSec2015::XSec(), genie::SlowRsclCharmDISPXSecLO::XSec(), genie::ReinDFRPXSec::XSec(), genie::AhrensNCELPXSec::XSec(), genie::StrumiaVissaniIBDPXSec::XSec(), genie::RosenbluthPXSec::XSec(), genie::QPMDISPXSec::XSec(), genie::AivazisCharmPXSecLO::XSec(), genie::IMDAnnihilationPXSec::XSec(), genie::P33PaschosLalakulichPXSec::XSec(), genie::AhrensDMELPXSec::XSec(), genie::QPMDMDISPXSec::XSec(), genie::LwlynSmithQELCCPXSec::XSec(), genie::BergerSehgalCOHPiPXSec2015::XSec(), genie::NuElectronPXSec::XSec(), genie::PaisQELLambdaPXSec::XSec(), genie::ReinSehgalCOHPiPXSec::XSec(), genie::KovalenkoQELCharmPXSec::XSec(), genie::BardinIMDRadCorPXSec::XSec(), genie::ReinSehgalRESPXSec::XSec(), genie::BSKLNBaseRESPXSec2014::XSec(), and genie::NievesQELCCPXSec::XSec().

143 {
144  // dT/dq0 x dC/dq3 - dT/dq3 x dC/dq0
145  double area = 0.0;
146  double insqrt = 0.0;
147  insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - ml * ml;
148  double numerator = dq3 / Enu;
149  if(insqrt < 0.0) {
150  area=0.0;
151  }
152  else {
153  area = numerator / TMath::Sqrt(insqrt);
154  }
155  return area;
156 }
ana::Sqrt
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double genie::utils::mec::Qvalue ( int  targetpdg,
int  nupdg 
)

Definition at line 158 of file MECUtils.cxx.

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

Referenced by TensorContraction().

159 {
160 // The had tensor calculation requires parameters that describe the
161 // nuclear density. For the Valencia model they are taken by
162 // C. Garcia-Recio, Nieves, Oset Nucl.Phys.A 547 (1992) 473-487 Table I
163 // and simplifies that the p and n density are not necessarily the same?
164 // Standard tables such as deVries et al are similar ~5%
165 // This is the only nuclear input on the hadron side of the Hadron Tensor.
166 // and is what makes PseudoPb and PseudoFe different than Ni56 and Rf208
167 //
168 
169  int nearpdg = targetpdg;
170 
171  if(nupdg > 0){
172  // get Z+1 for CC interaction e.g. 1000210400 from 1000200400
173  nearpdg = targetpdg + 10000;
174  } else {
175  // get Z-1 for CC interaction e.g. 1000210400 from 1000200400
176  nearpdg = targetpdg - 10000;
177  }
178 
179  TParticlePDG *partf = PDGLibrary::Instance()->Find(nearpdg);
180  TParticlePDG *parti = PDGLibrary::Instance()->Find(targetpdg);
181 
182  if(NULL == partf || NULL == parti){
183  // maybe not every valid nucleus in the table has Z+1 or Z-1
184  // for example, Ca40 did not.
185  LOG("MECUtils", pFATAL)
186  << "Can't get qvalue, nucleus " << targetpdg << " or " << nearpdg
187  << " is not in the table of nuclei in /data/evgen/pdg ";
188  exit(-1);
189  }
190 
191  double massf = partf->Mass();
192  double massi = parti->Mass();
193 
194  // keep this here as a reminder of what not to do
195  // +/- emass; // subtract electron from Z+1, add it to Z-1
196  // the lookup table gives the nucleus mass, not the atomic
197  // mass, so don't need this.
198 
199  return massf - massi; // in GeV.
200 }
pFATAL
#define pFATAL
Definition: Messenger.h:57
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
exit
exit(0)
double genie::utils::mec::TensorContraction ( const Interaction *  interaction,
int  tensor_pdg,
MECHadronTensor::MECHadronTensorType_t  tensor_type 
)

Definition at line 202 of file MECUtils.cxx.

References genie::Interaction::FSPrimLepton(), genie::Kinematics::GetKV(), genie::Interaction::InitState(), genie::Interaction::Kine(), genie::kKVctl, genie::kKVTl, genie::kRfLab, genie::Target::Pdg(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), and genie::InitialState::Tgt().

Referenced by genie::NievesSimoVacasMECPXSec2016::XSec().

206 {
207  int target_pdg = interaction->InitState().Tgt().Pdg();
208  int nu_pdg = interaction->InitState().ProbePdg();
209  double Enu = interaction->InitState().ProbeE(kRfLab);
210  double Ml = interaction->FSPrimLepton()->Mass();
211  double Tl = interaction->Kine().GetKV(kKVTl);
212  double costhl = interaction->Kine().GetKV(kKVctl);
213 
214  return genie::utils::mec::TensorContraction(
215  nu_pdg, target_pdg, Enu, Ml, Tl, costhl, tensor_pdg, tensor_type);
216 }
genie::Target::Pdg
int Pdg(void) const
Definition: Target.h:72
genie::Interaction::Kine
const Kinematics & Kine(void) const
Definition: Interaction.h:71
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:65
genie::Kinematics::GetKV
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:333
genie::Interaction::FSPrimLepton
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Definition: Interaction.cxx:124
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition: Interaction.h:69
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:67
genie::utils::mec::TensorContraction
double TensorContraction(const Interaction *interaction, int tensor_pdg, MECHadronTensor::MECHadronTensorType_t tensor_type)
Definition: MECUtils.cxx:202
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition: InitialState.cxx:365
double genie::utils::mec::TensorContraction ( int  nu_pdg,
int  target_pdg,
double  Enu,
double  M_l,
double  T_l,
double  costh_l,
int  tensor_pdg,
MECHadronTensor::MECHadronTensorType_t  tensor_type 
)

Definition at line 218 of file MECUtils.cxx.

References genie::units::cm2, E, MECModelEnuComparisons::i, genie::MECHadronTensor::Instance(), genie::constants::kGF2, genie::constants::kPi, LOG, pDEBUG, Qvalue(), std::sqrt(), ana::Sqrt(), genie::MECHadronTensor::TensorTable(), and xsec.

223 {
224  TLorentzVector v4lep;
225  TLorentzVector v4Nu(0,0,Enu,Enu); // assuming traveling along z:
226  TLorentzVector v4q;
227  double q0nucleus;
228  double facconv = 0.0389391289e15; // std::pow(0.19733,2)*1e15;
229 
230  double myQvalue = genie::utils::mec::Qvalue(targetpdg, nupdg);
231 
232  // Angles
233  double sinthl = 1. - costhl * costhl;
234  if(sinthl < 0.0) sinthl = 0.0;
235  else sinthl = TMath::Sqrt(sinthl);
236 
237  double Cosh = TMath::Cos(TMath::ACos(costhl)/2.);
238  double Sinh = TMath::Sin(TMath::ACos(costhl)/2.);
239 
240  // Lepton
241  v4lep.SetE( Tl + Ml );
242  // energy transfer from the lepton
243  double q0 = v4Nu.E() - v4lep.E();
244  // energy transfer that actually gets to the nucleons
245  q0nucleus = q0 - myQvalue;
246 
247  // Define some calculation placeholders
248  double part1, part2;
249  double modkprime ;
250  double w1, w2, w3, w4, w5;
251  double wtotd[5];
252 
253  double xsec = 0;
254  if (q0nucleus <= 0){
255  //nothing transfered to nucleus thus no 2p2h
256  xsec = 0.;
257  } else {
258  // energy was transfered to the nucleus
259  modkprime = std::sqrt(TMath::Power(v4lep.E(),2)-TMath::Power(Ml,2));
260  v4lep.SetX(modkprime*sinthl);
261  v4lep.SetY(0);
262  v4lep.SetZ(modkprime*costhl);
263 
264  //q: v4q = v4Nu - v4lep;
265  v4q.SetE(q0nucleus);
266  v4q.SetX(v4Nu.X() - v4lep.X());
267  v4q.SetY(v4Nu.Y() - v4lep.Y());
268  v4q.SetZ(v4Nu.Z() - v4lep.Z());
269 
270  MECHadronTensor * hadtensor = MECHadronTensor::Instance();
271  const vector <genie::BLI2DNonUnifGrid *> &
272  tensor_table = hadtensor->TensorTable(tensorpdg, tensor_type);
273 
274  for (int i=0 ; i < 5; i++){
275  wtotd[i] = tensor_table[i]->Evaluate(v4q.Vect().Mag(),v4q.E());
276  }
277 
278  // calculate hadron tensor components
279  // these are footnote 2 of Nieves PRC 70 055503
280  double W00 = wtotd[0];
281  double W0Z = wtotd[1];
282  double WXX = wtotd[2];
283  double WXY = wtotd[3];
284  double WZZ = wtotd[4];
285 
286  w1=WXX/2.;
287  w2=(W00+WXX+(q0*q0/(v4q.Vect().Mag()*v4q.Vect().Mag())
288  *(WZZ-WXX))
289  -(2.*q0/v4q.Vect().Mag()*W0Z))/2.;
290  w3=WXY/v4q.Vect().Mag();
291  w4=(WZZ-WXX)/(2.*v4q.Vect().Mag()*v4q.Vect().Mag());
292  w5=(W0Z-(q0/v4q.Vect().Mag()*(WZZ-WXX)))/v4q.Vect().Mag();
293  //w6 we have no need for w6, noted at the end of IIA.
294 
295  // adjust for anti neutrinos
296 
297  if (nupdg < 0) w3 = -1. * w3;
298 
299  // calculate cross section, in parts
300  double xw1 = w1*costhl;
301  double xw2 = w2/2.*costhl;
302  double xw3 = w3/2.*(v4lep.E()+modkprime-(v4lep.E()+v4Nu.E())*costhl);
303  double xw4 = w4/2.*(Ml*Ml*costhl+2.*v4lep.E()*(v4lep.E()+modkprime)*Sinh*Sinh);
304  double xw5 = w5*(v4lep.E()+modkprime)/2.;
305  part1 = xw1 - xw2 + xw3 + xw4 - xw5;
306 
307  double yw1 = 2.*w1*Sinh*Sinh;
308  double yw2 = w2*Cosh*Cosh;
309  double yw3 = w3*(v4lep.E()+v4Nu.E())*Sinh*Sinh;
310  double yw4 = Ml*Ml*part1/(v4lep.E()*(v4lep.E()+modkprime));
311  part2 = yw1 + yw2 - yw3 + yw4;
312 
313  xsec = modkprime*v4lep.E()*kGF2*2./kPi*part2*facconv;
314 
315  if( ! (xsec >= 0.0) ){
316  // Traps negative numbers and nan's.
317  // Sometimes costhl is just over 1.0 due to fluke or numerical
318  // precision, so sinthl would be undefined.
319  LOG("MECUtils", pDEBUG)
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 
327  xsec = 0;
328  }
329  }
330 
331 // LOG("MECUtils", pDEBUG)
332 // << "xsec(Enu = " << Enu << " GeV, Ml = " << Ml << " GeV; "
333 // << "Tl = " << Tl << " GeV, costhl = " << costhl << ") = "
334 // << xsec << " x 1E-41 cm^2";
335 
336  return (xsec * (1.0E-41 * units::cm2));
337 }
kPi
const double kPi
Definition: piN_asymmetric_decay.C:74
Munits::cm2
static constexpr Double_t cm2
Definition: Munits.h:141
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
genie::utils::mec::Qvalue
double Qvalue(int targetpdg, int nupdg)
Definition: MECUtils.cxx:158
genie::MECHadronTensor::TensorTable
const vector< genie::BLI2DNonUnifGrid * > & TensorTable(int targetpdg, MECHadronTensor::MECHadronTensorType_t type)
Definition: MECHadronTensor.cxx:125
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
E
Float_t E
Definition: plot.C:20
genie::MECHadronTensor
Singleton class to load and store MEC hadron tensor tables, to aid in the implementation (and improve...
Definition: MECHadronTensor.h:45
xsec
Double_t xsec[nknots]
Definition: testXsec.C:47
ana::Sqrt
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
genie::constants::kGF2
static const double kGF2
Definition: Constants.h:60
pDEBUG
#define pDEBUG
Definition: Messenger.h:64
Generated on Wed Feb 24 2021 14:19:14 for NOvA by   doxygen 1.8.11