GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
genie::COHKinematicsGenerator Class Reference

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

#include <COHKinematicsGenerator.h>

Inheritance diagram for genie::COHKinematicsGenerator:
Inheritance graph
[legend]
Collaboration diagram for genie::COHKinematicsGenerator:
Collaboration graph
[legend]

Public Member Functions

 COHKinematicsGenerator ()
 
 COHKinematicsGenerator (string config)
 
 ~COHKinematicsGenerator ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
void LoadConfig (void)
 
void CalculateKin_ReinSehgal (GHepRecord *event_rec) const
 
void CalculateKin_BergerSehgal (GHepRecord *event_rec) const
 
void CalculateKin_BergerSehgalFM (GHepRecord *event_rec) const
 
void CalculateKin_AlvarezRuso (GHepRecord *event_rec) const
 
void SetKinematics (const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction, Kinematics *kinematics) const
 
bool CheckKinematics (const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction) const
 
double ComputeMaxXSec (const Interaction *in) const
 
double MaxXSec_ReinSehgal (const Interaction *in) const
 
double MaxXSec_BergerSehgal (const Interaction *in) const
 
double MaxXSec_BergerSehgalFM (const Interaction *in) const
 
double MaxXSec_AlvarezRuso (const Interaction *in) const
 
double Energy (const Interaction *in) const
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Public Attributes

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

Private Member Functions

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

Private Attributes

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

Additional Inherited Members

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

Detailed Description

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

Author
Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool
Created:
October 03, 2004
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 30 of file COHKinematicsGenerator.h.

Constructor & Destructor Documentation

COHKinematicsGenerator::COHKinematicsGenerator ( )

Definition at line 42 of file COHKinematicsGenerator.cxx.

References fEnvelope.

42  :
43  KineGeneratorWithCache("genie::COHKinematicsGenerator")
44 {
45  fEnvelope = 0;
46 }
TF2 * fEnvelope
2-D envelope used for importance sampling
COHKinematicsGenerator::COHKinematicsGenerator ( string  config)

Definition at line 48 of file COHKinematicsGenerator.cxx.

References fEnvelope.

48  :
49  KineGeneratorWithCache("genie::COHKinematicsGenerator", config)
50 {
51  fEnvelope = 0;
52 }
TF2 * fEnvelope
2-D envelope used for importance sampling
COHKinematicsGenerator::~COHKinematicsGenerator ( )

Definition at line 54 of file COHKinematicsGenerator.cxx.

References fEnvelope.

55 {
56  if(fEnvelope) delete fEnvelope;
57 }
TF2 * fEnvelope
2-D envelope used for importance sampling

Member Function Documentation

void COHKinematicsGenerator::CalculateKin_AlvarezRuso ( GHepRecord event_rec) const

Definition at line 442 of file COHKinematicsGenerator.cxx.

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

Referenced by ProcessEventRecord().

443 {
444 
445  LOG("COHKinematics", pNOTICE) << "Using AlvarezRuso Model";
446  // Get the Primary Interacton object
447  Interaction * interaction = evrec->Summary();
448  interaction->SetBit(kISkipProcessChk);
449  interaction->SetBit(kISkipKinematicChk);
450 
451  // Initialise a random number generator
452  RandomGen * rnd = RandomGen::Instance();
453 
454  //-- For the subsequent kinematic selection with the rejection method:
455  // Calculate the max differential cross section or retrieve it from the
456  // cache. Throw an exception and quit the evg thread if a non-positive
457  // value is found.
458  // If the kinematics are generated uniformly over the allowed phase
459  // space the max xsec is irrelevant
460  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
461 
462  //Set up limits of integration variables
463  // Primary lepton energy
464  const double E_l_min = interaction->FSPrimLepton()->Mass();
465  const double E_l_max = interaction->InitStatePtr()->GetProbeP4(kRfLab)->E() - kPionMass;
466  // Primary lepton angle with respect to the beam axis
467  const double ctheta_l_min = 0.4;
468  const double ctheta_l_max = 1.0 - kASmallNum;
469  // Pion angle with respect to the beam axis
470  const double ctheta_pi_min = 0.4;
471  const double ctheta_pi_max = 1.0 - kASmallNum;
472  // Pion angle transverse to the beam axis
473  const double phi_min = kASmallNum;
474  const double phi_max = (2.0 * kPi) - kASmallNum;
475  //
476  const double d_E_l = E_l_max - E_l_min;
477  const double d_ctheta_l = ctheta_l_max - ctheta_l_min;
478  const double d_ctheta_pi = ctheta_pi_max - ctheta_pi_min;
479  const double d_phi = phi_max - phi_min;
480 
481  //------ Try to select a valid set of kinematics
482  unsigned int iter = 0;
483  bool accept=false;
484  double xsec=-1, g_E_l=-1, g_theta_l=-1, g_phi_l=-1, g_theta_pi=-1, g_phi_pi=-1;
485  double g_ctheta_l, g_ctheta_pi;
486 
487  while(1) {
488  iter++;
489  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
490 
491  //Select kinematic point
492  g_E_l = E_l_min + d_E_l * rnd->RndKine().Rndm();
493  g_ctheta_l = ctheta_l_min + d_ctheta_l * rnd->RndKine().Rndm();
494  g_ctheta_pi = ctheta_pi_min + d_ctheta_pi * rnd->RndKine().Rndm();
495  g_phi_l = phi_min + d_phi * rnd->RndKine().Rndm();
496  // random phi is relative to phi_l
497  g_phi_pi = g_phi_l + (phi_min + d_phi * rnd->RndKine().Rndm());
498  g_theta_l = TMath::ACos(g_ctheta_l);
499  g_theta_pi = TMath::ACos(g_ctheta_pi);
500 
501  LOG("COHKinematics", pINFO) << "Trying: Lep(" <<g_E_l << ", " <<
502  g_theta_l << ", " << g_phi_l << ") Pi(" <<
503  g_theta_pi << ", " << g_phi_pi << ")";
504 
505  this->SetKinematics(g_E_l, g_theta_l, g_phi_l, g_theta_pi, g_phi_pi,
506  interaction, interaction->KinePtr());
507 
508  // computing cross section for the current kinematics
509  xsec = fXSecModel->XSec(interaction,kPSElOlOpifE) / (1E-38 * units::cm2);
510 
511  if (!fGenerateUniformly) {
512  //-- decide whether to accept the current kinematics
513  double t = xsec_max * rnd->RndKine().Rndm();
514 
515  LOG("COHKinematics", pINFO) << "Got: xsec = " << xsec << ", t = " <<
516  t << " (max_xsec = " << xsec_max << ")";
517 
518  this->AssertXSecLimits(interaction, xsec, xsec_max);
519 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
520  LOG("COHKinematics", pDEBUG)
521  << "xsec= " << xsec << ", J= 1, Rnd= " << t;
522 #endif
523  accept = (t<xsec);
524  }
525  else {
526  accept = (xsec>0);
527  }
528 
529  //-- If the generated kinematics are accepted, finish-up module's job
530  if(accept) {
531  LOG("COHKinematics", pNOTICE) << "Selected: Lepton(" <<
532  g_E_l << ", " << g_theta_l << ", " <<
533  g_phi_l << ") Pion(" << g_theta_pi << ", " << g_phi_pi << ")";
534 
535  double E_l = g_E_l;
536  double theta_l = g_theta_l;
537  double theta_pi = g_theta_pi;
538  double phi_l = g_phi_l;
539  double phi_pi = g_phi_pi;
540  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
541  double E_nu = P4_nu.E();
542  double E_pi= E_nu-E_l;
543  double m_l = interaction->FSPrimLepton()->Mass();
544  double m_pi = this->pionMass(interaction);
545 
546  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
547  TVector3 lepton_3vector = TVector3(0,0,0);
548  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
549  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
550 
551  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
552  TVector3 pion_3vector = TVector3(0,0,0);
553  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
554  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
555 
556  TLorentzVector q = P4_nu - P4_lep;
557  double Q2 = -q.Mag2();
558  double x = Q2/(2*E_pi*constants::kNucleonMass);
559  double y = E_pi/E_nu;
560 
561  double t = TMath::Abs( (q - P4_pion).Mag2() );
562 
563  // for uniform kinematics, compute an event weight as
564  // wght = (phase space volume)*(differential xsec)/(event total xsec)
565  if(fGenerateUniformly) {
566  // Phase space volume needs checking
567  double vol = d_E_l*d_ctheta_l*d_phi*d_ctheta_pi*d_phi;
568  double totxsec = evrec->XSec();
569  double wght = (vol/totxsec)*xsec;
570  LOG("COHKinematics", pNOTICE) << "Kinematics wght = "<< wght;
571 
572  // apply computed weight to the current event weight
573  wght *= evrec->Weight();
574  LOG("COHKinematics", pNOTICE) << "Current event wght = " << wght;
575  evrec->SetWeight(wght);
576  }
577 
578  // reset bits
579  interaction->ResetBit(kISkipProcessChk);
580  interaction->ResetBit(kISkipKinematicChk);
581  // lock selected kinematics & clear running values
582  interaction->KinePtr()->Setx(x, true);
583  interaction->KinePtr()->Sety(y, true);
584  interaction->KinePtr()->Sett(t, true);
585  interaction->KinePtr()->SetW(kPionMass, true);
586  interaction->KinePtr()->SetQ2(2*kNucleonMass*x*y*E_nu, true);
587  interaction->KinePtr()->ClearRunningValues();
588  // set the cross section for the selected kinematics
589  evrec->SetDiffXSec(xsec,kPSElOlOpifE);
590  return;
591  }
592  }//while
593 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
static const double kNucleonMass
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
void SetKinematics(const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction, Kinematics *kinematics) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double cm2
Definition: Units.h:69
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:291
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
double pionMass(const Interaction *in) const
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void throwOnTooManyIterations(unsigned int iters, GHepRecord *evrec) const
#define pDEBUG
Definition: Messenger.h:63
void COHKinematicsGenerator::CalculateKin_BergerSehgal ( GHepRecord event_rec) const

Definition at line 86 of file COHKinematicsGenerator.cxx.

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

Referenced by ProcessEventRecord().

87 {
88  // Get the Primary Interacton object
89  Interaction * interaction = evrec->Summary();
90  interaction->SetBit(kISkipProcessChk);
91  interaction->SetBit(kISkipKinematicChk); // TODO: Why turn this off?
92 
93  // Initialise a random number generator
95 
96  //-- For the subsequent kinematic selection with the rejection method:
97  // Calculate the max differential cross section or retrieve it from the
98  // cache. Throw an exception and quit the evg thread if a non-positive
99  // value is found.
100  //
101  // TODO: We are not offering the "fGenerateUniformly" option here.
102  double xsec_max = this->MaxXSec(evrec);
103 
104  //-- Get the kinematical limits for the generated x,y
105  const KPhaseSpace & kps = interaction->PhaseSpace();
106  Range1D_t y = kps.YLim();
108  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
109 
110  const double ymin = y.min + kASmallNum;
111  const double ymax = y.max - kASmallNum;
112  const double dy = ymax - ymin;
113  const double Q2min = Q2.min + kASmallNum;
114  const double Q2max = Q2.max - kASmallNum;
115  const double dQ2 = Q2max - Q2min;
116 
117  unsigned int iter = 0;
118  bool accept=false;
119  double xsec=-1, gy=-1, gQ2=-1;
120 
121  while(1) {
122  iter++;
123  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
124 
125  //-- Select unweighted kinematics using importance sampling method.
126  // TODO: The importance sampling envelope is not used. Currently,
127  // we just employ a standard rejection-method approach.
128 
129  gy = ymin + dy * rnd->RndKine().Rndm();
130  gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
131 
132  LOG("COHKinematics", pINFO) <<
133  "Trying: Q^2 = " << gQ2 << ", y = " << gy; /* << ", t = " << gt; */
134 
135  interaction->KinePtr()->Sety(gy);
136  interaction->KinePtr()->SetQ2(gQ2);
137  kinematics::UpdateXFromQ2Y(interaction);
138 
139  // computing cross section for the current kinematics
140  xsec = fXSecModel->XSec(interaction, kPSQ2yfE);
141 
142  //-- decide whether to accept the current kinematics
143  accept = (xsec_max * rnd->RndKine().Rndm() < xsec);
144 
145  //-- If the generated kinematics are accepted, finish-up module's job
146  if(accept) {
147  LOG("COHKinematics", pNOTICE)
148  << "Selected: Q^2 = " << gQ2 << ", y = " << gy; /* << ", t = " << gt; */
149 
150  // the Berger-Sehgal COH cross section should be a triple differential cross section
151  // d^2xsec/dQ2dydt where t is the the square of the 4p transfer to the
152  // nucleus. The cross section used for kinematical selection should have
153  // the t-dependence integrated out. The t-dependence is of the form
154  // ~exp(-bt). Now that the x,y kinematical variables have been selected
155  // we can generate a t using the t-dependence as a PDF.
156  const InitialState & init_state = interaction->InitState();
157  double Ev = init_state.ProbeE(kRfLab);
158  double Epi = gy*Ev; // pion energy
159  double Epi2 = TMath::Power(Epi,2);
160  double pme2 = kPionMass2/Epi2;
161  double gx = interaction->KinePtr()->x(); // utils::kinematics::Q2YtoX(Ev,kNucleonMass,gQ2,gy); // gQ2 / ( 2. * gy * kNucleonMass * Ev);
162  double xME = kNucleonMass*gx/Epi;
163  double tA = 1. + xME - 0.5*pme2;
164  /* Range1D_t tl = kps.TLim(); // this becomes the bounds for t */
165  double tB = TMath::Sqrt(1.+ 2*xME) * TMath::Sqrt(1.-pme2);
166  double tmin = 2*Epi2 * (tA-tB);
167  double tmax = 2*Epi2 * (tA+tB);
168  double A = (double) init_state.Tgt().A();
169  double A13 = TMath::Power(A,1./3.);
170  double R = fRo * A13 * units::fermi; // nuclear radius
171  double R2 = TMath::Power(R,2.);
172  double b = 0.33333 * R2;
173  double tsum = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
174  double rt = tsum * rnd->RndKine().Rndm();
175  double gt = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/b;
176 
177  // TODO: If we re-install the fGenerateUniformly option, we
178  // would compute the event weight here.
179 
180  // reset bits
181  interaction->ResetBit(kISkipProcessChk);
182  interaction->ResetBit(kISkipKinematicChk);
183 
184  // lock selected kinematics & clear running values
185  interaction->KinePtr()->Setx(gx, true);
186  interaction->KinePtr()->Sety(gy, true);
187  interaction->KinePtr()->Sett(gt, true);
188  interaction->KinePtr()->SetW(this->pionMass(interaction), true);
189  interaction->KinePtr()->SetQ2(gQ2, true);
190  interaction->KinePtr()->ClearRunningValues();
191 
192  // set the cross section for the selected kinematics
193  evrec->SetDiffXSec(xsec * TMath::Exp(-b * gt) / tsum, kPSQ2yfE);
194 
195  return;
196  }
197  }// iterations
198 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
double fRo
nuclear scale parameter
static const double kNucleonMass
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
int A(void) const
Definition: Target.h:70
A simple [min,max] interval for doubles.
Definition: Range1.h:42
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
double x(bool selected=false) const
Definition: Kinematics.cxx:99
Range1D_t YLim(void) const
y limits
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
static constexpr double b
Definition: Units.h:78
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
double fQ2Min
lower bound of integration for Q^2 in Berger-Sehgal Model
Kinematical phase space.
Definition: KPhaseSpace.h:33
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:291
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
void UpdateXFromQ2Y(const Interaction *in)
Definition: KineUtils.cxx:1344
double pionMass(const Interaction *in) const
double gQ2
Definition: gtestDISSF.cxx:55
static constexpr double fermi
Definition: Units.h:55
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
double ProbeE(RefFrame_t rf) const
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void throwOnTooManyIterations(unsigned int iters, GHepRecord *evrec) const
Initial State information.
Definition: InitialState.h:48
void COHKinematicsGenerator::CalculateKin_BergerSehgalFM ( GHepRecord event_rec) const

Definition at line 200 of file COHKinematicsGenerator.cxx.

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

Referenced by ProcessEventRecord().

201 {
202  // Get the Primary Interacton object
203  Interaction * interaction = evrec->Summary();
204  interaction->SetBit(kISkipProcessChk);
205  interaction->SetBit(kISkipKinematicChk);
206 
207  // Initialise a random number generator
208  RandomGen * rnd = RandomGen::Instance();
209 
210  //-- For the subsequent kinematic selection with the rejection method:
211  // Calculate the max differential cross section or retrieve it from the
212  // cache. Throw an exception and quit the evg thread if a non-positive
213  // value is found.
214  //
215  // TODO: We are not offering the "fGenerateUniformly" option here.
216  double xsec_max = this->MaxXSec(evrec);
217 
218  //-- Get the kinematical limits for the generated x,y
219  const KPhaseSpace & kps = interaction->PhaseSpace();
220  Range1D_t y = kps.YLim();
222  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
223 
224  const double ymin = y.min + kASmallNum;
225  const double ymax = y.max - kASmallNum;
226  const double dy = ymax - ymin;
227  const double Q2min = Q2.min + kASmallNum;
228  const double Q2max = Q2.max - kASmallNum;
229  const double dQ2 = Q2max - Q2min;
230  const double tmin = kASmallNum;
231  const double tmax = fTMax - kASmallNum; // TODO: Choose realistic t bounds
232  const double dt = tmax - tmin;
233 
234  //-- Try to select a valid (Q^2,y,t) triple.
235 
236  unsigned int iter = 0;
237  bool accept=false;
238  double xsec=-1, gy=-1, gt=-1, gQ2=-1;
239 
240  while(1) {
241  iter++;
242  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
243 
244  //-- Select unweighted kinematics using importance sampling method.
245  // TODO: The importance sampling envelope is not used. Currently,
246  // we just employ a standard rejection-method approach.
247 
248  gy = ymin + dy * rnd->RndKine().Rndm();
249  gt = tmin + dt * rnd->RndKine().Rndm();
250  gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
251 
252  LOG("COHKinematics", pINFO) <<
253  "Trying: Q^2 = " << gQ2 << ", y = " << gy << ", t = " << gt;
254 
255  interaction->KinePtr()->Sety(gy);
256  interaction->KinePtr()->Sett(gt);
257  interaction->KinePtr()->SetQ2(gQ2);
258 
259  // computing cross section for the current kinematics
260  xsec = fXSecModel->XSec(interaction, kPSxyfE);
261 
262  //-- decide whether to accept the current kinematics
263  accept = (xsec_max * rnd->RndKine().Rndm() < xsec);
264 
265  //-- If the generated kinematics are accepted, finish-up module's job
266  if(accept) {
267  LOG("COHKinematics", pNOTICE)
268  << "Selected: Q^2 = " << gQ2 << ", y = " << gy << ", t = " << gt;
269 
270  // TODO: If we re-install the fGenerateUniformly option, we
271  // would compute the event weight here.
272 
273  // reset bits
274  interaction->ResetBit(kISkipProcessChk);
275  interaction->ResetBit(kISkipKinematicChk);
276 
277  // lock selected kinematics & clear running values
278  interaction->KinePtr()->SetQ2(gQ2, true);
279  interaction->KinePtr()->Sety(gy, true);
280  interaction->KinePtr()->Sett(gt, true);
281  interaction->KinePtr()->SetW(this->pionMass(interaction), true);
282  interaction->KinePtr()->ClearRunningValues();
283 
284  // set the cross section for the selected kinematics
285  evrec->SetDiffXSec(xsec, kPSxytfE);
286 
287  return;
288  }
289  }// iterations
290 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
A simple [min,max] interval for doubles.
Definition: Range1.h:42
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
Range1D_t YLim(void) const
y limits
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Summary information for an interaction.
Definition: Interaction.h:56
double fTMax
upper bound for t = (q - p_pi)^2
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fQ2Min
lower bound of integration for Q^2 in Berger-Sehgal Model
Kinematical phase space.
Definition: KPhaseSpace.h:33
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:291
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
double pionMass(const Interaction *in) const
double gQ2
Definition: gtestDISSF.cxx:55
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void throwOnTooManyIterations(unsigned int iters, GHepRecord *evrec) const
void COHKinematicsGenerator::CalculateKin_ReinSehgal ( GHepRecord event_rec) const

Definition at line 292 of file COHKinematicsGenerator.cxx.

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

Referenced by ProcessEventRecord().

293 {
294  // Get the Primary Interacton object
295  Interaction * interaction = evrec->Summary();
296  interaction->SetBit(kISkipProcessChk);
297  interaction->SetBit(kISkipKinematicChk);
298 
299  //-- Get the random number generators
300  RandomGen * rnd = RandomGen::Instance();
301 
302  //-- For the subsequent kinematic selection with the rejection method:
303  // Calculate the max differential cross section or retrieve it from the
304  // cache. Throw an exception and quit the evg thread if a non-positive
305  // value is found.
306  // If the kinematics are generated uniformly over the allowed phase
307  // space the max xsec is irrelevant
308  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
309 
310  //-- Get the kinematical limits for the generated x,y
311  const KPhaseSpace & kps = interaction->PhaseSpace();
312  Range1D_t y = kps.YLim();
313  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
314 
315  const double xmin = kASmallNum;
316  const double xmax = 1.- kASmallNum;
317  const double ymin = y.min + kASmallNum;
318  const double ymax = y.max - kASmallNum;
319  const double dx = xmax - xmin;
320  const double dy = ymax - ymin;
321 
322  //------ Try to select a valid x,y pair
323 
324  unsigned int iter = 0;
325  bool accept=false;
326  double xsec=-1, gx=-1, gy=-1;
327 
328  while(1) {
329  iter++;
330  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
331 
332  if(fGenerateUniformly) {
333  //-- Generate a x,y pair uniformly in the kinematically allowed range.
334  gx = xmin + dx * rnd->RndKine().Rndm();
335  gy = ymin + dy * rnd->RndKine().Rndm();
336 
337  } else {
338  //-- Select unweighted kinematics using importance sampling method.
339 
340  if(iter==1) {
341  LOG("COHKinematics", pNOTICE) << "Initializing the sampling envelope";
342  double Ev = interaction->InitState().ProbeE(kRfLab);
343  fEnvelope->SetRange(xmin,ymin,xmax,ymax);
344  fEnvelope->SetParameter(0, xsec_max);
345  fEnvelope->SetParameter(1, Ev);
346  }
347 
348  // Generate W,QD2 using the 2-D envelope as PDF
349  fEnvelope->GetRandom2(gx,gy);
350  }
351 
352  LOG("COHKinematics", pINFO) << "Trying: x = " << gx << ", y = " << gy;
353 
354  interaction->KinePtr()->Setx(gx);
355  interaction->KinePtr()->Sety(gy);
356 
357  // computing cross section for the current kinematics
358  xsec = fXSecModel->XSec(interaction, kPSxyfE);
359 
360  //-- decide whether to accept the current kinematics
361  if(!fGenerateUniformly) {
362  double max = fEnvelope->Eval(gx, gy);
363  double t = max * rnd->RndKine().Rndm();
364 
365  this->AssertXSecLimits(interaction, xsec, max);
366 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
367  LOG("COHKinematics", pDEBUG)
368  << "xsec= " << xsec << ", J= 1, Rnd= " << t;
369 #endif
370  accept = (t<xsec);
371  }
372  else {
373  accept = (xsec>0);
374  }
375 
376  //-- If the generated kinematics are accepted, finish-up module's job
377  if(accept) {
378  LOG("COHKinematics", pNOTICE) << "Selected: x = "<< gx << ", y = "<< gy;
379 
380  // the Rein-Sehgal COH cross section should be a triple differential cross section
381  // d^2xsec/dxdydt where t is the the square of the 4p transfer to the
382  // nucleus. The cross section used for kinematical selection should have
383  // the t-dependence integrated out. The t-dependence is of the form
384  // ~exp(-bt). Now that the x,y kinematical variables have been selected
385  // we can generate a t using the t-dependence as a PDF.
386  const InitialState & init_state = interaction->InitState();
387  double Ev = init_state.ProbeE(kRfLab);
388  double Epi = gy*Ev; // pion energy
389  double Epi2 = TMath::Power(Epi,2);
390  double pme2 = kPionMass2/Epi2;
391  double xME = kNucleonMass*gx/Epi;
392  double tA = 1. + xME - 0.5*pme2;
393  double tB = TMath::Sqrt(1.+ 2*xME) * TMath::Sqrt(1.-pme2);
394  double tmin = 2*Epi2 * (tA-tB);
395  double tmax = 2*Epi2 * (tA+tB);
396  double A = (double) init_state.Tgt().A();
397  double A13 = TMath::Power(A,1./3.);
398  double R = fRo * A13 * units::fermi; // nuclear radius
399  double R2 = TMath::Power(R,2.);
400  double b = 0.33333 * R2;
401  double tsum = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
402  double rt = tsum * rnd->RndKine().Rndm();
403  double gt = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/b;
404 
405  LOG("COHKinematics", pNOTICE)
406  << "Selected: t = "<< gt << ", from ["<< tmin << ", "<< tmax << "]";
407 
408  // for uniform kinematics, compute an event weight as
409  // wght = (phase space volume)*(differential xsec)/(event total xsec)
410  if(fGenerateUniformly) {
411  double vol = y.max-y.min; // dx=1, dt: irrelevant
412  double totxsec = evrec->XSec();
413  double wght = (vol/totxsec)*xsec;
414  LOG("COHKinematics", pNOTICE) << "Kinematics wght = "<< wght;
415 
416  // apply computed weight to the current event weight
417  wght *= evrec->Weight();
418  LOG("COHKinematics", pNOTICE) << "Current event wght = " << wght;
419  evrec->SetWeight(wght);
420  }
421 
422  // reset bits
423  interaction->ResetBit(kISkipProcessChk);
424  interaction->ResetBit(kISkipKinematicChk);
425 
426  // lock selected kinematics & clear running values
427  interaction->KinePtr()->Setx(gx, true);
428  interaction->KinePtr()->Sety(gy, true);
429  interaction->KinePtr()->Sett(gt, true);
430  interaction->KinePtr()->SetW(kPionMass, true);
431  interaction->KinePtr()->SetQ2(2*kNucleonMass*gx*gy*Ev, true);
432  interaction->KinePtr()->ClearRunningValues();
433 
434  // set the cross section for the selected kinematics
435  evrec->SetDiffXSec(xsec * TMath::Exp(-b * gt) / tsum, kPSxytfE);
436 
437  return;
438  }
439  }// iterations
440 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double fRo
nuclear scale parameter
static const double kNucleonMass
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
int A(void) const
Definition: Target.h:70
A simple [min,max] interval for doubles.
Definition: Range1.h:42
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
Range1D_t YLim(void) const
y limits
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
static constexpr double b
Definition: Units.h:78
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
Kinematical phase space.
Definition: KPhaseSpace.h:33
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:291
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
static constexpr double fermi
Definition: Units.h:55
TF2 * fEnvelope
2-D envelope used for importance sampling
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
double ProbeE(RefFrame_t rf) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void throwOnTooManyIterations(unsigned int iters, GHepRecord *evrec) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
bool COHKinematicsGenerator::CheckKinematics ( const double  E_l,
const double  theta_l,
const double  phi_l,
const double  theta_pi,
const double  phi_pi,
const Interaction interaction 
) const

Definition at line 641 of file COHKinematicsGenerator.cxx.

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

647 {
648  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
649  double E_nu = P4_nu.E();
650  double E_pi= E_nu-E_l;
651  double m_l = interaction->FSPrimLepton()->Mass();
652  double m_pi;
653  if ( interaction->ProcInfo().IsWeakCC() ) {
654  m_pi = constants::kPionMass;
655  }
656  else {
657  m_pi = constants::kPi0Mass;
658  }
659  if (E_l <= m_l) {
660  return false;
661  }
662  if (E_pi <= m_pi) {
663  return false;
664  }
665  return true;
666 }
bool IsWeakCC(void) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
double COHKinematicsGenerator::ComputeMaxXSec ( const Interaction in) const
virtual

Implements genie::KineGeneratorWithCache.

Definition at line 668 of file COHKinematicsGenerator.cxx.

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

669 {
670  // Computes the maximum differential cross section in the requested phase
671  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
672  // method and the value is cached at a circular cache branch for retrieval
673  // during subsequent event generation.
674 
675 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
676  SLOG("COHKinematics", pDEBUG)
677  << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
678 #endif
679  double max_xsec = 0.;
680  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
681  max_xsec = MaxXSec_ReinSehgal(in);
682  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015")) {
683  max_xsec = MaxXSec_BergerSehgal(in);
684  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015")) {
685  max_xsec = MaxXSec_BergerSehgalFM(in);
686  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
687  max_xsec = MaxXSec_AlvarezRuso(in);
688  }
689  else {
690  LOG("COHKinematicsGenerator",pFATAL) <<
691  "ComputeMaxXSec >> Cannot calculate max cross-section for " <<
692  fXSecModel->Id().Name();
693  }
694 
695  // Apply safety factor, since value retrieved from the cache might
696  // correspond to a slightly different energy.
697  max_xsec *= fSafetyFactor;
698 
699 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
700  SLOG("COHKinematics", pDEBUG) << in->AsString();
701  SLOG("COHKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
702  SLOG("COHKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
703 #endif
704 
705  return max_xsec;
706 }
double MaxXSec_ReinSehgal(const Interaction *in) const
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
#define pFATAL
Definition: Messenger.h:56
double MaxXSec_AlvarezRuso(const Interaction *in) const
string AsString(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string Name(void) const
Definition: AlgId.h:44
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
double MaxXSec_BergerSehgal(const Interaction *in) 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
double MaxXSec_BergerSehgalFM(const Interaction *in) const
#define pDEBUG
Definition: Messenger.h:63
void COHKinematicsGenerator::Configure ( const Registry config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 955 of file COHKinematicsGenerator.cxx.

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

956 {
957  Algorithm::Configure(config);
958  this->LoadConfig();
959 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void COHKinematicsGenerator::Configure ( string  config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 961 of file COHKinematicsGenerator.cxx.

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

962 {
963  Algorithm::Configure(config);
964  this->LoadConfig();
965 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
double COHKinematicsGenerator::Energy ( const Interaction in) const
virtual

Reimplemented from genie::KineGeneratorWithCache.

Definition at line 921 of file COHKinematicsGenerator.cxx.

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

922 {
923  // Override the base class Energy() method to cache the max xsec for the
924  // neutrino energy in the LAB rather than in the hit nucleon rest frame.
925 
926  const InitialState & init_state = interaction->InitState();
927  double E = init_state.ProbeE(kRfLab);
928  return E;
929 }
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
void COHKinematicsGenerator::LoadConfig ( void  )

Definition at line 967 of file COHKinematicsGenerator.cxx.

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

Referenced by Configure().

968 {
969  //-- COH model parameter Ro
970  GetParam( "COH-Ro", fRo );
971  //-- COH model parameter t_max for t = (q - p_pi)^2
972  GetParam( "COH-t-max", fTMax ) ;
973  //-- COH model bounds of integration for Q^2
974  GetParam( "COH-Q2-min", fQ2Min ) ;
975  GetParam( "COH-Q2-max", fQ2Max ) ;
976 
977  //-- max xsec safety factor (for rejection method) and min cached energy
978  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.6 ) ;
979  GetParamDef( "Cache-MinEnergy", fEMin, -1.0 ) ;
980 
981  //-- Generate kinematics uniformly over allowed phase space and compute
982  // an event weight?
983  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
984 
985  //-- Maximum allowed fractional cross section deviation from maxim cross
986  // section used in rejection method
987  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
988  assert(fMaxXSecDiffTolerance>=0);
989 
990  //-- Envelope employed when importance sampling is used
991  // (initialize with dummy range)
992  if(fEnvelope) delete fEnvelope;
993  fEnvelope = new TF2("CohKinEnvelope",
995  // stop ROOT from deleting this object of its own volition
996  gROOT->GetListOfFunctions()->Remove(fEnvelope);
997 }
double COHImportanceSamplingEnvelope(double *x, double *par)
Definition: KineUtils.cxx:1466
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double fRo
nuclear scale parameter
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec&gt;maxxsec
double fTMax
upper bound for t = (q - p_pi)^2
double fQ2Min
lower bound of integration for Q^2 in Berger-Sehgal Model
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
TF2 * fEnvelope
2-D envelope used for importance sampling
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
double COHKinematicsGenerator::MaxXSec_AlvarezRuso ( const Interaction in) const

Definition at line 851 of file COHKinematicsGenerator.cxx.

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

Referenced by ComputeMaxXSec().

852 {
853  // Computes the maximum differential cross section in the requested phase
854  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
855  // method and the value is cached at a circular cache branch for retrieval
856  // during subsequent event generation.
857 
858 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
859  SLOG("COHKinematics", pDEBUG)
860  << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
861 #endif
862  double max_xsec = 0.;
863  double Ev = in->InitState().ProbeE(kRfLab);
864 
865  const KPhaseSpace & kps = in->PhaseSpace();
866  Range1D_t y = kps.YLim();
867 
868  ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit2");
870  f.SetFactor(-1.); // Make it return negative of cross-section so we can minimize
871 
872  min->SetFunction( f );
873  min->SetMaxFunctionCalls(10000);
874  min->SetTolerance(0.05);
875 
876  const double min_el = in->FSPrimLepton()->Mass();
877  const double max_el = Ev - kPionMass;
878  const unsigned int n_el = 100;
879  const double d_el = (max_el - min_el) / double(n_el - 1);
880 
881  const double min_thetal = kASmallNum;
882  const double max_thetal = kPi / 4.0;
883  const unsigned int n_thetal = 10;
884  const double d_thetal = (max_thetal - min_thetal) / double(n_thetal - 1);
885 
886  const double min_thetapi = kASmallNum;
887  const double max_thetapi = kPi / 2.0;
888  const unsigned int n_thetapi = 10;
889  const double d_thetapi = (max_thetapi - min_thetapi) / double(n_thetapi - 1);
890 
891  //~ const double min_phipi = kPi;
892  //~ const double max_phipi = 0.5 * kPi;
893  const double min_phipi = kASmallNum;
894  const double max_phipi = 2*kPi-kASmallNum;
895  const unsigned int n_phipi = 10;
896  const double d_phipi = (max_phipi - min_phipi) / double(n_phipi - 1);
897 
898  min->SetLimitedVariable ( 0 ,"E_lep" , max_el -kASmallNum , d_el , min_el , max_el );
899  min->SetLimitedVariable ( 1 ,"theta_l" , min_thetal +kASmallNum , d_thetal , min_thetal , max_thetal );
900  min->SetLimitedVariable ( 2 ,"theta_pi" , min_thetapi+kASmallNum , d_thetapi , min_thetapi, max_thetapi );
901  min->SetLimitedVariable ( 3 ,"phi_pi" , min_phipi +kASmallNum , d_phipi , min_phipi , max_phipi );
902 
903  min->Minimize();
904  max_xsec = -min->MinValue(); //back to positive xsec
905 
906  // Apply safety factor, since value retrieved from the cache might
907  // correspond to a slightly different energy.
908  max_xsec *= fSafetyFactor;
909 
910 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
911  SLOG("COHKinematics", pDEBUG) << in->AsString();
912  SLOG("COHKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
913  SLOG("COHKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
914 #endif
915 
916  delete min;
917 
918  return max_xsec;
919 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t YLim(void) const
y limits
string AsString(void) const
Kinematical phase space.
Definition: KPhaseSpace.h:33
static const double kASmallNum
Definition: Controls.h:40
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
const InitialState & InitState(void) const
Definition: Interaction.h:69
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
double ProbeE(RefFrame_t rf) const
#define pDEBUG
Definition: Messenger.h:63
double COHKinematicsGenerator::MaxXSec_BergerSehgal ( const Interaction in) const

Definition at line 708 of file COHKinematicsGenerator.cxx.

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

Referenced by ComputeMaxXSec().

709 {
710  double max_xsec = 0;
711  const int NQ2 = 50;
712  const int Ny = 50;
713 
714  const KPhaseSpace & kps = in->PhaseSpace();
715  Range1D_t Q2r = kps.Q2Lim();
716  Q2r.max = fQ2Max;
717 
718  const double logQ2min = TMath::Log10(Q2r.min + kASmallNum);
719  const double logQ2max = TMath::Log10(Q2r.max);
720  const double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
721 
722  for(int i=0; i<NQ2; i++) {
723  double Q2 = TMath::Power(10, logQ2min + i * dlogQ2);
724  in->KinePtr()->SetQ2(Q2);
725 
726  Range1D_t yr = kps.YLim();
727  if ((yr.max < 0) || (yr.max < yr.min) ||
728  (yr.max > 1) || (yr.min < 0)) { // forbidden kinematics
729  continue;
730  }
731  const double logymin = TMath::Log10(yr.min);
732  const double logymax = TMath::Log10(yr.max);
733  const double dlogy = (logymax - logymin) /(Ny-1);
734 
735  for(int j=0; j<Ny; j++) {
736  double gy = TMath::Power(10, logymin + j * dlogy);
737  in->KinePtr()->Sety(gy);
738 
739  /* Range1D_t tl = kps.TLim(); // TESTING! - this becomes a loop over t */
741 
742  // Note: We're not stepping through log Q^2, log y - we "unpacked"
743  double xsec = fXSecModel->XSec(in, kPSQ2yfE);
744 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
745  LOG("COHKinematics", pDEBUG)
746  << "xsec(Q2= " << Q2 << ", y= " << gy << ", t = " << gt << ") = " << xsec;
747 #endif
748  max_xsec = TMath::Max(max_xsec, xsec);
749 
750  } // y
751  } // Q2
752  return max_xsec;
753 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
A simple [min,max] interval for doubles.
Definition: Range1.h:42
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
Kinematical phase space.
Definition: KPhaseSpace.h:33
static const double kASmallNum
Definition: Controls.h:40
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
double max
Definition: Range1.h:53
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
void UpdateXFromQ2Y(const Interaction *in)
Definition: KineUtils.cxx:1344
double min
Definition: Range1.h:52
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
#define pDEBUG
Definition: Messenger.h:63
double COHKinematicsGenerator::MaxXSec_BergerSehgalFM ( const Interaction in) const

Definition at line 755 of file COHKinematicsGenerator.cxx.

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

Referenced by ComputeMaxXSec().

756 {
757  double max_xsec = 0;
758  // How many sampling bins in each variable for max xsec calculation?
759  const int NQ2 = 50;
760  const int Ny = 50;
761  const int Nt = 50;
762 
763  const KPhaseSpace & kps = in->PhaseSpace();
764  Range1D_t Q2r = kps.Q2Lim();
765  Q2r.max = fQ2Max;
766 
767  const double logQ2min = TMath::Log10(Q2r.min + kASmallNum);
768  const double logQ2max = TMath::Log10(Q2r.max);
769  const double logtmin = TMath::Log10(kASmallNum);
770  const double logtmax = TMath::Log10(fTMax - kASmallNum);
771  const double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
772  const double dlogt = (logtmax - logtmin) /(Nt-1);
773 
774  for(int i=0; i<NQ2; i++) {
775  double Q2 = TMath::Power(10, logQ2min + i * dlogQ2);
776  in->KinePtr()->SetQ2(Q2);
777 
778  Range1D_t yr = kps.YLim();
779  if ((yr.max < 0) || (yr.max < yr.min) ||
780  (yr.max > 1) || (yr.min < 0)) { // forbidden kinematics
781  continue;
782  }
783  const double logymin = TMath::Log10(yr.min);
784  const double logymax = TMath::Log10(yr.max);
785  const double dlogy = (logymax - logymin) /(Ny-1);
786 
787  for(int j=0; j<Ny; j++) {
788  double gy = TMath::Power(10, logymin + j * dlogy);
789 
790  for(int k=0; k<Nt; k++) {
791  double gt = TMath::Power(10, logtmin + k * dlogt);
792 
793  in->KinePtr()->Sety(gy);
794  in->KinePtr()->Sett(gt);
795 
796  double xsec = fXSecModel->XSec(in, kPSxyfE);
797 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
798  LOG("COHKinematics", pDEBUG)
799  << "xsec(Q2= " << Q2 << ", y= " << gy << ", t = " << gt << ") = " << xsec;
800 #endif
801  max_xsec = TMath::Max(max_xsec, xsec);
802 
803  } // t
804  } // y
805  } // Q2
806  return max_xsec;
807 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t YLim(void) const
y limits
Range1D_t Q2Lim(void) const
Q2 limits.
double fTMax
upper bound for t = (q - p_pi)^2
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
Kinematical phase space.
Definition: KPhaseSpace.h:33
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:291
static const double kASmallNum
Definition: Controls.h:40
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
double max
Definition: Range1.h:53
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
double min
Definition: Range1.h:52
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
#define pDEBUG
Definition: Messenger.h:63
double COHKinematicsGenerator::MaxXSec_ReinSehgal ( const Interaction in) const

Definition at line 809 of file COHKinematicsGenerator.cxx.

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

Referenced by ComputeMaxXSec().

810 {
811  double max_xsec = 0;
812  double Ev = in->InitState().ProbeE(kRfLab);
813 
814  const int Nx = 50;
815  const int Ny = 50;
816 
817  const KPhaseSpace & kps = in->PhaseSpace();
818  Range1D_t y = kps.YLim();
819 
820  const double logxmin = TMath::Log10(1E-5);
821  const double logxmax = TMath::Log10(1.0);
822  const double logymin = TMath::Log10(y.min);
823  const double logymax = TMath::Log10(y.max);
824 
825  const double dlogx = (logxmax - logxmin) /(Nx-1);
826  const double dlogy = (logymax - logymin) /(Ny-1);
827 
828  for(int i=0; i<Nx; i++) {
829  double gx = TMath::Power(10, logxmin + i * dlogx);
830  for(int j=0; j<Ny; j++) {
831  double gy = TMath::Power(10, logymin + j * dlogy);
832 
833  double Q2 = 2*kNucleonMass*gx*gy*Ev;
834  if(Ev>1.0 && Q2>0.01) continue;
835 
836  in->KinePtr()->Setx(gx);
837  in->KinePtr()->Sety(gy);
838 
839  double xsec = fXSecModel->XSec(in, kPSxyfE);
840 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
841  LOG("COHKinematics", pDEBUG)
842  << "xsec(x= " << gx << ", y= " << gy << ") = " << xsec;
843 #endif
844  max_xsec = TMath::Max(max_xsec, xsec);
845 
846  }//y
847  }//x
848  return max_xsec;
849 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
static const double kNucleonMass
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t YLim(void) const
y limits
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
Kinematical phase space.
Definition: KPhaseSpace.h:33
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
double max
Definition: Range1.h:53
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
double ProbeE(RefFrame_t rf) const
#define pDEBUG
Definition: Messenger.h:63
double COHKinematicsGenerator::pionMass ( const Interaction in) const
private

Definition at line 931 of file COHKinematicsGenerator.cxx.

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

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

932 {
933  double m_pi = 0.0;
934  if ( in->ProcInfo().IsWeakCC() ) {
935  m_pi = constants::kPionMass;
936  } else {
937  m_pi = constants::kPi0Mass;
938  }
939  return m_pi;
940 }
bool IsWeakCC(void) const
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void COHKinematicsGenerator::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 59 of file COHKinematicsGenerator.cxx.

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

60 {
61  if(fGenerateUniformly) {
62  LOG("COHKinematics", pNOTICE)
63  << "Generating kinematics uniformly over the allowed phase space";
64  }
65 
66  //-- Access cross section algorithm for running thread
68  const EventGeneratorI * evg = rtinfo->RunningThread();
69  fXSecModel = evg->CrossSectionAlg();
70  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
72  } else if (fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015") {
74  } else if (fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015") {
76  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
78  }
79  else {
80  LOG("COHKinematicsGenerator",pFATAL) <<
81  "ProcessEventRecord >> Cannot calculate kinematics for " <<
82  fXSecModel->Id().Name();
83  }
84 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
void CalculateKin_BergerSehgalFM(GHepRecord *event_rec) const
#define pFATAL
Definition: Messenger.h:56
Defines the EventGeneratorI interface.
void CalculateKin_AlvarezRuso(GHepRecord *event_rec) const
void CalculateKin_ReinSehgal(GHepRecord *event_rec) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string Name(void) const
Definition: AlgId.h:44
void CalculateKin_BergerSehgal(GHepRecord *event_rec) const
static RunningThreadInfo * Instance(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
#define pNOTICE
Definition: Messenger.h:61
const EventGeneratorI * RunningThread(void)
Keep info on the event generation thread currently on charge. This is used so that event generation m...
void COHKinematicsGenerator::SetKinematics ( const double  E_l,
const double  theta_l,
const double  phi_l,
const double  theta_pi,
const double  phi_pi,
const Interaction interaction,
Kinematics kinematics 
) const

Definition at line 595 of file COHKinematicsGenerator.cxx.

References genie::Interaction::FSPrimLepton(), genie::InitialState::GetProbeP4(), genie::Interaction::InitStatePtr(), genie::ProcessInfo::IsWeakCC(), genie::constants::kNucleonMass, genie::constants::kPi0Mass, genie::constants::kPionMass, genie::kRfLab, genie::Interaction::ProcInfo(), genie::utils::kinematics::Q2(), genie::Kinematics::SetFSLeptonP4(), genie::Kinematics::SetHadSystP4(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), and genie::utils::kinematics::UpdateWQ2FromXY().

Referenced by CalculateKin_AlvarezRuso().

602 {
603  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
604  double E_nu = P4_nu.E();
605  double E_pi= E_nu-E_l;
606  double m_l = interaction->FSPrimLepton()->Mass();
607  double m_pi;
608  if ( interaction->ProcInfo().IsWeakCC() ) {
609  m_pi = constants::kPionMass;
610  } else {
611  m_pi = constants::kPi0Mass;
612  }
613  double p_l=0.0;
614  if (E_l > m_l) {
615  p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
616  }
617  TVector3 lepton_3vector = TVector3(0,0,0);
618  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
619  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
620 
621  double p_pi=0.0;
622  if (E_pi > m_pi) {
623  p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
624  }
625  TVector3 pion_3vector = TVector3(0,0,0);
626  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
627  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
628 
629  double Q2 = -(P4_nu-P4_lep).Mag2();
630  double x = Q2/(2*E_pi*constants::kNucleonMass);
631  double y = E_pi/E_nu;
632 
633  kinematics->Setx(x);
634  kinematics->Sety(y);
635  kinematics::UpdateWQ2FromXY(interaction);
636 
637  kinematics->SetFSLeptonP4(P4_lep );
638  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
639 }
bool IsWeakCC(void) const
static const double kNucleonMass
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1290
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
void COHKinematicsGenerator::throwOnTooManyIterations ( unsigned int  iters,
GHepRecord evrec 
) const
private

Definition at line 942 of file COHKinematicsGenerator.cxx.

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

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

944 {
945  LOG("COHKinematics", pWARN)
946  << "*** Could not select valid kinematics after "
947  << iters << " iterations";
948  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
950  exception.SetReason("Couldn't select kinematics");
951  exception.SwitchOnFastForward();
952  throw exception;
953 }
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pWARN
Definition: Messenger.h:60
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:117

Member Data Documentation

TF2* genie::COHKinematicsGenerator::fEnvelope
mutable

2-D envelope used for importance sampling

Definition at line 71 of file COHKinematicsGenerator.h.

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

double genie::COHKinematicsGenerator::fQ2Max
private

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

Definition at line 79 of file COHKinematicsGenerator.h.

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

double genie::COHKinematicsGenerator::fQ2Min
private

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

Definition at line 78 of file COHKinematicsGenerator.h.

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

double genie::COHKinematicsGenerator::fRo

nuclear scale parameter

Definition at line 72 of file COHKinematicsGenerator.h.

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

double genie::COHKinematicsGenerator::fTMax
private

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

Definition at line 80 of file COHKinematicsGenerator.h.

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


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