File indexing completed on 2024-09-11 04:33:23
0001
0002
0003 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0006
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014
0015 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0016 #include "DataFormats/DetId/interface/DetId.h"
0017 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0018 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0019 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0020 #include "DataFormats/Math/interface/deltaR.h"
0021 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0022
0023 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0024 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0025 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0026 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0027
0028 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0029 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0030
0031 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0032
0033 #include <map>
0034 #include <array>
0035 #include <string>
0036 #include <numeric>
0037 #include <cmath>
0038
0039 class HGCalHitCalibration : public DQMEDAnalyzer {
0040 public:
0041 explicit HGCalHitCalibration(const edm::ParameterSet&);
0042 ~HGCalHitCalibration() override;
0043
0044 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0045
0046 private:
0047 void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0048 void analyze(const edm::Event&, const edm::EventSetup&) override;
0049
0050 void fillWithRecHits(std::map<DetId, const HGCRecHit*>&, DetId, unsigned int, float, int&, float&);
0051
0052 edm::EDGetTokenT<HGCRecHitCollection> recHitsEE_;
0053 edm::EDGetTokenT<HGCRecHitCollection> recHitsFH_;
0054 edm::EDGetTokenT<HGCRecHitCollection> recHitsBH_;
0055 edm::EDGetTokenT<std::vector<CaloParticle> > caloParticles_;
0056 edm::EDGetTokenT<std::vector<reco::PFCluster> > hgcalMultiClusters_;
0057 edm::EDGetTokenT<std::vector<reco::GsfElectron> > electrons_;
0058 edm::EDGetTokenT<std::vector<reco::Photon> > photons_;
0059
0060 int algo_, depletion0_;
0061 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0062 bool rawRecHits_;
0063 int debug_;
0064 std::string folder_;
0065 hgcal::RecHitTools recHitTools_;
0066 static constexpr int depletion1_ = 200;
0067 static constexpr int depletion2_ = 300;
0068 static constexpr int scint_ = 400;
0069
0070 std::map<int, MonitorElement*> h_EoP_CPene_calib_fraction_;
0071 std::map<int, MonitorElement*> hgcal_EoP_CPene_calib_fraction_;
0072 std::map<int, MonitorElement*> hgcal_ele_EoP_CPene_calib_fraction_;
0073 std::map<int, MonitorElement*> hgcal_photon_EoP_CPene_calib_fraction_;
0074 MonitorElement* LayerOccupancy_;
0075
0076 static constexpr int layers_ = 60;
0077 std::array<float, layers_> Energy_layer_calib_;
0078 std::array<float, layers_> Energy_layer_calib_fraction_;
0079 };
0080
0081 HGCalHitCalibration::HGCalHitCalibration(const edm::ParameterSet& iConfig)
0082 : caloGeomToken_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0083 rawRecHits_(iConfig.getParameter<bool>("rawRecHits")),
0084 debug_(iConfig.getParameter<int>("debug")),
0085 folder_(iConfig.getParameter<std::string>("folder")) {
0086 auto detector = iConfig.getParameter<std::string>("detector");
0087 auto recHitsEE = iConfig.getParameter<edm::InputTag>("recHitsEE");
0088 auto recHitsFH = iConfig.getParameter<edm::InputTag>("recHitsFH");
0089 auto recHitsBH = iConfig.getParameter<edm::InputTag>("recHitsBH");
0090 auto caloParticles = iConfig.getParameter<edm::InputTag>("caloParticles");
0091 auto hgcalMultiClusters = iConfig.getParameter<edm::InputTag>("hgcalMultiClusters");
0092 auto electrons = iConfig.getParameter<edm::InputTag>("electrons");
0093 auto photons = iConfig.getParameter<edm::InputTag>("photons");
0094 if (detector == "all") {
0095 recHitsEE_ = consumes<HGCRecHitCollection>(recHitsEE);
0096 recHitsFH_ = consumes<HGCRecHitCollection>(recHitsFH);
0097 recHitsBH_ = consumes<HGCRecHitCollection>(recHitsBH);
0098 algo_ = 1;
0099 } else if (detector == "EM") {
0100 recHitsEE_ = consumes<HGCRecHitCollection>(recHitsEE);
0101 algo_ = 2;
0102 } else if (detector == "HAD") {
0103 recHitsFH_ = consumes<HGCRecHitCollection>(recHitsFH);
0104 recHitsBH_ = consumes<HGCRecHitCollection>(recHitsBH);
0105 algo_ = 3;
0106 }
0107 depletion0_ = iConfig.getParameter<int>("depletionFine");
0108 caloParticles_ = consumes<std::vector<CaloParticle> >(caloParticles);
0109 hgcalMultiClusters_ = consumes<std::vector<reco::PFCluster> >(hgcalMultiClusters);
0110 electrons_ = consumes<std::vector<reco::GsfElectron> >(electrons);
0111 photons_ = consumes<std::vector<reco::Photon> >(photons);
0112
0113
0114 Energy_layer_calib_.fill(0.);
0115 Energy_layer_calib_fraction_.fill(0.);
0116 }
0117
0118 HGCalHitCalibration::~HGCalHitCalibration() {
0119
0120
0121 }
0122
0123 void HGCalHitCalibration::bookHistograms(DQMStore::IBooker& ibooker,
0124 edm::Run const& iRun,
0125 edm::EventSetup const& ) {
0126 ibooker.cd();
0127 ibooker.setCurrentFolder(folder_);
0128 h_EoP_CPene_calib_fraction_[depletion0_] = ibooker.book1D("h_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0129 h_EoP_CPene_calib_fraction_[depletion1_] = ibooker.book1D("h_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0130 h_EoP_CPene_calib_fraction_[depletion2_] = ibooker.book1D("h_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0131 h_EoP_CPene_calib_fraction_[scint_] = ibooker.book1D("h_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0132 hgcal_EoP_CPene_calib_fraction_[depletion0_] =
0133 ibooker.book1D("hgcal_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0134 hgcal_EoP_CPene_calib_fraction_[depletion1_] =
0135 ibooker.book1D("hgcal_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0136 hgcal_EoP_CPene_calib_fraction_[depletion2_] =
0137 ibooker.book1D("hgcal_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0138 hgcal_EoP_CPene_calib_fraction_[scint_] = ibooker.book1D("hgcal_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0139 hgcal_ele_EoP_CPene_calib_fraction_[depletion0_] =
0140 ibooker.book1D("hgcal_ele_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0141 hgcal_ele_EoP_CPene_calib_fraction_[depletion1_] =
0142 ibooker.book1D("hgcal_ele_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0143 hgcal_ele_EoP_CPene_calib_fraction_[depletion2_] =
0144 ibooker.book1D("hgcal_ele_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0145 hgcal_ele_EoP_CPene_calib_fraction_[scint_] =
0146 ibooker.book1D("hgcal_ele_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0147 hgcal_photon_EoP_CPene_calib_fraction_[depletion0_] =
0148 ibooker.book1D("hgcal_photon_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0149 hgcal_photon_EoP_CPene_calib_fraction_[depletion1_] =
0150 ibooker.book1D("hgcal_photon_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0151 hgcal_photon_EoP_CPene_calib_fraction_[depletion2_] =
0152 ibooker.book1D("hgcal_photon_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0153 hgcal_photon_EoP_CPene_calib_fraction_[scint_] =
0154 ibooker.book1D("hgcal_photon_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0155 LayerOccupancy_ = ibooker.book1D("LayerOccupancy", "", layers_, 0., (float)layers_);
0156 }
0157
0158 void HGCalHitCalibration::fillWithRecHits(std::map<DetId, const HGCRecHit*>& hitmap,
0159 const DetId hitid,
0160 const unsigned int hitlayer,
0161 const float fraction,
0162 int& seedDet,
0163 float& seedEnergy) {
0164 if (hitmap.find(hitid) == hitmap.end()) {
0165
0166 IfLogTrace(debug_ > 0, "HGCalHitCalibration")
0167 << ">>> Failed to find detid " << std::hex << hitid.rawId() << std::dec << " Det " << hitid.det() << " Subdet "
0168 << hitid.subdetId() << std::endl;
0169 return;
0170 }
0171 if ((hitid.det() != DetId::Forward) && (hitid.det() != DetId::HGCalEE) && (hitid.det() != DetId::HGCalHSi) &&
0172 (hitid.det() != DetId::HGCalHSc)) {
0173 IfLogTrace(debug_ > 0, "HGCalHitCalibration")
0174 << ">>> Wrong type of detid " << std::hex << hitid.rawId() << std::dec << " Det " << hitid.det() << " Subdet "
0175 << hitid.subdetId() << std::endl;
0176 return;
0177 }
0178
0179 unsigned int layer = recHitTools_.getLayerWithOffset(hitid);
0180 assert(hitlayer == layer);
0181 Energy_layer_calib_fraction_[layer] += hitmap[hitid]->energy() * fraction;
0182 LayerOccupancy_->Fill(layer);
0183 if (seedEnergy < hitmap[hitid]->energy()) {
0184 seedEnergy = hitmap[hitid]->energy();
0185 seedDet = std::rint(recHitTools_.getSiThickness(hitid));
0186 if (hitid.det() == DetId::HGCalHSc) {
0187 seedDet = scint_;
0188 }
0189 }
0190 }
0191
0192 void HGCalHitCalibration::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0193 float constexpr max_dR2 = 0.0025;
0194
0195 const edm::ESHandle<CaloGeometry>& geom = iSetup.getHandle(caloGeomToken_);
0196 recHitTools_.setGeometry(*geom);
0197
0198 const edm::Handle<std::vector<CaloParticle> >& caloParticleHandle = iEvent.getHandle(caloParticles_);
0199 const std::vector<CaloParticle>& caloParticles = *caloParticleHandle;
0200
0201 const edm::Handle<std::vector<reco::PFCluster> >& hgcalMultiClustersHandle = iEvent.getHandle(hgcalMultiClusters_);
0202
0203 const edm::Handle<std::vector<reco::GsfElectron> >& PFElectronHandle = iEvent.getHandle(electrons_);
0204
0205 const edm::Handle<std::vector<reco::Photon> >& PFPhotonHandle = iEvent.getHandle(photons_);
0206
0207
0208 std::map<DetId, const HGCRecHit*> hitmap;
0209 switch (algo_) {
0210 case 1: {
0211 const auto& rechitsEE_handle = iEvent.getHandle(recHitsEE_);
0212 if (!rechitsEE_handle.isValid())
0213 return;
0214
0215 const auto& rechitsFH_handle = iEvent.getHandle(recHitsFH_);
0216 if (!rechitsFH_handle.isValid())
0217 return;
0218
0219 const auto& rechitsBH_handle = iEvent.getHandle(recHitsBH_);
0220 if (!rechitsBH_handle.isValid())
0221 return;
0222
0223 auto const& rechitsEE = *rechitsEE_handle;
0224 auto const& rechitsFH = *rechitsFH_handle;
0225 auto const& rechitsBH = *rechitsBH_handle;
0226 for (unsigned int i = 0; i < rechitsEE.size(); ++i) {
0227 hitmap[rechitsEE[i].detid()] = &rechitsEE[i];
0228 }
0229 for (unsigned int i = 0; i < rechitsFH.size(); ++i) {
0230 hitmap[rechitsFH[i].detid()] = &rechitsFH[i];
0231 }
0232 for (unsigned int i = 0; i < rechitsBH.size(); ++i) {
0233 hitmap[rechitsBH[i].detid()] = &rechitsBH[i];
0234 }
0235 break;
0236 }
0237 case 2: {
0238 const auto& rechitsEE_handle = iEvent.getHandle(recHitsEE_);
0239 if (!rechitsEE_handle.isValid())
0240 return;
0241
0242 auto const& rechitsEE = *rechitsEE_handle;
0243 for (unsigned int i = 0; i < rechitsEE.size(); i++) {
0244 hitmap[rechitsEE[i].detid()] = &rechitsEE[i];
0245 }
0246 break;
0247 }
0248 case 3: {
0249 const auto& rechitsFH_handle = iEvent.getHandle(recHitsFH_);
0250 if (!rechitsFH_handle.isValid())
0251 return;
0252
0253 const auto& rechitsBH_handle = iEvent.getHandle(recHitsBH_);
0254 if (!rechitsBH_handle.isValid())
0255 return;
0256
0257 auto const& rechitsFH = *rechitsFH_handle;
0258 auto const& rechitsBH = *rechitsBH_handle;
0259 for (unsigned int i = 0; i < rechitsFH.size(); i++) {
0260 hitmap[rechitsFH[i].detid()] = &rechitsFH[i];
0261 }
0262 for (unsigned int i = 0; i < rechitsBH.size(); i++) {
0263 hitmap[rechitsBH[i].detid()] = &rechitsBH[i];
0264 }
0265 break;
0266 }
0267 default:
0268 assert(false);
0269 break;
0270 }
0271
0272
0273 int seedDet = 0;
0274 float seedEnergy = 0.;
0275 IfLogTrace(debug_ > 0, "HGCalHitCalibration") << "Number of caloParticles: " << caloParticles.size() << std::endl;
0276 for (const auto& it_caloPart : caloParticles) {
0277 if (it_caloPart.g4Tracks()[0].eventId().event() != 0 or it_caloPart.g4Tracks()[0].eventId().bunchCrossing() != 0) {
0278 LogDebug("HGCalHitCalibration") << "Excluding CaloParticles from event: "
0279 << it_caloPart.g4Tracks()[0].eventId().event()
0280 << " with BX: " << it_caloPart.g4Tracks()[0].eventId().bunchCrossing()
0281 << std::endl;
0282 continue;
0283 }
0284
0285 const SimClusterRefVector& simClusterRefVector = it_caloPart.simClusters();
0286 Energy_layer_calib_.fill(0.);
0287 Energy_layer_calib_fraction_.fill(0.);
0288
0289 seedDet = 0;
0290 seedEnergy = 0.;
0291 for (const auto& it_sc : simClusterRefVector) {
0292 const SimCluster& simCluster = (*(it_sc));
0293 IfLogTrace(debug_ > 1, "HGCalHitCalibration")
0294 << ">>> SC.energy(): " << simCluster.energy() << " SC.simEnergy(): " << simCluster.simEnergy() << std::endl;
0295 const std::vector<std::pair<uint32_t, float> >& hits_and_fractions = simCluster.hits_and_fractions();
0296
0297
0298 for (const auto& it_haf : hits_and_fractions) {
0299 unsigned int hitlayer = recHitTools_.getLayerWithOffset(it_haf.first);
0300 DetId hitid = (it_haf.first);
0301
0302 if (rawRecHits_) {
0303 if ((hitid.det() == DetId::Forward) || (hitid.det() == DetId::HGCalEE) || (hitid.det() == DetId::HGCalHSi) ||
0304 (hitid.det() == DetId::HGCalHSc))
0305 fillWithRecHits(hitmap, hitid, hitlayer, it_haf.second, seedDet, seedEnergy);
0306 }
0307 }
0308 }
0309
0310 auto sumCalibRecHitCalib_fraction =
0311 std::accumulate(Energy_layer_calib_fraction_.begin(), Energy_layer_calib_fraction_.end(), 0.);
0312 IfLogTrace(debug_ > 0, "HGCalHitCalibration")
0313 << ">>> MC Energy: " << it_caloPart.energy() << " reco energy: " << sumCalibRecHitCalib_fraction << std::endl;
0314 if (h_EoP_CPene_calib_fraction_.count(seedDet))
0315 h_EoP_CPene_calib_fraction_[seedDet]->Fill(sumCalibRecHitCalib_fraction / it_caloPart.energy());
0316
0317
0318 if (hgcalMultiClustersHandle.isValid()) {
0319 const auto& clusters = *hgcalMultiClustersHandle;
0320 float total_energy = 0.;
0321 auto closest =
0322 std::min_element(clusters.begin(), clusters.end(), [&](const reco::PFCluster& a, const reco::PFCluster& b) {
0323 auto dR2_a = reco::deltaR2(it_caloPart, a);
0324 auto dR2_b = reco::deltaR2(it_caloPart, b);
0325 auto ERatio_a = a.correctedEnergy() / it_caloPart.energy();
0326 auto ERatio_b = b.correctedEnergy() / it_caloPart.energy();
0327
0328
0329 if (dR2_a < max_dR2 && dR2_b < max_dR2)
0330 return ERatio_a > ERatio_b;
0331 return dR2_a < dR2_b;
0332 });
0333 if (closest != clusters.end() && reco::deltaR2(*closest, it_caloPart) < 0.01) {
0334 total_energy = closest->correctedEnergy();
0335 seedDet = recHitTools_.getSiThickness(closest->seed());
0336 if (closest->seed().det() == DetId::HGCalHSc) {
0337 seedDet = scint_;
0338 }
0339 if (hgcal_EoP_CPene_calib_fraction_.count(seedDet)) {
0340 hgcal_EoP_CPene_calib_fraction_[seedDet]->Fill(total_energy / it_caloPart.energy());
0341 }
0342 }
0343 }
0344
0345 auto closest_fcn = [&](auto const& a, auto const& b) {
0346 auto dR2_a = reco::deltaR2(it_caloPart, a);
0347 auto dR2_b = reco::deltaR2(it_caloPart, b);
0348 auto ERatio_a = a.energy() / it_caloPart.energy();
0349 auto ERatio_b = b.energy() / it_caloPart.energy();
0350
0351
0352 if (dR2_a < max_dR2 && dR2_b < max_dR2)
0353 return ERatio_a > ERatio_b;
0354 return dR2_a < dR2_b;
0355 };
0356
0357
0358 if (PFElectronHandle.isValid()) {
0359 auto const& ele = (*PFElectronHandle);
0360 auto closest = std::min_element(ele.begin(), ele.end(), closest_fcn);
0361 if (closest != ele.end() &&
0362 ((closest->superCluster()->seed()->seed().det() == DetId::Forward) ||
0363 (closest->superCluster()->seed()->seed().det() == DetId::HGCalEE) ||
0364 (closest->superCluster()->seed()->seed().det() == DetId::HGCalHSi)) &&
0365 reco::deltaR2(*closest, it_caloPart) < 0.01) {
0366 seedDet = recHitTools_.getSiThickness(closest->superCluster()->seed()->seed());
0367 if (closest->superCluster()->seed()->seed().det() == DetId::HGCalHSc) {
0368 seedDet = scint_;
0369 }
0370 if (hgcal_ele_EoP_CPene_calib_fraction_.count(seedDet)) {
0371 hgcal_ele_EoP_CPene_calib_fraction_[seedDet]->Fill(closest->energy() / it_caloPart.energy());
0372 }
0373 }
0374 }
0375
0376
0377 if (PFPhotonHandle.isValid()) {
0378 auto const& photon = (*PFPhotonHandle);
0379 auto closest = std::min_element(photon.begin(), photon.end(), closest_fcn);
0380 if (closest != photon.end() &&
0381 ((closest->superCluster()->seed()->seed().det() == DetId::Forward) ||
0382 (closest->superCluster()->seed()->seed().det() == DetId::HGCalEE) ||
0383 (closest->superCluster()->seed()->seed().det() == DetId::HGCalHSi)) &&
0384 reco::deltaR2(*closest, it_caloPart) < 0.01) {
0385 seedDet = recHitTools_.getSiThickness(closest->superCluster()->seed()->seed());
0386 if (hgcal_photon_EoP_CPene_calib_fraction_.count(seedDet)) {
0387 hgcal_photon_EoP_CPene_calib_fraction_[seedDet]->Fill(closest->energy() / it_caloPart.energy());
0388 }
0389 }
0390 }
0391 }
0392 }
0393
0394
0395
0396 void HGCalHitCalibration::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0397 edm::ParameterSetDescription desc;
0398 desc.add<int>("debug", 0);
0399 desc.add<bool>("rawRecHits", true);
0400 desc.add<std::string>("folder", "HGCalHitCalibration");
0401 desc.add<std::string>("detector", "all");
0402 desc.add<int>("depletionFine", 120);
0403 desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0404 desc.add<edm::InputTag>("recHitsEE", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0405 desc.add<edm::InputTag>("recHitsFH", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0406 desc.add<edm::InputTag>("recHitsBH", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0407 desc.add<edm::InputTag>("hgcalMultiClusters", edm::InputTag("particleFlowClusterHGCal"));
0408 desc.add<edm::InputTag>("electrons", edm::InputTag("ecalDrivenGsfElectronsHGC"));
0409 desc.add<edm::InputTag>("photons", edm::InputTag("photonsHGC"));
0410 descriptions.add("hgcalHitCalibrationDefault", desc);
0411 }
0412
0413
0414 DEFINE_FWK_MODULE(HGCalHitCalibration);