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