diff --git a/common/framework/include/OMSim.hh b/common/framework/include/OMSim.hh index 4a1f4e1780..71d75502b4 100644 --- a/common/framework/include/OMSim.hh +++ b/common/framework/include/OMSim.hh @@ -19,7 +19,7 @@ #include "OMSimSteppingAction.hh" #include "OMSimUIinterface.hh" -#include +#include #include #include @@ -53,28 +53,24 @@ public: void startVisualisationIfRequested(); // OMSimDetectorConstruction* getDetectorConstruction(); - G4Navigator *getNavigator() { return mNavigator; }; + G4Navigator *getNavigator() { return mNavigator.get(); }; void extendOptions(po::options_description pNewOptions); po::options_description mGeneralOptions; private: void initialLoggerConfiguration(); + int determineNumberOfThreads(); po::variables_map parseArguments(int pArgumentCount, char *pArgumentVector[]); void setUserArgumentsToArgTable(po::variables_map pVariablesMap); void setGeneralOptions(); - G4RunManager *mRunManager = nullptr; - G4VisExecutive *mVisManager = nullptr; - G4Navigator *mNavigator = nullptr; - G4VUserPhysicsList *mPhysics = nullptr; - G4VUserPrimaryGeneratorAction *mGenAction = nullptr; - G4UserRunAction *mRunAction = nullptr; - G4UserEventAction *mEventAction = nullptr; - G4UserTrackingAction *mTracking = nullptr; - G4UserSteppingAction *mStepping = nullptr; - G4TouchableHistory *mHistory = nullptr; + std::unique_ptr mRunManager; + std::unique_ptr mVisManager; + std::unique_ptr mPhysics; + std::unique_ptr mHistory; + std::unique_ptr mNavigator; - G4double mStartingTime = 0; + std::chrono::high_resolution_clock::time_point mStartingTime; }; #endif // OMSIM_H diff --git a/common/framework/include/OMSimActionInitialization.hh b/common/framework/include/OMSimActionInitialization.hh new file mode 100644 index 0000000000..9dcd5620ba --- /dev/null +++ b/common/framework/include/OMSimActionInitialization.hh @@ -0,0 +1,19 @@ +#ifndef OMSimActionInitialization_h +#define OMSimActionInitialization_h 1 + +#include "G4VUserActionInitialization.hh" + +class OMSimActionInitialization : public G4VUserActionInitialization +{ +public: + OMSimActionInitialization(long pSeed); + virtual ~OMSimActionInitialization(); + + virtual void BuildForMaster() const; + virtual void Build() const; + +private: + long mMasterSeed; +}; + +#endif \ No newline at end of file diff --git a/common/framework/include/OMSimDataFileTypes.hh b/common/framework/include/OMSimDataFileTypes.hh index b7036e066d..b422342041 100644 --- a/common/framework/include/OMSimDataFileTypes.hh +++ b/common/framework/include/OMSimDataFileTypes.hh @@ -41,6 +41,7 @@ class abcDataFile { public: abcDataFile(G4String pFileName); + ~abcDataFile(); G4String mFileName; G4String mObjectName; @@ -66,7 +67,6 @@ public: abcMaterialData(G4String pFileName) : abcDataFile(pFileName){}; G4Material *mMaterial; G4MaterialPropertiesTable *mMPT; - G4NistManager *mMatDatBase; void createMaterial(); void extractAbsorptionLength(); diff --git a/common/framework/include/OMSimHitManager.hh b/common/framework/include/OMSimHitManager.hh index 802d28374f..eacc356e28 100644 --- a/common/framework/include/OMSimHitManager.hh +++ b/common/framework/include/OMSimHitManager.hh @@ -100,7 +100,6 @@ private: std::map moduleHits; }; G4ThreadLocal static ThreadLocalData* mThreadData; - G4bool mDataWasMerged = false; }; #endif diff --git a/common/framework/include/OMSimRunAction.hh b/common/framework/include/OMSimRunAction.hh index 1bbe56a298..72ee59bb90 100755 --- a/common/framework/include/OMSimRunAction.hh +++ b/common/framework/include/OMSimRunAction.hh @@ -2,7 +2,8 @@ #define OMSimRunAction_h 1 #include - +#include "OMSimHitManager.hh" +#include "OMSimLogger.hh" class G4Run; class OMSimRunAction : public G4UserRunAction @@ -12,8 +13,13 @@ public: ~OMSimRunAction(){}; public: - void BeginOfRunAction(const G4Run*){}; - void EndOfRunAction(const G4Run*){}; + void BeginOfRunAction(const G4Run *) {}; + void EndOfRunAction(const G4Run *run) override + { + log_debug("EndOfRunAction called"); + // Merge thread data here + OMSimHitManager::getInstance().mergeThreadData(); + } }; #endif diff --git a/common/framework/include/OMSimSensitiveDetector.hh b/common/framework/include/OMSimSensitiveDetector.hh index a34da5dcf1..8535b9d17e 100644 --- a/common/framework/include/OMSimSensitiveDetector.hh +++ b/common/framework/include/OMSimSensitiveDetector.hh @@ -51,7 +51,7 @@ class OMSimSensitiveDetector : public G4VSensitiveDetector { public: OMSimSensitiveDetector(G4String pName, DetectorType pDetectorType); - virtual ~OMSimSensitiveDetector() {}; + ~OMSimSensitiveDetector() {}; G4bool ProcessHits(G4Step *pStep, G4TouchableHistory *pTouchableHistory) override; void setPMTResponse(OMSimPMTResponse *pResponse); diff --git a/common/framework/include/OMSimUIinterface.hh b/common/framework/include/OMSimUIinterface.hh index 3e36279c41..c3010bc9ee 100644 --- a/common/framework/include/OMSimUIinterface.hh +++ b/common/framework/include/OMSimUIinterface.hh @@ -57,7 +57,7 @@ public: */ void runBeamOn(G4int pNumberOfEvents = -1) { - log_debug("Running beamOn command"); + log_trace("Running beamOn command"); G4int lNumEvents = pNumberOfEvents >= 0 ? pNumberOfEvents : OMSimCommandArgsTable::getInstance().get("numevents"); applyCommand("/run/beamOn ", lNumEvents); } diff --git a/common/framework/src/OMSim.cc b/common/framework/src/OMSim.cc index 99f9982acf..cb294bb068 100644 --- a/common/framework/src/OMSim.cc +++ b/common/framework/src/OMSim.cc @@ -10,12 +10,18 @@ #include "OMSim.hh" #include "OMSimLogger.hh" +#include "OMSimActionInitialization.hh" #include namespace po = boost::program_options; extern std::shared_ptr globalLogger; -OMSim::OMSim() : mStartingTime(clock() / CLOCKS_PER_SEC), mGeneralOptions("General options"), mRunManager(new G4RunManager()), mVisManager(new G4VisExecutive()), mNavigator(new G4Navigator()) +OMSim::OMSim() : +mStartingTime(std::chrono::high_resolution_clock::now()), +mGeneralOptions("General options"), +mRunManager(std::make_unique()), +mVisManager(std::make_unique()), +mNavigator(std::make_unique()) { setGeneralOptions(); initialLoggerConfiguration(); @@ -41,7 +47,8 @@ void OMSim::setGeneralOptions() ("glass", po::value()->default_value(0), "DEPRECATED. Index to select glass type [VITROVEX = 0, Chiba = 1, Kopp = 2, myVitroVex = 3, myChiba = 4, WOMQuartz = 5, fusedSilica = 6]") ("gel", po::value()->default_value(1), "DEPRECATED. Index to select gel type [Wacker = 0, Chiba = 1, IceCube = 2, Wacker_company = 3]") ("reflective_surface", po::value()->default_value(0), "DEPRECATED. Index to select reflective surface type [Refl_V95Gel = 0, Refl_V98Gel = 1, Refl_Aluminium = 2, Refl_Total98 = 3]") - ("pmt_model", po::value()->default_value(0), "DEPRECATED. R15458 (mDOM) = 0, R7081 (DOM) = 1, 4inch (LOM) = 2, R5912_20_100 (D-Egg)= 3"); + ("pmt_model", po::value()->default_value(0), "DEPRECATED. R15458 (mDOM) = 0, R7081 (DOM) = 1, 4inch (LOM) = 2, R5912_20_100 (D-Egg)= 3") + ("threads", po::value()->default_value(1), "number of threads to use."); } void OMSim::initialLoggerConfiguration() @@ -117,6 +124,21 @@ OMSimDetectorConstruction* OMSim::getDetectorConstruction() } */ +int OMSim::determineNumberOfThreads() +{ + int requestedThreads = OMSimCommandArgsTable::getInstance().get("threads"); + int availableThreads = G4Threading::G4GetNumberOfCores(); + + if (requestedThreads <= 0) { + log_info("Auto-detected {} available cores", availableThreads); + return availableThreads; + } else { + int threadsToUse = std::min(requestedThreads, availableThreads); + log_info("Using {} out of {} available cores", threadsToUse, availableThreads); + return threadsToUse; + } +} + /** * @brief Initialize the simulation constructing all Geant instances. */ @@ -131,41 +153,35 @@ void OMSim::initialiseSimulation(OMSimDetectorConstruction* pDetectorConstructio if (lArgs.get("save_args")) lArgs.writeToJson(lFileName); - CLHEP::HepRandom::setTheEngine(new CLHEP::RanluxEngine(lArgs.get("seed"), 3)); + //CLHEP::HepRandom::setTheEngine(new CLHEP::RanluxEngine(lArgs.get("seed"), 3)); + //CLHEP::HepRandom::setTheEngine(new CLHEP::MixMaxRng(lArgs.get("seed"))); + //G4Random::setTheEngine(new CLHEP::RanluxEngine(lArgs.get("seed"), 3)); + mRunManager->SetUserInitialization(pDetectorConstruction); - mPhysics = new OMSimPhysicsList; - mRunManager->SetUserInitialization(mPhysics); + mPhysics = std::make_unique(); + mRunManager->SetUserInitialization(mPhysics.get()); + mPhysics.release(); mVisManager->Initialize(); - mGenAction = new OMSimPrimaryGeneratorAction(); - mRunManager->SetUserAction(mGenAction); - - mRunAction = new OMSimRunAction(); - mRunManager->SetUserAction(mRunAction); - - mEventAction = new OMSimEventAction(); - mRunManager->SetUserAction(mEventAction); - - mTracking = new OMSimTrackingAction(); - mRunManager->SetUserAction(mTracking); + OMSimActionInitialization* actionInitialization = new OMSimActionInitialization(lArgs.get("seed")); + mRunManager->SetUserInitialization(actionInitialization); - mStepping = new OMSimSteppingAction(); - mRunManager->SetUserAction(mStepping); + // Set number of threads + int nThreads = determineNumberOfThreads(); + mRunManager->SetNumberOfThreads(nThreads); mRunManager->Initialize(); OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance(); lUIinterface.setUI(G4UImanager::GetUIpointer()); - mNavigator->SetWorldVolume(pDetectorConstruction->mWorldPhysical); - mNavigator->LocateGlobalPointAndSetup(G4ThreeVector(0., 0., 0.)); + mNavigator.get()->SetWorldVolume(pDetectorConstruction->mWorldPhysical); + mNavigator.get()->LocateGlobalPointAndSetup(G4ThreeVector(0., 0., 0.)); - mHistory = mNavigator->CreateTouchableHistory(); - - lUIinterface.applyCommand("/control/execute ", lArgs.get("visual")); + mHistory = std::unique_ptr(mNavigator->CreateTouchableHistory()); } /** @@ -238,23 +254,45 @@ bool OMSim::handleArguments(int pArgumentCount, char *pArgumentVector[]) OMSim::~OMSim() { log_trace("OMSim destructor"); - if (mRunManager) { - log_trace("Deleting RunManager"); - delete mRunManager; - mRunManager = nullptr; - } - - if (mVisManager) { - log_trace("Deleting VisManager"); - delete mVisManager; - mVisManager = nullptr; - } - - if (mNavigator) { - log_trace("Deleting Navigator"); - delete mNavigator; - mNavigator = nullptr; - } - double lFinishtime = clock() / CLOCKS_PER_SEC; - log_info("Computation time: {} {}", lFinishtime - mStartingTime, " seconds."); + // if (mRunManager) { + // log_trace("Deleting RunManager"); + // delete mRunManager; + // mRunManager = nullptr; + // } + + // if (mVisManager) { + // log_trace("Deleting VisManager"); + // delete mVisManager; + // mVisManager = nullptr; + // } + + // if (mNavigator) { + // log_trace("Deleting Navigator"); + // delete mNavigator; + // mNavigator = nullptr; + // } + + log_trace("OMSim destructor started"); + + log_trace("Resetting mHistory"); + mHistory.reset(); + + log_trace("Resetting mNavigator"); + mNavigator.reset(); + + log_trace("Resetting mPhysics"); + mPhysics.reset(); + + log_trace("Resetting mVisManager"); + mVisManager.reset(); + + log_trace("Resetting mRunManager"); + mRunManager.reset(); + + log_trace("OMSim destructor finished"); + + + std::chrono::high_resolution_clock::time_point lFinishtime = std::chrono::high_resolution_clock::now(); + const std::chrono::duration lDiff = lFinishtime - mStartingTime; + log_info("Computation time: {} {}", lDiff.count(), " seconds."); } \ No newline at end of file diff --git a/common/framework/src/OMSimActionInitialization.cc b/common/framework/src/OMSimActionInitialization.cc new file mode 100644 index 0000000000..117a08e6aa --- /dev/null +++ b/common/framework/src/OMSimActionInitialization.cc @@ -0,0 +1,36 @@ +#include "OMSimActionInitialization.hh" +#include "OMSimPrimaryGeneratorAction.hh" +#include "OMSimRunAction.hh" +#include "OMSimEventAction.hh" +#include "OMSimTrackingAction.hh" +#include "OMSimSteppingAction.hh" +#include + +OMSimActionInitialization::OMSimActionInitialization(long pSeed) + : G4VUserActionInitialization(), mMasterSeed(pSeed) +{ +} + +OMSimActionInitialization::~OMSimActionInitialization() +{ +} + +void OMSimActionInitialization::BuildForMaster() const +{ + SetUserAction(new OMSimRunAction); +} + +void OMSimActionInitialization::Build() const +{ + SetUserAction(new OMSimPrimaryGeneratorAction); + SetUserAction(new OMSimRunAction); + SetUserAction(new OMSimEventAction); + SetUserAction(new OMSimTrackingAction); + SetUserAction(new OMSimSteppingAction); + const long lPrime = 2147483647; + long lSeed = (mMasterSeed + (G4Threading::G4GetThreadId()+1) * lPrime); + lSeed = lSeed % std::numeric_limits::max(); + log_debug("Random engine of thread {} was assigned seed {}", G4Threading::G4GetThreadId(), lSeed); + G4Random::setTheSeed(lSeed); + G4Random::setTheEngine(new CLHEP::RanluxEngine); +} \ No newline at end of file diff --git a/common/framework/src/OMSimDataFileTypes.cc b/common/framework/src/OMSimDataFileTypes.cc index 7ebd85279a..0ada94422b 100644 --- a/common/framework/src/OMSimDataFileTypes.cc +++ b/common/framework/src/OMSimDataFileTypes.cc @@ -18,6 +18,11 @@ abcDataFile::abcDataFile(G4String pFileName) : mFileData(new ParameterTable()), { } +abcDataFile::~abcDataFile() +{ + //delete mFileData; +} + /** * @brief Sorts two vectors (pSortVector & pReferenceVector) based on the order of values in pReferenceVector. * @@ -70,7 +75,7 @@ void abcMaterialData::createMaterial() { mJsonTree = mFileData->appendAndReturnTree(mFileName); mMPT = new G4MaterialPropertiesTable(); - mMatDatBase = G4NistManager::Instance(); + G4NistManager *lMatDatBase = G4NistManager::Instance(); mObjectName = mJsonTree.get("jName"); log_trace("Creating new material from file {} with name {}", mFileName, mObjectName); @@ -90,10 +95,10 @@ void abcMaterialData::createMaterial() { std::string componentName = key.first; double componentFraction = key.second.get_value(); - mMaterial->AddMaterial(mMatDatBase->FindOrBuildMaterial(componentName), componentFraction); + mMaterial->AddMaterial(lMatDatBase->FindOrBuildMaterial(componentName), componentFraction); } - log_debug("New Material defined: {}", mMaterial->GetName()); + log_trace("New Material defined: {}", mMaterial->GetName()); } /** * @brief Extracts absorption length data from a json-file and adds it to the material property table. @@ -190,11 +195,11 @@ void IceCubeIce::extractInformation() { log_trace("Extracting ice properties from file {}", mFileName); mSpiceDepth_pos = OMSimCommandArgsTable::getInstance().get("depth_pos"); - + G4NistManager *lMatDatBase = G4NistManager::Instance(); createMaterial(); // creates IceCubeICE - G4Material *lIceMie = new G4Material("IceCubeICE_SPICE", mFileData->getValueWithUnit(mObjectName, "jDensity"), mMatDatBase->FindOrBuildMaterial("G4_WATER"), kStateSolid); // create IceCubeICE_SPICE - G4Material *lBubleColumnMie = new G4Material("Mat_BubColumn", mFileData->getValueWithUnit(mObjectName, "jDensity"), mMatDatBase->FindOrBuildMaterial("G4_WATER"), kStateSolid); // create IceCubeICE_SPICE + G4Material *lIceMie = new G4Material("IceCubeICE_SPICE", mFileData->getValueWithUnit(mObjectName, "jDensity"), lMatDatBase->FindOrBuildMaterial("G4_WATER"), kStateSolid); // create IceCubeICE_SPICE + G4Material *lBubleColumnMie = new G4Material("Mat_BubColumn", mFileData->getValueWithUnit(mObjectName, "jDensity"), lMatDatBase->FindOrBuildMaterial("G4_WATER"), kStateSolid); // create IceCubeICE_SPICE std::vector lMieScatteringLength; std::vector lMieScatteringLength_BubleColumn; std::vector lWavelength; @@ -345,13 +350,14 @@ void ReflectiveSurface::extractInformation() } mOpticalSurface->SetMaterialPropertiesTable(lMPT); - log_debug("New Optical Surface: {}", mObjectName); + log_trace("New Optical Surface: {}", mObjectName); } /* * %%%%%%%%%%%%%%%% Functions for scintillation properties %%%%%%%%%%%%%%%% */ + /** * @brief Extract sctintillation properties from json-file and adds them to existing material table */ @@ -366,7 +372,6 @@ void ScintillationProperties::extractInformation() if (lArgs.keyExists("temperature")) { - G4cout << lArgs.keyExists("temperature") << G4endl; lTemperature = lArgs.get("temperature"); } @@ -611,7 +616,7 @@ void CustomProperties::extractConstProperties() G4String lKey = item.first; G4double lValue = mFileData->getValueWithUnit(mJsonTree.get("jName"), "ConstProperties." + lKey); mMPT->AddConstProperty(lKey, lValue, true); - log_debug("Added {} constant property to {}.", lKey, mJsonTree.get("jMaterialName")); + log_trace("Added {} constant property to {}.", lKey, mJsonTree.get("jMaterialName")); } } @@ -638,7 +643,7 @@ void CustomProperties::extractProperties() mMPT->AddProperty(lKey, &lX[0], &lY[0], static_cast(lY.size()), true); G4String mssg = "Added " + lKey + " array property to " + mJsonTree.get("jName"); - log_debug(mssg); + log_trace(mssg); } } diff --git a/common/framework/src/OMSimHitManager.cc b/common/framework/src/OMSimHitManager.cc index 724e011c53..326f32bae1 100644 --- a/common/framework/src/OMSimHitManager.cc +++ b/common/framework/src/OMSimHitManager.cc @@ -3,7 +3,7 @@ #include "OMSimEventAction.hh" #include "OMSimLogger.hh" #include "G4EventManager.hh" -#include "G4Event.hh" +#include "G4Event.hh" #include G4Mutex OMSimHitManager::mMutex = G4Mutex(); @@ -16,6 +16,7 @@ G4ThreadLocal OMSimHitManager::ThreadLocalData *OMSimHitManager::mThreadData = n */ OMSimHitManager &OMSimHitManager::getInstance() { + if (!mInstance) { G4AutoLock lock(&mMutex); @@ -24,6 +25,7 @@ OMSimHitManager &OMSimHitManager::getInstance() mInstance = new OMSimHitManager(); } } + return *mInstance; } @@ -88,8 +90,10 @@ void OMSimHitManager::appendHitInfo( OMSimPMTResponse::PMTPulse pResponse, G4int pModuleNumber) { + if (!mThreadData) { + log_debug("Initialized mThreadData for thread {}", G4Threading::G4GetThreadId()); mThreadData = new ThreadLocalData(); } @@ -112,6 +116,7 @@ void OMSimHitManager::appendHitInfo( moduleHits.localPosition.push_back(pLocalPos); moduleHits.generationDetectionDistance.push_back(pDistance); moduleHits.PMTresponse.push_back(pResponse); + log_trace("Saved hit on module {} sensor {}", pModuleNumber, pPMTHitNumber); } /** @@ -155,8 +160,7 @@ void OMSimHitManager::reset() */ std::vector OMSimHitManager::countHits(int pModuleIndex) { - log_debug("Counting number of detected photons in module with index {}", pModuleIndex); - if (!mDataWasMerged){mergeThreadData();}; + log_trace("Counting number of detected photons in module with index {}", pModuleIndex); HitStats lHitsOfModule = mModuleHits[pModuleIndex]; G4int lNumberPMTs = mNumPMTs[pModuleIndex]; @@ -220,8 +224,7 @@ void OMSimHitManager::sortHitStatsByTime(HitStats &lHits) */ std::vector OMSimHitManager::calculateMultiplicity(const G4double pTimeWindow, int pModuleIndex) { - log_debug("Calculating multiplicity in time window {} for module with index", pTimeWindow, pModuleIndex); - if (!mDataWasMerged){mergeThreadData();}; + log_trace("Calculating multiplicity in time window {} for module with index", pTimeWindow, pModuleIndex); HitStats lHitsOfModule = mModuleHits[pModuleIndex]; G4int lNumberPMTs = mNumPMTs[pModuleIndex]; @@ -273,8 +276,12 @@ void OMSimHitManager::mergeThreadData() G4AutoLock lock(&mMutex); if (mThreadData) { + log_debug("Merging data of different threads"); for (const auto &[moduleIndex, hits] : mThreadData->moduleHits) { + log_debug("Thread ID: {} - Module Index: {} - Hit vector sizes: ={}", + G4Threading::G4GetThreadId(), moduleIndex, hits.eventId.size()); + auto &globalHits = mModuleHits[moduleIndex]; globalHits.eventId.insert(globalHits.eventId.end(), hits.eventId.begin(), hits.eventId.end()); globalHits.hitTime.insert(globalHits.hitTime.end(), hits.hitTime.begin(), hits.hitTime.end()); @@ -291,5 +298,4 @@ void OMSimHitManager::mergeThreadData() delete mThreadData; mThreadData = nullptr; } - mDataWasMerged = true; } \ No newline at end of file diff --git a/common/framework/src/OMSimInputData.cc b/common/framework/src/OMSimInputData.cc index 0bafbac854..dfe032e185 100644 --- a/common/framework/src/OMSimInputData.cc +++ b/common/framework/src/OMSimInputData.cc @@ -271,7 +271,7 @@ G4OpticalSurface *InputDataManager::getOpticalSurface(G4String pName) */ void InputDataManager::searchFolders() { - log_debug("Searching folders for data json files..."); + log_trace("Searching folders for data json files..."); std::vector directories = { "../common/data/Materials", "../common/data/PMTs", "../common/data/scintillation", // ... you can add more directories here ... diff --git a/common/framework/src/OMSimPMTResponse.cc b/common/framework/src/OMSimPMTResponse.cc index 24c91efade..a1acc32544 100644 --- a/common/framework/src/OMSimPMTResponse.cc +++ b/common/framework/src/OMSimPMTResponse.cc @@ -306,7 +306,7 @@ OMSimPMTResponse::PMTPulse OMSimPMTResponse::processPhotocathodeHit(G4double pX, mDOMPMTResponse::mDOMPMTResponse() { log_info("Using mDOM PMT response"); - log_debug("Opening mDOM photocathode scans data..."); + log_trace("Opening mDOM photocathode scans data..."); std::string lPath = "../common/data/PMT_scans/"; mRelativeDetectionEfficiencyInterp = new TGraph((lPath + "weightsVsR_vFit_220nm.txt").c_str()); @@ -324,7 +324,7 @@ mDOMPMTResponse::mDOMPMTResponse() mTransitTimeG2Dmap[lKey] = createHistogramFromData(lPath + "TransitTime_" + lWv + ".dat", ("TransitTime_" + lWv).c_str()); } - log_debug("Finished opening photocathode scans data..."); + log_trace("Finished opening photocathode scans data..."); } std::vector mDOMPMTResponse::getScannedWavelengths() @@ -431,7 +431,7 @@ LOMNNVTResponse::~LOMNNVTResponse() */ NoResponse::NoResponse() { - log_debug("OMSimResponse NoResponse initiated"); + log_trace("OMSimResponse NoResponse initiated"); } std::vector NoResponse::getScannedWavelengths() @@ -442,6 +442,7 @@ std::vector NoResponse::getScannedWavelengths() OMSimPMTResponse::PMTPulse NoResponse::processPhotocathodeHit(G4double pX, G4double pY, G4double pWavelength) { + OMSimPMTResponse::PMTPulse lPulse; lPulse.detectionProbability = -1; lPulse.PE = -1; diff --git a/common/framework/src/OMSimSensitiveDetector.cc b/common/framework/src/OMSimSensitiveDetector.cc index d5dd095947..7c820a2977 100644 --- a/common/framework/src/OMSimSensitiveDetector.cc +++ b/common/framework/src/OMSimSensitiveDetector.cc @@ -4,7 +4,6 @@ #include "G4LogicalVolume.hh" #include "G4ParticleDefinition.hh" #include "G4ParticleTypes.hh" -#include "G4SDManager.hh" #include "G4Step.hh" #include "G4TouchableHistory.hh" #include "G4Track.hh" diff --git a/common/geometry_construction/include/OMSimDetectorConstruction.hh b/common/geometry_construction/include/OMSimDetectorConstruction.hh index b897701e10..a2c21047a5 100644 --- a/common/geometry_construction/include/OMSimDetectorConstruction.hh +++ b/common/geometry_construction/include/OMSimDetectorConstruction.hh @@ -24,7 +24,8 @@ public: ~OMSimDetectorConstruction(); G4VPhysicalVolume *Construct(); - void setSensitiveDetector(G4LogicalVolume* logVol, G4VSensitiveDetector* aSD); + void ConstructSDandField() override; + void registerSensitiveDetector(G4LogicalVolume* logVol, G4VSensitiveDetector* aSD); G4VPhysicalVolume *mWorldPhysical; protected: @@ -33,6 +34,12 @@ protected: virtual void constructWorld() = 0; virtual void constructDetector() = 0; InputDataManager *mData; + + struct SDInfo { + G4LogicalVolume* logicalVolume; + G4VSensitiveDetector* sensitiveDetector; + }; + std::vector mSensitiveDetectors; }; #endif diff --git a/common/geometry_construction/include/OMSimPMTConstruction.hh b/common/geometry_construction/include/OMSimPMTConstruction.hh index 66305f9f0c..e0e836638b 100644 --- a/common/geometry_construction/include/OMSimPMTConstruction.hh +++ b/common/geometry_construction/include/OMSimPMTConstruction.hh @@ -70,7 +70,7 @@ private: G4double mMissingTubeLength; G4PVPlacement *mVacuumBackPhysical; - bool mCheckOverlaps = true; + bool mCheckOverlaps; // Variables from json files are saved in the following members G4double mTotalLenght; diff --git a/common/geometry_construction/src/OMSimDEGGHarness.cc b/common/geometry_construction/src/OMSimDEGGHarness.cc index eceb7e6e54..398ef0a8c8 100644 --- a/common/geometry_construction/src/OMSimDEGGHarness.cc +++ b/common/geometry_construction/src/OMSimDEGGHarness.cc @@ -26,8 +26,7 @@ void DEggHarness::construction() } else { - G4cout << "Penetrator not implemented for Geant4 native version, use -c if you need them" << G4endl; - G4cout << "Ropes not implemented for Geant4 native version, use -c if you need them" << G4endl; + log_warning("Penetrator and harness not implemented for Geant4 native version, use cad implementation if you need them"); G4VSolid *EggHarness = buildHarnessSolid(mRmin, mRmax, mSphi, mDphi, mStheta, mDtheta); G4LogicalVolume *lEggHarness = new G4LogicalVolume(EggHarness, mData->getMaterial("NoOptic_Reflector"), ""); appendComponent(EggHarness, lEggHarness, G4ThreeVector(0, 0, 0), G4RotationMatrix(), "dEGG_Harness"); diff --git a/common/geometry_construction/src/OMSimDetectorConstruction.cc b/common/geometry_construction/src/OMSimDetectorConstruction.cc index 89ff84891b..b8cebd6cfe 100644 --- a/common/geometry_construction/src/OMSimDetectorConstruction.cc +++ b/common/geometry_construction/src/OMSimDetectorConstruction.cc @@ -10,7 +10,6 @@ #include "G4SDManager.hh" - OMSimDetectorConstruction::OMSimDetectorConstruction() : mWorldSolid(0), mWorldLogical(0), mWorldPhysical(0) { @@ -18,7 +17,14 @@ OMSimDetectorConstruction::OMSimDetectorConstruction() OMSimDetectorConstruction::~OMSimDetectorConstruction() { - delete mData; + log_trace("OMSimDetectorConstruction destructor called"); + mSensitiveDetectors.clear(); + if (mData){ + delete mData; + mData = nullptr; + } + + log_trace("OMSimDetectorConstruction destructor finished"); } /** @@ -27,7 +33,7 @@ OMSimDetectorConstruction::~OMSimDetectorConstruction() */ G4VPhysicalVolume *OMSimDetectorConstruction::Construct() { - log_debug("Starting detector construction"); + log_trace("Starting detector construction"); mData = new InputDataManager(); mData->searchFolders(); constructWorld(); @@ -36,10 +42,26 @@ G4VPhysicalVolume *OMSimDetectorConstruction::Construct() return mWorldPhysical; } -void OMSimDetectorConstruction::setSensitiveDetector(G4LogicalVolume *pLogVol, G4VSensitiveDetector *pSD) -{ - log_debug("Setting logical volume {} as sensitive detector {}", pLogVol->GetName(), pSD->GetName()); - auto lSDManager = G4SDManager::GetSDMpointer(); - lSDManager->AddNewDetector(pSD); - SetSensitiveDetector(pLogVol, pSD); +void OMSimDetectorConstruction::ConstructSDandField() +{ + log_trace("ConstructSDandField started"); + + // G4SDManager* lSDManager = G4SDManager::GetSDMpointer(); + for (const auto& sdInfo : mSensitiveDetectors) { + G4String sdName = sdInfo.sensitiveDetector->GetName(); + // if (!lSDManager->FindSensitiveDetector(sdName)) { + // lSDManager->AddNewDetector(sdInfo.sensitiveDetector); + // log_trace("Added new sensitive detector: {}", sdName); + // } + SetSensitiveDetector(sdInfo.logicalVolume, sdInfo.sensitiveDetector); + log_trace("Set sensitive detector {} for logical volume {}", + sdName, sdInfo.logicalVolume->GetName()); + } + log_trace("ConstructSDandField finished"); +} + +void OMSimDetectorConstruction::registerSensitiveDetector(G4LogicalVolume *pLogVol, G4VSensitiveDetector *pSD) +{ + log_trace("Registering logical volume {} as sensitive detector {}", pLogVol->GetName(), pSD->GetName()); + mSensitiveDetectors.push_back({pLogVol, pSD}); } \ No newline at end of file diff --git a/common/geometry_construction/src/OMSimLOM18.cc b/common/geometry_construction/src/OMSimLOM18.cc index defe2e07e9..5b3de3de6d 100644 --- a/common/geometry_construction/src/OMSimLOM18.cc +++ b/common/geometry_construction/src/OMSimLOM18.cc @@ -275,7 +275,7 @@ void LOM18::placeCADSupportStructure(G4LogicalVolume* lInnerVolumeLogical) CADfile << "EverythingButPMTsGelpadsVesselPenetrator.obj"; //CADfile << "Internal.obj"; - G4cout << "using the following CAD file for support structure: " << CADfile.str() << G4endl; + log_debug("Using the following CAD file for support structure: {}", CADfile.str()); //load mesh auto mesh = CADMesh::TessellatedMesh::FromOBJ("../common/data/CADmeshes/LOM18/" + CADfile.str() ); diff --git a/common/geometry_construction/src/OMSimMDOM.cc b/common/geometry_construction/src/OMSimMDOM.cc index 06e74dae2e..71f0441e00 100644 --- a/common/geometry_construction/src/OMSimMDOM.cc +++ b/common/geometry_construction/src/OMSimMDOM.cc @@ -35,7 +35,7 @@ mDOM::mDOM(InputDataManager *pData, G4bool pPlaceHarness): OMSimOpticalModule(pD } construction(); - log_debug("Finished constructing mDOM"); + log_trace("Finished constructing mDOM"); } void mDOM::construction() diff --git a/common/geometry_construction/src/OMSimPMTConstruction.cc b/common/geometry_construction/src/OMSimPMTConstruction.cc index dbe7fb92ae..36b27cc05f 100644 --- a/common/geometry_construction/src/OMSimPMTConstruction.cc +++ b/common/geometry_construction/src/OMSimPMTConstruction.cc @@ -67,6 +67,7 @@ void OMSimPMTConstruction::construction() OMSimPMTResponse *OMSimPMTConstruction::getPMTResponseInstance() { + G4String jResponseData; try { @@ -103,9 +104,9 @@ OMSimPMTResponse *OMSimPMTConstruction::getPMTResponseInstance() void OMSimPMTConstruction::configureSensitiveVolume(OMSimDetectorConstruction *pDetConst, G4String pName) { - auto lSensitiveDetector = new OMSimSensitiveDetector(pName, DetectorType::PMT); + OMSimSensitiveDetector* lSensitiveDetector = new OMSimSensitiveDetector(pName, DetectorType::PMT); lSensitiveDetector->setPMTResponse(getPMTResponseInstance()); - pDetConst->setSensitiveDetector(mPhotocathodeLV, lSensitiveDetector); + pDetConst->registerSensitiveDetector(mPhotocathodeLV, lSensitiveDetector); } void OMSimPMTConstruction::constructHAcoating() @@ -394,11 +395,6 @@ void OMSimPMTConstruction::checkPhotocathodeThickness() G4double lEllipsePos_y = mData->getValueWithUnit(mSelectedPMT, lSide + ".jEllipsePos_y"); lSide = "jInnerShape"; - G4cout << "lOutRad " << (mData->getValueWithUnit(mSelectedPMT, lSide + ".jOutRad") - lOutRad) / nm << G4endl; - G4cout << "lEllipseXYaxis " << (mData->getValueWithUnit(mSelectedPMT, lSide + ".jEllipseXYaxis") - lEllipseXYaxis) / nm << G4endl; - G4cout << "lEllispesZaxis " << (mData->getValueWithUnit(mSelectedPMT, lSide + ".jEllipseZaxis") - lEllipseZaxis) / nm << G4endl; - G4cout << "lSpherePos_y " << (mData->getValueWithUnit(mSelectedPMT, lSide + ".jSpherePos_y") - lSpherePos_y) / nm << G4endl; - G4cout << "lEllipsePos_y " << (mData->getValueWithUnit(mSelectedPMT, lSide + ".jEllipsePos_y") - lEllipsePos_y) / nm << G4endl; } /** diff --git a/common/geometry_construction/src/abcDetectorComponent.cc b/common/geometry_construction/src/abcDetectorComponent.cc index de7630c3ac..58f2c99d20 100644 --- a/common/geometry_construction/src/abcDetectorComponent.cc +++ b/common/geometry_construction/src/abcDetectorComponent.cc @@ -122,7 +122,7 @@ void abcDetectorComponent::placeIt(G4ThreeVector pPosition, G4RotationMatrix pRo for (auto const &[key, Component] : mComponents) { G4String mssg = "Placing " + key + " in " + pMother->GetName() + "."; - log_debug(mssg); + log_trace(mssg); lTrans = getNewPosition(pPosition, pRotation, Component.Position, Component.Rotation); mLastPhysicals[key] = new G4PVPlacement(lTrans, Component.VLogical, Component.Name + pNameExtension, pMother, false, 0, mCheckOverlaps); } @@ -141,7 +141,7 @@ void abcDetectorComponent::placeIt(G4Transform3D pTrans, G4LogicalVolume *&pMoth G4Transform3D lTrans; for (auto const &[key, Component] : mComponents) { - log_debug("Placing {} in {}.", key, pMother->GetName()); + log_trace("Placing {} in {}.", key, pMother->GetName()); lTrans = getNewPosition(pTrans.getTranslation(), pTrans.getRotation(), Component.Position, Component.Rotation); mLastPhysicals[key] = new G4PVPlacement(lTrans, Component.VLogical, Component.Name + pNameExtension, pMother, false, 0, mCheckOverlaps); } @@ -186,7 +186,7 @@ G4SubtractionSolid *abcDetectorComponent::substractToVolume(G4VSolid *pInputVolu { lTrans = getNewPosition(pSubstractionPos, pSubstractionRot, Component.Position, Component.Rotation); G4String mssg = "Substracting " + key + " from " + pInputVolume->GetName() + "."; - log_debug(mssg); + log_trace(mssg); if (lCounter == 0) { lSubstractedVolume = new G4SubtractionSolid("SubstractedVolume", pInputVolume, Component.VSolid, lTrans); @@ -217,7 +217,7 @@ G4SubtractionSolid *abcDetectorComponent::substractToVolume(G4VSolid *pInputVolu for (auto const &[key, Component] : mComponents) { G4String mssg = "Substracting " + key + " from " + pInputVolume->GetName() + "."; - log_debug(mssg); + log_trace(mssg); if (lCounter == 0) { lSubstractedVolume = new G4SubtractionSolid("SubstractedVolume", pInputVolume, Component.VSolid, pTrans); diff --git a/doxygen_documentation/0common.h b/doxygen_documentation/0common.h index 3068fc31a4..c4c6d88cae 100644 --- a/doxygen_documentation/0common.h +++ b/doxygen_documentation/0common.h @@ -68,7 +68,7 @@ * @warning Only the mDOM PMT currently supports a detailed PMT response. * * For photon detection in both simple and complex geometries, the photons must be absorbed within the photocathode volume. The PMTs' photocathodes are made sensitive through the OMSimSensitiveDetector class, following Geant4's G4VSensitiveDetector pattern. This configuration is achieved by invoking `OMSimOpticalModule::configureSensitiveVolume` (or `OMSimPMTConstruction::configureSensitiveVolume` when simulating a single PMT). - * It is essential to invoke this method in the detector construction, as it needs the instance of `OMSimDetectorConstruction` to call `G4VUserDetectorConstruction::SetSensitiveDetector` for successful operation in Geant4 (refer to `OMSimDetectorConstruction::setSensitiveDetector`). + * It is essential to invoke this method in the detector construction, as it needs the instance of `OMSimDetectorConstruction` to call `G4VUserDetectorConstruction::SetSensitiveDetector` for successful operation in Geant4 (refer to `OMSimDetectorConstruction::registerSensitiveDetector`). * * Every step of a particle through the photocathode triggers the `OMSimSensitiveDetector::ProcessHits` method. It verifies if the particle is a photon and whether it was absorbed. For a deeper understanding of Geant4's philosophy concerning G4VSensitiveDetector, consult the Geant4 guide for application developers. * @@ -101,7 +101,7 @@ * OMSimSensitiveDetector* lSensitiveDetector = new OMSimSensitiveDetector("myDetector", DetectorType::GeneralPhotonDetector); * OMSimHitManager lHitManager = OMSimHitManager::getInstance(); * lHitManager.setNumberOfPMTs(1, lHitManager.getNextDetectorIndex()); - * setSensitiveDetector(lDetectorLV, lSensitiveDetector); + * registerSensitiveDetector(lDetectorLV, lSensitiveDetector); * ~~~~~~~~~~~~~ * In this case, `OMSimSensitiveDetector::ProcessHits` will store all absorbed photons. The number of photons absorbed will depend on the absorption length of the material connected to the logical volume. * If there's a need to make a volume sensitive to particles other than photons, add a new entry to the `DetectorType` enum (in `OMSimSensitiveDetector.hh`) and incorporate a new method that handles this scenario in `OMSimSensitiveDetector::ProcessHits`. diff --git a/effective_area/OMSim_effective_area.cc b/effective_area/OMSim_effective_area.cc index 3940ce5f97..fda4a6b44e 100755 --- a/effective_area/OMSim_effective_area.cc +++ b/effective_area/OMSim_effective_area.cc @@ -78,9 +78,10 @@ int main(int pArgumentCount, char *pArgumentVector[]) bool lContinue = lSimulation.handleArguments(pArgumentCount, pArgumentVector); if (!lContinue) return 0; - OMSimEffectiveAreaDetector* lDetectorConstruction = new OMSimEffectiveAreaDetector(); - lSimulation.initialiseSimulation(lDetectorConstruction); - + std::unique_ptr lDetectorConstruction = std::make_unique(); + lSimulation.initialiseSimulation(lDetectorConstruction.get()); + lDetectorConstruction.release(); + runEffectiveAreaSimulation(); lSimulation.startVisualisationIfRequested(); diff --git a/effective_area/include/OMSimPrimaryGeneratorAction.hh b/effective_area/include/OMSimPrimaryGeneratorAction.hh index ac99466fd8..b47bb05066 100755 --- a/effective_area/include/OMSimPrimaryGeneratorAction.hh +++ b/effective_area/include/OMSimPrimaryGeneratorAction.hh @@ -15,7 +15,7 @@ public: void GeneratePrimaries(G4Event* anEvent); private: - G4GeneralParticleSource* lParticleSource; + G4GeneralParticleSource* mParticleSource; }; #endif diff --git a/effective_area/src/OMSimEffectiveAreaDetector.cc b/effective_area/src/OMSimEffectiveAreaDetector.cc index 518a817478..68555ede7b 100644 --- a/effective_area/src/OMSimEffectiveAreaDetector.cc +++ b/effective_area/src/OMSimEffectiveAreaDetector.cc @@ -6,7 +6,6 @@ #include "OMSimMDOM.hh" #include "OMSimCommandArgsTable.hh" #include "OMSimHitManager.hh" -#include "G4SDManager.hh" #include #include "OMSimSensitiveDetector.hh" #include "G4LogicalVolume.hh" diff --git a/effective_area/src/OMSimPrimaryGeneratorAction.cc b/effective_area/src/OMSimPrimaryGeneratorAction.cc index 3c3e120310..c4bd561262 100755 --- a/effective_area/src/OMSimPrimaryGeneratorAction.cc +++ b/effective_area/src/OMSimPrimaryGeneratorAction.cc @@ -2,20 +2,21 @@ #include #include - +#include OMSimPrimaryGeneratorAction::OMSimPrimaryGeneratorAction() { - lParticleSource = new G4GeneralParticleSource (); + mParticleSource = new G4GeneralParticleSource (); } OMSimPrimaryGeneratorAction::~OMSimPrimaryGeneratorAction() { - delete lParticleSource; + delete mParticleSource; } void OMSimPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { - lParticleSource->GeneratePrimaryVertex(anEvent); + mParticleSource->SetParticlePolarization(G4RandomDirection()); + mParticleSource->GeneratePrimaryVertex(anEvent); } diff --git a/radioactive_decays/include/OMSimPrimaryGeneratorAction.hh b/radioactive_decays/include/OMSimPrimaryGeneratorAction.hh index 8bef764921..564994a3f5 100755 --- a/radioactive_decays/include/OMSimPrimaryGeneratorAction.hh +++ b/radioactive_decays/include/OMSimPrimaryGeneratorAction.hh @@ -26,7 +26,7 @@ public: void GeneratePrimaries(G4Event* anEvent); private: - G4GeneralParticleSource* lParticleSource; + G4GeneralParticleSource* mParticleSource; }; #endif diff --git a/radioactive_decays/src/OMSimDecaysGPS.cc b/radioactive_decays/src/OMSimDecaysGPS.cc index 068d6fafa2..c0eea02778 100644 --- a/radioactive_decays/src/OMSimDecaysGPS.cc +++ b/radioactive_decays/src/OMSimDecaysGPS.cc @@ -18,7 +18,7 @@ G4String OMSimDecaysGPS::getDecayTerminationNuclide() */ void OMSimDecaysGPS::generalGPS() { - log_debug("Configuring general GPS"); + log_trace("Configuring general GPS"); OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance(); OMSimCommandArgsTable &lArgs = OMSimCommandArgsTable::getInstance(); @@ -41,13 +41,13 @@ void OMSimDecaysGPS::generalGPS() if (lArgs.get("scint_off")) { - log_debug("Inactivating scintillation process"); + log_trace("Inactivating scintillation process"); lUIinterface.applyCommand("/process/inactivate Scintillation"); } if (lArgs.get("cherenkov_off")) { - log_debug("Inactivating Cerenkov process"); + log_trace("Inactivating Cerenkov process"); lUIinterface.applyCommand("/process/inactivate Cerenkov"); } } @@ -62,7 +62,7 @@ void OMSimDecaysGPS::generalGPS() */ void OMSimDecaysGPS::configureIsotopeGPS(G4String Isotope, G4String pVolumeName) { - log_debug("Configuring GPS for isotope {} in volume {}", Isotope, pVolumeName); + log_trace("Configuring GPS for isotope {} in volume {}", Isotope, pVolumeName); OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance(); generalGPS(); @@ -84,7 +84,7 @@ void OMSimDecaysGPS::configureIsotopeGPS(G4String Isotope, G4String pVolumeName) std::map OMSimDecaysGPS::calculateNumberOfDecays(G4MaterialPropertiesTable *pMPT, G4double pTimeWindow, G4double pMass) { std::map mNumberDecays; - log_debug("Calculating number of decays from material isotope activity..."); + log_trace("Calculating number of decays from material isotope activity..."); for (auto &pair : mIsotopeCommands) { G4String lIsotope = pair.first; @@ -101,7 +101,7 @@ std::map OMSimDecaysGPS::calculateNumberOfDecays(G4MaterialProp */ void OMSimDecaysGPS::simulateDecaysInPressureVessel(G4double pTimeWindow) { - log_debug("Simulating radioactive decays in pressure vessel in a time window of {} seconds", pTimeWindow); + log_trace("Simulating radioactive decays in pressure vessel in a time window of {} seconds", pTimeWindow); OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance(); @@ -120,7 +120,7 @@ void OMSimDecaysGPS::simulateDecaysInPressureVessel(G4double pTimeWindow) mNuclideStopName = mTerminationIsotopes[pair.first]; configureIsotopeGPS(lIsotope, lPressureVesselName); - log_debug("Simulating {} decays of {} in pressure vessel", lNrDecays, lIsotope); + log_trace("Simulating {} decays of {} in pressure vessel", lNrDecays, lIsotope); lUIinterface.runBeamOn(lNrDecays); } } @@ -131,7 +131,7 @@ void OMSimDecaysGPS::simulateDecaysInPressureVessel(G4double pTimeWindow) */ void OMSimDecaysGPS::simulateDecaysInPMTs(G4double pTimeWindow) { - log_debug("Simulating radioactive decays in glass of PMTs in a time window of {} seconds", pTimeWindow); + log_trace("Simulating radioactive decays in glass of PMTs in a time window of {} seconds", pTimeWindow); OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance(); G4double lMass = mOM->getPMTmanager()->getPMTGlassWeight(); G4LogicalVolume *lPVVolume = mOM->getPMTmanager()->getLogicalVolume(); @@ -149,7 +149,7 @@ void OMSimDecaysGPS::simulateDecaysInPMTs(G4double pTimeWindow) configureIsotopeGPS(lIsotope, "PMT"+std::to_string(pmt)); - log_debug("Simulating {} decays of {} in PMT {}", lNrDecays, lIsotope, pmt); + log_trace("Simulating {} decays of {} in PMT {}", lNrDecays, lIsotope, pmt); lUIinterface.runBeamOn(lNrDecays); } } diff --git a/radioactive_decays/src/OMSimPrimaryGeneratorAction.cc b/radioactive_decays/src/OMSimPrimaryGeneratorAction.cc index 3238cf630d..d038183b80 100755 --- a/radioactive_decays/src/OMSimPrimaryGeneratorAction.cc +++ b/radioactive_decays/src/OMSimPrimaryGeneratorAction.cc @@ -7,18 +7,18 @@ OMSimPrimaryGeneratorAction::OMSimPrimaryGeneratorAction() { - lParticleSource = new G4GeneralParticleSource (); - lParticleSource->SetParticleDefinition(G4GenericIon::GenericIonDefinition()); + mParticleSource = new G4GeneralParticleSource (); + mParticleSource->SetParticleDefinition(G4GenericIon::GenericIonDefinition()); } OMSimPrimaryGeneratorAction::~OMSimPrimaryGeneratorAction() { - delete lParticleSource; + delete mParticleSource; } void OMSimPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { - lParticleSource->GeneratePrimaryVertex(anEvent); + mParticleSource->GeneratePrimaryVertex(anEvent); } diff --git a/scintillation_yield_determination/src/OMSimYieldDetector.cc b/scintillation_yield_determination/src/OMSimYieldDetector.cc new file mode 100644 index 0000000000..5dd72c61f7 --- /dev/null +++ b/scintillation_yield_determination/src/OMSimYieldDetector.cc @@ -0,0 +1,129 @@ +#include "OMSimYieldDetector.hh" +#include "OMSimYieldSetup.hh" +#include "OMSimMDOM.hh" +#include "OMSimPDOM.hh" +#include "OMSimLOM16.hh" +#include "OMSimLOM18.hh" +#include "OMSimDEGG.hh" +#include "OMSimCommandArgsTable.hh" +#include "OMSimHitManager.hh" +#include "OMSimSensitiveDetector.hh" + +/** + * @brief Constructs the world volume (sphere). + */ +void OMSimYieldDetector::constructWorld() +{ + mWorldSolid = new G4Orb("World", OMSimCommandArgsTable::getInstance().get("world_radius") * m); + mWorldLogical = new G4LogicalVolume(mWorldSolid, mData->getMaterial("argWorld"), "World_log", 0, 0, 0); + mWorldPhysical = new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), mWorldLogical, "World_phys", 0, false, 0); + G4VisAttributes *World_vis = new G4VisAttributes(G4Colour(0.45, 0.5, 0.35, 0.)); + mWorldLogical->SetVisAttributes(World_vis); +} + +/** + * @brief Constructs the selected detector from the command line argument. + */ +void OMSimYieldDetector::constructDetector() +{ + OMSimHitManager &lHitManager = OMSimHitManager::getInstance(); + OMSimCommandArgsTable &lArgs = OMSimCommandArgsTable::getInstance(); + bool lPlaceHarness = lArgs.get("place_harness"); + + OMSimOpticalModule *lOpticalModule = nullptr; + + switch (lArgs.get("detector_type")) + { + + case 0: + { + log_info("Constructing Okamoto Cs-137 Source setup for electron yield"); + G4double lZrandom = 0; + G4double lYrandom = 0; + G4double lZrandomSource = 0; + G4double lYrandomSource = 0; + G4double lXSource = 0; + G4double lXSample = 0; + G4double lPMTrotY = 1*deg; + G4double lPMTrotX = 0*deg; + if (lArgs.get("systematics")) + { + lZrandom = G4RandGauss::shoot(0, 0.02) * cm; + lYrandom = G4RandGauss::shoot(0, 2) * mm; + lZrandomSource = G4RandGauss::shoot(0, 0.01) * cm; + lYrandomSource = G4RandGauss::shoot(0, 2) * mm; + lXSource = G4RandGauss::shoot(0, 2) * mm; + lXSample = G4RandGauss::shoot(0, 2) * mm; + lPMTrotY = G4RandGauss::shoot(1, 0.1) * deg; + lPMTrotX = G4RandGauss::shoot(0, 0.1) * deg; + } + + OMSimPMTConstruction *lPMTManager = new OMSimPMTConstruction(mData); + lPMTManager->selectPMT("argPMT"); + lPMTManager->construction(); + lPMTManager->placeIt(G4ThreeVector(0, 0, 0), G4RotationMatrix().rotateY(lPMTrotY).rotateX(lPMTrotX), mWorldLogical, "_0"); + lHitManager.setNumberOfPMTs(1, 0); + lPMTManager->configureSensitiveVolume(this, "/PMT/0"); + + OkamotoLargeSample *lOkamotoSample = new OkamotoLargeSample(mData); + + G4double lzSample = 4.46 * cm + lPMTManager->getDistancePMTCenterToTip() + lZrandom; + G4double lySample = 21.77 * mm - (40 * mm - 25 * mm)+lYrandom; + + lOkamotoSample->placeIt(G4ThreeVector(lXSample, -lySample, lzSample), G4RotationMatrix(), mWorldLogical, ""); + + G4double lzSource = lzSample + lOkamotoSample->getSampleThickness() + 0.96 * cm + lZrandomSource; + Cs137Source *lSource = new Cs137Source(mData); + lSource->placeIt(G4ThreeVector(lXSource, -lySample+lYrandomSource, lzSource), G4RotationMatrix(), mWorldLogical, ""); + mSource = lSource; + + break; + } + case 1: + { + log_info("Constructing single PMT"); + OMSimPMTConstruction *lPMTManager = new OMSimPMTConstruction(mData); + lPMTManager->selectPMT("argPMT"); + lPMTManager->construction(); + lPMTManager->placeIt(G4ThreeVector(0, 0, 0), G4RotationMatrix(), mWorldLogical, "_0"); + lHitManager.setNumberOfPMTs(1, 0); + lPMTManager->configureSensitiveVolume(this, "/PMT/0"); + break; + } + case 2: + { + lOpticalModule = new mDOM(mData, lPlaceHarness); + break; + } + case 3: + { + + lOpticalModule = new pDOM(mData, lPlaceHarness); + break; + } + case 4: + { + + lOpticalModule = new LOM16(mData, lPlaceHarness); + break; + } + case 5: + { + + lOpticalModule = new LOM18(mData, lPlaceHarness); + break; + } + case 6: + { + lOpticalModule = new DEGG(mData, lPlaceHarness); + break; + } + } + + if (lOpticalModule) + { + lOpticalModule->placeIt(G4ThreeVector(0, 0, 0), G4RotationMatrix(), mWorldLogical, ""); + lOpticalModule->configureSensitiveVolume(this); + mOpticalModule = lOpticalModule; + } +}