eic-opticks/u4/U4Genstep.cc

0001 #include "U4Genstep.h"
0002
0003 #include <string>
0004
0005 #include "G4Exception.hh"
0006 #include "G4MaterialPropertiesTable.hh"
0007 #include "G4ParticleDefinition.hh"
0008 #include "G4ProcessVector.hh"
0009 #include "G4Scintillation.hh"
0010 #include "G4SteppingManager.hh"
0011 #include "G4Track.hh"
0012
0013 U4Genstep::ProcessRange U4Genstep::Processes(G4SteppingManager* stepping_manager)
0014 {
0015     if (stepping_manager == nullptr)
0016         return {};
0017
0018     return Processes(
0019         stepping_manager->GetfStepStatus(),
0020         stepping_manager->GetfAtRestDoItVector(),
0021         stepping_manager->GetfPostStepDoItVector());
0022 }
0023
0024 U4Genstep::ProcessRange U4Genstep::Processes(
0025     G4StepStatus     status,
0026     G4ProcessVector* at_rest,
0027     G4ProcessVector* post_step)
0028 {
0029     G4ProcessVector* processes = status == fAtRestDoItProc ? at_rest : post_step;
0030     return {processes, processes ? processes->size() : 0};
0031 }
0032
0033 U4Genstep::ScintillationPlan U4Genstep::Scintillation(
0034     const G4Track*                   track,
0035     const G4Scintillation*           scintillation,
0036     const G4MaterialPropertiesTable* properties,
0037     G4int                            total_photons)
0038 {
0039     ScintillationPlan plan;
0040     if (scintillation == nullptr || properties == nullptr || total_photons <= 0)
0041         return plan;
0042
0043     const G4int component_keys[3] = {
0044         kSCINTILLATIONCOMPONENT1,
0045         kSCINTILLATIONCOMPONENT2,
0046         kSCINTILLATIONCOMPONENT3};
0047     plan.num_components = properties->GetProperty(component_keys[2])   ? 3
0048                           : properties->GetProperty(component_keys[1]) ? 2
0049                           : properties->GetProperty(component_keys[0]) ? 1
0050                                                                        : 0;
0051     if (plan.num_components == 0)
0052         return plan;
0053
0054     const char* standard_yield_keys[3] = {
0055         "SCINTILLATIONYIELD1",
0056         "SCINTILLATIONYIELD2",
0057         "SCINTILLATIONYIELD3"};
0058     const char* standard_time_keys[3] = {
0059         "SCINTILLATIONTIMECONSTANT1",
0060         "SCINTILLATIONTIMECONSTANT2",
0061         "SCINTILLATIONTIMECONSTANT3"};
0062
0063     G4String particle_prefix;
0064     if (scintillation->GetScintillationByParticleType())
0065     {
0066         if (track == nullptr)
0067             return {};
0068
0069         const G4ParticleDefinition* definition = track->GetParticleDefinition();
0070         const G4String&             name = definition->GetParticleName();
0071         particle_prefix = name == "proton" ? "PROTON" : name == "deuteron"                                            ? "DEUTERON"
0072                                                     : name == "triton"                                                ? "TRITON"
0073                                                     : name == "alpha"                                                 ? "ALPHA"
0074                                                     : name == "neutron" || definition->GetParticleType() == "nucleus" ? "ION"
0075                                                                                                                       : "ELECTRON";
0076     }
0077
0078     G4double yields[3] = {1., 0., 0.};
0079     for (G4int component = 0; component < plan.num_components; ++component)
0080     {
0081         const G4String suffix = std::to_string(component + 1);
0082         const G4String yield_key = particle_prefix.empty() ? G4String(standard_yield_keys[component])
0083                                                            : particle_prefix + "SCINTILLATIONYIELD" + suffix;
0084         const G4String time_key = particle_prefix.empty() ? G4String(standard_time_keys[component])
0085                                                           : particle_prefix + "SCINTILLATIONTIMECONSTANT" + suffix;
0086
0087         yields[component] = properties->ConstPropertyExists(yield_key) ? properties->GetConstProperty(yield_key)
0088                                                                        : (component == 0 ? 1. : 0.);
0089         plan.times[component] = properties->ConstPropertyExists(time_key) ? properties->GetConstProperty(time_key)
0090                                                                           : (properties->ConstPropertyExists(standard_time_keys[component]) ? properties->GetConstProperty(standard_time_keys[component])
0091                                                                                                                                             : 0.);
0092         plan.yield_sum += yields[component];
0093     }
0094
0095     if (plan.yield_sum <= 0.)
0096         return plan;
0097
0098     if (plan.num_components == 1)
0099     {
0100         plan.counts[0] = total_photons;
0101     }
0102     else if (plan.num_components == 2)
0103     {
0104         plan.counts[0] = G4int(yields[0] / plan.yield_sum * total_photons);
0105         plan.counts[1] = total_photons - plan.counts[0];
0106     }
0107     else
0108     {
0109         for (G4int component = 0; component < 3; ++component)
0110             plan.counts[component] =
0111                 G4int(yields[component] / plan.yield_sum * total_photons);
0112     }
0113
0114     for (G4int component = 0; component < plan.num_components; ++component)
0115     {
0116         if (properties->GetProperty(component_keys[component]) == nullptr)
0117             plan.counts[component] = 0;
0118     }
0119     return plan;
0120 }
0121
0122 bool U4Genstep::ValidateScintillation(
0123     const G4Scintillation* scintillation,
0124     const char*            origin)
0125 {
0126     if (scintillation == nullptr || !scintillation->GetFiniteRiseTime())
0127         return true;
0128
0129     G4ExceptionDescription description;
0130     description
0131         << "Finite scintillation rise time cannot be represented by "
0132         << "the current scintillation genstep. Disable "
0133         << "/process/optical/scintillation/setFiniteRiseTime before running.";
0134     G4Exception(
0135         origin ? origin : "U4Genstep::ValidateScintillation",
0136         "U4Genstep001",
0137         FatalException,
0138         description);
0139     return false;
0140 }