From 27ca8a3ccdd245ad2ac0105b99f752492e371681 Mon Sep 17 00:00:00 2001 From: Liam O'Sullivan Date: Tue, 17 Jun 2025 17:31:09 +0200 Subject: [PATCH] added changes from own edep-sim repo --- CMakeLists.txt | 4 +- app/edepSim.cc | 46 ++- io/EDepSimUnits.h | 326 ++++++++++++++++ io/TG4HitSegment.h | 7 +- src/CMakeLists.txt | 23 +- src/EDepSimBuilder.cc | 65 +++- src/EDepSimBuilder.hh | 32 +- src/EDepSimDetectorMessenger.cc | 123 +++++- src/EDepSimDetectorMessenger.hh | 9 +- src/EDepSimDokeBirks.cc | 32 +- src/EDepSimExtraPhysics.cc | 22 +- src/EDepSimHitSegment.cc | 106 +++--- src/EDepSimHitSegment.hh | 35 +- src/EDepSimPersistencyManager.cc | 21 +- src/EDepSimPhysicsList.cc | 1 + src/EDepSimRootGeometryManager.cc | 60 +-- src/EDepSimRootGeometryManager.hh | 20 +- src/EDepSimSegmentSD.cc | 5 +- src/EDepSimSegmentSD.hh | 24 +- src/EDepSimTrajectoryPoint.cc | 7 +- src/EDepSimTrajectoryPoint.hh | 18 + src/EDepSimUserDetectorConstruction.cc | 34 +- src/EDepSimUserPrimaryGeneratorMessenger.cc | 2 + src/NESTVersion098/G4S1Light.cc | 166 ++++---- src/NESTVersion098/G4S1Light.hh | 28 +- src/NESTVersion098/G4S2Light.cc | 1 - src/captain/CaptCryostatBuilder.cc | 100 +++-- src/captain/CaptExposedBuilder.cc | 16 +- src/captain/CaptWirePlaneBuilder.cc | 2 +- src/captain/CaptWorldBuilder.cc | 4 +- src/kinem/EDepSimGPSKinematicsGenerator.hh | 8 +- src/kinem/EDepSimHEPEVTKinematicsFactory.cc | 59 +++ src/kinem/EDepSimHEPEVTKinematicsFactory.hh | 63 +++ src/kinem/EDepSimHEPEVTKinematicsGenerator.cc | 360 ++++++++++++++++++ src/kinem/EDepSimHEPEVTKinematicsGenerator.hh | 69 ++++ .../EDepSimRooTrackerKinematicsGenerator.cc | 112 ++++-- .../EDepSimVConstrainedPositionGenerator.cc | 2 +- src/kinem/EDepSimVKinematicsFactory.hh | 6 +- src/kinem/EDepSimVPrimaryFactory.hh | 6 +- 39 files changed, 1637 insertions(+), 387 deletions(-) create mode 100644 io/EDepSimUnits.h create mode 100644 src/kinem/EDepSimHEPEVTKinematicsFactory.cc create mode 100644 src/kinem/EDepSimHEPEVTKinematicsFactory.hh create mode 100644 src/kinem/EDepSimHEPEVTKinematicsGenerator.cc create mode 100644 src/kinem/EDepSimHEPEVTKinematicsGenerator.hh diff --git a/CMakeLists.txt b/CMakeLists.txt index 3a68d57..c85f358 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,10 +9,12 @@ message("Energy Deposition Simulation -- ${VERSION}") # Define the options that can be set in the cache, or on the cmake # command line. set(CMAKE_BUILD_TYPE Debug) -set(EDEPSIM_DISPLAY TRUE CACHE BOOL +set(EDEPSIM_DISPLAY FALSE CACHE BOOL "If true, compile the edep-disp event display") set(EDEPSIM_READONLY FALSE CACHE BOOL "If true, then DO NOT use GEANT4") +set(EDEPSIM_USE_NEST TRUE CACHE BOOL + "If true, then make the NEST model available to be directly used") # Check to see if this is MACOS if(APPLE) diff --git a/app/edepSim.cc b/app/edepSim.cc index 881a2ed..30b0f76 100644 --- a/app/edepSim.cc +++ b/app/edepSim.cc @@ -36,10 +36,10 @@ void usage () { << std::endl; std::cout << " -v -- Increase the verbosity" << std::endl; std::cout << " -V =[quiet,log,info,verbose]" << std::endl - << " -- Change the named log level" + << " -- Change the named log level" << std::endl; std::cout << " -h -- This help message." << std::endl; - + exit(1); } @@ -66,10 +66,10 @@ int main(int argc,char** argv) { // integer or G4 will give an error. There could be error checking here, // but assume the user is "intelligent", and won't try to run "-e ten" // instead of "-e 10". - std::string beamOnCount = ""; - + std::string beamOnCount = ""; + if (argc<2) usage(); - + while (!errflg && ((c=getopt(argc,argv,"CdD:e:g:o:p:qsuUvV:h")) != -1)) { switch (c) { case 'C': { @@ -218,14 +218,14 @@ int main(int argc,char** argv) { EDepSim::LogManager::SetLogLevel(EDepSim::LogManager::VerboseLevel); EDepSimVerbose("Set log level to VerboseLevel"); } - - for (std::map::iterator i + + for (std::map::iterator i = namedLogLevel.begin(); i != namedLogLevel.end(); ++i) { EDepSim::LogManager::SetLogLevel(i->first.c_str(), i->second); } - + if (debugLevel == 1) { EDepSim::LogManager::SetDebugLevel(EDepSim::LogManager::WarnLevel); EDepSimWarn("Set debug level to WarnLevel"); @@ -239,32 +239,32 @@ int main(int argc,char** argv) { EDepSimTrace("Set debug level to TraceLevel"); } - for (std::map::iterator i + for (std::map::iterator i = namedDebugLevel.begin(); i != namedDebugLevel.end(); ++i) { EDepSim::LogManager::SetDebugLevel(i->first.c_str(), i->second); } - + // Set the mandatory initialization classes // Construct the default run manager G4RunManager* runManager = EDepSim::CreateRunManager(physicsList); - + // Create the persistency manager. The persistency manager must derive // from G4VPersistencyManager which will make this object available to the // G4RunManager as a singleton. There can only be one persistency manager // at a time. The Store methods will then be called by the run managers // Analyze methods. The persistency manager doesn't *have* to be derived - // from EDepSim::RootPersistencyManager, but a lot of the trajectory and hit - // handling functionality is handled by that class. + // from EDepSim::RootPersistencyManager, but it's probably the best way to + // save internal data structures. EDepSim::PersistencyManager* persistencyManager = NULL; persistencyManager = new EDepSim::RootPersistencyManager(); - + // Create a "no i/o" persistency manager. This doesn't actually save // anything, and it may be better to stop if there isn't a real // persistency manager declared. This is here as paranoia in case later - // gode gets clever about how the persistency manager is created. + // code gets clever about how the persistency manager is created. if (!persistencyManager) { persistencyManager = new EDepSim::PersistencyManager(); } @@ -278,7 +278,7 @@ int main(int argc,char** argv) { if (gdmlFilename != "") { UI->ApplyCommand("/edep/gdml/read "+gdmlFilename); } - + // Set the defaults for the simulation. UI->ApplyCommand("/edep/control edepsim-defaults 1.0"); @@ -286,7 +286,7 @@ int main(int argc,char** argv) { if (persistencyManager && ! outputFilename.empty()) { UI->ApplyCommand("/edep/db/open "+outputFilename); } - + std::signal(SIGILL, SIG_DFL); std::signal(SIGBUS, SIG_DFL); std::signal(SIGSEGV, SIG_DFL); @@ -295,7 +295,7 @@ int main(int argc,char** argv) { // generating the first event. This causes the executable to throw an // exception if the geometry has overlaps. if (validateGeometry) UI->ApplyCommand("/edep/validateGeometry"); - + // Set the random seed from the time. if (setSeed) UI->ApplyCommand("/edep/random/timeRandomSeed"); @@ -331,13 +331,19 @@ int main(int argc,char** argv) { UI->ApplyCommand("/run/beamOn " + beamOnCount); } } - + // If we have the persistency manager, then make sure it's closed. if (persistencyManager) { persistencyManager->Close(); delete persistencyManager; } delete runManager; - + return 0; } + +// Local Variables: +// mode:c++ +// c-basic-offset:4 +// compile-command:"$(git rev-parse --show-toplevel)/build/edep-build.sh force" +// End: diff --git a/io/EDepSimUnits.h b/io/EDepSimUnits.h new file mode 100644 index 0000000..83dd5ab --- /dev/null +++ b/io/EDepSimUnits.h @@ -0,0 +1,326 @@ +// -*- C++ -*- +// ---------------------------------------------------------------------- +// HEP coherent system of Units +// +// This file has been provided to CLHEP by Geant4 (simulation toolkit for +// HEP). And copied from the GEANT4 CLHEP into the EDepSim i/o library to +// specify the units being used in the i/o classes. These match the internal +// GEANT4 units. +// +// All of the definitions have been placed into the EDepSim namespace to avoid +// conflicts. The "unit" namespace can be used as an alternative by defining +// "USE_UNIT_NAMESPACE". +// +// The basic units are : +// millimeter (millimeter) +// nanosecond (nanosecond) +// Mega electron Volt (MeV) +// positron charge (eplus) +// degree Kelvin (kelvin) +// the amount of substance (mole) +// luminous intensity (candela) +// radian (radian) +// steradian (steradian) +// +// Below is a non exhaustive list of derived and pratical units +// (i.e. mostly the SI units). +// You can add your own units. +// +// The SI numerical value of the positron charge is defined here, +// as it is needed for conversion factor : positron charge = e_SI (coulomb) +// +// The others physical constants are defined in the header file : +// PhysicalConstants.h +// +// Authors: M.Maire, S.Giani +// +// History: +// +// 06.02.96 Created. +// 28.03.96 Added miscellaneous constants. +// 05.12.97 E.Tcherniaev: Redefined pascal (to avoid warnings on WinNT) +// 20.05.98 names: meter, second, gram, radian, degree +// (from Brian.Lasiuk@yale.edu (STAR)). Added luminous units. +// 05.08.98 angstrom, picobarn, microsecond, picosecond, petaelectronvolt +// 01.03.01 parsec +// 31.01.06 kilogray, milligray, microgray +// 29.04.08 use PDG 2006 value of e_SI +// 03.11.08 use PDG 2008 value of e_SI +// 19.08.15 added liter and its sub units (mma) +// 12.01.16 added symbols for microsecond (us) and picosecond (ps) (mma) + +#ifndef EDepSimUnits_h_Seen +#define EDepSimUnits_h_Seen + +#ifdef USE_UNIT_NAMESPACE +namespace unit { +#else +namespace EDepSim { +#endif + // + // + // + static constexpr double pi = 3.14159265358979323846; + static constexpr double twopi = 2*pi; + static constexpr double halfpi = pi/2; + static constexpr double pi2 = pi*pi; + + // + // Length [L] + // + static constexpr double millimeter = 1.; + static constexpr double millimeter2 = millimeter*millimeter; + static constexpr double millimeter3 = millimeter*millimeter*millimeter; + + static constexpr double centimeter = 10.*millimeter; + static constexpr double centimeter2 = centimeter*centimeter; + static constexpr double centimeter3 = centimeter*centimeter*centimeter; + + static constexpr double meter = 1000.*millimeter; + static constexpr double meter2 = meter*meter; + static constexpr double meter3 = meter*meter*meter; + + static constexpr double kilometer = 1000.*meter; + static constexpr double kilometer2 = kilometer*kilometer; + static constexpr double kilometer3 = kilometer*kilometer*kilometer; + + static constexpr double parsec = 3.0856775807e+16*meter; + + static constexpr double micrometer = 1.e-6 *meter; + static constexpr double nanometer = 1.e-9 *meter; + static constexpr double angstrom = 1.e-10*meter; + static constexpr double fermi = 1.e-15*meter; + + static constexpr double barn = 1.e-28*meter2; + static constexpr double millibarn = 1.e-3 *barn; + static constexpr double microbarn = 1.e-6 *barn; + static constexpr double nanobarn = 1.e-9 *barn; + static constexpr double picobarn = 1.e-12*barn; + + // symbols + static constexpr double nm = nanometer; + static constexpr double um = micrometer; + + static constexpr double mm = millimeter; + static constexpr double mm2 = millimeter2; + static constexpr double mm3 = millimeter3; + + static constexpr double cm = centimeter; + static constexpr double cm2 = centimeter2; + static constexpr double cm3 = centimeter3; + + static constexpr double liter = 1.e+3*cm3; + static constexpr double L = liter; + static constexpr double dL = 1.e-1*liter; + static constexpr double cL = 1.e-2*liter; + static constexpr double mL = 1.e-3*liter; + + static constexpr double m = meter; + static constexpr double m2 = meter2; + static constexpr double m3 = meter3; + + static constexpr double km = kilometer; + static constexpr double km2 = kilometer2; + static constexpr double km3 = kilometer3; + + static constexpr double pc = parsec; + + // + // Angle + // + static constexpr double radian = 1.; + static constexpr double milliradian = 1.e-3*radian; + static constexpr double degree = (pi/180.0)*radian; + + static constexpr double steradian = 1.; + + // symbols + static constexpr double rad = radian; + static constexpr double mrad = milliradian; + static constexpr double sr = steradian; + static constexpr double deg = degree; + + // + // Time [T] + // + static constexpr double nanosecond = 1.; + static constexpr double second = 1.e+9 *nanosecond; + static constexpr double millisecond = 1.e-3 *second; + static constexpr double microsecond = 1.e-6 *second; + static constexpr double picosecond = 1.e-12*second; + + static constexpr double hertz = 1./second; + static constexpr double kilohertz = 1.e+3*hertz; + static constexpr double megahertz = 1.e+6*hertz; + + // symbols + static constexpr double ns = nanosecond; + static constexpr double s = second; + static constexpr double ms = millisecond; + static constexpr double us = microsecond; + static constexpr double ps = picosecond; + + // + // Electric charge [Q] + // + static constexpr double eplus = 1. ;// positron charge + static constexpr double e_SI = 1.602176487e-19;// positron charge in coulomb + static constexpr double coulomb = eplus/e_SI;// coulomb = 6.24150 e+18 * eplus + + // + // Energy [E] + // + static constexpr double megaelectronvolt = 1. ; + static constexpr double electronvolt = 1.e-6*megaelectronvolt; + static constexpr double kiloelectronvolt = 1.e-3*megaelectronvolt; + static constexpr double gigaelectronvolt = 1.e+3*megaelectronvolt; + static constexpr double teraelectronvolt = 1.e+6*megaelectronvolt; + static constexpr double petaelectronvolt = 1.e+9*megaelectronvolt; + + static constexpr double joule = electronvolt/e_SI;// joule = 6.24150 e+12 * MeV + + // symbols + static constexpr double MeV = megaelectronvolt; + static constexpr double eV = electronvolt; + static constexpr double keV = kiloelectronvolt; + static constexpr double GeV = gigaelectronvolt; + static constexpr double TeV = teraelectronvolt; + static constexpr double PeV = petaelectronvolt; + + // + // Mass [E][T^2][L^-2] + // + static constexpr double kilogram = joule*second*second/(meter*meter); + static constexpr double gram = 1.e-3*kilogram; + static constexpr double milligram = 1.e-3*gram; + + // symbols + static constexpr double kg = kilogram; + static constexpr double g = gram; + static constexpr double mg = milligram; + + // + // Power [E][T^-1] + // + static constexpr double watt = joule/second;// watt = 6.24150 e+3 * MeV/ns + + // + // Force [E][L^-1] + // + static constexpr double newton = joule/meter;// newton = 6.24150 e+9 * MeV/mm + + // + // Pressure [E][L^-3] + // + static constexpr double pascal = newton/m2; // pascal = 6.24150 e+3 * MeV/mm3 + static constexpr double bar = 100000*pascal; // bar = 6.24150 e+8 * MeV/mm3 + static constexpr double atmosphere = 101325*pascal; // atm = 6.32420 e+8 * MeV/mm3 + + // + // Electric current [Q][T^-1] + // + static constexpr double ampere = coulomb/second; // ampere = 6.24150 e+9 * eplus/ns + static constexpr double milliampere = 1.e-3*ampere; + static constexpr double microampere = 1.e-6*ampere; + static constexpr double nanoampere = 1.e-9*ampere; + + // + // Electric potential [E][Q^-1] + // + static constexpr double megavolt = megaelectronvolt/eplus; + static constexpr double kilovolt = 1.e-3*megavolt; + static constexpr double volt = 1.e-6*megavolt; + + // + // Electric resistance [E][T][Q^-2] + // + static constexpr double ohm = volt/ampere;// ohm = 1.60217e-16*(MeV/eplus)/(eplus/ns) + + // + // Electric capacitance [Q^2][E^-1] + // + static constexpr double farad = coulomb/volt;// farad = 6.24150e+24 * eplus/Megavolt + static constexpr double millifarad = 1.e-3*farad; + static constexpr double microfarad = 1.e-6*farad; + static constexpr double nanofarad = 1.e-9*farad; + static constexpr double picofarad = 1.e-12*farad; + + // + // Magnetic Flux [T][E][Q^-1] + // + static constexpr double weber = volt*second;// weber = 1000*megavolt*ns + + // + // Magnetic Field [T][E][Q^-1][L^-2] + // + static constexpr double tesla = volt*second/meter2;// tesla =0.001*megavolt*ns/mm2 + + static constexpr double gauss = 1.e-4*tesla; + static constexpr double kilogauss = 1.e-1*tesla; + + // + // Inductance [T^2][E][Q^-2] + // + static constexpr double henry = weber/ampere;// henry = 1.60217e-7*MeV*(ns/eplus)**2 + + // + // Temperature + // + static constexpr double kelvin = 1.; + + // + // Amount of substance + // + static constexpr double mole = 1.; + + // + // Activity [T^-1] + // + static constexpr double becquerel = 1./second ; + static constexpr double curie = 3.7e+10 * becquerel; + static constexpr double kilobecquerel = 1.e+3*becquerel; + static constexpr double megabecquerel = 1.e+6*becquerel; + static constexpr double gigabecquerel = 1.e+9*becquerel; + static constexpr double millicurie = 1.e-3*curie; + static constexpr double microcurie = 1.e-6*curie; + static constexpr double Bq = becquerel; + static constexpr double kBq = kilobecquerel; + static constexpr double MBq = megabecquerel; + static constexpr double GBq = gigabecquerel; + static constexpr double Ci = curie; + static constexpr double mCi = millicurie; + static constexpr double uCi = microcurie; + + // + // Absorbed dose [L^2][T^-2] + // + static constexpr double gray = joule/kilogram ; + static constexpr double kilogray = 1.e+3*gray; + static constexpr double milligray = 1.e-3*gray; + static constexpr double microgray = 1.e-6*gray; + + // + // Luminous intensity [I] + // + static constexpr double candela = 1.; + + // + // Luminous flux [I] + // + static constexpr double lumen = candela*steradian; + + // + // Illuminance [I][L^-2] + // + static constexpr double lux = lumen/meter2; + + // + // Miscellaneous + // + static constexpr double perCent = 0.01 ; + static constexpr double perThousand = 0.001; + static constexpr double perMillion = 0.000001; + +} // namespace CLHEP + +#endif /* EDepSimUnits_h_Seen */ diff --git a/io/TG4HitSegment.h b/io/TG4HitSegment.h index eb42010..06bde3a 100644 --- a/io/TG4HitSegment.h +++ b/io/TG4HitSegment.h @@ -67,12 +67,15 @@ class TG4HitSegment : public TObject { /// The TrackId for each trajectory that contributed to this hit. This /// could contain the TrackId of the primary particle, but not /// necessarily. - Contributors Contrib; + const Contributors& GetContributors() const {return Contrib;} // The public fields are deprecated but still supported by default in the // current version. #define EDEPSIM_USE_PUBLIC_FIELDS - + + // Not hidden in private since it's the most commonly used interface. + Contributors Contrib; + #if defined(EDEPSIM_USE_PUBLIC_FIELDS)&&!defined(EDEPSIM_FORCE_PRIVATE_FIELDS)&&!defined(__CINT__) public: #ifdef EDEPSIM_WARN_PUBLIC_FIELDS diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e6f60e1..2accb25 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,8 +2,25 @@ message("Energy Deposition Simulation") # Suck up all of the .cc files for the source. This isn't good CMAKE # practice, but it makes maintaining this file easier. -file(GLOB_RECURSE source *.cc) -file(GLOB_RECURSE includes EDepSim*.hh) +# file(GLOB_RECURSE source *.cc) +# file(GLOB_RECURSE includes EDepSim*.hh) +file(GLOB edep_source *.cc) +file(GLOB edep_includes *.cc) +file(GLOB captain_source captain/*.cc) +file(GLOB captain_includes captain/EDepSim*.hh) +file(GLOB kinem_source kinem/*.cc) +file(GLOB kinem_includes kinem/EDepSim*.hh) + +if(EDEPSIM_USE_NEST) + file(GLOB nest_source NESTVersion098/*.cc) + file(GLOB nest_includes NESTVersion098/EDepSim*.hh) +else(EDEPSIM_USE_NEST) + message("Not including NEST in the compilation") + add_definitions(-DEDEPSIM_SKIP_NESTVersion098) +endif(EDEPSIM_USE_NEST) + +set(source ${edep_source} ${captain_source} ${kinem_source} ${nest_source}) +set(includes ${edep_includes} ${captain_includes} ${kinem_includes}) # Where the macro files can be found. add_definitions(-DDETSIM_INSTALL_PREFIX=\"${CMAKE_INSTALL_PREFIX}\") @@ -34,5 +51,3 @@ install(TARGETS edepsim # Install the header files. install(FILES ${includes} DESTINATION include/EDepSim) - - diff --git a/src/EDepSimBuilder.cc b/src/EDepSimBuilder.cc index 9c9ac1f..0417827 100644 --- a/src/EDepSimBuilder.cc +++ b/src/EDepSimBuilder.cc @@ -8,17 +8,17 @@ #include #include -EDepSim::Builder::Builder(G4String n, - EDepSim::UserDetectorConstruction* c) +EDepSim::Builder::Builder(G4String n, + EDepSim::UserDetectorConstruction* c) : fLocalName(n), fName(n), fConstruction(c), fParent(NULL), - fMessenger(NULL), fSensitiveDetector(NULL), + fMessenger(NULL), fSensitiveDetector(NULL), fOpacity(0.0), fCheck(false) { fMessenger = fConstruction->GetMessenger(); } -EDepSim::Builder::Builder(G4String n, EDepSim::Builder* p) +EDepSim::Builder::Builder(G4String n, EDepSim::Builder* p) : fLocalName(n), fName(n), fConstruction(NULL), fParent(p), - fMessenger(NULL), fSensitiveDetector(NULL), + fMessenger(NULL), fSensitiveDetector(NULL), fOpacity(0.0), fCheck(false) { fName = fParent->GetName() + "/" + fLocalName; fConstruction = fParent->GetConstruction(); @@ -63,7 +63,7 @@ void EDepSim::Builder::SetOpacity(double v) { /// Set the visibility of the constructed object. void EDepSim::Builder::SetDaughterOpacity(double v) { - for (std::map::iterator p + for (std::map::iterator p = fSubBuilders.begin(); p!=fSubBuilders.end(); ++p) { @@ -84,14 +84,14 @@ EDepSim::BuilderMessenger::BuilderMessenger(EDepSim::Builder* c, std::string guidance = "Commands for the " + c->GetName() + " Builder"; fDirectory->SetGuidance(guidance.c_str()); } - + fOpacityCMD = new G4UIcmdWithADouble(CommandName("opacity"),this); fOpacityCMD->SetGuidance("Set the log of the relative opacity." " Useful values are between [-0.5, 0.5], and" " +/-10 makes the object opague or transparent." " A value of 0 leaves the opacity unchanged."); fOpacityCMD->SetParameterName("opacity",false); - + fDaughterOpacityCMD = new G4UIcmdWithADouble( CommandName("daughterOpacity"),this); fDaughterOpacityCMD->SetGuidance( @@ -110,14 +110,20 @@ EDepSim::BuilderMessenger::BuilderMessenger(EDepSim::Builder* c, G4UIparameter* typeParam = new G4UIparameter('s'); typeParam->SetParameterName("Type"); fSensitiveCMD->SetParameter(typeParam); - - fMaximumHitSagittaCMD + + fMaximumHitSagittaCMD = new G4UIcmdWithADoubleAndUnit(CommandName("maxHitSagitta"),this); fMaximumHitSagittaCMD->SetGuidance("Set the maximum sagitta for a hit."); fMaximumHitSagittaCMD->SetParameterName("Sagitta",false); fMaximumHitSagittaCMD->SetUnitCategory("Length"); - - fMaximumHitLengthCMD + + fMaximumHitSeparationCMD + = new G4UIcmdWithADoubleAndUnit(CommandName("maxHitSeparation"),this); + fMaximumHitSeparationCMD->SetGuidance("Set the maximum separation for a hit."); + fMaximumHitSeparationCMD->SetParameterName("Separation",false); + fMaximumHitSeparationCMD->SetUnitCategory("Length"); + + fMaximumHitLengthCMD = new G4UIcmdWithADoubleAndUnit(CommandName("maxHitLength"),this); fMaximumHitLengthCMD->SetGuidance("Set the maximum length for a hit."); fMaximumHitLengthCMD->SetParameterName("HitLength",false); @@ -147,6 +153,11 @@ void EDepSim::BuilderMessenger::SetNewValue(G4UIcommand *cmd, G4String val) { SetMaximumHitSagitta( fMaximumHitSagittaCMD->GetNewDoubleValue(val)); } + else if (cmd==fMaximumHitSeparationCMD) { + fBuilder-> + SetMaximumHitSeparation( + fMaximumHitSeparationCMD->GetNewDoubleValue(val)); + } else if (cmd==fMaximumHitLengthCMD) { fBuilder-> SetMaximumHitLength( @@ -161,15 +172,16 @@ EDepSim::BuilderMessenger::~BuilderMessenger() { delete fCheckCMD; delete fSensitiveCMD; delete fMaximumHitSagittaCMD; + delete fMaximumHitSeparationCMD; delete fMaximumHitLengthCMD; } G4VisAttributes EDepSim::Builder::GetColor(G4Material* mat, double opacity) { EDepSim::RootGeometryManager* geoMan = EDepSim::RootGeometryManager::Get(); // Check if the opacity has been over-ruled by the macro... - if (fOpacity < -9.9) return G4VisAttributes::Invisible; + if (fOpacity < -9.9) return G4VisAttributes::GetInvisible(); opacity += fOpacity; - if (opacity < -9.9) return G4VisAttributes::Invisible; + if (opacity < -9.9) return G4VisAttributes::GetInvisible(); G4Color color = geoMan->GetG4Color(mat); double red = color.GetRed(); double green = color.GetGreen(); @@ -191,15 +203,16 @@ void EDepSim::Builder::SetSensitiveDetector(G4String name, G4String type) { EDepSim::SDFactory factory(type); SetSensitiveDetector(factory.MakeSD(name)); } - + void EDepSim::Builder::SetMaximumHitSagitta(double sagitta) { if (!fSensitiveDetector) { EDepSimError("Maximum hit sagitta must be set after the sensitive" - " detector is defined."); + " detector is defined."); EDepSimThrow("Builder does not have sensitive detector defined"); return; } - EDepSim::SegmentSD *segSD = dynamic_cast(fSensitiveDetector); + EDepSim::SegmentSD *segSD + = dynamic_cast(fSensitiveDetector); if (!segSD) { EDepSimThrow("Sensitive detector not derived from EDepSim::SegmentSD"); return; @@ -207,6 +220,22 @@ void EDepSim::Builder::SetMaximumHitSagitta(double sagitta) { segSD->SetMaximumHitSagitta(sagitta); } +void EDepSim::Builder::SetMaximumHitSeparation(double separation) { + if (!fSensitiveDetector) { + EDepSimError("Maximum hit separation must be set after the sensitive" + " detector is defined."); + EDepSimThrow("Builder does not have sensitive detector defined"); + return; + } + EDepSim::SegmentSD *segSD + = dynamic_cast(fSensitiveDetector); + if (!segSD) { + EDepSimThrow("Sensitive detector not derived from EDepSim::SegmentSD"); + return; + } + segSD->SetMaximumHitSeparation(separation); +} + void EDepSim::Builder::SetMaximumHitLength(double length) { if (!fSensitiveDetector) { EDepSimError("Maximum hit length must be set after the sensitive" @@ -221,5 +250,3 @@ void EDepSim::Builder::SetMaximumHitLength(double length) { } segSD->SetMaximumHitLength(length); } - - diff --git a/src/EDepSimBuilder.hh b/src/EDepSimBuilder.hh index fdc2dee..e8fa179 100644 --- a/src/EDepSimBuilder.hh +++ b/src/EDepSimBuilder.hh @@ -66,23 +66,23 @@ public: /// Set the relative opacity of the constructed object. void SetOpacity(double v); - + /// Get the relative opacity of the constructed object. double GetOpacity(void) const {return fOpacity;} /// Set the relative opacity of the object daughters void SetDaughterOpacity(double v); - + /// Return the detector construction that is building this piece. EDepSim::UserDetectorConstruction* GetConstruction(void) { return fConstruction; }; - + /// Set the sensitive detector for this component. virtual void SetSensitiveDetector(G4VSensitiveDetector* s) { fSensitiveDetector = s; } - + /// Get the sensitive detector for this component. virtual G4VSensitiveDetector* GetSensitiveDetector(void) { return fSensitiveDetector; @@ -90,12 +90,17 @@ public: /// Set the sensitive detector for this component by name. virtual void SetSensitiveDetector(G4String name, G4String type); - + /// Set the maximum sagitta for a track being added to a single hit - /// segment. + /// segment. This is passed to the sensitive detector. virtual void SetMaximumHitSagitta(double sagitta); - /// Set the maximum length of a single hit segment. + /// Set the maximum separation between deposits being added to a single + /// hit segment. This is passed to the sensitive detector. + virtual void SetMaximumHitSeparation(double separation); + + /// Set the maximum length of a single hit segment. This is passed to the + /// sensitive detector. virtual void SetMaximumHitLength(double length); /// Return the messenger for this constructor @@ -111,9 +116,9 @@ public: /// that will be used by the derived class should be added using this /// method. void AddBuilder(EDepSim::Builder* c) { - if (fSubBuilders.find(c->GetLocalName()) + if (fSubBuilders.find(c->GetLocalName()) != fSubBuilders.end()) { - std::cout << "Multiply defined constructor name " + std::cout << "Multiply defined constructor name " << c->GetName() << std::endl; EDepSimThrow("EDepSim::Builder::AddBuilder(): Multiple defines"); @@ -124,7 +129,7 @@ public: /// Get a sub constructor by name and do the dynamic cast. This will /// abort with an error message if you request an undefined name. template T& Get(G4String n) { - std::map::iterator p + std::map::iterator p = fSubBuilders.find(n); if (p==fSubBuilders.end()) { std::cout << "Error in " << GetName() << std::endl; @@ -143,7 +148,7 @@ public: T* c = dynamic_cast((*p).second); if (!c) { std::cout << "Error in " << GetName() << std::endl; - std::cout << " Cannot type cast " << n + std::cout << " Cannot type cast " << n << " to requested class" << std::endl; EDepSimThrow("EDepSim::Builder::Get<>:" " Bad typecast"); @@ -154,7 +159,7 @@ public: /// Find a sub constructor by name and do the dynamic cast. This returns /// a pointer that will be NULL if you request an undefined name. template T* Find(G4String n) { - std::map::iterator p + std::map::iterator p = fSubBuilders.find(n); if (p==fSubBuilders.end()) return NULL; T* c = dynamic_cast((*p).second); @@ -203,7 +208,7 @@ private: /// The parent of this constructor EDepSim::Builder* fParent; - /// The user interface messenger that will provide values for this class. + /// The user interface messenger that will provide values for this class. G4UImessenger* fMessenger; /// The sensitive detector for this tracking component. @@ -236,6 +241,7 @@ private: G4UIcmdWithABool* fCheckCMD; G4UIcommand* fSensitiveCMD; G4UIcmdWithADoubleAndUnit* fMaximumHitSagittaCMD; + G4UIcmdWithADoubleAndUnit* fMaximumHitSeparationCMD; G4UIcmdWithADoubleAndUnit* fMaximumHitLengthCMD; public: diff --git a/src/EDepSimDetectorMessenger.cc b/src/EDepSimDetectorMessenger.cc index d72cefb..9be7d28 100644 --- a/src/EDepSimDetectorMessenger.cc +++ b/src/EDepSimDetectorMessenger.cc @@ -55,13 +55,6 @@ EDepSim::DetectorMessenger::DetectorMessenger(EDepSim::UserDetectorConstruction* "generation, or reading of kinematics files."); fExportCmd->SetParameterName("RootFile", false); - fGDMLDir = new G4UIdirectory("/edep/gdml/"); - fGDMLDir->SetGuidance("Control over the gdml file."); - - fGDMLCmd = new G4UIcmdWithAString("/edep/gdml/read",this); - fGDMLCmd->SetGuidance("Define a GDML file to be used for the geometry."); - fGDMLCmd->AvailableForStates(G4State_PreInit); - fControlCmd = new G4UIcommand("/edep/control",this); fControlCmd->SetGuidance( "Set the run conditions for this simulation. This takes the name\n" @@ -97,6 +90,24 @@ EDepSim::DetectorMessenger::DetectorMessenger(EDepSim::UserDetectorConstruction* par = new G4UIparameter("Unit", 's', false); fHitSagittaCmd->SetParameter(par); + fHitSeparationCmd = new G4UIcommand("/edep/hitSeparation",this); + fHitSeparationCmd->SetGuidance( + "Set the maximum separation for hit segments."); + fHitSeparationCmd->AvailableForStates(G4State_PreInit); + + // The name of the sensitive detector. + par = new G4UIparameter("Sensitive", 's', false); + fHitSeparationCmd->SetParameter(par); + + // The name of the sensitive detector. + par = new G4UIparameter("Separation", 'd', false); + fHitSeparationCmd->SetParameter(par); + + // The name of the sensitive detector. + par = new G4UIparameter("Unit", 's', false); + fHitSeparationCmd->SetParameter(par); + + //////////// fHitLengthCmd = new G4UIcommand("/edep/hitLength",this); fHitLengthCmd->SetGuidance( "Set the length for hit segments to stop growing."); @@ -123,6 +134,39 @@ EDepSim::DetectorMessenger::DetectorMessenger(EDepSim::UserDetectorConstruction* par = new G4UIparameter("LogicalVolume", 's', false); fHitExcludedCmd->SetParameter(par); + ///////////////////////////////////////////////////////////// + // Add commands for the /edep/gdml/ sub-directory + fGDMLDir = new G4UIdirectory("/edep/gdml/"); + fGDMLDir->SetGuidance("Control over the gdml file."); + + fGDMLReadCmd = new G4UIcmdWithAString("/edep/gdml/read",this); + fGDMLReadCmd->SetGuidance("Set a GDML file to be used for the geometry."); + fGDMLReadCmd->AvailableForStates(G4State_PreInit); + + ///////////////////////////////////////////////////////////// + // Add commands for the /edep/material/ sub-directory + // + // These provide ways to inject information into the materials when it's + // not supported by the geometry creation tools. + fMaterialDir = new G4UIdirectory("/edep/material/"); + fMaterialDir->SetGuidance("Extra control over the material definitions."); + + fMaterialBirksCmd = new G4UIcommand("/edep/material/birksConstant",this); + fMaterialBirksCmd->SetGuidance("Override the Birks constant for material."); + fMaterialBirksCmd->AvailableForStates(G4State_PreInit); + + // The name of the material get a Birks constant. + par = new G4UIparameter("Material", 's', false); + fMaterialBirksCmd->SetParameter(par); + + // The value of the Birks constant. + par = new G4UIparameter("Birks", 'd', false); + fMaterialBirksCmd->SetParameter(par); + + // The unit for the Birks constant. + par = new G4UIparameter("Unit", 's', false); + fMaterialBirksCmd->SetParameter(par); + } @@ -132,13 +176,16 @@ EDepSim::DetectorMessenger::~DetectorMessenger() delete fPrintMassCmd; delete fValidateCmd; delete fExportCmd; - delete fGDMLCmd; delete fControlCmd; delete fHitSagittaCmd; + delete fHitSeparationCmd; delete fHitLengthCmd; delete fHitExcludedCmd; - delete fEDepSimDir; + delete fGDMLReadCmd; delete fGDMLDir; + delete fMaterialBirksCmd; + delete fMaterialDir; + delete fEDepSimDir; } @@ -157,7 +204,7 @@ void EDepSim::DetectorMessenger::SetNewValue(G4UIcommand * cmd, else if (cmd == fExportCmd) { EDepSim::RootGeometryManager::Get()->Export(newValue); } - else if (cmd == fGDMLCmd) { + else if (cmd == fGDMLReadCmd) { EDepSimLog("Read a gdml geometry from |" << newValue << "|"); G4GDMLParser* gdmlParser = fConstruction->GetGDMLParser(); if (gdmlParser) delete gdmlParser; @@ -209,6 +256,22 @@ void EDepSim::DetectorMessenger::SetNewValue(G4UIcommand * cmd, std::cout << "Invalid sensitive detector" << std::endl; } } + else if (cmd == fHitSeparationCmd) { + std::istringstream input((const char*)newValue); + std::string sdName; + double separation; + std::string unitName; + input >> sdName >> separation >> unitName; + separation *= G4UnitDefinition::GetValueOf(unitName); + SDFactory factory("segment"); + SegmentSD* sd = dynamic_cast(factory.MakeSD(sdName)); + if (sd) { + sd->SetMaximumHitSeparation(separation); + } + else { + std::cout << "Invalid sensitive detector" << std::endl; + } + } else if (cmd == fHitLengthCmd) { std::istringstream input((const char*)newValue); std::string sdName; @@ -231,5 +294,45 @@ void EDepSim::DetectorMessenger::SetNewValue(G4UIcommand * cmd, input >> logName; fConstruction->AddExcludedSensitiveDetector(logName); } + else if (cmd == fMaterialBirksCmd) { + std::istringstream input((const char*)newValue); + std::string matName; + double birksConstant; + std::string unitName; + input >> matName >> birksConstant >> unitName; + // Get the material that needs to be updated. + G4Material *material = G4Material::GetMaterial(matName.c_str()); + if (!material) { + EDepSimError("Unknown material for /edep/material/SetBirksConstant:" + << " " << matName); + throw std::runtime_error("Unknown material"); + } + + // parse the unitName which should be [length]/[energy] + std::size_t divide = unitName.find("/"); + if (divide == std::string::npos) { + EDepSimError("Illegal unit name for" + << " /edep/material/SetBirksConstant: " + << unitName); + throw std::runtime_error("Illegal input"); + } + std::string numerator = unitName.substr(0,divide); + std::string denominator = unitName.substr(divide+1); + birksConstant *= G4UnitDefinition::GetValueOf(numerator); + birksConstant /= G4UnitDefinition::GetValueOf(denominator); + + double bc = material->GetIonisation()->GetBirksConstant(); + if (bc > 1E-6) { + EDepSimError("Overriding Birks constant for " << matName + << " from nonzero value of" + << " " << bc/(mm/MeV) << " mm/MeV"); + } + + material->GetIonisation()->SetBirksConstant(birksConstant); + + EDepSimLog("Set Birks constant for" + << " " << matName + << " to " << birksConstant/(mm/MeV) << " mm/MeV"); + } } diff --git a/src/EDepSimDetectorMessenger.hh b/src/EDepSimDetectorMessenger.hh index 09878c2..11a6f29 100644 --- a/src/EDepSimDetectorMessenger.hh +++ b/src/EDepSimDetectorMessenger.hh @@ -25,18 +25,23 @@ private: EDepSim::UserDetectorConstruction* fConstruction; G4UIdirectory* fEDepSimDir; - G4UIdirectory* fGDMLDir; G4UIcmdWithoutParameter* fUpdateCmd; G4UIcmdWithAString* fPrintMassCmd; G4UIcmdWithoutParameter* fValidateCmd; G4UIcmdWithAString* fExportCmd; - G4UIcmdWithAString* fGDMLCmd; G4UIcommand* fControlCmd; G4UIcommand* fHitSagittaCmd; + G4UIcommand* fHitSeparationCmd; G4UIcommand* fHitLengthCmd; G4UIcommand* fHitExcludedCmd; + G4UIdirectory* fGDMLDir; + G4UIcmdWithAString* fGDMLReadCmd; + + G4UIdirectory* fMaterialDir; + G4UIcommand* fMaterialBirksCmd; + }; #endif diff --git a/src/EDepSimDokeBirks.cc b/src/EDepSimDokeBirks.cc index a9e2cc1..eea7ef0 100644 --- a/src/EDepSimDokeBirks.cc +++ b/src/EDepSimDokeBirks.cc @@ -21,7 +21,7 @@ EDepSim::DokeBirks::DokeBirks(const G4String& processName, fEmSaturation(nullptr) { SetProcessSubType(fScintillation); - + if (verboseLevel>0) { G4cout << GetProcessName() << " is created " << G4endl; } @@ -30,7 +30,7 @@ EDepSim::DokeBirks::DokeBirks(const G4String& processName, EDepSim::DokeBirks::~DokeBirks () {} //destructor needed to avoid linker error - + G4bool EDepSim::DokeBirks::IsApplicable(const G4ParticleDefinition& aParticleType) { @@ -60,7 +60,7 @@ G4VParticleChange* EDepSim::DokeBirks::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) { aParticleChange.Initialize(aTrack); - + const G4Material* aMaterial = aTrack.GetMaterial(); bool inLiquidArgon = true; @@ -81,21 +81,21 @@ EDepSim::DokeBirks::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) // Do Nothing inLiquidArgon = false; } - + const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); const G4ParticleDefinition *pDef = aParticle->GetDefinition(); G4String particleName = pDef->GetParticleName(); // G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint(); G4double totalEnergyDeposit = aStep.GetTotalEnergyDeposit(); - + if (totalEnergyDeposit <= 0.0) { - return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); } if (aStep.GetStepLength() <= 0.0) { - return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); } - + if (!inLiquidArgon) { if (!fEmSaturation) fEmSaturation = new G4EmSaturation(0); G4double nonIonizingEnergy @@ -120,7 +120,7 @@ EDepSim::DokeBirks::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) EDepSimError("LAr Volume " << aLogVolume->GetName() << " should have field!"); - return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); } G4ThreeVector thePosition = pPostStepPoint->GetPosition(); @@ -129,7 +129,7 @@ EDepSim::DokeBirks::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) thePosition.x(), thePosition.y(), thePosition.z(), theTime}; double bField[6]; aField->GetFieldValue(point,bField); - + double electricField = bField[3]*bField[3]; electricField += bField[4]*bField[4]; electricField += bField[5]*bField[5]; @@ -180,7 +180,7 @@ EDepSim::DokeBirks::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) LET *= ratio; dx /= ratio; } while (false); - + G4double recombProb = (dokeBirks[0]*LET)/(1+dokeBirks[1]*LET)+dokeBirks[2]; G4double ScintillationYield = 1 / (19.5*eV); @@ -196,7 +196,7 @@ EDepSim::DokeBirks::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) // Thomas-Imel fit the data. It's also because we are not interested in // extremely low energies, so Thomas-Imel makes a negligible correction. G4double nonIonizingEnergy = totalEnergyDeposit; - nonIonizingEnergy *= NumPhotons/NumQuanta; + nonIonizingEnergy *= NumPhotons/NumQuanta; aParticleChange.ProposeNonIonizingEnergyDeposit(nonIonizingEnergy); return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); } @@ -208,9 +208,9 @@ G4double EDepSim::DokeBirks::GetMeanFreePath(const G4Track&, G4ForceCondition* condition) { *condition = StronglyForced; - // Force the EDepSim::DokeBirks physics process to always happen, even if it's - // not the dominant process. In effect it is a meta-process on top of any - // and all other energy depositions which may occur, just like the + // Force the EDepSim::DokeBirks physics process to always happen, even if + // it's not the dominant process. In effect it is a meta-process on top + // of any and all other energy depositions which may occur, just like the // original G4Scintillation (disregard DBL_MAX, this function makes the // mean free path zero really, not infinite) @@ -224,6 +224,7 @@ G4double EDepSim::DokeBirks::GetMeanLifeTime(const G4Track&, { *condition = Forced; // This function and this condition has the same effect as the above + return DBL_MAX; } @@ -239,4 +240,3 @@ G4double EDepSim::DokeBirks::CalculateElectronLET ( G4double E) { else LET = 0; return LET; } - diff --git a/src/EDepSimExtraPhysics.cc b/src/EDepSimExtraPhysics.cc index e4aa7f4..2ea2594 100644 --- a/src/EDepSimExtraPhysics.cc +++ b/src/EDepSimExtraPhysics.cc @@ -1,9 +1,11 @@ #include "EDepSimExtraPhysics.hh" #include "EDepSimException.hh" +#ifndef EDEPSIM_SKIP_NESTVersion098 // Include NEST #include "NESTVersion098/G4S1Light.hh" #include "NESTVersion098/G4ThermalElectron.hh" +#endif #include "EDepSimDokeBirks.hh" #include "EDepSimLog.hh" @@ -20,13 +22,15 @@ #include #include -EDepSim::ExtraPhysics::ExtraPhysics() +EDepSim::ExtraPhysics::ExtraPhysics() : G4VPhysicsConstructor("EDepSimExtra"), fIonizationModel(1) { } EDepSim::ExtraPhysics::~ExtraPhysics() { } void EDepSim::ExtraPhysics::ConstructParticle() { +#ifndef EDEPSIM_SKIP_NESTVersion098 G4ThermalElectron::Definition(); +#endif } void EDepSim::ExtraPhysics::ConstructProcess() { @@ -44,20 +48,21 @@ void EDepSim::ExtraPhysics::ConstructProcess() { double charge = particle->GetPDGCharge(); if (!pman) { - EDepSimError("Particle " - << particleName + EDepSimError("Particle " + << particleName << " without a Process Manager."); EDepSimThrow("Particle without a Process Manager."); } - + // All charged particles should have a step limiter to make sure that // the steps do not get too long. if (std::abs(charge) > 0.1) { pman->AddDiscreteProcess(new G4StepLimiter("Step Limit")); } - + switch (fIonizationModel) { case 0: { +#ifndef EDEPSIM_SKIP_NESTVersion098 // Add nest to any applicable particle. G4S1Light* scintProcess = new G4S1Light(); scintProcess->SetScintillationYieldFactor(1.0); @@ -65,7 +70,10 @@ void EDepSim::ExtraPhysics::ConstructProcess() { pman->AddProcess(scintProcess,ordDefault, ordInActive,ordDefault); } - break; +#else + EDepSimThrow("NEST not enabled during compilation"); +#endif + break; } case 1: default: { // Add EDepSim::DokeBirks to any applicable particle. @@ -77,5 +85,5 @@ void EDepSim::ExtraPhysics::ConstructProcess() { break; } } - } + } } diff --git a/src/EDepSimHitSegment.cc b/src/EDepSimHitSegment.cc index d18d968..763db5a 100644 --- a/src/EDepSimHitSegment.cc +++ b/src/EDepSimHitSegment.cc @@ -33,20 +33,26 @@ G4Allocator edepHitSegmentAllocator; -EDepSim::HitSegment::HitSegment(double maxSagitta, double maxLength) - : fMaxSagitta(maxSagitta), fMaxLength(maxLength), - fPrimaryId(0), fEnergyDeposit(0), fSecondaryDeposit(0), fTrackLength(0), +EDepSim::HitSegment::HitSegment( + double maxSagitta, double maxSeparation, double maxLength) + : fMaxSagitta(maxSagitta), fMaxSeparation(maxSeparation), + fMaxLength(maxLength), + fPrimaryId(0), fEnergyDeposit(0), fSecondaryDeposit(0), fTrackLength(0), fStart(0,0,0,0), fStop(0,0,0,0) { - fPath.reserve(500); + fPath.reserve(50); + if (fMaxSeparation > fMaxSagitta) fMaxSeparation = fMaxSagitta; + if (fMaxSeparation > fMaxLength) fMaxSeparation = fMaxLength; + } EDepSim::HitSegment::HitSegment(const EDepSim::HitSegment& rhs) - : G4VHit(rhs), - fMaxSagitta(rhs.fMaxSagitta), fMaxLength(rhs.fMaxLength), + : G4VHit(rhs), + fMaxSagitta(rhs.fMaxSagitta), fMaxSeparation(rhs.fMaxSeparation), + fMaxLength(rhs.fMaxLength), fContributors(rhs.fContributors), fPrimaryId(rhs.fPrimaryId), fEnergyDeposit(rhs.fEnergyDeposit), fSecondaryDeposit(rhs.fSecondaryDeposit), - fTrackLength(rhs.fTrackLength), + fTrackLength(rhs.fTrackLength), fStart(rhs.fStart), fStop(rhs.fStop) {} @@ -54,7 +60,7 @@ EDepSim::HitSegment::~HitSegment() { } bool EDepSim::HitSegment::SameHit(G4Step* theStep) { // Check that the hit and new step are in the same volume - G4TouchableHandle touchable + G4TouchableHandle touchable = theStep->GetPreStepPoint()->GetTouchableHandle(); if (fHitVolume != touchable) { return false; @@ -86,10 +92,12 @@ bool EDepSim::HitSegment::SameHit(G4Step* theStep) { } } else { + // Short cut when separation is negative (or zero). + if (fMaxSeparation <= 0.0) return false; // This is not the same track that started this hit, but check to see // if it is a delta-ray that should be added to this segment. double separation = FindSeparation(theStep); - if (separation > fMaxSagitta) { + if (separation > fMaxSeparation) { return false; } } @@ -102,7 +110,7 @@ int EDepSim::HitSegment::FindPrimaryId(G4Track *theTrack) { } void EDepSim::HitSegment::AddStep(G4Step* theStep) { - G4TouchableHandle touchable + G4TouchableHandle touchable = theStep->GetPreStepPoint()->GetTouchableHandle(); G4ThreeVector prePos = theStep->GetPreStepPoint()->GetPosition(); G4ThreeVector postPos = theStep->GetPostStepPoint()->GetPosition(); @@ -112,26 +120,27 @@ void EDepSim::HitSegment::AddStep(G4Step* theStep) { double stepLength = (prePos-postPos).mag(); double trackLength = theStep->GetStepLength(); double nonIonizingDeposit = theStep->GetNonIonizingEnergyDeposit(); - + if (trackLength < 0.75*stepLength || trackLength < stepLength - 1*mm) { - EDepSimWarn("Track length shorter than step: " - << trackLength/mm << " mm " - << "<" << stepLength/mm << " mm"); + EDepSimWarn("Track length shorter than step: " + << trackLength/mm << " mm " + << "<" << stepLength/mm << " mm"); EDepSimWarn(" Volume: " - << theStep->GetTrack()->GetVolume()->GetName()); + << theStep->GetTrack()->GetVolume()->GetName()); EDepSimWarn(" Particle: " << particle->GetParticleName() - << " Depositing " << energyDeposit/MeV << " MeV"); - EDepSimWarn(" Total Energy: " - << theStep->GetTrack()->GetTotalEnergy()); + << " Depositing " << energyDeposit/MeV << " MeV"); + EDepSimWarn(" Total Energy: " + << theStep->GetTrack()->GetTotalEnergy()); } trackLength = std::max(trackLength,stepLength); - EDepSimTrace("Add Step with " << energyDeposit - << " in " << theStep->GetTrack()->GetVolume()->GetName()); + EDepSimTrace("Add Step with " << energyDeposit + << " in " << theStep->GetTrack()->GetVolume()->GetName()); if (energyDeposit <= 0.) { - EDepSimWarn("EDepSim::HitSegment:: No energy deposited: " << energyDeposit); + EDepSimWarn("EDepSim::HitSegment:: No energy deposited: " + << energyDeposit); } if (trackLength <= 0.) { @@ -148,7 +157,7 @@ void EDepSim::HitSegment::AddStep(G4Step* theStep) { stepLength = trackLength = std::min(0.5*mm,0.8*origStep); prePos = postPos - stepLength*mm*dir; EDepSimDebug("EDepSim::HitSegment:: " << particle->GetParticleName() - << " Deposited " << energyDeposit/MeV << " MeV"); + << " Deposited " << energyDeposit/MeV << " MeV"); EDepSimDebug(" Original step: " << origStep/mm << " mm"); EDepSimDebug(" New step: " << stepLength/mm << " mm"); } @@ -156,20 +165,21 @@ void EDepSim::HitSegment::AddStep(G4Step* theStep) { if (stepLength>fMaxLength || trackLength>fMaxLength) { G4Track* trk = theStep->GetTrack(); EDepSimWarn("EDepSim::HitSegment:: Long step in " - << trk->GetVolume()->GetName()); - EDepSimWarn(" Step Length is " - << stepLength/mm - << " mm and track length is " - << trackLength/mm << " mm"); + << trk->GetVolume()->GetName()); + EDepSimWarn(" Step Length is " + << stepLength/mm + << " mm and track length is " + << trackLength/mm << " mm"); EDepSimWarn(" PID: " << particle->GetParticleName() - << " E: " << trk->GetTotalEnergy()/MeV << " MeV" - << " (kin: " << trk->GetKineticEnergy()/MeV << " MeV" - << " Deposit: " - << energyDeposit/MeV << "MeV" - << " Status: " << theStep->GetPreStepPoint()->GetStepStatus() - << " -> " << theStep->GetPostStepPoint()->GetStepStatus()); - + << " E: " << trk->GetTotalEnergy()/MeV << " MeV" + << " (kin: " << trk->GetKineticEnergy()/MeV << " MeV" + << " Deposit: " + << energyDeposit/MeV << "MeV" + << " Status: " + << theStep->GetPreStepPoint()->GetStepStatus() + << " -> " << theStep->GetPostStepPoint()->GetStepStatus()); + const G4VProcess* proc = theStep->GetPostStepPoint() ->GetProcessDefinedStep(); if (!proc) { @@ -177,7 +187,8 @@ void EDepSim::HitSegment::AddStep(G4Step* theStep) { } else { EDepSimWarn(" Process: " << proc->GetProcessName() - <<"/"<< proc->GetProcessTypeName(proc->GetProcessType())); + << "/" + << proc->GetProcessTypeName(proc->GetProcessType())); } } @@ -189,7 +200,7 @@ void EDepSim::HitSegment::AddStep(G4Step* theStep) { theStep->GetPreStepPoint()->GetGlobalTime()); fPath.push_back(fStart.vect()); fStop.set(postPos.x(),postPos.y(),postPos.z(), - theStep->GetPostStepPoint()->GetGlobalTime()); + theStep->GetPostStepPoint()->GetGlobalTime()); fPath.push_back(fStop.vect()); fContributors.push_back(theStep->GetTrack()->GetTrackID()); } @@ -212,11 +223,11 @@ void EDepSim::HitSegment::AddStep(G4Step* theStep) { fTrackLength += trackLength; EDepSimDebug("EDepSim::HitSegment:: Deposit " << particle->GetParticleName() - << " adds " << energyDeposit/MeV << " MeV" - << " to " << fContributors.front() - << " w/ " << fEnergyDeposit - << " L " << trackLength - << " " << fTrackLength); + << " adds " << energyDeposit/MeV << " MeV" + << " to " << fContributors.front() + << " w/ " << fEnergyDeposit + << " L " << trackLength + << " " << fTrackLength); } double EDepSim::HitSegment::FindSagitta(G4Step* theStep) { @@ -228,9 +239,9 @@ double EDepSim::HitSegment::FindSagitta(G4Step* theStep) { if ((point-preStep).mag()>0.01*mm) return 10*m; const G4ThreeVector& postStep = theStep->GetPostStepPoint()->GetPosition(); - // The proposed new segment direction; + // The proposed new segment direction; G4ThreeVector newDir = (postStep-point).unit(); - + // Initialize the maximum sagitta found. double maxSagitta = 0.0; @@ -271,13 +282,12 @@ double EDepSim::HitSegment::FindSeparation(G4Step* theStep) { // the new step. double s1 = (frontDelta - cosDelta*dir).mag(); - // Find the distance from the segment center line to the final point of // the new step. G4ThreeVector rearDelta = postStep-front; cosDelta = rearDelta*dir; double s2 = (rearDelta - cosDelta*dir).mag(); - + return std::max(s1,s2); } @@ -295,9 +305,9 @@ void EDepSim::HitSegment::Draw(void) { } void EDepSim::HitSegment::Print(void) { - std::cout << "Hit Energy Deposit: " - << G4BestUnit(fEnergyDeposit,"Energy") - << std::endl; + std::cout << "Hit Energy Deposit: " + << G4BestUnit(fEnergyDeposit,"Energy") + << std::endl; } double EDepSim::HitSegment::GetLength() const { diff --git a/src/EDepSimHitSegment.hh b/src/EDepSimHitSegment.hh index 8916f29..ddb608e 100644 --- a/src/EDepSimHitSegment.hh +++ b/src/EDepSimHitSegment.hh @@ -2,7 +2,7 @@ // $Id: EDepSim::HitSegment.hh,v 1.5 2011/06/29 04:35:53 mcgrew Exp $ #ifndef EDepSim_HitSegment_h -#define EDepSim_HitSegment_h +#define EDepSim_HitSegment_h #include @@ -34,16 +34,17 @@ public: /// The default values are set so that normally, a scintillator element /// will only have a single hit for a through going track (& delta-rays). HitSegment(double maxSagitta = 1*CLHEP::mm, - double maxLength = 5*CLHEP::mm); - + double maxSeparation = 1*CLHEP::mm, + double maxLength = 5*CLHEP::mm); + HitSegment(const EDepSim::HitSegment& rhs); virtual ~HitSegment(); typedef G4THitsCollection HitSegmentCollection; - + inline void* operator new(size_t); inline void operator delete(void*); - + /// Add the effects of a part of a step to this hit. virtual void AddStep(G4Step* theStep); @@ -93,11 +94,11 @@ public: /// Get the total energy deposited in this hit. double GetEnergyDeposit(void) const {return fEnergyDeposit;} - + /// Get the secondary energy deposited in this hit (see the field /// documentation). double GetSecondaryDeposit(void) const {return fSecondaryDeposit;} - + /// Get the total charged track length in this hit. This includes all of /// the contributions from secondary particles that got lumped into this /// hit (e.g. the contributions from delta-rays). @@ -108,12 +109,12 @@ public: /// The position of the stopping point. const G4LorentzVector& GetStop() const {return fStop;} - + #ifdef BOGUS /// The X position of the hit starting point. Note that a hit by /// definition is in a single volume. If the hit is spread over two /// volumes, it's a result of round-off error (and is almost a bug). The - /// GeoNodeId should be defined by the average position of the hit. + /// GeoNodeId should be defined by the average position of the hit. double GetStartX(void) const {return fStartX;} /// The Y position of the hit starting point. Note that a hit by @@ -158,7 +159,7 @@ public: /// GeoNodeId should be defined by the average position of the hit. double GetStopT(void) const {return fStopT;} #endif - + /// Print the hit information. void ls(std::string = "") const; @@ -182,18 +183,24 @@ private: /// The sagitta tolerance for the segment. double fMaxSagitta; + /// The maximum distance between deposits that can be combined + /// into a segment. This sets the tolerance for combining things + /// like delta-rays into an existing segment. When this is zero, + /// then only one track id can be in a particular segment. + double fMaxSeparation; + /// The maximum length between the start and stop points of the segment. double fMaxLength; /// The TrackID for each trajectory that contributed to this hit. This /// could contain the TrackID of the primary particle, but not - /// necessarily. + /// necessarily. std::vector fContributors; /// The track id of the primary particle. int fPrimaryId; - /// The total energy deposit in this hit. + /// The total energy deposit in this hit. double fEnergyDeposit; /// The "secondary" energy deposit in this hit. This is used to help @@ -208,7 +215,7 @@ private: /// and N_e = N_q - N_ph. Thd fSecondaryDeposit value already includes /// the binomial fluctuation, so don't fluctuate N_ph or N_e. double fSecondaryDeposit; - + /// The total charged track length in this hit. This includes the /// contribution from all of the secondary particles (e.g. delta-rays) /// that are included in this hit. @@ -228,7 +235,7 @@ private: /// make sure that the current hit stays inside of it's allowed /// tolerances. std::vector fPath; - + }; extern G4Allocator edepHitSegmentAllocator; diff --git a/src/EDepSimPersistencyManager.cc b/src/EDepSimPersistencyManager.cc index 4e2e357..fb0fff3 100644 --- a/src/EDepSimPersistencyManager.cc +++ b/src/EDepSimPersistencyManager.cc @@ -432,9 +432,9 @@ void EDepSim::PersistencyManager::MarkTrajectories(const G4Event* event) { continue; } - // Make sure that the primary trajectory associated with this hit - // is saved. The primary trajectories are defined by - // EDepSim::TrajectoryMap::FindPrimaryId(). + // Explicitly save the primary. It will probably be marked again + // with the contributors, but that's OK. This catches some corner + // cases where the primary isn't what you would expect. int primaryId = g4Hit->GetPrimaryTrajectoryId(); EDepSim::Trajectory* ndTraj = dynamic_cast( @@ -445,6 +445,21 @@ void EDepSim::PersistencyManager::MarkTrajectories(const G4Event* event) { else { EDepSimError("Primary trajectory not found"); } + + // Make sure that all the contributors associated with this hit + // are saved. + for (int j = 0; j < g4Hit->GetContributorCount(); ++j) { + int contribId = g4Hit->GetContributor(j); + EDepSim::Trajectory* contribTraj + = dynamic_cast( + EDepSim::TrajectoryMap::Get(contribId)); + if (contribTraj) { + contribTraj->MarkTrajectory(false); + } + else { + EDepSimError("Contributor trajectory not found"); + } + } } } } diff --git a/src/EDepSimPhysicsList.cc b/src/EDepSimPhysicsList.cc index bfe46a3..0811db0 100644 --- a/src/EDepSimPhysicsList.cc +++ b/src/EDepSimPhysicsList.cc @@ -21,6 +21,7 @@ #include +#include #include #include diff --git a/src/EDepSimRootGeometryManager.cc b/src/EDepSimRootGeometryManager.cc index 9436ee3..e1448b7 100644 --- a/src/EDepSimRootGeometryManager.cc +++ b/src/EDepSimRootGeometryManager.cc @@ -22,6 +22,8 @@ #include #include #include +#include + #include #include @@ -46,11 +48,11 @@ #include #include #include +#include #include #include - #include #include #include @@ -82,7 +84,7 @@ int EDepSim::RootGeometryManager::GetNodeId(const G4ThreeVector& pos) { namespace { int CountVolumes(G4LogicalVolume* volume) { int count = 1; - for (std::size_t i = 0; iGetNoDaughters(); ++i) { + for (std::size_t i=0; i < (std::size_t)volume->GetNoDaughters(); ++i) { G4VPhysicalVolume* daughter = volume->GetDaughter(i); count += CountVolumes(daughter->GetLogicalVolume()); } @@ -221,8 +223,10 @@ void EDepSim::RootGeometryManager::Validate() { EDepSimLog("Geometry validated"); } -TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, - TGeoMatrix **returnMatrix) { +TGeoShape* EDepSim::RootGeometryManager::CreateShape( + const std::string& theName, + const G4VSolid* theSolid, + TGeoMatrix **returnMatrix) { const G4String geometryType = theSolid->GetEntityType(); TGeoShape* theShape = NULL; if (geometryType == "G4Box") { @@ -252,7 +256,8 @@ TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, double minPhiDeg = sphere->GetStartPhiAngle()/CLHEP::degree; double maxPhiDeg = minPhiDeg + sphere->GetDeltaPhiAngle()/CLHEP::degree; double minThetaDeg = sphere->GetStartThetaAngle()/CLHEP::degree; - double maxThetaDeg = minThetaDeg + sphere->GetDeltaThetaAngle()/CLHEP::degree; + double maxThetaDeg = minThetaDeg + + sphere->GetDeltaThetaAngle()/CLHEP::degree; theShape = new TGeoSphere(sphere->GetInnerRadius()/CLHEP::mm, sphere->GetOuterRadius()/CLHEP::mm, minThetaDeg, maxThetaDeg, @@ -306,7 +311,9 @@ TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, #else const G4PolyconeHistorical* param = polycone->GetOriginalParameters(); int numZ = param->Num_z_planes; - TGeoPcon* pcon = new TGeoPcon(phi/CLHEP::degree, dPhi/CLHEP::degree, numZ); + TGeoPcon* pcon = new TGeoPcon(phi/CLHEP::degree, + dPhi/CLHEP::degree, + numZ); // This depends on the older interface. It's not marked as // deprecated, but the documentation discourages it's use. for (int i = 0; i< numZ; ++i) { @@ -353,9 +360,9 @@ TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, const G4VSolid* solidB = sub->GetConstituentSolid(1); // solidA - solidB TGeoMatrix* matrixA = NULL; - TGeoShape* shapeA = CreateShape(solidA, &matrixA); + TGeoShape* shapeA = CreateShape(theName, solidA, &matrixA); TGeoMatrix* matrixB = NULL; - TGeoShape* shapeB = CreateShape(solidB, &matrixB); + TGeoShape* shapeB = CreateShape(theName, solidB, &matrixB); TGeoSubtraction* subtractNode = new TGeoSubtraction(shapeA,shapeB, matrixA, matrixB); theShape = new TGeoCompositeShape("name",subtractNode); @@ -366,7 +373,7 @@ TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, const G4VSolid* movedSolid = disp->GetConstituentMovedSolid(); G4RotationMatrix rotation = disp->GetObjectRotation(); G4ThreeVector displacement = disp->GetObjectTranslation(); - theShape = CreateShape(movedSolid); + theShape = CreateShape(theName, movedSolid); if (returnMatrix) { TGeoRotation* rotate = new TGeoRotation("rot", @@ -389,9 +396,9 @@ TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, const G4VSolid* solidB = sub->GetConstituentSolid(1); // solidA - solidB TGeoMatrix* matrixA = NULL; - TGeoShape* shapeA = CreateShape(solidA, &matrixA); + TGeoShape* shapeA = CreateShape(theName, solidA, &matrixA); TGeoMatrix* matrixB = NULL; - TGeoShape* shapeB = CreateShape(solidB, &matrixB); + TGeoShape* shapeB = CreateShape(theName, solidB, &matrixB); TGeoUnion* unionNode = new TGeoUnion(shapeA, shapeB, matrixA, matrixB); theShape = new TGeoCompositeShape("name",unionNode); @@ -403,14 +410,14 @@ TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, const G4VSolid* solidB = sub->GetConstituentSolid(1); // solidA - solidB TGeoMatrix* matrixA = NULL; - TGeoShape* shapeA = CreateShape(solidA, &matrixA); + TGeoShape* shapeA = CreateShape(theName, solidA, &matrixA); TGeoMatrix* matrixB = NULL; - TGeoShape* shapeB = CreateShape(solidB, &matrixB); - TGeoIntersection* intersectionNode = new TGeoIntersection(shapeA, shapeB, - matrixA, matrixB); + TGeoShape* shapeB = CreateShape(theName, solidB, &matrixB); + TGeoIntersection* intersectionNode + = new TGeoIntersection(shapeA, shapeB, + matrixA, matrixB); theShape = new TGeoCompositeShape("name",intersectionNode); } - else if (geometryType == "G4ExtrudedSolid"){ //This following only works when using the 'standard' //G4ExtrudedSolid Constructor. @@ -460,8 +467,16 @@ TGeoShape* EDepSim::RootGeometryManager::CreateShape(const G4VSolid* theSolid, //now assign 'theShape' to this complete extruded object. theShape = xtru; } + else if (geometryType == "G4EllipticalTube") { + const G4EllipticalTube* ellipticalTube + = dynamic_cast(theSolid); + theShape = new TGeoEltu(ellipticalTube->GetDx()/CLHEP::mm, + ellipticalTube->GetDy()/CLHEP::mm, + ellipticalTube->GetDz()/CLHEP::mm); + } else { - EDepSimThrow(geometryType+" --> shape not implemented"); + EDepSimThrow(theName + " :: " + geometryType + + " --> shape not implemented"); } return theShape; @@ -471,7 +486,7 @@ TGeoVolume* EDepSim::RootGeometryManager::CreateVolume(const G4VSolid* theSolid, std::string theName, TGeoMedium* theMedium) { - TGeoShape* theShape = CreateShape(theSolid); + TGeoShape* theShape = CreateShape(theName,theSolid); TGeoVolume* theVolume = new TGeoVolume(theName.c_str(), theShape, theMedium); @@ -481,7 +496,8 @@ EDepSim::RootGeometryManager::CreateVolume(const G4VSolid* theSolid, // Determine if a volume should copied to the ROOT geometry representation. // If this returns true, then the volume and all of it's children will not be // exported. -bool EDepSim::RootGeometryManager::IgnoreVolume(const G4VPhysicalVolume* theVol) { +bool EDepSim::RootGeometryManager::IgnoreVolume( + const G4VPhysicalVolume* theVol) { std::string theFullName = theVol->GetName(); std::string theShortName = theFullName; theShortName.erase(0,theShortName.rfind("/")+1); @@ -585,7 +601,7 @@ void EDepSim::RootGeometryManager::CreateMaterials( // Recurse through the children. for (std::size_t child = 0; - child < theLog->GetNoDaughters(); + child < (std::size_t) theLog->GetNoDaughters(); ++child) { G4VPhysicalVolume* theChild = theLog->GetDaughter(child); CreateMaterials(theChild); @@ -723,7 +739,7 @@ bool EDepSim::RootGeometryManager::CreateEnvelope( // Add the children to the daughter. double missingMass = 0.0; for (std::size_t child = 0; - child < theLog->GetNoDaughters(); + child < (std::size_t) theLog->GetNoDaughters(); ++child) { G4VPhysicalVolume* theChild = theLog->GetDaughter(child); if (theLog->GetNoDaughters() > 20000) { @@ -1010,7 +1026,7 @@ TGeoMedium* EDepSim::RootGeometryManager::AverageMaterial( stack.pop_back(); // Add all of the current children to the stack. for (std::size_t child = 0; - child < currentLog->GetNoDaughters(); + child < (std::size_t) currentLog->GetNoDaughters(); ++child) { stack.push_back(currentLog->GetDaughter(child)); } diff --git a/src/EDepSimRootGeometryManager.hh b/src/EDepSimRootGeometryManager.hh index 8c6d02b..b54078b 100644 --- a/src/EDepSimRootGeometryManager.hh +++ b/src/EDepSimRootGeometryManager.hh @@ -1,6 +1,6 @@ //////////////////////////////////////////////////////////// // $Id: EDepSim::RootGeometryManager.hh,v 1.13 2011/01/17 02:45:46 mcgrew Exp $ -// +// #ifndef EDepSim_RootGeometryManager_hh_Seen #define EDepSim_RootGeometryManager_hh_Seen @@ -51,10 +51,10 @@ public: /// Set the root geometry from a gdml file. void Update(std::string gdmlFile, bool validate); - + /// Make sure that the current geometry passes a bunch of tests. void Validate(); - + /// Get a volume ID base on the volume position. int GetNodeId(const G4ThreeVector& pos); @@ -83,7 +83,7 @@ protected: /// Get the opacity of the physical volume. This will return an opacity of /// zero (i.e. transparent) if the volume should be invisible. double GetOpacity(const G4VPhysicalVolume* vol); - + private: /// The pointer to the instantiation of this object. static EDepSim::RootGeometryManager* fThis; @@ -124,11 +124,13 @@ private: /// become impractially large. The decision is made based on traversing /// the geant geometry and counting the number of unique volumes. bool fCreateAllVolumes; - + /// Create a new ROOT shape object based on the G4 solid. This might be /// called recursively if a G4 boolean solid is found. - TGeoShape* CreateShape(const G4VSolid* theSolid, TGeoMatrix **mat = NULL); - + TGeoShape* CreateShape(const std::string& theName, + const G4VSolid* theSolid, + TGeoMatrix **mat = NULL); + /// Create a new ROOT volume object. TGeoVolume* CreateVolume(const G4VSolid* theSolid, std::string theName, @@ -137,7 +139,7 @@ private: /// Save the detector envelope. This is called recursively to fill in the /// entire detector. The G4 physical volume, theVol, is used to create a /// new root TGeoVolume which is added to the existing root TGeoVolume, - /// theMother, as a daughter. + /// theMother, as a daughter. /// /// RETURN VALUE: This returns true if one of the immediate daughter /// volumes was ignored. This tells the mother volume that it will need @@ -159,7 +161,7 @@ private: /// If this returns true, the volume and all of it's children will not be /// exported to ROOT. This allows the internal structure of low level /// objects (i.e. the internal "structure" of scintillating bars) to be - /// hidden. + /// hidden. virtual bool IgnoreVolume(const G4VPhysicalVolume* theVol); /// A method to flag that a volume mass should be printed. diff --git a/src/EDepSimSegmentSD.cc b/src/EDepSimSegmentSD.cc index 57d4fc5..9b06052 100644 --- a/src/EDepSimSegmentSD.cc +++ b/src/EDepSimSegmentSD.cc @@ -21,7 +21,9 @@ EDepSim::SegmentSD::SegmentSD(G4String name) :G4VSensitiveDetector(name), fHits(NULL), fHCID(-1), - fMaximumHitSagitta(1*CLHEP::mm), fMaximumHitLength(3*CLHEP::mm), + fMaximumHitSagitta(1*CLHEP::mm), + fMaximumHitSeparation(1*CLHEP::mm), + fMaximumHitLength(3*CLHEP::mm), fLastHit(0) { // In an unbelievably poor interface, the G4VSensitiveDetector class // exposes the protected field "std::vector collectionName" to @@ -73,6 +75,7 @@ G4bool EDepSim::SegmentSD::ProcessHits(G4Step* theStep, // If a hit wasn't found, create one. if (!currentHit) { currentHit = new EDepSim::HitSegment(fMaximumHitSagitta, + fMaximumHitSeparation, fMaximumHitLength); fLastHit = fHits->entries(); fHits->insert(currentHit); diff --git a/src/EDepSimSegmentSD.hh b/src/EDepSimSegmentSD.hh index af90f74..3533898 100644 --- a/src/EDepSimSegmentSD.hh +++ b/src/EDepSimSegmentSD.hh @@ -14,8 +14,7 @@ class G4Step; /// A sensitive detector to create EDepSim::HitSegment based hits. namespace EDepSim {class SegmentSD;} -class EDepSim::SegmentSD : public G4VSensitiveDetector -{ +class EDepSim::SegmentSD : public G4VSensitiveDetector { public: SegmentSD(G4String name); @@ -34,15 +33,25 @@ public: } double GetMaximumHitSagitta(void) {return fMaximumHitSagitta;} - /// Set the maximum length for the EDepSim::HitSegment objects created by this - /// sensitive detector. + /// Set the maximum separation between deposits for the + /// EDepSim::HitSegment objects created by this sensitive + /// detector. + void SetMaximumHitSeparation(double separation) { + EDepSimLog("Set max segment separation to " << separation + << " for " << GetName()); + fMaximumHitSeparation = separation; + } + double GetMaximumHitSeparation(void) {return fMaximumHitSeparation;} + + /// Set the maximum length for the EDepSim::HitSegment objects + /// created by this sensitive detector. void SetMaximumHitLength(double length) { EDepSimLog("Set max segment length to " << length - << " for " << GetName()); + << " for " << GetName()); fMaximumHitLength = length; } double GetMaximumHitLength(void) {return fMaximumHitLength;} - + private: /// The collection of hits that is being filled in the current event. It /// is constructed in Initialize, filled in ProcessHits, and added the the @@ -55,6 +64,9 @@ private: /// The maximum allowed sagitta; double fMaximumHitSagitta; + /// The maximum distance between deposits that can be combined. + double fMaximumHitSeparation; + /// The maximum allowed length; double fMaximumHitLength; diff --git a/src/EDepSimTrajectoryPoint.cc b/src/EDepSimTrajectoryPoint.cc index d035e3a..1f9e7fb 100644 --- a/src/EDepSimTrajectoryPoint.cc +++ b/src/EDepSimTrajectoryPoint.cc @@ -18,6 +18,11 @@ G4Allocator aTrajPointAllocator; +const G4String getTextForEnum(int enumVal) +{ + return EnumStrings[enumVal]; +} + EDepSim::TrajectoryPoint::TrajectoryPoint() : fTime(0.), fMomentum(0.,0.,0.), fStepStatus(fUndefined), @@ -128,7 +133,7 @@ std::vector* EDepSim::TrajectoryPoint::CreateAttValues() const { values->push_back(G4AttValue("Momentum", G4BestUnit(fMomentum,"Momentum"),"")); - values->push_back(G4AttValue("StepStatus",fStepStatus,"")); + values->push_back(G4AttValue("StepStatus",getTextForEnum(fStepStatus),"")); values->push_back(G4AttValue("PhysVolName",fPhysVolName,"")); diff --git a/src/EDepSimTrajectoryPoint.hh b/src/EDepSimTrajectoryPoint.hh index dbfce88..f40406e 100644 --- a/src/EDepSimTrajectoryPoint.hh +++ b/src/EDepSimTrajectoryPoint.hh @@ -92,6 +92,24 @@ public: G4ThreeVector fPrevPosition; }; +static const G4String EnumStrings[]{ + "fWorldBoundary", + // Step reached the world boundary + "fGeomBoundary", + // Step defined by a geometry boundary + "fAtRestDoItProc", + // Step defined by a PreStepDoItVector + "fAlongStepDoItProc", + // Step defined by a AlongStepDoItVector + "fPostStepDoItProc", + // Step defined by a PostStepDoItVector + "fUserDefinedLimit", + // Step defined by the user Step limit in the logical volume + "fExclusivelyForcedProc", + // Step defined by an exclusively forced PostStepDoIt process + "fUndefined" // Step not defined yet +}; + #if defined G4TRACKING_ALLOC_EXPORT extern G4DLLEXPORT G4Allocator aTrajPointAllocator; #else diff --git a/src/EDepSimUserDetectorConstruction.cc b/src/EDepSimUserDetectorConstruction.cc index 7d4e11a..ef64798 100644 --- a/src/EDepSimUserDetectorConstruction.cc +++ b/src/EDepSimUserDetectorConstruction.cc @@ -283,7 +283,7 @@ void EDepSim::UserDetectorConstruction::ConstructSDandField() { // handled later. G4LogicalVolume* logVolume = remainingVolumes.front(); remainingVolumes.pop(); - for (std::size_t i = 0; iGetNoDaughters(); ++i) { + for (std::size_t i=0; i<(std::size_t)logVolume->GetNoDaughters(); ++i) { G4VPhysicalVolume* daughter = logVolume->GetDaughter(i); remainingVolumes.push(daughter->GetLogicalVolume()); } @@ -415,6 +415,37 @@ void EDepSim::UserDetectorConstruction::ConstructSDandField() { } logVolume->SetFieldManager(manager,true); } + + /////////////////////////////////////////////////////////////// + // Handle custom material property values + // + // There are several material fields that cannot be directly set using a + // GDML interface. An example of one is be the Birks' constant. To work + // around this the GDML property element is used to set a + // MaterialPropertyVector that is then interpreted here to set the + // specific material field. + for (G4MaterialTable::iterator mat + = G4Material::GetMaterialTable()->begin(); + mat != G4Material::GetMaterialTable()->end(); ++mat) { + G4MaterialPropertiesTable* anMPT = (*mat)->GetMaterialPropertiesTable(); + if (!anMPT) continue; + // Check for a Birks constant value + if (anMPT->ConstPropertyExists("BIRKSCONSTANT")) { + double oldBC = (*mat)->GetIonisation()->GetBirksConstant(); + if (oldBC > 1E-6) { + EDepSimError("Overriding Birks constant for " + << (*mat)->GetName() + << " from nonzero value of" + << " " << oldBC/(mm/MeV) << " mm/MeV"); + } + + double bc = anMPT->GetConstProperty("BIRKSCONSTANT"); + EDepSimLog((*mat)->GetName() << " Birks constant set to" + << " " << bc/(mm/MeV) << " mm/MeV" + << " from BIRKSCONSTANT property"); + (*mat)->GetIonisation()->SetBirksConstant(bc); + } + } } void EDepSim::UserDetectorConstruction::DefineMaterials() { @@ -676,5 +707,6 @@ G4VPhysicalVolume* EDepSim::UserDetectorConstruction::ConstructDetector() { void EDepSim::UserDetectorConstruction::UpdateGeometry() { // Initialize the run manager. This sets everything in motion. + EDepSimLog("Update the geometry and initialize the G4RunManager"); G4RunManager::GetRunManager()->Initialize(); } diff --git a/src/EDepSimUserPrimaryGeneratorMessenger.cc b/src/EDepSimUserPrimaryGeneratorMessenger.cc index 347b41b..1a6b006 100644 --- a/src/EDepSimUserPrimaryGeneratorMessenger.cc +++ b/src/EDepSimUserPrimaryGeneratorMessenger.cc @@ -6,6 +6,7 @@ #include "kinem/EDepSimCombinationGenerator.hh" #include "kinem/EDepSimVKinematicsFactory.hh" #include "kinem/EDepSimGPSKinematicsFactory.hh" +#include "kinem/EDepSimHEPEVTKinematicsFactory.hh" #include "kinem/EDepSimRooTrackerKinematicsFactory.hh" #include "kinem/EDepSimNuMIRockKinematicsFactory.hh" @@ -65,6 +66,7 @@ EDepSim::UserPrimaryGeneratorMessenger::UserPrimaryGeneratorMessenger( ///////////////////////////////////////// AddKinematicsFactory(new EDepSim::GPSKinematicsFactory(this)); + AddKinematicsFactory(new EDepSim::HEPEVTKinematicsFactory(this)); AddKinematicsFactory(new EDepSim::NuMIRockKinematicsFactory(this)); AddKinematicsFactory(new EDepSim::RooTrackerKinematicsFactory(this)); diff --git a/src/NESTVersion098/G4S1Light.cc b/src/NESTVersion098/G4S1Light.cc index e9ee64e..1ae2706 100644 --- a/src/NESTVersion098/G4S1Light.cc +++ b/src/NESTVersion098/G4S1Light.cc @@ -85,10 +85,10 @@ G4S1Light::G4S1Light(const G4String& processName, { //luxManager = LUXSimManager::GetManager(); SetProcessSubType(fScintillation); - + fTrackSecondariesFirst = false; //particles die first, then scintillation is generated - + if (verboseLevel>0) { G4cout << GetProcessName() << " is created " << G4endl; } @@ -110,10 +110,10 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) // this is the most important function, where all light & charge yields happen! { aParticleChange.Initialize(aTrack); - + if ( !YieldFactor ) //set YF=0 when you want S1Light off in your sim return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); - + if( aTrack.GetParentID() == 0 && aTrack.GetCurrentStepNumber() == 1 ) { fExcitedNucleus = false; //an initialization or reset fVeryHighEnergy = false; //initializes or (later) resets this @@ -126,11 +126,11 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) G4String particleName = pDef->GetParticleName(); const G4Material* aMaterial = aStep.GetPreStepPoint()->GetMaterial(); const G4Material* bMaterial = aStep.GetPostStepPoint()->GetMaterial(); - + if((particleName == "neutron" || particleName == "antineutron") && aStep.GetTotalEnergyDeposit() <= 0) return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); - + // code for determining whether the present/next material is noble // element, or, in other words, for checking if either is a valid NEST // scintillating material, and save Z for later L calculation, or @@ -168,11 +168,11 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) j = 0; //no sites yet } //material properties initialized } //end of atomic number check - + if ( !NobleNow && !NobleLater ) return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); - - // retrieval of the particle's position, time, attributes at both the + + // retrieval of the particle's position, time, attributes at both the // beginning and the end of the current step along its track G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint(); G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint(); @@ -181,7 +181,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) G4double evtStrt = pPreStepPoint->GetGlobalTime(); G4double t0 = pPreStepPoint->GetLocalTime(); G4double t1 = pPostStepPoint->GetLocalTime(); - + // now check if we're entering a scintillating material (inside) or // leaving one (outside), in order to determine (later on in the code, // based on the booleans inside & outside) whether to add/subtract @@ -195,13 +195,13 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) ElementA = ElementB; aMaterialPropertiesTable = bMaterial->GetMaterialPropertiesTable(); } - if ( NobleNow && NobleLater && + if ( NobleNow && NobleLater && aMaterial->GetDensity() != bMaterial->GetDensity() ) InsAndOuts = true; - + // Get the LUXSimMaterials pointer //LUXSimMaterials *luxMaterials = LUXSimMaterials::GetMaterials(); - + // retrieve scintillation-related material properties G4double Density = aMaterial->GetDensity()/(CLHEP::g/CLHEP::cm3); G4double nDensity = Density*AVO; //molar mass factor applied below @@ -336,7 +336,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) R0 = 0.0*CLHEP::um; //all Doke/Birks interactions (except for alphas) G4double Townsend = (ElectricField/nDensity)*1e17; DokeBirks[0] = 0.0000; //all geminate (except at zero, low fields) - DokeBirks[2] = 0.1933*pow(Density,2.6199)+0.29754 - + DokeBirks[2] = 0.1933*pow(Density,2.6199)+0.29754 - (0.045439*pow(Density,2.4689)+0.066034)*log10(ElectricField); if ( ElectricField>6990 ) DokeBirks[2]=0.0; if ( ElectricField<1000 ) DokeBirks[2]=0.2; @@ -355,7 +355,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) tau1 = 3.5*CLHEP::ns; tau3 = 20.*CLHEP::ns; tauR = 40.*CLHEP::ns; } //solid Xe } - + // log present and running tally of energy deposition in this section G4double anExcitationEnergy = ((const G4Ions*)(pDef))-> GetExcitationEnergy(); //grab nuclear energy level @@ -363,13 +363,13 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) aMaterialPropertiesTable->GetConstProperty( "ENERGY_DEPOSIT_TOT" ); G4bool convert = false, annihil = false; //set up special cases for pair production and positron annihilation - if(pPreStepPoint->GetKineticEnergy()>=(2*CLHEP::electron_mass_c2) && - !pPostStepPoint->GetKineticEnergy() && + if(pPreStepPoint->GetKineticEnergy()>=(2*CLHEP::electron_mass_c2) && + !pPostStepPoint->GetKineticEnergy() && !aStep.GetTotalEnergyDeposit() && aParticle->GetPDGcode()==22) { convert = true; TotalEnergyDeposit = CLHEP::electron_mass_c2; } - if(pPreStepPoint->GetKineticEnergy() && - !pPostStepPoint->GetKineticEnergy() && + if(pPreStepPoint->GetKineticEnergy() && + !pPostStepPoint->GetKineticEnergy() && aParticle->GetPDGcode()==-11) { annihil = true; TotalEnergyDeposit += aStep.GetTotalEnergyDeposit(); } @@ -385,7 +385,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) AddConstProperty( "ENERGY_DEPOSIT_TOT", TotalEnergyDeposit ); //save current deposit for determining number of quanta produced now TotalEnergyDeposit = aStep.GetTotalEnergyDeposit(); - + // check what the current "goal" E is for dumping scintillation, // often the initial kinetic energy of the parent particle, and deal // with all other energy-related matters in this block of code @@ -426,7 +426,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) } } if ( InsAndOuts ) { - //G4double dribble = pPostStepPoint->GetKineticEnergy() - + //G4double dribble = pPostStepPoint->GetKineticEnergy() - //pPreStepPoint->GetKineticEnergy(); aMaterialPropertiesTable-> AddConstProperty("ENERGY_DEPOSIT_GOL",(-0.1*CLHEP::keV)+ @@ -461,7 +461,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL")==0 && aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_TOT")==0) return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); - + G4String procName; if ( aTrack.GetCreatorProcess() ) procName = aTrack.GetCreatorProcess()->GetProcessName(); @@ -469,7 +469,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) procName = "NULL"; if ( procName == "eBrem" && outside && !OutElectrons ) fMultipleScattering = true; - + // next 2 codeblocks deal with position-related things if ( fAlpha ) delta = 1000.*CLHEP::km; G4int i, k, counter = 0; G4double pos[3]; @@ -478,7 +478,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) fMultipleScattering = true; x1 = x0; //prevents generation of quanta outside active volume } //no scint. for e-'s that leave - + char xCoord[80]; char yCoord[80]; char zCoord[80]; G4bool exists = false; //for querying whether set-up of new site needed for(i=0;iAddConstProperty( xCoord, x1[0] ); @@ -505,14 +505,14 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) aMaterialPropertiesTable-> //save AddConstProperty( "TOTALNUM_INT_SITES", j ); } - + // this is where nuclear recoil "L" factor is handled: total yield is // reduced for nuclear recoil as per Lindhard theory - + //we assume you have a mono-elemental scintillator only //now, grab A's and Z's of current particle and of material (avg) G4double a1 = ElementA->GetA(); - z2 = pDef->GetAtomicNumber(); + z2 = pDef->GetAtomicNumber(); G4double a2 = (G4double)(pDef->GetAtomicMass()); if ( particleName == "alpha" || (z2 == 2 && a2 == 4) ) fAlpha = true; //used later to get S1 pulse shape correct for alpha @@ -522,7 +522,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) G4double gamma = 3.*pow(epsilon,0.15)+0.7*pow(epsilon,0.6)+epsilon; G4double kappa = 0.133*pow(z1,(2./3.))*pow(a2,(-1./2.))*(2./3.); //check if we are dealing with nuclear recoil (Z same as material) - if ( (z1 == z2 && pDef->GetParticleType() == "nucleus") || + if ( (z1 == z2 && pDef->GetParticleType() == "nucleus") || particleName == "neutron" || particleName == "antineutron" ) { YieldFactor=(kappa*gamma)/(1+kappa*gamma); //Lindhard factor if ( z1 == 18 && Phase == kStateLiquid ) @@ -534,11 +534,11 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) if ( z1 == 54 ) ThomasImel = 0.19; if ( z1 == 18 ) ThomasImel = 0.25; } //special TIB parameters for nuclear recoil only, in LXe and LAr - ExcitationRatio = 0.69337 + + ExcitationRatio = 0.69337 + 0.3065*exp(-0.008806*pow(ElectricField,0.76313)); } else YieldFactor = 1.000; //default - + // determine ultimate #quanta from current E-deposition (ph+e-) G4double MeanNumberOfQuanta = //total mean number of exc/ions ScintillationYield*TotalEnergyDeposit; @@ -555,15 +555,15 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) //less than zero because Gaussian fluctuated low, update to zero if(TotalEnergyDeposit < 1/ScintillationYield || NumQuanta < 0) NumQuanta = 0; - + // next section binomially assigns quanta to excitons and ions - G4int NumExcitons = + G4int NumExcitons = BinomFluct(NumQuanta,ExcitationRatio/(1+ExcitationRatio)); G4int NumIons = NumQuanta - NumExcitons; - + // this section calculates recombination following the modified Birks' // Law of Doke, deposition by deposition, and may be overridden later - // in code if a low enough energy necessitates switching to the + // in code if a low enough energy necessitates switching to the // Thomas-Imel box model for recombination instead (determined by site) G4double dE, dx=0, LET=0, recombProb; dE = TotalEnergyDeposit/CLHEP::MeV; @@ -582,7 +582,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) if(dx) LET = (dE/dx)*(1/Density); //lin. energy xfer (prop. to dE/dx) if ( LET > 0 && dE > 0 && dx > 0 ) { G4double ratio = CalculateElectronLET(dE*1e3,z1)/LET; - if ( j == 1 && ratio < 0.7 && !ThomasImelTail && + if ( j == 1 && ratio < 0.7 && !ThomasImelTail && particleName == "e-" ) { dx /= ratio; LET *= ratio; }} } @@ -615,7 +615,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) aParticleChange.ProposeNonIonizingEnergyDeposit(nonIonizingEnergy); return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); #endif - + // next section increments the numbers of excitons, ions, photons, and // electrons for the appropriate interaction site; it only appears to // be redundant by saving seemingly no longer needed exciton and ion @@ -636,7 +636,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) GetConstProperty( numEle ); aMaterialPropertiesTable->AddConstProperty( numPho, NumPhotons ); aMaterialPropertiesTable->AddConstProperty( numEle, NumElectrons ); - + // increment and save the total track length, and save interaction // times for later, when generating the scintillation quanta char trackL[80]; char time00[80]; char time01[80]; char energy[80]; @@ -665,29 +665,29 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) if (t1 > deltaTime) aMaterialPropertiesTable->AddConstProperty( time01, t1 ); } - + // begin the process of setting up creation of scint./ionization TotalEnergyDeposit=aMaterialPropertiesTable-> GetConstProperty("ENERGY_DEPOSIT_TOT"); //get the total E deposited InitialKinetEnergy=aMaterialPropertiesTable-> GetConstProperty("ENERGY_DEPOSIT_GOL"); //E that should have been - if(InitialKinetEnergy > HIENLIM && + if(InitialKinetEnergy > HIENLIM && abs(aParticle->GetPDGcode()) != 2112) fVeryHighEnergy=true; G4double safety; //margin of error for TotalE.. - InitialKinetEnergy if (fVeryHighEnergy && !fExcitedNucleus) safety = 0.2*CLHEP::keV; else safety = 2.*CLHEP::eV; - + //force a scintillation dump for NR and for full nuclear de-excitation - if( !anExcitationEnergy && pDef->GetParticleType() == "nucleus" && + if( !anExcitationEnergy && pDef->GetParticleType() == "nucleus" && aTrack.GetTrackStatus() != fAlive && !fAlpha ) InitialKinetEnergy = TotalEnergyDeposit; if ( particleName == "neutron" || particleName == "antineutron" ) InitialKinetEnergy = TotalEnergyDeposit; - + //force a dump of all saved scintillation under the following - //conditions: energy goal reached, and current particle dead, or an + //conditions: energy goal reached, and current particle dead, or an //error has occurred and total has exceeded goal (shouldn't happen) - if( fabs(TotalEnergyDeposit-InitialKinetEnergy)=InitialKinetEnergy ){ dx = 0; dE = 0; //calculate the total number of quanta from all sites and all @@ -711,10 +711,10 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) if (aTrack.GetTrackStatus() == fAlive ) aParticleChange.ProposeTrackStatus(fSuspend); } - + // begin the loop over all sites which generates all the quanta for(i=0;iGetConstProperty( energy ); t0 = aMaterialPropertiesTable->GetConstProperty( time00 ); t1 = aMaterialPropertiesTable->GetConstProperty( time01 ); - + //if site is small enough, override the Doke/Birks' model with - //Thomas-Imel, but not if we're dealing with super-high energy + //Thomas-Imel, but not if we're dealing with super-high energy //particles, and if it's NR force Thomas-Imel (though NR should be //already short enough in track even up to O(100) keV) if ( (delta < R0 && !fVeryHighEnergy) || z2 == z1 || fAlpha ) { @@ -785,18 +785,18 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) aMaterialPropertiesTable-> AddConstProperty( numEle, NumElectrons ); } - + // grab NumPhotons/NumElectrons, which come from Birks if // the Thomas-Imel block of code above was not executed NumPhotons = (G4int)aMaterialPropertiesTable-> GetConstProperty( numPho ); NumElectrons =(G4int)aMaterialPropertiesTable-> GetConstProperty( numEle ); - + // extra Fano factor caused by recomb. fluct. G4double FanoFactor =0; //ionization channel if(Phase == kStateLiquid && YieldFactor == 1) { - FanoFactor = + FanoFactor = 2575.9*pow((ElectricField+15.154),-0.64064)-1.4707; if ( fKr83m ) TotalEnergyDeposit = 4*CLHEP::keV; if ( (dE/CLHEP::keV) <= 100 && ElectricField >= 0 ) { @@ -810,7 +810,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) if ( Phase == kStateGas && Density>0.5 ) FanoFactor = 0.42857-4.7857*Density+7.8571*pow(Density,2.); if( FanoFactor <= 0 || fVeryHighEnergy ) FanoFactor = 0; - NumQuanta = NumPhotons + NumElectrons; + NumQuanta = NumPhotons + NumElectrons; if(z1==54 && FanoFactor) NumElectrons = G4int( floor(G4RandGauss::shoot(NumElectrons, sqrt(FanoFactor*NumElectrons))+0.5)); @@ -822,20 +822,20 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) NumPhotons = BinomFluct(NumPhotons,biExc); NumPhotons = BinomFluct(NumPhotons,QE_EFF); } NumElectrons = G4int(floor(NumElectrons*phe_per_e+0.5)); - + // new stuff to make Kr-83m work properly if(fKr83m || InitialKinetEnergy==9.4*CLHEP::keV) fKr83m += dE/CLHEP::keV; if(fKr83m > 41) fKr83m = 0; if ( SinglePhase ) //for a 1-phase det. don't propagate e-'s NumElectrons = 0; //saves simulation time - + // reset material properties numExc, numIon, numPho, numEle, as // their values have been used or stored elsewhere already aMaterialPropertiesTable->AddConstProperty( numExc, 0 ); aMaterialPropertiesTable->AddConstProperty( numIon, 0 ); aMaterialPropertiesTable->AddConstProperty( numPho, 0 ); aMaterialPropertiesTable->AddConstProperty( numEle, 0 ); - + // start particle creation loop if( InitialKinetEnergy < MAX_ENE && InitialKinetEnergy > MIN_ENE && !fMultipleScattering ) @@ -844,7 +844,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) for(k = 0; k < NumQuanta; k++) { G4double sampledEnergy; G4DynamicParticle* aQuantum; - + // Generate random direction G4double cost = 1. - 2.*G4UniformRand(); G4double sint = std::sqrt((1.-cost)*(1.+cost)); @@ -852,15 +852,15 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) G4double sinp = std::sin(phi); G4double cosp = std::cos(phi); G4double px = sint*cosp; G4double py = sint*sinp; G4double pz = cost; - + // Create momentum direction vector G4ParticleMomentum photonMomentum(px, py, pz); - + // case of photon-specific stuff if (k < NumPhotons) { - // Determine polarization of new photon + // Determine polarization of new photon G4double sx = cost*cosp; - G4double sy = cost*sinp; + G4double sy = cost*sinp; G4double sz = -sint; G4ThreeVector photonPolarization(sx, sy, sz); G4ThreeVector perp = photonMomentum.cross(photonPolarization); @@ -869,22 +869,22 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) cosp = std::cos(phi); photonPolarization = cosp * photonPolarization + sinp * perp; photonPolarization = photonPolarization.unit(); - + // Generate a new photon or electron: sampledEnergy = G4RandGauss::shoot(PhotMean,PhotWidth); - aQuantum = + aQuantum = new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(), photonMomentum); aQuantum->SetPolarization(photonPolarization.x(), photonPolarization.y(), photonPolarization.z()); } - + else { // this else statement is for ionization electrons if(ElectricField) { // point all electrons straight up, for drifting G4ParticleMomentum electronMomentum(0, 0, -FieldSign); - aQuantum = + aQuantum = new G4DynamicParticle(G4ThermalElectron::ThermalElectron(), electronMomentum); if ( Phase == kStateGas ) { @@ -898,7 +898,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) else { // use "photonMomentum" for the electrons in the case of zero // electric field, which is just randomized vector we made - aQuantum = + aQuantum = new G4DynamicParticle(G4ThermalElectron::ThermalElectron(), photonMomentum); sampledEnergy = 1.38e-23*(CLHEP::joule/CLHEP::kelvin)*Temperature; @@ -909,10 +909,10 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) aQuantum->SetKineticEnergy(sampledEnergy); if (verboseLevel>1) //verbosity stuff G4cout << "sampledEnergy = " << sampledEnergy << G4endl; - + // Generate new G4Track object: // emission time distribution - + // first an initial birth time is provided that is typically // <<1 ns after the initial interaction in the simulation, then // singlet, triplet lifetimes, and recombination time, are @@ -988,8 +988,8 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) if ( Phase == kStateLiquid ) aSecondaryTime -= tauTrap*CLHEP::ns*log(G4UniformRand()); } - - // emission position distribution -- + + // emission position distribution -- // Generate the position of a new photon or electron, with NO // stochastic variation because that could lead to particles // being mistakenly generated outside of your active region by @@ -1043,15 +1043,15 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) if(aSecondaryPosition[2] <= PMT && !GlobalFields) aSecondaryPosition[2] = PMT + R_TOL; } //end of electron diffusion code - + // GEANT4 business: stuff you need to make a new track if ( aSecondaryTime < 0 ) aSecondaryTime = 0; //no neg. time - G4Track * aSecondaryTrack = + G4Track * aSecondaryTrack = new G4Track(aQuantum,aSecondaryTime,aSecondaryPosition); if ( k < NumPhotons || radius < R_MAX ) aParticleChange.AddSecondary(aSecondaryTrack); } - + //reset bunch of things when done with an interaction site aMaterialPropertiesTable->AddConstProperty( xCoord, 999*CLHEP::km ); aMaterialPropertiesTable->AddConstProperty( yCoord, 999*CLHEP::km ); @@ -1060,7 +1060,7 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) aMaterialPropertiesTable->AddConstProperty( energy, 0*CLHEP::eV ); aMaterialPropertiesTable->AddConstProperty( time00, DBL_MAX ); aMaterialPropertiesTable->AddConstProperty( time01, -1*CLHEP::ns ); - + if (verboseLevel>0) { //more verbose stuff G4cout << "\n Exiting from G4S1Light::DoIt -- " << "NumberOfSecondaries = " @@ -1078,13 +1078,13 @@ G4S1Light::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) fExcitedNucleus = false; fAlpha = false; } - + //don't do anything when you're not ready to scintillate else { aParticleChange.SetNumberOfSecondaries(0); return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); } - + //the end (exiting) return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); } @@ -1115,7 +1115,7 @@ G4double G4S1Light::GetMeanLifeTime(const G4Track&, return DBL_MAX; } -G4double GetLiquidElectronDriftSpeed(G4double tempinput, G4double efieldinput, +G4double GetLiquidElectronDriftSpeed(G4double tempinput, G4double efieldinput, G4bool Miller, G4int Z) { if(efieldinput<0) efieldinput *= (-1); //Liquid equation one (165K) coefficients @@ -1150,7 +1150,7 @@ G4double GetLiquidElectronDriftSpeed(G4double tempinput, G4double efieldinput, thri=58.3785811395493, thrj=201512.080026704; G4double y1=0,y2=0,f1=0,f2=0,f3=0,edrift=0, - t1=0,t2=0,frac=0,slope=0,intercept=0; + t1=0,t2=0,slope=0,intercept=0; //Equations defined f1=onea/(1+exp(-(efieldinput-oneb)/onec))+oned/ @@ -1189,12 +1189,12 @@ G4double GetLiquidElectronDriftSpeed(G4double tempinput, G4double efieldinput, else if (tempinput == 200.0) edrift = f2; else if (tempinput == 230.0) edrift = f3; else { //Linear interpolation - frac=(tempinput-t1)/(t2-t1); + // frac=(tempinput-t1)/(t2-t1); slope = (y1-y2)/(t1-t2); intercept=y1-slope*t1; edrift=slope*tempinput+intercept; } - + if ( Miller ) { if ( efieldinput <= 40. ) edrift = -0.13274+0.041082*efieldinput-0.0006886*pow(efieldinput,2.)+ @@ -1248,7 +1248,7 @@ G4int BinomFluct ( G4int N0, G4double prob ) { G4int N1 = 0; if ( prob == 0.00 ) return N1; if ( prob == 1.00 ) return N0; - + if ( N0 < 10 ) { for(G4int i = 0; i < N0; i++) { if(G4UniformRand() < prob) N1++; @@ -1266,7 +1266,7 @@ void InitMatPropValues ( G4MaterialPropertiesTable *nobleElementMat ) { char xCoord[80]; char yCoord[80]; char zCoord[80]; char numExc[80]; char numIon[80]; char numPho[80]; char numEle[80]; char trackL[80]; char time00[80]; char time01[80]; char energy[80]; - + // for loop to initialize the interaction site mat'l properties for( G4int i=0; i<10000; i++ ) { sprintf(xCoord,"POS_X_%d",i); sprintf(yCoord,"POS_Y_%d",i); @@ -1287,7 +1287,7 @@ void InitMatPropValues ( G4MaterialPropertiesTable *nobleElementMat ) { nobleElementMat->AddConstProperty( time00, DBL_MAX ); nobleElementMat->AddConstProperty( time01,-1*CLHEP::ns ); } - + // we initialize the total number of interaction sites, a variable for // updating the amount of energy deposited thus far in the medium, and a // variable for storing the amount of energy expected to be deposited diff --git a/src/NESTVersion098/G4S1Light.hh b/src/NESTVersion098/G4S1Light.hh index 8b57a2d..6d6a6dd 100644 --- a/src/NESTVersion098/G4S1Light.hh +++ b/src/NESTVersion098/G4S1Light.hh @@ -32,7 +32,7 @@ #define phe_per_e 1 //S2 gain for quick studies // different field regions, for gamma-X studies -#define WIN 0*CLHEP::mm //top Cu block (also, quartz window) +#define WIN 0 //top Cu block (also, quartz window) #define TOP 0 //top grid wires #define ANE 0 //anode mesh #define SRF 0 //liquid-gas interface @@ -60,7 +60,7 @@ public: // methods, with descriptions G4double GetMeanFreePath(const G4Track& aTrack, G4double , G4ForceCondition* ); - // Returns infinity; i. e. the process does not limit the step, but + // Returns infinity; i. e. the process does not limit the step, but // sets the 'StronglyForced' condition for the DoIt to be invoked at // every step. @@ -69,15 +69,15 @@ public: // methods, with descriptions // Returns infinity; i. e. the process does not limit the time, but // sets the 'StronglyForced' condition for the DoIt to be invoked at // every step. - + // For in-flight particles losing energy (or those stopped) G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep); G4VParticleChange* AtRestDoIt ( const G4Track& aTrack, const G4Step& aStep); - + // These are the methods implementing the scintillation process. - + void SetTrackSecondariesFirst(const G4bool state); // If set, the primary particle tracking is interrupted and any // produced scintillation quanta are tracked next. When all have been @@ -85,14 +85,14 @@ public: // methods, with descriptions G4bool GetTrackSecondariesFirst() const; // Returns the boolean flag for tracking secondaries first. - + void SetScintillationYieldFactor(const G4double yieldfactor); // Called to set the scintillation quantum yield factor, useful for - // shutting off scintillation entirely, or for producing a universal + // shutting off scintillation entirely, or for producing a universal // re-scaling to for example represent detector effects. Internally is // used for Lindhard yield factor for NR. Default should be user-set // to be 1 (for ER) in your simulation -- see NEST readme - + G4double GetScintillationYieldFactor() const; // Returns the quantum (photon/electron) yield factor. See above. @@ -120,23 +120,23 @@ private: //////////////////// // Inline methods //////////////////// -inline +inline G4bool G4S1Light::IsApplicable(const G4ParticleDefinition& aParticleType) { if (aParticleType.GetParticleName() == "opticalphoton") return false; if (aParticleType.IsShortLived()) return false; if (aParticleType.GetParticleName() == "thermalelectron") return false; //if(abs(aParticleType.GetPDGEncoding())==2112 || //neutron (no E-dep.) - if(abs(aParticleType.GetPDGEncoding())==12 || //neutrinos (ditto) + if(abs(aParticleType.GetPDGEncoding())==12 || //neutrinos (ditto) abs(aParticleType.GetPDGEncoding())==14 || abs(aParticleType.GetPDGEncoding())==16) return false; - + return true; } -inline -void G4S1Light::SetTrackSecondariesFirst(const G4bool state) -{ +inline +void G4S1Light::SetTrackSecondariesFirst(const G4bool state) +{ fTrackSecondariesFirst = state; } diff --git a/src/NESTVersion098/G4S2Light.cc b/src/NESTVersion098/G4S2Light.cc index 608df59..cc563ed 100644 --- a/src/NESTVersion098/G4S2Light.cc +++ b/src/NESTVersion098/G4S2Light.cc @@ -296,7 +296,6 @@ G4double G4S2Light::GetMeanFreePath(const G4Track& aTrack, } else { *condition = NotForced; - const G4Material* aMaterial = aTrack.GetMaterial(); if ( aMaterial->GetDensity() == GRID_DENSITY ) return DBL_MAX; //pass electrons through grid wires return 0*CLHEP::nm; //kill elsewhere diff --git a/src/captain/CaptCryostatBuilder.cc b/src/captain/CaptCryostatBuilder.cc index b11e43c..a742bde 100644 --- a/src/captain/CaptCryostatBuilder.cc +++ b/src/captain/CaptCryostatBuilder.cc @@ -29,7 +29,7 @@ class CaptCryostatMessenger G4UIcmdWithAString* fVesselTypeCMD; public: - CaptCryostatMessenger(CaptCryostatBuilder* c) + CaptCryostatMessenger(CaptCryostatBuilder* c) : EDepSim::BuilderMessenger(c,"Control the driftRegion geometry."), fBuilder(c) { @@ -82,7 +82,7 @@ void CaptCryostatBuilder::Init(void) { SetMessenger(new CaptCryostatMessenger(this)); SetVesselType("CAPTAIN"); - + SetSensitiveDetector("cryo","segment"); AddBuilder(new CaptImmersedBuilder("Immersed",this)); @@ -107,21 +107,29 @@ void CaptCryostatBuilder::DefineCAPTAINVessel() { fInnerVessel.clear(); #include "captainInnerVessel.hxx" - + fOuterVessel.clear(); #include "captainOuterVessel.hxx" + double minZ = 1*CLHEP::km; + double maxZ = -1*CLHEP::km; + double maxOuter = 0.0; for (Shape::reverse_iterator p = fOuterVessel.rbegin(); p != fOuterVessel.rend(); ++p) { - fVesselEnvelope.push_back(*p); + minZ = std::min(p->fZ,minZ); + maxZ = std::max(p->fZ,maxZ); + maxOuter = std::max(p->fOuter,maxOuter); } for (Shape::reverse_iterator p = fInnerVessel.rbegin(); p != fInnerVessel.rend(); ++p) { - if (p->fZ >= fVesselEnvelope.back().fZ) continue; - fVesselEnvelope.push_back(*p); + minZ = std::min(p->fZ,minZ); + maxZ = std::max(p->fZ,maxZ); + maxOuter = std::max(p->fOuter,maxOuter); } - std::reverse(fVesselEnvelope.begin(), fVesselEnvelope.end()); - + maxOuter += 1*CLHEP::mm; + fVesselEnvelope.push_back(Point(minZ,0.0,maxOuter)); + fVesselEnvelope.push_back(Point(maxZ,0.0,maxOuter)); + } void CaptCryostatBuilder::DefineMiniCAPTAINVessel() { @@ -130,21 +138,29 @@ void CaptCryostatBuilder::DefineMiniCAPTAINVessel() { fInnerVessel.clear(); #include "miniCaptainInnerVessel.hxx" - + fOuterVessel.clear(); #include "miniCaptainOuterVessel.hxx" + double minZ = 1*CLHEP::km; + double maxZ = -1*CLHEP::km; + double maxOuter = 0.0; for (Shape::reverse_iterator p = fOuterVessel.rbegin(); p != fOuterVessel.rend(); ++p) { - fVesselEnvelope.push_back(*p); + minZ = std::min(p->fZ,minZ); + maxZ = std::max(p->fZ,maxZ); + maxOuter = std::max(p->fOuter,maxOuter); } for (Shape::reverse_iterator p = fInnerVessel.rbegin(); p != fInnerVessel.rend(); ++p) { - if (p->fZ >= fVesselEnvelope.back().fZ) continue; - fVesselEnvelope.push_back(*p); + minZ = std::min(p->fZ,minZ); + maxZ = std::max(p->fZ,maxZ); + maxOuter = std::max(p->fOuter,maxOuter); } - std::reverse(fVesselEnvelope.begin(), fVesselEnvelope.end()); - + maxOuter += 1*CLHEP::mm; + fVesselEnvelope.push_back(Point(minZ,0.0,maxOuter)); + fVesselEnvelope.push_back(Point(maxZ,0.0,maxOuter)); + } G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { @@ -171,20 +187,24 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { p != fVesselEnvelope.rend(); ++p) { conePlanes.push_back(- p->fZ); - coneMax.push_back(p->fOuter+10*CLHEP::cm); coneMin.push_back(0.0); + coneMax.push_back(p->fOuter+10*CLHEP::cm); } - - G4LogicalVolume* logVolume + + G4LogicalVolume* logVolume = new G4LogicalVolume( new G4Polycone(GetName(), - 0*CLHEP::degree, 360*CLHEP::degree, conePlanes.size(), + 0*CLHEP::degree, 360*CLHEP::degree, + conePlanes.size(), conePlanes.data(), - coneMin.data(), coneMax.data()), + coneMin.data(), + coneMax.data()), FindMaterial("Air"), GetName()); - logVolume->SetVisAttributes(G4VisAttributes::Invisible); - + logVolume->SetVisAttributes(G4VisAttributes::GetInvisible()); + +#define BUILD_OUTER_VESSEL +#ifdef BUILD_OUTER_VESSEL //////////////////////////////////////////////////////// // Define the outer vessel. //////////////////////////////////////////////////////// @@ -202,13 +222,15 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { G4LogicalVolume* logOuterVessel = new G4LogicalVolume( new G4Polycone(GetName()+"/OuterVessel", - 0*CLHEP::degree, 360*CLHEP::degree, conePlanes.size(), + 0*CLHEP::degree, 360*CLHEP::degree, + conePlanes.size(), conePlanes.data(), - coneMin.data(), coneMax.data()), + coneMin.data(), + coneMax.data()), FindMaterial("SS_304"), GetName()+"/OuterVessel"); logOuterVessel->SetVisAttributes(GetColor(logOuterVessel)); - + new G4PVPlacement(NULL, // rotation. G4ThreeVector(0,0,0), logOuterVessel, // logical volume @@ -217,8 +239,10 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { false, // (not used) 0, // Copy number (zero) Check()); // Check overlaps. +#endif - +#define BUILD_INNER_VESSEL +#ifdef BUILD_INNER_VESSEL //////////////////////////////////////////////////////// // Define the inner vessel. //////////////////////////////////////////////////////// @@ -236,13 +260,14 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { G4LogicalVolume* logInnerVessel = new G4LogicalVolume( new G4Polycone(GetName()+"/InnerVessel", - 0*CLHEP::degree, 360*CLHEP::degree, conePlanes.size(), + 0*CLHEP::degree, 360*CLHEP::degree, + conePlanes.size(), conePlanes.data(), coneMin.data(), coneMax.data()), FindMaterial("SS_304"), GetName()+"/InnerVessel"); logInnerVessel->SetVisAttributes(GetColor(logInnerVessel)); - + new G4PVPlacement(NULL, // rotation. G4ThreeVector(0,0,0), logInnerVessel, // logical volume @@ -271,17 +296,18 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { coneMin.push_back(0.0); coneMax.push_back(coneMax.back()); } - + fLiquidVolume = new G4LogicalVolume( new G4Polycone(GetName()+"/Liquid", - 0*CLHEP::degree, 360*CLHEP::degree, conePlanes.size(), + 0*CLHEP::degree, 360*CLHEP::degree, + conePlanes.size(), conePlanes.data(), coneMin.data(), coneMax.data()), FindMaterial("Argon_Liquid"), GetName()+"/Liquid"); fLiquidVolume->SetVisAttributes(GetColor(fLiquidVolume)); - + new G4PVPlacement(NULL, // rotation. G4ThreeVector(0,0,0), fLiquidVolume, // logical volume @@ -322,8 +348,9 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { std::cout << "Undefined immersed volume" << std::endl; std::exit(0); } +#endif - +#ifdef BUILD_ULLAGE //////////////////////////////////////////////////////// // Define the ullage volume. //////////////////////////////////////////////////////// @@ -343,17 +370,19 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { coneMin.push_back(0.0); coneMax.push_back(p->fInner); } - + fUllageVolume = new G4LogicalVolume( new G4Polycone(GetName()+"/Ullage", - 0*CLHEP::degree, 360*CLHEP::degree, conePlanes.size(), + 0*CLHEP::degree, 360*CLHEP::degree, + conePlanes.size(), conePlanes.data(), - coneMin.data(), coneMax.data()), + coneMin.data(), + coneMax.data()), FindMaterial("Argon_Gas"), GetName()+"/Ullage"); fUllageVolume->SetVisAttributes(GetColor(fUllageVolume)); - + new G4PVPlacement(NULL, // rotation. G4ThreeVector(0,0,0), fUllageVolume, // logical volume @@ -398,6 +427,7 @@ G4LogicalVolume *CaptCryostatBuilder::GetPiece(void) { std::cout << "Undefined exposed volume" << std::endl; std::exit(0); } +#endif return logVolume; } diff --git a/src/captain/CaptExposedBuilder.cc b/src/captain/CaptExposedBuilder.cc index e5bf61d..347d503 100644 --- a/src/captain/CaptExposedBuilder.cc +++ b/src/captain/CaptExposedBuilder.cc @@ -22,7 +22,7 @@ class CaptExposedMessenger CaptExposedBuilder* fBuilder; public: - CaptExposedMessenger(CaptExposedBuilder* c) + CaptExposedMessenger(CaptExposedBuilder* c) : EDepSim::BuilderMessenger(c,"Control the exposed geometry."), fBuilder(c) { @@ -44,27 +44,29 @@ void CaptExposedBuilder::Init(void) { CaptExposedBuilder::~CaptExposedBuilder() {} double CaptExposedBuilder::GetRadius() { - double radius = 2*CLHEP::cm; + double radius = 120*CLHEP::cm; return radius; } double CaptExposedBuilder::GetHeight() { - return 1*CLHEP::cm; + return 20*CLHEP::cm; } G4LogicalVolume *CaptExposedBuilder::GetPiece(void) { - G4LogicalVolume* logVolume + std::cout << " Radius " << GetRadius() << std::endl; + std::cout << " Height " << GetHeight() << std::endl; + G4LogicalVolume* logVolume = new G4LogicalVolume(new G4Tubs(GetName(), - 0.0, GetRadius(), GetHeight()/2, + 0.0, GetRadius(), GetHeight()/2, 0*CLHEP::degree, 360*CLHEP::degree), FindMaterial("Argon_Gas"), GetName()); logVolume->SetVisAttributes(GetColor(logVolume)); - G4ThreeVector center(0.0,0.0,-GetHeight()/2); + G4ThreeVector center(0.0,0.0,-GetHeight()/2-0.1*CLHEP::mm); fOffset = center; - + /// All the space above the drift region. center += G4ThreeVector(0.0,0.0,0.0); diff --git a/src/captain/CaptWirePlaneBuilder.cc b/src/captain/CaptWirePlaneBuilder.cc index b4b95bc..5e5888a 100644 --- a/src/captain/CaptWirePlaneBuilder.cc +++ b/src/captain/CaptWirePlaneBuilder.cc @@ -159,7 +159,7 @@ G4LogicalVolume *CaptWirePlaneBuilder::GetPiece(void) { GetHeight()/2), FindMaterial("Argon_Liquid"), GetName()+"/Wire"); - logWire->SetVisAttributes(G4VisAttributes::Invisible); + logWire->SetVisAttributes(G4VisAttributes::GetInvisible()); if (GetSensitiveDetector()) { logWire->SetSensitiveDetector(GetSensitiveDetector()); diff --git a/src/captain/CaptWorldBuilder.cc b/src/captain/CaptWorldBuilder.cc index 1fdd351..490b002 100644 --- a/src/captain/CaptWorldBuilder.cc +++ b/src/captain/CaptWorldBuilder.cc @@ -95,7 +95,7 @@ G4LogicalVolume *CaptWorldBuilder::GetPiece(void) { fHeight/2), FindMaterial("Air"), GetName()); - logVolume->SetVisAttributes(G4VisAttributes::Invisible); + logVolume->SetVisAttributes(G4VisAttributes::GetInvisible()); double floorThickness = 10*CLHEP::cm; G4LogicalVolume *logFloor @@ -105,7 +105,7 @@ G4LogicalVolume *CaptWorldBuilder::GetPiece(void) { floorThickness/2), FindMaterial("Cement"), GetName()+"/Floor"); - logFloor->SetVisAttributes(G4VisAttributes::Invisible); + logFloor->SetVisAttributes(G4VisAttributes::GetInvisible()); CaptCryostatBuilder& cryo = Get("Cryostat"); diff --git a/src/kinem/EDepSimGPSKinematicsGenerator.hh b/src/kinem/EDepSimGPSKinematicsGenerator.hh index 344ee1d..c35d8b1 100644 --- a/src/kinem/EDepSimGPSKinematicsGenerator.hh +++ b/src/kinem/EDepSimGPSKinematicsGenerator.hh @@ -9,17 +9,15 @@ class G4Event; class G4VPrimaryGenerator; -/// A base class for primary generator classes used in the EDepSim::. This is -/// used by the EDepSim::UserPrimaryGeneratorAction to generate particles which -/// will be tracked by the G4 simulation. +/// This is an interface to use the General Particle Source (GPS). namespace EDepSim {class GPSKinematicsGenerator;} class EDepSim::GPSKinematicsGenerator : public EDepSim::VKinematicsGenerator { public: - GPSKinematicsGenerator(const G4String& name, + GPSKinematicsGenerator(const G4String& name, G4VPrimaryGenerator* gps); virtual ~GPSKinematicsGenerator(); - /// Add a primary vertex to the event. + /// Add a primary vertex to the event. virtual GeneratorStatus GeneratePrimaryVertex(G4Event* evt, const G4LorentzVector& position); diff --git a/src/kinem/EDepSimHEPEVTKinematicsFactory.cc b/src/kinem/EDepSimHEPEVTKinematicsFactory.cc new file mode 100644 index 0000000..7a0ec3c --- /dev/null +++ b/src/kinem/EDepSimHEPEVTKinematicsFactory.cc @@ -0,0 +1,59 @@ +#include "kinem/EDepSimHEPEVTKinematicsFactory.hh" +#include "kinem/EDepSimHEPEVTKinematicsGenerator.hh" + +EDepSim::HEPEVTKinematicsFactory::HEPEVTKinematicsFactory( + EDepSim::UserPrimaryGeneratorMessenger* parent) + : EDepSim::VKinematicsFactory("hepevt",parent,false), + fInputFile("not-open"), + fFlavorName("pythia"), + fVerbosity(0), + fInputFileCMD(NULL), + fFlavorCMD(NULL), + fVerboseCMD(NULL) { + + fInputFileCMD = new G4UIcmdWithAString(CommandName("input"),this); + fInputFileCMD->SetGuidance("Set the input file."); + fInputFileCMD->SetParameterName("name",false); + + fFlavorCMD = new G4UIcmdWithAString(CommandName("flavor"),this); + fFlavorCMD->SetGuidance("Choose input format."); + fFlavorCMD->SetParameterName("flavorname",false); + fFlavorCMD->SetCandidates("pythia pbomb marley"); + + fVerboseCMD = new G4UIcmdWithAnInteger(CommandName("verbose"),this); + fVerboseCMD->SetGuidance("Set verbosity level (0 is default, 2 is max)."); + fVerboseCMD->SetParameterName("number",false); +} + +EDepSim::HEPEVTKinematicsFactory::~HEPEVTKinematicsFactory() { + if (fInputFileCMD) delete fInputFileCMD; + if (fFlavorCMD) delete fFlavorCMD; + if (fVerboseCMD) delete fVerboseCMD; +} + +EDepSim::VKinematicsGenerator* +EDepSim::HEPEVTKinematicsFactory::GetGenerator() { + EDepSim::VKinematicsGenerator* kine + = new EDepSim::HEPEVTKinematicsGenerator(GetName(), + GetInputFile(), + GetFlavor(), + GetVerbose()); + return kine; +} + +void EDepSim::HEPEVTKinematicsFactory::SetNewValue(G4UIcommand* command, + G4String newValue) { + if (command == fInputFileCMD) { + SetInputFile(newValue); + } + else if (command == fFlavorCMD) { + SetFlavor(newValue); + } + else if (command == fVerboseCMD) { + SetVerbose(fVerboseCMD->GetNewIntValue(newValue)); + } + else{ + EDepSimError("Nothing to set the value."); + EDepSimThrow("EDepSim::HEPKinematicsFactory::SetNewValue(): Error"); + } +} diff --git a/src/kinem/EDepSimHEPEVTKinematicsFactory.hh b/src/kinem/EDepSimHEPEVTKinematicsFactory.hh new file mode 100644 index 0000000..7e350f8 --- /dev/null +++ b/src/kinem/EDepSimHEPEVTKinematicsFactory.hh @@ -0,0 +1,63 @@ +#ifndef EDepSim_HEPEVTKinematicsFactory_hh_seen +#define EDepSim_HEPEVTKinematicsFactory_hh_seen + +#include "kinem/EDepSimVKinematicsFactory.hh" + +class G4VPrimaryGenerator; +namespace EDepSim {class HEPEVTKinematicsGenerator;} + +namespace EDepSim {class HEPEVTKinematicsFactory;} +class EDepSim::HEPEVTKinematicsFactory : public EDepSim::VKinematicsFactory { +public: + HEPEVTKinematicsFactory(EDepSim::UserPrimaryGeneratorMessenger* fParent); + virtual ~HEPEVTKinematicsFactory(); + + /// Return a new generator enclosing the current factory state. The new + /// generator method is pure virtual so it must be implemented by derived + /// classes. + virtual EDepSim::VKinematicsGenerator* GetGenerator(); + + /// Set the input file to read. + virtual void SetInputFile(const G4String& name) {fInputFile=name;} + + /// Get the input file to read. + virtual const G4String& GetInputFile() const {return fInputFile;} + + /// Set the amount of output from G4HEPEvtInterface. + virtual void SetVerbose(int v) {fVerbosity = v;} + + /// Get the expected level of output from G4HEPEvtInterface. + virtual int GetVerbose() {return fVerbosity;} + + /// Set the input HEPEVT format flavor + virtual void SetFlavor(const G4String& flavorname) {fFlavorName=flavorname;} + + /// Get the input HEPEVT format flavor + virtual const G4String& GetFlavor() const {return fFlavorName;} + + /// Handle the macro command inputs. + virtual void SetNewValue(G4UIcommand* command,G4String newValue); + +private: + + /// The text file of HEPEVT records. + G4String fInputFile; + + /// The flavor of HEPEVT input format + G4String fFlavorName; + + /// The verbosity of the event reading. + G4int fVerbosity; + + /// A command to get the file name to read for the HEPEVT records. + G4UIcmdWithAString* fInputFileCMD; + + /// A command to set the HEPEVT input "flavor", since there + /// apparently isn't a standardized format + G4UIcmdWithAString* fFlavorCMD; + + /// A command to set the HEPEVT verbosity. + G4UIcmdWithAnInteger* fVerboseCMD; + +}; +#endif diff --git a/src/kinem/EDepSimHEPEVTKinematicsGenerator.cc b/src/kinem/EDepSimHEPEVTKinematicsGenerator.cc new file mode 100644 index 0000000..7034a76 --- /dev/null +++ b/src/kinem/EDepSimHEPEVTKinematicsGenerator.cc @@ -0,0 +1,360 @@ +#include "kinem/EDepSimHEPEVTKinematicsGenerator.hh" +#include "EDepSimException.hh" +#include "EDepSimVertexInfo.hh" + +#include +#include +#include +#include +#include + +#ifdef USE_G4HEPEvtInterface +#include +#endif + +#include + +EDepSim::HEPEVTKinematicsGenerator::HEPEVTKinematicsGenerator( + const G4String& name, const G4String& fileName, + const G4String& flavor, int verbosity) + : EDepSim::VKinematicsGenerator(name), fGenerator(NULL), + fFileName(fileName), fFlavor(flavor), + fVerbosity(verbosity) {} + +EDepSim::HEPEVTKinematicsGenerator::~HEPEVTKinematicsGenerator() { +#ifdef USE_G4HEPEvtInterface + if (fGenerator) delete fGenerator; +#endif +} + +int EDepSim::HEPEVTKinematicsGenerator::GetTokens( + std::vector& tokens) { + + tokens.clear(); + if (!fInput.is_open() || fInput.eof()) { + EDepSimError("No more events: " << fFileName); + return 0; + } + + std::string line; + while (std::getline(fInput,line)) { + ++fCurrentLine; + // Strip out comments. + line = line.substr(0,line.find_first_of("#")); + // Break into tokens. + std::istringstream input(line); + std::string token; + tokens.clear(); + while (input >> token) tokens.push_back(token); + if (!tokens.empty()) return tokens.size(); + } + + EDepSimLog("No more data in " << fFileName); + return 0; +} + +int EDepSim::HEPEVTKinematicsGenerator::AsInteger(const std::string& token) { + std::istringstream input(token); + int value; + input >> value; + std::string trailer; + if (input >> trailer) { + EDepSimThrow(fFileName << ":" << fCurrentLine + << " -- Invalid integer: |" << token<< "|"); + } + return value; +} + +double EDepSim::HEPEVTKinematicsGenerator::AsReal(const std::string& token) { + std::istringstream input(token); + double value; + input >> value; + std::string trailer; + if (input >> trailer) { + EDepSimThrow(fFileName << ":" << fCurrentLine + << " -- Invalid real: |" << token << "|"); + } + return value; +} + +EDepSim::VKinematicsGenerator::GeneratorStatus +EDepSim::HEPEVTKinematicsGenerator::GeneratePrimaryVertex( + G4Event* evt, const G4LorentzVector& /* position */) { +#ifdef USE_G4HEPEvtInterface + if (!fGenerator) { + fGenerator = new G4HEPEvtInterface(fFileName,fVerbosity); + } + try { + fGenerator->GeneratePrimaryVertex(evt); + } + catch (...) { + EDepSimError("No more events."); + throw EDepSim::NoMoreEvents(); + } +#else + if (!fInput.is_open()) { + fInput.open(fFileName,std::ifstream::in); + fCurrentLine = 0; + if (!fInput.is_open()) { + EDepSimThrow("File not open: " << fFileName); + } + EDepSimNamedLog("HEPEVT", "Open HEPEVT file: " << fFileName); + } + + + if (fVerbosity>0) { + EDepSimNamedLog("HEPEVT", "Using HEPEVT input flavor: " << fFlavor); + } + + std::vector tokens; + int event_id = 0; + int vertex_id = 0; + int lines = 0; + double vtxX = -999; + double vtxY = -999; + double vtxZ = -999; + double vtxT = -999; + + // Store particles in a list before adding to vertex. This makes + // storing mother/daughter information possible. + std::vector HPlist; + + // Parse the header line for an interaction. Depending on the input + // flavor, this line may contain any number of the following: event ID, + // vertex ID, number of particles, and vertex (x,y,z,t). The number + // of particles is the only required field. Vertex positions should + // be given in cm, and the time should be in ns. + if (!GetTokens(tokens)) { + EDepSimLog("No more events in " << fFileName + << " after " << fCurrentLine << " lines"); + throw EDepSim::NoMoreEvents(); + } + switch (tokens.size()) { + // For pythia flavor input, the number of particle in the event + // is the only token. + case 1: + lines = AsInteger(tokens[0]); + break; + // The next two cases use the marley (15-number) input flavor + // For events with only a single vertex, separate each event by + // a line with two numbers: the eventID and number of particles + // in that event + case 2: + event_id = AsInteger(tokens[0]); + lines = AsInteger(tokens[1]); + break; + // For multi-vertex events, separate each vertex with a line + // containing three numbers: eventID, vertexID, and the number + // of particles in that vertex + case 3: + event_id = AsInteger(tokens[0]); + vertex_id = AsInteger(tokens[1]); + lines = AsInteger(tokens[2]); + evt->SetEventID(event_id); + break; + // Set vertex for pythia flavor + case 5: + vtxX = AsReal(tokens[1])*cm; + vtxY = AsReal(tokens[2])*cm; + vtxZ = AsReal(tokens[3])*cm; + vtxT = AsReal(tokens[4])*ns; + lines = AsInteger(tokens[0]); + break; + // Alternative header format for multi-vertex events: + // event ID, vertex ID, number of particles in that vertex, + // and vertex (x,y,z,t) + case 7: + event_id = AsInteger(tokens[0]); + vertex_id = AsInteger(tokens[1]); + lines = AsInteger(tokens[2]); + evt->SetEventID(event_id); + vtxX = AsReal(tokens[3])*cm; + vtxY = AsReal(tokens[4])*cm; + vtxZ = AsReal(tokens[5])*cm; + vtxT = AsReal(tokens[6])*ns; + break; + default: + EDepSimError("Syntax error at " << fFileName << ":" << fCurrentLine); + throw EDepSim::NoMoreEvents(); + } + + if (fVerbosity>0) { + EDepSimNamedLog("HEPEVT", "Reading " << lines << " particles"); + } + if (fVerbosity>1) { + EDepSim::LogManager::IncreaseIndentation(); + EDepSimNamedLog("HEPEVT", + "Vertex " + << vtxX << " mm, " + << vtxY << " mm, " + << vtxZ << " mm, " + << vtxT << " ns" ); + EDepSim::LogManager::DecreaseIndentation(); + } + + std::unique_ptr vertex( + new G4PrimaryVertex(vtxX,vtxY,vtxZ,vtxT)); + + // Add vertex ID + EDepSim::VertexInfo *vertexInfo = new EDepSim::VertexInfo; + vertex->SetUserInformation(vertexInfo); + vertexInfo->SetInteractionNumber(vertex_id); + + // Parse the particle lines for the vertex and check input flavor + // compatibility. This loses some of the context information for + // now and only adds the particles that should be tracked. + for (int line = 0; line < lines; ++line) { + if (!GetTokens(tokens)) { + EDepSimError("Input truncated in " << fFileName + << " after " << fCurrentLine << " lines"); + throw EDepSim::NoMoreEvents(); + } + if (tokens.size() != 8 && fFlavor == "pythia") { + EDepSimError("Invalid input at " << fFileName + << ":" << fCurrentLine + << " -- Wrong number of entries on line"); + throw EDepSim::NoMoreEvents(); + } + if (tokens.size() != 11 && fFlavor == "pbomb") { + EDepSimError("Invalid input at " << fFileName + << ":" << fCurrentLine + << " -- Wrong number of entries on line"); + throw EDepSim::NoMoreEvents(); + } + if (tokens.size() != 15 && fFlavor == "marley") { + EDepSimError("Invalid input at " << fFileName + << ":" << fCurrentLine + << " -- Wrong number of entries on line"); + throw EDepSim::NoMoreEvents(); + } + int status = AsInteger(tokens[0]); + if (status != 1) continue; + int pid = AsInteger(tokens[1]); + //int mother1, mother2 // ignored, here for documentation + int daughter1, daughter2; + double mass; + double momX, momY, momZ; + //double energy; + + // Read in pythia (8-number) input flavor + // This is the default used in Geant4; see + // https://apc.u-paris.fr/~franco/g4doxy4.10/html/class_g4_h_e_p_evt_interface.html + if (fFlavor == "pythia") { + //daughter1 = AsInteger(tokens[2]); // ignored, here for docu. + //daughter2 = AsInteger(tokens[3]); // ignored, here for docu. + momX = AsReal(tokens[4])*GeV; + momY = AsReal(tokens[5])*GeV; + momZ = AsReal(tokens[6])*GeV; + mass = AsReal(tokens[7])*GeV; + } + // Read in ParticleBomb flavor. Useful for multi-vertex events + // See http://deeplearnphysics.org/DLPGenerator/index.html + else if (fFlavor == "pbomb") { + //mother1 = AsInteger(tokens[2]); // ignored, here for docu. + //mother2 = AsInteger(tokens[3]); // ignored, here for docu. + daughter1 = AsInteger(tokens[4]); + daughter2 = AsInteger(tokens[5]); + momX = AsReal(tokens[6])*GeV; + momY = AsReal(tokens[7])*GeV; + momZ = AsReal(tokens[8])*GeV; + //energy = AsReal(tokens[9])*GeV; + mass = AsReal(tokens[10])*GeV; + } + // Read Marley (15-number) input flavor, also used by LArSoft's + // TextFileGen_module.cc; see + // https://internal.dunescience.org/doxygen/TextFileGen__module_8cc_source.html + else if (fFlavor == "marley") { + //mother1 = AsInteger(tokens[2]); // ignored, here for docu. + //mother2 = AsInteger(tokens[3]); // ignored, here for docu. + daughter1 = AsInteger(tokens[4]); + daughter2 = AsInteger(tokens[5]); + momX = AsReal(tokens[6])*GeV; + momY = AsReal(tokens[7])*GeV; + momZ = AsReal(tokens[8])*GeV; + //energy = AsReal(tokens[9])*GeV; // ignored, here for docu. + mass = AsReal(tokens[10])*GeV; + vtxX = AsReal(tokens[11])*cm; + vtxY = AsReal(tokens[12])*cm; + vtxZ = AsReal(tokens[13])*cm; + vtxT = AsReal(tokens[14])*ns; + vertex->SetPosition(vtxX, vtxY, vtxZ); + vertex->SetT0(vtxT); + if (fVerbosity>1) { + EDepSim::LogManager::IncreaseIndentation(); + EDepSimNamedLog("HEPEVT", + "Vertex " + << vtxX << " mm, " + << vtxY << " mm, " + << vtxZ << " mm, " + << vtxT << " ns" ); + EDepSim::LogManager::DecreaseIndentation(); + } + } + else { + EDepSimError("Unrecognized HEPEVT input flavor " << fFlavor); + throw EDepSim::NoMoreEvents(); + } + + std::unique_ptr part( + new G4PrimaryParticle(pid, momX, momY, momZ)); + part->SetMass(mass); + + // Use Geant's G4HEPEvtParticle class to store daughter information, + // similar to L120 here: + // https://apc.u-paris.fr/~franco/g4doxy/html/G4HEPEvtInterface_8cc-source.html#l00078 + std::unique_ptr hep_particle( + new G4HEPEvtParticle(part.release(), status, daughter1, daughter2)); + HPlist.push_back(hep_particle.release()); + //vertex->SetPrimary(hep_particle.release()->GetTheParticle()); + if (fVerbosity>1) { + EDepSim::LogManager::IncreaseIndentation(); + EDepSimNamedLog("HEPEVT", + "Primary " + << pid << " " << momX << " MeV/c, " + << momY << " MeV/c, " + << momZ << " MeV/c" ); + EDepSim::LogManager::DecreaseIndentation(); + } + } + + // Loop stored HEPParticles and set daughter information if applicable. + // Note that in Geant's notation, JDAHEP1 and JDAHEP2 correspond to + // daughter1 and daughter2, respectively + for(size_t i = 0; i < HPlist.size(); ++i) { + if (HPlist[i]->GetJDAHEP1() > 0) { + // Shift index to start from 0 + G4int jda1 = HPlist[i]->GetJDAHEP1()-1; + G4int jda2 = HPlist[i]->GetJDAHEP2()-1; + G4PrimaryParticle* mother = HPlist[i]->GetTheParticle(); + for( G4int j=jda1; j<=jda2; j++ ) { + G4PrimaryParticle* daughter = HPlist[j]->GetTheParticle(); + if(HPlist[j]->GetISTHEP()>0) { + mother->SetDaughter(daughter); + // Set daughter status to -1 for tracking + HPlist[j]->Done(); + } + } + } + } + + // Add initial particles to the vertex + for(size_t ii=0; iiGetISTHEP() > 0 ) { + G4PrimaryParticle* initialParticle = HPlist[ii]->GetTheParticle(); + vertex->SetPrimary(initialParticle); + } + } + + // Clear G4HEPEvtParticles + for(size_t iii=0;iiiAddPrimaryVertex(vertex.release()); +#endif + + return kSuccess; +} diff --git a/src/kinem/EDepSimHEPEVTKinematicsGenerator.hh b/src/kinem/EDepSimHEPEVTKinematicsGenerator.hh new file mode 100644 index 0000000..df44afd --- /dev/null +++ b/src/kinem/EDepSimHEPEVTKinematicsGenerator.hh @@ -0,0 +1,69 @@ +#ifndef EDepSim_HEPEVTKinematicsGenerator_hh_Seen +#define EDepSim_HEPEVTKinematicsGenerator_hh_Seen + +#include +#include + +#include "kinem/EDepSimVKinematicsGenerator.hh" + +#include +#include + +class G4Event; +class G4VPrimaryGenerator; + +/// A class to read a HEPEVT file. This normally uses an internal processor, +/// but can be forced (using #defines) to use the G4HEPEvtInterface. +namespace EDepSim {class HEPEVTKinematicsGenerator;} +class EDepSim::HEPEVTKinematicsGenerator : + public EDepSim::VKinematicsGenerator { +public: + HEPEVTKinematicsGenerator(const G4String& name, + const G4String& fileName, + const G4String& flavor, + int verbosity = 0); + virtual ~HEPEVTKinematicsGenerator(); + + /// Add a primary vertex to the event. + virtual GeneratorStatus + GeneratePrimaryVertex(G4Event* evt, + const G4LorentzVector& position); + +private: + /// The primary generator when G4HEPEvtInterface is used (usually NULL). + G4VPrimaryGenerator* fGenerator; + + /// The name of the input file (used for error messages) + std::string fFileName; + + ///The input flavor (type of HEPEVT format) to be read + std::string fFlavor; + + /// The number of lines read from the input file (used for error messages) + int fCurrentLine; + + /// The verbosity of the output + int fVerbosity; + + /// The input stream to be read. + std::ifstream fInput; + + /// Get the next line of "tokens" from the input file. This returns the + /// number of tokens that were read from the line, and zero when the file + /// is empty. Comments in the file are prefixed by "#", and tokens are + /// separated by white space. The expected format is + /// \code + /// token1 token2 token3 and so on # comments + /// \endcode + /// Lines that are empty, or only have comments are skipped, + int GetTokens(std::vector& tokens); + + /// Parse a token as an integer. This will throw an error if the string + /// is not a valid integer. + int AsInteger(const std::string& token); + + /// Parse a token as a real number This will throw an error if the string + /// is not a valid floating point number. + double AsReal(const std::string& token); +}; +#endif diff --git a/src/kinem/EDepSimRooTrackerKinematicsGenerator.cc b/src/kinem/EDepSimRooTrackerKinematicsGenerator.cc index 5f71e17..8ad9179 100644 --- a/src/kinem/EDepSimRooTrackerKinematicsGenerator.cc +++ b/src/kinem/EDepSimRooTrackerKinematicsGenerator.cc @@ -6,8 +6,6 @@ #include #include #include -#include -#include #include #include @@ -18,6 +16,7 @@ #include #include #include +#include #include #include @@ -61,42 +60,77 @@ EDepSim::RooTrackerKinematicsGenerator::RooTrackerKinematicsGenerator( std::string basename(filename,start_pos); fFilename = basename + ":" + treeName; + // Count fatal errors while connecting to the rooTracker tree. Some of + // the branches are not used, or can be ignored, but *some* absolutely + // must be correct or the output is wrong. + int fatalErrors = 0; + int branchResult; // Should be zero. + fEvtFlags = NULL; - fTree->SetBranchAddress("EvtFlags", &fEvtFlags); + branchResult = fTree->SetBranchAddress("EvtFlags", &fEvtFlags); fEvtCode = NULL; - fTree->SetBranchAddress("EvtCode", &fEvtCode); - fTree->SetBranchAddress("EvtNum", &fEvtNum); - fTree->SetBranchAddress("EvtXSec", &fEvtXSec); - fTree->SetBranchAddress("EvtDXSec", &fEvtDXSec); - fTree->SetBranchAddress("EvtWght", &fEvtWght); - fTree->SetBranchAddress("EvtProb", &fEvtProb); - fTree->SetBranchAddress("EvtVtx", fEvtVtx); - fTree->SetBranchAddress("StdHepN", &fStdHepN); - fTree->SetBranchAddress("StdHepPdg", fStdHepPdg); - fTree->SetBranchAddress("StdHepStatus", fStdHepStatus); - fTree->SetBranchAddress("StdHepX4", fStdHepX4); - fTree->SetBranchAddress("StdHepP4", fStdHepP4); - fTree->SetBranchAddress("StdHepPolz", fStdHepPolz); - fTree->SetBranchAddress("StdHepFd", fStdHepFd); - fTree->SetBranchAddress("StdHepLd", fStdHepLd); - fTree->SetBranchAddress("StdHepFm", fStdHepFm); - fTree->SetBranchAddress("StdHepLm", fStdHepLm); + branchResult = fTree->SetBranchAddress("EvtCode", &fEvtCode); + branchResult = fTree->SetBranchAddress("EvtNum", &fEvtNum); + if (branchResult) { + EDepSimError("Previous error causes a fatally malformed file"); + ++fatalErrors; + } + branchResult = fTree->SetBranchAddress("EvtXSec", &fEvtXSec); + branchResult = fTree->SetBranchAddress("EvtDXSec", &fEvtDXSec); + branchResult = fTree->SetBranchAddress("EvtWght", &fEvtWght); + branchResult = fTree->SetBranchAddress("EvtProb", &fEvtProb); + branchResult = fTree->SetBranchAddress("EvtVtx", fEvtVtx); + if (branchResult) { + EDepSimError("Previous error causes a fatally malformed file"); + ++fatalErrors; + } + branchResult = fTree->SetBranchAddress("StdHepN", &fStdHepN); + if (branchResult) { + EDepSimError("Previous error causes a fatally malformed file"); + ++fatalErrors; + } + branchResult = fTree->SetBranchAddress("StdHepPdg", fStdHepPdg); + if (branchResult) { + EDepSimError("Previous error causes a fatally malformed file"); + ++fatalErrors; + } + branchResult = fTree->SetBranchAddress("StdHepStatus", fStdHepStatus); + if (branchResult) { + EDepSimError("Previous error causes a fatally malformed file"); + ++fatalErrors; + } + branchResult = fTree->SetBranchAddress("StdHepX4", fStdHepX4); + branchResult = fTree->SetBranchAddress("StdHepP4", fStdHepP4); + if (branchResult) { + EDepSimError("Previous error causes a fatally malformed file"); + ++fatalErrors; + } + branchResult = fTree->SetBranchAddress("StdHepPolz", fStdHepPolz); + branchResult = fTree->SetBranchAddress("StdHepFd", fStdHepFd); + branchResult = fTree->SetBranchAddress("StdHepLd", fStdHepLd); + branchResult = fTree->SetBranchAddress("StdHepFm", fStdHepFm); + branchResult = fTree->SetBranchAddress("StdHepLm", fStdHepLm); #define PARENT_PARTICLE_PASS_THROUGH #ifdef PARENT_PARTICLE_PASS_THROUGH - fTree->SetBranchAddress("NuParentPdg", &fNuParentPdg); - fTree->SetBranchAddress("NuParentDecMode",&fNuParentDecMode); - fTree->SetBranchAddress("NuParentDecP4", fNuParentDecP4); - fTree->SetBranchAddress("NuParentDecX4", fNuParentDecX4); - fTree->SetBranchAddress("NuParentProP4", fNuParentProP4); - fTree->SetBranchAddress("NuParentProX4", fNuParentProX4); - fTree->SetBranchAddress("NuParentProNVtx",&fNuParentProNVtx); + branchResult = fTree->SetBranchAddress("NuParentPdg", &fNuParentPdg); + branchResult = fTree->SetBranchAddress("NuParentDecMode",&fNuParentDecMode); + branchResult = fTree->SetBranchAddress("NuParentDecP4", fNuParentDecP4); + branchResult = fTree->SetBranchAddress("NuParentDecX4", fNuParentDecX4); + branchResult = fTree->SetBranchAddress("NuParentProP4", fNuParentProP4); + branchResult = fTree->SetBranchAddress("NuParentProX4", fNuParentProX4); + branchResult = fTree->SetBranchAddress("NuParentProNVtx",&fNuParentProNVtx); #endif + if (fatalErrors) { + EDepSimError("Fatally malformed rooTracker file"); + throw std::runtime_error("Malformed rooTracker file"); + } + // Set the input tree to the current rootracker tree that this class is // using. EDepSim::KinemPassThrough::GetInstance()->AddInputTree(fTree, - filename, - GetName()); + filename, + GetName()); // Generate a vector of entries to be used by this generator. @@ -125,8 +159,10 @@ EDepSim::RooTrackerKinematicsGenerator::RooTrackerKinematicsGenerator( entry = (entry + stride)%entries; } if (order == "random") { - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - std::shuffle(fEntryVector.begin(), fEntryVector.end(),std::default_random_engine(seed)); + for (std::size_t i = fEntryVector.size()-1; i > 0; --i) { + int j = (i+1)*G4UniformRand(); + std::swap(fEntryVector[i],fEntryVector[j]); + } } if (firstEvent > 0) { @@ -175,6 +211,14 @@ EDepSim::RooTrackerKinematicsGenerator::GeneratePrimaryVertex( EDepSimVerbose("Use rooTracker event number " << fEvtNum << " (entry #" << entry << " in tree)"); + std::string eventCode("not-specified"); + if (fEvtCode) { + eventCode = fEvtCode->String(); + } + else { + EDepSimError("Malformed input file -- missing EvtCode"); + } + // Increment the next entry counter. ++fNextEntry; @@ -184,6 +228,7 @@ EDepSim::RooTrackerKinematicsGenerator::GeneratePrimaryVertex( // Check if this is an end-of-sequence marker. if (fStdHepN == 1 && fStdHepStatus[0]<0) { + EDepSimLog("End of event"); return kEndEvent; } @@ -204,7 +249,8 @@ EDepSim::RooTrackerKinematicsGenerator::GeneratePrimaryVertex( theVertex->SetUserInformation(vertexInfo); // Fill the information fields for this vertex. - vertexInfo->SetReaction(std::string(fEvtCode->String().Data())); + vertexInfo->SetReaction(eventCode); + // Set the file name for this event. std::ostringstream fs; fs << fFilename << ":" << entry; @@ -229,7 +275,7 @@ EDepSim::RooTrackerKinematicsGenerator::GeneratePrimaryVertex( // Add an information field to the vertex. EDepSim::VertexInfo *incomingVertexInfo = new EDepSim::VertexInfo; incomingVertexInfo->SetName("initial-state"); - incomingVertexInfo->SetReaction(std::string(fEvtCode->String().Data())); + incomingVertexInfo->SetReaction(eventCode); theIncomingVertex->SetUserInformation(incomingVertexInfo); // Fill the particles to be tracked (status ==1). These particles are diff --git a/src/kinem/EDepSimVConstrainedPositionGenerator.cc b/src/kinem/EDepSimVConstrainedPositionGenerator.cc index 31d6d6a..29a864b 100644 --- a/src/kinem/EDepSimVConstrainedPositionGenerator.cc +++ b/src/kinem/EDepSimVConstrainedPositionGenerator.cc @@ -268,7 +268,7 @@ void EDepSim::VConstrainedPositionGenerator::FindLimits() { } G4LogicalVolume* logVolume = phyVolume->GetLogicalVolume(); - for (std::size_t i=0; i< logVolume->GetNoDaughters(); ++i) { + for (std::size_t i=0; i<(std::size_t)logVolume->GetNoDaughters(); ++i) { volumes.push(QueueElement(logVolume->GetDaughter(i), translation, rotation)); diff --git a/src/kinem/EDepSimVKinematicsFactory.hh b/src/kinem/EDepSimVKinematicsFactory.hh index 7327837..5cf17a3 100644 --- a/src/kinem/EDepSimVKinematicsFactory.hh +++ b/src/kinem/EDepSimVKinematicsFactory.hh @@ -8,9 +8,9 @@ namespace EDepSim {class VKinematicsGenerator;} namespace EDepSim {class VKinematicsFactory;} class EDepSim::VKinematicsFactory : public EDepSim::VPrimaryFactory { public: - VKinematicsFactory(G4String name, - EDepSim::UserPrimaryGeneratorMessenger* fParent, - bool makeDirectory=true); + VKinematicsFactory(G4String name, + EDepSim::UserPrimaryGeneratorMessenger* fParent, + bool makeDirectory=true); virtual ~VKinematicsFactory(); /// Return a new generator enclosing the current factory state. The new diff --git a/src/kinem/EDepSimVPrimaryFactory.hh b/src/kinem/EDepSimVPrimaryFactory.hh index 9da7e30..b278ea5 100644 --- a/src/kinem/EDepSimVPrimaryFactory.hh +++ b/src/kinem/EDepSimVPrimaryFactory.hh @@ -20,9 +20,9 @@ namespace EDepSim {class VPrimaryFactory;} class EDepSim::VPrimaryFactory : public G4UImessenger { public: - VPrimaryFactory(G4String subdir, G4String name, - EDepSim::UserPrimaryGeneratorMessenger* parent, - bool makeDirectory); + VPrimaryFactory(G4String subdir, G4String name, + EDepSim::UserPrimaryGeneratorMessenger* parent, + bool makeDirectory); virtual ~VPrimaryFactory(); /// Return the full path for the factory. This returns the full name of