MECGenerator.cxx

Go to the documentation of this file.

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

 Costas Andreopoulos <c.andreopoulos \at cern.ch>
 University of Liverpool

 Steve Dytman <dytman+ \at pitt.edu>
 Pittsburgh University
*/
//____________________________________________________________________________

#include <TMath.h>
#include "Math/Minimizer.h"
#include "Math/Factory.h"

#include "Framework/Algorithm/AlgConfigPool.h"
#include "Framework/EventGen/XSecAlgorithmI.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Controls.h"
#include "Framework/EventGen/EVGThreadException.h"
#include "Framework/EventGen/RunningThreadInfo.h"
#include "Framework/EventGen/EventGeneratorI.h"
#include "Framework/GHEP/GHepStatus.h"
#include "Framework/GHEP/GHepFlags.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/Messenger/Messenger.h"
#include "Physics/Common/PrimaryLeptonUtils.h"
#include "Physics/Multinucleon/EventGen/MECGenerator.h"
#include "Physics/Multinucleon/XSection/MECUtils.h"
#include "Physics/Multinucleon/XSection/SuSAv2MECPXSec.h"

#include "Physics/NuclearState/NuclearModelI.h"
//#include "Physics/Multinucleon/XSection/MECHadronTensor.h"
#include "Physics/HadronTensors/HadronTensorI.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/KineUtils.h"
#include "Framework/Utils/PrintUtils.h"

using namespace genie;
using namespace genie::utils;
using namespace genie::constants;
using namespace genie::controls;

//___________________________________________________________________________
MECGenerator::MECGenerator() :
EventRecordVisitorI("genie::MECGenerator")
{

}
//___________________________________________________________________________
MECGenerator::MECGenerator(string config) :
EventRecordVisitorI("genie::MECGenerator", config)
{

}
//___________________________________________________________________________
MECGenerator::~MECGenerator()
{

}
//___________________________________________________________________________
void MECGenerator::ProcessEventRecord(GHepRecord * event) const
{
  //-- 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::AddTargetRemnant(GHepRecord * event) const
{
// 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::GenerateFermiMomentum(GHepRecord * event) const
{
// 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::SelectEmpiricalKinematics(GHepRecord * event) const
{
// 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::AddFinalStateLepton(GHepRecord * event) const
{
// 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 MECGenerator::RecoilNucleonCluster(GHepRecord * event) const
{
  // 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::DecayNucleonCluster(GHepRecord * event) const
{
// 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;
}
//___________________________________________________________________________
PDGCodeList MECGenerator::NucleonClusterConstituents(int pdgc) const
{
  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::SelectNSVLeptonKinematics (GHepRecord * event) const
{
  // -- 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
{
  // 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 ~~~";
}
//___________________________________________________________________________
void MECGenerator::GenerateNSVInitialHadrons(GHepRecord * event) const
{
    // 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);
}
//___________________________________________________________________________
void MECGenerator::Configure(const Registry & config)
{
    Algorithm::Configure(config);
    this->LoadConfig();
}
//___________________________________________________________________________
void MECGenerator::Configure(string config)
{
    Algorithm::Configure(config);
    this->LoadConfig();
}
//___________________________________________________________________________
void MECGenerator::LoadConfig(void)
{
    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. );
}
//___________________________________________________________________________
double MECGenerator::GetXSecMaxTlctl( const Interaction & in,
                                      const Range1D_t & Tl_range,
                                      const Range1D_t & ctl_range ) const {

  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 ;
}

//___________________________________________________________________________