Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-16 06:16:10

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   //In this way we can easily generalize to associations between other objects also.
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   density_ = consumes<Density>(edm::InputTag("hgcalLayerClusters"));
0055 
0056   simClusters_ = consumes<std::vector<SimCluster>>(pset.getParameter<edm::InputTag>("label_scl"));
0057 
0058   layerclusters_ = consumes<reco::CaloClusterCollection>(label_lcl);
0059 
0060   for (auto& itag : label_tst) {
0061     label_tstTokens.push_back(consumes<ticl::TracksterCollection>(itag));
0062   }
0063 
0064   simTracksters_ = consumes<ticl::TracksterCollection>(label_simTS);
0065   simTracksters_fromCPs_ = consumes<ticl::TracksterCollection>(label_simTSFromCP);
0066 
0067   associatorMapRtS = consumes<hgcal::RecoToSimCollection>(associator_);
0068   associatorMapStR = consumes<hgcal::SimToRecoCollection>(associator_);
0069 
0070   cpSelector = CaloParticleSelector(pset.getParameter<double>("ptMinCP"),
0071                                     pset.getParameter<double>("ptMaxCP"),
0072                                     pset.getParameter<double>("minRapidityCP"),
0073                                     pset.getParameter<double>("maxRapidityCP"),
0074                                     pset.getParameter<double>("lipCP"),
0075                                     pset.getParameter<double>("tipCP"),
0076                                     pset.getParameter<int>("minHitCP"),
0077                                     pset.getParameter<int>("maxSimClustersCP"),
0078                                     pset.getParameter<bool>("signalOnlyCP"),
0079                                     pset.getParameter<bool>("intimeOnlyCP"),
0080                                     pset.getParameter<bool>("chargedOnlyCP"),
0081                                     pset.getParameter<bool>("stableOnlyCP"),
0082                                     pset.getParameter<bool>("notConvertedOnlyCP"),
0083                                     pset.getParameter<std::vector<int>>("pdgIdCP"));
0084 
0085   tools_.reset(new hgcal::RecHitTools());
0086 
0087   particles_to_monitor_ = pset.getParameter<std::vector<int>>("pdgIdCP");
0088   totallayers_to_monitor_ = pset.getParameter<int>("totallayers_to_monitor");
0089   thicknesses_to_monitor_ = pset.getParameter<std::vector<int>>("thicknesses_to_monitor");
0090 
0091   //For the material budget file here
0092   std::ifstream fmb(cummatbudinxo_.fullPath().c_str());
0093   double thelay = 0.;
0094   double mbg = 0.;
0095   for (unsigned ilayer = 1; ilayer <= totallayers_to_monitor_; ++ilayer) {
0096     fmb >> thelay >> mbg;
0097     cummatbudg.insert(std::pair<double, double>(thelay, mbg));
0098   }
0099 
0100   fmb.close();
0101 
0102   ParameterSet psetForHistoProducerAlgo = pset.getParameter<ParameterSet>("histoProducerAlgoBlock");
0103   histoProducerAlgo_ = std::make_unique<HGVHistoProducerAlgo>(psetForHistoProducerAlgo);
0104 
0105   dirName_ = pset.getParameter<std::string>("dirName");
0106 }
0107 
0108 HGCalValidator::~HGCalValidator() {}
0109 
0110 void HGCalValidator::bookHistograms(DQMStore::IBooker& ibook,
0111                                     edm::Run const&,
0112                                     edm::EventSetup const& setup,
0113                                     Histograms& histograms) const {
0114   if (SaveGeneralInfo_) {
0115     ibook.cd();
0116     ibook.setCurrentFolder(dirName_ + "GeneralInfo");
0117     histoProducerAlgo_->bookInfo(ibook, histograms.histoProducerAlgo);
0118   }
0119 
0120   if (doCaloParticlePlots_) {
0121     ibook.cd();
0122 
0123     for (auto const particle : particles_to_monitor_) {
0124       ibook.setCurrentFolder(dirName_ + "SelectedCaloParticles/" + std::to_string(particle));
0125       histoProducerAlgo_->bookCaloParticleHistos(
0126           ibook, histograms.histoProducerAlgo, particle, totallayers_to_monitor_);
0127     }
0128     ibook.cd();
0129     ibook.setCurrentFolder(dirName_);
0130   }
0131 
0132   //Booking histograms concerning with simClusters
0133   if (doSimClustersPlots_) {
0134     ibook.cd();
0135     ibook.setCurrentFolder(dirName_ + label_SimClustersPlots_.label() + "/" + label_SimClustersLevel_.label());
0136     histoProducerAlgo_->bookSimClusterHistos(
0137         ibook, histograms.histoProducerAlgo, totallayers_to_monitor_, thicknesses_to_monitor_);
0138 
0139     for (unsigned int ws = 0; ws < label_clustersmask.size(); ws++) {
0140       ibook.cd();
0141       InputTag algo = label_clustersmask[ws];
0142       string dirName = dirName_ + label_SimClustersPlots_.label() + "/";
0143       if (!algo.process().empty())
0144         dirName += algo.process() + "_";
0145       LogDebug("HGCalValidator") << dirName << "\n";
0146       if (!algo.label().empty())
0147         dirName += algo.label() + "_";
0148       LogDebug("HGCalValidator") << dirName << "\n";
0149       if (!algo.instance().empty())
0150         dirName += algo.instance() + "_";
0151       LogDebug("HGCalValidator") << dirName << "\n";
0152 
0153       if (!dirName.empty()) {
0154         dirName.resize(dirName.size() - 1);
0155       }
0156 
0157       LogDebug("HGCalValidator") << dirName << "\n";
0158 
0159       ibook.setCurrentFolder(dirName);
0160 
0161       histoProducerAlgo_->bookSimClusterAssociationHistos(
0162           ibook, histograms.histoProducerAlgo, totallayers_to_monitor_, thicknesses_to_monitor_);
0163     }  //end of loop over masks
0164   }    //if for simCluster plots
0165 
0166   //Booking histograms concerning with hgcal layer clusters
0167   if (doLayerClustersPlots_) {
0168     ibook.cd();
0169     ibook.setCurrentFolder(dirName_ + label_layerClustersPlots_.label() + "/ClusterLevel");
0170     histoProducerAlgo_->bookClusterHistos_ClusterLevel(ibook,
0171                                                        histograms.histoProducerAlgo,
0172                                                        totallayers_to_monitor_,
0173                                                        thicknesses_to_monitor_,
0174                                                        cummatbudinxo_.fullPath());
0175     ibook.cd();
0176     ibook.setCurrentFolder(dirName_ + label_layerClustersPlots_.label() + "/" + label_LCToCPLinking_.label());
0177     histoProducerAlgo_->bookClusterHistos_LCtoCP_association(
0178         ibook, histograms.histoProducerAlgo, totallayers_to_monitor_);
0179 
0180     ibook.cd();
0181     ibook.setCurrentFolder(dirName_ + label_layerClustersPlots_.label() + "/CellLevel");
0182     histoProducerAlgo_->bookClusterHistos_CellLevel(
0183         ibook, histograms.histoProducerAlgo, totallayers_to_monitor_, thicknesses_to_monitor_);
0184   }
0185 
0186   //Booking histograms for Tracksters
0187   for (unsigned int www = 0; www < label_tst.size(); www++) {
0188     ibook.cd();
0189     InputTag algo = label_tst[www];
0190     string dirName = dirName_;
0191     if (!algo.process().empty())
0192       dirName += algo.process() + "_";
0193     LogDebug("HGCalValidator") << dirName << "\n";
0194     if (!algo.label().empty())
0195       dirName += algo.label() + "_";
0196     LogDebug("HGCalValidator") << dirName << "\n";
0197     if (!algo.instance().empty())
0198       dirName += algo.instance() + "_";
0199     LogDebug("HGCalValidator") << dirName << "\n";
0200 
0201     if (!dirName.empty()) {
0202       dirName.resize(dirName.size() - 1);
0203     }
0204 
0205     LogDebug("HGCalValidator") << dirName << "\n";
0206 
0207     ibook.setCurrentFolder(dirName);
0208 
0209     // Booking histograms concerning HGCal tracksters
0210     if (doTrackstersPlots_) {
0211       // Generic histos
0212       ibook.setCurrentFolder(dirName + "/" + label_TS_);
0213       histoProducerAlgo_->bookTracksterHistos(ibook, histograms.histoProducerAlgo, totallayers_to_monitor_);
0214       // CP Linking
0215       ibook.setCurrentFolder(dirName + "/" + label_TSToCPLinking_);
0216       histoProducerAlgo_->bookTracksterSTSHistos(
0217           ibook, histograms.histoProducerAlgo, HGVHistoProducerAlgo::validationType::Linking);
0218       // SimTracksters Pattern Recognition
0219       ibook.setCurrentFolder(dirName + "/" + label_TSToSTSPR_);
0220       histoProducerAlgo_->bookTracksterSTSHistos(
0221           ibook, histograms.histoProducerAlgo, HGVHistoProducerAlgo::validationType::PatternRecognition);
0222     }
0223   }  //end of booking Tracksters loop
0224 }
0225 
0226 void HGCalValidator::cpParametersAndSelection(const Histograms& histograms,
0227                                               std::vector<CaloParticle> const& cPeff,
0228                                               std::vector<SimVertex> const& simVertices,
0229                                               std::vector<size_t>& selected_cPeff,
0230                                               unsigned int layers,
0231                                               std::unordered_map<DetId, const HGCRecHit*> const& hitMap) const {
0232   selected_cPeff.reserve(cPeff.size());
0233 
0234   size_t j = 0;
0235   for (auto const& caloParticle : cPeff) {
0236     int id = caloParticle.pdgId();
0237 
0238     if (!doCaloParticleSelection_ || (doCaloParticleSelection_ && cpSelector(caloParticle, simVertices))) {
0239       selected_cPeff.push_back(j);
0240       if (doCaloParticlePlots_) {
0241         histoProducerAlgo_->fill_caloparticle_histos(
0242             histograms.histoProducerAlgo, id, caloParticle, simVertices, layers, hitMap);
0243       }
0244     }
0245     ++j;
0246   }  //end of loop over caloparticles
0247 }
0248 
0249 void HGCalValidator::dqmAnalyze(const edm::Event& event,
0250                                 const edm::EventSetup& setup,
0251                                 const Histograms& histograms) const {
0252   using namespace reco;
0253 
0254   LogDebug("HGCalValidator") << "\n===================================================="
0255                              << "\n"
0256                              << "Analyzing new event"
0257                              << "\n"
0258                              << "====================================================\n"
0259                              << "\n";
0260 
0261   edm::Handle<std::vector<SimVertex>> simVerticesHandle;
0262   event.getByToken(simVertices_, simVerticesHandle);
0263   std::vector<SimVertex> const& simVertices = *simVerticesHandle;
0264 
0265   edm::Handle<std::vector<CaloParticle>> caloParticleHandle;
0266   event.getByToken(label_cp_effic, caloParticleHandle);
0267   std::vector<CaloParticle> const& caloParticles = *caloParticleHandle;
0268 
0269   edm::Handle<ticl::TracksterCollection> simTracksterHandle;
0270   event.getByToken(simTracksters_, simTracksterHandle);
0271   ticl::TracksterCollection const& simTracksters = *simTracksterHandle;
0272 
0273   edm::Handle<ticl::TracksterCollection> simTracksterFromCPHandle;
0274   event.getByToken(simTracksters_fromCPs_, simTracksterFromCPHandle);
0275   ticl::TracksterCollection const& simTrackstersFromCPs = *simTracksterFromCPHandle;
0276 
0277   edm::Handle<std::map<uint, std::vector<uint>>> simTrackstersMapHandle;
0278   event.getByToken(simTrackstersMap_, simTrackstersMapHandle);
0279   const std::map<uint, std::vector<uint>> cpToSc_SimTrackstersMap = *simTrackstersMapHandle;
0280 
0281   edm::ESHandle<CaloGeometry> geom = setup.getHandle(caloGeomToken_);
0282   tools_->setGeometry(*geom);
0283   histoProducerAlgo_->setRecHitTools(tools_);
0284 
0285   edm::Handle<hgcal::SimToRecoCollection> simtorecoCollectionH;
0286   event.getByToken(associatorMapStR, simtorecoCollectionH);
0287   auto simRecColl = *simtorecoCollectionH;
0288   edm::Handle<hgcal::RecoToSimCollection> recotosimCollectionH;
0289   event.getByToken(associatorMapRtS, recotosimCollectionH);
0290   auto recSimColl = *recotosimCollectionH;
0291 
0292   edm::Handle<std::unordered_map<DetId, const HGCRecHit*>> hitMapHandle;
0293   event.getByToken(hitMap_, hitMapHandle);
0294   const std::unordered_map<DetId, const HGCRecHit*>* hitMap = &*hitMapHandle;
0295 
0296   //Some general info on layers etc.
0297   if (SaveGeneralInfo_) {
0298     histoProducerAlgo_->fill_info_histos(histograms.histoProducerAlgo, totallayers_to_monitor_);
0299   }
0300 
0301   std::vector<size_t> cPIndices;
0302   //Consider CaloParticles coming from the hard scatterer
0303   //excluding the PU contribution and save the indices.
0304   removeCPFromPU(caloParticles, cPIndices);
0305 
0306   // ##############################################
0307   // Fill caloparticles histograms
0308   // ##############################################
0309   // HGCRecHit are given to select the SimHits which are also reconstructed
0310   LogTrace("HGCalValidator") << "\n# of CaloParticles: " << caloParticles.size() << "\n" << std::endl;
0311   std::vector<size_t> selected_cPeff;
0312   cpParametersAndSelection(histograms, caloParticles, simVertices, selected_cPeff, totallayers_to_monitor_, *hitMap);
0313 
0314   //get collections from the event
0315   //simClusters
0316   edm::Handle<std::vector<SimCluster>> simClustersHandle;
0317   event.getByToken(simClusters_, simClustersHandle);
0318   std::vector<SimCluster> const& simClusters = *simClustersHandle;
0319 
0320   //Layer clusters
0321   edm::Handle<reco::CaloClusterCollection> clusterHandle;
0322   event.getByToken(layerclusters_, clusterHandle);
0323   const reco::CaloClusterCollection& clusters = *clusterHandle;
0324 
0325   //Density
0326   edm::Handle<Density> densityHandle;
0327   event.getByToken(density_, densityHandle);
0328   const Density& densities = *densityHandle;
0329 
0330   auto nSimClusters = simClusters.size();
0331   std::vector<size_t> sCIndices;
0332   //There shouldn't be any SimTracks from different crossings, but maybe they will be added later.
0333   //At the moment there should be one SimTrack in each SimCluster.
0334   for (unsigned int scId = 0; scId < nSimClusters; ++scId) {
0335     if (simClusters[scId].g4Tracks()[0].eventId().event() != 0 or
0336         simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() != 0) {
0337       LogDebug("HGCalValidator") << "Excluding SimClusters from event: "
0338                                  << simClusters[scId].g4Tracks()[0].eventId().event()
0339                                  << " with BX: " << simClusters[scId].g4Tracks()[0].eventId().bunchCrossing()
0340                                  << std::endl;
0341       continue;
0342     }
0343     sCIndices.emplace_back(scId);
0344   }
0345 
0346   // ##############################################
0347   // Fill simCluster histograms
0348   // ##############################################
0349   if (doSimClustersPlots_) {
0350     histoProducerAlgo_->fill_simCluster_histos(
0351         histograms.histoProducerAlgo, simClusters, totallayers_to_monitor_, thicknesses_to_monitor_);
0352 
0353     for (unsigned int ws = 0; ws < label_clustersmask.size(); ws++) {
0354       const auto& inputClusterMask = event.get(clustersMaskTokens_[ws]);
0355 
0356       edm::Handle<hgcal::SimToRecoCollectionWithSimClusters> simtorecoCollectionH;
0357       event.getByToken(associatorMapSimtR, simtorecoCollectionH);
0358       auto simRecColl = *simtorecoCollectionH;
0359       edm::Handle<hgcal::RecoToSimCollectionWithSimClusters> recotosimCollectionH;
0360       event.getByToken(associatorMapRtSim, recotosimCollectionH);
0361       auto recSimColl = *recotosimCollectionH;
0362 
0363       histoProducerAlgo_->fill_simClusterAssociation_histos(histograms.histoProducerAlgo,
0364                                                             ws,
0365                                                             clusterHandle,
0366                                                             clusters,
0367                                                             simClustersHandle,
0368                                                             simClusters,
0369                                                             sCIndices,
0370                                                             inputClusterMask,
0371                                                             *hitMap,
0372                                                             totallayers_to_monitor_,
0373                                                             recSimColl,
0374                                                             simRecColl);
0375 
0376       //General Info on simClusters
0377       LogTrace("HGCalValidator") << "\n# of SimClusters: " << nSimClusters
0378                                  << ", layerClusters mask label: " << label_clustersmask[ws].label() << "\n";
0379     }  //end of loop overs masks
0380   }
0381 
0382   // ##############################################
0383   // Fill layercluster histograms
0384   // ##############################################
0385   int w = 0;  //counter counting the number of sets of histograms
0386   if (doLayerClustersPlots_) {
0387     histoProducerAlgo_->fill_generic_cluster_histos(histograms.histoProducerAlgo,
0388                                                     w,
0389                                                     clusterHandle,
0390                                                     clusters,
0391                                                     densities,
0392                                                     caloParticleHandle,
0393                                                     caloParticles,
0394                                                     cPIndices,
0395                                                     selected_cPeff,
0396                                                     *hitMap,
0397                                                     cummatbudg,
0398                                                     totallayers_to_monitor_,
0399                                                     thicknesses_to_monitor_,
0400                                                     recSimColl,
0401                                                     simRecColl);
0402 
0403     for (unsigned int layerclusterIndex = 0; layerclusterIndex < clusters.size(); layerclusterIndex++) {
0404       histoProducerAlgo_->fill_cluster_histos(histograms.histoProducerAlgo, w, clusters[layerclusterIndex]);
0405     }
0406 
0407     //General Info on hgcalLayerClusters
0408     LogTrace("HGCalValidator") << "\n# of layer clusters with " << label_lcl.process() << ":" << label_lcl.label()
0409                                << ":" << label_lcl.instance() << ": " << clusters.size() << "\n";
0410   }
0411 
0412   // ##############################################
0413   // Fill Trackster histograms
0414   // ##############################################
0415   for (unsigned int wml = 0; wml < label_tstTokens.size(); wml++) {
0416     if (doTrackstersPlots_) {
0417       edm::Handle<ticl::TracksterCollection> tracksterHandle;
0418       event.getByToken(label_tstTokens[wml], tracksterHandle);
0419       const ticl::TracksterCollection& tracksters = *tracksterHandle;
0420 
0421       //General Info on Tracksters
0422       LogTrace("HGCalValidator") << "\n# of Tracksters from " << label_tst[wml].process() << ":"
0423                                  << label_tst[wml].label() << ":" << label_tst[wml].instance() << ": "
0424                                  << tracksters.size() << "\n"
0425                                  << std::endl;
0426 
0427       histoProducerAlgo_->fill_trackster_histos(histograms.histoProducerAlgo,
0428                                                 wml,
0429                                                 tracksters,
0430                                                 clusters,
0431                                                 simTracksters,
0432                                                 simTrackstersFromCPs,
0433                                                 cpToSc_SimTrackstersMap,
0434                                                 simClusters,
0435                                                 caloParticleHandle.id(),
0436                                                 caloParticles,
0437                                                 cPIndices,
0438                                                 selected_cPeff,
0439                                                 *hitMap,
0440                                                 totallayers_to_monitor_);
0441     }
0442   }  //end of loop over Trackster input labels
0443 }