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