genie::utils Namespace Reference

Root of GENIE utility namespaces. More...

Namespaces
	app_init
	Initialization code commonly occuring in GENIE apps, factored out from existing apps for convenience. Not generic GENIE initialization code.
	bwfunc
	Breit Wigner functions.
	config
	Simple functions for loading and reading nucleus dependent keys from config files.
	frgmfunc
	Fragmentation functions.
	geometry
	Geometry utilities.
	ghep
	GHEP event record utilities.
	gsl
	Simple utilities for integrating GSL in the GENIE framework.
	gui
	Simple utilities for GENIE Graphical User Interface widgets.
	hadxs
	Simple functions and data for computing hadron interaction xsecs.
	hnl
	Useful kinematic functions.
	intranuke
	intranuke2018
	kinematics
	Kinematical utilities.
	math
	Simple mathematical utilities not found in ROOT's TMath.
	mec
	MEC utilities.
	nnbar_osc
	nuclear
	Simple nuclear physics empirical formulas (densities, radii, ...) and empirical nuclear corrections.
	nucleon_decay
	phys
	Various physics formulas & utilities.
	prem
	Preliminary Earth Model.
	print
	Simple printing utilities.
	res
	Baryon Resonance utilities.
	str
	Utilities for string manipulation.
	style
	GENIE style!
	system
	System utilities.
	units
	Simple unit system utilities.
	xml

Classes
class	T2KEvGenMetaData
	Utility class to store MC job meta-data. More...
class	hnl
	Utilities for simulating the decay of Heavy Neutral Leptons. More...
class	neutron_osc
	Utilities for simulating neutron oscillation. More...
class	nucleon_decay
	Utilities for simulating nucleon decay. More...

Functions
ostream &	operator<< (ostream &stream, const T2KEvGenMetaData &md)
double	EnergyDeltaFunctionSolutionDMEL (const Interaction &inter)
DMELEvGen_BindingMode_t	StringToDMELBindingMode (const std::string &mode_str)
double	ComputeFullDMELPXSec (Interaction *interaction, const NuclearModelI *nucl_model, const XSecAlgorithmI *xsec_model, double cos_theta_0, double phi_0, double &Eb, DMELEvGen_BindingMode_t hitNucleonBindingMode, double min_angle_EM=0., bool bind_nucleon=true)
double	CosTheta0Max (const genie::Interaction &interaction)
void	BindHitNucleon (Interaction &interaction, const NuclearModelI &nucl_model, double &Eb, DMELEvGen_BindingMode_t hitNucleonBindingMode)
void	SetPrimaryLeptonPolarization (GHepRecord *ev)
double	EnergyDeltaFunctionSolutionQEL (const Interaction &inter)
QELEvGen_BindingMode_t	StringToQELBindingMode (const std::string &mode_str)
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)
void	BindHitNucleon (Interaction &interaction, const NuclearModelI &nucl_model, double &Eb, QELEvGen_BindingMode_t hitNucleonBindingMode)

Detailed Description

Root of GENIE utility namespaces.

Common functions used for handling generation of the primary lepton, regardless of whether the relevant class inherits from PrimaryLeptonGenerator or not.

Author

Steven Gardiner <gardiner fnal.gov> Fermi National Accelerator Laboratory

Created:

May 01, 2020

License:

Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Function Documentation

void genie::utils::BindHitNucleon	(	genie::Interaction &	interaction,
		const NuclearModelI &	nucl_model,
		double &	Eb,
		genie::DMELEvGen_BindingMode_t	hitNucleonBindingMode
	)

Definition at line 259 of file DMELUtils.cxx.

References genie::Target::A(), genie::PDGLibrary::Find(), genie::Target::HitNucP4Ptr(), genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCode(), genie::Target::IsNucleus(), genie::pdg::IsProton(), genie::kIAssumeFreeNucleon, genie::kOnShell, genie::kUseGroundStateRemnant, genie::kUseNuclearModel, genie::utils::res::Mass(), genie::Target::Mass(), genie::NuclearModelI::Momentum3(), genie::NuclearModelI::RemovalEnergy(), genie::InitialState::TgtPtr(), and genie::Target::Z().

Referenced by ComputeFullDMELPXSec(), ComputeFullQELPXSec(), genie::QELEventGenerator::ComputeMaxXSec(), genie::DMELEventGenerator::ComputeMaxXSec(), genie::QELEventGenerator::ProcessEventRecord(), and genie::DMELEventGenerator::ProcessEventRecord().

{
  genie::Target* tgt = interaction.InitState().TgtPtr();
  TLorentzVector* p4Ni = tgt->HitNucP4Ptr();

  // Initial nucleon 3-momentum (lab frame)
  TVector3 p3Ni = nucl_model.Momentum3();

  // Look up the (on-shell) mass of the initial nucleon
  TDatabasePDG* tb = TDatabasePDG::Instance();
  double mNi = tb->GetParticle( tgt->HitNucPdg() )->Mass();

  // Set the (possibly off-shell) initial nucleon energy based on
  // the selected binding energy mode. Always put the initial nucleon
  // on shell if it is not part of a composite nucleus
  double ENi = 0.;
  if ( tgt->IsNucleus() && hitNucleonBindingMode != genie::kOnShell ) {

    // For a nuclear target with a bound initial struck nucleon, take binding
    // energy effects and Pauli blocking into account when computing QE
    // differential cross sections
    interaction.ResetBit( kIAssumeFreeNucleon );

    // Initial nucleus mass
    double Mi = tgt->Mass();

    // Final nucleus mass
    double Mf = 0.;

    // If we use the removal energy reported by the nuclear
    // model, then it implies a certain value for the final
    // nucleus mass
    if ( hitNucleonBindingMode == genie::kUseNuclearModel ) {
      Eb = nucl_model.RemovalEnergy();
      // This equation is the definition that we assume
      // here for the "removal energy" (Eb) returned by the
      // nuclear model. It matches GENIE's convention for
      // the Bodek/Ritchie Fermi gas model.
      Mf = Mi + Eb - mNi;
    }
    // We can also assume that the final nucleus is in its
    // ground state. In this case, we can just look up its
    // mass from the standard table. This implies a particular
    // binding energy for the hit nucleon.
    else if ( hitNucleonBindingMode == genie::kUseGroundStateRemnant ) {
      // Determine the mass and proton numbers for the remnant nucleus
      int Af = tgt->A() - 1;
      int Zf = tgt->Z();
      if ( genie::pdg::IsProton( tgt->HitNucPdg()) ) --Zf;
      Mf = genie::PDGLibrary::Instance()->Find( genie::pdg::IonPdgCode(Af, Zf) )->Mass();

      // Deduce the binding energy from the final nucleus mass
      Eb = Mf - Mi + mNi;
    }

    // The (lab-frame) off-shell initial nucleon energy is the difference
    // between the lab frame total energies of the initial and remnant nuclei
    ENi = Mi - std::sqrt( Mf*Mf + p3Ni.Mag2() );
  }
  else {
    // Keep the struck nucleon on shell either because
    // hitNucleonBindingMode == kOnShell or because
    // the target is a single nucleon
    ENi = std::sqrt( p3Ni.Mag2() + std::pow(mNi, 2) );
    Eb = 0.;

    // If we're dealing with a nuclear target but using the on-shell
    // binding mode, set the "assume free nucleon" flag. This turns off
    // Pauli blocking and the requirement that q0 > 0 in the QE cross section
    // models (an on-shell nucleon *is* a free nucleon)
    if ( tgt->IsNucleus() ) interaction.SetBit( kIAssumeFreeNucleon );
  }

  // Update the initial nucleon lab-frame 4-momentum in the interaction with
  // its current components
  p4Ni->SetVect( p3Ni );
  p4Ni->SetE( ENi );

}

void genie::utils::BindHitNucleon	(	genie::Interaction &	interaction,
		const NuclearModelI &	nucl_model,
		double &	Eb,
		genie::QELEvGen_BindingMode_t	hitNucleonBindingMode
	)

Definition at line 261 of file QELUtils.cxx.

References genie::Target::A(), genie::PDGLibrary::Find(), genie::gAbortingInErr, genie::Target::HitNucP4Ptr(), genie::Target::HitNucPdg(), genie::Target::HitNucPosition(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCode(), genie::pdg::IsAntiNeutrino(), genie::ProcessInfo::IsEM(), genie::pdg::IsNeutrino(), genie::Target::IsNucleus(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), genie::kIAssumeFreeNucleon, genie::constants::kNucleonMass, genie::kOnShell, genie::kUseGroundStateRemnant, genie::kUseNuclearModel, genie::kValenciaStyleQValue, genie::NuclearModelI::LocalFermiMomentum(), LOG, genie::utils::res::Mass(), genie::Target::Mass(), genie::NuclearModelI::Momentum3(), pDEBUG, pFATAL, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), genie::utils::mec::Qvalue(), genie::Interaction::RecoilNucleonPdg(), genie::NuclearModelI::RemovalEnergy(), genie::InitialState::TgtPtr(), and genie::Target::Z().

{
  genie::Target* tgt = interaction.InitState().TgtPtr();
  TLorentzVector* p4Ni = tgt->HitNucP4Ptr();

  // Initial nucleon 3-momentum (lab frame)
  TVector3 p3Ni = nucl_model.Momentum3();

  // Look up the (on-shell) mass of the initial nucleon
  TDatabasePDG* tb = TDatabasePDG::Instance();
  double mNi = tb->GetParticle( tgt->HitNucPdg() )->Mass();

  // Set the (possibly off-shell) initial nucleon energy based on
  // the selected binding energy mode. Always put the initial nucleon
  // on shell if it is not part of a composite nucleus
  double ENi = 0.;
  if ( tgt->IsNucleus() && hitNucleonBindingMode != genie::kOnShell ) {

    // For a nuclear target with a bound initial struck nucleon, take binding
    // energy effects and Pauli blocking into account when computing QE
    // differential cross sections
    interaction.ResetBit( kIAssumeFreeNucleon );

    // Initial nucleus mass
    double Mi = tgt->Mass();

    // Final nucleus mass
    double Mf = 0.;

    // If we use the removal energy reported by the nuclear
    // model, then it implies a certain value for the final
    // nucleus mass
    if ( hitNucleonBindingMode == genie::kUseNuclearModel ) {
      Eb = nucl_model.RemovalEnergy();
      // This equation is the definition that we assume
      // here for the "removal energy" (Eb) returned by the
      // nuclear model. It matches GENIE's convention for
      // the Bodek/Ritchie Fermi gas model.
      Mf = Mi + Eb - mNi;
    }
    // We can also assume that the final nucleus is in its
    // ground state. In this case, we can just look up its
    // mass from the standard table. This implies a particular
    // binding energy for the hit nucleon.
    else if ( hitNucleonBindingMode == genie::kUseGroundStateRemnant ) {
      // Determine the mass and proton numbers for the remnant nucleus
      int Af = tgt->A() - 1;
      int Zf = tgt->Z();
      if ( genie::pdg::IsProton( tgt->HitNucPdg()) ) --Zf;
      Mf = genie::PDGLibrary::Instance()->Find( genie::pdg::IonPdgCode(Af, Zf) )->Mass();

      // Deduce the binding energy from the final nucleus mass
      Eb = Mf - Mi + mNi;
    }
    // A third option is to assign the binding energy and final nuclear
    // mass in a way that is equivalent to the original Valencia model
    // treatment
    else if ( hitNucleonBindingMode == genie::kValenciaStyleQValue ) {
      // Compute the Q-value needed for the ground-state-to-ground-state
      // transition without nucleon removal (see Sec. III B of
      // https://arxiv.org/abs/nucl-th/0408005). Note that this will be
      // zero for NC and EM scattering since the proton and nucleon numbers
      // are unchanged. Also get Q_LFG, the difference in Fermi energies
      // between the final and initial struck nucleon species. This will
      // likewise be zero for NC and EM interactions.
      const genie::ProcessInfo& info = interaction.ProcInfo();
      double Qvalue = 0.;
      double Q_LFG = 0.;
      if ( info.IsWeakCC() ) {
        // Without nucleon removal, the final nucleon number remains the same
        int Af = tgt->A();
        // For CC interactions, the proton number will change by one. Choose
        // the right change by checking whether we're working with a neutrino
        // or an antineutrino.
        int Zf = tgt->Z();
        int probe_pdg = interaction.InitState().ProbePdg();
        if ( genie::pdg::IsAntiNeutrino(probe_pdg) ) {
          --Zf;
        }
        else if ( genie::pdg::IsNeutrino(probe_pdg) ) {
          ++Zf;
        }
        else {
          LOG( "QELEvent", pFATAL ) << "Unhandled probe PDG code " << probe_pdg
            << " encountered for a CC interaction in"
            << " genie::utils::BindHitNucleon()";
          gAbortingInErr = true;
          std::exit( 1 );
        }

        // We have what we need to get the Q-value. Get the final nuclear
        // mass (without nucleon removal)
        double mf_keep_nucleon = genie::PDGLibrary::Instance()
          ->Find( genie::pdg::IonPdgCode(Af, Zf) )->Mass();

        Qvalue = mf_keep_nucleon - Mi;

        // Get the Fermi energies for the initial and final nucleons. Include
        // the radial dependence if using the LFG.
        double hit_nucleon_radius = tgt->HitNucPosition();

        // Average of the proton and neutron masses. It may actually be better
        // to use the exact on-shell masses here. However, the original paper
        // uses the same value for protons and neutrons when computing the
        // difference in Fermi energies. For consistency with the Valencia
        // model publication, I'll do the same.
        const double mN = genie::constants::kNucleonMass;

        double kF_Ni = nucl_model.LocalFermiMomentum( *tgt,
          tgt->HitNucPdg(), hit_nucleon_radius );
        double EFermi_Ni = std::sqrt( std::max(0., mN*mN + kF_Ni*kF_Ni) );

        double kF_Nf = nucl_model.LocalFermiMomentum( *tgt,
          interaction.RecoilNucleonPdg(), hit_nucleon_radius );
        // (On-shell) final nucleon mass
        //double mNf = interaction.RecoilNucleon()->Mass();
        double EFermi_Nf = std::sqrt( std::max(0., mN*mN + kF_Nf*kF_Nf) );

        // The difference in Fermi energies is Q_LFG
        Q_LFG = EFermi_Nf - EFermi_Ni;
      }
      // Fail here for interactions that aren't one of EM/NC/CC. This is
      // intended to help avoid silently doing the wrong thing in the future.
      else if ( !info.IsEM() && !info.IsWeakNC() ) {
        LOG( "QELEvent", pFATAL ) << "Unhandled process type \"" << info
          << "\" encountered in genie::utils::BindHitNucleon()";
        gAbortingInErr = true;
        std::exit( 1 );
      }

      // On-shell total energy of the initial struck nucleon
      double ENi_OnShell = std::sqrt( std::max(0., mNi*mNi + p3Ni.Mag2()) );

      // Total energy of the remnant nucleus (with the hit nucleon removed)
      double Ef = Mi - ENi_OnShell + Qvalue - Q_LFG;

      // Mass and kinetic energy of the remnant nucleus (with the hit nucleon
      // removed)
      Mf = std::sqrt( std::max(0., Ef*Ef - p3Ni.Mag2()) );

      // Deduce the binding energy from the final nucleus mass
      Eb = Mf - Mi + mNi;

      LOG( "QELEvent", pDEBUG ) << "Qvalue = " << Qvalue
        << ", Q_LFG = " << Q_LFG << " at radius = " << tgt->HitNucPosition();
    }

    // The (lab-frame) off-shell initial nucleon energy is the difference
    // between the lab frame total energies of the initial and remnant nuclei
    ENi = Mi - std::sqrt( Mf*Mf + p3Ni.Mag2() );
  }
  else {
    // Keep the struck nucleon on shell either because
    // hitNucleonBindingMode == kOnShell or because
    // the target is a single nucleon
    ENi = std::sqrt( p3Ni.Mag2() + std::pow(mNi, 2) );
    Eb = 0.;

    // If we're dealing with a nuclear target but using the on-shell
    // binding mode, set the "assume free nucleon" flag. This turns off
    // Pauli blocking and the requirement that q0 > 0 in the QE cross section
    // models (an on-shell nucleon *is* a free nucleon)
    if ( tgt->IsNucleus() ) interaction.SetBit( kIAssumeFreeNucleon );
  }

  LOG( "QELEvent", pDEBUG ) << "Eb = " << Eb << ", pNi = " << p3Ni.Mag()
    << ", ENi = " << ENi;

  // Update the initial nucleon lab-frame 4-momentum in the interaction with
  // its current components
  p4Ni->SetVect( p3Ni );
  p4Ni->SetE( ENi );

}

double genie::utils::ComputeFullDMELPXSec	(	genie::Interaction *	interaction,
		const NuclearModelI *	nucl_model,
		const XSecAlgorithmI *	xsec_model,
		double	cos_theta_0,
		double	phi_0,
		double &	Eb,
		genie::DMELEvGen_BindingMode_t	hitNucleonBindingMode,
		double	min_angle_EM = 0.,
		bool	bind_nucleon = true
	)

Definition at line 94 of file DMELUtils.cxx.

References BindHitNucleon(), genie::InitialState::CMEnergy(), genie::Interaction::FSPrimLepton(), genie::InitialState::GetProbeP4(), genie::Interaction::InitState(), genie::ProcessInfo::IsEM(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::constants::kPi, genie::kPSDMELEvGen, genie::kRfLab, genie::utils::res::Mass(), genie::Range1D_t::max, genie::Interaction::PhaseSpace(), genie::Interaction::ProcInfo(), genie::utils::kinematics::Q2(), genie::KPhaseSpace::Q2Lim(), genie::Interaction::RecoilNucleonPdg(), genie::units::s, genie::Kinematics::SetFSLeptonP4(), genie::Kinematics::SetHadSystP4(), genie::Kinematics::SetQ2(), and genie::XSecAlgorithmI::XSec().

Referenced by genie::DMELEventGenerator::ComputeMaxXSec(), and genie::DMELEventGenerator::ProcessEventRecord().

{
  // If requested, set the initial hit nucleon 4-momentum to be off-shell
  // according to the binding mode specified in the function call
  if ( bind_nucleon ) {
    genie::utils::BindHitNucleon(*interaction, *nucl_model, Eb,
      hitNucleonBindingMode);
  }

  // Mass of the outgoing lepton
  double lepMass = interaction->FSPrimLepton()->Mass();

  // Look up the (on-shell) mass of the final nucleon
  TDatabasePDG *tb = TDatabasePDG::Instance();
  double mNf = tb->GetParticle( interaction->RecoilNucleonPdg() )->Mass();

  // Mandelstam s for the probe/hit nucleon system
  double s = std::pow( interaction->InitState().CMEnergy(), 2 );

  // Return a differential cross section of zero if we're below threshold (and
  // therefore need to sample a new event)
  if ( std::sqrt(s) < lepMass + mNf ) return 0.;

  double outLeptonEnergy = ( s - mNf*mNf + lepMass*lepMass ) / (2 * std::sqrt(s));

  if (outLeptonEnergy*outLeptonEnergy - lepMass*lepMass < 0.) return 0.;
  double outMomentum = TMath::Sqrt(outLeptonEnergy*outLeptonEnergy - lepMass*lepMass);

  // Compute the boost vector for moving from the COM frame to the
  // lab frame, i.e., the velocity of the COM frame as measured
  // in the lab frame.
  TVector3 beta = COMframe2Lab( interaction->InitState() );

  // FullDifferentialXSec depends on theta_0 and phi_0, the lepton COM
  // frame angles with respect to the direction of the COM frame velocity
  // as measured in the lab frame. To generate the correct dependence
  // here, first set the lepton COM frame angles with respect to +z
  // (via TVector3::SetTheta() and TVector3::SetPhi()).
  TVector3 lepton3Mom(0., 0., outMomentum);
  lepton3Mom.SetTheta( TMath::ACos(cos_theta_0) );
  lepton3Mom.SetPhi( phi_0 );

  // Then rotate the lepton 3-momentum so that the old +z direction now
  // points along the COM frame velocity (beta)
  TVector3 zvec(0., 0., 1.);
  TVector3 rot = ( zvec.Cross(beta) ).Unit();
  double angle = beta.Angle( zvec );

  // Handle the edge case where beta is along -z, so the
  // cross product above vanishes
  if ( beta.Perp() == 0. && beta.Z() < 0. ) {
    rot = TVector3(0., 1., 0.);
    angle = genie::constants::kPi;
  }

  // Rotate if the rotation vector is not 0
  if ( rot.Mag() >= genie::controls::kASmallNum ) {
    lepton3Mom.Rotate(angle, rot);
  }

  // Construct the lepton 4-momentum in the COM frame
  TLorentzVector lepton(lepton3Mom, outLeptonEnergy);

  // The final state nucleon will have an equal and opposite 3-momentum
  // in the COM frame and will be on the mass shell
  TLorentzVector outNucleon(-1*lepton.Px(),-1*lepton.Py(),-1*lepton.Pz(),
    TMath::Sqrt(outMomentum*outMomentum + mNf*mNf));

  // Boost the 4-momenta for both particles into the lab frame
  lepton.Boost(beta);
  outNucleon.Boost(beta);

  // Check if event is at a low angle - if so return 0 and stop wasting time
  if (180 * lepton.Theta() / genie::constants::kPi < min_angle_EM && interaction->ProcInfo().IsEM()) {
    return 0;
  }

  TLorentzVector * nuP4 = interaction->InitState().GetProbeP4( genie::kRfLab );
  TLorentzVector qP4 = *nuP4 - lepton;
  delete nuP4;
  double Q2 = -1 * qP4.Mag2();

  interaction->KinePtr()->SetFSLeptonP4( lepton );
  interaction->KinePtr()->SetHadSystP4( outNucleon );
  interaction->KinePtr()->SetQ2( Q2 );

  // Check the Q2 range. If we're outside of it, don't bother
  // with the rest of the calculation.
  Range1D_t Q2lim = interaction->PhaseSpace().Q2Lim();
  if (Q2 < Q2lim.min || Q2 > Q2lim.max) return 0.;

  // Compute the QE cross section for the current kinematics
  double xsec = xsec_model->XSec(interaction, genie::kPSDMELEvGen);

  return xsec;
}

double genie::utils::ComputeFullQELPXSec	(	genie::Interaction *	interaction,
		const NuclearModelI *	nucl_model,
		const XSecAlgorithmI *	xsec_model,
		double	cos_theta_0,
		double	phi_0,
		double &	Eb,
		genie::QELEvGen_BindingMode_t	hitNucleonBindingMode,
		double	min_angle_EM = 0.,
		bool	bind_nucleon = true
	)

Definition at line 93 of file QELUtils.cxx.

References BindHitNucleon(), genie::InitialState::CMEnergy(), genie::Interaction::FSPrimLepton(), genie::InitialState::GetProbeP4(), genie::Interaction::InitState(), genie::ProcessInfo::IsEM(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::constants::kPi, genie::kPSQELEvGen, genie::kRfLab, genie::utils::res::Mass(), genie::Range1D_t::max, genie::Interaction::PhaseSpace(), genie::Interaction::ProcInfo(), genie::utils::kinematics::Q2(), genie::KPhaseSpace::Q2Lim(), genie::Interaction::RecoilNucleonPdg(), genie::units::s, genie::Kinematics::SetFSLeptonP4(), genie::Kinematics::SetHadSystP4(), genie::Kinematics::SetQ2(), and genie::XSecAlgorithmI::XSec().

Referenced by genie::QELEventGenerator::ComputeMaxXSec(), genie::utils::gsl::FullQELdXSec::DoEval(), and genie::QELEventGenerator::ProcessEventRecord().

{
  // If requested, set the initial hit nucleon 4-momentum to be off-shell
  // according to the binding mode specified in the function call
  if ( bind_nucleon ) {
    genie::utils::BindHitNucleon(*interaction, *nucl_model, Eb,
      hitNucleonBindingMode);
  }

  // Mass of the outgoing lepton
  double lepMass = interaction->FSPrimLepton()->Mass();

  // Look up the (on-shell) mass of the final nucleon
  TDatabasePDG *tb = TDatabasePDG::Instance();
  double mNf = tb->GetParticle( interaction->RecoilNucleonPdg() )->Mass();

  // Mandelstam s for the probe/hit nucleon system
  double s = std::pow( interaction->InitState().CMEnergy(), 2 );

  // Return a differential cross section of zero if we're below threshold (and
  // therefore need to sample a new event)
  if ( std::sqrt(s) < lepMass + mNf ) return 0.;

  double outLeptonEnergy = ( s - mNf*mNf + lepMass*lepMass ) / (2 * std::sqrt(s));

  if (outLeptonEnergy*outLeptonEnergy - lepMass*lepMass < 0.) return 0.;
  double outMomentum = TMath::Sqrt(outLeptonEnergy*outLeptonEnergy - lepMass*lepMass);

  // Compute the boost vector for moving from the COM frame to the
  // lab frame, i.e., the velocity of the COM frame as measured
  // in the lab frame.
  TVector3 beta = COMframe2Lab( interaction->InitState() );

  // FullDifferentialXSec depends on theta_0 and phi_0, the lepton COM
  // frame angles with respect to the direction of the COM frame velocity
  // as measured in the lab frame. To generate the correct dependence
  // here, first set the lepton COM frame angles with respect to +z
  // (via TVector3::SetTheta() and TVector3::SetPhi()).
  TVector3 lepton3Mom(0., 0., outMomentum);
  lepton3Mom.SetTheta( TMath::ACos(cos_theta_0) );
  lepton3Mom.SetPhi( phi_0 );

  // Then rotate the lepton 3-momentum so that the old +z direction now
  // points along the COM frame velocity (beta)
  TVector3 zvec(0., 0., 1.);
  TVector3 rot = ( zvec.Cross(beta) ).Unit();
  double angle = beta.Angle( zvec );

  // Handle the edge case where beta is along -z, so the
  // cross product above vanishes
  if ( beta.Perp() == 0. && beta.Z() < 0. ) {
    rot = TVector3(0., 1., 0.);
    angle = genie::constants::kPi;
  }

  // Rotate if the rotation vector is not 0
  if ( rot.Mag() >= genie::controls::kASmallNum ) {
    lepton3Mom.Rotate(angle, rot);
  }

  // Construct the lepton 4-momentum in the COM frame
  TLorentzVector lepton(lepton3Mom, outLeptonEnergy);

  // The final state nucleon will have an equal and opposite 3-momentum
  // in the COM frame and will be on the mass shell
  TLorentzVector outNucleon(-1*lepton.Px(),-1*lepton.Py(),-1*lepton.Pz(),
    TMath::Sqrt(outMomentum*outMomentum + mNf*mNf));

  // Boost the 4-momenta for both particles into the lab frame
  lepton.Boost(beta);
  outNucleon.Boost(beta);

  // Check if event is at a low angle - if so return 0 and stop wasting time
  if (180 * lepton.Theta() / genie::constants::kPi < min_angle_EM && interaction->ProcInfo().IsEM()) {
    return 0;
  }

  TLorentzVector * nuP4 = interaction->InitState().GetProbeP4( genie::kRfLab );
  TLorentzVector qP4 = *nuP4 - lepton;
  delete nuP4;
  double Q2 = -1 * qP4.Mag2();

  interaction->KinePtr()->SetFSLeptonP4( lepton );
  interaction->KinePtr()->SetHadSystP4( outNucleon );
  interaction->KinePtr()->SetQ2( Q2 );

  // Check the Q2 range. If we're outside of it, don't bother
  // with the rest of the calculation.
  Range1D_t Q2lim = interaction->PhaseSpace().Q2Lim();
  if (Q2 < Q2lim.min || Q2 > Q2lim.max) return 0.;

  // Compute the QE cross section for the current kinematics
  double xsec = xsec_model->XSec(interaction, genie::kPSQELEvGen);

  return xsec;
}

double genie::utils::CosTheta0Max

(

const genie::Interaction &

interaction

)

Definition at line 217 of file DMELUtils.cxx.

References genie::InitialState::CMEnergy(), genie::Interaction::FSPrimLepton(), genie::Target::HitNucMass(), genie::Target::HitNucP4(), genie::Interaction::InitState(), genie::Target::IsNucleus(), genie::kIAssumeFreeNucleon, genie::kRfLab, genie::InitialState::ProbeE(), genie::Interaction::RecoilNucleon(), and genie::InitialState::Tgt().

Referenced by genie::QELEventGenerator::ComputeMaxXSec(), genie::DMELEventGenerator::ComputeMaxXSec(), genie::QELEventGenerator::ProcessEventRecord(), and genie::DMELEventGenerator::ProcessEventRecord().

                                                                   {

  // q0 > 0 only needs to be enforced (indirectly via a calculation of
  // CosTheta0Max) for bound nucleons. The Q2 limits should take care of valid
  // kinematics for free nucleons.
  if ( !interaction.InitState().Tgt().IsNucleus()
    || interaction.TestBit(kIAssumeFreeNucleon) ) return 1.;

  double probe_E_lab = interaction.InitState().ProbeE( genie::kRfLab );

  TVector3 beta = COMframe2Lab( interaction.InitState() );
  double gamma = 1. / std::sqrt(1. - beta.Mag2());

  double sqrt_s = interaction.InitState().CMEnergy();
  double mNf = interaction.RecoilNucleon()->Mass();
  double ml = interaction.FSPrimLepton()->Mass();
  double lepton_E_COM = (sqrt_s*sqrt_s + ml*ml - mNf*mNf) / (2.*sqrt_s);

  // If there isn't enough available energy to create an on-shell
  // final lepton, then don't bother with the rest of the calculation
  // NOTE: C++11 would allow use to use lowest() here instead
  if ( lepton_E_COM <= ml ) return -std::numeric_limits<double>::max();

  double lepton_p_COM = std::sqrt( std::max(lepton_E_COM*lepton_E_COM
    - ml*ml, 0.) );

  // Possibly off-shell initial struck nucleon total energy
  // (BindHitNucleon() should have been called previously if needed)
  const TLorentzVector& p4Ni = interaction.InitState().Tgt().HitNucP4();
  double ENi = p4Ni.E();
  // On-shell mass of initial struck nucleon
  double mNi = interaction.InitState().Tgt().HitNucMass();
  // On-shell initial struck nucleon energy
  double ENi_on_shell = std::sqrt( mNi*mNi + p4Ni.Vect().Mag2() );
  // Energy needed to put initial nucleon on the mass shell
  double epsilon_B = ENi_on_shell - ENi;

  double cos_theta0_max = ( probe_E_lab - gamma*lepton_E_COM - epsilon_B )
    / ( gamma * lepton_p_COM * beta.Mag() );
  return cos_theta0_max;
}

double genie::utils::EnergyDeltaFunctionSolutionDMEL

(

const Interaction &

inter

)

Definition at line 51 of file DMELUtils.cxx.

References genie::Kinematics::FSLeptonP4(), genie::InitialState::GetProbeP4(), genie::Kinematics::HadSystP4(), genie::Target::HitNucP4(), genie::Interaction::InitStatePtr(), genie::Interaction::Kine(), genie::kRfLab, and genie::InitialState::TgtPtr().

{
  // Get the final-state lepton and struck nucleon 4-momenta in the lab frame
  const TLorentzVector& lepton = inter.Kine().FSLeptonP4();
  const TLorentzVector& outNucleon = inter.Kine().HadSystP4();

  // Get beta for a Lorentz boost from the lab frame to the COM frame and
  // vice-versa
  TLorentzVector* probe = inter.InitStatePtr()->GetProbeP4( kRfLab );
  const TLorentzVector& hit_nucleon = inter.InitStatePtr()->TgtPtr()
    ->HitNucP4();
  TLorentzVector total_p4 = (*probe) + hit_nucleon;
  TVector3 beta_COM_to_lab = total_p4.BoostVector();
  TVector3 beta_lab_to_COM = -beta_COM_to_lab;

  // Square of the Lorentz gamma factor for the boosts
  double gamma2 = std::pow(total_p4.Gamma(), 2);

  // Get the final-state lepton 4-momentum in the COM frame
  TLorentzVector leptonCOM = TLorentzVector( lepton );
  leptonCOM.Boost( beta_lab_to_COM );

  // Angle between COM outgoing lepton and lab-frame velocity of COM frame
  double theta0 = leptonCOM.Angle( beta_COM_to_lab );
  double cos_theta0_squared = std::pow(std::cos( theta0 ), 2);
  //double sin_theta0_squared = 1. - cos_theta0_squared; // sin^2(x) + cos^2(x) = 1

  // Difference in lab frame velocities between lepton and outgoing nucleon
  TVector3 lepton_vel = ( lepton.Vect() * (1. / lepton.E()) );
  TVector3 outNucleon_vel = ( outNucleon.Vect() * (1. / outNucleon.E()) );
  double vel_diff = ( lepton_vel - outNucleon_vel ).Mag();

  // Compute the factor in the kPSDMELEvGen phase space that arises from
  // eliminating the energy-conserving delta function
  double factor = std::sqrt( 1. + (1. - cos_theta0_squared)*(gamma2 - 1.) )
    * std::pow(leptonCOM.P(), 2) / vel_diff;

  delete probe;

  return factor;
}

double genie::utils::EnergyDeltaFunctionSolutionQEL

(

const Interaction &

inter

)

Definition at line 50 of file QELUtils.cxx.

References genie::Kinematics::FSLeptonP4(), genie::InitialState::GetProbeP4(), genie::Kinematics::HadSystP4(), genie::Target::HitNucP4(), genie::Interaction::InitStatePtr(), genie::Interaction::Kine(), genie::kRfLab, and genie::InitialState::TgtPtr().

Referenced by genie::LwlynSmithQELCCPXSec::FullDifferentialXSec(), and genie::NievesQELCCPXSec::XSec().

{
  // Get the final-state lepton and struck nucleon 4-momenta in the lab frame
  const TLorentzVector& lepton = inter.Kine().FSLeptonP4();
  const TLorentzVector& outNucleon = inter.Kine().HadSystP4();

  // Get beta for a Lorentz boost from the lab frame to the COM frame and
  // vice-versa
  TLorentzVector* probe = inter.InitStatePtr()->GetProbeP4( kRfLab );
  const TLorentzVector& hit_nucleon = inter.InitStatePtr()->TgtPtr()
    ->HitNucP4();
  TLorentzVector total_p4 = (*probe) + hit_nucleon;
  TVector3 beta_COM_to_lab = total_p4.BoostVector();
  TVector3 beta_lab_to_COM = -beta_COM_to_lab;

  // Square of the Lorentz gamma factor for the boosts
  double gamma2 = std::pow(total_p4.Gamma(), 2);

  // Get the final-state lepton 4-momentum in the COM frame
  TLorentzVector leptonCOM = TLorentzVector( lepton );
  leptonCOM.Boost( beta_lab_to_COM );

  // Angle between COM outgoing lepton and lab-frame velocity of COM frame
  double theta0 = leptonCOM.Angle( beta_COM_to_lab );
  double cos_theta0_squared = std::pow(std::cos( theta0 ), 2);
  //double sin_theta0_squared = 1. - cos_theta0_squared; // sin^2(x) + cos^2(x) = 1

  // Difference in lab frame velocities between lepton and outgoing nucleon
  TVector3 lepton_vel = ( lepton.Vect() * (1. / lepton.E()) );
  TVector3 outNucleon_vel = ( outNucleon.Vect() * (1. / outNucleon.E()) );
  double vel_diff = ( lepton_vel - outNucleon_vel ).Mag();

  // Compute the factor in the kPSQELEvGen phase space that arises from
  // eliminating the energy-conserving delta function
  double factor = std::sqrt( 1. + (1. - cos_theta0_squared)*(gamma2 - 1.) )
    * std::pow(leptonCOM.P(), 2) / vel_diff;

  delete probe;

  return factor;
}

ostream & genie::utils::operator<<	(	ostream &	stream,
		const T2KEvGenMetaData &	md
	)

Definition at line 22 of file T2KEvGenMetaData.cxx.

References genie::utils::T2KEvGenMetaData::Print().

  {
     md.Print(stream);
     return stream;
  }

void genie::utils::SetPrimaryLeptonPolarization

(

GHepRecord *

ev

)

Definition at line 23 of file PrimaryLeptonUtils.cxx.

References genie::GHepRecord::FinalStatePrimaryLepton(), genie::pdg::IsElectron(), genie::pdg::IsMuon(), genie::pdg::IsNeutrino(), genie::pdg::IsTau(), LOG, genie::GHepParticle::Name(), genie::GHepParticle::Pdg(), pERROR, pINFO, genie::GHepParticle::PolzAzimuthAngle(), genie::GHepParticle::PolzIsSet(), genie::GHepParticle::PolzPolarAngle(), genie::GHepParticle::Px(), genie::GHepParticle::Py(), genie::GHepParticle::Pz(), and genie::GHepParticle::SetPolarization().

Referenced by genie::MECGenerator::AddFinalStateLepton(), genie::QELEventGenerator::ProcessEventRecord(), genie::QELEventGeneratorSM::ProcessEventRecord(), genie::MECGenerator::SelectNSVLeptonKinematics(), and genie::PrimaryLeptonGenerator::SetPolarization().

{
// Moved out of the PrimaryLeptonGenerator class to make the same treatment
// accessible for generators that use a more unified approach (e.g.,
// QELEventGenerator and MECGenerator). -- S. Gardiner

// Set the final state lepton polarization. A mass-less lepton would be fully
// polarized. This would be exact for neutrinos and a very good approximation
// for electrons for the energies this generator is going to be used. This is
// not the case for muons and, mainly, for taus. I need to refine this later.
// How? See Kuzmin, Lyubushkin and Naumov, hep-ph/0312107

  // get the final state primary lepton
  GHepParticle * fsl = ev->FinalStatePrimaryLepton();
  if ( !fsl ) {
    LOG("LeptonicVertex", pERROR)
      << "Final state lepton not set yet! \n" << *ev;
    return;
  }

  // Get (px,py,pz) @ LAB
  TVector3 plab( fsl->Px(), fsl->Py(), fsl->Pz() );

  // In the limit m/E->0: leptons are left-handed and their anti-particles
  // are right-handed
  int pdgc = fsl->Pdg();
  if ( pdg::IsNeutrino(pdgc) || pdg::IsElectron(pdgc) ||
    pdg::IsMuon(pdgc) || pdg::IsTau(pdgc) )
  {
    plab *= -1; // left-handed
  }

  LOG("LeptonicVertex", pINFO)
    << "Setting polarization angles for particle: " << fsl->Name();

  fsl->SetPolarization( plab );

  if ( fsl->PolzIsSet() ) {
    LOG("LeptonicVertex", pINFO)
      << "Polarization (rad): Polar = "  << fsl->PolzPolarAngle()
      << ", Azimuthal = " << fsl->PolzAzimuthAngle();
  }
}

genie::DMELEvGen_BindingMode_t genie::utils::StringToDMELBindingMode

(

const std::string &

mode_str

)

Definition at line 195 of file DMELUtils.cxx.

References genie::gAbortingInErr, genie::kOnShell, genie::kUseGroundStateRemnant, genie::kUseNuclearModel, LOG, and pFATAL.

Referenced by genie::DMELEventGenerator::LoadConfig().

{
  // Translate the string setting the binding mode to the appropriate
  // enum value, or complain if one couldn't be found
  if ( binding_mode == "UseGroundStateRemnant" ) {
    return kUseGroundStateRemnant;
  }
  else if ( binding_mode == "UseNuclearModel" ) {
    return kUseNuclearModel;
  }
  else if ( binding_mode == "OnShell" ) {
    return kOnShell;
  }
  else {
    LOG("DMELEvent", pFATAL) << "Unrecognized setting \"" << binding_mode
      << "\" requested in genie::utils::StringToDMELBindingMode()";
    gAbortingInErr = true;
    std::exit(1);
  }
}

genie::QELEvGen_BindingMode_t genie::utils::StringToQELBindingMode

(

const std::string &

mode_str

)

Definition at line 194 of file QELUtils.cxx.

References genie::gAbortingInErr, genie::kOnShell, genie::kUseGroundStateRemnant, genie::kUseNuclearModel, genie::kValenciaStyleQValue, LOG, and pFATAL.

Referenced by genie::NewQELXSec::Integrate(), genie::QELEventGenerator::LoadConfig(), genie::LwlynSmithQELCCPXSec::LoadConfig(), and genie::NievesQELCCPXSec::LoadConfig().

{
  // Translate the string setting the binding mode to the appropriate
  // enum value, or complain if one couldn't be found
  if ( binding_mode == "UseGroundStateRemnant" ) {
    return kUseGroundStateRemnant;
  }
  else if ( binding_mode == "UseNuclearModel" ) {
    return kUseNuclearModel;
  }
  else if ( binding_mode == "OnShell" ) {
    return kOnShell;
  }
  else if ( binding_mode == "ValenciaStyleQValue" ) {
    return kValenciaStyleQValue;
  }
  else {
    LOG("QELEvent", pFATAL) << "Unrecognized setting \"" << binding_mode
      << "\" requested in genie::utils::StringToQELBindingMode()";
    gAbortingInErr = true;
    std::exit(1);
  }
}