Implement musun input generator

.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

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)

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<std::string> elements;
+  std::vector<double> 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

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

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

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;
+}

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

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

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

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 <filesystem>
+
+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 vec) override {};
+
+
+    void BeginOfRunAction(const G4Run*);
+    void EndOfRunAction(const G4Run*);
+
+  private:
+
+    void DefineCommands();
+    void SetMUSUNFile(G4String pathToFile);
+    void PrepareCopy(G4String pathToFile);
+
+    std::unique_ptr<G4ParticleGun> fGun = nullptr;
+    std::unique_ptr<G4GenericMessenger> 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

include/RMGMasterGenerator.hh

@@ -41,6 +41,7 @@ class RMGMasterGenerator : public G4VUserPrimaryGeneratorAction {
       kGPS,
       kBxDecay0,
       kCosmicMuons,
+      kMUSUNCosmicMuons,
       kUserDefined,
       kUndefined
     };