File indexing completed on 2023-06-02 00:51:03
0001 #include "Validation/HGCalValidation/interface/HGCalValidator.h"
0002
0003 #include "SimCalorimetry/HGCalAssociatorProducers/interface/AssociatorTools.h"
0004
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 using namespace std;
0009 using namespace edm;
0010
0011 HGCalValidator::HGCalValidator(const edm::ParameterSet& pset)
0012 : caloGeomToken_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0013 label_lcl(pset.getParameter<edm::InputTag>("label_lcl")),
0014 label_tst(pset.getParameter<std::vector<edm::InputTag>>("label_tst")),
0015 label_simTS(pset.getParameter<edm::InputTag>("label_simTS")),
0016 label_simTSFromCP(pset.getParameter<edm::InputTag>("label_simTSFromCP")),
0017 associator_(pset.getUntrackedParameter<edm::InputTag>("associator")),
0018 associatorSim_(pset.getUntrackedParameter<edm::InputTag>("associatorSim")),
0019 SaveGeneralInfo_(pset.getUntrackedParameter<bool>("SaveGeneralInfo")),
0020 doCaloParticlePlots_(pset.getUntrackedParameter<bool>("doCaloParticlePlots")),
0021 doCaloParticleSelection_(pset.getUntrackedParameter<bool>("doCaloParticleSelection")),
0022 doSimClustersPlots_(pset.getUntrackedParameter<bool>("doSimClustersPlots")),
0023 label_SimClustersPlots_(pset.getParameter<edm::InputTag>("label_SimClusters")),
0024 label_SimClustersLevel_(pset.getParameter<edm::InputTag>("label_SimClustersLevel")),
0025 doLayerClustersPlots_(pset.getUntrackedParameter<bool>("doLayerClustersPlots")),
0026 label_layerClustersPlots_(pset.getParameter<edm::InputTag>("label_layerClusterPlots")),
0027 label_LCToCPLinking_(pset.getParameter<edm::InputTag>("label_LCToCPLinking")),
0028 doTrackstersPlots_(pset.getUntrackedParameter<bool>("doTrackstersPlots")),
0029 label_TS_(pset.getParameter<std::string>("label_TS")),
0030 label_TSToCPLinking_(pset.getParameter<std::string>("label_TSToCPLinking")),
0031 label_TSToSTSPR_(pset.getParameter<std::string>("label_TSToSTSPR")),
0032 label_clustersmask(pset.getParameter<std::vector<edm::InputTag>>("LayerClustersInputMask")),
0033 cummatbudinxo_(pset.getParameter<edm::FileInPath>("cummatbudinxo")) {
0034
0035 const edm::InputTag& label_cp_effic_tag = pset.getParameter<edm::InputTag>("label_cp_effic");
0036 const edm::InputTag& label_cp_fake_tag = pset.getParameter<edm::InputTag>("label_cp_fake");
0037
0038 label_cp_effic = consumes<std::vector<CaloParticle>>(label_cp_effic_tag);
0039 label_cp_fake = consumes<std::vector<CaloParticle>>(label_cp_fake_tag);
0040
0041 simVertices_ = consumes<std::vector<SimVertex>>(pset.getParameter<edm::InputTag>("simVertices"));
0042
0043 for (auto& itag : label_clustersmask) {
0044 clustersMaskTokens_.push_back(consumes<std::vector<float>>(itag));
0045 }
0046
0047 associatorMapSimtR = consumes<hgcal::SimToRecoCollectionWithSimClusters>(associatorSim_);
0048 associatorMapRtSim = consumes<hgcal::RecoToSimCollectionWithSimClusters>(associatorSim_);
0049
0050 simTrackstersMap_ = consumes<std::map<uint, std::vector<uint>>>(edm::InputTag("ticlSimTracksters"));
0051
0052 hitMap_ = consumes<std::unordered_map<DetId, const HGCRecHit*>>(edm::InputTag("hgcalRecHitMapProducer"));
0053
0054 simClusters_ = consumes<std::vector<SimCluster>>(pset.getParameter<edm::InputTag>("label_scl"));
0055
0056 layerclusters_ = consumes<reco::CaloClusterCollection>(label_lcl);
0057
0058 for (auto& itag : label_tst) {
0059 label_tstTokens.push_back(consumes<ticl::TracksterCollection>(itag));
0060 }
0061
0062 simTracksters_ = consumes<ticl::TracksterCollection>(label_simTS);
0063 simTracksters_fromCPs_ = consumes<ticl::TracksterCollection>(label_simTSFromCP);
0064
0065 associatorMapRtS = consumes<hgcal::RecoToSimCollection>(associator_);
0066 associatorMapStR = consumes<hgcal::SimToRecoCollection>(associator_);
0067
0068 cpSelector = CaloParticleSelector(pset.getParameter<double>("ptMinCP"),
0069 pset.getParameter<double>("ptMaxCP"),
0070 pset.getParameter<double>("minRapidityCP"),
0071 pset.getParameter<double>("maxRapidityCP"),
0072 pset.getParameter<double>("lipCP"),
0073 pset.getParameter<double>("tipCP"),
0074 pset.getParameter<int>("minHitCP"),
0075 pset.getParameter<int>("maxSimClustersCP"),
0076 pset.getParameter<bool>("signalOnlyCP"),
0077 pset.getParameter<bool>("intimeOnlyCP"),
0078 pset.getParameter<bool>("chargedOnlyCP"),
0079 pset.getParameter<bool>("stableOnlyCP"),
0080 pset.getParameter<bool>("notConvertedOnlyCP"),
0081 pset.getParameter<std::vector<int>>("pdgIdCP"));
0082
0083 tools_.reset(new hgcal::RecHitTools());
0084
0085 particles_to_monitor_ = pset.getParameter<std::vector<int>>("pdgIdCP");
0086 totallayers_to_monitor_ = pset.getParameter<int>("totallayers_to_monitor");
0087 thicknesses_to_monitor_ = pset.getParameter<std::vector<int>>("thicknesses_to_monitor");
0088
0089
0090 std::ifstream fmb(cummatbudinxo_.fullPath().c_str());
0091 double thelay = 0.;
0092 double mbg = 0.;
0093 for (unsigned ilayer = 1; ilayer <= totallayers_to_monitor_; ++ilayer) {
0094 fmb >> thelay >> mbg;
0095 cummatbudg.insert(std::pair<double, double>(thelay, mbg));
0096 }
0097
0098 fmb.close();
0099
0100 ParameterSet psetForHistoProducerAlgo = pset.getParameter<ParameterSet>("histoProducerAlgoBlock");
0101 histoProducerAlgo_ = std::make_unique<HGVHistoProducerAlgo>(psetForHistoProducerAlgo);
0102
0103 dirName_ = pset.getParameter<std::string>("dirName");
0104 }
0105
0106 HGCalValidator::~HGCalValidator() {}
0107
0108 void HGCalValidator::bookHistograms(DQMStore::IBooker& ibook,
0109 edm::Run const&,
0110 edm::EventSetup const& setup,
0111 Histograms& histograms) const {
0112 if (SaveGeneralInfo_) {
0113 ibook.cd();
0114 ibook.setCurrentFolder(dirName_ + "GeneralInfo");
0115 histoProducerAlgo_->bookInfo(ibook, histograms.histoProducerAlgo);
0116 }
0117
0118 if (doCaloParticlePlots_) {
0119 ibook.cd();
0120
0121 for (auto const particle : particles_to_monitor_) {
0122 ibook.setCurrentFolder(dirName_ + "SelectedCaloParticles/" + std::to_string(particle));
0123 histoProducerAlgo_->bookCaloParticleHistos(
0124 ibook, histograms.histoProducerAlgo, particle, totallayers_to_monitor_);
0125 }
0126 ibook.cd();
0127 ibook.setCurrentFolder(dirName_);
0128 }
0129
0130
0131 if (doSimClustersPlots_) {
0132 ibook.cd();
0133 ibook.setCurrentFolder(dirName_ + label_SimClustersPlots_.label() + "/" + label_SimClustersLevel_.label());
0134 histoProducerAlgo_->bookSimClusterHistos(
0135 ibook, histograms.histoProducerAlgo, totallayers_to_monitor_, thicknesses_to_monitor_);
0136
0137 for (unsigned int ws = 0; ws < label_clustersmask.size(); ws++) {
0138 ibook.cd();
0139 InputTag algo = label_clustersmask[ws];
0140 string dirName = dirName_ + label_SimClustersPlots_.label() + "/";
0141 if (!algo.process().empty())
0142 dirName += algo.process() + "_";
0143 LogDebug("HGCalValidator") << dirName << "\n";
0144 if (!algo.label().empty())
0145 dirName += algo.label() + "_";
0146 LogDebug("HGCalValidator") << dirName << "\n";
0147 if (!algo.instance().empty())
0148 dirName += algo.instance() + "_";
0149 LogDebug("HGCalValidator") << dirName << "\n";
0150
0151 if (!dirName.empty()) {
0152 dirName.resize(dirName.size() - 1);
0153 }
0154
0155 LogDebug("HGCalValidator") << dirName << "\n";
0156
0157 ibook.setCurrentFolder(dirName);
0158
0159 histoProducerAlgo_->bookSimClusterAssociationHistos(
0160 ibook, histograms.histoProducerAlgo, totallayers_to_monitor_, thicknesses_to_monitor_);
0161 }
0162 }
0163
0164
0165 if (doLayerClustersPlots_) {
0166 ibook.cd();
0167 ibook.setCurrentFolder(dirName_ + label_layerClustersPlots_.label() + "/ClusterLevel");
0168 histoProducerAlgo_->bookClusterHistos_ClusterLevel(ibook,
0169 histograms.histoProducerAlgo,
0170 totallayers_to_monitor_,
0171 thicknesses_to_monitor_,
0172 cummatbudinxo_.fullPath());
0173 ibook.cd();
0174 ibook.setCurrentFolder(dirName_ + label_layerClustersPlots_.label() + "/" + label_LCToCPLinking_.label());
0175 histoProducerAlgo_->bookClusterHistos_LCtoCP_association(
0176 ibook, histograms.histoProducerAlgo, totallayers_to_monitor_);
0177
0178 ibook.cd();
0179 ibook.setCurrentFolder(dirName_ + label_layerClustersPlots_.label() + "/CellLevel");
0180 histoProducerAlgo_->bookClusterHistos_CellLevel(
0181 ibook, histograms.histoProducerAlgo, totallayers_to_monitor_, thicknesses_to_monitor_);
0182 }
0183
0184
0185 for (unsigned int www = 0; www < label_tst.size(); www++) {
0186 ibook.cd();
0187 InputTag algo = label_tst[www];
0188 string dirName = dirName_;
0189 if (!algo.process().empty())
0190 dirName += algo.process() + "_";
0191 LogDebug("HGCalValidator") << dirName << "\n";
0192 if (!algo.label().empty())
0193 dirName += algo.label() + "_";
0194 LogDebug("HGCalValidator") << dirName << "\n";
0195 if (!algo.instance().empty())
0196 dirName += algo.instance() + "_";
0197 LogDebug("HGCalValidator") << dirName << "\n";
0198
0199 if (!dirName.empty()) {
0200 dirName.resize(dirName.size() - 1);
0201 }
0202
0203 LogDebug("HGCalValidator") << dirName << "\n";
0204
0205 ibook.setCurrentFolder(dirName);
0206
0207
0208 if (doTrackstersPlots_) {
0209
0210 ibook.setCurrentFolder(dirName + "/" + label_TS_);
0211 histoProducerAlgo_->bookTracksterHistos(ibook, histograms.histoProducerAlgo, totallayers_to_monitor_);
0212
0213 ibook.setCurrentFolder(dirName + "/" + label_TSToCPLinking_);
0214 histoProducerAlgo_->bookTracksterSTSHistos(
0215 ibook, histograms.histoProducerAlgo, HGVHistoProducerAlgo::validationType::Linking);
0216
0217 ibook.setCurrentFolder(dirName + "/" + label_TSToSTSPR_);
0218 histoProducerAlgo_->bookTracksterSTSHistos(
0219 ibook, histograms.histoProducerAlgo, HGVHistoProducerAlgo::validationType::PatternRecognition);
0220 }
0221 }
0222 }
0223
0224 void HGCalValidator::cpParametersAndSelection(const Histograms& histograms,
0225 std::vector<CaloParticle> const& cPeff,
0226 std::vector<SimVertex> const& simVertices,
0227 std::vector<size_t>& selected_cPeff,
0228 unsigned int layers,
0229 std::unordered_map<DetId, const HGCRecHit*> const& hitMap) const {
0230 selected_cPeff.reserve(cPeff.size());
0231
0232 size_t j = 0;
0233 for (auto const& caloParticle : cPeff) {
0234 int id = caloParticle.pdgId();
0235
0236 if (!doCaloParticleSelection_ || (doCaloParticleSelection_ && cpSelector(caloParticle, simVertices))) {
0237 selected_cPeff.push_back(j);
0238 if (doCaloParticlePlots_) {
0239 histoProducerAlgo_->fill_caloparticle_histos(
0240 histograms.histoProducerAlgo, id, caloParticle, simVertices, layers, hitMap);
0241 }
0242 }
0243 ++j;
0244 }
0245 }
0246
0247 void HGCalValidator::dqmAnalyze(const edm::Event& event,
0248 const edm::EventSetup& setup,
0249 const Histograms& histograms) const {
0250 using namespace reco;
0251
0252 LogDebug("HGCalValidator") << "\n===================================================="
0253 << "\n"
0254 << "Analyzing new event"
0255 << "\n"
0256 << "====================================================\n"
0257 << "\n";
0258
0259 edm::Handle<std::vector<SimVertex>> simVerticesHandle;
0260 event.getByToken(simVertices_, simVerticesHandle);
0261 std::vector<SimVertex> const& simVertices = *simVerticesHandle;
0262
0263 edm::Handle<std::vector<CaloParticle>> caloParticleHandle;
0264 event.getByToken(label_cp_effic, caloParticleHandle);
0265 std::vector<CaloParticle> const& caloParticles = *caloParticleHandle;
0266
0267 edm::Handle<ticl::TracksterCollection> simTracksterHandle;
0268 event.getByToken(simTracksters_, simTracksterHandle);
0269 ticl::TracksterCollection const& simTracksters = *simTracksterHandle;
0270
0271 edm::Handle<ticl::TracksterCollection> simTracksterFromCPHandle;
0272 event.getByToken(simTracksters_fromCPs_, simTracksterFromCPHandle);
0273 ticl::TracksterCollection const& simTrackstersFromCPs = *simTracksterFromCPHandle;
0274
0275 edm::Handle<std::map<uint, std::vector<uint>>> simTrackstersMapHandle;
0276 event.getByToken(simTrackstersMap_, simTrackstersMapHandle);
0277 const std::map<uint, std::vector<uint>> cpToSc_SimTrackstersMap = *simTrackstersMapHandle;
0278
0279 edm::ESHandle<CaloGeometry> geom = setup.getHandle(caloGeomToken_);
0280 tools_->setGeometry(*geom);
0281 histoProducerAlgo_->setRecHitTools(tools_);
0282
0283 edm::Handle<hgcal::SimToRecoCollection> simtorecoCollectionH;
0284 event.getByToken(associatorMapStR, simtorecoCollectionH);
0285 auto simRecColl = *simtorecoCollectionH;
0286 edm::Handle<hgcal::RecoToSimCollection> recotosimCollectionH;
0287 event.getByToken(associatorMapRtS, recotosimCollectionH);
0288 auto recSimColl = *recotosimCollectionH;
0289
0290 edm::Handle<std::unordered_map<DetId, const HGCRecHit*>> hitMapHandle;
0291 event.getByToken(hitMap_, hitMapHandle);
0292 const std::unordered_map<DetId, const HGCRecHit*>* hitMap = &*hitMapHandle;
0293
0294
0295 if (SaveGeneralInfo_) {
0296 histoProducerAlgo_->fill_info_histos(histograms.histoProducerAlgo, totallayers_to_monitor_);
0297 }
0298
0299 std::vector<size_t> cPIndices;
0300
0301
0302 removeCPFromPU(caloParticles, cPIndices);
0303
0304
0305
0306
0307
0308 LogTrace("HGCalValidator") << "\n# of CaloParticles: " << caloParticles.size() << "\n" << std::endl;
0309 std::vector<size_t> selected_cPeff;
0310 cpParametersAndSelection(histograms, caloParticles, simVertices, selected_cPeff, totallayers_to_monitor_, *hitMap);
0311
0312
0313
0314 edm::Handle<std::vector<SimCluster>> simClustersHandle;
0315 event.getByToken(simClusters_, simClustersHandle);
0316 std::vector<SimCluster> const& simClusters = *simClustersHandle;
0317
0318
0319 edm::Handle<reco::CaloClusterCollection> clusterHandle;
0320 event.getByToken(layerclusters_, clusterHandle);
0321 const reco::CaloClusterCollection& clusters = *clusterHandle;
0322
0323 auto nSimClusters = simClusters.size();
0324 std::vector<size_t> sCIndices;
0325
0326
0327 for (unsigned int scId = 0; scId < nSimClusters; ++scId) {
0328 if (simClusters[scId].g4Tracks()[0].eventId().event() != 0 or
0329 simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() != 0) {
0330 LogDebug("HGCalValidator") << "Excluding SimClusters from event: "
0331 << simClusters[scId].g4Tracks()[0].eventId().event()
0332 << " with BX: " << simClusters[scId].g4Tracks()[0].eventId().bunchCrossing()
0333 << std::endl;
0334 continue;
0335 }
0336 sCIndices.emplace_back(scId);
0337 }
0338
0339
0340
0341
0342 if (doSimClustersPlots_) {
0343 histoProducerAlgo_->fill_simCluster_histos(
0344 histograms.histoProducerAlgo, simClusters, totallayers_to_monitor_, thicknesses_to_monitor_);
0345
0346 for (unsigned int ws = 0; ws < label_clustersmask.size(); ws++) {
0347 const auto& inputClusterMask = event.get(clustersMaskTokens_[ws]);
0348
0349 edm::Handle<hgcal::SimToRecoCollectionWithSimClusters> simtorecoCollectionH;
0350 event.getByToken(associatorMapSimtR, simtorecoCollectionH);
0351 auto simRecColl = *simtorecoCollectionH;
0352 edm::Handle<hgcal::RecoToSimCollectionWithSimClusters> recotosimCollectionH;
0353 event.getByToken(associatorMapRtSim, recotosimCollectionH);
0354 auto recSimColl = *recotosimCollectionH;
0355
0356 histoProducerAlgo_->fill_simClusterAssociation_histos(histograms.histoProducerAlgo,
0357 ws,
0358 clusterHandle,
0359 clusters,
0360 simClustersHandle,
0361 simClusters,
0362 sCIndices,
0363 inputClusterMask,
0364 *hitMap,
0365 totallayers_to_monitor_,
0366 recSimColl,
0367 simRecColl);
0368
0369
0370 LogTrace("HGCalValidator") << "\n# of SimClusters: " << nSimClusters
0371 << ", layerClusters mask label: " << label_clustersmask[ws].label() << "\n";
0372 }
0373 }
0374
0375
0376
0377
0378 int w = 0;
0379 if (doLayerClustersPlots_) {
0380 histoProducerAlgo_->fill_generic_cluster_histos(histograms.histoProducerAlgo,
0381 w,
0382 clusterHandle,
0383 clusters,
0384 caloParticleHandle,
0385 caloParticles,
0386 cPIndices,
0387 selected_cPeff,
0388 *hitMap,
0389 cummatbudg,
0390 totallayers_to_monitor_,
0391 thicknesses_to_monitor_,
0392 recSimColl,
0393 simRecColl);
0394
0395 for (unsigned int layerclusterIndex = 0; layerclusterIndex < clusters.size(); layerclusterIndex++) {
0396 histoProducerAlgo_->fill_cluster_histos(histograms.histoProducerAlgo, w, clusters[layerclusterIndex]);
0397 }
0398
0399
0400 LogTrace("HGCalValidator") << "\n# of layer clusters with " << label_lcl.process() << ":" << label_lcl.label()
0401 << ":" << label_lcl.instance() << ": " << clusters.size() << "\n";
0402 }
0403
0404
0405
0406
0407 for (unsigned int wml = 0; wml < label_tstTokens.size(); wml++) {
0408 if (doTrackstersPlots_) {
0409 edm::Handle<ticl::TracksterCollection> tracksterHandle;
0410 event.getByToken(label_tstTokens[wml], tracksterHandle);
0411 const ticl::TracksterCollection& tracksters = *tracksterHandle;
0412
0413
0414 LogTrace("HGCalValidator") << "\n# of Tracksters from " << label_tst[wml].process() << ":"
0415 << label_tst[wml].label() << ":" << label_tst[wml].instance() << ": "
0416 << tracksters.size() << "\n"
0417 << std::endl;
0418
0419 histoProducerAlgo_->fill_trackster_histos(histograms.histoProducerAlgo,
0420 wml,
0421 tracksters,
0422 clusters,
0423 simTracksters,
0424 simTrackstersFromCPs,
0425 cpToSc_SimTrackstersMap,
0426 simClusters,
0427 caloParticleHandle.id(),
0428 caloParticles,
0429 cPIndices,
0430 selected_cPeff,
0431 *hitMap,
0432 totallayers_to_monitor_);
0433 }
0434 }
0435 }