G4FastPathHadronicCrossSection.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4FastPathHadronicCrossSection.cc

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
#include "G4FastPathHadronicCrossSection.hh"
#include "G4ios.hh"
#include "G4DynamicParticle.hh"
#include "G4Material.hh"
#include "G4CrossSectionDataStore.hh"
#include <vector>
#if defined(WIN32)
  //Needed for M_LN10
  #define _USE_MATH_DEFINES // for C++
  #include <math.h>
#endif
#include <cmath>
#include <array>

#ifdef FPDEBUG
#define DBG( msg ) G4cout<< msg <<G4endl;
#define DUMP() G4cout<< <<G4endl;
#else
#define DBG(msg)
#define DUMP()
#endif

using namespace G4FastPathHadronicCrossSection;

//Utility functions used to perform fast-path calculations.
//See later for details
namespace {
  struct Point_t {
    double e;
    double xs;
  };
  int  simplify_function(G4double tolerance,
                       std::vector<Point_t> & raw_data,
                       std::vector<Point_t>  & simplified_data);
  void RemoveBias( std::vector <Point_t> &,
      std::vector <Point_t> &,
      std::vector <Point_t> &);
}

fastPathEntry::fastPathEntry(const G4ParticleDefinition* part, const G4Material* mat, G4double min) :
    particle(part),material(mat),min_cutoff(min),physicsVector(nullptr)
  {
  DBG("Initializing a fastPathEntry");
#ifdef FPDEBUG
  count = 0;
  slowpath_sum=0.;
  max_delta=0.;
  min_delta=0.;
  sum_delta=0.;
  sum_delta_square=0.;
#endif
}

fastPathEntry::~fastPathEntry()
{
  DBG("Deleting fastPathEntry");
  DBG("Dumping status for: "<<(particle?particle->GetParticleName():"PART_NONE")<<" "\
      <<(material?material->GetName():"MAT_NONE")<<" min_cutoff:"<<min_cutoff<<" "\
      <<" count:"<<count<<" slowpath_sum:"<<slowpath_sum<<" max_delta:"<<max_delta\
      <<" min_delta"<<min_delta<<" sum_delta"<<sum_delta<<" sum_delta_squared:"<<sum_delta_square);
  delete physicsVector;
}

//namespace {
//  static inline G4double exp10(G4double x) {
//    return std::exp( M_LN10*x);
//  }
//}

void fastPathEntry::Initialize(G4CrossSectionDataStore* xsds)
{
  //Check this method is called when G4CrossSectionDataStore is in the correct state:
  //  FastPath is enabled and we are indeed initializing
  assert( xsds->GetFastPathControlFlags().useFastPathIfAvailable &&
      xsds->GetFastPathControlFlags().initializationPhase );
  using std::log10;
  std::vector<Point_t> data_in;
  const fastPathParameters& params = xsds->GetFastPathParameters();
  G4double xs;

    //G4double max_query = params.queryMax;
    //G4int count = sampleCount;
    //G4double tol = dpTol;

  //Shift so max and min are >= 1.
  //Don't forget to shift back before computing XS
  G4double min = params.sampleMin;
  G4double max = params.sampleMax;
  G4double shift = 0.0;
  if(min < 1.0){
    shift = 1.0 - min;
    }
  min += shift;
  max += shift;

  G4double log_max = std::log10(params.sampleMax);
  G4double log_min = std::log10(params.sampleMin);
  G4double log_step = (log_max-log_min)/(1.0*params.sampleCount);

  G4double max_xs = 0.0;

  //Utility particle to calculate XS, with 0 kin energy by default
  static const G4ThreeVector constDirection(0.,0.,1.);
  G4DynamicParticle* probingParticle = new G4DynamicParticle( particle , constDirection , 0 );

  //add the cutoff energy
  probingParticle->SetKineticEnergy(min_cutoff);
  //Sample cross-section
  xs = xsds->GetCrossSection(probingParticle,material);
  data_in.push_back({min_cutoff,xs});

  G4double currEnergy = 0.0;
  //log results
  auto exp10 = [](G4double x){ return std::exp( M_LN10*x); };
  for(G4double log_currEnergy = log_min; log_currEnergy < log_max; log_currEnergy += log_step){
    currEnergy = exp10(log_currEnergy) - shift;
    if (currEnergy <  min_cutoff) continue;
    probingParticle->SetKineticEnergy(currEnergy);
    xs=xsds->GetCrossSection(probingParticle,material);
      //G4cout << "PRUTH: energy value " << currEnergy << ", XS value " << xs << G4endl;
      if (xs > max_xs) max_xs = xs;
      data_in.push_back({currEnergy,xs});
  } // --- end of loop i
  probingParticle->SetKineticEnergy(max-shift);
  xs = xsds->GetCrossSection(probingParticle,material);
  data_in.push_back({max-shift,xs});

  G4double tol = max_xs * 0.01;
  std::vector<Point_t> decimated_data;
  simplify_function(tol,  data_in,  decimated_data);
  std::vector<Point_t> debiased_data;
  RemoveBias( data_in,  decimated_data,  debiased_data);
  if ( physicsVector != nullptr ) delete physicsVector;
  physicsVector = new XSParam(decimated_data.size());
  G4int physicsVectorIndex = 0;
  for(size_t i = 0; i < decimated_data.size(); i++){
    physicsVector->PutValue(physicsVectorIndex++, decimated_data[i].e, decimated_data[i].xs);
  }
  //xsds->DumpFastPath(particle,material,G4cout);
}

cycleCountEntry::cycleCountEntry(const G4String& pname , const G4Material* mat) :
    particle(pname),material(mat),fastPath(nullptr),
    energy(-1.),crossSection(-1.)
{
  DBG("Initializing cache entry");
#ifdef FPDEBUG
  cacheHitCount = 0;
  initCyclesFastPath=0;
  invocationCountSlowPath=0;
  totalCyclesSlowPath=0;
  invocationCountFastPath=0;
  totalCyclesFastPath=0;
  invocationCountTriedOneLineCache=0;
  invocationCountOneLineCache=0;
#endif
}

cycleCountEntry::~cycleCountEntry()
{
  DBG("Deleting cache entry");
  DBG(particle<<" "<<material<<" ("<<(material?material->GetName():"MAT_NONE")<<") "<<" "\
      <<"fast path pointer:"<<fastPath<<" stored:"<<energy<<" "<<crossSection<<" "\
      <<cacheHitCount<<" "<<initCyclesFastPath<<" "<<invocationCountSlowPath<<" "\
      <<totalCyclesSlowPath<<" "<<invocationCountFastPath<<" "<<totalCyclesFastPath<<" "\
      <<invocationCountTriedOneLineCache<<" "<<invocationCountOneLineCache);
}

#ifdef FPDEBUG
namespace {
  static inline unsigned long long rdtsc() {
    unsigned hi=0,lo=0;
#if defined(__GNUC__) &&( defined(__i386__)|| defined(__x86_64__) )
    __asm__ __volatile__ ("rdtsc":"=a"(lo),"=d"(hi));
#endif
    return ((unsigned long long)lo) | ((unsigned long long)hi<<32 );
  }
}
void G4FastPathHadronicCrossSection::logStartCountCycles(timing& tm)
{
  tm.rdtsc_start=rdtsc();
}
void G4FastPathHadronicCrossSection::logStopCountCycles(timing& tm)
{
  tm.rdtsc_stop=rdtsc();
}
#endif
getCrossSectionCount::getCrossSectionCount() {
#ifdef FPDEBUG
  methodCalled = 0;
  hitOneLineCache=0;
  fastPath=0;
  slowPath=0;
  sampleZandA = 0;
#endif
}

namespace {
// Rob Fowler's simplify code

//  This is a curve simplification routine based on the Douglas-Peucker
//  algorithm.
//  Simplifying assumptions are that the input polyline is a piecewise
//  function with the x values monotonically increasing,  that the function
//  reaches an asymptote at the right (high energy) end.
//  Also, the correct error measure is the difference in y between the original
//  curve and the result.
//  In GEANT4 use, the assumption is that the calling program has identified
//  low- and high-energy cutoffs and that the vector passed in is restricted
//  to the region between the cutoffs.


// The raw_data vector comes in ordered left to right (small energy to large).
//  The simplified_data vector is initially empty.

//A.Dotti ( 16-July-2015): transform variable size C-array and use of size_t
//               to remove compilation warnings

int  simplify_function(G4double tolerance,
                       std::vector<Point_t> & raw_data,
                       std::vector<Point_t>  & simplified_data)
{

  int  gap_left, gap_right;  // indices of the current region

  G4double tolsq = tolerance*tolerance;  // Alternative to working with absolute values.

  std::vector<int> working_stack;
  //A stack of the points to the right of the current interval that
  // are known to be selected.

  gap_right = raw_data.size() - 1;  // index of the last element.

  gap_left = 0;

  DBG("First and last elements  " << gap_left <<"  " <<gap_right);

  simplified_data.push_back(raw_data[0]); //copy first element over.

  DBG("first point ( 0  "
          <<simplified_data[0].e  <<",  "<<simplified_data[0].xs <<" )");

  working_stack.push_back(gap_right); // 0th element on the stack.

  while ( !working_stack.empty() )
    {  G4double   a, slope, delta;
      G4double deltasq_max= tolsq;
      int i_max;

      gap_right = working_stack.back();  //get current TOS
      i_max = gap_right;


      if ( (gap_left +1) < gap_right )  // At least three points in the range.
        {
          // co-efficients for the left to right affine line segment
          slope =  (raw_data[gap_right].xs - raw_data[gap_left].xs) /
            (raw_data[gap_right].e - raw_data[gap_left].e);
          a = raw_data[gap_left].xs - slope * raw_data[gap_left].e;

          for ( int i = gap_left +1; i <gap_right;  i++) {
            delta = raw_data[i].xs - a - slope * raw_data[i].e;
            if ( delta * delta > deltasq_max){
              deltasq_max = delta * delta;
              i_max = i;
            }
          }
        }  else {
          DBG("      Less than 3 point interval at [ "<< gap_left <<", " <<gap_right<< " ]");
  }

      if(i_max < gap_right) {  //  Found a new point, push it on the stack
  working_stack.push_back(i_max);
  DBG("         pushing point " << i_max);
  gap_right = i_max;
      }
      else { // didn't find a new point betweek gap_left and gap_right.
  simplified_data.push_back(raw_data[gap_right]);
  DBG("inserting point ("
         <<gap_right  <<",  "<<raw_data[gap_right].e <<", "
         << raw_data[gap_right].xs <<" )");
  gap_left = gap_right;
  working_stack.pop_back();
  gap_right = working_stack.back();
  DBG("      new gap_right " << gap_right);
      }

    }
  DBG("Simplified curve size "<< simplified_data.size());
  return (simplified_data.size());
}

// Rob Fowler's debias code

//  This is a de-biasing routine applied after using a curve simplification
//  routine based on the Douglas-Peucker
//  algorithm.
//  Simplifying assumptions are that the input polyline is a piecewise
//  function with the x values monotonically increasing, and
//  The right error measure is the difference in y between the original
//  curve and the result.



void RemoveBias(std::vector<Point_t> &  original, std::vector<Point_t> & simplified,
    std::vector<Point_t> & result){


  const size_t originalSize = original.size();
  const size_t simplifiedSize = simplified.size();

  //Create index mapping array
  std::vector<G4int>  xindex(simplifiedSize,0);
  //G4int xindex[simplifiedSize];
  G4int lastmatch = 0;
  G4int  j = 0;

  DBG("   original and simplified vector sizes  " << originalSize <<"  "<<simplifiedSize);

  for (size_t k = 0; k <simplifiedSize; k++) {
    for (size_t i = lastmatch; i < originalSize; i++) {
      if (original[i].e == simplified[k].e) {
        xindex[j++] = i;
        lastmatch = i;
      }
    }
  }

  DBG("Matched  " << j << " values of the simplified vector");
  // Use short names here.
  G4int m = simplifiedSize;

  std::vector<G4double> GArea(m-1,0);
  //G4double GArea [m-1];
  G4double GAreatotal = 0;

  //Area of original simplified curve
  for(int i = 0; i < m-1; i++){
    G4double GAreatemp = 0;

    for(j = xindex[i]; j< xindex[i+1]; j++){
      G4double trap = (original[j+1].xs + original[j].xs) * (original[j+1].e - original[j].e)/2.0;
      GAreatemp = GAreatemp + trap;
    }
    GArea[i] = GAreatemp;
    GAreatotal = GAreatotal + GAreatemp;
  }

  DBG("   Area under the original curve " << GAreatotal);

  //aleph        Why is this not alpha?


  std::vector<G4double> aleph(m-1,0);
  //G4double aleph [m-1];
  for(int i = 0; i< m-1; i++){
    aleph[i] = (simplified[i+1].e - simplified[i].e)/2.0;
  }

  //solve for f
  std::vector<G4double> adjustedy(m-1,0);
  //G4double adjustedy [m];
  adjustedy[m-1] = simplified[m-1].xs;
  for(int i = 2; i < m+1; i++) {
    adjustedy[m-i] = (GArea[m-i]/aleph[m-i]) - adjustedy[m-i+1];
    if (adjustedy[m-i] <0.0) {
      adjustedy[m-i] = 0.0;
      DBG("   Fixing negative cross section at index " << (m-i));
    }
  }

  //error and difference tracking
  std::vector<G4double> difference(m,0.);
  //G4double difference [m];
  G4double maxdiff = 0;
  G4double adjustedarea = 0;
  G4double simplifiedarea = 0;
  for(int i = 0; i < m-1; i++){
    G4double trap;
    trap = (adjustedy[i+1]+adjustedy[i])*(simplified[i+1].e-simplified[i].e)/2.0;
    adjustedarea = adjustedarea+trap;
    trap = (simplified[i+1].xs+simplified[i].xs)*(simplified[i+1].e-simplified[i].e)/2.0;
    simplifiedarea = simplifiedarea + trap;
  }

  DBG("   Area:  Simplified curve = " <<simplifiedarea);
  DBG("    Area:  Debiased curve  = " << adjustedarea);

  for(int i = 0; i <m; i++) {
    difference[i] = simplified[i].xs-adjustedy[i];
  }
  for(int i = 0; i <m; i++){
    if(std::fabs(difference[i]) > maxdiff) {
      maxdiff = std::fabs(difference[i]);
    }
  }
  //  what is the significance of the loops above ?

  for(size_t i = 0; i < simplifiedSize; i++){
    result.push_back( {simplified[i].e , adjustedy[i] } );
  }

}


}