Skip to content

Commit

Permalink
Remove the kPSTnctnBnctl phase space since it is unique to
Browse files Browse the repository at this point in the history
QELEventGenerator and DMELEventGenerator and is never actually used for
sampling kinematics. Also remove DMELEventGenerator for now, as it appears
to have never been fully implemented. After fixes to QELEventGenerator have
been applied, DMELEventGenerator can be rewritten from a better
starting point.
  • Loading branch information
sjgardiner committed Mar 22, 2019
1 parent 1e237c6 commit 9ea3078
Show file tree
Hide file tree
Showing 13 changed files with 106 additions and 1,057 deletions.
44 changes: 0 additions & 44 deletions config/DMELEventGenerator.xml

This file was deleted.

1 change: 0 additions & 1 deletion config/master_config.xml
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@
<config alg="genie::VertexGenerator"> VertexGenerator.xml </config>
<config alg="genie::QELEventGenerator"> QELEventGenerator.xml </config>
<config alg="genie::QELEventGeneratorSM"> QELEventGeneratorSM.xml </config>
<config alg="genie::DMELEventGenerator"> DMELEventGenerator.xml </config>
<config alg="genie::QELPrimaryLeptonGenerator"> QELPrimaryLeptonGenerator.xml </config>
<config alg="genie::DISPrimaryLeptonGenerator"> DISPrimaryLeptonGenerator.xml </config>
<config alg="genie::RESPrimaryLeptonGenerator"> RESPrimaryLeptonGenerator.xml </config>
Expand Down
2 changes: 0 additions & 2 deletions src/Framework/Conventions/KinePhaseSpace.h
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,6 @@ typedef enum EKinePhaseSpace {
kPSElOlOpifE,
kPSElOlTpifE,
kPSTkTlctl,
kPSTnctnBnctl, // Nucleon momentum, angle, binding energy, lepton com angle
kPSQ2vfE,
kPSQELEvGen // Phase space used by genie::QELEventGenerator for sampling kinematic variables
// TODO: rename this value when the correct variables are identified
Expand Down Expand Up @@ -116,7 +115,6 @@ class KinePhaseSpace
case(kPSElOlOpifE) : return "<{El,Omega_l,Omega_pi}|E>"; break;
case(kPSElOlTpifE) : return "<{El,Omega_l,Theta_pi}|E>"; break;
case(kPSTkTlctl) : return "<{Tk,Tl,cos(theta_l)}|E>"; break;
case(kPSTnctnBnctl) : return "<centre-of-mass plep,Omega_lep, p_p, omega_p|E>"; break;
case(kPSQ2vfE) : return "<{Q2,v}|E>"; break;
// TODO: update this string when the appropriate kinematic variables are known
case(kPSQELEvGen) : return "<QELEvGen>"; break;
Expand Down
40 changes: 0 additions & 40 deletions src/Framework/Utils/KineUtils.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -251,46 +251,6 @@ double genie::utils::kinematics::Jacobian(
J = 2*constants::kPi;
}

// Transformation: centre-of-mass plep,Omega_lep, p_p, omega_p|E
// -> QELEvGen
// TODO: update this comment when kPSQELEvGen has been renamed
else if ( TransformMatched(fromps, tops, kPSTnctnBnctl, kPSQELEvGen, forward) )
{
// Get the final-state lepton and struck nucleon 4-momenta in the lab frame
const TLorentzVector& lepton = i->Kine().FSLeptonP4();
const TLorentzVector& outNucleon = i->Kine().HadSystP4();

// Get beta for a Lorentz boost from the lab frame to the COM frame and
// vice-versa
TLorentzVector* probe = i->InitStatePtr()->GetProbeP4(kRfLab);
const TLorentzVector& target = i->InitStatePtr()->TgtPtr()
->HitNucP4();
TLorentzVector total_p4 = (*probe) + target;
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 lab-frame muon velocity and velocity of COM frame
// as measured in lab frame
double theta_J = lepton.Angle( beta_COM_to_lab );
double cos_theta_J2 = std::pow(std::cos( theta_J ), 2);

// 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 Jacobian
J = 4. * kPi * leptonCOM.P() * leptonCOM.P()
* std::sqrt( 1. + (1. - cos_theta_J2)*(gamma2 - 1) ) / vel_diff;
}

else {
SLOG("KineLimits", pFATAL)
<< "*** Can not compute Jacobian for transforming: "
Expand Down
Loading

0 comments on commit 9ea3078

Please sign in to comment.