GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HEDISKinematicsGenerator.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  or see $GENIE/LICENSE
6 
7  Author: Alfonso Garcia <alfonsog \at nikhef.nl>
8  NIKHEF
9 
10  For the class documentation see the corresponding header file.
11 
12 */
13 //____________________________________________________________________________
14 
28 #include "Framework/Utils/Range1.h"
29 #include "Framework/Utils/RunOpt.h"
32 
33 #include <TMath.h>
34 
35 #include <string>
36 
37 using namespace genie;
38 using namespace genie::controls;
39 using namespace genie::utils;
40 
41 //___________________________________________________________________________
43 KineGeneratorWithCache("genie::HEDISKinematicsGenerator")
44 {
45 
46 }
47 //___________________________________________________________________________
49 KineGeneratorWithCache("genie::HEDISKinematicsGenerator", config)
50 {
51 
52 }
53 //___________________________________________________________________________
55 {
56 
57 }
58 //___________________________________________________________________________
60 {
61 
62  //-- Get the random number generators
64 
65  //-- Access cross section algorithm for running thread
67  const EventGeneratorI * evg = rtinfo->RunningThread();
68  fXSecModel = evg->CrossSectionAlg();
69 
70  //-- Get the interaction
71  Interaction * interaction = evrec->Summary();
72  interaction->SetBit(kISkipProcessChk);
73 
74  //-- Get neutrino energy and hit 'nucleon mass'
75  const InitialState & init_state = interaction->InitState();
76  double Ev = init_state.ProbeE(kRfLab);
77  double M = init_state.Tgt().HitNucP4().M(); // can be off m-shell
78 
79  //-- Get the physical W range
80  const KPhaseSpace & kps = interaction->PhaseSpace();
81  Range1D_t xl = kps.XLim();
82  Range1D_t Q2l = kps.Q2Lim();
83 
84  //-- x and y lower limit restrict by limits in SF tables
85  Q2l.min = TMath::Max(Q2l.min,fSFQ2min);
86  Q2l.max = TMath::Min(Q2l.max,fSFQ2max);
87  xl.min = TMath::Max(TMath::Max(xl.min,Q2l.min/2./M/Ev),fSFXmin);
88 
89  LOG("HEDISKinematics", pNOTICE) << "x: [" << xl.min << ", " << xl.max << "]";
90  LOG("HEDISKinematics", pNOTICE) << "log10Q2: [" << Q2l.min << ", " << Q2l.max << "]";
91 
92  //-- For the subsequent kinematic selection with the rejection method:
93 
94  //Scan through a wide region to find the maximum
95  Range1D_t xrange_wide(xl.min*fWideRange,xl.max/fWideRange);
96  Range1D_t Q2range_wide(Q2l.min*fWideRange,Q2l.max/fWideRange);
97  double x_wide = 0.;
98  double Q2_wide = 0.;
99  double xsec_wide = this->Scan(interaction,xrange_wide,Q2range_wide,fWideNKnotsX,fWideNKnotsQ2,2*M*Ev,x_wide,Q2_wide);
100 
101  //Scan through a fine region to find the maximum
102  Range1D_t xrange_fine(TMath::Max(x_wide/fFineRange,xrange_wide.min),TMath::Min(x_wide*fFineRange,xrange_wide.max));
103  Range1D_t Q2range_fine(TMath::Max(Q2_wide/fFineRange,Q2range_wide.min),TMath::Min(Q2_wide*fFineRange,Q2range_wide.max));
104  double x_fine = 0.;
105  double Q2_fine = 0.;
106  double xsec_fine = this->Scan(interaction,xrange_fine,Q2range_fine,fFineNKnotsX,fFineNKnotsQ2,2*M*Ev,x_fine,Q2_fine);
107 
108  //Apply safety factor
109  double xsec_max = fSafetyFactor * TMath::Max(xsec_wide,xsec_fine);
110 
111  //-- Try to select a valid (x,y) pair using the rejection method
112  double log10xmin = TMath::Log10(xl.min);
113  double log10xmax = TMath::Log10(xl.max);
114  double log10Q2min = TMath::Log10(Q2l.min);
115  double log10Q2max = TMath::Log10(Q2l.max);
116 
117  double dlog10x = log10xmax - log10xmin;
118  double dlog10Q2 = log10Q2max - log10Q2min;
119 
120  double gx=-1, gQ2=-1, xsec=-1;
121 
122  unsigned int iter = 0;
123  bool accept = false;
124  while(1) {
125  iter++;
126  if(iter > kRjMaxIterations) {
127  LOG("HEDISKinematics", pWARN)
128  << " Couldn't select kinematics after " << iter << " iterations";
129  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
131  exception.SetReason("Couldn't select kinematics");
132  exception.SwitchOnFastForward();
133  throw exception;
134  }
135 
136  gx = TMath::Power( 10., log10xmin + dlog10x * rnd->RndKine().Rndm() );
137  gQ2 = TMath::Power( 10., log10Q2min + dlog10Q2 * rnd->RndKine().Rndm() );
138 
139  interaction->KinePtr()->Setx(gx);
140  interaction->KinePtr()->SetQ2(gQ2);
141  kinematics::UpdateWYFromXQ2(interaction);
142 
143  LOG("HEDISKinematics", pDEBUG)
144  << "Trying: x = " << gx << ", Q2 = " << gQ2
145  << " (W = " << interaction->KinePtr()->W() << ","
146  << " y = " << interaction->KinePtr()->y() << ")";
147 
148  //-- compute the cross section for current kinematics
149  xsec = fXSecModel->XSec(interaction, kPSlog10xlog10Q2fE);
150 
151  //-- decide whether to accept the current kinematics
152  this->AssertXSecLimits(interaction, xsec, xsec_max);
153 
154  double t = xsec_max * rnd->RndKine().Rndm();
155 
156  LOG("HEDISKinematics", pDEBUG) << "xsec= " << xsec << ", Rnd= " << t;
157 
158  accept = (t < xsec);
159 
160  //-- If the generated kinematics are accepted, finish-up module's job
161  if(accept) {
162  LOG("HEDISKinematics", pNOTICE)
163  << "Selected: x = " << gx << ", Q2 = " << gQ2
164  << " (W = " << interaction->KinePtr()->W() << ","
165  << " (Y = " << interaction->KinePtr()->y() << ")";
166 
167  // reset trust bits
168  interaction->ResetBit(kISkipProcessChk);
169  interaction->ResetBit(kISkipKinematicChk);
170 
171  // set the cross section for the selected kinematics
172  evrec->SetDiffXSec(xsec,kPSxQ2fE);
173 
174  // lock selected kinematics & clear running values
175  interaction->KinePtr()->SetW (interaction->KinePtr()->W(), true);
176  interaction->KinePtr()->Sety (interaction->KinePtr()->y(), true);
177  interaction->KinePtr()->SetQ2(gQ2, true);
178  interaction->KinePtr()->Setx (gx, true);
179  interaction->KinePtr()->ClearRunningValues();
180  return;
181  }
182  } // iterations
183 }
184 //___________________________________________________________________________
185 double HEDISKinematicsGenerator::Scan(Interaction * interaction, Range1D_t xrange,Range1D_t Q2range, int NKnotsQ2, int NKnotsX, double ME2, double & x_scan, double & Q2_scan) const
186 {
187 
188  double xsec_scan = 0.;
189  Q2_scan = 0.;
190  x_scan = 0.;
191 
192  double stepQ2 = TMath::Power(Q2range.max/Q2range.min,1./(NKnotsQ2-1));
193 
194  for ( int iq=0; iq<NKnotsQ2; iq++) {
195  double Q2_aux = Q2range.min*TMath::Power(stepQ2,iq);
196  xrange.min = TMath::Max(xrange.min,Q2_aux/ME2);
197  double stepx = TMath::Power(xrange.max/xrange.min,1./(NKnotsX-1));
198  for ( int ix=0; ix<NKnotsX; ix++) {
199  double x_aux = xrange.min*TMath::Power(stepx,ix);
200  interaction->KinePtr()->Setx(x_aux);
201  interaction->KinePtr()->SetQ2(Q2_aux);
202  kinematics::UpdateWYFromXQ2(interaction);
203  double xsec_aux = fXSecModel->XSec(interaction, kPSlog10xlog10Q2fE);
204  LOG("HEDISKinematics", pDEBUG) << "x = " << x_aux << " , Q2 = " << Q2_aux << ", xsec = " << xsec_aux;
205  if (xsec_aux>xsec_scan) {
206  xsec_scan = xsec_aux;
207  Q2_scan = Q2_aux;
208  x_scan = x_aux;
209  }
210  }
211  }
212 
213  LOG("HEDISKinematics", pNOTICE) << "scan -> x = " << x_scan << " , Q2 = " << Q2_scan << ", xsec = " << xsec_scan;
214 
215  return xsec_scan;
216 
217 }
218 //___________________________________________________________________________
219 double HEDISKinematicsGenerator::ComputeMaxXSec(const Interaction * /* interaction */ ) const
220 {
221  return 0;
222 }
223 //___________________________________________________________________________
225 {
226  Algorithm::Configure(config);
227  this->LoadConfig();
228 }
229 //____________________________________________________________________________
231 {
232  Algorithm::Configure(config);
233  this->LoadConfig();
234 }
235 //____________________________________________________________________________
237 {
238 
239  //-- Safety factor for the maximum differential cross section
240  GetParamDef("MaxXSec-SafetyFactor", fSafetyFactor, 2. ) ;
241  //-- Maximum allowed fractional cross section deviation from maxim cross
242  // section used in rejection method
243  GetParamDef("MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
244  assert(fMaxXSecDiffTolerance>=0);
245 
246  GetParamDef("ScanWide-NKnotsX", fWideNKnotsX, 10 ) ;
247  GetParamDef("ScanWide-NKnotsQ2", fWideNKnotsQ2, 10 ) ;
248  GetParamDef("ScanWide-Range", fWideRange, 1.1 ) ;
249  GetParamDef("ScanFine-NKnotsX", fFineNKnotsX, 10 ) ;
250  GetParamDef("ScanFine-NKnotsQ2", fFineNKnotsQ2, 10 ) ;
251  GetParamDef("ScanFine-Range", fFineRange, 10. ) ;
252 
253  // Limits from the SF tables that are useful to reduce computation
254  // time of the max cross section
255  GetParam("XGrid-Min", fSFXmin ) ;
256  GetParam("Q2Grid-Min", fSFQ2min ) ;
257  GetParam("Q2Grid-Max", fSFQ2max ) ;
258 
259 }
double W(bool selected=false) const
Definition: Kinematics.cxx:157
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
double ComputeMaxXSec(const Interaction *interaction) const
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
double Scan(Interaction *interaction, Range1D_t xrange, Range1D_t Q2range, int NKnotsQ2, int NKnotsX, double ME2, double &x_scan, double &Q2_scan) const
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.
double fSFXmin
minimum value of x for which SF tables are computed
Range1D_t Q2Lim(void) const
Q2 limits.
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
void UpdateWYFromXQ2(const Interaction *in)
Definition: KineUtils.cxx:1326
double y(bool selected=false) const
Definition: Kinematics.cxx:112
void Configure(const Registry &config)
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
void ProcessEventRecord(GHepRecord *event_rec) const
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
double fSFQ2min
minimum value of Q2 for which SF tables are computed
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
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
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
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
double gQ2
Definition: gtestDISSF.cxx:55
Range1D_t XLim(void) const
x limits
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
double fSFQ2max
maximum value of Q2 for which SF tables are computed
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
const EventGeneratorI * RunningThread(void)
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