From a53ed41d3381b1eb781b655268f1508d192d054f Mon Sep 17 00:00:00 2001 From: mahmoudali2 Date: Tue, 20 Feb 2024 18:04:11 +0100 Subject: [PATCH 1/6] Adding the first draft for the Muon system digitizer --- CMakeLists.txt | 1 + MUONdigi/CMakeLists.txt | 43 +++++++ MUONdigi/include/MUONsimpleDigitizer.h | 75 ++++++++++++ MUONdigi/src/MUONsimpleDigitizer.cpp | 115 ++++++++++++++++++ MUONdigi/test/runMUONsimpleDigitizer.py | 150 ++++++++++++++++++++++++ 5 files changed, 384 insertions(+) create mode 100644 MUONdigi/CMakeLists.txt create mode 100644 MUONdigi/include/MUONsimpleDigitizer.h create mode 100644 MUONdigi/src/MUONsimpleDigitizer.cpp create mode 100644 MUONdigi/test/runMUONsimpleDigitizer.py diff --git a/CMakeLists.txt b/CMakeLists.txt index ffd4e85..535454c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,6 +51,7 @@ include(CTest) add_subdirectory(DCHdigi) add_subdirectory(ARCdigi) add_subdirectory(VTXdigi) +add_subdirectory(MUONdigi) set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH}) add_subdirectory(Tracking) diff --git a/MUONdigi/CMakeLists.txt b/MUONdigi/CMakeLists.txt new file mode 100644 index 0000000..81a4095 --- /dev/null +++ b/MUONdigi/CMakeLists.txt @@ -0,0 +1,43 @@ +project(MUONdigi) + +file(GLOB sources + ${PROJECT_SOURCE_DIR}/src/*.cpp +) + +file(GLOB headers + ${PROJECT_SOURCE_DIR}/include/*.h +) + +gaudi_add_module(MUONdigi + SOURCES ${sources} + LINK + Gaudi::GaudiKernel + EDM4HEP::edm4hep + k4FWCore::k4FWCore + DD4hep::DDRec +) + +target_include_directories(MUONdigi PUBLIC + $ + $ +) + +set_target_properties(MUONdigi PROPERTIES PUBLIC_HEADER "${headers}") + +file(GLOB scripts + ${PROJECT_SOURCE_DIR}/test/*.py +) + +file(COPY ${scripts} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/test) + +install(TARGETS MUONdigi + EXPORT ${CMAKE_PROJECT_NAME}Targets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/@{CMAKE_PROJECT_NAME}" COMPONENT dev +) + +install(FILES ${scripts} DESTINATION test) + +SET(test_name "test_MUONsimpleDigitizer") +ADD_TEST(NAME t_${test_name} COMMAND k4run test/runMUONsimpleDigitizer.py) diff --git a/MUONdigi/include/MUONsimpleDigitizer.h b/MUONdigi/include/MUONsimpleDigitizer.h new file mode 100644 index 0000000..14ea05f --- /dev/null +++ b/MUONdigi/include/MUONsimpleDigitizer.h @@ -0,0 +1,75 @@ +#pragma once + +// GAUDI +#include "Gaudi/Property.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/IRndmGenSvc.h" +#include "GaudiKernel/RndmGenerators.h" + +// K4FWCORE +#include "k4FWCore/DataHandle.h" +#include "k4Interface/IGeoSvc.h" + +// EDM4HEP +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/TrackerHitCollection.h" + +// DD4HEP +#include "DD4hep/Detector.h" // for dd4hep::VolumeManager +#include "DDSegmentation/BitFieldCoder.h" + +/** @class MUONsimpleDigitizer + * + * Algorithm for creating digitized Muon system hits (still based on edm4hep::TrackerHit) from edm4hep::SimTrackerHit. + * You have to specify the expected resolution in z and in xy (distance to the wire). The smearing is applied in the wire reference frame. + * + * @author Mahmoud Ali + * @date 2023-09 + * + */ + +class MUONsimpleDigitizer : public GaudiAlgorithm { +public: + explicit MUONsimpleDigitizer(const std::string&, ISvcLocator*); + virtual ~MUONsimpleDigitizer(); + /** Initialize. + * @return status code + */ + virtual StatusCode initialize() final; + /** Execute. + * @return status code + */ + virtual StatusCode execute() final; + /** Finalize. + * @return status code + */ + virtual StatusCode finalize() final; + +private: + // Input sim tracker hit collection name + DataHandle m_input_sim_hits{"inputSimHits", Gaudi::DataHandle::Reader, this}; + // Output digitized tracker hit collection name + DataHandle m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this}; + + // Detector readout name + Gaudi::Property m_readoutName{this, "readoutName", "MuonChamberBarrelReadout", "Name of the detector readout"}; + // Pointer to the geometry service + ServiceHandle m_geoSvc; + // Decoder for the cellID + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + // Volume manager to get the physical cell sensitive volume + dd4hep::VolumeManager m_volman; +/* + // z position resolution in mm + FloatProperty m_z_resolution{this, "zResolution", 1.0, + "Spatial resolution in the z direction (from reading out the wires at both sides) [mm]"}; + // xy resolution in mm + FloatProperty m_xy_resolution{this, "xyResolution", 0.1, "Spatial resolution in the xy direction [mm]"}; +*/ + // Random Number Service + IRndmGenSvc* m_randSvc; + // Gaussian random number generator used for the smearing of the z position + Rndm::Numbers m_gauss_z; + // Gaussian random number generator used for the smearing of the xy position + Rndm::Numbers m_gauss_xy; +}; diff --git a/MUONdigi/src/MUONsimpleDigitizer.cpp b/MUONdigi/src/MUONsimpleDigitizer.cpp new file mode 100644 index 0000000..c9f3142 --- /dev/null +++ b/MUONdigi/src/MUONsimpleDigitizer.cpp @@ -0,0 +1,115 @@ +#include "MUONsimpleDigitizer.h" + +// DD4hep +#include "DD4hep/Detector.h" +#include "DDRec/Vector3D.h" + +// ROOT +#include "Math/Cylindrical3D.h" + +DECLARE_COMPONENT(MUONsimpleDigitizer) + +MUONsimpleDigitizer::MUONsimpleDigitizer(const std::string& aName, ISvcLocator* aSvcLoc) + : GaudiAlgorithm(aName, aSvcLoc), m_geoSvc("GeoSvc", "MUONsimpleDigitizer") { + declareProperty("inputSimHits", m_input_sim_hits, "Input sim tracker hit collection name"); + declareProperty("outputDigiHits", m_output_digi_hits, "Output digitized tracker hit collection name"); +} + +MUONsimpleDigitizer::~MUONsimpleDigitizer() {} + +StatusCode MUONsimpleDigitizer::initialize() { + // Initialize random services + if (service("RndmGenSvc", m_randSvc).isFailure()) { + error() << "Couldn't get RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } +/* if (m_gauss_z.initialize(m_randSvc, Rndm::Gauss(0., m_z_resolution)).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } + if (m_gauss_xy.initialize(m_randSvc, Rndm::Gauss(0., m_xy_resolution)).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } +*/ + // check if readout exists + if (m_geoSvc->lcdd()->readouts().find(m_readoutName) == m_geoSvc->lcdd()->readouts().end()) { + error() << "Readout <<" << m_readoutName << ">> does not exist." << endmsg; + return StatusCode::FAILURE; + } + // set the cellID decoder + m_decoder = m_geoSvc->lcdd()->readout(m_readoutName).idSpec().decoder(); + // retrieve the volume manager + m_volman = m_geoSvc->lcdd()->volumeManager(); + + return StatusCode::SUCCESS; +} + +StatusCode MUONsimpleDigitizer::execute() { + // Get the input collection with Geant4 hits + const edm4hep::SimTrackerHitCollection* input_sim_hits = m_input_sim_hits.get(); + debug() << "Input Sim Hit collection size: " << input_sim_hits->size() << endmsg; + + // Digitize the sim hits + edm4hep::TrackerHitCollection* output_digi_hits = m_output_digi_hits.createAndPut(); + for (const auto& input_sim_hit : *input_sim_hits) { + auto output_digi_hit = output_digi_hits->create(); + // smear the hit position: need to go in the wire local frame to smear in the direction aligned/perpendicular with the wire for z/distanceToWire, taking e.g. stereo angle into account + // retrieve the cell detElement + dd4hep::DDSegmentation::CellID cellID = input_sim_hit.getCellID(); + auto cellDetElement = m_volman.lookupDetElement(cellID); + // retrieve the wire (in DD4hep 1.23 there is no easy way to access the volume daughters we have to pass by detElements, in later versions volumes can be used) + const std::string& wireDetElementName = + Form("superLayer_%d_layer_%d_phi_%d_wire", m_decoder->get(cellID, "superLayer"), + m_decoder->get(cellID, "layer"), m_decoder->get(cellID, "phi")); + dd4hep::DetElement wireDetElement = cellDetElement.child(wireDetElementName); + // get the transformation matrix used to place the wire + const auto& wireTransformMatrix = wireDetElement.nominal().worldTransformation(); + // Retrieve global position in mm and apply unit transformation (translation matrix is tored in cm) + double simHitGlobalPosition[3] = {input_sim_hit.getPosition().x * dd4hep::mm, + input_sim_hit.getPosition().y * dd4hep::mm, + input_sim_hit.getPosition().z * dd4hep::mm}; + double simHitLocalPosition[3] = {0, 0, 0}; + // get the simHit coordinate in cm in the wire reference frame to be able to apply smearing of radius perpendicular to the wire + wireTransformMatrix.MasterToLocal(simHitGlobalPosition, simHitLocalPosition); + debug() << "Cell ID string: " << m_decoder->valueString(cellID) << endmsg; + ; + debug() << "Global simHit x " << simHitGlobalPosition[0] << " --> Local simHit x " << simHitLocalPosition[0] + << " in cm" << endmsg; + debug() << "Global simHit y " << simHitGlobalPosition[1] << " --> Local simHit y " << simHitLocalPosition[1] + << " in cm" << endmsg; + debug() << "Global simHit z " << simHitGlobalPosition[2] << " --> Local simHit z " << simHitLocalPosition[2] + << " in cm" << endmsg; + // build a vector to easily apply smearing of distance to the wire + dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0], simHitLocalPosition[1], + simHitLocalPosition[2]); + // get the smeared distance to the wire (cylindrical coordinate as the smearing should be perpendicular to the wire) + double smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm; + // smear the z position (in local coordinate the z axis is aligned with the wire i.e. it take the stereo angle into account); + double smearedZ = simHitLocalPositionVector.z() + m_gauss_z.shoot() * dd4hep::mm; + // build the local position vector of the smeared hit using cylindrical coordinates. When we will have edm4hep::MUONHit there will be probably no need + ROOT::Math::Cylindrical3D digiHitLocalPositionVector(smearedDistanceToWire, smearedZ, + simHitLocalPositionVector.phi()); + debug() << "Local simHit distance to the wire: " << simHitLocalPositionVector.rho() + << " Local digiHit distance to the wire: " << smearedDistanceToWire << " in cm" << endmsg; + debug() << "Local simHit z: " << simHitLocalPositionVector.z() + << " Local digiHit distance to the wire: " << smearedZ << " in cm" << endmsg; + debug() << "Local simHit phi: " << simHitLocalPositionVector.phi() + << " Local digiHit distance to the wire: " << digiHitLocalPositionVector.Phi() << endmsg; + // go back to the global frame + double digiHitLocalPosition[3] = {digiHitLocalPositionVector.x(), digiHitLocalPositionVector.y(), + digiHitLocalPositionVector.z()}; + double digiHitGlobalPosition[3] = {0, 0, 0}; + wireTransformMatrix.LocalToMaster(digiHitLocalPosition, digiHitGlobalPosition); + // go back to mm + edm4hep::Vector3d digiHitGlobalPositionVector(digiHitGlobalPosition[0] / dd4hep::mm, + digiHitGlobalPosition[1] / dd4hep::mm, + digiHitGlobalPosition[2] / dd4hep::mm); + output_digi_hit.setPosition(digiHitGlobalPositionVector); + output_digi_hit.setCellID(cellID); + } + debug() << "Output Digi Hit collection size: " << output_digi_hits->size() << endmsg; + return StatusCode::SUCCESS; +} + +StatusCode MUONsimpleDigitizer::finalize() { return StatusCode::SUCCESS; } diff --git a/MUONdigi/test/runMUONsimpleDigitizer.py b/MUONdigi/test/runMUONsimpleDigitizer.py new file mode 100644 index 0000000..4aa179d --- /dev/null +++ b/MUONdigi/test/runMUONsimpleDigitizer.py @@ -0,0 +1,150 @@ +import os + +from Gaudi.Configuration import * + +from Configurables import FCCDataSvc +podioevent = FCCDataSvc("EventDataSvc") + +from GaudiKernel.SystemOfUnits import MeV, GeV, tesla +################## Particle gun setup +momentum = 5 # in GeV +#thetaMin = 90.25 # degrees +#thetaMax = 90.25 # degrees +thetaMin = 20 # degrees +thetaMax = 130 # degrees +pdgCode = 13 # 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ +magneticField = True +_pi = 3.14159 + +from Configurables import GenAlg +genAlg = GenAlg() +from Configurables import MomentumRangeParticleGun +pgun = MomentumRangeParticleGun("ParticleGun_Electron") +pgun.PdgCodes = [pdgCode] +pgun.MomentumMin = momentum * GeV +pgun.MomentumMax = momentum * GeV +pgun.PhiMin = 0 +pgun.PhiMax = 2 * _pi +pgun.ThetaMin = thetaMin * _pi / 180. +pgun.ThetaMax = thetaMax * _pi / 180. +genAlg.SignalProvider = pgun + +genAlg.hepmc.Path = "hepmc" + +from Configurables import HepMCToEDMConverter +hepmc_converter = HepMCToEDMConverter() +hepmc_converter.hepmc.Path="hepmc" +genParticlesOutputName = "genParticles" +hepmc_converter.GenParticles.Path = genParticlesOutputName +hepmc_converter.hepmcStatusList = [] + +################## Simulation setup +# Detector geometry +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +# if FCC_DETECTORS is empty, this should use relative path to working directory +path_to_detector = os.environ.get("FCCDETECTORS", "") +print(path_to_detector) +detectors_to_use=[ + 'Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v01.xml' + ] +# prefix all xmls with path_to_detector +geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use] +geoservice.OutputLevel = INFO + +# Geant4 service +# Configures the Geant simulation: geometry, physics list and user actions + +from Configurables import SimG4FullSimActions, SimG4Alg, SimG4PrimariesFromEdmTool, SimG4SaveParticleHistory +actions = SimG4FullSimActions() +# Uncomment if history from Geant4 decays is needed (e.g. to get the photons from pi0) and set actions=actions in SimG4Svc + uncomment saveHistTool in SimG4Alg +#actions.enableHistory=True +#actions.energyCut = 0.2 * GeV +#saveHistTool = SimG4SaveParticleHistory("saveHistory") + + +# Magnetic field +from Configurables import SimG4ConstantMagneticFieldTool +field = SimG4ConstantMagneticFieldTool("SimG4ConstantMagneticFieldTool", FieldComponentZ = -2 * tesla, FieldOn = magneticField, IntegratorStepper = "ClassicalRK4") + +from Configurables import SimG4Svc +geantservice = SimG4Svc("SimG4Svc", detector = 'SimG4DD4hepDetector', physicslist = "SimG4FtfpBert", actions = actions, magneticField = field) + +# Fixed seed to have reproducible results, change it for each job if you split one production into several jobs +# Mind that if you leave Gaudi handle random seed and some job start within the same second (very likely) you will have duplicates +geantservice.randomNumbersFromGaudi = False +geantservice.seedValue = 4242 + +# Range cut +geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] + +# Geant4 algorithm +# Translates EDM to G4Event, passes the event to G4, writes out outputs via tools +# and a tool that saves the calorimeter hits + +# next, create the G4 algorithm, giving the list of names of tools ("XX/YY") +from Configurables import SimG4PrimariesFromEdmTool +particle_converter = SimG4PrimariesFromEdmTool("EdmConverter") +particle_converter.GenParticles.Path = genParticlesOutputName + +from Configurables import SimG4SaveTrackerHits +savetrackertool = SimG4SaveTrackerHits("SimG4SaveTrackerHits", readoutNames=["MUONHits"]) +savetrackertool.SimTrackHits.Path = "DC_simTrackerHits" + + +from Configurables import SimG4Alg +geantsim = SimG4Alg("SimG4Alg", + outputs= [savetrackertool + #saveHistTool + ], + eventProvider=particle_converter, + OutputLevel=DEBUG) +# Digitize tracker hits +from Configurables import MUONsimpleDigitizer +muon_digitizer = MUONsimpleDigitizer("MUONsimpleDigitizer", + inputSimHits = savetrackertool.SimTrackHits.Path, + outputDigiHits = savetrackertool.SimTrackHits.Path.replace("sim", "digi"), + readoutName = "MUONHits", + xyResolution = 0.1, # mm + zResolution = 1, # mm + OutputLevel=DEBUG +) + +################ Output +from Configurables import PodioOutput +out = PodioOutput("out", + OutputLevel=INFO) +out.outputCommands = ["keep *"] + +import uuid +out.filename = "output_simplifiedDriftChamber_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+".root" + +#CPU information +from Configurables import AuditorSvc, ChronoAuditor +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +genAlg.AuditExecute = True +hepmc_converter.AuditExecute = True +geantsim.AuditExecute = True +out.AuditExecute = True + +from Configurables import EventCounter +event_counter = EventCounter('event_counter') +event_counter.Frequency = 1 + +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg = [ + event_counter, + genAlg, + hepmc_converter, + geantsim, + muon_digitizer, + out + ], + EvtSel = 'NONE', + EvtMax = 100, + ExtSvc = [geoservice, podioevent, geantservice, audsvc], + StopOnSignal = True, + ) From 01eab5375dcec683eaa61469b29451848aeaf817 Mon Sep 17 00:00:00 2001 From: mahmoudali2 Date: Mon, 18 Mar 2024 11:47:08 +0100 Subject: [PATCH 2/6] Fix muon digitizer --- MUONdigi/include/MUONsimpleDigitizer.h | 20 ++++++--- MUONdigi/src/MUONsimpleDigitizer.cpp | 55 +++++++++++++------------ MUONdigi/test/runMUONsimpleDigitizer.py | 31 +++++++++----- 3 files changed, 63 insertions(+), 43 deletions(-) diff --git a/MUONdigi/include/MUONsimpleDigitizer.h b/MUONdigi/include/MUONsimpleDigitizer.h index 14ea05f..c318ddc 100644 --- a/MUONdigi/include/MUONsimpleDigitizer.h +++ b/MUONdigi/include/MUONsimpleDigitizer.h @@ -12,7 +12,15 @@ // EDM4HEP #include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/TrackCollection.h" +#if __has_include("edm4hep/TrackerHit3DCollection.h") +#include "edm4hep/TrackerHit3DCollection.h" +#else #include "edm4hep/TrackerHitCollection.h" +namespace edm4hep { + using TrackerHit3DCollection = edm4hep::TrackerHitCollection; +} // namespace edm4hep +#endif // DD4HEP #include "DD4hep/Detector.h" // for dd4hep::VolumeManager @@ -20,8 +28,8 @@ /** @class MUONsimpleDigitizer * - * Algorithm for creating digitized Muon system hits (still based on edm4hep::TrackerHit) from edm4hep::SimTrackerHit. - * You have to specify the expected resolution in z and in xy (distance to the wire). The smearing is applied in the wire reference frame. + * Algorithm for creating digitized Muon system hits (still based on edm4hep::TrackerHit3D) from edm4hep::SimTrackerHit. + * You have to specify the expected resolution in z and in xy. * * @author Mahmoud Ali * @date 2023-09 @@ -49,7 +57,7 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { // Input sim tracker hit collection name DataHandle m_input_sim_hits{"inputSimHits", Gaudi::DataHandle::Reader, this}; // Output digitized tracker hit collection name - DataHandle m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this}; + DataHandle m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this}; // Detector readout name Gaudi::Property m_readoutName{this, "readoutName", "MuonChamberBarrelReadout", "Name of the detector readout"}; @@ -59,13 +67,13 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { dd4hep::DDSegmentation::BitFieldCoder* m_decoder; // Volume manager to get the physical cell sensitive volume dd4hep::VolumeManager m_volman; -/* + // z position resolution in mm FloatProperty m_z_resolution{this, "zResolution", 1.0, - "Spatial resolution in the z direction (from reading out the wires at both sides) [mm]"}; + "Spatial resolution in the z direction [mm]"}; // xy resolution in mm FloatProperty m_xy_resolution{this, "xyResolution", 0.1, "Spatial resolution in the xy direction [mm]"}; -*/ + // Random Number Service IRndmGenSvc* m_randSvc; // Gaussian random number generator used for the smearing of the z position diff --git a/MUONdigi/src/MUONsimpleDigitizer.cpp b/MUONdigi/src/MUONsimpleDigitizer.cpp index c9f3142..f0a128c 100644 --- a/MUONdigi/src/MUONsimpleDigitizer.cpp +++ b/MUONdigi/src/MUONsimpleDigitizer.cpp @@ -23,7 +23,7 @@ StatusCode MUONsimpleDigitizer::initialize() { error() << "Couldn't get RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } -/* if (m_gauss_z.initialize(m_randSvc, Rndm::Gauss(0., m_z_resolution)).isFailure()) { + if (m_gauss_z.initialize(m_randSvc, Rndm::Gauss(0., m_z_resolution)).isFailure()) { error() << "Couldn't initialize RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } @@ -31,16 +31,16 @@ StatusCode MUONsimpleDigitizer::initialize() { error() << "Couldn't initialize RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } -*/ + // check if readout exists - if (m_geoSvc->lcdd()->readouts().find(m_readoutName) == m_geoSvc->lcdd()->readouts().end()) { + if (m_geoSvc->getDetector()->readouts().find(m_readoutName) == m_geoSvc->getDetector()->readouts().end()) { error() << "Readout <<" << m_readoutName << ">> does not exist." << endmsg; return StatusCode::FAILURE; } // set the cellID decoder - m_decoder = m_geoSvc->lcdd()->readout(m_readoutName).idSpec().decoder(); + m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); // retrieve the volume manager - m_volman = m_geoSvc->lcdd()->volumeManager(); + m_volman = m_geoSvc->getDetector()->volumeManager(); return StatusCode::SUCCESS; } @@ -51,27 +51,30 @@ StatusCode MUONsimpleDigitizer::execute() { debug() << "Input Sim Hit collection size: " << input_sim_hits->size() << endmsg; // Digitize the sim hits - edm4hep::TrackerHitCollection* output_digi_hits = m_output_digi_hits.createAndPut(); + edm4hep::TrackerHit3DCollection* output_digi_hits = m_output_digi_hits.createAndPut(); for (const auto& input_sim_hit : *input_sim_hits) { auto output_digi_hit = output_digi_hits->create(); - // smear the hit position: need to go in the wire local frame to smear in the direction aligned/perpendicular with the wire for z/distanceToWire, taking e.g. stereo angle into account + // smear the hit position: // retrieve the cell detElement dd4hep::DDSegmentation::CellID cellID = input_sim_hit.getCellID(); + debug() << "Digitisation of " << m_readoutName << ", cellID: " << cellID << endmsg; auto cellDetElement = m_volman.lookupDetElement(cellID); - // retrieve the wire (in DD4hep 1.23 there is no easy way to access the volume daughters we have to pass by detElements, in later versions volumes can be used) - const std::string& wireDetElementName = - Form("superLayer_%d_layer_%d_phi_%d_wire", m_decoder->get(cellID, "superLayer"), - m_decoder->get(cellID, "layer"), m_decoder->get(cellID, "phi")); - dd4hep::DetElement wireDetElement = cellDetElement.child(wireDetElementName); - // get the transformation matrix used to place the wire - const auto& wireTransformMatrix = wireDetElement.nominal().worldTransformation(); + + /* const std::string& stripsDetElementName = + Form("system_%d_layer_%d_phi_%d", m_decoder->get(cellID, "layer"), m_decoder->get(cellID, "phi")); + dd4hep::DetElement stripsDetElement = cellDetElement.child(stripsDetElementName); + const auto& stripsTransformMatrix = stripsDetElement.nominal().worldTransformation(); + */ + + const auto& stripsTransformMatrix = m_volman.lookupVolumePlacement(cellID).matrix(); + // Retrieve global position in mm and apply unit transformation (translation matrix is tored in cm) double simHitGlobalPosition[3] = {input_sim_hit.getPosition().x * dd4hep::mm, input_sim_hit.getPosition().y * dd4hep::mm, input_sim_hit.getPosition().z * dd4hep::mm}; double simHitLocalPosition[3] = {0, 0, 0}; - // get the simHit coordinate in cm in the wire reference frame to be able to apply smearing of radius perpendicular to the wire - wireTransformMatrix.MasterToLocal(simHitGlobalPosition, simHitLocalPosition); + // get the simHit coordinate in cm in the strips reference frame + stripsTransformMatrix.MasterToLocal(simHitGlobalPosition, simHitLocalPosition); debug() << "Cell ID string: " << m_decoder->valueString(cellID) << endmsg; ; debug() << "Global simHit x " << simHitGlobalPosition[0] << " --> Local simHit x " << simHitLocalPosition[0] @@ -80,27 +83,27 @@ StatusCode MUONsimpleDigitizer::execute() { << " in cm" << endmsg; debug() << "Global simHit z " << simHitGlobalPosition[2] << " --> Local simHit z " << simHitLocalPosition[2] << " in cm" << endmsg; - // build a vector to easily apply smearing of distance to the wire + // build a vector to easily apply smearing of distance to the strips dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0], simHitLocalPosition[1], simHitLocalPosition[2]); - // get the smeared distance to the wire (cylindrical coordinate as the smearing should be perpendicular to the wire) - double smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm; - // smear the z position (in local coordinate the z axis is aligned with the wire i.e. it take the stereo angle into account); + // smear xy position + double smearedDistanceTostrips = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm; + // smear the z position double smearedZ = simHitLocalPositionVector.z() + m_gauss_z.shoot() * dd4hep::mm; // build the local position vector of the smeared hit using cylindrical coordinates. When we will have edm4hep::MUONHit there will be probably no need - ROOT::Math::Cylindrical3D digiHitLocalPositionVector(smearedDistanceToWire, smearedZ, + ROOT::Math::Cylindrical3D digiHitLocalPositionVector(smearedDistanceTostrips, smearedZ, simHitLocalPositionVector.phi()); - debug() << "Local simHit distance to the wire: " << simHitLocalPositionVector.rho() - << " Local digiHit distance to the wire: " << smearedDistanceToWire << " in cm" << endmsg; + debug() << "Local simHit distance to the strips: " << simHitLocalPositionVector.rho() + << " Local digiHit distance to the strips: " << smearedDistanceTostrips << " in cm" << endmsg; debug() << "Local simHit z: " << simHitLocalPositionVector.z() - << " Local digiHit distance to the wire: " << smearedZ << " in cm" << endmsg; + << " Local digiHit distance to the strips: " << smearedZ << " in cm" << endmsg; debug() << "Local simHit phi: " << simHitLocalPositionVector.phi() - << " Local digiHit distance to the wire: " << digiHitLocalPositionVector.Phi() << endmsg; + << " Local digiHit distance to the strips: " << digiHitLocalPositionVector.Phi() << endmsg; // go back to the global frame double digiHitLocalPosition[3] = {digiHitLocalPositionVector.x(), digiHitLocalPositionVector.y(), digiHitLocalPositionVector.z()}; double digiHitGlobalPosition[3] = {0, 0, 0}; - wireTransformMatrix.LocalToMaster(digiHitLocalPosition, digiHitGlobalPosition); + stripsTransformMatrix.LocalToMaster(digiHitLocalPosition, digiHitGlobalPosition); // go back to mm edm4hep::Vector3d digiHitGlobalPositionVector(digiHitGlobalPosition[0] / dd4hep::mm, digiHitGlobalPosition[1] / dd4hep::mm, diff --git a/MUONdigi/test/runMUONsimpleDigitizer.py b/MUONdigi/test/runMUONsimpleDigitizer.py index 4aa179d..a952a39 100644 --- a/MUONdigi/test/runMUONsimpleDigitizer.py +++ b/MUONdigi/test/runMUONsimpleDigitizer.py @@ -10,8 +10,8 @@ momentum = 5 # in GeV #thetaMin = 90.25 # degrees #thetaMax = 90.25 # degrees -thetaMin = 20 # degrees -thetaMax = 130 # degrees +thetaMin = 0 # degrees +thetaMax = 180 # degrees pdgCode = 13 # 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ magneticField = True _pi = 3.14159 @@ -46,7 +46,8 @@ path_to_detector = os.environ.get("FCCDETECTORS", "") print(path_to_detector) detectors_to_use=[ - 'Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v01.xml' + # 'Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v02.xml' + '/afs/cern.ch/work/m/maali/public/FCCDetectors/Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v02.xml' ] # prefix all xmls with path_to_detector geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use] @@ -88,23 +89,31 @@ particle_converter.GenParticles.Path = genParticlesOutputName from Configurables import SimG4SaveTrackerHits -savetrackertool = SimG4SaveTrackerHits("SimG4SaveTrackerHits", readoutNames=["MUONHits"]) -savetrackertool.SimTrackHits.Path = "DC_simTrackerHits" +#savetrackertool = SimG4SaveTrackerHits("SimG4SaveTrackerHits", readoutNames=["MuonChamberBarrelReadout"]) +#savetrackertool.SimTrackHits.Path = "MUON_simTrackerHits" +saveMuonBarrelTool = SimG4SaveTrackerHits("SimG4SaveMuonBarrelHits", readoutName="MuonChamberBarrelReadout") +saveMuonBarrelTool.SimTrackHits.Path = "muonBarrelSimHits" + +saveMuonPositiveEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonPositiveEndcapHits", readoutName="MuonChamberPositiveEndcapReadout") +saveMuonPositiveEndcapTool.SimTrackHits.Path = "muonPositiveEndcapSimHits" + +saveMuonNegativeEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonNegativeEndcapHits", readoutName="MuonChamberNegativeEndcapReadout") +saveMuonNegativeEndcapTool.SimTrackHits.Path = "muonNegativeEndcapSimHits" from Configurables import SimG4Alg geantsim = SimG4Alg("SimG4Alg", - outputs= [savetrackertool + outputs= [saveMuonBarrelTool, saveMuonPositiveEndcapTool, saveMuonNegativeEndcapTool #saveHistTool ], eventProvider=particle_converter, - OutputLevel=DEBUG) + OutputLevel=INFO) # Digitize tracker hits from Configurables import MUONsimpleDigitizer muon_digitizer = MUONsimpleDigitizer("MUONsimpleDigitizer", - inputSimHits = savetrackertool.SimTrackHits.Path, - outputDigiHits = savetrackertool.SimTrackHits.Path.replace("sim", "digi"), - readoutName = "MUONHits", + inputSimHits = saveMuonBarrelTool.SimTrackHits.Path, + outputDigiHits = saveMuonBarrelTool.SimTrackHits.Path.replace("sim", "digi"), + readoutName = "MuonChamberBarrelReadout", xyResolution = 0.1, # mm zResolution = 1, # mm OutputLevel=DEBUG @@ -117,7 +126,7 @@ out.outputCommands = ["keep *"] import uuid -out.filename = "output_simplifiedDriftChamber_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+".root" +out.filename = "output_simpleMuonSystem_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+".root" #CPU information from Configurables import AuditorSvc, ChronoAuditor From 36fc06324b54d3b5ef5accb23086e9bfebec33e5 Mon Sep 17 00:00:00 2001 From: mahmoudali2 Date: Wed, 3 Apr 2024 19:07:13 +0200 Subject: [PATCH 3/6] Changing Muon Hits collection into cartisian --- MUONdigi/CMakeLists.txt | 1 + MUONdigi/include/MUONsimpleDigitizer.h | 18 ++++++----- MUONdigi/src/MUONsimpleDigitizer.cpp | 41 +++++++++++++++---------- MUONdigi/test/runMUONsimpleDigitizer.py | 33 ++++++++++---------- 4 files changed, 53 insertions(+), 40 deletions(-) diff --git a/MUONdigi/CMakeLists.txt b/MUONdigi/CMakeLists.txt index 81a4095..4dcb36c 100644 --- a/MUONdigi/CMakeLists.txt +++ b/MUONdigi/CMakeLists.txt @@ -11,6 +11,7 @@ file(GLOB headers gaudi_add_module(MUONdigi SOURCES ${sources} LINK + Gaudi::GaudiAlgLib Gaudi::GaudiKernel EDM4HEP::edm4hep k4FWCore::k4FWCore diff --git a/MUONdigi/include/MUONsimpleDigitizer.h b/MUONdigi/include/MUONsimpleDigitizer.h index c318ddc..317c803 100644 --- a/MUONdigi/include/MUONsimpleDigitizer.h +++ b/MUONdigi/include/MUONsimpleDigitizer.h @@ -24,6 +24,8 @@ namespace edm4hep { // DD4HEP #include "DD4hep/Detector.h" // for dd4hep::VolumeManager +#include "DDRec/Vector3D.h" +#include "DDRec/SurfaceManager.h" #include "DDSegmentation/BitFieldCoder.h" /** @class MUONsimpleDigitizer @@ -60,7 +62,7 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { DataHandle m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this}; // Detector readout name - Gaudi::Property m_readoutName{this, "readoutName", "MuonChamberBarrelReadout", "Name of the detector readout"}; + Gaudi::Property m_readoutName{this, "readoutName", "mRWELLChamberReadout", "Name of the detector readout"}; // Pointer to the geometry service ServiceHandle m_geoSvc; // Decoder for the cellID @@ -69,15 +71,15 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { dd4hep::VolumeManager m_volman; // z position resolution in mm - FloatProperty m_z_resolution{this, "zResolution", 1.0, - "Spatial resolution in the z direction [mm]"}; + FloatProperty m_x_resolution{this, "xResolution", 1.0, + "Spatial resolution in the x direction [mm]"}; // xy resolution in mm - FloatProperty m_xy_resolution{this, "xyResolution", 0.1, "Spatial resolution in the xy direction [mm]"}; + FloatProperty m_y_resolution{this, "yResolution", 1.0, "Spatial resolution in the y direction [mm]"}; // Random Number Service IRndmGenSvc* m_randSvc; - // Gaussian random number generator used for the smearing of the z position - Rndm::Numbers m_gauss_z; - // Gaussian random number generator used for the smearing of the xy position - Rndm::Numbers m_gauss_xy; + // Gaussian random number generator used for the smearing of the x position + Rndm::Numbers m_gauss_x; + // Gaussian random number generator used for the smearing of the y position + Rndm::Numbers m_gauss_y; }; diff --git a/MUONdigi/src/MUONsimpleDigitizer.cpp b/MUONdigi/src/MUONsimpleDigitizer.cpp index f0a128c..351750f 100644 --- a/MUONdigi/src/MUONsimpleDigitizer.cpp +++ b/MUONdigi/src/MUONsimpleDigitizer.cpp @@ -5,7 +5,7 @@ #include "DDRec/Vector3D.h" // ROOT -#include "Math/Cylindrical3D.h" +//#include "Math/Cylindrical3D.h" DECLARE_COMPONENT(MUONsimpleDigitizer) @@ -23,11 +23,11 @@ StatusCode MUONsimpleDigitizer::initialize() { error() << "Couldn't get RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } - if (m_gauss_z.initialize(m_randSvc, Rndm::Gauss(0., m_z_resolution)).isFailure()) { + if (m_gauss_x.initialize(m_randSvc, Rndm::Gauss(0., m_x_resolution)).isFailure()) { error() << "Couldn't initialize RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } - if (m_gauss_xy.initialize(m_randSvc, Rndm::Gauss(0., m_xy_resolution)).isFailure()) { + if (m_gauss_y.initialize(m_randSvc, Rndm::Gauss(0., m_y_resolution)).isFailure()) { error() << "Couldn't initialize RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } @@ -58,7 +58,7 @@ StatusCode MUONsimpleDigitizer::execute() { // retrieve the cell detElement dd4hep::DDSegmentation::CellID cellID = input_sim_hit.getCellID(); debug() << "Digitisation of " << m_readoutName << ", cellID: " << cellID << endmsg; - auto cellDetElement = m_volman.lookupDetElement(cellID); + // auto cellDetElement = m_volman.lookupDetElement(cellID); /* const std::string& stripsDetElementName = Form("system_%d_layer_%d_phi_%d", m_decoder->get(cellID, "layer"), m_decoder->get(cellID, "phi")); @@ -72,7 +72,10 @@ StatusCode MUONsimpleDigitizer::execute() { double simHitGlobalPosition[3] = {input_sim_hit.getPosition().x * dd4hep::mm, input_sim_hit.getPosition().y * dd4hep::mm, input_sim_hit.getPosition().z * dd4hep::mm}; + dd4hep::rec::Vector3D simHitGlobalPositionVector(simHitGlobalPosition[0], simHitGlobalPosition[1], + simHitGlobalPosition[2]); double simHitLocalPosition[3] = {0, 0, 0}; + // get the simHit coordinate in cm in the strips reference frame stripsTransformMatrix.MasterToLocal(simHitGlobalPosition, simHitLocalPosition); debug() << "Cell ID string: " << m_decoder->valueString(cellID) << endmsg; @@ -87,21 +90,27 @@ StatusCode MUONsimpleDigitizer::execute() { dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0], simHitLocalPosition[1], simHitLocalPosition[2]); // smear xy position - double smearedDistanceTostrips = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm; + double smearedX = simHitLocalPositionVector.x() ; + double smearedY = simHitLocalPositionVector.y() + m_gauss_x.shoot() * dd4hep::mm; // smear the z position - double smearedZ = simHitLocalPositionVector.z() + m_gauss_z.shoot() * dd4hep::mm; - // build the local position vector of the smeared hit using cylindrical coordinates. When we will have edm4hep::MUONHit there will be probably no need - ROOT::Math::Cylindrical3D digiHitLocalPositionVector(smearedDistanceTostrips, smearedZ, - simHitLocalPositionVector.phi()); - debug() << "Local simHit distance to the strips: " << simHitLocalPositionVector.rho() - << " Local digiHit distance to the strips: " << smearedDistanceTostrips << " in cm" << endmsg; + double smearedZ = simHitLocalPositionVector.z() + m_gauss_y.shoot() * dd4hep::mm; + + double digiHitLocalPosition[3] = {smearedX, smearedY, smearedZ}; + + // build the local position vector of the smeared hit. + // dd4hep::rec::Vector3D digiHitLocalPositionVector(smearedX, smearedY, smearedZ); + debug() << "Local simHit x: " << simHitLocalPositionVector.x() + << " Local digiHit x: " << smearedX << " in cm" << endmsg; + debug() << "Local simHit y: " << simHitLocalPositionVector.y() + << " Local digiHit y: " << smearedY << " in cm" << endmsg; debug() << "Local simHit z: " << simHitLocalPositionVector.z() - << " Local digiHit distance to the strips: " << smearedZ << " in cm" << endmsg; - debug() << "Local simHit phi: " << simHitLocalPositionVector.phi() - << " Local digiHit distance to the strips: " << digiHitLocalPositionVector.Phi() << endmsg; + << " Local digiHit z: " << smearedZ << " in cm" << endmsg; + + + // go back to the global frame - double digiHitLocalPosition[3] = {digiHitLocalPositionVector.x(), digiHitLocalPositionVector.y(), - digiHitLocalPositionVector.z()}; + // double digiHitLocalPosition[3] = {digiHitLocalPositionVector.x(), digiHitLocalPositionVector.y(), + // digiHitLocalPositionVector.z()}; double digiHitGlobalPosition[3] = {0, 0, 0}; stripsTransformMatrix.LocalToMaster(digiHitLocalPosition, digiHitGlobalPosition); // go back to mm diff --git a/MUONdigi/test/runMUONsimpleDigitizer.py b/MUONdigi/test/runMUONsimpleDigitizer.py index a952a39..edffa0f 100644 --- a/MUONdigi/test/runMUONsimpleDigitizer.py +++ b/MUONdigi/test/runMUONsimpleDigitizer.py @@ -43,11 +43,12 @@ from Configurables import GeoSvc geoservice = GeoSvc("GeoSvc") # if FCC_DETECTORS is empty, this should use relative path to working directory -path_to_detector = os.environ.get("FCCDETECTORS", "") +path_to_detector = os.environ.get("K4GEO", "") print(path_to_detector) detectors_to_use=[ # 'Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v02.xml' - '/afs/cern.ch/work/m/maali/public/FCCDetectors/Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v02.xml' + '/afs/cern.ch/work/m/maali/public/k4geo/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml' + # '/afs/cern.ch/work/m/maali/public/FCCDetectors/Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v02.xml' ] # prefix all xmls with path_to_detector geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use] @@ -74,7 +75,7 @@ # Fixed seed to have reproducible results, change it for each job if you split one production into several jobs # Mind that if you leave Gaudi handle random seed and some job start within the same second (very likely) you will have duplicates geantservice.randomNumbersFromGaudi = False -geantservice.seedValue = 4242 +geantservice.seedValue = 3135 # Range cut geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] @@ -92,18 +93,18 @@ #savetrackertool = SimG4SaveTrackerHits("SimG4SaveTrackerHits", readoutNames=["MuonChamberBarrelReadout"]) #savetrackertool.SimTrackHits.Path = "MUON_simTrackerHits" -saveMuonBarrelTool = SimG4SaveTrackerHits("SimG4SaveMuonBarrelHits", readoutName="MuonChamberBarrelReadout") -saveMuonBarrelTool.SimTrackHits.Path = "muonBarrelSimHits" +SimG4SaveMuonHits = SimG4SaveTrackerHits("SimG4SaveMuonHits", readoutName="mRWELLChamberReadout") +SimG4SaveMuonHits.SimTrackHits.Path = "Muon_SimHits" -saveMuonPositiveEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonPositiveEndcapHits", readoutName="MuonChamberPositiveEndcapReadout") -saveMuonPositiveEndcapTool.SimTrackHits.Path = "muonPositiveEndcapSimHits" +#saveMuonPositiveEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonPositiveEndcapHits", readoutName="MuonChamberPositiveEndcapReadout") +#saveMuonPositiveEndcapTool.SimTrackHits.Path = "muonPositiveEndcapSimHits" -saveMuonNegativeEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonNegativeEndcapHits", readoutName="MuonChamberNegativeEndcapReadout") -saveMuonNegativeEndcapTool.SimTrackHits.Path = "muonNegativeEndcapSimHits" +#saveMuonNegativeEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonNegativeEndcapHits", readoutName="MuonChamberNegativeEndcapReadout") +#saveMuonNegativeEndcapTool.SimTrackHits.Path = "muonNegativeEndcapSimHits" from Configurables import SimG4Alg geantsim = SimG4Alg("SimG4Alg", - outputs= [saveMuonBarrelTool, saveMuonPositiveEndcapTool, saveMuonNegativeEndcapTool + outputs= [SimG4SaveMuonHits, #saveHistTool ], eventProvider=particle_converter, @@ -111,12 +112,12 @@ # Digitize tracker hits from Configurables import MUONsimpleDigitizer muon_digitizer = MUONsimpleDigitizer("MUONsimpleDigitizer", - inputSimHits = saveMuonBarrelTool.SimTrackHits.Path, - outputDigiHits = saveMuonBarrelTool.SimTrackHits.Path.replace("sim", "digi"), - readoutName = "MuonChamberBarrelReadout", - xyResolution = 0.1, # mm - zResolution = 1, # mm - OutputLevel=DEBUG + inputSimHits = SimG4SaveMuonHits.SimTrackHits.Path, + outputDigiHits = SimG4SaveMuonHits.SimTrackHits.Path.replace("sim", "digi"), + readoutName = "mRWELLChamberReadout", + xResolution = 0.4, # mm + yResolution = 0.4, # mm + OutputLevel=INFO ) ################ Output From d013b40a41aabdb5860005d56e61d2fe96617f84 Mon Sep 17 00:00:00 2001 From: mahmoudali2 Date: Tue, 9 Jul 2024 22:12:32 +0200 Subject: [PATCH 4/6] Updating muon-system readout --- MUONdigi/include/MUONsimpleDigitizer.h | 2 +- MUONdigi/test/runMUONsimpleDigitizer.py | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/MUONdigi/include/MUONsimpleDigitizer.h b/MUONdigi/include/MUONsimpleDigitizer.h index 317c803..2c27a8a 100644 --- a/MUONdigi/include/MUONsimpleDigitizer.h +++ b/MUONdigi/include/MUONsimpleDigitizer.h @@ -62,7 +62,7 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { DataHandle m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this}; // Detector readout name - Gaudi::Property m_readoutName{this, "readoutName", "mRWELLChamberReadout", "Name of the detector readout"}; + Gaudi::Property m_readoutName{this, "readoutName", "MuonSystemCollection", "Name of the detector readout"}; // Pointer to the geometry service ServiceHandle m_geoSvc; // Decoder for the cellID diff --git a/MUONdigi/test/runMUONsimpleDigitizer.py b/MUONdigi/test/runMUONsimpleDigitizer.py index edffa0f..e5f6985 100644 --- a/MUONdigi/test/runMUONsimpleDigitizer.py +++ b/MUONdigi/test/runMUONsimpleDigitizer.py @@ -47,7 +47,7 @@ print(path_to_detector) detectors_to_use=[ # 'Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v02.xml' - '/afs/cern.ch/work/m/maali/public/k4geo/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml' + '/afs/cern.ch/work/m/maali/public/k4geo/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml' # '/afs/cern.ch/work/m/maali/public/FCCDetectors/Detector/DetFCCeeIDEA/compact/IDEA_o1_v01/FCCee_DectMaster_v02.xml' ] # prefix all xmls with path_to_detector @@ -93,7 +93,7 @@ #savetrackertool = SimG4SaveTrackerHits("SimG4SaveTrackerHits", readoutNames=["MuonChamberBarrelReadout"]) #savetrackertool.SimTrackHits.Path = "MUON_simTrackerHits" -SimG4SaveMuonHits = SimG4SaveTrackerHits("SimG4SaveMuonHits", readoutName="mRWELLChamberReadout") +SimG4SaveMuonHits = SimG4SaveTrackerHits("SimG4SaveMuonHits", readoutName="MuonSystemCollection") SimG4SaveMuonHits.SimTrackHits.Path = "Muon_SimHits" #saveMuonPositiveEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonPositiveEndcapHits", readoutName="MuonChamberPositiveEndcapReadout") @@ -113,8 +113,8 @@ from Configurables import MUONsimpleDigitizer muon_digitizer = MUONsimpleDigitizer("MUONsimpleDigitizer", inputSimHits = SimG4SaveMuonHits.SimTrackHits.Path, - outputDigiHits = SimG4SaveMuonHits.SimTrackHits.Path.replace("sim", "digi"), - readoutName = "mRWELLChamberReadout", + outputDigiHits = SimG4SaveMuonHits.SimTrackHits.Path.replace("Sim", "Digi"), + readoutName = "MuonSystemCollection", xResolution = 0.4, # mm yResolution = 0.4, # mm OutputLevel=INFO @@ -154,7 +154,7 @@ out ], EvtSel = 'NONE', - EvtMax = 100, + EvtMax = 10000, ExtSvc = [geoservice, podioevent, geantservice, audsvc], StopOnSignal = True, ) From a5a4acecd3ed9dbc68f526360c1bd379b6fe5ab3 Mon Sep 17 00:00:00 2001 From: mahmoudali2 Date: Thu, 18 Jul 2024 15:40:52 +0200 Subject: [PATCH 5/6] Defining detector efficiency and adding Association extension for the MUONDigi --- MUONdigi/CMakeLists.txt | 44 +++++++++++++++++-- .../dataFormatExtension/muonSystemHit.yaml | 19 ++++++++ MUONdigi/include/MUONsimpleDigitizer.h | 18 ++++++-- MUONdigi/src/MUONsimpleDigitizer.cpp | 36 +++++++++------ MUONdigi/test/runMUONsimpleDigitizer.py | 27 ++++-------- 5 files changed, 106 insertions(+), 38 deletions(-) create mode 100644 MUONdigi/dataFormatExtension/muonSystemHit.yaml diff --git a/MUONdigi/CMakeLists.txt b/MUONdigi/CMakeLists.txt index 4dcb36c..5a56ae1 100644 --- a/MUONdigi/CMakeLists.txt +++ b/MUONdigi/CMakeLists.txt @@ -1,4 +1,38 @@ -project(MUONdigi) +set(PackageName MUONdigi) + +project(${PackageName}) + +# Build the extension data model and link it against the upstream model +PODIO_GENERATE_DATAMODEL(muon_extension dataFormatExtension/muonSystemHit.yaml ext_headers ext_sources + UPSTREAM_EDM edm4hep:${EDM4HEP_DATA_DIR}/edm4hep.yaml + IO_BACKEND_HANDLERS ${PODIO_IO_HANDLERS} + OUTPUT_FOLDER ${CMAKE_CURRENT_BINARY_DIR}) +PODIO_ADD_DATAMODEL_CORE_LIB(muon_extension "${ext_headers}" "${ext_sources}" + OUTPUT_FOLDER ${CMAKE_CURRENT_BINARY_DIR}) +target_link_libraries(muon_extension PUBLIC EDM4HEP::edm4hep) +PODIO_ADD_ROOT_IO_DICT(muon_extensionDict muon_extension "${ext_headers}" src/selection.xml + OUTPUT_FOLDER ${CMAKE_CURRENT_BINARY_DIR}) +#PODIO_ADD_SIO_IO_BLOCKS(muon_extension "${ext_headers}" "${ext_sources}") +add_library(muon_extension::muon_extensionDict ALIAS muon_extensionDict) +list(APPEND EXTENSION_INSTALL_LIBS muon_extension muon_extensionDict) +install(TARGETS ${EXTENSION_INSTALL_LIBS} + EXPORT ${PROJECT_NAME}Targets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/muon_extension" + COMPONENT dev) + +install(FILES + "${PROJECT_BINARY_DIR}/muon_extensionDictDict.rootmap" + DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT dev) + +install(FILES + dataFormatExtension/muonSystemHit.yaml + DESTINATION "${CMAKE_INSTALL_DATADIR}/muon_extension" COMPONENT dev) + +install(FILES + "${PROJECT_BINARY_DIR}/libmuon_extensionDict_rdict.pcm" + DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT dev) file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cpp @@ -8,19 +42,21 @@ file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.h ) -gaudi_add_module(MUONdigi +gaudi_add_module(${PackageName} SOURCES ${sources} LINK Gaudi::GaudiAlgLib Gaudi::GaudiKernel EDM4HEP::edm4hep + muon_extensionDict k4FWCore::k4FWCore DD4hep::DDRec ) -target_include_directories(MUONdigi PUBLIC +target_include_directories(${PackageName} PUBLIC $ $ + ${CMAKE_INSTALL_PREFIX}/muon_extension ) set_target_properties(MUONdigi PROPERTIES PUBLIC_HEADER "${headers}") @@ -41,4 +77,4 @@ install(TARGETS MUONdigi install(FILES ${scripts} DESTINATION test) SET(test_name "test_MUONsimpleDigitizer") -ADD_TEST(NAME t_${test_name} COMMAND k4run test/runMUONsimpleDigitizer.py) +ADD_TEST(NAME t_${test_name} COMMAND k4run test/runMUONsimpleDigitizer.py) \ No newline at end of file diff --git a/MUONdigi/dataFormatExtension/muonSystemHit.yaml b/MUONdigi/dataFormatExtension/muonSystemHit.yaml new file mode 100644 index 0000000..855a0bc --- /dev/null +++ b/MUONdigi/dataFormatExtension/muonSystemHit.yaml @@ -0,0 +1,19 @@ +--- +schema_version: 1 +options: + # should getters / setters be prefixed with get / set? + getSyntax: True + # should POD members be exposed with getters/setters in classes that have them as members? + exposePODMembers: False + includeSubfolder: True + +datatypes: + + extension::MCRecoMuonSystemDigiAssociation: + Description: "Association between a Digi hits and the corresponding simulated hit for the muon system" + Author: "B. Francois, CERN; Mahmoud Ali, INFN-Bo" + Members: + - float weight // weight of this association + OneToOneRelations: + - edm4hep::TrackerHit3D digi // reference to the digitized hit + - edm4hep::SimTrackerHit sim // reference to the simulated hit \ No newline at end of file diff --git a/MUONdigi/include/MUONsimpleDigitizer.h b/MUONdigi/include/MUONsimpleDigitizer.h index 2c27a8a..c70ff72 100644 --- a/MUONdigi/include/MUONsimpleDigitizer.h +++ b/MUONdigi/include/MUONsimpleDigitizer.h @@ -17,11 +17,16 @@ #include "edm4hep/TrackerHit3DCollection.h" #else #include "edm4hep/TrackerHitCollection.h" + +#include "podio/UserDataCollection.h" + namespace edm4hep { using TrackerHit3DCollection = edm4hep::TrackerHitCollection; } // namespace edm4hep #endif +#include "extension/MCRecoMuonSystemDigiAssociationCollection.h" + // DD4HEP #include "DD4hep/Detector.h" // for dd4hep::VolumeManager #include "DDRec/Vector3D.h" @@ -60,6 +65,8 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { DataHandle m_input_sim_hits{"inputSimHits", Gaudi::DataHandle::Reader, this}; // Output digitized tracker hit collection name DataHandle m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this}; + // Output association between digitized and simulated hit collections name + DataHandle m_output_sim_digi_association{"outputSimDigiAssociation", Gaudi::DataHandle::Writer, this}; // Detector readout name Gaudi::Property m_readoutName{this, "readoutName", "MuonSystemCollection", "Name of the detector readout"}; @@ -70,16 +77,21 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { // Volume manager to get the physical cell sensitive volume dd4hep::VolumeManager m_volman; - // z position resolution in mm + // x position resolution in mm FloatProperty m_x_resolution{this, "xResolution", 1.0, "Spatial resolution in the x direction [mm]"}; - // xy resolution in mm + // y resolution in mm FloatProperty m_y_resolution{this, "yResolution", 1.0, "Spatial resolution in the y direction [mm]"}; + // Detector efficiency + FloatProperty m_efficiency{this, "efficiency", 0.95, "Detector efficiency"}; + // Random Number Service IRndmGenSvc* m_randSvc; // Gaussian random number generator used for the smearing of the x position Rndm::Numbers m_gauss_x; // Gaussian random number generator used for the smearing of the y position Rndm::Numbers m_gauss_y; -}; + // Flat random number generator used for efficiency + Rndm::Numbers m_flat; +}; \ No newline at end of file diff --git a/MUONdigi/src/MUONsimpleDigitizer.cpp b/MUONdigi/src/MUONsimpleDigitizer.cpp index 351750f..d068940 100644 --- a/MUONdigi/src/MUONsimpleDigitizer.cpp +++ b/MUONdigi/src/MUONsimpleDigitizer.cpp @@ -4,6 +4,7 @@ #include "DD4hep/Detector.h" #include "DDRec/Vector3D.h" +#include "extension/MCRecoMuonSystemDigiAssociationCollection.h" // ROOT //#include "Math/Cylindrical3D.h" @@ -13,6 +14,8 @@ MUONsimpleDigitizer::MUONsimpleDigitizer(const std::string& aName, ISvcLocator* : GaudiAlgorithm(aName, aSvcLoc), m_geoSvc("GeoSvc", "MUONsimpleDigitizer") { declareProperty("inputSimHits", m_input_sim_hits, "Input sim tracker hit collection name"); declareProperty("outputDigiHits", m_output_digi_hits, "Output digitized tracker hit collection name"); + declareProperty("outputSimDigiAssociation", m_output_sim_digi_association, "Output name for the association between digitized and simulated hit collections"); + declareProperty("efficiency", m_efficiency, "Efficiency of the detector"); } MUONsimpleDigitizer::~MUONsimpleDigitizer() {} @@ -31,6 +34,10 @@ StatusCode MUONsimpleDigitizer::initialize() { error() << "Couldn't initialize RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } + if (m_flat.initialize(m_randSvc, Rndm::Flat(0., 1.)).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } // check if readout exists if (m_geoSvc->getDetector()->readouts().find(m_readoutName) == m_geoSvc->getDetector()->readouts().end()) { @@ -50,25 +57,27 @@ StatusCode MUONsimpleDigitizer::execute() { const edm4hep::SimTrackerHitCollection* input_sim_hits = m_input_sim_hits.get(); debug() << "Input Sim Hit collection size: " << input_sim_hits->size() << endmsg; + // Prepare a collection for the association between digitized and simulated hit, setting weights to 1 + extension::MCRecoMuonSystemDigiAssociationCollection* digi_sim_associations = m_output_sim_digi_association.createAndPut(); + // Digitize the sim hits edm4hep::TrackerHit3DCollection* output_digi_hits = m_output_digi_hits.createAndPut(); for (const auto& input_sim_hit : *input_sim_hits) { + // Apply efficiency + if (m_flat.shoot() > m_efficiency) { + continue; // Skip this hit + } + auto output_digi_hit = output_digi_hits->create(); // smear the hit position: // retrieve the cell detElement dd4hep::DDSegmentation::CellID cellID = input_sim_hit.getCellID(); debug() << "Digitisation of " << m_readoutName << ", cellID: " << cellID << endmsg; - // auto cellDetElement = m_volman.lookupDetElement(cellID); - - /* const std::string& stripsDetElementName = - Form("system_%d_layer_%d_phi_%d", m_decoder->get(cellID, "layer"), m_decoder->get(cellID, "phi")); - dd4hep::DetElement stripsDetElement = cellDetElement.child(stripsDetElementName); - const auto& stripsTransformMatrix = stripsDetElement.nominal().worldTransformation(); - */ + // auto cellDetElement = m_volman.lookupDetElement(cellID); const auto& stripsTransformMatrix = m_volman.lookupVolumePlacement(cellID).matrix(); - // Retrieve global position in mm and apply unit transformation (translation matrix is tored in cm) + // Retrieve global position in mm and apply unit transformation (translation matrix is stored in cm) double simHitGlobalPosition[3] = {input_sim_hit.getPosition().x * dd4hep::mm, input_sim_hit.getPosition().y * dd4hep::mm, input_sim_hit.getPosition().z * dd4hep::mm}; @@ -98,7 +107,7 @@ StatusCode MUONsimpleDigitizer::execute() { double digiHitLocalPosition[3] = {smearedX, smearedY, smearedZ}; // build the local position vector of the smeared hit. - // dd4hep::rec::Vector3D digiHitLocalPositionVector(smearedX, smearedY, smearedZ); + // dd4hep::rec::Vector3D digiHitLocalPositionVector(smearedX, smearedY, smearedZ); debug() << "Local simHit x: " << simHitLocalPositionVector.x() << " Local digiHit x: " << smearedX << " in cm" << endmsg; debug() << "Local simHit y: " << simHitLocalPositionVector.y() @@ -106,11 +115,12 @@ StatusCode MUONsimpleDigitizer::execute() { debug() << "Local simHit z: " << simHitLocalPositionVector.z() << " Local digiHit z: " << smearedZ << " in cm" << endmsg; - + // create the association between digitized and simulated hit + auto digi_sim_association = digi_sim_associations->create(); + digi_sim_association.setDigi(output_digi_hit); + digi_sim_association.setSim(input_sim_hit); // go back to the global frame - // double digiHitLocalPosition[3] = {digiHitLocalPositionVector.x(), digiHitLocalPositionVector.y(), - // digiHitLocalPositionVector.z()}; double digiHitGlobalPosition[3] = {0, 0, 0}; stripsTransformMatrix.LocalToMaster(digiHitLocalPosition, digiHitGlobalPosition); // go back to mm @@ -124,4 +134,4 @@ StatusCode MUONsimpleDigitizer::execute() { return StatusCode::SUCCESS; } -StatusCode MUONsimpleDigitizer::finalize() { return StatusCode::SUCCESS; } +StatusCode MUONsimpleDigitizer::finalize() { return StatusCode::SUCCESS; } \ No newline at end of file diff --git a/MUONdigi/test/runMUONsimpleDigitizer.py b/MUONdigi/test/runMUONsimpleDigitizer.py index e5f6985..627fc53 100644 --- a/MUONdigi/test/runMUONsimpleDigitizer.py +++ b/MUONdigi/test/runMUONsimpleDigitizer.py @@ -8,8 +8,6 @@ from GaudiKernel.SystemOfUnits import MeV, GeV, tesla ################## Particle gun setup momentum = 5 # in GeV -#thetaMin = 90.25 # degrees -#thetaMax = 90.25 # degrees thetaMin = 0 # degrees thetaMax = 180 # degrees pdgCode = 13 # 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ @@ -75,7 +73,7 @@ # Fixed seed to have reproducible results, change it for each job if you split one production into several jobs # Mind that if you leave Gaudi handle random seed and some job start within the same second (very likely) you will have duplicates geantservice.randomNumbersFromGaudi = False -geantservice.seedValue = 3135 +geantservice.seedValue = 4242 # Range cut geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] @@ -90,33 +88,25 @@ particle_converter.GenParticles.Path = genParticlesOutputName from Configurables import SimG4SaveTrackerHits -#savetrackertool = SimG4SaveTrackerHits("SimG4SaveTrackerHits", readoutNames=["MuonChamberBarrelReadout"]) -#savetrackertool.SimTrackHits.Path = "MUON_simTrackerHits" - SimG4SaveMuonHits = SimG4SaveTrackerHits("SimG4SaveMuonHits", readoutName="MuonSystemCollection") SimG4SaveMuonHits.SimTrackHits.Path = "Muon_SimHits" -#saveMuonPositiveEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonPositiveEndcapHits", readoutName="MuonChamberPositiveEndcapReadout") -#saveMuonPositiveEndcapTool.SimTrackHits.Path = "muonPositiveEndcapSimHits" - -#saveMuonNegativeEndcapTool = SimG4SaveTrackerHits("SimG4SaveMuonNegativeEndcapHits", readoutName="MuonChamberNegativeEndcapReadout") -#saveMuonNegativeEndcapTool.SimTrackHits.Path = "muonNegativeEndcapSimHits" - from Configurables import SimG4Alg geantsim = SimG4Alg("SimG4Alg", - outputs= [SimG4SaveMuonHits, - #saveHistTool - ], + outputs= [SimG4SaveMuonHits], eventProvider=particle_converter, OutputLevel=INFO) + # Digitize tracker hits from Configurables import MUONsimpleDigitizer muon_digitizer = MUONsimpleDigitizer("MUONsimpleDigitizer", inputSimHits = SimG4SaveMuonHits.SimTrackHits.Path, outputDigiHits = SimG4SaveMuonHits.SimTrackHits.Path.replace("Sim", "Digi"), + outputSimDigiAssociation = "MS_simDigiAssociation", readoutName = "MuonSystemCollection", xResolution = 0.4, # mm yResolution = 0.4, # mm + efficiency = 0.95, # efficiency of the detector (mRWELL chamber) OutputLevel=INFO ) @@ -127,7 +117,7 @@ out.outputCommands = ["keep *"] import uuid -out.filename = "output_simpleMuonSystem_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+".root" +out.filename = "output_MuonSystemDigi_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+".root" #CPU information from Configurables import AuditorSvc, ChronoAuditor @@ -137,6 +127,7 @@ genAlg.AuditExecute = True hepmc_converter.AuditExecute = True geantsim.AuditExecute = True +muon_digitizer.AuditExecute = True out.AuditExecute = True from Configurables import EventCounter @@ -154,7 +145,7 @@ out ], EvtSel = 'NONE', - EvtMax = 10000, + EvtMax = 100, ExtSvc = [geoservice, podioevent, geantservice, audsvc], StopOnSignal = True, - ) +) \ No newline at end of file From 7bdd29f16e01b06e84836d402649c2a3d3b646a0 Mon Sep 17 00:00:00 2001 From: mahmoudali2 Date: Thu, 25 Jul 2024 17:49:17 +0200 Subject: [PATCH 6/6] Adding validation for the MUON digitizer --- .../dataFormatExtension/muonSystemHit.yaml | 4 +- MUONdigi/include/MUONsimpleDigitizer.h | 16 ++++- MUONdigi/src/MUONsimpleDigitizer.cpp | 66 ++++++++++++++++--- MUONdigi/test/runMUONsimpleDigitizer.py | 7 +- 4 files changed, 76 insertions(+), 17 deletions(-) diff --git a/MUONdigi/dataFormatExtension/muonSystemHit.yaml b/MUONdigi/dataFormatExtension/muonSystemHit.yaml index 855a0bc..8ed466d 100644 --- a/MUONdigi/dataFormatExtension/muonSystemHit.yaml +++ b/MUONdigi/dataFormatExtension/muonSystemHit.yaml @@ -14,6 +14,8 @@ datatypes: Author: "B. Francois, CERN; Mahmoud Ali, INFN-Bo" Members: - float weight // weight of this association + - edm4hep::Vector3d positionDifference // the difference between simHit and DigiHit positions. OneToOneRelations: - edm4hep::TrackerHit3D digi // reference to the digitized hit - - edm4hep::SimTrackerHit sim // reference to the simulated hit \ No newline at end of file + - edm4hep::SimTrackerHit sim // reference to the simulated hit + \ No newline at end of file diff --git a/MUONdigi/include/MUONsimpleDigitizer.h b/MUONdigi/include/MUONsimpleDigitizer.h index c70ff72..9e1c2a6 100644 --- a/MUONdigi/include/MUONsimpleDigitizer.h +++ b/MUONdigi/include/MUONsimpleDigitizer.h @@ -6,9 +6,11 @@ #include "GaudiKernel/IRndmGenSvc.h" #include "GaudiKernel/RndmGenerators.h" -// K4FWCORE +// K4FWCORE & podio #include "k4FWCore/DataHandle.h" #include "k4Interface/IGeoSvc.h" +#include "podio/UserDataCollection.h" +#include "edm4hep/Vector3d.h" // EDM4HEP #include "edm4hep/SimTrackerHitCollection.h" @@ -18,8 +20,6 @@ #else #include "edm4hep/TrackerHitCollection.h" -#include "podio/UserDataCollection.h" - namespace edm4hep { using TrackerHit3DCollection = edm4hep::TrackerHitCollection; } // namespace edm4hep @@ -83,15 +83,25 @@ class MUONsimpleDigitizer : public GaudiAlgorithm { // y resolution in mm FloatProperty m_y_resolution{this, "yResolution", 1.0, "Spatial resolution in the y direction [mm]"}; + // z resolution in mm + FloatProperty m_z_resolution{this, "zResolution", 1.0, "Spatial resolution in the z direction [mm]"}; + // Detector efficiency FloatProperty m_efficiency{this, "efficiency", 0.95, "Detector efficiency"}; + // Declaration of validation distribution + DataHandle> m_simDigiDifferenceX{"simDigiDifferenceX", Gaudi::DataHandle::Writer, this}; // mm + DataHandle> m_simDigiDifferenceY{"simDigiDifferenceY", Gaudi::DataHandle::Writer, this}; // mm + DataHandle> m_simDigiDifferenceZ{"simDigiDifferenceZ", Gaudi::DataHandle::Writer, this}; // mm + // Random Number Service IRndmGenSvc* m_randSvc; // Gaussian random number generator used for the smearing of the x position Rndm::Numbers m_gauss_x; // Gaussian random number generator used for the smearing of the y position Rndm::Numbers m_gauss_y; + // Gaussian random number generator used for the smearing of the z position + Rndm::Numbers m_gauss_z; // Flat random number generator used for efficiency Rndm::Numbers m_flat; }; \ No newline at end of file diff --git a/MUONdigi/src/MUONsimpleDigitizer.cpp b/MUONdigi/src/MUONsimpleDigitizer.cpp index d068940..4f0f131 100644 --- a/MUONdigi/src/MUONsimpleDigitizer.cpp +++ b/MUONdigi/src/MUONsimpleDigitizer.cpp @@ -3,6 +3,8 @@ // DD4hep #include "DD4hep/Detector.h" #include "DDRec/Vector3D.h" +#include "podio/UserDataCollection.h" +#include "edm4hep/Vector3d.h" #include "extension/MCRecoMuonSystemDigiAssociationCollection.h" // ROOT @@ -34,6 +36,10 @@ StatusCode MUONsimpleDigitizer::initialize() { error() << "Couldn't initialize RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } + if (m_gauss_z.initialize(m_randSvc, Rndm::Gauss(0., m_z_resolution)).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } if (m_flat.initialize(m_randSvc, Rndm::Flat(0., 1.)).isFailure()) { error() << "Couldn't initialize RndmGenSvc!" << endmsg; return StatusCode::FAILURE; @@ -57,11 +63,17 @@ StatusCode MUONsimpleDigitizer::execute() { const edm4hep::SimTrackerHitCollection* input_sim_hits = m_input_sim_hits.get(); debug() << "Input Sim Hit collection size: " << input_sim_hits->size() << endmsg; + // Digitize the sim hits + edm4hep::TrackerHit3DCollection* output_digi_hits = m_output_digi_hits.createAndPut(); + // Prepare a collection for the association between digitized and simulated hit, setting weights to 1 extension::MCRecoMuonSystemDigiAssociationCollection* digi_sim_associations = m_output_sim_digi_association.createAndPut(); - // Digitize the sim hits - edm4hep::TrackerHit3DCollection* output_digi_hits = m_output_digi_hits.createAndPut(); + // Prepare collections for the validation distributions + auto simDigiDifferenceX = m_simDigiDifferenceX.createAndPut(); + auto simDigiDifferenceY = m_simDigiDifferenceY.createAndPut(); + auto simDigiDifferenceZ = m_simDigiDifferenceZ.createAndPut(); + for (const auto& input_sim_hit : *input_sim_hits) { // Apply efficiency if (m_flat.shoot() > m_efficiency) { @@ -99,10 +111,10 @@ StatusCode MUONsimpleDigitizer::execute() { dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0], simHitLocalPosition[1], simHitLocalPosition[2]); // smear xy position - double smearedX = simHitLocalPositionVector.x() ; - double smearedY = simHitLocalPositionVector.y() + m_gauss_x.shoot() * dd4hep::mm; + double smearedX = simHitLocalPositionVector.x() + m_gauss_x.shoot() * dd4hep::mm; + double smearedY = simHitLocalPositionVector.y() + m_gauss_y.shoot() * dd4hep::mm; // smear the z position - double smearedZ = simHitLocalPositionVector.z() + m_gauss_y.shoot() * dd4hep::mm; + double smearedZ = simHitLocalPositionVector.z() + m_gauss_z.shoot() * dd4hep::mm; double digiHitLocalPosition[3] = {smearedX, smearedY, smearedZ}; @@ -114,11 +126,6 @@ StatusCode MUONsimpleDigitizer::execute() { << " Local digiHit y: " << smearedY << " in cm" << endmsg; debug() << "Local simHit z: " << simHitLocalPositionVector.z() << " Local digiHit z: " << smearedZ << " in cm" << endmsg; - - // create the association between digitized and simulated hit - auto digi_sim_association = digi_sim_associations->create(); - digi_sim_association.setDigi(output_digi_hit); - digi_sim_association.setSim(input_sim_hit); // go back to the global frame double digiHitGlobalPosition[3] = {0, 0, 0}; @@ -129,9 +136,48 @@ StatusCode MUONsimpleDigitizer::execute() { digiHitGlobalPosition[2] / dd4hep::mm); output_digi_hit.setPosition(digiHitGlobalPositionVector); output_digi_hit.setCellID(cellID); + + // create the association between digitized and simulated hit + auto digi_sim_association = digi_sim_associations->create(); + digi_sim_association.setDigi(output_digi_hit); + digi_sim_association.setSim(input_sim_hit); + + // Validating if the digitization; Calculate the position difference in global coordinates + double dx = (digiHitGlobalPosition[0] - simHitGlobalPosition[0]) / dd4hep::mm; + double dy = (digiHitGlobalPosition[1] - simHitGlobalPosition[1]) / dd4hep::mm; + double dz = (digiHitGlobalPosition[2] - simHitGlobalPosition[2]) / dd4hep::mm; + + edm4hep::Vector3d simDigiHitsPositionsDifference(dx / dd4hep::mm, + dy / dd4hep::mm, + dz / dd4hep::mm); + simDigiDifferenceX->push_back(dx); + simDigiDifferenceY->push_back(dy); + simDigiDifferenceZ->push_back(dz); + //digi_sim_association.setPositionDifference(simDigiHitsPositionsDifference); } debug() << "Output Digi Hit collection size: " << output_digi_hits->size() << endmsg; return StatusCode::SUCCESS; +/* + // Validating if the digitization is working + for (const auto& association : *digi_sim_associations) { + const auto simHit = association.getSim(); + const auto digiHit = association.getDigi(); + + const auto simPosition = simHit.getPosition(); + const auto digiPosition = digiHit.getPosition(); + + double dx = digiPosition.x - simPosition.x; + double dy = digiPosition.y - simPosition.y; + double dz = digiPosition.z - simPosition.z; + + edm4hep::Vector3d simDigiHitsPositionsDifference(dx / dd4hep::mm, + dy / dd4hep::mm, + dz / dd4hep::mm); + association.setPositionDifference(simDigiHitsPositionsDifference); + + debug() << "Hit position Difference for cellID " << simHit.getCellID() << ": dx = " << dx << " mm, dy = " << dy << " mm, dz = " << dz << " mm" << endmsg; + } +*/ } StatusCode MUONsimpleDigitizer::finalize() { return StatusCode::SUCCESS; } \ No newline at end of file diff --git a/MUONdigi/test/runMUONsimpleDigitizer.py b/MUONdigi/test/runMUONsimpleDigitizer.py index 627fc53..4b33c8a 100644 --- a/MUONdigi/test/runMUONsimpleDigitizer.py +++ b/MUONdigi/test/runMUONsimpleDigitizer.py @@ -106,8 +106,9 @@ readoutName = "MuonSystemCollection", xResolution = 0.4, # mm yResolution = 0.4, # mm + zResolution = 0.4, # mm efficiency = 0.95, # efficiency of the detector (mRWELL chamber) - OutputLevel=INFO + OutputLevel=DEBUG ) ################ Output @@ -117,7 +118,7 @@ out.outputCommands = ["keep *"] import uuid -out.filename = "output_MuonSystemDigi_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+".root" +out.filename = "output_MuonSystemDigi"+"_pMin_"+str(momentum)+"_GeV"+"_pdgId_"+str(pdgCode)+".root" #CPU information from Configurables import AuditorSvc, ChronoAuditor @@ -145,7 +146,7 @@ out ], EvtSel = 'NONE', - EvtMax = 100, + EvtMax = 10000, ExtSvc = [geoservice, podioevent, geantservice, audsvc], StopOnSignal = True, ) \ No newline at end of file