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.

	intranuke

	intranuke2018

	kinematics
	Kinematical utilities.

	math
	Simple mathematical utilities not found in ROOT's TMath.

	mec
	MEC utilities.

	nhl

	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	neutron_osc
	Utilities for simulating neutron oscillation. More...

class	nhl
	Utilities for simulating the decay of Neutral Heavy Leptons. More...

class	nucleon_decay
	Utilities for simulating nucleon decay. More...

class	T2KEvGenMetaData
	Utility class to store MC job meta-data. 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

May 01, 2020

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.

 {
   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.

 {
   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.

 {
   // 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.

 {
   // 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.

                                                                    {
 
   // 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.

 {
   // 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.

 {
   // 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.

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

void genie::utils::SetPrimaryLeptonPolarization ( GHepRecord * ev )

Definition at line 23 of file PrimaryLeptonUtils.cxx.

 {
 // 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.

 {
   // 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.

 {
   // 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);
   }
 }

Namespaces

Classes

Functions

Detailed Description

Function Documentation