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