diff --git a/VTXdigi/include/VTXdigitizer.h b/VTXdigi/include/VTXdigitizer.h index e0cf0de..84515e5 100644 --- a/VTXdigi/include/VTXdigitizer.h +++ b/VTXdigi/include/VTXdigitizer.h @@ -30,11 +30,13 @@ namespace edm4hep { #include "DDSegmentation/BitFieldCoder.h" +#include + /** @class VTXdigitizer * * Algorithm for creating digitized (meaning 'reconstructed' for now) vertex detector hits (edm4hep::TrackerHit3D) from Geant4 hits (edm4hep::SimTrackerHit). * - * @author Brieuc Francois + * @author Brieuc Francois, Armin Ilg * @date 2023-03 * */ @@ -76,13 +78,13 @@ class VTXdigitizer : public Gaudi::Algorithm { dd4hep::VolumeManager m_volman; // x resolution in um - FloatProperty m_x_resolution{this, "xResolution", 0.1, "Spatial resolution in the x direction [um] (r-phi direction in barrel, z direction in disks)"}; + Gaudi::Property> m_x_resolution{this, "xResolution", {0.1}, "Spatial resolutions in the x direction per layer [um] (r-phi direction in barrel, z direction in disks)"}; // y resolution in um - FloatProperty m_y_resolution{this, "yResolution", 0.1, "Spatial resolution in the y direction [um] (r direction in barrel, r-phi direction in disks)"}; + Gaudi::Property> m_y_resolution{this, "yResolution", {0.1}, "Spatial resolutions in the y direction per layer [um] (r direction in barrel, r-phi direction in disks)"}; // t resolution in ns - FloatProperty m_t_resolution{this, "tResolution", 0.1, "Time resolution [ns]"}; + Gaudi::Property> m_t_resolution{this, "tResolution", {0.1}, "Time resolutions per layer [ns]"}; // Surface manager used to project hits onto sensitive surface with forceHitsOntoSurface argument mutable const dd4hep::rec::SurfaceMap* _map; @@ -94,7 +96,7 @@ class VTXdigitizer : public Gaudi::Algorithm { IRndmGenSvc* m_randSvc; // Gaussian random number generator used for smearing - Rndm::Numbers m_gauss_x; - Rndm::Numbers m_gauss_y; - Rndm::Numbers m_gauss_time; + std::vector m_gauss_x_vec; + std::vector m_gauss_y_vec; + std::vector m_gauss_t_vec; }; diff --git a/VTXdigi/src/VTXdigitizer.cpp b/VTXdigi/src/VTXdigitizer.cpp index 94b04d0..ff6d35a 100644 --- a/VTXdigi/src/VTXdigitizer.cpp +++ b/VTXdigi/src/VTXdigitizer.cpp @@ -17,17 +17,29 @@ StatusCode VTXdigitizer::initialize() { error() << "Couldn't get RndmGenSvc!" << endmsg; return StatusCode::FAILURE; } - if (m_gauss_x.initialize(m_randSvc, Rndm::Gauss(0., m_x_resolution)).isFailure()) { - error() << "Couldn't initialize RndmGenSvc!" << endmsg; - return StatusCode::FAILURE; + + m_gauss_x_vec.resize(m_x_resolution.size()); + for (size_t i = 0; i < m_x_resolution.size(); ++i) { + if (m_gauss_x_vec[i].initialize(m_randSvc, Rndm::Gauss(0., m_x_resolution[i])).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } } - if (m_gauss_y.initialize(m_randSvc, Rndm::Gauss(0., m_y_resolution)).isFailure()) { - error() << "Couldn't initialize RndmGenSvc!" << endmsg; - return StatusCode::FAILURE; + + m_gauss_y_vec.resize(m_y_resolution.size()); + for (size_t i = 0; i < m_y_resolution.size(); ++i) { + if (m_gauss_y_vec[i].initialize(m_randSvc, Rndm::Gauss(0., m_y_resolution[i])).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } } - if (m_gauss_time.initialize(m_randSvc, Rndm::Gauss(0., m_t_resolution)).isFailure()) { - error() << "Couldn't initialize RndmGenSvc!" << endmsg; - return StatusCode::FAILURE; + + m_gauss_t_vec.resize(m_t_resolution.size()); + for (size_t i = 0; i < m_t_resolution.size(); ++i) { + if (m_gauss_t_vec[i].initialize(m_randSvc, Rndm::Gauss(0., m_t_resolution[i])).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } } // check if readout exists @@ -39,6 +51,13 @@ StatusCode VTXdigitizer::initialize() { // set the cellID decoder m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); // Can be used to access e.g. layer index: m_decoder->get(cellID, "layer"), + if (m_decoder->fieldDescription().find("layer") == std::string::npos){ + error() + << " Readout " << m_readoutName << " does not contain layer id!" + << endmsg; + return StatusCode::FAILURE; + } + // retrieve the volume manager m_volman = m_geoSvc->getDetector()->volumeManager(); @@ -126,22 +145,21 @@ StatusCode VTXdigitizer::execute(const EventContext&) const { // Smear the hit in the local sensor coordinates double digiHitLocalPosition[3]; - if (m_readoutName == "VTXIBCollection" || - m_readoutName == "VTXOBCollection" || - m_readoutName == "VertexBarrelCollection" || + int iLayer = m_decoder->get(cellID, "layer"); + debug() << "readout: " << m_readoutName << ", layer id: " << iLayer << endmsg; + if (m_readoutName == "VertexBarrelCollection" || m_readoutName == "SiWrBCollection") { // In barrel, the sensor box is along y-z digiHitLocalPosition[0] = simHitLocalPositionVector.x(); - digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_x.shoot() * dd4hep::mm; - digiHitLocalPosition[2] = simHitLocalPositionVector.z() + m_gauss_y.shoot() * dd4hep::mm; - } else if (m_readoutName == "VTXDCollection" || - m_readoutName == "VertexEndcapCollection" || + digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_x_vec[iLayer].shoot() * dd4hep::mm; + digiHitLocalPosition[2] = simHitLocalPositionVector.z() + m_gauss_y_vec[iLayer].shoot() * dd4hep::mm; + } else if (m_readoutName == "VertexEndcapCollection" || m_readoutName == "SiWrDCollection") { // In the disks, the sensor box is already in x-y - digiHitLocalPosition[0] = simHitLocalPositionVector.x() + m_gauss_x.shoot() * dd4hep::mm; - digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_y.shoot() * dd4hep::mm; + digiHitLocalPosition[0] = simHitLocalPositionVector.x() + m_gauss_x_vec[iLayer].shoot() * dd4hep::mm; + digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_y_vec[iLayer].shoot() * dd4hep::mm; digiHitLocalPosition[2] = simHitLocalPositionVector.z(); } else { error() - << "VTX readout name (m_readoutName) unknown!" + << "VTX readout name (m_readoutName) unknown or xResolution/yResolution/tResolution not defining all detector layer resolutions!" << endmsg; return StatusCode::FAILURE; } @@ -167,7 +185,7 @@ StatusCode VTXdigitizer::execute(const EventContext&) const { output_digi_hit.setPosition(digiHitGlobalPositionVector); // Apply time smearing - output_digi_hit.setTime(input_sim_hit.getTime() + m_gauss_time.shoot()); + output_digi_hit.setTime(input_sim_hit.getTime() + m_gauss_t_vec[iLayer].shoot()); output_digi_hit.setCellID(cellID); diff --git a/VTXdigi/test/runVTXdigitizer.py b/VTXdigi/test/runVTXdigitizer.py index 943542e..bbf65da 100644 --- a/VTXdigi/test/runVTXdigitizer.py +++ b/VTXdigi/test/runVTXdigitizer.py @@ -16,22 +16,22 @@ ################## Vertex sensor resolutions # IDEA -innerVertexResolution_x = 0.003 # [mm], assume 5 µm resolution for ARCADIA sensor -innerVertexResolution_y = 0.003 # [mm], assume 5 µm resolution for ARCADIA sensor +innerVertexResolution_x = 0.003 # [mm], assume 3 µm resolution for ARCADIA sensor +innerVertexResolution_y = 0.003 # [mm], assume 3 µm resolution for ARCADIA sensor innerVertexResolution_t = 1000 # [ns] outerVertexResolution_x = 0.050/math.sqrt(12) # [mm], assume ATLASPix3 sensor with 50 µm pitch outerVertexResolution_y = 0.150/math.sqrt(12) # [mm], assume ATLASPix3 sensor with 150 µm pitch outerVertexResolution_t = 1000 # [ns] # CLD -vertexBarrelResolution_x = 0.003 # [mm], assume 5 µm resolution -vertexBarrelResolution_y = 0.003 # [mm], assume 5 µm resolution +vertexBarrelResolution_x = 0.003 # [mm], assume 3 µm resolution +vertexBarrelResolution_y = 0.003 # [mm], assume 3 µm resolution vertexBarrelResolution_t = 1000 # [ns] -vertexEndcapResolution_x = 0.003 # [mm], assume 5 µm resolution -vertexEndcapResolution_y = 0.003 # [mm], assume 5 µm resolution +vertexEndcapResolution_x = 0.003 # [mm], assume 3 µm resolution +vertexEndcapResolution_y = 0.003 # [mm], assume 3 µm resolution vertexEndcapResolution_t = 1000 # [ns] - +# IDEA silicon wrapper siWrapperResolution_x = 0.050/math.sqrt(12) # [mm] siWrapperResolution_y = 1.0/math.sqrt(12) # [mm] siWrapperResolution_t = 0.040 # [ns], assume 40 ps timing resolution for a single layer -> Should lead to <30 ps resolution when >1 hit @@ -113,21 +113,18 @@ from Configurables import SimG4SaveTrackerHits ### CLD -SimG4SaveTrackerHitsB = SimG4SaveTrackerHits("SimG4SaveTrackerHitsB", readoutName="VertexBarrelCollection") -SimG4SaveTrackerHitsB.SimTrackHits.Path = "VTXB_simTrackerHits" +# SimG4SaveTrackerHitsB = SimG4SaveTrackerHits("SimG4SaveTrackerHitsB", readoutName="VertexBarrelCollection") +# SimG4SaveTrackerHitsB.SimTrackHits.Path = "VTXB_simTrackerHits" -SimG4SaveTrackerHitsE = SimG4SaveTrackerHits("SimG4SaveTrackerHitsE", readoutName="VertexEndcapCollection") -SimG4SaveTrackerHitsE.SimTrackHits.Path = "VTXE_simTrackerHits" +# SimG4SaveTrackerHitsE = SimG4SaveTrackerHits("SimG4SaveTrackerHitsE", readoutName="VertexEndcapCollection") +# SimG4SaveTrackerHitsE.SimTrackHits.Path = "VTXE_simTrackerHits" ### IDEA -SimG4SaveTrackerHitsIB = SimG4SaveTrackerHits("SimG4SaveTrackerHitsIB", readoutName="VTXIBCollection") -SimG4SaveTrackerHitsIB.SimTrackHits.Path = "VTXIB_simTrackerHits" - -SimG4SaveTrackerHitsOB = SimG4SaveTrackerHits("SimG4SaveTrackerHitsOB", readoutName="VTXOBCollection") -SimG4SaveTrackerHitsOB.SimTrackHits.Path = "VTXOB_simTrackerHits" +SimG4SaveTrackerHitsB = SimG4SaveTrackerHits("SimG4SaveTrackerHitsB", readoutName="VertexBarrelCollection") +SimG4SaveTrackerHitsB.SimTrackHits.Path = "VTXB_simTrackerHits" -SimG4SaveTrackerHitsD = SimG4SaveTrackerHits("SimG4SaveTrackerHitsD", readoutName="VTXDCollection") +SimG4SaveTrackerHitsD = SimG4SaveTrackerHits("SimG4SaveTrackerHitsD", readoutName="VertexEndcapCollection") SimG4SaveTrackerHitsD.SimTrackHits.Path = "VTXD_simTrackerHits" SimG4SaveTrackerHitsSiWrB = SimG4SaveTrackerHits("SimG4SaveTrackerHitsSiWrB", readoutName="SiWrBCollection") @@ -148,8 +145,8 @@ # IDEA geantsim = SimG4Alg("SimG4Alg", - outputs= [SimG4SaveTrackerHitsIB, SimG4SaveTrackerHitsOB, SimG4SaveTrackerHitsD, - #SimG4SaveTrackerHitsSiWrB, SimG4SaveTrackerHitsSiWrD, + outputs= [SimG4SaveTrackerHitsB, SimG4SaveTrackerHitsD, + SimG4SaveTrackerHitsSiWrB, SimG4SaveTrackerHitsSiWrD, #saveHistTool ], eventProvider=particle_converter, @@ -159,97 +156,86 @@ from Configurables import VTXdigitizer ### For CLD. Not working yet, SimG4 doesn't produce hits in CLD vertex yet -vtxb_digitizer = VTXdigitizer("VTXBdigitizer", +# cld_vtxb_digitizer = VTXdigitizer("VTXBdigitizer", +# inputSimHits = SimG4SaveTrackerHitsB.SimTrackHits.Path, +# outputDigiHits = SimG4SaveTrackerHitsB.SimTrackHits.Path.replace("sim", "digi"), +# outputSimDigiAssociation = SimG4SaveTrackerHitsB.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), +# detectorName = "Vertex", +# readoutName = "VertexBarrelCollection", +# xResolutions = [vertexBarrelResolution_x, vertexBarrelResolution_x, vertexBarrelResolution_x, vertexBarrelResolution_x, vertexBarrelResolution_x, vertexBarrelResolution_x], +# yResolutions = [vertexBarrelResolution_y, vertexBarrelResolution_y, vertexBarrelResolution_y, vertexBarrelResolution_y, vertexBarrelResolution_y, vertexBarrelResolution_y], +# tResolutions = [vertexBarrelResolution_t, vertexBarrelResolution_t, vertexBarrelResolution_t, vertexBarrelResolution_t, vertexBarrelResolution_t, vertexBarrelResolution_t], +# forceHitsOntoSurface = False, +# OutputLevel = INFO +# ) + +# cld_vtxd_digitizer = VTXdigitizer("VTXDdigitizer", +# inputSimHits = SimG4SaveTrackerHitsE.SimTrackHits.Path, +# outputDigiHits = SimG4SaveTrackerHitsE.SimTrackHits.Path.replace("sim", "digi"), +# outputSimDigiAssociation = SimG4SaveTrackerHitsE.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), +# detectorName = "Vertex", +# readoutName = "VertexEndcapCollection", +# xResolutions = [vertexEndcapResolution_x, vertexEndcapResolution_x, vertexEndcapResolution_x, vertexEndcapResolution_x, vertexEndcapResolution_x, vertexEndcapResolution_x], +# yResolutions = [vertexEndcapResolution_y, vertexEndcapResolution_y, vertexEndcapResolution_y, vertexEndcapResolution_y, vertexEndcapResolution_y, vertexEndcapResolution_y], +# tResolutions = [vertexEndcapResolution_t, vertexEndcapResolution_t, vertexEndcapResolution_t, vertexEndcapResolution_t, vertexEndcapResolution_t, vertexEndcapResolution_t], +# forceHitsOntoSurface = False, +# OutputLevel = INFO +# ) + + +### For IDEA +idea_vtxb_digitizer = VTXdigitizer("VTXBdigitizer", inputSimHits = SimG4SaveTrackerHitsB.SimTrackHits.Path, outputDigiHits = SimG4SaveTrackerHitsB.SimTrackHits.Path.replace("sim", "digi"), outputSimDigiAssociation = SimG4SaveTrackerHitsB.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), detectorName = "Vertex", readoutName = "VertexBarrelCollection", - xResolution = vertexBarrelResolution_x, - yResolution = vertexBarrelResolution_y, - tResolution = vertexBarrelResolution_t, + xResolution = [innerVertexResolution_x, innerVertexResolution_x, innerVertexResolution_x, outerVertexResolution_x, outerVertexResolution_x], # mm, r-phi direction + yResolution = [innerVertexResolution_y, innerVertexResolution_y, innerVertexResolution_y, outerVertexResolution_y, outerVertexResolution_y], # mm, z direction + tResolution = [innerVertexResolution_t, innerVertexResolution_t, innerVertexResolution_t, outerVertexResolution_t, outerVertexResolution_t], # ns forceHitsOntoSurface = False, OutputLevel = INFO ) -vtxe_digitizer = VTXdigitizer("VTXEdigitizer", - inputSimHits = SimG4SaveTrackerHitsE.SimTrackHits.Path, - outputDigiHits = SimG4SaveTrackerHitsE.SimTrackHits.Path.replace("sim", "digi"), - outputSimDigiAssociation = SimG4SaveTrackerHitsE.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), +idea_vtxd_digitizer = VTXdigitizer("VTXDdigitizer", + inputSimHits = SimG4SaveTrackerHitsD.SimTrackHits.Path, + outputDigiHits = SimG4SaveTrackerHitsD.SimTrackHits.Path.replace("sim", "digi"), + outputSimDigiAssociation = SimG4SaveTrackerHitsD.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), detectorName = "Vertex", readoutName = "VertexEndcapCollection", - xResolution = vertexEndcapResolution_x, - yResolution = vertexEndcapResolution_y, - tResolution = vertexEndcapResolution_t, + xResolution = [outerVertexResolution_x, outerVertexResolution_x, outerVertexResolution_x], # mm, r direction + yResolution = [outerVertexResolution_y, outerVertexResolution_y, outerVertexResolution_y], # mm, phi direction + tResolution = [outerVertexResolution_t, outerVertexResolution_t, outerVertexResolution_t], # ns forceHitsOntoSurface = False, OutputLevel = INFO ) - -### For IDEA -vtxib_digitizer = VTXdigitizer("VTXIBdigitizer", - inputSimHits = SimG4SaveTrackerHitsIB.SimTrackHits.Path, - outputDigiHits = SimG4SaveTrackerHitsIB.SimTrackHits.Path.replace("sim", "digi"), - outputSimDigiAssociation = SimG4SaveTrackerHitsIB.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), - detectorName = "Vertex", - readoutName = "VTXIBCollection", - xResolution = innerVertexResolution_x, # mm, r-phi direction - yResolution = innerVertexResolution_y, # mm, z direction - tResolution = innerVertexResolution_t, - forceHitsOntoSurface = False, - OutputLevel = INFO -) - -vtxob_digitizer = VTXdigitizer("VTXOBdigitizer", - inputSimHits = SimG4SaveTrackerHitsOB.SimTrackHits.Path, - outputDigiHits = SimG4SaveTrackerHitsOB.SimTrackHits.Path.replace("sim", "digi"), - outputSimDigiAssociation = SimG4SaveTrackerHitsOB.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), - detectorName = "Vertex", - readoutName = "VTXOBCollection", - xResolution = outerVertexResolution_x, # mm, r-phi direction - yResolution = outerVertexResolution_y, # mm, z direction - tResolution = outerVertexResolution_t, # ns - forceHitsOntoSurface = False, - OutputLevel = INFO +idea_siwrb_digitizer = VTXdigitizer("SiWrBdigitizer", + inputSimHits = SimG4SaveTrackerHitsSiWrB.SimTrackHits.Path, + outputDigiHits = SimG4SaveTrackerHitsSiWrB.SimTrackHits.Path.replace("sim", "digi"), + outputSimDigiAssociation = SimG4SaveTrackerHitsSiWrB.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), + detectorName = "SiWrB", + readoutName = "SiWrBCollection", + xResolution = [siWrapperResolution_x, siWrapperResolution_x], # mm, r-phi direction + yResolution = [siWrapperResolution_y, siWrapperResolution_y], # mm, z direction + tResolution = [siWrapperResolution_t, siWrapperResolution_t], # ns + forceHitsOntoSurface = False, + OutputLevel = INFO ) -vtxd_digitizer = VTXdigitizer("VTXDdigitizer", - inputSimHits = SimG4SaveTrackerHitsD.SimTrackHits.Path, - outputDigiHits = SimG4SaveTrackerHitsD.SimTrackHits.Path.replace("sim", "digi"), - outputSimDigiAssociation = SimG4SaveTrackerHitsD.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), - detectorName = "Vertex", - readoutName = "VTXDCollection", - xResolution = outerVertexResolution_x, # mm, r direction - yResolution = outerVertexResolution_y, # mm, phi direction - tResolution = outerVertexResolution_t, # ns - forceHitsOntoSurface = False, - OutputLevel = INFO +idea_siwrd_digitizer = VTXdigitizer("SiWrDdigitizer", + inputSimHits = SimG4SaveTrackerHitsSiWrD.SimTrackHits.Path, + outputDigiHits = SimG4SaveTrackerHitsSiWrD.SimTrackHits.Path.replace("sim", "digi"), + outputSimDigiAssociation = SimG4SaveTrackerHitsSiWrD.SimTrackHits.Path.replace("simTrackerHits", "simDigiAssociation"), + detectorName = "SiWrD", + readoutName = "SiWrDCollection", + xResolution = [siWrapperResolution_x, siWrapperResolution_x], # mm, r direction + yResolution = [siWrapperResolution_y, siWrapperResolution_y], # mm, phi direction + tResolution = [siWrapperResolution_t, siWrapperResolution_t], # ns + forceHitsOntoSurface = False, + OutputLevel = INFO ) -#siwrb_digitizer = VTXdigitizer("SiWrBdigitizer", -# inputSimHits = SimG4SaveTrackerHitsSiWrB.SimTrackHits.Path, -# outputDigiHits = SimG4SaveTrackerHitsSiWrB.SimTrackHits.Path.replace("sim", "digi"), -# detectorName = "SiliconWrapper", -# readoutName = "SiWrBCollection", -# xResolution = siWrapperResolution_x, # mm, r direction -# yResolution = siWrapperResolution_y, # mm, phi direction -# tResolution = siWrapperResolution_t, # ns -# forceHitsOntoSurface = False, -# OutputLevel = INFO -#) -# -#siwrd_digitizer = VTXdigitizer("SiWrDdigitizer", -# inputSimHits = SimG4SaveTrackerHitsSiWrD.SimTrackHits.Path, -# outputDigiHits = SimG4SaveTrackerHitsSiWrD.SimTrackHits.Path.replace("sim", "digi"), -# detectorName = "SiliconWrapper", -# readoutName = "SiWrDCollection", -# xResolution = siWrapperResolution_x, # mm, r direction -# yResolution = siWrapperResolution_y, # mm, phi direction -# tResolution = siWrapperResolution_t, # ns -# forceHitsOntoSurface = False, -# OutputLevel = INFO -#) - # run the genfit tracking # from Configurables import GenFitter # genfitter = GenFitter("GenFitter", inputHits = savetrackertool.SimTrackHits.Path.replace("sim", "digi"), outputTracks = "genfit_tracks") @@ -281,7 +267,7 @@ # genAlg, # hepmc_converter, # geantsim, -# vtxb_digitizer,vtxe_digitizer, +# cld_vtxb_digitizer, cld_vtxd_digitizer, # out # ], # EvtSel = 'NONE', @@ -296,8 +282,8 @@ genAlg, hepmc_converter, geantsim, - vtxib_digitizer, vtxob_digitizer, vtxd_digitizer, - #siwrb_digitizer, siwrd_digitizer, + idea_vtxb_digitizer, idea_vtxd_digitizer, + idea_siwrb_digitizer, idea_siwrd_digitizer, out ], EvtSel = 'NONE',