From a0f18bdc9ac0ed360d85cd039ab48c852acc1c95 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Tue, 21 Jan 2025 15:10:49 +0100 Subject: [PATCH 1/9] more thread-safety in RMGAnalysisReader --- include/RMGAnalysisReader.hh | 40 ++++++++++++++++++++++++++++++++++-- include/RMGVertexFromFile.hh | 1 - src/RMGAnalysisReader.cc | 22 +++++++++++--------- src/RMGVertexFromFile.cc | 14 +++++++------ 4 files changed, 58 insertions(+), 19 deletions(-) diff --git a/include/RMGAnalysisReader.hh b/include/RMGAnalysisReader.hh index b4d49a8..3e55c9a 100644 --- a/include/RMGAnalysisReader.hh +++ b/include/RMGAnalysisReader.hh @@ -23,7 +23,21 @@ #include "RMGVVertexGenerator.hh" -class RMGAnalysisReader { +/** + * @brief wrapper around \ref(G4VAnalasisReader) instances with special handling for LH5 files. + * + * @details notes for threadsafe use: + * - opening/closing can only be performed on the master thread. + * - in a multithreaded application, all function calls need to by guarded by a mutex. Worker + * threads can use the reader + * - instance after opening, but only one worker can use the reader at a time. + * - the instance obtained by \ref(GetReader) should only be bound once to variables that are of same + * storage duration as the \ref(RMGAnalysisReader) instance. Example: If you use a static instance + * in a class, only bind to static class fields (of the same class) to read the values into. + * - it is EXTREMELY important to always use overloads with the ntuple id parameter, and to always + * pass the result of \ref(GetNtupleID) to them. This affects ntuple reading and column registration. + */ +class RMGAnalysisReader final { public: @@ -35,11 +49,33 @@ class RMGAnalysisReader { RMGAnalysisReader(RMGAnalysisReader&&) = delete; RMGAnalysisReader& operator=(RMGAnalysisReader&&) = delete; - bool OpenFile(std::string& name, std::string ntuple_dir_name, std::string ntuple_name); + /** + * @brief open an input file for reading of one specific ntuple. + * + * @details This function can only be used on the master thread. + * + * @param file_name the input file name. the file format is determined from the file extension. + * @param ntuple_dir_name the first part of the input table name. For the table addressed by + * "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". + */ + bool OpenFile(std::string& file_name, std::string ntuple_dir_name, std::string ntuple_name); + /** + * @brief if any file is open for reading, close the reader. + * + * @details this invalidates all reader instances obtained via \ref(GetReader). This function + * can only be used on the master thread. */ void CloseFile(); + /** + * @brief get a pointer to the current underlying G4VAnalysisReader. */ [[nodiscard]] inline auto GetReader() { return fReader; } + /** + * @brief return the ntuple id of the current tuple. If the value is negative, no ntuple is opened. */ [[nodiscard]] inline auto GetNtupleId() const { return fNtupleId; } + /** + * @brief get the file name of the current open file, or an empty string. */ [[nodiscard]] inline auto& GetFileName() const { return fFileName; } private: diff --git a/include/RMGVertexFromFile.hh b/include/RMGVertexFromFile.hh index d645684..afa8889 100644 --- a/include/RMGVertexFromFile.hh +++ b/include/RMGVertexFromFile.hh @@ -25,7 +25,6 @@ #include "RMGAnalysisReader.hh" #include "RMGVVertexGenerator.hh" -class G4VAnalysisReader; class RMGVertexFromFile : public RMGVVertexGenerator { public: diff --git a/src/RMGAnalysisReader.cc b/src/RMGAnalysisReader.cc index 53f7de5..1cfd9f5 100644 --- a/src/RMGAnalysisReader.cc +++ b/src/RMGAnalysisReader.cc @@ -35,7 +35,7 @@ namespace fs = std::filesystem; #include "RMGLog.hh" -bool RMGAnalysisReader::OpenFile(std::string& name, std::string ntuple_dir_name, +bool RMGAnalysisReader::OpenFile(std::string& file_name, std::string ntuple_dir_name, std::string ntuple_name) { // reader initialization should only happen on the master thread (otherwise it will fail). @@ -45,13 +45,13 @@ bool RMGAnalysisReader::OpenFile(std::string& name, std::string ntuple_dir_name, fFileIsTemp = false; - if (!fs::exists(fs::path(name))) { - RMGLog::Out(RMGLog::error, "input file ", name, " does not exist"); + if (!fs::exists(fs::path(file_name))) { + RMGLog::Out(RMGLog::error, "input file ", file_name, " does not exist"); return false; } // NOTE: GetExtension() returns a default extension if there is no file extension - auto ext = G4Analysis::GetExtension(name); + auto ext = G4Analysis::GetExtension(file_name); if (ext == "root") { fReader = G4RootAnalysisReader::Instance(); } else if (ext == "hdf5") { @@ -65,13 +65,13 @@ bool RMGAnalysisReader::OpenFile(std::string& name, std::string ntuple_dir_name, #if RMG_HAS_HDF5 std::uniform_int_distribution dist(10000, 99999); std::random_device rd; - auto path = fs::path(name); + auto path = fs::path(file_name); std::string new_fn = ".rmg-vtx-" + std::to_string(dist(rd)) + "." + path.stem().string() + ".hdf5"; std::error_code ec; if (!fs::copy_file(path, fs::path(new_fn)), ec) { - RMGLog::Out(RMGLog::error, "copy of input file ", name, " failed. Does it exist? (", + RMGLog::Out(RMGLog::error, "copy of input file ", file_name, " failed. Does it exist? (", ec.message(), ")"); return false; } @@ -82,7 +82,7 @@ bool RMGAnalysisReader::OpenFile(std::string& name, std::string ntuple_dir_name, return false; } - name = new_fn; + file_name = new_fn; fReader = G4Hdf5AnalysisReader::Instance(); fFileIsTemp = true; #else @@ -99,10 +99,9 @@ bool RMGAnalysisReader::OpenFile(std::string& name, std::string ntuple_dir_name, if (RMGLog::GetLogLevel() <= RMGLog::debug) fReader->SetVerboseLevel(10); - fReader->SetFileName(name); - fFileName = name; + fFileName = file_name; - fNtupleId = fReader->GetNtuple(ntuple_name); + fNtupleId = fReader->GetNtuple(ntuple_name, file_name); if (fNtupleId < 0) { RMGLog::Out(RMGLog::error, "Ntuple named '", ntuple_name, "' could not be found in input file!"); return false; @@ -115,6 +114,8 @@ void RMGAnalysisReader::CloseFile() { if (!fReader) return; + // fReader is a thread-local singleton. Do not delete it here, otherwise geant4 will to delete it + // again. also do not close files, as there is no way to close a specific file only. fReader = nullptr; if (fFileIsTemp) { std::error_code ec; @@ -122,4 +123,5 @@ void RMGAnalysisReader::CloseFile() { } fFileName = ""; fFileIsTemp = false; + fNtupleId = -1; } diff --git a/src/RMGVertexFromFile.cc b/src/RMGVertexFromFile.cc index bb30d6f..5e51b36 100644 --- a/src/RMGVertexFromFile.cc +++ b/src/RMGVertexFromFile.cc @@ -35,11 +35,12 @@ void RMGVertexFromFile::OpenFile(std::string& name) { if (!fReader->OpenFile(name, fNtupleDirectoryName, "pos")) return; - if (fReader->GetNtupleId() >= 0) { + auto nt = fReader->GetNtupleId(); + if (nt >= 0) { // bind the static variables once here, and only use them later. - fReader->GetReader()->SetNtupleDColumn("xloc_in_m", fXpos); - fReader->GetReader()->SetNtupleDColumn("yloc_in_m", fYpos); - fReader->GetReader()->SetNtupleDColumn("zloc_in_m", fZpos); + fReader->GetReader()->SetNtupleDColumn(nt, "xloc_in_m", fXpos); + fReader->GetReader()->SetNtupleDColumn(nt, "yloc_in_m", fYpos); + fReader->GetReader()->SetNtupleDColumn(nt, "zloc_in_m", fZpos); } } @@ -48,9 +49,10 @@ bool RMGVertexFromFile::GenerateVertex(G4ThreeVector& vertex) { G4AutoLock lock(&fMutex); if (fReader->GetReader() && fReader->GetNtupleId() >= 0) { - fXpos = fYpos = fZpos = NAN; + fXpos = fYpos = fZpos = NAN; // initialize sentinel values. - if (fReader->GetReader()->GetNtupleRow()) { + auto nt = fReader->GetNtupleId(); + if (fReader->GetReader()->GetNtupleRow(nt)) { // check for NaN sentinel values - i.e. non-existing columns (there is no error message). if (std::isnan(fXpos) || std::isnan(fYpos) || std::isnan(fZpos)) { RMGLog::Out(RMGLog::error, "At least one of the columns does not exist"); From 78af7adc64cf2cf2f1a2740c023c4a31b7db6299 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Tue, 21 Jan 2025 15:12:28 +0100 Subject: [PATCH 2/9] do set the input ntuple directory correctly --- docs/rmg-commands.md | 2 +- src/RMGAnalysisReader.cc | 2 +- src/RMGConvertLH5.cc | 4 ---- src/RMGVertexFromFile.cc | 2 +- tests/vertex/macros/vert.hdf5 | Bin 19896 -> 19896 bytes 5 files changed, 3 insertions(+), 7 deletions(-) diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index 9b24129..d82980e 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -1308,7 +1308,7 @@ Set name of the file containing vertex positions for the next run. See the docum Change the default input directory/group for ntuples. :::{note} -this option only has an effect for LH5 input files. +this option only has an effect for LH5 or HDF5 input files. ::: * **Parameter** – `nt_directory` diff --git a/src/RMGAnalysisReader.cc b/src/RMGAnalysisReader.cc index 1cfd9f5..f9a45ad 100644 --- a/src/RMGAnalysisReader.cc +++ b/src/RMGAnalysisReader.cc @@ -101,7 +101,7 @@ bool RMGAnalysisReader::OpenFile(std::string& file_name, std::string ntuple_dir_ fFileName = file_name; - fNtupleId = fReader->GetNtuple(ntuple_name, file_name); + fNtupleId = fReader->GetNtuple(ntuple_name, file_name, ntuple_dir_name); if (fNtupleId < 0) { RMGLog::Out(RMGLog::error, "Ntuple named '", ntuple_name, "' could not be found in input file!"); return false; diff --git a/src/RMGConvertLH5.cc b/src/RMGConvertLH5.cc index d5e3153..9fcee71 100644 --- a/src/RMGConvertLH5.cc +++ b/src/RMGConvertLH5.cc @@ -454,10 +454,6 @@ bool RMGConvertLH5::ConvertFromLH5Internal() { if (ntuples_group.attrExists("datatype")) ntuples_group.removeAttr("datatype"); - if (ntuple_group_name != "default_ntuples") { - hfile.moveLink(ntuple_group_name, "default_ntuples"); - } - hfile.close(); LH5Log(RMGLog::summary, "Done updating HDF5 file ", fHdf5FileName, " to LH5"); diff --git a/src/RMGVertexFromFile.cc b/src/RMGVertexFromFile.cc index 5e51b36..0f266ed 100644 --- a/src/RMGVertexFromFile.cc +++ b/src/RMGVertexFromFile.cc @@ -105,7 +105,7 @@ void RMGVertexFromFile::DefineCommands() { fMessenger->DeclareProperty("NtupleDirectory", fNtupleDirectoryName) .SetGuidance("Change the default input directory/group for ntuples.") - .SetGuidance("note: this option only has an effect for LH5 input files.") + .SetGuidance("note: this option only has an effect for LH5 or HDF5 input files.") .SetParameterName("nt_directory", false) .SetDefaultValue(fNtupleDirectoryName) .SetStates(G4State_PreInit, G4State_Idle); diff --git a/tests/vertex/macros/vert.hdf5 b/tests/vertex/macros/vert.hdf5 index 7b3670f54b7c3bab5cc1bb518d5838f80bb9b5cb..c53b6aebb93c36adcda2277b70f04ab478898ce2 100644 GIT binary patch delta 81 zcmdlnn{mf%#t9o3EjDgsXJXWtypzdX(SZR11fXks= dfr(LLvL&1PWC1pfiQ7+1ejv!P*?_&n4FI|{7z_Xa From 4deeaf19f7e03ea48b9ec1f39468c9ab66e7198a Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Tue, 21 Jan 2025 16:57:32 +0100 Subject: [PATCH 3/9] tests: generate vertex input as fixture --- tests/vertex/CMakeLists.txt | 15 +++++-- tests/vertex/create-test-input.py | 42 ++++++++++++++++++++ tests/vertex/macros/create-test-vert-lh5.py | 16 -------- tests/vertex/macros/vert-hdf5.mac | 2 +- tests/vertex/macros/vert-lh5.mac | 2 +- tests/vertex/macros/vert.hdf5 | Bin 19896 -> 0 bytes tests/vertex/macros/vert.lh5 | Bin 24576 -> 0 bytes 7 files changed, 56 insertions(+), 21 deletions(-) create mode 100755 tests/vertex/create-test-input.py delete mode 100755 tests/vertex/macros/create-test-vert-lh5.py delete mode 100644 tests/vertex/macros/vert.hdf5 delete mode 100644 tests/vertex/macros/vert.lh5 diff --git a/tests/vertex/CMakeLists.txt b/tests/vertex/CMakeLists.txt index 4b16e00..b89e5c0 100644 --- a/tests/vertex/CMakeLists.txt +++ b/tests/vertex/CMakeLists.txt @@ -2,19 +2,28 @@ file( GLOB _aux RELATIVE ${PROJECT_SOURCE_DIR} - macros/*.mac macros/*.hdf5 macros/*.lh5 gdml/*.gdml gdml/*.xml) + macros/*.mac *.py gdml/*.gdml gdml/*.xml) + +set(_remage_from_lh5 "false") +if(RMG_HAS_HDF5) + set(_remage_from_lh5 $) +endif() # copy them to the build area foreach(_file ${_aux}) configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) endforeach() +add_test(NAME vertex/gen-input COMMAND ${PYTHONPATH} create-test-input.py ${_remage_from_lh5}) +set_tests_properties(vertex/gen-input PROPERTIES LABELS extra FIXTURES_SETUP vertex-input-fixture) + set(_macros vert-hdf5.mac vert-lh5.mac) foreach(_mac ${_macros}) add_test(NAME vertex/${_mac} COMMAND ${REMAGE_PYEXE} -g gdml/geometry.gdml -- macros/${_mac}) - set_tests_properties(vertex/${_mac} PROPERTIES LABELS extra) + set_tests_properties(vertex/${_mac} PROPERTIES LABELS extra FIXTURES_SETUP vertex-input-fixture) add_test(NAME vertex-mt/${_mac} COMMAND ${REMAGE_PYEXE} -g gdml/geometry.gdml -t 2 macros/${_mac}) - set_tests_properties(vertex-mt/${_mac} PROPERTIES LABELS "mt extra") + set_tests_properties(vertex-mt/${_mac} PROPERTIES LABELS "mt extra" FIXTURES_SETUP + vertex-input-fixture) endforeach() diff --git a/tests/vertex/create-test-input.py b/tests/vertex/create-test-input.py new file mode 100755 index 0000000..a2bf22c --- /dev/null +++ b/tests/vertex/create-test-input.py @@ -0,0 +1,42 @@ +#!/bin/env python3 +from __future__ import annotations + +import shutil +import subprocess +import sys +from pathlib import Path + +import numpy as np +from lgdo import Array, Table, lh5 + +rmg_from_lh5 = sys.argv[1] + + +def _convert_lh5_to_hdf5(lh5_fn): + if rmg_from_lh5 == "false": + return + + orig_path = Path(lh5_fn) + new_path = orig_path.with_suffix(".hdf5") + shutil.copy(orig_path, new_path) + ret = subprocess.run( + [rmg_from_lh5, "--ntuple-group", "vtx", "--", new_path], check=False + ) + assert ret.returncode == 0 + + +# ================================================== +# position input + +xloc = Array(np.array([0, 1, 2] * 100, dtype=np.float64), attrs={"units": "m"}) +yloc = Array(np.array([3, 4, 5] * 100, dtype=np.float64), attrs={"units": "m"}) +zloc = Array(np.array([6, 7, 8] * 100, dtype=np.float64), attrs={"units": "m"}) + +vtx_pos_file = "macros/vtx-pos.lh5" +lh5.write( + Table({"xloc": xloc, "yloc": yloc, "zloc": zloc}), + "vtx/pos", + lh5_file=vtx_pos_file, + wo_mode="of", +) +_convert_lh5_to_hdf5(vtx_pos_file) diff --git a/tests/vertex/macros/create-test-vert-lh5.py b/tests/vertex/macros/create-test-vert-lh5.py deleted file mode 100755 index d6ab8fe..0000000 --- a/tests/vertex/macros/create-test-vert-lh5.py +++ /dev/null @@ -1,16 +0,0 @@ -#!/bin/env python3 -from __future__ import annotations - -import numpy as np -from lgdo import Array, Table, lh5 - -xloc = Array(np.array([0, 1, 2] * 10, dtype=np.float64), attrs={"units": "m"}) -yloc = Array(np.array([3, 4, 5] * 10, dtype=np.float64), attrs={"units": "m"}) -zloc = Array(np.array([6, 7, 8] * 10, dtype=np.float64), attrs={"units": "m"}) - -lh5.write( - Table({"xloc": xloc, "yloc": yloc, "zloc": zloc}), - "vtx/pos", - lh5_file="vert.lh5", - wo_mode="of", -) diff --git a/tests/vertex/macros/vert-hdf5.mac b/tests/vertex/macros/vert-hdf5.mac index 8fba466..8128f49 100644 --- a/tests/vertex/macros/vert-hdf5.mac +++ b/tests/vertex/macros/vert-hdf5.mac @@ -2,7 +2,7 @@ /RMG/Generator/Confine FromFile -/RMG/Generator/Confinement/FromFile/FileName macros/vert.hdf5 +/RMG/Generator/Confinement/FromFile/FileName macros/vtx-pos.hdf5 /RMG/Generator/Select GPS /gps/particle ion diff --git a/tests/vertex/macros/vert-lh5.mac b/tests/vertex/macros/vert-lh5.mac index ba9bb5a..5acef89 100644 --- a/tests/vertex/macros/vert-lh5.mac +++ b/tests/vertex/macros/vert-lh5.mac @@ -2,7 +2,7 @@ /RMG/Generator/Confine FromFile -/RMG/Generator/Confinement/FromFile/FileName macros/vert.lh5 +/RMG/Generator/Confinement/FromFile/FileName macros/vtx-pos.lh5 /RMG/Generator/Select GPS /gps/particle ion diff --git a/tests/vertex/macros/vert.hdf5 b/tests/vertex/macros/vert.hdf5 deleted file mode 100644 index c53b6aebb93c36adcda2277b70f04ab478898ce2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 19896 zcmeHPO>7!R6n@LvKodK{Y0{slZ8&nIl`C1aRjCwG`j@0Xhui?e+E4*Y0lT=TUMkT; zrCcicNTo*drIL?6R;nWT*kg`XPOj8Lj=81t=Di2mC4jMAnRsU{_RYL`J3I63dvD&% z&id!`PSdjD=M(VHBXJ{}H@au1sVd%17 z6I4VGwXzHI1oZ1$w{KWr_r3Za>#f}9fb&2j07$#yHr-U@vYapb{O~?sCvuNDA6q;@ zqi66^8Mpj6EC6{T3nB`;D16_nw9)^(tojgEKf%T6R4c10zW`}F63Ip zvQsPMHl1dx;?{L}Og)G>or>$T)H9&=bNy0>@fZCueo4b;P9qQ)<(J~DM{=3yo#>VD zk4wd#{QhSFlDbL@=CwVww;B#I(4MeK(vDK4=@i>;a|;XAfj$n)Wc-$Ky32k@+C3@G z+>Yt$)V<9SLxzcZQc*3&xy5Rsg&q*cW&DYJ+jXm*RH|HBdpDI*Q}3o8z47`a zUmw+Ruf)2e;#`^5kD7n)yrQdvdqiN8A3J=D1kJe4fCY?*9~1v2ablFmH-V#3Sa(pb zoqD@jQ7O^uvN)(6nNIlcn2GJr;Ayxbiw&VMpqPTX|wdlJnF z79C|zkA36K|7B19@$KyFDMGOe4_Fqjj9C^jWO+~Cv)Za0OP1mv@|-IT%NJq!DX;tE zECI_yc1Txn5#l<})u}X=Tln6sT4+1E3^gyx_q9WzWE>f^jMq2M%X~Miw77VC9Y7)* zEqiJwio-(v_aKc=$5$;4e8~Q1ymW9%ah`?GS@cFIPsfIU$S1rk(L2G{TWWF!0fT@+ zz#w1{FbJFn1UgnRpAKLD_qX>Wr7y$bM-Qw&Hs)?27pzN>t_e9g1H=4#BZ$)}e+l zF${x%LBJqj5HJXwcLZiI4o!bRpwC1c+WUOG21p!wm>Mq^w9;wg&|gd(8Z;F#{WS;} z1PlTO0fT^A{|xK(&}e@>lp-t6>M!*%9aBFlNo&;l`|m8;e;4JQY_^h;Jyj?;ypGw^ z*U$QC8eE?JE`AsnLd{lR&*Rp`VLWzb9STqT1Q%LF$Mti+zR5lHFB^sfb!4m56}$RH(zwC8?|~1>V=yA&S%YS>VFL* t_ILM``#jC_!#$LvpWk2n{Jc9YU$o~Zco_fZRWyT(7Z1r9Tm%U({SU843-|y4 diff --git a/tests/vertex/macros/vert.lh5 b/tests/vertex/macros/vert.lh5 deleted file mode 100644 index d37a6e42b6c806dbb759c2b0fc72606cc2e41714..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 24576 zcmeI1OKVd>6vxjbedN^_Rje)~E9s((2whYZjrL&{Vv7oPReEg`B19X~RMHlqxO3N~ zOFu%vFCgN^ogcx43qe5_3hwHhnR9A#Y3Rd%B>CTvd*+@w^ZMO$&&<8^GB-0fkhqYD z#A2jJMia|EQIzZVZ@zvu#6O~O(Ve3-ov~@-)lNGsd9*JF=S#UU$r~YWwJXkK8lofv4m`8~jfL}f1y&QTkg*mNX)Lt=fFCX+O|FY&(DYB@jR z!x$xL%3T{-KSf-=q4nZMp-{?89o!EtXJ|EFC~c z`{Kq2-!wk-ed}NS*x?rV?0Sa!?cc7$eRDSLi8pKifBJFO_R}6M_hWu$TFEuEWOSWA z>FM_L{GIw|O>PglIAklU29HYlZJvcDcnX@}`RBRS>ockLPo!Z9Go+7OP^ko)*Dhb* zuH=J~EvRuLTu#T*tX!YeA1)`;RnzvQ7^R0a1qlKm00JNY0w4eaAOHd&00JNY0^LV| z_bK?ho__aHIob?GrZGdm8y$S`Ys|(K2!H?xfB*=9KnDi3;~L%)6wK6I+Bo#|}PP~~ls-hGQ-@v7{t zIpb7`tBDV7GCp(@1V8`;KmY_l00ck)1V8`;KmY_lfas+A|D)LdhueSv2!H?xfB*=9 w00@8p2!H?xfI!<3IO+TUCA|N)-91K!K>!3m00ck)1V8`;KmY_l00cnb56fT^vH$=8 From a9d52cdb0c7973d78e48c89630d8e6e97f1a8647 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Tue, 21 Jan 2025 16:54:36 +0100 Subject: [PATCH 4/9] add generator to read momenta and energies from file --- docs/rmg-commands.md | 34 ++++- include/RMGGeneratorFromFile.hh | 83 ++++++++++ include/RMGMasterGenerator.hh | 1 + src/CMakeLists.txt | 2 + src/RMGGeneratorFromFile.cc | 144 ++++++++++++++++++ src/RMGMasterGenerator.cc | 2 + src/remage-doc-dump.cc | 2 + tests/basics/CMakeLists.txt | 2 +- tests/vertex/CMakeLists.txt | 13 +- tests/vertex/create-test-input.py | 34 ++++- tests/vertex/macros/_init.mac | 5 + tests/vertex/macros/kin-hdf5.mac | 6 + tests/vertex/macros/kin-lh5.mac | 6 + .../macros/{vert-hdf5.mac => pos-hdf5.mac} | 2 +- tests/vertex/macros/pos-kin-lh5.mac | 9 ++ .../macros/{vert-lh5.mac => pos-lh5.mac} | 2 +- tests/vertex/run-vtx.py | 83 ++++++++++ 17 files changed, 416 insertions(+), 14 deletions(-) create mode 100644 include/RMGGeneratorFromFile.hh create mode 100644 src/RMGGeneratorFromFile.cc create mode 100644 tests/vertex/macros/kin-hdf5.mac create mode 100644 tests/vertex/macros/kin-lh5.mac rename tests/vertex/macros/{vert-hdf5.mac => pos-hdf5.mac} (93%) create mode 100644 tests/vertex/macros/pos-kin-lh5.mac rename tests/vertex/macros/{vert-lh5.mac => pos-lh5.mac} (93%) create mode 100755 tests/vertex/run-vtx.py diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index d82980e..4d30d77 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -700,6 +700,7 @@ Commands for controlling generators * `/RMG/Generator/MUSUNCosmicMuons/` – Commands for controlling the MUSUN µ generator * `/RMG/Generator/CosmicMuons/` – Commands for controlling the µ generator +* `/RMG/Generator/FromFile/` – Commands for controlling reading event data from file * `/RMG/Generator/Confinement/` – Commands for controlling primary confinement **Commands:** @@ -723,7 +724,7 @@ Select event generator * **Parameter** – `generator` * **Parameter type** – `s` * **Omittable** – `False` - * **Candidates** – `G4gun GPS BxDecay0 CosmicMuons MUSUNCosmicMuons UserDefined Undefined` + * **Candidates** – `G4gun GPS BxDecay0 FromFile CosmicMuons MUSUNCosmicMuons UserDefined Undefined` ## `/RMG/Generator/MUSUNCosmicMuons/` @@ -940,6 +941,37 @@ Maximum zenith angle of the generated muon position on the sphere * **Default value** – `deg` * **Candidates** – `rad mrad deg radian milliradian degree` +## `/RMG/Generator/FromFile/` + +Commands for controlling reading event data from file + + +**Commands:** + +* `FileName` – Set name of the file containing vertex kinetics for the next run. See the documentation for a specification of the format. +* `NtupleDirectory` – Change the default input directory/group for ntuples. + +### `/RMG/Generator/FromFile/FileName` + +Set name of the file containing vertex kinetics for the next run. See the documentation for a specification of the format. + +* **Parameter** – `filename` + * **Parameter type** – `s` + * **Omittable** – `False` + +### `/RMG/Generator/FromFile/NtupleDirectory` + +Change the default input directory/group for ntuples. + +:::{note} +this option only has an effect for LH5 or HDF5 input files. +::: + +* **Parameter** – `nt_directory` + * **Parameter type** – `s` + * **Omittable** – `False` + * **Default value** – `vtx` + ## `/RMG/Generator/Confinement/` Commands for controlling primary confinement diff --git a/include/RMGGeneratorFromFile.hh b/include/RMGGeneratorFromFile.hh new file mode 100644 index 0000000..ea8b501 --- /dev/null +++ b/include/RMGGeneratorFromFile.hh @@ -0,0 +1,83 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// This program is free software: you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this program. If not, see . + +#ifndef _RMG_GENERATOR_FROM_FILE_HH_ +#define _RMG_GENERATOR_FROM_FILE_HH_ + +#include + +#include "CLHEP/Units/SystemOfUnits.h" +#include "G4GenericMessenger.hh" +#include "G4ParticleGun.hh" + +#include "RMGAnalysisReader.hh" +#include "RMGVGenerator.hh" + +namespace u = CLHEP; + +class G4Event; +class RMGGeneratorFromFile : public RMGVGenerator { + + public: + + RMGGeneratorFromFile(); + ~RMGGeneratorFromFile() = default; + + RMGGeneratorFromFile(RMGGeneratorFromFile const&) = delete; + RMGGeneratorFromFile& operator=(RMGGeneratorFromFile const&) = delete; + RMGGeneratorFromFile(RMGGeneratorFromFile&&) = delete; + RMGGeneratorFromFile& operator=(RMGGeneratorFromFile&&) = delete; + + void GeneratePrimaries(G4Event* event) override; + void SetParticlePosition(G4ThreeVector pos) override { fParticlePosition = pos; } + + void BeginOfRunAction(const G4Run*) override; + void EndOfRunAction(const G4Run*) override; + + void OpenFile(std::string& name); + + private: + + struct RowData { + int fG4Pid = -1; + double fEkin = NAN; + double fPx = NAN; + double fPy = NAN; + double fPz = NAN; + RowData() {}; // NOLINT(modernize-use-equals-default) + + [[nodiscard]] bool IsValid() const { + return fG4Pid != -1 && !std::isnan(fEkin) && !std::isnan(fPx) && !std::isnan(fPy) && + !std::isnan(fPz); + } + }; + + static G4Mutex fMutex; + static RMGAnalysisReader* fReader; + inline static RowData fRowData{}; + + std::unique_ptr fMessenger = nullptr; + void DefineCommands(); + + std::string fNtupleDirectoryName = "vtx"; + + std::unique_ptr fGun = nullptr; + + G4ThreeVector fParticlePosition; +}; + +#endif + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/include/RMGMasterGenerator.hh b/include/RMGMasterGenerator.hh index ffb556b..2ccf461 100644 --- a/include/RMGMasterGenerator.hh +++ b/include/RMGMasterGenerator.hh @@ -39,6 +39,7 @@ class RMGMasterGenerator : public G4VUserPrimaryGeneratorAction { kG4gun, kGPS, kBxDecay0, + kFromFile, kCosmicMuons, kMUSUNCosmicMuons, kUserDefined, diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 607761b..4c941d2 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -9,6 +9,7 @@ set(PROJECT_PUBLIC_HEADERS ${_root}/include/RMGExceptionHandler.hh ${_root}/include/RMGEventAction.hh ${_root}/include/RMGGeneratorCosmicMuons.hh + ${_root}/include/RMGGeneratorFromFile.hh ${_root}/include/RMGGeneratorMUSUNCosmicMuons.hh ${_root}/include/RMGGeneratorG4Gun.hh ${_root}/include/RMGGeneratorGPS.hh @@ -52,6 +53,7 @@ set(PROJECT_SOURCES ${_root}/src/RMGHardwareMessenger.cc ${_root}/src/RMGEventAction.cc ${_root}/src/RMGGeneratorCosmicMuons.cc + ${_root}/src/RMGGeneratorFromFile.cc ${_root}/src/RMGGeneratorMUSUNCosmicMuons.cc ${_root}/src/RMGGeneratorUtil.cc ${_root}/src/RMGGermaniumDetector.cc diff --git a/src/RMGGeneratorFromFile.cc b/src/RMGGeneratorFromFile.cc new file mode 100644 index 0000000..3ab12e6 --- /dev/null +++ b/src/RMGGeneratorFromFile.cc @@ -0,0 +1,144 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// This program is free software: you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this program. If not, see . + +#include "RMGGeneratorFromFile.hh" + +#include + +#include "G4ParticleGun.hh" +#include "G4ParticleMomentum.hh" +#include "G4ParticleTypes.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; + +G4Mutex RMGGeneratorFromFile::fMutex = G4MUTEX_INITIALIZER; + +RMGAnalysisReader* RMGGeneratorFromFile::fReader = new RMGAnalysisReader(); + +RMGGeneratorFromFile::RMGGeneratorFromFile() : RMGVGenerator("FromFile") { + this->DefineCommands(); + fGun = std::make_unique(); +} + +void RMGGeneratorFromFile::OpenFile(std::string& name) { + + // reader initialization should only happen on the master thread (otherwise it will fail). + if (!G4Threading::IsMasterThread()) return; + G4AutoLock lock(&fMutex); + + if (!fReader->OpenFile(name, fNtupleDirectoryName, "kin")) return; + + auto nt = fReader->GetNtupleId(); + if (nt >= 0) { + // bind the static variables once here, and only use them later. + fReader->GetReader()->SetNtupleIColumn(nt, "g4_pid", fRowData.fG4Pid); + fReader->GetReader()->SetNtupleDColumn(nt, "ekin_in_MeV", fRowData.fEkin); + fReader->GetReader()->SetNtupleDColumn(nt, "px_in_MeV", fRowData.fPx); + fReader->GetReader()->SetNtupleDColumn(nt, "py_in_MeV", fRowData.fPy); + fReader->GetReader()->SetNtupleDColumn(nt, "pz_in_MeV", fRowData.fPz); + } +} + +void RMGGeneratorFromFile::BeginOfRunAction(const G4Run*) { + + if (!G4Threading::IsMasterThread()) return; + G4AutoLock lock(&fMutex); + + if (!fReader->GetReader()) { + RMGLog::Out(RMGLog::fatal, "vertex file '", fReader->GetFileName(), + "' not found or in wrong format"); + } +} + +void RMGGeneratorFromFile::EndOfRunAction(const G4Run*) { + + if (!G4Threading::IsMasterThread()) return; + G4AutoLock lock(&fMutex); + + fReader->CloseFile(); +} + + +void RMGGeneratorFromFile::GeneratePrimaries(G4Event* event) { + + G4AutoLock lock(&fMutex); + + if (!fReader->GetReader() || fReader->GetNtupleId() < 0) { + RMGLog::Out(RMGLog::error, "Ntuple named 'pos' could not be found in input file!"); + return; + } + + fRowData = RowData(); // initialize sentinel values. + + auto nt = fReader->GetNtupleId(); + if (!fReader->GetReader()->GetNtupleRow(nt)) { + RMGLog::Out(RMGLog::error, "No more vertices available in input file!"); + return; + } + + // check for NaN sentinel values - i.e. non-existing columns (there is no error message). + if (!fRowData.IsValid()) { + RMGLog::Out(RMGLog::error, "At least one of the columns does not exist"); + return; + } + + auto particle = G4ParticleTable::GetParticleTable()->FindParticle(fRowData.fG4Pid); + if (!particle) { + RMGLog::Out(RMGLog::error, "invalid particle PDG id ", fRowData.fG4Pid); + return; + } + + RMGLog::OutFormat(RMGLog::debug, + "particle {:d} (px,py,pz) = ({:.4g}, {:.4g}, {:.4g}) MeV ; Ekin = {:.4g} MeV", + fRowData.fG4Pid, fRowData.fPx, fRowData.fPy, fRowData.fPz, fRowData.fEkin); + + G4ThreeVector momentum{fRowData.fPx, fRowData.fPy, fRowData.fPz}; + + fGun->SetParticleDefinition(particle); + fGun->SetParticlePosition(fParticlePosition); + fGun->SetParticleMomentumDirection(momentum * u::MeV); + fGun->SetParticleEnergy(fRowData.fEkin * u::MeV); + + fGun->GeneratePrimaryVertex(event); +} + +void RMGGeneratorFromFile::DefineCommands() { + + fMessenger = std::make_unique(this, "/RMG/Generator/FromFile/", + "Commands for controlling reading event data from file"); + + fMessenger->DeclareMethod("FileName", &RMGGeneratorFromFile::OpenFile) + .SetGuidance("Set name of the file containing vertex kinetics for the next run. See the " + "documentation for a specification of the format.") + .SetParameterName("filename", false) + .SetStates(G4State_PreInit, G4State_Idle); + + fMessenger->DeclareProperty("NtupleDirectory", fNtupleDirectoryName) + .SetGuidance("Change the default input directory/group for ntuples.") + .SetGuidance("note: this option only has an effect for LH5 or HDF5 input files.") + .SetParameterName("nt_directory", false) + .SetDefaultValue(fNtupleDirectoryName) + .SetStates(G4State_PreInit, G4State_Idle); +} + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/src/RMGMasterGenerator.cc b/src/RMGMasterGenerator.cc index 0a79c24..5b71735 100644 --- a/src/RMGMasterGenerator.cc +++ b/src/RMGMasterGenerator.cc @@ -23,6 +23,7 @@ #if RMG_HAS_BXDECAY0 #include "RMGGeneratorDecay0.hh" #endif +#include "RMGGeneratorFromFile.hh" #include "RMGGeneratorG4Gun.hh" #include "RMGGeneratorGPS.hh" #include "RMGGeneratorMUSUNCosmicMuons.hh" @@ -113,6 +114,7 @@ void RMGMasterGenerator::SetGenerator(RMGMasterGenerator::Generator gen) { case Generator::kMUSUNCosmicMuons: fGeneratorObj = std::make_unique(); break; + case Generator::kFromFile: fGeneratorObj = std::make_unique(); break; case Generator::kUndefined: case Generator::kUserDefined: break; default: diff --git a/src/remage-doc-dump.cc b/src/remage-doc-dump.cc index 32e3c70..b1db846 100644 --- a/src/remage-doc-dump.cc +++ b/src/remage-doc-dump.cc @@ -21,6 +21,7 @@ #include "G4UImanager.hh" #include "RMGGeneratorCosmicMuons.hh" +#include "RMGGeneratorFromFile.hh" #include "RMGGeneratorMUSUNCosmicMuons.hh" #include "RMGGermaniumOutputScheme.hh" #include "RMGIsotopeFilterOutputScheme.hh" @@ -50,6 +51,7 @@ void init_extra() { // generators new RMGGeneratorMUSUNCosmicMuons(); new RMGGeneratorCosmicMuons(); + new RMGGeneratorFromFile(); // confinments new RMGVertexConfinement(); new RMGVertexFromFile(); diff --git a/tests/basics/CMakeLists.txt b/tests/basics/CMakeLists.txt index b701cff..2293695 100644 --- a/tests/basics/CMakeLists.txt +++ b/tests/basics/CMakeLists.txt @@ -30,7 +30,7 @@ endif() if(BxDecay0_THREADSAFE) foreach(_mac ${_macros_extra}) add_test(NAME basics-mt/${_mac} COMMAND ${REMAGE_PYEXE} -g gdml/main.gdml -t 2 macros/${_mac}) - set_tests_properties(basics-mt/${_mac} PROPERTIES LABELS "mt extra") + set_tests_properties(basics-mt/${_mac} PROPERTIES LABELS "mt;extra") endforeach() # Geant4 11.0 has navigation bugs leading to flaky CI test runs. diff --git a/tests/vertex/CMakeLists.txt b/tests/vertex/CMakeLists.txt index b89e5c0..dfc81b3 100644 --- a/tests/vertex/CMakeLists.txt +++ b/tests/vertex/CMakeLists.txt @@ -17,13 +17,14 @@ endforeach() add_test(NAME vertex/gen-input COMMAND ${PYTHONPATH} create-test-input.py ${_remage_from_lh5}) set_tests_properties(vertex/gen-input PROPERTIES LABELS extra FIXTURES_SETUP vertex-input-fixture) -set(_macros vert-hdf5.mac vert-lh5.mac) +set(_macros pos-hdf5 pos-lh5 kin-lh5 kin-hdf5 pos-kin-lh5) foreach(_mac ${_macros}) - add_test(NAME vertex/${_mac} COMMAND ${REMAGE_PYEXE} -g gdml/geometry.gdml -- macros/${_mac}) - set_tests_properties(vertex/${_mac} PROPERTIES LABELS extra FIXTURES_SETUP vertex-input-fixture) - add_test(NAME vertex-mt/${_mac} COMMAND ${REMAGE_PYEXE} -g gdml/geometry.gdml -t 2 - macros/${_mac}) - set_tests_properties(vertex-mt/${_mac} PROPERTIES LABELS "mt extra" FIXTURES_SETUP + add_test(NAME vertex/${_mac} COMMAND ${PYTHONPATH} run-vtx.py ${REMAGE_PYEXE} ${_mac} 1) + set_tests_properties(vertex/${_mac} PROPERTIES LABELS extra FIXTURES_REQUIRED vertex-input-fixture) + + # note: only 4 threads available on GitHub actions. + add_test(NAME vertex-mt/${_mac} COMMAND ${PYTHONPATH} run-vtx.py ${REMAGE_PYEXE} ${_mac} 4) + set_tests_properties(vertex-mt/${_mac} PROPERTIES LABELS "mt;extra" FIXTURES_REQUIRED vertex-input-fixture) endforeach() diff --git a/tests/vertex/create-test-input.py b/tests/vertex/create-test-input.py index a2bf22c..eafed34 100755 --- a/tests/vertex/create-test-input.py +++ b/tests/vertex/create-test-input.py @@ -20,17 +20,20 @@ def _convert_lh5_to_hdf5(lh5_fn): new_path = orig_path.with_suffix(".hdf5") shutil.copy(orig_path, new_path) ret = subprocess.run( - [rmg_from_lh5, "--ntuple-group", "vtx", "--", new_path], check=False + [rmg_from_lh5, "--ntuple-group", "vtx", "--", str(new_path)], check=False ) assert ret.returncode == 0 +INPUT_FILE_ROWS = int(1e6) +rng = np.random.default_rng(123456) + # ================================================== # position input -xloc = Array(np.array([0, 1, 2] * 100, dtype=np.float64), attrs={"units": "m"}) -yloc = Array(np.array([3, 4, 5] * 100, dtype=np.float64), attrs={"units": "m"}) -zloc = Array(np.array([6, 7, 8] * 100, dtype=np.float64), attrs={"units": "m"}) +xloc = Array(np.linspace(0, 1000, INPUT_FILE_ROWS), attrs={"units": "m"}) +yloc = Array(rng.uniform(-1, 1, INPUT_FILE_ROWS), attrs={"units": "m"}) +zloc = Array(rng.uniform(-1, 1, INPUT_FILE_ROWS), attrs={"units": "m"}) vtx_pos_file = "macros/vtx-pos.lh5" lh5.write( @@ -40,3 +43,26 @@ def _convert_lh5_to_hdf5(lh5_fn): wo_mode="of", ) _convert_lh5_to_hdf5(vtx_pos_file) + + +# ================================================== +# kinetcis input. + +p_theta = np.acos(rng.uniform(-1, 1, INPUT_FILE_ROWS)) +p_phi = rng.uniform(0, 2 * np.pi, INPUT_FILE_ROWS) + +px = Array(np.cos(p_phi) * np.sin(p_theta)) +py = Array(np.sin(p_phi) * np.sin(p_theta)) +pz = Array(np.cos(p_theta)) + +pdg = Array(np.ones(INPUT_FILE_ROWS, dtype=np.int64) * 11) +ekin = Array(np.linspace(1, 10, INPUT_FILE_ROWS), attrs={"units": "MeV"}) + +vtx_kin_file = "macros/vtx-kin.lh5" +lh5.write( + Table({"g4_pid": pdg, "ekin": ekin, "px": px, "py": py, "pz": pz}), + "vtx/kin", + lh5_file=vtx_kin_file, + wo_mode="of", +) +_convert_lh5_to_hdf5(vtx_kin_file) diff --git a/tests/vertex/macros/_init.mac b/tests/vertex/macros/_init.mac index 301cf8b..a1e497f 100644 --- a/tests/vertex/macros/_init.mac +++ b/tests/vertex/macros/_init.mac @@ -4,6 +4,11 @@ /RMG/Geometry/GDMLDisableOverlapCheck true +/RMG/Geometry/RegisterDetector Germanium det1 11 +/RMG/Geometry/RegisterDetector Germanium det2 12 + /run/initialize /RMG/Manager/Logging/LogLevel detail + +/RMG/Output/Vertex/StorePrimaryParticleInformation true diff --git a/tests/vertex/macros/kin-hdf5.mac b/tests/vertex/macros/kin-hdf5.mac new file mode 100644 index 0000000..05e3927 --- /dev/null +++ b/tests/vertex/macros/kin-hdf5.mac @@ -0,0 +1,6 @@ +/control/execute macros/_init.mac + +/RMG/Generator/Select FromFile +/RMG/Generator/FromFile/FileName macros/vtx-kin.hdf5 + +/run/beamOn 100000 diff --git a/tests/vertex/macros/kin-lh5.mac b/tests/vertex/macros/kin-lh5.mac new file mode 100644 index 0000000..5a0ebe5 --- /dev/null +++ b/tests/vertex/macros/kin-lh5.mac @@ -0,0 +1,6 @@ +/control/execute macros/_init.mac + +/RMG/Generator/Select FromFile +/RMG/Generator/FromFile/FileName macros/vtx-kin.lh5 + +/run/beamOn 100000 diff --git a/tests/vertex/macros/vert-hdf5.mac b/tests/vertex/macros/pos-hdf5.mac similarity index 93% rename from tests/vertex/macros/vert-hdf5.mac rename to tests/vertex/macros/pos-hdf5.mac index 8128f49..a60dfb0 100644 --- a/tests/vertex/macros/vert-hdf5.mac +++ b/tests/vertex/macros/pos-hdf5.mac @@ -10,4 +10,4 @@ /gps/energy 0 keV /gps/ang/type iso -/run/beamOn 10 +/run/beamOn 100000 diff --git a/tests/vertex/macros/pos-kin-lh5.mac b/tests/vertex/macros/pos-kin-lh5.mac new file mode 100644 index 0000000..d736a53 --- /dev/null +++ b/tests/vertex/macros/pos-kin-lh5.mac @@ -0,0 +1,9 @@ +/control/execute macros/_init.mac + +/RMG/Generator/Confine FromFile +/RMG/Generator/Confinement/FromFile/FileName macros/vtx-pos.lh5 + +/RMG/Generator/Select FromFile +/RMG/Generator/FromFile/FileName macros/vtx-kin.lh5 + +/run/beamOn 100000 diff --git a/tests/vertex/macros/vert-lh5.mac b/tests/vertex/macros/pos-lh5.mac similarity index 93% rename from tests/vertex/macros/vert-lh5.mac rename to tests/vertex/macros/pos-lh5.mac index 5acef89..54448fa 100644 --- a/tests/vertex/macros/vert-lh5.mac +++ b/tests/vertex/macros/pos-lh5.mac @@ -10,4 +10,4 @@ /gps/energy 0 keV /gps/ang/type iso -/run/beamOn 10 +/run/beamOn 100000 diff --git a/tests/vertex/run-vtx.py b/tests/vertex/run-vtx.py new file mode 100755 index 0000000..b8874f9 --- /dev/null +++ b/tests/vertex/run-vtx.py @@ -0,0 +1,83 @@ +#!/bin/env python3 + +from __future__ import annotations + +import subprocess +import sys +from pathlib import Path + +import numpy as np +from lgdo import lh5 + +rmg = sys.argv[1] +macro = sys.argv[2] +nthread = int(sys.argv[3]) + +output_stem = macro + ("-mt" if nthread > 1 else "") +output_lh5 = f"{output_stem}.lh5" + +# clean up files from last run(s). +for old_f in Path().glob(f"{output_stem}*.lh5"): + Path(old_f).unlink() + +# run remage, produce lh5 output. +cmd = [ + rmg, + "-g", + "gdml/geometry.gdml", + "-o", + output_lh5, + "-t", + str(nthread), + "-l", + "summary", + "--", + f"macros/{macro}.mac", +] +proc = subprocess.run(cmd, check=False) + +if proc.returncode not in [0, 2]: + print("remage subprocess failed") + sys.exit(proc.returncode) + +# list of output files. +files = output_lh5 +if nthread > 1: + files = [f"{output_stem}_t{i}.lh5" for i in range(nthread)] + +# validate output against input files. +if "pos" in macro: + input_pos = lh5.read("vtx/pos", lh5_file="macros/vtx-pos.lh5").view_as("pd") + output_pos = lh5.read("stp/vertices", lh5_file=files).view_as("pd") + + output_pos = output_pos.sort_values("xloc") # sort by linear column. + output_pos = output_pos[["xloc", "yloc", "zloc"]] + uniq, cnt = np.unique(output_pos["xloc"], return_counts=True) + if not np.all(cnt <= 1): + raise ValueError(uniq[cnt > 1]) + + input_pos = input_pos.iloc[0 : len(output_pos)][["xloc", "yloc", "zloc"]] + + assert np.all(np.isclose(output_pos.to_numpy(), input_pos.to_numpy())) + +if "kin" in macro: + input_kin = lh5.read("vtx/kin", lh5_file="macros/vtx-kin.lh5").view_as("pd") + output_kin = lh5.read("stp/particles", lh5_file=files).view_as("pd") + + output_kin = output_kin.sort_values("ekin") # sort by linear column. + output_kin = output_kin[["px", "py", "pz", "ekin", "particle"]] + uniq, cnt = np.unique(output_kin["px"], return_counts=True) + if not np.all(cnt <= 1): + raise ValueError(uniq[cnt > 1]) + + output_p_scale = np.sqrt( + output_kin["px"] ** 2 + output_kin["py"] ** 2 + output_kin["pz"] ** 2 + ) + output_kin["px"] /= output_p_scale + output_kin["py"] /= output_p_scale + output_kin["pz"] /= output_p_scale + + input_kin = input_kin.iloc[0 : len(output_kin)] + input_kin = input_kin[["px", "py", "pz", "ekin", "g4_pid"]] + + assert np.all(np.isclose(output_kin.to_numpy(), input_kin.to_numpy())) From 2aaab8f89242e207b64ce8ddfc6731e1773c10e5 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Thu, 23 Jan 2025 01:47:53 +0100 Subject: [PATCH 5/9] refactor RMGAnalysisReader --- include/RMGAnalysisReader.hh | 103 +++++++++++++++++++++++++------- include/RMGGeneratorFromFile.hh | 1 - include/RMGVertexFromFile.hh | 1 - src/RMGAnalysisReader.cc | 46 ++++++++++---- src/RMGGeneratorFromFile.cc | 64 +++++++++----------- src/RMGVertexFromFile.cc | 29 ++++----- tests/vertex/CMakeLists.txt | 3 +- 7 files changed, 157 insertions(+), 90 deletions(-) diff --git a/include/RMGAnalysisReader.hh b/include/RMGAnalysisReader.hh index 3e55c9a..964e8cb 100644 --- a/include/RMGAnalysisReader.hh +++ b/include/RMGAnalysisReader.hh @@ -18,29 +18,86 @@ #include -#include "G4ThreeVector.hh" +#include "G4AutoLock.hh" #include "G4VAnalysisReader.hh" -#include "RMGVVertexGenerator.hh" +#include "RMGConfig.hh" /** - * @brief wrapper around \ref(G4VAnalasisReader) instances with special handling for LH5 files. + * @brief wrapper around @ref G4VAnalysisReader instances with special handling for LH5 files. * * @details notes for threadsafe use: * - opening/closing can only be performed on the master thread. - * - in a multithreaded application, all function calls need to by guarded by a mutex. Worker - * threads can use the reader - * - instance after opening, but only one worker can use the reader at a time. - * - the instance obtained by \ref(GetReader) should only be bound once to variables that are of same - * storage duration as the \ref(RMGAnalysisReader) instance. Example: If you use a static instance - * in a class, only bind to static class fields (of the same class) to read the values into. - * - it is EXTREMELY important to always use overloads with the ntuple id parameter, and to always - * pass the result of \ref(GetNtupleID) to them. This affects ntuple reading and column registration. + * - in a multithreaded application, all function calls are guarded by a mutex. Worker + * threads can use the reader instance after opening, but only one worker can use the reader at a + * time. + * - the reader access handles are generally thread-safe, if no other thread uses any + * @ref G4VAnalysisReader directly. + * - the reader should only be bound once to variables that are of static + * storage duration. Example: only bind to static class fields to read the values into. Do only + * unloack the reader access handle after reading/checking the read data. */ class RMGAnalysisReader final { public: + /** + * @brief thread-safe access handle to the underlying reader. */ + class Access final { + friend class RMGAnalysisReader; + + public: + + ~Access() { unlock(); } + + Access(Access const&) = delete; + Access& operator=(Access&) = delete; + Access& operator=(Access const&) = delete; + Access& operator=(Access&&) = delete; + + /** + * @brief unlock this access handle before it exits the scope. */ + void unlock() { + fReader = nullptr; + fNtupleId = -1; + if (fLock) { fLock.unlock(); } + } + + /** + * @brief wraps @ref G4VAnalysisReader::GetNtupleRow. */ + [[nodiscard]] auto GetNtupleRow() { return fReader->GetNtupleRow(fNtupleId); } + /** + * @brief wraps @ref G4VAnalysisReader::SetNtupleDColumn. */ + auto SetNtupleDColumn(const std::string& name, G4double& value) { + return fReader->SetNtupleDColumn(fNtupleId, name, value); + } + /** + * @brief wraps @ref G4VAnalysisReader::SetNtupleFColumn. */ + auto SetNtupleFColumn(const std::string& name, G4float& value) { + return fReader->SetNtupleFColumn(fNtupleId, name, value); + } + /** + * @brief wraps @ref G4VAnalysisReader::SetNtupleIColumn. */ + auto SetNtupleIColumn(const std::string& name, G4int& value) { + return fReader->SetNtupleIColumn(fNtupleId, name, value); + } + + /** + * @brief check whether this access handle is still valid. */ + operator bool() const { return fReader != nullptr && fNtupleId >= 0 && fLock; } + + private: + + // only allow creation or moving in parent. + inline Access(G4AutoLock lock, G4VAnalysisReader* reader, int nt) + : fLock(std::move(lock)), fReader(reader), fNtupleId(nt) {}; + Access(Access&&) = default; + + G4VAnalysisReader* fReader = nullptr; + int fNtupleId = -1; + G4AutoLock fLock; + }; + RMGAnalysisReader() = default; ~RMGAnalysisReader() = default; @@ -52,7 +109,8 @@ class RMGAnalysisReader final { /** * @brief open an input file for reading of one specific ntuple. * - * @details This function can only be used on the master thread. + * @details This function can only be used on the master thread. This operation acquires a + * global lock to avoid problems. * * @param file_name the input file name. the file format is determined from the file extension. * @param ntuple_dir_name the first part of the input table name. For the table addressed by @@ -60,26 +118,31 @@ class RMGAnalysisReader final { * @param ntuple_name the first part of the input table name. For the table addressed by * "dir/table" this is "table". */ - bool OpenFile(std::string& file_name, std::string ntuple_dir_name, std::string ntuple_name); + [[nodiscard]] Access OpenFile(std::string& file_name, std::string ntuple_dir_name, + std::string ntuple_name); + /** * @brief if any file is open for reading, close the reader. * - * @details this invalidates all reader instances obtained via \ref(GetReader). This function - * can only be used on the master thread. */ + * @details this invalidates all readers obtained via @ref RMGAnalysisReader::GetLockedReader. + * This function can only be used on the master thread. This operation acquires a global lock to + * avoid problems. */ void CloseFile(); /** - * @brief get a pointer to the current underlying G4VAnalysisReader. */ - [[nodiscard]] inline auto GetReader() { return fReader; } - /** - * @brief return the ntuple id of the current tuple. If the value is negative, no ntuple is opened. */ - [[nodiscard]] inline auto GetNtupleId() const { return fNtupleId; } + * @brief get an access handle to the current underlying G4VAnalysisReader. + * + * @details this acquires a global lock to avoid problems. The lock is held until the access handle is discarded. */ + [[nodiscard]] Access GetLockedReader() const; + /** * @brief get the file name of the current open file, or an empty string. */ [[nodiscard]] inline auto& GetFileName() const { return fFileName; } private: + static G4Mutex fMutex; + G4VAnalysisReader* fReader = nullptr; int fNtupleId = -1; diff --git a/include/RMGGeneratorFromFile.hh b/include/RMGGeneratorFromFile.hh index ea8b501..5de1e69 100644 --- a/include/RMGGeneratorFromFile.hh +++ b/include/RMGGeneratorFromFile.hh @@ -64,7 +64,6 @@ class RMGGeneratorFromFile : public RMGVGenerator { } }; - static G4Mutex fMutex; static RMGAnalysisReader* fReader; inline static RowData fRowData{}; diff --git a/include/RMGVertexFromFile.hh b/include/RMGVertexFromFile.hh index afa8889..0cffa63 100644 --- a/include/RMGVertexFromFile.hh +++ b/include/RMGVertexFromFile.hh @@ -46,7 +46,6 @@ class RMGVertexFromFile : public RMGVVertexGenerator { private: - static G4Mutex fMutex; static RMGAnalysisReader* fReader; inline static double fXpos = NAN, fYpos = NAN, fZpos = NAN; diff --git a/src/RMGAnalysisReader.cc b/src/RMGAnalysisReader.cc index f9a45ad..369429e 100644 --- a/src/RMGAnalysisReader.cc +++ b/src/RMGAnalysisReader.cc @@ -35,19 +35,31 @@ namespace fs = std::filesystem; #include "RMGLog.hh" -bool RMGAnalysisReader::OpenFile(std::string& file_name, std::string ntuple_dir_name, - std::string ntuple_name) { +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); // reader initialization should only happen on the master thread (otherwise it will fail). - if (!G4Threading::IsMasterThread()) return false; + if (!G4Threading::IsMasterThread()) { + RMGLog::OutDev(RMGLog::fatal, "can only be used on the master thread"); + std::abort(); + } - if (fReader) return false; + G4AutoLock lock(&fMutex); + + if (fReader) return std::move(invalid_access); fFileIsTemp = false; if (!fs::exists(fs::path(file_name))) { RMGLog::Out(RMGLog::error, "input file ", file_name, " does not exist"); - return false; + return std::move(invalid_access); } // NOTE: GetExtension() returns a default extension if there is no file extension @@ -73,13 +85,13 @@ bool RMGAnalysisReader::OpenFile(std::string& file_name, std::string ntuple_dir_ if (!fs::copy_file(path, fs::path(new_fn)), ec) { RMGLog::Out(RMGLog::error, "copy of input file ", file_name, " failed. Does it exist? (", ec.message(), ")"); - return false; + return std::move(invalid_access); } if (!RMGConvertLH5::ConvertFromLH5(new_fn, ntuple_dir_name, false)) { RMGLog::Out(RMGLog::error, "Conversion of input file ", new_fn, " to LH5 failed. Data is potentially corrupted."); - return false; + return std::move(invalid_access); } file_name = new_fn; @@ -94,7 +106,7 @@ bool RMGAnalysisReader::OpenFile(std::string& file_name, std::string ntuple_dir_ } else if (ext == "xml") { fReader = G4XmlAnalysisReader::Instance(); } else { - RMGLog::OutFormat(RMGLog::fatal, "File Extension '.{}' not recognized!"); + RMGLog::OutFormat(RMGLog::fatal, "File Extension '.{}' not recognized!", ext); } if (RMGLog::GetLogLevel() <= RMGLog::debug) fReader->SetVerboseLevel(10); @@ -104,13 +116,18 @@ bool RMGAnalysisReader::OpenFile(std::string& file_name, std::string ntuple_dir_ fNtupleId = fReader->GetNtuple(ntuple_name, file_name, ntuple_dir_name); if (fNtupleId < 0) { RMGLog::Out(RMGLog::error, "Ntuple named '", ntuple_name, "' could not be found in input file!"); - return false; + return std::move(invalid_access); } - return true; + return {std::move(lock), fReader, fNtupleId}; } void RMGAnalysisReader::CloseFile() { - if (!G4Threading::IsMasterThread()) return; + if (!G4Threading::IsMasterThread()) { + RMGLog::OutDev(RMGLog::fatal, "can only be used on the master thread"); + return; + } + + G4AutoLock lock(&fMutex); if (!fReader) return; @@ -125,3 +142,10 @@ void RMGAnalysisReader::CloseFile() { fFileIsTemp = false; fNtupleId = -1; } + +RMGAnalysisReader::Access RMGAnalysisReader::GetLockedReader() const { + + G4AutoLock lock(&fMutex); + + return {std::move(lock), fReader, fNtupleId}; +} diff --git a/src/RMGGeneratorFromFile.cc b/src/RMGGeneratorFromFile.cc index 3ab12e6..8f68ba1 100644 --- a/src/RMGGeneratorFromFile.cc +++ b/src/RMGGeneratorFromFile.cc @@ -18,21 +18,13 @@ #include #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; -G4Mutex RMGGeneratorFromFile::fMutex = G4MUTEX_INITIALIZER; - RMGAnalysisReader* RMGGeneratorFromFile::fReader = new RMGAnalysisReader(); RMGGeneratorFromFile::RMGGeneratorFromFile() : RMGVGenerator("FromFile") { @@ -44,27 +36,23 @@ void RMGGeneratorFromFile::OpenFile(std::string& name) { // reader initialization should only happen on the master thread (otherwise it will fail). if (!G4Threading::IsMasterThread()) return; - G4AutoLock lock(&fMutex); - - if (!fReader->OpenFile(name, fNtupleDirectoryName, "kin")) return; - - auto nt = fReader->GetNtupleId(); - if (nt >= 0) { - // bind the static variables once here, and only use them later. - fReader->GetReader()->SetNtupleIColumn(nt, "g4_pid", fRowData.fG4Pid); - fReader->GetReader()->SetNtupleDColumn(nt, "ekin_in_MeV", fRowData.fEkin); - fReader->GetReader()->SetNtupleDColumn(nt, "px_in_MeV", fRowData.fPx); - fReader->GetReader()->SetNtupleDColumn(nt, "py_in_MeV", fRowData.fPy); - fReader->GetReader()->SetNtupleDColumn(nt, "pz_in_MeV", fRowData.fPz); - } + + auto reader = fReader->OpenFile(name, fNtupleDirectoryName, "kin"); + if (!reader) return; + + // bind the static variables once here, and only use them later. + reader.SetNtupleIColumn("g4_pid", fRowData.fG4Pid); + reader.SetNtupleDColumn("ekin_in_MeV", fRowData.fEkin); + reader.SetNtupleDColumn("px", fRowData.fPx); + reader.SetNtupleDColumn("py", fRowData.fPy); + reader.SetNtupleDColumn("pz", fRowData.fPz); } void RMGGeneratorFromFile::BeginOfRunAction(const G4Run*) { if (!G4Threading::IsMasterThread()) return; - G4AutoLock lock(&fMutex); - if (!fReader->GetReader()) { + if (!fReader->GetLockedReader()) { RMGLog::Out(RMGLog::fatal, "vertex file '", fReader->GetFileName(), "' not found or in wrong format"); } @@ -73,7 +61,6 @@ void RMGGeneratorFromFile::BeginOfRunAction(const G4Run*) { void RMGGeneratorFromFile::EndOfRunAction(const G4Run*) { if (!G4Threading::IsMasterThread()) return; - G4AutoLock lock(&fMutex); fReader->CloseFile(); } @@ -81,43 +68,46 @@ void RMGGeneratorFromFile::EndOfRunAction(const G4Run*) { void RMGGeneratorFromFile::GeneratePrimaries(G4Event* event) { - G4AutoLock lock(&fMutex); + auto locked_reader = fReader->GetLockedReader(); - if (!fReader->GetReader() || fReader->GetNtupleId() < 0) { + if (!locked_reader) { RMGLog::Out(RMGLog::error, "Ntuple named 'pos' could not be found in input file!"); return; } fRowData = RowData(); // initialize sentinel values. - auto nt = fReader->GetNtupleId(); - if (!fReader->GetReader()->GetNtupleRow(nt)) { + if (!locked_reader.GetNtupleRow()) { RMGLog::Out(RMGLog::error, "No more vertices available in input file!"); return; } + // make copy of data and exit critical section. + auto row_data = fRowData; + locked_reader.unlock(); + // check for NaN sentinel values - i.e. non-existing columns (there is no error message). - if (!fRowData.IsValid()) { + if (!row_data.IsValid()) { RMGLog::Out(RMGLog::error, "At least one of the columns does not exist"); return; } - auto particle = G4ParticleTable::GetParticleTable()->FindParticle(fRowData.fG4Pid); + auto particle = G4ParticleTable::GetParticleTable()->FindParticle(row_data.fG4Pid); if (!particle) { - RMGLog::Out(RMGLog::error, "invalid particle PDG id ", fRowData.fG4Pid); + RMGLog::Out(RMGLog::error, "invalid particle PDG id ", row_data.fG4Pid); return; } RMGLog::OutFormat(RMGLog::debug, - "particle {:d} (px,py,pz) = ({:.4g}, {:.4g}, {:.4g}) MeV ; Ekin = {:.4g} MeV", - fRowData.fG4Pid, fRowData.fPx, fRowData.fPy, fRowData.fPz, fRowData.fEkin); + "particle {:d} (px,py,pz) = ({:.4g}, {:.4g}, {:.4g}); Ekin = {:.4g} MeV", row_data.fG4Pid, + row_data.fPx, row_data.fPy, row_data.fPz, row_data.fEkin); - G4ThreeVector momentum{fRowData.fPx, fRowData.fPy, fRowData.fPz}; + G4ThreeVector momentum{row_data.fPx, row_data.fPy, row_data.fPz}; fGun->SetParticleDefinition(particle); fGun->SetParticlePosition(fParticlePosition); - fGun->SetParticleMomentumDirection(momentum * u::MeV); - fGun->SetParticleEnergy(fRowData.fEkin * u::MeV); + fGun->SetParticleMomentumDirection(momentum); + fGun->SetParticleEnergy(row_data.fEkin * u::MeV); fGun->GeneratePrimaryVertex(event); } diff --git a/src/RMGVertexFromFile.cc b/src/RMGVertexFromFile.cc index 0f266ed..4054e91 100644 --- a/src/RMGVertexFromFile.cc +++ b/src/RMGVertexFromFile.cc @@ -16,13 +16,10 @@ #include "RMGVertexFromFile.hh" #include "CLHEP/Units/SystemOfUnits.h" -#include "G4AutoLock.hh" #include "G4Threading.hh" #include "RMGLog.hh" -G4Mutex RMGVertexFromFile::fMutex = G4MUTEX_INITIALIZER; - RMGAnalysisReader* RMGVertexFromFile::fReader = new RMGAnalysisReader(); RMGVertexFromFile::RMGVertexFromFile() : RMGVVertexGenerator("FromFile") { this->DefineCommands(); } @@ -31,28 +28,24 @@ void RMGVertexFromFile::OpenFile(std::string& name) { // reader initialization should only happen on the master thread (otherwise it will fail). if (!G4Threading::IsMasterThread()) return; - G4AutoLock lock(&fMutex); - if (!fReader->OpenFile(name, fNtupleDirectoryName, "pos")) return; + auto reader = fReader->OpenFile(name, fNtupleDirectoryName, "pos"); + if (!reader) return; - auto nt = fReader->GetNtupleId(); - if (nt >= 0) { - // bind the static variables once here, and only use them later. - fReader->GetReader()->SetNtupleDColumn(nt, "xloc_in_m", fXpos); - fReader->GetReader()->SetNtupleDColumn(nt, "yloc_in_m", fYpos); - fReader->GetReader()->SetNtupleDColumn(nt, "zloc_in_m", fZpos); - } + // bind the static variables once here, and only use them later. + reader.SetNtupleDColumn("xloc_in_m", fXpos); + reader.SetNtupleDColumn("yloc_in_m", fYpos); + reader.SetNtupleDColumn("zloc_in_m", fZpos); } bool RMGVertexFromFile::GenerateVertex(G4ThreeVector& vertex) { - G4AutoLock lock(&fMutex); + auto reader = fReader->GetLockedReader(); - if (fReader->GetReader() && fReader->GetNtupleId() >= 0) { + if (reader) { fXpos = fYpos = fZpos = NAN; // initialize sentinel values. - auto nt = fReader->GetNtupleId(); - if (fReader->GetReader()->GetNtupleRow(nt)) { + if (reader.GetNtupleRow()) { // check for NaN sentinel values - i.e. non-existing columns (there is no error message). if (std::isnan(fXpos) || std::isnan(fYpos) || std::isnan(fZpos)) { RMGLog::Out(RMGLog::error, "At least one of the columns does not exist"); @@ -76,9 +69,8 @@ bool RMGVertexFromFile::GenerateVertex(G4ThreeVector& vertex) { void RMGVertexFromFile::BeginOfRunAction(const G4Run*) { if (!G4Threading::IsMasterThread()) return; - G4AutoLock lock(&fMutex); - if (!fReader->GetReader()) { + if (!fReader->GetLockedReader()) { RMGLog::Out(RMGLog::fatal, "vertex file '", fReader->GetFileName(), "' not found or in wrong format"); } @@ -87,7 +79,6 @@ void RMGVertexFromFile::BeginOfRunAction(const G4Run*) { void RMGVertexFromFile::EndOfRunAction(const G4Run*) { if (!G4Threading::IsMasterThread()) return; - G4AutoLock lock(&fMutex); fReader->CloseFile(); } diff --git a/tests/vertex/CMakeLists.txt b/tests/vertex/CMakeLists.txt index dfc81b3..6aefe66 100644 --- a/tests/vertex/CMakeLists.txt +++ b/tests/vertex/CMakeLists.txt @@ -21,7 +21,8 @@ set(_macros pos-hdf5 pos-lh5 kin-lh5 kin-hdf5 pos-kin-lh5) foreach(_mac ${_macros}) add_test(NAME vertex/${_mac} COMMAND ${PYTHONPATH} run-vtx.py ${REMAGE_PYEXE} ${_mac} 1) - set_tests_properties(vertex/${_mac} PROPERTIES LABELS extra FIXTURES_REQUIRED vertex-input-fixture) + set_tests_properties(vertex/${_mac} PROPERTIES LABELS extra FIXTURES_REQUIRED + vertex-input-fixture) # note: only 4 threads available on GitHub actions. add_test(NAME vertex-mt/${_mac} COMMAND ${PYTHONPATH} run-vtx.py ${REMAGE_PYEXE} ${_mac} 4) From 38ced948b701ca4e1699e15f46eadb2eea4f6dd6 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Fri, 24 Jan 2025 15:45:53 +0100 Subject: [PATCH 6/9] add basic unit support to RMGAnalysisReader --- include/RMGAnalysisReader.hh | 30 ++++++++++++++--- include/RMGConvertLH5.hh | 8 +++-- src/RMGAnalysisReader.cc | 54 ++++++++++++++++++++++++------ src/RMGConvertLH5.cc | 18 ++++++---- src/RMGGeneratorFromFile.cc | 16 +++++---- src/RMGVertexFromFile.cc | 19 ++++++++--- src/remage-from-lh5.cc | 4 ++- tests/vertex/CMakeLists.txt | 2 +- tests/vertex/create-test-input.py | 30 ++++++++++------- tests/vertex/macros/pos-lh5-mm.mac | 13 +++++++ tests/vertex/run-vtx.py | 32 +++++++++++++++--- 11 files changed, 172 insertions(+), 54 deletions(-) create mode 100644 tests/vertex/macros/pos-lh5-mm.mac diff --git a/include/RMGAnalysisReader.hh b/include/RMGAnalysisReader.hh index 964e8cb..820026d 100644 --- a/include/RMGAnalysisReader.hh +++ b/include/RMGAnalysisReader.hh @@ -16,6 +16,7 @@ #ifndef _RMG_ANALYSIS_READER_HH_ #define _RMG_ANALYSIS_READER_HH_ +#include #include #include "G4AutoLock.hh" @@ -60,6 +61,7 @@ class RMGAnalysisReader final { void unlock() { fReader = nullptr; fNtupleId = -1; + fUnits = nullptr; if (fLock) { fLock.unlock(); } } @@ -68,20 +70,31 @@ class RMGAnalysisReader final { [[nodiscard]] auto GetNtupleRow() { return fReader->GetNtupleRow(fNtupleId); } /** * @brief wraps @ref G4VAnalysisReader::SetNtupleDColumn. */ - auto SetNtupleDColumn(const std::string& name, G4double& value) { + auto SetNtupleDColumn(const std::string& name, G4double& value, + const std::vector& allowed_units = {}) { + AssertUnit(name, allowed_units); return fReader->SetNtupleDColumn(fNtupleId, name, value); } /** * @brief wraps @ref G4VAnalysisReader::SetNtupleFColumn. */ - auto SetNtupleFColumn(const std::string& name, G4float& value) { + auto SetNtupleFColumn(const std::string& name, G4float& value, + const std::vector& allowed_units = {}) { + AssertUnit(name, allowed_units); return fReader->SetNtupleFColumn(fNtupleId, name, value); } /** * @brief wraps @ref G4VAnalysisReader::SetNtupleIColumn. */ - auto SetNtupleIColumn(const std::string& name, G4int& value) { + auto SetNtupleIColumn(const std::string& name, G4int& value, + const std::vector& allowed_units = {}) { + AssertUnit(name, allowed_units); return fReader->SetNtupleIColumn(fNtupleId, name, value); } + /** + * @brief get unit information for the column. an empty string means either no unit + * attached or no support by the file format. */ + [[nodiscard]] std::string GetUnit(const std::string& name) const; + /** * @brief check whether this access handle is still valid. */ operator bool() const { return fReader != nullptr && fNtupleId >= 0 && fLock; } @@ -89,12 +102,16 @@ class RMGAnalysisReader final { private: // only allow creation or moving in parent. - inline Access(G4AutoLock lock, G4VAnalysisReader* reader, int nt) - : fLock(std::move(lock)), fReader(reader), fNtupleId(nt) {}; + inline Access(G4AutoLock lock, G4VAnalysisReader* reader, int nt, + const std::map* u) + : fLock(std::move(lock)), fReader(reader), fNtupleId(nt), fUnits(u) {}; Access(Access&&) = default; + void AssertUnit(const std::string& name, const std::vector& allowed_units) const; + G4VAnalysisReader* fReader = nullptr; int fNtupleId = -1; + const std::map* fUnits; G4AutoLock fLock; }; @@ -146,6 +163,9 @@ class RMGAnalysisReader final { G4VAnalysisReader* fReader = nullptr; int fNtupleId = -1; + std::map fUnits{}; + bool fHasUnits = false; + std::string fFileName; bool fFileIsTemp = false; }; diff --git a/include/RMGConvertLH5.hh b/include/RMGConvertLH5.hh index 0516fa5..cda94f8 100644 --- a/include/RMGConvertLH5.hh +++ b/include/RMGConvertLH5.hh @@ -18,6 +18,7 @@ #define _RMG_CONVERT_LH5_HH #include +#include #include #include #include @@ -33,7 +34,8 @@ class RMGConvertLH5 { public: static bool ConvertToLH5(std::string, std::string, bool, bool part_of_batch = false); - static bool ConvertFromLH5(std::string, std::string, bool, bool part_of_batch = false); + static bool ConvertFromLH5(std::string, std::string, bool, bool part_of_batch, + std::map>&); inline static bool fIsStandalone = false; @@ -77,12 +79,12 @@ class RMGConvertLH5 { //////////////////////////////////////////////////////////////////////////////////////////// // LH5 -> HDF5 (input files): - bool ConvertFromLH5Internal(); + bool ConvertFromLH5Internal(std::map>&); static inline const std::regex table_dtype_re = std::regex("^table\\{.*\\}$"); std::string HDFDataTypeToForm(H5::DataType); - bool ConvertTableToNTuple(H5::Group&); + bool ConvertTableToNTuple(H5::Group&, std::map&); //////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/RMGAnalysisReader.cc b/src/RMGAnalysisReader.cc index 369429e..e795a49 100644 --- a/src/RMGAnalysisReader.cc +++ b/src/RMGAnalysisReader.cc @@ -13,26 +13,24 @@ // You should have received a copy of the GNU Lesser General Public License // along with this program. If not, see . +#include "RMGAnalysisReader.hh" + #include #include - -#include "RMGVertexFromFile.hh" namespace fs = std::filesystem; #include "G4AnalysisUtilities.hh" -#include "G4RootAnalysisReader.hh" -#include "G4Threading.hh" -#include "G4VAnalysisReader.hh" +#include "G4CsvAnalysisReader.hh" #if RMG_HAS_HDF5 #include "G4Hdf5AnalysisReader.hh" #endif -#include "G4CsvAnalysisReader.hh" +#include "G4RootAnalysisReader.hh" +#include "G4Threading.hh" #include "G4XmlAnalysisReader.hh" #if RMG_HAS_HDF5 #include "RMGConvertLH5.hh" #endif - #include "RMGLog.hh" G4Mutex RMGAnalysisReader::fMutex = G4MUTEX_INITIALIZER; @@ -43,7 +41,7 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_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); + auto invalid_access = RMGAnalysisReader::Access(std::move(invalid_lock), nullptr, -1, nullptr); // reader initialization should only happen on the master thread (otherwise it will fail). if (!G4Threading::IsMasterThread()) { @@ -62,6 +60,8 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, 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); if (ext == "root") { @@ -88,12 +88,13 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, return std::move(invalid_access); } - if (!RMGConvertLH5::ConvertFromLH5(new_fn, ntuple_dir_name, false)) { + if (!RMGConvertLH5::ConvertFromLH5(new_fn, ntuple_dir_name, false, false, units_map)) { RMGLog::Out(RMGLog::error, "Conversion of input file ", new_fn, " to LH5 failed. Data is potentially corrupted."); return std::move(invalid_access); } + fHasUnits = true; file_name = new_fn; fReader = G4Hdf5AnalysisReader::Instance(); fFileIsTemp = true; @@ -118,7 +119,8 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(std::string& file_name, RMGLog::Out(RMGLog::error, "Ntuple named '", ntuple_name, "' could not be found in input file!"); return std::move(invalid_access); } - return {std::move(lock), fReader, fNtupleId}; + fUnits = units_map[ntuple_name]; + return {std::move(lock), fReader, fNtupleId, fHasUnits ? &fUnits : nullptr}; } void RMGAnalysisReader::CloseFile() { @@ -141,11 +143,41 @@ void RMGAnalysisReader::CloseFile() { fFileName = ""; fFileIsTemp = false; fNtupleId = -1; + fHasUnits = false; + fHasUnits = {}; } RMGAnalysisReader::Access RMGAnalysisReader::GetLockedReader() const { G4AutoLock lock(&fMutex); - return {std::move(lock), fReader, fNtupleId}; + return {std::move(lock), fReader, fNtupleId, fHasUnits ? &fUnits : nullptr}; +} + +G4AutoLock RMGAnalysisReader::GetLock() const { + + // 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(); + } + + G4AutoLock lock(&fMutex); + + return std::move(lock); +} + +/* ========================================================================================== */ + +std::string RMGAnalysisReader::Access::GetUnit(const std::string& name) const { + if (!fUnits) return ""; + return fUnits->at(name); +} + +void RMGAnalysisReader::Access::AssertUnit(const std::string& name, + const std::vector& allowed_units) const { + if (!fUnits) return; + if (std::find(allowed_units.begin(), allowed_units.end(), GetUnit(name)) == allowed_units.end()) { + RMGLog::Out(RMGLog::fatal, "invalid unit '", GetUnit(name), "' for column ", name); + } } diff --git a/src/RMGConvertLH5.cc b/src/RMGConvertLH5.cc index 9fcee71..2dd8ebd 100644 --- a/src/RMGConvertLH5.cc +++ b/src/RMGConvertLH5.cc @@ -345,7 +345,8 @@ std::string RMGConvertLH5::HDFDataTypeToForm(H5::DataType dtype) { } } -bool RMGConvertLH5::ConvertTableToNTuple(H5::Group& det_group) { +bool RMGConvertLH5::ConvertTableToNTuple(H5::Group& det_group, + std::map& units_map) { const std::string ntuple_name = det_group.getObjName(); const std::string ntuple_log_prefix = "ntuple " + ntuple_name + " - "; LH5Log(RMGLog::detail, ntuple_log_prefix, "visiting"); @@ -381,7 +382,9 @@ bool RMGConvertLH5::ConvertTableToNTuple(H5::Group& det_group) { auto dset_column = det_group.openDataSet(lgdo_name); if (dset_column.attrExists("units")) { - column += "_in_" + GetStringAttribute(dset_column, "units").value(); + units_map[lgdo_name] = GetStringAttribute(dset_column, "units").value(); + } else { + units_map[lgdo_name] = ""; } names_string += column + std::string("\0", 1); forms_string += HDFDataTypeToForm(dset_column.getDataType()) + std::string("\0", 1); @@ -426,7 +429,8 @@ bool RMGConvertLH5::ConvertTableToNTuple(H5::Group& det_group) { return true; } -bool RMGConvertLH5::ConvertFromLH5Internal() { +bool RMGConvertLH5::ConvertFromLH5Internal( + std::map>& units_map) { // using the core driver with no backing storage will allow to change the file purely in-memory. // warning: this will internally allocate approx. the full file size! H5::FileAccPropList fapl; @@ -449,7 +453,8 @@ bool RMGConvertLH5::ConvertFromLH5Internal() { bool ntuple_success = true; for (auto& ntuple : ntuples) { auto det_group = ntuples_group.openGroup(ntuple); - ntuple_success &= ConvertTableToNTuple(det_group); + units_map[ntuple] = {}; + ntuple_success &= ConvertTableToNTuple(det_group, units_map[ntuple]); } if (ntuples_group.attrExists("datatype")) ntuples_group.removeAttr("datatype"); @@ -462,10 +467,11 @@ bool RMGConvertLH5::ConvertFromLH5Internal() { } bool RMGConvertLH5::ConvertFromLH5(std::string lh5_file_name, std::string ntuple_group_name, - bool dry_run, bool part_of_batch) { + bool dry_run, bool part_of_batch, + std::map>& units_map) { auto conv = RMGConvertLH5(lh5_file_name, ntuple_group_name, dry_run, part_of_batch); try { - return conv.ConvertFromLH5Internal(); + return conv.ConvertFromLH5Internal(units_map); } catch (const H5::Exception& e) { conv.LH5Log(RMGLog::error, e.getDetailMsg()); return false; diff --git a/src/RMGGeneratorFromFile.cc b/src/RMGGeneratorFromFile.cc index 8f68ba1..d8794d0 100644 --- a/src/RMGGeneratorFromFile.cc +++ b/src/RMGGeneratorFromFile.cc @@ -41,11 +41,11 @@ void RMGGeneratorFromFile::OpenFile(std::string& name) { if (!reader) return; // bind the static variables once here, and only use them later. - reader.SetNtupleIColumn("g4_pid", fRowData.fG4Pid); - reader.SetNtupleDColumn("ekin_in_MeV", fRowData.fEkin); - reader.SetNtupleDColumn("px", fRowData.fPx); - reader.SetNtupleDColumn("py", fRowData.fPy); - reader.SetNtupleDColumn("pz", fRowData.fPz); + reader.SetNtupleIColumn("g4_pid", fRowData.fG4Pid, {""}); + reader.SetNtupleDColumn("ekin", fRowData.fEkin, {"eV", "keV", "MeV", "GeV"}); + reader.SetNtupleDColumn("px", fRowData.fPx, {""}); + reader.SetNtupleDColumn("py", fRowData.fPy, {""}); + reader.SetNtupleDColumn("pz", fRowData.fPz, {""}); } void RMGGeneratorFromFile::BeginOfRunAction(const G4Run*) { @@ -84,6 +84,7 @@ void RMGGeneratorFromFile::GeneratePrimaries(G4Event* event) { // make copy of data and exit critical section. auto row_data = fRowData; + auto unit_name = locked_reader.GetUnit("ekin"); locked_reader.unlock(); // check for NaN sentinel values - i.e. non-existing columns (there is no error message). @@ -104,10 +105,13 @@ void RMGGeneratorFromFile::GeneratePrimaries(G4Event* event) { G4ThreeVector momentum{row_data.fPx, row_data.fPy, row_data.fPz}; + const std::map units = {{"", u::MeV}, {"eV", u::eV}, {"keV", u::keV}, + {"MeV", u::MeV}, {"GeV", u::GeV}}; + fGun->SetParticleDefinition(particle); fGun->SetParticlePosition(fParticlePosition); fGun->SetParticleMomentumDirection(momentum); - fGun->SetParticleEnergy(row_data.fEkin * u::MeV); + fGun->SetParticleEnergy(row_data.fEkin * units.at(unit_name)); fGun->GeneratePrimaryVertex(event); } diff --git a/src/RMGVertexFromFile.cc b/src/RMGVertexFromFile.cc index 4054e91..18c829b 100644 --- a/src/RMGVertexFromFile.cc +++ b/src/RMGVertexFromFile.cc @@ -33,9 +33,16 @@ void RMGVertexFromFile::OpenFile(std::string& name) { if (!reader) return; // bind the static variables once here, and only use them later. - reader.SetNtupleDColumn("xloc_in_m", fXpos); - reader.SetNtupleDColumn("yloc_in_m", fYpos); - reader.SetNtupleDColumn("zloc_in_m", fZpos); + reader.SetNtupleDColumn("xloc", fXpos, {"nm", "um", "mm", "cm", "m"}); + reader.SetNtupleDColumn("yloc", fYpos, {"nm", "um", "mm", "cm", "m"}); + reader.SetNtupleDColumn("zloc", fZpos, {"nm", "um", "mm", "cm", "m"}); + + const auto xunit = reader.GetUnit("xloc"); + const auto yunit = reader.GetUnit("yloc"); + const auto zunit = reader.GetUnit("zloc"); + if (xunit != yunit || xunit != zunit) { + RMGLog::OutFormat(RMGLog::fatal, " ('{}', '{}', '{}')", xunit, yunit, zunit); + } } bool RMGVertexFromFile::GenerateVertex(G4ThreeVector& vertex) { @@ -53,7 +60,11 @@ bool RMGVertexFromFile::GenerateVertex(G4ThreeVector& vertex) { return false; } - vertex = G4ThreeVector{fXpos, fYpos, fZpos} * CLHEP::m; + const auto unit_name = reader.GetUnit("xloc"); + const std::map units = {{"", CLHEP::m}, {"nm", CLHEP::nm}, + {"um", CLHEP::um}, {"mm", CLHEP::mm}, {"cm", CLHEP::cm}, {"m", CLHEP::m}}; + + vertex = G4ThreeVector{fXpos, fYpos, fZpos} * units.at(unit_name); return true; } diff --git a/src/remage-from-lh5.cc b/src/remage-from-lh5.cc index d2de39e..9ef9470 100644 --- a/src/remage-from-lh5.cc +++ b/src/remage-from-lh5.cc @@ -50,7 +50,9 @@ int main(int argc, char** argv) { RMGLog::OutFormat(RMGLog::error, "{} does not exist", file_name); continue; } - RMGConvertLH5::ConvertFromLH5(file_name, ntuple_group_name, dry_run, !file_names.empty()); + std::map> units_map{}; + RMGConvertLH5::ConvertFromLH5(file_name, ntuple_group_name, dry_run, !file_names.empty(), + units_map); } struct rusage usage{}; diff --git a/tests/vertex/CMakeLists.txt b/tests/vertex/CMakeLists.txt index 6aefe66..6689ff6 100644 --- a/tests/vertex/CMakeLists.txt +++ b/tests/vertex/CMakeLists.txt @@ -17,7 +17,7 @@ endforeach() add_test(NAME vertex/gen-input COMMAND ${PYTHONPATH} create-test-input.py ${_remage_from_lh5}) set_tests_properties(vertex/gen-input PROPERTIES LABELS extra FIXTURES_SETUP vertex-input-fixture) -set(_macros pos-hdf5 pos-lh5 kin-lh5 kin-hdf5 pos-kin-lh5) +set(_macros pos-hdf5 pos-lh5 pos-lh5-mm kin-lh5 kin-hdf5 pos-kin-lh5) foreach(_mac ${_macros}) add_test(NAME vertex/${_mac} COMMAND ${PYTHONPATH} run-vtx.py ${REMAGE_PYEXE} ${_mac} 1) diff --git a/tests/vertex/create-test-input.py b/tests/vertex/create-test-input.py index eafed34..0b554b3 100755 --- a/tests/vertex/create-test-input.py +++ b/tests/vertex/create-test-input.py @@ -25,25 +25,31 @@ def _convert_lh5_to_hdf5(lh5_fn): assert ret.returncode == 0 -INPUT_FILE_ROWS = int(1e6) +INPUT_FILE_ROWS = int(1e5) rng = np.random.default_rng(123456) # ================================================== # position input -xloc = Array(np.linspace(0, 1000, INPUT_FILE_ROWS), attrs={"units": "m"}) -yloc = Array(rng.uniform(-1, 1, INPUT_FILE_ROWS), attrs={"units": "m"}) -zloc = Array(rng.uniform(-1, 1, INPUT_FILE_ROWS), attrs={"units": "m"}) -vtx_pos_file = "macros/vtx-pos.lh5" -lh5.write( - Table({"xloc": xloc, "yloc": yloc, "zloc": zloc}), - "vtx/pos", - lh5_file=vtx_pos_file, - wo_mode="of", -) -_convert_lh5_to_hdf5(vtx_pos_file) +def _create_pos_input(dtype, vtx_pos_file, units): + attrs = {"units": units} + xloc = Array(np.linspace(0, 1000, INPUT_FILE_ROWS).astype(dtype), attrs=attrs) + yloc = Array(rng.uniform(-1, 1, INPUT_FILE_ROWS).astype(dtype), attrs=attrs) + zloc = Array(rng.uniform(-1, 1, INPUT_FILE_ROWS).astype(dtype), attrs=attrs) + + lh5.write( + Table({"xloc": xloc, "yloc": yloc, "zloc": zloc}), + "vtx/pos", + lh5_file=vtx_pos_file, + wo_mode="of", + ) + _convert_lh5_to_hdf5(vtx_pos_file) + +_create_pos_input(np.float64, "macros/vtx-pos.lh5", "m") +# _create_pos_input(np.float32, "macros/vtx-pos-float32.lh5", "m") +_create_pos_input(np.float64, "macros/vtx-pos-mm.lh5", "mm") # ================================================== # kinetcis input. diff --git a/tests/vertex/macros/pos-lh5-mm.mac b/tests/vertex/macros/pos-lh5-mm.mac new file mode 100644 index 0000000..5e0ef1c --- /dev/null +++ b/tests/vertex/macros/pos-lh5-mm.mac @@ -0,0 +1,13 @@ +/control/execute macros/_init.mac + +/RMG/Generator/Confine FromFile + +/RMG/Generator/Confinement/FromFile/FileName macros/vtx-pos-mm.lh5 + +/RMG/Generator/Select GPS +/gps/particle ion +/gps/ion 81 208 +/gps/energy 0 keV +/gps/ang/type iso + +/run/beamOn 100000 diff --git a/tests/vertex/run-vtx.py b/tests/vertex/run-vtx.py index b8874f9..f4c3010 100755 --- a/tests/vertex/run-vtx.py +++ b/tests/vertex/run-vtx.py @@ -2,6 +2,7 @@ from __future__ import annotations +import re import subprocess import sys from pathlib import Path @@ -20,6 +21,18 @@ for old_f in Path().glob(f"{output_stem}*.lh5"): Path(old_f).unlink() +# "parse" the macro. +pos_input_file = None +kin_input_file = None +with Path(f"macros/{macro}.mac").open("r") as macro_file: + for line in macro_file: + m = re.match(r"/RMG/Generator/FromFile/FileName (.*)\n$", line) + if m: + kin_input_file = m.group(1) + m = re.match(r"/RMG/Generator/Confinement/FromFile/FileName (.*)\n$", line) + if m: + pos_input_file = m.group(1) + # run remage, produce lh5 output. cmd = [ rmg, @@ -47,28 +60,37 @@ # validate output against input files. if "pos" in macro: - input_pos = lh5.read("vtx/pos", lh5_file="macros/vtx-pos.lh5").view_as("pd") + pos_input_file = pos_input_file.replace(".hdf5", ".lh5") + input_pos = lh5.read("vtx/pos", lh5_file=pos_input_file).view_as("pd") output_pos = lh5.read("stp/vertices", lh5_file=files).view_as("pd") output_pos = output_pos.sort_values("xloc") # sort by linear column. output_pos = output_pos[["xloc", "yloc", "zloc"]] uniq, cnt = np.unique(output_pos["xloc"], return_counts=True) if not np.all(cnt <= 1): - raise ValueError(uniq[cnt > 1]) + msg = f"non-unique pos 'indices' {uniq[cnt > 1]}" + raise ValueError(msg) + + input_pos = input_pos.iloc[0 : len(output_pos)] + input_pos = input_pos[["xloc", "yloc", "zloc"]] - input_pos = input_pos.iloc[0 : len(output_pos)][["xloc", "yloc", "zloc"]] + # re-scale to accommodate for different units. + if "vtx-pos-mm.lh5" in pos_input_file: + input_pos /= 1000 assert np.all(np.isclose(output_pos.to_numpy(), input_pos.to_numpy())) if "kin" in macro: - input_kin = lh5.read("vtx/kin", lh5_file="macros/vtx-kin.lh5").view_as("pd") + kin_input_file = kin_input_file.replace(".hdf5", ".lh5") + input_kin = lh5.read("vtx/kin", lh5_file=kin_input_file).view_as("pd") output_kin = lh5.read("stp/particles", lh5_file=files).view_as("pd") output_kin = output_kin.sort_values("ekin") # sort by linear column. output_kin = output_kin[["px", "py", "pz", "ekin", "particle"]] uniq, cnt = np.unique(output_kin["px"], return_counts=True) if not np.all(cnt <= 1): - raise ValueError(uniq[cnt > 1]) + msg = f"non-unique kin 'indices' {uniq[cnt > 1]}" + raise ValueError(msg) output_p_scale = np.sqrt( output_kin["px"] ** 2 + output_kin["py"] ** 2 + output_kin["pz"] ** 2 From 63e759cedf191a198c7283e25fad3fbd514a6150 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Sat, 25 Jan 2025 01:14:50 +0100 Subject: [PATCH 7/9] 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 820026d..a50d380 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 cdfdd37..63896e4 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 e795a49..78f8f93 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 039f67a..95bc776 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); } From 3f91a8c4a692322c307e32df1a0aaf0d8f402493 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Sun, 26 Jan 2025 13:04:08 +0100 Subject: [PATCH 8/9] add one more check --- include/RMGAnalysisReader.hh | 32 +++++++++++++++++++++----------- src/RMGAnalysisReader.cc | 13 ++++++++++--- 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/include/RMGAnalysisReader.hh b/include/RMGAnalysisReader.hh index a50d380..8e0525e 100644 --- a/include/RMGAnalysisReader.hh +++ b/include/RMGAnalysisReader.hh @@ -43,7 +43,8 @@ class RMGAnalysisReader final { public: /** - * @brief thread-safe access handle to the underlying reader. */ + * @brief thread-safe access handle to the underlying reader. This handle can be used to set-up + * ntuple reading (in setup mode) or to read rows from the ntuple. */ class Access final { friend class RMGAnalysisReader; @@ -67,11 +68,15 @@ class RMGAnalysisReader final { /** * @brief wraps @ref G4VAnalysisReader::GetNtupleRow. */ - [[nodiscard]] auto GetNtupleRow() { return fReader->GetNtupleRow(fNtupleId); } + [[nodiscard]] auto GetNtupleRow() { + AssertSetup(false); + return fReader->GetNtupleRow(fNtupleId); + } /** * @brief wraps @ref G4VAnalysisReader::SetNtupleDColumn. */ auto SetNtupleDColumn(const std::string& name, G4double& value, const std::vector& allowed_units = {}) { + AssertSetup(true); AssertUnit(name, allowed_units); return fReader->SetNtupleDColumn(fNtupleId, name, value); } @@ -79,6 +84,7 @@ class RMGAnalysisReader final { * @brief wraps @ref G4VAnalysisReader::SetNtupleFColumn. */ auto SetNtupleFColumn(const std::string& name, G4float& value, const std::vector& allowed_units = {}) { + AssertSetup(true); AssertUnit(name, allowed_units); return fReader->SetNtupleFColumn(fNtupleId, name, value); } @@ -86,6 +92,7 @@ class RMGAnalysisReader final { * @brief wraps @ref G4VAnalysisReader::SetNtupleIColumn. */ auto SetNtupleIColumn(const std::string& name, G4int& value, const std::vector& allowed_units = {}) { + AssertSetup(true); AssertUnit(name, allowed_units); return fReader->SetNtupleIColumn(fNtupleId, name, value); } @@ -103,16 +110,18 @@ class RMGAnalysisReader final { // only allow creation or moving in parent. inline Access(G4AutoLock lock, G4VAnalysisReader* reader, int nt, - const std::map* u) - : fLock(std::move(lock)), fReader(reader), fNtupleId(nt), fUnits(u) {}; + const std::map* u, bool setup) + : fLock(std::move(lock)), fReader(reader), fNtupleId(nt), fUnits(u), fCanSetup(setup) {}; Access(Access&&) = default; void AssertUnit(const std::string& name, const std::vector& allowed_units) const; + void AssertSetup(bool setup) const; G4VAnalysisReader* fReader = nullptr; int fNtupleId = -1; const std::map* fUnits; G4AutoLock fLock; + bool fCanSetup = false; }; RMGAnalysisReader() = default; @@ -124,10 +133,11 @@ class RMGAnalysisReader final { RMGAnalysisReader& operator=(RMGAnalysisReader&&) = delete; /** - * @brief open an input file for reading of one specific ntuple. + * @brief open an input file for reading of one specific ntuple. The return access handle can be + * used to connect to-be-read column to variables. * * @details This function can only be used on the master thread. This operation acquires a - * global lock to avoid problems. + * global lock across all readers that will be held until the access handle is discarded. * * @param file_name the input file name. the file format is determined from the file extension. * @param ntuple_dir_name the first part of the input table name. For the table addressed by @@ -146,17 +156,17 @@ class RMGAnalysisReader final { std::string ntuple_name, std::string force_ext = ""); /** - * @brief if any file is open for reading, close the reader. + * @brief if any file is open for reading, close the reader. Also clean-up temporary files. * - * @details this invalidates all readers obtained via @ref RMGAnalysisReader::GetLockedReader. - * This function can only be used on the master thread. This operation acquires a global lock to - * avoid problems. */ + * @details This function can only be used on the master thread. This operation acquires a global + * across all readers. This function will not actually free resources allocated for the reader by Geant4. */ void CloseFile(); /** * @brief get an access handle to the current underlying G4VAnalysisReader. * - * @details this acquires a global lock to avoid problems. The lock is held until the access handle is discarded. */ + * @details The return access handle can be used to read row(s) from the ntuple. This function + * acquires a global across all readers that will be held until the access handle is discarded. */ [[nodiscard]] Access GetLockedReader() const; /** diff --git a/src/RMGAnalysisReader.cc b/src/RMGAnalysisReader.cc index 78f8f93..b73f989 100644 --- a/src/RMGAnalysisReader.cc +++ b/src/RMGAnalysisReader.cc @@ -66,7 +66,8 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(const std::string& file_na // 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); + auto invalid_access = + RMGAnalysisReader::Access(std::move(invalid_lock), nullptr, -1, nullptr, false); if (fReader) return std::move(invalid_access); @@ -135,7 +136,7 @@ RMGAnalysisReader::Access RMGAnalysisReader::OpenFile(const std::string& file_na return std::move(invalid_access); } fUnits = units_map[ntuple_name]; - return {std::move(lock), fReader, fNtupleId, fHasUnits ? &fUnits : nullptr}; + return {std::move(lock), fReader, fNtupleId, fHasUnits ? &fUnits : nullptr, true}; } void RMGAnalysisReader::CloseFile() { @@ -166,7 +167,7 @@ RMGAnalysisReader::Access RMGAnalysisReader::GetLockedReader() const { G4AutoLock lock(&fMutex); - return {std::move(lock), fReader, fNtupleId, fHasUnits ? &fUnits : nullptr}; + return {std::move(lock), fReader, fNtupleId, fHasUnits ? &fUnits : nullptr, false}; } G4AutoLock RMGAnalysisReader::GetLock() const { @@ -196,3 +197,9 @@ void RMGAnalysisReader::Access::AssertUnit(const std::string& name, RMGLog::Out(RMGLog::fatal, "invalid unit '", GetUnit(name), "' for column ", name); } } + +void RMGAnalysisReader::Access::AssertSetup(bool setup) const { + if (setup != fCanSetup) { + RMGLog::Out(RMGLog::fatal, "invalid operation in mode fCanSetup=", fCanSetup); + } +} From 9186f465737856013e485bc5a316e7d9fe436ebf Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Tue, 28 Jan 2025 17:25:12 +0100 Subject: [PATCH 9/9] tests: add MUSUN example to tests --- examples/05-MUSUN/CMakeLists.txt | 11 +++++++++++ examples/05-MUSUN/main.cc | 4 +++- examples/05-MUSUN/run.mac | 6 ++---- tests/CMakeLists.txt | 4 ++++ tests/examples/CMakeLists.txt | 5 +++++ 5 files changed, 25 insertions(+), 5 deletions(-) create mode 100644 tests/examples/CMakeLists.txt diff --git a/examples/05-MUSUN/CMakeLists.txt b/examples/05-MUSUN/CMakeLists.txt index 48eef28..519909b 100644 --- a/examples/05-MUSUN/CMakeLists.txt +++ b/examples/05-MUSUN/CMakeLists.txt @@ -10,3 +10,14 @@ endif() add_executable(05-MUSUN main.cc HPGeTestStand.hh HPGeTestStand.cc) target_link_libraries(05-MUSUN PUBLIC RMG::remage) + +# collect auxiliary files +file( + GLOB _aux + RELATIVE ${PROJECT_SOURCE_DIR} + *.csv *.mac) + +# copy them to the build area +foreach(_file ${_aux}) + configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) +endforeach() diff --git a/examples/05-MUSUN/main.cc b/examples/05-MUSUN/main.cc index 7632fda..95324cc 100644 --- a/examples/05-MUSUN/main.cc +++ b/examples/05-MUSUN/main.cc @@ -11,6 +11,7 @@ int main(int argc, char** argv) { manager.SetUserInit(new HPGeTestStand()); manager.SetNumberOfThreads(2); manager.EnablePersistency(); + manager.SetOutputOverwriteFiles(true); manager.GetDetectorConstruction()->RegisterDetector(kGermanium, "HPGe1", 0); manager.GetDetectorConstruction()->RegisterDetector(kGermanium, "HPGe2", 1); manager.GetDetectorConstruction()->RegisterDetector(kGermanium, "HPGe3", 2); @@ -22,5 +23,6 @@ int main(int argc, char** argv) { manager.Initialize(); manager.Run(); - return 0; + if (manager.HadError()) return 1; + return manager.HadWarning() ? 2 : 0; } diff --git a/examples/05-MUSUN/run.mac b/examples/05-MUSUN/run.mac index a59b48a..198951b 100644 --- a/examples/05-MUSUN/run.mac +++ b/examples/05-MUSUN/run.mac @@ -6,11 +6,9 @@ /RMG/Output/FileName output.root - - /RMG/Generator/Confine UnConfined /RMG/Generator/Select MUSUNCosmicMuons -/RMG/Generator/MUSUNCosmicMuons/SetMUSUNFile MUSUN_10k_events.csv +/RMG/Generator/MUSUNCosmicMuons/SetMUSUNFile MUSUN_100_events.csv -/run/beamOn 1000 +/run/beamOn 100 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 67bfe4c..86a41a7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -15,3 +15,7 @@ add_subdirectory(internals) add_subdirectory(output) add_subdirectory(python) add_subdirectory(vertex) + +if(RMG_BUILD_EXAMPLES) + add_subdirectory(examples) +endif() diff --git a/tests/examples/CMakeLists.txt b/tests/examples/CMakeLists.txt new file mode 100644 index 0000000..1ec935d --- /dev/null +++ b/tests/examples/CMakeLists.txt @@ -0,0 +1,5 @@ +add_test( + NAME examples/05-MUSUN + WORKING_DIRECTORY $ + COMMAND 05-MUSUN run.mac) +set_tests_properties(examples/05-MUSUN PROPERTIES LABELS extra)