genie::KPhaseSpace Class Reference

Kinematical phase space. More...

#include "/cvmfs/nova.opensciencegrid.org/externals/genie/v3_00_06_p01/Linux64bit+3.10-2.17-e17-debug/GENIE-Generator/src/Framework/Interaction/KPhaseSpace.h"

Inheritance diagram for genie::KPhaseSpace:

Public Member Functions
	KPhaseSpace (void)
	KPhaseSpace (const Interaction *in)
	~KPhaseSpace (void)
void	UseInteraction (const Interaction *in)
double	Threshold (void) const
	Energy threshold. More...
bool	IsAboveThreshold (void) const
	Checks whether the interaction is above the energy threshold. More...
bool	IsAllowed (void) const
	Check whether the current kinematics is in the allowed phase space. More...
Range1D_t	Limits (KineVar_t kvar) const
	Return the kinematical variable limits. More...
double	Minimum (KineVar_t kvar) const
double	Maximum (KineVar_t kvar) const
Range1D_t	WLim (void) const
	W limits. More...
Range1D_t	Q2Lim_W (void) const
	Q2 limits @ fixed W. More...
Range1D_t	q2Lim_W (void) const
	q2 limits @ fixed W More...
Range1D_t	Q2Lim (void) const
	Q2 limits. More...
Range1D_t	q2Lim (void) const
	q2 limits More...
Range1D_t	XLim (void) const
	x limits More...
Range1D_t	YLim (void) const
	y limits More...
Range1D_t	YLim_X (void) const
	y limits @ fixed x More...
Range1D_t	YLim (double xsi) const
	y limits (COH) More...
Range1D_t	YLim_X (double xsi) const
	y limits @ fixed x (COH) More...
Range1D_t	TLim (void) const
	t limits More...

Static Public Member Functions
static double	GetTMaxDFR ()

Private Member Functions
void	Init (void)

Private Attributes
const Interaction *	fInteraction

Detailed Description

Kinematical phase space.

Author

Costas Andreopoulos <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

May 06, 2004

Copyright (c) 2003-2019, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 34 of file KPhaseSpace.h.

Constructor & Destructor Documentation

genie::KPhaseSpace::KPhaseSpace

(

void

)

KPhaseSpace::KPhaseSpace

(

const Interaction *

in

)

Definition at line 47 of file KPhaseSpace.cxx.

References UseInteraction().

                                               :
TObject(), fInteraction(NULL)
{
  this->UseInteraction(in);
}

KPhaseSpace::~KPhaseSpace

(

void

)

Definition at line 53 of file KPhaseSpace.cxx.

{

}

Member Function Documentation

double KPhaseSpace::GetTMaxDFR ( )

static

Definition at line 58 of file KPhaseSpace.cxx.

References genie::AlgConfigPool::CommonList(), genie::Registry::GetDouble(), genie::AlgConfigPool::Instance(), and r().

Referenced by genie::DFRKinematicsGenerator::ProcessEventRecord(), and TLim().

{
  static bool tMaxLoaded = false;
  static double DFR_tMax = -1;

  if (!tMaxLoaded)
  {
    AlgConfigPool * confp = AlgConfigPool::Instance();
    const Registry * r = confp->CommonList( "Param", "Diffractive" ) ;
    double tmax = r->GetDouble("DFR-t-max");
    DFR_tMax = tmax;
    tMaxLoaded = true;
  }

  return DFR_tMax;

}

void genie::KPhaseSpace::Init ( void  )

private

bool KPhaseSpace::IsAboveThreshold

(

void

)

const

Checks whether the interaction is above the energy threshold.

Definition at line 225 of file KPhaseSpace.cxx.

References E, fInteraction, genie::Interaction::InitState(), genie::ProcessInfo::IsAMNuGamma(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDarkMatterElastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::ProcessInfo::IsIMDAnnihilation(), genie::ProcessInfo::IsInverseBetaDecay(), genie::ProcessInfo::IsInverseMuDecay(), genie::ProcessInfo::IsMEC(), genie::ProcessInfo::IsNuElectronElastic(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::ProcessInfo::IsSingleKaon(), genie::kRfHitNucRest, genie::kRfLab, LOG, pDEBUG, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), and Threshold().

Referenced by genie::COHXSecAR::Integrate(), genie::IMDXSec::Integrate(), genie::COHXSec::Integrate(), genie::RESXSec::Integrate(), genie::DISXSec::Integrate(), genie::QELXSec::Integrate(), genie::DMDISXSec::Integrate(), genie::DMELXSec::Integrate(), genie::MECXSec::Integrate(), genie::AlamSimoAtharVacasSKXSec::Integrate(), genie::DFRXSec::Integrate(), genie::NuElectronXSec::Integrate(), genie::SmithMonizQELCCXSec::Integrate(), genie::ReinSehgalRESXSec::Integrate(), genie::ReinSehgalRESXSecFast::Integrate(), genie::ReinSehgalSPPXSec::Integrate(), genie::FermiMover::KickHitNucleon(), genie::XSecAlgorithmI::ValidKinematics(), genie::StrumiaVissaniIBDPXSec::ValidKinematics(), and genie::KLVOxygenIBDPXSec::ValidKinematics().

{
  double E    = 0.;
  double Ethr = this->Threshold();

  const ProcessInfo &  pi         = fInteraction->ProcInfo();
  const InitialState & init_state = fInteraction->InitState();

  if (pi.IsCoherent()       ||
      pi.IsInverseMuDecay() ||
      pi.IsIMDAnnihilation() ||
      pi.IsNuElectronElastic() ||
      pi.IsMEC())
  {
      E = init_state.ProbeE(kRfLab);
  }

  if(pi.IsQuasiElastic()            ||
     pi.IsDarkMatterElastic()       ||
     pi.IsInverseBetaDecay()        ||
     pi.IsResonant()                ||
     pi.IsDeepInelastic()           ||
     pi.IsDarkMatterDeepInelastic() ||
     pi.IsDiffractive()             ||
     pi.IsSingleKaon()              ||
     pi.IsAMNuGamma())
  {
      E = init_state.ProbeE(kRfHitNucRest);
  }

  LOG("KPhaseSpace", pDEBUG) << "E = " << E << ", Ethr = " << Ethr;
  return (E>Ethr);
}

bool KPhaseSpace::IsAllowed

(

void

)

const

Check whether the current kinematics is in the allowed phase space.

Definition at line 259 of file KPhaseSpace.cxx.

References fInteraction, genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDarkMatterElastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::ProcessInfo::IsIMDAnnihilation(), genie::ProcessInfo::IsInverseBetaDecay(), genie::ProcessInfo::IsInverseMuDecay(), genie::ProcessInfo::IsMEC(), genie::ProcessInfo::IsNuElectronElastic(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::ProcessInfo::IsSingleKaon(), genie::utils::math::IsWithinLimits(), genie::Interaction::Kine(), LOG, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, python.hepunit::pi, genie::Interaction::ProcInfo(), genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), Q2Lim(), Q2Lim_W(), genie::Kinematics::t(), confusionMatrixTree::t, moon_position_table_new3::tl(), TLim(), genie::utils::kinematics::UpdateWQ2FromXY(), W, genie::Kinematics::W(), WLim(), genie::Kinematics::x(), submit_syst::x, XLim(), genie::Kinematics::y(), submit_syst::y, and YLim().

Referenced by genie::XSecAlgorithmI::ValidKinematics().

{
  const ProcessInfo & pi   = fInteraction->ProcInfo();
  const Kinematics &  kine = fInteraction->Kine();

  // ASK single kaon:
  // XSec code returns zero when kinematics are not allowed
  // Here just let kinematics always be allowed
  if(pi.IsSingleKaon()) {
    return true;
  }

  // QEL:
  //  Check the running Q2 vs the Q2 limits
  if(pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic()) {
    Range1D_t Q2l = this->Q2Lim();
    double    Q2  = kine.Q2();
    bool in_phys = math::IsWithinLimits(Q2, Q2l);
    bool allowed = in_phys;
    return allowed;
  }

  // RES
  //   Check the running W vs the W limits
  //   & the running Q2 vs Q2 limits for the given W
  if(pi.IsResonant()) {
    Range1D_t Wl  = this->WLim();
    Range1D_t Q2l = this->Q2Lim_W();
    double    W   = kine.W();
    double    Q2  = kine.Q2();
    bool in_phys = (math::IsWithinLimits(Q2, Q2l) && math::IsWithinLimits(W, Wl));
    bool allowed = in_phys;
    return allowed;
  }

  // DIS
  if(pi.IsDeepInelastic() || pi.IsDarkMatterDeepInelastic()) {
    Range1D_t Wl  = this->WLim();
    Range1D_t Q2l = this->Q2Lim_W();
    double    W   = kine.W();
    double    Q2  = kine.Q2();
    bool in_phys = (math::IsWithinLimits(Q2, Q2l) && math::IsWithinLimits(W, Wl));
    bool allowed = in_phys;
    return allowed;
  }

  //IMD
  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation() || pi.IsNuElectronElastic()) {
    Range1D_t yl = this->YLim();
    double    y  = kine.y();
    bool in_phys = math::IsWithinLimits(y, yl);
    bool allowed = in_phys;
    return allowed;
  }

  //COH
  if (pi.IsCoherent()) {
    Range1D_t xl = this->XLim();
    Range1D_t yl = this->YLim();
    double    x  = kine.x();
    double    y  = kine.y();
    bool in_phys = (math::IsWithinLimits(x, xl) && math::IsWithinLimits(y, yl));
    bool allowed = in_phys;
    return allowed;
  }

  // DFR
  if (pi.IsDiffractive()) {
    // first two checks are the same as RES & DIS
    Range1D_t Wl  = this->WLim();
    Range1D_t Q2l = this->Q2Lim_W();

    kinematics::UpdateWQ2FromXY(fInteraction);
    double    W   = kine.W();
    double    Q2  = kine.Q2();

    LOG("KPhaseSpace", pDEBUG) << " W = " << W << ", limits = [" << Wl.min << "," << Wl.max << "];";
    LOG("KPhaseSpace", pDEBUG) << " Q2 = " << Q2 << ", limits = [" << Q2l.min << "," << Q2l.max << "];";
    bool in_phys = math::IsWithinLimits(W, Wl);
    in_phys = in_phys && math::IsWithinLimits(Q2, Q2l);

    // extra check: there's a t minimum.
    // but only check if W, Q2 is reasonable
    // (otherwise get NaNs in tmin)
    if (in_phys)
    {
      double    t   = kine.t();
      Range1D_t tl  = this->TLim();
      LOG("KPhaseSpace", pDEBUG) << " t = " << t << ", limits = [" << tl.min << "," << tl.max << "];";
      in_phys = in_phys && math::IsWithinLimits(t, tl);
    }
    LOG("KPhaseSpace", pDEBUG) << " phase space point is " << ( in_phys ? "ALLOWED" : "NOT ALLOWED");


    bool allowed = in_phys;
    return allowed;
  }

  // was MECTensor
  if (pi.IsMEC()){
    Range1D_t Q2l = this->Q2Lim();
    double    Q2  = kine.Q2();
    bool in_phys = math::IsWithinLimits(Q2, Q2l);
    bool allowed = in_phys;
    return allowed;
  }

  return false;
}

Range1D_t KPhaseSpace::Limits

(

KineVar_t

kvar

)

const

Return the kinematical variable limits.

Definition at line 191 of file KPhaseSpace.cxx.

References ana::assert(), genie::KineVar::AsString(), fInteraction, genie::kKVQ2, genie::kKVq2, genie::kKVt, genie::kKVW, genie::kKVx, genie::kKVy, LOG, pERROR, Q2Lim(), q2Lim(), R, TLim(), WLim(), XLim(), and YLim().

Referenced by genie::ReinSehgalRESXSecWithCache::CacheResExcitationXSec(), genie::IBDKinematicsGenerator::ComputeMaxXSec(), genie::NuEKinematicsGenerator::ComputeMaxXSec(), genie::QELKinematicsGenerator::ComputeMaxXSec(), genie::DISKinematicsGenerator::ComputeMaxXSec(), genie::DMELKinematicsGenerator::ComputeMaxXSec(), genie::DMDISKinematicsGenerator::ComputeMaxXSec(), genie::IMDXSec::Integrate(), genie::COHXSecAR::Integrate(), genie::COHXSec::Integrate(), genie::RESXSec::Integrate(), genie::QELXSec::Integrate(), genie::DMELXSec::Integrate(), genie::DFRXSec::Integrate(), genie::NuElectronXSec::Integrate(), genie::SmithMonizQELCCXSec::Integrate(), genie::ReinSehgalRESXSec::Integrate(), Maximum(), Minimum(), genie::utils::kinematics::PhaseSpaceVolume(), PrintLimits(), genie::DFRKinematicsGenerator::ProcessEventRecord(), genie::NuEKinematicsGenerator::ProcessEventRecord(), genie::QELKinematicsGenerator::ProcessEventRecord(), genie::RESKinematicsGenerator::ProcessEventRecord(), genie::IBDKinematicsGenerator::ProcessEventRecord(), genie::DMELKinematicsGenerator::ProcessEventRecord(), genie::DISKinematicsGenerator::ProcessEventRecord(), genie::DMDISKinematicsGenerator::ProcessEventRecord(), genie::QELKinematicsGenerator::SpectralFuncExperimentalCode(), and genie::DMELKinematicsGenerator::SpectralFuncExperimentalCode().

{
  // Compute limits for the input kinematic variable irrespective of any other
  // relevant kinematical variable
  //
  assert(fInteraction);

  switch(kvar) {
  case(kKVW)  : return this->WLim();  break;
  case(kKVQ2) : return this->Q2Lim(); break;
  case(kKVq2) : return this->q2Lim(); break;
  case(kKVx)  : return this->XLim();  break;
  case(kKVy)  : return this->YLim();  break;
  case(kKVt)  : return this->TLim();  break;
  default:
    LOG("KPhaseSpace", pERROR)
      << "Couldn't compute limits for " << KineVar::AsString(kvar);
    Range1D_t R(-1.,-1);
    return R;
  }
}

double KPhaseSpace::Maximum

(

KineVar_t

kvar

)

const

Definition at line 219 of file KPhaseSpace.cxx.

References Limits(), and genie::Range1D_t::max.

{
  Range1D_t lim = this->Limits(kvar);
  return lim.max;
}

double KPhaseSpace::Minimum

(

KineVar_t

kvar

)

const

Definition at line 213 of file KPhaseSpace.cxx.

References Limits(), and genie::Range1D_t::min.

{
  Range1D_t lim = this->Limits(kvar);
  return lim.min;
}

Range1D_t KPhaseSpace::Q2Lim

(

void

)

const

Q2 limits.

Definition at line 500 of file KPhaseSpace.cxx.

References genie::XclsTag::CharmHadronPdg(), genie::utils::kinematics::CohQ2Lim(), genie::utils::kinematics::DarkQ2Lim(), genie::utils::kinematics::DarkQ2Lim_W(), genie::Interaction::ExclTag(), genie::PDGLibrary::Find(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4Ptr(), genie::utils::kinematics::InelQ2Lim(), genie::utils::kinematics::InelQ2Lim_W(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::XclsTag::IsCharmEvent(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDarkMatterElastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsEM(), genie::ProcessInfo::IsInverseBetaDecay(), genie::ProcessInfo::IsMEC(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::XclsTag::IsStrangeEvent(), genie::ProcessInfo::IsWeakCC(), genie::controls::kMinQ2Limit_VLE, kPi0Mass, genie::constants::kPionMass, genie::kRfHitNucRest, genie::Range1D_t::max, genie::Range1D_t::min, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::Interaction::RecoilNucleon(), genie::XclsTag::StrangeHadronPdg(), genie::InitialState::Tgt(), and W.

Referenced by genie::DISXSec::CacheFreeNucleonXSec(), genie::DMDISXSec::CacheFreeNucleonXSec(), genie::utils::ComputeFullQELPXSec(), genie::DISXSec::Integrate(), genie::DMDISXSec::Integrate(), IsAllowed(), Limits(), genie::COHKinematicsGenerator::MaxXSec_BergerSehgal(), genie::COHKinematicsGenerator::MaxXSec_BergerSehgalFM(), q2Lim(), and Q2Lim_W().

{
  // Computes momentum transfer (Q2>0) limits irrespective of the invariant mass
  // For QEL this is identical to Q2Lim_W (since W is fixed)
  // For RES & DIS, the calculation proceeds as in kinematics::InelQ2Lim().
  //
  Range1D_t Q2l;
  Q2l.min = -1;
  Q2l.max = -1;

  const ProcessInfo & pi = fInteraction->ProcInfo();

  bool is_em = pi.IsEM();
  bool is_qel   = pi.IsQuasiElastic()  || pi.IsInverseBetaDecay();
  bool is_inel  = pi.IsDeepInelastic() || pi.IsResonant();
  bool is_coh   = pi.IsCoherent();
  bool is_dme   = pi.IsDarkMatterElastic();
  bool is_dmdis = pi.IsDarkMatterDeepInelastic();

  if(!is_qel && !is_inel && !is_coh && !is_dme && !is_dmdis) return Q2l;

  const InitialState & init_state = fInteraction->InitState();
  double Ev  = init_state.ProbeE(kRfHitNucRest);
  double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
  double ml  = fInteraction->FSPrimLepton()->Mass();

  if(is_coh) {
    bool pionIsCharged = pi.IsWeakCC();
    double mpi = pionIsCharged ? kPionMass : kPi0Mass;
    Q2l = kinematics::CohQ2Lim(M, mpi, ml, Ev);
    return Q2l;
  }

  const XclsTag & xcls = fInteraction->ExclTag();

  // quasi-elastic
  if(is_qel) {
    double W = fInteraction->RecoilNucleon()->Mass();
    if(xcls.IsCharmEvent()) {
      int charm_pdgc = xcls.CharmHadronPdg();
      W = PDGLibrary::Instance()->Find(charm_pdgc)->Mass();
    }  else if(xcls.IsStrangeEvent()) {
      int strange_pdgc = xcls.StrangeHadronPdg();
      W = PDGLibrary::Instance()->Find(strange_pdgc)->Mass();
    }
    if (pi.IsInverseBetaDecay()) {
      Q2l = kinematics::InelQ2Lim_W(Ev,M,ml,W,controls::kMinQ2Limit_VLE);
    } else {
     Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W); 
    }

    return Q2l;
  }

    // dark mattter elastic
  if(is_dme) {
    double W = fInteraction->RecoilNucleon()->Mass();
    if(xcls.IsCharmEvent()) {
      int charm_pdgc = xcls.CharmHadronPdg();
      W = PDGLibrary::Instance()->Find(charm_pdgc)->Mass();
    }  else if(xcls.IsStrangeEvent()) {
      int strange_pdgc = xcls.StrangeHadronPdg();
      W = PDGLibrary::Instance()->Find(strange_pdgc)->Mass();
    }
    if (pi.IsInverseBetaDecay()) {
      Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W,controls::kMinQ2Limit_VLE);
    } else {
      Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W);
    }


    return Q2l;
  }

  // was MECTensor
  // TODO: Q2maxConfig
  if (pi.IsMEC()){
    double W = fInteraction->RecoilNucleon()->Mass();
    Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W); 
    double Q2maxConfig = 1.44; // need to pull from config file somehow?
    if (Q2l.max > Q2maxConfig) Q2l.max = Q2maxConfig;
    return Q2l;
  }

  if (is_dmdis) {
    Q2l = kinematics::DarkQ2Lim(Ev,M,ml);
    return Q2l;
  }

  // inelastic
  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim(Ev,ml,M) : kinematics::InelQ2Lim(Ev,M,ml); 
  return Q2l;
}

Range1D_t KPhaseSpace::q2Lim

(

void

)

const

q2 limits

Definition at line 594 of file KPhaseSpace.cxx.

References genie::Range1D_t::max, genie::Range1D_t::min, genie::utils::kinematics::Q2(), q2, and Q2Lim().

Referenced by Limits().

{
// As Q2Lim(void) but with reversed sign (Q2 -> q2)
//
  Range1D_t Q2 = this->Q2Lim();
  Range1D_t q2;
  q2.min = - Q2.max;
  q2.max = - Q2.min;
  return q2;
}

Range1D_t KPhaseSpace::Q2Lim_W

(

void

)

const

Q2 limits @ fixed W.

Definition at line 439 of file KPhaseSpace.cxx.

References genie::utils::kinematics::DarkQ2Lim_W(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4Ptr(), genie::utils::kinematics::InelQ2Lim_W(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDarkMatterElastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::ProcessInfo::IsEM(), genie::ProcessInfo::IsInverseBetaDecay(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::controls::kMinQ2Limit_VLE, genie::kRfHitNucRest, genie::Range1D_t::max, genie::Range1D_t::min, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), Q2Lim(), genie::Interaction::RecoilNucleon(), genie::InitialState::Tgt(), and W.

Referenced by genie::DFRKinematicsGenerator::ComputeMaxXSec(), genie::RESKinematicsGenerator::ComputeMaxXSec(), genie::utils::gsl::d2XSecRESFast_dWQ2_E::DoEval(), IsAllowed(), genie::utils::kinematics::PhaseSpaceVolume(), genie::RESKinematicsGenerator::ProcessEventRecord(), and q2Lim_W().

{
  // Computes momentum transfer (Q2>0) limits @ the input invariant mass
  // The calculation proceeds as in kinematics::InelQ2Lim_W().
  // For QEL, W is set to the recoil nucleon mass
  //
  // TODO: For now, choosing to handle Q2 at fixed W for coherent in the
  // same way as for the general Q2 limits... but shouldn't we just use
  // W = m_pi? - which we do in Q2Lim() anyway... seems like there are
  // cleanup opportunities here.

  Range1D_t Q2l;
  Q2l.min = -1;
  Q2l.max = -1;

  const ProcessInfo & pi = fInteraction->ProcInfo();

  bool is_em = pi.IsEM();
  bool is_qel   = pi.IsQuasiElastic()  || pi.IsInverseBetaDecay();
  bool is_inel  = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive();
  bool is_coh   = pi.IsCoherent();
  bool is_dme   = pi.IsDarkMatterElastic();
  bool is_dmdis = pi.IsDarkMatterDeepInelastic();

  if(!is_qel && !is_inel && !is_coh && !is_dme && !is_dmdis) return Q2l;

  if(is_coh) {
    return Q2Lim();
  }

  const InitialState & init_state = fInteraction->InitState();
  double Ev  = init_state.ProbeE(kRfHitNucRest);
  double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
  double ml  = fInteraction->FSPrimLepton()->Mass();

  double W = 0;
  if(is_qel || is_dme) W = fInteraction->RecoilNucleon()->Mass();
  else       W = kinematics::W(fInteraction);

  if (pi.IsInverseBetaDecay()) {
     Q2l = kinematics::InelQ2Lim_W(Ev,M,ml,W, controls::kMinQ2Limit_VLE);
  } else if (is_dme || is_dmdis) {
    Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W);
  } else {
     Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : Q2l = kinematics::InelQ2Lim_W(Ev,M,ml,W); 
  }

  return Q2l;
}

Range1D_t KPhaseSpace::q2Lim_W

(

void

)

const

q2 limits @ fixed W

Definition at line 489 of file KPhaseSpace.cxx.

References genie::Range1D_t::max, genie::Range1D_t::min, genie::utils::kinematics::Q2(), q2, and Q2Lim_W().

{
// As Q2Lim_W(void) but with reversed sign (Q2 -> q2)
//
  Range1D_t Q2 = this->Q2Lim_W();
  Range1D_t q2;
  q2.min = - Q2.max;
  q2.max = - Q2.min;
  return q2;
}

double KPhaseSpace::Threshold

(

void

)

const

Energy threshold.

Definition at line 81 of file KPhaseSpace.cxx.

References ana::assert(), genie::XclsTag::CharmHadronPdg(), genie::Interaction::ExclTag(), exit(), genie::PDGLibrary::Find(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Target::HitNucIsSet(), genie::Target::HitNucP4Ptr(), genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::ProcessInfo::IsAMNuGamma(), genie::XclsTag::IsCharmEvent(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDarkMatterElastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::ProcessInfo::IsGlashowResonance(), genie::ProcessInfo::IsIMDAnnihilation(), genie::XclsTag::IsInclusiveCharm(), genie::ProcessInfo::IsInverseBetaDecay(), genie::ProcessInfo::IsInverseMuDecay(), genie::ProcessInfo::IsMEC(), genie::ProcessInfo::IsNuElectronElastic(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::ProcessInfo::IsSingleKaon(), genie::ProcessInfo::IsWeakCC(), genie::controls::kASmallNum, genie::constants::kElectronMass, genie::constants::kElectronMass2, genie::constants::kLightestChmHad, genie::constants::kMuonMass2, genie::constants::kNeutronMass, genie::constants::kNucleonMass, genie::constants::kPhotontest, kPi0Mass, genie::constants::kPionMass, genie::constants::kProtonMass, m, Munits::m2, genie::XclsTag::NProtons(), genie::Target::Pdg(), pERROR, python.hepunit::pi, genie::Interaction::ProcInfo(), genie::Interaction::RecoilNucleon(), SLOG, genie::XclsTag::StrangeHadronPdg(), genie::pdg::SwitchProtonNeutron(), genie::InitialState::Tgt(), and Wmin.

Referenced by genie::DISXSec::CacheFreeNucleonXSec(), genie::DMDISXSec::CacheFreeNucleonXSec(), genie::ReinSehgalRESXSecWithCache::CacheResExcitationXSec(), genie::ReinSehgalRESXSecWithCacheFast::CacheResExcitationXSec(), genie::XSecSplineList::CreateSpline(), IsAboveThreshold(), and genie::FermiMover::KickHitNucleon().

{
  const ProcessInfo &  pi         = fInteraction->ProcInfo();
  const InitialState & init_state = fInteraction->InitState();
  const XclsTag &      xcls       = fInteraction->ExclTag();
  const Target &       tgt        = init_state.Tgt();

  double ml = fInteraction->FSPrimLepton()->Mass();

  if (pi.IsSingleKaon()) {
    int kaon_pdgc = xcls.StrangeHadronPdg();
    double Mi   = tgt.HitNucP4Ptr()->M(); // initial nucleon mass
    // Final nucleon can be different for K0 interaction
    double Mf = (xcls.NProtons()==1) ? kProtonMass : kNeutronMass;
    double mk   = PDGLibrary::Instance()->Find(kaon_pdgc)->Mass();
  //double ml   = PDGLibrary::Instance()->Find(fInteraction->FSPrimLeptonPdg())->Mass();
    double mtot = ml + mk + Mf; // total mass of FS particles
    double Ethresh = (mtot*mtot - Mi*Mi)/(2. * Mf);
    return Ethresh;
  }

  if (pi.IsCoherent()) {
    int tgtpdgc = tgt.Pdg(); // nuclear target PDG code (10LZZZAAAI)
    double mpi  = pi.IsWeakCC() ? kPionMass : kPi0Mass;
    double MA   = PDGLibrary::Instance()->Find(tgtpdgc)->Mass();
    double m    = ml + mpi;
    double m2   = TMath::Power(m,2);
    double Ethr = m + 0.5*m2/MA;
    return TMath::Max(0.,Ethr);
  }

  if(pi.IsQuasiElastic()            ||
     pi.IsDarkMatterElastic()       ||
     pi.IsInverseBetaDecay()        ||
     pi.IsResonant()                ||
     pi.IsDeepInelastic()           ||
     pi.IsDarkMatterDeepInelastic() ||
     pi.IsDiffractive())
  {
    assert(tgt.HitNucIsSet());
    double Mn   = tgt.HitNucP4Ptr()->M();
    double Mn2  = TMath::Power(Mn,2);
    double Wmin = kNucleonMass + kPionMass;
    if ( pi.IsQuasiElastic() || pi.IsDarkMatterElastic() || pi.IsInverseBetaDecay() ) {
      int finalNucPDG = tgt.HitNucPdg();
      if ( pi.IsWeakCC() ) finalNucPDG = pdg::SwitchProtonNeutron( finalNucPDG );
      Wmin = PDGLibrary::Instance()->Find( finalNucPDG )->Mass();
    }
    if (pi.IsResonant()) {
        Wmin = kNucleonMass + kPhotontest;
    }

    if(xcls.IsCharmEvent()) {
       if(xcls.IsInclusiveCharm()) {
          Wmin = kNucleonMass+kLightestChmHad;
       } else {
          int cpdg = xcls.CharmHadronPdg();
          double mchm = PDGLibrary::Instance()->Find(cpdg)->Mass();
          if(pi.IsQuasiElastic() || pi.IsInverseBetaDecay()) {
            Wmin = mchm + controls::kASmallNum;
          }
          else {
            Wmin = kNeutronMass + mchm + controls::kASmallNum;
          }
       }//incl.?
    }//charm?

    double smin = TMath::Power(Wmin+ml,2.);
    double Ethr = 0.5*(smin-Mn2)/Mn;
    // threshold is different for dark matter case
    if(pi.IsDarkMatterElastic() || pi.IsDarkMatterDeepInelastic()) {
      Ethr = TMath::Max(0.5*(smin-Mn2-ml*ml)/Mn,ml);
    }

    return TMath::Max(0.,Ethr);
  }

  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation()) {
    double Ethr = 0.5 * (kMuonMass2-kElectronMass2)/kElectronMass;
    return TMath::Max(0.,Ethr);
  }

  if(pi.IsNuElectronElastic() || pi.IsGlashowResonance() ) {
    return 0;
  }
  if(pi.IsAMNuGamma()) {
    return 0;
  }
  if (pi.IsMEC()) {
    if (tgt.HitNucIsSet()) {
        double Mn   = tgt.HitNucP4Ptr()->M();
        double Mn2  = TMath::Power(Mn,2);
        double Wmin = fInteraction->RecoilNucleon()->Mass(); // mass of the recoil nucleon cluster
        double smin = TMath::Power(Wmin+ml,2.);
        double Ethr = 0.5*(smin-Mn2)/Mn;
        return TMath::Max(0.,Ethr);
    }
    else {
        // this was ... if (pi.IsMECTensor())
        return ml;
    }
  }

  SLOG("KPhaseSpace", pERROR)
         << "Can't compute threshold for \n" << *fInteraction;
  exit(1);

  return 99999999;
}

Range1D_t KPhaseSpace::TLim

(

void

)

const

t limits

Definition at line 809 of file KPhaseSpace.cxx.

References std::abs(), fInteraction, GetTMaxDFR(), genie::Target::HitNucMass(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDiffractive(), std::isnan(), genie::ProcessInfo::IsWeakCC(), genie::Interaction::Kine(), kPi0Mass, genie::constants::kPionMass, genie::kRfHitNucRest, LOG, genie::Range1D_t::max, genie::Range1D_t::min, pERROR, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), pWARN, genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), std::sqrt(), genie::InitialState::Tgt(), moon_position_table_new3::tl(), genie::utils::kinematics::UpdateWQ2FromXY(), and genie::Kinematics::y().

Referenced by genie::DFRKinematicsGenerator::ComputeMaxXSec(), IsAllowed(), and Limits().

{
  // t limits for Coherent pion production from
  //   Kartavtsev, Paschos, and Gounaris, PRD 74 054007, and
  //   Paschos and Schalla, PRD 80, 03305
  // TODO: Attempt to assign t bounds for other reactions?
  Range1D_t tl;
  tl.min = -1;
  tl.max = -1;

  const InitialState & init_state = fInteraction->InitState();
  const ProcessInfo & pi = fInteraction->ProcInfo();
  const Kinematics & kine = fInteraction->Kine();
  kinematics::UpdateWQ2FromXY(fInteraction);
  double Ev = init_state.ProbeE(kRfHitNucRest);
  double Q2 = kine.Q2();
  double nu = Ev * kine.y();
  bool pionIsCharged = pi.IsWeakCC();
  double mpi = pionIsCharged ? kPionMass : kPi0Mass;
  double mpi2 = mpi*mpi;

  //COH
  if(pi.IsCoherent()) {
    tl.min = 1.0 * (Q2 + mpi2)/(2.0 * nu) * (Q2 + mpi2)/(2.0 * nu);
    tl.max = 0.05;
    return tl;
  }
  // DFR
  else if (pi.IsDiffractive()) {
    // diffractive tmin from Nucl.Phys.B278,61 (1986), eq. 12
    double M = init_state.Tgt().HitNucMass();
    double M2 = M*M;
    double nuSqPlusQ2 = nu*nu + Q2;
    double nuOverM = nu / M;
    double mpiQ2term = mpi2 - Q2 - 2*nu*nu;
    double A1 = 1 + 2*nuOverM + nuOverM*nuOverM - nuSqPlusQ2/M2;
    double A2 = (1+nuOverM) * mpiQ2term + 2*nuOverM*nuSqPlusQ2;
    double A3 = mpiQ2term*mpiQ2term - 4*nuSqPlusQ2*(nu*nu - mpi2);

    tl.min = std::abs( (A2 + sqrt(A2*A2 - A1*A3)) / A1 );  // GENIE's convention is that t is positive
    bool tminIsNaN;
    // use std::isnan when C++11 is around
#if __cplusplus >= 201103L
      tminIsNaN = std::isnan(tl.min);
#else
      // this the old-fashioned way to check for NaN:
      // NaN's aren't equal to anything, including themselves
      tminIsNaN = tl.min != tl.min;
#endif
    if (tminIsNaN)
    {
      LOG("KPhaseSpace", pERROR)
        << "tmin for diffractive scattering is NaN "
        << "( Enu = " << Ev << ", Q2 = " << Q2 << ", nu = " << nu << ")";
      throw genie::exceptions::InteractionException("NaN tmin for diffractive scattering");
    }
    tl.max = this->GetTMaxDFR();

    return tl;
  }

  // RES+DIS
  // IMD
  LOG("KPhaseSpace", pWARN) << "It is not sensible to ask for t limits for events that are not coherent or diffractive.";
  return tl;
}

void KPhaseSpace::UseInteraction

(

const Interaction *

in

)

Definition at line 76 of file KPhaseSpace.cxx.

References fInteraction, and in.

Referenced by KPhaseSpace().

{
  fInteraction = in;
}

Range1D_t KPhaseSpace::WLim

(

void

)

const

W limits.

Definition at line 369 of file KPhaseSpace.cxx.

References genie::utils::kinematics::DarkWLim(), genie::Interaction::ExclTag(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4Ptr(), genie::utils::kinematics::InelWLim(), genie::Interaction::InitState(), genie::XclsTag::IsCharmEvent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDarkMatterElastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::ProcessInfo::IsEM(), genie::ProcessInfo::IsInverseBetaDecay(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::constants::kLightestChmHad, genie::constants::kNeutronMass, genie::constants::kPionMass, genie::kRfHitNucRest, LOG, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::Interaction::RecoilNucleon(), and genie::InitialState::Tgt().

Referenced by genie::DISXSec::CacheFreeNucleonXSec(), genie::DMDISXSec::CacheFreeNucleonXSec(), genie::DFRKinematicsGenerator::ComputeMaxXSec(), genie::RESKinematicsGenerator::ComputeMaxXSec(), genie::utils::gsl::d2XSecRESFast_dWQ2_E::d2XSecRESFast_dWQ2_E(), genie::DISXSec::Integrate(), genie::DMDISXSec::Integrate(), IsAllowed(), and Limits().

{
// Computes hadronic invariant mass limits.
// For QEL the range reduces to the recoil nucleon mass.
// For DIS & RES the calculation proceeds as in kinematics::InelWLim().
// It is not computed for other interactions
//
  Range1D_t Wl;
  Wl.min = -1;
  Wl.max = -1;

  const ProcessInfo & pi = fInteraction->ProcInfo();

  bool is_em = pi.IsEM();
  bool is_qel  = pi.IsQuasiElastic()  || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic();
  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive();
  bool is_dmdis = pi.IsDarkMatterDeepInelastic();

  if(is_qel) {
    double MR = fInteraction->RecoilNucleon()->Mass();
    Wl.min = MR;
    Wl.max = MR;
    return Wl;
  }
  if(is_inel) {
    const InitialState & init_state = fInteraction->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M(); //can be off m/shell
    double ml = fInteraction->FSPrimLepton()->Mass();

    Wl = is_em ? kinematics::electromagnetic::InelWLim(Ev,ml,M) : kinematics::InelWLim(Ev,M,ml); 

    if(fInteraction->ExclTag().IsCharmEvent()) {
      //Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass+kLightestChmHad);
      Wl.min = TMath::Max(Wl.min, kNeutronMass+kLightestChmHad);
    }
    else if (fInteraction->ProcInfo().IsDiffractive())
      Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass);
    else if ( fInteraction->ProcInfo().IsDeepInelastic() ) {
      Wl.min = TMath::Max(Wl.min, kNeutronMass + kPionMass );
    }

    // sanity check
    if(Wl.min>Wl.max) {Wl.min=-1; Wl.max=-1;}

    return Wl;
  }
  if(is_dmdis) {
    const InitialState & init_state = fInteraction->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M(); //can be off m/shell
    double ml = fInteraction->FSPrimLepton()->Mass();
    Wl = kinematics::DarkWLim(Ev,M,ml);
    if(fInteraction->ExclTag().IsCharmEvent()) {
      //Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass+kLightestChmHad);
      Wl.min = TMath::Max(Wl.min, kNeutronMass+kLightestChmHad);
    }
    else if (fInteraction->ProcInfo().IsDiffractive())
      Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass);

    LOG("KPhaseSpace", pDEBUG) << "Found nominal limits: " << Wl.min << ", " << Wl.max;

    // sanity check
    if(Wl.min>Wl.max) {Wl.min=-1; Wl.max=-1;}

    return Wl;
  }
  return Wl;
}

Range1D_t KPhaseSpace::XLim

(

void

)

const

x limits

Definition at line 605 of file KPhaseSpace.cxx.

References genie::utils::kinematics::CohXLim(), genie::utils::kinematics::DarkXLim(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4Ptr(), genie::utils::kinematics::InelXLim(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDarkMatterElastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::ProcessInfo::IsEM(), genie::ProcessInfo::IsInverseBetaDecay(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsResonant(), genie::controls::kASmallNum, genie::kRfHitNucRest, genie::Range1D_t::max, genie::Range1D_t::min, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by genie::DFRKinematicsGenerator::ComputeMaxXSec(), genie::utils::gsl::d5XSecAR::DoEval(), genie::utils::gsl::d5Xsec_dEldOmegaldOmegapi::DoEval(), genie::utils::gsl::d4Xsec_dEldThetaldOmegapi::DoEval(), genie::utils::gsl::d3Xsec_dOmegaldThetapi::DoEval(), IsAllowed(), and Limits().

{
  // Computes x-limits;

  Range1D_t xl;
  xl.min = -1;
  xl.max = -1;

  const ProcessInfo & pi = fInteraction->ProcInfo();
  bool is_em = pi.IsEM();

  //RES+DIS
  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
  if(is_inel) {
    const InitialState & init_state  = fInteraction->InitState();
    double Ev  = init_state.ProbeE(kRfHitNucRest);
    double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
    double ml  = fInteraction->FSPrimLepton()->Mass();
    xl = is_em ? kinematics::electromagnetic::InelXLim(Ev,ml,M) : kinematics::InelXLim(Ev,M,ml); 
    return xl;
  }
  //DMDIS
  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
  if(is_dmdis) {
    const InitialState & init_state  = fInteraction->InitState();
    double Ev  = init_state.ProbeE(kRfHitNucRest);
    double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
    double ml  = fInteraction->FSPrimLepton()->Mass();
    xl = kinematics::DarkXLim(Ev,M,ml);
    return xl;
  }
  //COH
  bool is_coh = pi.IsCoherent();
  if(is_coh) {
    xl = kinematics::CohXLim();
    return xl;
  }
  //QEL
  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic();
  if(is_qel) {
    xl.min = 1;
    xl.max = 1;
    return xl;
  }
  bool is_dfr = pi.IsDiffractive();
  if(is_dfr) {
    xl.min =      controls::kASmallNum;
    xl.max = 1. - controls::kASmallNum;
    return xl;
  }

  return xl;
}

Range1D_t KPhaseSpace::YLim

(

void

)

const

y limits

Definition at line 659 of file KPhaseSpace.cxx.

References genie::utils::kinematics::CohYLim(), genie::utils::kinematics::DarkYLim(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4Ptr(), genie::utils::kinematics::InelYLim(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsDiffractive(), genie::ProcessInfo::IsEM(), genie::ProcessInfo::IsIMDAnnihilation(), genie::ProcessInfo::IsInverseMuDecay(), genie::ProcessInfo::IsNuElectronElastic(), genie::ProcessInfo::IsResonant(), genie::controls::kASmallNum, genie::constants::kElectronMass, genie::constants::kPionMass, genie::kRfHitNucRest, genie::kRfLab, genie::Range1D_t::max, genie::Range1D_t::min, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by genie::COHKinematicsGenerator::CalculateKin_BergerSehgal(), genie::COHKinematicsGenerator::CalculateKin_BergerSehgalFM(), genie::COHKinematicsGenerator::CalculateKin_ReinSehgal(), genie::DFRKinematicsGenerator::ComputeMaxXSec(), genie::COHElKinematicsGenerator::ComputeMaxXSec(), IsAllowed(), Limits(), genie::COHKinematicsGenerator::MaxXSec_AlvarezRuso(), genie::COHKinematicsGenerator::MaxXSec_BergerSehgal(), genie::COHKinematicsGenerator::MaxXSec_BergerSehgalFM(), genie::COHKinematicsGenerator::MaxXSec_ReinSehgal(), genie::COHElKinematicsGenerator::ProcessEventRecord(), YLim(), and YLim_X().

{
  Range1D_t yl;
  yl.min = -1;
  yl.max = -1;

  const ProcessInfo & pi = fInteraction->ProcInfo();
  bool is_em = pi.IsEM();

  //RES+DIS
  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
  if(is_inel) {
    const InitialState & init_state = fInteraction->InitState();
    double Ev  = init_state.ProbeE(kRfHitNucRest);
    double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
    double ml  = fInteraction->FSPrimLepton()->Mass();
    yl = is_em ? kinematics::electromagnetic::InelYLim(Ev,ml,M) : kinematics::InelYLim(Ev,M,ml); 
    return yl;
  }
  //DMDIS
  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
  if(is_dmdis) {
    const InitialState & init_state = fInteraction->InitState();
    double Ev  = init_state.ProbeE(kRfHitNucRest);
    double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
    double ml  = fInteraction->FSPrimLepton()->Mass();
    yl = kinematics::DarkYLim(Ev,M,ml);
    return yl;
  }
  //COH
  bool is_coh = pi.IsCoherent();
  if(is_coh) {
    const InitialState & init_state = fInteraction->InitState();
    double EvL = init_state.ProbeE(kRfLab);
    double ml  = fInteraction->FSPrimLepton()->Mass();
    yl = kinematics::CohYLim(EvL,ml);
    return yl;
  }
  // IMD
  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation() || pi.IsNuElectronElastic()) {
    const InitialState & init_state = fInteraction->InitState();
    double Ev = init_state.ProbeE(kRfLab);
    double ml = fInteraction->FSPrimLepton()->Mass();
    double me = kElectronMass;
    yl.min = controls::kASmallNum;
    yl.max = 1 - (ml*ml + me*me)/(2*me*Ev) - controls::kASmallNum;
    return yl;
  }
  bool is_dfr = pi.IsDiffractive();
  if(is_dfr) {
    const InitialState & init_state = fInteraction -> InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double ml = fInteraction->FSPrimLepton()->Mass();
    yl.min = kPionMass/Ev + controls::kASmallNum;
    yl.max = 1. -ml/Ev - controls::kASmallNum;
    return yl;
  }
  return yl;
}

Range1D_t KPhaseSpace::YLim

(

double

xsi

)

const

y limits (COH)

Definition at line 764 of file KPhaseSpace.cxx.

References genie::utils::kinematics::CohYLim(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsWeakCC(), genie::Interaction::Kine(), kPi0Mass, genie::constants::kPionMass, genie::kRfHitNucRest, genie::Target::Mass(), genie::Range1D_t::max, genie::Range1D_t::min, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), genie::InitialState::Tgt(), and YLim().

{
  // Paschos-Schalla xsi parameter for y-limits in COH
  // From PRD 80, 033005 (2009)

  Range1D_t yl;
  yl.min = -1;
  yl.max = -1;

  const ProcessInfo & pi = fInteraction->ProcInfo();

  //COH
  bool is_coh = pi.IsCoherent();
  if(is_coh) {
    const InitialState & init_state = fInteraction->InitState();
    const Kinematics & kine = fInteraction->Kine();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double Q2 = kine.Q2();
    bool pionIsCharged = pi.IsWeakCC();
    double Mn = init_state.Tgt().Mass();
    double mpi = pionIsCharged ? kPionMass : kPi0Mass;
    double mlep = fInteraction->FSPrimLepton()->Mass();
    yl = kinematics::CohYLim(Mn, mpi, mlep, Ev, Q2, xsi);
    return yl;
  } else {
    return this->YLim();
  }
}

Range1D_t KPhaseSpace::YLim_X

(

void

)

const

y limits @ fixed x

Definition at line 719 of file KPhaseSpace.cxx.

References genie::utils::kinematics::CohYLim(), genie::utils::kinematics::DarkYLim_X(), fInteraction, genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4Ptr(), genie::utils::kinematics::InelYLim_X(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherent(), genie::ProcessInfo::IsDarkMatterDeepInelastic(), genie::ProcessInfo::IsDeepInelastic(), genie::ProcessInfo::IsEM(), genie::ProcessInfo::IsResonant(), genie::Interaction::Kine(), genie::kRfHitNucRest, genie::kRfLab, genie::Range1D_t::max, genie::Range1D_t::min, python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::InitialState::Tgt(), genie::Kinematics::x(), and submit_syst::x.

Referenced by YLim_X().

{
// Computes kinematical limits for y @ the input x

  Range1D_t yl;
  yl.min = -1;
  yl.max = -1;

  const ProcessInfo & pi = fInteraction->ProcInfo();
  bool is_em = pi.IsEM();

  //RES+DIS
  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
  if(is_inel) {
    const InitialState & init_state = fInteraction->InitState();
    double Ev  = init_state.ProbeE(kRfHitNucRest);
    double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
    double ml  = fInteraction->FSPrimLepton()->Mass();
    double x   = fInteraction->Kine().x();
    yl = is_em ? kinematics::electromagnetic::InelYLim_X(Ev,ml,M,x) : kinematics::InelYLim_X(Ev,M,ml,x); 
    return yl;
  }
  //DMDIS
  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
  if(is_dmdis) {
    const InitialState & init_state = fInteraction->InitState();
    double Ev  = init_state.ProbeE(kRfHitNucRest);
    double M   = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
    double ml  = fInteraction->FSPrimLepton()->Mass();
    double x   = fInteraction->Kine().x();
    yl = kinematics::DarkYLim_X(Ev,M,ml,x);
    return yl;
  }
  //COH
  bool is_coh = pi.IsCoherent();
  if(is_coh) {
    const InitialState & init_state = fInteraction->InitState();
    double EvL = init_state.ProbeE(kRfLab);
    double ml  = fInteraction->FSPrimLepton()->Mass();
    yl = kinematics::CohYLim(EvL,ml);
    return yl;
  }
  return yl;
}

Range1D_t KPhaseSpace::YLim_X

(

double

xsi

)

const

y limits @ fixed x (COH)

Definition at line 793 of file KPhaseSpace.cxx.

References fInteraction, genie::ProcessInfo::IsCoherent(), python.hepunit::pi, genie::Interaction::ProcInfo(), YLim(), and YLim_X().

{
  // Paschos-Schalla xsi parameter for y-limits in COH
  // From PRD 80, 033005 (2009)

  const ProcessInfo & pi = fInteraction->ProcInfo();

  //COH
  bool is_coh = pi.IsCoherent();
  if(is_coh) {
    return this->YLim(xsi);
  } else {
    return this->YLim_X();
  }
}

Member Data Documentation

const Interaction* genie::KPhaseSpace::fInteraction

private

Definition at line 74 of file KPhaseSpace.h.

Referenced by IsAboveThreshold(), IsAllowed(), Limits(), Q2Lim(), Q2Lim_W(), Threshold(), TLim(), UseInteraction(), WLim(), XLim(), YLim(), and YLim_X().

---

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

• Framework/Interaction/KPhaseSpace.h
• Framework/Interaction/KPhaseSpace.cxx