SmithMonizUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Igor Kakorin <kakorin@jinr.ru>, Joint Institute for Nuclear Research
8  adapted from fortran code provided by
9  Konstantin Kuzmin <kkuzmin@theor.jinr.ru>,
10  Joint Institute for Nuclear Research, Institute for Theoretical and Experimental Physics
11  Vladimir Lyubushkin,
12  Joint Institute for Nuclear Research
13  Vadim Naumov <vnaumov@theor.jinr.ru>,
14  Joint Institute for Nuclear Research
15  based on code of Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
16  University of Liverpool & STFC Rutherford Appleton Lab
17 
18  For the class documentation see the corresponding header file.
19 
20 
21 */
22 //____________________________________________________________________________
23 
24 #include <sstream>
25 
26 #include <TMath.h>
27 #include <Math/IFunction.h>
28 #include <Math/RootFinder.h>
29 
39 #include "Framework/Utils/Range1.h"
41 
42 using namespace genie;
43 using namespace genie::constants;
44 using std::ostringstream;
45 
46 //____________________________________________________________________________
48 Algorithm("genie::SmithMonizUtils")
49 {
50 
51 }
52 //____________________________________________________________________________
54 Algorithm("genie::SmithMonizUtils", config)
55 {
56 
57 }
58 //____________________________________________________________________________
60 {
61 
62 }
63 //____________________________________________________________________________
65 {
66  Algorithm::Configure(config);
67  this->LoadConfig();
68 }
69 //____________________________________________________________________________
71 {
72  Algorithm::Configure(config);
73  this->LoadConfig();
74 }
75 //____________________________________________________________________________
77 {
78 
79 
80  GetParam( "FermiMomentumTable", fKFTable);
81  GetParam( "RFG-UseParametrization", fUseParametrization);
82 
83 
84  // Load removal energy for specific nuclei from either the algorithm's
85  // configuration file or the UserPhysicsOptions file.
86  // If none is used use Wapstra's semi-empirical formula.
87  const std::string keyStart = "RFG-NucRemovalE@Pdg=";
88 
89  RgIMap entries = GetConfig().GetItemMap();
90 
91  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it){
92  const std::string& key = it->first;
93  int pdg = 0;
94  int Z = 0;
95  if (0 == key.compare(0, keyStart.size(), keyStart.c_str())) {
96  pdg = atoi(key.c_str() + keyStart.size());
97  Z = pdg::IonPdgCodeToZ(pdg);
98  }
99  if (0 != pdg && 0 != Z) {
100  ostringstream key_ss ;
101  key_ss << keyStart << pdg;
102  RgKey rgkey = key_ss.str();
103  double eb ;
104  GetParam( rgkey, eb ) ;
105  eb = TMath::Max(eb, 0.);
106  fNucRmvE.insert(map<int,double>::value_type(Z,eb));
107  }
108  }
109 
110 }
111 //____________________________________________________________________________
112 // Set the variables necessary for further calculations
114 {
115  fInteraction = interaction;
116  // Get kinematics & init-state parameters
117  // unused // const Kinematics & kinematics = interaction -> Kine();
118  const InitialState & init_state = interaction -> InitState();
119  const Target & target = init_state.Tgt();
120  PDGLibrary * pdglib = PDGLibrary::Instance();
121 
122  E_nu = interaction->InitState().ProbeE(kRfLab); // Neutrino energy (GeV)
123 
124  assert(target.HitNucIsSet());
125  // get lepton&nuclear masses (init & final state nucleus)
126  m_lep = interaction->FSPrimLepton()->Mass(); // Mass of final charged lepton (GeV)
127  mm_lep = TMath::Power(m_lep, 2);
128  int nucl_pdg_ini = target.HitNucPdg();
129  m_ini = target.HitNucMass();
130  mm_ini = TMath::Power(m_ini, 2);
131  int nucl_pdg_fin = genie::pdg::SwitchProtonNeutron(nucl_pdg_ini);
132  TParticlePDG * nucl_fin = pdglib->Find( nucl_pdg_fin );
133  m_fin = nucl_fin -> Mass(); // Mass of final hadron or hadron system (GeV)
134  mm_fin = TMath::Power(m_fin, 2);
135  m_tar = target.Mass(); // Mass of target nucleus (GeV)
136  mm_tar = TMath::Power(m_tar, 2);
137 
138  // For hydrogen and deuterium RFG is not applied
139  if (target.A()<3)
140  {
141  E_BIN = P_Fermi = m_rnu = mm_rnu = 0;
142  return;
143  }
144 
145  bool is_p = pdg::IsProton(nucl_pdg_ini);
146  int Zi = target.Z();
147  int Ai = target.A();
148  int Zf = (is_p) ? Zi-1 : Zi;
149  int Af = Ai-1;
150  TParticlePDG * nucl_f = pdglib->Find( pdg::IonPdgCode(Af, Zf) );
151  if(!nucl_f) {
152  LOG("SmithMoniz", pFATAL)
153  << "Unknwown nuclear target! No target with code: "
154  << pdg::IonPdgCode(Af, Zf) << " in PDGLibrary!";
155  exit(1);
156  }
157 
158  m_rnu = nucl_f -> Mass(); // Mass of residual nucleus (GeV)
159  mm_rnu = TMath::Power(m_rnu, 2);
160 
161  int Z = target.Z();
162  int A = target.A();
163  int N = A-Z;
164 
165  // Maximum value of Fermi momentum of target nucleon (GeV)
166  if (A<6 || !fUseParametrization)
167  {
168  //-- look up the Fermi momentum for this Target
170  const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
171  P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucl_pdg_ini);
172  }
173  else
174  {
175  //-- define the Fermi momentum for this Target
176  //
178  //-- correct the Fermi momentum for the struck nucleon
179  if(is_p) P_Fermi *= TMath::Power( 2.*Z/A, 1./3);
180  else
181  P_Fermi *= TMath::Power( 2.*N/A, 1./3);
182  }
183 
184  // Neutrino binding energy (GeV)
185  if (target.A()<6 || !fUseParametrization)
186  {
187  map<int,double>::const_iterator it = fNucRmvE.find(Z);
188  if(it != fNucRmvE.end()) E_BIN = it->second;
190  }
191  else
193 
194 
195 }
196 //____________________________________________________________________________
197 // Get the neutrino energy threshold
199 {
200 
202 
203 
204  const int MFC = 10000; // Maximum of function call
205  const double EPSABS = 0;
206  const double EPSREL = 1.0e-08;
207  double E_min= ((m_lep+m_rnu+m_fin)*(m_lep+m_rnu+m_fin)-mm_tar)/(2*m_tar);
208  double Enu_2= 5.0e+00;
209  double Enu_rf;
210  if (QEL_EnuMin_SM(E_min)>0)
211  {
212  // C++ analog of fortran function Enu_rf= DZEROX(E_min,Enu_2,EPS,MFC,QEL_EnuMin_SM,1)
213  ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
214  rfgb.Solve(QEL_EnuMin_SM_, E_min, Enu_2, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
215  Enu_rf = rfgb.Root();
216 
217  }
218  else
219  {
220  Enu_rf=-1.0e+01;
221  }
222  E_min = TMath::Max(E_min,Enu_rf);
223  if(E_min<0)
224  {
225  E_min = 0;
226  LOG("SmithMoniz", pDEBUG) << "E_min = " << E_min;
227  }
228  return E_min;
229 
230 }
231 //____________________________________________________________________________
232 // The auxiliary function for determining energy threshold
233 double SmithMonizUtils::QEL_EnuMin_SM(double Enu) const
234 {
235  const double EPS = 1.0e-06;
236  const double Delta= 1.0e-14;
237  const double Precision = std::numeric_limits<double>::epsilon();
238  Enu_in = Enu;
239  double s = 2*Enu*m_tar+mm_tar;
240  double W2 = (m_rnu+m_fin)*(m_rnu+m_fin);
241  double c = 0.5*(W2+mm_lep-mm_tar*(W2-mm_lep)/s);
242  double sqrtD = TMath::Sqrt(TMath::Max(Precision,LambdaFUNCTION(1.0, mm_lep/s, W2/s)));
243  double tmp = 0.5*(s-mm_tar);
244  double Q2_lim1= tmp*(1.0-sqrtD)-c;
245  double Q2_lim2= tmp*(1.0+sqrtD)-c;
246  // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_lim1,Q2_lim2,EPS,Delta,Q2_0,F_MIN,LLM)
248  double Q2_0,F_MIN;
249  bool LLM;
250  DMINFC(Q2lim1_SM_,Q2_lim1,Q2_lim2,EPS,Delta,Q2_0,F_MIN,LLM);
251  return F_MIN;
252 }
253 //____________________________________________________________________________
254 // The auxiliary function for determining Q2-range
255 double SmithMonizUtils::Q2lim1_SM(double Q2) const
256 {
257  double s = 2*Enu_in*m_tar+mm_tar;
258  double nu_max = (s-mm_tar-mm_lep*(s-mm_tar)/(Q2+mm_lep)-mm_tar*(Q2+mm_lep)/(s-mm_tar))/(2*m_tar);
259  double E = sqrt(P_Fermi*P_Fermi+mm_ini);
260  double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
261  double a = 0.5*(Q2+mm_fin-b);
262  double nu_1 = (a*(E-E_BIN)-P_Fermi*TMath::Sqrt(a*a+Q2*b))/b;
263  return nu_1-nu_max;
264 
265 }
266 //____________________________________________________________________________
267 // The auxiliary function for determining Q2-range
268 double SmithMonizUtils::Q2lim2_SM(double Q2) const
269 {
270  double nu_min = ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/(2*m_tar);
271  double E = sqrt(P_Fermi*P_Fermi+mm_ini);
272  double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
273  double a = (Q2+mm_fin-b)*0.5;
274  double nu_2 = (a*(E-E_BIN)+P_Fermi*TMath::Sqrt(a*a+Q2*b))/b;
275  return nu_min-nu_2;
276 
277 }
278 //____________________________________________________________________________
279 // Return allowed Q2-range
281 {
282 
283 
284  const int MFC = 1000; // Maximum of function call
285  const double EPS = 1.0e-08;
286  const double Delta= 1.0e-14;
287  const double EPSABS = 0;
288  const double EPSREL = 1.0e-08;
289 
290  Enu_in = E_nu;
291  double s = 2*E_nu*m_tar+mm_tar;
292  double W2 = (m_rnu+m_fin)*(m_rnu+m_fin);
293  double c = 0.5*(W2+mm_lep-mm_tar*(W2-mm_lep)/s);
294  double sqrtD = TMath::Sqrt(LambdaFUNCTION(1.0, mm_lep/s, W2/s));
295  double tmp = 0.5*(s-mm_tar);
296  double Q2_min = tmp*(1.0-sqrtD)-c;
297  double Q2_max = tmp*(1.0+sqrtD)-c;
298  // if the nucleus is hydrogen or deuterium then there is no need in further calculation
299  if (E_BIN == 0 && P_Fermi == 0)
300  {
301  Q2_min= TMath::Max(Q2_min,0.0);
302  Range1D_t R(Q2_min,Q2_max);
303  return R;
304  }
305  double F_MIN, Q2_0;
306  bool LLM;
307  // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM)
309  DMINFC(Q2lim1_SM_,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM);
310  if (F_MIN>0)
311  {
312  LOG("SmithMoniz", pFATAL)
313  << "No overlapped area for energy " << E_nu << "\n" <<
314  "Q2_min=" << Q2_min << " Q2_max=" << Q2_max << "\n" <<
315  "Q2_0=" << Q2_0 << " F_MIN=" << F_MIN;
316  exit(1);
317  }
318  if (Q2lim1_SM(Q2_min)>0)
319  {
320  //C++ analog of fortran function Q2_RF=DZEROX(Q2_min,Q2_0,EPS,MFC,Q2lim1_SM,1)
321  ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
322  rfgb.Solve(Q2lim1_SM_, Q2_min, Q2_0, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
323  double Q2_RF = rfgb.Root();
324  Q2_min= TMath::Max(Q2_min,Q2_RF);
325  }
326  if(Q2lim1_SM(Q2_max)>0)
327  {
328  // C++ analog of fortran function Q2_RF=DZEROX(Q2_0,Q2_max,Eps,MFC,Q2lim1_SM,1)
329  ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
330  rfgb.Solve(Q2lim1_SM_, Q2_0, Q2_max, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
331  double Q2_RF = rfgb.Root();
332  Q2_max= TMath::Min(Q2_max,Q2_RF);
333  }
335  if (Q2lim2_SM(Q2_min)>0)
336  {
337  if(Q2lim2_SM(Q2_max)>0)
338  {
339  LOG("SmithMoniz", pFATAL) << "No overlapped area for energy " << E_nu;
340  exit(1);
341  }
342  else
343  {
344  // C++ analog of fortran function Q2_RF = DZEROX(Q2_min,Q2_max,Eps,MFC,Q2lim2_SM,1)
345  ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
346  rfgb.Solve(Q2lim2_SM_, Q2_min,Q2_max, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
347  double Q2_RF = rfgb.Root();
348  Q2_min= TMath::Max(Q2_min,Q2_RF);
349  }
350  }
351 
352  Q2_min= TMath::Max(Q2_min,0.0);
353  Range1D_t R(Q2_min,Q2_max);
354  return R;
355 
356 }
357 //____________________________________________________________________________
358 // Return allowed v-range for given Q2
360 {
361  double s = 2*E_nu*m_tar+mm_tar;
362  double nu_min= ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/(2*m_tar);
363  double nu_max= (s-mm_tar-mm_lep*(s-mm_tar)/(Q2+mm_lep)-mm_tar*(Q2+mm_lep)/(s-mm_tar))/(2*m_tar);
364  double E = TMath::Sqrt(P_Fermi*P_Fermi+mm_ini);
365  double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
366  double a = (Q2+mm_fin-b)*0.5;
367  double tmp1 = a*(E-E_BIN);
368  double tmp2 = P_Fermi*TMath::Sqrt(a*a+Q2*b);
369  double nu_1 = (tmp1-tmp2)/b;
370  double nu_2 = (tmp1+tmp2)/b;
371  nu_min= TMath::Max(nu_min,nu_1);
372  nu_max= TMath::Min(nu_max,nu_2);
373  Range1D_t R;
374  if (E_BIN == 0 && P_Fermi == 0)
375  nu_max=nu_min;
376  if (nu_min<=nu_max)
377  R = Range1D_t(nu_min,nu_max);
378  else
379  R = Range1D_t(0.5*(nu_min+nu_max),0.5*(nu_min+nu_max));
380  return R;
381 
382 }
383 //____________________________________________________________________________
384 // Return allowed Fermi momentum range for given Q2 and v
386 {
387  double qv = TMath::Sqrt(nu*nu+Q2);
388  double c_f = (nu-E_BIN)/qv;
389  double d_f = (E_BIN*E_BIN-2*nu*E_BIN-Q2+mm_ini-mm_fin)/(2*qv*m_ini);
390  double Ef_min= m_ini*(c_f*d_f+TMath::Sqrt(1.0-c_f*c_f+d_f*d_f))/(1.0-c_f*c_f);
391  double kF_min= TMath::Sqrt(TMath::Max(Ef_min*Ef_min-mm_ini,0.0));
392  double kF_max= P_Fermi;
393  Range1D_t R;
394  if (kF_min<=kF_max)
395  R = Range1D_t(kF_min,kF_max);
396  else
397  R = Range1D_t(0.5*(kF_min+kF_max),0.5*(kF_min+kF_max));
398  return R;
399 
400 }
401 //____________________________________________________________________________
402 double SmithMonizUtils::LambdaFUNCTION(double a, double b, double c) const
403 {
404  return a*a+b*b+c*c-2*(a*b+b*c+a*c);
405 }
406 //____________________________________________________________________________
407 // C++ implementation of DMINC function from CERN library
408 void SmithMonizUtils::DMINFC(Func1D<SmithMonizUtils> &F, double A,double B, double EPS, double DELTA, double &X, double &Y, bool &LLM) const
409 {
410  const double W5 = TMath::Sqrt(5);
411  const double HV = (3-W5)/2;
412  const double HW = (W5-1)/2;
413  const double R1 = 1.0;
414  const double HF = R1/2;
415 
416  int N = -1;
417  if(A!=B) N = TMath::Nint(2.08*TMath::Log(TMath::Abs((A-B)/EPS)));
418  double C = A;
419  double D = B;
420  if(A>B)
421  {
422  C = B;
423  D = A;
424  }
425  bool LLT = true;
426  bool LGE = true;
427  double V, FV, W, FW, H;
428  while(true)
429  {
430  H = D-C;
431  if(N<0)
432  {
433  X = HF*(C+D);
434  Y = F(X);
435  LLM = TMath::Abs(X-A)>DELTA && TMath::Abs(X-B)>DELTA;
436  return;
437  }
438  if(LLT)
439  {
440  V = C+HV*H;
441  FV = F(V);
442  }
443  if(LGE)
444  {
445  W = C+HW*H;
446  FW = F(W);
447  }
448  if(FV<FW)
449  {
450  LLT = true;
451  LGE = false;
452  D = W;
453  W = V;
454  FW = FV;
455  }
456  else
457  {
458  LLT = false;
459  LGE = true;
460  C = V;
461  V = W;
462  FV = FW;
463  }
464  N = N-1;
465  }
466 }
467 //____________________________________________________________________________
468 // Density of Fermi gas, for case T_Fermi=0 is a step functio, which is blurred at T_Fermi!=0
469 double SmithMonizUtils::rho(double P_Fermi, double T_Fermi, double p)
470 {
471 
472  if (T_Fermi==0) //Pure Fermi gaz with T_Fermi=0
473  if(p<=P_Fermi) return 1.0;
474  else return 0.0;
475  else //Fermi-Dirac distribution
476  return 1.0/(1.0 + TMath::Exp(-(P_Fermi-p)/T_Fermi));
477 
478 
479 }
480 //____________________________________________________________________________
482 {
483  return E_BIN;
484 }
485 //____________________________________________________________________________
487 {
488  return P_Fermi;
489 }
490 //____________________________________________________________________________
491 double SmithMonizUtils::GetTheta_k(double v, double qv) const
492 {
493  return TMath::ACos((v + (mm_lep-v*v+qv*qv)/2/E_nu)/qv);
494 }
495 //____________________________________________________________________________
496 double SmithMonizUtils::GetTheta_p(double pv, double v, double qv, double &E_p) const
497 {
498  E_p = TMath::Sqrt(mm_ini+pv*pv)-E_BIN;
499  if (pv != 0)
500  return TMath::ACos(((v-E_BIN)*(2*E_p+v+E_BIN)-qv*qv+mm_ini-mm_fin)/(2*pv*qv));
501  else
502  return 0;
503 }
504 //____________________________________________________________________________
506 {
507  return vQES_SM_lim(Q2).min;
508 }
509 //____________________________________________________________________________
511 {
512  return -1.0*vQES_SM_lim(Q2).max;
513 }
514 //____________________________________________________________________________
516 {
517  double vol = 0.0;
518  if (ps == kPSQ2fE)
519  {
520  Range1D_t rQ2 = Q2QES_SM_lim();
521  vol = rQ2.max - rQ2.min;
522  return vol;
523  }
524  else if (ps == kPSQ2vfE)
525  {
526  Range1D_t rQ2 = Q2QES_SM_lim();
527  const int kNQ2 = 101;
528  double dQ2 = (rQ2.max-rQ2.min)/(kNQ2-1);
529  for(int iq2=0; iq2<kNQ2; iq2++) {
530  double Q2 = (rQ2.min + iq2*dQ2);
531  Range1D_t rvlims = vQES_SM_lim(Q2);
532  double dv = (rvlims.max-rvlims.min);
533  vol += (dQ2*dv);
534  }
535  return vol;
536  }
537  else
538  return 0;
539 }
void SetInteraction(const Interaction *i)
const float DELTA
Range1D_t Q2QES_SM_lim(void) const
double mm_lep
Squared mass of final charged lepton (GeV)
#define F(x, y, z)
Basic constants.
const XML_Char * target
Definition: expat.h:268
double m_ini
Mass of initial hadron or hadron system (GeV)
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
double E_BIN
Binding energy (GeV)
int HitNucPdg(void) const
Definition: Target.cxx:321
double E_nu_thr_SM(void) const
int A(void) const
Definition: Target.h:71
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double HitNucMass(void) const
Definition: Target.cxx:250
double m_lep
Mass of final charged lepton (GeV)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
#define pFATAL
Definition: Messenger.h:57
double mm_fin
Squared mass of final hadron or hadron system (GeV)
double vQES_SM_low_bound(double Q2) const
static FermiMomentumTablePool * Instance(void)
int SwitchProtonNeutron(int pdgc)
Definition: PDGUtils.cxx:319
double E_nu
Neutrino energy (GeV)
map< int, double > fNucRmvE
Algorithm abstract base class.
Definition: Algorithm.h:54
double mm_ini
Sqared mass of initial hadron or hadron system (GeV)
double Mass(Resonance_t res)
resonance mass (GeV)
A table of Fermi momentum constants.
Definition: config.py:1
double GetTheta_k(double v, double qv) const
Float_t tmp
Definition: plot.C:36
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Mass(void) const
Definition: Target.cxx:241
Range1D_t vQES_SM_lim(double Q2) const
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
double GetBindingEnergy(void) const
Float_t Y
Definition: plot.C:38
Float_t Z
Definition: plot.C:38
double Q2lim1_SM(double Q2) const
double mm_rnu
Squared mass of residual nucleus (GeV)
const double C
const XML_Char * s
Definition: expat.h:262
Summary information for an interaction.
Definition: Interaction.h:56
const Interaction * fInteraction
Range1D_t kFQES_SM_lim(double nu, double Q2) const
double BindEnergyPerNucleon(const Target &target)
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const RgIMap & GetItemMap(void) const
Definition: Registry.h:162
Float_t E
Definition: plot.C:20
const FermiMomentumTable * GetTable(string name)
double Q2lim2_SM(double Q2) const
const double a
double P_Fermi
Maximum value of Fermi momentum of target nucleon (GeV)
double BindEnergyPerNucleonParametrization(const Target &target)
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
double m_rnu
Mass of residual nucleus (GeV)
#define R(x)
int Z(void) const
Definition: Target.h:69
double m_tar
Mass of target nucleus (GeV)
double vQES_SM_upper_bound(double Q2) const
double QEL_EnuMin_SM(double E_nu) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
double Enu_in
Running neutrino energy (GeV)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
static const double A
Definition: Units.h:82
double max
Definition: Range1.h:54
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static double rho(double P_Fermi, double T_Fermi, double p)
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
bool HitNucIsSet(void) const
Definition: Target.cxx:300
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double PhaseSpaceVolume(KinePhaseSpace_t ps) const
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
double epsilon
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
exit(0)
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:53
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:53
double GetTheta_p(double pv, double v, double qv, double &E_p) const
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
double mm_tar
Squared mass of target nucleus (GeV)
void Configure(const Registry &config)
Float_t X
Definition: plot.C:38
#define W(x)
double ProbeE(RefFrame_t rf) const
double LambdaFUNCTION(double a, double b, double c) const
void DMINFC(Func1D< SmithMonizUtils > &F, double A, double B, double EPS, double DELTA, double &X, double &Y, bool &LLM) const
double m_fin
Mass of final hadron or hadron system (GeV)
double GetFermiMomentum(void) const
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:46