DISKinematicsGenerator.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
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Mar 03, 2009 - CA
14  Moved into the new DIS package from its previous location (EVGModules).
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 #include <cfloat>
22 
23 #include <TMath.h>
24 
41 
42 using namespace genie;
43 using namespace genie::controls;
44 using namespace genie::utils;
45 
46 //___________________________________________________________________________
48 KineGeneratorWithCache("genie::DISKinematicsGenerator")
49 {
50 
51 }
52 //___________________________________________________________________________
54 KineGeneratorWithCache("genie::DISKinematicsGenerator", config)
55 {
56 
57 }
58 //___________________________________________________________________________
60 {
61 
62 }
63 //___________________________________________________________________________
65 {
66  if(fGenerateUniformly) {
67  LOG("DISKinematics", pNOTICE)
68  << "Generating kinematics uniformly over the allowed phase space";
69  }
70 
71  //-- Get the random number generators
73 
74  //-- Access cross section algorithm for running thread
76  const EventGeneratorI * evg = rtinfo->RunningThread();
77  fXSecModel = evg->CrossSectionAlg();
78 
79  //-- Get the interaction
80  Interaction * interaction = evrec->Summary();
81  interaction->SetBit(kISkipProcessChk);
82 
83  //-- Get neutrino energy and hit 'nucleon mass'
84  const InitialState & init_state = interaction->InitState();
85  double Ev = init_state.ProbeE(kRfHitNucRest);
86  double M = init_state.Tgt().HitNucP4().M(); // can be off m-shell
87 
88  //-- Get the physical W range
89  const KPhaseSpace & kps = interaction->PhaseSpace();
90  Range1D_t W = kps.Limits(kKVW);
91  if(W.max <=0 || W.min>=W.max) {
92  LOG("DISKinematics", pWARN) << "No available phase space";
93  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
95  exception.SetReason("No available phase space");
96  exception.SwitchOnFastForward();
97  throw exception;
98  }
99 
100  Range1D_t xl = kps.Limits(kKVx);
101  Range1D_t yl = kps.Limits(kKVy);
102 
103  LOG("DISKinematics", pNOTICE) << "x: [" << xl.min << ", " << xl.max << "]";
104  LOG("DISKinematics", pNOTICE) << "y: [" << yl.min << ", " << yl.max << "]";
105 
106  assert(xl.min>0 && yl.min>0);
107 
108  //-- For the subsequent kinematic selection with the rejection method:
109  // Calculate the max differential cross section or retrieve it from the
110  // cache. Throw an exception and quit the evg thread if a non-positive
111  // value is found.
112  // If the kinematics are generated uniformly over the allowed phase
113  // space the max xsec is irrelevant
114  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
115 
116  //-- Try to select a valid (x,y) pair using the rejection method
117 
118  double dx = xl.max - xl.min;
119  double dy = yl.max - yl.min;
120  double gx=-1, gy=-1, gW=-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("DISKinematics", 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  //-- random x,y
137  gx = xl.min + dx * rnd->RndKine().Rndm();
138  gy = yl.min + dy * rnd->RndKine().Rndm();
139  interaction->KinePtr()->Setx(gx);
140  interaction->KinePtr()->Sety(gy);
141  kinematics::UpdateWQ2FromXY(interaction);
142 
143  LOG("DISKinematics", pNOTICE)
144  << "Trying: x = " << gx << ", y = " << gy
145  << " (W = " << interaction->KinePtr()->W() << ","
146  << " (Q2 = " << interaction->KinePtr()->Q2() << ")";
147 
148  //-- compute the cross section for current kinematics
149  xsec = fXSecModel->XSec(interaction, kPSxyfE);
150 
151  //-- decide whether to accept the current kinematics
152  if(!fGenerateUniformly) {
153  this->AssertXSecLimits(interaction, xsec, xsec_max);
154  double t = xsec_max * rnd->RndKine().Rndm();
155  double J = 1;
156 
157 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
158  LOG("DISKinematics", pDEBUG)
159  << "xsec= " << xsec << ", J= " << J << ", Rnd= " << t;
160 #endif
161  accept = (t < J*xsec);
162  }
163  else {
164  accept = (xsec>0);
165  }
166 
167  //-- If the generated kinematics are accepted, finish-up module's job
168  if(accept) {
169  LOG("DISKinematics", pNOTICE)
170  << "Selected: x = " << gx << ", y = " << gy
171  << " (W = " << interaction->KinePtr()->W() << ","
172  << " (Q2 = " << interaction->KinePtr()->Q2() << ")";
173 
174  // reset trust bits
175  interaction->ResetBit(kISkipProcessChk);
176  interaction->ResetBit(kISkipKinematicChk);
177 
178  // set the cross section for the selected kinematics
179  evrec->SetDiffXSec(xsec,kPSxyfE);
180 
181  // for uniform kinematics, compute an event weight as
182  // wght = (phase space volume)*(differential xsec)/(event total xsec)
183  if(fGenerateUniformly) {
184  double vol = kinematics::PhaseSpaceVolume(interaction,kPSxyfE);
185  double totxsec = evrec->XSec();
186  double wght = (vol/totxsec)*xsec;
187  LOG("DISKinematics", pNOTICE) << "Kinematics wght = "<< wght;
188 
189  // apply computed weight to the current event weight
190  wght *= evrec->Weight();
191  LOG("DISKinematics", pNOTICE) << "Current event wght = " << wght;
192  evrec->SetWeight(wght);
193  }
194 
195  // compute W,Q2 for selected x,y
196  bool is_em = interaction->ProcInfo().IsEM();
197  kinematics::XYtoWQ2(Ev,M,gW,gQ2,gx,gy);
198 
199  LOG("DISKinematics", pNOTICE)
200  << "Selected x,y => W = " << gW << ", Q2 = " << gQ2;
201 
202  // lock selected kinematics & clear running values
203  interaction->KinePtr()->SetW (gW, true);
204  interaction->KinePtr()->SetQ2(gQ2, true);
205  interaction->KinePtr()->Setx (gx, true);
206  interaction->KinePtr()->Sety (gy, 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. ) ;
241 
242  //-- Generate kinematics uniformly over allowed phase space and compute
243  // an event weight?
244  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
245 
246 }
247 //____________________________________________________________________________
249  const Interaction * interaction) const
250 {
251 // Computes the maximum differential cross section in the requested phase
252 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
253 // method and the value is cached at a circular cache branch for retrieval
254 // during subsequent event generation.
255 // The computed max differential cross section does not need to be the exact
256 // maximum. The number used in the rejection method will be scaled up by a
257 // safety factor. But this needs to be fast - do not use a very fine grid.
258 
259 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
260  LOG("DISKinematics", pDEBUG)<< "Computing max xsec in allowed phase space";
261 #endif
262  double max_xsec = 0.0;
263 
264  const InitialState & init_state = interaction->InitState();
265  //double Ev = init_state.ProbeE(kRfHitNucRest);
266  //const ProcessInfo & proc = interaction->ProcInfo();
267  const Target & tgt = init_state.Tgt();
268 
269  int Ny = 20;
270  int Nx = 40;
271  int Nxb = 3;
272 
273  double xpeak = .2;
274  double xwindow = .2;
275  double ypeak = .5;
276  double ywindow = .5;
277 
278  if(tgt.HitQrkIsSet()) {
279  if(tgt.HitSeaQrk()) {
280  xpeak = .1;
281  xwindow = .1;
282  ypeak = .7;
283  ywindow = .3;
284  } else {
285  xpeak = .2;
286  xwindow = .2;
287  ypeak = .7;
288  ywindow = .3;
289  }
290  }
291 
292  const KPhaseSpace & kps = interaction->PhaseSpace();
293  Range1D_t xl = kps.Limits(kKVx);
294  Range1D_t yl = kps.Limits(kKVy);
295 
296  double xmin = TMath::Max(xpeak-xwindow, TMath::Max(xl.min, 5E-3));
297  double xmax = TMath::Min(xpeak+xwindow, xl.max);
298  double ymin = TMath::Max(ypeak-ywindow, TMath::Max(yl.min, 2E-3));
299  double ymax = TMath::Min(ypeak+ywindow, yl.max);
300  double dx = (xmax-xmin)/(Nx-1);
301  double dy = (ymax-ymin)/(Ny-1);
302 
303 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
304  LOG("DISKinematics", pDEBUG)
305  << "Searching max. in x [" << xmin << ", " << xmax << "], y [" << ymin << ", " << ymax << "]";
306 #endif
307  double xseclast_y = -1;
308  bool increasing_y;
309 
310  for(int i=0; i<Ny; i++) {
311  double gy = ymin + i*dy;
312  //double gy = TMath::Power(10., logymin + i*dlogy);
313  interaction->KinePtr()->Sety(gy);
314 
315 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
316  LOG("DISKinematics", pDEBUG) << "y = " << gy;
317 #endif
318  double xseclast_x = -1;
319  bool increasing_x;
320 
321  for(int j=0; j<Nx; j++) {
322  double gx = xmin + j*dx;
323  //double gx = TMath::Power(10., logxmin + j*dlogx);
324  interaction->KinePtr()->Setx(gx);
325  kinematics::UpdateWQ2FromXY(interaction);
326 
327  double xsec = fXSecModel->XSec(interaction, kPSxyfE);
328 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
329  LOG("DISKinematics", pINFO)
330  << "xsec(y=" << gy << ", x=" << gx << ") = " << xsec;
331 #endif
332  // update maximum xsec
333  max_xsec = TMath::Max(xsec, max_xsec);
334 
335  increasing_x = xsec-xseclast_x>=0;
336  xseclast_x = xsec;
337 
338  // once the cross section stops increasing, I reduce the step size and
339  // step backwards a little bit to handle cases that the max cross section
340  // is grossly underestimated (very peaky distribution & large step)
341  if(!increasing_x) {
342 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
343  LOG("DISKinematics", pDEBUG)
344  << "d2xsec/dxdy|x stopped increasing. Stepping back & exiting x loop";
345 #endif
346  //double dlogxn = dlogx/(Nxb+1);
347  double dxn = dx/(Nxb+1);
348  for(int ik=0; ik<Nxb; ik++) {
349  //gx = TMath::Exp(TMath::Log(gx) - dlogxn);
350  gx = gx - dxn;
351  interaction->KinePtr()->Setx(gx);
352  kinematics::UpdateWQ2FromXY(interaction);
353  xsec = fXSecModel->XSec(interaction, kPSxyfE);
354 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
355  LOG("DISKinematics", pINFO)
356  << "xsec(y=" << gy << ", x=" << gx << ") = " << xsec;
357 #endif
358  }
359  break;
360  } // stepping back within last bin
361  } // x
362  increasing_y = max_xsec-xseclast_y>=0;
363  xseclast_y = max_xsec;
364  if(!increasing_y) {
365 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
366  LOG("DISKinematics", pDEBUG)
367  << "d2xsec/dxdy stopped increasing. Exiting y loop";
368 #endif
369  break;
370  }
371  }// y
372 
373  // Apply safety factor, since value retrieved from the cache might
374  // correspond to a slightly different energy
375  // max_xsec *= fSafetyFactor;
376  //max_xsec *= ( (Ev<3.0) ? 2.5 : fSafetyFactor);
377  max_xsec *= 3;
378 
379 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
380  SLOG("DISKinematics", pDEBUG) << interaction->AsString();
381  SLOG("DISKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
382  SLOG("DISKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel;
383 #endif
384 
385  return max_xsec;
386 }
387 //___________________________________________________________________________
388 
double ComputeMaxXSec(const Interaction *interaction) const
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
virtual void SetWeight(double wght)
Definition: GHepRecord.h:131
bool HitSeaQrk(void) const
Definition: Target.cxx:316
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 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
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
Definition: config.py:1
void Configure(const Registry &config)
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
double PhaseSpaceVolume(const Interaction *const i, KinePhaseSpace_t ps)
Definition: KineUtils.cxx:34
Double_t ymax
Definition: plot.C:25
Summary information for an interaction.
Definition: Interaction.h:56
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
Float_t E
Definition: plot.C:20
void XYtoWQ2(double Ev, double M, double &W, double &Q2, double x, double y)
Definition: KineUtils.cxx:1069
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
Kinematical phase space.
Definition: KPhaseSpace.h:34
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
bool IsEM(void) const
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
bool HitQrkIsSet(void) const
Definition: Target.cxx:309
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.
void ProcessEventRecord(GHepRecord *event_rec) const
assert(nhit_max >=nhit_nbins)
virtual double XSec(void) const
Definition: GHepRecord.h:127
double gQ2
Definition: gtestDISSF.cxx:55
const InitialState & InitState(void) const
Definition: Interaction.h:69
Double_t ymin
Definition: plot.C:24
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
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
const Target & Tgt(void) const
Definition: InitialState.h:67
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
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
#define W(x)
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.
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