Drift chamber digitized hit in global coordinate (#19)

* DCHdigi in global coordinates

* Make the DCHdigi in global coordinates the default

* [DCHdigi] write down assumptions

* Add an association between DCHdigi and simHits

* Add debug collections and optional local coordinate collection

* Fixes

---------

Co-authored-by: jmcarcell <[email protected]>
Co-authored-by: Juan Miguel Carceller <[email protected]>

DCHdigi/dataFormatExtension/driftChamberHit.yaml

@@ -9,8 +9,8 @@ options:
 
 datatypes:
 
-  extension::DriftChamberDigi:
-    Description: "Drift chamber digitized hit (before tracking)"
+  extension::DriftChamberDigiLocal:
+    Description: "Drift chamber digitized hit (before tracking) in local coordinates"
     Author: "B. Francois, CERN"
     Members:
       - uint64_t cellID               // ID of the wire that created this hit
@@ -20,3 +20,24 @@ datatypes:
       - float eDep                    // energy deposited on the hit [GeV].
       - float eDepError               // error measured on eDep [GeV].
       - uint32_t clusterCount         // number of clusters associated to this hit
+
+  extension::DriftChamberDigi:
+    Description: "Drift chamber digitized hit (before tracking) in global coordinates. Assumes that the hits are radially in the middle of the cells"
+    Author: "B. Francois, CERN"
+    Members:
+      - uint64_t cellID               // ID of the wire that created this hit
+      - edm4hep::Vector3d leftPosition // position of the hit assuming it was on the left side of the wire, radially in the middle of the cell [mm]
+      - edm4hep::Vector3d rightPosition // position of the hit assuming it was on the right side of the wire, radially in the middle of the cell [mm]
+      - float time                    // time of the hit [ns].
+      - float eDep                    // energy deposited on the hit [GeV].
+      - float eDepError               // error measured on eDep [GeV].
+      - uint32_t clusterCount         // number of clusters associated to this hit
+
+  extension::MCRecoDriftChamberDigiAssociation:
+    Description: "Association between a DriftChamberDigi and the corresponding simulated hit"
+    Author: "B. Francois, CERN"
+    Members:
+      - float weight              // weight of this association
+    OneToOneRelations:
+      - extension::DriftChamberDigi digi    // reference to the digitized hit
+      - edm4hep::SimTrackerHit sim // reference to the simulated hit

DCHdigi/include/DCHsimpleDigitizerExtendedEdm.h

@@ -10,11 +10,14 @@
 #include "k4FWCore/DataHandle.h"
 #include "k4Interface/IGeoSvc.h"
 
-// EDM4HEP
+// EDM4HEP & PODIO
 #include "edm4hep/SimTrackerHitCollection.h"
+#include "podio/UserDataCollection.h"
 
 // EDM4HEP extension
 #include "extension/DriftChamberDigiCollection.h"
+#include "extension/DriftChamberDigiLocalCollection.h"
+#include "extension/MCRecoDriftChamberDigiAssociationCollection.h"
 
 // DD4HEP
 #include "DD4hep/Detector.h"  // for dd4hep::VolumeManager
@@ -52,6 +55,10 @@ class DCHsimpleDigitizerExtendedEdm : public GaudiAlgorithm {
   DataHandle<edm4hep::SimTrackerHitCollection> m_input_sim_hits{"inputSimHits", Gaudi::DataHandle::Reader, this};
   // Output digitized tracker hit collection name
   DataHandle<extension::DriftChamberDigiCollection> m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this};
+  // Output association between digitized and simulated hit collections
+  DataHandle<extension::MCRecoDriftChamberDigiAssociationCollection> m_output_sim_digi_association{"outputSimDigiAssociation", Gaudi::DataHandle::Writer, this};
+  // Output digitized tracker hit in local coordinates collection name. Only filled in debug mode
+  DataHandle<extension::DriftChamberDigiLocalCollection> m_output_digi_local_hits{"outputDigiLocalHits", Gaudi::DataHandle::Writer, this};
 
   // Detector readout name
   Gaudi::Property<std::string> m_readoutName{this, "readoutName", "CDCHHits", "Name of the detector readout"};
@@ -67,6 +74,13 @@ class DCHsimpleDigitizerExtendedEdm : public GaudiAlgorithm {
                                "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]"};
+  // Flag to produce debugging distributions
+  Gaudi::Property<bool> m_debugMode{this, "debugMode", false, "Flag to produce debugging distributions"};
+  // Declaration of debugging distributions
+  DataHandle<podio::UserDataCollection<double>> m_leftHitSimHitDeltaDistToWire{"leftHitSimHitDeltaDistToWire", Gaudi::DataHandle::Writer, this}; // mm
+  DataHandle<podio::UserDataCollection<double>> m_leftHitSimHitDeltaLocalZ{"leftHitSimHitDeltaLocalZ", Gaudi::DataHandle::Writer, this}; // mm
+  DataHandle<podio::UserDataCollection<double>> m_rightHitSimHitDeltaDistToWire{"rightHitSimHitDeltaDistToWire", Gaudi::DataHandle::Writer, this}; // mm
+  DataHandle<podio::UserDataCollection<double>> m_rightHitSimHitDeltaLocalZ{"rightHitSimHitDeltaLocalZ", Gaudi::DataHandle::Writer, this}; // mm
 
   // Random Number Service
   IRndmGenSvc* m_randSvc;

DCHdigi/src/DCHsimpleDigitizerExtendedEdm.cpp

@@ -13,6 +13,7 @@ DCHsimpleDigitizerExtendedEdm::DCHsimpleDigitizerExtendedEdm(const std::string&
     : GaudiAlgorithm(aName, aSvcLoc), m_geoSvc("GeoSvc", "DCHsimpleDigitizerExtendedEdm") {
   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");
 }
 
 DCHsimpleDigitizerExtendedEdm::~DCHsimpleDigitizerExtendedEdm() {}
@@ -50,8 +51,20 @@ StatusCode DCHsimpleDigitizerExtendedEdm::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
+  // Prepare a collection for digitized hits in local coordinate (only filled in debug mode)
   extension::DriftChamberDigiCollection* output_digi_hits = m_output_digi_hits.createAndPut();
+  extension::DriftChamberDigiLocalCollection* output_digi_local_hits = m_output_digi_local_hits.createAndPut();
+
+  // Prepare a collection for the association between digitized and simulated hit, setting weights to 1
+  extension::MCRecoDriftChamberDigiAssociationCollection* digi_sim_associations = m_output_sim_digi_association.createAndPut();
+
+  // Prepare collections for the debugging distributions
+  auto leftHitSimHitDeltaDistToWire = m_leftHitSimHitDeltaDistToWire.createAndPut();
+  auto leftHitSimHitDeltaLocalZ = m_leftHitSimHitDeltaLocalZ.createAndPut();
+  auto rightHitSimHitDeltaDistToWire = m_rightHitSimHitDeltaDistToWire.createAndPut();
+  auto rightHitSimHitDeltaLocalZ = m_rightHitSimHitDeltaLocalZ.createAndPut();
+
+  // Digitize the sim hits
   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
@@ -63,9 +76,9 @@ StatusCode DCHsimpleDigitizerExtendedEdm::execute() {
         Form("superLayer_%ld_layer_%ld_phi_%ld_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
+    // get the transformation matrix used to place the wire (DD4hep works with cm)
     const auto& wireTransformMatrix = wireDetElement.nominal().worldTransformation();
-    // Retrieve global position in mm and apply unit transformation (translation matrix is stored in cm)
+    // Retrieve global position in mm and transform it to cm because the DD4hep 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};
@@ -79,23 +92,71 @@ StatusCode DCHsimpleDigitizerExtendedEdm::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
-    dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0], simHitLocalPosition[1],
-                                                    simHitLocalPosition[2]);
+    // build a vector to easily apply smearing of distance to the wire, going back to mm
+    dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0] / dd4hep::mm, simHitLocalPosition[1] / dd4hep::mm,
+                                                    simHitLocalPosition[2] / dd4hep::mm);
     // get the smeared distance to the wire (cylindrical coordinate as the smearing should be perpendicular to the wire)
-    debug() << "Original distance to wire: " << simHitLocalPositionVector.rho() << endmsg;
-    double smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm;
+    debug() << "Original distance to wire: " << simHitLocalPositionVector.rho() << " mm" << endmsg;
+    double smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot();
     while(smearedDistanceToWire < 0){
       debug() << "Negative smearedDistanceToWire (" << smearedDistanceToWire << ") shooting another random number" << endmsg;
-      smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm;
+      smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot();
     }
-    debug() << "Smeared distance to wire: " << smearedDistanceToWire << endmsg;
+    debug() << "Smeared distance to wire: " << smearedDistanceToWire << " mm " << endmsg;
     // 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;
+    double smearedZ = simHitLocalPositionVector.z() + m_gauss_z.shoot();
+    // NB: here we assume the hit is radially in the middle of the cell
+    // create the left and right hit local position (in cm again because of the transform matrix)
+    double leftHitLocalPosition[3]  = {-1 * smearedDistanceToWire  * dd4hep::mm, 0, smearedZ * dd4hep::mm};
+    double rightHitLocalPosition[3]  = {smearedDistanceToWire * dd4hep::mm, 0, smearedZ * dd4hep::mm};
+    // transform the left and right hit local position in global coordinate (still cm here)
+    double leftHitGlobalPosition[3]  = {0, 0, 0};
+    double rightHitGlobalPosition[3]  = {0, 0, 0};
+    wireTransformMatrix.LocalToMaster(leftHitLocalPosition, leftHitGlobalPosition);
+    wireTransformMatrix.LocalToMaster(rightHitLocalPosition, rightHitGlobalPosition);
+    //std::cout << (leftHitGlobalPosition[2]==rightHitGlobalPosition[2]) <<  std::endl; // FIXME why are left and right global z coordinates the same?
+    debug() << "Global leftHit x " << leftHitGlobalPosition[0] << " --> Local leftHit x " << leftHitLocalPosition[0]
+            << " in cm" << endmsg;
+    debug() << "Global leftHit y " << leftHitGlobalPosition[1] << " --> Local leftHit y " << leftHitLocalPosition[1]
+            << " in cm" << endmsg;
+    debug() << "Global leftHit z " << leftHitGlobalPosition[2] << " --> Local leftHit z " << leftHitLocalPosition[2]
+            << " in cm" << endmsg;
+    debug() << "Global rightHit x " << rightHitGlobalPosition[0] << " --> Local rightHit x " << rightHitLocalPosition[0]
+            << " in cm" << endmsg;
+    debug() << "Global rightHit y " << rightHitGlobalPosition[1] << " --> Local rightHit y " << rightHitLocalPosition[1]
+            << " in cm" << endmsg;
+    debug() << "Global rightHit z " << rightHitGlobalPosition[2] << " --> Local rightHit z " << rightHitLocalPosition[2]
+            << " in cm" << endmsg;
     // fill the output DriftChamberDigi (making sure we are back in mm)
     output_digi_hit.setCellID(cellID);
-    output_digi_hit.setDistanceToWire(smearedDistanceToWire / dd4hep::mm);
-    output_digi_hit.setZPositionAlongWire(smearedZ / dd4hep::mm);
+    edm4hep::Vector3d leftHitGlobalPositionVector = edm4hep::Vector3d(leftHitGlobalPosition[0] / dd4hep::mm, leftHitGlobalPosition[1] / dd4hep::mm, leftHitGlobalPosition[2] / dd4hep::mm);
+    edm4hep::Vector3d rightHitGlobalPositionVector = edm4hep::Vector3d(rightHitGlobalPosition[0] / dd4hep::mm, rightHitGlobalPosition[1] / dd4hep::mm, rightHitGlobalPosition[2] / dd4hep::mm);
+    output_digi_hit.setLeftPosition(leftHitGlobalPositionVector);
+    output_digi_hit.setRightPosition(rightHitGlobalPositionVector);
+    output_digi_hit.setTime(input_sim_hit.getTime()); // will apply smearing when we know more from R&D teams
+    //output_digi_hit.setEDep(input_sim_hit.getEDep()); // will enable this when we know more from R&D teams
+
+
+    // 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);
+
+    // if required, populate debugging distributions
+    if(m_debugMode) {
+      // Fill the local coordinate digi hit
+      auto output_digi_local_hit = output_digi_local_hits->create();
+      output_digi_local_hit.setCellID(cellID);
+      output_digi_local_hit.setDistanceToWire(smearedDistanceToWire);
+      output_digi_local_hit.setZPositionAlongWire(smearedZ);
+      // produce a vector instead of using directly smearedDistanceToWire or smearedZ for a more thorough testing
+      dd4hep::rec::Vector3D leftHitLocalPositionVector(leftHitLocalPosition[0] / dd4hep::mm, leftHitLocalPosition[1] / dd4hep::mm, leftHitLocalPosition[2] / dd4hep::mm);
+      dd4hep::rec::Vector3D rightHitLocalPositionVector(rightHitLocalPosition[0] / dd4hep::mm, rightHitLocalPosition[1] / dd4hep::mm, rightHitLocalPosition[2] / dd4hep::mm);
+      leftHitSimHitDeltaDistToWire->push_back(leftHitLocalPositionVector.rho() - simHitLocalPositionVector.rho());
+      leftHitSimHitDeltaLocalZ->push_back(leftHitLocalPositionVector.z() - simHitLocalPositionVector.z());
+      rightHitSimHitDeltaDistToWire->push_back(rightHitLocalPositionVector.rho() - simHitLocalPositionVector.rho());
+      rightHitSimHitDeltaLocalZ->push_back(rightHitLocalPositionVector.z() - simHitLocalPositionVector.z());
+    }
   }
   debug() << "Output Digi Hit collection size: " << output_digi_hits->size() << endmsg;
   return StatusCode::SUCCESS;

DCHdigi/test/runDCHsimpleDigitizerExtendedEdm.py

@@ -7,7 +7,7 @@
 
 from GaudiKernel.SystemOfUnits import MeV, GeV, tesla
 ################## Particle gun setup
-momentum = 5 # in GeV
+momentum = 10 # in GeV
 #thetaMin = 90.25 # degrees
 #thetaMax = 90.25 # degrees
 thetaMin = 20 # degrees
@@ -60,7 +60,7 @@
 regiontool.volumeNames = ["CDCH"]
 regiontool.OutputLevel = DEBUG
 from GaudiKernel.SystemOfUnits import mm
-regiontool.maxStep = 0.4*mm
+#regiontool.maxStep = 0.4*mm
 
 from Configurables import SimG4UserLimitPhysicsList
 # create overlay on top of FTFP_BERT physics list, attaching fast sim/parametrization process
@@ -122,20 +122,35 @@
 dch_digitizer = DCHsimpleDigitizerExtendedEdm("DCHsimpleDigitizerExtendedEdm",
     inputSimHits = savetrackertool.SimTrackHits.Path,
     outputDigiHits = savetrackertool.SimTrackHits.Path.replace("sim", "digi"),
+    outputSimDigiAssociation = "DC_simDigiAssociation",
     readoutName = "CDCHHits",
     xyResolution = 0.1, # mm
     zResolution = 1, # mm
-    OutputLevel=INFO
+    debugMode = False,
+    OutputLevel = INFO
 )
 
+# Derive performance quantities
+#from Configurables import DCHsimpleDigitizerExtendedEdmPerformance
+#dch_perf = DCHsimpleDigitizerExtendedEdmPerformance("DCHsimpleDigitizerExtendedEdmPerformance",
+#        inputSimHits = savetrackertool.SimTrackHits.Path,
+#        inputDigiHits = dch_digitizer.outputDigiHits,
+#        inputSimDigiAssociation = dch_digitizer.outputSimDigiAssociation
+#)
+
+
+
 ################ Output
 from Configurables import PodioOutput
 out = PodioOutput("out",
                   OutputLevel=INFO)
 out.outputCommands = ["keep *"]
+if not dch_digitizer.debugMode:
+    out.outputCommands.append("drop *HitSimHitDelta*")
+    out.outputCommands.append("drop outputDigiLocalHits")
 
 import uuid
-out.filename = "output_simplifiedDriftChamber_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+"_stepLength_"+str(regiontool.maxStep)+".root"
+out.filename = "output_simplifiedDriftChamber_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+"_stepLength_default.root"
 
 #CPU information
 from Configurables import AuditorSvc, ChronoAuditor
@@ -159,10 +174,11 @@
               hepmc_converter,
               geantsim,
               dch_digitizer,
+              #dch_perf,
               out
               ],
     EvtSel = 'NONE',
-    EvtMax   = 10,
+    EvtMax   = 100,
     ExtSvc = [geoservice, podioevent, geantservice, audsvc],
     StopOnSignal = True,
  )