40 using namespace genie;
67 <<
"Generating kinematics uniformly over the allowed phase space";
98 TVector3
beta = pnuc4.BoostVector();
101 double enu = P4_nu.E();
107 const double Tkmax = enu - mk - ml;
108 const double Tlmax = enu - mk - ml;
112 LOG(
"SKKinematics",
pWARN) <<
"No available phase space";
115 exception.
SetReason(
"No available phase space");
120 const double Tkmin = 0.0;
121 const double Tlmin = 0.0;
124 const double xmax = 0.69314718056;
125 const double phikqmin = 0.0;
126 const double phikqmax = 2.0 *
kPi;
127 const double dtk = Tkmax - Tkmin;
128 const double dtl = Tlmax - Tlmin;
129 const double dx = xmax -
xmin;
130 const double dphikq = phikqmax - phikqmin;
134 unsigned int iter = 0;
139 double costhetal = -1;
146 <<
"*** Could not select a valid (tk, tl, costhetal) triplet after " 147 << iter <<
" iterations";
150 exception.
SetReason(
"Couldn't select kinematics");
156 tk = Tkmin + dtk * rnd->
RndKine().Rndm();
157 tl = Tlmin + dtl * rnd->
RndKine().Rndm();
158 double x = xmin + dx * rnd->
RndKine().Rndm();
159 costhetal = 1.0 - TMath::Exp(x);
160 phikq = phikqmin + dphikq * rnd->
RndKine().Rndm();
162 tk = Tkmin + dtk * rnd->
RndKine().Rndm();
163 tl = Tlmin + dtl * rnd->
RndKine().Rndm();
164 double x = xmin + dx * rnd->
RndKine().Rndm();
165 costhetal = 1.0 - TMath::Exp(x);
166 phikq = phikqmin + dphikq * rnd->
RndKine().Rndm();
169 LOG(
"SKKinematics",
pDEBUG) <<
"Trying: Tk = " << tk <<
", Tl = " << tl <<
", cosThetal = " << costhetal <<
", phikq = " << phikq;
181 TVector3 lepton_3vector = TVector3(0,0,0);
182 lepton_3vector.SetMagThetaPhi(pl,TMath::ACos(costhetal),0.0);
183 TLorentzVector P4_lep( lepton_3vector, tl+ml );
184 TLorentzVector q = P4_nu - P4_lep;
185 double Q2 = -q.Mag2();
186 double xbj = Q2/(2*M*q.E());
187 double y = q.E()/P4_nu.E();
188 double W2 = (pnuc4+q).
Mag2();
197 double max = xsec_max;
198 double t = max * rnd->
RndKine().Rndm();
199 double J = TMath::Abs(1. - costhetal);
203 if( xsec*J > xsec_max ) {
205 <<
"!!!!!!XSEC ABOVE MAX!!!!! xsec= " << xsec <<
", J= " << J <<
", xsec*J = " << xsec*J <<
" max= " << xsec_max;
208 accept = (t< J*
xsec);
224 LOG(
"SKKinematics",
pNOTICE) <<
"Current event wght = " << wght;
227 LOG(
"SKKinematics",
pWARN) <<
"\nLepton energy (rest frame) = " << el <<
" kaon = " << tl + mk;
258 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 260 <<
"Scanning the allowed phase space {K} for the max(dxsec/d{K})";
270 const int Nctl = 100;
279 const double Tkmax = enu - mk - ml;
280 const double Tlmax = enu - mk - ml;
281 const double Tkmin = 0.0;
282 const double Tlmin = 0.0;
286 const double xmax = 0.69314718056;
287 const double dtk = (Tkmax - Tkmin)/Ntk;
288 const double dtl = (Tlmax - Tlmin)/Ntl;
289 const double dx = (xmax -
xmin)/Nctl;
292 for(
int i = 1;
i < Ntk;
i++) {
293 double tk = Tkmin + dtk*
i;
294 for(
int j = 1;
j < Ntl;
j++) {
295 double tl = Tlmin + dtl*
j;
296 for(
int k = 0; k < Nctl; k++) {
297 double log_1_minus_cosine_theta_lepton = xmin + dx*k;
299 double ctl = 1.0 - TMath::Exp(log_1_minus_cosine_theta_lepton);
317 if( xsec > max_xsec ) {
327 LOG(
"SKKinematics",
pINFO) <<
"Max XSec is " << max_xsec <<
" for enu = " << enu <<
" tk = " << max_tk <<
" tl = " << max_tl <<
" cosine theta = " << max_ctl;
333 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 335 SLOG(
"SKKinematics",
pDEBUG) <<
"Max xsec in phase space = " << max_xsec;
virtual double MaxXSec(GHepRecord *evrec) const
virtual void SetWeight(double wght)
std::map< std::string, double > xmax
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
THE MAIN GENIE PROJECT NAMESPACE
void Configure(const Registry &config)
double Q2(const Interaction *const i)
virtual Interaction * Summary(void) const
static RandomGen * Instance()
Access instance.
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
Defines the EventGeneratorI interface.
virtual double Weight(void) const
::xsd::cxx::tree::exception< char > exception
Abstract class. Provides a data caching mechanism for for concrete implementations of the EventRecord...
A singleton holding random number generator classes. All random number generation in GENIE should tak...
string AsString(void) const
const XSecAlgorithmI * fXSecModel
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
Summary information for an interaction.
const TLorentzVector & HitNucP4(void) const
double ComputeMaxXSec(const Interaction *in) const
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
int StrangeHadronPdg(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
virtual void Configure(const Registry &config)
void ProcessEventRecord(GHepRecord *event_rec) const
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Misc GENIE control constants.
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
void Setx(double x, bool selected=false)
void SetKV(KineVar_t kv, double value)
static RunningThreadInfo * Instance(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
void SetW(double W, bool selected=false)
void SwitchOnFastForward(void)
static PDGLibrary * Instance(void)
void SetReason(string reason)
virtual TBits * EventFlags(void) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
A registry. Provides the container for algorithm configuration parameters.
void Sety(double y, bool selected=false)
const XclsTag & ExclTag(void) const
TParticlePDG * Find(int pdgc)
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
assert(nhit_max >=nhit_nbins)
InitialState * InitStatePtr(void) const
const InitialState & InitState(void) const
double fMinLog1MinusCosTheta
double Energy(const Interaction *in) const
void ClearRunningValues(void)
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
static const unsigned int kRjMaxIterations
void CalculateKin_AtharSingleKaon(GHepRecord *event_rec) const
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const EventGeneratorI * RunningThread(void)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
double ProbeE(RefFrame_t rf) const
GENIE's GHEP MC event record.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Keep info on the event generation thread currently on charge. This is used so that event generation m...
Root of GENIE utility namespaces.
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Initial State information.