eic-opticks/u4/U4Physics.cc

0001 /**
0002 U4Physics.cc
0003 ==============
0004
0005 Boundary class changes need to match in all the below::
0006
0007     U4OpBoundaryProcess.h
0008     U4Physics.cc
0009     U4Recorder.cc
0010     U4StepPoint.cc
0011
0012 **/
0013
0014
0015 #include <iomanip>
0016 #include "ssys.h"
0017 #include "sstr.h"
0018 #include "U4Physics.hh"
0019 #include "U4OpBoundaryProcess.h"
0020 #include "G4ProcessManager.hh"
0021 #include "G4FastSimulationManagerProcess.hh"
0022
0023
0024 #if defined(WITH_CUSTOM4) && defined(WITH_PMTSIM)
0025 #include "G4OpBoundaryProcess.hh"
0026 #include "C4OpBoundaryProcess.hh"
0027 #include "PMTSimParamSvc/PMTAccessor.h"
0028 #elif defined(WITH_CUSTOM4) && !defined(WITH_PMTSIM)
0029 #include "G4OpBoundaryProcess.hh"
0030 #include "C4OpBoundaryProcess.hh"
0031 #include "SPMTAccessor.h"
0032 #elif defined(WITH_INSTRUMENTED_DEBUG)
0033 #include "InstrumentedG4OpBoundaryProcess.hh"
0034 #else
0035 #include "G4OpBoundaryProcess.hh"
0036 #endif
0037
0038
0039 #include "SLOG.hh"
0040 const plog::Severity U4Physics::LEVEL = SLOG::EnvLevel("U4Physics", "DEBUG") ;
0041
0042 U4Physics::U4Physics()
0043     :
0044     fCerenkov(nullptr),
0045     fScintillation(nullptr),
0046     fAbsorption(nullptr),
0047     fRayleigh(nullptr),
0048     fBoundary(nullptr),
0049     fFastSim(nullptr)
0050 {
0051     Cerenkov_DISABLE = EInt(_Cerenkov_DISABLE, "0") ;
0052     Scintillation_DISABLE = EInt(_Scintillation_DISABLE, "0" );
0053     OpAbsorption_DISABLE = EInt(_OpAbsorption_DISABLE, "0") ;
0054     OpRayleigh_DISABLE = EInt(_OpRayleigh_DISABLE, "0") ;
0055     OpBoundaryProcess_DISABLE = EInt(_OpBoundaryProcess_DISABLE, "0") ;
0056     OpBoundaryProcess_LASTPOST = EInt(_OpBoundaryProcess_LASTPOST, "0") ;
0057     FastSim_ENABLE = EInt(_FastSim_ENABLE, "0") ;
0058 }
0059
0060
0061
0062
0063 #include "G4BosonConstructor.hh"
0064 #include "G4LeptonConstructor.hh"
0065 #include "G4MesonConstructor.hh"
0066 #include "G4BaryonConstructor.hh"
0067 #include "G4IonConstructor.hh"
0068
0069
0070 void U4Physics::ConstructParticle()
0071 {
0072     G4BosonConstructor::ConstructParticle();
0073     G4LeptonConstructor::ConstructParticle();
0074     G4MesonConstructor::ConstructParticle();
0075     G4BaryonConstructor::ConstructParticle();
0076     G4IonConstructor::ConstructParticle();
0077 }
0078
0079 void U4Physics::ConstructProcess()
0080 {
0081     AddTransportation();
0082     ConstructEM();
0083     ConstructOp();
0084 }
0085
0086 // from OpNovicePhysicsList::ConstructEM
0087
0088 #include "G4ComptonScattering.hh"
0089 #include "G4GammaConversion.hh"
0090 #include "G4PhotoElectricEffect.hh"
0091
0092 #include "G4eMultipleScattering.hh"
0093 #include "G4MuMultipleScattering.hh"
0094 #include "G4hMultipleScattering.hh"
0095
0096 #include "G4eIonisation.hh"
0097 #include "G4eBremsstrahlung.hh"
0098 #include "G4eplusAnnihilation.hh"
0099
0100 #include "G4MuIonisation.hh"
0101 #include "G4MuBremsstrahlung.hh"
0102 #include "G4MuPairProduction.hh"
0103
0104 #include "G4hIonisation.hh"
0105
0106 void U4Physics::ConstructEM()
0107 {
0108     G4int em_verbosity = 0 ;
0109     G4EmParameters* empar = G4EmParameters::Instance() ;
0110     empar->SetVerbose(em_verbosity);
0111     empar->SetWorkerVerbose(em_verbosity);
0112
0113   auto particleIterator=GetParticleIterator();
0114   particleIterator->reset();
0115   while( (*particleIterator)() )
0116   {
0117     G4ParticleDefinition* particle = particleIterator->value();
0118
0119     G4ProcessManager* pmanager = particle->GetProcessManager();
0120     G4String particleName = particle->GetParticleName();
0121
0122     if (particleName == "gamma") {
0123     // gamma
0124       // Construct processes for gamma
0125       pmanager->AddDiscreteProcess(new G4GammaConversion());
0126       pmanager->AddDiscreteProcess(new G4ComptonScattering());
0127       pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());
0128
0129     } else if (particleName == "e-") {
0130     //electron
0131       // Construct processes for electron
0132       pmanager->AddProcess(new G4eMultipleScattering(),-1, 1, 1);
0133       pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
0134       pmanager->AddProcess(new G4eBremsstrahlung(),   -1, 3, 3);
0135
0136     } else if (particleName == "e+") {
0137     //positron
0138       // Construct processes for positron
0139       pmanager->AddProcess(new G4eMultipleScattering(),-1, 1, 1);
0140       pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
0141       pmanager->AddProcess(new G4eBremsstrahlung(),   -1, 3, 3);
0142       pmanager->AddProcess(new G4eplusAnnihilation(),  0,-1, 4);
0143
0144     } else if( particleName == "mu+" ||
0145                particleName == "mu-"    ) {
0146     //muon
0147      // Construct processes for muon
0148      pmanager->AddProcess(new G4MuMultipleScattering(),-1, 1, 1);
0149      pmanager->AddProcess(new G4MuIonisation(),      -1, 2, 2);
0150      pmanager->AddProcess(new G4MuBremsstrahlung(),  -1, 3, 3);
0151      pmanager->AddProcess(new G4MuPairProduction(),  -1, 4, 4);
0152
0153     } else {
0154       if ((particle->GetPDGCharge() != 0.0) &&
0155           (particle->GetParticleName() != "chargedgeantino") &&
0156           !particle->IsShortLived()) {
0157        // all others charged particles except geantino
0158        pmanager->AddProcess(new G4hMultipleScattering(),-1,1,1);
0159        pmanager->AddProcess(new G4hIonisation(),       -1,2,2);
0160      }
0161     }
0162   }
0163 }
0164
0165 #include "Local_G4Cerenkov_modified.hh"
0166 #include "Local_DsG4Scintillation.hh"
0167
0168 #include "ShimG4OpAbsorption.hh"
0169 #include "ShimG4OpRayleigh.hh"
0170
0171
0172 std::string U4Physics::desc() const
0173 {
0174     std::stringstream ss ;
0175     ss
0176         << "U4Physics::desc" << "\n"
0177         << std::setw(60) << _Cerenkov_DISABLE           << " : " << Cerenkov_DISABLE << "\n"
0178         << std::setw(60) << _Scintillation_DISABLE      << " : " << Scintillation_DISABLE << "\n"
0179         << std::setw(60) << _OpAbsorption_DISABLE       << " : " << OpAbsorption_DISABLE << "\n"
0180         << std::setw(60) << _OpRayleigh_DISABLE         << " : " << OpRayleigh_DISABLE << "\n"
0181         << std::setw(60) << _OpBoundaryProcess_DISABLE  << " : " << OpBoundaryProcess_DISABLE << "\n"
0182         << std::setw(60) << _OpBoundaryProcess_LASTPOST << " : " << OpBoundaryProcess_LASTPOST << "\n"
0183         << std::setw(60) << _FastSim_ENABLE             << " : " << FastSim_ENABLE << "\n"
0184         ;
0185     std::string str = ss.str();
0186     return str ;
0187 }
0188
0189
0190 std::string U4Physics::Desc()  // static
0191 {
0192     std::stringstream ss ;
0193 #ifdef DEBUG_TAG
0194     ss << ( ShimG4OpAbsorption::FLOAT ? "ShimG4OpAbsorption_FLOAT" : "ShimG4OpAbsorption_ORIGINAL" ) ;
0195     ss << "_" ;
0196     ss << ( ShimG4OpRayleigh::FLOAT ? "ShimG4OpRayleigh_FLOAT" : "ShimG4OpRayleigh_ORIGINAL" ) ;
0197 #endif
0198     std::string str = ss.str();
0199     return str ;
0200 }
0201
0202
0203
0204
0205 std::string U4Physics::Switches()  // static
0206 {
0207     std::stringstream ss ;
0208     ss << "U4Physics::Switches" << std::endl ;
0209 #if defined(WITH_CUSTOM4)
0210     ss << "WITH_CUSTOM4" << std::endl ;
0211 #else
0212     ss << "NOT:WITH_CUSTOM4" << std::endl ;
0213 #endif
0214 #if defined(WITH_PMTSIM)
0215     ss << "WITH_PMTSIM" << std::endl ;
0216 #else
0217     ss << "NOT:WITH_PMTSIM" << std::endl ;
0218 #endif
0219 #if defined(WITH_CUSTOM4) && defined(WITH_PMTSIM)
0220     ss << "WITH_CUSTOM4_AND_WITH_PMTSIM" << std::endl ;
0221 #else
0222     ss << "NOT:WITH_CUSTOM4_AND_WITH_PMTSIM" << std::endl ;
0223 #endif
0224
0225 #if defined(WITH_CUSTOM4) && !defined(WITH_PMTSIM)
0226     ss << "WITH_CUSTOM4_AND_NOT_WITH_PMTSIM" << std::endl ;
0227 #else
0228     ss << "NOT:WITH_CUSTOM4_AND_NOT_WITH_PMTSIM" << std::endl ;
0229 #endif
0230
0231 #if defined(DEBUG_TAG)
0232     ss << "DEBUG_TAG" << std::endl ;
0233 #else
0234     ss << "NOT:DEBUG_TAG" << std::endl ;
0235 #endif
0236     std::string str = ss.str();
0237     return str ;
0238 }
0239
0240
0241
0242 int U4Physics::EInt(const char* key, const char* fallback)  // static
0243 {
0244     const char* val_ = getenv(key) ;
0245     int val =  std::atoi(val_ ? val_ : fallback) ;
0246     return val ;
0247 }
0248
0249
0250 /**
0251 U4Physics::ConstructOp
0252 -----------------------
0253
0254 Scintillation needs to come after absorption for reemission
0255 to sometimes happen for fStopAndKill
0256
0257 But suspect coming after boundary may be causing
0258 the need for UseGivenVelocity_KLUDGE to get velocity and times correct see::
0259
0260     ~/o/notes/issues/Geant4_UseGivenVelocity_after_refraction_is_there_a_better_way_than_the_kludge_fix.rst
0261
0262
0263 **/
0264
0265
0266 void U4Physics::ConstructOp()
0267 {
0268     LOG(info) << desc() ;
0269
0270     if(Cerenkov_DISABLE == 0)
0271     {
0272         fCerenkov = new Local_G4Cerenkov_modified ;
0273         fCerenkov->SetMaxNumPhotonsPerStep(10000);
0274         fCerenkov->SetMaxBetaChangePerStep(10.0);
0275         fCerenkov->SetTrackSecondariesFirst(true);
0276         fCerenkov->SetVerboseLevel(EInt("Local_G4Cerenkov_modified_verboseLevel", "0"));
0277     }
0278
0279     if(Scintillation_DISABLE == 0)
0280     {
0281         fScintillation = new Local_DsG4Scintillation(EInt("Local_DsG4Scintillation_opticksMode","0")) ;
0282         fScintillation->SetTrackSecondariesFirst(true);
0283     }
0284
0285     if(FastSim_ENABLE == 1 )
0286     {
0287         fFastSim  = new G4FastSimulationManagerProcess("fast_sim_man");
0288     }
0289
0290     if(OpAbsorption_DISABLE == 0)
0291     {
0292 #ifdef DEBUG_TAG
0293         fAbsorption = new ShimG4OpAbsorption();
0294 #else
0295         fAbsorption = new G4OpAbsorption();
0296 #endif
0297     }
0298
0299     if(OpRayleigh_DISABLE == 0)
0300     {
0301 #ifdef DEBUG_TAG
0302         fRayleigh = new ShimG4OpRayleigh();
0303 #else
0304         fRayleigh = new G4OpRayleigh();
0305 #endif
0306     }
0307
0308     if(OpBoundaryProcess_DISABLE == 0)
0309     {
0310         fBoundary = CreateBoundaryProcess();
0311         LOG(info) << " fBoundary " << fBoundary ;
0312     }
0313
0314
0315
0316   auto particleIterator=GetParticleIterator();
0317   particleIterator->reset();
0318   while( (*particleIterator)() )
0319   {
0320         G4ParticleDefinition* particle = particleIterator->value();
0321         G4ProcessManager* pmanager = particle->GetProcessManager();
0322         G4String particleName = particle->GetParticleName();
0323
0324         if ( fCerenkov && fCerenkov->IsApplicable(*particle))
0325         {
0326             pmanager->AddProcess(fCerenkov);
0327             pmanager->SetProcessOrdering(fCerenkov,idxPostStep);
0328         }
0329
0330         if ( fScintillation && fScintillation->IsApplicable(*particle) && particleName != "opticalphoton")
0331         {
0332             pmanager->AddProcess(fScintillation);
0333             pmanager->SetProcessOrderingToLast(fScintillation, idxAtRest);
0334             pmanager->SetProcessOrderingToLast(fScintillation, idxPostStep);
0335         }
0336
0337         if (particleName == "opticalphoton")
0338         {
0339             ConstructOp_opticalphoton(pmanager, particleName);
0340         }
0341     }
0342 }
0343
0344
0345 /**
0346 U4Physics::ConstructOp_opticalphoton
0347 ------------------------------------
0348
0349 TO AVOID THE UseGivenVelocity KLUDGE scintillation process needs to be after absorption
0350 BUT UNFORTUNATELY PUTTING Scintillation AFTER Absorption prevents REEMISSION from happening
0351
0352 * ~/o/notes/issues/Geant4_UseGivenVelocity_KLUDGE_may_be_avoided_by_doing_PostStepDoIt_for_boundary_after_scintillation
0353 * ~/o/notes/issues/G4CXTest_GEOM_shakedown.rst
0354
0355 **/
0356
0357 void U4Physics::ConstructOp_opticalphoton(G4ProcessManager* pmanager, const G4String& particleName)
0358 {
0359     assert( particleName == "opticalphoton" );
0360
0361     if(fScintillation)
0362     {
0363         pmanager->AddProcess(fScintillation);
0364         pmanager->SetProcessOrderingToLast(fScintillation, idxAtRest);
0365         pmanager->SetProcessOrderingToLast(fScintillation, idxPostStep);
0366     }
0367     if(fAbsorption)    pmanager->AddDiscreteProcess(fAbsorption);
0368     if(fRayleigh)      pmanager->AddDiscreteProcess(fRayleigh);
0369     if(fBoundary)      pmanager->AddDiscreteProcess(fBoundary);
0370     if(fFastSim)       pmanager->AddDiscreteProcess(fFastSim);
0371 }
0372
0373
0374
0375
0376
0377
0378 /**
0379 U4Physics::CreateBoundaryProcess
0380 ---------------------------------
0381
0382 Looks like this needs updating now that it
0383 is normal to use WITH_CUSTOM4 within junosw+opticks
0384 without using WITH_PMTSIM
0385
0386 * NB : BOUNDARY CLASS CHANGES HERE MUST PARALLEL THOSE IN U4OpBoundaryProcess.h
0387
0388   * OTHERWISE GET "UNEXPECTED BoundaryFlag ZERO "
0389
0390 **/
0391
0392 G4VProcess* U4Physics::CreateBoundaryProcess()  // static
0393 {
0394     G4VProcess* proc = nullptr ;
0395
0396 #if defined(WITH_PMTSIM) && defined(WITH_CUSTOM4)
0397     const char* path = "$PMTSimParamData_BASE" ;  // directory with PMTSimParamData subfolder
0398     const PMTSimParamData* data = PMTAccessor::LoadData(path) ;
0399     LOG(LEVEL) << "load path "  << path << " giving PMTSimParamData.data: " << ( data ? "YES" : "NO" ) ;
0400     //LOG_IF(LEVEL, data != nullptr ) << *data ;
0401
0402     const PMTAccessor* pmt = PMTAccessor::Create(data) ;
0403     const C4IPMTAccessor* ipmt = pmt ;
0404     proc = new C4OpBoundaryProcess(ipmt);
0405
0406     LOG(LEVEL) << "create C4OpBoundaryProcess :  WITH_CUSTOM4 WITH_PMTSIM " ;
0407
0408 #elif defined(WITH_CUSTOM4)
0409     const char* jpmt = spath::Resolve("$CFBaseFromGEOM/CSGFoundry/SSim/extra/jpmt");
0410     const SPMTAccessor* pmt = SPMTAccessor::Load(jpmt) ;
0411     const char* geom = ssys::getenvvar("GEOM", "no-GEOM") ;
0412     LOG_IF(fatal, pmt == nullptr )
0413         << " FAILED TO SPMTAccessor::Load from [" << jpmt << "]"
0414         << " GEOM " << ( geom ? geom : "-" )
0415         ;
0416
0417     assert(pmt) ;  // trying to get C4 to work without the PMT info, just assert when really need PMT info
0418     const C4IPMTAccessor* ipmt = pmt ;
0419     proc = new C4OpBoundaryProcess(ipmt);
0420     LOG(LEVEL) << "create C4OpBoundaryProcess :  WITH_CUSTOM4 NOT:WITH_PMTSIM " ;
0421
0422 #elif defined(WITH_INSTRUMENTED_DEBUG)
0423     proc = new InstrumentedG4OpBoundaryProcess();
0424     LOG(LEVEL) << "create InstrumentedG4OpBoundaryProcess : NOT (WITH_PMTSIM and WITH_CUSTOM4) " ;
0425 #else
0426     proc = new G4OpBoundaryProcess();
0427     //LOG(LEVEL) << "create G4OpBoundaryProcess : NOT (WITH_PMTSIM and WITH_CUSTOM4) " ;
0428 #endif
0429     return proc ;
0430 }
0431