QELUtils.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: Steven Gardiner <gardiner \at fnal.gov>
8  Fermi National Accelerator Laboratory
9 */
10 //____________________________________________________________________________
11 
12 // standard library includes
13 #include <cmath>
14 #include <limits>
15 #include <sstream>
16 
17 // ROOT includes
18 #include "TDatabasePDG.h"
19 #include "TLorentzVector.h"
20 
21 // GENIE includes
31 
32 
33 // Helper functions local to this source file
34 namespace {
35 
36  TVector3 COMframe2Lab(const genie::InitialState& initialState)
37  {
38  TLorentzVector* k4 = initialState.GetProbeP4( genie::kRfLab );
39  TLorentzVector* p4 = initialState.TgtPtr()->HitNucP4Ptr();
40  TLorentzVector totMom = *k4 + *p4;
41 
42  TVector3 beta = totMom.BoostVector();
43 
44  delete k4;
45 
46  return beta;
47  }
48 
49 }
50 
52  const genie::Interaction& inter)
53 {
54  // Get the final-state lepton and struck nucleon 4-momenta in the lab frame
55  const TLorentzVector& lepton = inter.Kine().FSLeptonP4();
56  const TLorentzVector& outNucleon = inter.Kine().HadSystP4();
57 
58  // Get beta for a Lorentz boost from the lab frame to the COM frame and
59  // vice-versa
60  TLorentzVector* probe = inter.InitStatePtr()->GetProbeP4( kRfLab );
61  const TLorentzVector& hit_nucleon = inter.InitStatePtr()->TgtPtr()
62  ->HitNucP4();
63  TLorentzVector total_p4 = (*probe) + hit_nucleon;
64  TVector3 beta_COM_to_lab = total_p4.BoostVector();
65  TVector3 beta_lab_to_COM = -beta_COM_to_lab;
66 
67  // Square of the Lorentz gamma factor for the boosts
68  double gamma2 = std::pow(total_p4.Gamma(), 2);
69 
70  // Get the final-state lepton 4-momentum in the COM frame
71  TLorentzVector leptonCOM = TLorentzVector( lepton );
72  leptonCOM.Boost( beta_lab_to_COM );
73 
74  // Angle between COM outgoing lepton and lab-frame velocity of COM frame
75  double theta0 = leptonCOM.Angle( beta_COM_to_lab );
76  double cos_theta0_squared = std::pow(std::cos( theta0 ), 2);
77  //double sin_theta0_squared = 1. - cos_theta0_squared; // sin^2(x) + cos^2(x) = 1
78 
79  // Difference in lab frame velocities between lepton and outgoing nucleon
80  TVector3 lepton_vel = ( lepton.Vect() * (1. / lepton.E()) );
81  TVector3 outNucleon_vel = ( outNucleon.Vect() * (1. / outNucleon.E()) );
82  double vel_diff = ( lepton_vel - outNucleon_vel ).Mag();
83 
84  // Compute the factor in the kPSQELEvGen phase space that arises from
85  // eliminating the energy-conserving delta function
86  double factor = std::sqrt( 1. + (1. - cos_theta0_squared)*(gamma2 - 1.) )
87  * std::pow(leptonCOM.P(), 2) / vel_diff;
88 
89  delete probe;
90 
91  return factor;
92 }
93 
95  const genie::NuclearModelI* nucl_model, const genie::XSecAlgorithmI* xsec_model,
96  double cos_theta_0, double phi_0, double& Eb,
97  genie::QELEvGen_BindingMode_t hitNucleonBindingMode, double min_angle_EM,
98  bool bind_nucleon)
99 {
100  // If requested, set the initial hit nucleon 4-momentum to be off-shell
101  // according to the binding mode specified in the function call
102  if ( bind_nucleon ) {
103  genie::utils::BindHitNucleon(*interaction, *nucl_model, Eb,
104  hitNucleonBindingMode);
105  }
106 
107  // Mass of the outgoing lepton
108  double lepMass = interaction->FSPrimLepton()->Mass();
109 
110  // Look up the (on-shell) mass of the final nucleon
111  TDatabasePDG *tb = TDatabasePDG::Instance();
112  double mNf = tb->GetParticle( interaction->RecoilNucleonPdg() )->Mass();
113 
114  // Mandelstam s for the probe/hit nucleon system
115  double s = std::pow( interaction->InitState().CMEnergy(), 2 );
116 
117  // Return a differential cross section of zero if we're below threshold (and
118  // therefore need to sample a new event)
119  if ( std::sqrt(s) < lepMass + mNf ) return 0.;
120 
121  double outLeptonEnergy = ( s - mNf*mNf + lepMass*lepMass ) / (2 * std::sqrt(s));
122 
123  if (outLeptonEnergy*outLeptonEnergy - lepMass*lepMass < 0.) return 0.;
124  double outMomentum = TMath::Sqrt(outLeptonEnergy*outLeptonEnergy - lepMass*lepMass);
125 
126  // Compute the boost vector for moving from the COM frame to the
127  // lab frame, i.e., the velocity of the COM frame as measured
128  // in the lab frame.
129  TVector3 beta = COMframe2Lab( interaction->InitState() );
130 
131  // FullDifferentialXSec depends on theta_0 and phi_0, the lepton COM
132  // frame angles with respect to the direction of the COM frame velocity
133  // as measured in the lab frame. To generate the correct dependence
134  // here, first set the lepton COM frame angles with respect to +z
135  // (via TVector3::SetTheta() and TVector3::SetPhi()).
136  TVector3 lepton3Mom(0., 0., outMomentum);
137  lepton3Mom.SetTheta( TMath::ACos(cos_theta_0) );
138  lepton3Mom.SetPhi( phi_0 );
139 
140  // Then rotate the lepton 3-momentum so that the old +z direction now
141  // points along the COM frame velocity (beta)
142  TVector3 zvec(0., 0., 1.);
143  TVector3 rot = ( zvec.Cross(beta) ).Unit();
144  double angle = beta.Angle( zvec );
145 
146  // Handle the edge case where beta is along -z, so the
147  // cross product above vanishes
148  if ( beta.Perp() == 0. && beta.Z() < 0. ) {
149  rot = TVector3(0., 1., 0.);
150  angle = genie::constants::kPi;
151  }
152 
153  // Rotate if the rotation vector is not 0
154  if ( rot.Mag() >= genie::controls::kASmallNum ) {
155  lepton3Mom.Rotate(angle, rot);
156  }
157 
158  // Construct the lepton 4-momentum in the COM frame
159  TLorentzVector lepton(lepton3Mom, outLeptonEnergy);
160 
161  // The final state nucleon will have an equal and opposite 3-momentum
162  // in the COM frame and will be on the mass shell
163  TLorentzVector outNucleon(-1*lepton.Px(),-1*lepton.Py(),-1*lepton.Pz(),
164  TMath::Sqrt(outMomentum*outMomentum + mNf*mNf));
165 
166  // Boost the 4-momenta for both particles into the lab frame
167  lepton.Boost(beta);
168  outNucleon.Boost(beta);
169 
170  // Check if event is at a low angle - if so return 0 and stop wasting time
171  if (180 * lepton.Theta() / genie::constants::kPi < min_angle_EM && interaction->ProcInfo().IsEM()) {
172  return 0;
173  }
174 
175  TLorentzVector * nuP4 = interaction->InitState().GetProbeP4( genie::kRfLab );
176  TLorentzVector qP4 = *nuP4 - lepton;
177  delete nuP4;
178  double Q2 = -1 * qP4.Mag2();
179 
180  interaction->KinePtr()->SetFSLeptonP4( lepton );
181  interaction->KinePtr()->SetHadSystP4( outNucleon );
182  interaction->KinePtr()->SetQ2( Q2 );
183 
184  // Check the Q2 range. If we're outside of it, don't bother
185  // with the rest of the calculation.
186  Range1D_t Q2lim = interaction->PhaseSpace().Q2Lim();
187  if (Q2 < Q2lim.min || Q2 > Q2lim.max) return 0.;
188 
189  // Compute the QE cross section for the current kinematics
190  double xsec = xsec_model->XSec(interaction, genie::kPSQELEvGen);
191 
192  return xsec;
193 }
194 
196  const std::string& binding_mode)
197 {
198  // Translate the string setting the binding mode to the appropriate
199  // enum value, or complain if one couldn't be found
200  if ( binding_mode == "UseGroundStateRemnant" ) {
201  return kUseGroundStateRemnant;
202  }
203  else if ( binding_mode == "UseNuclearModel" ) {
204  return kUseNuclearModel;
205  }
206  else if ( binding_mode == "OnShell" ) {
207  return kOnShell;
208  }
209  else {
210  LOG("QELEvent", pFATAL) << "Unrecognized setting \"" << binding_mode
211  << "\" requested in genie::utils::StringToQELBindingMode()";
212  gAbortingInErr = true;
213  std::exit(1);
214  }
215 }
216 
217 double genie::utils::CosTheta0Max(const genie::Interaction& interaction) {
218 
219  // q0 > 0 only needs to be enforced (indirectly via a calculation of
220  // CosTheta0Max) for bound nucleons. The Q2 limits should take care of valid
221  // kinematics for free nucleons.
222  if ( !interaction.InitState().Tgt().IsNucleus()
223  || interaction.TestBit(kIAssumeFreeNucleon) ) return 1.;
224 
225  double probe_E_lab = interaction.InitState().ProbeE( genie::kRfLab );
226 
227  TVector3 beta = COMframe2Lab( interaction.InitState() );
228  double gamma = 1. / std::sqrt(1. - beta.Mag2());
229 
230  double sqrt_s = interaction.InitState().CMEnergy();
231  double mNf = interaction.RecoilNucleon()->Mass();
232  double ml = interaction.FSPrimLepton()->Mass();
233  double lepton_E_COM = (sqrt_s*sqrt_s + ml*ml - mNf*mNf) / (2.*sqrt_s);
234 
235  // If there isn't enough available energy to create an on-shell
236  // final lepton, then don't bother with the rest of the calculation
237  // NOTE: C++11 would allow use to use lowest() here instead
238  if ( lepton_E_COM <= ml ) return -std::numeric_limits<double>::max();
239 
240  double lepton_p_COM = std::sqrt( std::max(lepton_E_COM*lepton_E_COM
241  - ml*ml, 0.) );
242 
243  // Possibly off-shell initial struck nucleon total energy
244  // (BindHitNucleon() should have been called previously if needed)
245  const TLorentzVector& p4Ni = interaction.InitState().Tgt().HitNucP4();
246  double ENi = p4Ni.E();
247  // On-shell mass of initial struck nucleon
248  double mNi = interaction.InitState().Tgt().HitNucMass();
249  // On-shell initial struck nucleon energy
250  double ENi_on_shell = std::sqrt( mNi*mNi + p4Ni.Vect().Mag2() );
251  // Energy needed to put initial nucleon on the mass shell
252  double epsilon_B = ENi_on_shell - ENi;
253 
254  double cos_theta0_max = ( probe_E_lab - gamma*lepton_E_COM - epsilon_B )
255  / ( gamma * lepton_p_COM * beta.Mag() );
256  return cos_theta0_max;
257 }
258 
260  const genie::NuclearModelI& nucl_model, double& Eb,
261  genie::QELEvGen_BindingMode_t hitNucleonBindingMode)
262 {
263  genie::Target* tgt = interaction.InitState().TgtPtr();
264  TLorentzVector* p4Ni = tgt->HitNucP4Ptr();
265 
266  // Initial nucleon 3-momentum (lab frame)
267  TVector3 p3Ni = nucl_model.Momentum3();
268 
269  // Look up the (on-shell) mass of the initial nucleon
270  TDatabasePDG* tb = TDatabasePDG::Instance();
271  double mNi = tb->GetParticle( tgt->HitNucPdg() )->Mass();
272 
273  // Set the (possibly off-shell) initial nucleon energy based on
274  // the selected binding energy mode. Always put the initial nucleon
275  // on shell if it is not part of a composite nucleus
276  double ENi = 0.;
277  if ( tgt->IsNucleus() && hitNucleonBindingMode != genie::kOnShell ) {
278 
279  // For a nuclear target with a bound initial struck nucleon, take binding
280  // energy effects and Pauli blocking into account when computing QE
281  // differential cross sections
282  interaction.ResetBit( kIAssumeFreeNucleon );
283 
284  // Initial nucleus mass
285  double Mi = tgt->Mass();
286 
287  // Final nucleus mass
288  double Mf = 0.;
289 
290  // If we use the removal energy reported by the nuclear
291  // model, then it implies a certain value for the final
292  // nucleus mass
293  if ( hitNucleonBindingMode == genie::kUseNuclearModel ) {
294  Eb = nucl_model.RemovalEnergy();
295  // This equation is the definition that we assume
296  // here for the "removal energy" (Eb) returned by the
297  // nuclear model. It matches GENIE's convention for
298  // the Bodek/Ritchie Fermi gas model.
299  Mf = Mi + Eb - mNi;
300  }
301  // We can also assume that the final nucleus is in its
302  // ground state. In this case, we can just look up its
303  // mass from the standard table. This implies a particular
304  // binding energy for the hit nucleon.
305  else if ( hitNucleonBindingMode == genie::kUseGroundStateRemnant ) {
306  // Determine the mass and proton numbers for the remnant nucleus
307  int Af = tgt->A() - 1;
308  int Zf = tgt->Z();
309  if ( genie::pdg::IsProton( tgt->HitNucPdg()) ) --Zf;
310  Mf = genie::PDGLibrary::Instance()->Find( genie::pdg::IonPdgCode(Af, Zf) )->Mass();
311 
312  // Deduce the binding energy from the final nucleus mass
313  Eb = Mf - Mi + mNi;
314  }
315 
316  // The (lab-frame) off-shell initial nucleon energy is the difference
317  // between the lab frame total energies of the initial and remnant nuclei
318  ENi = Mi - std::sqrt( Mf*Mf + p3Ni.Mag2() );
319  }
320  else {
321  // Keep the struck nucleon on shell either because
322  // hitNucleonBindingMode == kOnShell or because
323  // the target is a single nucleon
324  ENi = std::sqrt( p3Ni.Mag2() + std::pow(mNi, 2) );
325  Eb = 0.;
326 
327  // If we're dealing with a nuclear target but using the on-shell
328  // binding mode, set the "assume free nucleon" flag. This turns off
329  // Pauli blocking and the requirement that q0 > 0 in the QE cross section
330  // models (an on-shell nucleon *is* a free nucleon)
331  if ( tgt->IsNucleus() ) interaction.SetBit( kIAssumeFreeNucleon );
332  }
333 
334  // Update the initial nucleon lab-frame 4-momentum in the interaction with
335  // its current components
336  p4Ni->SetVect( p3Ni );
337  p4Ni->SetE( ENi );
338 
339 }
Cross Section Calculation Interface.
T max(const caf::Proxy< T > &a, T b)
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Double_t angle
Definition: plot.C:86
double RemovalEnergy(void) const
Definition: NuclearModelI.h:57
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int HitNucPdg(void) const
Definition: Target.cxx:321
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
int RecoilNucleonPdg(void) const
recoil nucleon pdg
int A(void) const
Definition: Target.h:71
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double HitNucMass(void) const
Definition: Target.cxx:250
T sqrt(T number)
Definition: d0nt_math.hpp:156
#define pFATAL
Definition: Messenger.h:57
bool IsNucleus(void) const
Definition: Target.cxx:289
constexpr T pow(T x)
Definition: pow.h:72
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:67
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
double Mass(Resonance_t res)
resonance mass (GeV)
Double_t beta
double ComputeFullQELPXSec(Interaction *interaction, const NuclearModelI *nucl_model, const XSecAlgorithmI *xsec_model, double cos_theta_0, double phi_0, double &Eb, QELEvGen_BindingMode_t hitNucleonBindingMode, double min_angle_EM=0., bool bind_nucleon=true)
Definition: QELUtils.cxx:94
Range1D_t Q2Lim(void) const
Q2 limits.
void BindHitNucleon(Interaction &interaction, const NuclearModelI &nucl_model, double &Eb, QELEvGen_BindingMode_t hitNucleonBindingMode)
Definition: QELUtils.cxx:259
double Mass(void) const
Definition: Target.cxx:241
enum genie::EQELEvGenBindingMode QELEvGen_BindingMode_t
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
const XML_Char * s
Definition: expat.h:262
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:92
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:66
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
Definition: QELUtils.cxx:51
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
int Z(void) const
Definition: Target.h:69
static const double kASmallNum
Definition: Controls.h:40
TVector3 Unit() const
Double_t xsec[nknots]
Definition: testXsec.C:47
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
double max
Definition: Range1.h:54
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double CosTheta0Max(const genie::Interaction &interaction)
Definition: QELUtils.cxx:217
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:195
double CMEnergy() const
centre-of-mass energy (sqrt s)
Target * TgtPtr(void) const
Definition: InitialState.h:68
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
float Mag() const
exit(0)
T cos(T number)
Definition: d0nt_math.hpp:78
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:317
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
const Target & Tgt(void) const
Definition: InitialState.h:67
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
bool gAbortingInErr
Definition: Messenger.cxx:56
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:38
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
Initial State information.
Definition: InitialState.h:49
enum BeamMode string