GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DISKinematicsGenerator.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"
31 
32 using namespace genie;
33 using namespace genie::controls;
34 using namespace genie::utils;
35 
36 //___________________________________________________________________________
38 KineGeneratorWithCache("genie::DISKinematicsGenerator")
39 {
40 
41 }
42 //___________________________________________________________________________
44 KineGeneratorWithCache("genie::DISKinematicsGenerator", config)
45 {
46 
47 }
48 //___________________________________________________________________________
50 {
51 
52 }
53 //___________________________________________________________________________
55 {
56  if(fGenerateUniformly) {
57  LOG("DISKinematics", pNOTICE)
58  << "Generating kinematics uniformly over the allowed phase space";
59  }
60 
61  //-- Get the random number generators
63 
64  //-- Access cross section algorithm for running thread
66  const EventGeneratorI * evg = rtinfo->RunningThread();
67  fXSecModel = evg->CrossSectionAlg();
68 
69  //-- Get the interaction
70  Interaction * interaction = evrec->Summary();
71  interaction->SetBit(kISkipProcessChk);
72 
73  //-- Get neutrino energy and hit 'nucleon mass'
74  const InitialState & init_state = interaction->InitState();
75  double Ev = init_state.ProbeE(kRfHitNucRest);
76  double M = init_state.Tgt().HitNucP4().M(); // can be off m-shell
77 
78  //-- Get the physical W range
79  const KPhaseSpace & kps = interaction->PhaseSpace();
80  Range1D_t W = kps.Limits(kKVW);
81  if(W.max <=0 || W.min>=W.max) {
82  LOG("DISKinematics", pWARN) << "No available phase space";
83  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
85  exception.SetReason("No available phase space");
86  exception.SwitchOnFastForward();
87  throw exception;
88  }
89 
90  Range1D_t xl = kps.Limits(kKVx);
91  Range1D_t yl = kps.Limits(kKVy);
92 
93  LOG("DISKinematics", pNOTICE) << "x: [" << xl.min << ", " << xl.max << "]";
94  LOG("DISKinematics", pNOTICE) << "y: [" << yl.min << ", " << yl.max << "]";
95 
96  assert(xl.min>0 && yl.min>0);
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 (x,y) pair using the rejection method
107 
108  double dx = xl.max - xl.min;
109  double dy = yl.max - yl.min;
110  double gx=-1, gy=-1, gW=-1, gQ2=-1, xsec=-1;
111 
112  unsigned int iter = 0;
113  bool accept = false;
114  while(1) {
115  iter++;
116  if(iter > kRjMaxIterations) {
117  LOG("DISKinematics", pWARN)
118  << " Couldn't select kinematics after " << iter << " iterations";
119  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
121  exception.SetReason("Couldn't select kinematics");
122  exception.SwitchOnFastForward();
123  throw exception;
124  }
125 
126  //-- random x,y
127  gx = xl.min + dx * rnd->RndKine().Rndm();
128  gy = yl.min + dy * rnd->RndKine().Rndm();
129  interaction->KinePtr()->Setx(gx);
130  interaction->KinePtr()->Sety(gy);
131  kinematics::UpdateWQ2FromXY(interaction);
132 
133  LOG("DISKinematics", pNOTICE)
134  << "Trying: x = " << gx << ", y = " << gy
135  << " (W = " << interaction->KinePtr()->W() << ","
136  << " (Q2 = " << interaction->KinePtr()->Q2() << ")";
137 
138  //-- compute the cross section for current kinematics
139  xsec = fXSecModel->XSec(interaction, kPSxyfE);
140 
141  //-- decide whether to accept the current kinematics
142  if(!fGenerateUniformly) {
143  this->AssertXSecLimits(interaction, xsec, xsec_max);
144  double t = xsec_max * rnd->RndKine().Rndm();
145  double J = 1;
146 
147 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
148  LOG("DISKinematics", pDEBUG)
149  << "xsec= " << xsec << ", J= " << J << ", Rnd= " << t;
150 #endif
151  accept = (t < J*xsec);
152  }
153  else {
154  accept = (xsec>0);
155  }
156 
157  //-- If the generated kinematics are accepted, finish-up module's job
158  if(accept) {
159  LOG("DISKinematics", pNOTICE)
160  << "Selected: x = " << gx << ", y = " << gy
161  << " (W = " << interaction->KinePtr()->W() << ","
162  << " (Q2 = " << interaction->KinePtr()->Q2() << ")";
163 
164  // reset trust bits
165  interaction->ResetBit(kISkipProcessChk);
166  interaction->ResetBit(kISkipKinematicChk);
167 
168  // set the cross section for the selected kinematics
169  evrec->SetDiffXSec(xsec,kPSxyfE);
170 
171  // for uniform kinematics, compute an event weight as
172  // wght = (phase space volume)*(differential xsec)/(event total xsec)
173  if(fGenerateUniformly) {
174  double vol = kinematics::PhaseSpaceVolume(interaction,kPSxyfE);
175  double totxsec = evrec->XSec();
176  double wght = (vol/totxsec)*xsec;
177  LOG("DISKinematics", pNOTICE) << "Kinematics wght = "<< wght;
178 
179  // apply computed weight to the current event weight
180  wght *= evrec->Weight();
181  LOG("DISKinematics", pNOTICE) << "Current event wght = " << wght;
182  evrec->SetWeight(wght);
183  }
184 
185  // compute W,Q2 for selected x,y
186  //bool is_em = interaction->ProcInfo().IsEM();
187  kinematics::XYtoWQ2(Ev,M,gW,gQ2,gx,gy);
188 
189  LOG("DISKinematics", pNOTICE)
190  << "Selected x,y => W = " << gW << ", Q2 = " << gQ2;
191 
192  // lock selected kinematics & clear running values
193  interaction->KinePtr()->SetW (gW, true);
194  interaction->KinePtr()->SetQ2(gQ2, true);
195  interaction->KinePtr()->Setx (gx, true);
196  interaction->KinePtr()->Sety (gy, true);
197  interaction->KinePtr()->ClearRunningValues();
198  return;
199  }
200  } // iterations
201 }
202 //___________________________________________________________________________
204 {
205  Algorithm::Configure(config);
206  this->LoadConfig();
207 }
208 //____________________________________________________________________________
210 {
211  Algorithm::Configure(config);
212  this->LoadConfig();
213 }
214 //____________________________________________________________________________
216 {
217 // Reads its configuration data from its configuration Registry and loads them
218 // in private data members to avoid looking up at the Registry all the time.
219 
220  //-- Safety factor for the maximum differential cross section
221  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.25 ) ;
222 
223  //-- Minimum energy for which max xsec would be cached, forcing explicit
224  // calculation for lower eneries
225  GetParamDef( "Cache-MinEnergy", fEMin, 0.8 ) ;
226 
227  //-- Maximum allowed fractional cross section deviation from maxim cross
228  // section used in rejection method
229  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
230  assert(fMaxXSecDiffTolerance>=0);
231 
232  //-- Generate kinematics uniformly over allowed phase space and compute
233  // an event weight?
234  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
235 
236 }
237 //____________________________________________________________________________
239  const Interaction * interaction) const
240 {
241 // Computes the maximum differential cross section in the requested phase
242 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
243 // method and the value is cached at a circular cache branch for retrieval
244 // during subsequent event generation.
245 // The computed max differential cross section does not need to be the exact
246 // maximum. The number used in the rejection method will be scaled up by a
247 // safety factor. But this needs to be fast - do not use a very fine grid.
248 
249 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
250  LOG("DISKinematics", pDEBUG)<< "Computing max xsec in allowed phase space";
251 #endif
252  double max_xsec = 0.0;
253 
254  const InitialState & init_state = interaction->InitState();
255  //double Ev = init_state.ProbeE(kRfHitNucRest);
256  //const ProcessInfo & proc = interaction->ProcInfo();
257  const Target & tgt = init_state.Tgt();
258 
259  int Ny = 20;
260  int Nx = 40;
261  int Nxb = 3;
262 
263  double xpeak = .2;
264  double xwindow = .2;
265  double ypeak = .5;
266  double ywindow = .5;
267 
268  if(tgt.HitQrkIsSet()) {
269  if(tgt.HitSeaQrk()) {
270  xpeak = .1;
271  xwindow = .1;
272  ypeak = .7;
273  ywindow = .3;
274  } else {
275  xpeak = .2;
276  xwindow = .2;
277  ypeak = .7;
278  ywindow = .3;
279  }
280  }
281 
282  const KPhaseSpace & kps = interaction->PhaseSpace();
283  Range1D_t xl = kps.Limits(kKVx);
284  Range1D_t yl = kps.Limits(kKVy);
285 
286  double xmin = TMath::Max(xpeak-xwindow, TMath::Max(xl.min, 5E-3));
287  double xmax = TMath::Min(xpeak+xwindow, xl.max);
288  double ymin = TMath::Max(ypeak-ywindow, TMath::Max(yl.min, 2E-3));
289  double ymax = TMath::Min(ypeak+ywindow, yl.max);
290  double dx = (xmax-xmin)/(Nx-1);
291  double dy = (ymax-ymin)/(Ny-1);
292 
293 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
294  LOG("DISKinematics", pDEBUG)
295  << "Searching max. in x [" << xmin << ", " << xmax << "], y [" << ymin << ", " << ymax << "]";
296 #endif
297  double xseclast_y = -1;
298  bool increasing_y;
299 
300  for(int i=0; i<Ny; i++) {
301  double gy = ymin + i*dy;
302  //double gy = TMath::Power(10., logymin + i*dlogy);
303  interaction->KinePtr()->Sety(gy);
304 
305 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
306  LOG("DISKinematics", pDEBUG) << "y = " << gy;
307 #endif
308  double xseclast_x = -1;
309  bool increasing_x;
310 
311  for(int j=0; j<Nx; j++) {
312  double gx = xmin + j*dx;
313  //double gx = TMath::Power(10., logxmin + j*dlogx);
314  interaction->KinePtr()->Setx(gx);
315  kinematics::UpdateWQ2FromXY(interaction);
316 
317  double xsec = fXSecModel->XSec(interaction, kPSxyfE);
318 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
319  LOG("DISKinematics", pINFO)
320  << "xsec(y=" << gy << ", x=" << gx << ") = " << xsec;
321 #endif
322  // update maximum xsec
323  max_xsec = TMath::Max(xsec, max_xsec);
324 
325  increasing_x = xsec-xseclast_x>=0;
326  xseclast_x = xsec;
327 
328  // once the cross section stops increasing, I reduce the step size and
329  // step backwards a little bit to handle cases that the max cross section
330  // is grossly underestimated (very peaky distribution & large step)
331  if(!increasing_x) {
332 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
333  LOG("DISKinematics", pDEBUG)
334  << "d2xsec/dxdy|x stopped increasing. Stepping back & exiting x loop";
335 #endif
336  //double dlogxn = dlogx/(Nxb+1);
337  double dxn = dx/(Nxb+1);
338  for(int ik=0; ik<Nxb; ik++) {
339  //gx = TMath::Exp(TMath::Log(gx) - dlogxn);
340  gx = gx - dxn;
341  interaction->KinePtr()->Setx(gx);
342  kinematics::UpdateWQ2FromXY(interaction);
343  xsec = fXSecModel->XSec(interaction, kPSxyfE);
344 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
345  LOG("DISKinematics", pINFO)
346  << "xsec(y=" << gy << ", x=" << gx << ") = " << xsec;
347 #endif
348  }
349  break;
350  } // stepping back within last bin
351  } // x
352  increasing_y = max_xsec-xseclast_y>=0;
353  xseclast_y = max_xsec;
354  if(!increasing_y) {
355 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
356  LOG("DISKinematics", pDEBUG)
357  << "d2xsec/dxdy stopped increasing. Exiting y loop";
358 #endif
359  break;
360  }
361  }// y
362 
363  // Apply safety factor, since value retrieved from the cache might
364  // correspond to a slightly different energy
365  // max_xsec *= fSafetyFactor;
366  //max_xsec *= ( (Ev<3.0) ? 2.5 : fSafetyFactor);
367  max_xsec *= 3;
368 
369 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
370  SLOG("DISKinematics", pDEBUG) << interaction->AsString();
371  SLOG("DISKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
372  SLOG("DISKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel;
373 #endif
374 
375  return max_xsec;
376 }
377 //___________________________________________________________________________
double ComputeMaxXSec(const Interaction *interaction) const
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 HitSeaQrk(void) const
Definition: Target.cxx:299
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 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
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
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:29
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:36
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
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
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
#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
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
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
void ProcessEventRecord(GHepRecord *event_rec) const
virtual double XSec(void) const
Definition: GHepRecord.h:126
double gQ2
Definition: gtestDISSF.cxx:55
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
const Target & Tgt(void) const
Definition: InitialState.h:66
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
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...
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