diff --git a/.gitignore b/.gitignore index f5c4ebc..6a3549d 100644 --- a/.gitignore +++ b/.gitignore @@ -6,7 +6,9 @@ /docs/api /docs/extra build +install .local +.vscode # G4 outputs *.root @@ -25,6 +27,7 @@ compile_commands.json *.swo .mypy_cache __pycache__ +*.dat # python .venv diff --git a/examples/05-MUSUN/CMakeLists.txt b/examples/05-MUSUN/CMakeLists.txt new file mode 100644 index 0000000..32c7529 --- /dev/null +++ b/examples/05-MUSUN/CMakeLists.txt @@ -0,0 +1,9 @@ +cmake_minimum_required(VERSION 3.8) +project(05-MUSUN) + +set(CMAKE_BUILD_TYPE RelWithDebInfo) + +find_package(remage) + +add_executable(05-MUSUN main.cc HPGeTestStand.hh HPGeTestStand.cc) +target_link_libraries(05-MUSUN PUBLIC RMG::remage) diff --git a/examples/05-MUSUN/HPGeTestStand.cc b/examples/05-MUSUN/HPGeTestStand.cc new file mode 100644 index 0000000..fccc53c --- /dev/null +++ b/examples/05-MUSUN/HPGeTestStand.cc @@ -0,0 +1,81 @@ +#include "HPGeTestStand.hh" + +#include "G4Box.hh" +#include "G4LogicalVolume.hh" +#include "G4Material.hh" +#include "G4NistManager.hh" +#include "G4PVPlacement.hh" +#include "G4UnitsTable.hh" + +namespace u = CLHEP; + +G4VPhysicalVolume* HPGeTestStand::DefineGeometry() { + + G4State state; + std::string name, symbol; + std::vector elements; + std::vector mass_fraction; + double density; + double temperature; + double pressure; + double abundance; + int n_isotopes; + int n_components; + int n_atoms; + + auto man = G4NistManager::Instance(); + + // define enriched germanium + auto Ge70 = new G4Isotope(name = "Ge70", 32, 70, 69.92 * u::g / u::mole); + auto Ge72 = new G4Isotope(name = "Ge72", 32, 72, 71.92 * u::g / u::mole); + auto Ge73 = new G4Isotope(name = "Ge73", 32, 73, 73.00 * u::g / u::mole); + auto Ge74 = new G4Isotope(name = "Ge74", 32, 74, 74.00 * u::g / u::mole); + auto Ge76 = new G4Isotope(name = "Ge76", 32, 76, 76.00 * u::g / u::mole); + + auto el_enr_ge = new G4Element(name = "EnrichedGermanium", symbol = "EnrGe", n_isotopes = 5); + el_enr_ge->AddIsotope(Ge70, abundance = 0.0 * u::perCent); + el_enr_ge->AddIsotope(Ge72, abundance = 0.1 * u::perCent); + el_enr_ge->AddIsotope(Ge73, abundance = 0.2 * u::perCent); + el_enr_ge->AddIsotope(Ge74, abundance = 13.1 * u::perCent); + el_enr_ge->AddIsotope(Ge76, abundance = 86.6 * u::perCent); + + auto LAr = man->FindOrBuildMaterial("G4_lAr"); + auto worldMaterial = man->FindOrBuildMaterial("G4_Galactic"); + + auto mat_enr_ge = new G4Material("CryogenicEnrichedGermanium", + density = 5.56 * u::g / (u::cm * u::cm * u::cm), n_components = 1, state = G4State::kStateSolid, + temperature = LAr->GetTemperature(), pressure = LAr->GetPressure()); + + mat_enr_ge->AddElement(el_enr_ge, n_atoms = 1); + + auto world_s = new G4Box("WorldLAr", 20 * u::m, 20 * u::m, 20 * u::m); + + auto world_l = new G4LogicalVolume(world_s, worldMaterial, "WorldLAr"); + + auto world_p = new G4PVPlacement(nullptr, G4ThreeVector(), world_l, "World", 0, false, 0); + + auto hpge_s = new G4Box("HPGe", 5 * u::m, 5 * u::m, 5 * u::cm); + + auto hpge1_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe1"); + auto hpge2_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe2"); + auto hpge3_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe3"); + auto hpge4_l = + new G4LogicalVolume(hpge_s, G4Material::GetMaterial("CryogenicEnrichedGermanium"), "HPGe4"); + + auto spacing = 5.5 * u::m; + new G4PVPlacement(nullptr, G4ThreeVector(+spacing, +spacing, 0), hpge1_l, "HPGe1", world_l, false, + 0); + new G4PVPlacement(nullptr, G4ThreeVector(+spacing, -spacing, 0), hpge2_l, "HPGe2", world_l, false, + 0); + new G4PVPlacement(nullptr, G4ThreeVector(-spacing, +spacing, 0), hpge3_l, "HPGe3", world_l, false, + 0); + new G4PVPlacement(nullptr, G4ThreeVector(-spacing, -spacing, 0), hpge4_l, "HPGe4", world_l, false, + 0); + + return world_p; +} + +// vim: tabstop=2 textwidth=2 expandtab diff --git a/examples/05-MUSUN/HPGeTestStand.hh b/examples/05-MUSUN/HPGeTestStand.hh new file mode 100644 index 0000000..6bd3bd1 --- /dev/null +++ b/examples/05-MUSUN/HPGeTestStand.hh @@ -0,0 +1,13 @@ +#ifndef _HPGE_TEST_STAND_HH_ +#define _HPGE_TEST_STAND_HH_ + +#include "RMGHardware.hh" + +class HPGeTestStand : public RMGHardware { + + public: + + G4VPhysicalVolume* DefineGeometry(); +}; + +#endif diff --git a/examples/05-MUSUN/MUSUN_100_events.csv b/examples/05-MUSUN/MUSUN_100_events.csv new file mode 100644 index 0000000..8f28042 --- /dev/null +++ b/examples/05-MUSUN/MUSUN_100_events.csv @@ -0,0 +1,100 @@ + 1 10 526.0 800.0 768.5 -45.8 1.0479 3.2079 + 2 10 57.2 800.0 -432.6 276.2 0.9065 2.2608 + 3 11 525.2 -739.0 -662.9 1000.0 0.4317 2.5394 + 4 11 77.2 -373.1 -263.0 1000.0 0.6689 1.2178 + 5 10 113.9 571.8 800.0 942.6 0.6654 4.7304 + 6 11 153.6 -800.0 45.5 -655.6 1.1282 5.9799 + 7 10 27.6 127.2 -800.0 807.0 0.4609 2.3555 + 8 11 32.9 -221.2 -747.5 1000.0 0.1594 6.2084 + 9 10 34.6 339.8 -800.0 748.2 0.7838 2.7358 + 10 11 128.7 788.9 800.0 886.0 0.5197 4.5260 + 11 10 63.8 167.6 185.5 1000.0 0.2098 1.7774 + 12 11 448.5 -607.0 724.4 1000.0 0.8797 0.6395 + 13 11 135.7 -788.2 -603.2 1000.0 0.3784 0.6812 + 14 10 189.7 800.0 -523.8 -865.0 0.8261 3.5386 + 15 10 783.2 -217.9 -800.0 832.6 0.9215 1.8035 + 16 11 211.4 -800.0 -607.7 914.8 0.8207 5.7886 + 17 11 548.3 -517.5 720.9 1000.0 0.3760 0.8090 + 18 11 1077.3 193.2 495.9 1000.0 0.2487 1.3030 + 19 10 209.6 723.1 -100.0 1000.0 0.9276 0.2391 + 20 11 244.6 -123.4 -800.0 -487.7 0.9416 2.6353 + 21 11 84.7 -336.9 -61.9 1000.0 0.5787 0.8478 + 22 10 31.8 -136.6 -353.4 1000.0 0.2543 1.9878 + 23 10 19.6 -800.0 -407.0 762.3 0.3635 0.1139 + 24 11 110.6 106.8 562.6 1000.0 0.1193 1.8106 + 25 10 208.4 -83.3 -368.8 1000.0 0.9640 6.2737 + 26 10 102.6 -370.4 476.9 1000.0 0.3535 4.2869 + 27 11 116.8 -287.3 -4.3 1000.0 0.3400 5.3986 + 28 10 133.8 -800.0 -718.0 150.2 0.2875 6.1861 + 29 11 117.2 -237.3 -800.0 -419.4 0.5466 0.9147 + 30 10 116.8 -121.9 61.4 1000.0 0.6917 1.8001 + 31 11 182.1 -588.4 -800.0 851.4 0.8999 1.1730 + 32 10 242.0 126.2 800.0 423.3 1.0170 5.8627 + 33 11 115.2 411.8 -10.0 1000.0 0.5997 0.4238 + 34 10 105.1 800.0 -380.0 -0.2 1.0274 2.7817 + 35 11 28.7 -522.7 769.1 1000.0 0.5153 3.3020 + 36 11 804.6 800.0 351.8 671.7 1.2103 2.1982 + 37 11 158.9 -245.7 -800.0 182.8 0.8884 0.7329 + 38 11 156.7 174.0 -448.2 1000.0 0.5366 0.6615 + 39 10 125.4 800.0 -323.6 421.2 0.7321 3.3997 + 40 11 154.6 -587.6 -383.9 1000.0 0.4975 1.0033 + 41 10 100.3 -518.0 800.0 5.4 0.5021 5.8716 + 42 11 38.6 -539.5 -754.6 1000.0 0.6915 3.1957 + 43 10 164.4 800.0 390.0 125.9 0.3870 3.8795 + 44 10 497.4 -800.0 666.1 -441.6 0.8870 0.5854 + 45 11 538.9 -479.8 -800.0 953.7 0.7508 2.3873 + 46 10 113.8 800.0 -699.1 587.9 0.4458 2.5153 + 47 10 89.5 800.0 -604.8 -674.5 0.6067 3.7204 + 48 11 170.8 590.7 -321.3 1000.0 0.6105 1.8563 + 49 11 32.2 800.0 -390.0 -431.0 0.4343 3.5999 + 50 10 50.1 -433.0 -800.0 639.5 0.8544 0.1427 + 51 10 152.6 800.0 -376.8 187.2 0.7371 3.9973 + 52 10 244.2 -391.0 800.0 391.9 0.5853 3.5183 + 53 10 1207.8 800.0 583.4 863.1 1.3995 2.7492 + 54 11 60.6 800.0 420.7 210.7 0.8514 3.6544 + 55 11 196.0 800.0 -564.6 -677.7 0.6648 3.3774 + 56 10 29.4 -222.5 -613.6 1000.0 0.8477 1.4582 + 57 10 509.4 -800.0 -200.8 574.5 0.7651 4.9589 + 58 11 398.2 800.0 -530.7 482.6 0.8029 2.4177 + 59 11 155.8 -364.4 -139.0 1000.0 0.8622 2.0692 + 60 10 129.5 -403.3 -642.0 1000.0 0.4503 0.8428 + 61 10 279.8 -284.6 -800.0 403.4 0.6222 2.5170 + 62 11 990.4 -79.9 85.4 1000.0 0.7738 6.2719 + 63 11 26.9 -124.3 -800.0 -622.6 1.2285 1.8099 + 64 11 112.0 -800.0 110.5 -416.8 0.9586 0.5349 + 65 10 29.1 253.0 392.0 1000.0 0.6087 0.7518 + 66 10 123.4 -731.5 800.0 672.0 1.2164 3.3188 + 67 10 135.7 -800.0 -162.3 291.2 0.6749 5.5593 + 68 10 378.1 -800.0 542.7 -9.0 0.4532 6.2200 + 69 11 158.6 40.6 -636.5 1000.0 0.4328 0.0397 + 70 10 1421.5 -108.5 -800.0 -718.5 0.7949 1.3889 + 71 11 385.8 198.9 -666.0 1000.0 1.0735 2.6832 + 72 10 238.9 -58.9 214.0 1000.0 0.6053 1.4336 + 73 10 83.8 -800.0 545.4 162.7 0.7204 0.0742 + 74 10 1.6 786.2 380.7 1000.0 0.7064 2.6703 + 75 10 258.2 -800.0 15.4 437.1 0.8364 0.9567 + 76 10 9.8 -800.0 -322.2 64.8 0.2215 5.7997 + 77 10 962.0 -800.0 -55.1 413.8 0.5464 5.9295 + 78 10 169.3 -145.6 -180.1 1000.0 0.5768 0.1291 + 79 10 65.4 -535.7 -800.0 -89.1 0.8314 1.4677 + 80 10 73.7 705.3 184.9 1000.0 0.7741 4.8536 + 81 11 127.4 -800.0 -3.2 674.7 0.7153 0.4078 + 82 11 223.5 800.0 478.0 946.1 1.1002 3.5350 + 83 10 54.7 -629.3 321.3 1000.0 0.6905 3.9056 + 84 11 284.9 514.3 -800.0 -778.5 0.7618 1.1300 + 85 10 71.4 373.8 -326.4 1000.0 0.5017 1.3038 + 86 11 277.3 -90.7 722.4 1000.0 0.2517 4.1904 + 87 10 94.5 226.0 -601.3 1000.0 0.4027 5.7292 + 88 10 59.4 411.4 -800.0 786.0 0.6348 2.4486 + 89 10 234.6 -656.9 305.7 1000.0 0.5538 0.7453 + 90 10 2.6 795.7 800.0 432.0 0.7128 4.6691 + 91 10 5.0 800.0 -118.2 95.9 1.2710 3.1173 + 92 10 126.5 646.3 -800.0 868.6 0.9698 2.7141 + 93 11 573.0 -800.0 -226.7 -726.8 0.4192 5.2811 + 94 10 370.5 -800.0 -163.3 -201.6 0.4569 0.7669 + 95 11 113.7 -630.1 -190.5 1000.0 0.2919 1.3824 + 96 10 194.1 -800.0 -414.2 355.5 0.4413 5.5602 + 97 10 13.9 200.4 750.0 1000.0 0.1852 4.8597 + 98 10 151.2 800.0 -650.1 19.4 0.9356 3.2461 + 99 10 82.2 -179.1 574.5 1000.0 0.3053 1.6857 + 100 10 50.1 -459.1 403.4 1000.0 0.4999 2.7358 diff --git a/examples/05-MUSUN/main.cc b/examples/05-MUSUN/main.cc new file mode 100644 index 0000000..18b7cbe --- /dev/null +++ b/examples/05-MUSUN/main.cc @@ -0,0 +1,26 @@ +#include "RMGLog.hh" +#include "RMGManager.hh" + +#include "HPGeTestStand.hh" + +int main(int argc, char** argv) { + + // RMGLog::SetLogLevel(RMGLog::debug); + + RMGManager manager("05-MUSUN", argc, argv); + manager.SetUserInit(new HPGeTestStand()); + manager.SetNumberOfThreads(2); + manager.EnablePersistency(); + manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe1", 0); + manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe2", 1); + manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe3", 2); + manager.GetDetectorConstruction()->RegisterDetector(RMGHardware::kGermanium, "HPGe4", 3); + + std::string macro = argc > 1 ? argv[1] : ""; + if (!macro.empty()) manager.IncludeMacroFile(macro); + + manager.Initialize(); + manager.Run(); + + return 0; +} diff --git a/examples/05-MUSUN/run.mac b/examples/05-MUSUN/run.mac new file mode 100644 index 0000000..a59b48a --- /dev/null +++ b/examples/05-MUSUN/run.mac @@ -0,0 +1,16 @@ +# /RMG/Manager/Logging/LogLevel debug + +#/tracking/verbose 1 + +/run/initialize + +/RMG/Output/FileName output.root + + + +/RMG/Generator/Confine UnConfined +/RMG/Generator/Select MUSUNCosmicMuons + +/RMG/Generator/MUSUNCosmicMuons/SetMUSUNFile MUSUN_10k_events.csv + +/run/beamOn 1000 diff --git a/examples/05-MUSUN/vis-cutaway.mac b/examples/05-MUSUN/vis-cutaway.mac new file mode 100644 index 0000000..822dda7 --- /dev/null +++ b/examples/05-MUSUN/vis-cutaway.mac @@ -0,0 +1,14 @@ +/run/initialize + +/vis/open OIQt +/vis/scene/create + +/vis/viewer/set/style surface +/vis/viewer/set/background white + +/vis/viewer/addCutawayPlane 0 0 0 m 1 0 0 + +/vis/drawVolume + +/vis/sceneHandler/attach +/vis/viewer/flush diff --git a/examples/05-MUSUN/vis-traj.mac b/examples/05-MUSUN/vis-traj.mac new file mode 100644 index 0000000..96cb0f6 --- /dev/null +++ b/examples/05-MUSUN/vis-traj.mac @@ -0,0 +1,22 @@ +/RMG/Manager/Logging/LogLevel debug + +/run/initialize + +/vis/open OGL #OIQt +/vis/scene/create +/vis/sceneHandler/attach + +/vis/drawVolume + +/vis/viewer/set/defaultColour black +/vis/viewer/set/background white + +/vis/scene/add/trajectories smooth +/vis/scene/add/hits +/vis/scene/endOfEventAction accumulate + +/RMG/Generator/Select MUSUNCosmicMuons + +/RMG/Generator/MUSUNCosmicMuons/SetMUSUNFile MUSUN_10k_events.dat + +/run/beamOn 10 diff --git a/include/RMGGeneratorMUSUNCosmicMuons.hh b/include/RMGGeneratorMUSUNCosmicMuons.hh new file mode 100644 index 0000000..9dec8e0 --- /dev/null +++ b/include/RMGGeneratorMUSUNCosmicMuons.hh @@ -0,0 +1,71 @@ +#ifndef _RMG_GENERATOR_MUSUN_COSMIC_MUONS_HH_ +#define _RMG_GENERATOR_MUSUN_COSMIC_MUONS_HH_ + +#include "RMGVGenerator.hh" +#include "RMGVVertexGenerator.hh" + +#include "CLHEP/Units/SystemOfUnits.h" +#include "G4GenericMessenger.hh" +#include "G4ParticleGun.hh" +#include "G4VUserPrimaryGeneratorAction.hh" + +#include "G4CsvAnalysisReader.hh" +#include + +namespace u = CLHEP; + +struct RMGGeneratorMUSUNCosmicMuons_Data { + G4int fID; + G4int fType; + G4double fEkin; + G4double fX; + G4double fY; + G4double fZ; + G4double fTheta; + G4double fPhi; + G4double fPx; + G4double fPy; + G4double fPz; +}; + + +class G4Event; +class RMGGeneratorMUSUNCosmicMuons : public RMGVGenerator { + + public: + + RMGGeneratorMUSUNCosmicMuons(); + ~RMGGeneratorMUSUNCosmicMuons() = default; + + RMGGeneratorMUSUNCosmicMuons(RMGGeneratorMUSUNCosmicMuons const&) = delete; + RMGGeneratorMUSUNCosmicMuons& operator=(RMGGeneratorMUSUNCosmicMuons const&) = delete; + RMGGeneratorMUSUNCosmicMuons(RMGGeneratorMUSUNCosmicMuons&&) = delete; + RMGGeneratorMUSUNCosmicMuons& operator=(RMGGeneratorMUSUNCosmicMuons&&) = delete; + + void GeneratePrimaries(G4Event* event); + virtual void SetParticlePosition(G4ThreeVector) override{}; + + + void BeginOfRunAction(const G4Run*); + void EndOfRunAction(const G4Run*); + + private: + + void DefineCommands(); + void SetMUSUNFile(G4String pathToFile); + void PrepareCopy(G4String pathToFile); + + std::unique_ptr fGun = nullptr; + std::unique_ptr fMessenger = nullptr; + G4String fPathToFile = ""; + std::filesystem::path fPathToTmpFolder; + G4String fPathToTmpFile = ""; + + static G4CsvAnalysisReader* fAnalysisReader; + + static RMGGeneratorMUSUNCosmicMuons_Data* input_data; +}; + +#endif + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/include/RMGMasterGenerator.hh b/include/RMGMasterGenerator.hh index 898c5eb..8777869 100644 --- a/include/RMGMasterGenerator.hh +++ b/include/RMGMasterGenerator.hh @@ -41,6 +41,7 @@ class RMGMasterGenerator : public G4VUserPrimaryGeneratorAction { kGPS, kBxDecay0, kCosmicMuons, + kMUSUNCosmicMuons, kUserDefined, kUndefined }; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b01a70a..7e5e6c7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,7 @@ set(PROJECT_PUBLIC_HEADERS ${_root}/include/RMGExceptionHandler.hh ${_root}/include/RMGEventAction.hh ${_root}/include/RMGGeneratorCosmicMuons.hh + ${_root}/include/RMGGeneratorMUSUNCosmicMuons.hh ${_root}/include/RMGGeneratorG4Gun.hh ${_root}/include/RMGGeneratorGPS.hh ${_root}/include/RMGGeneratorUtil.hh @@ -38,6 +39,7 @@ set(PROJECT_SOURCES ${_root}/src/RMGHardwareMessenger.cc ${_root}/src/RMGEventAction.cc ${_root}/src/RMGGeneratorCosmicMuons.cc + ${_root}/src/RMGGeneratorMUSUNCosmicMuons.cc ${_root}/src/RMGGeneratorUtil.cc ${_root}/src/RMGGermaniumDetector.cc ${_root}/src/RMGGermaniumOutputScheme.cc diff --git a/src/RMGGeneratorMUSUNCosmicMuons.cc b/src/RMGGeneratorMUSUNCosmicMuons.cc new file mode 100644 index 0000000..1a73d2c --- /dev/null +++ b/src/RMGGeneratorMUSUNCosmicMuons.cc @@ -0,0 +1,204 @@ +#include "RMGGeneratorMUSUNCosmicMuons.hh" + +#include +#include +#include + +#include "G4GenericMessenger.hh" +#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" + +#include "G4AutoLock.hh" +#include "G4CsvAnalysisReader.hh" + +namespace u = CLHEP; +namespace { + G4Mutex RMGGeneratorMUSUNCosmicMuonsDestrMutex = G4MUTEX_INITIALIZER; + G4Mutex RMGGeneratorMUSUNCosmicMuonsMutex = G4MUTEX_INITIALIZER; +} // namespace + +G4CsvAnalysisReader* RMGGeneratorMUSUNCosmicMuons::fAnalysisReader = 0; +RMGGeneratorMUSUNCosmicMuons_Data* RMGGeneratorMUSUNCosmicMuons::input_data = 0; + +RMGGeneratorMUSUNCosmicMuons::RMGGeneratorMUSUNCosmicMuons() : RMGVGenerator("MUSUNCosmicMuons") { + this->DefineCommands(); + fGun = std::make_unique(); + fPathToTmpFolder = std::filesystem::temp_directory_path(); +} + +void RMGGeneratorMUSUNCosmicMuons::PrepareCopy(G4String 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. + To do this, we copy the file to /var/tmp with the appropriate header. + To determine the header format, we need to determine the number of columns. + */ + + // Define fPathToTmpFile + std::filesystem::path originalFilePath((std::string)pathToFile); + std::filesystem::path fileName = originalFilePath.filename(); + fPathToTmpFile = (G4String)(fPathToTmpFolder / + fileName); //.substr(0, fileName.find_last_of(".")) + "_nt_MUSUN.csv"; + + // 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 (std::filesystem::exists((std::string)fPathToTmpFile)) { + RMGLog::Out(RMGLog::warning, "Temporary file already exists. Deleting it."); + std::filesystem::remove((std::string)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; + return; + } + std::istringstream iss(firstLine); + std::vector tokens(std::istream_iterator{iss}, + std::istream_iterator()); + int numColumns = tokens.size(); + + // Define header template + std::string header_template = "#class tools::wcsv::ntuple\n" + "#title MUSUN\n" + "#separator 11\n" + "#vector_separator 10\n" + "#column int ID\n" + "#column int type\n" + "#column double Ekin\n" + "#column double x\n" + "#column double y\n" + "#column double z\n"; + + // Based on the number of columns, add additional columns + if (numColumns == 8) { + header_template += "#column double theta\n" + "#column double phi\n"; + } else if (numColumns == 9) { + header_template += "#column double px\n" + "#column double py\n" + "#column double pz\n"; + } else + RMGLog::Out(RMGLog::fatal, + "MUSUN format not identified! It has " + to_string(numColumns) + " columns. Exit."); + + + // Create a temporary file and write the header + std::ofstream tmpFile(fPathToTmpFile); + if (!tmpFile) RMGLog::Out(RMGLog::fatal, "Unable to create temporary file! Exit."); + + tmpFile << header_template; + + // Copy contents from original file to temporary file + tmpFile << firstLine << std::endl; + tmpFile << originalFile.rdbuf(); + + // Close files + originalFile.close(); + tmpFile.close(); +} + +void RMGGeneratorMUSUNCosmicMuons::BeginOfRunAction(const G4Run*) { + + G4AutoLock lock(&RMGGeneratorMUSUNCosmicMuonsMutex); + if (!fAnalysisReader) { + 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)); + } + 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; + } +} + + +void RMGGeneratorMUSUNCosmicMuons::GeneratePrimaries(G4Event* event) { + + G4AutoLock lock(&RMGGeneratorMUSUNCosmicMuonsMutex); + + fAnalysisReader->GetNtupleRow(); + + G4ParticleTable* 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}); + + 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.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); + } else { + 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); + } + + RMGLog::OutFormat(RMGLog::debug, "...energy {:.4g} GeV", input_data->fEkin); + fGun->SetParticleEnergy(input_data->fEkin * u::GeV); + + fGun->GeneratePrimaryVertex(event); +} + +void RMGGeneratorMUSUNCosmicMuons::SetMUSUNFile(G4String pathToFile) { fPathToFile = pathToFile; } + +void RMGGeneratorMUSUNCosmicMuons::DefineCommands() { + + // NOTE: SetUnit(Category) is not thread-safe + + fMessenger = std::make_unique(this, "/RMG/Generator/MUSUNCosmicMuons/", + "Commands for controlling the MUSUN µ generator"); + + fMessenger->DeclareMethod("SetMUSUNFile", &RMGGeneratorMUSUNCosmicMuons::SetMUSUNFile) + .SetGuidance("Set the MUSUN input file") + .SetParameterName("MUSUNFileName", false) + .SetToBeBroadcasted(true) + .SetStates(G4State_PreInit, G4State_Idle); +} + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/src/RMGMasterGenerator.cc b/src/RMGMasterGenerator.cc index 04e08ec..6409149 100644 --- a/src/RMGMasterGenerator.cc +++ b/src/RMGMasterGenerator.cc @@ -26,6 +26,7 @@ #include "RMGGeneratorDecay0.hh" #endif #include "RMGGeneratorCosmicMuons.hh" +#include "RMGGeneratorMUSUNCosmicMuons.hh" #include "RMGLog.hh" #include "RMGVertexFromFile.hh" @@ -110,6 +111,9 @@ void RMGMasterGenerator::SetGenerator(RMGMasterGenerator::Generator gen) { case Generator::kCosmicMuons: fGeneratorObj = std::make_unique(); break; + case Generator::kMUSUNCosmicMuons: + fGeneratorObj = std::make_unique(); + break; case Generator::kUndefined: case Generator::kUserDefined: break; default: