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