KineGeneratorWithCache.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2020, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
7  University of Liverpool & STFC Rutherford Appleton Laboratory
8 */
9 //____________________________________________________________________________
10 
11 #include <sstream>
12 #include <cstdlib>
13 #include <map>
14 
15 //#include <TSQLResult.h>
16 //#include <TSQLRow.h>
17 #include <TMath.h>
18 
24 #include "Framework/Utils/Cache.h"
27 
28 using std::ostringstream;
29 using std::map;
30 
31 using namespace genie;
32 
33 //___________________________________________________________________________
36 {
37 
38 }
39 //___________________________________________________________________________
42 {
43 
44 }
45 //___________________________________________________________________________
47 EventRecordVisitorI(name, config)
48 {
49 
50 }
51 //___________________________________________________________________________
53 {
54 
55 }
56 //___________________________________________________________________________
58 {
59  LOG("Kinematics", pINFO)
60  << "Getting max. differential xsec for the rejection method";
61 
62  double xsec_max = -1;
63  Interaction * interaction = event_rec->Summary();
64 
65  LOG("Kinematics", pINFO)
66  << "Attempting to find a cached max{dxsec/dK} value";
67  xsec_max = this->FindMaxXSec(interaction);
68  if(xsec_max>0) return xsec_max;
69 
70  LOG("Kinematics", pINFO)
71  << "Attempting to compute the max{dxsec/dK} value";
72  xsec_max = this->ComputeMaxXSec(interaction);
73  if(xsec_max>0) {
74  LOG("Kinematics", pINFO) << "max{dxsec/dK} = " << xsec_max;
75  this->CacheMaxXSec(interaction, xsec_max);
76  return xsec_max;
77  }
78 
79  LOG("Kinematics", pNOTICE)
80  << "Can not generate event kinematics {K} (max_xsec({K};E)<=0)";
81  // xsec for selected kinematics = 0
82  event_rec->SetDiffXSec(0,kPSNull);
83  // switch on error flag
84  event_rec->EventFlags()->SetBitNumber(kKineGenErr, true);
85  // reset 'trust' bits
86  interaction->ResetBit(kISkipProcessChk);
87  interaction->ResetBit(kISkipKinematicChk);
88  // throw exception
90  exception.SetReason("kinematics generation: max_xsec({K};E)<=0");
91  exception.SwitchOnFastForward();
92  throw exception;
93 
94  return 0;
95 }
96 //___________________________________________________________________________
98  const Interaction * interaction) const
99 {
100 // Find a cached max xsec for the specified xsec algorithm & interaction and
101 // close to the specified energy
102 
103  // get neutrino energy
104  double E = this->Energy(interaction);
105  LOG("Kinematics", pINFO) << "E = " << E;
106 
107  if(E < fEMin) {
108  LOG("Kinematics", pINFO)
109  << "Below minimum energy - Forcing explicit calculation";
110  return -1.;
111  }
112 
113  // access the the cache branch
114  CacheBranchFx * cb = this->AccessCacheBranch(interaction);
115 
116  // if there are enough points stored in the cache buffer to build a
117  // spline, then intepolate
118  if( cb->Spl() ) {
119  if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax()) {
120  double spl_max_xsec = cb->Spl()->Evaluate(E);
121  LOG("Kinematics", pINFO)
122  << "\nInterpolated: max xsec (E=" << E << ") = " << spl_max_xsec;
123  return spl_max_xsec;
124  }
125  LOG("Kinematics", pINFO)
126  << "Outside spline boundaries - Forcing explicit calculation";
127  return -1.;
128  }
129 
130  // if there are not enough points at the cache buffer to have a spline,
131  // look whether there is another point that is sufficiently close
132  double dE = TMath::Min(0.25, 0.05*E);
133  const map<double,double> & fmap = cb->Map();
134  map<double,double>::const_iterator iter = fmap.lower_bound(E);
135  if(iter != fmap.end()) {
136  if(TMath::Abs(E - iter->first) < dE) return iter->second;
137  }
138 
139  return -1;
140 
141 /*
142  // build the search rule
143  double dE = TMath::Min(0.25, 0.05*E);
144  ostringstream search;
145  search << "(x-" << E << " < " << dE << ") && (x>=" << E << ")";
146 
147  // query for all the entries at a window around the current energy
148  TSQLResult * result = cb->Ntuple()->Query("x:y", search.str().c_str());
149  int nrows = result->GetRowCount();
150  LOG("Kinematics", pDEBUG)
151  << "Found " << nrows << " rows with " << search.str();
152  if(nrows <= 0) {
153  delete result;
154  return -1;
155  }
156 
157  // and now select the entry with the closest energy
158  double max_xsec = -1.0;
159  double Ep = 0;
160  double dEmin = 999;
161  TSQLRow * row = 0;
162  while( (row = result->Next()) ) {
163  double cE = atof( row->GetField(0) );
164  double cxsec = atof( row->GetField(1) );
165  double dE = TMath::Abs(E-cE);
166  if(dE < dEmin) {
167  max_xsec = cxsec;
168  Ep = cE;
169  dEmin = TMath::Min(dE,dEmin);
170  }
171  delete row;
172  }
173  delete result;
174 
175  LOG("Kinematics", pINFO)
176  << "\nRetrieved: max xsec = " << max_xsec << " cached at E = " << Ep;
177 
178  return max_xsec;
179 */
180 }
181 //___________________________________________________________________________
183  const Interaction * interaction, double max_xsec) const
184 {
185  LOG("Kinematics", pINFO)
186  << "Adding the computed max{dxsec/dK} value to cache";
187  CacheBranchFx * cb = this->AccessCacheBranch(interaction);
188 
189  double E = this->Energy(interaction);
190  if(max_xsec>0) cb->AddValues(E,max_xsec);
191 
192  if(! cb->Spl() ) {
193  if( cb->Map().size() > 40 ) cb->CreateSpline();
194  }
195 
196  if( cb->Spl() ) {
197  if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() ) {
198  cb->CreateSpline();
199  }
200  }
201 }
202 //___________________________________________________________________________
204 {
205 // Returns the neutrino energy at the struck nucleon rest frame. Kinematic
206 // generators should override this method if they need to cache the max-xsec
207 // values for another energy value (eg kinematic generators for IMD or COH)
208 
209  const InitialState & init_state = interaction->InitState();
210  double E = init_state.ProbeE(kRfHitNucRest);
211  return E;
212 }
213 //___________________________________________________________________________
215  const Interaction * interaction) const
216 {
217 // Returns the cache branch for this algorithm and this interaction. If no
218 // branch is found then one is created.
219 
220  Cache * cache = Cache::Instance();
221 
222  // build the cache branch key as: namespace::algorithm/config/interaction
223  string algkey = this->Id().Key();
224  string intkey = interaction->AsString();
225  string key = cache->CacheBranchKey(algkey, intkey);
226 
227  CacheBranchFx * cache_branch =
228  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
229  if(!cache_branch) {
230  //-- create the cache branch at the first pass
231  LOG("Kinematics", pINFO) << "No Max d^nXSec/d{K}^n cache branch found";
232  LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key;
233 
234  cache_branch = new CacheBranchFx("max[d^nXSec/d^n{K}] over phase space");
235  cache->AddCacheBranch(key, cache_branch);
236  }
237  assert(cache_branch);
238 
239  return cache_branch;
240 }
241 //___________________________________________________________________________
243  const Interaction * interaction, double xsec, double xsec_max) const
244 {
245  // check the computed cross section for the current kinematics against the
246  // maximum cross section used in the rejection MC method for the current
247  // interaction at the current energy.
248  if(xsec>xsec_max) {
249  double f = 200*(xsec-xsec_max)/(xsec_max+xsec);
250  if(f>fMaxXSecDiffTolerance) {
251  LOG("Kinematics", pFATAL)
252  << "xsec: (curr) = " << xsec
253  << " > (max) = " << xsec_max << "\n for " << *interaction;
254  LOG("Kinematics", pFATAL)
255  << "*** Exceeding estimated maximum differential cross section";
256  std::terminate();
257  } else {
258  LOG("Kinematics", pWARN)
259  << "xsec: (curr) = " << xsec
260  << " > (max) = " << xsec_max << "\n for " << *interaction;
261  LOG("Kinematics", pWARN)
262  << "*** The fractional deviation of " << f << " % was allowed";
263  }
264  }
265 
266  // this should never happen - print an error mesg just in case...
267  if(xsec<0) {
268  LOG("Kinematics", pERROR)
269  << "Negative cross section for current kinematics!! \n" << *interaction;
270  }
271 }
272 //___________________________________________________________________________
static QCString name
Definition: declinfo.cpp:673
virtual double MaxXSec(GHepRecord *evrec) const
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
#define pERROR
Definition: Messenger.h:59
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
#define pFATAL
Definition: Messenger.h:56
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
intermediate_table::const_iterator const_iterator
double Evaluate(double x) const
Definition: Spline.cxx:361
string AsString(void) const
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:88
Summary information for an interaction.
Definition: Interaction.h:56
void AddValues(double x, double y)
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
Spline * Spl(void) const
Definition: CacheBranchFx.h:47
def key(type, name=None)
Definition: graph.py:13
static Config * config
Definition: config.cpp:1054
string CacheBranchKey(string k0, string k1="", string k2="") const
Definition: Cache.cxx:93
virtual double ComputeMaxXSec(const Interaction *in) const =0
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:80
GENIE Cache Memory.
Definition: Cache.h:38
double XMax(void) const
Definition: Spline.h:77
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:97
const map< double, double > & Map(void) const
Definition: CacheBranchFx.h:46
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:117
E
Definition: 018_def.c:13
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
virtual double Energy(const Interaction *in) const
virtual CacheBranchFx * AccessCacheBranch(const Interaction *in) const
virtual void CacheMaxXSec(const Interaction *in, double xsec) const
const InitialState & InitState(void) const
Definition: Interaction.h:69
#define pNOTICE
Definition: Messenger.h:61
static Cache * Instance(void)
Definition: Cache.cxx:67
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:37
double XMin(void) const
Definition: Spline.h:76
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
string Key(void) const
Definition: AlgId.h:46
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
virtual double FindMaxXSec(const Interaction *in) const
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