BergerSehgalFMCOHPiPXSec2015.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  For the class documentation see the corresponding header file.
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
25 
26 using namespace genie;
27 using namespace genie::constants;
28 using namespace genie::utils;
29 
30 //____________________________________________________________________________
32  XSecAlgorithmI("genie::BergerSehgalFMCOHPiPXSec2015")
33 {
34 
35 }
36 //____________________________________________________________________________
38  XSecAlgorithmI("genie::BergerSehgalFMCOHPiPXSec2015", config)
39 {
40 
41 }
42 //____________________________________________________________________________
44 {
45 
46 }
47 //____________________________________________________________________________
49  const Interaction * interaction, KinePhaseSpace_t /*kps*/) const
50 {
51  // Here we are following PRD 79, 053003 (2009) by Berger and Sehgal
52  // This method computes the differential cross section represented
53  // in Eq.'s 6 (CC) and 7 (NC) from that paper.
54 
55  // We have additionally modified the formulae to account for a
56  // non-infinite mass for the target nucleus
57 
58  if(! this -> ValidProcess (interaction) ) return 0.;
59  if(! this -> ValidKinematics (interaction) ) return 0.;
60 
61  const Kinematics & kinematics = interaction -> Kine();
62  const InitialState & init_state = interaction -> InitState();
63 
64  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
65  double xsec = 0.0;
66 
67  double E = init_state.ProbeE(kRfLab); // nu E
68  double Q2 = kinematics.Q2();
69  double y = kinematics.y(); // inelasticity
70  double t = kinematics.t(); // fun exists?
71  assert(E > 0.);
72  assert(y > 0.);
73  assert(y < 1.);
74  double ppistar = PionCOMAbsMomentum(interaction); // |Center of Mass Momentum|
75  if (ppistar <= 0.0) {
76  LOG("BergerSehgalFMCohPi", pDEBUG) <<
77  "Pion COM momentum negative for Q2 = " << Q2 <<
78  " y = " << y;
79  return 0.0;
80  }
81  double front = ExactKinematicTerm(interaction);
82  if (front <= 0.0) {
83  LOG("BergerSehgalFMCohPi", pDEBUG) << "Exact kin. form = " << front <<
84  " E = " << E << " Q2 = " << Q2 << " y = " << y;
85  return 0.0;
86  }
87 
88  double A = (double) init_state.Tgt().A(); // mass number
89  double A2 = TMath::Power(A, 2.);
90  double A_3 = TMath::Power(A, 1./3.);
91  double M = init_state.Tgt().Mass();
92  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
93  double M_pi2 = M_pi * M_pi;
94  double Epi = y * E - t / (2 * M);
95  double Epi2 = Epi * Epi;
96  double ma2 = fMa * fMa;
97  double Ga = ma2 / (ma2 + Q2);
98  double Ga2 = Ga * Ga;
99  double Ro2 = TMath::Power(fRo * units::fermi, 2.);
100  double ppi2 = Epi2 - M_pi2;
101  double ppi = ppi2 > 0.0 ? sqrt(ppi2) : 0.0;
102  // double fp = 0.93 * kPionMass; // unused // pion decay constant
103 
104  double costheta = (t - Q2 - M_pi * M_pi) / (2 * ( (y *E - Epi) * Epi -
105  ppi * sqrt(TMath::Power(y * E - Epi, 2.) + t) ) );
106 
107  if ((costheta > 1.0) || (costheta < -1.0)) return 0.0;
108 
109  // tot. pi+N xsec
110  double sTot =
111  utils::hadxs::berger::TotalPionNucleonXSec(Epi, pionIsCharged);
112  double sTot2 = sTot * sTot;
113  // inel. pi+N xsec
114  double sInel =
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 * sInel);
122 
123  // A_RS for BS version of RS, and/or Tpi>1.0
124  //double RS_factor = (A2 * Fabs) / (16.0 * kPi) * (sTot2);
125  double R = fRo * A_3 * units::fermi; // nuclear radius
126  double R2 = R * R; //
127  double b = 0.33333 * R2; // Eq. 12 in BS
128  double expbt = TMath::Exp( -b * t );
129  double dsigEldt = sTot2 / (16. * kPi); // Eq. 11 in BS
130  double dsigpiNdt = A2 * dsigEldt * expbt * Fabs; // Eq. 10 in BS
131 
132  double tpilow = 0.0;
133  double siglow = 0.0;
134  double tpihigh = 0.0;
135  double sighigh = 0.0;
136  double dsigdt = 0.0;
137  double tpi = 0.0;
138  int xsec_stat = 0;
139 
140  // differential cross section for pion-nucleus in t (ds/dt term from
141  // Eq. 7 in BS. we will initially set the value to a "Rein-Sehgal style"
142  // computation and update to use the Berger-Sehgal pion-nucleus cross
143  // section where appropriate.
144  double edep_dsigpiNdt = dsigpiNdt;
145 
146  // c.o.m.
147  tpi = Epi - M_pi;
148 
149  if (tpi <= 1.0 && fRSPionXSec == false) {
150  // use the Berger-Sehgal pion-nucleus cross section. note we're only
151  // checking on the pion energy and the conditional flag - is it really
152  // reasonable to ever use this value for non-Carbon targets?
153  xsec_stat =
155  tpi, ppistar, t, A,
156  tpilow, siglow,
157  tpihigh, sighigh);
158  if (xsec_stat != 0)
159  LOG("BergerSehgalFMCohPi", pWARN) <<
160  "Unable to retrieve pion-nucleus cross section with A = " <<
161  A << ", t_pi = " << tpi;
162  dsigdt = siglow + (sighigh - siglow) * (tpi - tpilow) / (tpihigh - tpilow);
163  dsigdt = dsigdt / (2.0 * ppistar * ppistar) * units::mb;
164  edep_dsigpiNdt = dsigdt;
165  }
166 
167  // complete calculation of Eq. 7 in BS paper
168  xsec = front * Ga2 * edep_dsigpiNdt;
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("BergerSehgalFMCohPi", pDEBUG)
194  << "\n momentum transfer .............. Q2 = " << Q2
195  << "\n mass number .................... A = " << A
196  << "\n pion energy .................... Epi = " << Epi
197  << "\n propagator term ................ propg = " << propg
198  << "\n Re/Im of fwd pion scat. ampl. .. Re/Im = " << fReIm
199  << "\n total pi+N cross section ....... sigT = " << sTot
200  << "\n inelastic pi+N cross section ... sigI = " << sInel
201  << "\n nuclear size scale ............. Ro = " << fRo
202  << "\n pion absorption factor ......... Fabs = " << Fabs
203  << "\n t integration factor ........... tint = " << tint;
204  LOG("BergerSehgalFMCohPi", pINFO)
205  << "d3xsec/dQ2dydt[COHPi] (x= " << x << ", y="
206  << y << ", E=" << E << ") = "<< xsec;
207 #endif
208 
209  //----- The algorithm computes d^3xsec/dQ^2dydt
210  // Check whether Jacobian tranformation is needed...
211 
212  return xsec;
213 }
214 //____________________________________________________________________________
216  const Interaction * interaction) const
217 {
218  // This function is a bit inefficient but is being encapsulated as
219  // such in order to possibly migrate into a general kinematics check.
220  const Kinematics & kinematics = interaction -> Kine();
221  const InitialState & init_state = interaction -> InitState();
222 
223  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
224  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
225  double E = init_state.ProbeE(kRfLab); // nu E
226  double Q2 = kinematics.Q2();
227  double y = kinematics.y(); // inelasticity
228  double fp2 = (0.93 * M_pi)*(0.93 * M_pi);
229 
230  double term = ((kGF2 * fp2) / (4.0 * kPi2)) *
231  ((E * (1.0 - y)) / sqrt(y*E * y*E + Q2)) *
232  (1.0 - Q2 / (4.0 * E*E * (1.0 - y)));
233  return term;
234 }
235 //____________________________________________________________________________
237  const Interaction * interaction) const
238 {
239  // This function is a bit inefficient but is being encapsulated as
240  // such in order to possibly migrate into a general kinematics check.
241  const Kinematics & kinematics = interaction -> Kine();
242  const InitialState & init_state = interaction -> InitState();
243 
244  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
245  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
246  double E = init_state.ProbeE(kRfLab); // nu E
247  double Q2 = kinematics.Q2();
248  double y = kinematics.y(); // inelasticity
249  double MT = init_state.Tgt().Mass();
250 
251  double W2 = MT * MT - Q2 + 2.0 * y * E * MT;
252  double arg = (2.0 * MT * (y * E - M_pi) - Q2 - M_pi * M_pi) *
253  (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 BergerSehgalFMCOHPiPXSec2015::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.IsCoherent()) 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"));
315 }
316 //____________________________________________________________________________
317 
Cross Section Calculation Interface.
const double kPi
Basic constants.
bool IsWeak(void) const
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
const XML_Char * target
Definition: expat.h:268
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
static const double fermi
Definition: Units.h:63
int A(void) const
Definition: Target.h:71
T sqrt(T number)
Definition: d0nt_math.hpp:156
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:99
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Mass(void) const
Definition: Target.cxx:241
double y(bool selected=false) const
Definition: Kinematics.cxx:122
double ExactKinematicTerm(const Interaction *i) const
const double C
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 "style" pion-nucleon xsecs.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
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
double PionCOMAbsMomentum(const Interaction *i) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
#define R(x)
#define pINFO
Definition: Messenger.h:63
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
double fRo
nuclear size scale parameter
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static const double A
Definition: Units.h:82
static const double kPionMass
Definition: Constants.h:74
double Integral(const Interaction *i) const
bool HitNucIsSet(void) const
Definition: Target.cxx:300
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const Var kPi0Mass
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
const hit & b
Definition: hits.cxx:21
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:91
assert(nhit_max >=nhit_nbins)
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:180
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
static const double kGF2
Definition: Constants.h:60
static const double mb
Definition: Units.h:87
void kinematics()
Definition: kinematics.C:10
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
static const double kPi2
Definition: Constants.h:39
Root of GENIE utility namespaces.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
int PionNucleusXSec(double tpi, double ppistar, double t_new, double A, double &tpilow, double &siglow, double &tpihigh, double &sighigh)
Definition: HadXSUtils.cxx:214
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97
static const double kPionMass2
Definition: Constants.h:87
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353