GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NuEKinematicsGenerator.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 
23 
24 using namespace genie;
25 using namespace genie::controls;
26 using namespace genie::utils;
27 
28 //___________________________________________________________________________
30 KineGeneratorWithCache("genie::NuEKinematicsGenerator")
31 {
32 
33 }
34 //___________________________________________________________________________
36 KineGeneratorWithCache("genie::NuEKinematicsGenerator", config)
37 {
38 
39 }
40 //___________________________________________________________________________
42 {
43 
44 }
45 //___________________________________________________________________________
47 {
48  if(fGenerateUniformly) {
49  LOG("NuEKinematics", pNOTICE)
50  << "Generating kinematics uniformly over the allowed phase space";
51  }
52 
53  //-- Get the random number generators
55 
56  //-- Access cross section algorithm for running thread
58  const EventGeneratorI * evg = rtinfo->RunningThread();
59  fXSecModel = evg->CrossSectionAlg();
60 
61  //-- For the subsequent kinematic selection with the rejection method:
62  // Calculate the max differential cross section or retrieve it from the
63  // cache. Throw an exception and quit the evg thread if a non-positive
64  // value is found.
65  // If the kinematics are generated uniformly over the allowed phase
66  // space the max xsec is irrelevant
67  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
68 
69  //-- y range
70  const KPhaseSpace & kps = evrec->Summary()->PhaseSpace();
71  Range1D_t yl = kps.Limits(kKVy);
72  double ymin = yl.min;
73  double ymax = yl.max;
74  double dy = ymax-ymin;
75 
76  double xsec = -1;
77  Interaction * interaction = evrec->Summary();
78 
79  //-- Try to select a valid inelastisity y
80  unsigned int iter = 0;
81  bool accept = false;
82  while(1) {
83  iter++;
84  if(iter > kRjMaxIterations) {
85  LOG("NuEKinematics", pWARN)
86  << "*** Could not select a valid y after "
87  << iter << " iterations";
88  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
90  exception.SetReason("Couldn't select kinematics");
91  exception.SwitchOnFastForward();
92  throw exception;
93  }
94 
95  double y = ymin + dy * rnd->RndKine().Rndm();
96  interaction->KinePtr()->Sety(y);
97 
98  LOG("NuEKinematics", pINFO) << "Trying: y = " << y;
99 
100  //-- computing cross section for the current kinematics
101  xsec = fXSecModel->XSec(interaction, kPSyfE);
102 
103  //-- decide whether to accept the current kinematics
104  if(!fGenerateUniformly) {
105  this->AssertXSecLimits(interaction, xsec, xsec_max);
106 
107  double t = xsec_max * rnd->RndKine().Rndm();
108  LOG("NuEKinematics", pDEBUG) << "xsec= "<< xsec<< ", J= 1, Rnd= "<< t;
109 
110  accept = (t<xsec);
111  } else {
112  accept = (xsec>0);
113  }
114 
115  //-- If the generated kinematics are accepted, finish-up module's job
116  if(accept) {
117  LOG("NuEKinematics", pINFO) << "Selected: y = " << y;
118 
119  // set the cross section for the selected kinematics
120  evrec->SetDiffXSec(xsec,kPSyfE);
121 
122  // for uniform kinematics, compute an event weight as
123  // wght = (phase space volume)*(differential xsec)/(event total xsec)
124  if(fGenerateUniformly) {
125  double vol = kinematics::PhaseSpaceVolume(interaction,kPSyfE);
126  double totxsec = evrec->XSec();
127  double wght = (vol/totxsec)*xsec;
128  LOG("NuEKinematics", pNOTICE) << "Kinematics wght = "<< wght;
129 
130  // apply computed weight to the current event weight
131  wght *= evrec->Weight();
132  LOG("NuEKinematics", pNOTICE) << "Current event wght = " << wght;
133  evrec->SetWeight(wght);
134  }
135 
136  // lock selected kinematics & clear running values
137  interaction->KinePtr()->Sety(y, true);
138  interaction->KinePtr()->ClearRunningValues();
139 
140  return;
141  }
142  }// iterations
143 }
144 //___________________________________________________________________________
146  const Interaction * interaction) const
147 {
148 // Computes the maximum differential cross section in the requested phase
149 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
150 // method and the value is cached at a circular cache branch for retrieval
151 // during subsequent event generation.
152 // The computed max differential cross section does not need to be the exact
153 // maximum. The number used in the rejection method will be scaled up by a
154 // safety factor. But it needs to be fast - do not use a very small y step.
155 
156  const int N = 40;
157 //const int Nb = 6;
158 
159  const KPhaseSpace & kps = interaction->PhaseSpace();
160  Range1D_t yl = kps.Limits(kKVy);
161  const double ymin = yl.min;
162  const double ymax = yl.max;
163 
164  double max_xsec = -1.0;
165 
166  double dy = (ymax-ymin)/(N-1);
167 //double xseclast = -1;
168 //bool increasing;
169 
170  for(int i=0; i<N; i++) {
171  double y = ymin + i * dy;
172  interaction->KinePtr()->Sety(y);
173  double xsec = fXSecModel->XSec(interaction, kPSyfE);
174 
175  SLOG("NuEKinematics", pDEBUG) << "xsec(y = " << y << ") = " << xsec;
176  max_xsec = TMath::Max(xsec, max_xsec);
177 /*
178  increasing = xsec-xseclast>=0;
179  xseclast = xsec;
180 
181  // once the cross section stops increasing, I reduce the step size and
182  // step backwards a little bit to handle cases that the max cross section
183  // is grossly underestimated (very peaky distribution & large step)
184  if(!increasing) {
185  dy/=(Nb+1);
186  for(int ib=0; ib<Nb; ib++) {
187  y = y-dy;
188  if(y<ymin) break;
189  interaction->KinePtr()->Sety(y);
190  xsec = fXSecModel->XSec(interaction, kPSyfE);
191  SLOG("NuEKinematics", pDEBUG) << "xsec(y = " << y << ") = " << xsec;
192  max_xsec = TMath::Max(xsec, max_xsec);
193  }
194  break;
195  }
196 */
197  }//y
198 
199  // Apply safety factor, since value retrieved from the cache might
200  // correspond to a slightly different energy.
201  max_xsec *= fSafetyFactor;
202 
203  SLOG("NuEKinematics", pDEBUG) << interaction->AsString();
204  SLOG("NuEKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
205  SLOG("NuEKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel;
206 
207  return max_xsec;
208 }
209 //___________________________________________________________________________
210 double NuEKinematicsGenerator::Energy(const Interaction * interaction) const
211 {
212 // Override the base class Energy() method to cache the max xsec for the
213 // neutrino energy in the LAB rather than in the hit nucleon rest frame.
214 
215  const InitialState & init_state = interaction->InitState();
216  double E = init_state.ProbeE(kRfLab);
217  return E;
218 }
219 //___________________________________________________________________________
221 {
222  Algorithm::Configure(config);
223  this->LoadConfig();
224 }
225 //____________________________________________________________________________
227 {
228  Algorithm::Configure(config);
229  this->LoadConfig();
230 }
231 //____________________________________________________________________________
233 {
234  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 2.00 ) ;
235  GetParamDef( "Cache-MinEnergy", fEMin, 1.00 ) ;
236 
237  GetParamDef("MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 0. ) ;
238  assert(fMaxXSecDiffTolerance>=0);
239 
240  //-- Generate kinematics uniformly over allowed phase space and compute
241  // an event weight?
242  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
243 
244 }
245 //____________________________________________________________________________
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?
void Configure(const Registry &config)
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
double ComputeMaxXSec(const Interaction *in) const
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.
void ProcessEventRecord(GHepRecord *event_rec) const
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
virtual double Weight(void) const
Definition: GHepRecord.h:124
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
Summary information for an interaction.
Definition: Interaction.h:56
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
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
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
static RunningThreadInfo * Instance(void)
double max
Definition: Range1.h:53
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
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
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...
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Definition: GHepRecord.h:133
double Energy(const Interaction *in) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63