GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ReinSehgalCOHPiPXSec.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
15 #include "Framework/Conventions/GBuild.h"
25 
26 using namespace genie;
27 using namespace genie::constants;
28 using namespace genie::utils;
29 
30 //____________________________________________________________________________
32 XSecAlgorithmI("genie::ReinSehgalCOHPiPXSec")
33 {
34 
35 }
36 //____________________________________________________________________________
38 XSecAlgorithmI("genie::ReinSehgalCOHPiPXSec", config)
39 {
40 
41 }
42 //____________________________________________________________________________
44 {
45 
46 }
47 //____________________________________________________________________________
49  const Interaction * interaction, KinePhaseSpace_t kps) const
50 {
51  if(! this -> ValidProcess (interaction) ) return 0.;
52  if(! this -> ValidKinematics (interaction) ) return 0.;
53 
54  const Kinematics & kinematics = interaction -> Kine();
55  const InitialState & init_state = interaction -> InitState();
56 
57  //----- Compute the coherent NC pi0 production d2xsec/dxdy
58  // see page 34 in Nucl.Phys.B223:29-144 (1983)
59  double E = init_state.ProbeE(kRfLab); // neutrino energy
60  double x = kinematics.x(); // bjorken x
61  double y = kinematics.y(); // inelasticity y
62  double Q2 = 2.*x*y*kNucleonMass*E; // momentum transfer Q2>0
63  double A = (double) init_state.Tgt().A(); // mass number
64  double A2 = TMath::Power(A,2.);
65  double A_3 = TMath::Power(A,1./3.);
66  double Gf = kGF2 * kNucleonMass / (32 * kPi3);
67  double fp = 0.93 * kPionMass; // pion decay constant
68  double fp2 = TMath::Power(fp,2.);
69  double Epi = y*E; // pion energy
70  double ma2 = TMath::Power(fMa,2);
71  double propg = TMath::Power(ma2/(ma2+Q2),2.); // propagator term
72  double r2 = TMath::Power(fReIm,2.);
73  double sTot = utils::hadxs::TotalPionNucleonXSec(Epi); // tot. pi+N xsec
74  double sTot2 = TMath::Power(sTot,2.);
75  double sInel = utils::hadxs::InelasticPionNucleonXSec(Epi); // inel. pi+N xsec
76  double Ro2 = TMath::Power(fRo*units::fermi,2.);
77 
78  // effect of pion absorption in the nucleus
79  double Fabs = TMath::Exp( -9.*A_3*sInel / (16.*kPi*Ro2) );
80 
81  // the xsec in Nucl.Phys.B223:29-144 (1983) is d^3xsec/dxdydt but the only
82  // t-dependent factor is an exp(-bt) so it can be integrated analyticaly
83  double Epi2 = TMath::Power(Epi,2.);
84  double R = fRo * A_3 * units::fermi; // nuclear radius
85  double R2 = TMath::Power(R,2.);
86  double b = 0.33333 * R2;
87  double MxEpi = kNucleonMass*x/Epi;
88  double mEpi2 = kPionMass2/Epi2;
89  double tA = 1. + MxEpi - 0.5*mEpi2;
90  double tB = TMath::Sqrt(1. + 2*MxEpi) * TMath::Sqrt(1.-mEpi2);
91  double tmin = 2*Epi2 * (tA-tB);
92  double tmax = 2*Epi2 * (tA+tB);
93  double tint = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b; // t integral
94 
95  double xsec = Gf*fp2 * A2 * E*(1-y) * sTot2 * (1+r2)*propg * Fabs*tint;
96 
97 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
98  LOG("ReinSehgalCohPi", pDEBUG)
99  << "\n momentum transfer .............. Q2 = " << Q2
100  << "\n mass number .................... A = " << A
101  << "\n pion energy .................... Epi = " << Epi
102  << "\n propagator term ................ propg = " << propg
103  << "\n Re/Im of fwd pion scat. ampl. .. Re/Im = " << fReIm
104  << "\n total pi+N cross section ....... sigT = " << sTot
105  << "\n inelastic pi+N cross section ... sigI = " << sInel
106  << "\n nuclear size scale ............. Ro = " << fRo
107  << "\n pion absorption factor ......... Fabs = " << Fabs
108  << "\n t integration range ............ [" << tmin << "," << tmax << "]"
109  << "\n t integration factor ........... tint = " << tint;
110 #endif
111 
112  // compute the cross section for the CC case
113 
114  if(interaction->ProcInfo().IsWeakCC()) {
115  // Check whether a modification to Adler's PCAC theorem is applied for
116  // including the effect of the muon mass.
117  // See Rein and Sehgal, PCAC and the Deficit of Forward Muons in pi+
118  // Production by Neutrinos, hep-ph/0606185
119  double C = 1.;
120  if(fModPCAC) {
121  double ml = interaction->FSPrimLepton()->Mass();
122  double ml2 = TMath::Power(ml,2);
123  double Q2min = ml2 * y/(1-y);
124  if(Q2>Q2min) {
125  double C1 = TMath::Power(1-0.5*Q2min/(Q2+kPionMass2), 2);
126  double C2 = 0.25*y*Q2min*(Q2-Q2min)/ TMath::Power(Q2+kPionMass2,2);
127  C = C1+C2;
128  } else {
129  C = 0.;
130  }
131  }
132  xsec *= (2.*C);
133  }
134 
135 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
136  LOG("ReinSehgalCohPi", pINFO)
137  << "d2xsec/dxdy[COHPi] (x= " << x << ", y="
138  << y << ", E=" << E << ") = "<< xsec;
139 #endif
140 
141  //----- The algorithm computes d^2xsec/dxdy
142  // Check whether variable tranformation is needed
143  if(kps!=kPSxyfE) {
144  double J = utils::kinematics::Jacobian(interaction,kPSxyfE,kps);
145  xsec *= J;
146  }
147 
148  return xsec;
149 }
150 //____________________________________________________________________________
151 double ReinSehgalCOHPiPXSec::Integral(const Interaction * interaction) const
152 {
153  double xsec = fXSecIntegrator->Integrate(this,interaction);
154  return xsec;
155 }
156 //____________________________________________________________________________
157 bool ReinSehgalCOHPiPXSec::ValidProcess(const Interaction * interaction) const
158 {
159  if(interaction->TestBit(kISkipProcessChk)) return true;
160 
161  const InitialState & init_state = interaction->InitState();
162  const ProcessInfo & proc_info = interaction->ProcInfo();
163  const Target & target = init_state.Tgt();
164 
165  int nu = init_state.ProbePdg();
166 
167  if (!proc_info.IsCoherentProduction()) return false;
168  if (!proc_info.IsWeak()) return false;
169  if (target.HitNucIsSet()) return false;
170  if (!(target.A()>1)) return false;
171  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
172 
173  return true;
174 }
175 //____________________________________________________________________________
177 {
178  Algorithm::Configure(config);
179  this->LoadConfig();
180 }
181 //____________________________________________________________________________
183 {
184  Algorithm::Configure(config);
185  this->LoadConfig();
186 }
187 //____________________________________________________________________________
189 {
190 
191  GetParam( "COH-Ma", fMa ) ;
192  GetParam( "COH-Ro", fRo ) ;
193  GetParam( "COH-ReImAmpl", fReIm ) ;
194  GetParam( "COH-UseModifiedPCAC", fModPCAC ) ;
195 
196  //-- load the differential cross section integrator
198  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
199  assert(fXSecIntegrator);
200 }
201 //____________________________________________________________________________
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
Cross Section Integrator Interface.
static const double kNucleonMass
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 x(bool selected=false) const
Definition: Kinematics.cxx:99
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
bool IsCoherentProduction(void) const
enum genie::EKinePhaseSpace KinePhaseSpace_t
bool fModPCAC
use modified PCAC (including f/s lepton mass)
double fReIm
Re/Im {forward pion scattering amplitude}.
double y(bool selected=false) const
Definition: Kinematics.cxx:112
static constexpr double b
Definition: Units.h:78
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
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?
double Integral(const Interaction *i) const
#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
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
const XSecIntegratorI * fXSecIntegrator
#define pINFO
Definition: Messenger.h:62
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:20
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
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 fRo
nuclear size scale parameter
void Configure(const Registry &config)
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 TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:51
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345