Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Correction of MTD ETL Validation: hits accumulation on a pixel-basis #43928

Closed
wants to merge 7 commits into from
Closed

Correction of MTD ETL Validation: hits accumulation on a pixel-basis #43928

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
618dc98
fix of ETL simHit accumulation on a pixel basis
Dec 8, 2023
0fb5456
cleaning
Dec 8, 2023
d880143
fix of ETL digiHits accumulation on a pixel basis
Dec 8, 2023
81f1b48
fix of ETL localReco Validation: resolution on a pixel basis
Dec 11, 2023
cdffeca
cleaning
Dec 11, 2023
2c29f78
fix in track validation>: ETL hits handled on a pixel basis
Dec 12, 2023
2c671c5
code format
Feb 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 30 additions & 4 deletions Validation/MtdValidation/plugins/EtlDigiHitsValidation.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@
#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
#include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"

#include "Geometry/MTDGeometryBuilder/interface/MTDGeomUtil.h"

class EtlDigiHitsValidation : public DQMEDAnalyzer {
public:
explicit EtlDigiHitsValidation(const edm::ParameterSet&);
Expand Down Expand Up @@ -78,8 +80,6 @@ class EtlDigiHitsValidation : public DQMEDAnalyzer {
MonitorElement* meHitQvsEta_[4];
MonitorElement* meHitTvsPhi_[4];
MonitorElement* meHitTvsEta_[4];

std::array<std::unordered_map<uint32_t, uint32_t>, 4> ndigiPerLGAD_;
};

// ------------ constructor and destructor --------------
Expand All @@ -97,13 +97,37 @@ EtlDigiHitsValidation::~EtlDigiHitsValidation() {}
void EtlDigiHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
using namespace edm;

using namespace mtd;
using namespace std;

// Struct to identify the pixel
struct ETLPixelId {
const uint32_t detid_;
const uint8_t row_;
const uint8_t col_;
ETLPixelId() : detid_(0), row_(0), col_(0) {}
ETLPixelId(const ETLDetId& id, uint8_t row, uint8_t col) : detid_(id.rawId()), row_(row), col_(col) {}
bool operator==(const ETLPixelId& other) const {
return detid_ == other.detid_ && row_ == other.row_ && col_ == other.col_;
}
};

struct PixelKey_hash {
size_t operator()(const ETLPixelId& key) const {
return std::hash<uint32_t>{}(key.detid_) ^ std::hash<uint8_t>{}(key.row_) ^ std::hash<uint8_t>{}(key.col_);
}
};

auto geometryHandle = iSetup.getTransientHandle(mtdgeoToken_);
const MTDGeometry* geom = geometryHandle.product();

auto etlDigiHitsHandle = makeValid(iEvent.getHandle(etlDigiHitsToken_));

// --- Loop over the ETL DIGI hits

std::array<std::unordered_map<ETLPixelId, uint32_t, PixelKey_hash>, 4> ndigiPerLGAD_;
MTDGeomUtil geomUtil;

unsigned int n_digi_etl[4] = {0, 0, 0, 0};
for (size_t i = 0; i < 4; i++) {
ndigiPerLGAD_[i].clear();
Expand Down Expand Up @@ -181,9 +205,11 @@ void EtlDigiHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSe

n_digi_etl[idet]++;
size_t ncount(0);
ndigiPerLGAD_[idet].emplace(detId.rawId(), ncount);
ndigiPerLGAD_[idet].at(detId.rawId())++;

std::pair<uint8_t, uint8_t> pixel = geomUtil.pixelInModule(detId, local_point);
ETLPixelId pixelId(detId.rawId(), pixel.first, pixel.second);
ndigiPerLGAD_[idet].emplace(pixelId, ncount);
ndigiPerLGAD_[idet].at(pixelId)++;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ctarricone the update of this class is not needed, because of the meaning of this counter: number of DIGIs per LGAD. Each pixel has its own separate DIGI, and there is one DIGI per pixel by construction, see https://github.com/cms-sw/cmssw/blob/master/SimFastTiming/FastTimingCommon/src/ETLDeviceSim.cc#L130
and the definition of MTDCellId.

Counting DIGIs per LGAD means counting how many pixels have fired in an LGAD for a given event, and this is what we are interested to monitor here. In the case of your class, you might see just either zero or one.

} // dataFrame loop

for (int i = 0; i < 4; i++) {
Expand Down
97 changes: 69 additions & 28 deletions Validation/MtdValidation/plugins/EtlLocalRecoValidation.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
#include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
#include "Geometry/MTDGeometryBuilder/interface/MTDGeomUtil.h"

#include "RecoLocalFastTime/Records/interface/MTDCPERecord.h"
#include "RecoLocalFastTime/FTLClusterizer/interface/MTDClusterParameterEstimator.h"
Expand Down Expand Up @@ -175,6 +176,27 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
using namespace edm;
using namespace std;
using namespace geant_units::operators;
using namespace mtd;

// Struct to identify the pixel
struct ETLPixelId {
const uint32_t detid_;
const uint8_t row_;
const uint8_t col_;
ETLPixelId() : detid_(0), row_(0), col_(0) {}
ETLPixelId(const ETLDetId& id, uint8_t row, uint8_t col) : detid_(id.rawId()), row_(row), col_(col) {}
bool operator==(const ETLPixelId& other) const {
return detid_ == other.detid_ && row_ == other.row_ && col_ == other.col_;
}
};

struct PixelKey_hash {
size_t operator()(const ETLPixelId& key) const {
return std::hash<uint32_t>{}(key.detid_) ^ std::hash<uint8_t>{}(key.row_) ^ std::hash<uint8_t>{}(key.col_);
}
};

MTDGeomUtil geomUtil;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the use of MTDGeomUtil is possible if geometry and topology are properly initialised with the set methods:

https://github.com/cms-sw/cmssw/blob/master/Geometry/MTDGeometryBuilder/interface/MTDGeomUtil.h#L28

For ETL at least the geometry needs to be set

https://github.com/cms-sw/cmssw/blob/master/Geometry/MTDGeometryBuilder/src/MTDGeomUtil.cc#L132

otherwise the code will happily crash. The comment is valid for all classes affected by this PR.


auto geometryHandle = iSetup.getTransientHandle(mtdgeoToken_);
const MTDGeometry* geom = geometryHandle.product();
Expand Down Expand Up @@ -202,7 +224,8 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
#endif

// --- Loop over the ETL SIM hits
std::unordered_map<uint32_t, MTDHit> m_etlSimHits[4];
//std::unordered_map<uint32_t, MTDHit> m_etlSimHits[4];
std::unordered_map<ETLPixelId, MTDHit, PixelKey_hash> m_etlSimHits[4];
for (auto const& simHit : etlSimHits) {
// --- Use only hits compatible with the in-time bunch-crossing
if (simHit.tof() < 0 || simHit.tof() > 25.)
Expand All @@ -225,7 +248,10 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
continue;
}

auto simHitIt = m_etlSimHits[idet].emplace(id.rawId(), MTDHit()).first;
//auto simHitIt = m_etlSimHits[idet].emplace(id.rawId(), MTDHit()).first;
std::pair<uint8_t, uint8_t> pixel = geomUtil.pixelInModule(id, simHit.localPosition());
ETLPixelId pixelId(id.rawId(), pixel.first, pixel.second);
auto simHitIt = m_etlSimHits[idet].emplace(pixelId, MTDHit()).first;

// --- Accumulate the energy (in MeV) of SIM hits in the same detector cell
(simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
Expand Down Expand Up @@ -312,16 +338,19 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
meHitTvsEta_[idet]->Fill(global_point.eta(), recHit.time());

// Resolution histograms
if (m_etlSimHits[idet].count(detId.rawId()) == 1) {
if (m_etlSimHits[idet][detId.rawId()].energy > hitMinEnergy2Dis_) {
float time_res = recHit.time() - m_etlSimHits[idet][detId.rawId()].time;
float energy_res = recHit.energy() - m_etlSimHits[idet][detId.rawId()].energy;
std::pair<uint8_t, uint8_t> pixel = geomUtil.pixelInModule(detId, local_point);
ETLPixelId pixelId(detId.rawId(), pixel.first, pixel.second);

if (m_etlSimHits[idet].count(pixelId) == 1) {
if (m_etlSimHits[idet][pixelId].energy > hitMinEnergy2Dis_) {
float time_res = recHit.time() - m_etlSimHits[idet][pixelId].time;
float energy_res = recHit.energy() - m_etlSimHits[idet][pixelId].energy;

meTimeRes_->Fill(time_res);
meEnergyRes_->Fill(energy_res);

meTPullvsEta_->Fill(std::abs(global_point.eta()), time_res / recHit.timeError());
meTPullvsE_->Fill(m_etlSimHits[idet][detId.rawId()].energy, time_res / recHit.timeError());
meTPullvsE_->Fill(m_etlSimHits[idet][pixelId].energy, time_res / recHit.timeError());
}
}

Expand Down Expand Up @@ -401,13 +430,23 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS

// Match the RECO hit to the corresponding SIM hit
for (const auto& recHit : *etlRecHitsHandle) {
ETLDetId hitId(recHit.id().rawId());
ETLDetId detId(recHit.id().rawId());

DetId geoId = detId.geographicalId();
const MTDGeomDet* thedet = geom->idToDet(geoId);
const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());

if (m_etlSimHits[idet].count(hitId.rawId()) == 0)
Local3DPoint local_point(topo.localX(recHit.row()), topo.localY(recHit.column()), 0.);

std::pair<uint8_t, uint8_t> pixel = geomUtil.pixelInModule(detId, local_point);
ETLPixelId pixelId(detId.rawId(), pixel.first, pixel.second);

if (m_etlSimHits[idet].count(pixelId) == 0)
continue;

// Check the hit position
if (hitId.zside() != cluId.zside() || hitId.mtdRR() != cluId.mtdRR() || hitId.module() != cluId.module() ||
if (detId.zside() != cluId.zside() || detId.mtdRR() != cluId.mtdRR() || detId.module() != cluId.module() ||
recHit.row() != hit_row || recHit.column() != hit_col)
continue;

Expand All @@ -416,18 +455,18 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
continue;

// SIM hit's position in the module reference frame
Local3DPoint local_point_sim(convertMmToCm(m_etlSimHits[idet][recHit.id().rawId()].x),
convertMmToCm(m_etlSimHits[idet][recHit.id().rawId()].y),
convertMmToCm(m_etlSimHits[idet][recHit.id().rawId()].z));
Local3DPoint local_point_sim(convertMmToCm(m_etlSimHits[idet][pixelId].x),
convertMmToCm(m_etlSimHits[idet][pixelId].y),
convertMmToCm(m_etlSimHits[idet][pixelId].z));

// Calculate the SIM cluster's position in the module reference frame
cluLocXSIM += local_point_sim.x() * m_etlSimHits[idet][recHit.id().rawId()].energy;
cluLocYSIM += local_point_sim.y() * m_etlSimHits[idet][recHit.id().rawId()].energy;
cluLocZSIM += local_point_sim.z() * m_etlSimHits[idet][recHit.id().rawId()].energy;
cluLocXSIM += local_point_sim.x() * m_etlSimHits[idet][pixelId].energy;
cluLocYSIM += local_point_sim.y() * m_etlSimHits[idet][pixelId].energy;
cluLocZSIM += local_point_sim.z() * m_etlSimHits[idet][pixelId].energy;

// Calculate the SIM cluster energy and time
cluEneSIM += m_etlSimHits[idet][recHit.id().rawId()].energy;
cluTimeSIM += m_etlSimHits[idet][recHit.id().rawId()].time * m_etlSimHits[idet][recHit.id().rawId()].energy;
cluEneSIM += m_etlSimHits[idet][pixelId].energy;
cluTimeSIM += m_etlSimHits[idet][pixelId].time * m_etlSimHits[idet][pixelId].energy;

break;

Expand Down Expand Up @@ -500,31 +539,33 @@ void EtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS

for (const auto& uRecHit : *etlUncalibRecHitsHandle) {
ETLDetId detId = uRecHit.id();

int idet = detId.zside() + detId.nDisc();

// --- Skip UncalibratedRecHits not matched to SimHits
if (m_etlSimHits[idet].count(detId.rawId()) != 1)
continue;

DetId geoId = detId.geographicalId();
const MTDGeomDet* thedet = geom->idToDet(geoId);
if (thedet == nullptr)
throw cms::Exception("EtlLocalRecoValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
<< detId.rawId() << ") is invalid!" << std::dec << std::endl;

const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());

Local3DPoint local_point(topo.localX(uRecHit.row()), topo.localY(uRecHit.column()), 0.);
const auto& global_point = thedet->toGlobal(local_point);

std::pair<uint8_t, uint8_t> pixel = geomUtil.pixelInModule(detId, local_point);
ETLPixelId pixelId(detId.rawId(), pixel.first, pixel.second);

// --- Skip UncalibratedRecHits not matched to SimHits
if (m_etlSimHits[idet].count(pixelId) == 0)
continue;

if (thedet == nullptr)
throw cms::Exception("EtlLocalRecoValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
<< detId.rawId() << ") is invalid!" << std::dec << std::endl;

// --- Fill the histograms

if (uRecHit.amplitude().first < hitMinAmplitude_)
continue;

float time_res = uRecHit.time().first - m_etlSimHits[idet][detId.rawId()].time;
float time_res = uRecHit.time().first - m_etlSimHits[idet][pixelId].time;

int iside = (detId.zside() == -1 ? 0 : 1);

Expand Down
40 changes: 34 additions & 6 deletions Validation/MtdValidation/plugins/EtlSimHitsValidation.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"

#include "Geometry/MTDGeometryBuilder/interface/MTDGeomUtil.h"

#include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
Expand Down Expand Up @@ -101,14 +103,37 @@ void EtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet
using namespace edm;
using namespace geant_units::operators;

using namespace mtd;
using namespace std;

// Struct to identify the pixel
struct ETLPixelId {
const uint32_t detid_;
const uint8_t row_;
const uint8_t col_;
ETLPixelId() : detid_(0), row_(0), col_(0) {}
ETLPixelId(const ETLDetId& id, uint8_t row, uint8_t col) : detid_(id.rawId()), row_(row), col_(col) {}
bool operator==(const ETLPixelId& other) const {
return detid_ == other.detid_ && row_ == other.row_ && col_ == other.col_;
}
};

struct PixelKey_hash {
size_t operator()(const ETLPixelId& key) const {
return std::hash<uint32_t>{}(key.detid_) ^ std::hash<uint8_t>{}(key.row_) ^ std::hash<uint8_t>{}(key.col_);
}
};

auto geometryHandle = iSetup.getTransientHandle(mtdgeoToken_);
const MTDGeometry* geom = geometryHandle.product();

MTDGeomUtil geomUtil;

auto etlSimHitsHandle = makeValid(iEvent.getHandle(etlSimHitsToken_));
MixCollection<PSimHit> etlSimHits(etlSimHitsHandle.product());

std::unordered_map<uint32_t, MTDHit> m_etlHits[4];
std::unordered_map<uint32_t, std::set<int> > m_etlTrkPerCell[4];
std::unordered_map<ETLPixelId, MTDHit, PixelKey_hash> m_etlHits[4];
std::unordered_map<ETLPixelId, std::set<int>, PixelKey_hash> m_etlTrkPerCell[4];

// --- Loop over the ETL SIM hits

Expand Down Expand Up @@ -139,9 +164,10 @@ void EtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet
continue;
}

m_etlTrkPerCell[idet][id.rawId()].insert(simHit.trackId());

auto simHitIt = m_etlHits[idet].emplace(id.rawId(), MTDHit()).first;
std::pair<uint8_t, uint8_t> pixel = geomUtil.pixelInModule(id, simHit.localPosition());
ETLPixelId pixelId(id.rawId(), pixel.first, pixel.second);
m_etlTrkPerCell[idet][pixelId].insert(simHit.trackId());
auto simHitIt = m_etlHits[idet].emplace(pixelId, MTDHit()).first;

// --- Accumulate the energy (in MeV) of SIM hits in the same detector cell
(simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
Expand Down Expand Up @@ -176,7 +202,9 @@ void EtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSet
if ((hit.second).energy < hitMinEnergy2Dis_)
continue;
// --- Get the SIM hit global position
ETLDetId detId(hit.first);

ETLDetId detId;
detId = hit.first.detid_;
DetId geoId = detId.geographicalId();
const MTDGeomDet* thedet = geom->idToDet(geoId);
if (thedet == nullptr)
Expand Down
35 changes: 30 additions & 5 deletions Validation/MtdValidation/plugins/MtdTracksValidation.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
#include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
#include "Geometry/MTDGeometryBuilder/interface/MTDGeomUtil.h"

#include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h"
#include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h"
Expand Down Expand Up @@ -296,14 +297,35 @@ void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetu
using namespace edm;
using namespace geant_units::operators;
using namespace std;
using namespace mtd;

// Struct to identify the pixel
struct ETLPixelId {
const uint32_t detid_;
const uint8_t row_;
const uint8_t col_;
ETLPixelId() : detid_(0), row_(0), col_(0) {}
ETLPixelId(const ETLDetId& id, uint8_t row, uint8_t col) : detid_(id.rawId()), row_(row), col_(col) {}
bool operator==(const ETLPixelId& other) const {
return detid_ == other.detid_ && row_ == other.row_ && col_ == other.col_;
}
};

struct PixelKey_hash {
size_t operator()(const ETLPixelId& key) const {
return std::hash<uint32_t>{}(key.detid_) ^ std::hash<uint8_t>{}(key.row_) ^ std::hash<uint8_t>{}(key.col_);
}
};

MTDGeomUtil geomUtil;

auto GenRecTrackHandle = makeValid(iEvent.getHandle(GenRecTrackToken_));
auto RecVertexHandle = makeValid(iEvent.getHandle(RecVertexToken_));

std::unordered_map<uint32_t, MTDHit> m_btlHits;
std::unordered_map<uint32_t, MTDHit> m_etlHits;
std::unordered_map<ETLPixelId, MTDHit, PixelKey_hash> m_etlHits;
std::unordered_map<uint32_t, std::set<unsigned long int>> m_btlTrkPerCell;
std::unordered_map<uint32_t, std::set<unsigned long int>> m_etlTrkPerCell;
std::unordered_map<ETLPixelId, std::set<unsigned long int>, PixelKey_hash> m_etlTrkPerCell;
std::map<TrackingParticleRef, std::vector<uint32_t>> m_tp2detid;

const auto& tMtd = iEvent.get(tmtdToken_);
Expand Down Expand Up @@ -358,8 +380,11 @@ void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetu
}
DetId id = simHit.detUnitId();
auto const thisHId = uniqueId(simHit.trackId(), simHit.eventId());
m_etlTrkPerCell[id.rawId()].insert(thisHId);
auto simHitIt = m_etlHits.emplace(id.rawId(), MTDHit()).first;
std::pair<uint8_t, uint8_t> pixel = geomUtil.pixelInModule(id, simHit.localPosition());
ETLPixelId pixelId(id.rawId(), pixel.first, pixel.second);
m_etlTrkPerCell[pixelId].insert(thisHId);
auto simHitIt = m_etlHits.emplace(pixelId, MTDHit()).first;

// --- Accumulate the energy (in MeV) of SIM hits in the same detector cell
(simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
}
Expand Down Expand Up @@ -393,7 +418,7 @@ void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetu
}
for (auto const& simtrack : cell.second) {
if (thisTId == simtrack) {
m_tp2detid[tpref].emplace_back(cell.first);
m_tp2detid[tpref].emplace_back(cell.first.detid_);
break;
}
}
Expand Down