G4FastPathHadronicCrossSection.hh

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

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
#ifndef G4FastPathHadronicCrossSection_hh
#define G4FastPathHadronicCrossSection_hh

#include "G4PhysicsFreeVector.hh"
#include "G4ParticleDefinition.hh"
#include "G4Material.hh"
#include <functional>
#include <utility>
#include <unordered_map>
#include <iostream>
#include <set>
#include <stdint.h>

class G4DynamicParticle;
class G4Material;
class G4CrossSectionDataStore;

//To measure performances and debug info on fast cross-section enable this
//#define FPDEBUG

//TODO: Move all logging and debug functionality to separate header
namespace G4FastPathHadronicCrossSection {
  //This data type contains the simplified representation of the
  //cross-section, by default it is a G4PhysicsVector type
  using XSParam=G4PhysicsFreeVector;
  //The key used to search in the cache.
  using G4CrossSectionDataStore_Key=std::pair<const G4ParticleDefinition*,const G4Material*>;
  //This represents the fast XS implementation.
  struct fastPathEntry{
    //fastPathEntry();
    fastPathEntry(const G4ParticleDefinition *par,const G4Material* mat,G4double min_cutoff);
    ~fastPathEntry();
    inline G4double GetCrossSection(G4double ene) const { return physicsVector->Value(ene); }
    void Initialize(G4CrossSectionDataStore* );
    const G4ParticleDefinition * const particle;
    const G4Material * const material;
    const G4double min_cutoff;

    XSParam *physicsVector;
#       ifdef FPDEBUG
    //stats for debug
    G4int count;
    G4double slowpath_sum; //sum of all slowpath xs
    G4double max_delta;
    G4double min_delta;
    G4double sum_delta;
    G4double sum_delta_square;
#   endif
  };

  //A cache entry.
  struct cycleCountEntry{
    cycleCountEntry(const G4String& pname , const G4Material* mat);
    ~cycleCountEntry();
    const G4String& particle;
    const G4Material * const material;

    //optional fastPathEntry
    fastPathEntry* fastPath;

    //cache per element of material test
    G4double energy;
    G4double crossSection;
#   ifdef FPDEBUG
    uint64_t cacheHitCount;//
    uint64_t initCyclesFastPath;
    uint64_t invocationCountSlowPath;
    uint64_t totalCyclesSlowPath;
    uint64_t invocationCountFastPath;
    uint64_t totalCyclesFastPath;
    uint64_t invocationCountTriedOneLineCache;//
    uint64_t invocationCountOneLineCache;//
#   endif
  };

  struct timing {
    unsigned long long rdtsc_start;
    unsigned long long rdtsc_stop;
  };

  struct getCrossSectionCount {
    getCrossSectionCount();
    inline void MethodCalled();
    inline void HitOneLine();
    inline void FastPath();
    inline void SlowPath();
    inline void SampleZandA();
#ifdef FPDEBUG
    uint64_t methodCalled;
    uint64_t hitOneLineCache;
    uint64_t fastPath;
    uint64_t slowPath;
    uint64_t sampleZandA;
#endif
  };

  //Hashing the key
  struct G4CrossSectionDataStore_Key_Hash {
    std::hash<uint64_t> hash_uint64_t;
    inline size_t operator()(const G4CrossSectionDataStore_Key& x) const throw() {
      return hash_uint64_t(hash_uint64_t( ((uint64_t)(x.first)) ) +  hash_uint64_t(((uint64_t)(x.second))));
    }
  };
  //Equality for two key elements
  struct G4CrossSectionDataStore_Key_EqualTo {
    inline bool operator()(const G4CrossSectionDataStore_Key& lhs, const G4CrossSectionDataStore_Key& rhs ) const {
      //TODO: Verify this: particles are singletons, materials use operator==
            //TODO: in ref-10, G4Material::operator== becomes deleted, investigating why
      return (lhs.first==rhs.first)&&(lhs.second == rhs.second);
    }
  };
//  The cache itself
  using G4CrossSectionDataStore_Cache=std::unordered_map<G4CrossSectionDataStore_Key,cycleCountEntry*,
      G4CrossSectionDataStore_Key_Hash,G4CrossSectionDataStore_Key_EqualTo>;

  struct fastPathRequestConfig_t {
    G4CrossSectionDataStore_Key part_mat;
    G4double min_cutoff;
  };
  //Two of the elements are identical if the part_mat part is
  struct fastPathRequestConfig_Less {
    std::less<G4CrossSectionDataStore_Key> less;
    inline bool operator()(const fastPathRequestConfig_t& lhs,const fastPathRequestConfig_t& rhs ) const {
      return less(lhs.part_mat,rhs.part_mat);
    }
  };
  using G4CrossSectionDataStore_Requests=std::set<fastPathRequestConfig_t,fastPathRequestConfig_Less>;

  //Configure the caching mechanism
  struct controlFlag {
    G4bool prevCalcUsedFastPath;
    G4bool useFastPathIfAvailable;
    G4bool initializationPhase;
    controlFlag() : prevCalcUsedFastPath(false),useFastPathIfAvailable(false),initializationPhase(false) {}
  };
  //Parameters to control sampling
  struct fastPathParameters {
    fastPathParameters() {
      //default
      //TODO: are these ok?
        queryMax = 10000;
        sampleMin = 0.0001;
        sampleMax = 10000;
        sampleCount = 200000;
        dpTol = 0.01;
    }
    //PRUTH vars for sampling and surragate model
    G4double queryMax;
    G4double sampleMin;
    G4double sampleMax;
    G4int sampleCount;
    G4double dpTol;
  };

  //Logging functionalities, disabled if not in FPDEBUG mode
  static inline void logInvocationTriedOneLine( cycleCountEntry* );
  static inline void logInvocationOneLine( cycleCountEntry* );
  static inline void logHit(cycleCountEntry*);
  static inline void logInvocationCountFastPath( cycleCountEntry* );
  static inline void logInvocationCountSlowPAth( cycleCountEntry* );
#ifdef FPDEBUG
  void logStartCountCycles( timing& );
  void logStopCountCycles( timing& );
#else
  inline void logStartCountCycles(timing&) {}
  inline void logStopCountCycles(timing&) {}
#endif
  static inline void logInitCyclesFastPath( cycleCountEntry* , timing& );
  static inline void logTotalCyclesFastPath( cycleCountEntry* , timing& );
  static inline void logTotalCyclesSlowPath( cycleCountEntry* , timing& );
  static inline void logTiming( cycleCountEntry* , fastPathEntry* , timing& );
}

inline std::ostream& operator<<(std::ostream& os, const G4FastPathHadronicCrossSection::fastPathEntry& fp);

//Implementation of inline functions. Note the ifdef

namespace G4FastPathHadronicCrossSection {

#ifdef FPDEBUG
inline void logInvocationTriedOneLine(cycleCountEntry* cl ) {
  if ( cl != nullptr ) ++(cl->invocationCountTriedOneLineCache);
}
inline void logInvocationOneLine( cycleCountEntry* cl ) {
  if ( cl != nullptr ) ++(cl->invocationCountOneLineCache);
}
inline void logHit(cycleCountEntry* cl) {
  if ( cl != nullptr ) ++(cl->cacheHitCount);
}
inline void logInvocationCountFastPath( cycleCountEntry* cl )
{
  if ( cl != nullptr ) ++(cl->invocationCountFastPath);
}
inline void logInvocationCountSlowPAth( cycleCountEntry* cl)
{
  if ( cl != nullptr ) ++(cl->invocationCountSlowPath);
}

inline void logInitCyclesFastPath(cycleCountEntry* cl,timing& tm)
{
  if ( cl != nullptr ) cl->initCyclesFastPath = tm.rdtsc_stop - tm.rdtsc_start;
}
inline void logTotalCyclesFastPath( cycleCountEntry* cl,timing& tm)
{
  if ( cl!=nullptr ) cl->totalCyclesFastPath = tm.rdtsc_stop - tm.rdtsc_start;
}
inline void logTotalCyclesSlowPath( cycleCountEntry* cl,timing& tm)
{
  if ( cl!=nullptr ) cl->totalCyclesSlowPath = tm.rdtsc_stop - tm.rdtsc_start;
}
inline void logTiming( cycleCountEntry* entry , fastPathEntry* fast_entry, timing& timing)
{
  if (fast_entry != nullptr ) {
    if ( entry->invocationCountFastPath == 0 ) {
      //PRUTH style initialization
      G4FastPathHadronicCrossSection::logInitCyclesFastPath(entry,timing);
      G4FastPathHadronicCrossSection::logInvocationCountFastPath(entry);
      } else {
      //PRUTH comment to understand:
      //the first one includes the initialization... don't count it for now
      G4FastPathHadronicCrossSection::logTotalCyclesFastPath(entry,timing);
      G4FastPathHadronicCrossSection::logInvocationCountFastPath(entry);
      }
       } else {
         G4FastPathHadronicCrossSection::logInvocationCountSlowPAth(entry);
         G4FastPathHadronicCrossSection::logTotalCyclesSlowPath(entry,timing);
       }
}
#else
inline void logInvocationTriedOneLine(cycleCountEntry*){}
inline void logInvocationOneLine( cycleCountEntry*){}
inline void logHit(cycleCountEntry*){}
inline void logInvocationCountFastPath( cycleCountEntry*){}
inline void logInvocationCountSlowPAth( cycleCountEntry*){}
inline void logInitCyclesFastPath( cycleCountEntry* , timing& ){}
inline void logTotalCyclesFastPath( cycleCountEntry* , timing& ){}
inline void logTotalCyclesSlowPath( cycleCountEntry* , timing& ){}
inline void logTiming( cycleCountEntry* , fastPathEntry* , timing& ) {}
#endif

inline void getCrossSectionCount::MethodCalled() {
#ifdef FPDEBUG
  ++methodCalled;
#endif
}

inline void getCrossSectionCount::HitOneLine() {
#ifdef FPDEBUG
  ++hitOneLineCache;
#endif
}

inline void getCrossSectionCount::FastPath() {
#ifdef FPDEBUG
  ++fastPath;
#endif
}

inline void getCrossSectionCount::SlowPath() {
#ifdef FPDEBUG
  ++slowPath;
#endif
}

inline void getCrossSectionCount::SampleZandA() {
#ifdef FPDEBUG
  ++sampleZandA;
#endif
}
}//namespace

inline std::ostream& operator<<(std::ostream& os, const G4FastPathHadronicCrossSection::fastPathEntry& fp) {
  using CLHEP::MeV;
  os<<"#Particle: "<<(fp.particle!=nullptr?fp.particle->GetParticleName():"UNDEFINED")<<"\n";
  os<<"#Material: "<<(fp.material!=nullptr?fp.material->GetName():"UNDEFINED")<<"\n";
  os<<"#min_cutoff(MeV): "<<fp.min_cutoff/MeV<<"\n";
#ifdef FPDEBUG
  os<<"#DEBUG COUNTERS: count="<<fp.count<<" slowpath_sum="<<fp.slowpath_sum<<" max_delta="<<fp.max_delta;
  os<<" min_delta="<<fp.min_delta<<" sum_delta="<<fp.sum_delta<<" sum_delta_square="<<fp.sum_delta_square<<"\n";
#endif
  os<<*(fp.physicsVector)<<"\n";
  return os;
}

#endif //G4FastPathHadronicCrossSection_hh