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

Generates resonance event (v+N->l+Resonance) kinematics. Is a concrete implementation of the EventRecordVisitorI interface. More...

#include <RESKinematicsGenerator.h>

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

Public Member Functions

 RESKinematicsGenerator ()
 
 RESKinematicsGenerator (string config)
 
 ~RESKinematicsGenerator ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
- 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...
 

Private Member Functions

void LoadConfig (void)
 
double ComputeMaxXSec (const Interaction *interaction) const
 

Private Attributes

TF2 * fEnvelope
 2-D envelope used for importance sampling More...
 
double fWcut
 Wcut parameter in DIS/RES join scheme. 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 double Energy (const Interaction *in) 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 resonance event (v+N->l+Resonance) kinematics. Is a concrete implementation of the EventRecordVisitorI interface.

Author
Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool
Created:
November 18, 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 29 of file RESKinematicsGenerator.h.

Constructor & Destructor Documentation

RESKinematicsGenerator::RESKinematicsGenerator ( )

Definition at line 38 of file RESKinematicsGenerator.cxx.

References fEnvelope.

38  :
39 KineGeneratorWithCache("genie::RESKinematicsGenerator")
40 {
41  fEnvelope = 0;
42 }
TF2 * fEnvelope
2-D envelope used for importance sampling
RESKinematicsGenerator::RESKinematicsGenerator ( string  config)

Definition at line 44 of file RESKinematicsGenerator.cxx.

References fEnvelope.

44  :
45 KineGeneratorWithCache("genie::RESKinematicsGenerator", config)
46 {
47  fEnvelope = 0;
48 }
TF2 * fEnvelope
2-D envelope used for importance sampling
RESKinematicsGenerator::~RESKinematicsGenerator ( )

Definition at line 50 of file RESKinematicsGenerator.cxx.

References fEnvelope.

51 {
52  if(fEnvelope) delete fEnvelope;
53 }
TF2 * fEnvelope
2-D envelope used for importance sampling

Member Function Documentation

double RESKinematicsGenerator::ComputeMaxXSec ( const Interaction interaction) const
privatevirtual

Implements genie::KineGeneratorWithCache.

Definition at line 269 of file RESKinematicsGenerator.cxx.

References genie::Interaction::ExclTag(), genie::KineGeneratorWithCache::fSafetyFactor, fWcut, genie::KineGeneratorWithCache::fXSecModel, genie::ProcessInfo::IsEM(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::controls::kMinQ2Limit, genie::utils::kinematics::electromagnetic::kMinQ2Limit, genie::XclsTag::KnownResonance(), genie::kPSWQD2fE, genie::kRfHitNucRest, LOG, genie::utils::res::Mass(), genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, genie::Interaction::PhaseSpace(), genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::utils::kinematics::Q2(), genie::KPhaseSpace::Q2Lim_W(), genie::XclsTag::Resonance(), genie::Kinematics::SetQ2(), genie::Kinematics::SetW(), genie::utils::kinematics::W(), genie::KPhaseSpace::WLim(), and genie::XSecAlgorithmI::XSec().

271 {
272 // Computes the maximum differential cross section in the requested phase
273 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
274 // method and the value is cached at a circular cache branch for retrieval
275 // during subsequent event generation.
276 // The computed max differential cross section does not need to be the exact
277 // maximum. The number used in the rejection method will be scaled up by a
278 // safety factor. But this needs to be fast - do not use a very fine grid.
279 
280  double max_xsec = 0.;
281 
282  const InitialState & init_state = interaction -> InitState();
283  double E = init_state.ProbeE(kRfHitNucRest);
284  bool is_em = interaction->ProcInfo().IsEM();
286 
287  double md;
288  if(!interaction->ExclTag().KnownResonance()) md=1.23;
289  else {
290  Resonance_t res = interaction->ExclTag().Resonance();
291  md=res::Mass(res);
292  }
293 
294  // ** 2-D Scan
295  const KPhaseSpace & kps = interaction->PhaseSpace();
296  Range1D_t rW = kps.WLim();
297 
298  int NW = 20;
299  double Wmin = rW.min + kASmallNum;
300  double Wmax = rW.max - kASmallNum;
301 
302  Wmax = TMath::Min(Wmax,fWcut);
303 
304  Wmin = TMath::Max(Wmin, md-.3);
305  Wmax = TMath::Min(Wmax, md+.3);
306 
307  if(Wmax-Wmin<0.05) { NW=1; Wmin=Wmax; }
308 
309  double dW = (NW>1) ? (Wmax-Wmin)/(NW-1) : 0.;
310 
311  for(int iw=0; iw<NW; iw++)
312  {
313  double W = Wmin + iw*dW;
314  interaction->KinePtr()->SetW(W);
315 
316  int NQ2 = 25;
317  int NQ2b = 4;
318 
319  Range1D_t rQ2 = kps.Q2Lim_W();
320  if( rQ2.max < Q2Thres || rQ2.min <=0 ) continue;
321  if( rQ2.max-rQ2.min<0.02 ) {NQ2=5; NQ2b=3;}
322 
323  double logQ2min = TMath::Log(rQ2.min+kASmallNum);
324  double logQ2max = TMath::Log(rQ2.max-kASmallNum);
325  double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
326  double xseclast = -1;
327  bool increasing = true;
328 
329  for(int iq2=0; iq2<NQ2; iq2++)
330  {
331  double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2);
332  interaction->KinePtr()->SetQ2(Q2);
333  double xsec = fXSecModel->XSec(interaction, kPSWQD2fE);
334  LOG("RESKinematics", pDEBUG)
335  << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec;
336  max_xsec = TMath::Max(xsec, max_xsec);
337  increasing = xsec-xseclast>=0;
338  xseclast=xsec;
339 
340  // once the cross section stops increasing, I reduce the step size and
341  // step backwards a little bit to handle cases that the max cross section
342  // is grossly underestimated (very peaky distribution & large step)
343  if(!increasing)
344  {
345  dlogQ2/=NQ2b;
346  for(int iq2b=0; iq2b<NQ2b; iq2b++)
347  {
348  Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
349  if(Q2 < rQ2.min) continue;
350  interaction->KinePtr()->SetQ2(Q2);
351  xsec = fXSecModel->XSec(interaction, kPSWQD2fE);
352  LOG("RESKinematics", pDEBUG)
353  << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec;
354  max_xsec = TMath::Max(xsec, max_xsec);
355  }
356  break;
357  }
358  } // Q2
359  }//W
360 
361  // Apply safety factor, since value retrieved from the cache might
362  // correspond to a slightly different energy
363  // Apply larger safety factor for smaller energies.
364  max_xsec *= ( (E<md) ? 2. : fSafetyFactor);
365 
366  return max_xsec;
367 }
double fWcut
Wcut parameter in DIS/RES join scheme.
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
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
A simple [min,max] interval for doubles.
Definition: Range1.h:42
bool KnownResonance(void) const
Definition: XclsTag.h:68
double Mass(Resonance_t res)
resonance mass (GeV)
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
enum genie::EResonance Resonance_t
static const double kMinQ2Limit
Definition: Controls.h:41
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
#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
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
bool IsEM(void) const
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
double ProbeE(RefFrame_t rf) const
Range1D_t WLim(void) const
W limits.
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
void RESKinematicsGenerator::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 227 of file RESKinematicsGenerator.cxx.

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

228 {
229  Algorithm::Configure(config);
230  this->LoadConfig();
231 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void RESKinematicsGenerator::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 233 of file RESKinematicsGenerator.cxx.

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

234 {
235  Algorithm::Configure(config);
236  this->LoadConfig();
237 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void RESKinematicsGenerator::LoadConfig ( void  )
private

Definition at line 239 of file RESKinematicsGenerator.cxx.

References genie::KineGeneratorWithCache::fEMin, fEnvelope, genie::KineGeneratorWithCache::fGenerateUniformly, genie::KineGeneratorWithCache::fMaxXSecDiffTolerance, genie::KineGeneratorWithCache::fSafetyFactor, fWcut, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), and genie::utils::kinematics::RESImportanceSamplingEnvelope().

Referenced by Configure().

240 {
241  // Safety factor for the maximum differential cross section
242  this->GetParamDef("MaxXSec-SafetyFactor", fSafetyFactor, 1.25);
243 
244  // Minimum energy for which max xsec would be cached, forcing explicit
245  // calculation for lower eneries
246  this->GetParamDef("Cache-MinEnergy", fEMin, 0.5);
247 
248  // Load Wcut used in DIS/RES join scheme
249  this->GetParam("Wcut", fWcut);
250 
251  // Maximum allowed fractional cross section deviation from maxim cross
252  // section used in rejection method
253  this->GetParamDef("MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999.);
254  assert(fMaxXSecDiffTolerance>=0);
255 
256  // Generate kinematics uniformly over allowed phase space and compute
257  // an event weight?
258  this->GetParamDef("UniformOverPhaseSpace", fGenerateUniformly, false);
259 
260  // Envelope employed when importance sampling is used
261  // (initialize with dummy range)
262  if(fEnvelope) delete fEnvelope;
263  fEnvelope = new TF2("res-envelope",
265  // stop ROOT from deleting this object of its own volition
266  gROOT->GetListOfFunctions()->Remove(fEnvelope);
267 }
double fWcut
Wcut parameter in DIS/RES join scheme.
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec&gt;maxxsec
TF2 * fEnvelope
2-D envelope used for importance sampling
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
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 RESImportanceSamplingEnvelope(double *x, double *par)
Definition: KineUtils.cxx:1372
void RESKinematicsGenerator::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 55 of file RESKinematicsGenerator.cxx.

References genie::KineGeneratorWithCache::AssertXSecLimits(), genie::Kinematics::ClearRunningValues(), genie::EventGeneratorI::CrossSectionAlg(), genie::GHepRecord::EventFlags(), genie::KineGeneratorWithCache::fGenerateUniformly, genie::KineGeneratorWithCache::fXSecModel, gQ2, genie::Target::HitNucP4(), genie::RunningThreadInfo::Instance(), genie::RandomGen::Instance(), genie::ProcessInfo::IsEM(), genie::controls::kASmallNum, genie::Interaction::KinePtr(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::kKineGenErr, genie::kKVW, genie::kPSWQ2fE, genie::kPSWQD2fE, genie::kRfHitNucRest, genie::controls::kRjMaxIterations, genie::KPhaseSpace::Limits(), LOG, genie::Range1D_t::max, genie::KineGeneratorWithCache::MaxXSec(), genie::Range1D_t::min, genie::Interaction::PhaseSpace(), genie::utils::kinematics::PhaseSpaceVolume(), pINFO, pNOTICE, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), pWARN, genie::utils::kinematics::Q2(), genie::KPhaseSpace::Q2Lim_W(), genie::utils::kinematics::Q2toQD2(), genie::utils::kinematics::QD2toQ2(), genie::RandomGen::RndKine(), genie::RunningThreadInfo::RunningThread(), genie::GHepRecord::SetDiffXSec(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::Kinematics::SetW(), genie::GHepRecord::SetWeight(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::GHepRecord::Summary(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), genie::InitialState::Tgt(), genie::utils::kinematics::W(), genie::GHepRecord::Weight(), genie::utils::kinematics::WQ2toXY(), genie::XSecAlgorithmI::XSec(), and genie::GHepRecord::XSec().

56 {
57  if(fGenerateUniformly) {
58  LOG("RESKinematics", pNOTICE)
59  << "Generating kinematics uniformly over the allowed phase space";
60  }
61 
62  //-- Get the random number generators
64 
65  //-- Access cross section algorithm for running thread
67  const EventGeneratorI * evg = rtinfo->RunningThread();
68  fXSecModel = evg->CrossSectionAlg();
69 
70  //-- Get the interaction from the GHEP record
71  Interaction * interaction = evrec->Summary();
72  interaction->SetBit(kISkipProcessChk);
73 
74  //-- EM or CC/NC process (different importance sampling methods)
75  bool is_em = interaction->ProcInfo().IsEM();
76 
77  //-- Compute the W limits
78  // (the physically allowed W's, unless an external cut is imposed)
79  const KPhaseSpace & kps = interaction->PhaseSpace();
80  Range1D_t W = kps.Limits(kKVW);
81 
82  if(W.max <=0 || W.min>=W.max) {
83  LOG("RESKinematics", pWARN) << "No available phase space";
84  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
86  exception.SetReason("No available phase space");
87  exception.SwitchOnFastForward();
88  throw exception;
89  }
90 
91  const InitialState & init_state = interaction -> InitState();
92  double E = init_state.ProbeE(kRfHitNucRest);
93 
94  //-- For the subsequent kinematic selection with the rejection method:
95  // Calculate the max differential cross section or retrieve it from the
96  // cache. Throw an exception and quit the evg thread if a non-positive
97  // value is found.
98  // If the kinematics are generated uniformly over the allowed phase
99  // space the max xsec is irrelevant
100  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
101 
102  //-- Try to select a valid W, Q2 pair using the rejection method
103  double dW = W.max - W.min;
104  double xsec = -1;
105 
106  unsigned int iter = 0;
107  bool accept = false;
108  while(1) {
109  iter++;
110  if(iter > kRjMaxIterations) {
111  LOG("RESKinematics", pWARN)
112  << "*** Could not select a valid (W,Q^2) pair after "
113  << iter << " iterations";
114  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
116  exception.SetReason("Couldn't select kinematics");
117  exception.SwitchOnFastForward();
118  throw exception;
119  }
120 
121  double gW = 0; // current hadronic invariant mass
122  double gQ2 = 0; // current momentum transfer
123  double gQD2 = 0; // tranformed Q2 to take out dipole form
124 
125  if(fGenerateUniformly)
126  {
127  //-- Generate a W uniformly in the kinematically allowed range.
128  // For the generated W, compute the Q2 range and generate a value
129  // uniformly over that range
130  gW = W.min + dW * rnd->RndKine().Rndm();
131  interaction->KinePtr()->SetW(gW);
132  Range1D_t Q2 = kps.Q2Lim_W();
133  if(Q2.max<=0. || Q2.min>=Q2.max) continue;
134  gQ2 = Q2.min + (Q2.max-Q2.min) * rnd->RndKine().Rndm();
135  interaction->SetBit(kISkipKinematicChk);
136  }
137  else
138  {
139  // neutrino scattering
140  // Selecting unweighted event kinematics using an importance sampling
141  // method. Q2 with be transformed to QD2 to take out the dipole form.
142  interaction->KinePtr()->SetW(W.min);
143  Range1D_t Q2 = kps.Q2Lim_W();
144  double Q2min = -99.;
145  if (is_em)
146  Q2min = Q2.min + kASmallNum;
147  else
148  Q2min = 0 + kASmallNum;
149  double Q2max = Q2.max - kASmallNum;
150 
151  // In unweighted mode - use transform that takes out the dipole form
152  double QD2min = utils::kinematics::Q2toQD2(Q2max);
153  double QD2max = utils::kinematics::Q2toQD2(Q2min);
154 
155  gW = W.min + dW * rnd->RndKine().Rndm();
156  gQD2 = QD2min + (QD2max - QD2min) * rnd->RndKine().Rndm();
157 
158  // QD2 -> Q2
159  gQ2 = utils::kinematics::QD2toQ2(gQD2);
160  } // uniformly over phase space?
161 
162  LOG("RESKinematics", pINFO) << "Trying: W = " << gW << ", Q2 = " << gQ2;
163 
164  //-- Set kinematics for current trial
165  interaction->KinePtr()->SetW(gW);
166  interaction->KinePtr()->SetQ2(gQ2);
167 
168  //-- Computing cross section for the current kinematics
169  xsec = fXSecModel->XSec(interaction, kPSWQD2fE);
170 
171  //-- Decide whether to accept the current kinematics
172  if(!fGenerateUniformly)
173  {
174  // unified neutrino / electron scattering
175  double t = xsec_max * rnd->RndKine().Rndm();
176  this->AssertXSecLimits(interaction, xsec, xsec_max);
177  accept = (t < xsec);
178  } // charged lepton or neutrino scattering?
179  else
180  {
181  accept = (xsec>0);
182  } // uniformly over phase space
183 
184  //-- If the generated kinematics are accepted, finish-up module's job
185  if(accept) {
186  LOG("RESKinematics", pINFO)
187  << "Selected: W = " << gW << ", Q2 = " << gQ2;
188  // reset 'trust' bits
189  interaction->ResetBit(kISkipProcessChk);
190  interaction->ResetBit(kISkipKinematicChk);
191 
192  // compute x,y for selected W,Q2
193  // note: hit nucleon can be off the mass-shell
194  double gx=-1, gy=-1;
195  double M = init_state.Tgt().HitNucP4().M();
196  kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy);
197 
198  // set the cross section for the selected kinematics
199  evrec->SetDiffXSec(xsec,kPSWQ2fE);
200 
201  // for uniform kinematics, compute an event weight as
202  // wght = (phase space volume)*(differential xsec)/(event total xsec)
203  if(fGenerateUniformly) {
204  double vol = kinematics::PhaseSpaceVolume(interaction,kPSWQ2fE);
205  double totxsec = evrec->XSec();
206  double wght = (vol/totxsec)*xsec;
207  LOG("RESKinematics", pNOTICE) << "Kinematics wght = "<< wght;
208 
209  // apply computed weight to the current event weight
210  wght *= evrec->Weight();
211  LOG("RESKinematics", pNOTICE) << "Current event wght = " << wght;
212  evrec->SetWeight(wght);
213  }
214 
215  // lock selected kinematics & clear running values
216  interaction->KinePtr()->SetQ2(gQ2, true);
217  interaction->KinePtr()->SetW (gW, true);
218  interaction->KinePtr()->Setx (gx, true);
219  interaction->KinePtr()->Sety (gy, true);
220  interaction->KinePtr()->ClearRunningValues();
221 
222  return;
223  } // accept
224  } // iterations
225 }
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
bool fGenerateUniformly
uniform over allowed phase space + event weight?
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 void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
Defines the EventGeneratorI interface.
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
double PhaseSpaceVolume(const Interaction *const i, KinePhaseSpace_t ps)
Definition: KineUtils.cxx:36
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
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
double QD2toQ2(double QD2)
Definition: KineUtils.cxx:1071
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1132
Kinematical phase space.
Definition: KPhaseSpace.h:33
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
double Q2toQD2(double Q2)
Definition: KineUtils.cxx:1061
#define pWARN
Definition: Messenger.h:60
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.
bool IsEM(void) const
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
static RunningThreadInfo * Instance(void)
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
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
double gQ2
Definition: gtestDISSF.cxx:55
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
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
const EventGeneratorI * RunningThread(void)
double ProbeE(RefFrame_t rf) const
Keep info on the event generation thread currently on charge. This is used so that event generation m...
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48

Member Data Documentation

TF2* genie::RESKinematicsGenerator::fEnvelope
mutableprivate

2-D envelope used for importance sampling

Definition at line 48 of file RESKinematicsGenerator.h.

Referenced by LoadConfig(), RESKinematicsGenerator(), and ~RESKinematicsGenerator().

double genie::RESKinematicsGenerator::fWcut
private

Wcut parameter in DIS/RES join scheme.

Definition at line 49 of file RESKinematicsGenerator.h.

Referenced by ComputeMaxXSec(), and LoadConfig().


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