BergerSehgalCOHPiPXSec2015.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 STFC, Rutherford Appleton Laboratory - March 11, 2005
9 
10 For the class documentation see the corresponding header file.
11 */
12 //____________________________________________________________________________
13 
14 #include <TMath.h>
15 
28 
29 using namespace genie;
30 using namespace genie::constants;
31 using namespace genie::utils;
32 
33 //____________________________________________________________________________
35  XSecAlgorithmI("genie::BergerSehgalCOHPiPXSec2015")
36 {
37 
38 }
39 //____________________________________________________________________________
41  XSecAlgorithmI("genie::BergerSehgalCOHPiPXSec2015", config)
42 {
43 
44 }
45 //____________________________________________________________________________
47 {
48 
49 }
50 //____________________________________________________________________________
52  const Interaction * interaction, KinePhaseSpace_t kps) const
53 {
54  // Here we are following PRD 79, 053003 (2009) by Berger and Sehgal
55  // This method computes the differential cross section represented
56  // in Eq.'s 6 (CC) and 7 (NC) from that paper.
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 x = kinematics.x();
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("BergerSehgalCohPi", pDEBUG) << "Pion COM momentum negative for Q2 = " << Q2 <<
77  " x = " << x << " y = " << y;
78  return 0.0;
79  }
80  double front = ExactKinematicTerm(interaction);
81  if (front <= 0.0) {
82  LOG("BergerSehgalCohPi", pDEBUG) << "Exact kin. form = " << front <<
83  " E = " << E << " Q2 = " << Q2 << " y = " << y << " x = " << x;
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 Epi = y*E; // ~pion energy
93  double ma2 = TMath::Power(fMa, 2); // "axial mass" squared
94  double Ga = ma2 / (ma2 + Q2);
95  double Ga2 = TMath::Power(Ga, 2.); // propagator term
96  double Ro2 = TMath::Power(fRo * units::fermi, 2.);
97 
98  // the xsec is d^3xsec/dQ^2dydt but the only t-dependent factor
99  // is an exp(-bt) so it can be integrated analyticaly
100  double Epi2 = TMath::Power(Epi, 2.);
101  double R = fRo * A_3 * units::fermi; // nuclear radius
102  double R2 = TMath::Power(R, 2.);
103  double b = 0.33333 * R2;
104  double MxEpi = M * x / Epi;
105  double mEpi2 = (M_pi * M_pi) / Epi2;
106  double tA = 1. + MxEpi - 0.5 * mEpi2;
107  double tB = TMath::Sqrt(1.0 + 2 * MxEpi) * TMath::Sqrt(1.0 - mEpi2);
108  double tmin = 2 * Epi2 * (tA - tB);
109  double tmax = 2 * Epi2 * (tA + tB);
110  if (tmin < 1.0e-8) {
111  tmin = 1.0e-8;
112  }
113 
114  /* const KPhaseSpace & kphase = interaction->PhaseSpace(); */
115  /* Range1D_t tl = kphase.TLim(); // TESTING! */
116 
117  double sigtot_pin = utils::hadxs::berger::PionNucleonXSec(Epi, /* get_total = */ true, pionIsCharged);
118  double sigel_pin = utils::hadxs::berger::PionNucleonXSec(Epi, /* get_total = */ false, pionIsCharged);
119  double siginel_pin = sigtot_pin - sigel_pin;
120 
121  // fabs (F_{abs}) describes the average attenuation of a pion emerging
122  // from a sphere of nuclear matter with radius = R_0 A^{1/3}. it is
123  // Eq. 13 in Berger-Sehgal PRD 79, 053003
124  double fabs_input = (9.0 * A_3) / (16.0 * kPi * Ro2);
125  double fabs = TMath::Exp( -1.0 * fabs_input * siginel_pin);
126 
127  // my old hackery to get things to work, A. Mislivec provided a better alt.
128  // double factor = 0.1; // to go from 10^-37 cm^2 -> 10^-38 cm^2
129  // double RS_factor = (units::mb*A2*fabs)/(16.0*kPi) * (sigtot_pin*sigtot_pin);
130 
131  // A_RS for BS version of RS, and/or Tpi>1.0
132  double RS_factor = (A2 * fabs) / (16.0 * kPi) * (sigtot_pin * sigtot_pin);
133 
134  // get the pion-nucleus cross section on carbon, fold it into differential cross section
135  double tpi = (E * y) - M_pi - ((Q2 + M_pi * M_pi) / (2 * M));
136  double tpilow = 0.0;
137  double siglow = 0.0;
138  double tpihigh = 0.0;
139  double sighigh = 0.0;
140  double dsigdzfit = 0.0;
141  double dsigdtfit = 0.0;
142  int xsec_stat = 0;
143  double dsig = 0.0;
144  double tstep = 100;
145  double logt_step = TMath::Abs(log(tmax) - log(tmin)) / tstep;
146  double logt = log(tmin) - logt_step/2.0;
147  double t_itt = TMath::Exp(logt);
148  double t_width = 0.0;
149 
150  for (double t_step = 0; t_step<tstep; t_step++) {
151 
152  logt = logt + logt_step;
153  t_itt = TMath::Exp(logt);
154  t_width = t_itt*logt_step;
155 
156  if (tpi <= 1.0 && fRSPionXSec == false) {
157  xsec_stat = utils::hadxs::berger::PionNucleusXSec(tpi, ppistar, t_itt, A, tpilow, siglow, tpihigh, sighigh);
158  if(xsec_stat){
159  LOG("BergerSehgalCohPi", pERROR) << "Call to PionNucleusXSec code failed - return xsec of 0.0";
160  return 0.0;
161  }
162  dsigdzfit = siglow + (sighigh - siglow) * (tpi - tpilow) / (tpihigh - tpilow);
163  dsigdtfit = dsigdzfit / (2.0 * ppistar * ppistar);
164  // we are handed a cross section in mb, need to convert it to GeV^{-2}
165  dsig += 1.0 * front * Ga2 * t_width * dsigdtfit * units::mb;
166  }
167  else {
168  dsig += /*factor **/ front * Ga2 * t_width * RS_factor * exp(-1.0*b*t_itt);
169  }
170 
171  }
172  xsec = dsig;
173 
174  // Correction for finite final state lepton mass.
175  // Lepton mass modification is part of Berger-Sehgal and is not optional.
176  if (pionIsCharged) {
177  double C = 1.;
178  // First, we need to remove the leading G_{A}^2 which is required for NC.
179  xsec /= Ga2;
180  // Next, \cos^2 \theta_{Cabibbo} appears in the CC xsec, but not the NC.
181  xsec *= fCos8c2;
182  double ml = interaction->FSPrimLepton()->Mass();
183  double ml2 = TMath::Power(ml,2);
184  double Q2min = ml2 * y/(1-y);
185  if(Q2 > Q2min) {
186  double C1 = TMath::Power(Ga - 0.5 * Q2min / (Q2 + kPionMass2), 2);
187  double C2 = 0.25 * y * Q2min * (Q2 - Q2min) /
188  TMath::Power(Q2 + kPionMass2, 2);
189  C = C1 + C2;
190  } else {
191  C = 0.;
192  }
193  xsec *= (2. * C); // *2 is for CC vs NC in BS
194  }
195 
196 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
197  LOG("BergerSehgalCohPi", pDEBUG)
198  << "\n momentum transfer .............. Q2 = " << Q2
199  << "\n mass number .................... A = " << A
200  << "\n pion energy .................... Epi = " << Epi
201  << "\n propagator term ................ Ga2 = " << Ga2
202  << "\n total pi+N cross section ....... sigT = " << sigtot_pin
203  << "\n inelastic pi+N cross section ... sigI = " << siginel_pin
204  << "\n nuclear size scale ............. Ro = " << fRo
205  << "\n pion absorption factor ......... Fabs = " << fabs
206  << "\n t integration range ............ [" << tmin << "," << tmax << "]"
207  LOG("BergerSehgalCohPi", pINFO)
208  << "d2xsec/dQ2dy[COHPi] (Q2= " << Q2 << ", y="
209  << y << ", E=" << E << ") = "<< xsec;
210 #endif
211 
212  //----- The algorithm computes d^2xsec/dQ2dy
213  // Check whether variable tranformation is needed? May be working with logs.
214  // kPSlogQ2logyfE is possible - all others will not succeed
215  if(kps != kPSQ2yfE) {
216  double J = utils::kinematics::Jacobian(interaction,kPSQ2yfE, kps);
217  xsec *= J;
218  }
219  return xsec;
220 }
221 //____________________________________________________________________________
223 {
224  // This function is a bit inefficient but is being encapsulated as
225  // such in order to possibly migrate into a general kinematics check.
226  const Kinematics & kinematics = interaction -> Kine();
227  const InitialState & init_state = interaction -> InitState();
228 
229  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
230  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
231  double E = init_state.ProbeE(kRfLab); // nu E
232  double Q2 = kinematics.Q2();
233  double y = kinematics.y(); // inelasticity
234  double fp2 = (0.93 * M_pi)*(0.93 * M_pi);
235 
236  double term = ((kGF2 * fp2) / (4.0 * kPi2)) *
237  ((E * (1.0 - y)) / sqrt(y*E * y*E + Q2)) *
238  (1.0 - Q2 / (4.0 * E*E * (1.0 - y)));
239  return term;
240 }
241 //____________________________________________________________________________
243 {
244  // This function is a bit inefficient but is being encapsulated as
245  // such in order to possibly migrate into a general kinematics check.
246  const Kinematics & kinematics = interaction -> Kine();
247  const InitialState & init_state = interaction -> InitState();
248 
249  bool pionIsCharged = interaction->ProcInfo().IsWeakCC();
250  double M_pi = pionIsCharged ? kPionMass : kPi0Mass;
251  double E = init_state.ProbeE(kRfLab); // nu E
252  double Q2 = kinematics.Q2();
253  double y = kinematics.y(); // inelasticity
254  double MT = init_state.Tgt().Mass();
255 
256  double W2 = MT*MT - Q2 + 2.0 * y * E * MT;
257  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);
258  if (arg < 0) return arg;
259  double ppistar = TMath::Sqrt(arg) / 2.0 / TMath::Sqrt(W2);
260 
261  return ppistar;
262 }
263 //____________________________________________________________________________
264 double BergerSehgalCOHPiPXSec2015::Integral(const Interaction * interaction) const
265 {
266  double xsec = fXSecIntegrator->Integrate(this,interaction);
267  return xsec;
268 }
269 //____________________________________________________________________________
271 {
272  if(interaction->TestBit(kISkipProcessChk)) return true;
273 
274  const InitialState & init_state = interaction->InitState();
275  const ProcessInfo & proc_info = interaction->ProcInfo();
276  const Target & target = init_state.Tgt();
277 
278  int nu = init_state.ProbePdg();
279 
280  if (!proc_info.IsCoherent()) return false;
281  if (!proc_info.IsWeak()) return false;
282  if (target.HitNucIsSet()) return false;
283  if (!(target.A()>1)) return false;
284  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
285 
286  return true;
287 }
288 //____________________________________________________________________________
290 {
291  Algorithm::Configure(config);
292  this->LoadConfig();
293 }
294 //____________________________________________________________________________
296 {
297  Algorithm::Configure(config);
298  this->LoadConfig();
299 }
300 //____________________________________________________________________________
302 {
303  GetParam( "COH-Ma",fMa ) ;
304  GetParam( "COH-Ro", fRo ) ;
305 
306  double thc ;
307  GetParam( "CabibboAngle", thc ) ;
308  fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
309 
310  // fRSPionXSec => Do not use the pion-nucleus cross section from Table 1 in PRD 79, 053003
311  // Instead, use the Rein-Sehgal "style" pion-nucleon cross section and scale by A
312  // for all pion energies.
313  GetParam( "COH-UseRSPionXSec", fRSPionXSec ) ;
314 
315  //-- load the differential cross section integrator
317  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
319 }
320 //____________________________________________________________________________
321 
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
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
Cross Section Integrator Interface.
bool fRSPionXSec
Use Rein-Sehgal "style" pion-nucleon xsecs.
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
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
Definition: HadXSUtils.cxx:105
double fRo
nuclear size scale parameter
double x(bool selected=false) const
Definition: Kinematics.cxx:109
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Mass(void) const
Definition: Target.cxx:241
double PionCOMAbsMomentum(const Interaction *i) const
double y(bool selected=false) const
Definition: Kinematics.cxx:122
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?
#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
#define R(x)
#define pINFO
Definition: Messenger.h:63
double ExactKinematicTerm(const Interaction *i) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Double_t xsec[nknots]
Definition: testXsec.C:47
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static const double A
Definition: Units.h:82
static const double kPionMass
Definition: Constants.h:74
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
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 ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
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
Float_t e
Definition: plot.C:35
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
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: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