Skip to content

Commit

Permalink
Add debug collections and optional local coordinate collection
Browse files Browse the repository at this point in the history
  • Loading branch information
@BrieucF
BrieucF committed Feb 23, 2024
1 parent fba9898 commit b34b20d
Show file tree
Hide file tree
Showing 3 changed files with 76 additions and 16 deletions.
13 changes: 12 additions & 1 deletion DCHdigi/include/DCHsimpleDigitizerExtendedEdm.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,11 +10,13 @@
#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
Expand Down Expand Up @@ -55,6 +57,8 @@ class DCHsimpleDigitizerExtendedEdm : public GaudiAlgorithm {
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"};
Expand All @@ -70,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;
Expand Down
73 changes: 59 additions & 14 deletions DCHdigi/src/DCHsimpleDigitizerExtendedEdm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,11 +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;

// 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
extension::DriftChamberDigiCollection* 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 wire local frame to smear in the direction aligned/perpendicular with the wire for z/distanceToWire, taking e.g. stereo angle into account
Expand All @@ -67,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};
Expand All @@ -83,35 +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;
debug() << "Original distance to wire: " << simHitLocalPositionVector.rho() << " mm" << endmsg;
double smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm;
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
double leftHitLocalPosition[3] = {-1 * smearedDistanceToWire, 0, smearedZ};
double rightHitLocalPosition[3] = {smearedDistanceToWire, 0, smearedZ};
// 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.setLeftPosition(edm4hep::Vector3d({leftHitGlobalPosition[0] / dd4hep::mm, leftHitGlobalPosition[1] / dd4hep::mm, leftHitGlobalPosition[2] / dd4hep::mm}));
output_digi_hit.setRightPosition(edm4hep::Vector3d({rightHitGlobalPosition[0] / dd4hep::mm, rightHitGlobalPosition[1] / dd4hep::mm, rightHitGlobalPosition[2] / 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;
Expand Down
6 changes: 5 additions & 1 deletion DCHdigi/test/runDCHsimpleDigitizerExtendedEdm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,8 @@
readoutName = "CDCHHits",
xyResolution = 0.1, # mm
zResolution = 1, # mm
OutputLevel=INFO
debugMode = False,
OutputLevel = INFO
)

# Derive performance quantities
Expand All @@ -144,6 +145,9 @@
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_default.root"
Expand Down

0 comments on commit b34b20d

Please sign in to comment.