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