Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FastPathHadronicCrossSection.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
27 #include "G4ios.hh"
28 #include "G4DynamicParticle.hh"
29 #include "G4Material.hh"
31 #include <vector>
32 #if defined(WIN32)
33  //Needed for M_LN10
34  #define _USE_MATH_DEFINES // for C++
35  #include <math.h>
36 #endif
37 #include <cmath>
38 #include <array>
39 
40 #ifdef FPDEBUG
41 #define DBG( msg ) G4cout<< msg <<G4endl;
42 #define DUMP() G4cout<< <<G4endl;
43 #else
44 #define DBG(msg)
45 #define DUMP()
46 #endif
47 
48 using namespace G4FastPathHadronicCrossSection;
49 
50 //Utility functions used to perform fast-path calculations.
51 //See later for details
52 namespace {
53  struct Point_t {
54  double e;
55  double xs;
56  };
57  int simplify_function(G4double tolerance,
58  std::vector<Point_t> & raw_data,
59  std::vector<Point_t> & simplified_data);
60  void RemoveBias( std::vector <Point_t> &,
61  std::vector <Point_t> &,
62  std::vector <Point_t> &);
63 }
64 
66  particle(part),material(mat),min_cutoff(min),physicsVector(nullptr)
67  {
68  DBG("Initializing a fastPathEntry");
69 #ifdef FPDEBUG
70  count = 0;
71  slowpath_sum=0.;
72  max_delta=0.;
73  min_delta=0.;
74  sum_delta=0.;
75  sum_delta_square=0.;
76 #endif
77 }
78 
80 {
81  DBG("Deleting fastPathEntry");
82  DBG("Dumping status for: "<<(particle?particle->GetParticleName():"PART_NONE")<<" "\
83  <<(material?material->GetName():"MAT_NONE")<<" min_cutoff:"<<min_cutoff<<" "\
84  <<" count:"<<count<<" slowpath_sum:"<<slowpath_sum<<" max_delta:"<<max_delta\
85  <<" min_delta"<<min_delta<<" sum_delta"<<sum_delta<<" sum_delta_squared:"<<sum_delta_square);
86  delete physicsVector;
87 }
88 
89 //namespace {
90 // static inline G4double exp10(G4double x) {
91 // return std::exp( M_LN10*x);
92 // }
93 //}
94 
96 {
97  //Check this method is called when G4CrossSectionDataStore is in the correct state:
98  // FastPath is enabled and we are indeed initializing
101  using std::log10;
102  std::vector<Point_t> data_in;
103  const fastPathParameters& params = xsds->GetFastPathParameters();
104  G4double xs;
105 
106  //G4double max_query = params.queryMax;
107  //G4int count = sampleCount;
108  //G4double tol = dpTol;
109 
110  //Shift so max and min are >= 1.
111  //Don't forget to shift back before computing XS
112  G4double min = params.sampleMin;
113  G4double max = params.sampleMax;
114  G4double shift = 0.0;
115  if(min < 1.0){
116  shift = 1.0 - min;
117  }
118  min += shift;
119  max += shift;
120 
121  G4double log_max = std::log10(params.sampleMax);
122  G4double log_min = std::log10(params.sampleMin);
123  G4double log_step = (log_max-log_min)/(1.0*params.sampleCount);
124 
125  G4double max_xs = 0.0;
126 
127  //Utility particle to calculate XS, with 0 kin energy by default
128  static const G4ThreeVector constDirection(0.,0.,1.);
129  G4DynamicParticle* probingParticle = new G4DynamicParticle( particle , constDirection , 0 );
130 
131  //add the cutoff energy
132  probingParticle->SetKineticEnergy(min_cutoff);
133  //Sample cross-section
134  xs = xsds->GetCrossSection(probingParticle,material);
135  data_in.push_back({min_cutoff,xs});
136 
137  G4double currEnergy = 0.0;
138  //log results
139  auto exp10 = [](G4double x){ return std::exp( M_LN10*x); };
140  for(G4double log_currEnergy = log_min; log_currEnergy < log_max; log_currEnergy += log_step){
141  currEnergy = exp10(log_currEnergy) - shift;
142  if (currEnergy < min_cutoff) continue;
143  probingParticle->SetKineticEnergy(currEnergy);
144  xs=xsds->GetCrossSection(probingParticle,material);
145  //G4cout << "PRUTH: energy value " << currEnergy << ", XS value " << xs << G4endl;
146  if (xs > max_xs) max_xs = xs;
147  data_in.push_back({currEnergy,xs});
148  } // --- end of loop i
149  probingParticle->SetKineticEnergy(max-shift);
150  xs = xsds->GetCrossSection(probingParticle,material);
151  data_in.push_back({max-shift,xs});
152 
153  G4double tol = max_xs * 0.01;
154  std::vector<Point_t> decimated_data;
155  simplify_function(tol, data_in, decimated_data);
156  std::vector<Point_t> debiased_data;
157  RemoveBias( data_in, decimated_data, debiased_data);
158  if ( physicsVector != nullptr ) delete physicsVector;
159  physicsVector = new XSParam(decimated_data.size());
160  G4int physicsVectorIndex = 0;
161  for(size_t i = 0; i < decimated_data.size(); i++){
162  physicsVector->PutValue(physicsVectorIndex++, decimated_data[i].e, decimated_data[i].xs);
163  }
164  //xsds->DumpFastPath(particle,material,G4cout);
165 }
166 
168  particle(pname),material(mat),fastPath(nullptr),
169  energy(-1.),crossSection(-1.)
170 {
171  DBG("Initializing cache entry");
172 #ifdef FPDEBUG
173  cacheHitCount = 0;
174  initCyclesFastPath=0;
175  invocationCountSlowPath=0;
176  totalCyclesSlowPath=0;
177  invocationCountFastPath=0;
178  totalCyclesFastPath=0;
179  invocationCountTriedOneLineCache=0;
180  invocationCountOneLineCache=0;
181 #endif
182 }
183 
185 {
186  DBG("Deleting cache entry");
187  DBG(particle<<" "<<material<<" ("<<(material?material->GetName():"MAT_NONE")<<") "<<" "\
188  <<"fast path pointer:"<<fastPath<<" stored:"<<energy<<" "<<crossSection<<" "\
189  <<cacheHitCount<<" "<<initCyclesFastPath<<" "<<invocationCountSlowPath<<" "\
190  <<totalCyclesSlowPath<<" "<<invocationCountFastPath<<" "<<totalCyclesFastPath<<" "\
191  <<invocationCountTriedOneLineCache<<" "<<invocationCountOneLineCache);
192 }
193 
194 #ifdef FPDEBUG
195 namespace {
196  static inline unsigned long long rdtsc() {
197  unsigned hi=0,lo=0;
198 #if defined(__GNUC__) &&( defined(__i386__)|| defined(__x86_64__) )
199  __asm__ __volatile__ ("rdtsc":"=a"(lo),"=d"(hi));
200 #endif
201  return ((unsigned long long)lo) | ((unsigned long long)hi<<32 );
202  }
203 }
205 {
206  tm.rdtsc_start=rdtsc();
207 }
209 {
210  tm.rdtsc_stop=rdtsc();
211 }
212 #endif
214 #ifdef FPDEBUG
215  methodCalled = 0;
216  hitOneLineCache=0;
217  fastPath=0;
218  slowPath=0;
219  sampleZandA = 0;
220 #endif
221 }
222 
223 namespace {
224 // Rob Fowler's simplify code
225 
226 // This is a curve simplification routine based on the Douglas-Peucker
227 // algorithm.
228 // Simplifying assumptions are that the input polyline is a piecewise
229 // function with the x values monotonically increasing, that the function
230 // reaches an asymptote at the right (high energy) end.
231 // Also, the correct error measure is the difference in y between the original
232 // curve and the result.
233 // In GEANT4 use, the assumption is that the calling program has identified
234 // low- and high-energy cutoffs and that the vector passed in is restricted
235 // to the region between the cutoffs.
236 
237 
238 // The raw_data vector comes in ordered left to right (small energy to large).
239 // The simplified_data vector is initially empty.
240 
241 //A.Dotti ( 16-July-2015): transform variable size C-array and use of size_t
242 // to remove compilation warnings
243 
244 int simplify_function(G4double tolerance,
245  std::vector<Point_t> & raw_data,
246  std::vector<Point_t> & simplified_data)
247 {
248 
249  int gap_left, gap_right; // indices of the current region
250 
251  G4double tolsq = tolerance*tolerance; // Alternative to working with absolute values.
252 
253  std::vector<int> working_stack;
254  //A stack of the points to the right of the current interval that
255  // are known to be selected.
256 
257  gap_right = raw_data.size() - 1; // index of the last element.
258 
259  gap_left = 0;
260 
261  DBG("First and last elements " << gap_left <<" " <<gap_right);
262 
263  simplified_data.push_back(raw_data[0]); //copy first element over.
264 
265  DBG("first point ( 0 "
266  <<simplified_data[0].e <<", "<<simplified_data[0].xs <<" )");
267 
268  working_stack.push_back(gap_right); // 0th element on the stack.
269 
270  while ( !working_stack.empty() )
271  { G4double a, slope, delta;
272  G4double deltasq_max= tolsq;
273  int i_max;
274 
275  gap_right = working_stack.back(); //get current TOS
276  i_max = gap_right;
277 
278 
279  if ( (gap_left +1) < gap_right ) // At least three points in the range.
280  {
281  // co-efficients for the left to right affine line segment
282  slope = (raw_data[gap_right].xs - raw_data[gap_left].xs) /
283  (raw_data[gap_right].e - raw_data[gap_left].e);
284  a = raw_data[gap_left].xs - slope * raw_data[gap_left].e;
285 
286  for ( int i = gap_left +1; i <gap_right; i++) {
287  delta = raw_data[i].xs - a - slope * raw_data[i].e;
288  if ( delta * delta > deltasq_max){
289  deltasq_max = delta * delta;
290  i_max = i;
291  }
292  }
293  } else {
294  DBG(" Less than 3 point interval at [ "<< gap_left <<", " <<gap_right<< " ]");
295  }
296 
297  if(i_max < gap_right) { // Found a new point, push it on the stack
298  working_stack.push_back(i_max);
299  DBG(" pushing point " << i_max);
300  gap_right = i_max;
301  }
302  else { // didn't find a new point betweek gap_left and gap_right.
303  simplified_data.push_back(raw_data[gap_right]);
304  DBG("inserting point ("
305  <<gap_right <<", "<<raw_data[gap_right].e <<", "
306  << raw_data[gap_right].xs <<" )");
307  gap_left = gap_right;
308  working_stack.pop_back();
309  gap_right = working_stack.back();
310  DBG(" new gap_right " << gap_right);
311  }
312 
313  }
314  DBG("Simplified curve size "<< simplified_data.size());
315  return (simplified_data.size());
316 }
317 
318 // Rob Fowler's debias code
319 
320 // This is a de-biasing routine applied after using a curve simplification
321 // routine based on the Douglas-Peucker
322 // algorithm.
323 // Simplifying assumptions are that the input polyline is a piecewise
324 // function with the x values monotonically increasing, and
325 // The right error measure is the difference in y between the original
326 // curve and the result.
327 
328 
329 
330 void RemoveBias(std::vector<Point_t> & original, std::vector<Point_t> & simplified,
331  std::vector<Point_t> & result){
332 
333 
334  const size_t originalSize = original.size();
335  const size_t simplifiedSize = simplified.size();
336 
337  //Create index mapping array
338  std::vector<G4int> xindex(simplifiedSize,0);
339  //G4int xindex[simplifiedSize];
340  G4int lastmatch = 0;
341  G4int j = 0;
342 
343  DBG(" original and simplified vector sizes " << originalSize <<" "<<simplifiedSize);
344 
345  for (size_t k = 0; k <simplifiedSize; k++) {
346  for (size_t i = lastmatch; i < originalSize; i++) {
347  if (original[i].e == simplified[k].e) {
348  xindex[j++] = i;
349  lastmatch = i;
350  }
351  }
352  }
353 
354  DBG("Matched " << j << " values of the simplified vector");
355  // Use short names here.
356  G4int m = simplifiedSize;
357 
358  std::vector<G4double> GArea(m-1,0);
359  //G4double GArea [m-1];
360  G4double GAreatotal = 0;
361 
362  //Area of original simplified curve
363  for(int i = 0; i < m-1; i++){
364  G4double GAreatemp = 0;
365 
366  for(j = xindex[i]; j< xindex[i+1]; j++){
367  G4double trap = (original[j+1].xs + original[j].xs) * (original[j+1].e - original[j].e)/2.0;
368  GAreatemp = GAreatemp + trap;
369  }
370  GArea[i] = GAreatemp;
371  GAreatotal = GAreatotal + GAreatemp;
372  }
373 
374  DBG(" Area under the original curve " << GAreatotal);
375 
376  //aleph Why is this not alpha?
377 
378 
379  std::vector<G4double> aleph(m-1,0);
380  //G4double aleph [m-1];
381  for(int i = 0; i< m-1; i++){
382  aleph[i] = (simplified[i+1].e - simplified[i].e)/2.0;
383  }
384 
385  //solve for f
386  std::vector<G4double> adjustedy(m-1,0);
387  //G4double adjustedy [m];
388  adjustedy[m-1] = simplified[m-1].xs;
389  for(int i = 2; i < m+1; i++) {
390  adjustedy[m-i] = (GArea[m-i]/aleph[m-i]) - adjustedy[m-i+1];
391  if (adjustedy[m-i] <0.0) {
392  adjustedy[m-i] = 0.0;
393  DBG(" Fixing negative cross section at index " << (m-i));
394  }
395  }
396 
397  //error and difference tracking
398  std::vector<G4double> difference(m,0.);
399  //G4double difference [m];
400  G4double maxdiff = 0;
401  G4double adjustedarea = 0;
402  G4double simplifiedarea = 0;
403  for(int i = 0; i < m-1; i++){
404  G4double trap;
405  trap = (adjustedy[i+1]+adjustedy[i])*(simplified[i+1].e-simplified[i].e)/2.0;
406  adjustedarea = adjustedarea+trap;
407  trap = (simplified[i+1].xs+simplified[i].xs)*(simplified[i+1].e-simplified[i].e)/2.0;
408  simplifiedarea = simplifiedarea + trap;
409  }
410 
411  DBG(" Area: Simplified curve = " <<simplifiedarea);
412  DBG(" Area: Debiased curve = " << adjustedarea);
413 
414  for(int i = 0; i <m; i++) {
415  difference[i] = simplified[i].xs-adjustedy[i];
416  }
417  for(int i = 0; i <m; i++){
418  if(std::fabs(difference[i]) > maxdiff) {
419  maxdiff = std::fabs(difference[i]);
420  }
421  }
422  // what is the significance of the loops above ?
423 
424  for(size_t i = 0; i < simplifiedSize; i++){
425  result.push_back( {simplified[i].e , adjustedy[i] } );
426  }
427 
428 }
429 
430 
431 }
G4double G4ParticleHPJENDLHEData::G4double result
const G4FastPathHadronicCrossSection::controlFlag & GetFastPathControlFlags() const
fastPathEntry(const G4ParticleDefinition *par, const G4Material *mat, G4double min_cutoff)
void Initialize(G4CrossSectionDataStore *)
#define DBG(msg)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const G4String & particle
XSParam * physicsVector
const G4String & GetName() const
Definition: G4Material.hh:178
void PutValue(size_t index, G4double energy, G4double dataValue)
~fastPathEntry()
const G4ParticleDefinition *const particle
tuple x
Definition: test.py:50
G4double energy
const G4double min_cutoff
const G4Material *const material
int G4int
Definition: G4Types.hh:78
const G4Material *const material
const G4String & GetParticleName() const
const G4FastPathHadronicCrossSection::fastPathParameters & GetFastPathParameters() const
string material
Definition: eplot.py:19
static constexpr double m
Definition: G4SIunits.hh:129
cycleCountEntry(const G4String &pname, const G4Material *mat)
G4double GetCrossSection(const G4DynamicParticle *, const G4Material *)
void SetKineticEnergy(G4double aEnergy)
string pname
Definition: eplot.py:33
fastPathEntry * fastPath
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
G4double crossSection
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
~cycleCountEntry()
double G4double
Definition: G4Types.hh:76