From 480427dce06867e73397fb29acdb1508a45fcf03 Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Wed, 15 May 2024 13:05:23 +0200 Subject: [PATCH 01/14] transfer over of GenieGenerator --- SimCore/CMakeLists.txt | 18 +- SimCore/include/SimCore/Event/GHepParticle.h | 62 +++ SimCore/include/SimCore/Event/GTruth.h | 154 ++++++ .../SimCore/Generators/GenieGenerator.h | 111 +++++ SimCore/include/SimCore/PrimaryGenerator.h | 4 + SimCore/include/SimCore/Simulator.h | 13 + .../include/SimCore/UserEventInformation.h | 17 +- SimCore/python/generators.py | 72 +++ SimCore/src/SimCore/Event/GHepParticle.cxx | 14 + SimCore/src/SimCore/Event/GTruth.cxx | 14 + .../SimCore/G4User/PrimaryGeneratorAction.cxx | 7 +- .../src/SimCore/Generators/GenieGenerator.cxx | 447 ++++++++++++++++++ SimCore/src/SimCore/Simulator.cxx | 163 +++++++ cmake/BuildMacros.cmake | 165 ++++++- 14 files changed, 1249 insertions(+), 12 deletions(-) create mode 100644 SimCore/include/SimCore/Event/GHepParticle.h create mode 100644 SimCore/include/SimCore/Event/GTruth.h create mode 100644 SimCore/include/SimCore/Generators/GenieGenerator.h create mode 100644 SimCore/src/SimCore/Event/GHepParticle.cxx create mode 100644 SimCore/src/SimCore/Event/GTruth.cxx create mode 100644 SimCore/src/SimCore/Generators/GenieGenerator.cxx diff --git a/SimCore/CMakeLists.txt b/SimCore/CMakeLists.txt index c80435e43..4cfc110d5 100644 --- a/SimCore/CMakeLists.txt +++ b/SimCore/CMakeLists.txt @@ -15,6 +15,8 @@ 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) @@ -24,11 +26,16 @@ if(BUILD_EVENT_ONLY) 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 "GTruth") + register_event_object(module_path "SimCore/Event" namespace "ldmx" + class "GHepParticle" 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} + register_target + include ) return() @@ -37,6 +44,9 @@ endif() # Configure Geant4 setup_geant4_target() +#Configure Genie +setup_genie_target() + # include dark brem simulation add_subdirectory(G4DarkBreM) @@ -64,12 +74,14 @@ file(GLOB SRC_FILES CONFIGURE_DEPDENDS ${PROJECT_SOURCE_DIR}/src/SimCore/Geo/[a- setup_library(module SimCore dependencies Geant4::Interface ROOT::Physics + ${GENIE_LIBS} Framework::Configure Framework::Framework SimCore::G4User G4DarkBreM DetDescr::DetDescr Boost::log + log4cpp "${registered_targets}" sources ${SRC_FILES}) @@ -81,7 +93,7 @@ setup_library(module SimCore # Primary Generators library setup_library(module SimCore name Generators - dependencies SimCore::SimCore SimCore::LHE) + dependencies SimCore::SimCore SimCore::LHE ${GENIE_LIBS}) # Sensitive Detectors library setup_library(module SimCore diff --git a/SimCore/include/SimCore/Event/GHepParticle.h b/SimCore/include/SimCore/Event/GHepParticle.h new file mode 100644 index 000000000..b15606d9e --- /dev/null +++ b/SimCore/include/SimCore/Event/GHepParticle.h @@ -0,0 +1,62 @@ +//////////////////////////////////////////////////////////////////////// +/// \file GHepParticle.h +/// \Class to hold the particle information needed to recreate a genie::EventRecord +/// \author wketchum@fnal.gov +/// +/// See https://github.com/GENIE-MC/Generator/blob/master/src/Framework/GHEP/GHepParticle.h +/// +/// Adopted for LDMX-SW by wketchum@fnal.gov +//////////////////////////////////////////////////////////////////////// + +#ifndef SIMCORE_EVENT_GHEPPARTICLE_H +#define SIMCORE_EVENT_GHEPPARTICLE_H + +#include "TObject.h" +#include + +namespace ldmx { + + class GHepParticle { + + public: + + GHepParticle() {} + + void Clear() {} + void Print() const; + + int fPosition; + + int fPdgCode; + int fStatus; + int fRescatterCode; + int fFirstMother; + int fLastMother; + int fFirstDaugher; + int fLastDaughter; + + double fP_x; + double fP_y; + double fP_z; + double fP_t; + + double fX_x; + double fX_y; + double fX_z; + double fX_t; + + double fPolzTheta; + double fPolzPhi; + + double fRemovalEnergy; + bool fIsBound; + + public: + + ClassDef(GHepParticle, 1); + + }; + +} // end simb namespace + +#endif // SIMCORE_EVENT_GTRUTH_H diff --git a/SimCore/include/SimCore/Event/GTruth.h b/SimCore/include/SimCore/Event/GTruth.h new file mode 100644 index 000000000..aa63e9f57 --- /dev/null +++ b/SimCore/include/SimCore/Event/GTruth.h @@ -0,0 +1,154 @@ +//////////////////////////////////////////////////////////////////////// +/// \file GTruth.h +/// \Class to hold the additional information needed to recreate a genie::EventRecord +/// \author nathan.mayer@tufts.edu +/// +/// See https://github.com/NuSoftHEP/nusimdata/blob/develop/nusimdata/SimulationBase/GTruth.h +/// +/// Adopted for LDMX-SW by wketchum@fnal.gov +//////////////////////////////////////////////////////////////////////// + +/// This class stores/retrieves the additional information needed +/// (and not in MCTruth) to recreate a genie::EventRecord +/// for genie based event reweighting. + +#ifndef SIMCORE_EVENT__GTRUTH_H +#define SIMCORE_EVENT__GTRUTH_H + +#include "TObject.h" + +#include +#include + +namespace ldmx { + + class GTruth { + + public: + + GTruth() {} + + void Clear() {} + void Print() const; + + // genie::GHepRecord info + // holds a genie::Interaction + + //TLorentzVector fVertex; + double fVertex_x; + double fVertex_y; + double fVertex_z; + double fVertex_t; + + + // skipping TBits data members EventFlags and EventMask + double fweight; ///< event interaction weight (genie internal) + double fprobability; ///< interaction probability + double fXsec; ///< cross section of interaction + double fDiffXsec; ///< differential cross section of interaction + int fGPhaseSpace; ///< phase space system of DiffXSec + + // genie::Interaction + // container for InitialState, ProcessInfo, XclsTag, Kinematics, KPhaseSpace) + // holds no fundamental type info + + // genie:::InitialState info (sub-object to genie::Interactions) + int fProbePDG; + // holds a genie::Target + + //TLorentzVector fProbeP4; + double fProbe_px; + double fProbe_py; + double fProbe_pz; + double fProbe_e; + + //TLorentzVector fTgtP4; // added version 13 + double fTgt_px; + double fTgt_py; + double fTgt_pz; + double fTgt_e; + + // genie::Target info (sub-object to genie::InitialState) + int ftgtZ; + int ftgtA; + int ftgtPDG; ///< PDG of Target Nucleus, nucleon only if free + int fHitNucPDG; ///< hit nucleon PDG code // added version 13 + int fHitQrkPDG; ///< hit quark PDG code // added version 13 + bool fIsSeaQuark; + //TLorentzVector fHitNucP4; + double fHitNuc_px; + double fHitNuc_py; + double fHitNuc_pz; + double fHitNuc_e; + + double fHitNucPos; // added version 12 + + // genie::ProcessInfo (sub-object to genie::Interactions) + int fGscatter; ///< neutrino scattering code + int fGint; ///< interaction code + + // genie::Kinematics info (sub-object to genie::Interactions) + + ///< these are for the internal (on shell) genie kinematics + ///< this list might be an incomplete transcription of map + //double fgQ2; + //double fgq2; + //double fgW; + //double fgT; + //double fgX; + //double fgY; + std::map fKV; + + ///< a common running variable to be recorded + //double fgWrun; + + //TLorentzVector fFSleptonP4; ///< generated final state primary lepton (LAB frame) // added version 13 + double fFSlepton_px; + double fFSlepton_py; + double fFSlepton_pz; + double fFSlepton_e; + + //TLorentzVector fFShadSystP4; ///< generated final state hadronic system (LAB frame) + double fFShadSyst_px; + double fFShadSyst_py; + double fFShadSyst_pz; + double fFShadSyst_e; + + // genie::XclsTag info (sub-object to genie::Interactions) + bool fIsCharm; ///< did the interaction produce a charmed hadron? + int fCharmHadronPdg; // added version 13 + bool fIsStrange; ///< strange production // added version 13 + int fStrangeHadronPdg; // added version 13 + int fNumProton; ///< number of protons after reaction, before FSI + int fNumNeutron; ///< number of neutrons after reaction, before FSI + int fNumPi0; ///< number of pi0 after reaction, before FSI + int fNumPiPlus; ///< number of pi pluses after reaction, before FSI + int fNumPiMinus; ///< number of pi minuses after reaction, before FSI + int fNumSingleGammas; ///< number of gammas after reaction, before FSI + int fNumRho0; ///< number of pi0 after reaction, before FSI + int fNumRhoPlus; ///< number of pi pluses after reaction, before FSI + int fNumRhoMinus; ///< number of pi minuses after reaction, before FSI + int fResNum; ///< resonance number + int fDecayMode; // added version 13 + bool fIsFinalQuarkEvent; + int fFinalQuarkPdg; + bool fIsFinalLeptonEvent; + int fFinalLeptonPdg; + + + // genie::KPhaseSpace (sub-object to genie::Interactions) + // has no relevant private data + + // Flag for values that might not have been set + static constexpr double kUndefinedValue = -99999; + + public: + //friend std::ostream& operator<< (std::ostream& output, const simb::GTruth >ruth); + + ClassDef(GTruth, 1); + + }; + +} // end simb namespace + +#endif // SIMCORE_EVENT__GTRUTH_H diff --git a/SimCore/include/SimCore/Generators/GenieGenerator.h b/SimCore/include/SimCore/Generators/GenieGenerator.h new file mode 100644 index 000000000..8bd928ceb --- /dev/null +++ b/SimCore/include/SimCore/Generators/GenieGenerator.h @@ -0,0 +1,111 @@ +/** + * @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" + +//------------// +// LDMX // +//------------// +#include "SimCore/PrimaryGenerator.h" + +#include +#include + +// 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 + */ + genie::GEVGDriver evg_driver_; + + 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..6b70dca7a 100644 --- a/SimCore/include/SimCore/PrimaryGenerator.h +++ b/SimCore/include/SimCore/PrimaryGenerator.h @@ -23,6 +23,8 @@ #include "Framework/RunHeader.h" #include "SimCore/Factory.h" +#include "UserEventInformation.h" + // Forward Declarations class G4Event; @@ -69,6 +71,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/Simulator.h b/SimCore/include/SimCore/Simulator.h index 7a6b7d416..ecf763bbb 100644 --- a/SimCore/include/SimCore/Simulator.h +++ b/SimCore/include/SimCore/Simulator.h @@ -25,6 +25,13 @@ #include "SimCore/DetectorConstruction.h" #include "SimCore/RunManager.h" #include "SimCore/SimulatorBase.h" +#include "SimCore/Event/GTruth.h" +#include "SimCore/Event/GHepParticle.h" + +namespace genie { + class EventRecord; + class GHepParticle; +} class G4UImanager; class G4UIsession; @@ -116,6 +123,12 @@ class Simulator : public SimulatorBase { */ void setSeeds(std::vector seeds); + void FillGTruth(const genie::EventRecord* record, + ldmx::GTruth& truth); + + void FillGHepParticles(const genie::EventRecord* record, + std::vector& ghep_particle_list); + private: /// Number of events started int numEventsBegan_{0}; diff --git a/SimCore/include/SimCore/UserEventInformation.h b/SimCore/include/SimCore/UserEventInformation.h index 65e9e0080..9d38bfb78 100644 --- a/SimCore/include/SimCore/UserEventInformation.h +++ b/SimCore/include/SimCore/UserEventInformation.h @@ -2,6 +2,9 @@ #define SIMCORE_USEREVENTINFORMATION_H #include "G4VUserEventInformation.hh" + +#include "GENIE/Framework/EventGen/EventRecord.h" + namespace simcore { /** @@ -13,7 +16,8 @@ class UserEventInformation : public G4VUserEventInformation { UserEventInformation() = default; /// Destructor - virtual ~UserEventInformation() = default; + virtual ~UserEventInformation() + { if(genie_event_) delete genie_event_;} /// Print the information associated with the track void Print() const final override; @@ -116,7 +120,10 @@ class UserEventInformation : public G4VUserEventInformation { * @returns true if it was */ bool wasLastStepEN() const { return last_step_en_; } - + + void setGENIEEventRecord(genie::EventRecord *event) { genie_event_ = event; } + genie::EventRecord* getGENIEEventRecord() { return genie_event_; } + private: /// Total number of brem candidates in the event int bremCandidateCount_{0}; @@ -165,6 +172,12 @@ class UserEventInformation : public G4VUserEventInformation { * dark brem did not occur within the event in question. */ double db_material_z_{-1.}; + + /** + * Pointer for a GENIE event record. + */ + genie::EventRecord* genie_event_{NULL}; + }; } // namespace simcore diff --git a/SimCore/python/generators.py b/SimCore/python/generators.py index 3ab1f6710..099bc77ea 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/src/SimCore/Event/GHepParticle.cxx b/SimCore/src/SimCore/Event/GHepParticle.cxx new file mode 100644 index 000000000..b016b353f --- /dev/null +++ b/SimCore/src/SimCore/Event/GHepParticle.cxx @@ -0,0 +1,14 @@ +#include "SimCore/Event/GHepParticle.h" + +#include + +ClassImp(ldmx::GHepParticle) + +namespace ldmx{ + + void GHepParticle::Print() const { + std::cout << "GHEPPARTICLE! TODO!" << std::endl ; + + } + +} diff --git a/SimCore/src/SimCore/Event/GTruth.cxx b/SimCore/src/SimCore/Event/GTruth.cxx new file mode 100644 index 000000000..6c996f447 --- /dev/null +++ b/SimCore/src/SimCore/Event/GTruth.cxx @@ -0,0 +1,14 @@ +#include "SimCore/Event/GTruth.h" + +#include + +ClassImp(ldmx::GTruth) + +namespace ldmx{ + + void GTruth::Print() const { + std::cout << "GTRUTH! TODO!" << std::endl ; + + } + +} diff --git a/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx b/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx index bc925a265..341782377 100644 --- a/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx +++ b/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx @@ -70,12 +70,11 @@ 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); - + PrimaryGenerator::Factory::get().apply([event](const auto& generator) { - generator->GeneratePrimaryVertex(event); + generator->GeneratePrimaryVertex(event); }); // smear all primary vertices (if activated) diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx new file mode 100644 index 000000000..2b85cb870 --- /dev/null +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -0,0 +1,447 @@ +/** + * @file GenieGenerator.cxx + * @brief Simple GENIE event generator. + * @author Wesley Ketchum, FNAL + */ + +#include "SimCore/Generators/GenieGenerator.h" +#include "SimCore/UserPrimaryParticleInformation.h" +#include "SimCore/Event/GTruth.h" + +// GENIE +#include "Framework/Utils/AppInit.h" +#include "GENIE/Framework/Utils/RunOpt.h" +#include "Framework/Utils/XSecSplineList.h" +#include "GENIE/Framework/Interaction/InitialState.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepParticle.h" + + +#include "Framework/Conventions/XmlParserStatus.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/Conventions/Controls.h" +#include "Framework/Conventions/Units.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/EventGen/GFluxI.h" +#include "Framework/EventGen/GEVGDriver.h" +#include "Framework/EventGen/GMCJDriver.h" +#include "Framework/EventGen/GMCJMonitor.h" +#include "Framework/EventGen/InteractionList.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Ntuple/NtpWriter.h" +#include "Framework/Ntuple/NtpMCFormat.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/RunOpt.h" +#include "Framework/Utils/XSecSplineList.h" +#include "Framework/Utils/StringUtils.h" +#include "Framework/Utils/PrintUtils.h" +#include "Framework/Utils/SystemUtils.h" +#include "Framework/Utils/CmdLnArgParser.h" + +// Geant4 +#include "G4Event.hh" +#include "G4PhysicalConstants.hh" +#include "Randomize.hh" + + +// ROOT +#include +#include +#include + +// standard +#include +#include +#include + +namespace simcore { +namespace generators { + +void GenieGenerator::fillConfig(const framework::config::Parameters& p) +{ + + energy_ = p.getParameter("energy");// * GeV; + + targets_ = p.getParameter< std::vector >("targets"); + abundances_ = p.getParameter< std::vector >("abundances"); + + time_ = p.getParameter("time");// * ns; + position_ = p.getParameter< std::vector >("position"); // mm + beam_size_ = p.getParameter< std::vector >("beam_size"); // mm + direction_ = p.getParameter< std::vector >("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; i_a0) + 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; i_d(""), + 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_tgetSeed(); + 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_ - (double)initial_e.GetMass()*(double)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 genie_info = new UserEventInformation; + genie_info->setGENIEEventRecord(genie_event); + event->SetUserInformation(genie_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_pStatus()!=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"); +} + +} // namespace generators +} // namespace simcore + +DECLARE_GENERATOR(simcore::generators::GenieGenerator) diff --git a/SimCore/src/SimCore/Simulator.cxx b/SimCore/src/SimCore/Simulator.cxx index d8d1c2b61..bab21838c 100644 --- a/SimCore/src/SimCore/Simulator.cxx +++ b/SimCore/src/SimCore/Simulator.cxx @@ -15,6 +15,12 @@ #include "Framework/RandomNumberSeedService.h" #include "Framework/Version.h" //for LDMX_INSTALL path + +/*~~~~~~~~~~~~~*/ +/* GENIE */ +/*~~~~~~~~~~~~~*/ +#include "GENIE/Framework/GHEP/GHepParticle.h" + /*~~~~~~~~~~~~~*/ /* SimCore */ /*~~~~~~~~~~~~~*/ @@ -27,6 +33,8 @@ #include "SimCore/SensitiveDetector.h" #include "SimCore/UserEventInformation.h" #include "SimCore/XsecBiasingOperator.h" +#include "SimCore/Event/GTruth.h" +#include "SimCore/Event/GHepParticle.h" /*~~~~~~~~~~~~~~*/ /* Geant4 */ @@ -128,6 +136,141 @@ void Simulator::onNewRun(const ldmx::RunHeader& runHeader) { run_ = runHeader.getRunNumber(); } +void Simulator::FillGHepParticles(const genie::EventRecord* record, + std::vector& ghep_particle_list) +{ + for(int position=0; position< record->GetEntries(); ++position) + { + const genie::GHepParticle* particle = record->Particle(position); + + ldmx::GHepParticle my_particle; + + my_particle.fPosition = position; + + my_particle.fPdgCode = particle->Pdg(); + my_particle.fStatus = particle->Status(); + my_particle.fRescatterCode = particle->RescatterCode(); + my_particle.fFirstMother = particle->FirstMother(); + my_particle.fLastMother = particle->LastMother(); + my_particle.fFirstDaugher = particle->FirstDaughter(); + my_particle.fLastDaughter = particle->LastDaughter(); + + my_particle.fP_x = particle->P4()->Px(); + my_particle.fP_y = particle->P4()->Py(); + my_particle.fP_z = particle->P4()->Pz(); + my_particle.fP_t = particle->P4()->E(); + + my_particle.fX_x = particle->X4()->X(); + my_particle.fX_y = particle->X4()->Y(); + my_particle.fX_z = particle->X4()->Z(); + my_particle.fX_t = particle->X4()->T(); + + my_particle.fPolzTheta = particle->PolzPolarAngle(); + my_particle.fPolzPhi = particle->PolzAzimuthAngle(); + + my_particle.fRemovalEnergy = particle->RemovalEnergy(); + my_particle.fIsBound = particle->IsBound(); + + ghep_particle_list.push_back(my_particle); + } +} + +void Simulator::FillGTruth(const genie::EventRecord* record, + ldmx::GTruth& truth) { + + //interactions info + genie::Interaction *inter = record->Summary(); + const genie::ProcessInfo &procInfo = inter->ProcInfo(); + truth.fGint = (int)procInfo.InteractionTypeId(); + truth.fGscatter = (int)procInfo.ScatteringTypeId(); + + //Event info + truth.fweight = record->Weight(); + truth.fprobability = record->Probability(); + truth.fXsec = record->XSec(); + truth.fDiffXsec = record->DiffXSec(); + truth.fGPhaseSpace = (int)record->DiffXSecVars(); + + //Initial State info + const genie::InitialState &initState = inter->InitState(); + truth.fProbePDG = initState.ProbePdg(); + truth.fProbe_px = initState.GetProbeP4()->Px(); + truth.fProbe_py = initState.GetProbeP4()->Py(); + truth.fProbe_pz = initState.GetProbeP4()->Pz(); + truth.fProbe_e = initState.GetProbeP4()->E(); + + truth.fTgt_px = initState.GetTgtP4()->Px(); + truth.fTgt_py = initState.GetTgtP4()->Py(); + truth.fTgt_pz = initState.GetTgtP4()->Pz(); + truth.fTgt_e = initState.GetTgtP4()->E(); + + //Target info + const genie::Target &tgt = initState.Tgt(); + truth.ftgtZ = tgt.Z(); + truth.ftgtA = tgt.A(); + truth.ftgtPDG = tgt.Pdg(); + truth.fHitNucPDG = tgt.HitNucPdg(); + truth.fHitQrkPDG = tgt.HitQrkPdg(); + truth.fIsSeaQuark = tgt.HitSeaQrk(); + truth.fHitNuc_px = tgt.HitNucP4Ptr()->Px(); + truth.fHitNuc_py = tgt.HitNucP4Ptr()->Py(); + truth.fHitNuc_pz = tgt.HitNucP4Ptr()->Pz(); + truth.fHitNuc_e = tgt.HitNucP4Ptr()->E(); + truth.fHitNucPos = tgt.HitNucPosition(); + + + truth.fVertex_x = record->Vertex()->X(); + truth.fVertex_y = record->Vertex()->Y(); + truth.fVertex_z = record->Vertex()->Z(); + truth.fVertex_t = record->Vertex()->T(); + + //true reaction information and byproducts + //(PRE FSI) + const genie::XclsTag &exclTag = inter->ExclTag(); + truth.fIsCharm = exclTag.IsCharmEvent(); + truth.fCharmHadronPdg = exclTag.CharmHadronPdg(); + truth.fIsStrange = exclTag.IsStrangeEvent(); + truth.fStrangeHadronPdg = exclTag.StrangeHadronPdg(); + truth.fResNum = (int)exclTag.Resonance(); + truth.fDecayMode = exclTag.DecayMode(); + + truth.fNumPiPlus = exclTag.NPiPlus(); + truth.fNumPiMinus = exclTag.NPiMinus(); + truth.fNumPi0 = exclTag.NPi0(); + truth.fNumProton = exclTag.NProtons(); + truth.fNumNeutron = exclTag.NNeutrons(); + truth.fNumSingleGammas = exclTag.NSingleGammas(); + + truth.fNumSingleGammas = exclTag.NSingleGammas(); + truth.fNumRho0 = exclTag.NRho0(); + truth.fNumRhoPlus = exclTag.NRhoPlus(); + truth.fNumRhoMinus = exclTag.NRhoMinus(); + + truth.fIsFinalQuarkEvent = exclTag.IsFinalQuarkEvent(); + truth.fFinalQuarkPdg = exclTag.FinalQuarkPdg(); + truth.fIsFinalLeptonEvent = exclTag.IsFinalLeptonEvent(); + truth.fFinalLeptonPdg = exclTag.FinalLeptonPdg(); + + // Get the GENIE kinematics info + const genie::Kinematics &kine = inter->Kine(); + for(int kvar=genie::KineVar_t::kKVNull; kvar!=genie::KineVar_t::kNumOfKineVar; ++kvar) + truth.fKV[kvar] = kine.GetKV((genie::KineVar_t)kvar); + + truth.fFSlepton_px = kine.FSLeptonP4().Px(); + truth.fFSlepton_py = kine.FSLeptonP4().Py(); + truth.fFSlepton_pz = kine.FSLeptonP4().Pz(); + truth.fFSlepton_e = kine.FSLeptonP4().E(); + + truth.fFShadSyst_px = kine.HadSystP4().Px(); + truth.fFShadSyst_py = kine.HadSystP4().Py(); + truth.fFShadSyst_pz = kine.HadSystP4().Pz(); + truth.fFShadSyst_e = kine.HadSystP4().E(); + + return; + +} + + void Simulator::produce(framework::Event& event) { // Generate and process a Geant4 event. numEventsBegan_++; @@ -157,6 +300,26 @@ 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()); + if(event_info->getGENIEEventRecord()){ + //event_info->getGENIEEventRecord()->Print(); + + ldmx::GTruth gtruth; + std::vector ghep_particles; + + FillGTruth(event_info->getGENIEEventRecord(),gtruth); + FillGHepParticles(event_info->getGENIEEventRecord(),ghep_particles); + + event.add("SimGTruth", gtruth); + event.add("SimGHepParticles", ghep_particles); + } saveTracks(event); saveSDHits(event); diff --git a/cmake/BuildMacros.cmake b/cmake/BuildMacros.cmake index 3716712e4..fb7102ba6 100644 --- a/cmake/BuildMacros.cmake +++ b/cmake/BuildMacros.cmake @@ -67,6 +67,160 @@ 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 + ROOT::MathCore + ROOT::MathMore + 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 + Genie::GRwFwk + Genie::GRwClc + Genie::GRwIO + 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 +251,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 +285,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}) From 11f36ec082f8253036559fa1ba7fc38d8ddf2ccc Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Wed, 15 May 2024 22:21:15 +0200 Subject: [PATCH 02/14] remove genie reweight libraries for now --- cmake/BuildMacros.cmake | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cmake/BuildMacros.cmake b/cmake/BuildMacros.cmake index fb7102ba6..023d57e76 100644 --- a/cmake/BuildMacros.cmake +++ b/cmake/BuildMacros.cmake @@ -213,9 +213,9 @@ SET( GENIE_LIBS Genie::GPhHEDISEG Genie::GTlGeo Genie::GTlFlx - Genie::GRwFwk - Genie::GRwClc - Genie::GRwIO +# Genie::GRwFwk +# Genie::GRwClc +# Genie::GRwIO blas) message(STATUS "Setting GENIE_LIBS to ${GENIE_LIBS}") endif() From 7a54c71c30a6cb55fda1c159b4d6c9a2d51c2369 Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Wed, 17 Jul 2024 08:30:35 +0200 Subject: [PATCH 03/14] transfer from wketchum/Genie_Generator_HepMC3Refactor branch in standalone SimCore branch --- SimCore/CMakeLists.txt | 21 ++- SimCore/include/SimCore/Event/EventWeights.h | 82 +++++++++ SimCore/include/SimCore/Event/GHepParticle.h | 62 ------- SimCore/include/SimCore/Event/GTruth.h | 154 ----------------- .../include/SimCore/Event/HepMC3GenEvent.h | 33 ++++ .../SimCore/Generators/GenieGenerator.h | 4 +- .../SimCore/Reweight/GenieReweightProducer.h | 70 ++++++++ SimCore/include/SimCore/Simulator.h | 15 +- .../include/SimCore/UserEventInformation.h | 16 +- SimCore/python/genie_reweight.py | 12 ++ SimCore/src/SimCore/Event/EventWeights.cxx | 25 +++ SimCore/src/SimCore/Event/GHepParticle.cxx | 14 -- SimCore/src/SimCore/Event/GTruth.cxx | 14 -- SimCore/src/SimCore/Event/HepMC3GenEvent.cxx | 26 +++ .../src/SimCore/Generators/GenieGenerator.cxx | 12 +- .../Reweight/GenieReweightProducer.cxx | 135 +++++++++++++++ SimCore/src/SimCore/Simulator.cxx | 158 +----------------- 17 files changed, 427 insertions(+), 426 deletions(-) create mode 100644 SimCore/include/SimCore/Event/EventWeights.h delete mode 100644 SimCore/include/SimCore/Event/GHepParticle.h delete mode 100644 SimCore/include/SimCore/Event/GTruth.h create mode 100644 SimCore/include/SimCore/Event/HepMC3GenEvent.h create mode 100644 SimCore/include/SimCore/Reweight/GenieReweightProducer.h create mode 100644 SimCore/python/genie_reweight.py create mode 100644 SimCore/src/SimCore/Event/EventWeights.cxx delete mode 100644 SimCore/src/SimCore/Event/GHepParticle.cxx delete mode 100644 SimCore/src/SimCore/Event/GTruth.cxx create mode 100644 SimCore/src/SimCore/Event/HepMC3GenEvent.cxx create mode 100644 SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx diff --git a/SimCore/CMakeLists.txt b/SimCore/CMakeLists.txt index 4cfc110d5..b25211661 100644 --- a/SimCore/CMakeLists.txt +++ b/SimCore/CMakeLists.txt @@ -20,6 +20,8 @@ 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" @@ -27,13 +29,14 @@ if(BUILD_EVENT_ONLY) 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 "GTruth") + class "EventWeights") + register_event_object(module_path "SimCore/Event" namespace "ldmx" - class "GHepParticle" type "collection") + class "HepMC3GenEvent" type "collection") # Generate the files needed to build the event classes. setup_library(module SimCore name Event - dependencies ROOT::Core ${GENIE_LIBS} + dependencies ROOT::Core ${GENIE_LIBS} HepMC3 HepMC3rootIO register_target include ) @@ -74,14 +77,15 @@ file(GLOB SRC_FILES CONFIGURE_DEPDENDS ${PROJECT_SOURCE_DIR}/src/SimCore/Geo/[a- setup_library(module SimCore dependencies Geant4::Interface ROOT::Physics - ${GENIE_LIBS} + HepMC3 + ${GENIE_LIBS} Framework::Configure Framework::Framework SimCore::G4User G4DarkBreM DetDescr::DetDescr Boost::log - log4cpp + log4cpp "${registered_targets}" sources ${SRC_FILES}) @@ -93,7 +97,12 @@ setup_library(module SimCore # Primary Generators library setup_library(module SimCore name Generators - dependencies SimCore::SimCore SimCore::LHE ${GENIE_LIBS}) + dependencies SimCore::SimCore SimCore::LHE ROOT::MathMore ${GENIE_LIBS} HepMC3) + +# Reweight Library +setup_library(module SimCore + name Reweight + dependencies SimCore::SimCore SimCore::Generators 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..9cc2ea122 --- /dev/null +++ b/SimCore/include/SimCore/Event/EventWeights.h @@ -0,0 +1,82 @@ +// +// Created by Wesley Ketchum on 4/27/24. +// + +#ifndef SIM_CORE_EventWeights_H +#define SIM_CORE_EventWeights_H + +#include "TObject.h" +#include +#include + +namespace ldmx { + + class EventWeights { + public: + + //Enum to define + enum VariationType { + kINVALID = -1, + kUNKNOWN = 0, + kGENIE_GENERIC = 1000 + }; + + EventWeights() = default; + virtual ~EventWeights() = default; + + //constructor where variation map is declared + EventWeights( std::map< VariationType, std::vector > var_map ) + : variations_map_(var_map) {} + + void Clear(); + void Print(); + + std::vector getWeights() const { return weights_; } + std::map< VariationType, std::vector > 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< VariationType, double > 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< VariationType, std::vector > 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"; + } + 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; + + return VariationType::kUNKNOWN; + } + + }; +} + +#endif //SIM_CORE_EventWeights_H diff --git a/SimCore/include/SimCore/Event/GHepParticle.h b/SimCore/include/SimCore/Event/GHepParticle.h deleted file mode 100644 index b15606d9e..000000000 --- a/SimCore/include/SimCore/Event/GHepParticle.h +++ /dev/null @@ -1,62 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -/// \file GHepParticle.h -/// \Class to hold the particle information needed to recreate a genie::EventRecord -/// \author wketchum@fnal.gov -/// -/// See https://github.com/GENIE-MC/Generator/blob/master/src/Framework/GHEP/GHepParticle.h -/// -/// Adopted for LDMX-SW by wketchum@fnal.gov -//////////////////////////////////////////////////////////////////////// - -#ifndef SIMCORE_EVENT_GHEPPARTICLE_H -#define SIMCORE_EVENT_GHEPPARTICLE_H - -#include "TObject.h" -#include - -namespace ldmx { - - class GHepParticle { - - public: - - GHepParticle() {} - - void Clear() {} - void Print() const; - - int fPosition; - - int fPdgCode; - int fStatus; - int fRescatterCode; - int fFirstMother; - int fLastMother; - int fFirstDaugher; - int fLastDaughter; - - double fP_x; - double fP_y; - double fP_z; - double fP_t; - - double fX_x; - double fX_y; - double fX_z; - double fX_t; - - double fPolzTheta; - double fPolzPhi; - - double fRemovalEnergy; - bool fIsBound; - - public: - - ClassDef(GHepParticle, 1); - - }; - -} // end simb namespace - -#endif // SIMCORE_EVENT_GTRUTH_H diff --git a/SimCore/include/SimCore/Event/GTruth.h b/SimCore/include/SimCore/Event/GTruth.h deleted file mode 100644 index aa63e9f57..000000000 --- a/SimCore/include/SimCore/Event/GTruth.h +++ /dev/null @@ -1,154 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -/// \file GTruth.h -/// \Class to hold the additional information needed to recreate a genie::EventRecord -/// \author nathan.mayer@tufts.edu -/// -/// See https://github.com/NuSoftHEP/nusimdata/blob/develop/nusimdata/SimulationBase/GTruth.h -/// -/// Adopted for LDMX-SW by wketchum@fnal.gov -//////////////////////////////////////////////////////////////////////// - -/// This class stores/retrieves the additional information needed -/// (and not in MCTruth) to recreate a genie::EventRecord -/// for genie based event reweighting. - -#ifndef SIMCORE_EVENT__GTRUTH_H -#define SIMCORE_EVENT__GTRUTH_H - -#include "TObject.h" - -#include -#include - -namespace ldmx { - - class GTruth { - - public: - - GTruth() {} - - void Clear() {} - void Print() const; - - // genie::GHepRecord info - // holds a genie::Interaction - - //TLorentzVector fVertex; - double fVertex_x; - double fVertex_y; - double fVertex_z; - double fVertex_t; - - - // skipping TBits data members EventFlags and EventMask - double fweight; ///< event interaction weight (genie internal) - double fprobability; ///< interaction probability - double fXsec; ///< cross section of interaction - double fDiffXsec; ///< differential cross section of interaction - int fGPhaseSpace; ///< phase space system of DiffXSec - - // genie::Interaction - // container for InitialState, ProcessInfo, XclsTag, Kinematics, KPhaseSpace) - // holds no fundamental type info - - // genie:::InitialState info (sub-object to genie::Interactions) - int fProbePDG; - // holds a genie::Target - - //TLorentzVector fProbeP4; - double fProbe_px; - double fProbe_py; - double fProbe_pz; - double fProbe_e; - - //TLorentzVector fTgtP4; // added version 13 - double fTgt_px; - double fTgt_py; - double fTgt_pz; - double fTgt_e; - - // genie::Target info (sub-object to genie::InitialState) - int ftgtZ; - int ftgtA; - int ftgtPDG; ///< PDG of Target Nucleus, nucleon only if free - int fHitNucPDG; ///< hit nucleon PDG code // added version 13 - int fHitQrkPDG; ///< hit quark PDG code // added version 13 - bool fIsSeaQuark; - //TLorentzVector fHitNucP4; - double fHitNuc_px; - double fHitNuc_py; - double fHitNuc_pz; - double fHitNuc_e; - - double fHitNucPos; // added version 12 - - // genie::ProcessInfo (sub-object to genie::Interactions) - int fGscatter; ///< neutrino scattering code - int fGint; ///< interaction code - - // genie::Kinematics info (sub-object to genie::Interactions) - - ///< these are for the internal (on shell) genie kinematics - ///< this list might be an incomplete transcription of map - //double fgQ2; - //double fgq2; - //double fgW; - //double fgT; - //double fgX; - //double fgY; - std::map fKV; - - ///< a common running variable to be recorded - //double fgWrun; - - //TLorentzVector fFSleptonP4; ///< generated final state primary lepton (LAB frame) // added version 13 - double fFSlepton_px; - double fFSlepton_py; - double fFSlepton_pz; - double fFSlepton_e; - - //TLorentzVector fFShadSystP4; ///< generated final state hadronic system (LAB frame) - double fFShadSyst_px; - double fFShadSyst_py; - double fFShadSyst_pz; - double fFShadSyst_e; - - // genie::XclsTag info (sub-object to genie::Interactions) - bool fIsCharm; ///< did the interaction produce a charmed hadron? - int fCharmHadronPdg; // added version 13 - bool fIsStrange; ///< strange production // added version 13 - int fStrangeHadronPdg; // added version 13 - int fNumProton; ///< number of protons after reaction, before FSI - int fNumNeutron; ///< number of neutrons after reaction, before FSI - int fNumPi0; ///< number of pi0 after reaction, before FSI - int fNumPiPlus; ///< number of pi pluses after reaction, before FSI - int fNumPiMinus; ///< number of pi minuses after reaction, before FSI - int fNumSingleGammas; ///< number of gammas after reaction, before FSI - int fNumRho0; ///< number of pi0 after reaction, before FSI - int fNumRhoPlus; ///< number of pi pluses after reaction, before FSI - int fNumRhoMinus; ///< number of pi minuses after reaction, before FSI - int fResNum; ///< resonance number - int fDecayMode; // added version 13 - bool fIsFinalQuarkEvent; - int fFinalQuarkPdg; - bool fIsFinalLeptonEvent; - int fFinalLeptonPdg; - - - // genie::KPhaseSpace (sub-object to genie::Interactions) - // has no relevant private data - - // Flag for values that might not have been set - static constexpr double kUndefinedValue = -99999; - - public: - //friend std::ostream& operator<< (std::ostream& output, const simb::GTruth >ruth); - - ClassDef(GTruth, 1); - - }; - -} // end simb namespace - -#endif // SIMCORE_EVENT__GTRUTH_H diff --git a/SimCore/include/SimCore/Event/HepMC3GenEvent.h b/SimCore/include/SimCore/Event/HepMC3GenEvent.h new file mode 100644 index 000000000..630cdba09 --- /dev/null +++ b/SimCore/include/SimCore/Event/HepMC3GenEvent.h @@ -0,0 +1,33 @@ +// +// 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 "TObject.h" +#include "HepMC3/Data/GenEventData.h" + +namespace ldmx { + + class HepMC3GenEvent : public HepMC3::GenEventData { + + public: + + //HepMC3GenEvent() : HepMC3::GenEventData() {} + //HepMC3GenEvent(const HepMC3::GenEvent& event) : HepMC3::GenEvent(event) {} + + void Clear(); + void Print() const; + + public: + ClassDef(HepMC3GenEvent, 1); + + }; + +} + +#endif //SIM_CORE_HEPMC3GENEVENT_H diff --git a/SimCore/include/SimCore/Generators/GenieGenerator.h b/SimCore/include/SimCore/Generators/GenieGenerator.h index 8bd928ceb..bb72df7ca 100644 --- a/SimCore/include/SimCore/Generators/GenieGenerator.h +++ b/SimCore/include/SimCore/Generators/GenieGenerator.h @@ -18,6 +18,7 @@ // GENIE // //------------// #include "GENIE/Framework/EventGen/GEVGDriver.h" +#include "GENIE/Framework/EventGen/HepMC3Converter.h" //------------// // LDMX // @@ -71,9 +72,10 @@ class GenieGenerator : public simcore::PrimaryGenerator { private: /** - * The GENIE event generator driver + * The GENIE event generator driver and convertor to HepMC3GenEvent */ genie::GEVGDriver evg_driver_; + genie::HepMC3Converter hepMC3Converter_; int verbosity_; double energy_; diff --git a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h new file mode 100644 index 000000000..9ef6d90f4 --- /dev/null +++ b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h @@ -0,0 +1,70 @@ +// +// Created by Wesley Ketchum on 4/29/24. +// + +#ifndef SIMCORE_GENIEREWEIGHTPRODUCER_H +#define SIMCORE_GENIEREWEIGHTPRODUCER_H + +#include "SimCore/Event/EventWeights.h" + +#include "Framework/EventProcessor.h" + +#include +#include + +namespace genie { + class Interaction; + class HepMC3Converter; +} + +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&); + + //produce on the event + virtual void produce(framework::Event& event); + + private: + + //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_; + + //variations to run + std::map< ldmx::EventWeights::VariationType, std::vector > variation_map_; + + //hepmc3 convertor + genie::HepMC3Converter* hepMC3Converter_; + + genie::rew::GReWeight* genie_rw_; + + }; + +} + +#endif //SIMCORE_GENIEREWEIGHTPRODUCER_H diff --git a/SimCore/include/SimCore/Simulator.h b/SimCore/include/SimCore/Simulator.h index fc09f9938..5b713d7c3 100644 --- a/SimCore/include/SimCore/Simulator.h +++ b/SimCore/include/SimCore/Simulator.h @@ -25,12 +25,13 @@ #include "SimCore/DetectorConstruction.h" #include "SimCore/RunManager.h" #include "SimCore/SimulatorBase.h" -#include "SimCore/Event/GTruth.h" -#include "SimCore/Event/GHepParticle.h" namespace genie { - class EventRecord; - class GHepParticle; + class Interaction; +} + +namespace ldmx { + class HepMC3GenEvent; } class G4UImanager; @@ -115,12 +116,6 @@ class Simulator : public SimulatorBase { */ void setSeeds(std::vector seeds); - void FillGTruth(const genie::EventRecord* record, - ldmx::GTruth& truth); - - void FillGHepParticles(const genie::EventRecord* record, - std::vector& ghep_particle_list); - private: /// Number of events started int numEventsBegan_{0}; diff --git a/SimCore/include/SimCore/UserEventInformation.h b/SimCore/include/SimCore/UserEventInformation.h index 9d38bfb78..09c5316da 100644 --- a/SimCore/include/SimCore/UserEventInformation.h +++ b/SimCore/include/SimCore/UserEventInformation.h @@ -3,7 +3,8 @@ #include "G4VUserEventInformation.hh" -#include "GENIE/Framework/EventGen/EventRecord.h" +#include "SimCore/Event/HepMC3GenEvent.h" +#include namespace simcore { @@ -16,8 +17,7 @@ class UserEventInformation : public G4VUserEventInformation { UserEventInformation() = default; /// Destructor - virtual ~UserEventInformation() - { if(genie_event_) delete genie_event_;} + virtual ~UserEventInformation() {} /// Print the information associated with the track void Print() const final override; @@ -121,9 +121,9 @@ class UserEventInformation : public G4VUserEventInformation { */ bool wasLastStepEN() const { return last_step_en_; } - void setGENIEEventRecord(genie::EventRecord *event) { genie_event_ = event; } - genie::EventRecord* getGENIEEventRecord() { return genie_event_; } - + 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}; @@ -174,9 +174,9 @@ class UserEventInformation : public G4VUserEventInformation { double db_material_z_{-1.}; /** - * Pointer for a GENIE event record. + * a collection of HepMC3 event records. */ - genie::EventRecord* genie_event_{NULL}; + std::vector< ldmx::HepMC3GenEvent > hepmc3_events_; }; } // namespace simcore diff --git a/SimCore/python/genie_reweight.py b/SimCore/python/genie_reweight.py new file mode 100644 index 000000000..0982a7ca5 --- /dev/null +++ b/SimCore/python/genie_reweight.py @@ -0,0 +1,12 @@ +from LDMX.Framework import ldmxcfg + +class GenieReweightProducer(ldmxcfg.Producer) : + + def __init__(self,name='genieEventWeights'): + super().__init__(name,"simcore::GenieReweightProducer","SimCore::Reweight") + + 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..c05e6c535 --- /dev/null +++ b/SimCore/src/SimCore/Event/EventWeights.cxx @@ -0,0 +1,25 @@ +// +// 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 - -ClassImp(ldmx::GHepParticle) - -namespace ldmx{ - - void GHepParticle::Print() const { - std::cout << "GHEPPARTICLE! TODO!" << std::endl ; - - } - -} diff --git a/SimCore/src/SimCore/Event/GTruth.cxx b/SimCore/src/SimCore/Event/GTruth.cxx deleted file mode 100644 index 6c996f447..000000000 --- a/SimCore/src/SimCore/Event/GTruth.cxx +++ /dev/null @@ -1,14 +0,0 @@ -#include "SimCore/Event/GTruth.h" - -#include - -ClassImp(ldmx::GTruth) - -namespace ldmx{ - - void GTruth::Print() const { - std::cout << "GTRUTH! TODO!" << std::endl ; - - } - -} diff --git a/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx b/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx new file mode 100644 index 000000000..d08cb53a9 --- /dev/null +++ b/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx @@ -0,0 +1,26 @@ +// +// Created by Wesley Ketchum on 4/26/24. +// + +#include "SimCore/Event/HepMC3GenEvent.h" +#include "HepMC3/GenEvent.h" +#include "HepMC3/Print.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 + } + +} diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx index 2b85cb870..eff8416d5 100644 --- a/SimCore/src/SimCore/Generators/GenieGenerator.cxx +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -6,7 +6,6 @@ #include "SimCore/Generators/GenieGenerator.h" #include "SimCore/UserPrimaryParticleInformation.h" -#include "SimCore/Event/GTruth.h" // GENIE #include "Framework/Utils/AppInit.h" @@ -15,7 +14,7 @@ #include "GENIE/Framework/Interaction/InitialState.h" #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepParticle.h" - +#include "HepMC3/GenEvent.h" #include "Framework/Conventions/XmlParserStatus.h" #include "Framework/Conventions/GBuild.h" @@ -383,9 +382,12 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* event) while(!genie_event) genie_event = evg_driver_.GenerateEvent(e_p4); - auto genie_info = new UserEventInformation; - genie_info->setGENIEEventRecord(genie_event); - event->SetUserInformation(genie_info); + 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 diff --git a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx new file mode 100644 index 000000000..c1a65f791 --- /dev/null +++ b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx @@ -0,0 +1,135 @@ +// +// Created by Wesley Ketchum on 4/29/24. +// + +#include "SimCore/Reweight/GenieReweightProducer.h" +#include "SimCore/Event/HepMC3GenEvent.h" +#include "SimCore/Event/EventWeights.h" + +#include "Framework/EventGen/HepMC3Converter.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/Interaction/Interaction.h" + +#include "RwFramework/GReWeightI.h" +#include "RwFramework/GSystSet.h" +#include "RwFramework/GSyst.h" +#include "RwFramework/GReWeight.h" +#include "RwCalculators/GReWeightNuXSecNCEL.h" +#include "RwCalculators/GReWeightNuXSecCCQE.h" +#include "RwCalculators/GReWeightNuXSecCCRES.h" +#include "RwCalculators/GReWeightNuXSecCOH.h" +#include "RwCalculators/GReWeightNonResonanceBkg.h" +#include "RwCalculators/GReWeightFGM.h" +#include "RwCalculators/GReWeightDISNuclMod.h" +#include "RwCalculators/GReWeightResonanceDecay.h" +#include "RwCalculators/GReWeightFZone.h" +#include "RwCalculators/GReWeightINuke.h" +#include "RwCalculators/GReWeightAGKY.h" +#include "RwCalculators/GReWeightNuXSecCCQEaxial.h" +#include "RwCalculators/GReWeightNuXSecCCQEvec.h" +#include "RwCalculators/GReWeightNuXSecNCRES.h" +#include "RwCalculators/GReWeightNuXSecDIS.h" + +#include "RwCalculators/GReWeightINukeParams.h" +#include "RwCalculators/GReWeightNuXSecNC.h" +#include "RwCalculators/GReWeightXSecEmpiricalMEC.h" +#include "RwCalculators/GReWeightXSecMEC.h" +#include "RwCalculators/GReWeightDeltaradAngle.h" + + +#include "HepMC3/GenEvent.h" +#include "HepMC3/Print.h" + +#include +#include + + + +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"); + seed_ = ps.getParameter("seed"); + n_weights_ = (size_t)(ps.getParameter("n_weights")); + auto var_types_strings = ps.getParameter< std::vector >("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_wAdoptWghtCalc( "hadro_fzone", new genie::rew::GReWeightFZone ); + genie_rw_->AdoptWghtCalc( "hadro_intranuke", new genie::rew::GReWeightINuke ); + auto & syst = genie_rw_->Systematics(); + syst.Init(genie::rew::GSyst::FromString("FrInel_pi")); + } + + void GenieReweightProducer::produce(framework::Event &event) { + + //grab the input hepmc3 event collection + auto hepmc3_col = event.getObject< std::vector >(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; + + auto & syst = genie_rw_->Systematics(); + auto var_value = variation_map_[ldmx::EventWeights::kGENIE_GENERIC][i_w]; + syst.Set(genie::rew::GSyst::FromString("FrInel_pi"), var_value); + genie_rw_->Reconfigure(); + + //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) 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) 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 + + ev_weights.Print(); + + event.add(eventWeightsCollName_, ev_weights); + + } +} + +DECLARE_PRODUCER_NS(simcore, GenieReweightProducer); diff --git a/SimCore/src/SimCore/Simulator.cxx b/SimCore/src/SimCore/Simulator.cxx index e1f76ac38..302c52032 100644 --- a/SimCore/src/SimCore/Simulator.cxx +++ b/SimCore/src/SimCore/Simulator.cxx @@ -16,11 +16,6 @@ #include "Framework/Version.h" //for LDMX_INSTALL path -/*~~~~~~~~~~~~~*/ -/* GENIE */ -/*~~~~~~~~~~~~~*/ -#include "GENIE/Framework/GHEP/GHepParticle.h" - /*~~~~~~~~~~~~~*/ /* SimCore */ /*~~~~~~~~~~~~~*/ @@ -33,8 +28,7 @@ #include "SimCore/SensitiveDetector.h" #include "SimCore/UserEventInformation.h" #include "SimCore/XsecBiasingOperator.h" -#include "SimCore/Event/GTruth.h" -#include "SimCore/Event/GHepParticle.h" +#include "SimCore/Event/HepMC3GenEvent.h" /*~~~~~~~~~~~~~~*/ /* Geant4 */ @@ -136,141 +130,6 @@ void Simulator::onNewRun(const ldmx::RunHeader& runHeader) { run_ = runHeader.getRunNumber(); } -void Simulator::FillGHepParticles(const genie::EventRecord* record, - std::vector& ghep_particle_list) -{ - for(int position=0; position< record->GetEntries(); ++position) - { - const genie::GHepParticle* particle = record->Particle(position); - - ldmx::GHepParticle my_particle; - - my_particle.fPosition = position; - - my_particle.fPdgCode = particle->Pdg(); - my_particle.fStatus = particle->Status(); - my_particle.fRescatterCode = particle->RescatterCode(); - my_particle.fFirstMother = particle->FirstMother(); - my_particle.fLastMother = particle->LastMother(); - my_particle.fFirstDaugher = particle->FirstDaughter(); - my_particle.fLastDaughter = particle->LastDaughter(); - - my_particle.fP_x = particle->P4()->Px(); - my_particle.fP_y = particle->P4()->Py(); - my_particle.fP_z = particle->P4()->Pz(); - my_particle.fP_t = particle->P4()->E(); - - my_particle.fX_x = particle->X4()->X(); - my_particle.fX_y = particle->X4()->Y(); - my_particle.fX_z = particle->X4()->Z(); - my_particle.fX_t = particle->X4()->T(); - - my_particle.fPolzTheta = particle->PolzPolarAngle(); - my_particle.fPolzPhi = particle->PolzAzimuthAngle(); - - my_particle.fRemovalEnergy = particle->RemovalEnergy(); - my_particle.fIsBound = particle->IsBound(); - - ghep_particle_list.push_back(my_particle); - } -} - -void Simulator::FillGTruth(const genie::EventRecord* record, - ldmx::GTruth& truth) { - - //interactions info - genie::Interaction *inter = record->Summary(); - const genie::ProcessInfo &procInfo = inter->ProcInfo(); - truth.fGint = (int)procInfo.InteractionTypeId(); - truth.fGscatter = (int)procInfo.ScatteringTypeId(); - - //Event info - truth.fweight = record->Weight(); - truth.fprobability = record->Probability(); - truth.fXsec = record->XSec(); - truth.fDiffXsec = record->DiffXSec(); - truth.fGPhaseSpace = (int)record->DiffXSecVars(); - - //Initial State info - const genie::InitialState &initState = inter->InitState(); - truth.fProbePDG = initState.ProbePdg(); - truth.fProbe_px = initState.GetProbeP4()->Px(); - truth.fProbe_py = initState.GetProbeP4()->Py(); - truth.fProbe_pz = initState.GetProbeP4()->Pz(); - truth.fProbe_e = initState.GetProbeP4()->E(); - - truth.fTgt_px = initState.GetTgtP4()->Px(); - truth.fTgt_py = initState.GetTgtP4()->Py(); - truth.fTgt_pz = initState.GetTgtP4()->Pz(); - truth.fTgt_e = initState.GetTgtP4()->E(); - - //Target info - const genie::Target &tgt = initState.Tgt(); - truth.ftgtZ = tgt.Z(); - truth.ftgtA = tgt.A(); - truth.ftgtPDG = tgt.Pdg(); - truth.fHitNucPDG = tgt.HitNucPdg(); - truth.fHitQrkPDG = tgt.HitQrkPdg(); - truth.fIsSeaQuark = tgt.HitSeaQrk(); - truth.fHitNuc_px = tgt.HitNucP4Ptr()->Px(); - truth.fHitNuc_py = tgt.HitNucP4Ptr()->Py(); - truth.fHitNuc_pz = tgt.HitNucP4Ptr()->Pz(); - truth.fHitNuc_e = tgt.HitNucP4Ptr()->E(); - truth.fHitNucPos = tgt.HitNucPosition(); - - - truth.fVertex_x = record->Vertex()->X(); - truth.fVertex_y = record->Vertex()->Y(); - truth.fVertex_z = record->Vertex()->Z(); - truth.fVertex_t = record->Vertex()->T(); - - //true reaction information and byproducts - //(PRE FSI) - const genie::XclsTag &exclTag = inter->ExclTag(); - truth.fIsCharm = exclTag.IsCharmEvent(); - truth.fCharmHadronPdg = exclTag.CharmHadronPdg(); - truth.fIsStrange = exclTag.IsStrangeEvent(); - truth.fStrangeHadronPdg = exclTag.StrangeHadronPdg(); - truth.fResNum = (int)exclTag.Resonance(); - truth.fDecayMode = exclTag.DecayMode(); - - truth.fNumPiPlus = exclTag.NPiPlus(); - truth.fNumPiMinus = exclTag.NPiMinus(); - truth.fNumPi0 = exclTag.NPi0(); - truth.fNumProton = exclTag.NProtons(); - truth.fNumNeutron = exclTag.NNeutrons(); - truth.fNumSingleGammas = exclTag.NSingleGammas(); - - truth.fNumSingleGammas = exclTag.NSingleGammas(); - truth.fNumRho0 = exclTag.NRho0(); - truth.fNumRhoPlus = exclTag.NRhoPlus(); - truth.fNumRhoMinus = exclTag.NRhoMinus(); - - truth.fIsFinalQuarkEvent = exclTag.IsFinalQuarkEvent(); - truth.fFinalQuarkPdg = exclTag.FinalQuarkPdg(); - truth.fIsFinalLeptonEvent = exclTag.IsFinalLeptonEvent(); - truth.fFinalLeptonPdg = exclTag.FinalLeptonPdg(); - - // Get the GENIE kinematics info - const genie::Kinematics &kine = inter->Kine(); - for(int kvar=genie::KineVar_t::kKVNull; kvar!=genie::KineVar_t::kNumOfKineVar; ++kvar) - truth.fKV[kvar] = kine.GetKV((genie::KineVar_t)kvar); - - truth.fFSlepton_px = kine.FSLeptonP4().Px(); - truth.fFSlepton_py = kine.FSLeptonP4().Py(); - truth.fFSlepton_pz = kine.FSLeptonP4().Pz(); - truth.fFSlepton_e = kine.FSLeptonP4().E(); - - truth.fFShadSyst_px = kine.HadSystP4().Px(); - truth.fFShadSyst_py = kine.HadSystP4().Py(); - truth.fFShadSyst_pz = kine.HadSystP4().Pz(); - truth.fFShadSyst_e = kine.HadSystP4().E(); - - return; - -} - - void Simulator::produce(framework::Event& event) { // Generate and process a Geant4 event. numEventsBegan_++; @@ -308,18 +167,13 @@ void Simulator::produce(framework::Event& event) { auto event_info = static_cast( runManager_->GetCurrentEvent()->GetUserInformation()); - if(event_info->getGENIEEventRecord()){ - //event_info->getGENIEEventRecord()->Print(); - - ldmx::GTruth gtruth; - std::vector ghep_particles; - FillGTruth(event_info->getGENIEEventRecord(),gtruth); - FillGHepParticles(event_info->getGENIEEventRecord(),ghep_particles); - - event.add("SimGTruth", gtruth); - event.add("SimGHepParticles", ghep_particles); + 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); From 4feda33fac9fb288364e037790c999ce5e5cf522 Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Wed, 17 Jul 2024 08:35:12 +0200 Subject: [PATCH 04/14] reenable GENIE Reweight libraries --- cmake/BuildMacros.cmake | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cmake/BuildMacros.cmake b/cmake/BuildMacros.cmake index 023d57e76..fb7102ba6 100644 --- a/cmake/BuildMacros.cmake +++ b/cmake/BuildMacros.cmake @@ -213,9 +213,9 @@ SET( GENIE_LIBS Genie::GPhHEDISEG Genie::GTlGeo Genie::GTlFlx -# Genie::GRwFwk -# Genie::GRwClc -# Genie::GRwIO + Genie::GRwFwk + Genie::GRwClc + Genie::GRwIO blas) message(STATUS "Setting GENIE_LIBS to ${GENIE_LIBS}") endif() From 5d4039c722f5d8e556e3295dfc7916c95060b4f4 Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Thu, 18 Jul 2024 13:18:16 +0200 Subject: [PATCH 05/14] fix the linking orders, and get the reweight producer to be working --- SimCore/CMakeLists.txt | 6 +++--- .../SimCore/Reweight/GenieReweightProducer.h | 3 +++ SimCore/python/genie_reweight.py | 1 + SimCore/src/SimCore/Generators/GenieGenerator.cxx | 1 + .../src/SimCore/Reweight/GenieReweightProducer.cxx | 9 +++++++++ cmake/BuildMacros.cmake | 14 ++++++++------ 6 files changed, 25 insertions(+), 9 deletions(-) diff --git a/SimCore/CMakeLists.txt b/SimCore/CMakeLists.txt index b25211661..6460a1489 100644 --- a/SimCore/CMakeLists.txt +++ b/SimCore/CMakeLists.txt @@ -77,8 +77,8 @@ 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 - ${GENIE_LIBS} Framework::Configure Framework::Framework SimCore::G4User @@ -97,12 +97,12 @@ setup_library(module SimCore # Primary Generators library setup_library(module SimCore name Generators - dependencies SimCore::SimCore SimCore::LHE ROOT::MathMore ${GENIE_LIBS} HepMC3) + dependencies SimCore::SimCore SimCore::LHE SimCore::Event ${GENIE_LIBS} HepMC3) # Reweight Library setup_library(module SimCore name Reweight - dependencies SimCore::SimCore SimCore::Generators Framework::Framework ${GENIE_LIBS} HepMC3) + 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/Reweight/GenieReweightProducer.h b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h index 9ef6d90f4..c42805a76 100644 --- a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h +++ b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h @@ -55,6 +55,9 @@ namespace simcore { // number of weights to be calculated per event size_t n_weights_; + // GENIE tune + std::string tune_; + //variations to run std::map< ldmx::EventWeights::VariationType, std::vector > variation_map_; diff --git a/SimCore/python/genie_reweight.py b/SimCore/python/genie_reweight.py index 0982a7ca5..5ea739f21 100644 --- a/SimCore/python/genie_reweight.py +++ b/SimCore/python/genie_reweight.py @@ -8,5 +8,6 @@ def __init__(self,name='genieEventWeights'): self.seed = 10 self.n_weights = 100 self.var_types = ["GENIE_GENERIC"] + self.tune = "G18_02a_02_11b" self.eventWeightsCollName = "genieEventWeights" \ No newline at end of file diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx index eff8416d5..2f3677009 100644 --- a/SimCore/src/SimCore/Generators/GenieGenerator.cxx +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -52,6 +52,7 @@ #include #include #include +#include "Math/Interpolator.h" // standard #include diff --git a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx index c1a65f791..84b3dcdc6 100644 --- a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx +++ b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx @@ -9,6 +9,7 @@ #include "Framework/EventGen/HepMC3Converter.h" #include "Framework/EventGen/EventRecord.h" #include "Framework/Interaction/Interaction.h" +#include "Framework/Utils/RunOpt.h" #include "RwFramework/GReWeightI.h" #include "RwFramework/GSystSet.h" @@ -64,6 +65,7 @@ namespace simcore { eventWeightsCollName_ = ps.getParameter("eventWeightsCollName"); seed_ = ps.getParameter("seed"); n_weights_ = (size_t)(ps.getParameter("n_weights")); + tune_ = ps.getParameter("tune"); auto var_types_strings = ps.getParameter< std::vector >("var_types"); std::default_random_engine generator(seed_); @@ -74,6 +76,13 @@ namespace simcore { for(size_t i_w=0; i_wSetTuneName(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(); diff --git a/cmake/BuildMacros.cmake b/cmake/BuildMacros.cmake index fb7102ba6..347107b30 100644 --- a/cmake/BuildMacros.cmake +++ b/cmake/BuildMacros.cmake @@ -159,9 +159,12 @@ SET( GENIE_LIBS # Genie::GFwMsg # Genie::GFwReg # Genie::GFwAlg - ROOT::MathCore - ROOT::MathMore - Genie::GFwInt + + Genie::GRwFwk + Genie::GRwClc + Genie::GRwIO + + Genie::GFwInt # Genie::GFwGHEP Genie::GFwNum Genie::GFwUtl @@ -213,9 +216,8 @@ SET( GENIE_LIBS Genie::GPhHEDISEG Genie::GTlGeo Genie::GTlFlx - Genie::GRwFwk - Genie::GRwClc - Genie::GRwIO + ROOT::MathCore + ROOT::MathMore blas) message(STATUS "Setting GENIE_LIBS to ${GENIE_LIBS}") endif() From b2761b349a46469264680326693dd63b6b9b0916 Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Fri, 19 Jul 2024 12:22:46 +0200 Subject: [PATCH 06/14] updates to put the GENIE tune information into the run header, and extract that in the GenieReweightProducer --- .../SimCore/Reweight/GenieReweightProducer.h | 5 ++ SimCore/python/genie_reweight.py | 1 - .../src/SimCore/Generators/GenieGenerator.cxx | 1 + .../Reweight/GenieReweightProducer.cxx | 46 ++++++++++++++----- 4 files changed, 41 insertions(+), 12 deletions(-) diff --git a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h index c42805a76..df713f0a9 100644 --- a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h +++ b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h @@ -35,6 +35,9 @@ namespace simcore { //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); @@ -66,6 +69,8 @@ namespace simcore { genie::rew::GReWeight* genie_rw_; + void reconfigureGenieReweight(); + }; } diff --git a/SimCore/python/genie_reweight.py b/SimCore/python/genie_reweight.py index 5ea739f21..0982a7ca5 100644 --- a/SimCore/python/genie_reweight.py +++ b/SimCore/python/genie_reweight.py @@ -8,6 +8,5 @@ def __init__(self,name='genieEventWeights'): self.seed = 10 self.n_weights = 100 self.var_types = ["GENIE_GENERIC"] - self.tune = "G18_02a_02_11b" self.eventWeightsCollName = "genieEventWeights" \ No newline at end of file diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx index 2f3677009..4bd80e789 100644 --- a/SimCore/src/SimCore/Generators/GenieGenerator.cxx +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -442,6 +442,7 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* 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 diff --git a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx index 84b3dcdc6..cd34ef71c 100644 --- a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx +++ b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx @@ -37,7 +37,6 @@ #include "RwCalculators/GReWeightXSecMEC.h" #include "RwCalculators/GReWeightDeltaradAngle.h" - #include "HepMC3/GenEvent.h" #include "HepMC3/Print.h" @@ -65,7 +64,6 @@ namespace simcore { eventWeightsCollName_ = ps.getParameter("eventWeightsCollName"); seed_ = ps.getParameter("seed"); n_weights_ = (size_t)(ps.getParameter("n_weights")); - tune_ = ps.getParameter("tune"); auto var_types_strings = ps.getParameter< std::vector >("var_types"); std::default_random_engine generator(seed_); @@ -77,16 +75,42 @@ namespace simcore { variation_map_[vtype].push_back(normal_distribution(generator)); } - genie::RunOpt::Instance()->SetTuneName(tune_); - if ( ! genie::RunOpt::Instance()->Tune() ) { - EXCEPTION_RAISE("ConfigurationException","No TuneId in RunOption."); + } + + void GenieReweightProducer::reconfigureGenieReweight() + { + 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(); + syst.Init(genie::rew::GSyst::FromString("FrInel_pi")); + } + + 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; } - 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(); - syst.Init(genie::rew::GSyst::FromString("FrInel_pi")); + if (tune_ != new_tune) { + std::cout << "Found new tune " << new_tune << " (used to be " << tune_ << ")" << std::endl; + tune_ = new_tune; + reconfigureGenieReweight(); + } } void GenieReweightProducer::produce(framework::Event &event) { @@ -134,7 +158,7 @@ namespace simcore { } //end loop over weights - ev_weights.Print(); + //ev_weights.Print(); event.add(eventWeightsCollName_, ev_weights); From 1e085f242fa0446e1ce93a4406087a001dd10205 Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Mon, 22 Jul 2024 07:43:54 +0200 Subject: [PATCH 07/14] flesh out the genie reweighting to take in FSI calculators --- SimCore/include/SimCore/Event/EventWeights.h | 73 ++++++++++++++++++- .../SimCore/Reweight/GenieReweightProducer.h | 39 +++++++++- SimCore/python/genie_reweight.py | 1 + .../Reweight/GenieReweightProducer.cxx | 52 ++++++++++--- 4 files changed, 151 insertions(+), 14 deletions(-) diff --git a/SimCore/include/SimCore/Event/EventWeights.h b/SimCore/include/SimCore/Event/EventWeights.h index 9cc2ea122..e4e5687df 100644 --- a/SimCore/include/SimCore/Event/EventWeights.h +++ b/SimCore/include/SimCore/Event/EventWeights.h @@ -18,7 +18,22 @@ namespace ldmx { enum VariationType { kINVALID = -1, kUNKNOWN = 0, - kGENIE_GENERIC = 1000 + + 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; @@ -57,10 +72,36 @@ namespace ldmx { switch (vtype) { case VariationType::kINVALID: return "INVALID"; - case VariationType ::kUNKNOWN: + case VariationType::kUNKNOWN: return "UNKNOWN"; - case VariationType ::kGENIE_GENERIC: + 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"; } @@ -72,6 +113,32 @@ namespace ldmx { 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; } diff --git a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h index df713f0a9..855f4a9fb 100644 --- a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h +++ b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h @@ -9,6 +9,8 @@ #include "Framework/EventProcessor.h" +#include "RwFramework/GSyst.h" + #include #include @@ -24,6 +26,7 @@ namespace genie::rew { namespace simcore { class GenieReweightProducer : public framework::Producer { + public: //constructor @@ -43,6 +46,9 @@ namespace simcore { private: + //seed to use + int verbosity_; + //input hepmc3 collection name std::string hepmc3CollName_; @@ -69,10 +75,39 @@ namespace simcore { genie::rew::GReWeight* genie_rw_; - void reconfigureGenieReweight(); + 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; + } + return genie::rew::EGSyst::kNullSystematic; + } }; - } #endif //SIMCORE_GENIEREWEIGHTPRODUCER_H diff --git a/SimCore/python/genie_reweight.py b/SimCore/python/genie_reweight.py index 0982a7ca5..9703c8358 100644 --- a/SimCore/python/genie_reweight.py +++ b/SimCore/python/genie_reweight.py @@ -5,6 +5,7 @@ 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"] diff --git a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx index cd34ef71c..cf4951496 100644 --- a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx +++ b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx @@ -62,6 +62,7 @@ namespace simcore { 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< std::vector >("var_types"); @@ -77,7 +78,7 @@ namespace simcore { } - void GenieReweightProducer::reconfigureGenieReweight() + void GenieReweightProducer::reinitializeGenieReweight() { delete genie_rw_; genie_rw_ = new genie::rew::GReWeight; @@ -90,8 +91,10 @@ namespace simcore { genie_rw_->AdoptWghtCalc( "hadro_fzone", new genie::rew::GReWeightFZone ); genie_rw_->AdoptWghtCalc( "hadro_intranuke", new genie::rew::GReWeightINuke ); + auto & syst = genie_rw_->Systematics(); - syst.Init(genie::rew::GSyst::FromString("FrInel_pi")); + for (auto var : variation_map_) + syst.Init(variation_type_to_genie_dial(var.first)); } void GenieReweightProducer::onNewRun(const ldmx::RunHeader &runHeader) { @@ -107,10 +110,43 @@ namespace simcore { } if (tune_ != new_tune) { - std::cout << "Found new tune " << new_tune << " (used to be " << tune_ << ")" << std::endl; + if(verbosity_>=0) + std::cout << "Found new tune " << new_tune << " (used to be " << tune_ << ")" << std::endl; tune_ = new_tune; - reconfigureGenieReweight(); + 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) { @@ -125,10 +161,7 @@ namespace simcore { double running_weight=1; - auto & syst = genie_rw_->Systematics(); - auto var_value = variation_map_[ldmx::EventWeights::kGENIE_GENERIC][i_w]; - syst.Set(genie::rew::GSyst::FromString("FrInel_pi"), var_value); - genie_rw_->Reconfigure(); + 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) { @@ -158,7 +191,8 @@ namespace simcore { } //end loop over weights - //ev_weights.Print(); + if(verbosity_>=1) + ev_weights.Print(); event.add(eventWeightsCollName_, ev_weights); From 3a77fe37d689c32f303879f146b4f3052eb0966a Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Thu, 25 Jul 2024 15:09:28 +0200 Subject: [PATCH 08/14] add in a string writer so we can more easily transfer back and forth in python --- .../include/SimCore/Event/HepMC3GenEvent.h | 5 +++++ SimCore/src/SimCore/Event/HepMC3GenEvent.cxx | 20 ++++++++++++++++++- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/SimCore/include/SimCore/Event/HepMC3GenEvent.h b/SimCore/include/SimCore/Event/HepMC3GenEvent.h index 630cdba09..57c99b8cf 100644 --- a/SimCore/include/SimCore/Event/HepMC3GenEvent.h +++ b/SimCore/include/SimCore/Event/HepMC3GenEvent.h @@ -9,6 +9,7 @@ #define SIM_CORE_HEPMC3GENEVENT_H #include "TObject.h" +#include "HepMC3/GenEvent.h" #include "HepMC3/Data/GenEventData.h" namespace ldmx { @@ -23,6 +24,10 @@ namespace ldmx { void Clear(); void Print() const; + HepMC3::GenEvent getHepMCGenEvent() const; + + std::string get_as_string() const; + public: ClassDef(HepMC3GenEvent, 1); diff --git a/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx b/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx index d08cb53a9..91fe53fcd 100644 --- a/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx +++ b/SimCore/src/SimCore/Event/HepMC3GenEvent.cxx @@ -3,8 +3,11 @@ // #include "SimCore/Event/HepMC3GenEvent.h" -#include "HepMC3/GenEvent.h" #include "HepMC3/Print.h" +#include "HepMC3/WriterAscii.h" + +#include +#include namespace ldmx{ @@ -23,4 +26,19 @@ namespace ldmx{ 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(); + + } } From 1c0236788158de8852f6e0088090a6b16a8f7a9e Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Thu, 1 Aug 2024 17:56:10 +0200 Subject: [PATCH 09/14] add a GENIE truth DQM --- DQM/include/DQM/GenieTruthDQM.h | 58 ++++++++++ DQM/python/dqm.py | 20 +++- DQM/src/DQM/GenieTruthDQM.cxx | 190 ++++++++++++++++++++++++++++++++ 3 files changed, 266 insertions(+), 2 deletions(-) create mode 100644 DQM/include/DQM/GenieTruthDQM.h create mode 100644 DQM/src/DQM/GenieTruthDQM.cxx diff --git a/DQM/include/DQM/GenieTruthDQM.h b/DQM/include/DQM/GenieTruthDQM.h new file mode 100644 index 000000000..c4e63dd98 --- /dev/null +++ b/DQM/include/DQM/GenieTruthDQM.h @@ -0,0 +1,58 @@ +// +// 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 caa78936e..0c4ee2599 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -623,7 +623,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(), @@ -672,5 +689,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..d296496cd --- /dev/null +++ b/DQM/src/DQM/GenieTruthDQM.cxx @@ -0,0 +1,190 @@ +// +// Created by Wesley Ketchum on 8/1/24. +// + +#include "DQM/GenieTruthDQM.h" + +#include "SimCore/Event/HepMC3GenEvent.h" + +#include "GENIE/Framework/Conventions/KineVar.h" + +#include + +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< std::vector >(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("hitnuc_px", hitnuc_4vec[0]); + ntuple_.setVar("hitnuc_py", hitnuc_4vec[1]); + ntuple_.setVar("hitnuc_pz", hitnuc_4vec[2]); + ntuple_.setVar("hitnuc_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; + } + + +} + +DECLARE_ANALYZER_NS(dqm, GenieTruthDQM); From 4665e26f3ae449e665ab67cef8231a1faffbc3ca Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Thu, 12 Sep 2024 13:27:56 +0200 Subject: [PATCH 10/14] fix some typos, and decrease verbosity --- DQM/src/DQM/GenieTruthDQM.cxx | 8 ++++---- SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/DQM/src/DQM/GenieTruthDQM.cxx b/DQM/src/DQM/GenieTruthDQM.cxx index d296496cd..3d8b95cea 100644 --- a/DQM/src/DQM/GenieTruthDQM.cxx +++ b/DQM/src/DQM/GenieTruthDQM.cxx @@ -159,10 +159,10 @@ namespace dqm { auto hitnuc_4vec_ptr = hepmc3_ev.attribute("GENIE.Interaction.HitNucleonP4"); if(hitnuc_4vec_ptr) { auto hitnuc_4vec = hitnuc_4vec_ptr->value(); - ntuple_.setVar("hitnuc_px", hitnuc_4vec[0]); - ntuple_.setVar("hitnuc_py", hitnuc_4vec[1]); - ntuple_.setVar("hitnuc_pz", hitnuc_4vec[2]); - ntuple_.setVar("hitnuc_e", hitnuc_4vec[3]); + 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! diff --git a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx index cf4951496..57b7e81aa 100644 --- a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx +++ b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx @@ -172,13 +172,13 @@ namespace simcore { hepmc3_genev.read_data(hepmc3_ev); //print it out to check it ... - if(i_w==0) HepMC3::Print::line(hepmc3_genev, true); //print attributes + 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) genie_ev_record_ptr->Print(std::cout); + 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); From 6a062ccdc1c503a4480dfd01189fdda48f45c1ef Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Wed, 16 Oct 2024 12:39:20 +0000 Subject: [PATCH 11/14] Apply clang-format --- DQM/include/DQM/GenieTruthDQM.h | 3 +- DQM/src/DQM/GenieTruthDQM.cxx | 337 +++++------ SimCore/include/SimCore/Event/EventWeights.h | 278 ++++----- .../include/SimCore/Event/HepMC3GenEvent.h | 33 +- .../SimCore/Generators/GenieGenerator.h | 23 +- SimCore/include/SimCore/PrimaryGenerator.h | 3 +- .../SimCore/Reweight/GenieReweightProducer.h | 185 +++--- SimCore/include/SimCore/Simulator.h | 4 +- .../include/SimCore/UserEventInformation.h | 17 +- SimCore/src/SimCore/Event/EventWeights.cxx | 30 +- SimCore/src/SimCore/Event/HepMC3GenEvent.cxx | 59 +- .../SimCore/G4User/PrimaryGeneratorAction.cxx | 4 +- .../src/SimCore/Generators/GenieGenerator.cxx | 533 ++++++++---------- .../Reweight/GenieReweightProducer.cxx | 307 +++++----- SimCore/src/SimCore/Simulator.cxx | 7 +- 15 files changed, 900 insertions(+), 923 deletions(-) diff --git a/DQM/include/DQM/GenieTruthDQM.h b/DQM/include/DQM/GenieTruthDQM.h index c4e63dd98..127e887f9 100644 --- a/DQM/include/DQM/GenieTruthDQM.h +++ b/DQM/include/DQM/GenieTruthDQM.h @@ -38,7 +38,7 @@ class GenieTruthDQM : public framework::Analyzer { /** * Grab the run number... */ - virtual void onNewRun(const ldmx::RunHeader &runHeader); + virtual void onNewRun(const ldmx::RunHeader& runHeader); /** * Fills histograms/ntuples @@ -54,5 +54,4 @@ class GenieTruthDQM : public framework::Analyzer { }; } // namespace dqm - #endif // DQM_GENIETRUTH_H diff --git a/DQM/src/DQM/GenieTruthDQM.cxx b/DQM/src/DQM/GenieTruthDQM.cxx index 3d8b95cea..7aa3cd9ca 100644 --- a/DQM/src/DQM/GenieTruthDQM.cxx +++ b/DQM/src/DQM/GenieTruthDQM.cxx @@ -4,187 +4,200 @@ #include "DQM/GenieTruthDQM.h" -#include "SimCore/Event/HepMC3GenEvent.h" +#include #include "GENIE/Framework/Conventions/KineVar.h" - -#include +#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"); +void GenieTruthDQM::configure(framework::config::Parameters& ps) { + hepmc3CollName_ = ps.getParameter("hepmc3CollName"); + hepmc3PassName_ = ps.getParameter("hepmc3PassName"); + return; +} - 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"); +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"); +} - 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"); +void GenieTruthDQM::onNewRun(const ldmx::RunHeader& runHeader) { + run_number_ = runHeader.getRunNumber(); +} - 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"); +void GenieTruthDQM::analyze(const framework::Event& event) { + getHistoDirectory(); - 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"); + ntuple_.clear(); + ntuple_.setVar("run", run_number_); + ntuple_.setVar("event", event.getEventNumber()); - } + auto hepmc3_col = event.getObject >( + hepmc3CollName_, hepmc3PassName_); - void GenieTruthDQM::onNewRun(const ldmx::RunHeader &runHeader) { - run_number_ = runHeader.getRunNumber(); + if (hepmc3_col.size() < 1) { + ntuple_.fill(); + return; } - 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< std::vector >(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]); + 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]); } + } - //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(); + // 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]); + } - return; + // 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/SimCore/include/SimCore/Event/EventWeights.h b/SimCore/include/SimCore/Event/EventWeights.h index e4e5687df..35d3f78a6 100644 --- a/SimCore/include/SimCore/Event/EventWeights.h +++ b/SimCore/include/SimCore/Event/EventWeights.h @@ -5,145 +5,147 @@ #ifndef SIM_CORE_EventWeights_H #define SIM_CORE_EventWeights_H -#include "TObject.h" -#include #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< VariationType, std::vector > var_map ) - : variations_map_(var_map) {} - - void Clear(); - void Print(); - - std::vector getWeights() const { return weights_; } - std::map< VariationType, std::vector > 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< VariationType, double > 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< VariationType, std::vector > 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; - } - - }; -} - -#endif //SIM_CORE_EventWeights_H +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 index 57c99b8cf..deac2cfd2 100644 --- a/SimCore/include/SimCore/Event/HepMC3GenEvent.h +++ b/SimCore/include/SimCore/Event/HepMC3GenEvent.h @@ -8,31 +8,28 @@ #ifndef SIM_CORE_HEPMC3GENEVENT_H #define SIM_CORE_HEPMC3GENEVENT_H -#include "TObject.h" -#include "HepMC3/GenEvent.h" #include "HepMC3/Data/GenEventData.h" +#include "HepMC3/GenEvent.h" +#include "TObject.h" namespace ldmx { - class HepMC3GenEvent : public HepMC3::GenEventData { - - public: - - //HepMC3GenEvent() : HepMC3::GenEventData() {} - //HepMC3GenEvent(const HepMC3::GenEvent& event) : HepMC3::GenEvent(event) {} - - void Clear(); - void Print() const; +class HepMC3GenEvent : public HepMC3::GenEventData { + public: + // HepMC3GenEvent() : HepMC3::GenEventData() {} + // HepMC3GenEvent(const HepMC3::GenEvent& event) : HepMC3::GenEvent(event) {} - HepMC3::GenEvent getHepMCGenEvent() const; + void Clear(); + void Print() const; - std::string get_as_string() const; + HepMC3::GenEvent getHepMCGenEvent() const; - public: - ClassDef(HepMC3GenEvent, 1); + std::string get_as_string() const; - }; + public: + ClassDef(HepMC3GenEvent, 1); +}; -} +} // namespace ldmx -#endif //SIM_CORE_HEPMC3GENEVENT_H +#endif // SIM_CORE_HEPMC3GENEVENT_H diff --git a/SimCore/include/SimCore/Generators/GenieGenerator.h b/SimCore/include/SimCore/Generators/GenieGenerator.h index bb72df7ca..5c933856b 100644 --- a/SimCore/include/SimCore/Generators/GenieGenerator.h +++ b/SimCore/include/SimCore/Generators/GenieGenerator.h @@ -13,7 +13,6 @@ #include "TRandom.h" #include "TRandomGen.h" - //------------// // GENIE // //------------// @@ -23,10 +22,10 @@ //------------// // LDMX // //------------// -#include "SimCore/PrimaryGenerator.h" - -#include #include +#include + +#include "SimCore/PrimaryGenerator.h" // Forward declarations class G4Event; @@ -56,7 +55,8 @@ class GenieGenerator : public simcore::PrimaryGenerator { * tune : name of GENIE tune * seed : seed for random generator */ - GenieGenerator(const std::string& name, const framework::config::Parameters& parameters); + GenieGenerator(const std::string& name, + const framework::config::Parameters& parameters); /// Destructor ~GenieGenerator(); @@ -98,13 +98,14 @@ class GenieGenerator : public simcore::PrimaryGenerator { 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 - + 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 diff --git a/SimCore/include/SimCore/PrimaryGenerator.h b/SimCore/include/SimCore/PrimaryGenerator.h index 6b70dca7a..536ab3c9a 100644 --- a/SimCore/include/SimCore/PrimaryGenerator.h +++ b/SimCore/include/SimCore/PrimaryGenerator.h @@ -22,7 +22,6 @@ #include "Framework/Configure/Parameters.h" #include "Framework/RunHeader.h" #include "SimCore/Factory.h" - #include "UserEventInformation.h" // Forward Declarations @@ -72,7 +71,7 @@ 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 index 855f4a9fb..964bfeee5 100644 --- a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h +++ b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h @@ -5,109 +5,104 @@ #ifndef SIMCORE_GENIEREWEIGHTPRODUCER_H #define SIMCORE_GENIEREWEIGHTPRODUCER_H -#include "SimCore/Event/EventWeights.h" +#include +#include #include "Framework/EventProcessor.h" - #include "RwFramework/GSyst.h" - -#include -#include +#include "SimCore/Event/EventWeights.h" namespace genie { - class Interaction; - class HepMC3Converter; -} +class Interaction; +class HepMC3Converter; +} // namespace genie namespace genie::rew { - class GReWeight; +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< ldmx::EventWeights::VariationType, std::vector > 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; - } - return genie::rew::EGSyst::kNullSystematic; - } - - }; -} - -#endif //SIMCORE_GENIEREWEIGHTPRODUCER_H +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; + } + 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 de2489233..0a0b3d23c 100644 --- a/SimCore/include/SimCore/Simulator.h +++ b/SimCore/include/SimCore/Simulator.h @@ -27,11 +27,11 @@ #include "SimCore/SimulatorBase.h" namespace genie { - class Interaction; +class Interaction; } namespace ldmx { - class HepMC3GenEvent; +class HepMC3GenEvent; } class G4UImanager; diff --git a/SimCore/include/SimCore/UserEventInformation.h b/SimCore/include/SimCore/UserEventInformation.h index 3a1a8e799..7fbe28f99 100644 --- a/SimCore/include/SimCore/UserEventInformation.h +++ b/SimCore/include/SimCore/UserEventInformation.h @@ -1,10 +1,10 @@ #ifndef SIMCORE_USEREVENTINFORMATION_H #define SIMCORE_USEREVENTINFORMATION_H -#include "G4VUserEventInformation.hh" +#include +#include "G4VUserEventInformation.hh" #include "SimCore/Event/HepMC3GenEvent.h" -#include namespace simcore { @@ -120,9 +120,13 @@ class UserEventInformation : public G4VUserEventInformation { * @returns true if it was */ bool wasLastStepEN() const { return last_step_en_; } - - void addHepMC3GenEvent(ldmx::HepMC3GenEvent event) { hepmc3_events_.push_back(event); } - std::vector getHepMC3GenEvents() { return hepmc3_events_; } + + 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 @@ -176,8 +180,7 @@ class UserEventInformation : public G4VUserEventInformation { /** * a collection of HepMC3 event records. */ - std::vector< ldmx::HepMC3GenEvent > hepmc3_events_; - + std::vector hepmc3_events_; }; } // namespace simcore diff --git a/SimCore/src/SimCore/Event/EventWeights.cxx b/SimCore/src/SimCore/Event/EventWeights.cxx index c05e6c535..9c524142c 100644 --- a/SimCore/src/SimCore/Event/EventWeights.cxx +++ b/SimCore/src/SimCore/Event/EventWeights.cxx @@ -3,23 +3,25 @@ // #include "SimCore/Event/EventWeights.h" + #include namespace ldmx { - void EventWeights::Clear() { - weights_.clear(); - variations_map_.clear(); - } +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 #include +#include -namespace ldmx{ +#include "HepMC3/Print.h" +#include "HepMC3/WriterAscii.h" - 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(); - } +namespace ldmx { - void HepMC3GenEvent::Print() const { - HepMC3::GenEvent ev; ev.read_data(*this); - HepMC3::Print::line(ev,true); //print attributes - } +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(); +} - HepMC3::GenEvent HepMC3GenEvent::getHepMCGenEvent() const { - HepMC3::GenEvent ev; ev.read_data(*this); - return ev; - } +void HepMC3GenEvent::Print() const { + HepMC3::GenEvent ev; + ev.read_data(*this); + HepMC3::Print::line(ev, true); // print attributes +} - std::string HepMC3GenEvent::get_as_string() const { - HepMC3::GenEvent ev; ev.read_data(*this); +HepMC3::GenEvent HepMC3GenEvent::getHepMCGenEvent() const { + HepMC3::GenEvent ev; + ev.read_data(*this); + return ev; +} - std::stringstream ss; - HepMC3::WriterAscii writer(ss); +std::string HepMC3GenEvent::get_as_string() const { + HepMC3::GenEvent ev; + ev.read_data(*this); - writer.write_event(ev); - return ss.str(); + 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 341782377..eae6d9963 100644 --- a/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx +++ b/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx @@ -72,9 +72,9 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event) { // Make our information container and give it to geant4 // G4Event owns the event information and will delete it auto event_info = new UserEventInformation; - + PrimaryGenerator::Factory::get().apply([event](const auto& generator) { - generator->GeneratePrimaryVertex(event); + generator->GeneratePrimaryVertex(event); }); // smear all primary vertices (if activated) diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx index 4bd80e789..07e2be530 100644 --- a/SimCore/src/SimCore/Generators/GenieGenerator.cxx +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -5,211 +5,190 @@ */ #include "SimCore/Generators/GenieGenerator.h" + #include "SimCore/UserPrimaryParticleInformation.h" // GENIE -#include "Framework/Utils/AppInit.h" -#include "GENIE/Framework/Utils/RunOpt.h" -#include "Framework/Utils/XSecSplineList.h" -#include "GENIE/Framework/Interaction/InitialState.h" -#include "Framework/EventGen/EventRecord.h" -#include "Framework/GHEP/GHepParticle.h" -#include "HepMC3/GenEvent.h" - -#include "Framework/Conventions/XmlParserStatus.h" -#include "Framework/Conventions/GBuild.h" #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/GFluxI.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/NtpWriter.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/XSecSplineList.h" #include "Framework/Utils/StringUtils.h" -#include "Framework/Utils/PrintUtils.h" #include "Framework/Utils/SystemUtils.h" -#include "Framework/Utils/CmdLnArgParser.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 + #include "Math/Interpolator.h" // standard +#include #include #include -#include namespace simcore { namespace generators { -void GenieGenerator::fillConfig(const framework::config::Parameters& p) -{ +void GenieGenerator::fillConfig(const framework::config::Parameters& p) { + energy_ = p.getParameter("energy"); // * GeV; - energy_ = p.getParameter("energy");// * GeV; + targets_ = p.getParameter >("targets"); + abundances_ = p.getParameter >("abundances"); - targets_ = p.getParameter< std::vector >("targets"); - abundances_ = p.getParameter< std::vector >("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 - time_ = p.getParameter("time");// * ns; - position_ = p.getParameter< std::vector >("position"); // mm - beam_size_ = p.getParameter< std::vector >("beam_size"); // mm - direction_ = p.getParameter< std::vector >("direction"); - target_thickness_ = p.getParameter("target_thickness"); //mm - - tune_ = p.getParameter("tune"); + tune_ = p.getParameter("tune"); spline_file_ = p.getParameter("spline_file"); - message_threshold_file_ = p.getParameter("message_threshold_file"); + message_threshold_file_ = + p.getParameter("message_threshold_file"); - verbosity_ = p.getParameter("verbosity"); - + 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; - } +bool GenieGenerator::validateConfig() { + bool ret = true; - 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 (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(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 (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(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]); + 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; + // 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; i_a0) + 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; 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; + << ", 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; i_d(""), - const_cast("--event-generator-list"), - const_cast("EM")}; - genie::RunOpt::Instance()->ReadFromCommandLine(3,inarr); + char* inarr[3] = {const_cast(""), + const_cast("--event-generator-list"), + const_cast("EM")}; + genie::RunOpt::Instance()->ReadFromCommandLine(3, inarr); } - //set message thresholds + // set message thresholds genie::utils::app_init::MesgThresholds(message_threshold_file_); - //set tune info + // set tune info genie::RunOpt::Instance()->SetTuneName(tune_); - if ( ! genie::RunOpt::Instance()->Tune() ) { - EXCEPTION_RAISE("ConfigurationException","No TuneId in RunOption."); + if (!genie::RunOpt::Instance()->Tune()) { + EXCEPTION_RAISE("ConfigurationException", "No TuneId in RunOption."); } genie::RunOpt::Instance()->BuildTune(); @@ -220,168 +199,156 @@ void GenieGenerator::initializeGENIE() std::cout << "Initializing GENIE with seed " << seed << std::endl; genie::utils::app_init::RandGen(seed); */ - - //give it the splint file and require it + + // give it the splint file and require it genie::utils::app_init::XSecTable(spline_file_, true); - //set GHEP print level (needed?) + // set GHEP print level (needed?) genie::GHepRecord::SetPrintLevel(0); - //setup for event driver - evg_driver_.SetEventGeneratorList(genie::RunOpt::Instance()->EventGeneratorList()); + // setup for event driver + evg_driver_.SetEventGeneratorList( + genie::RunOpt::Instance()->EventGeneratorList()); evg_driver_.SetUnphysEventMask(*genie::RunOpt::Instance()->UnphysEventMask()); - } -void GenieGenerator::calculateTotalXS() -{ - //initializing... +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_tgetSeed(); - if(verbosity_>=1) + 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; - - } + if (targets_.size() > 0) { + double rand_uniform = G4Random::getTheGenerator()->flat(); - 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) + 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); + << "(" << 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 + // setup the initial election TParticle initial_e; initial_e.SetPdgCode(11); - double elec_i_p = std::sqrt(energy_*energy_ - (double)initial_e.GetMass()*(double)initial_e.GetMass()); - initial_e.SetMomentum(elec_i_p*direction_[0], - elec_i_p*direction_[1], - elec_i_p*direction_[2], - energy_); + double elec_i_p = + std::sqrt(energy_ * energy_ - + (double)initial_e.GetMass() * (double)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) + 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; + << "(" << 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) + // 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); + 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); @@ -389,60 +356,56 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* event) hepmc3_genie->write_data(hepmc3_ldmx_genie); ev_info->addHepMC3GenEvent(hepmc3_ldmx_genie); event->SetUserInformation(ev_info); - - //setup the primary vertex now - + + // setup the primary vertex now + G4PrimaryVertex* vertex = new G4PrimaryVertex(); - vertex->SetPosition(x_pos, - y_pos, - z_pos); + vertex->SetPosition(x_pos, y_pos, z_pos); vertex->SetWeight(genie_event->Weight()); - - //loop over the entries and add to the G4Event + + // loop over the entries and add to the G4Event int nentries = genie_event->GetEntries(); - if(verbosity_>=1){ + if (verbosity_ >= 1) { std::cout << "---------- " - << "Generated Event " << n_events_generated_+1 - << " ----------" - << std::endl; + << "Generated Event " << n_events_generated_ + 1 << " ----------" + << std::endl; } - - for(int i_p=0; i_pStatus()!=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 + 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; - + // 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_); + rh.setStringParameter(id + " Class", "simcore::generators::GenieGenerator"); + rh.setStringParameter(id + "GenieTune", tune_); } } // namespace generators diff --git a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx index 57b7e81aa..9e412b65e 100644 --- a/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx +++ b/SimCore/src/SimCore/Reweight/GenieReweightProducer.cxx @@ -3,200 +3,201 @@ // #include "SimCore/Reweight/GenieReweightProducer.h" -#include "SimCore/Event/HepMC3GenEvent.h" -#include "SimCore/Event/EventWeights.h" -#include "Framework/EventGen/HepMC3Converter.h" +#include +#include + #include "Framework/EventGen/EventRecord.h" +#include "Framework/EventGen/HepMC3Converter.h" #include "Framework/Interaction/Interaction.h" #include "Framework/Utils/RunOpt.h" - -#include "RwFramework/GReWeightI.h" -#include "RwFramework/GSystSet.h" -#include "RwFramework/GSyst.h" -#include "RwFramework/GReWeight.h" -#include "RwCalculators/GReWeightNuXSecNCEL.h" -#include "RwCalculators/GReWeightNuXSecCCQE.h" -#include "RwCalculators/GReWeightNuXSecCCRES.h" -#include "RwCalculators/GReWeightNuXSecCOH.h" -#include "RwCalculators/GReWeightNonResonanceBkg.h" -#include "RwCalculators/GReWeightFGM.h" +#include "HepMC3/GenEvent.h" +#include "HepMC3/Print.h" +#include "RwCalculators/GReWeightAGKY.h" #include "RwCalculators/GReWeightDISNuclMod.h" -#include "RwCalculators/GReWeightResonanceDecay.h" +#include "RwCalculators/GReWeightDeltaradAngle.h" +#include "RwCalculators/GReWeightFGM.h" #include "RwCalculators/GReWeightFZone.h" #include "RwCalculators/GReWeightINuke.h" -#include "RwCalculators/GReWeightAGKY.h" +#include "RwCalculators/GReWeightINukeParams.h" +#include "RwCalculators/GReWeightNonResonanceBkg.h" +#include "RwCalculators/GReWeightNuXSecCCQE.h" #include "RwCalculators/GReWeightNuXSecCCQEaxial.h" #include "RwCalculators/GReWeightNuXSecCCQEvec.h" -#include "RwCalculators/GReWeightNuXSecNCRES.h" +#include "RwCalculators/GReWeightNuXSecCCRES.h" +#include "RwCalculators/GReWeightNuXSecCOH.h" #include "RwCalculators/GReWeightNuXSecDIS.h" - -#include "RwCalculators/GReWeightINukeParams.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 "RwCalculators/GReWeightDeltaradAngle.h" - -#include "HepMC3/GenEvent.h" -#include "HepMC3/Print.h" - -#include -#include - - +#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< std::vector >("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("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)); + } +} - genie::RunOpt::Instance()->SetTuneName(tune_); - if ( ! genie::RunOpt::Instance()->Tune() ) { - EXCEPTION_RAISE("ConfigurationException","No TuneId in RunOption."); - } - genie::RunOpt::Instance()->BuildTune(); +void GenieReweightProducer::reinitializeGenieReweight() { + delete genie_rw_; + genie_rw_ = new genie::rew::GReWeight; - genie_rw_->AdoptWghtCalc( "hadro_fzone", new genie::rew::GReWeightFZone ); - genie_rw_->AdoptWghtCalc( "hadro_intranuke", new genie::rew::GReWeightINuke ); + genie::RunOpt::Instance()->SetTuneName(tune_); + if (!genie::RunOpt::Instance()->Tune()) { + EXCEPTION_RAISE("ConfigurationException", "No TuneId in RunOption."); + } + genie::RunOpt::Instance()->BuildTune(); - auto & syst = genie_rw_->Systematics(); - for (auto var : variation_map_) - syst.Init(variation_type_to_genie_dial(var.first)); - } + genie_rw_->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); + genie_rw_->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); - 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(); + 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; } - void GenieReweightProducer::produce(framework::Event &event) { - - //grab the input hepmc3 event collection - auto hepmc3_col = event.getObject< std::vector >(hepmc3CollName_, hepmc3PassName_); - - //create an output weights - ldmx::EventWeights ev_weights(variation_map_); + if (tune_ != new_tune) { + if (verbosity_ >= 0) + std::cout << "Found new tune " << new_tune << " (used to be " << tune_ + << ")" << std::endl; + tune_ = new_tune; + reinitializeGenieReweight(); + } +} - for(size_t i_w=0; i_w < n_weights_; ++i_w) { +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(); +} - double running_weight=1; +void GenieReweightProducer::produce(framework::Event& event) { + // grab the input hepmc3 event collection + auto hepmc3_col = event.getObject >( + hepmc3CollName_, hepmc3PassName_); - reconfigureGenieReweight(i_w); + // create an output weights + ldmx::EventWeights ev_weights(variation_map_); - //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) { + for (size_t i_w = 0; i_w < n_weights_; ++i_w) { + double running_weight = 1; - 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); + reconfigureGenieReweight(i_w); - //print it out to check it ... - if(i_w==0 && verbosity_>=1) HepMC3::Print::line(hepmc3_genev, true); //print attributes + // 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); - //now convert to genie event record - auto genie_ev_record_ptr = hepMC3Converter_->RetrieveGHEP(hepmc3_genev); + // print it out to check it ... + if (i_w == 0 && verbosity_ >= 1) + HepMC3::Print::line(hepmc3_genev, true); // print attributes - //print that out too ... - if(i_w==0 && verbosity_>=1) genie_ev_record_ptr->Print(std::cout); + // now convert to genie event record + auto genie_ev_record_ptr = hepMC3Converter_->RetrieveGHEP(hepmc3_genev); - //auto this_weight = 1.0 + var_value*0.05; - auto this_weight = genie_rw_->CalcWeight(*genie_ev_record_ptr); + // print that out too ... + if (i_w == 0 && verbosity_ >= 1) genie_ev_record_ptr->Print(std::cout); - running_weight = running_weight * this_weight; + // auto this_weight = 1.0 + var_value*0.05; + auto this_weight = genie_rw_->CalcWeight(*genie_ev_record_ptr); - }//end loop over interactions in event + running_weight = running_weight * this_weight; - ev_weights.addWeight(running_weight); + } // end loop over interactions in event - } //end loop over weights + ev_weights.addWeight(running_weight); - if(verbosity_>=1) - ev_weights.Print(); + } // end loop over weights - event.add(eventWeightsCollName_, ev_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 302c52032..26466fe1d 100644 --- a/SimCore/src/SimCore/Simulator.cxx +++ b/SimCore/src/SimCore/Simulator.cxx @@ -15,12 +15,12 @@ #include "Framework/RandomNumberSeedService.h" #include "Framework/Version.h" //for LDMX_INSTALL path - /*~~~~~~~~~~~~~*/ /* SimCore */ /*~~~~~~~~~~~~~*/ #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" @@ -28,7 +28,6 @@ #include "SimCore/SensitiveDetector.h" #include "SimCore/UserEventInformation.h" #include "SimCore/XsecBiasingOperator.h" -#include "SimCore/Event/HepMC3GenEvent.h" /*~~~~~~~~~~~~~~*/ /* Geant4 */ @@ -164,12 +163,12 @@ void Simulator::produce(framework::Event& event) { 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){ + for (auto& hepmc3ev : hepmc3_events) { hepmc3ev.event_number = event.getEventHeader().getEventNumber(); } event.add("SimHepMC3Events", hepmc3_events); From 96f698064f435dddc2b6f1290bf81aab5b492c9d Mon Sep 17 00:00:00 2001 From: Sanjit <68086841+SanjitMasanam@users.noreply.github.com> Date: Thu, 17 Oct 2024 11:35:11 -0700 Subject: [PATCH 12/14] Introduce `ldmx-reduced-v2` (#1489) Implemented reduced ldmx v2 geometry + modified dependent scripts --- DetDescr/python/EcalGeometry.py | 19 +- DetDescr/python/HcalGeometry.py | 2 +- Detectors/data/ldmx-reduced-v2/constants.gdml | 471 ++++++ Detectors/data/ldmx-reduced-v2/detector.gdml | 215 +++ Detectors/data/ldmx-reduced-v2/ecal.gdml | 1297 ++++++++++++++ .../ecal_motherboard5_assembly.gdml | 286 ++++ .../ecal_motherboard6_assembly.gdml | 372 ++++ .../ldmx-reduced-v2/ecal_support_box.gdml | 1491 +++++++++++++++++ Detectors/data/ldmx-reduced-v2/hcal.gdml | 448 +++++ Detectors/data/ldmx-reduced-v2/materials.gdml | 45 + Detectors/data/ldmx-reduced-v2/recoil.gdml | 140 ++ .../data/ldmx-reduced-v2/scoring_planes.gdml | 363 ++++ Detectors/data/ldmx-reduced-v2/tagger.gdml | 89 + Detectors/data/ldmx-reduced-v2/target.gdml | 112 ++ .../data/ldmx-reduced-v2/trig_scint.gdml | 109 ++ Ecal/python/digi.py | 12 + 16 files changed, 5467 insertions(+), 4 deletions(-) create mode 100644 Detectors/data/ldmx-reduced-v2/constants.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/detector.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/ecal.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/ecal_motherboard5_assembly.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/ecal_motherboard6_assembly.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/ecal_support_box.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/hcal.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/materials.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/recoil.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/scoring_planes.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/tagger.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/target.gdml create mode 100644 Detectors/data/ldmx-reduced-v2/trig_scint.gdml 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 + ] From d5757c090577817dbbe55dbdccdba2e878eed1f7 Mon Sep 17 00:00:00 2001 From: Wesley Ketchum Date: Tue, 22 Oct 2024 15:38:19 +0200 Subject: [PATCH 13/14] cleanup some compile warnings and fix failing Hcal and Run_SimCore_basic tests --- SimCore/include/SimCore/Event/HepMC3GenEvent.h | 2 +- SimCore/include/SimCore/Reweight/GenieReweightProducer.h | 2 ++ SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx | 1 + SimCore/src/SimCore/Generators/GenieGenerator.cxx | 6 +++--- 4 files changed, 7 insertions(+), 4 deletions(-) diff --git a/SimCore/include/SimCore/Event/HepMC3GenEvent.h b/SimCore/include/SimCore/Event/HepMC3GenEvent.h index deac2cfd2..ef65a22d6 100644 --- a/SimCore/include/SimCore/Event/HepMC3GenEvent.h +++ b/SimCore/include/SimCore/Event/HepMC3GenEvent.h @@ -16,7 +16,7 @@ namespace ldmx { class HepMC3GenEvent : public HepMC3::GenEventData { public: - // HepMC3GenEvent() : HepMC3::GenEventData() {} + // ~HepMC3GenEvent() {} // HepMC3GenEvent(const HepMC3::GenEvent& event) : HepMC3::GenEvent(event) {} void Clear(); diff --git a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h index 964bfeee5..57b863003 100644 --- a/SimCore/include/SimCore/Reweight/GenieReweightProducer.h +++ b/SimCore/include/SimCore/Reweight/GenieReweightProducer.h @@ -99,6 +99,8 @@ class GenieReweightProducer : public framework::Producer { 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; } diff --git a/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx b/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx index eae6d9963..a60f168bd 100644 --- a/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx +++ b/SimCore/src/SimCore/G4User/PrimaryGeneratorAction.cxx @@ -72,6 +72,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event) { // Make our information container and give it to geant4 // G4Event owns the event information and will delete it auto event_info = new UserEventInformation; + event->SetUserInformation(event_info); PrimaryGenerator::Factory::get().apply([event](const auto& generator) { generator->GeneratePrimaryVertex(event); diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx index 07e2be530..daa08207a 100644 --- a/SimCore/src/SimCore/Generators/GenieGenerator.cxx +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -147,7 +147,7 @@ bool GenieGenerator::validateConfig() { std::cout << "Will renormalize abundances to unity!" << std::endl; } - for (size_t i_a; i_a < abundances_.size(); ++i_a) { + for (size_t i_a=0; i_a < abundances_.size(); ++i_a) { abundances_[i_a] = abundances_[i_a] / abundance_sum; if (verbosity_ > 0) @@ -164,7 +164,7 @@ bool GenieGenerator::validateConfig() { << direction_[2] << ")" << std::endl; ret = false; } - for (size_t i_d; i_d < direction_.size(); ++i_d) + 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.); @@ -330,7 +330,7 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* event) { initial_e.SetPdgCode(11); double elec_i_p = std::sqrt(energy_ * energy_ - - (double)initial_e.GetMass() * (double)initial_e.GetMass()); + 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; From d0b995e41c9ed51f0d29c51eed01f8eed1562101 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 22 Oct 2024 13:41:11 +0000 Subject: [PATCH 14/14] Apply clang-format --- SimCore/src/SimCore/Generators/GenieGenerator.cxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/SimCore/src/SimCore/Generators/GenieGenerator.cxx b/SimCore/src/SimCore/Generators/GenieGenerator.cxx index daa08207a..dc550616f 100644 --- a/SimCore/src/SimCore/Generators/GenieGenerator.cxx +++ b/SimCore/src/SimCore/Generators/GenieGenerator.cxx @@ -147,7 +147,7 @@ bool GenieGenerator::validateConfig() { std::cout << "Will renormalize abundances to unity!" << std::endl; } - for (size_t i_a=0; i_a < abundances_.size(); ++i_a) { + for (size_t i_a = 0; i_a < abundances_.size(); ++i_a) { abundances_[i_a] = abundances_[i_a] / abundance_sum; if (verbosity_ > 0) @@ -164,7 +164,7 @@ bool GenieGenerator::validateConfig() { << direction_[2] << ")" << std::endl; ret = false; } - for (size_t i_d=0; i_d < direction_.size(); ++i_d) + 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.); @@ -329,8 +329,7 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* event) { TParticle initial_e; initial_e.SetPdgCode(11); double elec_i_p = - std::sqrt(energy_ * energy_ - - initial_e.GetMass() * initial_e.GetMass()); + 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;