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:
genie::KineGeneratorWithCache genie::EventRecordVisitorI genie::Algorithm

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)
 
- Protected Member Functions inherited from genie::KineGeneratorWithCache
 KineGeneratorWithCache ()
 
 KineGeneratorWithCache (string name)
 
 KineGeneratorWithCache (string name, string config)
 
 ~KineGeneratorWithCache ()
 
virtual double MaxXSec (GHepRecord *evrec) const
 
virtual double FindMaxXSec (const Interaction *in) const
 
virtual void CacheMaxXSec (const Interaction *in, double xsec) const
 
virtual double Energy (const Interaction *in) const
 
virtual CacheBranchFxAccessCacheBranch (const Interaction *in) 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
 
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
 maxxsec -> maxxsec * safety_factor 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< boolfOwnerships
 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 <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

November 18, 2004

Copyright (c) 2003-2020, 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.

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.

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

Definition at line 50 of file RESKinematicsGenerator.cxx.

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 313 of file RESKinematicsGenerator.cxx.

315 {
316 // Computes the maximum differential cross section in the requested phase
317 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
318 // method and the value is cached at a circular cache branch for retrieval
319 // during subsequent event generation.
320 // The computed max differential cross section does not need to be the exact
321 // maximum. The number used in the rejection method will be scaled up by a
322 // safety factor. But this needs to be fast - do not use a very fine grid.
323 
324  double max_xsec = 0.;
325 
326  const InitialState & init_state = interaction -> InitState();
327  double E = init_state.ProbeE(kRfHitNucRest);
328  bool is_em = interaction->ProcInfo().IsEM();
330 
331 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
332  LOG("RESKinematics", pDEBUG) << "Scanning phase space for E= " << E;
333 #endif
334  //bool scan1d = (E>1.0);
335  bool scan1d = false;
336 
337  double md;
338  if(!interaction->ExclTag().KnownResonance()) md=1.23;
339  else {
340  Resonance_t res = interaction->ExclTag().Resonance();
341  md=res::Mass(res);
342  }
343 
344  if(scan1d) {
345  // ** 1-D Scan
346  //
347 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
348  LOG("RESKinematics", pDEBUG)
349  << "Will search for max{xsec} along W(=MRes) = " << md;
350 #endif
351  // Set W around the value where d^2xsec/dWdQ^2 peaks
352  interaction->KinePtr()->SetW(md);
353 
354  const KPhaseSpace & kps = interaction->PhaseSpace();
355  Range1D_t rQ2 = kps.Q2Lim_W();
356  if( rQ2.max < Q2Thres || rQ2.min <=0 ) return 0.;
357 
358  int NQ2 = 25;
359  int NQ2b = 5;
360  double logQ2min = TMath::Log(rQ2.min+kASmallNum);
361  double logQ2max = TMath::Log(rQ2.max-kASmallNum);
362  double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
363 
364  double xseclast = -1;
365  bool increasing = true;
366 
367  for(int iq2=0; iq2<NQ2; iq2++) {
368  double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2);
369  interaction->KinePtr()->SetQ2(Q2);
370  double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
371 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
372  LOG("RESKinematics", pDEBUG)
373  << "xsec(W= " << md << ", Q2= " << Q2 << ") = " << xsec;
374 #endif
375  max_xsec = TMath::Max(xsec, max_xsec);
376  increasing = xsec-xseclast>=0;
377  xseclast=xsec;
378 
379  // once the cross section stops increasing, I reduce the step size and
380  // step backwards a little bit to handle cases that the max cross section
381  // is grossly underestimated (very peaky distribution & large step)
382  if(!increasing) {
383  dlogQ2/=NQ2b;
384  for(int iq2b=0; iq2b<NQ2b; iq2b++) {
385  Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
386  if(Q2 < rQ2.min) continue;
387  interaction->KinePtr()->SetQ2(Q2);
388  xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
389 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
390  LOG("RESKinematics", pDEBUG)
391  << "xsec(W= " << md << ", Q2= " << Q2 << ") = " << xsec;
392 #endif
393  max_xsec = TMath::Max(xsec, max_xsec);
394  }
395  break;
396  }
397  } // Q2
398 
399  } else {
400 
401  // ** 2-D Scan
402  //
403  const KPhaseSpace & kps = interaction->PhaseSpace();
404  Range1D_t rW = kps.WLim();
405 
406  int NW = 20;
407  double Wmin = rW.min + kASmallNum;
408  double Wmax = rW.max - kASmallNum;
409 
410  Wmax = TMath::Min(Wmax,fWcut);
411 
412  Wmin = TMath::Max(Wmin, md-.3);
413  Wmax = TMath::Min(Wmax, md+.3);
414 
415  if(Wmax-Wmin<0.05) { NW=1; Wmin=Wmax; }
416 
417  double dW = (NW>1) ? (Wmax-Wmin)/(NW-1) : 0.;
418 
419  for(int iw=0; iw<NW; iw++) {
420  double W = Wmin + iw*dW;
421  interaction->KinePtr()->SetW(W);
422 
423  int NQ2 = 25;
424  int NQ2b = 4;
425 
426  Range1D_t rQ2 = kps.Q2Lim_W();
427  if( rQ2.max < Q2Thres || rQ2.min <=0 ) continue;
428  if( rQ2.max-rQ2.min<0.02 ) {NQ2=5; NQ2b=3;}
429 
430  double logQ2min = TMath::Log(rQ2.min+kASmallNum);
431  double logQ2max = TMath::Log(rQ2.max-kASmallNum);
432  double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
433  double xseclast = -1;
434  bool increasing = true;
435 
436  for(int iq2=0; iq2<NQ2; iq2++) {
437  double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2);
438  interaction->KinePtr()->SetQ2(Q2);
439  double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
440  LOG("RESKinematics", pDEBUG)
441  << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec;
442  max_xsec = TMath::Max(xsec, max_xsec);
443  increasing = xsec-xseclast>=0;
444  xseclast=xsec;
445 
446  // once the cross section stops increasing, I reduce the step size and
447  // step backwards a little bit to handle cases that the max cross section
448  // is grossly underestimated (very peaky distribution & large step)
449  if(!increasing) {
450  dlogQ2/=NQ2b;
451  for(int iq2b=0; iq2b<NQ2b; iq2b++) {
452  Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
453  if(Q2 < rQ2.min) continue;
454  interaction->KinePtr()->SetQ2(Q2);
455  xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
456  LOG("RESKinematics", pDEBUG)
457  << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec;
458  max_xsec = TMath::Max(xsec, max_xsec);
459  }
460  break;
461  }
462  } // Q2
463  }//W
464  }//2d scan
465 
466  // Apply safety factor, since value retrieved from the cache might
467  // correspond to a slightly different energy
468  // Apply larger safety factor for smaller energies.
469  max_xsec *= ( (E<md) ? 2. : fSafetyFactor);
470 
471 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
472  LOG("RESKinematics", pDEBUG) << interaction->AsString();
473  LOG("RESKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
474  LOG("RESKinematics", pDEBUG) << "Computed using " << fXSecModel->Id();
475 #endif
476 
477  return max_xsec;
478 }
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:1064
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
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.
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
string AsString(void) const
static const double kMinQ2Limit
Definition: Controls.h:41
#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
bool IsEM(void) const
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:97
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
E
Definition: 018_def.c:13
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 271 of file RESKinematicsGenerator.cxx.

272 {
273  Algorithm::Configure(config);
274  this->LoadConfig();
275 }
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 277 of file RESKinematicsGenerator.cxx.

278 {
280  this->LoadConfig();
281 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void RESKinematicsGenerator::LoadConfig ( void  )
private

Definition at line 283 of file RESKinematicsGenerator.cxx.

284 {
285  // Safety factor for the maximum differential cross section
286  this->GetParamDef("MaxXSec-SafetyFactor", fSafetyFactor, 1.25);
287 
288  // Minimum energy for which max xsec would be cached, forcing explicit
289  // calculation for lower eneries
290  this->GetParamDef("Cache-MinEnergy", fEMin, 0.5);
291 
292  // Load Wcut used in DIS/RES join scheme
293  this->GetParam("Wcut", fWcut);
294 
295  // Maximum allowed fractional cross section deviation from maxim cross
296  // section used in rejection method
297  this->GetParamDef("MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999.);
298  assert(fMaxXSecDiffTolerance>=0);
299 
300  // Generate kinematics uniformly over allowed phase space and compute
301  // an event weight?
302  this->GetParamDef("UniformOverPhaseSpace", fGenerateUniformly, false);
303 
304  // Envelope employed when importance sampling is used
305  // (initialize with dummy range)
306  if(fEnvelope) delete fEnvelope;
307  fEnvelope = new TF2("res-envelope",
309  // stop ROOT from deleting this object of its own volition
310  gROOT->GetListOfFunctions()->Remove(fEnvelope);
311 }
double fWcut
Wcut parameter in DIS/RES join scheme.
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>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:1359
void RESKinematicsGenerator::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 55 of file RESKinematicsGenerator.cxx.

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  //costas says this is bad if(W.max>1.7) W.max=1.7;
82 //assert(W.min>=0. && W.min<W.max);
83 
84  if(W.max <=0 || W.min>=W.max) {
85  LOG("RESKinematics", pWARN) << "No available phase space";
86  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
88  exception.SetReason("No available phase space");
89  exception.SwitchOnFastForward();
90  throw exception;
91  }
92 
93  const InitialState & init_state = interaction -> InitState();
94  double E = init_state.ProbeE(kRfHitNucRest);
95  // double M = init_state.Tgt().HitNucP4().M();
96  // double ml = interaction->FSPrimLepton()->Mass();
97 
98  //-- For the subsequent kinematic selection with the rejection method:
99  // Calculate the max differential cross section or retrieve it from the
100  // cache. Throw an exception and quit the evg thread if a non-positive
101  // value is found.
102  // If the kinematics are generated uniformly over the allowed phase
103  // space the max xsec is irrelevant
104  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
105 
106  //-- Try to select a valid W, Q2 pair using the rejection method
107  double dW = W.max - W.min;
108  double xsec = -1;
109 
110  unsigned int iter = 0;
111  bool accept = false;
112  while(1) {
113  iter++;
114  if(iter > kRjMaxIterations) {
115  LOG("RESKinematics", pWARN)
116  << "*** Could not select a valid (W,Q^2) pair after "
117  << iter << " iterations";
118  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
120  exception.SetReason("Couldn't select kinematics");
121  exception.SwitchOnFastForward();
122  throw exception;
123  }
124 
125  double gW = 0; // current hadronic invariant mass
126  double gQ2 = 0; // current momentum transfer
127  double gQD2 = 0; // tranformed Q2 to take out dipole form
128 
129  if(fGenerateUniformly) {
130  //-- Generate a W uniformly in the kinematically allowed range.
131  // For the generated W, compute the Q2 range and generate a value
132  // uniformly over that range
133  gW = W.min + dW * rnd->RndKine().Rndm();
134  Range1D_t Q2 = kps.Q2Lim_W();
135  if(Q2.max<=0. || Q2.min>=Q2.max) continue;
136  gQ2 = Q2.min + (Q2.max-Q2.min) * rnd->RndKine().Rndm();
137 
138  interaction->SetBit(kISkipKinematicChk);
139 
140  } else {
141 
142 
143  // > neutrino scattering
144  // Selecting unweighted event kinematics using an importance sampling
145  // method. Q2 with be transformed to QD2 to take out the dipole form.
146  // An importance sampling envelope will be constructed for W.
147  // first pass, configure the sampling envelope
148  if(iter==1) {
149  LOG("RESKinematics", pINFO) << "Initializing the sampling envelope";
150  if(!fEnvelope) {
151  LOG("RESKinematics", pFATAL) << "Null sampling envelope!";
152  exit(1);
153  }
154  interaction->KinePtr()->SetW(W.min);
155  Range1D_t Q2 = kps.Q2Lim_W();
156  double Q2min = -99.;
157  if (is_em) { Q2min = Q2.min + kASmallNum; }
158  else { Q2min = 0 + kASmallNum; }
159  double Q2max = Q2.max - kASmallNum;
160 
161  // In unweighted mode - use transform that takes out the dipole form
162  double QD2min = utils::kinematics::Q2toQD2(Q2max);
163  double QD2max = utils::kinematics::Q2toQD2(Q2min);
164 
165 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
166  LOG("RESKinematics", pDEBUG)
167  << "Q^2: [" << Q2min << ", " << Q2max << "] => "
168  << "QD^2: [" << QD2min << ", " << QD2max << "]";
169 #endif
170  double mR, gR;
171  if(!interaction->ExclTag().KnownResonance()) {
172  mR = 1.2;
173  gR = 0.6;
174  } else {
175  Resonance_t res = interaction->ExclTag().Resonance();
176  mR = res::Mass(res);
177  gR = (E>mR) ? 0.220 : 0.400;
178  }
179 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
180  LOG("RESKinematics", pDEBUG)
181  << "(m,g) = (" << mR << ", " << gR
182  << "), max(xsec,W) = (" << xsec_max << ", " << W.max << ")";
183 #endif
184  fEnvelope->SetRange(QD2min,W.min,QD2max,W.max); // range
185  fEnvelope->SetParameter(0, mR); // resonance mass
186  fEnvelope->SetParameter(1, gR); // resonance width
187  fEnvelope->SetParameter(2, xsec_max); // max differential xsec
188  fEnvelope->SetParameter(3, W.max); // kinematically allowed Wmax
189  }// first pass
190 
191  // Generate W,QD2 using the 2-D envelope as PDF
192  fEnvelope->GetRandom2(gQD2,gW);
193 
194  // QD2 -> Q2
195  gQ2 = utils::kinematics::QD2toQ2(gQD2);
196  } // uniformly over phase space?
197 
198  LOG("RESKinematics", pINFO) << "Trying: W = " << gW << ", Q2 = " << gQ2;
199 
200  //-- Set kinematics for current trial
201  interaction->KinePtr()->SetW(gW);
202  interaction->KinePtr()->SetQ2(gQ2);
203 
204  //-- Computing cross section for the current kinematics
205  xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
206 
207  //-- Decide whether to accept the current kinematics
208  if(!fGenerateUniformly) {
209 
210  // unified neutrino / electron scattering
211  double max = fEnvelope->Eval(gQD2, gW);
212  double t = max * rnd->RndKine().Rndm();
213  double J = kinematics::Jacobian(interaction,kPSWQ2fE,kPSWQD2fE);
214 
215  this->AssertXSecLimits(interaction, xsec, max);
216 
217 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
218  LOG("RESKinematics", pDEBUG)
219  << "xsec= " << xsec << ", J= " << J << ", Rnd= " << t;
220 #endif
221  accept = (t < J*xsec);
222  } // charged lepton or neutrino scattering?
223  else {
224  accept = (xsec>0);
225  } // uniformly over phase space
226 
227  //-- If the generated kinematics are accepted, finish-up module's job
228  if(accept) {
229  LOG("RESKinematics", pINFO)
230  << "Selected: W = " << gW << ", Q2 = " << gQ2;
231  // reset 'trust' bits
232  interaction->ResetBit(kISkipProcessChk);
233  interaction->ResetBit(kISkipKinematicChk);
234 
235  // compute x,y for selected W,Q2
236  // note: hit nucleon can be off the mass-shell
237  double gx=-1, gy=-1;
238  double M = init_state.Tgt().HitNucP4().M();
239  //double M = init_state.Tgt().HitNucMass();
240  kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy);
241 
242  // set the cross section for the selected kinematics
243  evrec->SetDiffXSec(xsec,kPSWQ2fE);
244 
245  // for uniform kinematics, compute an event weight as
246  // wght = (phase space volume)*(differential xsec)/(event total xsec)
247  if(fGenerateUniformly) {
248  double vol = kinematics::PhaseSpaceVolume(interaction,kPSWQ2fE);
249  double totxsec = evrec->XSec();
250  double wght = (vol/totxsec)*xsec;
251  LOG("RESKinematics", pNOTICE) << "Kinematics wght = "<< wght;
252 
253  // apply computed weight to the current event weight
254  wght *= evrec->Weight();
255  LOG("RESKinematics", pNOTICE) << "Current event wght = " << wght;
256  evrec->SetWeight(wght);
257  }
258 
259  // lock selected kinematics & clear running values
260  interaction->KinePtr()->SetQ2(gQ2, true);
261  interaction->KinePtr()->SetW (gW, true);
262  interaction->KinePtr()->Setx (gx, true);
263  interaction->KinePtr()->Sety (gy, true);
264  interaction->KinePtr()->ClearRunningValues();
265 
266  return;
267  } // accept
268  } // iterations
269 }
virtual double MaxXSec(GHepRecord *evrec) const
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
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 double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
A simple [min,max] interval for doubles.
Definition: Range1.h:42
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
bool KnownResonance(void) const
Definition: XclsTag.h:68
#define pFATAL
Definition: Messenger.h:56
Defines the EventGeneratorI interface.
double Mass(Resonance_t res)
resonance mass (GeV)
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
enum genie::EResonance Resonance_t
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
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:1058
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1119
Kinematical phase space.
Definition: KPhaseSpace.h:33
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
static int max(int a, int b)
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
TF2 * fEnvelope
2-D envelope used for importance sampling
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
double Q2toQD2(double Q2)
Definition: KineUtils.cxx:1048
#define pWARN
Definition: Messenger.h:60
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
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
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
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...
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

TF2* genie::RESKinematicsGenerator::fEnvelope
mutableprivate

2-D envelope used for importance sampling

Definition at line 48 of file RESKinematicsGenerator.h.

double genie::RESKinematicsGenerator::fWcut
private

Wcut parameter in DIS/RES join scheme.

Definition at line 49 of file RESKinematicsGenerator.h.


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