File indexing completed on 2024-12-12 23:19:34
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <string>
0015
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020
0021 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0022 #include "DQMServices/Core/interface/DQMStore.h"
0023
0024 #include "DataFormats/Common/interface/ValidHandle.h"
0025 #include "DataFormats/Math/interface/GeantUnits.h"
0026 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0027
0028 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0029 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0030 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0031
0032 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0033 #include "Geometry/Records/interface/MTDTopologyRcd.h"
0034 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0035 #include "Geometry/MTDGeometryBuilder/interface/MTDTopology.h"
0036
0037 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0038 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0039
0040 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0041
0042 #include "MTDHit.h"
0043
0044 class BtlSimHitsValidation : public DQMEDAnalyzer {
0045 public:
0046 explicit BtlSimHitsValidation(const edm::ParameterSet&);
0047 ~BtlSimHitsValidation() override;
0048
0049 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 private:
0052 void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0053
0054 void analyze(const edm::Event&, const edm::EventSetup&) override;
0055
0056
0057
0058 const std::string folder_;
0059 const float hitMinEnergy_;
0060
0061 static constexpr float BXTime_ = 25.;
0062
0063 edm::EDGetTokenT<CrossingFrame<PSimHit> > btlSimHitsToken_;
0064
0065 edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0066 edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0067
0068
0069
0070 MonitorElement* meNevents_;
0071
0072 MonitorElement* meNhits_;
0073 MonitorElement* meNtrkPerCell_;
0074
0075 MonitorElement* meHitEnergy_;
0076 MonitorElement* meHitLogEnergy_;
0077 MonitorElement* meHitTime_;
0078
0079 MonitorElement* meHitXlocal_;
0080 MonitorElement* meHitYlocal_;
0081 MonitorElement* meHitZlocal_;
0082
0083 MonitorElement* meOccupancy_;
0084
0085 MonitorElement* meHitX_;
0086 MonitorElement* meHitY_;
0087 MonitorElement* meHitZ_;
0088 MonitorElement* meHitPhi_;
0089 MonitorElement* meHitEta_;
0090
0091 MonitorElement* meHitTvsE_;
0092 MonitorElement* meHitEvsPhi_;
0093 MonitorElement* meHitEvsEta_;
0094 MonitorElement* meHitEvsZ_;
0095 MonitorElement* meHitTvsPhi_;
0096 MonitorElement* meHitTvsEta_;
0097 MonitorElement* meHitTvsZ_;
0098
0099 MonitorElement* meHitTrkID1_;
0100 MonitorElement* meHitTrkID2_;
0101 MonitorElement* meHitTrkID3_;
0102 MonitorElement* meHitTrkID4_;
0103
0104 MonitorElement* meHitDeltaT2_;
0105 MonitorElement* meHitDeltaT3_;
0106 MonitorElement* meHitDeltaT4_;
0107
0108 MonitorElement* meHitDeltaE12vsE1_;
0109
0110 static constexpr float cellEneCut_ = 10.;
0111
0112 MonitorElement* meHitE1overEcellBulk1_;
0113 MonitorElement* meHitE1overEcellBulk2_;
0114 MonitorElement* meHitE1overEcellBulk3_;
0115 MonitorElement* meHitE1overEcellBulk4_;
0116 MonitorElement* meHitE1overEcellTail1_;
0117 MonitorElement* meHitE1overEcellTail2_;
0118 MonitorElement* meHitE1overEcellTail3_;
0119 MonitorElement* meHitE1overEcellTail4_;
0120 };
0121
0122
0123 BtlSimHitsValidation::BtlSimHitsValidation(const edm::ParameterSet& iConfig)
0124 : folder_(iConfig.getParameter<std::string>("folder")),
0125 hitMinEnergy_(iConfig.getParameter<double>("hitMinimumEnergy")) {
0126 btlSimHitsToken_ = consumes<CrossingFrame<PSimHit> >(iConfig.getParameter<edm::InputTag>("inputTag"));
0127 mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>();
0128 mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>();
0129 }
0130
0131 BtlSimHitsValidation::~BtlSimHitsValidation() {}
0132
0133
0134 void BtlSimHitsValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0135 using namespace edm;
0136 using namespace geant_units::operators;
0137
0138 auto geometryHandle = iSetup.getTransientHandle(mtdgeoToken_);
0139 const MTDGeometry* geom = geometryHandle.product();
0140
0141 auto topologyHandle = iSetup.getTransientHandle(mtdtopoToken_);
0142 const MTDTopology* topology = topologyHandle.product();
0143
0144 auto btlSimHitsHandle = makeValid(iEvent.getHandle(btlSimHitsToken_));
0145 MixCollection<PSimHit> btlSimHits(btlSimHitsHandle.product());
0146
0147 std::unordered_map<uint32_t, std::unordered_map<uint64_t, MTDHit> > m_btlHitsPerCell;
0148
0149
0150 for (auto const& simHit : btlSimHits) {
0151
0152 if (simHit.tof() < 0. || simHit.tof() > BXTime_) {
0153 continue;
0154 }
0155
0156 DetId id = simHit.detUnitId();
0157
0158
0159 uint64_t globalTrkID = ((uint64_t)simHit.eventId().rawId() << 32) | simHit.trackId();
0160
0161
0162 m_btlHitsPerCell[id.rawId()][globalTrkID].energy += convertGeVToMeV(simHit.energyLoss());
0163
0164
0165 if (m_btlHitsPerCell[id.rawId()][globalTrkID].time == 0. ||
0166 simHit.tof() < m_btlHitsPerCell[id.rawId()][globalTrkID].time) {
0167 m_btlHitsPerCell[id.rawId()][globalTrkID].time = simHit.tof();
0168
0169 auto hit_pos = simHit.localPosition();
0170 m_btlHitsPerCell[id.rawId()][globalTrkID].x = hit_pos.x();
0171 m_btlHitsPerCell[id.rawId()][globalTrkID].y = hit_pos.y();
0172 m_btlHitsPerCell[id.rawId()][globalTrkID].z = hit_pos.z();
0173 }
0174
0175 }
0176
0177
0178
0179
0180
0181 if (!m_btlHitsPerCell.empty()) {
0182 meNhits_->Fill(log10(m_btlHitsPerCell.size()));
0183 }
0184
0185
0186 for (auto const& cell : m_btlHitsPerCell) {
0187
0188 const auto& m_hits = cell.second;
0189
0190
0191 if (m_hits.empty()) {
0192 continue;
0193 }
0194
0195
0196 std::vector<uint64_t> v_hitID;
0197 float ene_tot_cell = 0.;
0198
0199 for (auto const& hit : m_hits) {
0200 ene_tot_cell += hit.second.energy;
0201 v_hitID.push_back(hit.first);
0202 }
0203
0204 meHitLogEnergy_->Fill(log10(ene_tot_cell));
0205
0206
0207 if (ene_tot_cell < hitMinEnergy_) {
0208 continue;
0209 }
0210
0211
0212 bool swapped;
0213 for (unsigned int ihit = 0; ihit < v_hitID.size() - 1; ++ihit) {
0214 swapped = false;
0215 for (unsigned int jhit = 0; jhit < v_hitID.size() - ihit - 1; ++jhit) {
0216 if (m_hits.at(v_hitID[jhit]).time > m_hits.at(v_hitID[jhit + 1]).time) {
0217 std::swap(v_hitID[jhit], v_hitID[jhit + 1]);
0218 swapped = true;
0219 }
0220 }
0221 if (swapped == false) {
0222 break;
0223 }
0224 }
0225
0226
0227 float deltaT_max = 0.;
0228 for (unsigned int ihit = 0; ihit < v_hitID.size() - 1; ++ihit) {
0229 float deltaT = m_hits.at(v_hitID[ihit + 1]).time - m_hits.at(v_hitID[ihit]).time;
0230
0231 if (deltaT > deltaT_max) {
0232 deltaT_max = deltaT;
0233 }
0234 }
0235
0236
0237 BTLDetId detId(cell.first);
0238 DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode()));
0239 const MTDGeomDet* thedet = geom->idToDet(geoId);
0240 if (thedet == nullptr) {
0241 throw cms::Exception("BtlSimHitsValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0242 << detId.rawId() << ") is invalid!" << std::dec << std::endl;
0243 }
0244 const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0245 const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0246
0247 Local3DPoint local_point(convertMmToCm(m_hits.at(v_hitID[0]).x),
0248 convertMmToCm(m_hits.at(v_hitID[0]).y),
0249 convertMmToCm(m_hits.at(v_hitID[0]).z));
0250
0251 local_point = topo.pixelToModuleLocalPoint(local_point, detId.row(topo.nrows()), detId.column(topo.nrows()));
0252 const auto& global_point = thedet->toGlobal(local_point);
0253
0254
0255
0256 meHitEnergy_->Fill(ene_tot_cell);
0257 meHitTime_->Fill(m_hits.at(v_hitID[0]).time);
0258
0259 meHitXlocal_->Fill(m_hits.at(v_hitID[0]).x);
0260 meHitYlocal_->Fill(m_hits.at(v_hitID[0]).y);
0261 meHitZlocal_->Fill(m_hits.at(v_hitID[0]).z);
0262
0263 meOccupancy_->Fill(global_point.z(), global_point.phi());
0264
0265 meHitX_->Fill(global_point.x());
0266 meHitY_->Fill(global_point.y());
0267 meHitZ_->Fill(global_point.z());
0268 meHitPhi_->Fill(global_point.phi());
0269 meHitEta_->Fill(global_point.eta());
0270
0271 meHitTvsE_->Fill(ene_tot_cell, m_hits.at(v_hitID[0]).time);
0272 meHitEvsPhi_->Fill(global_point.phi(), ene_tot_cell);
0273 meHitEvsEta_->Fill(global_point.eta(), ene_tot_cell);
0274 meHitEvsZ_->Fill(global_point.z(), ene_tot_cell);
0275 meHitTvsPhi_->Fill(global_point.phi(), m_hits.at(v_hitID[0]).time);
0276 meHitTvsEta_->Fill(global_point.eta(), m_hits.at(v_hitID[0]).time);
0277 meHitTvsZ_->Fill(global_point.z(), m_hits.at(v_hitID[0]).time);
0278
0279 meNtrkPerCell_->Fill(v_hitID.size());
0280
0281
0282 int trackID = (int)(v_hitID[0] & 0xFFFFFFFF) / 100000000;
0283 meHitTrkID1_->Fill(trackID);
0284
0285
0286 if (ene_tot_cell < cellEneCut_) {
0287 meHitE1overEcellBulk1_->Fill(m_hits.at(v_hitID[0]).energy / ene_tot_cell);
0288 }
0289
0290 else {
0291 meHitE1overEcellTail1_->Fill(m_hits.at(v_hitID[0]).energy / ene_tot_cell);
0292 }
0293
0294
0295 if (v_hitID.size() == 2) {
0296 trackID = (int)(v_hitID[1] & 0xFFFFFFFF) / 100000000;
0297 meHitTrkID2_->Fill(trackID);
0298
0299 meHitDeltaT2_->Fill(deltaT_max);
0300
0301
0302 if (ene_tot_cell < cellEneCut_) {
0303 meHitE1overEcellBulk2_->Fill(m_hits.at(v_hitID[1]).energy / ene_tot_cell);
0304 }
0305
0306 else {
0307 meHitE1overEcellTail2_->Fill(m_hits.at(v_hitID[1]).energy / ene_tot_cell);
0308 }
0309
0310 meHitDeltaE12vsE1_->Fill(m_hits.at(v_hitID[0]).energy,
0311 m_hits.at(v_hitID[0]).energy - m_hits.at(v_hitID[1]).energy);
0312
0313 }
0314
0315 else if (v_hitID.size() == 3) {
0316 trackID = (int)(v_hitID[2] & 0xFFFFFFFF) / 100000000;
0317 meHitTrkID3_->Fill(trackID);
0318
0319 meHitDeltaT3_->Fill(deltaT_max);
0320
0321
0322 if (ene_tot_cell < cellEneCut_) {
0323 meHitE1overEcellBulk3_->Fill(m_hits.at(v_hitID[2]).energy / ene_tot_cell);
0324 }
0325
0326 else {
0327 meHitE1overEcellTail3_->Fill(m_hits.at(v_hitID[2]).energy / ene_tot_cell);
0328 }
0329
0330 }
0331
0332 else if (v_hitID.size() >= 4) {
0333 for (unsigned int ihit = 3; ihit < v_hitID.size(); ++ihit) {
0334 trackID = (int)(v_hitID[ihit] & 0xFFFFFFFF) / 100000000;
0335 meHitTrkID4_->Fill(trackID);
0336
0337
0338 if (ene_tot_cell < cellEneCut_) {
0339 meHitE1overEcellBulk4_->Fill(m_hits.at(v_hitID[ihit]).energy / ene_tot_cell);
0340 }
0341
0342 else {
0343 meHitE1overEcellTail4_->Fill(m_hits.at(v_hitID[ihit]).energy / ene_tot_cell);
0344 }
0345 }
0346
0347 meHitDeltaT4_->Fill(deltaT_max);
0348 }
0349
0350 }
0351
0352
0353 meNevents_->Fill(0.5);
0354 }
0355
0356
0357 void BtlSimHitsValidation::bookHistograms(DQMStore::IBooker& ibook,
0358 edm::Run const& run,
0359 edm::EventSetup const& iSetup) {
0360 ibook.setCurrentFolder(folder_);
0361
0362
0363
0364 meNevents_ = ibook.book1D("BtlNevents", "Number of events", 1, 0., 1.);
0365
0366 meNhits_ = ibook.book1D("BtlNhits", "Number of BTL cells with SIM hits;log_{10}(N_{BTL cells})", 100, 0., 5.25);
0367 meNtrkPerCell_ = ibook.book1D("BtlNtrkPerCell", "Number of tracks per BTL cell;N_{trk}", 10, 0., 10.);
0368
0369 meHitEnergy_ = ibook.book1D("BtlHitEnergy", "BTL SIM hits energy;E_{SIM} [MeV]", 100, 0., 20.);
0370 meHitLogEnergy_ = ibook.book1D("BtlHitLogEnergy", "BTL SIM hits energy;log_{10}(E_{SIM} [MeV])", 200, -6., 3.);
0371 meHitTime_ = ibook.book1D("BtlHitTime", "BTL SIM hits ToA;ToA_{SIM} [ns]", 100, 0., 25.);
0372
0373 meHitXlocal_ = ibook.book1D("BtlHitXlocal", "BTL SIM local X;X_{SIM}^{LOC} [mm]", 100, -30., 30.);
0374 meHitYlocal_ = ibook.book1D("BtlHitYlocal", "BTL SIM local Y;Y_{SIM}^{LOC} [mm]", 100, -1.65, 1.65);
0375 meHitZlocal_ = ibook.book1D("BtlHitZlocal", "BTL SIM local Z;Z_{SIM}^{LOC} [mm]", 100, -2., 2.);
0376
0377 meOccupancy_ = ibook.book2D(
0378 "BtlOccupancy", "BTL SIM hits occupancy;z_{SIM} [cm];#phi_{SIM} [rad]", 130, -260., 260., 200, -3.15, 3.15);
0379
0380 meHitX_ = ibook.book1D("BtlHitX", "BTL SIM hits X;X_{SIM} [cm]", 100, -120., 120.);
0381 meHitY_ = ibook.book1D("BtlHitY", "BTL SIM hits Y;Y_{SIM} [cm]", 100, -120., 120.);
0382 meHitZ_ = ibook.book1D("BtlHitZ", "BTL SIM hits Z;Z_{SIM} [cm]", 100, -260., 260.);
0383 meHitPhi_ = ibook.book1D("BtlHitPhi", "BTL SIM hits #phi;#phi_{SIM} [rad]", 200, -3.15, 3.15);
0384 meHitEta_ = ibook.book1D("BtlHitEta", "BTL SIM hits #eta;#eta_{SIM}", 100, -1.55, 1.55);
0385
0386 meHitTvsE_ =
0387 ibook.bookProfile("BtlHitTvsE", "BTL SIM time vs energy;E_{SIM} [MeV];T_{SIM} [ns]", 50, 0., 20., 0., 100.);
0388 meHitEvsPhi_ = ibook.bookProfile(
0389 "BtlHitEvsPhi", "BTL SIM energy vs #phi;#phi_{SIM} [rad];E_{SIM} [MeV]", 50, -3.15, 3.15, 0., 100.);
0390 meHitEvsEta_ =
0391 ibook.bookProfile("BtlHitEvsEta", "BTL SIM energy vs #eta;#eta_{SIM};E_{SIM} [MeV]", 50, -1.55, 1.55, 0., 100.);
0392 meHitEvsZ_ =
0393 ibook.bookProfile("BtlHitEvsZ", "BTL SIM energy vs Z;Z_{SIM} [cm];E_{SIM} [MeV]", 50, -260., 260., 0., 100.);
0394 meHitTvsPhi_ = ibook.bookProfile(
0395 "BtlHitTvsPhi", "BTL SIM time vs #phi;#phi_{SIM} [rad];T_{SIM} [ns]", 50, -3.15, 3.15, 0., 100.);
0396 meHitTvsEta_ =
0397 ibook.bookProfile("BtlHitTvsEta", "BTL SIM time vs #eta;#eta_{SIM};T_{SIM} [ns]", 50, -1.55, 1.55, 0., 100.);
0398 meHitTvsZ_ =
0399 ibook.bookProfile("BtlHitTvsZ", "BTL SIM time vs Z;Z_{SIM} [cm];T_{SIM} [ns]", 50, -260., 260., 0., 100.);
0400
0401 meHitTrkID1_ = ibook.book1I("BtlHitTrkID1", "Category of the 1^{st} hit in time;Hit category", 8, 0, 8);
0402 meHitTrkID2_ = ibook.book1I("BtlHitTrkID2", "Category of the 2^{nd} hit in time;Hit category", 8, 0, 8);
0403 meHitTrkID3_ = ibook.book1I("BtlHitTrkID3", "Category of the 3^{rd} hit in time;Hit category", 8, 0, 8);
0404 meHitTrkID4_ = ibook.book1I("BtlHitTrkID4", "Category of the #geq4^{th} hit in time;Hit category", 8, 0, 8);
0405
0406 meHitDeltaT2_ = ibook.book1D(
0407 "BtlHitDeltaT2", "Time interval between hits in the same cell (2 hits);#DeltaT_{2} [ns]", 100., 0., 25.);
0408 meHitDeltaT3_ = ibook.book1D(
0409 "BtlHitDeltaT3", "Max time interval between hits in the same cell (3 hits);#DeltaT_{3} [ns]", 100., 0., 25.);
0410 meHitDeltaT4_ = ibook.book1D("BtlHitDeltaT4",
0411 "Max time interval between hits in the same cell (#geq4 hits);#DeltaT_{#geq4} [ns]",
0412 100.,
0413 0.,
0414 25.);
0415
0416 meHitDeltaE12vsE1_ = ibook.bookProfile(
0417 "BtlHitDeltaE12vsE", "E_{1}-E_{2} vs E_{1};E_{1} [MeV];E_{1}-E_{2} [MeV]", 100, 0., 100., 0., 100.);
0418
0419 meHitE1overEcellBulk1_ = ibook.book1D("BtlHitE1overEtotBulk1",
0420 "Energy fraction of the 1^{st} hit in time (E_{cell}<10 MeV);E_{1} / E_{cell}",
0421 100,
0422 0.,
0423 1.);
0424 meHitE1overEcellBulk2_ = ibook.book1D("BtlHitE1overEtotBulk2",
0425 "Energy fraction of the 2^{nd} hit in time (E_{cell}<10 MeV);E_{2} / E_{cell}",
0426 100,
0427 0.,
0428 1.);
0429 meHitE1overEcellBulk3_ = ibook.book1D("BtlHitE1overEtotBulk3",
0430 "Energy fraction of the 3^{rd} hit in time (E_{cell}<10 MeV);E_{3} / E_{cell}",
0431 100,
0432 0.,
0433 1.);
0434 meHitE1overEcellBulk4_ =
0435 ibook.book1D("BtlHitE1overEtotBulk4",
0436 "Energy fraction of the #geq4^{th} hit in time (E_{cell}<10 MeV);E_{#geq4} / E_{cell}",
0437 100,
0438 0.,
0439 1.);
0440 meHitE1overEcellTail1_ = ibook.book1D("BtlHitE1overEtotTail1",
0441 "Energy fraction of the 1^{st} hit in time (E_{cell}>10 MeV);E_{1} / E_{cell}",
0442 100,
0443 0.,
0444 1.);
0445 meHitE1overEcellTail2_ = ibook.book1D("BtlHitE1overEtotTail2",
0446 "Energy fraction of the 2^{nd} hit in time (E_{cell}>10 MeV);E_{2} / E_{cell}",
0447 100,
0448 0.,
0449 1.);
0450 meHitE1overEcellTail3_ = ibook.book1D("BtlHitE1overEtotTail3",
0451 "Energy fraction of the 3^{rd} hit in time (E_{cell}>10 MeV);E_{3} / E_{cell}",
0452 100,
0453 0.,
0454 1.);
0455 meHitE1overEcellTail4_ =
0456 ibook.book1D("BtlHitE1overEtotTail4",
0457 "Energy fraction of the #geq4^{th} hit in time (E_{cell}>10 MeV);E_{#geq4} / E_{cell}",
0458 100,
0459 0.,
0460 1.);
0461 }
0462
0463
0464 void BtlSimHitsValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0465 edm::ParameterSetDescription desc;
0466
0467 desc.add<std::string>("folder", "MTD/BTL/SimHits");
0468 desc.add<edm::InputTag>("inputTag", edm::InputTag("mix", "g4SimHitsFastTimerHitsBarrel"));
0469 desc.add<double>("hitMinimumEnergy", 1.);
0470
0471 descriptions.add("btlSimHitsValid", desc);
0472 }
0473
0474 DEFINE_FWK_MODULE(BtlSimHitsValidation);