Generates values for the kinematic variables describing coherent neutrino-nucleus pion production events. Is a concrete implementation of the EventRecordVisitorI interface. More...

#include <COHKinematicsGenerator.h>

Inheritance diagram for genie::COHKinematicsGenerator:

[legend]

Collaboration diagram for genie::COHKinematicsGenerator:

[legend]

Public Member Functions
	COHKinematicsGenerator ()
	COHKinematicsGenerator (string config)
	~COHKinematicsGenerator ()
void	ProcessEventRecord (GHepRecord *event_rec) const
void	Configure (const Registry &config)
void	Configure (string config)
void	LoadConfig (void)
void	CalculateKin_ReinSehgal (GHepRecord *event_rec) const
void	CalculateKin_BergerSehgal (GHepRecord *event_rec) const
void	CalculateKin_BergerSehgalFM (GHepRecord *event_rec) const
void	CalculateKin_AlvarezRuso (GHepRecord *event_rec) const
void	SetKinematics (const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction interaction, Kinematics kinematics) const
bool	CheckKinematics (const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction) const
double	ComputeMaxXSec (const Interaction *in) const
double	MaxXSec_ReinSehgal (const Interaction *in) const
double	MaxXSec_BergerSehgal (const Interaction *in) const
double	MaxXSec_BergerSehgalFM (const Interaction *in) const
double	MaxXSec_AlvarezRuso (const Interaction *in) const
double	Energy (const Interaction *in) const
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.
virtual AlgStatus_t	GetStatus (void) const
	Get algorithm status.
virtual bool	AllowReconfig (void) const
virtual AlgCmp_t	Compare (const Algorithm *alg) const
	Compare with input algorithm.
virtual void	SetId (const AlgId &id)
	Set algorithm ID.
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.

Public Attributes
TF2 *	fEnvelope
	2-D envelope used for importance sampling
double	fRo
	nuclear scale parameter

Private Member Functions
double	pionMass (const Interaction *in) const
void	throwOnTooManyIterations (unsigned int iters, GHepRecord *evrec) const

Private Attributes
double	fQ2Min
	lower bound of integration for Q^2 in Berger-Sehgal Model
double	fQ2Max
	upper bound of integration for Q^2 in Berger-Sehgal Model
double	fTMax
	upper bound for t = (q - p_pi)^2

Additional Inherited Members
Static Public Member Functions inherited from genie::Algorithm
static string	BuildParamVectKey (const std::string &comm_name, unsigned int i)
static string	BuildParamVectSizeKey (const std::string &comm_name)
static string	BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
static string	BuildParamMatRowSizeKey (const std::string &comm_name)
static string	BuildParamMatColSizeKey (const std::string &comm_name)
Protected Member Functions inherited from genie::KineGeneratorWithCache
	KineGeneratorWithCache ()
	KineGeneratorWithCache (string name)
	KineGeneratorWithCache (string name, string config)
	~KineGeneratorWithCache ()
virtual double	ComputeMaxXSec (const Interaction *in, const int nkey) const
virtual double	MaxXSec (GHepRecord *evrec, const int nkey=0) const
virtual double	FindMaxXSec (const Interaction *in, const int nkey=0) const
virtual void	CacheMaxXSec (const Interaction *in, double xsec, const int nkey=0) const
virtual CacheBranchFx *	AccessCacheBranch (const Interaction *in, const int nkey=0) const
virtual void	AssertXSecLimits (const Interaction *in, double xsec, double xsec_max) const
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.
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.
int	GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
template<class T>
int	GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
	Handle to load matrix of parameters.
template<class T>
int	GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
int	GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
int	AddTopRegistry (Registry *rp, bool owns=true)
	add registry with top priority, also update ownership
int	AddLowRegistry (Registry *rp, bool owns=true)
	add registry with lowest priority, also update ownership
int	MergeTopRegistry (const Registry &r)
int	AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
	Add registries with top priority, also udated Ownerships.
Protected Attributes inherited from genie::KineGeneratorWithCache
const XSecAlgorithmI *	fXSecModel
double	fSafetyFactor
	ComputeMaxXSec -> ComputeMaxXSec * fSafetyFactor.
std::vector< double >	vSafetyFactors
	MaxXSec -> MaxXSec * fSafetyFactors[nkey].
int	fNumOfSafetyFactors
	Number of given safety factors.
std::vector< string >	vInterpolatorTypes
	Type of interpolator for each key in a branch.
int	fNumOfInterpolatorTypes
	Number of given interpolators types.
double	fMaxXSecDiffTolerance
	max{100(xsec-maxxsec)/.5(xsec+maxxsec)} if xsec>maxxsec
double	fEMin
	min E for which maxxsec is cached - forcing explicit calc.
bool	fGenerateUniformly
	uniform over allowed phase space + event weight?
Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig
bool	fOwnsSubstruc
	true if it owns its substructure (sub-algs,...)
AlgId	fID
	algorithm name and configuration set
vector< Registry * >	fConfVect
vector< bool >	fOwnerships
	ownership for every registry in fConfVect
AlgStatus_t	fStatus
	algorithm execution status
AlgMap *	fOwnedSubAlgMp
	local pool for owned sub-algs (taken out of the factory pool)

Detailed Description

Generates values for the kinematic variables describing coherent neutrino-nucleus pion production events. Is a concrete implementation of the EventRecordVisitorI interface.

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

Created:\n October 03, 2004

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

Definition at line 30 of file COHKinematicsGenerator.h.

Constructor & Destructor Documentation

◆ COHKinematicsGenerator() [1/2]

COHKinematicsGenerator::COHKinematicsGenerator ( )

Definition at line 42 of file COHKinematicsGenerator.cxx.

                                               :
  KineGeneratorWithCache("genie::COHKinematicsGenerator")
{
  fEnvelope = 0;
}

References fEnvelope, and genie::KineGeneratorWithCache::KineGeneratorWithCache().

◆ COHKinematicsGenerator() [2/2]

COHKinematicsGenerator::COHKinematicsGenerator ( string config )

Definition at line 48 of file COHKinematicsGenerator.cxx.

                                                            :
  KineGeneratorWithCache("genie::COHKinematicsGenerator", config)
{
  fEnvelope = 0;
}

References fEnvelope, and genie::KineGeneratorWithCache::KineGeneratorWithCache().

◆ ~COHKinematicsGenerator()

COHKinematicsGenerator::~COHKinematicsGenerator ( )

Definition at line 54 of file COHKinematicsGenerator.cxx.

{
  if(fEnvelope) delete fEnvelope;
}

References fEnvelope.

Member Function Documentation

◆ CalculateKin_AlvarezRuso()

void COHKinematicsGenerator::CalculateKin_AlvarezRuso ( GHepRecord * event_rec ) const

Definition at line 442 of file COHKinematicsGenerator.cxx.

{
 
  LOG("COHKinematics", pNOTICE) << "Using AlvarezRuso Model";
  // Get the Primary Interacton object
  Interaction * interaction = evrec->Summary();
  interaction->SetBit(kISkipProcessChk);
  interaction->SetBit(kISkipKinematicChk);
 
  // Initialise a random number generator
  RandomGen * rnd = RandomGen::Instance();
 
  //-- For the subsequent kinematic selection with the rejection method:
  //   Calculate the max differential cross section or retrieve it from the
  //   cache. Throw an exception and quit the evg thread if a non-positive
  //   value is found.
  //   If the kinematics are generated uniformly over the allowed phase
  //   space the max xsec is irrelevant
  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
 
  //Set up limits of integration variables
  // Primary lepton energy
  const double E_l_min = interaction->FSPrimLepton()->Mass();
  const double E_l_max = interaction->InitStatePtr()->GetProbeP4(kRfLab)->E() - kPionMass;
  // Primary lepton angle with respect to the beam axis
  const double ctheta_l_min = 0.4;
  const double ctheta_l_max = 1.0 - kASmallNum;
  // Pion angle with respect to the beam axis
  const double ctheta_pi_min = 0.4;
  const double ctheta_pi_max = 1.0 - kASmallNum;
  // Pion angle transverse to the beam axis
  const double phi_min = kASmallNum;
  const double phi_max = (2.0 * kPi) - kASmallNum;
  //
  const double d_E_l = E_l_max - E_l_min;
  const double d_ctheta_l  = ctheta_l_max  - ctheta_l_min;
  const double d_ctheta_pi = ctheta_pi_max - ctheta_pi_min;
  const double d_phi = phi_max - phi_min;
 
  //------ Try to select a valid set of kinematics
  unsigned int iter = 0;
  bool accept=false;
  double xsec=-1, g_E_l=-1, g_theta_l=-1, g_phi_l=-1, g_theta_pi=-1, g_phi_pi=-1;
  double g_ctheta_l, g_ctheta_pi;
 
  while(1) {
    iter++;
    if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
 
    //Select kinematic point
    g_E_l = E_l_min + d_E_l * rnd->RndKine().Rndm();
    g_ctheta_l  = ctheta_l_min  + d_ctheta_l  * rnd->RndKine().Rndm();
    g_ctheta_pi = ctheta_pi_min + d_ctheta_pi * rnd->RndKine().Rndm();
    g_phi_l = phi_min + d_phi * rnd->RndKine().Rndm();
    // random phi is relative to phi_l
    g_phi_pi = g_phi_l + (phi_min + d_phi * rnd->RndKine().Rndm());
    g_theta_l = TMath::ACos(g_ctheta_l);
    g_theta_pi = TMath::ACos(g_ctheta_pi);
 
    LOG("COHKinematics", pINFO) << "Trying: Lep(" <<g_E_l << ", " <<
      g_theta_l << ", " << g_phi_l << ") Pi(" <<
      g_theta_pi << ",     " << g_phi_pi << ")";
 
    this->SetKinematics(g_E_l, g_theta_l, g_phi_l, g_theta_pi, g_phi_pi,
                        interaction, interaction->KinePtr());
 
    // computing cross section for the current kinematics
    xsec = fXSecModel->XSec(interaction,kPSElOlOpifE) / (1E-38 * units::cm2);
 
    if (!fGenerateUniformly) {
      //-- decide whether to accept the current kinematics
      double t   = xsec_max * rnd->RndKine().Rndm();
 
      LOG("COHKinematics", pINFO) << "Got: xsec = " << xsec << ", t = " <<
        t << " (max_xsec = " << xsec_max << ")";
 
      this->AssertXSecLimits(interaction, xsec, xsec_max);
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
      LOG("COHKinematics", pDEBUG)
        << "xsec= " << xsec << ", J= 1, Rnd= " << t;
#endif
      accept = (t<xsec);
    }
    else {
      accept = (xsec>0);
    }
 
    //-- If the generated kinematics are accepted, finish-up module's job
    if(accept) {
      LOG("COHKinematics", pNOTICE) << "Selected: Lepton(" <<
        g_E_l << ", " << g_theta_l << ", " <<
        g_phi_l << ") Pion(" << g_theta_pi << ", " << g_phi_pi << ")";
 
      double E_l = g_E_l;
      double theta_l = g_theta_l;
      double theta_pi = g_theta_pi;
      double phi_l = g_phi_l;
      double phi_pi = g_phi_pi;
      const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
      double E_nu       = P4_nu.E();
      double E_pi= E_nu-E_l;
      double m_l = interaction->FSPrimLepton()->Mass();
      double m_pi = this->pionMass(interaction);
 
      double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
      TVector3 lepton_3vector = TVector3(0,0,0);
      lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
      TLorentzVector P4_lep    = TLorentzVector(lepton_3vector , E_l );
 
      double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
      TVector3 pion_3vector = TVector3(0,0,0);
      pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
      TLorentzVector P4_pion   = TLorentzVector(pion_3vector   , E_pi);
 
      TLorentzVector q = P4_nu - P4_lep;
      double Q2 = -q.Mag2();
      double x = Q2/(2*E_pi*constants::kNucleonMass);
      double y = E_pi/E_nu;
 
      double t = TMath::Abs( (q - P4_pion).Mag2() );
 
      // for uniform kinematics, compute an event weight as
      // wght = (phase space volume)*(differential xsec)/(event total xsec)
      if(fGenerateUniformly) {
        // Phase space volume needs checking
        double vol     = d_E_l*d_ctheta_l*d_phi*d_ctheta_pi*d_phi;
        double totxsec = evrec->XSec();
        double wght    = (vol/totxsec)*xsec;
        LOG("COHKinematics", pNOTICE)  << "Kinematics wght = "<< wght;
 
        // apply computed weight to the current event weight
        wght *= evrec->Weight();
        LOG("COHKinematics", pNOTICE) << "Current event wght = " << wght;
        evrec->SetWeight(wght);
      }
 
      // reset bits
      interaction->ResetBit(kISkipProcessChk);
      interaction->ResetBit(kISkipKinematicChk);
      // lock selected kinematics & clear running values
      interaction->KinePtr()->Setx(x, true);
      interaction->KinePtr()->Sety(y, true);
      interaction->KinePtr()->Sett(t, true);
      interaction->KinePtr()->SetW(kPionMass, true);
      interaction->KinePtr()->SetQ2(2*kNucleonMass*x*y*E_nu, true);
      interaction->KinePtr()->ClearRunningValues();
      // set the cross section for the selected kinematics
      evrec->SetDiffXSec(xsec,kPSElOlOpifE);
      return;
    }
  }//while
}

References genie::KineGeneratorWithCache::AssertXSecLimits(), genie::Kinematics::ClearRunningValues(), genie::units::cm2, genie::KineGeneratorWithCache::fGenerateUniformly, genie::Interaction::FSPrimLepton(), genie::KineGeneratorWithCache::fXSecModel, genie::InitialState::GetProbeP4(), genie::Interaction::InitStatePtr(), genie::RandomGen::Instance(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::constants::kNucleonMass, genie::constants::kPi, genie::constants::kPionMass, genie::kPSElOlOpifE, genie::kRfLab, genie::controls::kRjMaxIterations, LOG, genie::KineGeneratorWithCache::MaxXSec(), pDEBUG, pINFO, pionMass(), pNOTICE, genie::RandomGen::RndKine(), genie::GHepRecord::SetDiffXSec(), SetKinematics(), genie::Kinematics::SetQ2(), genie::Kinematics::Sett(), genie::Kinematics::SetW(), genie::GHepRecord::SetWeight(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::GHepRecord::Summary(), throwOnTooManyIterations(), genie::GHepRecord::Weight(), and genie::GHepRecord::XSec().

Referenced by ProcessEventRecord().

◆ CalculateKin_BergerSehgal()

void COHKinematicsGenerator::CalculateKin_BergerSehgal ( GHepRecord * event_rec ) const

Definition at line 86 of file COHKinematicsGenerator.cxx.

{
  // Get the Primary Interacton object
  Interaction * interaction = evrec->Summary();
  interaction->SetBit(kISkipProcessChk);
  interaction->SetBit(kISkipKinematicChk);    // TODO: Why turn this off?
 
  // Initialise a random number generator
  RandomGen * rnd = RandomGen::Instance();
 
  //-- For the subsequent kinematic selection with the rejection method:
  //   Calculate the max differential cross section or retrieve it from the
  //   cache. Throw an exception and quit the evg thread if a non-positive
  //   value is found.
  //
  //   TODO: We are not offering the "fGenerateUniformly" option here.
  double xsec_max = this->MaxXSec(evrec);
 
  //-- Get the kinematical limits for the generated x,y
  const KPhaseSpace & kps = interaction->PhaseSpace();
  Range1D_t y  = kps.YLim();
  Range1D_t Q2(fQ2Min,fQ2Max);
  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
 
  const double ymin  = y.min + kASmallNum;
  const double ymax  = y.max - kASmallNum;
  const double dy    = ymax - ymin;
  const double Q2min = Q2.min + kASmallNum;
  const double Q2max = Q2.max - kASmallNum;
  const double dQ2   = Q2max - Q2min;
 
  unsigned int iter = 0;
  bool accept=false;
  double xsec=-1, gy=-1, gQ2=-1;
 
  while(1) {
    iter++;
    if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
 
    //-- Select unweighted kinematics using importance sampling method.
    // TODO: The importance sampling envelope is not used. Currently,
    // we just employ a standard rejection-method approach.
 
    gy  = ymin  + dy  * rnd->RndKine().Rndm();
    gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
 
    LOG("COHKinematics", pINFO) <<
      "Trying: Q^2 = " << gQ2 << ", y = " << gy; /* << ", t = " << gt; */
 
    interaction->KinePtr()->Sety(gy);
    interaction->KinePtr()->SetQ2(gQ2);
    kinematics::UpdateXFromQ2Y(interaction);
 
    // computing cross section for the current kinematics
    xsec = fXSecModel->XSec(interaction, kPSQ2yfE);
 
    //-- decide whether to accept the current kinematics
    accept = (xsec_max * rnd->RndKine().Rndm() < xsec);
 
    //-- If the generated kinematics are accepted, finish-up module's job
    if(accept) {
      LOG("COHKinematics", pNOTICE)
        << "Selected: Q^2 = " << gQ2 << ", y = " << gy; /* << ", t = " << gt; */
 
      // the Berger-Sehgal COH cross section should be a triple differential cross section
      // d^2xsec/dQ2dydt where t is the the square of the 4p transfer to the
      // nucleus. The cross section used for kinematical selection should have
      // the t-dependence integrated out. The t-dependence is of the form
      // ~exp(-bt). Now that the x,y kinematical variables have been selected
      // we can generate a t using the t-dependence as a PDF.
      const InitialState & init_state = interaction->InitState();
      double Ev    = init_state.ProbeE(kRfLab);
      double Epi   = gy*Ev; // pion energy
      double Epi2  = TMath::Power(Epi,2);
      double pme2  = kPionMass2/Epi2;
      double gx    = interaction->KinePtr()->x(); // utils::kinematics::Q2YtoX(Ev,kNucleonMass,gQ2,gy); // gQ2 / ( 2. * gy * kNucleonMass * Ev);
      double xME   = kNucleonMass*gx/Epi;
      double tA    = 1. + xME - 0.5*pme2;
      /* Range1D_t tl = kps.TLim();   // this becomes the bounds for t */
      double tB    = TMath::Sqrt(1.+ 2*xME) * TMath::Sqrt(1.-pme2);
      double tmin  = 2*Epi2 * (tA-tB);
      double tmax  = 2*Epi2 * (tA+tB);
      double A     = (double) init_state.Tgt().A();
      double A13   = TMath::Power(A,1./3.);
      double R     = fRo * A13 * units::fermi; // nuclear radius
      double R2    = TMath::Power(R,2.);
      double b     = 0.33333 * R2;
      double tsum  = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
      double rt    = tsum * rnd->RndKine().Rndm();
      double gt    = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/b;
 
      // TODO: If we re-install the fGenerateUniformly option, we
      // would compute the event weight here.
 
      // reset bits
      interaction->ResetBit(kISkipProcessChk);
      interaction->ResetBit(kISkipKinematicChk);
 
      // lock selected kinematics & clear running values
      interaction->KinePtr()->Setx(gx, true);
      interaction->KinePtr()->Sety(gy, true);
      interaction->KinePtr()->Sett(gt, true);
      interaction->KinePtr()->SetW(this->pionMass(interaction), true);
      interaction->KinePtr()->SetQ2(gQ2, true);
      interaction->KinePtr()->ClearRunningValues();
 
      // set the cross section for the selected kinematics
      evrec->SetDiffXSec(xsec * TMath::Exp(-b * gt) / tsum, kPSQ2yfE);
 
      return;
    }
  }// iterations
}

References genie::Target::A(), genie::Kinematics::ClearRunningValues(), genie::units::fermi, fQ2Max, fQ2Min, fRo, genie::KineGeneratorWithCache::fXSecModel, gQ2, genie::Interaction::InitState(), genie::RandomGen::Instance(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::constants::kNucleonMass, genie::constants::kPionMass2, genie::kPSQ2yfE, genie::kRfLab, genie::controls::kRjMaxIterations, LOG, genie::Range1D_t::max, genie::KineGeneratorWithCache::MaxXSec(), genie::Range1D_t::min, genie::Interaction::PhaseSpace(), pINFO, pionMass(), pNOTICE, genie::InitialState::ProbeE(), genie::RandomGen::RndKine(), genie::GHepRecord::SetDiffXSec(), genie::Kinematics::SetQ2(), genie::Kinematics::Sett(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::GHepRecord::Summary(), genie::InitialState::Tgt(), throwOnTooManyIterations(), genie::utils::kinematics::UpdateXFromQ2Y(), genie::Kinematics::x(), and genie::KPhaseSpace::YLim().

Referenced by ProcessEventRecord().

◆ CalculateKin_BergerSehgalFM()

void COHKinematicsGenerator::CalculateKin_BergerSehgalFM ( GHepRecord * event_rec ) const

Definition at line 200 of file COHKinematicsGenerator.cxx.

{
  // Get the Primary Interacton object
  Interaction * interaction = evrec->Summary();
  interaction->SetBit(kISkipProcessChk);
  interaction->SetBit(kISkipKinematicChk);
 
  // Initialise a random number generator
  RandomGen * rnd = RandomGen::Instance();
 
  //-- For the subsequent kinematic selection with the rejection method:
  //   Calculate the max differential cross section or retrieve it from the
  //   cache. Throw an exception and quit the evg thread if a non-positive
  //   value is found.
  //
  //   TODO: We are not offering the "fGenerateUniformly" option here.
  double xsec_max = this->MaxXSec(evrec);
 
  //-- Get the kinematical limits for the generated x,y
  const KPhaseSpace & kps = interaction->PhaseSpace();
  Range1D_t y  = kps.YLim();
  Range1D_t Q2(fQ2Min,fQ2Max);
  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
 
  const double ymin  = y.min + kASmallNum;
  const double ymax  = y.max - kASmallNum;
  const double dy    = ymax - ymin;
  const double Q2min = Q2.min + kASmallNum;
  const double Q2max = Q2.max - kASmallNum;
  const double dQ2   = Q2max - Q2min;
  const double tmin  = kASmallNum;
  const double tmax  = fTMax - kASmallNum; // TODO: Choose realistic t bounds
  const double dt    = tmax - tmin;
 
  //-- Try to select a valid (Q^2,y,t) triple.
 
  unsigned int iter = 0;
  bool accept=false;
  double xsec=-1, gy=-1, gt=-1, gQ2=-1;
 
  while(1) {
    iter++;
    if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
 
    //-- Select unweighted kinematics using importance sampling method.
    // TODO: The importance sampling envelope is not used. Currently,
    // we just employ a standard rejection-method approach.
 
    gy  = ymin  + dy  * rnd->RndKine().Rndm();
    gt  = tmin  + dt  * rnd->RndKine().Rndm();
    gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
 
    LOG("COHKinematics", pINFO) <<
      "Trying: Q^2 = " << gQ2 << ", y = " << gy << ", t = " << gt;
 
    interaction->KinePtr()->Sety(gy);
    interaction->KinePtr()->Sett(gt);
    interaction->KinePtr()->SetQ2(gQ2);
 
    // computing cross section for the current kinematics
    xsec = fXSecModel->XSec(interaction, kPSxyfE);
 
    //-- decide whether to accept the current kinematics
    accept = (xsec_max * rnd->RndKine().Rndm() < xsec);
 
    //-- If the generated kinematics are accepted, finish-up module's job
    if(accept) {
      LOG("COHKinematics", pNOTICE)
        << "Selected: Q^2 = " << gQ2 << ", y = " << gy << ", t = " << gt;
 
      // TODO: If we re-install the fGenerateUniformly option, we
      // would compute the event weight here.
 
      // reset bits
      interaction->ResetBit(kISkipProcessChk);
      interaction->ResetBit(kISkipKinematicChk);
 
      // lock selected kinematics & clear running values
      interaction->KinePtr()->SetQ2(gQ2, true);
      interaction->KinePtr()->Sety(gy, true);
      interaction->KinePtr()->Sett(gt, true);
      interaction->KinePtr()->SetW(this->pionMass(interaction), true);
      interaction->KinePtr()->ClearRunningValues();
 
      // set the cross section for the selected kinematics
      evrec->SetDiffXSec(xsec, kPSxytfE);
 
      return;
    }
  }// iterations
}

References genie::Kinematics::ClearRunningValues(), fQ2Max, fQ2Min, fTMax, genie::KineGeneratorWithCache::fXSecModel, gQ2, genie::RandomGen::Instance(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::kPSxyfE, genie::kPSxytfE, genie::controls::kRjMaxIterations, LOG, genie::Range1D_t::max, genie::KineGeneratorWithCache::MaxXSec(), genie::Range1D_t::min, genie::Interaction::PhaseSpace(), pINFO, pionMass(), pNOTICE, genie::RandomGen::RndKine(), genie::GHepRecord::SetDiffXSec(), genie::Kinematics::SetQ2(), genie::Kinematics::Sett(), genie::Kinematics::SetW(), genie::Kinematics::Sety(), genie::GHepRecord::Summary(), throwOnTooManyIterations(), and genie::KPhaseSpace::YLim().

Referenced by ProcessEventRecord().

◆ CalculateKin_ReinSehgal()

void COHKinematicsGenerator::CalculateKin_ReinSehgal ( GHepRecord * event_rec ) const

Definition at line 292 of file COHKinematicsGenerator.cxx.

{
  // Get the Primary Interacton object
  Interaction * interaction = evrec->Summary();
  interaction->SetBit(kISkipProcessChk);
  interaction->SetBit(kISkipKinematicChk);
 
  //-- Get the random number generators
  RandomGen * rnd = RandomGen::Instance();
 
  //-- For the subsequent kinematic selection with the rejection method:
  //   Calculate the max differential cross section or retrieve it from the
  //   cache. Throw an exception and quit the evg thread if a non-positive
  //   value is found.
  //   If the kinematics are generated uniformly over the allowed phase
  //   space the max xsec is irrelevant
  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
 
  //-- Get the kinematical limits for the generated x,y
  const KPhaseSpace & kps = interaction->PhaseSpace();
  Range1D_t y = kps.YLim();
  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
 
  const double xmin = kASmallNum;
  const double xmax = 1.- kASmallNum;
  const double ymin = y.min + kASmallNum;
  const double ymax = y.max - kASmallNum;
  const double dx   = xmax - xmin;
  const double dy   = ymax - ymin;
 
  //------ Try to select a valid x,y pair
 
  unsigned int iter = 0;
  bool accept=false;
  double xsec=-1, gx=-1, gy=-1;
 
  while(1) {
    iter++;
    if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
 
    if(fGenerateUniformly) {
      //-- Generate a x,y pair uniformly in the kinematically allowed range.
      gx = xmin + dx * rnd->RndKine().Rndm();
      gy = ymin + dy * rnd->RndKine().Rndm();
 
    } else {
      //-- Select unweighted kinematics using importance sampling method.
 
      if(iter==1) {
        LOG("COHKinematics", pNOTICE) << "Initializing the sampling envelope";
        double Ev = interaction->InitState().ProbeE(kRfLab);
        fEnvelope->SetRange(xmin,ymin,xmax,ymax);
        fEnvelope->SetParameter(0, xsec_max);
        fEnvelope->SetParameter(1, Ev);
      }
 
      // Generate W,QD2 using the 2-D envelope as PDF
      fEnvelope->GetRandom2(gx,gy);
    }
 
    LOG("COHKinematics", pINFO) << "Trying: x = " << gx << ", y = " << gy;
 
    interaction->KinePtr()->Setx(gx);
    interaction->KinePtr()->Sety(gy);
 
    // computing cross section for the current kinematics
    xsec = fXSecModel->XSec(interaction, kPSxyfE);
 
    //-- decide whether to accept the current kinematics
    if(!fGenerateUniformly) {
      double max = fEnvelope->Eval(gx, gy);
      double t   = max * rnd->RndKine().Rndm();
 
      this->AssertXSecLimits(interaction, xsec, max);
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
      LOG("COHKinematics", pDEBUG)
        << "xsec= " << xsec << ", J= 1, Rnd= " << t;
#endif
      accept = (t<xsec);
    }
    else {
      accept = (xsec>0);
    }
 
    //-- If the generated kinematics are accepted, finish-up module's job
    if(accept) {
      LOG("COHKinematics", pNOTICE) << "Selected: x = "<< gx << ", y = "<< gy;
 
      // the Rein-Sehgal COH cross section should be a triple differential cross section
      // d^2xsec/dxdydt where t is the the square of the 4p transfer to the
      // nucleus. The cross section used for kinematical selection should have
      // the t-dependence integrated out. The t-dependence is of the form
      // ~exp(-bt). Now that the x,y kinematical variables have been selected
      // we can generate a t using the t-dependence as a PDF.
      const InitialState & init_state = interaction->InitState();
      double Ev    = init_state.ProbeE(kRfLab);
      double Epi   = gy*Ev; // pion energy
      double Epi2  = TMath::Power(Epi,2);
      double pme2  = kPionMass2/Epi2;
      double xME   = kNucleonMass*gx/Epi;
      double tA    = 1. + xME - 0.5*pme2;
      double tB    = TMath::Sqrt(1.+ 2*xME) * TMath::Sqrt(1.-pme2);
      double tmin  = 2*Epi2 * (tA-tB);
      double tmax  = 2*Epi2 * (tA+tB);
      double A     = (double) init_state.Tgt().A();
      double A13   = TMath::Power(A,1./3.);
      double R     = fRo * A13 * units::fermi; // nuclear radius
      double R2    = TMath::Power(R,2.);
      double b     = 0.33333 * R2;
      double tsum  = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
      double rt    = tsum * rnd->RndKine().Rndm();
      double gt    = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/b;
 
      LOG("COHKinematics", pNOTICE)
        << "Selected: t = "<< gt << ", from ["<< tmin << ", "<< tmax << "]";
 
      // for uniform kinematics, compute an event weight as
      // wght = (phase space volume)*(differential xsec)/(event total xsec)
      if(fGenerateUniformly) {
        double vol     = y.max-y.min; // dx=1, dt: irrelevant
        double totxsec = evrec->XSec();
        double wght    = (vol/totxsec)*xsec;
        LOG("COHKinematics", pNOTICE)  << "Kinematics wght = "<< wght;
 
        // apply computed weight to the current event weight
        wght *= evrec->Weight();
        LOG("COHKinematics", pNOTICE) << "Current event wght = " << wght;
        evrec->SetWeight(wght);
      }
 
      // reset bits
      interaction->ResetBit(kISkipProcessChk);
      interaction->ResetBit(kISkipKinematicChk);
 
      // lock selected kinematics & clear running values
      interaction->KinePtr()->Setx(gx, true);
      interaction->KinePtr()->Sety(gy, true);
      interaction->KinePtr()->Sett(gt, true);
      interaction->KinePtr()->SetW(kPionMass, true);
      interaction->KinePtr()->SetQ2(2*kNucleonMass*gx*gy*Ev, true);
      interaction->KinePtr()->ClearRunningValues();
 
      // set the cross section for the selected kinematics
      evrec->SetDiffXSec(xsec * TMath::Exp(-b * gt) / tsum, kPSxytfE);
 
      return;
    }
  }// iterations
}

References genie::Target::A(), genie::KineGeneratorWithCache::AssertXSecLimits(), genie::Kinematics::ClearRunningValues(), fEnvelope, genie::units::fermi, genie::KineGeneratorWithCache::fGenerateUniformly, fRo, genie::KineGeneratorWithCache::fXSecModel, genie::Interaction::InitState(), genie::RandomGen::Instance(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::constants::kNucleonMass, genie::constants::kPionMass, genie::constants::kPionMass2, genie::kPSxyfE, genie::kPSxytfE, genie::kRfLab, genie::controls::kRjMaxIterations, LOG, genie::Range1D_t::max, genie::KineGeneratorWithCache::MaxXSec(), genie::Range1D_t::min, pDEBUG, genie::Interaction::PhaseSpace(), pINFO, pNOTICE, genie::InitialState::ProbeE(), genie::RandomGen::RndKine(), genie::GHepRecord::SetDiffXSec(), genie::Kinematics::SetQ2(), genie::Kinematics::Sett(), genie::Kinematics::SetW(), genie::GHepRecord::SetWeight(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::GHepRecord::Summary(), genie::InitialState::Tgt(), throwOnTooManyIterations(), genie::GHepRecord::Weight(), genie::GHepRecord::XSec(), and genie::KPhaseSpace::YLim().

Referenced by ProcessEventRecord().

◆ CheckKinematics()

bool COHKinematicsGenerator::CheckKinematics	(	const double	E_l,
		const double	theta_l,
		const double	phi_l,
		const double	theta_pi,
		const double	phi_pi,
		const Interaction *	interaction ) const

Definition at line 641 of file COHKinematicsGenerator.cxx.

{
  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
  double E_nu       = P4_nu.E();
  double E_pi= E_nu-E_l;
  double m_l = interaction->FSPrimLepton()->Mass();
  double m_pi;
  if ( interaction->ProcInfo().IsWeakCC() ) {
    m_pi = constants::kPionMass;
  }
  else {
    m_pi = constants::kPi0Mass;
  }
  if (E_l <= m_l) {
    return false;
  }
  if (E_pi <= m_pi) {
    return false;
  }
  return true;
}

References genie::Interaction::FSPrimLepton(), genie::Interaction::InitStatePtr(), genie::ProcessInfo::IsWeakCC(), genie::constants::kPi0Mass, genie::constants::kPionMass, genie::kRfLab, and genie::Interaction::ProcInfo().

◆ ComputeMaxXSec()

double COHKinematicsGenerator::ComputeMaxXSec ( const Interaction * in ) const

virtual

Implements genie::KineGeneratorWithCache.

Definition at line 668 of file COHKinematicsGenerator.cxx.

{
  // Computes the maximum differential cross section in the requested phase
  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
  // method and the value is cached at a circular cache branch for retrieval
  // during subsequent event generation.
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  SLOG("COHKinematics", pDEBUG)
    << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
#endif
  double max_xsec = 0.;
  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
    max_xsec = MaxXSec_ReinSehgal(in);
  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015")) {
    max_xsec = MaxXSec_BergerSehgal(in);
  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015")) {
    max_xsec = MaxXSec_BergerSehgalFM(in);
  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
    max_xsec = MaxXSec_AlvarezRuso(in);
  }
  else {
    LOG("COHKinematicsGenerator",pFATAL) <<
      "ComputeMaxXSec >> Cannot calculate max cross-section for " <<
      fXSecModel->Id().Name();
  }
 
  // Apply safety factor, since value retrieved from the cache might
  // correspond to a slightly different energy.
  max_xsec *= fSafetyFactor;
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  SLOG("COHKinematics", pDEBUG) << in->AsString();
  SLOG("COHKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
  SLOG("COHKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
#endif
 
  return max_xsec;
}

References genie::Interaction::AsString(), genie::KineGeneratorWithCache::fSafetyFactor, genie::KineGeneratorWithCache::fXSecModel, LOG, MaxXSec_AlvarezRuso(), MaxXSec_BergerSehgal(), MaxXSec_BergerSehgalFM(), MaxXSec_ReinSehgal(), pDEBUG, pFATAL, and SLOG.

◆ Configure() [1/2]

void COHKinematicsGenerator::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 955 of file COHKinematicsGenerator.cxx.

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

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

◆ Configure() [2/2]

void COHKinematicsGenerator::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 961 of file COHKinematicsGenerator.cxx.

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

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

◆ Energy()

double COHKinematicsGenerator::Energy ( const Interaction * in ) const

virtual

Reimplemented from genie::KineGeneratorWithCache.

Definition at line 921 of file COHKinematicsGenerator.cxx.

{
  // Override the base class Energy() method to cache the max xsec for the
  // neutrino energy in the LAB rather than in the hit nucleon rest frame.
 
  const InitialState & init_state = interaction->InitState();
  double E = init_state.ProbeE(kRfLab);
  return E;
}

References genie::Interaction::InitState(), genie::kRfLab, and genie::InitialState::ProbeE().

◆ LoadConfig()

void COHKinematicsGenerator::LoadConfig ( void )

Definition at line 967 of file COHKinematicsGenerator.cxx.

{
  //-- COH model parameter Ro
  GetParam( "COH-Ro", fRo );
  //-- COH model parameter t_max for t = (q - p_pi)^2
  GetParam( "COH-t-max", fTMax ) ;
  //-- COH model bounds of integration for Q^2
  GetParam( "COH-Q2-min", fQ2Min ) ;
  GetParam( "COH-Q2-max", fQ2Max ) ;
 
  //-- max xsec safety factor (for rejection method) and min cached energy
  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.6 ) ;
  GetParamDef( "Cache-MinEnergy", fEMin,  -1.0 ) ;
 
  //-- Generate kinematics uniformly over allowed phase space and compute
  //   an event weight?
  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
 
  //-- Maximum allowed fractional cross section deviation from maxim cross
  //   section used in rejection method
  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
    assert(fMaxXSecDiffTolerance>=0);
 
  //-- Envelope employed when importance sampling is used
  //   (initialize with dummy range)
  if(fEnvelope) delete fEnvelope;
  fEnvelope = new TF2("CohKinEnvelope",
                      kinematics::COHImportanceSamplingEnvelope,0.,1,0.,1,2);
  // stop ROOT from deleting this object of its own volition
  gROOT->GetListOfFunctions()->Remove(fEnvelope);
}

References genie::utils::kinematics::COHImportanceSamplingEnvelope(), genie::KineGeneratorWithCache::fEMin, fEnvelope, genie::KineGeneratorWithCache::fGenerateUniformly, genie::KineGeneratorWithCache::fMaxXSecDiffTolerance, fQ2Max, fQ2Min, fRo, genie::KineGeneratorWithCache::fSafetyFactor, fTMax, genie::Algorithm::GetParam(), and genie::Algorithm::GetParamDef().

Referenced by Configure(), and Configure().

◆ MaxXSec_AlvarezRuso()

double COHKinematicsGenerator::MaxXSec_AlvarezRuso ( const Interaction * in ) const

Definition at line 851 of file COHKinematicsGenerator.cxx.

{
  // Computes the maximum differential cross section in the requested phase
  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
  // method and the value is cached at a circular cache branch for retrieval
  // during subsequent event generation.
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  SLOG("COHKinematics", pDEBUG)
    << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
#endif
  double max_xsec = 0.;
  double Ev = in->InitState().ProbeE(kRfLab);
 
  const KPhaseSpace & kps = in->PhaseSpace();
  Range1D_t y = kps.YLim();
 
  ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit2");
  gsl::d4Xsec_dEldThetaldOmegapi f(fXSecModel,in);
  f.SetFactor(-1.); // Make it return negative of cross-section so we can minimize
 
  min->SetFunction( f );
  min->SetMaxFunctionCalls(10000);
  min->SetTolerance(0.05);
 
  const double min_el = in->FSPrimLepton()->Mass();
  const double max_el = Ev - kPionMass;
  const unsigned int n_el = 100;
  const double d_el = (max_el - min_el) / double(n_el - 1);
 
  const double min_thetal = kASmallNum;
  const double max_thetal = kPi / 4.0;
  const unsigned int n_thetal = 10;
  const double d_thetal = (max_thetal - min_thetal) / double(n_thetal - 1);
 
  const double min_thetapi = kASmallNum;
  const double max_thetapi = kPi / 2.0;
  const unsigned int n_thetapi = 10;
  const double d_thetapi = (max_thetapi - min_thetapi) / double(n_thetapi - 1);
 
  //~ const double min_phipi = kPi;
  //~ const double max_phipi = 0.5 * kPi;
  const double min_phipi = kASmallNum;
  const double max_phipi = 2*kPi-kASmallNum;
  const unsigned int n_phipi = 10;
  const double d_phipi = (max_phipi - min_phipi) / double(n_phipi - 1);
 
  min->SetLimitedVariable ( 0 ,"E_lep"    , max_el     -kASmallNum , d_el      , min_el     , max_el      );
  min->SetLimitedVariable ( 1 ,"theta_l"  , min_thetal +kASmallNum , d_thetal  , min_thetal , max_thetal  );
  min->SetLimitedVariable ( 2 ,"theta_pi" , min_thetapi+kASmallNum , d_thetapi , min_thetapi, max_thetapi );
  min->SetLimitedVariable ( 3 ,"phi_pi"   , min_phipi  +kASmallNum , d_phipi   , min_phipi  , max_phipi   );
 
  min->Minimize();
  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;
 
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  SLOG("COHKinematics", pDEBUG) << in->AsString();
  SLOG("COHKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
  SLOG("COHKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
#endif
 
  delete min;
 
  return max_xsec;
}

References genie::Interaction::AsString(), genie::KineGeneratorWithCache::fSafetyFactor, genie::Interaction::FSPrimLepton(), genie::KineGeneratorWithCache::fXSecModel, genie::Interaction::InitState(), genie::controls::kASmallNum, genie::constants::kPi, genie::constants::kPionMass, genie::kRfLab, pDEBUG, genie::Interaction::PhaseSpace(), genie::InitialState::ProbeE(), genie::utils::gsl::d4Xsec_dEldThetaldOmegapi::SetFactor(), SLOG, and genie::KPhaseSpace::YLim().

Referenced by ComputeMaxXSec().

◆ MaxXSec_BergerSehgal()

double COHKinematicsGenerator::MaxXSec_BergerSehgal ( const Interaction * in ) const

Definition at line 708 of file COHKinematicsGenerator.cxx.

{
  double max_xsec = 0;
  const int NQ2   = 50;
  const int Ny    = 50;
 
  const KPhaseSpace & kps = in->PhaseSpace();
  Range1D_t Q2r = kps.Q2Lim();
  Q2r.max = fQ2Max;
 
  const double logQ2min = TMath::Log10(Q2r.min + kASmallNum);
  const double logQ2max = TMath::Log10(Q2r.max);
  const double dlogQ2   = (logQ2max - logQ2min) /(NQ2-1);
 
  for(int i=0; i<NQ2; i++) {
    double Q2 = TMath::Power(10, logQ2min + i * dlogQ2);
    in->KinePtr()->SetQ2(Q2);
 
    Range1D_t yr = kps.YLim();
    if ((yr.max < 0) || (yr.max < yr.min) ||
        (yr.max > 1) || (yr.min < 0)) { // forbidden kinematics
      continue;
    }
    const double logymin  = TMath::Log10(yr.min);
    const double logymax  = TMath::Log10(yr.max);
    const double dlogy    = (logymax - logymin) /(Ny-1);
 
    for(int j=0; j<Ny; j++) {
      double gy = TMath::Power(10, logymin + j * dlogy);
      in->KinePtr()->Sety(gy);
 
      /* Range1D_t tl = kps.TLim();   // TESTING! - this becomes a loop over t */
      kinematics::UpdateXFromQ2Y(in);
 
      // Note: We're not stepping through log Q^2, log y - we "unpacked"
      double xsec = fXSecModel->XSec(in, kPSQ2yfE);
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
      LOG("COHKinematics", pDEBUG)
        << "xsec(Q2= " << Q2 << ", y= " << gy << ", t = " << gt << ") = " << xsec;
#endif
      max_xsec = TMath::Max(max_xsec, xsec);
 
    } // y
  } // Q2
  return max_xsec;
}

References fQ2Max, genie::KineGeneratorWithCache::fXSecModel, genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::kPSQ2yfE, LOG, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, genie::Interaction::PhaseSpace(), genie::KPhaseSpace::Q2Lim(), genie::Kinematics::SetQ2(), genie::Kinematics::Sety(), genie::utils::kinematics::UpdateXFromQ2Y(), and genie::KPhaseSpace::YLim().

Referenced by ComputeMaxXSec().

◆ MaxXSec_BergerSehgalFM()

double COHKinematicsGenerator::MaxXSec_BergerSehgalFM ( const Interaction * in ) const

Definition at line 755 of file COHKinematicsGenerator.cxx.

{
  double max_xsec = 0;
  // How many sampling bins in each variable for max xsec calculation?
  const int NQ2   = 50;
  const int Ny    = 50;
  const int Nt    = 50;
 
  const KPhaseSpace & kps = in->PhaseSpace();
  Range1D_t Q2r = kps.Q2Lim();
  Q2r.max = fQ2Max;
 
  const double logQ2min = TMath::Log10(Q2r.min + kASmallNum);
  const double logQ2max = TMath::Log10(Q2r.max);
  const double logtmin  = TMath::Log10(kASmallNum);
  const double logtmax  = TMath::Log10(fTMax - kASmallNum);
  const double dlogQ2   = (logQ2max - logQ2min) /(NQ2-1);
  const double dlogt    = (logtmax - logtmin) /(Nt-1);
 
  for(int i=0; i<NQ2; i++) {
    double Q2 = TMath::Power(10, logQ2min + i * dlogQ2);
    in->KinePtr()->SetQ2(Q2);
 
    Range1D_t yr = kps.YLim();
    if ((yr.max < 0) || (yr.max < yr.min) ||
        (yr.max > 1) || (yr.min < 0)) { // forbidden kinematics
      continue;
    }
    const double logymin  = TMath::Log10(yr.min);
    const double logymax  = TMath::Log10(yr.max);
    const double dlogy    = (logymax - logymin) /(Ny-1);
 
    for(int j=0; j<Ny; j++) {
      double gy = TMath::Power(10, logymin + j * dlogy);
 
      for(int k=0; k<Nt; k++) {
        double gt = TMath::Power(10, logtmin + k * dlogt);
 
        in->KinePtr()->Sety(gy);
        in->KinePtr()->Sett(gt);
 
        double xsec = fXSecModel->XSec(in, kPSxyfE);
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
        LOG("COHKinematics", pDEBUG)
          << "xsec(Q2= " << Q2 << ", y= " << gy << ", t = " << gt << ") = " << xsec;
#endif
        max_xsec = TMath::Max(max_xsec, xsec);
 
      } // t
    } // y
  } // Q2
  return max_xsec;
}

References fQ2Max, fTMax, genie::KineGeneratorWithCache::fXSecModel, genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::kPSxyfE, LOG, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, genie::Interaction::PhaseSpace(), genie::KPhaseSpace::Q2Lim(), genie::Kinematics::SetQ2(), genie::Kinematics::Sett(), genie::Kinematics::Sety(), and genie::KPhaseSpace::YLim().

Referenced by ComputeMaxXSec().

◆ MaxXSec_ReinSehgal()

double COHKinematicsGenerator::MaxXSec_ReinSehgal ( const Interaction * in ) const

Definition at line 809 of file COHKinematicsGenerator.cxx.

{
  double max_xsec = 0;
  double Ev = in->InitState().ProbeE(kRfLab);
 
  const int Nx = 50;
  const int Ny = 50;
 
  const KPhaseSpace & kps = in->PhaseSpace();
  Range1D_t y = kps.YLim();
 
  const double logxmin = TMath::Log10(1E-5);
  const double logxmax = TMath::Log10(1.0);
  const double logymin = TMath::Log10(y.min);
  const double logymax = TMath::Log10(y.max);
 
  const double dlogx   = (logxmax - logxmin) /(Nx-1);
  const double dlogy   = (logymax - logymin) /(Ny-1);
 
  for(int i=0; i<Nx; i++) {
    double gx = TMath::Power(10, logxmin + i * dlogx);
    for(int j=0; j<Ny; j++) {
      double gy = TMath::Power(10, logymin + j * dlogy);
 
      double Q2 = 2*kNucleonMass*gx*gy*Ev;
      if(Ev>1.0 && Q2>0.01) continue;
 
      in->KinePtr()->Setx(gx);
      in->KinePtr()->Sety(gy);
 
      double xsec = fXSecModel->XSec(in, kPSxyfE);
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
      LOG("COHKinematics", pDEBUG)
        << "xsec(x= " << gx << ", y= " << gy << ") = " << xsec;
#endif
      max_xsec = TMath::Max(max_xsec, xsec);
 
    }//y
  }//x
  return max_xsec;
}

References genie::KineGeneratorWithCache::fXSecModel, genie::Interaction::InitState(), genie::Interaction::KinePtr(), genie::constants::kNucleonMass, genie::kPSxyfE, genie::kRfLab, LOG, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, genie::Interaction::PhaseSpace(), genie::InitialState::ProbeE(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), and genie::KPhaseSpace::YLim().

Referenced by ComputeMaxXSec().

◆ pionMass()

double COHKinematicsGenerator::pionMass ( const Interaction * in ) const

private

Definition at line 931 of file COHKinematicsGenerator.cxx.

{
  double m_pi = 0.0;
  if ( in->ProcInfo().IsWeakCC() ) {
    m_pi = constants::kPionMass;
  } else {
    m_pi = constants::kPi0Mass;
  }
  return m_pi;
}

References genie::ProcessInfo::IsWeakCC(), genie::constants::kPi0Mass, genie::constants::kPionMass, and genie::Interaction::ProcInfo().

Referenced by CalculateKin_AlvarezRuso(), CalculateKin_BergerSehgal(), and CalculateKin_BergerSehgalFM().

◆ ProcessEventRecord()

void COHKinematicsGenerator::ProcessEventRecord ( GHepRecord * event_rec ) const

virtual

Implements genie::EventRecordVisitorI.

Definition at line 59 of file COHKinematicsGenerator.cxx.

{
  if(fGenerateUniformly) {
    LOG("COHKinematics", pNOTICE)
      << "Generating kinematics uniformly over the allowed phase space";
  }
 
  //-- Access cross section algorithm for running thread
  RunningThreadInfo * rtinfo = RunningThreadInfo::Instance();
  const EventGeneratorI * evg = rtinfo->RunningThread();
  fXSecModel = evg->CrossSectionAlg();
  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
    CalculateKin_ReinSehgal(evrec);
  } else if (fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015") {
    CalculateKin_BergerSehgal(evrec);
  } else if (fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015") {
    CalculateKin_BergerSehgalFM(evrec);
  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
    CalculateKin_AlvarezRuso(evrec);
  }
  else {
    LOG("COHKinematicsGenerator",pFATAL) <<
      "ProcessEventRecord >> Cannot calculate kinematics for " <<
      fXSecModel->Id().Name();
  }
}

References CalculateKin_AlvarezRuso(), CalculateKin_BergerSehgal(), CalculateKin_BergerSehgalFM(), CalculateKin_ReinSehgal(), genie::EventGeneratorI::CrossSectionAlg(), genie::KineGeneratorWithCache::fGenerateUniformly, genie::KineGeneratorWithCache::fXSecModel, genie::RunningThreadInfo::Instance(), LOG, pFATAL, pNOTICE, and genie::RunningThreadInfo::RunningThread().

◆ SetKinematics()

void COHKinematicsGenerator::SetKinematics	(	const double	E_l,
		const double	theta_l,
		const double	phi_l,
		const double	theta_pi,
		const double	phi_pi,
		const Interaction *	interaction,
		Kinematics *	kinematics ) const

Definition at line 595 of file COHKinematicsGenerator.cxx.

{
  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
  double E_nu       = P4_nu.E();
  double E_pi= E_nu-E_l;
  double m_l = interaction->FSPrimLepton()->Mass();
  double m_pi;
  if ( interaction->ProcInfo().IsWeakCC() ) {
    m_pi = constants::kPionMass;
  } else {
    m_pi = constants::kPi0Mass;
  }
  double p_l=0.0;
  if (E_l > m_l) {
    p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
  }
  TVector3 lepton_3vector = TVector3(0,0,0);
  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
  TLorentzVector P4_lep    = TLorentzVector(lepton_3vector , E_l );
 
  double p_pi=0.0;
  if (E_pi > m_pi) {
    p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
  }
  TVector3 pion_3vector = TVector3(0,0,0);
  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
  TLorentzVector P4_pion   = TLorentzVector(pion_3vector   , E_pi);
 
  double Q2 = -(P4_nu-P4_lep).Mag2();
  double x = Q2/(2*E_pi*constants::kNucleonMass);
  double y = E_pi/E_nu;
 
  kinematics->Setx(x);
  kinematics->Sety(y);
  kinematics::UpdateWQ2FromXY(interaction);
 
  kinematics->SetFSLeptonP4(P4_lep );
  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
}

References genie::Interaction::FSPrimLepton(), genie::Interaction::InitStatePtr(), genie::ProcessInfo::IsWeakCC(), genie::constants::kNucleonMass, genie::constants::kPi0Mass, genie::constants::kPionMass, genie::kRfLab, genie::Interaction::ProcInfo(), and genie::utils::kinematics::UpdateWQ2FromXY().

Referenced by CalculateKin_AlvarezRuso().

◆ throwOnTooManyIterations()

void COHKinematicsGenerator::throwOnTooManyIterations	(	unsigned int	iters,
		GHepRecord *	evrec ) const

private

Definition at line 942 of file COHKinematicsGenerator.cxx.

{
  LOG("COHKinematics", pWARN)
    << "*** Could not select valid kinematics after "
    << iters << " iterations";
  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
  genie::exceptions::EVGThreadException exception;
  exception.SetReason("Couldn't select kinematics");
  exception.SwitchOnFastForward();
  throw exception;
}

References genie::GHepRecord::EventFlags(), genie::kKineGenErr, LOG, pWARN, genie::exceptions::EVGThreadException::SetReason(), and genie::exceptions::EVGThreadException::SwitchOnFastForward().

Referenced by CalculateKin_AlvarezRuso(), CalculateKin_BergerSehgal(), CalculateKin_BergerSehgalFM(), and CalculateKin_ReinSehgal().

Member Data Documentation

◆ fEnvelope

TF2* genie::COHKinematicsGenerator::fEnvelope

mutable

2-D envelope used for importance sampling

Definition at line 71 of file COHKinematicsGenerator.h.

Referenced by CalculateKin_ReinSehgal(), COHKinematicsGenerator(), COHKinematicsGenerator(), LoadConfig(), and ~COHKinematicsGenerator().

◆ fQ2Max

double genie::COHKinematicsGenerator::fQ2Max

private

upper bound of integration for Q^2 in Berger-Sehgal Model

Definition at line 79 of file COHKinematicsGenerator.h.

Referenced by CalculateKin_BergerSehgal(), CalculateKin_BergerSehgalFM(), LoadConfig(), MaxXSec_BergerSehgal(), and MaxXSec_BergerSehgalFM().

◆ fQ2Min

double genie::COHKinematicsGenerator::fQ2Min

private

lower bound of integration for Q^2 in Berger-Sehgal Model

Definition at line 78 of file COHKinematicsGenerator.h.

Referenced by CalculateKin_BergerSehgal(), CalculateKin_BergerSehgalFM(), and LoadConfig().

◆ fRo

double genie::COHKinematicsGenerator::fRo

nuclear scale parameter

Definition at line 72 of file COHKinematicsGenerator.h.

Referenced by CalculateKin_BergerSehgal(), CalculateKin_ReinSehgal(), and LoadConfig().

◆ fTMax

double genie::COHKinematicsGenerator::fTMax

private

upper bound for t = (q - p_pi)^2

Definition at line 80 of file COHKinematicsGenerator.h.

Referenced by CalculateKin_BergerSehgalFM(), LoadConfig(), and MaxXSec_BergerSehgalFM().

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

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ COHKinematicsGenerator() [1/2]

◆ COHKinematicsGenerator() [2/2]

◆ ~COHKinematicsGenerator()

Member Function Documentation

◆ CalculateKin_AlvarezRuso()

◆ CalculateKin_BergerSehgal()

◆ CalculateKin_BergerSehgalFM()

◆ CalculateKin_ReinSehgal()

◆ CheckKinematics()

◆ ComputeMaxXSec()

◆ Configure() [1/2]

◆ Configure() [2/2]

◆ Energy()

◆ LoadConfig()

◆ MaxXSec_AlvarezRuso()

◆ MaxXSec_BergerSehgal()

◆ MaxXSec_BergerSehgalFM()

◆ MaxXSec_ReinSehgal()

◆ pionMass()

◆ ProcessEventRecord()

◆ SetKinematics()

◆ throwOnTooManyIterations()

Member Data Documentation

◆ fEnvelope

◆ fQ2Max

◆ fQ2Min

◆ fRo

◆ fTMax