eic-opticks/u4/ShimG4OpAbsorption.cc

0001 #include <csignal>
0002
0003 #include "G4SystemOfUnits.hh"
0004
0005 #include "ShimG4OpAbsorption.hh"
0006 #include "SEvt.hh"
0007 #include "SLOG.hh"
0008 #include "ssys.h"
0009 #include "U4UniformRand.h"
0010 #include "U4Stack.h"
0011
0012
0013 const plog::Severity ShimG4OpAbsorption::LEVEL = SLOG::EnvLevel("ShimG4OpAbsorption", "DEBUG") ;
0014
0015
0016 ShimG4OpAbsorption::ShimG4OpAbsorption(const G4String& processName, G4ProcessType type )
0017     :
0018     G4OpAbsorption(processName, type)
0019 {
0020 }
0021
0022 ShimG4OpAbsorption::~ShimG4OpAbsorption(){}
0023
0024
0025
0026
0027 //#ifdef DEBUG_TAG
0028 /**
0029 ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft
0030 --------------------------------------------------------
0031
0032 Shim makes process classname appear in SBacktrace.h enabling U4Random::flat/U4Stack::Classify
0033
0034 **/
0035
0036 const bool ShimG4OpAbsorption::FLOAT = ssys::getenvbool("ShimG4OpAbsorption_FLOAT") ;
0037 const int  ShimG4OpAbsorption::PIDX  = ssys::getenvint("PIDX",-1);
0038 const bool ShimG4OpAbsorption::PIDX_ENABLED = ssys::getenvbool("ShimG4OpAbsorption__PIDX_ENABLED") ;
0039
0040
0041 // void ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft(){ G4VProcess::ResetNumberOfInteractionLengthLeft(); }
0042  void ShimG4OpAbsorption::ResetNumberOfInteractionLengthLeft()
0043 {
0044     G4double u = G4UniformRand() ;
0045
0046     LOG(LEVEL)
0047         << U4UniformRand::Desc(u, SEvt::UU )
0048         ;
0049
0050 #ifndef PRODUCTION
0051     SEvt::AddTag( 1, U4Stack_AbsorptionDiscreteReset, u );
0052 #endif
0053
0054     if(FLOAT)
0055     {
0056         float f = -1.f*std::log( float(u) ) ;
0057         theNumberOfInteractionLengthLeft = f ;
0058     }
0059     else
0060     {
0061         theNumberOfInteractionLengthLeft = -1.*G4Log(u) ;
0062     }
0063     theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
0064
0065 }
0066
0067  G4double ShimG4OpAbsorption::GetMeanFreePath(const G4Track& aTrack,  G4double, G4ForceCondition* )
0068 {
0069      G4double len = G4OpAbsorption::GetMeanFreePath( aTrack, 0, 0 );
0070      //std::cout << "ShimG4OpAbsorption::GetMeanFreePath " << len << std::endl ;
0071
0072      //std::raise(SIGINT);
0073
0074      return len ;
0075 }
0076
0077
0078  G4double ShimG4OpAbsorption::PostStepGetPhysicalInteractionLength(const G4Track& track, G4double   previousStepSize, G4ForceCondition* condition)
0079 {
0080   if ( (previousStepSize < 0.0) || (theNumberOfInteractionLengthLeft<=0.0)) {
0081     // beggining of tracking (or just after DoIt of this process)
0082     ResetNumberOfInteractionLengthLeft();
0083   } else if ( previousStepSize > 0.0) {
0084     // subtract NumberOfInteractionLengthLeft
0085     SubtractNumberOfInteractionLengthLeft(previousStepSize);
0086   } else {
0087     // zero step
0088     //  DO NOTHING
0089   }
0090
0091   // condition is set to "Not Forced"
0092   *condition = NotForced;
0093
0094   // get mean free path
0095   currentInteractionLength = GetMeanFreePath(track, previousStepSize, condition);
0096
0097   G4double value;
0098   if (currentInteractionLength <DBL_MAX)
0099   {
0100      if( FLOAT )
0101      {
0102           float fvalue = float(theNumberOfInteractionLengthLeft) * float(currentInteractionLength) ;
0103           value = fvalue ;
0104
0105      }
0106      else
0107      {
0108           value = theNumberOfInteractionLengthLeft * currentInteractionLength ;
0109      }
0110
0111
0112      //std::cout << "ShimG4OpAbsorption::PostStepGetPhysicalInteractionLength " << track.GetTrackID() << std::endl ;
0113
0114       if(track.GetTrackID() - 1 == PIDX && PIDX_ENABLED)
0115       {
0116            std::cout
0117                << "ShimG4OpAbsorption::PostStepGetPhysicalInteractionLength"
0118                << " PIDX " << PIDX
0119                << " currentInteractionLength " << std::setw(10) << std::fixed << std::setprecision(7) << currentInteractionLength
0120                << " theNumberOfInteractionLengthLeft " << std::setw(10) << std::fixed << std::setprecision(7) << theNumberOfInteractionLengthLeft
0121                << " value " << std::setw(10) << std::fixed << std::setprecision(7) << value
0122                << std::endl
0123                ;
0124
0125       }
0126
0127
0128
0129
0130   } else {
0131     value = DBL_MAX;
0132   }
0133 #ifdef G4VERBOSE
0134   if (verboseLevel>1){
0135     G4cout << "G4VDiscreteProcess::PostStepGetPhysicalInteractionLength ";
0136     G4cout << "[ " << GetProcessName() << "]" <<G4endl;
0137     track.GetDynamicParticle()->DumpInfo();
0138     G4cout << " in Material  " <<  track.GetMaterial()->GetName() <<G4endl;
0139     G4cout << "InteractionLength= " << value/cm <<"[cm] " <<G4endl;
0140   }
0141 #endif
0142   return value;
0143
0144 }
0145
0146 //#endif
0147