IBDKinematicsGenerator.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: Corey Reed <cjreed \at nikhef.nl>
8  using code from the QELKinematicGenerator written by
9  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
10 
11  For the class documentation see the corresponding header file.
12 
13 */
14 //____________________________________________________________________________
15 
16 #include <TMath.h>
17 
34 
35 using namespace genie;
36 using namespace genie::controls;
37 using namespace genie::constants;
38 using namespace genie::utils;
39 
40 //___________________________________________________________________________
42 KineGeneratorWithCache("genie::IBDKinematicsGenerator")
43 {
44 
45 }
46 //___________________________________________________________________________
48 KineGeneratorWithCache("genie::IBDKinematicsGenerator", config)
49 {
50 
51 }
52 //___________________________________________________________________________
54 {
55 
56 }
57 //___________________________________________________________________________
59 {
60  if(fGenerateUniformly) {
61  LOG("IBD", 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 and set the 'trust' bits
74  Interaction * interaction = evrec->Summary();
75  interaction->SetBit(kISkipProcessChk);
76  interaction->SetBit(kISkipKinematicChk);
77 
78  //-- Note: The kinematic generator would be using the free nucleon cross
79  // section (even for nuclear targets) so as not to double-count nuclear
80  // suppression. This assumes that a) the nuclear suppression was turned
81  // on when computing the cross sections for selecting the current event
82  // and that b) if the event turns out to be unphysical (Pauli-blocked)
83  // the next attempted event will be forced to QEL again.
84  // (discussion with Hugh - GENIE/NeuGEN integration workshop - 07APR2006
85  interaction->SetBit(kIAssumeFreeNucleon);
86 
87  //-- Get the limits for the generated Q2
88  const KPhaseSpace & kps = interaction->PhaseSpace();
89  const Range1D_t Q2 = kps.Limits(kKVQ2);
90 
91  if(Q2.max <=0 || Q2.min>=Q2.max) {
92  LOG("IBD", pWARN) << "No available phase space";
93  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
95  exception.SetReason("No available phase space");
96  exception.SwitchOnFastForward();
97  throw exception;
98  }
99 
100  //-- For the subsequent kinematic selection with the rejection method:
101  // Calculate the max differential cross section or retrieve it from the
102  // cache. Throw an exception and quit the evg thread if a non-positive
103  // value is found.
104  // If the kinematics are generated uniformly over the allowed phase
105  // space the max xsec is irrelevant
106  const double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
107 
108  //-- Try to select a valid Q2 using the rejection method
109 
110  // kinematical limits
111  const double Q2min = Q2.min;
112  const double Q2max = Q2.max;
113  double xsec = -1.;
114  double gQ2 = 0.;
115 
116  unsigned int iter = 0;
117  bool accept = false;
118  while(1) {
119  iter++;
120  if(iter > kRjMaxIterations) {
121  LOG("IBD", pWARN)
122  << "Couldn't select a valid Q^2 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  //-- Generate a Q2 value within the allowed phase space
131  gQ2 = Q2min + (Q2max-Q2min) * rnd->RndKine().Rndm();
132  interaction->KinePtr()->SetQ2(gQ2);
133  LOG("IBD", pINFO) << "Trying: Q^2 = " << gQ2;
134 
135  //-- Computing cross section for the current kinematics
136  xsec = fXSecModel->XSec(interaction, kPSQ2fE);
137 
138  //-- Decide whether to accept the current kinematics
139  if(!fGenerateUniformly) {
140  this->AssertXSecLimits(interaction, xsec, xsec_max);
141  const double t = xsec_max * rnd->RndKine().Rndm();
142 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
143  LOG("IBD", pDEBUG)
144  << "dxsec/dQ2 = " << xsec << ", rnd = " << t;
145 #endif
146  accept = (t < xsec);
147  } else {
148  accept = (xsec>0);
149  }
150 
151  //-- If the generated kinematics are accepted, finish-up module's job
152  if(accept) {
153  LOG("IBD", pINFO) << "Selected: Q^2 = " << gQ2;
154 
155  // reset bits
156  interaction->ResetBit(kISkipProcessChk);
157  interaction->ResetBit(kISkipKinematicChk);
158  interaction->ResetBit(kIAssumeFreeNucleon);
159 
160  // compute the rest of the kinematical variables
161 
162  // get neutrino energy at struck nucleon rest frame and the
163  // struck nucleon mass (can be off the mass shell)
164  const InitialState & init_state = interaction->InitState();
165  const double E = init_state.ProbeE(kRfHitNucRest);
166  const double M = init_state.Tgt().HitNucP4().M();
167 
168  LOG("IBD", pNOTICE) << "E = " << E << ", M = "<< M;
169 
170  // The hadronic inv. mass is equal to the recoil nucleon on-shell mass.
171  const int rpdgc = interaction->RecoilNucleonPdg();
172  assert(rpdgc);
173  const double gW = PDGLibrary::Instance()->Find(rpdgc)->Mass();
174 
175  LOG("IBD", pNOTICE) << "Selected: W = "<< gW;
176 
177  // (W,Q2) -> (x,y)
178  double gx=0, gy=0;
179  kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy);
180 
181  // set the cross section for the selected kinematics
182  evrec->SetDiffXSec(xsec,kPSQ2fE);
183 
184  // for uniform kinematics, compute an event weight as
185  // wght = (phase space volume)*(differential xsec)/(event total xsec)
186  if(fGenerateUniformly) {
187  const double vol = kinematics::PhaseSpaceVolume(interaction,kPSQ2fE);
188  const double totxsec = evrec->XSec();
189  double wght = (vol/totxsec)*xsec;
190  LOG("IBD", pNOTICE) << "Kinematics wght = "<< wght;
191 
192  // apply computed weight to the current event weight
193  wght *= evrec->Weight();
194  LOG("IBD", pNOTICE) << "Current event wght = " << wght;
195  evrec->SetWeight(wght);
196  }
197 
198  // lock selected kinematics & clear running values
199  interaction->KinePtr()->SetQ2(gQ2, true);
200  interaction->KinePtr()->SetW (gW, true);
201  interaction->KinePtr()->Setx (gx, true);
202  interaction->KinePtr()->Sety (gy, true);
203  interaction->KinePtr()->ClearRunningValues();
204 
205  return;
206  }
207  }// iterations
208 }
209 //___________________________________________________________________________
211 {
212  Algorithm::Configure(config);
213  this->LoadConfig();
214 }
215 //____________________________________________________________________________
217 {
218  Algorithm::Configure(config);
219  this->LoadConfig();
220 }
221 //____________________________________________________________________________
223 {
224  //-- Safety factor for the maximum differential cross section
225  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.25 ) ;
226 
227  //-- Minimum energy for which max xsec would be cached, forcing explicit
228  // calculation for lower eneries
229  GetParamDef( "Cache-MinEnergy", fEMin, 1.00 ) ;
230 
231  //-- Maximum allowed fractional cross section deviation from maxim cross
232  // section used in rejection method
233  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
235 
236  //-- Generate kinematics uniformly over allowed phase space and compute
237  // an event weight?
238  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
239 
240 }
241 //____________________________________________________________________________
243  const Interaction * interaction) const
244 {
245 // Computes the maximum differential cross section in the requested phase
246 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
247 // method and the value is cached at a circular cache branch for retrieval
248 // during subsequent event generation.
249 // The computed max differential cross section does not need to be the exact
250 // maximum. The number used in the rejection method will be scaled up by a
251 // safety factor. But it needs to be fast - do not use a very small dQ2 step.
252 
253  double max_xsec = 0.0;
254 
255  const KPhaseSpace & kps = interaction->PhaseSpace();
256  Range1D_t rQ2 = kps.Limits(kKVQ2);
257  if(rQ2.min <=0 || rQ2.max <= rQ2.min) return 0.;
258 
259  //const double logQ2min = TMath::Log(rQ2.min + kASmallNum);
260  //const double logQ2max = TMath::Log(rQ2.max - kASmallNum);
261  const double logQ2min = TMath::Log(rQ2.min);
262  const double logQ2max = TMath::Log(rQ2.max);
263 
264  const int N = 15;
265  const int Nb = 10;
266 
267  double dlogQ2 = (logQ2max - logQ2min) /(N-1);
268  double xseclast = -1;
269  bool increasing;
270 
271  for(int i=0; i<N; i++) {
272  double Q2 = TMath::Exp(logQ2min + i * dlogQ2);
273  interaction->KinePtr()->SetQ2(Q2);
274  double xsec = fXSecModel->XSec(interaction, kPSQ2fE);
275 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
276  LOG("IBD", pDEBUG) << "xsec(Q2= " << Q2 << ") = " << xsec;
277 #endif
278  max_xsec = TMath::Max(xsec, max_xsec);
279  increasing = xsec-xseclast>=0;
280  xseclast = xsec;
281 
282  // once the cross section stops increasing, I reduce the step size and
283  // step backwards a little bit to handle cases that the max cross section
284  // is grossly underestimated (very peaky distribution & large step)
285  if(!increasing) {
286  dlogQ2/=(Nb+1);
287  for(int ib=0; ib<Nb; ib++) {
288  Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
289  if(Q2 < rQ2.min) continue;
290  interaction->KinePtr()->SetQ2(Q2);
291  xsec = fXSecModel->XSec(interaction, kPSQ2fE);
292 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
293  LOG("IBD", pDEBUG) << "xsec(Q2= " << Q2 << ") = " << xsec;
294 #endif
295  max_xsec = TMath::Max(xsec, max_xsec);
296  }
297  break;
298  }
299  }//Q^2
300 
301  // Apply safety factor, since value retrieved from the cache might
302  // correspond to a slightly different energy
303  max_xsec *= fSafetyFactor;
304 
305 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
306  SLOG("IBD", pDEBUG) << interaction->AsString();
307  SLOG("IBD", pDEBUG) << "Max xsec in phase space = " << max_xsec;
308  SLOG("IBD", pDEBUG) << "Computed using alg = " << *fXSecModel;
309 #endif
310 
311  return max_xsec;
312 }
313 //___________________________________________________________________________
314 
virtual double MaxXSec(GHepRecord *evrec) const
void Configure(const Registry &config)
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
virtual void SetWeight(double wght)
Definition: GHepRecord.h:131
bool fGenerateUniformly
uniform over allowed phase space + event weight?
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
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.
int RecoilNucleonPdg(void) const
recoil nucleon pdg
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 ComputeMaxXSec(const Interaction *in) const
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...
void ProcessEventRecord(GHepRecord *event_rec) const
#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
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
Kinematical phase space.
Definition: KPhaseSpace.h:34
#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)
double max
Definition: Range1.h:54
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
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
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
assert(nhit_max >=nhit_nbins)
virtual double XSec(void) const
Definition: GHepRecord.h:127
double gQ2
Definition: gtestDISSF.cxx:55
const int Nb
const InitialState & InitState(void) const
Definition: Interaction.h:69
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
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
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