GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BergerSehgalFMCOHPiPXSec2015.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::BergerSehgalFMCOHPiPXSec2015")
32 {
33 
34 }
35 //____________________________________________________________________________
37  XSecAlgorithmI("genie::BergerSehgalFMCOHPiPXSec2015", 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  // We have additionally modified the formulae to account for a
55  // non-infinite mass for the target nucleus
56 
57  if(! this -> ValidProcess (interaction) ) return 0.;
58  if(! this -> ValidKinematics (interaction) ) return 0.;
59 
60  const Kinematics & kinematics = interaction -> Kine();
61  const InitialState & init_state = interaction -> InitState();
62 
63  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
64  double xsec = 0.0;
65 
66  double E = init_state.ProbeE(kRfLab); // nu E
67  double Q2 = kinematics.Q2();
68  double y = kinematics.y(); // inelasticity
69  double t = kinematics.t(); // fun exists?
70  assert(E > 0.);
71  assert(y > 0.);
72  assert(y < 1.);
73  double ppistar = PionCOMAbsMomentum(interaction); // |Center of Mass Momentum|
74  if (ppistar <= 0.0) {
75  LOG("BergerSehgalFMCohPi", pDEBUG) <<
76  "Pion COM momentum negative for Q2 = " << Q2 <<
77  " y = " << y;
78  return 0.0;
79  }
80  double front = ExactKinematicTerm(interaction);
81  if (front <= 0.0) {
82  LOG("BergerSehgalFMCohPi", pDEBUG) << "Exact kin. form = " << front <<
83  " E = " << E << " Q2 = " << Q2 << " y = " << y;
84  return 0.0;
85  }
86 
87  double A = (double) init_state.Tgt().A(); // mass number
88  double A2 = TMath::Power(A, 2.);
89  double A_3 = TMath::Power(A, 1./3.);
90  double M = init_state.Tgt().Mass();
91  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
92  double M_pi2 = M_pi * M_pi;
93  double Epi = y * E - t / (2 * M);
94  double Epi2 = Epi * Epi;
95  double ma2 = fMa * fMa;
96  double Ga = ma2 / (ma2 + Q2);
97  double Ga2 = Ga * Ga;
98  double Ro2 = TMath::Power(fRo * units::fermi, 2.);
99  double ppi2 = Epi2 - M_pi2;
100  double ppi = ppi2 > 0.0 ? sqrt(ppi2) : 0.0;
101  // double fp = 0.93 * kPionMass; // unused // pion decay constant
102 
103  double costheta = (t - Q2 - M_pi * M_pi) / (2 * ( (y *E - Epi) * Epi -
104  ppi * sqrt(TMath::Power(y * E - Epi, 2.) + t) ) );
105 
106  if ((costheta > 1.0) || (costheta < -1.0)) return 0.0;
107 
108  // tot. pi+N xsec
109  double sTot =
110  utils::hadxs::berger::TotalPionNucleonXSec(Epi, pionIsCharged);
111  double sTot2 = sTot * sTot;
112  // inel. pi+N xsec
113  double sInel =
115 
116  // Fabs (F_{abs}) describes the average attenuation of a pion emerging
117  // from a sphere of nuclear matter with radius = R_0 A^{1/3}. it is
118  // Eq. 13 in Berger-Sehgal PRD 79, 053003
119  double Fabs_input = (9.0 * A_3) / (16.0 * kPi * Ro2);
120  double Fabs = TMath::Exp( -1.0 * Fabs_input * sInel);
121 
122  // A_RS for BS version of RS, and/or Tpi>1.0
123  //double RS_factor = (A2 * Fabs) / (16.0 * kPi) * (sTot2);
124  double R = fRo * A_3 * units::fermi; // nuclear radius
125  double R2 = R * R; //
126  double b = 0.33333 * R2; // Eq. 12 in BS
127  double expbt = TMath::Exp( -b * t );
128  double dsigEldt = sTot2 / (16. * kPi); // Eq. 11 in BS
129  double dsigpiNdt = A2 * dsigEldt * expbt * Fabs; // Eq. 10 in BS
130 
131  double tpilow = 0.0;
132  double siglow = 0.0;
133  double tpihigh = 0.0;
134  double sighigh = 0.0;
135  double dsigdt = 0.0;
136  double tpi = 0.0;
137  int xsec_stat = 0;
138 
139  // differential cross section for pion-nucleus in t (ds/dt term from
140  // Eq. 7 in BS. we will initially set the value to a "Rein-Sehgal style"
141  // computation and update to use the Berger-Sehgal pion-nucleus cross
142  // section where appropriate.
143  double edep_dsigpiNdt = dsigpiNdt;
144 
145  // c.o.m.
146  tpi = Epi - M_pi;
147 
148  if (tpi <= 1.0 && fRSPionXSec == false) {
149  // use the Berger-Sehgal pion-nucleus cross section. note we're only
150  // checking on the pion energy and the conditional flag - is it really
151  // reasonable to ever use this value for non-Carbon targets?
152  xsec_stat =
154  tpi, ppistar, t, A,
155  tpilow, siglow,
156  tpihigh, sighigh);
157  if (xsec_stat != 0)
158  LOG("BergerSehgalFMCohPi", pWARN) <<
159  "Unable to retrieve pion-nucleus cross section with A = " <<
160  A << ", t_pi = " << tpi;
161  dsigdt = siglow + (sighigh - siglow) * (tpi - tpilow) / (tpihigh - tpilow);
162  dsigdt = dsigdt / (2.0 * ppistar * ppistar) * units::mb;
163  edep_dsigpiNdt = dsigdt;
164  }
165 
166  // complete calculation of Eq. 7 in BS paper
167  xsec = front * Ga2 * edep_dsigpiNdt;
168 
169  // Correction for finite final state lepton mass.
170  // Lepton mass modification is part of Berger-Sehgal and is not optional.
171  if (pionIsCharged) {
172  double C = 1.;
173  // First, we need to remove the leading G_{A}^2 which is required for NC.
174  xsec /= Ga2;
175  // Next, \cos^2 \theta_{Cabibbo} appears in the CC xsec, but not the NC.
176  xsec *= fCos8c2;
177  double ml = interaction->FSPrimLepton()->Mass();
178  double ml2 = TMath::Power(ml,2);
179  double Q2min = ml2 * y/(1-y);
180  if(Q2 > Q2min) {
181  double C1 = TMath::Power(Ga - 0.5 * Q2min / (Q2 + kPionMass2), 2);
182  double C2 = 0.25 * y * Q2min * (Q2 - Q2min) /
183  TMath::Power(Q2 + kPionMass2, 2);
184  C = C1 + C2;
185  } else {
186  C = 0.;
187  }
188  xsec *= (2. * C); // *2 is for CC vs NC in BS
189  }
190 
191 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
192  LOG("BergerSehgalFMCohPi", pDEBUG)
193  << "\n momentum transfer .............. Q2 = " << Q2
194  << "\n mass number .................... A = " << A
195  << "\n pion energy .................... Epi = " << Epi
196  << "\n propagator term ................ propg = " << propg
197  << "\n Re/Im of fwd pion scat. ampl. .. Re/Im = " << fReIm
198  << "\n total pi+N cross section ....... sigT = " << sTot
199  << "\n inelastic pi+N cross section ... sigI = " << sInel
200  << "\n nuclear size scale ............. Ro = " << fRo
201  << "\n pion absorption factor ......... Fabs = " << Fabs
202  << "\n t integration factor ........... tint = " << tint;
203  LOG("BergerSehgalFMCohPi", pINFO)
204  << "d3xsec/dQ2dydt[COHPi] (x= " << x << ", y="
205  << y << ", E=" << E << ") = "<< xsec;
206 #endif
207 
208  //----- The algorithm computes d^3xsec/dQ^2dydt
209  // Check whether Jacobian tranformation is needed...
210 
211  return xsec;
212 }
213 //____________________________________________________________________________
215  const Interaction * interaction) const
216 {
217  // This function is a bit inefficient but is being encapsulated as
218  // such in order to possibly migrate into a general kinematics check.
219  const Kinematics & kinematics = interaction -> Kine();
220  const InitialState & init_state = interaction -> InitState();
221 
222  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
223  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
224  double E = init_state.ProbeE(kRfLab); // nu E
225  double Q2 = kinematics.Q2();
226  double y = kinematics.y(); // inelasticity
227  double fp2 = (0.93 * M_pi)*(0.93 * M_pi);
228 
229  double term = ((kGF2 * fp2) / (4.0 * kPi2)) *
230  ((E * (1.0 - y)) / sqrt(y*E * y*E + Q2)) *
231  (1.0 - Q2 / (4.0 * E*E * (1.0 - y)));
232  return term;
233 }
234 //____________________________________________________________________________
236  const Interaction * interaction) const
237 {
238  // This function is a bit inefficient but is being encapsulated as
239  // such in order to possibly migrate into a general kinematics check.
240  const Kinematics & kinematics = interaction -> Kine();
241  const InitialState & init_state = interaction -> InitState();
242 
243  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
244  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
245  double E = init_state.ProbeE(kRfLab); // nu E
246  double Q2 = kinematics.Q2();
247  double y = kinematics.y(); // inelasticity
248  double MT = init_state.Tgt().Mass();
249 
250  double W2 = MT * MT - Q2 + 2.0 * y * E * MT;
251  double arg = (2.0 * MT * (y * E - M_pi) - Q2 - M_pi * M_pi) *
252  (2.0 * MT * (y * E + M_pi) - Q2 - M_pi * M_pi);
253  if (arg < 0) return arg;
254  double ppistar = TMath::Sqrt(arg) / 2.0 / TMath::Sqrt(W2);
255 
256  return ppistar;
257 }
258 //____________________________________________________________________________
259 double BergerSehgalFMCOHPiPXSec2015::Integral(const Interaction * interaction) const
260 {
261  double xsec = fXSecIntegrator->Integrate(this,interaction);
262  return xsec;
263 }
264 //____________________________________________________________________________
266 {
267  if(interaction->TestBit(kISkipProcessChk)) return true;
268 
269  const InitialState & init_state = interaction->InitState();
270  const ProcessInfo & proc_info = interaction->ProcInfo();
271  const Target & target = init_state.Tgt();
272 
273  int nu = init_state.ProbePdg();
274 
275  if (!proc_info.IsCoherentProduction()) return false;
276  if (!proc_info.IsWeak()) return false;
277  if (target.HitNucIsSet()) return false;
278  if (!(target.A()>1)) return false;
279  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
280 
281  return true;
282 }
283 //____________________________________________________________________________
285 {
286  Algorithm::Configure(config);
287  this->LoadConfig();
288 }
289 //____________________________________________________________________________
291 {
292  Algorithm::Configure(config);
293  this->LoadConfig();
294 }
295 //____________________________________________________________________________
297 {
298  GetParam( "COH-Ma",fMa ) ;
299  GetParam( "COH-Ro", fRo ) ;
300 
301  double thc ;
302  GetParam( "CabibboAngle", thc ) ;
303  fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
304 
305  // fRSPionXSec => Do not use the pion-nucleus cross section from Table 1 in PRD 79, 053003
306  // Instead, use the Rein-Sehgal "style" pion-nucleon cross section and scale by A
307  // for all pion energies.
308  GetParam( "COH-UseRSPionXSec", fRSPionXSec ) ;
309 
310  //-- load the differential cross section integrator
312  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
313  assert(fXSecIntegrator);
314 }
315 //____________________________________________________________________________
Cross Section Calculation Interface.
bool IsWeak(void) const
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int A(void) const
Definition: Target.h:70
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:90
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
bool IsCoherentProduction(void) const
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Mass(void) const
Definition: Target.cxx:224
double y(bool selected=false) const
Definition: Kinematics.cxx:112
static constexpr double b
Definition: Units.h:78
double ExactKinematicTerm(const Interaction *i) const
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?
bool fRSPionXSec
Use Rein-Sehgal &quot;style&quot; pion-nucleon xsecs.
#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
double PionCOMAbsMomentum(const Interaction *i) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
double fRo
nuclear size scale parameter
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
double Integral(const Interaction *i) const
bool HitNucIsSet(void) const
Definition: Target.cxx:283
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:82
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 t(bool selected=false) const
Definition: Kinematics.cxx:170
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
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