RosenbluthPXSec.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
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Sep 19, 2009 - CA
14  First included in v2.5.1.
15  @ Nov 26, 2009 - CA
16  Fix mistake in convertion from dsigma/dOmega --> dsigma/dQ2 uncovered at
17  the first comparison against electron QE data.
18 
19 */
20 //____________________________________________________________________________
21 
22 #include <TMath.h>
23 
36 
37 using namespace genie;
38 using namespace genie::utils;
39 using namespace genie::constants;
40 
41 //____________________________________________________________________________
43 XSecAlgorithmI("genie::RosenbluthPXSec")
44 {
45 
46 }
47 //____________________________________________________________________________
49 XSecAlgorithmI("genie::RosenbluthPXSec", config)
50 {
51 
52 }
53 //____________________________________________________________________________
55 {
56  if (fCleanUpfElFFModel) {
57  delete fElFFModel;
58  }
59 }
60 //____________________________________________________________________________
62  const Interaction * interaction, KinePhaseSpace_t kps) const
63 {
64  if(! this -> ValidProcess (interaction) ) return 0.;
65  if(! this -> ValidKinematics (interaction) ) return 0.;
66 
67  // Get interaction information
68  const InitialState & init_state = interaction -> InitState();
69  const Kinematics & kinematics = interaction -> Kine();
70  const Target & target = init_state.Tgt();
71 
72  int nucpdgc = target.HitNucPdg();
73  double E = init_state.ProbeE(kRfHitNucRest);
74  double Q2 = kinematics.Q2();
75  double M = target.HitNucMass();
76 
77  double E2 = E*E;
78  double E3 = E*E2;
79  double M2 = M*M;
80 
81  // Calculate scattering angle
82  //
83  // Q^2 = 4 * E^2 * sin^2 (theta/2) / ( 1 + 2 * (E/M) * sin^2(theta/2) ) =>
84  // sin^2 (theta/2) = MQ^2 / (4ME^2 - 2EQ^2)
85 
86  double sin2_halftheta = M*Q2 / (4*M*E2 - 2*E*Q2);
87  double sin4_halftheta = TMath::Power(sin2_halftheta, 2.);
88  double cos2_halftheta = 1.-sin2_halftheta;
89  //unused double cos_halftheta = TMath::Sqrt(cos2_halftheta);
90  double tan2_halftheta = sin2_halftheta/cos2_halftheta;
91 
92  // Calculate the elastic nucleon form factors
93  fELFF.Calculate(interaction);
94  double Gm = pdg::IsProton(nucpdgc) ? fELFF.Gmp() : fELFF.Gmn();
95  double Ge = pdg::IsProton(nucpdgc) ? fELFF.Gep() : fELFF.Gen();
96  double Ge2 = Ge*Ge;
97  double Gm2 = Gm*Gm;
98 
99  // Calculate tau and the virtual photon polarization (epsilon)
100  double tau = Q2/(4*M2);
101  double epsilon = 1. / (1. + 2.*(1.+tau)*tan2_halftheta);
102 
103  // Calculate the scattered lepton energy
104  double Ep = E / (1. + 2.*(E/M)*sin2_halftheta);
105  //double Ep2 = Ep*Ep; // unused variable
106 
107  // Calculate the Mott cross section dsigma/dOmega
108  double xsec_mott = (0.25 * kAem2 * Ep / E3) * (cos2_halftheta/sin4_halftheta);
109 
110  // Calculate the electron-nucleon elastic cross section dsigma/dOmega
111  double xsec = xsec_mott * (Ge2 + (tau/epsilon)*Gm2) / (1+tau);
112 
113  // Convert dsigma/dOmega --> dsigma/dQ2
114  xsec *= (kPi/(Ep*E));
115 
116  // The algorithm computes dxsec/dQ2
117  // Check whether variable tranformation is needed
118  if(kps!=kPSQ2fE) {
119  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
120 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
121  LOG("Rosenbluth", pDEBUG)
122  << "Jacobian for transformation to: "
123  << KinePhaseSpace::AsString(kps) << ", J = " << J;
124 #endif
125  xsec *= J;
126  }
127 
128  // If requested, return the free nucleon xsec even for input nuclear tgt
129  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
130 
131  // Take into account the number of nucleons/tgt
132  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
133  xsec *= NNucl;
134 
135  // Compute & apply nuclear suppression factor
136  // (R(Q2) is adapted from NeuGEN - see comments therein)
137  double R = nuclear::NuclQELXSecSuppression("Default", 0.5, interaction);
138  xsec *= R;
139 
140  return xsec;
141 }
142 //____________________________________________________________________________
143 double RosenbluthPXSec::Integral(const Interaction * interaction) const
144 {
145  double xsec = fXSecIntegrator->Integrate(this,interaction);
146  return xsec;
147 }
148 //____________________________________________________________________________
149 bool RosenbluthPXSec::ValidProcess(const Interaction * interaction) const
150 {
151  if(interaction->TestBit(kISkipProcessChk)) return true;
152 
153  const ProcessInfo & proc_info = interaction->ProcInfo();
154 
155  if(!proc_info.IsQuasiElastic()) return false;
156  if(!proc_info.IsEM()) return false;
157 
158  const InitialState & init_state = interaction->InitState();
159 
160  int hitnuc = init_state.Tgt().HitNucPdg();
161  bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc));
162  if (!is_pn) return false;
163 
164  int probe = init_state.ProbePdg();
165  bool is_chgl = pdg::IsChargedLepton(probe);
166  if (!is_chgl) return false;
167 
168  return true;
169 }
170 //____________________________________________________________________________
172 {
173  Algorithm::Configure(config);
174  this->LoadConfig();
175 }
176 //____________________________________________________________________________
178 {
179  Algorithm::Configure(config);
180  this->LoadConfig();
181 }
182 //____________________________________________________________________________
184 {
185  fElFFModel = 0;
186 
187  // load elastic form factors model
188  fElFFModel = dynamic_cast<const ELFormFactorsModelI *> ( this -> SubAlg("ElasticFormFactorsModel" ) ) ;
189 
191 
192  fCleanUpfElFFModel = false;
193  bool useFFTE = false ;
194  GetParam( "UseElFFTransverseEnhancement", useFFTE ) ;
195  if( useFFTE ) {
196  const ELFormFactorsModelI* sub_alg = fElFFModel;
197  fElFFModel = dynamic_cast<const ELFormFactorsModelI *> ( this -> SubAlg("TransverseEnhancement") ) ;
198  dynamic_cast<const TransverseEnhancementFFModel*>(fElFFModel)->SetElFFBaseModel( sub_alg );
199  fCleanUpfElFFModel = true;
200  }
202 
203  // load XSec Integrator
205  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
207 }
208 //____________________________________________________________________________
Cross Section Calculation Interface.
const double kPi
Basic constants.
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
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.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int HitNucPdg(void) const
Definition: Target.cxx:321
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
double HitNucMass(void) const
Definition: Target.cxx:250
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:99
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Gen(void) const
Get the computed form factor Gen.
Definition: ELFormFactors.h:57
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
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?
const XSecIntegratorI * fXSecIntegrator
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Pure abstract base class. Defines the ELFormFactorsModelI interface to be implemented by any algorith...
Float_t E
Definition: plot.C:20
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
static const double kAem2
Definition: Constants.h:58
static string AsString(KinePhaseSpace_t kps)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
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)
int Z(void) const
Definition: Target.h:69
Double_t xsec[nknots]
Definition: testXsec.C:47
void Calculate(const Interaction *interaction)
Calculate the form factors for the input interaction using the attached algorithm.
bool IsEM(void) const
void Configure(const Registry &config)
double Gmn(void) const
Get the computed form factor Gmn.
Definition: ELFormFactors.h:60
int N(void) const
Definition: Target.h:70
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
double epsilon
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double Gmp(void) const
Get the computed form factor Gmp.
Definition: ELFormFactors.h:54
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
double Gep(void) const
Get the computed form factor Gep.
Definition: ELFormFactors.h:51
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
void SetModel(const ELFormFactorsModelI *model)
Attach an algorithm.
Modification of magnetic form factors to match observed enhancement in transverse cross section of th...
void kinematics()
Definition: kinematics.C:10
const ELFormFactorsModelI * fElFFModel
double ProbeE(RefFrame_t rf) const
double Integral(const Interaction *i) const
Root of GENIE utility namespaces.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353