First version of VTX digitisation for IDEA vertex detector

Makefile

@@ -3,7 +3,6 @@ MAKEFLAGS += --no-print-directory
 make:
 	@ mkdir -p build install ; \
 	cd build ; \
-	source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh ; \
 	cmake .. -DCMAKE_INSTALL_PREFIX=../install ; \
 	make install -j4 ; \
 	cd .. ; \

Tracking/test/runGenFitTrackingOnSimplifiedDriftChamber.py

@@ -91,32 +91,42 @@
 saveDCHsimHitTool = SimG4SaveTrackerHits("saveDCHsimHitTool", readoutNames=["SimplifiedDriftChamberCollection"])
 saveDCHsimHitTool.SimTrackHits.Path = "DC_simTrackerHits"
 
-saveVTXBsimHitTool = SimG4SaveTrackerHits("saveVTXBsimHitTool", readoutNames=["VertexBarrelCollection"])
-saveVTXBsimHitTool.SimTrackHits.Path = "VTXB_simTrackerHits"
+saveVTXIBsimHitTool = SimG4SaveTrackerHits("saveVTXIBsimHitTool", readoutNames=["VTXIBCollection"])
+saveVTXIBsimHitTool.SimTrackHits.Path = "VTXIB_simTrackerHits"
 
-saveVTXEsimHitTool = SimG4SaveTrackerHits("saveVTXEsimHitTool", readoutNames=["VertexEndcapCollection"])
-saveVTXEsimHitTool.SimTrackHits.Path = "VTXE_simTrackerHits"
+saveVTXOBsimHitTool = SimG4SaveTrackerHits("saveVTXOBsimHitTool", readoutNames=["VTXOBCollection"])
+saveVTXOBsimHitTool.SimTrackHits.Path = "VTXOB_simTrackerHits"
+
+saveVTXDsimHitTool = SimG4SaveTrackerHits("saveVTXDsimHitTool", readoutNames=["VTXDCollection"])
+saveVTXDsimHitTool.SimTrackHits.Path = "VTXD_simTrackerHits"
 
 from Configurables import SimG4Alg
 geantsim = SimG4Alg("SimG4Alg",
-                       outputs= [saveVTXBsimHitTool, saveVTXEsimHitTool, saveDCHsimHitTool,
+                       outputs= [saveVTXIBsimHitTool, saveVTXOBsimHitTool, saveVTXDsimHitTool, saveDCHsimHitTool,
                                  #saveHistTool
                        ],
                        eventProvider = particle_converter,
                        OutputLevel = INFO)
 
 # Digitize tracker hits (for now digitization is reconstruction)
-vtxb_reco_hit_name = saveVTXBsimHitTool.SimTrackHits.Path.replace("sim", "reco")
+vtxib_reco_hit_name = saveVTXIBsimHitTool.SimTrackHits.Path.replace("sim", "reco")
+from Configurables import VTXdigitizer
+vtxib_digitizer = VTXdigitizer("VTXIBdigitizer",
+    inputSimHits = saveVTXIBsimHitTool.SimTrackHits.Path,
+    outputDigiHits = vtxib_reco_hit_name
+)
+
+vtxob_reco_hit_name = saveVTXOBsimHitTool.SimTrackHits.Path.replace("sim", "reco")
 from Configurables import VTXdigitizer
-vtxb_digitizer = VTXdigitizer("VTXBdigitizer",
-    inputSimHits = saveVTXBsimHitTool.SimTrackHits.Path,
-    outputDigiHits = vtxb_reco_hit_name
+vtxob_digitizer = VTXdigitizer("VTXOBdigitizer",
+    inputSimHits = saveVTXOBsimHitTool.SimTrackHits.Path,
+    outputDigiHits = vtxob_reco_hit_name
 )
 
-vtxe_reco_hit_name = saveVTXEsimHitTool.SimTrackHits.Path.replace("sim", "reco")
-vtxe_digitizer = VTXdigitizer("VTXEdigitizer",
-    inputSimHits = saveVTXEsimHitTool.SimTrackHits.Path,
-    outputDigiHits = vtxe_reco_hit_name
+vtxd_reco_hit_name = saveVTXDsimHitTool.SimTrackHits.Path.replace("sim", "reco")
+vtxd_digitizer = VTXdigitizer("VTXDdigitizer",
+    inputSimHits = saveVTXDsimHitTool.SimTrackHits.Path,
+    outputDigiHits = vtxd_reco_hit_name
 )
 
 dch_reco_hit_name = saveDCHsimHitTool.SimTrackHits.Path.replace("sim", "reco")
@@ -161,7 +171,7 @@
               hepmc_converter,
               geantsim,
               vtxb_digitizer,
-              vtxe_digitizer,
+              vtxd_digitizer,
               dch_digitizer,
               genfitter,
               out

VTXdigi/CMakeLists.txt

@@ -16,6 +16,7 @@ gaudi_add_module(${PackageName}
   Gaudi::GaudiKernel
   EDM4HEP::edm4hep
   k4FWCore::k4FWCore
+  DD4hep::DDRec
 )
 
 target_include_directories(${PackageName} PUBLIC
@@ -39,3 +40,6 @@ install(TARGETS ${PackageName}
 )
 
 install(FILES ${scripts} DESTINATION test)
+
+SET(test_name "test_runVTXdigitizer")
+ADD_TEST(NAME t_${test_name} COMMAND k4run test/runVTXdigitizer.py)

VTXdigi/include/VTXdigitizer.h

@@ -3,14 +3,24 @@
 // 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 "DDRec/Vector3D.h"
+#include "DDRec/SurfaceManager.h"
+
+#include "DDSegmentation/BitFieldCoder.h"
+
 /** @class VTXdigitizer
  *
  *  Algorithm for creating digitized (meaning 'reconstructed' for now) vertex detector hits (edm4hep::TrackerHit) from Geant4 hits (edm4hep::SimTrackerHit).
@@ -42,4 +52,38 @@ class VTXdigitizer : public GaudiAlgorithm {
   DataHandle<edm4hep::SimTrackerHitCollection> m_input_sim_hits{"inputSimHits", Gaudi::DataHandle::Reader, this};
   // Output digitized vertex hit collection name
   DataHandle<edm4hep::TrackerHitCollection> m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this};
+
+  // Detector name
+  Gaudi::Property<std::string> m_detectorName{this, "detectorName", "Vertex", "Name of the detector (default: Vertex)"};
+  // Detector readout names
+  Gaudi::Property<std::string> m_readoutName{this, "readoutName", "VertexBarrelCollection", "Name of the detector readout"};
+  // Pointer to the geometry service
+  ServiceHandle<IGeoSvc> m_geoSvc;
+  // Decoder for the cellID
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+  // Volume manager to get the physical cell sensitive volume
+  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)"};
+
+  // 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)"};
+
+  // t resolution in ns
+  FloatProperty m_t_resolution{this, "tResolution", 0.1, "Time resolution [ns]"};
+
+  // Surface manager used to project hits onto sensitive surface with forceHitsOntoSurface argument
+  const dd4hep::rec::SurfaceMap* _map;
+
+  // Option to force hits onto sensitive surface
+  BooleanProperty m_forceHitsOntoSurface{this, "forceHitsOntoSurface", false, "Project hits onto the surface in case they are not yet on the surface (default: false"};
+
+  // Random Number Service
+  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;
 };

VTXdigi/src/VTXdigitizer.cpp

@@ -2,26 +2,175 @@
 
 DECLARE_COMPONENT(VTXdigitizer)
 
-VTXdigitizer::VTXdigitizer(const std::string& aName, ISvcLocator* aSvcLoc) : GaudiAlgorithm(aName, aSvcLoc) {
+VTXdigitizer::VTXdigitizer(const std::string& aName, ISvcLocator* aSvcLoc)
+    : GaudiAlgorithm(aName, aSvcLoc), m_geoSvc("GeoSvc", "VTXdigitizer") {
   declareProperty("inputSimHits", m_input_sim_hits, "Input sim vertex hit collection name");
   declareProperty("outputDigiHits", m_output_digi_hits, "Output digitized vertex hit collection name");
 }
 
 VTXdigitizer::~VTXdigitizer() {}
 
-StatusCode VTXdigitizer::initialize() { return StatusCode::SUCCESS; }
+StatusCode VTXdigitizer::initialize() {
+  // Initialize random services
+  if (service("RndmGenSvc", m_randSvc).isFailure()) {
+    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;
+  }
+  if (m_gauss_y.initialize(m_randSvc, Rndm::Gauss(0., m_y_resolution)).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;
+  }
+
+  // check if readout exists
+  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->getDetector()->readout(m_readoutName).idSpec().decoder(); // Can be used to access e.g. layer index: m_decoder->get(cellID, "layer"),
+
+  // retrieve the volume manager
+  m_volman = m_geoSvc->getDetector()->volumeManager();
+
+  return StatusCode::SUCCESS;
+}
 
 StatusCode VTXdigitizer::execute() {
   // Get the input collection with Geant4 hits
   const edm4hep::SimTrackerHitCollection* input_sim_hits = m_input_sim_hits.get();
   verbose() << "Input Sim Hit collection size: " << input_sim_hits->size() << endmsg;
 
+  unsigned nDismissedHits=0;
+
   // 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 local frame of the silicon sensor to smear in the direction along/perpendicular to the stave
+
+    // retrieve the cell detElement
+    dd4hep::DDSegmentation::CellID cellID = input_sim_hit.getCellID();
+    debug() << "Digitisation of " << m_readoutName << ", cellID: " << cellID << endmsg;
+
+    // Get transformation matrix of sensor
+    const auto& sensorTransformMatrix = m_volman.lookupVolumePlacement(cellID).matrix();
+
+    // 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};
+    dd4hep::rec::Vector3D simHitGlobalPositionVector(simHitGlobalPosition[0], simHitGlobalPosition[1],
+                                                    simHitGlobalPosition[2]);
+    double simHitLocalPosition[3]  = {0, 0, 0};
+
+    if(m_forceHitsOntoSurface){
+      dd4hep::Detector& theDetector = dd4hep::Detector::getInstance();
+      dd4hep::rec::SurfaceManager& surfMan = *theDetector.extension<dd4hep::rec::SurfaceManager>() ;
+      dd4hep::DetElement det = m_geoSvc->getDetector()->detector(m_detectorName) ;
+      _map = surfMan.map( det.name() ) ;
+
+      if( ! _map )
+        error() << " Could not find surface map for detector " << det.name() << " in SurfaceManager " << endmsg;
+
+      dd4hep::rec::SurfaceMap::const_iterator sI = _map->find(cellID) ;
+      if( sI == _map->end() )
+        error() << " VTXdigitizer: no surface found for cellID " << m_decoder->valueString(cellID) << std::endl << endmsg;
+
+      dd4hep::rec::Vector3D newPos ;
+      const dd4hep::rec::ISurface* surf = sI->second ;    
+
+      // Check if Hit is inside sensitive 
+      if ( ! surf->insideBounds( simHitGlobalPositionVector ) ) {
+
+        info() << "Hit at " << simHitGlobalPositionVector << " is not on surface " << *surf  
+                << ". Distance: " << surf->distance(simHitGlobalPositionVector )
+                << std::endl << endmsg;        
+
+        if( m_forceHitsOntoSurface ){
+          dd4hep::rec::Vector2D lv = surf->globalToLocal(simHitGlobalPositionVector) ;
+          dd4hep::rec::Vector3D oldPosOnSurf = surf->localToGlobal( lv ) ; 
+
+          info() << "Moved hit to " << oldPosOnSurf << ", distance " << (oldPosOnSurf-simHitGlobalPositionVector).r()
+                  << std::endl << endmsg;
+
+          simHitGlobalPositionVector = oldPosOnSurf ;
+        } 
+        else
+          ++nDismissedHits;
+      }
+    }
+
+    // get the simHit coordinate in cm in the sensor reference frame to be able to apply smearing
+    sensorTransformMatrix.MasterToLocal(simHitGlobalPosition, simHitLocalPosition);
+    debug() << "Cell ID string: " << m_decoder->valueString(cellID) << endmsg;
+    ;
+    debug() << "Global simHit x " << simHitGlobalPosition[0] << " [mm] --> Local simHit x " << simHitLocalPosition[0]
+            << " [cm]" << endmsg;
+    debug() << "Global simHit y " << simHitGlobalPosition[1] << " [mm] --> Local simHit y " << simHitLocalPosition[1]
+            << " [cm]" << endmsg;
+    debug() << "Global simHit z " << simHitGlobalPosition[2] << " [mm] --> Local simHit z " << simHitLocalPosition[2]
+            << " [cm]" << endmsg;
+
+    // build a vector to easily apply smearing of distance to the wire
+    dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0], simHitLocalPosition[1],
+                                                    simHitLocalPosition[2]);
+
+    // Smear the hit in the local sensor coordinates
+    double digiHitLocalPosition[3];
+    if (m_readoutName == "VTXIBCollection" ||
+        m_readoutName == "VTXOBCollection" ||
+        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" ||
+               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[2] = simHitLocalPositionVector.z();
+    } else {
+      error()
+          << "VTX readout name (m_readoutName) unknown!"
+          << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    // go back to the global frame
+    double digiHitGlobalPosition[3] = {0, 0, 0};
+    sensorTransformMatrix.LocalToMaster(digiHitLocalPosition, digiHitGlobalPosition);
+
+    // go back to mm
+    edm4hep::Vector3d digiHitGlobalPositionVector(digiHitGlobalPosition[0] / dd4hep::mm,
+                                                  digiHitGlobalPosition[1] / dd4hep::mm,
+                                                  digiHitGlobalPosition[2] / dd4hep::mm);
+
+    debug() << "Global digiHit x " << digiHitGlobalPositionVector[0] << " [mm] --> Local digiHit x "
+            << digiHitLocalPosition[0] << " [cm]" << endmsg;
+    debug() << "Global digiHit y " << digiHitGlobalPositionVector[1] << " [mm] --> Local digiHit y "
+            << digiHitLocalPosition[1] << " [cm]" << endmsg;
+    debug() << "Global digiHit z " << digiHitGlobalPositionVector[2] << " [mm] --> Local digiHit z "
+            << digiHitLocalPosition[2] << " [cm]" << endmsg;
+    debug() << "Moving to next hit... " << std::endl << endmsg;
+
     output_digi_hit.setEDep(input_sim_hit.getEDep());
-    output_digi_hit.setPosition(input_sim_hit.getPosition());
+    output_digi_hit.setPosition(digiHitGlobalPositionVector);
+
+    // Apply time smearing
+    output_digi_hit.setTime(input_sim_hit.getTime() + m_gauss_time.shoot());
+
+    output_digi_hit.setCellID(cellID);
   }
   return StatusCode::SUCCESS;
 }