KPhaseSpace.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  Changes required to implement the GENIE Boosted Dark Matter module
11  were installed by Josh Berger (Univ. of Wisconsin)
12 */
13 //____________________________________________________________________________
14 
15 #include <cmath>
16 #include <cstdlib>
17 
18 #include <TMath.h>
19 
21 
33 
34 using namespace genie;
35 using namespace genie::utils;
36 using namespace genie::constants;
37 
39 
40 //____________________________________________________________________________
41 KPhaseSpace::KPhaseSpace(void) :
42 TObject(), fInteraction(NULL)
43 {
44  this->UseInteraction(0);
45 }
46 //___________________________________________________________________________
48 TObject(), fInteraction(NULL)
49 {
50  this->UseInteraction(in);
51 }
52 //___________________________________________________________________________
54 {
55 
56 }
57 //___________________________________________________________________________
59 {
60  static bool tMaxLoaded = false;
61  static double DFR_tMax = -1;
62 
63  if (!tMaxLoaded)
64  {
66  const Registry * r = confp->CommonList( "Param", "Diffractive" ) ;
67  double tmax = r->GetDouble("DFR-t-max");
68  DFR_tMax = tmax;
69  tMaxLoaded = true;
70  }
71 
72  return DFR_tMax;
73 
74 }
75 //___________________________________________________________________________
77 {
78  fInteraction = in;
79 }
80 //___________________________________________________________________________
81 double KPhaseSpace::Threshold(void) const
82 {
83  const ProcessInfo & pi = fInteraction->ProcInfo();
84  const InitialState & init_state = fInteraction->InitState();
85  const XclsTag & xcls = fInteraction->ExclTag();
86  const Target & tgt = init_state.Tgt();
87 
88  double ml = fInteraction->FSPrimLepton()->Mass();
89 
90  if (pi.IsSingleKaon()) {
91  int kaon_pdgc = xcls.StrangeHadronPdg();
92  double Mi = tgt.HitNucP4Ptr()->M(); // initial nucleon mass
93  // Final nucleon can be different for K0 interaction
94  double Mf = (xcls.NProtons()==1) ? kProtonMass : kNeutronMass;
95  double mk = PDGLibrary::Instance()->Find(kaon_pdgc)->Mass();
96  //double ml = PDGLibrary::Instance()->Find(fInteraction->FSPrimLeptonPdg())->Mass();
97  double mtot = ml + mk + Mf; // total mass of FS particles
98  double Ethresh = (mtot*mtot - Mi*Mi)/(2. * Mf);
99  return Ethresh;
100  }
101 
102  if (pi.IsCoherent()) {
103  int tgtpdgc = tgt.Pdg(); // nuclear target PDG code (10LZZZAAAI)
104  double mpi = pi.IsWeakCC() ? kPionMass : kPi0Mass;
105  double MA = PDGLibrary::Instance()->Find(tgtpdgc)->Mass();
106  double m = ml + mpi;
107  double m2 = TMath::Power(m,2);
108  double Ethr = m + 0.5*m2/MA;
109  return TMath::Max(0.,Ethr);
110  }
111 
112  if(pi.IsQuasiElastic() ||
113  pi.IsDarkMatterElastic() ||
114  pi.IsInverseBetaDecay() ||
115  pi.IsResonant() ||
116  pi.IsDeepInelastic() ||
118  pi.IsDiffractive())
119  {
120  assert(tgt.HitNucIsSet());
121  double Mn = tgt.HitNucP4Ptr()->M();
122  double Mn2 = TMath::Power(Mn,2);
123  double Wmin = kNucleonMass + kPionMass;
124  if ( pi.IsQuasiElastic() || pi.IsDarkMatterElastic() || pi.IsInverseBetaDecay() ) {
125  int finalNucPDG = tgt.HitNucPdg();
126  if ( pi.IsWeakCC() ) finalNucPDG = pdg::SwitchProtonNeutron( finalNucPDG );
127  Wmin = PDGLibrary::Instance()->Find( finalNucPDG )->Mass();
128  }
129  if (pi.IsResonant()) {
130  Wmin = kNucleonMass + kPhotontest;
131  }
132 
133  if(xcls.IsCharmEvent()) {
134  if(xcls.IsInclusiveCharm()) {
136  } else {
137  int cpdg = xcls.CharmHadronPdg();
138  double mchm = PDGLibrary::Instance()->Find(cpdg)->Mass();
139  if(pi.IsQuasiElastic() || pi.IsInverseBetaDecay()) {
140  Wmin = mchm + controls::kASmallNum;
141  }
142  else {
143  Wmin = kNeutronMass + mchm + controls::kASmallNum;
144  }
145  }//incl.?
146  }//charm?
147 
148  double smin = TMath::Power(Wmin+ml,2.);
149  double Ethr = 0.5*(smin-Mn2)/Mn;
150  // threshold is different for dark matter case
152  Ethr = TMath::Max(0.5*(smin-Mn2-ml*ml)/Mn,ml);
153  }
154 
155  return TMath::Max(0.,Ethr);
156  }
157 
158  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation()) {
159  double Ethr = 0.5 * (kMuonMass2-kElectronMass2)/kElectronMass;
160  return TMath::Max(0.,Ethr);
161  }
162 
163  if(pi.IsNuElectronElastic() || pi.IsGlashowResonance() ) {
164  return 0;
165  }
166  if(pi.IsAMNuGamma()) {
167  return 0;
168  }
169  if (pi.IsMEC()) {
170  if (tgt.HitNucIsSet()) {
171  double Mn = tgt.HitNucP4Ptr()->M();
172  double Mn2 = TMath::Power(Mn,2);
173  double Wmin = fInteraction->RecoilNucleon()->Mass(); // mass of the recoil nucleon cluster
174  double smin = TMath::Power(Wmin+ml,2.);
175  double Ethr = 0.5*(smin-Mn2)/Mn;
176  return TMath::Max(0.,Ethr);
177  }
178  else {
179  // this was ... if (pi.IsMECTensor())
180  return ml;
181  }
182  }
183 
184  SLOG("KPhaseSpace", pERROR)
185  << "Can't compute threshold for \n" << *fInteraction;
186  exit(1);
187 
188  return 99999999;
189 }
190 //___________________________________________________________________________
192 {
193  // Compute limits for the input kinematic variable irrespective of any other
194  // relevant kinematical variable
195  //
197 
198  switch(kvar) {
199  case(kKVW) : return this->WLim(); break;
200  case(kKVQ2) : return this->Q2Lim(); break;
201  case(kKVq2) : return this->q2Lim(); break;
202  case(kKVx) : return this->XLim(); break;
203  case(kKVy) : return this->YLim(); break;
204  case(kKVt) : return this->TLim(); break;
205  default:
206  LOG("KPhaseSpace", pERROR)
207  << "Couldn't compute limits for " << KineVar::AsString(kvar);
208  Range1D_t R(-1.,-1);
209  return R;
210  }
211 }
212 //____________________________________________________________________________
213 double KPhaseSpace::Minimum(KineVar_t kvar) const
214 {
215  Range1D_t lim = this->Limits(kvar);
216  return lim.min;
217 }
218 //___________________________________________________________________________
219 double KPhaseSpace::Maximum(KineVar_t kvar) const
220 {
221  Range1D_t lim = this->Limits(kvar);
222  return lim.max;
223 }
224 //___________________________________________________________________________
226 {
227  double E = 0.;
228  double Ethr = this->Threshold();
229 
230  const ProcessInfo & pi = fInteraction->ProcInfo();
231  const InitialState & init_state = fInteraction->InitState();
232 
233  if (pi.IsCoherent() ||
234  pi.IsInverseMuDecay() ||
235  pi.IsIMDAnnihilation() ||
236  pi.IsNuElectronElastic() ||
237  pi.IsMEC())
238  {
239  E = init_state.ProbeE(kRfLab);
240  }
241 
242  if(pi.IsQuasiElastic() ||
243  pi.IsDarkMatterElastic() ||
244  pi.IsInverseBetaDecay() ||
245  pi.IsResonant() ||
246  pi.IsDeepInelastic() ||
248  pi.IsDiffractive() ||
249  pi.IsSingleKaon() ||
250  pi.IsAMNuGamma())
251  {
252  E = init_state.ProbeE(kRfHitNucRest);
253  }
254 
255  LOG("KPhaseSpace", pDEBUG) << "E = " << E << ", Ethr = " << Ethr;
256  return (E>Ethr);
257 }
258 //___________________________________________________________________________
259 bool KPhaseSpace::IsAllowed(void) const
260 {
261  const ProcessInfo & pi = fInteraction->ProcInfo();
262  const Kinematics & kine = fInteraction->Kine();
263 
264  // ASK single kaon:
265  // XSec code returns zero when kinematics are not allowed
266  // Here just let kinematics always be allowed
267  if(pi.IsSingleKaon()) {
268  return true;
269  }
270 
271  // QEL:
272  // Check the running Q2 vs the Q2 limits
273  if(pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic()) {
274  Range1D_t Q2l = this->Q2Lim();
275  double Q2 = kine.Q2();
276  bool in_phys = math::IsWithinLimits(Q2, Q2l);
277  bool allowed = in_phys;
278  return allowed;
279  }
280 
281  // RES
282  // Check the running W vs the W limits
283  // & the running Q2 vs Q2 limits for the given W
284  if(pi.IsResonant()) {
285  Range1D_t Wl = this->WLim();
286  Range1D_t Q2l = this->Q2Lim_W();
287  double W = kine.W();
288  double Q2 = kine.Q2();
289  bool in_phys = (math::IsWithinLimits(Q2, Q2l) && math::IsWithinLimits(W, Wl));
290  bool allowed = in_phys;
291  return allowed;
292  }
293 
294  // DIS
295  if(pi.IsDeepInelastic() || pi.IsDarkMatterDeepInelastic()) {
296  Range1D_t Wl = this->WLim();
297  Range1D_t Q2l = this->Q2Lim_W();
298  double W = kine.W();
299  double Q2 = kine.Q2();
300  bool in_phys = (math::IsWithinLimits(Q2, Q2l) && math::IsWithinLimits(W, Wl));
301  bool allowed = in_phys;
302  return allowed;
303  }
304 
305  //IMD
306  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation() || pi.IsNuElectronElastic()) {
307  Range1D_t yl = this->YLim();
308  double y = kine.y();
309  bool in_phys = math::IsWithinLimits(y, yl);
310  bool allowed = in_phys;
311  return allowed;
312  }
313 
314  //COH
315  if (pi.IsCoherent()) {
316  Range1D_t xl = this->XLim();
317  Range1D_t yl = this->YLim();
318  double x = kine.x();
319  double y = kine.y();
320  bool in_phys = (math::IsWithinLimits(x, xl) && math::IsWithinLimits(y, yl));
321  bool allowed = in_phys;
322  return allowed;
323  }
324 
325  // DFR
326  if (pi.IsDiffractive()) {
327  // first two checks are the same as RES & DIS
328  Range1D_t Wl = this->WLim();
329  Range1D_t Q2l = this->Q2Lim_W();
330 
332  double W = kine.W();
333  double Q2 = kine.Q2();
334 
335  LOG("KPhaseSpace", pDEBUG) << " W = " << W << ", limits = [" << Wl.min << "," << Wl.max << "];";
336  LOG("KPhaseSpace", pDEBUG) << " Q2 = " << Q2 << ", limits = [" << Q2l.min << "," << Q2l.max << "];";
337  bool in_phys = math::IsWithinLimits(W, Wl);
338  in_phys = in_phys && math::IsWithinLimits(Q2, Q2l);
339 
340  // extra check: there's a t minimum.
341  // but only check if W, Q2 is reasonable
342  // (otherwise get NaNs in tmin)
343  if (in_phys)
344  {
345  double t = kine.t();
346  Range1D_t tl = this->TLim();
347  LOG("KPhaseSpace", pDEBUG) << " t = " << t << ", limits = [" << tl.min << "," << tl.max << "];";
348  in_phys = in_phys && math::IsWithinLimits(t, tl);
349  }
350  LOG("KPhaseSpace", pDEBUG) << " phase space point is " << ( in_phys ? "ALLOWED" : "NOT ALLOWED");
351 
352 
353  bool allowed = in_phys;
354  return allowed;
355  }
356 
357  // was MECTensor
358  if (pi.IsMEC()){
359  Range1D_t Q2l = this->Q2Lim();
360  double Q2 = kine.Q2();
361  bool in_phys = math::IsWithinLimits(Q2, Q2l);
362  bool allowed = in_phys;
363  return allowed;
364  }
365 
366  return false;
367 }
368 //___________________________________________________________________________
370 {
371 // Computes hadronic invariant mass limits.
372 // For QEL the range reduces to the recoil nucleon mass.
373 // For DIS & RES the calculation proceeds as in kinematics::InelWLim().
374 // It is not computed for other interactions
375 //
376  Range1D_t Wl;
377  Wl.min = -1;
378  Wl.max = -1;
379 
380  const ProcessInfo & pi = fInteraction->ProcInfo();
381 
382  bool is_em = pi.IsEM();
383  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic();
384  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive();
385  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
386 
387  if(is_qel) {
388  double MR = fInteraction->RecoilNucleon()->Mass();
389  Wl.min = MR;
390  Wl.max = MR;
391  return Wl;
392  }
393  if(is_inel) {
394  const InitialState & init_state = fInteraction->InitState();
395  double Ev = init_state.ProbeE(kRfHitNucRest);
396  double M = init_state.Tgt().HitNucP4Ptr()->M(); //can be off m/shell
397  double ml = fInteraction->FSPrimLepton()->Mass();
398 
399  Wl = is_em ? kinematics::electromagnetic::InelWLim(Ev,ml,M) : kinematics::InelWLim(Ev,M,ml);
400 
402  //Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass+kLightestChmHad);
403  Wl.min = TMath::Max(Wl.min, kNeutronMass+kLightestChmHad);
404  }
405  else if (fInteraction->ProcInfo().IsDiffractive())
406  Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass);
407  else if ( fInteraction->ProcInfo().IsDeepInelastic() ) {
408  Wl.min = TMath::Max(Wl.min, kNeutronMass + kPionMass );
409  }
410 
411  // sanity check
412  if(Wl.min>Wl.max) {Wl.min=-1; Wl.max=-1;}
413 
414  return Wl;
415  }
416  if(is_dmdis) {
417  const InitialState & init_state = fInteraction->InitState();
418  double Ev = init_state.ProbeE(kRfHitNucRest);
419  double M = init_state.Tgt().HitNucP4Ptr()->M(); //can be off m/shell
420  double ml = fInteraction->FSPrimLepton()->Mass();
421  Wl = kinematics::DarkWLim(Ev,M,ml);
423  //Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass+kLightestChmHad);
424  Wl.min = TMath::Max(Wl.min, kNeutronMass+kLightestChmHad);
425  }
426  else if (fInteraction->ProcInfo().IsDiffractive())
427  Wl.min = TMath::Max(Wl.min, kNeutronMass+kPionMass);
428 
429  LOG("KPhaseSpace", pDEBUG) << "Found nominal limits: " << Wl.min << ", " << Wl.max;
430 
431  // sanity check
432  if(Wl.min>Wl.max) {Wl.min=-1; Wl.max=-1;}
433 
434  return Wl;
435  }
436  return Wl;
437 }
438 //____________________________________________________________________________
440 {
441  // Computes momentum transfer (Q2>0) limits @ the input invariant mass
442  // The calculation proceeds as in kinematics::InelQ2Lim_W().
443  // For QEL, W is set to the recoil nucleon mass
444  //
445  // TODO: For now, choosing to handle Q2 at fixed W for coherent in the
446  // same way as for the general Q2 limits... but shouldn't we just use
447  // W = m_pi? - which we do in Q2Lim() anyway... seems like there are
448  // cleanup opportunities here.
449 
450  Range1D_t Q2l;
451  Q2l.min = -1;
452  Q2l.max = -1;
453 
454  const ProcessInfo & pi = fInteraction->ProcInfo();
455 
456  bool is_em = pi.IsEM();
457  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay();
458  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive();
459  bool is_coh = pi.IsCoherent();
460  bool is_dme = pi.IsDarkMatterElastic();
461  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
462 
463  if(!is_qel && !is_inel && !is_coh && !is_dme && !is_dmdis) return Q2l;
464 
465  if(is_coh) {
466  return Q2Lim();
467  }
468 
469  const InitialState & init_state = fInteraction->InitState();
470  double Ev = init_state.ProbeE(kRfHitNucRest);
471  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
472  double ml = fInteraction->FSPrimLepton()->Mass();
473 
474  double W = 0;
475  if(is_qel || is_dme) W = fInteraction->RecoilNucleon()->Mass();
476  else W = kinematics::W(fInteraction);
477 
478  if (pi.IsInverseBetaDecay()) {
480  } else if (is_dme || is_dmdis) {
481  Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W);
482  } else {
483  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : Q2l = kinematics::InelQ2Lim_W(Ev,M,ml,W);
484  }
485 
486  return Q2l;
487 }
488 //____________________________________________________________________________
490 {
491 // As Q2Lim_W(void) but with reversed sign (Q2 -> q2)
492 //
493  Range1D_t Q2 = this->Q2Lim_W();
494  Range1D_t q2;
495  q2.min = - Q2.max;
496  q2.max = - Q2.min;
497  return q2;
498 }
499 //____________________________________________________________________________
501 {
502  // Computes momentum transfer (Q2>0) limits irrespective of the invariant mass
503  // For QEL this is identical to Q2Lim_W (since W is fixed)
504  // For RES & DIS, the calculation proceeds as in kinematics::InelQ2Lim().
505  //
506  Range1D_t Q2l;
507  Q2l.min = -1;
508  Q2l.max = -1;
509 
510  const ProcessInfo & pi = fInteraction->ProcInfo();
511 
512  bool is_em = pi.IsEM();
513  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay();
514  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
515  bool is_coh = pi.IsCoherent();
516  bool is_dme = pi.IsDarkMatterElastic();
517  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
518 
519  if(!is_qel && !is_inel && !is_coh && !is_dme && !is_dmdis) return Q2l;
520 
521  const InitialState & init_state = fInteraction->InitState();
522  double Ev = init_state.ProbeE(kRfHitNucRest);
523  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
524  double ml = fInteraction->FSPrimLepton()->Mass();
525 
526  if(is_coh) {
527  bool pionIsCharged = pi.IsWeakCC();
528  double mpi = pionIsCharged ? kPionMass : kPi0Mass;
529  Q2l = kinematics::CohQ2Lim(M, mpi, ml, Ev);
530  return Q2l;
531  }
532 
533  const XclsTag & xcls = fInteraction->ExclTag();
534 
535  // quasi-elastic
536  if(is_qel) {
537  double W = fInteraction->RecoilNucleon()->Mass();
538  if(xcls.IsCharmEvent()) {
539  int charm_pdgc = xcls.CharmHadronPdg();
540  W = PDGLibrary::Instance()->Find(charm_pdgc)->Mass();
541  } else if(xcls.IsStrangeEvent()) {
542  int strange_pdgc = xcls.StrangeHadronPdg();
543  W = PDGLibrary::Instance()->Find(strange_pdgc)->Mass();
544  }
545  if (pi.IsInverseBetaDecay()) {
547  } else {
548  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W);
549  }
550 
551  return Q2l;
552  }
553 
554  // dark mattter elastic
555  if(is_dme) {
556  double W = fInteraction->RecoilNucleon()->Mass();
557  if(xcls.IsCharmEvent()) {
558  int charm_pdgc = xcls.CharmHadronPdg();
559  W = PDGLibrary::Instance()->Find(charm_pdgc)->Mass();
560  } else if(xcls.IsStrangeEvent()) {
561  int strange_pdgc = xcls.StrangeHadronPdg();
562  W = PDGLibrary::Instance()->Find(strange_pdgc)->Mass();
563  }
564  if (pi.IsInverseBetaDecay()) {
566  } else {
567  Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W);
568  }
569 
570 
571  return Q2l;
572  }
573 
574  // was MECTensor
575  // TODO: Q2maxConfig
576  if (pi.IsMEC()){
577  double W = fInteraction->RecoilNucleon()->Mass();
578  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W);
579  double Q2maxConfig = 1.44; // need to pull from config file somehow?
580  if (Q2l.max > Q2maxConfig) Q2l.max = Q2maxConfig;
581  return Q2l;
582  }
583 
584  if (is_dmdis) {
585  Q2l = kinematics::DarkQ2Lim(Ev,M,ml);
586  return Q2l;
587  }
588 
589  // inelastic
590  Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim(Ev,ml,M) : kinematics::InelQ2Lim(Ev,M,ml);
591  return Q2l;
592 }
593 //____________________________________________________________________________
595 {
596 // As Q2Lim(void) but with reversed sign (Q2 -> q2)
597 //
598  Range1D_t Q2 = this->Q2Lim();
599  Range1D_t q2;
600  q2.min = - Q2.max;
601  q2.max = - Q2.min;
602  return q2;
603 }
604 //____________________________________________________________________________
606 {
607  // Computes x-limits;
608 
609  Range1D_t xl;
610  xl.min = -1;
611  xl.max = -1;
612 
613  const ProcessInfo & pi = fInteraction->ProcInfo();
614  bool is_em = pi.IsEM();
615 
616  //RES+DIS
617  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
618  if(is_inel) {
619  const InitialState & init_state = fInteraction->InitState();
620  double Ev = init_state.ProbeE(kRfHitNucRest);
621  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
622  double ml = fInteraction->FSPrimLepton()->Mass();
623  xl = is_em ? kinematics::electromagnetic::InelXLim(Ev,ml,M) : kinematics::InelXLim(Ev,M,ml);
624  return xl;
625  }
626  //DMDIS
627  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
628  if(is_dmdis) {
629  const InitialState & init_state = fInteraction->InitState();
630  double Ev = init_state.ProbeE(kRfHitNucRest);
631  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
632  double ml = fInteraction->FSPrimLepton()->Mass();
633  xl = kinematics::DarkXLim(Ev,M,ml);
634  return xl;
635  }
636  //COH
637  bool is_coh = pi.IsCoherent();
638  if(is_coh) {
639  xl = kinematics::CohXLim();
640  return xl;
641  }
642  //QEL
643  bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic();
644  if(is_qel) {
645  xl.min = 1;
646  xl.max = 1;
647  return xl;
648  }
649  bool is_dfr = pi.IsDiffractive();
650  if(is_dfr) {
652  xl.max = 1. - controls::kASmallNum;
653  return xl;
654  }
655 
656  return xl;
657 }
658 //____________________________________________________________________________
660 {
661  Range1D_t yl;
662  yl.min = -1;
663  yl.max = -1;
664 
665  const ProcessInfo & pi = fInteraction->ProcInfo();
666  bool is_em = pi.IsEM();
667 
668  //RES+DIS
669  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
670  if(is_inel) {
671  const InitialState & init_state = fInteraction->InitState();
672  double Ev = init_state.ProbeE(kRfHitNucRest);
673  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
674  double ml = fInteraction->FSPrimLepton()->Mass();
675  yl = is_em ? kinematics::electromagnetic::InelYLim(Ev,ml,M) : kinematics::InelYLim(Ev,M,ml);
676  return yl;
677  }
678  //DMDIS
679  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
680  if(is_dmdis) {
681  const InitialState & init_state = fInteraction->InitState();
682  double Ev = init_state.ProbeE(kRfHitNucRest);
683  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
684  double ml = fInteraction->FSPrimLepton()->Mass();
685  yl = kinematics::DarkYLim(Ev,M,ml);
686  return yl;
687  }
688  //COH
689  bool is_coh = pi.IsCoherent();
690  if(is_coh) {
691  const InitialState & init_state = fInteraction->InitState();
692  double EvL = init_state.ProbeE(kRfLab);
693  double ml = fInteraction->FSPrimLepton()->Mass();
694  yl = kinematics::CohYLim(EvL,ml);
695  return yl;
696  }
697  // IMD
698  if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation() || pi.IsNuElectronElastic()) {
699  const InitialState & init_state = fInteraction->InitState();
700  double Ev = init_state.ProbeE(kRfLab);
701  double ml = fInteraction->FSPrimLepton()->Mass();
702  double me = kElectronMass;
704  yl.max = 1 - (ml*ml + me*me)/(2*me*Ev) - controls::kASmallNum;
705  return yl;
706  }
707  bool is_dfr = pi.IsDiffractive();
708  if(is_dfr) {
709  const InitialState & init_state = fInteraction -> InitState();
710  double Ev = init_state.ProbeE(kRfHitNucRest);
711  double ml = fInteraction->FSPrimLepton()->Mass();
713  yl.max = 1. -ml/Ev - controls::kASmallNum;
714  return yl;
715  }
716  return yl;
717 }
718 //____________________________________________________________________________
720 {
721 // Computes kinematical limits for y @ the input x
722 
723  Range1D_t yl;
724  yl.min = -1;
725  yl.max = -1;
726 
727  const ProcessInfo & pi = fInteraction->ProcInfo();
728  bool is_em = pi.IsEM();
729 
730  //RES+DIS
731  bool is_inel = pi.IsDeepInelastic() || pi.IsResonant();
732  if(is_inel) {
733  const InitialState & init_state = fInteraction->InitState();
734  double Ev = init_state.ProbeE(kRfHitNucRest);
735  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
736  double ml = fInteraction->FSPrimLepton()->Mass();
737  double x = fInteraction->Kine().x();
738  yl = is_em ? kinematics::electromagnetic::InelYLim_X(Ev,ml,M,x) : kinematics::InelYLim_X(Ev,M,ml,x);
739  return yl;
740  }
741  //DMDIS
742  bool is_dmdis = pi.IsDarkMatterDeepInelastic();
743  if(is_dmdis) {
744  const InitialState & init_state = fInteraction->InitState();
745  double Ev = init_state.ProbeE(kRfHitNucRest);
746  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
747  double ml = fInteraction->FSPrimLepton()->Mass();
748  double x = fInteraction->Kine().x();
749  yl = kinematics::DarkYLim_X(Ev,M,ml,x);
750  return yl;
751  }
752  //COH
753  bool is_coh = pi.IsCoherent();
754  if(is_coh) {
755  const InitialState & init_state = fInteraction->InitState();
756  double EvL = init_state.ProbeE(kRfLab);
757  double ml = fInteraction->FSPrimLepton()->Mass();
758  yl = kinematics::CohYLim(EvL,ml);
759  return yl;
760  }
761  return yl;
762 }
763 //____________________________________________________________________________
764 Range1D_t KPhaseSpace::YLim(double xsi) const
765 {
766  // Paschos-Schalla xsi parameter for y-limits in COH
767  // From PRD 80, 033005 (2009)
768 
769  Range1D_t yl;
770  yl.min = -1;
771  yl.max = -1;
772 
773  const ProcessInfo & pi = fInteraction->ProcInfo();
774 
775  //COH
776  bool is_coh = pi.IsCoherent();
777  if(is_coh) {
778  const InitialState & init_state = fInteraction->InitState();
779  const Kinematics & kine = fInteraction->Kine();
780  double Ev = init_state.ProbeE(kRfHitNucRest);
781  double Q2 = kine.Q2();
782  bool pionIsCharged = pi.IsWeakCC();
783  double Mn = init_state.Tgt().Mass();
784  double mpi = pionIsCharged ? kPionMass : kPi0Mass;
785  double mlep = fInteraction->FSPrimLepton()->Mass();
786  yl = kinematics::CohYLim(Mn, mpi, mlep, Ev, Q2, xsi);
787  return yl;
788  } else {
789  return this->YLim();
790  }
791 }
792 //____________________________________________________________________________
794 {
795  // Paschos-Schalla xsi parameter for y-limits in COH
796  // From PRD 80, 033005 (2009)
797 
798  const ProcessInfo & pi = fInteraction->ProcInfo();
799 
800  //COH
801  bool is_coh = pi.IsCoherent();
802  if(is_coh) {
803  return this->YLim(xsi);
804  } else {
805  return this->YLim_X();
806  }
807 }
808 //____________________________________________________________________________
810 {
811  // t limits for Coherent pion production from
812  // Kartavtsev, Paschos, and Gounaris, PRD 74 054007, and
813  // Paschos and Schalla, PRD 80, 03305
814  // TODO: Attempt to assign t bounds for other reactions?
815  Range1D_t tl;
816  tl.min = -1;
817  tl.max = -1;
818 
819  const InitialState & init_state = fInteraction->InitState();
820  const ProcessInfo & pi = fInteraction->ProcInfo();
821  const Kinematics & kine = fInteraction->Kine();
823  double Ev = init_state.ProbeE(kRfHitNucRest);
824  double Q2 = kine.Q2();
825  double nu = Ev * kine.y();
826  bool pionIsCharged = pi.IsWeakCC();
827  double mpi = pionIsCharged ? kPionMass : kPi0Mass;
828  double mpi2 = mpi*mpi;
829 
830  //COH
831  if(pi.IsCoherent()) {
832  tl.min = 1.0 * (Q2 + mpi2)/(2.0 * nu) * (Q2 + mpi2)/(2.0 * nu);
833  tl.max = 0.05;
834  return tl;
835  }
836  // DFR
837  else if (pi.IsDiffractive()) {
838  // diffractive tmin from Nucl.Phys.B278,61 (1986), eq. 12
839  double M = init_state.Tgt().HitNucMass();
840  double M2 = M*M;
841  double nuSqPlusQ2 = nu*nu + Q2;
842  double nuOverM = nu / M;
843  double mpiQ2term = mpi2 - Q2 - 2*nu*nu;
844  double A1 = 1 + 2*nuOverM + nuOverM*nuOverM - nuSqPlusQ2/M2;
845  double A2 = (1+nuOverM) * mpiQ2term + 2*nuOverM*nuSqPlusQ2;
846  double A3 = mpiQ2term*mpiQ2term - 4*nuSqPlusQ2*(nu*nu - mpi2);
847 
848  tl.min = std::abs( (A2 + sqrt(A2*A2 - A1*A3)) / A1 ); // GENIE's convention is that t is positive
849  bool tminIsNaN;
850  // use std::isnan when C++11 is around
851 #if __cplusplus >= 201103L
852  tminIsNaN = std::isnan(tl.min);
853 #else
854  // this the old-fashioned way to check for NaN:
855  // NaN's aren't equal to anything, including themselves
856  tminIsNaN = tl.min != tl.min;
857 #endif
858  if (tminIsNaN)
859  {
860  LOG("KPhaseSpace", pERROR)
861  << "tmin for diffractive scattering is NaN "
862  << "( Enu = " << Ev << ", Q2 = " << Q2 << ", nu = " << nu << ")";
863  throw genie::exceptions::InteractionException("NaN tmin for diffractive scattering");
864  }
865  tl.max = this->GetTMaxDFR();
866 
867  return tl;
868  }
869 
870  // RES+DIS
871  // IMD
872  LOG("KPhaseSpace", pWARN) << "It is not sensible to ask for t limits for events that are not coherent or diffractive.";
873  return tl;
874 }
875 //____________________________________________________________________________
bool IsResonant(void) const
Definition: ProcessInfo.cxx:92
double W(bool selected=false) const
Definition: Kinematics.cxx:167
Basic constants.
const Var kPi0Mass
bool IsWeakCC(void) const
static const double kMinQ2Limit_VLE
Definition: Controls.h:42
Range1D_t YLim_X(void) const
y limits @ fixed x
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int HitNucPdg(void) const
Definition: Target.cxx:321
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
Range1D_t InelWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:288
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double Threshold(void) const
Energy threshold.
Definition: KPhaseSpace.cxx:81
bool IsInverseMuDecay(void) const
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
Range1D_t q2Lim(void) const
q2 limits
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:41
double HitNucMass(void) const
Definition: Target.cxx:250
int CharmHadronPdg(void) const
Definition: XclsTag.h:50
T sqrt(T number)
Definition: d0nt_math.hpp:156
bool IsStrangeEvent(void) const
Definition: XclsTag.h:51
int Pdg(void) const
Definition: Target.h:72
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:829
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:265
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
static const double kPhotontest
Definition: Constants.h:80
int SwitchProtonNeutron(int pdgc)
Definition: PDGUtils.cxx:319
Range1D_t CohYLim(double Mn, double mpi, double mlep, double Ev, double Q2, double xsi)
Definition: KineUtils.cxx:762
bool IsInverseBetaDecay(void) const
double x(bool selected=false) const
Definition: Kinematics.cxx:109
static const double kLightestChmHad
Definition: Constants.h:79
Range1D_t YLim(void) const
y limits
bool IsDiffractive(void) const
bool IsIMDAnnihilation(void) const
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:489
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
Range1D_t InelQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:359
Range1D_t Q2Lim(void) const
Q2 limits.
double Mass(void) const
Definition: Target.cxx:241
static const double kElectronMass
Definition: Constants.h:71
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
float abs(float number)
Definition: d0nt_math.hpp:39
Registry * CommonList(const string &file_id, const string &set_name) const
static double GetTMaxDFR()
Definition: KPhaseSpace.cxx:58
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:37
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:309
double y(bool selected=false) const
Definition: Kinematics.cxx:122
bool IsCharmEvent(void) const
Definition: XclsTag.h:48
bool IsSingleKaon(void) const
Definition: ProcessInfo.cxx:77
Double_t q2[12][num]
Definition: f2_nu.C:137
Summary information for an interaction.
Definition: Interaction.h:56
Range1D_t InelXLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:387
Range1D_t DarkQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:887
Definition: NueSkimmer.h:24
int StrangeHadronPdg(void) const
Definition: XclsTag.h:53
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const Interaction * fInteraction
Definition: KPhaseSpace.h:74
static const double kElectronMass2
Definition: Constants.h:84
Exception used inside Interaction classes.
bool IsNuElectronElastic(void) const
Float_t E
Definition: plot.C:20
ClassImp(KPhaseSpace) KPhaseSpace
Definition: KPhaseSpace.cxx:38
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
bool IsAMNuGamma(void) const
static const double kMuonMass2
Definition: Constants.h:85
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
Kinematical phase space.
Definition: KPhaseSpace.h:34
static const double kNeutronMass
Definition: Constants.h:77
Range1D_t CohQ2Lim(double Mn, double mpi, double mlep, double Ev)
Definition: KineUtils.cxx:673
#define R(x)
Range1D_t CohXLim(void)
Definition: KineUtils.cxx:665
static const double kASmallNum
Definition: Controls.h:40
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:72
double Minimum(KineVar_t kvar) const
#define pWARN
Definition: Messenger.h:61
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsMEC(void) const
bool IsEM(void) const
enum genie::EKineVar KineVar_t
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:82
double Maximum(KineVar_t kvar) const
void UseInteraction(const Interaction *in)
Definition: KPhaseSpace.cxx:76
static constexpr Double_t m2
Definition: Munits.h:145
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1185
double max
Definition: Range1.h:54
ifstream in
Definition: comparison.C:7
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static const double kPionMass
Definition: Constants.h:74
bool IsAllowed(void) const
Check whether the current kinematics is in the allowed phase space.
bool HitNucIsSet(void) const
Definition: Target.cxx:300
Range1D_t InelYLim_X(double Ev, double M, double ml, double x)
Definition: KineUtils.cxx:442
Range1D_t DarkYLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:935
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
static string AsString(KineVar_t kv)
Definition: KineVar.h:68
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:63
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
Range1D_t XLim(void) const
x limits
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double Wmin[kNWBins]
double min
Definition: Range1.h:53
double t(bool selected=false) const
Definition: Kinematics.cxx:180
Range1D_t DarkYLim_X(double Ev, double M, double ml, double x)
Definition: KineUtils.cxx:951
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
const Target & Tgt(void) const
Definition: InitialState.h:67
Range1D_t TLim(void) const
t limits
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:85
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:87
#define W(x)
Range1D_t InelYLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:407
double ProbeE(RefFrame_t rf) const
Range1D_t DarkXLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:915
Range1D_t q2Lim_W(void) const
q2 limits @ fixed W
static const double kProtonMass
Definition: Constants.h:76
Range1D_t DarkWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:806
Root of GENIE utility namespaces.
Range1D_t WLim(void) const
W limits.
static AlgConfigPool * Instance()
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97
int NProtons(void) const
Definition: XclsTag.h:54