Simulate the primary MEC interaction. More...

#include <MECGenerator.h>

Inheritance diagram for genie::MECGenerator:

Collaboration diagram for genie::MECGenerator:

Public Member Functions
	MECGenerator ()

	MECGenerator (string config)

	~MECGenerator ()

void	ProcessEventRecord (GHepRecord *event) const

void	Configure (const Registry &config)

void	Configure (string config)

Public Member Functions inherited from genie::EventRecordVisitorI
virtual	~EventRecordVisitorI ()

Public Member Functions inherited from genie::Algorithm
virtual	~Algorithm ()

virtual void	FindConfig (void)

virtual const Registry &	GetConfig (void) const

Registry *	GetOwnedConfig (void)

virtual const AlgId &	Id (void) const
	Get algorithm ID. More...

virtual AlgStatus_t	GetStatus (void) const
	Get algorithm status. More...

virtual bool	AllowReconfig (void) const

virtual AlgCmp_t	Compare (const Algorithm *alg) const
	Compare with input algorithm. More...

virtual void	SetId (const AlgId &id)
	Set algorithm ID. More...

virtual void	SetId (string name, string config)

const Algorithm *	SubAlg (const RgKey &registry_key) const

void	AdoptConfig (void)

void	AdoptSubstructure (void)

virtual void	Print (ostream &stream) const
	Print algorithm info. More...

Private Member Functions
void	LoadConfig (void)

void	AddNucleonCluster (GHepRecord *event) const

void	AddTargetRemnant (GHepRecord *event) const

void	GenerateFermiMomentum (GHepRecord *event) const

void	SelectEmpiricalKinematics (GHepRecord *event) const

void	AddFinalStateLepton (GHepRecord *event) const

void	RecoilNucleonCluster (GHepRecord *event) const

void	DecayNucleonCluster (GHepRecord *event) const

void	SelectNSVLeptonKinematics (GHepRecord *event) const

void	SelectSuSALeptonKinematics (GHepRecord *event) const

void	GenerateNSVInitialHadrons (GHepRecord *event) const

PDGCodeList	NucleonClusterConstituents (int pdgc) const

double	GetXSecMaxTlctl (const Interaction &inter, const Range1D_t &Tl_range, const Range1D_t &ctl_range) const

Private Attributes
const XSecAlgorithmI *	fXSecModel

TGenPhaseSpace	fPhaseSpaceGenerator

const NuclearModelI *	fNuclModel

double	fSafetyFactor

int	fFunctionCalls

double	fRelTolerance

int	fMinScanPointsTmu

int	fMinScanPointsCosth

double	fQ3Max

double	fSuSAMaxXSecDiffTolerance

Additional Inherited Members
Static Public Member Functions inherited from genie::Algorithm
static string	BuildParamVectKey (const std::string &comm_name, unsigned int i)

static string	BuildParamVectSizeKey (const std::string &comm_name)

static string	BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)

static string	BuildParamMatRowSizeKey (const std::string &comm_name)

static string	BuildParamMatColSizeKey (const std::string &comm_name)

Protected Member Functions inherited from genie::EventRecordVisitorI
	EventRecordVisitorI ()

	EventRecordVisitorI (string name)

	EventRecordVisitorI (string name, string config)

Protected Member Functions inherited from genie::Algorithm
	Algorithm ()

	Algorithm (string name)

	Algorithm (string name, string config)

void	Initialize (void)

void	DeleteConfig (void)

void	DeleteSubstructure (void)

Registry *	ExtractLocalConfig (const Registry &in) const

Registry *	ExtractLowerConfig (const Registry &in, const string &alg_key) const
	Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...

template<class T >
bool	GetParam (const RgKey &name, T &p, bool is_top_call=true) const

template<class T >
bool	GetParamDef (const RgKey &name, T &p, const T &def) const

template<class T >
int	GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
	Handle to load vectors of parameters. More...

int	GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

template<class T >
int	GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
	Handle to load matrix of parameters. More...

template<class T >
int	GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const

int	GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

int	AddTopRegistry (Registry *rp, bool owns=true)
	add registry with top priority, also update ownership More...

int	AddLowRegistry (Registry *rp, bool owns=true)
	add registry with lowest priority, also update ownership More...

int	MergeTopRegistry (const Registry &r)

int	AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
	Add registries with top priority, also udated Ownerships. More...

Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig

bool	fOwnsSubstruc
	true if it owns its substructure (sub-algs,...) More...

AlgId	fID
	algorithm name and configuration set More...

vector< Registry * >	fConfVect

vector< bool >	fOwnerships
	ownership for every registry in fConfVect More...

AlgStatus_t	fStatus
	algorithm execution status More...

AlgMap *	fOwnedSubAlgMp
	local pool for owned sub-algs (taken out of the factory pool) More...

Detailed Description

Simulate the primary MEC interaction.

Author: Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool

Steve Dytman <dytman+ pitt.edu> Pittsburgh University

Created:: Sep. 22, 2008

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

Definition at line 36 of file MECGenerator.h.

Constructor & Destructor Documentation

MECGenerator::MECGenerator ( )

Definition at line 51 of file MECGenerator.cxx.

                            :
 EventRecordVisitorI("genie::MECGenerator")
 {
 
 }

MECGenerator::MECGenerator ( string config )

Definition at line 57 of file MECGenerator.cxx.

                                         :
 EventRecordVisitorI("genie::MECGenerator", config)
 {
 
 }

MECGenerator::~MECGenerator ( )

Definition at line 63 of file MECGenerator.cxx.

 {
 
 }

Member Function Documentation

void MECGenerator::AddFinalStateLepton ( GHepRecord * event ) const

private

Definition at line 326 of file MECGenerator.cxx.

References genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4(), genie::Interaction::InitState(), genie::RandomGen::Instance(), genie::Interaction::Kine(), genie::kIStStableFinalState, genie::constants::kPi, genie::kRfHitNucRest, LOG, pNOTICE, genie::GHepRecord::Probe(), genie::InitialState::ProbeE(), genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), genie::RandomGen::RndLep(), genie::utils::SetPrimaryLeptonPolarization(), genie::InitialState::Tgt(), genie::GHepParticle::X4(), and genie::Kinematics::y().

Referenced by ProcessEventRecord().

 {
 // Add the final-state primary lepton in the event record.
 // Compute its 4-momentum based on the selected interaction kinematics.
 //
   Interaction * interaction = event->Summary();
 
   // The boost back to the lab frame was missing, that is included now with the introduction of the beta factor
   const InitialState & init_state = interaction->InitState();
   const TLorentzVector & pnuc4 = init_state.Tgt().HitNucP4(); //[@LAB]
   TVector3 beta = pnuc4.BoostVector();
 
   // Boosting the incoming neutrino to the NN-cluster rest frame
   // Neutrino 4p
   TLorentzVector * p4v = event->Probe()->GetP4(); // v 4p @ LAB
   p4v->Boost(-1.*beta);                           // v 4p @ NN-cluster rest frame
 
   // Look-up selected kinematics
   double Q2 = interaction->Kine().Q2(true);
   double y  = interaction->Kine().y(true);
 
   // Auxiliary params
   double Ev  = interaction->InitState().ProbeE(kRfHitNucRest);
   LOG("MEC", pNOTICE) << "neutrino energy = " << Ev;
   double ml  = interaction->FSPrimLepton()->Mass();
   double ml2 = TMath::Power(ml,2);
 
   // Compute the final state primary lepton energy and momentum components
   // along and perpendicular the neutrino direction
   double El  = (1-y)*Ev;
   double plp = El - 0.5*(Q2+ml2)/Ev;                          // p(//)
   double plt = TMath::Sqrt(TMath::Max(0.,El*El-plp*plp-ml2)); // p(-|)
 
   LOG("MEC", pNOTICE)
           << "fsl: E = " << El << ", |p//| = " << plp << ", |pT| = " << plt;
 
   // Randomize transverse components
   RandomGen * rnd = RandomGen::Instance();
   double phi  = 2*kPi * rnd->RndLep().Rndm();
   double pltx = plt * TMath::Cos(phi);
   double plty = plt * TMath::Sin(phi);
 
   // Take a unit vector along the neutrino direction
   // WE NEED THE UNIT VECTOR ALONG THE NEUTRINO DIRECTION IN THE NN-CLUSTER REST FRAME, NOT IN THE LAB FRAME
  TVector3 unit_nudir = p4v->Vect().Unit();      //We use this one, which is in the NN-cluster rest frame
   // Rotate lepton momentum vector from the reference frame (x'y'z') where
   // {z':(neutrino direction), z'x':(theta plane)} to the LAB
   TVector3 p3l(pltx,plty,plp);
   p3l.RotateUz(unit_nudir);
 
   // Lepton 4-momentum in LAB
   TLorentzVector p4l(p3l,El);
 
   // Boost final state primary lepton to the lab frame
   p4l.Boost(beta); // active Lorentz transform
 
   // Figure out the final-state primary lepton PDG code
   int pdgc = interaction->FSPrimLepton()->PdgCode();
 
   // Lepton 4-position (= interacton vtx)
   TLorentzVector v4(*event->Probe()->X4());
 
   // Add the final-state lepton to the event record
   int momidx = event->ProbePosition();
   event->AddParticle(
     pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
 
   // Set its polarization
   utils::SetPrimaryLeptonPolarization( event );
 }

void genie::MECGenerator::AddNucleonCluster ( GHepRecord * event ) const

private

void MECGenerator::AddTargetRemnant ( GHepRecord * event ) const

private

Definition at line 117 of file MECGenerator.cxx.

References genie::units::A, genie::GHepParticle::A(), genie::pdg::IonPdgCode(), genie::kIStStableFinalState, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::GHepParticle::P4(), genie::GHepParticle::Pdg(), and genie::GHepParticle::Z().

Referenced by ProcessEventRecord().

 {
 // Add the remnant nucleus (= initial nucleus - nucleon cluster) in the
 // event record.
 
   GHepParticle * target  = event->TargetNucleus();
   GHepParticle * cluster = event->HitNucleon();
 
   int Z = target->Z();
   int A = target->A();
 
   if(cluster->Pdg() == kPdgClusterNN) { A-=2; ;     }
   if(cluster->Pdg() == kPdgClusterNP) { A-=2; Z-=1; }
   if(cluster->Pdg() == kPdgClusterPP) { A-=2; Z-=2; }
 
   int ipdgc = pdg::IonPdgCode(A, Z);
 
   const TLorentzVector & p4cluster = *(cluster->P4());
   const TLorentzVector & p4tgt     = *(target ->P4());
 
   const TLorentzVector p4 = p4tgt - p4cluster;
   const TLorentzVector v4(0.,0.,0., 0.);
 
   int momidx = event->TargetNucleusPosition();
 
   /*
   if( A == 2 && Z == 2){
     // residual nucleus was just two protons, not a nucleus we know.
     // this might not conserve energy...
     event->AddParticle(kPdgProton,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
     // v4,v4 because I'm lazy, give the four momentum to one of the protons, not the other
     event->AddParticle(kPdgProton,kIStStableFinalState, momidx,-1,-1,-1, v4,v4);
   } else if ( A == 2 && Z == 0){
     // residual nucleus was just two neutrons, not a nucleus we know.
     // this might not conserve energy...
     event->AddParticle(kPdgNeutron,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
     // v4,v4 because I'm lazy, give the four momentum to one of the protons, not the other
     event->AddParticle(kPdgNeutron,kIStStableFinalState, momidx,-1,-1,-1, v4,v4);
   } else {
     // regular nucleus, including deuterium
     event->AddParticle(ipdgc,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
   }
   */
 
   event->AddParticle(ipdgc,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
 
 }

void MECGenerator::Configure ( const Registry & config )

virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 1321 of file MECGenerator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

 {
     Algorithm::Configure(config);
     this->LoadConfig();
 }

void MECGenerator::Configure ( string config )

virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 1327 of file MECGenerator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

 {
     Algorithm::Configure(config);
     this->LoadConfig();
 }

void MECGenerator::DecayNucleonCluster ( GHepRecord * event ) const

private

accept_decay

Definition at line 436 of file MECGenerator.cxx.

References genie::PDGLibrary::Find(), fPhaseSpaceGenerator, genie::GHepParticle::GetP4(), genie::GHepParticle::GetX4(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::kHadroSysGenErr, genie::kIStHadronInTheNucleus, genie::controls::kMaxUnweightDecayIterations, LOG, genie::units::m, NucleonClusterConstituents(), genie::utils::print::P4AsString(), genie::GHepParticle::Pdg(), pERROR, pINFO, pNOTICE, pWARN, genie::RandomGen::RndDec(), genie::exceptions::EVGThreadException::SetReason(), genie::exceptions::EVGThreadException::SetReturnStep(), and genie::exceptions::EVGThreadException::SwitchOnStepBack().

Referenced by ProcessEventRecord().

 {
 // Perform a phase-space decay of the nucleon cluster and add its decay
 // products in the event record
 //
   LOG("MEC", pINFO) << "Decaying nucleon cluster...";
 
   // get di-nucleon cluster
   int nucleon_cluster_id = 5;
   GHepParticle * nucleon_cluster = event->Particle(nucleon_cluster_id);
   assert(nucleon_cluster);
 
   // get decay products
   PDGCodeList pdgv = this->NucleonClusterConstituents(nucleon_cluster->Pdg());
   LOG("MEC", pINFO) << "Decay product IDs: " << pdgv;
 
   // Get the decay product masses
   vector<int>::const_iterator pdg_iter;
   int i = 0;
   double * mass = new double[pdgv.size()];
   double   sum  = 0;
   for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
     int pdgc = *pdg_iter;
     double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
     mass[i++] = m;
     sum += m;
   }
 
   LOG("MEC", pINFO)
     << "Performing a phase space decay to "
     << pdgv.size() << " particles / total mass = " << sum;
 
   TLorentzVector * p4d = nucleon_cluster->GetP4();
   TLorentzVector * v4d = nucleon_cluster->GetX4();
 
   LOG("MEC", pINFO)
     << "Decaying system p4 = " << utils::print::P4AsString(p4d);
 
   // Set the decay
   bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
   if(!permitted) {
      LOG("MEC", pERROR)
        << " *** Phase space decay is not permitted \n"
        << " Total particle mass = " << sum << "\n"
        << " Decaying system p4 = " << utils::print::P4AsString(p4d);
      // clean-up
      delete [] mass;
      delete p4d;
      delete v4d;
      // throw exception
      event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
      genie::exceptions::EVGThreadException exception;
      exception.SetReason("Decay not permitted kinematically");
      exception.SwitchOnStepBack();
      exception.SetReturnStep(0);
      throw exception;
   }
 
   // Get the maximum weight
   double wmax = -1;
   for(int idec=0; idec<200; idec++) {
      double w = fPhaseSpaceGenerator.Generate();
      wmax = TMath::Max(wmax,w);
   }
   assert(wmax>0);
   wmax *= 2;
 
   LOG("MEC", pNOTICE)
      << "Max phase space gen. weight = " << wmax;
 
   RandomGen * rnd = RandomGen::Instance();
   bool accept_decay=false;
   unsigned int itry=0;
   while(!accept_decay)
   {
      itry++;
 
      if(itry > controls::kMaxUnweightDecayIterations) {
        // report, clean-up and return
        LOG("MEC", pWARN)
            << "Couldn't generate an unweighted phase space decay after "
            << itry << " attempts";
        // clean up
        delete [] mass;
        delete p4d;
        delete v4d;
        // throw exception
        event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
        genie::exceptions::EVGThreadException exception;
        exception.SetReason("Couldn't select decay after N attempts");
        exception.SwitchOnStepBack();
        exception.SetReturnStep(0);
        throw exception;
      }
      double w  = fPhaseSpaceGenerator.Generate();
      if(w > wmax) {
         LOG("MEC", pWARN)
            << "Decay weight = " << w << " > max decay weight = " << wmax;
      }
      double gw = wmax * rnd->RndDec().Rndm();
      accept_decay = (gw<=w);
 
      LOG("MEC", pINFO)
         << "Decay weight = " << w << " / R = " << gw
         << " - accepted: " << accept_decay;
 
   } //!accept_decay
 
   // Insert the decay products in the event record
   TLorentzVector v4(*v4d);
   GHepStatus_t ist = kIStHadronInTheNucleus;
   int idp = 0;
   for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
      int pdgc = *pdg_iter;
      TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
      event->AddParticle(pdgc, ist, nucleon_cluster_id,-1,-1,-1, *p4fin, v4);
      idp++;
   }
 
   // Clean-up
   delete [] mass;
   delete p4d;
   delete v4d;
 }

void MECGenerator::GenerateFermiMomentum ( GHepRecord * event ) const

private

Definition at line 165 of file MECGenerator.cxx.

References genie::PDGLibrary::Find(), fNuclModel, genie::NuclearModelI::GenerateNucleon(), genie::PDGLibrary::Instance(), LOG, genie::utils::res::Mass(), genie::NuclearModelI::Momentum3(), NucleonClusterConstituents(), genie::GHepParticle::Pdg(), pINFO, and genie::GHepParticle::SetMomentum().

Referenced by ProcessEventRecord().

 {
 // Generate the initial state di-nucleon cluster 4-momentum.
 // Draw Fermi momenta for each of the two nucleons.
 // Sum the two Fermi momenta to calculate the di-nucleon momentum.
 // For simplicity, keep the di-nucleon cluster on the mass shell.
 //
   GHepParticle * target_nucleus = event->TargetNucleus();
   assert(target_nucleus);
   GHepParticle * nucleon_cluster = event->HitNucleon();
   assert(nucleon_cluster);
   GHepParticle * remnant_nucleus = event->RemnantNucleus();
   assert(remnant_nucleus);
 
   // generate a Fermi momentum for each nucleon
 
   Target tgt(target_nucleus->Pdg());
   PDGCodeList pdgv = this->NucleonClusterConstituents(nucleon_cluster->Pdg());
   assert(pdgv.size()==2);
   tgt.SetHitNucPdg(pdgv[0]);
   fNuclModel->GenerateNucleon(tgt);
   TVector3 p3a = fNuclModel->Momentum3();
   tgt.SetHitNucPdg(pdgv[1]);
   fNuclModel->GenerateNucleon(tgt);
   TVector3 p3b = fNuclModel->Momentum3();
 
   LOG("FermiMover", pINFO)
      << "1st nucleon (code = " << pdgv[0] << ") generated momentum: ("
      << p3a.Px() << ", " << p3a.Py() << ", " << p3a.Pz() << "), "
      << "|p| = " << p3a.Mag();
   LOG("FermiMover", pINFO)
      << "2nd nucleon (code = " << pdgv[1] << ") generated momentum: ("
      << p3b.Px() << ", " << p3b.Py() << ", " << p3b.Pz() << "), "
      << "|p| = " << p3b.Mag();
 
   // calcute nucleon cluster momentum
 
   TVector3 p3 = p3a + p3b;
 
   LOG("FermiMover", pINFO)
      << "di-nucleon cluster momentum: ("
      << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
      << "|p| = " << p3.Mag();
 
   // target (initial) nucleus and nucleon cluster mass
 
   double Mi  = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass(); // initial nucleus mass
   double M2n = PDGLibrary::Instance()->Find(nucleon_cluster->Pdg())-> Mass(); // nucleon cluster mass
 
   // nucleon cluster energy
 
   double EN = TMath::Sqrt(p3.Mag2() + M2n*M2n);
 
   // set the remnant nucleus and nucleon cluster 4-momenta at GHEP record
 
   TLorentzVector p4nclust   (   p3.Px(),    p3.Py(),    p3.Pz(),  EN   );
   TLorentzVector p4remnant   (-1*p3.Px(), -1*p3.Py(), -1*p3.Pz(), Mi-EN);
 
   nucleon_cluster->SetMomentum(p4nclust);
   remnant_nucleus->SetMomentum(p4remnant);
 
   // set the nucleon cluster 4-momentum at the interaction summary
 
   event->Summary()->InitStatePtr()->TgtPtr()->SetHitNucP4(p4nclust);
 }

void MECGenerator::GenerateNSVInitialHadrons ( GHepRecord * event ) const

private

Definition at line 1126 of file MECGenerator.cxx.

References genie::PDGLibrary::Find(), fNuclModel, genie::NuclearModelI::GenerateNucleon(), genie::PDGLibrary::Instance(), genie::Interaction::KinePtr(), genie::kIStDecayedState, genie::kKineGenErr, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::controls::kRjMaxIterations, LOG, genie::utils::res::Mass(), genie::NuclearModelI::Momentum3(), NucleonClusterConstituents(), genie::GHepParticle::P4(), pDEBUG, genie::GHepParticle::Pdg(), pERROR, pWARN, genie::NuclearModelI::RemovalEnergy(), genie::Kinematics::SetHadSystP4(), genie::GHepParticle::SetMomentum(), genie::exceptions::EVGThreadException::SetReason(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), and genie::GHepParticle::X4().

Referenced by ProcessEventRecord().

 {
     // We need a kinematic limits accept/reject loop here, so generating the
     // initial hadrons is combined with generating the recoil hadrons...
 
     LOG("MEC",pDEBUG) << "Generate Initial Hadrons - Start";
 
     Interaction * interaction = event->Summary();
     GHepParticle * neutrino = event->Probe();
     assert(neutrino);
     TLorentzVector p4nu(*neutrino->P4());
 
     // get final state primary lepton & its 4-momentum
     GHepParticle * fsl = event->FinalStatePrimaryLepton();
     assert(fsl);
     TLorentzVector p4l(*fsl->P4());
 
     // calculate 4-momentum transfer from these
     TLorentzVector Q4 = p4nu - p4l;
 
     // get the target two-nucleon cluster and nucleus.
     // the remnant nucleus is apparently set, except for its momentum.
     GHepParticle * target_nucleus = event->TargetNucleus();
     assert(target_nucleus);
     GHepParticle * initial_nucleon_cluster = event->HitNucleon();
     assert(initial_nucleon_cluster);
     GHepParticle * remnant_nucleus = event->RemnantNucleus();
     assert(remnant_nucleus);
 
     // -- make a two-nucleon system, then give them some momenta.
 
     // instantiate an empty local target nucleus, so I can use existing methods
     // to get a momentum from the prevailing Fermi-motion distribution
     Target tgt(target_nucleus->Pdg());
 
     // NucleonClusterConstituents is an implementation within this class, called with this
     // It using the nucleon cluster from the earlier tests for a pn state,
     // the method returns a vector of pdgs, which hopefully will be of size two.
 
     PDGCodeList pdgv = this->NucleonClusterConstituents(initial_nucleon_cluster->Pdg());
     assert(pdgv.size()==2);
 
     // These things need to be saved through to the end of the accept loop.
     bool accept = false;
     TVector3 p31i;
     TVector3 p32i;
     unsigned int iter = 0;
 
     int initial_nucleon_cluster_pdg = initial_nucleon_cluster->Pdg();
     int final_nucleon_cluster_pdg = 0;
 
     //e-scat case:
     if (neutrino->Pdg() == 11) {
       final_nucleon_cluster_pdg = initial_nucleon_cluster->Pdg();
     }
     //neutrino case
     else if (neutrino->Pdg() > 0) {
         if (initial_nucleon_cluster->Pdg() == kPdgClusterNP) {
             final_nucleon_cluster_pdg = kPdgClusterPP;
         }
         else if (initial_nucleon_cluster->Pdg() == kPdgClusterNN) {
             final_nucleon_cluster_pdg = kPdgClusterNP;
         }
         else {
             LOG("MEC", pERROR) << "Wrong pdg for a CC neutrino MEC interaction"
                 << initial_nucleon_cluster->Pdg();
         }
     }
     else if (neutrino->Pdg() < 0) {
         if (initial_nucleon_cluster->Pdg() == kPdgClusterNP) {
             final_nucleon_cluster_pdg = kPdgClusterNN;
         }
         else if (initial_nucleon_cluster->Pdg() == kPdgClusterPP) {
             final_nucleon_cluster_pdg = kPdgClusterNP;
         }
         else {
             LOG("MEC", pERROR) << "Wrong pdg for a CC anti-neutrino MEC interaction"
                 << initial_nucleon_cluster->Pdg();
         }
     }
 
     TLorentzVector p4initial_cluster;
     TLorentzVector p4final_cluster;
     TLorentzVector p4remnant_nucleus;
     double removalenergy1;
     double removalenergy2;
 
     //===========================================================================
     // Choose two nucleons from the prevailing fermi-motion distribution.
     // Some produce kinematically unallowed combinations initial cluster and Q2
     // Find out, and if so choose them again with this accept/reject loop.
     // Some kinematics are especially tough
     while(!accept){
         iter++;
         if(iter > kRjMaxIterations*1000) {
             // error if try too many times
             LOG("MEC", pWARN)
                 << "Couldn't select a valid W, Q^2 pair after "
                 << iter << " iterations";
             event->EventFlags()->SetBitNumber(kKineGenErr, true);
             genie::exceptions::EVGThreadException exception;
             exception.SetReason("Couldn't select initial hadron kinematics");
             exception.SwitchOnFastForward();
             throw exception;
         }
 
         // generate nucleons
         // tgt is a Target object for local use, just waiting to be filled.
         // this sets the pdg of each nucleon and its momentum from a Fermi-gas
         // Nieves et al. would use a local Fermi gas here, not this, but ok.
         // so momentum from global Fermi gas, local Fermi gas, or spectral function
         // and removal energy ~0.025 GeV, correlated with density, or from SF distribution
         tgt.SetHitNucPdg(pdgv[0]);
         fNuclModel->GenerateNucleon(tgt);
         p31i = fNuclModel->Momentum3();
         removalenergy1 = fNuclModel->RemovalEnergy();
         tgt.SetHitNucPdg(pdgv[1]);
         fNuclModel->GenerateNucleon(tgt);
         p32i = fNuclModel->Momentum3();
         removalenergy2 = fNuclModel->RemovalEnergy();
 
         // not sure -- could give option to use Nieves q-value here.
 
         // Now write down the initial cluster four-vector for this choice
         TVector3 p3i = p31i + p32i;
         double mass2 = PDGLibrary::Instance()->Find( initial_nucleon_cluster_pdg )->Mass();
         mass2 *= mass2;
         double energy = TMath::Sqrt(p3i.Mag2() + mass2);
         p4initial_cluster.SetPxPyPzE(p3i.Px(),p3i.Py(),p3i.Pz(),energy);
 
         // cast the removal energy as the energy component of a 4-vector for later.
         TLorentzVector tLVebind(0., 0., 0., -1.0 * (removalenergy1 + removalenergy2));
 
         // RIK: You might ask why is this the right place to subtract ebind?
         // It is okay. Physically, I'm subtracting it from q. The energy
         // transfer to the nucleon is 50 MeV less. the energy cost to put this
         // cluster on-shell. What Jan says he does in PRC.86.015504 is this:
         //   The nucleons are assumed to be in a potential well
         //     V = Efermi + 8 MeV.
         //   The Fermi energy is subtracted from each initial-state nucleon
         //   (I guess he does it dynamically based on Ef = p2/2M or so which
         //   is what we are doing above, on average). Then after the reaction,
         // another 8 MeV is subtracted at that point; a small adjustment to the
         // momentum is needed to keep the nucleon on shell.
 
         // calculate recoil nucleon cluster 4-momentum (tLVebind is negative)
         p4final_cluster = p4initial_cluster + Q4 + tLVebind;
 
         // Test if the resulting four-vector corresponds to a high-enough invariant mass.
         // Fail the accept if we couldn't put this thing on-shell.
         if (p4final_cluster.M() <
                 PDGLibrary::Instance()->Find(final_nucleon_cluster_pdg )->Mass()) {
             accept = false;
         } else {
             accept = true;
         }
 
     }  // end accept loop
 
     // we got here if we accepted the final state two-nucleon system
     // so now we need to write everything to ghep
 
     // First the initial state nucleons.
     initial_nucleon_cluster->SetMomentum(p4initial_cluster);
 
     // and the remnant nucleus
     double Mi  = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass();
     remnant_nucleus->SetMomentum(-1.0*p4initial_cluster.Px(),
             -1.0*p4initial_cluster.Py(),
             -1.0*p4initial_cluster.Pz(),
             Mi - p4initial_cluster.E() + removalenergy1 + removalenergy2);
 
     // Now the final nucleon cluster.
 
     // Getting this v4 is a position, i.e. a position within the nucleus (?)
     // possibly it takes on a non-trivial value only for local Fermi gas
     // or for sophisticated treatments of intranuclear rescattering.
     TLorentzVector v4(*neutrino->X4());
 
     // Now write the (undecayed) final two-nucleon system
     GHepParticle p1(final_nucleon_cluster_pdg, kIStDecayedState,
             2, -1, -1, -1, p4final_cluster, v4);
 
     //p1.SetRemovalEnergy(removalenergy1 + removalenergy2);
     // The "bound particle" concept applies only to p or n.
     // Instead, add this directly to the remnant nucleon a few lines above.
 
     // actually, this is not an status1 particle, so it is not picked up
     // by the aggregator. anyway, the aggregator does not run until the very end.
 
     event->AddParticle(p1);
 
     interaction->KinePtr()->SetHadSystP4(p4final_cluster);
 }

double MECGenerator::GetXSecMaxTlctl	(	const Interaction &	inter,
		const Range1D_t &	Tl_range,
		const Range1D_t &	ctl_range
	)		const

private

Definition at line 1354 of file MECGenerator.cxx.

References fFunctionCalls, fMinScanPointsCosth, fMinScanPointsTmu, fSafetyFactor, genie::Interaction::FSPrimLepton(), fXSecModel, genie::Interaction::InitState(), genie::kRfHitNucRest, genie::Range1D_t::min, and genie::InitialState::ProbeE().

Referenced by SelectNSVLeptonKinematics().

                                                                           {
 
   ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit2");
 
   double Enu = in.InitState().ProbeE(kRfHitNucRest);
   double LepMass = in.FSPrimLepton()->Mass();
 
   genie::utils::mec::gsl::d2Xsec_dTCosth f(fXSecModel,in, Enu, LepMass, -1.) ;
 
   std::array<string,2> names = { "Tl", "CosThetal" } ;
   std::array<Range1D_t,2> ranges = { Tl_range, ctl_range } ;
 
   std::array<double,2> start, steps, temp_point ;
   steps[0] = ( ranges[0].max - ranges[0].min ) / ( fMinScanPointsTmu + 1 ) ;
   steps[1] = ( ranges[1].max - ranges[1].min ) / ( fMinScanPointsCosth + 1 ) ;
 
   double xsec = 0 ;
 
   // preliimnary scan
   for ( unsigned int i = 1 ; i <= (unsigned int) fMinScanPointsTmu ; ++i ) {
     temp_point[0] = ranges[0].min + steps[0]*i ;
 
     for ( unsigned int j = 1 ; j <= (unsigned int) fMinScanPointsCosth ; ++j ) {
       temp_point[1] = ranges[1].min + steps[1]*j ;
 
       double temp_xsec = - f( temp_point.data() ) ;
       if ( temp_xsec > xsec ) {
         start = temp_point ;
         xsec = temp_xsec ;
       }
     }
   }
 
   // Set minimizer function and absolute tolerance :
   min->SetFunction( f );
   min->SetMaxFunctionCalls(fFunctionCalls);
   //  min->SetTolerance( fRelTolerance * xsec );
 
   for ( unsigned int i = 0 ; i < ranges.size() ; ++i ) {
     min -> SetLimitedVariable( i, names[i], start[i], steps[i], ranges[i].min, ranges[i].max ) ;
   }
 
   min->Minimize();
 
   double max_xsec = -min->MinValue(); //back to positive xsec
 
   // Apply safety factor, since value retrieved from the cache might
   // correspond to a slightly different energy.
   max_xsec *= fSafetyFactor;
 
   return max_xsec ;
 }

void MECGenerator::LoadConfig ( void )

private

Definition at line 1333 of file MECGenerator.cxx.

References fFunctionCalls, fMinScanPointsCosth, fMinScanPointsTmu, fNuclModel, fQ3Max, fRelTolerance, fSafetyFactor, fSuSAMaxXSecDiffTolerance, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), and genie::Algorithm::SubAlg().

Referenced by Configure().

 {
     fNuclModel = 0;
     RgKey nuclkey = "NuclearModel";
     fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
     assert(fNuclModel);
 
     GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.6 ) ;
     GetParam( "MaxXSec-FunctionCalls", fFunctionCalls ) ;
     GetParam( "MaxXSec-RelativeTolerance", fRelTolerance ) ;
     GetParam( "MaxXSec-MinScanPointsTmu", fMinScanPointsTmu ) ;
     GetParam( "MaxXSec-MinScanPointsCosth", fMinScanPointsCosth ) ;
 
     GetParam( "NSV-Q3Max", fQ3Max ) ;
 
     // Maximum allowed percentage deviation from the maximum cross section used
     // in the accept/reject loop for selecting lepton kinematics for SuSAv2.
     // Similar to the tolerance used by QELEventGenerator.
     GetParamDef( "SuSA-MaxXSec-DiffTolerance", fSuSAMaxXSecDiffTolerance, 999999. );
 }

PDGCodeList MECGenerator::NucleonClusterConstituents ( int pdgc ) const

private

Definition at line 561 of file MECGenerator.cxx.

References genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::kPdgNeutron, genie::kPdgProton, LOG, pERROR, and genie::PDGCodeList::push_back().

Referenced by DecayNucleonCluster(), GenerateFermiMomentum(), and GenerateNSVInitialHadrons().

 {
   bool allowdup = true;
   PDGCodeList pdgv(allowdup);
 
   if(pdgc == kPdgClusterNN) {
      pdgv.push_back(kPdgNeutron);
      pdgv.push_back(kPdgNeutron);
   }
   else
   if(pdgc == kPdgClusterNP) {
      pdgv.push_back(kPdgNeutron);
      pdgv.push_back(kPdgProton);
   }
   else
   if(pdgc == kPdgClusterPP) {
      pdgv.push_back(kPdgProton);
      pdgv.push_back(kPdgProton);
   }
   else
   {
      LOG("MEC", pERROR)
         << "Unknown di-nucleon cluster PDG code (" << pdgc << ")";
   }
 
   return pdgv;
 }

void MECGenerator::ProcessEventRecord ( GHepRecord * event ) const

virtual

shortly, this will be handled by the InitialStateAppender module

Implements genie::EventRecordVisitorI.

Definition at line 68 of file MECGenerator.cxx.

References AddFinalStateLepton(), AddTargetRemnant(), genie::EventGeneratorI::CrossSectionAlg(), DecayNucleonCluster(), fXSecModel, GenerateFermiMomentum(), GenerateNSVInitialHadrons(), genie::Algorithm::Id(), genie::RunningThreadInfo::Instance(), LOG, genie::AlgId::Name(), pFATAL, RecoilNucleonCluster(), genie::RunningThreadInfo::RunningThread(), SelectEmpiricalKinematics(), SelectNSVLeptonKinematics(), and SelectSuSALeptonKinematics().

 {
   //-- Access cross section algorithm for running thread
   RunningThreadInfo * rtinfo = RunningThreadInfo::Instance();
   const EventGeneratorI * evg = rtinfo->RunningThread();
   fXSecModel = evg->CrossSectionAlg();
   if (fXSecModel->Id().Name() == "genie::EmpiricalMECPXSec2015") {
       this -> AddTargetRemnant      (event); /// shortly, this will be handled by the InitialStateAppender module
       this -> GenerateFermiMomentum(event);
       this -> SelectEmpiricalKinematics(event);
       // TODO: test removing `AddFinalStateLepton` and replacing it with
       //    PrimaryLeptonGenerator::ProcessEventRecord(evrec);
       this -> AddFinalStateLepton(event);
       // TODO: maybe put `RecoilNucleonCluster` in a `MECHadronicSystemGenerator` class,
       // name it `GenerateEmpiricalDiNucleonCluster` or something...
       this -> RecoilNucleonCluster(event);
       // TODO: `DecayNucleonCluster` should probably be in `MECHadronicSystemGenerator`,
       // if we make that...
       this -> DecayNucleonCluster(event);
   } else if (fXSecModel->Id().Name() == "genie::NievesSimoVacasMECPXSec2016") {
       this -> SelectNSVLeptonKinematics(event);
       this -> AddTargetRemnant(event);
       this -> GenerateNSVInitialHadrons(event);
       // Note: this method in `MECTensor/MECTensorGenerator.cxx` appeared to be a straight
       // copy of an earlier version of the `DecayNucleonCluster` method here - but, watch
       // for this...
       this -> DecayNucleonCluster(event);
   }  else if (fXSecModel->Id().Name() == "genie::SuSAv2MECPXSec") {
       event->Print();
       this -> SelectSuSALeptonKinematics(event);
       event->Print();
       this -> AddTargetRemnant(event);
       event->Print();
       this -> GenerateNSVInitialHadrons(event);
       event->Print();
       // Note: this method in `MECTensor/MECTensorGenerator.cxx` appeared to be a straight
       // copy of an earlier version of the `DecayNucleonCluster` method here - but, watch
       // for this...
       this -> DecayNucleonCluster(event);
   }
   else {
       LOG("MECGenerator",pFATAL) <<
           "ProcessEventRecord >> Cannot calculate kinematics for " <<
           fXSecModel->Id().Name();
   }
 
 
 }

void MECGenerator::RecoilNucleonCluster ( GHepRecord * event ) const

private

Definition at line 397 of file MECGenerator.cxx.

References genie::GHepParticle::GetP4(), genie::kIStDecayedState, LOG, genie::GHepParticle::P4(), pINFO, and genie::GHepParticle::X4().

Referenced by ProcessEventRecord().

 {
   // get di-nucleon cluster & its 4-momentum
   GHepParticle * nucleon_cluster = event->HitNucleon();
   assert(nucleon_cluster);
   TLorentzVector * tmp=nucleon_cluster->GetP4();
   TLorentzVector p4cluster(*tmp);
   delete tmp;
 
   // get neutrino & its 4-momentum
   GHepParticle * neutrino = event->Probe();
   assert(neutrino);
   TLorentzVector p4v(*neutrino->P4());
 
   // get final state primary lepton & its 4-momentum
   GHepParticle * fsl = event->FinalStatePrimaryLepton();
   assert(fsl);
   TLorentzVector p4l(*fsl->P4());
 
   // calculate 4-momentum transfer
   TLorentzVector q = p4v - p4l;
 
   // calculate recoil nucleon cluster 4-momentum
   TLorentzVector p4cluster_recoil = p4cluster + q;
 
   // work-out recoil nucleon cluster code
   LOG("MEC", pINFO) << "Interaction summary";
   LOG("MEC", pINFO) << *event->Summary();
   int recoil_nucleon_cluster_pdg = event->Summary()->RecoilNucleonPdg();
 
   // vtx position in nucleus coord system
   TLorentzVector v4(*neutrino->X4());
 
   // add to the event record
   event->AddParticle(
     recoil_nucleon_cluster_pdg, kIStDecayedState,
     2, -1, -1, -1, p4cluster_recoil, v4);
 }

void MECGenerator::SelectEmpiricalKinematics ( GHepRecord * event ) const

private

Definition at line 231 of file MECGenerator.cxx.

References genie::Kinematics::ClearRunningValues(), genie::EventGeneratorI::CrossSectionAlg(), genie::PDGLibrary::Find(), fXSecModel, gQ2, genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::RunningThreadInfo::Instance(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::utils::mec::J(), genie::Interaction::KinePtr(), genie::kKineGenErr, genie::kPSWQ2fE, genie::kRfHitNucRest, genie::controls::kRjMaxIterations, LOG, pINFO, pNOTICE, genie::InitialState::ProbeE(), pWARN, genie::utils::kinematics::Q2(), genie::RandomGen::RndKine(), genie::RunningThreadInfo::RunningThread(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), genie::InitialState::Tgt(), genie::utils::kinematics::W(), genie::utils::kinematics::WQ2toXY(), and genie::XSecAlgorithmI::XSec().

Referenced by ProcessEventRecord().

 {
 // Select interaction kinematics using the rejection method.
 //
 
   // Access cross section algorithm for running thread
   RunningThreadInfo * rtinfo = RunningThreadInfo::Instance();
   const EventGeneratorI * evg = rtinfo->RunningThread();
   fXSecModel = evg->CrossSectionAlg();
 
   Interaction * interaction = event->Summary();
   double Ev = interaction->InitState().ProbeE(kRfHitNucRest);
 
   // **** NOTE / TODO:
   // **** Hardcode bogus limits for the time-being
   // **** Should be able to get limits via Interaction::KPhaseSpace
   double Q2min =  0.01;
   double Q2max =  8.00;
   double Wmin  =  1.88;
   double Wmax  =  3.00;
 
   // Scan phase-space for the maximum differential cross section
   // at the current neutrino energy
   const int nq=30;
   const int nw=20;
   double dQ2 = (Q2max-Q2min) / (nq-1);
   double dW  = (Wmax-Wmin )  / (nw-1);
   double xsec_max =  0;
   for(int iw=0; iw<nw; iw++) {
     for(int iq=0; iq<nq; iq++) {
       double Q2 = Q2min + iq*dQ2;
       double W  = Wmin  + iw*dW;
       interaction->KinePtr()->SetQ2(Q2);
       interaction->KinePtr()->SetW (W);
       double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
       xsec_max = TMath::Max(xsec, xsec_max);
     }
   }
   LOG("MEC", pNOTICE) << "xsec_max (E = " << Ev << " GeV) = " << xsec_max;
 
   // Select kinematics
   RandomGen * rnd = RandomGen::Instance();
   unsigned int iter = 0;
   bool accept = false;
   while(1) {
      iter++;
      if(iter > kRjMaxIterations) {
         LOG("MEC", pWARN)
            << "Couldn't select a valid W, Q^2 pair after "
            << iter << " iterations";
         event->EventFlags()->SetBitNumber(kKineGenErr, true);
         genie::exceptions::EVGThreadException exception;
         exception.SetReason("Couldn't select kinematics");
         exception.SwitchOnFastForward();
         throw exception;
      }
 
      // Generate next pair
      double gQ2 = Q2min + (Q2max-Q2min) * rnd->RndKine().Rndm();
      double gW  = Wmin  + (Wmax -Wmin ) * rnd->RndKine().Rndm();
 
      // Calculate d2sigma/dQ2dW
      interaction->KinePtr()->SetQ2(gQ2);
      interaction->KinePtr()->SetW (gW);
      double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
 
      // Decide whether to accept the current kinematics
      double t = xsec_max * rnd->RndKine().Rndm();
      double J = 1; // jacobean
      accept = (t < J*xsec);
 
      // If the generated kinematics are accepted, finish-up module's job
      if(accept) {
         LOG("MEC", pINFO) << "Selected: Q^2 = " << gQ2 << ", W = " << gW;
         double gx = 0;
         double gy = 0;
         //  More accurate calculation of the mass of the cluster than 2*Mnucl
         int nucleon_cluster_pdg = interaction->InitState().Tgt().HitNucPdg();
         double M2n = PDGLibrary::Instance()->Find(nucleon_cluster_pdg)->Mass();
         //bool is_em = interaction->ProcInfo().IsEM();
         kinematics::WQ2toXY(Ev,M2n,gW,gQ2,gx,gy);
 
         LOG("MEC", pINFO) << "x = " << gx << ", y = " << gy;
         // lock selected kinematics & clear running values
         interaction->KinePtr()->SetQ2(gQ2, true);
         interaction->KinePtr()->SetW (gW,  true);
         interaction->KinePtr()->Setx (gx,  true);
         interaction->KinePtr()->Sety (gy,  true);
         interaction->KinePtr()->ClearRunningValues();
 
         return;
      }//accepted?
   }//iter
 }

void MECGenerator::SelectNSVLeptonKinematics ( GHepRecord * event ) const

private

Definition at line 589 of file MECGenerator.cxx.

References genie::Interaction::ExclTagPtr(), fQ3Max, genie::Interaction::FSPrimLepton(), fXSecModel, genie::utils::mec::Getq0q3FromTlCostl(), GetXSecMaxTlctl(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::RandomGen::Instance(), genie::Interaction::KinePtr(), genie::kIStNucleonTarget, genie::kIStStableFinalState, genie::kKineGenErr, genie::kKVctl, genie::kKVQ0, genie::kKVQ3, genie::kKVTl, genie::kNoResonance, genie::constants::kNucleonMass, genie::kP33_1232, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::constants::kPi, genie::kPSTlctl, genie::kRfHitNucRest, genie::controls::kRjMaxIterations, LOG, pDEBUG, pERROR, pINFO, genie::GHepRecord::Probe(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), pWARN, genie::utils::kinematics::Q2(), genie::utils::kinematics::Q2YtoX(), genie::RandomGen::RndKine(), genie::RandomGen::RndLep(), genie::Kinematics::SetFSLeptonP4(), genie::Target::SetHitNucPdg(), genie::Kinematics::SetKV(), genie::utils::SetPrimaryLeptonPolarization(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::XclsTag::SetResonance(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), genie::InitialState::TgtPdg(), genie::InitialState::TgtPtr(), genie::GHepParticle::X4(), genie::XSecAlgorithmI::XSec(), and genie::utils::kinematics::XYtoW().

Referenced by ProcessEventRecord().

 {
   // -- implementation -- //
   // The IFIC Valencia model can provide three different hadron tensors.
   // The user probably wants all CC QE-like 2p2h events
   // But could in principle get the no-delta component if they want (deactivated incode)
   int FullDeltaNodelta = 1;  // 1:  full, 2:  only delta, 3:  zero delta
 
   // -- Event Properties -----------------------------//
   Interaction * interaction = event->Summary();
   Kinematics * kinematics = interaction->KinePtr();
 
   double Enu = interaction->InitState().ProbeE(kRfHitNucRest);
 
   int NuPDG = interaction->InitState().ProbePdg();
   int TgtPDG = interaction->InitState().TgtPdg();
   // interacton vtx
   TLorentzVector v4(*event->Probe()->X4());
   TLorentzVector tempp4(0.,0.,0.,0.);
 
   // -- Lepton Kinematic Limits ----------------------------------------- //
   double Costh = 0.0; // lepton angle
   double CosthMax = 1.0;
   double CosthMin = -1.0;
 
   double T = 0.0;  // lepton kinetic energy
   double TMax = std::numeric_limits<double>::max();
   double TMin = 0.0;
 
   double Plep = 0.0; // lepton 3 momentum
   double Elep = 0.0; // lepton energy
   double LepMass = interaction->FSPrimLepton()->Mass();
 
   double Q0 = 0.0; // energy component of q four vector
   double Q3 = 0.0; // magnitude of transfered 3 momentum
   double Q2 = 0.0; // properly Q^2 (Q squared) - transfered 4 momentum.
 
   // Set lepton KE TMax for for throwing rndm in the accept/reject loop.
   // We can accidentally set it too high, because the xsec will return zero.
   // This way if someone reuses this code, they are not tripped up by it.
   TMax = Enu - LepMass;
 
   // Warn if fQ3Max value is suspect 
   // (i.e. above the Valencia model's rage of validity)
   // This is just to warn users who have swapped out a MEC model
   // without remembering to change fQ3Max.
   if(fQ3Max>1.21){
     LOG("MEC", pWARN)
         << "fQ3 max is larger than expected for Valencia MEC: "
         << fQ3Max << ". Are you sure this is correct?";
   }
 
   // Set Tmin for throwing rndm in the accept/reject loop
   // the hadron tensors we expect will be limited in q3
   // therefore also the outgoing lepton KE can't be too low or costheta too backward
   // make the accept/reject loop more efficient by using Min values.
   if(Enu < fQ3Max){
     TMin = 0 ;
     CosthMin = -1 ;
   } else {
     TMin = TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu - fQ3Max), 2)) - LepMass;
     CosthMin = TMath::Sqrt(1 - TMath::Power((fQ3Max / Enu ), 2));
   }
 
   // Compute the maximum xsec value:
   Range1D_t Tl_range ( TMin, TMax ) ;
   Range1D_t ctl_range ( CosthMin, CosthMax ) ;
   double XSecMax = GetXSecMaxTlctl( *interaction, Tl_range, ctl_range ) ;
 
   // -- Generate and Test the Kinematics----------------------------------//
 
   RandomGen * rnd = RandomGen::Instance();
   bool accept = false;
   unsigned int iter = 0;
 
   // loop over different (randomly) selected T and Costh
   while (!accept) {
       iter++;
       if(iter > kRjMaxIterations) {
           // error if try too many times
           LOG("MEC", pWARN)
               << "Couldn't select a valid Tmu, CosTheta pair after "
               << iter << " iterations";
           event->EventFlags()->SetBitNumber(kKineGenErr, true);
           genie::exceptions::EVGThreadException exception;
           exception.SetReason("Couldn't select lepton kinematics");
           exception.SwitchOnFastForward();
           throw exception;
       }
 
       // generate random kinetic energy T and Costh
       T = TMin + (TMax-TMin)*rnd->RndKine().Rndm();
       Costh = CosthMin + (CosthMax-CosthMin)*rnd->RndKine().Rndm();
 
       // Calculate useful values for judging this choice
       genie::utils::mec::Getq0q3FromTlCostl(T, Costh, Enu, LepMass, Q0, Q3);
 
       // Don't bother doing hard work if the selected Q3 is greater than Q3Max
       if (Q3 > fQ3Max) continue ;
 
       Plep = TMath::Sqrt( T * (T + (2.0 * LepMass)));  // ok is sqrt(E2 - m2)
       kinematics->SetKV(kKVTl, T);
       kinematics->SetKV(kKVctl, Costh);
       kinematics->SetKV( kKVQ0, Q0 ) ;
       kinematics->SetKV( kKVQ3, Q3 ) ;
 
       // decide whether to accept or reject these kinematics
       // AND set the chosen two-nucleon system
 
       if (FullDeltaNodelta == 1){
         // this block for the user who wants all CC QE-like 2p2h events
         // We need four different cross sections. Right now, pursue the
         // inelegant method of calling XSec four times - there is
         // definitely some runtime inefficiency here, but it is not awful
 
         // first, get delta-less all
         if (NuPDG > 0) {
           interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNN);
         }
         else {
           interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterPP);
         }
         double XSec = fXSecModel->XSec(interaction, kPSTlctl);
         // now get all with delta
         interaction->ExclTagPtr()->SetResonance(genie::kP33_1232);
         double XSecDelta = fXSecModel->XSec(interaction, kPSTlctl);
         // get PN with delta
         interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNP);
         double XSecDeltaPN = fXSecModel->XSec(interaction, kPSTlctl);
         // now get delta-less PN
         interaction->ExclTagPtr()->SetResonance(genie::kNoResonance);
         double XSecPN = fXSecModel->XSec(interaction, kPSTlctl);
 
         if (XSec > XSecMax) {
           LOG("MEC", pERROR) << "XSec is > XSecMax for nucleus " << TgtPDG << " "
                              << XSec << " > " << XSecMax
                              << " don't let this happen.";
         }
         assert(XSec <= XSecMax);
         accept = XSec > XSecMax*rnd->RndKine().Rndm();
         LOG("MEC", pINFO) << "Xsec, Max, Accept: " << XSec << ", "
                           << XSecMax << ", " << accept;
 
         if(accept){
           // If it passes the All cross section we still need to do two things:
           // * Was the initial state pn or not?
           // * Do we assign the reaction to have had a Delta on the inside?
 
           // PDD means from the part of the XSec with an internal Delta line
           // that (at the diagram level) did not produce a pion in the final state.
 
           bool isPDD = false;
 
           // Find out if we should use a pn initial state
           double myrand = rnd->RndKine().Rndm();
           double pnFraction = XSecPN / XSec;
           LOG("MEC", pDEBUG) << "Test for pn: xsec_pn = " << XSecPN
                              << "; xsec = " << XSec
                              << "; pn_fraction = " << pnFraction
                              << "; random number val = " << myrand;
 
           if (myrand <= pnFraction) {
             // yes it is, add a PN initial state to event record
             event->AddParticle(kPdgClusterNP, kIStNucleonTarget,
                                1, -1, -1, -1, tempp4, v4);
             interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNP);
 
             // Its a pn, so test for Delta by comparing DeltaPN/PN
             if (rnd->RndKine().Rndm() <= XSecDeltaPN / XSecPN) {
               isPDD = true;
             }
           }
           else {
             // no it is not a PN, add either NN or PP initial state to event record.
             if (NuPDG > 0) {
               event->AddParticle(kPdgClusterNN, kIStNucleonTarget,
                                  1, -1, -1, -1, tempp4, v4);
               interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNN);
             }
             else {
               event->AddParticle(kPdgClusterPP, kIStNucleonTarget,
                                  1, -1, -1, -1, tempp4, v4);
               interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterPP);
               }
             // its not pn, so test for Delta (XSecDelta-XSecDeltaPN)/(XSec-XSecPN)
             // right, both numerator and denominator are total not pn.
             if (rnd->RndKine().Rndm() <=
                 (XSecDelta - XSecDeltaPN) / (XSec - XSecPN)) {
               isPDD = true;
             }
           }
 
           // now test whether we tagged this as a pion event
           // and assign that fact to the Exclusive State tag
           // later, we can query const XclsTag & xcls = interaction->ExclTag()
           if (isPDD){
             interaction->ExclTagPtr()->SetResonance(genie::kP33_1232);
           }
 
 
           } // end if accept
       } // end if delta == 1
 
         /* One can make simpler versions of the above for the
            FullDeltaNodelta == 2 (only delta)
            or
            FullDeltaNodelta == 3 (set Delta FF = 1, lose interference effect).
            but I (Rik) don't see what the use-case is for these, genratorly speaking.
         */
 
   } // end while
 
   // -- finish lepton kinematics
   // If the code got here, then we accepted some kinematics
   // and we can proceed to generate the final state.
 
   // define coordinate system wrt neutrino: z along neutrino, xy perp
 
   // Cos theta gives us z, the rest in xy:
   double PlepZ = Plep * Costh;
   double PlepXY = Plep * TMath::Sqrt(1. - TMath::Power(Costh,2));
 
   // random rotation about unit vector for phi direction
   double phi= 2 * kPi * rnd->RndLep().Rndm();
   // now fill x and y from PlepXY
   double PlepX = PlepXY * TMath::Cos(phi);
   double PlepY = PlepXY * TMath::Sin(phi);
 
   // Rotate lepton momentum vector from the reference frame (x'y'z') where
   // {z':(neutrino direction), z'x':(theta plane)} to the LAB
   TVector3 unit_nudir = event->Probe()->P4()->Vect().Unit();
   TVector3 p3l(PlepX, PlepY, PlepZ);
   p3l.RotateUz(unit_nudir);
 
   // Lepton 4-momentum in LAB
   Elep = TMath::Sqrt(LepMass*LepMass + PlepX*PlepX + PlepY*PlepY + PlepZ*PlepZ);
   TLorentzVector p4l(p3l,Elep);
 
   // Figure out the final-state primary lepton PDG code
   int pdgc = interaction->FSPrimLepton()->PdgCode();
   int momidx = event->ProbePosition();
 
   // -- Store Values ------------------------------------------//
   // -- Interaction: Q2
   Q2 = Q3*Q3 - Q0*Q0;
   double gy = Q0 / Enu;
   double gx = kinematics::Q2YtoX(Enu, 2 * kNucleonMass, Q2, gy);
   double gW = kinematics::XYtoW(Enu, 2 * kNucleonMass, gx, gy);
 
   interaction->KinePtr()->SetQ2(Q2, true);
   interaction->KinePtr()->Sety(gy, true);
   interaction->KinePtr()->Setx(gx, true);
   interaction->KinePtr()->SetW(gW, true);
   interaction->KinePtr()->SetFSLeptonP4(p4l);
   // in later methods
   // will also set the four-momentum and W^2 of the hadron system.
 
   // -- Lepton
   event->AddParticle( pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
 
   // Set the final-state lepton polarization
   utils::SetPrimaryLeptonPolarization( event );
 
   LOG("MEC",pDEBUG) << "~~~ LEPTON DONE ~~~";
 }

void MECGenerator::SelectSuSALeptonKinematics ( GHepRecord * event ) const

private

Definition at line 855 of file MECGenerator.cxx.

References fQ3Max, genie::Interaction::FSPrimLepton(), fSuSAMaxXSecDiffTolerance, fXSecModel, genie::utils::mec::GetMaxXSecTlctl(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::RandomGen::Instance(), genie::ProcessInfo::IsEM(), genie::Interaction::KinePtr(), genie::kIStNucleonTarget, genie::kIStStableFinalState, genie::kKineGenErr, genie::kKVctl, genie::kKVTl, genie::controls::kMinQ2Limit, genie::constants::kNucleonMass, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::constants::kPi, genie::kPSTlctl, genie::kRfLab, genie::controls::kRjMaxIterations, LOG, pDEBUG, pERROR, pFATAL, pINFO, genie::GHepRecord::Probe(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), pWARN, genie::utils::kinematics::Q2(), genie::utils::kinematics::Q2YtoX(), genie::RandomGen::RndKine(), genie::RandomGen::RndLep(), genie::Kinematics::SetFSLeptonP4(), genie::Target::SetHitNucPdg(), genie::Kinematics::SetKV(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), genie::InitialState::TgtPdg(), genie::InitialState::TgtPtr(), genie::GHepParticle::X4(), genie::XSecAlgorithmI::XSec(), and genie::utils::kinematics::XYtoW().

Referenced by ProcessEventRecord().

 {
   // Event Properties
   Interaction* interaction = event->Summary();
   Kinematics* kinematics = interaction->KinePtr();
 
   // Choose the appropriate minimum Q^2 value based on the interaction
   // mode (this is important for EM interactions since the differential
   // cross section blows up as Q^2 --> 0)
   double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit
   if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics
     ::electromagnetic::kMinQ2Limit; // EM limit
 
   LOG("MEC", pDEBUG) << "Q2min = " << Q2min;
 
   double Enu = interaction->InitState().ProbeE( kRfLab );
 
   int NuPDG = interaction->InitState().ProbePdg();
   int TgtPDG = interaction->InitState().TgtPdg();
 
   // Interacton vtx
   TLorentzVector v4( *event->Probe()->X4() );
   TLorentzVector tempp4( 0., 0., 0., 0. );
 
   // Lepton Kinematic Limits
   double Costh = 0.0; // lepton angle
   double CosthMax = 1.0;
   double CosthMin = -1.0;
 
   double T = 0.0;  // lepton kinetic energy
   double TMax = std::numeric_limits<double>::max();
   double TMin = 0.0;
 
   double Plep = 0.0; // lepton 3 momentum
   double Elep = 0.0; // lepton energy
   double LepMass = interaction->FSPrimLepton()->Mass();
 
   double Q0 = 0.0; // energy component of q four vector
   double Q3 = 0.0; // magnitude of transfered 3 momentum
   double Q2 = 0.0; // properly Q^2 (Q squared) - transfered 4 momentum.
 
   // Set lepton KE TMax for for throwing rndm in the accept/reject loop.
   // We can accidentally set it too high, because the xsec will return zero.
   // This way if someone reuses this code, they are not tripped up by it.
   TMax = Enu - LepMass;
 
   // Warn if fQ3Max value is suspect 
   // (i.e. below the SuSA model's rage of validity)
   // This is just to warn users who have swapped out a MEC model
   // without remembering to change fQ3Max.
   if(fQ3Max<1.99){
     LOG("MEC", pWARN)
         << "fQ3 max is smaller than expected for SuSAv2 MEC: "
         << fQ3Max << ". Are you sure this is correct?";
   }
 
   // TODO: double-check the limits below
 
   // Set Tmin for throwing rndm in the accept/reject loop
   // the hadron tensors we expect will be limited in q3
   // therefore also the outgoing lepton KE can't be too low or costheta too backward
   // make the accept/reject loop more efficient by using Min values.
   if ( Enu < fQ3Max ) {
     TMin = 0;
     CosthMin = -1;
   } else {
     TMin = TMath::Sqrt( TMath::Power(LepMass, 2) + TMath::Power(Enu - fQ3Max, 2) ) - LepMass;
     CosthMin = TMath::Sqrt( 1. - TMath::Power(fQ3Max / Enu, 2) );
   }
 
   // Generate and Test the Kinematics
 
   RandomGen* rnd = RandomGen::Instance();
   bool accept = false;
   unsigned int iter = 0;
   unsigned int maxIter = kRjMaxIterations;
 
   // TODO: revisit this
   // e-scat xsecs blow up close to theta=0, MC methods won't work ...
   if ( NuPDG == 11 ) maxIter *= 100000;
 
   // Scan the accessible phase space to find the maximum differential cross
   // section to throw against
   double XSecMax = utils::mec::GetMaxXSecTlctl( *fXSecModel, *interaction );
 
   // loop over different (randomly) selected T and Costh
   while ( !accept ) {
     ++iter;
     if ( iter > maxIter ) {
       // error if try too many times
       LOG("MEC", pWARN)
         << "Couldn't select a valid Tmu, CosTheta pair after "
         << iter << " iterations";
       event->EventFlags()->SetBitNumber( kKineGenErr, true );
       genie::exceptions::EVGThreadException exception;
       exception.SetReason( "Couldn't select lepton kinematics" );
       exception.SwitchOnFastForward();
       throw exception;
     }
 
     // generate random kinetic energy T and Costh
     T = TMin + (TMax-TMin)*rnd->RndKine().Rndm();
     Costh = CosthMin + (CosthMax-CosthMin)*rnd->RndKine().Rndm();
 
     // Calculate useful values for judging this choice
     Plep = TMath::Sqrt( T * (T + (2.0 * LepMass)));  // ok is sqrt(E2 - m2)
     Q3 = TMath::Sqrt(Plep*Plep + Enu*Enu - 2.0 * Plep * Enu * Costh);
 
     // TODO: implement this more cleanly (throw Costh from restricted range)
     Q0 = Enu - (T + LepMass);
     Q2 = Q3*Q3 - Q0*Q0;
 
     LOG("MEC", pDEBUG) << "T = " << T << ", Costh = " << Costh
       << ", Q2 = " << Q2;
 
     // Don't bother doing hard work if the selected Q3 is greater than Q3Max
     // or if Q2 falls below the minimum allowed Q^2 value
     if ( Q3 < fQ3Max && Q2 >= Q2min ) {
 
       kinematics->SetKV(kKVTl, T);
       kinematics->SetKV(kKVctl, Costh);
 
       // decide whether to accept or reject these kinematics
       // AND set the chosen two-nucleon system
 
       LOG("MEC", pDEBUG) << " T, Costh: " << T << ", " << Costh ;
 
       // Get total xsec (nn+np)
       double XSec = fXSecModel->XSec( interaction, kPSTlctl );
 
       if ( XSec > XSecMax ) {
         LOG("MEC", pERROR) << "XSec is > XSecMax for nucleus " << TgtPDG << " "
           << XSec << " > " << XSecMax << " don't let this happen.";
 
         double percent_deviation = 200. * ( XSec - XSecMax ) / ( XSecMax + XSec );
 
         if ( percent_deviation > fSuSAMaxXSecDiffTolerance ) {
           LOG( "Kinematics", pFATAL ) << "xsec: (curr) = " << XSec
             << " > (max) = " << XSecMax << "\n for " << *interaction;
           LOG( "Kinematics", pFATAL )
              << "*** Exceeding estimated maximum differential cross section";
           std::terminate();
         }
         else {
           LOG( "Kinematics", pWARN ) << "xsec: (curr) = " << XSec
             << " > (max) = " << XSecMax << "\n for " << *interaction;
           LOG("Kinematics", pWARN) << "*** The fractional deviation of "
             << percent_deviation << " % was allowed";
         }
       }
 
       accept = XSec > XSecMax*rnd->RndKine().Rndm();
       LOG("MEC", pINFO) << "Xsec, Max, Accept: " << XSec << ", "
         << XSecMax << ", " << accept;
 
       if ( accept ) {
         // Now that we've selected kinematics, we also need to choose the
         // isospin of the initial hit nucleon pair
 
         // Find out if we should use a pn initial state
         double myrand_pn = rnd->RndKine().Rndm();
         double pnFraction = dynamic_cast< const SuSAv2MECPXSec* >( fXSecModel )
           ->PairRatio( interaction );
 
         LOG("MEC", pINFO) << "Test for pn: "
           << "; xsec = " << XSec << "; pn_fraction = " << pnFraction
           << "; random number val = " << myrand_pn;
 
         double myrand_pp = rnd->RndKine().Rndm();
         double ppFraction = 0 ;
 
         if ( interaction->ProcInfo().IsEM() ) {
       // calculate ppFraction in the EM case
           ppFraction = dynamic_cast< const SuSAv2MECPXSec* >( fXSecModel )
           ->PairRatio( interaction ,"ppFraction");
 
       LOG("MEC", pINFO) << "Test for pp: "
                         << "; xsec = " << XSec << "; pp_fraction = " << ppFraction
                         << "; random number val = " << myrand_pp;
         }
 
         if ( myrand_pn <= pnFraction ) {
           // yes it is, add a PN initial state to event record
           event->AddParticle(kPdgClusterNP, kIStNucleonTarget,
             1, -1, -1, -1, tempp4, v4);
           interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg( kPdgClusterNP );
         }
         else {
           // no it is not a PN, add either NN or PP initial state to event record (EM case).
           if ( interaction->ProcInfo().IsEM() ) {
             if ( myrand_pp <= ppFraction/(1. - pnFraction) ) {
               // record a PP pair:
               event->AddParticle(kPdgClusterPP, kIStNucleonTarget,
                                  1, -1, -1, -1, tempp4, v4);
               interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg( kPdgClusterPP );
             } else {
               // record a NN pair:
               event->AddParticle(kPdgClusterNN, kIStNucleonTarget,
                                  1, -1, -1, -1, tempp4, v4);
               interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg( kPdgClusterNN );
             }
           } else {
             // no it is not a PN, add either NN or PP initial state to event record (CC cases).
             if ( NuPDG > 0 ) {
               event->AddParticle(kPdgClusterNN, kIStNucleonTarget,
                                  1, -1, -1, -1, tempp4, v4);
               interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg( kPdgClusterNN );
             }
             else {
               event->AddParticle(kPdgClusterPP, kIStNucleonTarget,
                                  1, -1, -1, -1, tempp4, v4);
               interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg( kPdgClusterPP );
             }
           }
         }
       } // end if accept
     } // end if passes q3 test
   } // end while
 
   // -- finish lepton kinematics
   // If the code got here, then we accepted some kinematics
   // and we can proceed to generate the final state.
 
   // define coordinate system wrt neutrino: z along neutrino, xy perp
 
   // Cos theta gives us z, the rest in xy:
   double PlepZ = Plep * Costh;
   double PlepXY = Plep * TMath::Sqrt( 1. - TMath::Power(Costh,2) );
 
   // random rotation about unit vector for phi direction
   double phi = 2. * kPi * rnd->RndLep().Rndm();
   // now fill x and y from PlepXY
   double PlepX = PlepXY * TMath::Cos(phi);
   double PlepY = PlepXY * TMath::Sin(phi);
 
   // Rotate lepton momentum vector from the reference frame (x'y'z') where
   // {z':(neutrino direction), z'x':(theta plane)} to the LAB
   TVector3 unit_nudir = event->Probe()->P4()->Vect().Unit();
   TVector3 p3l( PlepX, PlepY, PlepZ );
   p3l.RotateUz( unit_nudir );
 
   // Lepton 4-momentum in LAB
   Elep = TMath::Sqrt( LepMass*LepMass + PlepX*PlepX + PlepY*PlepY + PlepZ*PlepZ );
   TLorentzVector p4l( p3l, Elep );
 
   // Figure out the final-state primary lepton PDG code
   int pdgc = interaction->FSPrimLepton()->PdgCode();
   int momidx = event->ProbePosition();
 
   // -- Store Values ------------------------------------------//
   // -- Interaction: Q2
   Q0 = Enu - Elep;
   Q2 = Q3*Q3 - Q0*Q0;
   double gy = Q0 / Enu;
   double gx = kinematics::Q2YtoX(Enu, 2 * kNucleonMass, Q2, gy);
   double gW = kinematics::XYtoW(Enu, 2 * kNucleonMass, gx, gy);
 
   interaction->KinePtr()->SetQ2(Q2, true);
   interaction->KinePtr()->Sety(gy, true);
   interaction->KinePtr()->Setx(gx, true);
   interaction->KinePtr()->SetW(gW, true);
   interaction->KinePtr()->SetFSLeptonP4(p4l);
   // in later methods
   // will also set the four-momentum and W^2 of the hadron system.
 
   // -- Lepton
   event->AddParticle( pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4 );
 
   LOG("MEC", pDEBUG) << "~~~ LEPTON DONE ~~~";
 }