From 63e759cedf191a198c7283e25fad3fbd514a6150 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Sat, 25 Jan 2025 01:14:50 +0100 Subject: [PATCH] migrate musun reader to RMGAnalysisReader --- include/RMGAnalysisReader.hh | 15 ++- include/RMGGeneratorMUSUNCosmicMuons.hh | 10 +- src/RMGAnalysisReader.cc | 45 +++++--- src/RMGGeneratorMUSUNCosmicMuons.cc | 130 ++++++++++++------------ 4 files changed, 113 insertions(+), 87 deletions(-) diff --git a/include/RMGAnalysisReader.hh b/include/RMGAnalysisReader.hh index 820026da..a50d380b 100644 --- a/include/RMGAnalysisReader.hh +++ b/include/RMGAnalysisReader.hh @@ -134,9 +134,16 @@ class RMGAnalysisReader final { * "dir/table" this is "dir". * @param ntuple_name the first part of the input table name. For the table addressed by * "dir/table" this is "table". + * @param lock a lock instance obtained by @ref RMGAnalysisReader::GetLock - optional. + * @param force_ext force a file extension/reader type. should be a lowercase file extension + * like @c lh5, @c hdf5, @c csv, @c root. */ - [[nodiscard]] Access OpenFile(std::string& file_name, std::string ntuple_dir_name, - std::string ntuple_name); + [[nodiscard]] Access OpenFile(const std::string& file_name, std::string ntuple_dir_name, + std::string ntuple_name, G4AutoLock lock, std::string force_ext = ""); + /** + * @overload */ + [[nodiscard]] Access OpenFile(const std::string& file_name, std::string ntuple_dir_name, + std::string ntuple_name, std::string force_ext = ""); /** * @brief if any file is open for reading, close the reader. @@ -152,6 +159,10 @@ class RMGAnalysisReader final { * @details this acquires a global lock to avoid problems. The lock is held until the access handle is discarded. */ [[nodiscard]] Access GetLockedReader() const; + /** + * @brief acquires a global lock to the analysis reader mutex. */ + [[nodiscard]] G4AutoLock GetLock() const; + /** * @brief get the file name of the current open file, or an empty string. */ [[nodiscard]] inline auto& GetFileName() const { return fFileName; } diff --git a/include/RMGGeneratorMUSUNCosmicMuons.hh b/include/RMGGeneratorMUSUNCosmicMuons.hh index cdfdd376..63896e4f 100644 --- a/include/RMGGeneratorMUSUNCosmicMuons.hh +++ b/include/RMGGeneratorMUSUNCosmicMuons.hh @@ -4,10 +4,10 @@ #include #include "CLHEP/Units/SystemOfUnits.h" -#include "G4CsvAnalysisReader.hh" #include "G4GenericMessenger.hh" #include "G4ParticleGun.hh" +#include "RMGAnalysisReader.hh" #include "RMGVGenerator.hh" #include "RMGVVertexGenerator.hh" @@ -51,17 +51,17 @@ class RMGGeneratorMUSUNCosmicMuons : public RMGVGenerator { void DefineCommands(); void SetMUSUNFile(G4String pathToFile); - void PrepareCopy(G4String pathToFile); + void PrepareCopy(std::string pathToFile); std::unique_ptr fGun = nullptr; std::unique_ptr fMessenger = nullptr; G4String fPathToFile = ""; std::filesystem::path fPathToTmpFolder; - G4String fPathToTmpFile = ""; + std::filesystem::path fPathToTmpFile; - static G4CsvAnalysisReader* fAnalysisReader; + static RMGAnalysisReader* fAnalysisReader; - static RMGGeneratorMUSUNCosmicMuons_Data* input_data; + static RMGGeneratorMUSUNCosmicMuons_Data* fInputData; }; #endif diff --git a/src/RMGAnalysisReader.cc b/src/RMGAnalysisReader.cc index e795a495..78f8f93d 100644 --- a/src/RMGAnalysisReader.cc +++ b/src/RMGAnalysisReader.cc @@ -35,13 +35,8 @@ namespace fs = std::filesystem; G4Mutex RMGAnalysisReader::fMutex = G4MUTEX_INITIALIZER; -RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, - std::string ntuple_dir_name, std::string ntuple_name) { - - // create an invalid lock to return on error. - G4AutoLock invalid_lock(&fMutex, std::defer_lock); - invalid_lock.release(); - auto invalid_access = RMGAnalysisReader::Access(std::move(invalid_lock), nullptr, -1, nullptr); +RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(const std::string& file_name, + std::string ntuple_dir_name, std::string ntuple_name, std::string force_ext) { // reader initialization should only happen on the master thread (otherwise it will fail). if (!G4Threading::IsMasterThread()) { @@ -51,19 +46,42 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, G4AutoLock lock(&fMutex); + return OpenFile(file_name, ntuple_dir_name, ntuple_name, std::move(lock), force_ext); +} + +RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(const std::string& file_name, + std::string ntuple_dir_name, std::string ntuple_name, G4AutoLock lock, std::string force_ext) { + + if (lock.mutex() != &fMutex || !lock) { + RMGLog::OutDev(RMGLog::fatal, "used lock protecting the wrong mutex or unheld lock"); + std::abort(); + } + + // reader initialization should only happen on the master thread (otherwise it will fail). + if (!G4Threading::IsMasterThread()) { + RMGLog::OutDev(RMGLog::fatal, "can only be used on the master thread"); + std::abort(); + } + + // create an invalid lock to return on error. + G4AutoLock invalid_lock(&fMutex, std::defer_lock); + invalid_lock.release(); + auto invalid_access = RMGAnalysisReader::Access(std::move(invalid_lock), nullptr, -1, nullptr); + if (fReader) return std::move(invalid_access); fFileIsTemp = false; + fFileName = file_name; + auto path = fs::path(file_name); - if (!fs::exists(fs::path(file_name))) { + if (!fs::exists(path)) { RMGLog::Out(RMGLog::error, "input file ", file_name, " does not exist"); return std::move(invalid_access); } std::map> units_map{}; - // NOTE: GetExtension() returns a default extension if there is no file extension - auto ext = G4Analysis::GetExtension(file_name); + auto ext = !force_ext.empty() ? force_ext : G4Analysis::GetExtension(file_name); if (ext == "root") { fReader = G4RootAnalysisReader::Instance(); } else if (ext == "hdf5") { @@ -77,7 +95,6 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, #if RMG_HAS_HDF5 std::uniform_int_distribution dist(10000, 99999); std::random_device rd; - auto path = fs::path(file_name); std::string new_fn = ".rmg-vtx-" + std::to_string(dist(rd)) + "." + path.stem().string() + ".hdf5"; @@ -95,7 +112,7 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, } fHasUnits = true; - file_name = new_fn; + fFileName = new_fn; fReader = G4Hdf5AnalysisReader::Instance(); fFileIsTemp = true; #else @@ -112,9 +129,7 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, if (RMGLog::GetLogLevel() <= RMGLog::debug) fReader->SetVerboseLevel(10); - fFileName = file_name; - - fNtupleId = fReader->GetNtuple(ntuple_name, file_name, ntuple_dir_name); + fNtupleId = fReader->GetNtuple(ntuple_name, fFileName, ntuple_dir_name); if (fNtupleId < 0) { RMGLog::Out(RMGLog::error, "Ntuple named '", ntuple_name, "' could not be found in input file!"); return std::move(invalid_access); diff --git a/src/RMGGeneratorMUSUNCosmicMuons.cc b/src/RMGGeneratorMUSUNCosmicMuons.cc index 039f67a3..95bc776e 100644 --- a/src/RMGGeneratorMUSUNCosmicMuons.cc +++ b/src/RMGGeneratorMUSUNCosmicMuons.cc @@ -4,29 +4,20 @@ #include #include -#include "G4AutoLock.hh" -#include "G4CsvAnalysisReader.hh" #include "G4GenericMessenger.hh" #include "G4ParticleGun.hh" -#include "G4ParticleMomentum.hh" -#include "G4ParticleTypes.hh" +#include "G4ParticleTable.hh" #include "G4ThreeVector.hh" #include "Randomize.hh" #include "RMGHardware.hh" #include "RMGLog.hh" -#include "RMGManager.hh" -#include "RMGTools.hh" #include "RMGVGenerator.hh" namespace u = CLHEP; -namespace { - G4Mutex RMGGeneratorMUSUNCosmicMuonsDestrMutex = G4MUTEX_INITIALIZER; - G4Mutex RMGGeneratorMUSUNCosmicMuonsMutex = G4MUTEX_INITIALIZER; -} // namespace -G4CsvAnalysisReader* RMGGeneratorMUSUNCosmicMuons::fAnalysisReader = nullptr; -RMGGeneratorMUSUNCosmicMuons_Data* RMGGeneratorMUSUNCosmicMuons::input_data = nullptr; +RMGAnalysisReader* RMGGeneratorMUSUNCosmicMuons::fAnalysisReader = new RMGAnalysisReader(); +RMGGeneratorMUSUNCosmicMuons_Data* RMGGeneratorMUSUNCosmicMuons::fInputData = nullptr; RMGGeneratorMUSUNCosmicMuons::RMGGeneratorMUSUNCosmicMuons() : RMGVGenerator("MUSUNCosmicMuons") { this->DefineCommands(); @@ -34,7 +25,7 @@ RMGGeneratorMUSUNCosmicMuons::RMGGeneratorMUSUNCosmicMuons() : RMGVGenerator("MU fPathToTmpFolder = std::filesystem::temp_directory_path(); } -void RMGGeneratorMUSUNCosmicMuons::PrepareCopy(G4String pathToFile) { +void RMGGeneratorMUSUNCosmicMuons::PrepareCopy(std::string pathToFile) { /* The working assumption is that the user uses the output directly from MUSUN, i.e. there is no header. To allow proper multiprocessing, we want the file to be read using G4CsvAnalysisReader. @@ -43,24 +34,23 @@ void RMGGeneratorMUSUNCosmicMuons::PrepareCopy(G4String pathToFile) { */ // Define fPathToTmpFile - std::filesystem::path originalFilePath((std::string)pathToFile); + std::filesystem::path originalFilePath(pathToFile); std::filesystem::path fileName = originalFilePath.filename(); - fPathToTmpFile = (G4String)(fPathToTmpFolder / - fileName); //.substr(0, fileName.find_last_of(".")) + "_nt_MUSUN.csv"; + fPathToTmpFile = fPathToTmpFolder / fileName; // Check if the original file exists / the tmp file does not exist std::ifstream originalFile(pathToFile); - if (!originalFile) RMGLog::Out(RMGLog::fatal, "MUSUN file not found! Exit."); + if (!originalFile) RMGLog::Out(RMGLog::fatal, "MUSUN file ", pathToFile, " not found! Exit."); - if (std::filesystem::exists((std::string)fPathToTmpFile)) { + if (std::filesystem::exists(fPathToTmpFile)) { RMGLog::Out(RMGLog::warning, "Temporary file already exists. Deleting it."); - std::filesystem::remove((std::string)fPathToTmpFile); + std::filesystem::remove(fPathToTmpFile); } // Counting the number of columns to identify which header to use std::string firstLine; if (!std::getline(originalFile, firstLine)) { - std::cerr << "Error: File is empty" << std::endl; + RMGLog::Out(RMGLog::error, "Error: File is empty"); return; } std::istringstream iss(firstLine); @@ -110,76 +100,86 @@ void RMGGeneratorMUSUNCosmicMuons::PrepareCopy(G4String pathToFile) { void RMGGeneratorMUSUNCosmicMuons::BeginOfRunAction(const G4Run*) { - G4AutoLock lock(&RMGGeneratorMUSUNCosmicMuonsMutex); - if (!fAnalysisReader) { + if (!fInputData) { + if (!G4Threading::IsMasterThread()) { + RMGLog::Out(RMGLog::fatal, "on worker, but data not initialized"); + std::abort(); + return; + } + + // include in lock + auto lock = fAnalysisReader->GetLock(); PrepareCopy(fPathToFile); - using G4AnalysisReader = G4CsvAnalysisReader; - fAnalysisReader = G4AnalysisReader::Instance(); - fAnalysisReader->SetVerboseLevel(1); - fAnalysisReader->SetFileName(fPathToTmpFile); - G4int ntupleId = fAnalysisReader->GetNtuple("MUSUN", fPathToTmpFile); - if (ntupleId < 0) RMGLog::Out(RMGLog::fatal, "Temp MUSUN file not found! Exit."); - - input_data = new RMGGeneratorMUSUNCosmicMuons_Data; - fAnalysisReader->SetNtupleIColumn(0, "ID", (input_data->fID)); - fAnalysisReader->SetNtupleIColumn(0, "type", (input_data->fType)); - fAnalysisReader->SetNtupleDColumn(0, "Ekin", (input_data->fEkin)); - fAnalysisReader->SetNtupleDColumn(0, "x", (input_data->fX)); - fAnalysisReader->SetNtupleDColumn(0, "y", (input_data->fY)); - fAnalysisReader->SetNtupleDColumn(0, "z", (input_data->fZ)); - fAnalysisReader->SetNtupleDColumn(0, "theta", (input_data->fTheta)); - fAnalysisReader->SetNtupleDColumn(0, "phi", (input_data->fPhi)); - fAnalysisReader->SetNtupleDColumn(0, "px", (input_data->fPx)); - fAnalysisReader->SetNtupleDColumn(0, "py", (input_data->fPy)); - fAnalysisReader->SetNtupleDColumn(0, "pz", (input_data->fPz)); + + auto reader = fAnalysisReader->OpenFile(fPathToTmpFile, "", "MUSUN", std::move(lock), "csv"); + if (!reader) RMGLog::Out(RMGLog::fatal, "Temp MUSUN file not found! Exit."); + + fInputData = new RMGGeneratorMUSUNCosmicMuons_Data; + reader.SetNtupleIColumn("ID", (fInputData->fID)); + reader.SetNtupleIColumn("type", (fInputData->fType)); + reader.SetNtupleDColumn("Ekin", (fInputData->fEkin)); + reader.SetNtupleDColumn("x", (fInputData->fX)); + reader.SetNtupleDColumn("y", (fInputData->fY)); + reader.SetNtupleDColumn("z", (fInputData->fZ)); + reader.SetNtupleDColumn("theta", (fInputData->fTheta)); + reader.SetNtupleDColumn("phi", (fInputData->fPhi)); + reader.SetNtupleDColumn("px", (fInputData->fPx)); + reader.SetNtupleDColumn("py", (fInputData->fPy)); + reader.SetNtupleDColumn("pz", (fInputData->fPz)); } - lock.unlock(); } void RMGGeneratorMUSUNCosmicMuons::EndOfRunAction(const G4Run*) { - G4AutoLock lock(&RMGGeneratorMUSUNCosmicMuonsDestrMutex); - - if (fAnalysisReader) { - std::filesystem::remove((std::string)fPathToTmpFile); - // delete fAnalysisReader; - // fAnalysisReader = 0; - // delete input_data; - // input_data = 0; - } + if (!G4Threading::IsMasterThread()) return; + + auto reader = fAnalysisReader->GetLockedReader(); + if (reader) { std::filesystem::remove(fPathToTmpFile); } } void RMGGeneratorMUSUNCosmicMuons::GeneratePrimaries(G4Event* event) { - G4AutoLock lock(&RMGGeneratorMUSUNCosmicMuonsMutex); + auto reader = fAnalysisReader->GetLockedReader(); + + if (!reader) { + RMGLog::Out(RMGLog::error, "MUSUN ntuple could not be found in input file!"); + return; + } + + if (!reader.GetNtupleRow()) { + RMGLog::Out(RMGLog::error, "no more input rows in MUSUN file!"); + return; + } - fAnalysisReader->GetNtupleRow(); + // copy data and end critical section. + RMGGeneratorMUSUNCosmicMuons_Data input_data = *fInputData; + reader.unlock(); - G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable(); - if (input_data->fType == 10) fGun->SetParticleDefinition(theParticleTable->FindParticle("mu-")); + auto theParticleTable = G4ParticleTable::GetParticleTable(); + if (input_data.fType == 10) fGun->SetParticleDefinition(theParticleTable->FindParticle("mu-")); else fGun->SetParticleDefinition(theParticleTable->FindParticle("mu+")); RMGLog::OutFormat(RMGLog::debug, "...origin ({:.4g}, {:.4g}, {:.4g}) m", - input_data->fX * u::cm / u::m, input_data->fY * u::cm / u::m, input_data->fZ * u::cm / u::m); - fGun->SetParticlePosition({input_data->fX * u::cm, input_data->fY * u::cm, input_data->fZ * u::cm}); + input_data.fX * u::cm / u::m, input_data.fY * u::cm / u::m, input_data.fZ * u::cm / u::m); + fGun->SetParticlePosition({input_data.fX * u::cm, input_data.fY * u::cm, input_data.fZ * u::cm}); - if (input_data->fTheta != 0 && input_data->fPhi != 0) { + if (input_data.fTheta != 0 && input_data.fPhi != 0) { G4ThreeVector d_cart(1, 1, 1); - d_cart.setTheta(input_data->fTheta); // in rad - d_cart.setPhi(input_data->fPhi); // in rad + d_cart.setTheta(input_data.fTheta); // in rad + d_cart.setPhi(input_data.fPhi); // in rad d_cart.setMag(1 * u::m); fGun->SetParticleMomentumDirection(d_cart); RMGLog::OutFormat(RMGLog::debug, "...direction (θ,φ) = ({:.4g}, {:.4g}) deg", - input_data->fTheta / u::deg, input_data->fPhi / u::deg); + input_data.fTheta / u::deg, input_data.fPhi / u::deg); } else { - G4ThreeVector d_cart(input_data->fPx, input_data->fPy, input_data->fPz); + G4ThreeVector d_cart(input_data.fPx, input_data.fPy, input_data.fPz); fGun->SetParticleMomentumDirection(d_cart); RMGLog::OutFormat(RMGLog::debug, "...direction (px,py,pz) = ({:.4g}, {:.4g}, {:.4g}) deg", - input_data->fPx, input_data->fPy, input_data->fPz); + input_data.fPx, input_data.fPy, input_data.fPz); } - RMGLog::OutFormat(RMGLog::debug, "...energy {:.4g} GeV", input_data->fEkin); - fGun->SetParticleEnergy(input_data->fEkin * u::GeV); + RMGLog::OutFormat(RMGLog::debug, "...energy {:.4g} GeV", input_data.fEkin); + fGun->SetParticleEnergy(input_data.fEkin * u::GeV); fGun->GeneratePrimaryVertex(event); }