Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::KPhaseSpace Class Reference

Kinematical phase space. More...

#include <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 InteractionfInteraction
 

Detailed Description

Kinematical phase space.

Author
Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

May 06, 2004

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

Definition at line 33 of file KPhaseSpace.h.

Constructor & Destructor Documentation

genie::KPhaseSpace::KPhaseSpace ( void  )
KPhaseSpace::KPhaseSpace ( const Interaction in)

Definition at line 46 of file KPhaseSpace.cxx.

46  :
47 TObject(), fInteraction(NULL)
48 {
49  this->UseInteraction(in);
50 }
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
void UseInteraction(const Interaction *in)
Definition: KPhaseSpace.cxx:75
KPhaseSpace::~KPhaseSpace ( void  )

Definition at line 52 of file KPhaseSpace.cxx.

53 {
54 
55 }

Member Function Documentation

double KPhaseSpace::GetTMaxDFR ( )
static

Definition at line 57 of file KPhaseSpace.cxx.

58 {
59  static bool tMaxLoaded = false;
60  static double DFR_tMax = -1;
61 
62  if (!tMaxLoaded)
63  {
65  const Registry * r = confp->CommonList( "Param", "Diffractive" ) ;
66  double tmax = r->GetDouble("DFR-t-max");
67  DFR_tMax = tmax;
68  tMaxLoaded = true;
69  }
70 
71  return DFR_tMax;
72 
73 }
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:40
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:474
Registry * CommonList(const string &file_id, const string &set_name) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
static AlgConfigPool * Instance()
void genie::KPhaseSpace::Init ( void  )
private
bool KPhaseSpace::IsAboveThreshold ( void  ) const

Checks whether the interaction is above the energy threshold.

Definition at line 249 of file KPhaseSpace.cxx.

250 {
251  double E = 0.;
252  double Ethr = this->Threshold();
253 
254  const ProcessInfo & pi = fInteraction->ProcInfo();
255  const InitialState & init_state = fInteraction->InitState();
256 
257  if (pi.IsCoherentElastic() ||
258  pi.IsCoherentProduction() ||
259  pi.IsInverseMuDecay() ||
260  pi.IsIMDAnnihilation() ||
261  pi.IsNuElectronElastic() ||
263  pi.IsMEC() ||
264  pi.IsGlashowResonance())
265  {
266  E = init_state.ProbeE(kRfLab);
267  }
268 
269  if(pi.IsQuasiElastic() ||
270  pi.IsDarkMatterElastic() ||
271  pi.IsInverseBetaDecay() ||
272  pi.IsResonant() ||
273  pi.IsDeepInelastic() ||
275  pi.IsDiffractive() ||
276  pi.IsSingleKaon() ||
277  pi.IsAMNuGamma())
278  {
279  E = init_state.ProbeE(kRfHitNucRest);
280  }
281 
282  LOG("KPhaseSpace", pDEBUG) << "E = " << E << ", Ethr = " << Ethr;
283  return (E>Ethr);
284 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
bool IsDarkMatterElectronElastic(void) const
double Threshold(void) const
Energy threshold.
Definition: KPhaseSpace.cxx:80
bool IsInverseMuDecay(void) const
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
bool IsInverseBetaDecay(void) const
bool IsDiffractive(void) const
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
bool IsIMDAnnihilation(void) const
bool IsSingleKaon(void) const
Definition: ProcessInfo.cxx:79
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsCoherentElastic(void) const
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
bool IsNuElectronElastic(void) const
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsAMNuGamma(void) const
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:74
bool IsMEC(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
float pi
Definition: units.py:11
bool IsGlashowResonance(void) const
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
bool KPhaseSpace::IsAllowed ( void  ) const

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

Definition at line 286 of file KPhaseSpace.cxx.

287 {
288  const ProcessInfo & pi = fInteraction->ProcInfo();
289  const Kinematics & kine = fInteraction->Kine();
290 
291  // ASK single kaon:
292  // XSec code returns zero when kinematics are not allowed
293  // Here just let kinematics always be allowed
294  if(pi.IsSingleKaon()) {
295  return true;
296  }
297 
298  // QEL:
299  // Check the running Q2 vs the Q2 limits
300  if(pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic()) {
301  Range1D_t Q2l = this->Q2Lim();
302  double Q2 = kine.Q2();
303  bool in_phys = math::IsWithinLimits(Q2, Q2l);
304  bool allowed = in_phys;
305  return allowed;
306  }
307 
308  // RES
309  // Check the running W vs the W limits
310  // & the running Q2 vs Q2 limits for the given W
311  if(pi.IsResonant()) {
312  Range1D_t Wl = this->WLim();
313  Range1D_t Q2l = this->Q2Lim_W();
314  double W = kine.W();
315  double Q2 = kine.Q2();
316  bool in_phys = (math::IsWithinLimits(Q2, Q2l) && math::IsWithinLimits(W, Wl));
317  bool allowed = in_phys;
318  return allowed;
319  }
320 
321  // DIS
322  if(pi.IsDeepInelastic() || pi.IsDarkMatterDeepInelastic()) {
323  Range1D_t Wl = this->WLim();
324  Range1D_t Q2l = this->Q2Lim_W();
325  double W = kine.W();
326  double Q2 = kine.Q2();
327  bool in_phys = (math::IsWithinLimits(Q2, Q2l) && math::IsWithinLimits(W, Wl));
328  bool allowed = in_phys;
329  return allowed;
330  }
331 
332  //IMD
334  Range1D_t yl = this->YLim();
335  double y = kine.y();
336  bool in_phys = math::IsWithinLimits(y, yl);
337  bool allowed = in_phys;
338  return allowed;
339  }
340 
341  //COH
342  if (pi.IsCoherentProduction()) {
343  Range1D_t xl = this->XLim();
344  Range1D_t yl = this->YLim();
345  double x = kine.x();
346  double y = kine.y();
347  bool in_phys = (math::IsWithinLimits(x, xl) && math::IsWithinLimits(y, yl));
348  bool allowed = in_phys;
349  return allowed;
350  }
351 
352  // CEvNS
353  if (pi.IsCoherentElastic()) {
354  double Q2 = kine.Q2();
355  bool allowed (Q2 > 0);
356  return allowed;
357  }
358 
359  // DFR
360  if (pi.IsDiffractive()) {
361  // first two checks are the same as RES & DIS
362  Range1D_t Wl = this->WLim();
363  Range1D_t Q2l = this->Q2Lim_W();
364 
366  double W = kine.W();
367  double Q2 = kine.Q2();
368 
369  LOG("KPhaseSpace", pDEBUG) << " W = " << W << ", limits = [" << Wl.min << "," << Wl.max << "];";
370  LOG("KPhaseSpace", pDEBUG) << " Q2 = " << Q2 << ", limits = [" << Q2l.min << "," << Q2l.max << "];";
371  bool in_phys = math::IsWithinLimits(W, Wl);
372  in_phys = in_phys && math::IsWithinLimits(Q2, Q2l);
373 
374  // extra check: there's a t minimum.
375  // but only check if W, Q2 is reasonable
376  // (otherwise get NaNs in tmin)
377  if (in_phys)
378  {
379  double t = kine.t();
380  Range1D_t tl = this->TLim();
381  LOG("KPhaseSpace", pDEBUG) << " t = " << t << ", limits = [" << tl.min << "," << tl.max << "];";
382  in_phys = in_phys && math::IsWithinLimits(t, tl);
383  }
384  LOG("KPhaseSpace", pDEBUG) << " phase space point is " << ( in_phys ? "ALLOWED" : "NOT ALLOWED");
385 
386 
387  bool allowed = in_phys;
388  return allowed;
389  }
390 
391  // was MECTensor
392  if (pi.IsMEC()){
393  Range1D_t Q2l = this->Q2Lim();
394  double Q2 = kine.Q2();
395  bool in_phys = math::IsWithinLimits(Q2, Q2l);
396  bool allowed = in_phys;
397  return allowed;
398  }
399 
400  return false;
401 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
double W(bool selected=false) const
Definition: Kinematics.cxx:157
bool IsDarkMatterElectronElastic(void) const
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
A simple [min,max] interval for doubles.
Definition: Range1.h:42
bool IsInverseMuDecay(void) const
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:258
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
bool IsInverseBetaDecay(void) const
double x(bool selected=false) const
Definition: Kinematics.cxx:99
Range1D_t YLim(void) const
y limits
bool IsDiffractive(void) const
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
bool IsIMDAnnihilation(void) const
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
Range1D_t Q2Lim(void) const
Q2 limits.
double y(bool selected=false) const
Definition: Kinematics.cxx:112
bool IsSingleKaon(void) const
Definition: ProcessInfo.cxx:79
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsCoherentElastic(void) const
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
bool IsNuElectronElastic(void) const
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const Kinematics & Kine(void) const
Definition: Interaction.h:71
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:74
bool IsMEC(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1277
double max
Definition: Range1.h:53
Range1D_t XLim(void) const
x limits
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
list x
Definition: train.py:276
double min
Definition: Range1.h:52
double t(bool selected=false) const
Definition: Kinematics.cxx:170
float pi
Definition: units.py:11
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
Range1D_t TLim(void) const
t limits
bool IsGlashowResonance(void) const
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
Range1D_t WLim(void) const
W limits.
#define pDEBUG
Definition: Messenger.h:63
Range1D_t KPhaseSpace::Limits ( KineVar_t  kvar) const

Return the kinematical variable limits.

Definition at line 215 of file KPhaseSpace.cxx.

216 {
217  // Compute limits for the input kinematic variable irrespective of any other
218  // relevant kinematical variable
219  //
220  assert(fInteraction);
221 
222  switch(kvar) {
223  case(kKVW) : return this->WLim(); break;
224  case(kKVQ2) : return this->Q2Lim(); break;
225  case(kKVq2) : return this->q2Lim(); break;
226  case(kKVx) : return this->XLim(); break;
227  case(kKVy) : return this->YLim(); break;
228  case(kKVt) : return this->TLim(); break;
229  default:
230  LOG("KPhaseSpace", pERROR)
231  << "Couldn't compute limits for " << KineVar::AsString(kvar);
232  Range1D_t R(-1.,-1);
233  return R;
234  }
235 }
#define pERROR
Definition: Messenger.h:59
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t q2Lim(void) const
q2 limits
Range1D_t YLim(void) const
y limits
Range1D_t Q2Lim(void) const
Q2 limits.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
static string AsString(KineVar_t kv)
Definition: KineVar.h:71
Range1D_t XLim(void) const
x limits
Range1D_t TLim(void) const
t limits
Range1D_t WLim(void) const
W limits.
double KPhaseSpace::Maximum ( KineVar_t  kvar) const

Definition at line 243 of file KPhaseSpace.cxx.

244 {
245  Range1D_t lim = this->Limits(kvar);
246  return lim.max;
247 }
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
double max
Definition: Range1.h:53
double KPhaseSpace::Minimum ( KineVar_t  kvar) const

Definition at line 237 of file KPhaseSpace.cxx.

238 {
239  Range1D_t lim = this->Limits(kvar);
240  return lim.min;
241 }
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
double min
Definition: Range1.h:52
Range1D_t KPhaseSpace::Q2Lim ( void  ) const

Q2 limits.

Definition at line 534 of file KPhaseSpace.cxx.

535 {
536  // Computes momentum transfer (Q2>0) limits irrespective of the invariant mass
537  // For QEL this is identical to Q2Lim_W (since W is fixed)
538  // For RES & DIS, the calculation proceeds as in kinematics::InelQ2Lim().
539  //
540  Range1D_t Q2l;
541  Q2l.min = -1;
542  Q2l.max = -1;
543 
544  const ProcessInfo & pi = fInteraction->ProcInfo();
545 
546  bool is_em = pi.IsEM();
547  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay();
548  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
549  bool is_coh = pi.IsCoherentProduction();
550  bool is_cevns = pi.IsCoherentElastic();
551  bool is_dme = pi.IsDarkMatterElastic();
552  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
553 
554  if(!is_qel && !is_inel && !is_coh && !is_cevns && !is_dme && !is_dmdis) return Q2l;
555 
556  const InitialState & init_state = fInteraction->InitState();
557  double Ev = init_state.ProbeE(kRfHitNucRest);
558  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
559  double ml = fInteraction->FSPrimLepton()->Mass();
560 
561  if(is_cevns) {
562  double Ev_lab = init_state.ProbeE(kRfLab);
563  Q2l = kinematics::CEvNSQ2Lim(Ev_lab);
564  return Q2l;
565  }
566 
567  const XclsTag & xcls = fInteraction->ExclTag();
568 
569  if(is_coh) {
570 
571  double m_other = controls::kASmallNum ;
572  // as a default the mass of hadronic system is the mass of the photon.
573  // which is assumed to be a small number to avoid divergences
574 
575  if ( xcls.NPions() > 0 ) {
576  bool pionIsCharged = pi.IsWeakCC();
577  m_other = pionIsCharged ? kPionMass : kPi0Mass;
578  }
579 
580  Q2l = kinematics::CohQ2Lim(M, m_other, ml, Ev);
581  return Q2l;
582  }
583 
584  // quasi-elastic
585  if(is_qel) {
586  double W = fInteraction->RecoilNucleon()->Mass();
587  if(xcls.IsCharmEvent()) {
588  int charm_pdgc = xcls.CharmHadronPdg();
589  W = PDGLibrary::Instance()->Find(charm_pdgc)->Mass();
590  } else if(xcls.IsStrangeEvent()) {
591  int strange_pdgc = xcls.StrangeHadronPdg();
592  W = PDGLibrary::Instance()->Find(strange_pdgc)->Mass();
593  }
594  if (pi.IsInverseBetaDecay()) {
596  } else {
597  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W);
598  }
599 
600  return Q2l;
601  }
602 
603  // dark mattter elastic
604  if(is_dme) {
605  double W = fInteraction->RecoilNucleon()->Mass();
606  if(xcls.IsCharmEvent()) {
607  int charm_pdgc = xcls.CharmHadronPdg();
608  W = PDGLibrary::Instance()->Find(charm_pdgc)->Mass();
609  } else if(xcls.IsStrangeEvent()) {
610  int strange_pdgc = xcls.StrangeHadronPdg();
611  W = PDGLibrary::Instance()->Find(strange_pdgc)->Mass();
612  }
613  if (pi.IsInverseBetaDecay()) {
615  } else {
616  Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W);
617  }
618 
619 
620  return Q2l;
621  }
622 
623  // was MECTensor
624  // TODO: Q2maxConfig
625  if (pi.IsMEC()){
626  double W = fInteraction->RecoilNucleon()->Mass();
627  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W);
628  double Q2maxConfig = 1.44; // need to pull from config file somehow?
629  if (Q2l.max > Q2maxConfig) Q2l.max = Q2maxConfig;
630  return Q2l;
631  }
632 
633  if (is_dmdis) {
634  Q2l = kinematics::DarkQ2Lim(Ev,M,ml);
635  return Q2l;
636  }
637 
638  // inelastic
639  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim(Ev,ml,M) : kinematics::InelQ2Lim(Ev,M,ml);
640  return Q2l;
641 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
bool IsWeakCC(void) const
static const double kMinQ2Limit_VLE
Definition: Controls.h:42
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
A simple [min,max] interval for doubles.
Definition: Range1.h:42
int NPions(void) const
Definition: XclsTag.h:62
static const double kPi0Mass
Definition: Constants.h:74
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
int CharmHadronPdg(void) const
Definition: XclsTag.h:52
bool IsStrangeEvent(void) const
Definition: XclsTag.h:53
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:902
bool IsInverseBetaDecay(void) const
Range1D_t CEvNSQ2Lim(double Ev)
Definition: KineUtils.cxx:873
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
Range1D_t InelQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:416
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:366
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
Range1D_t DarkQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:960
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
bool IsCoherentElastic(void) const
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
static const double kASmallNum
Definition: Controls.h:40
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:74
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsMEC(void) const
bool IsEM(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
double max
Definition: Range1.h:53
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
static const double kPionMass
Definition: Constants.h:73
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
Range1D_t CohQ2Lim(double Mn, double m_produced, double mlep, double Ev)
Definition: KineUtils.cxx:730
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
float pi
Definition: units.py:11
const Target & Tgt(void) const
Definition: InitialState.h:66
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
Range1D_t KPhaseSpace::q2Lim ( void  ) const

q2 limits

Definition at line 643 of file KPhaseSpace.cxx.

644 {
645 // As Q2Lim(void) but with reversed sign (Q2 -> q2)
646 //
647  Range1D_t Q2 = this->Q2Lim();
648  Range1D_t q2;
649  q2.min = - Q2.max;
650  q2.max = - Q2.min;
651  return q2;
652 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t Q2Lim(void) const
Q2 limits.
double max
Definition: Range1.h:53
double min
Definition: Range1.h:52
Range1D_t KPhaseSpace::Q2Lim_W ( void  ) const

Q2 limits @ fixed W.

Definition at line 473 of file KPhaseSpace.cxx.

474 {
475  // Computes momentum transfer (Q2>0) limits @ the input invariant mass
476  // The calculation proceeds as in kinematics::InelQ2Lim_W().
477  // For QEL, W is set to the recoil nucleon mass
478  //
479  // TODO: For now, choosing to handle Q2 at fixed W for coherent in the
480  // same way as for the general Q2 limits... but shouldn't we just use
481  // W = m_pi? - which we do in Q2Lim() anyway... seems like there are
482  // cleanup opportunities here.
483 
484  Range1D_t Q2l;
485  Q2l.min = -1;
486  Q2l.max = -1;
487 
488  const ProcessInfo & pi = fInteraction->ProcInfo();
489 
490  bool is_em = pi.IsEM();
491  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay();
492  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive();
493  bool is_coh = pi.IsCoherentProduction();
494  bool is_dme = pi.IsDarkMatterElastic();
495  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
496 
497  if(!is_qel && !is_inel && !is_coh && !is_dme && !is_dmdis) return Q2l;
498 
499  if(is_coh) {
500  return Q2Lim();
501  }
502 
503  const InitialState & init_state = fInteraction->InitState();
504  double Ev = init_state.ProbeE(kRfHitNucRest);
505  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
506  double ml = fInteraction->FSPrimLepton()->Mass();
507 
508  double W = 0;
509  if(is_qel || is_dme) W = fInteraction->RecoilNucleon()->Mass();
510  else W = kinematics::W(fInteraction);
511 
512  if (pi.IsInverseBetaDecay()) {
514  } else if (is_dme || is_dmdis) {
515  Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W);
516  } else {
517  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W);
518  }
519 
520  return Q2l;
521 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
static const double kMinQ2Limit_VLE
Definition: Controls.h:42
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
A simple [min,max] interval for doubles.
Definition: Range1.h:42
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:902
bool IsInverseBetaDecay(void) const
bool IsDiffractive(void) const
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
Range1D_t Q2Lim(void) const
Q2 limits.
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:366
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:74
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
double max
Definition: Range1.h:53
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
float pi
Definition: units.py:11
const Target & Tgt(void) const
Definition: InitialState.h:66
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
Range1D_t KPhaseSpace::q2Lim_W ( void  ) const

q2 limits @ fixed W

Definition at line 523 of file KPhaseSpace.cxx.

524 {
525 // As Q2Lim_W(void) but with reversed sign (Q2 -> q2)
526 //
527  Range1D_t Q2 = this->Q2Lim_W();
528  Range1D_t q2;
529  q2.min = - Q2.max;
530  q2.max = - Q2.min;
531  return q2;
532 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
double max
Definition: Range1.h:53
double min
Definition: Range1.h:52
double KPhaseSpace::Threshold ( void  ) const

Energy threshold.

Definition at line 80 of file KPhaseSpace.cxx.

81 {
82  const ProcessInfo & pi = fInteraction->ProcInfo();
83  const InitialState & init_state = fInteraction->InitState();
84  const XclsTag & xcls = fInteraction->ExclTag();
85  const Target & tgt = init_state.Tgt();
86 
87  double ml = fInteraction->FSPrimLepton()->Mass();
88 
89  if( ! pi.IsKnown() ) return 0;
90 
91  if (pi.IsSingleKaon()) {
92  int kaon_pdgc = xcls.StrangeHadronPdg();
93  double Mi = tgt.HitNucP4Ptr()->M(); // initial nucleon mass
94  // Final nucleon can be different for K0 interaction
95  double Mf = (xcls.NProtons()==1) ? kProtonMass : kNeutronMass;
96  double mk = PDGLibrary::Instance()->Find(kaon_pdgc)->Mass();
97  //double ml = PDGLibrary::Instance()->Find(fInteraction->FSPrimLeptonPdg())->Mass();
98  double mtot = ml + mk + Mf; // total mass of FS particles
99  double Ethresh = (mtot*mtot - Mi*Mi)/(2. * Mf);
100  return Ethresh;
101  }
102 
103  if(pi.IsCoherentElastic()) {
104  return ml + 0.5*ml*ml/tgt.Mass();
105  }
106 
107  if (pi.IsCoherentProduction()) {
108 
109  int tgtpdgc = tgt.Pdg(); // nuclear target PDG code (10LZZZAAAI)
110  double MA = PDGLibrary::Instance()->Find(tgtpdgc)->Mass();
111 
112  double m_other = controls::kASmallNum ;
113  // as a default the mass of hadronic system is the mass of the photon.
114  // which is assumed to be a small number to avoid divergences
115 
116  if ( xcls.NPions() > 0 ) {
117  m_other = pi.IsWeakCC() ? kPionMass : kPi0Mass;
118  }
119 
120  double m = ml + m_other ;
121  double m2 = TMath::Power(m,2);
122  double Ethr = m + 0.5*m2/MA;
123 
124  return TMath::Max(0.,Ethr);
125  }
126 
127  if(pi.IsQuasiElastic() ||
128  pi.IsDarkMatterElastic() ||
129  pi.IsInverseBetaDecay() ||
130  pi.IsResonant() ||
131  pi.IsDeepInelastic() ||
133  pi.IsDiffractive())
134  {
135  assert(tgt.HitNucIsSet());
136  double Mn = tgt.HitNucP4Ptr()->M();
137  double Mn2 = TMath::Power(Mn,2);
138  double Wmin = kNucleonMass + kPionMass;
139  if ( pi.IsQuasiElastic() || pi.IsDarkMatterElastic() || pi.IsInverseBetaDecay() ) {
140  int finalNucPDG = tgt.HitNucPdg();
141  if ( pi.IsWeakCC() ) finalNucPDG = pdg::SwitchProtonNeutron( finalNucPDG );
142  Wmin = PDGLibrary::Instance()->Find( finalNucPDG )->Mass();
143  }
144  if (pi.IsResonant()) {
145  Wmin = kNucleonMass + kPhotontest;
146  }
147 
148  if(xcls.IsCharmEvent()) {
149  if(xcls.IsInclusiveCharm()) {
151  } else {
152  int cpdg = xcls.CharmHadronPdg();
153  double mchm = PDGLibrary::Instance()->Find(cpdg)->Mass();
154  if(pi.IsQuasiElastic() || pi.IsInverseBetaDecay()) {
155  Wmin = mchm + controls::kASmallNum;
156  }
157  else {
158  Wmin = kNeutronMass + mchm + controls::kASmallNum;
159  }
160  }//incl.?
161  }//charm?
162 
163  double smin = TMath::Power(Wmin+ml,2.);
164  double Ethr = 0.5*(smin-Mn2)/Mn;
165  // threshold is different for dark matter case
167  // Correction to minimum kinematic variables
168  Wmin = Mn;
169  smin = TMath::Power(Wmin+ml,2);
170  Ethr = TMath::Max(0.5*(smin-Mn2-ml*ml)/Mn,ml);
171  }
172 
173  return TMath::Max(0.,Ethr);
174  }
175 
176  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation()) {
177  double Ethr = 0.5 * (kMuonMass2-kElectronMass2)/kElectronMass;
178  return TMath::Max(0.,Ethr);
179  }
180 
182  return 0;
183  }
184  if(pi.IsAMNuGamma()) {
185  return 0;
186  }
187  if (pi.IsMEC()) {
188  if (tgt.HitNucIsSet()) {
189  double Mn = tgt.HitNucP4Ptr()->M();
190  double Mn2 = TMath::Power(Mn,2);
191  double Wmin = fInteraction->RecoilNucleon()->Mass(); // mass of the recoil nucleon cluster
192  double smin = TMath::Power(Wmin+ml,2.);
193  double Ethr = 0.5*(smin-Mn2)/Mn;
194  return TMath::Max(0.,Ethr);
195  }
196  else {
197  // this was ... if (pi.IsMECTensor())
198  return ml;
199  }
200  }
201  if(pi.IsGlashowResonance()) {
202  double Ethr = 0.5 * (ml*ml-kElectronMass2)/kElectronMass;
203  return TMath::Max(0.,Ethr);
204  }
205 
206 
207  SLOG("KPhaseSpace", pERROR)
208  << "Can't compute threshold for \n" << *fInteraction;
209  throw genie::exceptions::InteractionException("Can't compute threshold");
210  //exit(1);
211 
212  return 99999999;
213 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
bool IsWeakCC(void) const
#define pERROR
Definition: Messenger.h:59
bool IsDarkMatterElectronElastic(void) const
static const double kNucleonMass
Definition: Constants.h:77
int HitNucPdg(void) const
Definition: Target.cxx:304
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
int NPions(void) const
Definition: XclsTag.h:62
bool IsInverseMuDecay(void) const
static const double kPi0Mass
Definition: Constants.h:74
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
int CharmHadronPdg(void) const
Definition: XclsTag.h:52
int Pdg(void) const
Definition: Target.h:71
static const double kPhotontest
Definition: Constants.h:79
int SwitchProtonNeutron(int pdgc)
Definition: PDGUtils.cxx:353
cpdg
Definition: tracks.py:338
bool IsInverseBetaDecay(void) const
static const double kLightestChmHad
Definition: Constants.h:78
bool IsDiffractive(void) const
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
bool IsIMDAnnihilation(void) const
double Mass(void) const
Definition: Target.cxx:224
static const double kElectronMass
Definition: Constants.h:70
bool IsKnown(void) const
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
bool IsSingleKaon(void) const
Definition: ProcessInfo.cxx:79
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
bool IsCoherentElastic(void) const
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
static const double kElectronMass2
Definition: Constants.h:83
Exception used inside Interaction classes.
bool IsNuElectronElastic(void) const
static constexpr double m2
Definition: Units.h:72
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsAMNuGamma(void) const
static const double kMuonMass2
Definition: Constants.h:84
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
static const double kNeutronMass
Definition: Constants.h:76
static const double kASmallNum
Definition: Controls.h:40
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:74
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsMEC(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
static const double kPionMass
Definition: Constants.h:73
bool HitNucIsSet(void) const
Definition: Target.cxx:283
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:54
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
float pi
Definition: units.py:11
const Target & Tgt(void) const
Definition: InitialState.h:66
bool IsGlashowResonance(void) const
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
static const double kProtonMass
Definition: Constants.h:75
Initial State information.
Definition: InitialState.h:48
int NProtons(void) const
Definition: XclsTag.h:56
Range1D_t KPhaseSpace::TLim ( void  ) const

t limits

Definition at line 888 of file KPhaseSpace.cxx.

889 {
890  // t limits for Coherent pion production from
891  // Kartavtsev, Paschos, and Gounaris, PRD 74 054007, and
892  // Paschos and Schalla, PRD 80, 03305
893  // TODO: Attempt to assign t bounds for other reactions?
894  Range1D_t tl;
895  tl.min = -1;
896  tl.max = -1;
897 
898  const InitialState & init_state = fInteraction->InitState();
899  const ProcessInfo & pi = fInteraction->ProcInfo();
900  const Kinematics & kine = fInteraction->Kine();
902  double Ev = init_state.ProbeE(kRfHitNucRest);
903  double Q2 = kine.Q2();
904  double nu = Ev * kine.y();
905 
906  //COH
907  if(pi.IsCoherentProduction()) {
908 
909  double m_other = controls::kASmallNum ;
910  // as a default the mass of hadronic system is the mass of the photon.
911  // which is assumed to be a small number to avoid divergences
912 
913  const XclsTag & xcls = fInteraction -> ExclTag() ;
914 
915  if ( xcls.NPions() > 0 ) {
916  bool pionIsCharged = pi.IsWeakCC();
917  m_other = pionIsCharged ? kPionMass : kPi0Mass;
918  }
919 
920  double m_other2 = m_other * m_other ;
921 
922  tl.min = 1.0 * (Q2 + m_other2)/(2.0 * nu) * (Q2 + m_other2)/(2.0 * nu);
923  tl.max = 0.05;
924  return tl;
925  }
926  // DFR
927  else if (pi.IsDiffractive()) {
928 
929  // diffractive tmin from Nucl.Phys.B278,61 (1986), eq. 12
930 
931  bool pionIsCharged = pi.IsWeakCC();
932  double mpi = pionIsCharged ? kPionMass : kPi0Mass;
933  double mpi2 = mpi*mpi;
934 
935  double M = init_state.Tgt().HitNucMass();
936  double M2 = M*M;
937  double nuSqPlusQ2 = nu*nu + Q2;
938  double nuOverM = nu / M;
939  double mpiQ2term = mpi2 - Q2 - 2*nu*nu;
940  double A1 = 1 + 2*nuOverM + nuOverM*nuOverM - nuSqPlusQ2/M2;
941  double A2 = (1+nuOverM) * mpiQ2term + 2*nuOverM*nuSqPlusQ2;
942  double A3 = mpiQ2term*mpiQ2term - 4*nuSqPlusQ2*(nu*nu - mpi2);
943 
944  tl.min = std::abs( (A2 + sqrt(A2*A2 - A1*A3)) / A1 ); // GENIE's convention is that t is positive
945  bool tminIsNaN;
946  // use std::isnan when C++11 is around
947 #if __cplusplus >= 201103L
948  tminIsNaN = std::isnan(tl.min);
949 #else
950  // this the old-fashioned way to check for NaN:
951  // NaN's aren't equal to anything, including themselves
952  tminIsNaN = tl.min != tl.min;
953 #endif
954  if (tminIsNaN)
955  {
956  LOG("KPhaseSpace", pERROR)
957  << "tmin for diffractive scattering is NaN "
958  << "( Enu = " << Ev << ", Q2 = " << Q2 << ", nu = " << nu << ")";
959  throw genie::exceptions::InteractionException("NaN tmin for diffractive scattering");
960  }
961  tl.max = this->GetTMaxDFR();
962 
963  return tl;
964  }
965 
966  // RES+DIS
967  // IMD
968  LOG("KPhaseSpace", pWARN) << "It is not sensible to ask for t limits for events that are not coherent or diffractive.";
969  return tl;
970 }
bool IsWeakCC(void) const
#define pERROR
Definition: Messenger.h:59
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
A simple [min,max] interval for doubles.
Definition: Range1.h:42
int NPions(void) const
Definition: XclsTag.h:62
static const double kPi0Mass
Definition: Constants.h:74
double HitNucMass(void) const
Definition: Target.cxx:233
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
bool IsDiffractive(void) const
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
static double GetTMaxDFR()
Definition: KPhaseSpace.cxx:57
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
double y(bool selected=false) const
Definition: Kinematics.cxx:112
T abs(T value)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
Exception used inside Interaction classes.
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const Kinematics & Kine(void) const
Definition: Interaction.h:71
static const double kASmallNum
Definition: Controls.h:40
#define pWARN
Definition: Messenger.h:60
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1277
double max
Definition: Range1.h:53
static const double kPionMass
Definition: Constants.h:73
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
float pi
Definition: units.py:11
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
const Target & Tgt(void) const
Definition: InitialState.h:66
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
void KPhaseSpace::UseInteraction ( const Interaction in)

Definition at line 75 of file KPhaseSpace.cxx.

76 {
77  fInteraction = in;
78 }
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
Range1D_t KPhaseSpace::WLim ( void  ) const

W limits.

Definition at line 403 of file KPhaseSpace.cxx.

404 {
405 // Computes hadronic invariant mass limits.
406 // For QEL the range reduces to the recoil nucleon mass.
407 // For DIS & RES the calculation proceeds as in kinematics::InelWLim().
408 // It is not computed for other interactions
409 //
410  Range1D_t Wl;
411  Wl.min = -1;
412  Wl.max = -1;
413 
414  const ProcessInfo & pi = fInteraction->ProcInfo();
415 
416  bool is_em = pi.IsEM();
417  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic();
418  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive();
419  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
420 
421  if(is_qel) {
422  double MR = fInteraction->RecoilNucleon()->Mass();
423  Wl.min = MR;
424  Wl.max = MR;
425  return Wl;
426  }
427  if(is_inel) {
428  const InitialState & init_state = fInteraction->InitState();
429  double Ev = init_state.ProbeE(kRfHitNucRest);
430  double M = init_state.Tgt().HitNucP4Ptr()->M(); //can be off m/shell
431  double ml = fInteraction->FSPrimLepton()->Mass();
432 
433  Wl = is_em ? kinematics::electromagnetic::InelWLim(Ev,ml,M) : kinematics::InelWLim(Ev,M,ml);
434 
436  //Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass+kLightestChmHad);
437  Wl.min = TMath::Max(Wl.min, kNeutronMass+kLightestChmHad);
438  }
439  else if (fInteraction->ProcInfo().IsDiffractive())
440  Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass);
441  else if ( fInteraction->ProcInfo().IsDeepInelastic() ) {
442  Wl.min = TMath::Max(Wl.min, kNeutronMass + kPionMass );
443  }
444 
445  // sanity check
446  if(Wl.min>Wl.max) {Wl.min=-1; Wl.max=-1;}
447 
448  return Wl;
449  }
450  if(is_dmdis) {
451  const InitialState & init_state = fInteraction->InitState();
452  double Ev = init_state.ProbeE(kRfHitNucRest);
453  double M = init_state.Tgt().HitNucP4Ptr()->M(); //can be off m/shell
454  double ml = fInteraction->FSPrimLepton()->Mass();
455  Wl = kinematics::DarkWLim(Ev,M,ml);
457  //Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass+kLightestChmHad);
458  Wl.min = TMath::Max(Wl.min, kNeutronMass+kLightestChmHad);
459  }
460  else if (fInteraction->ProcInfo().IsDiffractive())
461  Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass);
462 
463  LOG("KPhaseSpace", pDEBUG) << "Found nominal limits: " << Wl.min << ", " << Wl.max;
464 
465  // sanity check
466  if(Wl.min>Wl.max) {Wl.min=-1; Wl.max=-1;}
467 
468  return Wl;
469  }
470  return Wl;
471 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
Range1D_t InelWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:345
A simple [min,max] interval for doubles.
Definition: Range1.h:42
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
bool IsInverseBetaDecay(void) const
static const double kLightestChmHad
Definition: Constants.h:78
bool IsDiffractive(void) const
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
static const double kNeutronMass
Definition: Constants.h:76
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:74
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
double max
Definition: Range1.h:53
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static const double kPionMass
Definition: Constants.h:73
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
float pi
Definition: units.py:11
const Target & Tgt(void) const
Definition: InitialState.h:66
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
double ProbeE(RefFrame_t rf) const
Range1D_t DarkWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:879
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
Range1D_t KPhaseSpace::XLim ( void  ) const

x limits

Definition at line 654 of file KPhaseSpace.cxx.

655 {
656  // Computes x-limits;
657 
658  Range1D_t xl;
659  xl.min = -1;
660  xl.max = -1;
661 
662  const ProcessInfo & pi = fInteraction->ProcInfo();
663  bool is_em = pi.IsEM();
664 
665  //RES+DIS
666  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
667  if(is_inel) {
668  const InitialState & init_state = fInteraction->InitState();
669  double Ev = init_state.ProbeE(kRfHitNucRest);
670  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
671  double ml = fInteraction->FSPrimLepton()->Mass();
672  xl = is_em ? kinematics::electromagnetic::InelXLim(Ev,ml,M) : kinematics::InelXLim(Ev,M,ml);
673  return xl;
674  }
675  //DMDIS
676  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
677  if(is_dmdis) {
678  const InitialState & init_state = fInteraction->InitState();
679  double Ev = init_state.ProbeE(kRfHitNucRest);
680  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
681  double ml = fInteraction->FSPrimLepton()->Mass();
682  xl = kinematics::DarkXLim(Ev,M,ml);
683  return xl;
684  }
685  //COH
686  bool is_coh = pi.IsCoherentProduction();
687  if(is_coh) {
688  xl = kinematics::CohXLim();
689  return xl;
690  }
691  //QEL
692  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic();
693  if(is_qel) {
694  xl.min = 1;
695  xl.max = 1;
696  return xl;
697  }
698  bool is_dfr = pi.IsDiffractive();
699  if(is_dfr) {
701  xl.max = 1. - controls::kASmallNum;
702  return xl;
703  }
704 
705  return xl;
706 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
A simple [min,max] interval for doubles.
Definition: Range1.h:42
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
bool IsInverseBetaDecay(void) const
bool IsDiffractive(void) const
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
Range1D_t InelXLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:444
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
Range1D_t CohXLim(void)
Definition: KineUtils.cxx:722
static const double kASmallNum
Definition: Controls.h:40
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:74
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
double max
Definition: Range1.h:53
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
float pi
Definition: units.py:11
const Target & Tgt(void) const
Definition: InitialState.h:66
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
double ProbeE(RefFrame_t rf) const
Range1D_t DarkXLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:988
Initial State information.
Definition: InitialState.h:48
Range1D_t KPhaseSpace::YLim ( void  ) const

y limits

Definition at line 708 of file KPhaseSpace.cxx.

709 {
710  Range1D_t yl;
711  yl.min = -1;
712  yl.max = -1;
713 
714  const ProcessInfo & pi = fInteraction->ProcInfo();
715  bool is_em = pi.IsEM();
716 
717  //RES+DIS
718  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
719  if(is_inel) {
720  const InitialState & init_state = fInteraction->InitState();
721  double Ev = init_state.ProbeE(kRfHitNucRest);
722  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
723  double ml = fInteraction->FSPrimLepton()->Mass();
724  yl = is_em ? kinematics::electromagnetic::InelYLim(Ev,ml,M) : kinematics::InelYLim(Ev,M,ml);
725  return yl;
726  }
727  //DMDIS
728  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
729  if(is_dmdis) {
730  const InitialState & init_state = fInteraction->InitState();
731  double Ev = init_state.ProbeE(kRfHitNucRest);
732  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
733  double ml = fInteraction->FSPrimLepton()->Mass();
734  yl = kinematics::DarkYLim(Ev,M,ml);
735  return yl;
736  }
737  //COH
738  bool is_coh = pi.IsCoherentProduction();
739  if(is_coh) {
740  const InitialState & init_state = fInteraction->InitState();
741  double EvL = init_state.ProbeE(kRfLab);
742  double ml = fInteraction->FSPrimLepton()->Mass();
743  yl = kinematics::CohYLim(EvL,ml);
744  return yl;
745  }
746  // IMD
747  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation() || pi.IsNuElectronElastic()) {
748  const InitialState & init_state = fInteraction->InitState();
749  double Ev = init_state.ProbeE(kRfLab);
750  double ml = fInteraction->FSPrimLepton()->Mass();
751  double me = kElectronMass;
753  yl.max = 1 - (ml*ml + me*me)/(2*me*Ev) - controls::kASmallNum;
754  return yl;
755  }
756  // EDIT: y limits are different for massive probe
757  if(pi.IsDarkMatterElectronElastic()) {
758  const InitialState & init_state = fInteraction->InitState();
759  double Ev = init_state.ProbeE(kRfLab);
760  double ml = fInteraction->FSPrimLepton()->Mass();
761  double me = kElectronMass;
762  yl.min = (Ev*me*me + ml*ml*(Ev + 2.0*me)) / (Ev * (2.0*Ev*me + me*me + ml*ml)) + controls::kASmallNum;
763  yl.max = 1.0 - controls::kASmallNum;
764  return yl;
765  }
766  bool is_dfr = pi.IsDiffractive();
767  if(is_dfr) {
768  const InitialState & init_state = fInteraction -> InitState();
769  double Ev = init_state.ProbeE(kRfHitNucRest);
770  double ml = fInteraction->FSPrimLepton()->Mass();
772  yl.max = 1. -ml/Ev - controls::kASmallNum;
773  return yl;
774  }
775  // GLRES
776  if(pi.IsGlashowResonance()) {
777  const InitialState & init_state = fInteraction->InitState();
778  double Ev = init_state.ProbeE(kRfLab);
779  double ml = fInteraction->FSPrimLepton()->Mass();
780  double me = kElectronMass;
781  yl.min = (ml*ml+me*me)/2/Ev/me + controls::kASmallNum;
782  yl.max = (4*Ev*(Ev+me) + (ml*ml+me*me))/2/Ev/(me+2*Ev) - controls::kASmallNum;
783  return yl;
784  }
785  return yl;
786 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
bool IsDarkMatterElectronElastic(void) const
A simple [min,max] interval for doubles.
Definition: Range1.h:42
bool IsInverseMuDecay(void) const
bool IsDiffractive(void) const
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
bool IsIMDAnnihilation(void) const
static const double kElectronMass
Definition: Constants.h:70
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
bool IsNuElectronElastic(void) const
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
static const double kASmallNum
Definition: Controls.h:40
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
double max
Definition: Range1.h:53
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static const double kPionMass
Definition: Constants.h:73
Range1D_t CohYLim(double Mn, double m_produced, double mlep, double Ev, double Q2, double xsi)
Definition: KineUtils.cxx:827
Range1D_t DarkYLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:1008
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
float pi
Definition: units.py:11
const Target & Tgt(void) const
Definition: InitialState.h:66
bool IsGlashowResonance(void) const
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
Range1D_t InelYLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:464
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
Range1D_t KPhaseSpace::YLim ( double  xsi) const

y limits (COH)

Definition at line 833 of file KPhaseSpace.cxx.

834 {
835  // Paschos-Schalla xsi parameter for y-limits in COH
836  // From PRD 80, 033005 (2009)
837 
838  Range1D_t yl;
839  yl.min = -1;
840  yl.max = -1;
841 
842  const ProcessInfo & pi = fInteraction->ProcInfo();
843 
844  //COH
845  bool is_coh = pi.IsCoherentProduction();
846  if(is_coh) {
847  const InitialState & init_state = fInteraction->InitState();
848  const Kinematics & kine = fInteraction->Kine();
849  double Ev = init_state.ProbeE(kRfHitNucRest);
850  double Q2 = kine.Q2();
851  double Mn = init_state.Tgt().Mass();
852  double mlep = fInteraction->FSPrimLepton()->Mass();
853 
854  double m_other = controls::kASmallNum ;
855  // as a default the mass of hadronic system is the mass of the photon.
856  // which is assumed to be a small number to avoid divergences
857 
858  const XclsTag & xcls = fInteraction -> ExclTag() ;
859 
860  if ( xcls.NPions() > 0 ) {
861  bool pionIsCharged = pi.IsWeakCC();
862  m_other = pionIsCharged ? kPionMass : kPi0Mass;
863  }
864 
865  yl = kinematics::CohYLim(Mn, m_other, mlep, Ev, Q2, xsi);
866  return yl;
867  } else {
868  return this->YLim();
869  }
870 }
bool IsWeakCC(void) const
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
A simple [min,max] interval for doubles.
Definition: Range1.h:42
int NPions(void) const
Definition: XclsTag.h:62
static const double kPi0Mass
Definition: Constants.h:74
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
Range1D_t YLim(void) const
y limits
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
double Mass(void) const
Definition: Target.cxx:224
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const Kinematics & Kine(void) const
Definition: Interaction.h:71
static const double kASmallNum
Definition: Controls.h:40
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
double max
Definition: Range1.h:53
static const double kPionMass
Definition: Constants.h:73
Range1D_t CohYLim(double Mn, double m_produced, double mlep, double Ev, double Q2, double xsi)
Definition: KineUtils.cxx:827
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
float pi
Definition: units.py:11
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
const Target & Tgt(void) const
Definition: InitialState.h:66
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
Range1D_t KPhaseSpace::YLim_X ( void  ) const

y limits @ fixed x

Definition at line 788 of file KPhaseSpace.cxx.

789 {
790 // Computes kinematical limits for y @ the input x
791 
792  Range1D_t yl;
793  yl.min = -1;
794  yl.max = -1;
795 
796  const ProcessInfo & pi = fInteraction->ProcInfo();
797  bool is_em = pi.IsEM();
798 
799  //RES+DIS
800  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
801  if(is_inel) {
802  const InitialState & init_state = fInteraction->InitState();
803  double Ev = init_state.ProbeE(kRfHitNucRest);
804  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
805  double ml = fInteraction->FSPrimLepton()->Mass();
806  double x = fInteraction->Kine().x();
807  yl = is_em ? kinematics::electromagnetic::InelYLim_X(Ev,ml,M,x) : kinematics::InelYLim_X(Ev,M,ml,x);
808  return yl;
809  }
810  //DMDIS
811  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
812  if(is_dmdis) {
813  const InitialState & init_state = fInteraction->InitState();
814  double Ev = init_state.ProbeE(kRfHitNucRest);
815  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
816  double ml = fInteraction->FSPrimLepton()->Mass();
817  double x = fInteraction->Kine().x();
818  yl = kinematics::DarkYLim_X(Ev,M,ml,x);
819  return yl;
820  }
821  //COH
822  bool is_coh = pi.IsCoherentProduction();
823  if(is_coh) {
824  const InitialState & init_state = fInteraction->InitState();
825  double EvL = init_state.ProbeE(kRfLab);
826  double ml = fInteraction->FSPrimLepton()->Mass();
827  yl = kinematics::CohYLim(EvL,ml);
828  return yl;
829  }
830  return yl;
831 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:94
A simple [min,max] interval for doubles.
Definition: Range1.h:42
double x(bool selected=false) const
Definition: Kinematics.cxx:99
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const Kinematics & Kine(void) const
Definition: Interaction.h:71
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
double max
Definition: Range1.h:53
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
Range1D_t InelYLim_X(double Ev, double M, double ml, double x)
Definition: KineUtils.cxx:499
Range1D_t CohYLim(double Mn, double m_produced, double mlep, double Ev, double Q2, double xsi)
Definition: KineUtils.cxx:827
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
list x
Definition: train.py:276
double min
Definition: Range1.h:52
Range1D_t DarkYLim_X(double Ev, double M, double ml, double x)
Definition: KineUtils.cxx:1024
float pi
Definition: units.py:11
const Target & Tgt(void) const
Definition: InitialState.h:66
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
Range1D_t KPhaseSpace::YLim_X ( double  xsi) const

y limits @ fixed x (COH)

Definition at line 872 of file KPhaseSpace.cxx.

873 {
874  // Paschos-Schalla xsi parameter for y-limits in COH
875  // From PRD 80, 033005 (2009)
876 
877  const ProcessInfo & pi = fInteraction->ProcInfo();
878 
879  //COH
880  bool is_coh = pi.IsCoherentProduction();
881  if(is_coh) {
882  return this->YLim(xsi);
883  } else {
884  return this->YLim_X();
885  }
886 }
Range1D_t YLim_X(void) const
y limits @ fixed x
Range1D_t YLim(void) const
y limits
bool IsCoherentProduction(void) const
Definition: ProcessInfo.cxx:99
const Interaction * fInteraction
Definition: KPhaseSpace.h:73
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
float pi
Definition: units.py:11

Member Data Documentation

const Interaction* genie::KPhaseSpace::fInteraction
private

Definition at line 73 of file KPhaseSpace.h.


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