COHElKinematicsGenerator.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 
14 */
15 //____________________________________________________________________________
16 
17 #include <cstdlib>
18 
19 #include <TMath.h>
20 
35 
36 using namespace genie;
37 using namespace genie::constants;
38 using namespace genie::controls;
39 using namespace genie::utils;
40 
41 //___________________________________________________________________________
43 KineGeneratorWithCache("genie::COHElKinematicsGenerator")
44 {
45 
46 }
47 //___________________________________________________________________________
49 KineGeneratorWithCache("genie::COHElKinematicsGenerator", config)
50 {
51 
52 }
53 //___________________________________________________________________________
55 {
56 
57 }
58 //___________________________________________________________________________
60 {
61  if(fGenerateUniformly) {
62  LOG("COHElKinematics", pNOTICE)
63  << "Generating kinematics uniformly over the allowed phase space";
64  }
65 
66  Interaction * interaction = evrec->Summary();
67  interaction->SetBit(kISkipProcessChk);
68  interaction->SetBit(kISkipKinematicChk);
69 
70  //-- Get the random number generators
72 
73  //-- Access cross section algorithm for running thread
75  const EventGeneratorI * evg = rtinfo->RunningThread();
76  fXSecModel = evg->CrossSectionAlg();
77 
78  //-- For the subsequent kinematic selection with the rejection method:
79  // Calculate the max differential cross section or retrieve it from the
80  // cache. Throw an exception and quit the evg thread if a non-positive
81  // value is found.
82  // If the kinematics are generated uniformly over the allowed phase
83  // space the max xsec is irrelevant
84  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
85 
86  //-- Get the kinematical limits for the generated x,y
87  const KPhaseSpace & kps = interaction->PhaseSpace();
88  Range1D_t y = kps.YLim();
89  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
90  const double ymin = y.min + kASmallNum;
91  const double ymax = y.max - kASmallNum;
92  const double dy = ymax - ymin;
93 
94  //------ Try to select a valid y
95  unsigned int iter = 0;
96  bool accept=false;
97  double xsec=-1, gy=-1;
98 
99  while(1) {
100  iter++;
101  if(iter > kRjMaxIterations) {
102  LOG("COHElKinematics", pWARN)
103  << "*** Could not select a valid y after " << iter << " iterations";
104  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
106  exception.SetReason("Couldn't select kinematics");
107  exception.SwitchOnFastForward();
108  throw exception;
109  }
110 
111  gy = ymin + dy * rnd->RndKine().Rndm();
112  LOG("COHElKinematics", pINFO) << "Trying: y = " << gy;
113 
114  interaction->KinePtr()->Sety(gy);
115 
116  // computing cross section for the current kinematics
117  xsec = fXSecModel->XSec(interaction, kPSyfE);
118 
119  //-- decide whether to accept the current kinematics
120  if(!fGenerateUniformly) {
121  double t = xsec_max * rnd->RndKine().Rndm();
122  this->AssertXSecLimits(interaction, xsec, xsec_max);
123  LOG("COHElKinematics", pINFO) << "xsec= " << xsec << ", J= 1, Rnd= " << t;
124  accept = (t<xsec);
125  }
126  else {
127  accept = (xsec>0);
128  }
129 
130  //-- If the generated kinematics are accepted, finish-up module's job
131  if(accept) {
132  LOG("COHElKinematics", pNOTICE) << "Selected: y = " << gy;
133 
134  // for uniform kinematics, compute an event weight as
135  // wght = (phase space volume)*(differential xsec)/(event total xsec)
136  if(fGenerateUniformly) {
137  double vol = y.max-y.min;
138  double totxsec = evrec->XSec();
139  double wght = (vol/totxsec)*xsec;
140  LOG("COHElKinematics", pNOTICE) << "Kinematics wght = "<< wght;
141 
142  // apply computed weight to the current event weight
143  wght *= evrec->Weight();
144  LOG("COHElKinematics", pNOTICE) << "Current event wght = " << wght;
145  evrec->SetWeight(wght);
146  }
147 
148  // reset bits
149  interaction->ResetBit(kISkipProcessChk);
150  interaction->ResetBit(kISkipKinematicChk);
151 
152  // lock selected kinematics & clear running values
153  double Ev = interaction->InitState().ProbeE(kRfLab);
154  double x = 1;
155  double Q2 = 2*kNucleonMass*x*gy*Ev;
156  interaction->KinePtr()->Setx(x, true);
157  interaction->KinePtr()->Sety(gy, true);
158  interaction->KinePtr()->Sett(0, true);
159  interaction->KinePtr()->SetW(0, true);
160  interaction->KinePtr()->SetQ2(Q2, true);
161  interaction->KinePtr()->ClearRunningValues();
162 
163  // set the cross section for the selected kinematics
164  // TODO: this is kPSxyfE on the dev branch, but I think kPSyfE is right
165  evrec->SetDiffXSec(xsec,kPSyfE);
166 
167  return;
168  }
169  }// iterations
170 }
171 //___________________________________________________________________________
173 {
174 // Computes the maximum differential cross section in the requested phase
175 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
176 // method and the value is cached at a circular cache branch for retrieval
177 // during subsequent event generation.
178 
179  SLOG("COHElKinematics", pDEBUG)
180  << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
181 
182  double max_xsec = 0.;
183  const int N = 50;
184 
185  const KPhaseSpace & kps = in->PhaseSpace();
186  Range1D_t yr = kps.YLim();
187 
188  const double logymin = TMath::Log10(yr.min);
189  const double logymax = TMath::Log10(yr.max);
190  const double dlogy = (logymax - logymin) /(N-1);
191 
192  for(int i=0; i<N; i++) {
193  double y = TMath::Power(10, logymin+i*dlogy);
194  in->KinePtr()->Sety(y);
195 
196  double xsec = fXSecModel->XSec(in, kPSyfE);
197  LOG("COHElKinematics", pDEBUG) << "xsec(y= " << y << ") = " << xsec;
198  max_xsec = TMath::Max(max_xsec, xsec);
199  }//y
200 
201  // Apply safety factor, since value retrieved from the cache might
202  // correspond to a slightly different energy.
203  max_xsec *= fSafetyFactor;
204 
205  SLOG("COHElKinematics", pDEBUG) << in->AsString();
206  SLOG("COHElKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
207  SLOG("COHElKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
208 
209  return max_xsec;
210 }
211 //___________________________________________________________________________
212 double COHElKinematicsGenerator::Energy(const Interaction * interaction) const
213 {
214 // Override the base class Energy() method to cache the max xsec for the
215 // neutrino energy in the LAB rather than in the hit nucleon rest frame.
216 
217  const InitialState & init_state = interaction->InitState();
218  double E = init_state.ProbeE(kRfLab);
219  return E;
220 }
221 //___________________________________________________________________________
223 {
224  Algorithm::Configure(config);
225  this->LoadConfig();
226 }
227 //____________________________________________________________________________
229 {
230  Algorithm::Configure(config);
231  this->LoadConfig();
232 }
233 //____________________________________________________________________________
235 {
236  //-- max xsec safety factor (for rejection method) and min cached energy
237  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.6 ) ;
238  GetParamDef( "Cache-MinEnergy", fEMin, -1.0 ) ;
239 
240  //-- Generate kinematics uniformly over allowed phase space and compute
241  // an event weight?
242  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
243 
244  //-- Maximum allowed fractional cross section deviation from maxim cross
245  // section used in rejection method
246  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
248 }
249 //____________________________________________________________________________
250 
virtual double MaxXSec(GHepRecord *evrec) const
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
virtual void SetWeight(double wght)
Definition: GHepRecord.h:131
void ProcessEventRecord(GHepRecord *event_rec) const
bool fGenerateUniformly
uniform over allowed phase space + event weight?
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
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
Range1D_t YLim(void) const
y limits
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
string AsString(void) const
double ComputeMaxXSec(const Interaction *in) const
Double_t ymax
Definition: plot.C:25
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 dy[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
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
Kinematical phase space.
Definition: KPhaseSpace.h:34
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:301
static const double kASmallNum
Definition: Controls.h:40
#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
static RunningThreadInfo * Instance(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
double max
Definition: Range1.h:54
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
ifstream in
Definition: comparison.C:7
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
void Configure(const Registry &config)
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
double Energy(const Interaction *in) const
assert(nhit_max >=nhit_nbins)
virtual double XSec(void) const
Definition: GHepRecord.h:127
const InitialState & InitState(void) const
Definition: Interaction.h:69
Double_t ymin
Definition: plot.C:24
double min
Definition: Range1.h:53
#define pNOTICE
Definition: Messenger.h:62
void ClearRunningValues(void)
Definition: Kinematics.cxx:357
bool GetParamDef(const RgKey &name, T &p, const T &def) const
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
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