diff --git a/DQM/include/DQM/GenieTruthDQM.h b/DQM/include/DQM/GenieTruthDQM.h new file mode 100644 index 000000000..127e887f9 --- /dev/null +++ b/DQM/include/DQM/GenieTruthDQM.h @@ -0,0 +1,57 @@ +// +// Created by Wesley Ketchum on 8/1/24. +// + +#ifndef DQM_GENIETRUTH_H +#define DQM_GENIETRUTH_H + +// LDMX Framework +#include "Framework/Configure/Parameters.h" // Needed to import parameters from configuration file +#include "Framework/EventProcessor.h" //Needed to declare processor + +namespace dqm { + +/** + * @class GenieTruthDQM + * @brief Generate histograms/ntuple to extract genie output info + */ +class GenieTruthDQM : public framework::Analyzer { + public: + /** + * Constructor + * + * Blank Analyzer constructor + */ + GenieTruthDQM(const std::string& name, framework::Process& process) + : framework::Analyzer(name, process) {} + + /** + * Input python configuration parameters + */ + virtual void configure(framework::config::Parameters& ps); + + /** + * Construct histograms/ntuples + */ + virtual void onProcessStart(); + + /** + * Grab the run number... + */ + virtual void onNewRun(const ldmx::RunHeader& runHeader); + + /** + * Fills histograms/ntuples + */ + virtual void analyze(const framework::Event& event); + + private: + /// Pass Name for genie objects + std::string hepmc3CollName_; + std::string hepmc3PassName_; + + int run_number_; +}; +} // namespace dqm + +#endif // DQM_GENIETRUTH_H diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index 3d66bc76d..381fefef5 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -680,7 +680,24 @@ def __init__(self, name='SampleValidation') : # daughters of hard brem histograms self.build1DHistogram("pdgid_hardbremdaughters", "ID of hard brem daughters", 20, 0, 20) self.build1DHistogram("startZ_hardbremdaughters", "Start z position of hard brem daughters [mm]", 200, -1000, 1000) - + +class GenieTruthDQM(ldmxcfg.Analyzer) : + """Configured GenieTruthDQM python object + + Contains an instance of GenieTruthDQM that + has already been configured. + + Examples + -------- + from LDMX.DQM import dqm + p.sequence.append( dqm.GenieTruthDQM() ) + """ + + def __init__(self,name='GenieTruthDQM',coll_name="",pass_name="") : + super().__init__(name,'dqm::GenieTruthDQM','DQM') + + self.hepmc3CollName = coll_name + self.hepmc3PassName = pass_name ecal_dqm = [ EcalDigiVerify(), @@ -731,5 +748,4 @@ def __init__(self, name='SampleValidation') : Trigger() ] - all_dqm = ecal_dqm + hcal_dqm + recoil_dqm + trigScint_dqm + trigger_dqm diff --git a/DQM/src/DQM/GenieTruthDQM.cxx b/DQM/src/DQM/GenieTruthDQM.cxx new file mode 100644 index 000000000..7aa3cd9ca --- /dev/null +++ b/DQM/src/DQM/GenieTruthDQM.cxx @@ -0,0 +1,203 @@ +// +// Created by Wesley Ketchum on 8/1/24. +// + +#include "DQM/GenieTruthDQM.h" + +#include + +#include "GENIE/Framework/Conventions/KineVar.h" +#include "SimCore/Event/HepMC3GenEvent.h" + +namespace dqm { + +void GenieTruthDQM::configure(framework::config::Parameters& ps) { + hepmc3CollName_ = ps.getParameter("hepmc3CollName"); + hepmc3PassName_ = ps.getParameter("hepmc3PassName"); + return; +} + +void GenieTruthDQM::onProcessStart() { + getHistoDirectory(); + + ntuple_.create("genie_truth"); + + // event info + ntuple_.addVar("genie_truth", "run"); + ntuple_.addVar("genie_truth", "event"); + + ntuple_.addVar("genie_truth", "interaction_type"); + ntuple_.addVar("genie_truth", "scattering_type"); + ntuple_.addVar("genie_truth", "rescatter_code"); + + ntuple_.addVar("genie_truth", "x_bj"); + ntuple_.addVar("genie_truth", "y_inel"); + ntuple_.addVar("genie_truth", "Q2"); + ntuple_.addVar("genie_truth", "W"); + + ntuple_.addVar("genie_truth", "lep_pdg"); + + ntuple_.addVar("genie_truth", "lep_i_px"); + ntuple_.addVar("genie_truth", "lep_i_py"); + ntuple_.addVar("genie_truth", "lep_i_pz"); + ntuple_.addVar("genie_truth", "lep_i_e"); + + ntuple_.addVar("genie_truth", "lep_f_px"); + ntuple_.addVar("genie_truth", "lep_f_py"); + ntuple_.addVar("genie_truth", "lep_f_pz"); + ntuple_.addVar("genie_truth", "lep_f_e"); + + ntuple_.addVar("genie_truth", "tgt_pdg"); + ntuple_.addVar("genie_truth", "tgt_px"); + ntuple_.addVar("genie_truth", "tgt_py"); + ntuple_.addVar("genie_truth", "tgt_pz"); + ntuple_.addVar("genie_truth", "tgt_e"); + + ntuple_.addVar("genie_truth", "hnuc_pdg"); + ntuple_.addVar("genie_truth", "hnuc_px"); + ntuple_.addVar("genie_truth", "hnuc_py"); + ntuple_.addVar("genie_truth", "hnuc_pz"); + ntuple_.addVar("genie_truth", "hnuc_e"); + + ntuple_.addVar("genie_truth", "hqrk_pdg"); + ntuple_.addVar("genie_truth", "hqrk_sea"); + ntuple_.addVar("genie_truth", "hadsys_px"); + ntuple_.addVar("genie_truth", "hadsys_py"); + ntuple_.addVar("genie_truth", "hadsys_pz"); + ntuple_.addVar("genie_truth", "hadsys_e"); +} + +void GenieTruthDQM::onNewRun(const ldmx::RunHeader& runHeader) { + run_number_ = runHeader.getRunNumber(); +} + +void GenieTruthDQM::analyze(const framework::Event& event) { + getHistoDirectory(); + + ntuple_.clear(); + ntuple_.setVar("run", run_number_); + ntuple_.setVar("event", event.getEventNumber()); + + auto hepmc3_col = event.getObject >( + hepmc3CollName_, hepmc3PassName_); + + if (hepmc3_col.size() < 1) { + ntuple_.fill(); + return; + } + + auto const& hepmc3_ev = hepmc3_col.at(0).getHepMCGenEvent(); + + // set interaction/scattering codes + auto interaction_type_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.InteractionType"); + if (interaction_type_ptr) + ntuple_.setVar("interaction_type", interaction_type_ptr->value()); + auto scattering_type_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.ScatteringType"); + if (scattering_type_ptr) + ntuple_.setVar("scattering_type", scattering_type_ptr->value()); + auto rescatter_code_ptr = + hepmc3_ev.attribute("GENIE.RescatterCode"); + if (rescatter_code_ptr) + ntuple_.setVar("rescatter_code", scattering_type_ptr->value()); + + // get kinematic vars + auto kvar_labels_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.KineVarLabels"); + auto kvar_values_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.KineVarValues"); + if (kvar_labels_ptr && kvar_values_ptr) { + auto kvar_labels = kvar_labels_ptr->value(); + auto kvar_values = kvar_values_ptr->value(); + for (size_t i = 0; i < kvar_labels.size(); ++i) { + if (kvar_labels[i] == genie::EKineVar::kKVSelx) + ntuple_.setVar("x_bj", kvar_values[i]); + else if (kvar_labels[i] == genie::EKineVar::kKVSely) + ntuple_.setVar("y_inel", kvar_values[i]); + else if (kvar_labels[i] == genie::EKineVar::kKVSelQ2) + ntuple_.setVar("Q2", kvar_values[i]); + else if (kvar_labels[i] == genie::EKineVar::kKVSelW) + ntuple_.setVar("W", kvar_values[i]); + } + } + + // electron info + auto lep_pdg_ptr = + hepmc3_ev.attribute("GENIE.Interaction.ProbePDG"); + if (lep_pdg_ptr) ntuple_.setVar("lep_pdg", lep_pdg_ptr->value()); + + auto lep_i_4vec_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.ProbeP4"); + if (lep_i_4vec_ptr) { + auto lep_i_4vec = lep_i_4vec_ptr->value(); + ntuple_.setVar("lep_i_px", lep_i_4vec[0]); + ntuple_.setVar("lep_i_py", lep_i_4vec[1]); + ntuple_.setVar("lep_i_pz", lep_i_4vec[2]); + ntuple_.setVar("lep_i_e", lep_i_4vec[3]); + } + auto lep_f_4vec_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.FSLeptonP4"); + if (lep_f_4vec_ptr) { + auto lep_f_4vec = lep_f_4vec_ptr->value(); + ntuple_.setVar("lep_f_px", lep_f_4vec[0]); + ntuple_.setVar("lep_f_py", lep_f_4vec[1]); + ntuple_.setVar("lep_f_pz", lep_f_4vec[2]); + ntuple_.setVar("lep_f_e", lep_f_4vec[3]); + } + + // target info + auto tgt_pdg_ptr = + hepmc3_ev.attribute("GENIE.Interaction.TargetPDG"); + if (tgt_pdg_ptr) ntuple_.setVar("tgt_pdg", tgt_pdg_ptr->value()); + + auto tgt_4vec_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.TargetP4"); + if (tgt_4vec_ptr) { + auto tgt_4vec = tgt_4vec_ptr->value(); + ntuple_.setVar("tgt_px", tgt_4vec[0]); + ntuple_.setVar("tgt_py", tgt_4vec[1]); + ntuple_.setVar("tgt_pz", tgt_4vec[2]); + ntuple_.setVar("tgt_e", tgt_4vec[3]); + } + + // hit nucleon info + auto hnuc_pdg_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.HitNucleonPDG"); + if (hnuc_pdg_ptr) ntuple_.setVar("hnuc_pdg", hnuc_pdg_ptr->value()); + + auto hitnuc_4vec_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.HitNucleonP4"); + if (hitnuc_4vec_ptr) { + auto hitnuc_4vec = hitnuc_4vec_ptr->value(); + ntuple_.setVar("hnuc_px", hitnuc_4vec[0]); + ntuple_.setVar("hnuc_py", hitnuc_4vec[1]); + ntuple_.setVar("hnuc_pz", hitnuc_4vec[2]); + ntuple_.setVar("hnuc_e", hitnuc_4vec[3]); + } + // hit quark info + // note: it's only there for some interaction types! + auto hqrkpdg_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.HitQuarkPDG"); + if (hqrkpdg_ptr) ntuple_.setVar("hqrk_pdg", hqrkpdg_ptr->value()); + auto hqrksea_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.HitSeaQuark"); + if (hqrksea_ptr) ntuple_.setVar("hqrk_sea", hqrkpdg_ptr->value()); + + auto hadsys_4vec_ptr = hepmc3_ev.attribute( + "GENIE.Interaction.HadSystP4"); + if (hadsys_4vec_ptr) { + auto hadsys_4vec = hadsys_4vec_ptr->value(); + ntuple_.setVar("hadsys_px", hadsys_4vec[0]); + ntuple_.setVar("hadsys_py", hadsys_4vec[1]); + ntuple_.setVar("hadsys_pz", hadsys_4vec[2]); + ntuple_.setVar("hadsys_e", hadsys_4vec[3]); + } + ntuple_.fill(); + + return; +} + +} // namespace dqm + +DECLARE_ANALYZER_NS(dqm, GenieTruthDQM); diff --git a/DetDescr/python/EcalGeometry.py b/DetDescr/python/EcalGeometry.py index 589e434f1..ddddcb0ba 100644 --- a/DetDescr/python/EcalGeometry.py +++ b/DetDescr/python/EcalGeometry.py @@ -124,9 +124,8 @@ def v14() : eg.layer_shift_x = 2*eg.moduleMinR / eg.nCellRHeight return eg - def reduced() : - eg = EcalGeometry(detectors_valid = ["ldmx-reduced","ldmx-reduced[.].*"], + eg = EcalGeometry(detectors_valid = ["ldmx-reduced"], gap = 1.5, layerZPositions = [ 7.932, 14.532, 32.146, 40.746, 58.110, 67.710 @@ -139,5 +138,19 @@ def reduced() : eg.layer_shift_x = 2*eg.moduleMinR / eg.nCellRHeight return eg + def reduced_v2() : + eg = EcalGeometry(detectors_valid = ["ldmx-reduced-v2"], + gap = 1.5, + layerZPositions = [ + 9.932, 22.412, 40.81, 54.556, 71.954, 85.67 + ], + ecalFrontZ = 240.0, + cornersSideUp = True, + layer_shift_odd = True, + ) + # shift by a single cell diameter + eg.layer_shift_x = 2*eg.moduleMinR / eg.nCellRHeight + return eg + def geometries() : - return [EcalGeometry.v9(), EcalGeometry.v12(), EcalGeometry.v13(), EcalGeometry.v14(), EcalGeometry.reduced()] + return [EcalGeometry.v9(), EcalGeometry.v12(), EcalGeometry.v13(), EcalGeometry.v14(), EcalGeometry.reduced(), EcalGeometry.reduced_v2()] diff --git a/DetDescr/python/HcalGeometry.py b/DetDescr/python/HcalGeometry.py index 59aff9f7d..51c4b8c7a 100644 --- a/DetDescr/python/HcalGeometry.py +++ b/DetDescr/python/HcalGeometry.py @@ -487,5 +487,5 @@ def make_v14(self): ] # added the reduced geometry temporarily, for the final geometry # we should have a new function "reduced()" with the prototype geom - self.v14.detectors_valid = ["ldmx-det-v14", "ldmx-det-v14.*", "ldmx-reduced","ldmx-lyso-r1-v14", "ldmx-lyso-r1-v14.*"] + self.v14.detectors_valid = ["ldmx-det-v14", "ldmx-det-v14.*", "ldmx-reduced", "ldmx-reduced-v2","ldmx-lyso-r1-v14", "ldmx-lyso-r1-v14.*"] self.v14.y_offset = 19.05 diff --git a/Detectors/data/ldmx-reduced-v2/constants.gdml b/Detectors/data/ldmx-reduced-v2/constants.gdml new file mode 100644 index 000000000..ce4ced272 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/constants.gdml @@ -0,0 +1,471 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/detector.gdml b/Detectors/data/ldmx-reduced-v2/detector.gdml new file mode 100644 index 000000000..558557180 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/detector.gdml @@ -0,0 +1,215 @@ + + + +]> + + + + + &constants; + + + + + + + + + + + + + + + + + + + + + + + + + &materials; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/ecal.gdml b/Detectors/data/ldmx-reduced-v2/ecal.gdml new file mode 100644 index 000000000..93b9276c3 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/ecal.gdml @@ -0,0 +1,1297 @@ + + +]> + + + + + + + &constants; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + + + + + + + + + + + + + + + + + +
+
+ + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + + + + + + + + + + + + + + + + + +
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/ecal_motherboard5_assembly.gdml b/Detectors/data/ldmx-reduced-v2/ecal_motherboard5_assembly.gdml new file mode 100644 index 000000000..f06730b67 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/ecal_motherboard5_assembly.gdml @@ -0,0 +1,286 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/ecal_motherboard6_assembly.gdml b/Detectors/data/ldmx-reduced-v2/ecal_motherboard6_assembly.gdml new file mode 100644 index 000000000..edde53477 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/ecal_motherboard6_assembly.gdml @@ -0,0 +1,372 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/ecal_support_box.gdml b/Detectors/data/ldmx-reduced-v2/ecal_support_box.gdml new file mode 100644 index 000000000..c8a4d7c1f --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/ecal_support_box.gdml @@ -0,0 +1,1491 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/hcal.gdml b/Detectors/data/ldmx-reduced-v2/hcal.gdml new file mode 100644 index 000000000..a5e565aa9 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/hcal.gdml @@ -0,0 +1,448 @@ + + +]> + + + + + &constants; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/materials.gdml b/Detectors/data/ldmx-reduced-v2/materials.gdml new file mode 100644 index 000000000..2403d69f1 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/materials.gdml @@ -0,0 +1,45 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/recoil.gdml b/Detectors/data/ldmx-reduced-v2/recoil.gdml new file mode 100644 index 000000000..4afc3a90f --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/recoil.gdml @@ -0,0 +1,140 @@ + + + +]> + + + &constants; + + + + + + + + + + + + + + + + + + + + + + &materials; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/scoring_planes.gdml b/Detectors/data/ldmx-reduced-v2/scoring_planes.gdml new file mode 100644 index 000000000..3acac00b0 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/scoring_planes.gdml @@ -0,0 +1,363 @@ + + +]> + + + + + &constants; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/tagger.gdml b/Detectors/data/ldmx-reduced-v2/tagger.gdml new file mode 100644 index 000000000..1ae8f801b --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/tagger.gdml @@ -0,0 +1,89 @@ + + + + +]> + + + &constants; + + + + + + + + + + + + + + + + + + + + + &materials; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/target.gdml b/Detectors/data/ldmx-reduced-v2/target.gdml new file mode 100644 index 000000000..f6a86d6cb --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/target.gdml @@ -0,0 +1,112 @@ + + +]> + + + + + &constants; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detectors/data/ldmx-reduced-v2/trig_scint.gdml b/Detectors/data/ldmx-reduced-v2/trig_scint.gdml new file mode 100644 index 000000000..1c1baf4c6 --- /dev/null +++ b/Detectors/data/ldmx-reduced-v2/trig_scint.gdml @@ -0,0 +1,109 @@ + + +]> + + + + + &constants; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Ecal/python/digi.py b/Ecal/python/digi.py index c0b9bb72c..bc41dd461 100644 --- a/Ecal/python/digi.py +++ b/Ecal/python/digi.py @@ -211,3 +211,15 @@ def v14(self) : 10.915, 10.915, 14.783, 18.539, 18.539, 18.539, 18.539, 18.539, 18.539, 18.539, 18.539, 18.539, 9.938 ] + + def reduced_v2(self) : + """Generated for the reduced v2 geometry + + TODO: The secondOrderEnergyCorrection for this geometry has yet to be calculated, + so the v14 secondOrderEnergyCorrection is being used as a placeholder. + """ + + self.secondOrderEnergyCorrection = 4000. / 3940.5; + self.layerWeights = [ + 2.312, 5.417, 9.837, 11.910, 11.910, 11.910 + ] diff --git a/SimCore/CMakeLists.txt b/SimCore/CMakeLists.txt index c80435e43..6460a1489 100644 --- a/SimCore/CMakeLists.txt +++ b/SimCore/CMakeLists.txt @@ -15,20 +15,30 @@ find_package(ROOT CONFIG REQUIRED) find_package(Boost REQUIRED COMPONENTS log) +include_directories( /usr/local/include/GENIE ) + option(BUILD_EVENT_ONLY "Build the event library." ON) if(BUILD_EVENT_ONLY) + install(FILES /usr/local/lib/libHepMC3rootIO_rdict.pcm DESTINATION lib) + register_event_object(module_path "SimCore/Event" namespace "ldmx" class "SimCalorimeterHit" type "collection") register_event_object(module_path "SimCore/Event" namespace "ldmx" class "SimTrackerHit" type "collection") register_event_object(module_path "SimCore/Event" namespace "ldmx" class "SimParticle" type "map" key "int") + register_event_object(module_path "SimCore/Event" namespace "ldmx" + class "EventWeights") + + register_event_object(module_path "SimCore/Event" namespace "ldmx" + class "HepMC3GenEvent" type "collection") # Generate the files needed to build the event classes. setup_library(module SimCore name Event - dependencies ROOT::Core - register_target) + dependencies ROOT::Core ${GENIE_LIBS} HepMC3 HepMC3rootIO + register_target + include ) return() @@ -37,6 +47,9 @@ endif() # Configure Geant4 setup_geant4_target() +#Configure Genie +setup_genie_target() + # include dark brem simulation add_subdirectory(G4DarkBreM) @@ -64,12 +77,15 @@ file(GLOB SRC_FILES CONFIGURE_DEPDENDS ${PROJECT_SOURCE_DIR}/src/SimCore/Geo/[a- setup_library(module SimCore dependencies Geant4::Interface ROOT::Physics + ROOT::MathMore + HepMC3 Framework::Configure Framework::Framework SimCore::G4User G4DarkBreM DetDescr::DetDescr Boost::log + log4cpp "${registered_targets}" sources ${SRC_FILES}) @@ -81,7 +97,12 @@ setup_library(module SimCore # Primary Generators library setup_library(module SimCore name Generators - dependencies SimCore::SimCore SimCore::LHE) + dependencies SimCore::SimCore SimCore::LHE SimCore::Event ${GENIE_LIBS} HepMC3) + +# Reweight Library +setup_library(module SimCore + name Reweight + dependencies SimCore::SimCore SimCore::Generators SimCore::Event Framework::Framework ${GENIE_LIBS} HepMC3) # Sensitive Detectors library setup_library(module SimCore diff --git a/SimCore/include/SimCore/Event/EventWeights.h b/SimCore/include/SimCore/Event/EventWeights.h new file mode 100644 index 000000000..35d3f78a6 --- /dev/null +++ b/SimCore/include/SimCore/Event/EventWeights.h @@ -0,0 +1,151 @@ +// +// Created by Wesley Ketchum on 4/27/24. +// + +#ifndef SIM_CORE_EventWeights_H +#define SIM_CORE_EventWeights_H + +#include +#include + +#include "TObject.h" + +namespace ldmx { + +class EventWeights { + public: + // Enum to define + enum VariationType { + kINVALID = -1, + kUNKNOWN = 0, + + kGENIE_GENERIC = 1000, + kGENIE_INukeTwkDial = kGENIE_GENERIC + 100, + kGENIE_INukeTwkDial_MFP_pi = kGENIE_INukeTwkDial + 1, + kGENIE_INukeTwkDial_MFP_N = kGENIE_INukeTwkDial + 2, + kGENIE_INukeTwkDial_FrCEx_pi = kGENIE_INukeTwkDial + 3, + kGENIE_INukeTwkDial_FrInel_pi = kGENIE_INukeTwkDial + 4, + kGENIE_INukeTwkDial_FrAbs_pi = kGENIE_INukeTwkDial + 5, + kGENIE_INukeTwkDial_FrPiProd_pi = kGENIE_INukeTwkDial + 6, + kGENIE_INukeTwkDial_FrCEx_N = kGENIE_INukeTwkDial + 7, + kGENIE_INukeTwkDial_FrInel_N = kGENIE_INukeTwkDial + 8, + kGENIE_INukeTwkDial_FrAbs_N = kGENIE_INukeTwkDial + 9, + kGENIE_INukeTwkDial_FrPiProd_N = kGENIE_INukeTwkDial + 10, + + kGENIE_HadrNuclTwkDial = kGENIE_GENERIC + 150, + kGENIE_HadrNuclTwkDial_FormZone = kGENIE_HadrNuclTwkDial + 1 + }; + + EventWeights() = default; + virtual ~EventWeights() = default; + + // constructor where variation map is declared + EventWeights(std::map > var_map) + : variations_map_(var_map) {} + + void Clear(); + void Print(); + + std::vector getWeights() const { return weights_; } + std::map > getVariations() const { + return variations_map_; + } + + size_t getNumWeights() const { return weights_.size(); } + size_t getNumVariationTypes() const { return variations_map_.size(); } + + double getNthWeight(size_t i_w) const { return weights_.at(i_w); } + + std::map getVariationsNthWeight(size_t i_w) const; + + void addWeight(double w) { weights_.push_back(w); } + void setWeight(size_t i_w, double w) { weights_[i_w] = w; } + + private: + // set of event weights (n weights per event) + std::vector weights_; + + // map of the variations used: key is variation type, value is list of + // variation values (n per event) + std::map > variations_map_; + + public: + inline static std::string variation_type_to_string( + const VariationType& vtype) { + switch (vtype) { + case VariationType::kINVALID: + return "INVALID"; + case VariationType::kUNKNOWN: + return "UNKNOWN"; + case VariationType::kGENIE_GENERIC: + return "GENIE_GENERIC"; + case VariationType::kGENIE_INukeTwkDial: + return "GENIE_INukeTwkDial"; + case VariationType::kGENIE_INukeTwkDial_MFP_pi: + return "GENIE_INukeTwkDial_MFP_pi"; + case VariationType::kGENIE_INukeTwkDial_MFP_N: + return "GENIE_INukeTwkDial_MFP_N"; + case VariationType::kGENIE_INukeTwkDial_FrCEx_pi: + return "GENIE_INukeTwkDial_FrCEx_pi"; + case VariationType::kGENIE_INukeTwkDial_FrInel_pi: + return "GENIE_INukeTwkDial_FrInel_pi"; + case VariationType::kGENIE_INukeTwkDial_FrAbs_pi: + return "GENIE_INukeTwkDial_FrAbs_pi"; + case VariationType::kGENIE_INukeTwkDial_FrPiProd_pi: + return "GENIE_INukeTwkDial_FrPiProd_pi"; + case VariationType::kGENIE_INukeTwkDial_FrCEx_N: + return "GENIE_INukeTwkDial_FrCEx_N"; + case VariationType::kGENIE_INukeTwkDial_FrInel_N: + return "GENIE_INukeTwkDial_FrInel_N"; + case VariationType::kGENIE_INukeTwkDial_FrAbs_N: + return "GENIE_INukeTwkDial_FrAbs_N"; + case VariationType::kGENIE_INukeTwkDial_FrPiProd_N: + return "GENIE_INukeTwkDial_FrPiProd_N"; + case VariationType ::kGENIE_HadrNuclTwkDial: + return "GENIE_HadrNuclTwkDial"; + case VariationType ::kGENIE_HadrNuclTwkDial_FormZone: + return "GENIE_HadrNuclTwkDial_FormZone"; + } + return "UNKNOWN"; + } + inline static VariationType string_to_variation_type( + const std::string& typestring) { + if (typestring == "UNKNOWN") + return VariationType::kUNKNOWN; + else if (typestring == "INVALID") + return VariationType::kINVALID; + else if (typestring == "GENIE_GENERIC") + return VariationType::kGENIE_GENERIC; + else if (typestring == "GENIE_INukeTwkDial") + return VariationType::kGENIE_INukeTwkDial; + else if (typestring == "GENIE_INukeTwkDial_MFP_pi") + return VariationType::kGENIE_INukeTwkDial_MFP_pi; + else if (typestring == "GENIE_INukeTwkDial_MFP_N") + return VariationType::kGENIE_INukeTwkDial_MFP_N; + else if (typestring == "GENIE_INukeTwkDial_FrCEx_pi") + return VariationType::kGENIE_INukeTwkDial_FrCEx_pi; + else if (typestring == "GENIE_INukeTwkDial_FrInel_pi") + return VariationType::kGENIE_INukeTwkDial_FrInel_pi; + else if (typestring == "GENIE_INukeTwkDial_FrAbs_pi") + return VariationType::kGENIE_INukeTwkDial_FrAbs_pi; + else if (typestring == "GENIE_INukeTwkDial_FrPiProd_pi") + return VariationType::kGENIE_INukeTwkDial_FrPiProd_pi; + else if (typestring == "GENIE_INukeTwkDial_FrCEx_N") + return VariationType::kGENIE_INukeTwkDial_FrCEx_N; + else if (typestring == "GENIE_INukeTwkDial_FrInel_N") + return VariationType::kGENIE_INukeTwkDial_FrInel_N; + else if (typestring == "GENIE_INukeTwkDial_FrAbs_N") + return VariationType::kGENIE_INukeTwkDial_FrAbs_N; + else if (typestring == "GENIE_INukeTwkDial_FrPiProd_N") + return VariationType::kGENIE_INukeTwkDial_FrPiProd_N; + else if (typestring == "GENIE_HadrNuclTwkDial") + return VariationType::kGENIE_HadrNuclTwkDial; + else if (typestring == "GENIE_HadrNuclTwkDial_FormZone") + return VariationType::kGENIE_HadrNuclTwkDial_FormZone; + + return VariationType::kUNKNOWN; + } +}; +} // namespace ldmx + +#endif // SIM_CORE_EventWeights_H diff --git a/SimCore/include/SimCore/Event/HepMC3GenEvent.h b/SimCore/include/SimCore/Event/HepMC3GenEvent.h new file mode 100644 index 000000000..ef65a22d6 --- /dev/null +++ b/SimCore/include/SimCore/Event/HepMC3GenEvent.h @@ -0,0 +1,35 @@ +// +// Created by Wesley Ketchum on 4/26/24. +// +// This is just a simple extension of the HepMC3::GenEvent +// for use on the ldmx Event Bus. +// + +#ifndef SIM_CORE_HEPMC3GENEVENT_H +#define SIM_CORE_HEPMC3GENEVENT_H + +#include "HepMC3/Data/GenEventData.h" +#include "HepMC3/GenEvent.h" +#include "TObject.h" + +namespace ldmx { + +class HepMC3GenEvent : public HepMC3::GenEventData { + public: + // ~HepMC3GenEvent() {} + // HepMC3GenEvent(const HepMC3::GenEvent& event) : HepMC3::GenEvent(event) {} + + void Clear(); + void Print() const; + + HepMC3::GenEvent getHepMCGenEvent() const; + + std::string get_as_string() const; + + public: + ClassDef(HepMC3GenEvent, 1); +}; + +} // namespace ldmx + +#endif // SIM_CORE_HEPMC3GENEVENT_H diff --git a/SimCore/include/SimCore/Generators/GenieGenerator.h b/SimCore/include/SimCore/Generators/GenieGenerator.h new file mode 100644 index 000000000..5c933856b --- /dev/null +++ b/SimCore/include/SimCore/Generators/GenieGenerator.h @@ -0,0 +1,114 @@ +/** + * @file GenieGenerator.h + * @brief Simple GENIE event generator. + * @author Wesley Ketchum, FNAL + */ + +#ifndef SIMCORE_GENIE_GENERATOR_H +#define SIMCORE_GENIE_GENERATOR_H + +//----------// +// ROOT // +//----------// +#include "TRandom.h" +#include "TRandomGen.h" + +//------------// +// GENIE // +//------------// +#include "GENIE/Framework/EventGen/GEVGDriver.h" +#include "GENIE/Framework/EventGen/HepMC3Converter.h" + +//------------// +// LDMX // +//------------// +#include +#include + +#include "SimCore/PrimaryGenerator.h" + +// Forward declarations +class G4Event; + +namespace simcore { +namespace generators { + +/** + * @class GenieGenerator + * @brief Class that uses GENIE's GEVGDriver to generator eN interactions. + */ +class GenieGenerator : public simcore::PrimaryGenerator { + public: + /** + * Constructor. + * + * @param parameters Parameters used to configure GENIE generator. + * + * Parameters: + * verbosity : > 1 means print configuration + * energy : energy of initial electron (GeV) + * targets : list of ten-digit 10LZZZAAAI target codes + * abundances: list of relative abundances for the given targets + * position : position of interaction from (mm three-vector) + * time : time to shoot at (ns) + * direction : direction to shoot in (unitless three-vector) + * tune : name of GENIE tune + * seed : seed for random generator + */ + GenieGenerator(const std::string& name, + const framework::config::Parameters& parameters); + + /// Destructor + ~GenieGenerator(); + + /** + * Generate the primary vertices in the Geant4 event. + * + * @param event The Geant4 event. + */ + void GeneratePrimaryVertex(G4Event* event) final override; + + void RecordConfig(const std::string& id, ldmx::RunHeader& rh) final override; + + private: + /** + * The GENIE event generator driver and convertor to HepMC3GenEvent + */ + genie::GEVGDriver evg_driver_; + genie::HepMC3Converter hepMC3Converter_; + + int verbosity_; + double energy_; + std::vector targets_; + std::vector abundances_; + std::vector position_; + std::vector beam_size_; + double target_thickness_; + double time_; + std::vector direction_; + + std::string tune_; + std::string spline_file_; + + std::string message_threshold_file_; + + std::vector ev_weighting_integral_; + size_t n_events_generated_; + std::vector n_events_by_target_; + std::vector xsec_by_target_; + + double xsec_total_; + + void fillConfig( + const framework::config::Parameters&); /// fill the configuration + bool validateConfig(); /// simple validation check on configuration params + + void initializeGENIE(); /// GENIE initialization + void calculateTotalXS(); /// GENIE initialization + +}; // ParticleGun + +} // namespace generators +} // namespace simcore + +#endif // SIMCORE_PARTICLE_GUN_H diff --git a/SimCore/include/SimCore/PrimaryGenerator.h b/SimCore/include/SimCore/PrimaryGenerator.h index b6e99cf2f..536ab3c9a 100644 --- a/SimCore/include/SimCore/PrimaryGenerator.h +++ b/SimCore/include/SimCore/PrimaryGenerator.h @@ -22,6 +22,7 @@ #include "Framework/Configure/Parameters.h" #include "Framework/RunHeader.h" #include "SimCore/Factory.h" +#include "UserEventInformation.h" // Forward Declarations class G4Event; @@ -69,6 +70,8 @@ class PrimaryGenerator : public G4VPrimaryGenerator { */ virtual void RecordConfig(const std::string& id, ldmx::RunHeader& rh) = 0; + std::string Name() { return name_; } + protected: /// Name of the PrimaryGenerator std::string name_{""}; diff --git a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h new file mode 100644 index 000000000..57b863003 --- /dev/null +++ b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h @@ -0,0 +1,110 @@ +// +// Created by Wesley Ketchum on 4/29/24. +// + +#ifndef SIMCORE_GENIEREWEIGHTPRODUCER_H +#define SIMCORE_GENIEREWEIGHTPRODUCER_H + +#include +#include + +#include "Framework/EventProcessor.h" +#include "RwFramework/GSyst.h" +#include "SimCore/Event/EventWeights.h" + +namespace genie { +class Interaction; +class HepMC3Converter; +} // namespace genie + +namespace genie::rew { +class GReWeight; +} + +namespace simcore { + +class GenieReweightProducer : public framework::Producer { + public: + // constructor + GenieReweightProducer(const std::string& name, framework::Process& process); + + // default destructor + virtual ~GenieReweightProducer(); + + // configuration + virtual void configure(framework::config::Parameters&); + + // on new run + virtual void onNewRun(const ldmx::RunHeader& runHeader); + + // produce on the event + virtual void produce(framework::Event& event); + + private: + // seed to use + int verbosity_; + + // input hepmc3 collection name + std::string hepmc3CollName_; + + // input hepmc3 pass name + std::string hepmc3PassName_; + + // output EventWeights collection name + std::string eventWeightsCollName_; + + // seed to use + int seed_; + + // number of weights to be calculated per event + size_t n_weights_; + + // GENIE tune + std::string tune_; + + // variations to run + std::map > + variation_map_; + + // hepmc3 convertor + genie::HepMC3Converter* hepMC3Converter_; + + genie::rew::GReWeight* genie_rw_; + + void reinitializeGenieReweight(); + void reconfigureGenieReweight(size_t); + + inline static genie::rew::EGSyst variation_type_to_genie_dial( + const ldmx::EventWeights::VariationType& vtype) { + switch (vtype) { + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_MFP_pi: + return genie::rew::EGSyst::kINukeTwkDial_MFP_pi; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_MFP_N: + return genie::rew::EGSyst::kINukeTwkDial_MFP_N; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrCEx_pi: + return genie::rew::EGSyst::kINukeTwkDial_FrCEx_pi; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrInel_pi: + return genie::rew::EGSyst::kINukeTwkDial_FrInel_pi; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrAbs_pi: + return genie::rew::EGSyst::kINukeTwkDial_FrAbs_pi; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrPiProd_pi: + return genie::rew::EGSyst::kINukeTwkDial_FrPiProd_pi; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrCEx_N: + return genie::rew::EGSyst::kINukeTwkDial_FrCEx_N; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrInel_N: + return genie::rew::EGSyst::kINukeTwkDial_FrInel_N; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrAbs_N: + return genie::rew::EGSyst::kINukeTwkDial_FrAbs_N; + case ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrPiProd_N: + return genie::rew::EGSyst::kINukeTwkDial_FrPiProd_N; + case ldmx::EventWeights::VariationType ::kGENIE_HadrNuclTwkDial_FormZone: + return genie::rew::EGSyst::kHadrNuclTwkDial_FormZone; + default: + return genie::rew::EGSyst::kNullSystematic; + } + return genie::rew::EGSyst::kNullSystematic; + } +}; +} // namespace simcore + +#endif // SIMCORE_GENIEREWEIGHTPRODUCER_H diff --git a/SimCore/include/SimCore/Simulator.h b/SimCore/include/SimCore/Simulator.h index 21fc0071a..0a0b3d23c 100644 --- a/SimCore/include/SimCore/Simulator.h +++ b/SimCore/include/SimCore/Simulator.h @@ -26,6 +26,14 @@ #include "SimCore/RunManager.h" #include "SimCore/SimulatorBase.h" +namespace genie { +class Interaction; +} + +namespace ldmx { +class HepMC3GenEvent; +} + class G4UImanager; class G4UIsession; class G4RunManager; diff --git a/SimCore/include/SimCore/UserEventInformation.h b/SimCore/include/SimCore/UserEventInformation.h index 162ab33ed..7fbe28f99 100644 --- a/SimCore/include/SimCore/UserEventInformation.h +++ b/SimCore/include/SimCore/UserEventInformation.h @@ -1,7 +1,11 @@ #ifndef SIMCORE_USEREVENTINFORMATION_H #define SIMCORE_USEREVENTINFORMATION_H +#include + #include "G4VUserEventInformation.hh" +#include "SimCore/Event/HepMC3GenEvent.h" + namespace simcore { /** @@ -13,7 +17,7 @@ class UserEventInformation : public G4VUserEventInformation { UserEventInformation() = default; /// Destructor - virtual ~UserEventInformation() = default; + virtual ~UserEventInformation() {} /// Print the information associated with the track void Print() const override; @@ -117,6 +121,13 @@ class UserEventInformation : public G4VUserEventInformation { */ bool wasLastStepEN() const { return last_step_en_; } + void addHepMC3GenEvent(ldmx::HepMC3GenEvent event) { + hepmc3_events_.push_back(event); + } + std::vector getHepMC3GenEvents() { + return hepmc3_events_; + } + private: /// Total number of brem candidates in the event int bremCandidateCount_{0}; @@ -165,6 +176,11 @@ class UserEventInformation : public G4VUserEventInformation { * dark brem did not occur within the event in question. */ double db_material_z_{-1.}; + + /** + * a collection of HepMC3 event records. + */ + std::vector hepmc3_events_; }; } // namespace simcore diff --git a/SimCore/python/generators.py b/SimCore/python/generators.py index 95fe09b1f..04397143d 100644 --- a/SimCore/python/generators.py +++ b/SimCore/python/generators.py @@ -193,6 +193,78 @@ def __init__(self,name,initCommands) : super().__init__(name,'simcore::generators::GeneralParticleSource') self.initCommands = initCommands +class genie(simcfg.PrimaryGenerator) : + """Simple GENIE generator + + Parameters + ---------- + name : str + name of new primary generator + + Attributes + ---------- + verbosity : int, optional + Verbosity flag for this generator + energy : float + Energy of particle to shoot [GeV] + targets : list of int + List of target nuclei PDG codes + abundances : list of float + List of abundances for target nuclei + time : double + Time to shoot from [ns] + position : list of double + Position of interaction from [mm] + direction : list of double + Unit vector direction of incoming electron + tune: str + Name of GENIE tune to use + target_thickness : double + Thickness of target to use for generation + beam_size : list of double + uniform beam size width to use + + Examples + -------- + myGenie = genie( name='myGenie', + energy = 4.0, + targets = [ 1000220480 ], + abundances = [ 1.0 ], + time = 0.0, + position = [0.,0.,0.], + direction = [0.,0.,1.], + tune='G18_02a_00_000', + verbosity=True) + """ + + def __init__(self,name, + energy=4.0, + targets = [], + target_thickness = 0.3504, + abundances = [], + time = 0.0, + position = [ 0.0, 0.0, 0.0 ], + beam_size = [ 0.0, 0.0 ], + direction = [ 0.0, 0.0, 1.0 ], + tune = 'default', + spline_file = '', + message_threshold_file = "/usr/local/GENIE/Generator/config/Messenger.xml", + verbosity=0 ) : + super().__init__( name , "simcore::generators::GenieGenerator" ) + + self.energy = energy + self.targets = targets + self.target_thickness = target_thickness + self.abundances = abundances + self.time = time + self.position = position + self.beam_size = beam_size + self.direction = direction + self.tune = tune + self.spline_file = spline_file + self.message_threshold_file = message_threshold_file + self.verbosity = verbosity + def _single_e_upstream_tagger(position, momentum, energy): """Internal helper function for creating electron beam guns upstream of tagger diff --git a/SimCore/python/genie_reweight.py b/SimCore/python/genie_reweight.py new file mode 100644 index 000000000..9703c8358 --- /dev/null +++ b/SimCore/python/genie_reweight.py @@ -0,0 +1,13 @@ +from LDMX.Framework import ldmxcfg + +class GenieReweightProducer(ldmxcfg.Producer) : + + def __init__(self,name='genieEventWeights'): + super().__init__(name,"simcore::GenieReweightProducer","SimCore::Reweight") + + self.verbosity = 0 + self.seed = 10 + self.n_weights = 100 + self.var_types = ["GENIE_GENERIC"] + + self.eventWeightsCollName = "genieEventWeights" \ No newline at end of file diff --git a/SimCore/src/SimCore/Event/EventWeights.cxx b/SimCore/src/SimCore/Event/EventWeights.cxx new file mode 100644 index 000000000..9c524142c --- /dev/null +++ b/SimCore/src/SimCore/Event/EventWeights.cxx @@ -0,0 +1,27 @@ +// +// Created by Wesley Ketchum on 4/28/24. +// + +#include "SimCore/Event/EventWeights.h" + +#include + +namespace ldmx { + +void EventWeights::Clear() { + weights_.clear(); + variations_map_.clear(); +} + +void EventWeights::Print() { + std::cout << "Num weights: " << getNumWeights() << std::endl; + for (size_t i_w = 0; i_w < weights_.size(); ++i_w) { + std::cout << "\t" << i_w << ": weight=" << weights_[i_w]; + for (auto const& v : variations_map_) { + std::cout << "\n\t var_type=" << v.first + << ", var_value=" << v.second[i_w]; + } + std::cout << std::endl; + } +} +} // namespace ldmx \ No newline at end of file diff --git a/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx b/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx new file mode 100644 index 000000000..8e1a77aae --- /dev/null +++ b/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx @@ -0,0 +1,47 @@ +// +// Created by Wesley Ketchum on 4/26/24. +// + +#include "SimCore/Event/HepMC3GenEvent.h" + +#include +#include + +#include "HepMC3/Print.h" +#include "HepMC3/WriterAscii.h" + +namespace ldmx { + +void HepMC3GenEvent::Clear() { + this->particles.clear(); + this->vertices.clear(); + this->links1.clear(); + this->links2.clear(); + this->attribute_id.clear(); + this->attribute_name.clear(); + this->attribute_string.clear(); +} + +void HepMC3GenEvent::Print() const { + HepMC3::GenEvent ev; + ev.read_data(*this); + HepMC3::Print::line(ev, true); // print attributes +} + +HepMC3::GenEvent HepMC3GenEvent::getHepMCGenEvent() const { + HepMC3::GenEvent ev; + ev.read_data(*this); + return ev; +} + +std::string HepMC3GenEvent::get_as_string() const { + HepMC3::GenEvent ev; + ev.read_data(*this); + + std::stringstream ss; + HepMC3::WriterAscii writer(ss); + + writer.write_event(ev); + return ss.str(); +} +} // namespace ldmx diff --git a/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx b/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx index bc925a265..a60f168bd 100644 --- a/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx +++ b/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx @@ -70,7 +70,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event) { } // Make our information container and give it to geant4 - // G4Event owns the event information and will delete it + // G4Event owns the event information and will delete it auto event_info = new UserEventInformation; event->SetUserInformation(event_info); diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx new file mode 100644 index 000000000..dc550616f --- /dev/null +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -0,0 +1,413 @@ +/** + * @file GenieGenerator.cxx + * @brief Simple GENIE event generator. + * @author Wesley Ketchum, FNAL + */ + +#include "SimCore/Generators/GenieGenerator.h" + +#include "SimCore/UserPrimaryParticleInformation.h" + +// GENIE +#include "Framework/Conventions/Controls.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Conventions/XmlParserStatus.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/EventGen/GEVGDriver.h" +#include "Framework/EventGen/GFluxI.h" +#include "Framework/EventGen/GMCJDriver.h" +#include "Framework/EventGen/GMCJMonitor.h" +#include "Framework/EventGen/InteractionList.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Ntuple/NtpMCFormat.h" +#include "Framework/Ntuple/NtpWriter.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/Numerical/Spline.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/Utils/AppInit.h" +#include "Framework/Utils/CmdLnArgParser.h" +#include "Framework/Utils/PrintUtils.h" +#include "Framework/Utils/RunOpt.h" +#include "Framework/Utils/StringUtils.h" +#include "Framework/Utils/SystemUtils.h" +#include "Framework/Utils/XSecSplineList.h" +#include "GENIE/Framework/Interaction/InitialState.h" +#include "GENIE/Framework/Utils/RunOpt.h" +#include "HepMC3/GenEvent.h" + +// Geant4 +#include "G4Event.hh" +#include "G4PhysicalConstants.hh" +#include "Randomize.hh" + +// ROOT +#include +#include +#include + +#include "Math/Interpolator.h" + +// standard +#include +#include +#include + +namespace simcore { +namespace generators { + +void GenieGenerator::fillConfig(const framework::config::Parameters& p) { + energy_ = p.getParameter("energy"); // * GeV; + + targets_ = p.getParameter >("targets"); + abundances_ = p.getParameter >("abundances"); + + time_ = p.getParameter("time"); // * ns; + position_ = p.getParameter >("position"); // mm + beam_size_ = p.getParameter >("beam_size"); // mm + direction_ = p.getParameter >("direction"); + target_thickness_ = p.getParameter("target_thickness"); // mm + + tune_ = p.getParameter("tune"); + spline_file_ = p.getParameter("spline_file"); + + message_threshold_file_ = + p.getParameter("message_threshold_file"); + + verbosity_ = p.getParameter("verbosity"); +} + +bool GenieGenerator::validateConfig() { + bool ret = true; + + if (targets_.size() == 0 || abundances_.size() == 0) { + std::cout << "targets and/or abundances sizes are zero." + << " " << targets_.size() << ", " << abundances_.size() + << std::endl; + ret = false; + } + if (targets_.size() != abundances_.size()) { + std::cout << "targets and abundances sizes unequal." + << " " << targets_.size() << " != " << abundances_.size() + << std::endl; + ret = false; + } + + if (position_.size() != 3 || direction_.size() != 3) { + std::cout << "position and/or direction sizes are not 3." + << " " << position_.size() << ", " << direction_.size() + << std::endl; + ret = false; + } + + if (target_thickness_ < 0) { + std::cout << "target thickness cannot be less than 0. " << target_thickness_ + << std::endl; + std::cout << "Taking absolute value." << std::endl; + target_thickness_ = std::abs(target_thickness_); + } + + if (beam_size_.size() != 2) { + if (beam_size_.size() == 0) { + std::cout << "beam size not set. Using zero." << std::endl; + beam_size_.resize(2); + beam_size_[0] = 0.0; + beam_size_[1] = 0.0; + } else { + std::cout << "beam size is set, but does not have size 2." + << " " << beam_size_.size() << std::endl; + ret = false; + } + } else if (beam_size_[0] < 0 || beam_size_[1] < 0) { + std::cout << "Beam size set as negative value? " + << "(" << beam_size_[0] << "," << beam_size_[1] << ")" + << std::endl; + std::cout << "Changing to positive." << std::endl; + beam_size_[0] = std::abs(beam_size_[0]); + beam_size_[1] = std::abs(beam_size_[1]); + } + + // normalize abundances + double abundance_sum = 0; + for (auto a : abundances_) { + abundance_sum += a; + } + + if (std::abs(abundance_sum) < 1e-6) { + std::cout << "abundances list sums to zero? " << abundance_sum << std::endl; + ret = false; + } + + if (std::abs(abundance_sum - 1.0) > 2e-2) { + std::cout << "abundances list sums is not unity (" << abundance_sum + << " instead.)" << std::endl; + std::cout << "Will renormalize abundances to unity!" << std::endl; + } + + for (size_t i_a = 0; i_a < abundances_.size(); ++i_a) { + abundances_[i_a] = abundances_[i_a] / abundance_sum; + + if (verbosity_ > 0) + std::cout << "Target=" << targets_[i_a] + << ", Abundance=" << abundances_[i_a] << std::endl; + } + + double dir_total_sq = 0; + for (auto d : direction_) dir_total_sq += d * d; + + if (dir_total_sq < 1e-6) { + std::cout << "direction vector is zero or negative? " + << "(" << direction_[0] << "," << direction_[1] << "," + << direction_[2] << ")" << std::endl; + ret = false; + } + for (size_t i_d = 0; i_d < direction_.size(); ++i_d) + direction_[i_d] = direction_[i_d] / std::sqrt(dir_total_sq); + + xsec_by_target_.resize(targets_.size(), -999.); + n_events_by_target_.resize(targets_.size(), 0); + + return ret; +} + +void GenieGenerator::initializeGENIE() { + // initialize some RunOpt by hacking the command line interface + { + char* inarr[3] = {const_cast(""), + const_cast("--event-generator-list"), + const_cast("EM")}; + genie::RunOpt::Instance()->ReadFromCommandLine(3, inarr); + } + + // set message thresholds + genie::utils::app_init::MesgThresholds(message_threshold_file_); + + // set tune info + genie::RunOpt::Instance()->SetTuneName(tune_); + if (!genie::RunOpt::Instance()->Tune()) { + EXCEPTION_RAISE("ConfigurationException", "No TuneId in RunOption."); + } + genie::RunOpt::Instance()->BuildTune(); + + /* + //set random seed + auto seed = G4Random::getTheEngine()->getSeed(); + if(verbosity_>=1) + std::cout << "Initializing GENIE with seed " << seed << std::endl; + genie::utils::app_init::RandGen(seed); + */ + + // give it the splint file and require it + genie::utils::app_init::XSecTable(spline_file_, true); + + // set GHEP print level (needed?) + genie::GHepRecord::SetPrintLevel(0); + + // setup for event driver + evg_driver_.SetEventGeneratorList( + genie::RunOpt::Instance()->EventGeneratorList()); + evg_driver_.SetUnphysEventMask(*genie::RunOpt::Instance()->UnphysEventMask()); +} + +void GenieGenerator::calculateTotalXS() { + // initializing... + xsec_total_ = 0; + ev_weighting_integral_.resize(targets_.size(), 0.0); + + // calculate the total xsec per target... + for (size_t i_t = 0; i_t < targets_.size(); ++i_t) { + genie::InitialState initial_state(targets_[i_t], 11); + evg_driver_.Configure(initial_state); + evg_driver_.UseSplines(); + + // setup the initial election + TParticle initial_e; + initial_e.SetPdgCode(11); + auto elec_i_p = std::sqrt(energy_ * energy_ - + initial_e.GetMass() * initial_e.GetMass()); + initial_e.SetMomentum(elec_i_p * direction_[0], elec_i_p * direction_[1], + elec_i_p * direction_[2], energy_); + TLorentzVector e_p4; + initial_e.Momentum(e_p4); + + xsec_by_target_[i_t] = evg_driver_.XSecSum(e_p4); + xsec_total_ += xsec_by_target_[i_t] * abundances_[i_t]; + + ev_weighting_integral_[i_t] = xsec_total_; // running sum + + // print... + std::cout << "Target=" << targets_[i_t] + << "\tAbundance=" << abundances_[i_t] + << "\tXSEC=" << xsec_by_target_[i_t] / genie::units::millibarn + << " mb" << std::endl; + } + std::cout << "Total XSEC = " << xsec_total_ / genie::units::millibarn << " mb" + << std::endl; + + // renormalize our weighting integral + for (size_t i_t = 0; i_t < ev_weighting_integral_.size(); ++i_t) + ev_weighting_integral_[i_t] = ev_weighting_integral_[i_t] / xsec_total_; +} + +GenieGenerator::GenieGenerator(const std::string& name, + const framework::config::Parameters& p) + : PrimaryGenerator(name, p) { + fillConfig(p); + + if (!validateConfig()) + EXCEPTION_RAISE("ConfigurationException", "Configuration not valid."); + + n_events_generated_ = 0; + initializeGENIE(); + calculateTotalXS(); +} + +GenieGenerator::~GenieGenerator() { + std::cout << "--- GENIE Generation Summary BEGIN ---" << std::endl; + double total_xsec = 0; + for (size_t i_t = 0; i_t < targets_.size(); ++i_t) { + std::cout << "Target=" << targets_[i_t] + << "\tAbundance=" << abundances_[i_t] + << "\tXSEC=" << xsec_by_target_[i_t] / genie::units::millibarn + << " mb" + << "\tEvents=" << n_events_by_target_[i_t] << std::endl; + if (n_events_by_target_[i_t] > 0) + total_xsec += xsec_by_target_[i_t] * abundances_[i_t]; + } + + std::cout << "Total events generated = " << n_events_generated_ << std::endl; + std::cout << "Total XSEC = " << total_xsec / genie::units::millibarn << " mb" + << std::endl; + + std::cout << "--- GENIE Generation Summary *END* ---" << std::endl; +} + +void GenieGenerator::GeneratePrimaryVertex(G4Event* event) { + if (n_events_generated_ == 0) { + // set random seed + // have to do this here since seeds aren't properly set until we know the + // run number + auto seed = G4Random::getTheEngine()->getSeed(); + if (verbosity_ >= 1) + std::cout << "Initializing GENIE with seed " << seed << std::endl; + genie::utils::app_init::RandGen(seed); + } + + auto nucl_target_i = 0; + + if (targets_.size() > 0) { + double rand_uniform = G4Random::getTheGenerator()->flat(); + + nucl_target_i = std::distance( + ev_weighting_integral_.begin(), + std::lower_bound(ev_weighting_integral_.begin(), + ev_weighting_integral_.end(), rand_uniform)); + if (verbosity_ >= 1) + std::cout << "Random number = " << rand_uniform << ", target picked " + << targets_.at(nucl_target_i) << std::endl; + } + + auto x_pos = position_[0] + + (G4Random::getTheGenerator()->flat() - 0.5) * beam_size_[0]; + auto y_pos = position_[1] + + (G4Random::getTheGenerator()->flat() - 0.5) * beam_size_[1]; + auto z_pos = position_[2] + + (G4Random::getTheGenerator()->flat() - 0.5) * target_thickness_; + if (verbosity_ >= 1) + std::cout << "Generating interaction at (x,y,z)=" + << "(" << x_pos << "," << y_pos << "," << z_pos << ")" + << std::endl; + + genie::InitialState initial_state(targets_.at(nucl_target_i), 11); + evg_driver_.Configure(initial_state); + evg_driver_.UseSplines(); + + // setup the initial election + TParticle initial_e; + initial_e.SetPdgCode(11); + double elec_i_p = + std::sqrt(energy_ * energy_ - initial_e.GetMass() * initial_e.GetMass()); + initial_e.SetMomentum(elec_i_p * direction_[0], elec_i_p * direction_[1], + elec_i_p * direction_[2], energy_); + TLorentzVector e_p4; + initial_e.Momentum(e_p4); + if (verbosity_ >= 1) + std::cout << "Generating interation with (px,py,pz,e)=" + << "(" << e_p4.Px() << "," << e_p4.Py() << "," << e_p4.Pz() << "," + << e_p4.E() << ")" << std::endl; + + // calculate total xsec + if (n_events_by_target_[nucl_target_i] == 0) + xsec_by_target_[nucl_target_i] = evg_driver_.XSecSum(e_p4); + + n_events_by_target_[nucl_target_i] += 1; + + // GENIE magic + genie::EventRecord* genie_event = NULL; + while (!genie_event) genie_event = evg_driver_.GenerateEvent(e_p4); + + auto ev_info = new UserEventInformation; + auto hepmc3_genie = hepMC3Converter_.ConvertToHepMC3(*genie_event); + ldmx::HepMC3GenEvent hepmc3_ldmx_genie; + hepmc3_genie->write_data(hepmc3_ldmx_genie); + ev_info->addHepMC3GenEvent(hepmc3_ldmx_genie); + event->SetUserInformation(ev_info); + + // setup the primary vertex now + + G4PrimaryVertex* vertex = new G4PrimaryVertex(); + vertex->SetPosition(x_pos, y_pos, z_pos); + vertex->SetWeight(genie_event->Weight()); + + // loop over the entries and add to the G4Event + int nentries = genie_event->GetEntries(); + + if (verbosity_ >= 1) { + std::cout << "---------- " + << "Generated Event " << n_events_generated_ + 1 << " ----------" + << std::endl; + } + + for (int i_p = 0; i_p < nentries; ++i_p) { + genie::GHepParticle* p = (genie::GHepParticle*)(*genie_event)[i_p]; + + // make sure it's a final state particle + if (p->Status() != 1) continue; + + if (verbosity_ >= 1) + std::cout << "\tAdding particle " << p->Pdg() << " with status " + << p->Status() << " energy " << p->E() << " ..." << std::endl; + + G4PrimaryParticle* primary = new G4PrimaryParticle(); + primary->SetPDGcode(p->Pdg()); + primary->Set4Momentum(p->Px() * CLHEP::GeV, p->Py() * CLHEP::GeV, + p->Pz() * CLHEP::GeV, p->E() * CLHEP::GeV); + primary->SetProperTime(time_ * CLHEP::ns); + + UserPrimaryParticleInformation* primaryInfo = + new UserPrimaryParticleInformation(); + primaryInfo->setHepEvtStatus(1); + primary->SetUserInformation(primaryInfo); + + vertex->SetPrimary(primary); + } + + // add the vertex to the event + event->AddPrimaryVertex(vertex); + + ++n_events_generated_; + // delete genie_event; +} + +void GenieGenerator::RecordConfig(const std::string& id, ldmx::RunHeader& rh) { + rh.setStringParameter(id + " Class", "simcore::generators::GenieGenerator"); + rh.setStringParameter(id + "GenieTune", tune_); +} + +} // namespace generators +} // namespace simcore + +DECLARE_GENERATOR(simcore::generators::GenieGenerator) diff --git a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx new file mode 100644 index 000000000..9e412b65e --- /dev/null +++ b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx @@ -0,0 +1,203 @@ +// +// Created by Wesley Ketchum on 4/29/24. +// + +#include "SimCore/Reweight/GenieReweightProducer.h" + +#include +#include + +#include "Framework/EventGen/EventRecord.h" +#include "Framework/EventGen/HepMC3Converter.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Utils/RunOpt.h" +#include "HepMC3/GenEvent.h" +#include "HepMC3/Print.h" +#include "RwCalculators/GReWeightAGKY.h" +#include "RwCalculators/GReWeightDISNuclMod.h" +#include "RwCalculators/GReWeightDeltaradAngle.h" +#include "RwCalculators/GReWeightFGM.h" +#include "RwCalculators/GReWeightFZone.h" +#include "RwCalculators/GReWeightINuke.h" +#include "RwCalculators/GReWeightINukeParams.h" +#include "RwCalculators/GReWeightNonResonanceBkg.h" +#include "RwCalculators/GReWeightNuXSecCCQE.h" +#include "RwCalculators/GReWeightNuXSecCCQEaxial.h" +#include "RwCalculators/GReWeightNuXSecCCQEvec.h" +#include "RwCalculators/GReWeightNuXSecCCRES.h" +#include "RwCalculators/GReWeightNuXSecCOH.h" +#include "RwCalculators/GReWeightNuXSecDIS.h" +#include "RwCalculators/GReWeightNuXSecNC.h" +#include "RwCalculators/GReWeightNuXSecNCEL.h" +#include "RwCalculators/GReWeightNuXSecNCRES.h" +#include "RwCalculators/GReWeightResonanceDecay.h" +#include "RwCalculators/GReWeightXSecEmpiricalMEC.h" +#include "RwCalculators/GReWeightXSecMEC.h" +#include "RwFramework/GReWeight.h" +#include "RwFramework/GReWeightI.h" +#include "RwFramework/GSyst.h" +#include "RwFramework/GSystSet.h" +#include "SimCore/Event/EventWeights.h" +#include "SimCore/Event/HepMC3GenEvent.h" + +namespace simcore { + +GenieReweightProducer::GenieReweightProducer(const std::string& name, + framework::Process& process) + : Producer(name, process) { + hepMC3Converter_ = new genie::HepMC3Converter; + genie_rw_ = new genie::rew::GReWeight; +} + +GenieReweightProducer::~GenieReweightProducer() { + delete hepMC3Converter_; + delete genie_rw_; +} + +void GenieReweightProducer::configure(framework::config::Parameters& ps) { + hepmc3CollName_ = ps.getParameter("hepmc3CollName"); + hepmc3PassName_ = ps.getParameter("hepmc3PassName"); + eventWeightsCollName_ = ps.getParameter("eventWeightsCollName"); + verbosity_ = ps.getParameter("verbosity"); + seed_ = ps.getParameter("seed"); + n_weights_ = (size_t)(ps.getParameter("n_weights")); + auto var_types_strings = + ps.getParameter >("var_types"); + + std::default_random_engine generator(seed_); + std::normal_distribution normal_distribution; + + for (auto const& vt_str : var_types_strings) { + auto vtype = ldmx::EventWeights::string_to_variation_type(vt_str); + for (size_t i_w = 0; i_w < n_weights_; ++i_w) + variation_map_[vtype].push_back(normal_distribution(generator)); + } +} + +void GenieReweightProducer::reinitializeGenieReweight() { + delete genie_rw_; + genie_rw_ = new genie::rew::GReWeight; + + genie::RunOpt::Instance()->SetTuneName(tune_); + if (!genie::RunOpt::Instance()->Tune()) { + EXCEPTION_RAISE("ConfigurationException", "No TuneId in RunOption."); + } + genie::RunOpt::Instance()->BuildTune(); + + genie_rw_->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); + genie_rw_->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); + + auto& syst = genie_rw_->Systematics(); + for (auto var : variation_map_) + syst.Init(variation_type_to_genie_dial(var.first)); +} + +void GenieReweightProducer::onNewRun(const ldmx::RunHeader& runHeader) { + const std::string genie_tune_par = "GenieTune"; + std::string new_tune; + for (auto par : runHeader.getStringParameters()) + if (par.first.size() >= genie_tune_par.size() && + par.first.compare(par.first.size() - genie_tune_par.size(), + genie_tune_par.size(), genie_tune_par) == 0) { + new_tune = par.second; + break; + } + + if (tune_ != new_tune) { + if (verbosity_ >= 0) + std::cout << "Found new tune " << new_tune << " (used to be " << tune_ + << ")" << std::endl; + tune_ = new_tune; + reinitializeGenieReweight(); + } +} + +void GenieReweightProducer::reconfigureGenieReweight(size_t i_w) { + auto& syst = genie_rw_->Systematics(); + for (auto var : variation_map_) { + if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_MFP_pi) + syst.Set(genie::rew::GSyst::FromString("MFP_pi"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_MFP_N) + syst.Set(genie::rew::GSyst::FromString("MFP_N"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrCEx_pi) + syst.Set(genie::rew::GSyst::FromString("FrCEx_pi"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrInel_pi) + syst.Set(genie::rew::GSyst::FromString("FrInel_pi"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrAbs_pi) + syst.Set(genie::rew::GSyst::FromString("FrAbs_pi"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrPiProd_pi) + syst.Set(genie::rew::GSyst::FromString("FrPiProd_pi"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrCEx_N) + syst.Set(genie::rew::GSyst::FromString("FrCEx_N"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrInel_N) + syst.Set(genie::rew::GSyst::FromString("FrInel_N"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrAbs_N) + syst.Set(genie::rew::GSyst::FromString("FrAbs_N"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_INukeTwkDial_FrPiProd_N) + syst.Set(genie::rew::GSyst::FromString("FrPiProd_N"), var.second[i_w]); + else if (var.first == + ldmx::EventWeights::VariationType::kGENIE_HadrNuclTwkDial_FormZone) + syst.Set(genie::rew::GSyst::FromString("FormZone"), var.second[i_w]); + } + genie_rw_->Reconfigure(); +} + +void GenieReweightProducer::produce(framework::Event& event) { + // grab the input hepmc3 event collection + auto hepmc3_col = event.getObject >( + hepmc3CollName_, hepmc3PassName_); + + // create an output weights + ldmx::EventWeights ev_weights(variation_map_); + + for (size_t i_w = 0; i_w < n_weights_; ++i_w) { + double running_weight = 1; + + reconfigureGenieReweight(i_w); + + // setup a loop here ... but we're going to force only looping over one + // interaction if it exists for now. + for (size_t i_ev = 0; i_ev < 1; ++i_ev) { + auto const& hepmc3_ev = hepmc3_col.at(i_ev); + // fill our event data into a HepMC3GenEvent + HepMC3::GenEvent hepmc3_genev; + hepmc3_genev.read_data(hepmc3_ev); + + // print it out to check it ... + if (i_w == 0 && verbosity_ >= 1) + HepMC3::Print::line(hepmc3_genev, true); // print attributes + + // now convert to genie event record + auto genie_ev_record_ptr = hepMC3Converter_->RetrieveGHEP(hepmc3_genev); + + // print that out too ... + if (i_w == 0 && verbosity_ >= 1) genie_ev_record_ptr->Print(std::cout); + + // auto this_weight = 1.0 + var_value*0.05; + auto this_weight = genie_rw_->CalcWeight(*genie_ev_record_ptr); + + running_weight = running_weight * this_weight; + + } // end loop over interactions in event + + ev_weights.addWeight(running_weight); + + } // end loop over weights + + if (verbosity_ >= 1) ev_weights.Print(); + + event.add(eventWeightsCollName_, ev_weights); +} +} // namespace simcore + +DECLARE_PRODUCER_NS(simcore, GenieReweightProducer); diff --git a/SimCore/src/SimCore/Simulator.cxx b/SimCore/src/SimCore/Simulator.cxx index 43300c234..26466fe1d 100644 --- a/SimCore/src/SimCore/Simulator.cxx +++ b/SimCore/src/SimCore/Simulator.cxx @@ -20,6 +20,7 @@ /*~~~~~~~~~~~~~*/ #include "SimCore/APrimePhysics.h" #include "SimCore/DetectorConstruction.h" +#include "SimCore/Event/HepMC3GenEvent.h" #include "SimCore/G4Session.h" #include "SimCore/G4User/TrackingAction.h" #include "SimCore/Geo/ParserFactory.h" @@ -157,6 +158,21 @@ void Simulator::produce(framework::Event& event) { event_header.setStringParameter("eventSeed", stream.str()); + /* + PrimaryGenerator::Factory::get().apply([](auto gen){ + std::cout << gen->Name() << std::endl; + }); + */ + + auto event_info = static_cast( + runManager_->GetCurrentEvent()->GetUserInformation()); + + auto hepmc3_events = event_info->getHepMC3GenEvents(); + for (auto& hepmc3ev : hepmc3_events) { + hepmc3ev.event_number = event.getEventHeader().getEventNumber(); + } + event.add("SimHepMC3Events", hepmc3_events); + saveTracks(event); saveSDHits(event); diff --git a/cmake/BuildMacros.cmake b/cmake/BuildMacros.cmake index 3716712e4..347107b30 100644 --- a/cmake/BuildMacros.cmake +++ b/cmake/BuildMacros.cmake @@ -67,6 +67,162 @@ macro(setup_geant4_target) endmacro() + +function(add_genie_target gname) + +#message(STATUS "Adding library Genie::${gname}") + +add_library(Genie::${gname} SHARED IMPORTED GLOBAL) # or STATIC instead of SHARED +set_target_properties(Genie::${gname} PROPERTIES + IMPORTED_LOCATION "/usr/local/lib/lib${gname}.so" + INTERFACE_INCLUDE_DIRECTORIES "/usr/local/include/GENIE" +) +endfunction() + +macro(setup_genie_target) +if(NOT DEFINED GENIE_LIBS) + +find_package(ROOT CONFIG REQUIRED) + +add_genie_target(GFwMsg) +add_genie_target(GFwReg) +add_genie_target(GFwAlg) +add_genie_target(GFwInt) +add_genie_target(GFwGHEP) +add_genie_target(GFwNum) +add_genie_target(GFwUtl) +add_genie_target(GFwParDat) +add_genie_target(GFwEG) +add_genie_target(GFwNtp) + +add_genie_target(GPhCmn) +add_genie_target(GPhDcy) +add_genie_target(GPhHadTens) +add_genie_target(GPhMEL) +add_genie_target(GPhPDF) +add_genie_target(GPhXSIg) +add_genie_target(GPhNuclSt) +add_genie_target(GPhDeEx) +add_genie_target(GPhHadnz) +add_genie_target(GPhHadTransp) +add_genie_target(GPhAMNGXS) +add_genie_target(GPhAMNGEG) +add_genie_target(GPhChmXS) +add_genie_target(GPhCohXS) +add_genie_target(GPhCohEG) +add_genie_target(GPhDISXS) +add_genie_target(GPhDISEG) +add_genie_target(GPhDfrcXS) +add_genie_target(GPhDfrcEG) +add_genie_target(GPhHELptnXS) +add_genie_target(GPhHELptnEG) +add_genie_target(GPhIBDXS) +add_genie_target(GPhIBDEG) +add_genie_target(GPhMNucXS) +add_genie_target(GPhMNucEG) +add_genie_target(GPhNuElXS) +add_genie_target(GPhNuElEG) +add_genie_target(GPhQELXS) +add_genie_target(GPhQELEG) +add_genie_target(GPhResXS) +add_genie_target(GPhResEG) +add_genie_target(GPhStrXS) +add_genie_target(GPhStrEG) +add_genie_target(GPhHEDISXS) +add_genie_target(GPhHEDISEG) +add_genie_target(GTlFlx) +add_genie_target(GTlGeo) + +add_genie_target(GRwFwk) +add_genie_target(GRwIO) +add_genie_target(GRwClc) + +add_library(Pythia6 SHARED IMPORTED GLOBAL) # or STATIC instead of SHARED +set_target_properties(Pythia6 PROPERTIES + IMPORTED_LOCATION "/usr/local/pythia6/libPythia6.so" +) + +add_library(blas SHARED IMPORTED GLOBAL) # or STATIC instead of SHARED +set_target_properties(blas PROPERTIES + IMPORTED_LOCATION "/usr/lib/x86_64-linux-gnu/libblas.so.3" +) + +#add_library(Genie::Interface INTERFACE IMPORTED GLOBAL) +#set_property(TARGET Genie::Interface PROPERTY +# INTERFACE_LINK_LIBRARIES +SET( GENIE_LIBS +# xml2 + log4cpp +# LHAPDF +# Pythia6 + gsl +# Genie::GFwMsg +# Genie::GFwReg +# Genie::GFwAlg + + Genie::GRwFwk + Genie::GRwClc + Genie::GRwIO + + Genie::GFwInt +# Genie::GFwGHEP + Genie::GFwNum + Genie::GFwUtl + Genie::GFwParDat + Genie::GFwEG + Genie::GFwAlg + Genie::GFwGHEP + Genie::GFwMsg + xml2 + Genie::GFwReg + ROOT::EG + ROOT::EGPythia6 + Pythia6 + Genie::GFwNtp + Genie::GPhXSIg + Genie::GPhPDF + Genie::GPhNuclSt + Genie::GPhCmn + Genie::GPhDcy + Genie::GPhHadTransp + Genie::GPhHadnz + Genie::GPhHadTens + Genie::GPhDeEx + Genie::GPhAMNGXS + Genie::GPhAMNGEG + Genie::GPhChmXS + Genie::GPhCohXS + Genie::GPhCohEG + Genie::GPhDISXS + Genie::GPhDISEG + Genie::GPhDfrcXS + Genie::GPhDfrcEG + Genie::GPhHELptnXS + Genie::GPhHELptnEG + Genie::GPhIBDXS + Genie::GPhIBDEG + Genie::GPhMNucXS + Genie::GPhMNucEG + Genie::GPhMEL + Genie::GPhNuElXS + Genie::GPhNuElEG + Genie::GPhQELXS + Genie::GPhQELEG + Genie::GPhResXS + Genie::GPhResEG + Genie::GPhStrXS + Genie::GPhStrEG + Genie::GPhHEDISXS + Genie::GPhHEDISEG + Genie::GTlGeo + Genie::GTlFlx + ROOT::MathCore + ROOT::MathMore + blas) +message(STATUS "Setting GENIE_LIBS to ${GENIE_LIBS}") +endif() +endmacro() + macro(setup_lcio_target) # If it doesn't exists, create an imported target for LCIO @@ -97,7 +253,7 @@ endmacro() macro(setup_library) - set(options interface register_target) + set(options interface register_target include) set(oneValueArgs module name) set(multiValueArgs dependencies sources) cmake_parse_arguments(setup_library "${options}" "${oneValueArgs}" @@ -131,14 +287,19 @@ macro(setup_library) # Setup the include directories if(setup_library_interface) target_include_directories(${library_name} - INTERFACE ${PROJECT_SOURCE_DIR}/include) + INTERFACE ${PROJECT_SOURCE_DIR}/include) else() target_include_directories(${library_name} PUBLIC ${PROJECT_SOURCE_DIR}/include) endif() # Setup the targets to link against - target_link_libraries(${library_name} PUBLIC ${setup_library_dependencies}) + if(setup_library_interface) + target_link_libraries(${library_name} INTERFACE ${setup_library_dependencies}) + else() + target_link_libraries(${library_name} PUBLIC ${setup_library_dependencies}) + endif() + enable_sanitizers(${library_name}) enable_compiler_warnings(${library_name}) enable_ipo(${library_name})