GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DFRKinematicsGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <cfloat>
12 
13 #include <TMath.h>
14 
16 #include "Framework/Conventions/GBuild.h"
32 
33 using namespace genie;
34 using namespace genie::controls;
35 using namespace genie::constants;
36 using namespace genie::utils;
37 
38 //___________________________________________________________________________
40 KineGeneratorWithCache("genie::DFRKinematicsGenerator")
41 {
42 
43 }
44 //___________________________________________________________________________
46 KineGeneratorWithCache("genie::DFRKinematicsGenerator", config)
47 {
48 
49 }
50 //___________________________________________________________________________
52 {
53 
54 }
55 //___________________________________________________________________________
57 {
58  if(fGenerateUniformly) {
59  LOG("DFRKinematics", pNOTICE)
60  << "Generating kinematics uniformly over the allowed phase space";
61  }
62 
63  //-- Get the random number generators
65 
66  //-- Access cross section algorithm for running thread
68  const EventGeneratorI * evg = rtinfo->RunningThread();
69  fXSecModel = evg->CrossSectionAlg();
70 
71  //-- Get the interaction
72  Interaction * interaction = evrec->Summary();
73  interaction->SetBit(kISkipProcessChk);
74 
75  //-- Get neutrino energy and hit 'nucleon mass'
76  const InitialState & init_state = interaction->InitState();
77  double Ev = init_state.ProbeE(kRfHitNucRest);
78  double M = init_state.Tgt().HitNucP4().M(); // can be off m-shell
79 
80  //-- Get the physical W range
81  const KPhaseSpace & kps = interaction->PhaseSpace();
82 
83  Range1D_t xl = kps.Limits(kKVx);
84  Range1D_t yl = kps.Limits(kKVy);
85 
86  // Note that the t limits used here are NOT the same as what you get from KPhaseSpace::TLim().
87  // We use a larger range here because the limits from KPhaseSpace::TLim() depend on x and y,
88  // and using the rejection method in a region that's changing depending on some of the
89  // values guarantees that some regions will be oversampled. We want to avoid that.
90  Range1D_t tl;
91  tl.min = 0;
93 
94  LOG("DFRKinematics", pNOTICE) << "x: [" << xl.min << ", " << xl.max << "]";
95  LOG("DFRKinematics", pNOTICE) << "y: [" << yl.min << ", " << yl.max << "]";
96  LOG("DFRKinematics", pNOTICE) << "t: [" << tl.min << ", " << tl.max << "]";
97 
98 /*
99  Range1D_t W = kps.Limits(kKVW);
100  if(W.max <=0 || W.min>=W.max) {
101  LOG("DFRKinematics", pWARN) << "No available phase space";
102  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
103  genie::exceptions::EVGThreadException exception;
104  exception.SetReason("No available phase space");
105  exception.SwitchOnFastForward();
106  throw exception;
107  }
108 */
109 
110  assert(xl.min>0 && yl.min>0);
111 
112  //-- For the subsequent kinematic selection with the rejection method:
113  // Calculate the max differential cross section or retrieve it from the
114  // cache. Throw an exception and quit the evg thread if a non-positive
115  // value is found.
116  // If the kinematics are generated uniformly over the allowed phase
117  // space the max xsec is irrelevant
118  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
119 
120  //-- Try to select a valid (x,y,t) triplet using the rejection method
121 
122  double dx = xl.max - xl.min;
123  double dy = yl.max - yl.min;
124  double dt = tl.max - tl.min;
125  double gx=-1, gy=-1, gt=-1, gW=-1, gQ2=-1, xsec=-1;
126 
127  unsigned int iter = 0;
128  bool accept = false;
129  while(true) {
130  iter++;
131  if(iter > kRjMaxIterations) {
132  LOG("DFRKinematics", pWARN)
133  << " Couldn't select kinematics after " << iter << " iterations";
134  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
136  exception.SetReason("Couldn't select kinematics");
137  exception.SwitchOnFastForward();
138  throw exception;
139  }
140 
141  //-- random x,y,t
142  gx = xl.min + dx * rnd->RndKine().Rndm();
143  gy = yl.min + dy * rnd->RndKine().Rndm();
144  gt = tl.min + dt * rnd->RndKine().Rndm();
145 
146  interaction->KinePtr()->Setx(gx);
147  interaction->KinePtr()->Sety(gy);
148  interaction->KinePtr()->Sett(gt);
149 
150  LOG("DFRKinematics", pDEBUG)
151  << "Trying: x = " << gx << ", y = " << gy << ", t = " << gt;
152 
153  //-- compute the cross section for current kinematics
154  xsec = fXSecModel->XSec(interaction, kPSxytfE);
155 
156  //-- decide whether to accept the current kinematics
157  if(!fGenerateUniformly) {
158  this->AssertXSecLimits(interaction, xsec, xsec_max);
159  double n = xsec_max * rnd->RndKine().Rndm();
160  double J = 1;
161 
162 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
163  LOG("DFRKinematics", pDEBUG)
164  << "xsec= " << xsec << ", J= " << J << ", Rnd= " << n;
165 #endif
166  accept = (n < J*xsec);
167  }
168  else {
169  accept = (xsec>0);
170  }
171 
172  //-- If the generated kinematics are accepted, finish-up module's job
173  if(accept) {
174  // reset trust bits
175  interaction->ResetBit(kISkipProcessChk);
176  interaction->ResetBit(kISkipKinematicChk);
177 
178  // for uniform kinematics, compute an event weight as
179  // wght = (phase space volume)*(differential xsec)/(event total xsec)
180  if(fGenerateUniformly) {
181  double vol = kinematics::PhaseSpaceVolume(interaction,kPSxytfE);
182  double totxsec = evrec->XSec();
183  double wght = (vol/totxsec)*xsec;
184  LOG("DFRKinematics", pNOTICE) << "Kinematics wght = "<< wght;
185 
186  // apply computed weight to the current event weight
187  wght *= evrec->Weight();
188  LOG("DFRKinematics", pNOTICE) << "Current event wght = " << wght;
189  evrec->SetWeight(wght);
190  }
191 
192  // compute W,Q2 for selected x,y
193  kinematics::XYtoWQ2(Ev,M,gW,gQ2,gx,gy);
194 
195  LOG("DFRKinematics", pNOTICE)
196  << "Selected x,y => W = " << gW << ", Q2 = " << gQ2;
197 
198  // set the cross section for the selected kinematics
199  evrec->SetDiffXSec(xsec, kPSxytfE);
200 
201  // lock selected kinematics & clear running values
202  interaction->KinePtr()->SetW (gW, true);
203  interaction->KinePtr()->SetQ2(gQ2, true);
204  interaction->KinePtr()->Setx (gx, true);
205  interaction->KinePtr()->Sety (gy, true);
206  interaction->KinePtr()->Sett (gt, true);
207  interaction->KinePtr()->ClearRunningValues();
208  return;
209  }
210  } // iterations
211 }
212 //___________________________________________________________________________
214 {
215  Algorithm::Configure(config);
216  this->LoadConfig();
217 }
218 //____________________________________________________________________________
220 {
221  Algorithm::Configure(config);
222  this->LoadConfig();
223 }
224 //____________________________________________________________________________
226 {
227 // Reads its configuration data from its configuration Registry and loads them
228 // in private data members to avoid looking up at the Registry all the time.
229 
230  //-- Safety factor for the maximum differential cross section
231  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.25 ) ;
232 
233  //-- Minimum energy for which max xsec would be cached, forcing explicit
234  // calculation for lower eneries
235  GetParamDef( "Cache-MinEnergy", fEMin, 0.8 ) ;
236 
237  //-- Maximum allowed fractional cross section deviation from maxim cross
238  // section used in rejection method
239  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
240  assert(fMaxXSecDiffTolerance>=0);
241 
242  //-- Generate kinematics uniformly over allowed phase space and compute
243  // an event weight?
244  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
245 
246  GetParam( "DFR-Beta", fBeta ) ;
247 
248 }
249 //____________________________________________________________________________
251  const Interaction * interaction) const
252 {
253 // Computes the maximum differential cross section in the requested phase
254 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
255 // method and the value is cached at a circular cache branch for retrieval
256 // during subsequent event generation.
257 // The computed max differential cross section does not need to be the exact
258 // maximum. The number used in the rejection method will be scaled up by a
259 // safety factor. But this needs to be fast - do not use a very fine grid.
260 
261 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
262  LOG("DFRKinematics", pDEBUG)
263  << "Computing max xsec in allowed phase space";
264 #endif
265  double max_xsec = 0.0;
266 
267  const KPhaseSpace & kps = interaction->PhaseSpace();
268  Range1D_t xl = kps.XLim();
269  Range1D_t yl = kps.YLim();
270  Range1D_t Wl = kps.WLim();
271 
272  int Ny = 20;
273  int Nx = 20;
274  int Nt = 10;
275  double xmin = xl.min;
276  double xmax = xl.max;
277  double ymin = yl.min;
278  double ymax = yl.max;
279  double dx = (xmax-xmin)/(Nx-1);
280  double dy = (ymax-ymin)/(Ny-1);
281 
282 
283 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
284  LOG("DFRKinematics", pDEBUG)
285  << "Searching max. in x [" << xmin << ", " << xmax
286  << "], y [" << ymin << ", " << ymax << "], z [" << zmin << ", " << zmax << "]";
287 #endif
288  double xseclast_y = -1;
289  bool increasing_y;
290 
291  for(int i=0; i<Ny; i++) {
292  double gy = ymin + i*dy;
293  interaction->KinePtr()->Sety(gy);
294 
295 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
296  LOG("DFRKinematics", pDEBUG) << "y = " << gy;
297 #endif
298  for(int j=0; j<Nx; j++) {
299  double gx = xmin + j*dx;
300  interaction->KinePtr()->Setx(gx);
301  kinematics::UpdateWQ2FromXY(interaction);
302  Range1D_t Q2l = kps.Q2Lim_W();
303 
304  // the t range depends on x and y,
305  // and you get NaNs if you try to calculate t
306  // for an unphysical (x,y) combination
307  bool in_phys = math::IsWithinLimits(interaction->KinePtr()->W(), Wl);
308  in_phys = in_phys && math::IsWithinLimits(interaction->KinePtr()->Q2(), Q2l);
309  if (!in_phys)
310  continue;
311 
312  // t range depends on x and y
313  Range1D_t tl = kps.TLim();
314  double tmin = tl.min;
315  double tmax = tl.max;
316  double dt = (tmax-tmin)/(Nt-1);
317  for(int k=0; k<Nt; k++) {
318  double gt = tmin + k*dt;
319  interaction->KinePtr()->Sett(gt);
320 
321  double xsec = fXSecModel->XSec(interaction, kPSxytfE);
322 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
323  LOG("DFRKinematics", pINFO)
324  << "xsec(y=" << gy << ", x=" << gx << ", t=" << gt << ") = " << xsec;
325 #endif
326  // update maximum xsec
327  max_xsec = TMath::Max(xsec, max_xsec);
328  } // t
329  } // x
330  increasing_y = max_xsec-xseclast_y>=0;
331  xseclast_y = max_xsec;
332  if(!increasing_y) {
333 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
334  LOG("DFRKinematics", pDEBUG)
335  << "d2xsec/dxdy stopped increasing. Exiting y loop";
336 #endif
337  break;
338  }
339  }// y
340 
341  // Apply safety factor, since value retrieved from the cache might
342  // correspond to a slightly different energy
343  max_xsec *= 3;
344 
345 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
346  SLOG("DFRKinematics", pDEBUG) << interaction->AsString();
347  SLOG("DFRKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
348  SLOG("DFRKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel;
349 #endif
350 
351  return max_xsec;
352 }
353 //___________________________________________________________________________
double W(bool selected=false) const
Definition: Kinematics.cxx:157
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
virtual void SetWeight(double wght)
Definition: GHepRecord.h:130
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
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 ComputeMaxXSec(const Interaction *interaction) const
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
A simple [min,max] interval for doubles.
Definition: Range1.h:42
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:258
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec&gt;maxxsec
Defines the EventGeneratorI interface.
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
virtual double Weight(void) const
Definition: GHepRecord.h:124
Range1D_t YLim(void) const
y limits
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
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...
Definition: RandomGen.h:29
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
string AsString(void) const
static double GetTMaxDFR()
Definition: KPhaseSpace.cxx:59
double PhaseSpaceVolume(const Interaction *const i, KinePhaseSpace_t ps)
Definition: KineUtils.cxx:36
void Configure(const Registry &config)
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
void XYtoWQ2(double Ev, double M, double &W, double &Q2, double x, double y)
Definition: KineUtils.cxx:1155
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:291
#define pINFO
Definition: Messenger.h:62
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 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
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1290
static RunningThreadInfo * Instance(void)
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:117
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
virtual double XSec(void) const
Definition: GHepRecord.h:126
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 Q2(bool selected=false) const
Definition: Kinematics.cxx:125
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
Range1D_t TLim(void) const
t limits
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
const EventGeneratorI * RunningThread(void)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
void ProcessEventRecord(GHepRecord *event_rec) const
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
Keep info on the event generation thread currently on charge. This is used so that event generation m...
Range1D_t WLim(void) const
W limits.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Definition: GHepRecord.h:133
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63