Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::HEDISKinematicsGenerator Class Reference

Generates values for the kinematic variables describing HEDIS v interaction events. Is a concrete implementation of the EventRecordVisitorI interface. More...

#include <HEDISKinematicsGenerator.h>

Inheritance diagram for genie::HEDISKinematicsGenerator:
genie::KineGeneratorWithCache genie::EventRecordVisitorI genie::Algorithm

Public Member Functions

 HEDISKinematicsGenerator ()
 
 HEDISKinematicsGenerator (string config)
 
 ~HEDISKinematicsGenerator ()
 
double ComputeMaxXSec (const Interaction *interaction) const
 
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

double Scan (Interaction *interaction, Range1D_t xrange, Range1D_t Q2range, int NKnotsQ2, int NKnotsX, double ME2, double &x_scan, double &Q2_scan) const
 
void LoadConfig (void)
 

Private Attributes

int fWideNKnotsX
 
int fWideNKnotsQ2
 
double fWideRange
 
int fFineNKnotsX
 
int fFineNKnotsQ2
 
double fFineRange
 
double fSFXmin
 minimum value of x for which SF tables are computed More...
 
double fSFQ2min
 minimum value of Q2 for which SF tables are computed More...
 
double fSFQ2max
 maximum value of Q2 for which SF tables are computed 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 values for the kinematic variables describing HEDIS v interaction events. Is a concrete implementation of the EventRecordVisitorI interface.

Max Xsec are precomputed as the total xsec and they are stored in ascii files. This saves a lot of computational time.

Author
Alfonso Garcia <alfonsog nikhef.nl> NIKHEF

August 28, 2019

Copyright (c) 2003-2018, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 31 of file HEDISKinematicsGenerator.h.

Constructor & Destructor Documentation

HEDISKinematicsGenerator::HEDISKinematicsGenerator ( )

Definition at line 42 of file HEDISKinematicsGenerator.cxx.

42  :
43 KineGeneratorWithCache("genie::HEDISKinematicsGenerator")
44 {
45 
46 }
HEDISKinematicsGenerator::HEDISKinematicsGenerator ( string  config)

Definition at line 48 of file HEDISKinematicsGenerator.cxx.

48  :
49 KineGeneratorWithCache("genie::HEDISKinematicsGenerator", config)
50 {
51 
52 }
static Config * config
Definition: config.cpp:1054
HEDISKinematicsGenerator::~HEDISKinematicsGenerator ( )

Definition at line 54 of file HEDISKinematicsGenerator.cxx.

55 {
56 
57 }

Member Function Documentation

double HEDISKinematicsGenerator::ComputeMaxXSec ( const Interaction interaction) const
virtual

Implements genie::KineGeneratorWithCache.

Definition at line 219 of file HEDISKinematicsGenerator.cxx.

220 {
221  return 0;
222 }
void HEDISKinematicsGenerator::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 224 of file HEDISKinematicsGenerator.cxx.

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

231 {
233  this->LoadConfig();
234 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void HEDISKinematicsGenerator::LoadConfig ( void  )
private

Definition at line 236 of file HEDISKinematicsGenerator.cxx.

237 {
238 
239  //-- Safety factor for the maximum differential cross section
240  GetParamDef("MaxXSec-SafetyFactor", fSafetyFactor, 2. ) ;
241  //-- Maximum allowed fractional cross section deviation from maxim cross
242  // section used in rejection method
243  GetParamDef("MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
244  assert(fMaxXSecDiffTolerance>=0);
245 
246  GetParamDef("ScanWide-NKnotsX", fWideNKnotsX, 10 ) ;
247  GetParamDef("ScanWide-NKnotsQ2", fWideNKnotsQ2, 10 ) ;
248  GetParamDef("ScanWide-Range", fWideRange, 1.1 ) ;
249  GetParamDef("ScanFine-NKnotsX", fFineNKnotsX, 10 ) ;
250  GetParamDef("ScanFine-NKnotsQ2", fFineNKnotsQ2, 10 ) ;
251  GetParamDef("ScanFine-Range", fFineRange, 10. ) ;
252 
253  // Limits from the SF tables that are useful to reduce computation
254  // time of the max cross section
255  GetParam("XGrid-Min", fSFXmin ) ;
256  GetParam("Q2Grid-Min", fSFQ2min ) ;
257  GetParam("Q2Grid-Max", fSFQ2max ) ;
258 
259 }
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
double fSFXmin
minimum value of x for which SF tables are computed
double fSFQ2min
minimum value of Q2 for which SF tables are computed
double fSFQ2max
maximum value of Q2 for which SF tables are computed
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
void HEDISKinematicsGenerator::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 59 of file HEDISKinematicsGenerator.cxx.

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
71  Interaction * interaction = evrec->Summary();
72  interaction->SetBit(kISkipProcessChk);
73 
74  //-- Get neutrino energy and hit 'nucleon mass'
75  const InitialState & init_state = interaction->InitState();
76  double Ev = init_state.ProbeE(kRfLab);
77  double M = init_state.Tgt().HitNucP4().M(); // can be off m-shell
78 
79  //-- Get the physical W range
80  const KPhaseSpace & kps = interaction->PhaseSpace();
81  Range1D_t xl = kps.XLim();
82  Range1D_t Q2l = kps.Q2Lim();
83 
84  //-- x and y lower limit restrict by limits in SF tables
85  Q2l.min = TMath::Max(Q2l.min,fSFQ2min);
86  Q2l.max = TMath::Min(Q2l.max,fSFQ2max);
87  xl.min = TMath::Max(TMath::Max(xl.min,Q2l.min/2./M/Ev),fSFXmin);
88 
89  LOG("HEDISKinematics", pNOTICE) << "x: [" << xl.min << ", " << xl.max << "]";
90  LOG("HEDISKinematics", pNOTICE) << "log10Q2: [" << Q2l.min << ", " << Q2l.max << "]";
91 
92  //-- For the subsequent kinematic selection with the rejection method:
93 
94  //Scan through a wide region to find the maximum
95  Range1D_t xrange_wide(xl.min*fWideRange,xl.max/fWideRange);
96  Range1D_t Q2range_wide(Q2l.min*fWideRange,Q2l.max/fWideRange);
97  double x_wide = 0.;
98  double Q2_wide = 0.;
99  double xsec_wide = this->Scan(interaction,xrange_wide,Q2range_wide,fWideNKnotsX,fWideNKnotsQ2,2*M*Ev,x_wide,Q2_wide);
100 
101  //Scan through a fine region to find the maximum
102  Range1D_t xrange_fine(TMath::Max(x_wide/fFineRange,xrange_wide.min),TMath::Min(x_wide*fFineRange,xrange_wide.max));
103  Range1D_t Q2range_fine(TMath::Max(Q2_wide/fFineRange,Q2range_wide.min),TMath::Min(Q2_wide*fFineRange,Q2range_wide.max));
104  double x_fine = 0.;
105  double Q2_fine = 0.;
106  double xsec_fine = this->Scan(interaction,xrange_fine,Q2range_fine,fFineNKnotsX,fFineNKnotsQ2,2*M*Ev,x_fine,Q2_fine);
107 
108  //Apply safety factor
109  double xsec_max = fSafetyFactor * TMath::Max(xsec_wide,xsec_fine);
110 
111  //-- Try to select a valid (x,y) pair using the rejection method
112  double log10xmin = TMath::Log10(xl.min);
113  double log10xmax = TMath::Log10(xl.max);
114  double log10Q2min = TMath::Log10(Q2l.min);
115  double log10Q2max = TMath::Log10(Q2l.max);
116 
117  double dlog10x = log10xmax - log10xmin;
118  double dlog10Q2 = log10Q2max - log10Q2min;
119 
120  double gx=-1, gQ2=-1, xsec=-1;
121 
122  unsigned int iter = 0;
123  bool accept = false;
124  while(1) {
125  iter++;
126  if(iter > kRjMaxIterations) {
127  LOG("HEDISKinematics", pWARN)
128  << " Couldn't select kinematics after " << iter << " iterations";
129  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
131  exception.SetReason("Couldn't select kinematics");
132  exception.SwitchOnFastForward();
133  throw exception;
134  }
135 
136  gx = TMath::Power( 10., log10xmin + dlog10x * rnd->RndKine().Rndm() );
137  gQ2 = TMath::Power( 10., log10Q2min + dlog10Q2 * rnd->RndKine().Rndm() );
138 
139  interaction->KinePtr()->Setx(gx);
140  interaction->KinePtr()->SetQ2(gQ2);
141  kinematics::UpdateWYFromXQ2(interaction);
142 
143  LOG("HEDISKinematics", pDEBUG)
144  << "Trying: x = " << gx << ", Q2 = " << gQ2
145  << " (W = " << interaction->KinePtr()->W() << ","
146  << " y = " << interaction->KinePtr()->y() << ")";
147 
148  //-- compute the cross section for current kinematics
149  xsec = fXSecModel->XSec(interaction, kPSlog10xlog10Q2fE);
150 
151  //-- decide whether to accept the current kinematics
152  this->AssertXSecLimits(interaction, xsec, xsec_max);
153 
154  double t = xsec_max * rnd->RndKine().Rndm();
155 
156  LOG("HEDISKinematics", pDEBUG) << "xsec= " << xsec << ", Rnd= " << t;
157 
158  accept = (t < xsec);
159 
160  //-- If the generated kinematics are accepted, finish-up module's job
161  if(accept) {
162  LOG("HEDISKinematics", pNOTICE)
163  << "Selected: x = " << gx << ", Q2 = " << gQ2
164  << " (W = " << interaction->KinePtr()->W() << ","
165  << " (Y = " << interaction->KinePtr()->y() << ")";
166 
167  // reset trust bits
168  interaction->ResetBit(kISkipProcessChk);
169  interaction->ResetBit(kISkipKinematicChk);
170 
171  // set the cross section for the selected kinematics
172  evrec->SetDiffXSec(xsec,kPSxQ2fE);
173 
174  // lock selected kinematics & clear running values
175  interaction->KinePtr()->SetW (interaction->KinePtr()->W(), true);
176  interaction->KinePtr()->Sety (interaction->KinePtr()->y(), true);
177  interaction->KinePtr()->SetQ2(gQ2, true);
178  interaction->KinePtr()->Setx (gx, true);
179  interaction->KinePtr()->ClearRunningValues();
180  return;
181  }
182  } // iterations
183 }
double W(bool selected=false) const
Definition: Kinematics.cxx:157
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
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
double Scan(Interaction *interaction, Range1D_t xrange, Range1D_t Q2range, int NKnotsQ2, int NKnotsX, double ME2, double &x_scan, double &Q2_scan) 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.
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.
double fSFXmin
minimum value of x for which SF tables are computed
Range1D_t Q2Lim(void) const
Q2 limits.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
void UpdateWYFromXQ2(const Interaction *in)
Definition: KineUtils.cxx:1313
double y(bool selected=false) const
Definition: Kinematics.cxx:112
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 fSFQ2min
minimum value of Q2 for which SF tables are computed
Kinematical phase space.
Definition: KPhaseSpace.h:33
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
#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
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
double gQ2
Definition: gtestDISSF.cxx:55
Range1D_t XLim(void) const
x limits
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
double fSFQ2max
maximum value of Q2 for which SF tables are computed
#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
double HEDISKinematicsGenerator::Scan ( Interaction interaction,
Range1D_t  xrange,
Range1D_t  Q2range,
int  NKnotsQ2,
int  NKnotsX,
double  ME2,
double &  x_scan,
double &  Q2_scan 
) const
private

Definition at line 185 of file HEDISKinematicsGenerator.cxx.

186 {
187 
188  double xsec_scan = 0.;
189  Q2_scan = 0.;
190  x_scan = 0.;
191 
192  double stepQ2 = TMath::Power(Q2range.max/Q2range.min,1./(NKnotsQ2-1));
193 
194  for ( int iq=0; iq<NKnotsQ2; iq++) {
195  double Q2_aux = Q2range.min*TMath::Power(stepQ2,iq);
196  xrange.min = TMath::Max(xrange.min,Q2_aux/ME2);
197  double stepx = TMath::Power(xrange.max/xrange.min,1./(NKnotsX-1));
198  for ( int ix=0; ix<NKnotsX; ix++) {
199  double x_aux = xrange.min*TMath::Power(stepx,ix);
200  interaction->KinePtr()->Setx(x_aux);
201  interaction->KinePtr()->SetQ2(Q2_aux);
202  kinematics::UpdateWYFromXQ2(interaction);
203  double xsec_aux = fXSecModel->XSec(interaction, kPSlog10xlog10Q2fE);
204  LOG("HEDISKinematics", pDEBUG) << "x = " << x_aux << " , Q2 = " << Q2_aux << ", xsec = " << xsec_aux;
205  if (xsec_aux>xsec_scan) {
206  xsec_scan = xsec_aux;
207  Q2_scan = Q2_aux;
208  x_scan = x_aux;
209  }
210  }
211  }
212 
213  LOG("HEDISKinematics", pNOTICE) << "scan -> x = " << x_scan << " , Q2 = " << Q2_scan << ", xsec = " << xsec_scan;
214 
215  return xsec_scan;
216 
217 }
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.
void UpdateWYFromXQ2(const Interaction *in)
Definition: KineUtils.cxx:1313
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
double max
Definition: Range1.h:53
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

int genie::HEDISKinematicsGenerator::fFineNKnotsQ2
private

Definition at line 58 of file HEDISKinematicsGenerator.h.

int genie::HEDISKinematicsGenerator::fFineNKnotsX
private

Definition at line 57 of file HEDISKinematicsGenerator.h.

double genie::HEDISKinematicsGenerator::fFineRange
private

Definition at line 59 of file HEDISKinematicsGenerator.h.

double genie::HEDISKinematicsGenerator::fSFQ2max
private

maximum value of Q2 for which SF tables are computed

Definition at line 63 of file HEDISKinematicsGenerator.h.

double genie::HEDISKinematicsGenerator::fSFQ2min
private

minimum value of Q2 for which SF tables are computed

Definition at line 62 of file HEDISKinematicsGenerator.h.

double genie::HEDISKinematicsGenerator::fSFXmin
private

minimum value of x for which SF tables are computed

Definition at line 61 of file HEDISKinematicsGenerator.h.

int genie::HEDISKinematicsGenerator::fWideNKnotsQ2
private

Definition at line 55 of file HEDISKinematicsGenerator.h.

int genie::HEDISKinematicsGenerator::fWideNKnotsX
private

Definition at line 54 of file HEDISKinematicsGenerator.h.

double genie::HEDISKinematicsGenerator::fWideRange
private

Definition at line 56 of file HEDISKinematicsGenerator.h.


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