GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BergerSehgalCOHPiPXSec2015.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  G. Perdue, H. Gallagher, D. Cherdack
7 */
8 //____________________________________________________________________________
9 
10 #include <TMath.h>
11 
14 #include "Framework/Conventions/GBuild.h"
24 
25 using namespace genie;
26 using namespace genie::constants;
27 using namespace genie::utils;
28 
29 //____________________________________________________________________________
31  XSecAlgorithmI("genie::BergerSehgalCOHPiPXSec2015")
32 {
33 
34 }
35 //____________________________________________________________________________
37  XSecAlgorithmI("genie::BergerSehgalCOHPiPXSec2015", config)
38 {
39 
40 }
41 //____________________________________________________________________________
43 {
44 
45 }
46 //____________________________________________________________________________
48  const Interaction * interaction, KinePhaseSpace_t kps) const
49 {
50  // Here we are following PRD 79, 053003 (2009) by Berger and Sehgal
51  // This method computes the differential cross section represented
52  // in Eq.'s 6 (CC) and 7 (NC) from that paper.
53 
54  if(! this -> ValidProcess (interaction) ) return 0.;
55  if(! this -> ValidKinematics (interaction) ) return 0.;
56 
57  const Kinematics & kinematics = interaction -> Kine();
58  const InitialState & init_state = interaction -> InitState();
59 
60  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
61  double xsec = 0.0;
62 
63  double E = init_state.ProbeE(kRfLab); // nu E
64  double Q2 = kinematics.Q2();
65  double y = kinematics.y(); // inelasticity
66  double x = kinematics.x();
67  assert(E > 0.);
68  assert(y > 0.);
69  assert(y < 1.);
70  double ppistar = PionCOMAbsMomentum(interaction); // |Center of Mass Momentum|
71  if (ppistar <= 0.0) {
72  LOG("BergerSehgalCohPi", pDEBUG) << "Pion COM momentum negative for Q2 = " << Q2 <<
73  " x = " << x << " y = " << y;
74  return 0.0;
75  }
76  double front = ExactKinematicTerm(interaction);
77  if (front <= 0.0) {
78  LOG("BergerSehgalCohPi", pDEBUG) << "Exact kin. form = " << front <<
79  " E = " << E << " Q2 = " << Q2 << " y = " << y << " x = " << x;
80  return 0.0;
81  }
82 
83  double A = (double) init_state.Tgt().A(); // mass number
84  double A2 = TMath::Power(A, 2.);
85  double A_3 = TMath::Power(A, 1./3.);
86  double M = init_state.Tgt().Mass();
87  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
88  double Epi = y*E; // ~pion energy
89  double ma2 = TMath::Power(fMa, 2); // "axial mass" squared
90  double Ga = ma2 / (ma2 + Q2);
91  double Ga2 = TMath::Power(Ga, 2.); // propagator term
92  double Ro2 = TMath::Power(fRo * units::fermi, 2.);
93 
94  // the xsec is d^3xsec/dQ^2dydt but the only t-dependent factor
95  // is an exp(-bt) so it can be integrated analyticaly
96  double Epi2 = TMath::Power(Epi, 2.);
97  double R = fRo * A_3 * units::fermi; // nuclear radius
98  double R2 = TMath::Power(R, 2.);
99  double b = 0.33333 * R2;
100  double MxEpi = M * x / Epi;
101  double mEpi2 = (M_pi * M_pi) / Epi2;
102  double tA = 1. + MxEpi - 0.5 * mEpi2;
103  double tB = TMath::Sqrt(1.0 + 2 * MxEpi) * TMath::Sqrt(1.0 - mEpi2);
104  double tmin = 2 * Epi2 * (tA - tB);
105  double tmax = 2 * Epi2 * (tA + tB);
106  if (tmin < 1.0e-8) {
107  tmin = 1.0e-8;
108  }
109 
110  /* const KPhaseSpace & kphase = interaction->PhaseSpace(); */
111  /* Range1D_t tl = kphase.TLim(); // TESTING! */
112 
113  double sigtot_pin = utils::hadxs::berger::PionNucleonXSec(Epi, /* get_total = */ true, pionIsCharged);
114  double sigel_pin = utils::hadxs::berger::PionNucleonXSec(Epi, /* get_total = */ false, pionIsCharged);
115  double siginel_pin = sigtot_pin - sigel_pin;
116 
117  // fabs (F_{abs}) describes the average attenuation of a pion emerging
118  // from a sphere of nuclear matter with radius = R_0 A^{1/3}. it is
119  // Eq. 13 in Berger-Sehgal PRD 79, 053003
120  double fabs_input = (9.0 * A_3) / (16.0 * kPi * Ro2);
121  double fabs = TMath::Exp( -1.0 * fabs_input * siginel_pin);
122 
123  // my old hackery to get things to work, A. Mislivec provided a better alt.
124  // double factor = 0.1; // to go from 10^-37 cm^2 -> 10^-38 cm^2
125  // double RS_factor = (units::mb*A2*fabs)/(16.0*kPi) * (sigtot_pin*sigtot_pin);
126 
127  // A_RS for BS version of RS, and/or Tpi>1.0
128  double RS_factor = (A2 * fabs) / (16.0 * kPi) * (sigtot_pin * sigtot_pin);
129 
130  // get the pion-nucleus cross section on carbon, fold it into differential cross section
131  double tpi = (E * y) - M_pi - ((Q2 + M_pi * M_pi) / (2 * M));
132  double tpilow = 0.0;
133  double siglow = 0.0;
134  double tpihigh = 0.0;
135  double sighigh = 0.0;
136  double dsigdzfit = 0.0;
137  double dsigdtfit = 0.0;
138  int xsec_stat = 0;
139  double dsig = 0.0;
140  double tstep = 100;
141  double logt_step = TMath::Abs(log(tmax) - log(tmin)) / tstep;
142  double logt = log(tmin) - logt_step/2.0;
143  double t_itt = TMath::Exp(logt);
144  double t_width = 0.0;
145 
146  for (double t_step = 0; t_step<tstep; t_step++) {
147 
148  logt = logt + logt_step;
149  t_itt = TMath::Exp(logt);
150  t_width = t_itt*logt_step;
151 
152  if (tpi <= 1.0 && fRSPionXSec == false) {
153  xsec_stat = utils::hadxs::berger::PionNucleusXSec(tpi, ppistar, t_itt, A, tpilow, siglow, tpihigh, sighigh);
154  if(xsec_stat){
155  LOG("BergerSehgalCohPi", pERROR) << "Call to PionNucleusXSec code failed - return xsec of 0.0";
156  return 0.0;
157  }
158  dsigdzfit = siglow + (sighigh - siglow) * (tpi - tpilow) / (tpihigh - tpilow);
159  dsigdtfit = dsigdzfit / (2.0 * ppistar * ppistar);
160  // we are handed a cross section in mb, need to convert it to GeV^{-2}
161  dsig += 1.0 * front * Ga2 * t_width * dsigdtfit * units::mb;
162  }
163  else {
164  dsig += /*factor **/ front * Ga2 * t_width * RS_factor * exp(-1.0*b*t_itt);
165  }
166 
167  }
168  xsec = dsig;
169 
170  // Correction for finite final state lepton mass.
171  // Lepton mass modification is part of Berger-Sehgal and is not optional.
172  if (pionIsCharged) {
173  double C = 1.;
174  // First, we need to remove the leading G_{A}^2 which is required for NC.
175  xsec /= Ga2;
176  // Next, \cos^2 \theta_{Cabibbo} appears in the CC xsec, but not the NC.
177  xsec *= fCos8c2;
178  double ml = interaction->FSPrimLepton()->Mass();
179  double ml2 = TMath::Power(ml,2);
180  double Q2min = ml2 * y/(1-y);
181  if(Q2 > Q2min) {
182  double C1 = TMath::Power(Ga - 0.5 * Q2min / (Q2 + kPionMass2), 2);
183  double C2 = 0.25 * y * Q2min * (Q2 - Q2min) /
184  TMath::Power(Q2 + kPionMass2, 2);
185  C = C1 + C2;
186  } else {
187  C = 0.;
188  }
189  xsec *= (2. * C); // *2 is for CC vs NC in BS
190  }
191 
192 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
193  LOG("BergerSehgalCohPi", pDEBUG)
194  << "\n momentum transfer .............. Q2 = " << Q2
195  << "\n mass number .................... A = " << A
196  << "\n pion energy .................... Epi = " << Epi
197  << "\n propagator term ................ Ga2 = " << Ga2
198  << "\n total pi+N cross section ....... sigT = " << sigtot_pin
199  << "\n inelastic pi+N cross section ... sigI = " << siginel_pin
200  << "\n nuclear size scale ............. Ro = " << fRo
201  << "\n pion absorption factor ......... Fabs = " << fabs
202  << "\n t integration range ............ [" << tmin << "," << tmax << "]"
203  LOG("BergerSehgalCohPi", pINFO)
204  << "d2xsec/dQ2dy[COHPi] (Q2= " << Q2 << ", y="
205  << y << ", E=" << E << ") = "<< xsec;
206 #endif
207 
208  //----- The algorithm computes d^2xsec/dQ2dy
209  // Check whether variable tranformation is needed? May be working with logs.
210  // kPSlogQ2logyfE is possible - all others will not succeed
211  if(kps != kPSQ2yfE) {
212  double J = utils::kinematics::Jacobian(interaction,kPSQ2yfE, kps);
213  xsec *= J;
214  }
215  return xsec;
216 }
217 //____________________________________________________________________________
219 {
220  // This function is a bit inefficient but is being encapsulated as
221  // such in order to possibly migrate into a general kinematics check.
222  const Kinematics & kinematics = interaction -> Kine();
223  const InitialState & init_state = interaction -> InitState();
224 
225  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
226  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
227  double E = init_state.ProbeE(kRfLab); // nu E
228  double Q2 = kinematics.Q2();
229  double y = kinematics.y(); // inelasticity
230  double fp2 = (0.93 * M_pi)*(0.93 * M_pi);
231 
232  double term = ((kGF2 * fp2) / (4.0 * kPi2)) *
233  ((E * (1.0 - y)) / sqrt(y*E * y*E + Q2)) *
234  (1.0 - Q2 / (4.0 * E*E * (1.0 - y)));
235  return term;
236 }
237 //____________________________________________________________________________
239 {
240  // This function is a bit inefficient but is being encapsulated as
241  // such in order to possibly migrate into a general kinematics check.
242  const Kinematics & kinematics = interaction -> Kine();
243  const InitialState & init_state = interaction -> InitState();
244 
245  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
246  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
247  double E = init_state.ProbeE(kRfLab); // nu E
248  double Q2 = kinematics.Q2();
249  double y = kinematics.y(); // inelasticity
250  double MT = init_state.Tgt().Mass();
251 
252  double W2 = MT*MT - Q2 + 2.0 * y * E * MT;
253  double arg = (2.0*MT*(y*E - M_pi) - Q2 - M_pi*M_pi)*(2.0*MT*(y*E + M_pi) - Q2 - M_pi*M_pi);
254  if (arg < 0) return arg;
255  double ppistar = TMath::Sqrt(arg) / 2.0 / TMath::Sqrt(W2);
256 
257  return ppistar;
258 }
259 //____________________________________________________________________________
260 double BergerSehgalCOHPiPXSec2015::Integral(const Interaction * interaction) const
261 {
262  double xsec = fXSecIntegrator->Integrate(this,interaction);
263  return xsec;
264 }
265 //____________________________________________________________________________
267 {
268  if(interaction->TestBit(kISkipProcessChk)) return true;
269 
270  const InitialState & init_state = interaction->InitState();
271  const ProcessInfo & proc_info = interaction->ProcInfo();
272  const Target & target = init_state.Tgt();
273 
274  int nu = init_state.ProbePdg();
275 
276  if (!proc_info.IsCoherentProduction()) return false;
277  if (!proc_info.IsWeak()) return false;
278  if (target.HitNucIsSet()) return false;
279  if (!(target.A()>1)) return false;
280  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
281 
282  return true;
283 }
284 //____________________________________________________________________________
286 {
287  Algorithm::Configure(config);
288  this->LoadConfig();
289 }
290 //____________________________________________________________________________
292 {
293  Algorithm::Configure(config);
294  this->LoadConfig();
295 }
296 //____________________________________________________________________________
298 {
299  GetParam( "COH-Ma",fMa ) ;
300  GetParam( "COH-Ro", fRo ) ;
301 
302  double thc ;
303  GetParam( "CabibboAngle", thc ) ;
304  fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
305 
306  // fRSPionXSec => Do not use the pion-nucleus cross section from Table 1 in PRD 79, 053003
307  // Instead, use the Rein-Sehgal "style" pion-nucleon cross section and scale by A
308  // for all pion energies.
309  GetParam( "COH-UseRSPionXSec", fRSPionXSec ) ;
310 
311  //-- load the differential cross section integrator
313  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
314  assert(fXSecIntegrator);
315 }
316 //____________________________________________________________________________
Cross Section Calculation Interface.
bool IsWeak(void) const
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
#define pERROR
Definition: Messenger.h:59
Cross Section Integrator Interface.
bool fRSPionXSec
Use Rein-Sehgal &quot;style&quot; pion-nucleon xsecs.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int A(void) const
Definition: Target.h:70
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
Definition: HadXSUtils.cxx:96
double fRo
nuclear size scale parameter
double x(bool selected=false) const
Definition: Kinematics.cxx:99
bool IsCoherentProduction(void) const
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Mass(void) const
Definition: Target.cxx:224
double PionCOMAbsMomentum(const Interaction *i) const
double y(bool selected=false) const
Definition: Kinematics.cxx:112
static constexpr double b
Definition: Units.h:78
Summary information for an interaction.
Definition: Interaction.h:56
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
const double e
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
static constexpr double mb
Definition: Units.h:79
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
#define pINFO
Definition: Messenger.h:62
double ExactKinematicTerm(const Interaction *i) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool HitNucIsSet(void) const
Definition: Target.cxx:283
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
static constexpr double fermi
Definition: Units.h:55
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double Integral(const Interaction *i) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
int PionNucleusXSec(double tpi, double ppistar, double t_new, double A, double &tpilow, double &siglow, double &tpihigh, double &sighigh)
Definition: HadXSUtils.cxx:205
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345