Adding validation for the MUON digitizer

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
+      - edm4hep::SimTrackerHit sim // reference to the simulated hit
+

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<podio::UserDataCollection<double>> m_simDigiDifferenceX{"simDigiDifferenceX", Gaudi::DataHandle::Writer, this}; // mm
+  DataHandle<podio::UserDataCollection<double>> m_simDigiDifferenceY{"simDigiDifferenceY", Gaudi::DataHandle::Writer, this}; // mm
+  DataHandle<podio::UserDataCollection<double>> 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;
 };

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; }

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,
 )