ReinSehgalRESPXSec.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: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Oct 05, 2009 - CA
14  Modified code to handle charged lepton scattering too.
15  Also, the helicity amplitude code now returns a `const RSHelicityAmpl &'.
16  @ July 23, 2010 - CA
17  BaryonResParams, and BreitWignerI, BaryonResDataSetI implementations are
18  now redundant. Get resonance parameters from BaryonResUtils and use the
19  Breit-Weigner functions from utils::bwfunc.
20  @ May 01, 2012 - CA
21  Pick nutau/nutaubar scaling factors from new location.
22  @ May 01, 2016 - Libo Jiang
23  Add W dependence to Delta->N gamma
24 
25 */
26 //____________________________________________________________________________
27 
28 #include <TMath.h>
29 #include <TSystem.h>
30 
45 #include "Framework/Utils/Range1.h"
46 #include "Framework/Utils/BWFunc.h"
54 
55 using namespace genie;
56 using namespace genie::constants;
57 
58 //____________________________________________________________________________
60 XSecAlgorithmI("genie::ReinSehgalRESPXSec")
61 {
62  fNuTauRdSpl = 0;
63  fNuTauBarRdSpl = 0;
64 }
65 //____________________________________________________________________________
67 XSecAlgorithmI("genie::ReinSehgalRESPXSec", config)
68 {
69  fNuTauRdSpl = 0;
70  fNuTauBarRdSpl = 0;
71 }
72 //____________________________________________________________________________
74 {
75  if(fNuTauRdSpl) delete fNuTauRdSpl;
76  if(fNuTauBarRdSpl) delete fNuTauBarRdSpl;
77 }
78 //____________________________________________________________________________
80  const Interaction * interaction, KinePhaseSpace_t kps) const
81 {
82  if(! this -> ValidProcess (interaction) ) return 0.;
83  if(! this -> ValidKinematics (interaction) ) return 0.;
84 
85  const InitialState & init_state = interaction -> InitState();
86  const ProcessInfo & proc_info = interaction -> ProcInfo();
87  const Target & target = init_state.Tgt();
88 
89  // Get kinematical parameters
90  const Kinematics & kinematics = interaction -> Kine();
91  double W = kinematics.W();
92  double q2 = kinematics.q2();
93 
94  // Under the DIS/RES joining scheme, xsec(RES)=0 for W>=Wcut
95  if(fUsingDisResJoin) {
96  if(W>=fWcut) {
97 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
98  LOG("ReinSehgalRes", pDEBUG)
99  << "RES/DIS Join Scheme: XSec[RES, W=" << W
100  << " >= Wcut=" << fWcut << "] = 0";
101 #endif
102  return 0;
103  }
104  }
105 
106  // Get the input baryon resonance
107  Resonance_t resonance = interaction->ExclTag().Resonance();
108  string resname = utils::res::AsString(resonance);
109  bool is_delta = utils::res::IsDelta (resonance);
110 
111  // Get the neutrino, hit nucleon & weak current
112  int nucpdgc = target.HitNucPdg();
113  int probepdgc = init_state.ProbePdg();
114  bool is_nu = pdg::IsNeutrino (probepdgc);
115  bool is_nubar = pdg::IsAntiNeutrino (probepdgc);
116  bool is_lplus = pdg::IsPosChargedLepton (probepdgc);
117  bool is_lminus = pdg::IsNegChargedLepton (probepdgc);
118  bool is_p = pdg::IsProton (nucpdgc);
119  bool is_n = pdg::IsNeutron (nucpdgc);
120  bool is_CC = proc_info.IsWeakCC();
121  bool is_NC = proc_info.IsWeakNC();
122  bool is_EM = proc_info.IsEM();
123 
124  if(is_CC && !is_delta) {
125  if((is_nu && is_p) || (is_nubar && is_n)) return 0;
126  }
127 
128  // Get baryon resonance parameters
129  int IR = utils::res::ResonanceIndex (resonance);
130  int LR = utils::res::OrbitalAngularMom (resonance);
131  double MR = utils::res::Mass (resonance);
132  double WR = utils::res::Width (resonance);
134 
135  // Following NeuGEN, avoid problems with underlying unphysical
136  // model assumptions by restricting the allowed W phase space
137  // around the resonance peak
138  if (fNormBW) {
139  if (W > MR + fN0ResMaxNWidths * WR && IR==0) return 0.;
140  else if (W > MR + fN2ResMaxNWidths * WR && IR==2) return 0.;
141  else if (W > MR + fGnResMaxNWidths * WR) return 0.;
142  }
143 
144  // Compute auxiliary & kinematical factors
145  double E = init_state.ProbeE(kRfHitNucRest);
146  double Mnuc = target.HitNucMass();
147  double W2 = TMath::Power(W, 2);
148  double Mnuc2 = TMath::Power(Mnuc, 2);
149  double k = 0.5 * (W2 - Mnuc2)/Mnuc;
150  double v = k - 0.5 * q2/Mnuc;
151  double v2 = TMath::Power(v, 2);
152  double Q2 = v2 - q2;
153  double Q = TMath::Sqrt(Q2);
154  double Eprime = E - v;
155  double U = 0.5 * (E + Eprime + Q) / E;
156  double V = 0.5 * (E + Eprime - Q) / E;
157  double U2 = TMath::Power(U, 2);
158  double V2 = TMath::Power(V, 2);
159  double UV = U*V;
160 
161 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
162  LOG("ReinSehgalRes", pDEBUG)
163  << "Kinematical params V = " << V << ", U = " << U;
164 #endif
165 
166  // Calculate the Feynman-Kislinger-Ravndall parameters
167 
168  double Go = TMath::Power(1 - 0.25 * q2/Mnuc2, 0.5-IR);
169  double GV = Go * TMath::Power( 1./(1-q2/fMv2), 2);
170  double GA = Go * TMath::Power( 1./(1-q2/fMa2), 2);
171 
172  if(is_EM) {
173  GA = 0.; // zero the axial term for EM scattering
174  }
175 
176  double d = TMath::Power(W+Mnuc,2.) - q2;
177  double sq2omg = TMath::Sqrt(2./fOmega);
178  double nomg = IR * fOmega;
179  double mq_w = Mnuc*Q/W;
180 
181  fFKR.Lamda = sq2omg * mq_w;
182  fFKR.Tv = GV / (3.*W*sq2omg);
183  fFKR.Rv = kSqrt2 * mq_w*(W+Mnuc)*GV / d;
184  fFKR.S = (-q2/Q2) * (3*W*Mnuc + q2 - Mnuc2) * GV / (6*Mnuc2);
185  fFKR.Ta = (2./3.) * (fZeta/sq2omg) * mq_w * GA / d;
186  fFKR.Ra = (kSqrt2/6.) * fZeta * (GA/W) * (W+Mnuc + 2*nomg*W/d );
187  fFKR.B = fZeta/(3.*W*sq2omg) * (1 + (W2-Mnuc2+q2)/ d) * GA;
188  fFKR.C = fZeta/(6.*Q) * (W2 - Mnuc2 + nomg*(W2-Mnuc2+q2)/d) * (GA/Mnuc);
189  fFKR.R = fFKR.Rv;
190  fFKR.Rplus = - (fFKR.Rv + fFKR.Ra);
191  fFKR.Rminus = - (fFKR.Rv - fFKR.Ra);
192  fFKR.T = fFKR.Tv;
193  fFKR.Tplus = - (fFKR.Tv + fFKR.Ta);
194  fFKR.Tminus = - (fFKR.Tv - fFKR.Ta);
195 
196 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
197  LOG("FKR", pDEBUG)
198  << "FKR params for RES = " << resname << " : " << fFKR;
199 #endif
200 
201  // Calculate the Rein-Sehgal Helicity Amplitudes
202 
203  const RSHelicityAmplModelI * hamplmod = 0;
204  if(is_CC) {
205  hamplmod = fHAmplModelCC;
206  }
207  else
208  if(is_NC) {
209  if (is_p) { hamplmod = fHAmplModelNCp;}
210  else { hamplmod = fHAmplModelNCn;}
211  }
212  else
213  if(is_EM) {
214  if (is_p) { hamplmod = fHAmplModelEMp;}
215  else { hamplmod = fHAmplModelEMn;}
216  }
217  assert(hamplmod);
218 
219  const RSHelicityAmpl & hampl = hamplmod->Compute(resonance, fFKR);
220 
221 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
222  LOG("RSHAmpl", pDEBUG)
223  << "Helicity Amplitudes for RES = " << resname << " : " << hampl;
224 #endif
225 
226  double g2 = kGF2;
227  if(is_CC) g2 = kGF2*fVud2;
228  // For EM interaction replace G_{Fermi} with :
229  // a_{em} * pi / ( sqrt(2) * sin^2(theta_weinberg) * Mass_{W}^2 }
230  // See C.Quigg, Gauge Theories of the Strong, Weak and E/M Interactions,
231  // ISBN 0-8053-6021-2, p.112 (6.3.57)
232  // Also, take int account that the photon propagator is 1/p^2 but the
233  // W propagator is 1/(p^2-Mass_{W}^2), so weight the EM case with
234  // Mass_{W}^4 / q^4
235  // So, overall:
236  // G_{Fermi}^2 --> a_{em}^2 * pi^2 / (2 * sin^4(theta_weinberg) * q^{4})
237  //
238  if(is_EM) {
239  double q4 = q2*q2;
240  g2 = kAem2 * kPi2 / (2.0 * fSin48w * q4);
241  }
242 
243  // Compute the cross section
244 
245  double sig0 = 0.125*(g2/kPi)*(-q2/Q2)*(W/Mnuc);
246  double scLR = W/Mnuc;
247  double scS = (Mnuc/W)*(-Q2/q2);
248  double sigL = scLR* (hampl.Amp2Plus3 () + hampl.Amp2Plus1 ());
249  double sigR = scLR* (hampl.Amp2Minus3() + hampl.Amp2Minus1());
250  double sigS = scS * (hampl.Amp20Plus () + hampl.Amp20Minus());
251 
252 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
253  LOG("ReinSehgalRes", pDEBUG) << "sig_{0} = " << sig0;
254  LOG("ReinSehgalRes", pDEBUG) << "sig_{L} = " << sigL;
255  LOG("ReinSehgalRes", pDEBUG) << "sig_{R} = " << sigR;
256  LOG("ReinSehgalRes", pDEBUG) << "sig_{S} = " << sigS;
257 #endif
258 
259  double xsec = 0.0;
260  if (is_nu || is_lminus) {
261  xsec = sig0*(V2*sigR + U2*sigL + 2*UV*sigS);
262  }
263  else
264  if (is_nubar || is_lplus) {
265  xsec = sig0*(U2*sigR + V2*sigL + 2*UV*sigS);
266  }
267  xsec = TMath::Max(0.,xsec);
268 
269  double mult = 1.0;
270  if(is_CC && is_delta) {
271  if((is_nu && is_p) || (is_nubar && is_n)) mult=3.0;
272  }
273  xsec *= mult;
274 
275  // Check whether the cross section is to be weighted with a
276  // Breit-Wigner distribution (default: true)
277  double bw = 1.0;
278  if(fWghtBW) {
279  //different Delta photon decay branch
280  if(is_delta){
281  bw = utils::bwfunc::BreitWignerLGamma(W,LR,MR,WR,NR);
282  }
283  else{
284  bw = utils::bwfunc::BreitWignerL(W,LR,MR,WR,NR);
285  }
286  }
287 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
288  LOG("ReinSehgalRes", pDEBUG)
289  << "BreitWigner(RES=" << resname << ", W=" << W << ") = " << bw;
290 #endif
291  xsec *= bw;
292 
293  // Apply NeuGEN nutau cross section reduction factors
294  double rf = 1.0;
295  Spline * spl = 0;
296  if (is_CC && fUsingNuTauScaling) {
297  if (pdg::IsNuTau (probepdgc)) spl = fNuTauRdSpl;
298  else if (pdg::IsAntiNuTau(probepdgc)) spl = fNuTauBarRdSpl;
299 
300  if(spl) {
301  if(E <spl->XMax()) rf = spl->Evaluate(E);
302  }
303  }
304  xsec *= rf;
305 
306  // Apply given scaling factor
307  double xsec_scale = 1.;
308  if (is_CC) { xsec_scale = fXSecScaleCC; }
309  else if (is_NC) { xsec_scale = fXSecScaleNC; }
310  xsec *= xsec_scale;
311 
312 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
313  LOG("ReinSehgalRes", pINFO)
314  << "\n d2xsec/dQ2dW" << "[" << interaction->AsString()
315  << "](W=" << W << ", q2=" << q2 << ", E=" << E << ") = " << xsec;
316 #endif
317 
318  // The algorithm computes d^2xsec/dWdQ2
319  // Check whether variable tranformation is needed
320  if(kps!=kPSWQ2fE) {
321  double J = utils::kinematics::Jacobian(interaction,kPSWQ2fE,kps);
322  xsec *= J;
323  }
324 
325  // If requested return the free nucleon xsec even for input nuclear tgt
326  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
327 
328 
329  int Z = target.Z();
330  int A = target.A();
331  int N = A-Z;
332 
333  // Take into account the number of scattering centers in the target
334  int NNucl = (is_p) ? Z : N;
335 
336  xsec*=NNucl; // nuclear xsec (no nuclear suppression factor)
337 
338  if (fUsePauliBlocking && A!=1)
339  {
340  // Calculation of Pauli blocking according references:
341  //
342  // [1] S.L. Adler, S. Nussinov, and E.A. Paschos, "Nuclear
343  // charge exchange corrections to leptonic pion production
344  // in the (3,3) resonance region," Phys. Rev. D 9 (1974)
345  // 2125-2143 [Erratum Phys. Rev. D 10 (1974) 1669].
346  // [2] J.Y. Yu, "Neutrino interactions and nuclear effects in
347  // oscillation experiments and the nonperturbative disper-
348  // sive sector in strong (quasi-)abelian fields," Ph. D.
349  // Thesis, Dortmund U., Dortmund, 2002 (unpublished).
350  // [3] E.A. Paschos, J.Y. Yu, and M. Sakuda, "Neutrino pro-
351  // duction of resonances," Phys. Rev. D 69 (2004) 014013
352  // [arXiv: hep-ph/0308130].
353 
354  double P_Fermi = 0.0;
355 
356  // Maximum value of Fermi momentum of target nucleon (GeV)
357  if (A<6 || !fUseRFGParametrization)
358  {
359  // Look up the Fermi momentum for this target
361  const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
362  P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucpdgc);
363  }
364  else {
365  // Define the Fermi momentum for this target
367  // Correct the Fermi momentum for the struck nucleon
368  if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); }
369  else { P_Fermi *= TMath::Power( 2.*N/A, 1./3); }
370  }
371 
372  double FactorPauli_RES = 1.0;
373 
374  double k0 = 0., q = 0., q0 = 0.;
375 
376  if (P_Fermi > 0.)
377  {
378  k0 = (W2-Mnuc2-Q2)/(2*W);
379  k = TMath::Sqrt(k0*k0+Q2); // previous value of k is overridden
380  q0 = (W2-Mnuc2+kPionMass2)/(2*W);
381  q = TMath::Sqrt(q0*q0-kPionMass2);
382  }
383 
384  if (2*P_Fermi < k-q)
385  FactorPauli_RES = 1.0;
386  if (2*P_Fermi >= k+q)
387  FactorPauli_RES = ((3*k*k+q*q)/(2*P_Fermi)-(5*TMath::Power(k,4)+TMath::Power(q,4)+10*k*k*q*q)/(40*TMath::Power(P_Fermi,3)))/(2*k);
388  if (2*P_Fermi >= k-q && 2*P_Fermi <= k+q)
389  FactorPauli_RES = ((q+k)*(q+k)-4*P_Fermi*P_Fermi/5-TMath::Power(k-q, 3)/(2*P_Fermi)+TMath::Power(k-q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*q*k);
390 
391  xsec *= FactorPauli_RES;
392  }
393 
394  return xsec;
395 }
396 //____________________________________________________________________________
397 double ReinSehgalRESPXSec::Integral(const Interaction * interaction) const
398 {
399  double xsec = fXSecIntegrator->Integrate(this,interaction);
400  return xsec;
401 }
402 //____________________________________________________________________________
403 bool ReinSehgalRESPXSec::ValidProcess(const Interaction * interaction) const
404 {
405  if(interaction->TestBit(kISkipProcessChk)) return true;
406 
407  const InitialState & init_state = interaction->InitState();
408  const ProcessInfo & proc_info = interaction->ProcInfo();
409  const XclsTag & xcls = interaction->ExclTag();
410 
411  if(!proc_info.IsResonant()) return false;
412  if(!xcls.KnownResonance()) return false;
413 
414  int hitnuc = init_state.Tgt().HitNucPdg();
415  bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc));
416 
417  if (!is_pn) return false;
418 
419  int probe = init_state.ProbePdg();
420  bool is_weak = proc_info.IsWeak();
421  bool is_em = proc_info.IsEM();
422  bool nu_weak = (pdg::IsNeutralLepton(probe) && is_weak);
423  bool l_em = (pdg::IsChargedLepton(probe) && is_em );
424 
425  if (!nu_weak && !l_em) return false;
426 
427  return true;
428 }
429 //____________________________________________________________________________
431 {
432  Algorithm::Configure(config);
433  this->LoadConfig();
434 }
435 //____________________________________________________________________________
437 {
438  Algorithm::Configure(config);
439  this->LoadConfig();
440 }
441 //____________________________________________________________________________
443 {
444  // Cross section scaling factors
445  this->GetParam( "RES-CC-XSecScale", fXSecScaleCC ) ;
446  this->GetParam( "RES-NC-XSecScale", fXSecScaleNC ) ;
447 
448  this->GetParam( "RES-Zeta", fZeta ) ;
449  this->GetParam( "RES-Omega", fOmega ) ;
450 
451  double ma, mv ;
452  this->GetParam( "RES-Ma", ma ) ;
453  this->GetParam( "RES-Mv", mv ) ;
454  fMa2 = TMath::Power(ma,2);
455  fMv2 = TMath::Power(mv,2);
456 
457  this->GetParamDef( "BreitWignerWeight", fWghtBW, true ) ;
458  this->GetParamDef( "BreitWignerNorm", fNormBW, true);
459 
460  double thw ;
461  this->GetParam( "WeinbergAngle", thw ) ;
462  fSin48w = TMath::Power( TMath::Sin(thw), 4 );
463  double Vud;
464  this->GetParam("CKM-Vud", Vud );
465  fVud2 = TMath::Power( Vud, 2 );
466  this->GetParam("FermiMomentumTable", fKFTable);
467  this->GetParam("RFG-UseParametrization", fUseRFGParametrization);
468  this->GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
469 
470  // Load all the sub-algorithms needed
471 
472  fHAmplModelCC = 0;
473  fHAmplModelNCp = 0;
474  fHAmplModelNCn = 0;
475  fHAmplModelEMp = 0;
476  fHAmplModelEMn = 0;
477 
478  AlgFactory * algf = AlgFactory::Instance();
479 
480  fHAmplModelCC = dynamic_cast<const RSHelicityAmplModelI *> (
481  algf->GetAlgorithm("genie::RSHelicityAmplModelCC","Default"));
482  fHAmplModelNCp = dynamic_cast<const RSHelicityAmplModelI *> (
483  algf->GetAlgorithm("genie::RSHelicityAmplModelNCp","Default"));
484  fHAmplModelNCn = dynamic_cast<const RSHelicityAmplModelI *> (
485  algf->GetAlgorithm("genie::RSHelicityAmplModelNCn","Default"));
486  fHAmplModelEMp = dynamic_cast<const RSHelicityAmplModelI *> (
487  algf->GetAlgorithm("genie::RSHelicityAmplModelEMp","Default"));
488  fHAmplModelEMn = dynamic_cast<const RSHelicityAmplModelI *> (
489  algf->GetAlgorithm("genie::RSHelicityAmplModelEMn","Default"));
490 
496 
497  // Use algorithm within a DIS/RES join scheme. If yes get Wcut
498  this->GetParam( "UseDRJoinScheme", fUsingDisResJoin ) ;
499  fWcut = 999999;
500  if(fUsingDisResJoin) {
501  this->GetParam( "Wcut", fWcut ) ;
502  }
503 
504  // NeuGEN limits in the allowed resonance phase space:
505  // W < min{ Wmin(physical), (res mass) + x * (res width) }
506  // It limits the integration area around the peak and avoids the
507  // problem with huge xsec increase at low Q2 and high W.
508  // In correspondence with Hugh, Rein said that the underlying problem
509  // are unphysical assumptions in the model.
510  this->GetParamDef( "MaxNWidthForN2Res", fN2ResMaxNWidths, 2.0 ) ;
511  this->GetParamDef( "MaxNWidthForN0Res", fN0ResMaxNWidths, 6.0 ) ;
512  this->GetParamDef( "MaxNWidthForGNRes", fGnResMaxNWidths, 4.0 ) ;
513 
514  // NeuGEN reduction factors for nu_tau: a gross estimate of the effect of
515  // neglected form factors in the R/S model
516  this->GetParamDef( "UseNuTauScalingFactors", fUsingNuTauScaling, true ) ;
517  if(fUsingNuTauScaling) {
518  if(fNuTauRdSpl) delete fNuTauRdSpl;
519  if(fNuTauBarRdSpl) delete fNuTauBarRdSpl;
520 
521  assert( std::getenv( "GENIE") );
522  string base = std::getenv( "GENIE") ;
523 
524  string filename = base + "/data/evgen/rein_sehgal/res/nutau_xsec_scaling_factors.dat";
525  LOG("ReinSehgalRes", pNOTICE)
526  << "Loading nu_tau xsec reduction spline from: " << filename;
527  fNuTauRdSpl = new Spline(filename);
528 
529  filename = base + "/data/evgen/rein_sehgal/res/nutaubar_xsec_scaling_factors.dat";
530  LOG("ReinSehgalRes", pNOTICE)
531  << "Loading bar{nu_tau} xsec reduction spline from: " << filename;
532  fNuTauBarRdSpl = new Spline(filename);
533  }
534 
535  // load the differential cross section integrator
537  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
539 }
540 //____________________________________________________________________________
541 
bool IsDelta(Resonance_t res)
is it a Delta resonance?
bool IsResonant(void) const
Definition: ProcessInfo.cxx:92
Cross Section Calculation Interface.
virtual const RSHelicityAmpl & Compute(Resonance_t res, const FKR &fkr) const =0
const double kPi
double W(bool selected=false) const
Definition: Kinematics.cxx:167
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
bool IsNuTau(int pdgc)
Definition: PDGUtils.cxx:160
Basic constants.
bool fNormBW
normalize resonance breit-wigner to 1?
bool IsWeak(void) const
bool IsWeakCC(void) const
static const double kSqrt2
Definition: Constants.h:116
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
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
double Rminus
Definition: FKR.h:51
Spline * fNuTauBarRdSpl
xsec reduction spline for nu_tau_bar
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
bool fWghtBW
weight with resonance breit-wigner?
int HitNucPdg(void) const
Definition: Target.cxx:321
const RSHelicityAmplModelI * fHAmplModelEMn
const RSHelicityAmplModelI * fHAmplModelNCp
double Ra
Definition: FKR.h:43
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
int A(void) const
Definition: Target.h:71
bool KnownResonance(void) const
Definition: XclsTag.h:61
double HitNucMass(void) const
Definition: Target.cxx:250
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
static FermiMomentumTablePool * Instance(void)
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
bool IsAntiNuTau(int pdgc)
Definition: PDGUtils.cxx:175
double Lamda
Definition: FKR.h:38
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:99
double Mass(Resonance_t res)
resonance mass (GeV)
double R
Definition: FKR.h:46
A table of Fermi momentum constants.
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition: config.py:1
double Width(Resonance_t res)
resonance width (GeV)
string filename
Definition: shutoffs.py:106
double fMv2
(vector mass)^2
double BreitWignerLGamma(double W, int L, double mass, double width0, double norm)
Definition: BWFunc.cxx:30
double Evaluate(double x) const
Definition: Spline.cxx:362
double BreitWignerL(double W, int L, double mass, double width0, double norm)
Definition: BWFunc.cxx:107
enum genie::EKinePhaseSpace KinePhaseSpace_t
const RSHelicityAmplModelI * fHAmplModelCC
double BWNorm(Resonance_t res, double N0ResMaxNWidths=6, double N2ResMaxNWidths=2, double GnResMaxNWidths=4)
breit-wigner normalization factor
enum genie::EResonance Resonance_t
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
const RSHelicityAmplModelI * fHAmplModelNCn
string AsString(void) const
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:37
double Integral(const Interaction *i) const
const RSHelicityAmplModelI * fHAmplModelEMp
Float_t Z
Definition: plot.C:38
Double_t q2[12][num]
Definition: f2_nu.C:137
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
double fZeta
FKR parameter Zeta.
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:140
Summary information for an interaction.
Definition: Interaction.h:56
double Tv
Definition: FKR.h:39
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
double q2(bool selected=false) const
Definition: Kinematics.cxx:151
A class holding the Rein-Sehgal&#39;s helicity amplitudes.
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
bool IsWeakNC(void) const
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 double NR
Float_t E
Definition: plot.C:20
const FermiMomentumTable * GetTable(string name)
static const double kAem2
Definition: Constants.h:58
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
double fWcut
apply DIS/RES joining scheme < Wcut
std::string getenv(std::string const &name)
double T
Definition: FKR.h:47
double Rv
Definition: FKR.h:40
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
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
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
Float_t d
Definition: plot.C:236
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
int Z(void) const
Definition: Target.h:69
#define pINFO
Definition: Messenger.h:63
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
Resonance_t Resonance(void) const
Definition: XclsTag.h:62
Double_t xsec[nknots]
Definition: testXsec.C:47
bool IsEM(void) const
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:93
double C
Definition: FKR.h:45
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
static const double A
Definition: Units.h:82
double Tplus
Definition: FKR.h:48
const XSecIntegratorI * fXSecIntegrator
double fXSecScaleCC
external CC xsec scaling factor
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
double B
Definition: FKR.h:44
double fOmega
FKR parameter Omega.
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double Rplus
Definition: FKR.h:50
bool fUseRFGParametrization
use parametrization for fermi momentum insted of table?
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool fUsingDisResJoin
use a DIS/RES joining scheme?
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
double Tminus
Definition: FKR.h:49
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double fSin48w
sin^4(Weingberg angle)
double fN0ResMaxNWidths
limits allowed phase space for n=0 res
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
const InitialState & InitState(void) const
Definition: Interaction.h:69
const char * AsString(Resonance_t res)
resonance id -> string
double fGnResMaxNWidths
limits allowed phase space for other res
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
void Configure(const Registry &config)
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
double fXSecScaleNC
external NC xsec scaling factor
double fMa2
(axial mass)^2
bool fUsingNuTauScaling
use NeuGEN nutau xsec reduction factors?
static const double kGF2
Definition: Constants.h:60
string fKFTable
table of Fermi momentum (kF) constants for various nuclei
void kinematics()
Definition: kinematics.C:10
#define W(x)
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:131
static const double kPi2
Definition: Constants.h:39
double S
Definition: FKR.h:41
double Ta
Definition: FKR.h:42
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
int ResonanceIndex(Resonance_t res)
resonance idx, quark model / SU(6)
bool fUsePauliBlocking
account for Pauli blocking?
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double fN2ResMaxNWidths
limits allowed phase space for n=2 res
Spline * fNuTauRdSpl
xsec reduction spline for nu_tau
static const double kPionMass2
Definition: Constants.h:87
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353