genie::MECGenerator Class Reference

Simulate the primary MEC interaction. More...

#include <MECGenerator.h>

Inheritance 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)
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
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 <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

Steve Dytman <dytman+ pitt.edu> Pittsburgh University

Sep. 22, 2008

Copyright (c) 2003-2020, 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.

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

{
// 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 1273 of file MECGenerator.cxx.

{
    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 1279 of file MECGenerator.cxx.

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

void MECGenerator::DecayNucleonCluster ( GHepRecord *  event ) const

private

accept_decay

Definition at line 436 of file MECGenerator.cxx.

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

{
// 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 1078 of file MECGenerator.cxx.

{
    // 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 1306 of file MECGenerator.cxx.

                                                                          {

  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 1285 of file MECGenerator.cxx.

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

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

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

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

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

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

  // 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 845 of file MECGenerator.cxx.

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

  // 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 = 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;

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

Member Data Documentation

int genie::MECGenerator::fFunctionCalls

private

Definition at line 75 of file MECGenerator.h.

int genie::MECGenerator::fMinScanPointsCosth

private

Definition at line 78 of file MECGenerator.h.

int genie::MECGenerator::fMinScanPointsTmu

private

Definition at line 77 of file MECGenerator.h.

const NuclearModelI* genie::MECGenerator::fNuclModel

private

Definition at line 72 of file MECGenerator.h.

TGenPhaseSpace genie::MECGenerator::fPhaseSpaceGenerator

mutableprivate

Definition at line 71 of file MECGenerator.h.

double genie::MECGenerator::fQ3Max

private

Definition at line 80 of file MECGenerator.h.

double genie::MECGenerator::fRelTolerance

private

Definition at line 76 of file MECGenerator.h.

double genie::MECGenerator::fSafetyFactor

private

Definition at line 74 of file MECGenerator.h.

double genie::MECGenerator::fSuSAMaxXSecDiffTolerance

private

Definition at line 84 of file MECGenerator.h.

const XSecAlgorithmI* genie::MECGenerator::fXSecModel

mutableprivate

Definition at line 70 of file MECGenerator.h.

---

The documentation for this class was generated from the following files:

• Generator/src/Physics/Multinucleon/EventGen/MECGenerator.h
• Generator/src/Physics/Multinucleon/EventGen/MECGenerator.cxx