ReinSehgalCOHPiPXSec.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  University of Liverpool & STFC Rutherford Appleton Lab - March 11, 2005
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Nov 21, 2007 - CA
14  Was renamed to ReinSehgalCOHPiPXSec (from ReinSehgalCOHPXSec)
15  @ Mar 31, 2009 - CA
16  Fixed a minor bug in the C2 term controlling the forward mu- suppression
17  predicted by Adler's PCAC (mpi -> mpi^2)
18 */
19 //____________________________________________________________________________
20 
21 #include <TMath.h>
22 
35 
36 using namespace genie;
37 using namespace genie::constants;
38 using namespace genie::utils;
39 
40 //____________________________________________________________________________
42 XSecAlgorithmI("genie::ReinSehgalCOHPiPXSec")
43 {
44 
45 }
46 //____________________________________________________________________________
48 XSecAlgorithmI("genie::ReinSehgalCOHPiPXSec", config)
49 {
50 
51 }
52 //____________________________________________________________________________
54 {
55 
56 }
57 //____________________________________________________________________________
59  const Interaction * interaction, KinePhaseSpace_t kps) const
60 {
61  if(! this -> ValidProcess (interaction) ) return 0.;
62  if(! this -> ValidKinematics (interaction) ) return 0.;
63 
64  const Kinematics & kinematics = interaction -> Kine();
65  const InitialState & init_state = interaction -> InitState();
66 
67  //----- Compute the coherent NC pi0 production d2xsec/dxdy
68  // see page 34 in Nucl.Phys.B223:29-144 (1983)
69  double E = init_state.ProbeE(kRfLab); // neutrino energy
70  double x = kinematics.x(); // bjorken x
71  double y = kinematics.y(); // inelasticity y
72  double Q2 = 2.*x*y*kNucleonMass*E; // momentum transfer Q2>0
73  double A = (double) init_state.Tgt().A(); // mass number
74  double A2 = TMath::Power(A,2.);
75  double A_3 = TMath::Power(A,1./3.);
76  double Gf = kGF2 * kNucleonMass / (32 * kPi3);
77  double fp = 0.93 * kPionMass; // pion decay constant
78  double fp2 = TMath::Power(fp,2.);
79  double Epi = y*E; // pion energy
80  double ma2 = TMath::Power(fMa,2);
81  double propg = TMath::Power(ma2/(ma2+Q2),2.); // propagator term
82  double r2 = TMath::Power(fReIm,2.);
83  double sTot = utils::hadxs::TotalPionNucleonXSec(Epi); // tot. pi+N xsec
84  double sTot2 = TMath::Power(sTot,2.);
85  double sInel = utils::hadxs::InelasticPionNucleonXSec(Epi); // inel. pi+N xsec
86  double Ro2 = TMath::Power(fRo*units::fermi,2.);
87 
88  // effect of pion absorption in the nucleus
89  double Fabs = TMath::Exp( -9.*A_3*sInel / (16.*kPi*Ro2) );
90 
91  // the xsec in Nucl.Phys.B223:29-144 (1983) is d^3xsec/dxdydt but the only
92  // t-dependent factor is an exp(-bt) so it can be integrated analyticaly
93  double Epi2 = TMath::Power(Epi,2.);
94  double R = fRo * A_3 * units::fermi; // nuclear radius
95  double R2 = TMath::Power(R,2.);
96  double b = 0.33333 * R2;
97  double MxEpi = kNucleonMass*x/Epi;
98  double mEpi2 = kPionMass2/Epi2;
99  double tA = 1. + MxEpi - 0.5*mEpi2;
100  double tB = TMath::Sqrt(1. + 2*MxEpi) * TMath::Sqrt(1.-mEpi2);
101  double tmin = 2*Epi2 * (tA-tB);
102  double tmax = 2*Epi2 * (tA+tB);
103  double tint = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b; // t integral
104 
105  double xsec = Gf*fp2 * A2 * E*(1-y) * sTot2 * (1+r2)*propg * Fabs*tint;
106 
107 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
108  LOG("ReinSehgalCohPi", pDEBUG)
109  << "\n momentum transfer .............. Q2 = " << Q2
110  << "\n mass number .................... A = " << A
111  << "\n pion energy .................... Epi = " << Epi
112  << "\n propagator term ................ propg = " << propg
113  << "\n Re/Im of fwd pion scat. ampl. .. Re/Im = " << fReIm
114  << "\n total pi+N cross section ....... sigT = " << sTot
115  << "\n inelastic pi+N cross section ... sigI = " << sInel
116  << "\n nuclear size scale ............. Ro = " << fRo
117  << "\n pion absorption factor ......... Fabs = " << Fabs
118  << "\n t integration range ............ [" << tmin << "," << tmax << "]"
119  << "\n t integration factor ........... tint = " << tint;
120 #endif
121 
122  // compute the cross section for the CC case
123 
124  if(interaction->ProcInfo().IsWeakCC()) {
125  // Check whether a modification to Adler's PCAC theorem is applied for
126  // including the effect of the muon mass.
127  // See Rein and Sehgal, PCAC and the Deficit of Forward Muons in pi+
128  // Production by Neutrinos, hep-ph/0606185
129  double C = 1.;
130  if(fModPCAC) {
131  double ml = interaction->FSPrimLepton()->Mass();
132  double ml2 = TMath::Power(ml,2);
133  double Q2min = ml2 * y/(1-y);
134  if(Q2>Q2min) {
135  double C1 = TMath::Power(1-0.5*Q2min/(Q2+kPionMass2), 2);
136  double C2 = 0.25*y*Q2min*(Q2-Q2min)/ TMath::Power(Q2+kPionMass2,2);
137  C = C1+C2;
138  } else {
139  C = 0.;
140  }
141  }
142  xsec *= (2.*C);
143  }
144 
145 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
146  LOG("ReinSehgalCohPi", pINFO)
147  << "d2xsec/dxdy[COHPi] (x= " << x << ", y="
148  << y << ", E=" << E << ") = "<< xsec;
149 #endif
150 
151  //----- The algorithm computes d^2xsec/dxdy
152  // Check whether variable tranformation is needed
153  if(kps!=kPSxyfE) {
154  double J = utils::kinematics::Jacobian(interaction,kPSxyfE,kps);
155  xsec *= J;
156  }
157 
158  return xsec;
159 }
160 //____________________________________________________________________________
161 double ReinSehgalCOHPiPXSec::Integral(const Interaction * interaction) const
162 {
163  double xsec = fXSecIntegrator->Integrate(this,interaction);
164  return xsec;
165 }
166 //____________________________________________________________________________
167 bool ReinSehgalCOHPiPXSec::ValidProcess(const Interaction * interaction) const
168 {
169  if(interaction->TestBit(kISkipProcessChk)) return true;
170 
171  const InitialState & init_state = interaction->InitState();
172  const ProcessInfo & proc_info = interaction->ProcInfo();
173  const Target & target = init_state.Tgt();
174 
175  int nu = init_state.ProbePdg();
176 
177  if (!proc_info.IsCoherent()) return false;
178  if (!proc_info.IsWeak()) return false;
179  if (target.HitNucIsSet()) return false;
180  if (!(target.A()>1)) return false;
181  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
182 
183  return true;
184 }
185 //____________________________________________________________________________
187 {
188  Algorithm::Configure(config);
189  this->LoadConfig();
190 }
191 //____________________________________________________________________________
193 {
194  Algorithm::Configure(config);
195  this->LoadConfig();
196 }
197 //____________________________________________________________________________
199 {
200 
201  GetParam( "COH-Ma", fMa ) ;
202  GetParam( "COH-Ro", fRo ) ;
203  GetParam( "COH-ReImAmpl", fReIm ) ;
204  GetParam( "COH-UseModifiedPCAC", fModPCAC ) ;
205 
206  //-- load the differential cross section integrator
208  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
210 }
211 //____________________________________________________________________________
212 
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
Cross Section Integrator Interface.
static const double kNucleonMass
Definition: Constants.h:78
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
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double x(bool selected=false) const
Definition: Kinematics.cxx:109
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: config.py:1
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:122
const double C
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
static const double kPi3
Definition: Constants.h:40
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: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)
const XSecIntegratorI * fXSecIntegrator
#define pINFO
Definition: Messenger.h:63
Double_t xsec[nknots]
Definition: testXsec.C:47
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:29
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
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
FILE * fp
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
void kinematics()
Definition: kinematics.C:10
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...
Root of GENIE utility namespaces.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:60
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
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