Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:27:49 

0001 // user include files
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 
0038 class HGCalHitCalibration : public DQMEDAnalyzer {
0039 public:
0040   explicit HGCalHitCalibration(const edm::ParameterSet&);
0041   ~HGCalHitCalibration() override;
0042 
0043   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0044 
0045 private:
0046   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0047   void analyze(const edm::Event&, const edm::EventSetup&) override;
0048 
0049   void fillWithRecHits(std::map<DetId, const HGCRecHit*>&, DetId, unsigned int, float, int&, float&);
0050 
0051   edm::EDGetTokenT<HGCRecHitCollection> recHitsEE_;
0052   edm::EDGetTokenT<HGCRecHitCollection> recHitsFH_;
0053   edm::EDGetTokenT<HGCRecHitCollection> recHitsBH_;
0054   edm::EDGetTokenT<std::vector<CaloParticle> > caloParticles_;
0055   edm::EDGetTokenT<std::vector<reco::PFCluster> > hgcalMultiClusters_;
0056   edm::EDGetTokenT<std::vector<reco::GsfElectron> > electrons_;
0057   edm::EDGetTokenT<std::vector<reco::Photon> > photons_;
0058 
0059   int algo_, depletion0_;
0060   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0061   bool rawRecHits_;
0062   int debug_;
0063   hgcal::RecHitTools recHitTools_;
0064   static constexpr int depletion1_ = 200;
0065   static constexpr int depletion2_ = 300;
0066   static constexpr int scint_ = 400;
0067 
0068   std::map<int, MonitorElement*> h_EoP_CPene_calib_fraction_;
0069   std::map<int, MonitorElement*> hgcal_EoP_CPene_calib_fraction_;
0070   std::map<int, MonitorElement*> hgcal_ele_EoP_CPene_calib_fraction_;
0071   std::map<int, MonitorElement*> hgcal_photon_EoP_CPene_calib_fraction_;
0072   MonitorElement* LayerOccupancy_;
0073 
0074   static constexpr int layers_ = 60;
0075   std::array<float, layers_> Energy_layer_calib_;
0076   std::array<float, layers_> Energy_layer_calib_fraction_;
0077 };
0078 
0079 HGCalHitCalibration::HGCalHitCalibration(const edm::ParameterSet& iConfig)
0080     : caloGeomToken_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0081       rawRecHits_(iConfig.getParameter<bool>("rawRecHits")),
0082       debug_(iConfig.getParameter<int>("debug")) {
0083   auto detector = iConfig.getParameter<std::string>("detector");
0084   auto recHitsEE = iConfig.getParameter<edm::InputTag>("recHitsEE");
0085   auto recHitsFH = iConfig.getParameter<edm::InputTag>("recHitsFH");
0086   auto recHitsBH = iConfig.getParameter<edm::InputTag>("recHitsBH");
0087   auto caloParticles = iConfig.getParameter<edm::InputTag>("caloParticles");
0088   auto hgcalMultiClusters = iConfig.getParameter<edm::InputTag>("hgcalMultiClusters");
0089   auto electrons = iConfig.getParameter<edm::InputTag>("electrons");
0090   auto photons = iConfig.getParameter<edm::InputTag>("photons");
0091   if (detector == "all") {
0092     recHitsEE_ = consumes<HGCRecHitCollection>(recHitsEE);
0093     recHitsFH_ = consumes<HGCRecHitCollection>(recHitsFH);
0094     recHitsBH_ = consumes<HGCRecHitCollection>(recHitsBH);
0095     algo_ = 1;
0096   } else if (detector == "EM") {
0097     recHitsEE_ = consumes<HGCRecHitCollection>(recHitsEE);
0098     algo_ = 2;
0099   } else if (detector == "HAD") {
0100     recHitsFH_ = consumes<HGCRecHitCollection>(recHitsFH);
0101     recHitsBH_ = consumes<HGCRecHitCollection>(recHitsBH);
0102     algo_ = 3;
0103   }
0104   depletion0_ = iConfig.getParameter<int>("depletionFine");
0105   caloParticles_ = consumes<std::vector<CaloParticle> >(caloParticles);
0106   hgcalMultiClusters_ = consumes<std::vector<reco::PFCluster> >(hgcalMultiClusters);
0107   electrons_ = consumes<std::vector<reco::GsfElectron> >(electrons);
0108   photons_ = consumes<std::vector<reco::Photon> >(photons);
0109 
0110   // size should be HGC layers 52 is enough
0111   Energy_layer_calib_.fill(0.);
0112   Energy_layer_calib_fraction_.fill(0.);
0113 }
0114 
0115 HGCalHitCalibration::~HGCalHitCalibration() {
0116   // do anything here that needs to be done at desctruction time
0117   // (e.g. close files, deallocate resources etc.)
0118 }
0119 
0120 void HGCalHitCalibration::bookHistograms(DQMStore::IBooker& ibooker,
0121                                          edm::Run const& iRun,
0122                                          edm::EventSetup const& /* iSetup */) {
0123   ibooker.cd();
0124   ibooker.setCurrentFolder("HGCalHitCalibration");
0125   h_EoP_CPene_calib_fraction_[depletion0_] = ibooker.book1D("h_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0126   h_EoP_CPene_calib_fraction_[depletion1_] = ibooker.book1D("h_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0127   h_EoP_CPene_calib_fraction_[depletion2_] = ibooker.book1D("h_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0128   h_EoP_CPene_calib_fraction_[scint_] = ibooker.book1D("h_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0129   hgcal_EoP_CPene_calib_fraction_[depletion0_] =
0130       ibooker.book1D("hgcal_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0131   hgcal_EoP_CPene_calib_fraction_[depletion1_] =
0132       ibooker.book1D("hgcal_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0133   hgcal_EoP_CPene_calib_fraction_[depletion2_] =
0134       ibooker.book1D("hgcal_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0135   hgcal_EoP_CPene_calib_fraction_[scint_] = ibooker.book1D("hgcal_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0136   hgcal_ele_EoP_CPene_calib_fraction_[depletion0_] =
0137       ibooker.book1D("hgcal_ele_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0138   hgcal_ele_EoP_CPene_calib_fraction_[depletion1_] =
0139       ibooker.book1D("hgcal_ele_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0140   hgcal_ele_EoP_CPene_calib_fraction_[depletion2_] =
0141       ibooker.book1D("hgcal_ele_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0142   hgcal_ele_EoP_CPene_calib_fraction_[scint_] =
0143       ibooker.book1D("hgcal_ele_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0144   hgcal_photon_EoP_CPene_calib_fraction_[depletion0_] =
0145       ibooker.book1D("hgcal_photon_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
0146   hgcal_photon_EoP_CPene_calib_fraction_[depletion1_] =
0147       ibooker.book1D("hgcal_photon_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
0148   hgcal_photon_EoP_CPene_calib_fraction_[depletion2_] =
0149       ibooker.book1D("hgcal_photon_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
0150   hgcal_photon_EoP_CPene_calib_fraction_[scint_] =
0151       ibooker.book1D("hgcal_photon_EoP_CPene_scint_calib_fraction", "", 1000, -0.5, 2.5);
0152   LayerOccupancy_ = ibooker.book1D("LayerOccupancy", "", layers_, 0., (float)layers_);
0153 }
0154 
0155 void HGCalHitCalibration::fillWithRecHits(std::map<DetId, const HGCRecHit*>& hitmap,
0156                                           const DetId hitid,
0157                                           const unsigned int hitlayer,
0158                                           const float fraction,
0159                                           int& seedDet,
0160                                           float& seedEnergy) {
0161   if (hitmap.find(hitid) == hitmap.end()) {
0162     // Hit was not reconstructed
0163     IfLogTrace(debug_ > 0, "HGCalHitCalibration")
0164         << ">>> Failed to find detid " << std::hex << hitid.rawId() << std::dec << " Det " << hitid.det() << " Subdet "
0165         << hitid.subdetId() << std::endl;
0166     return;
0167   }
0168   if ((hitid.det() != DetId::Forward) && (hitid.det() != DetId::HGCalEE) && (hitid.det() != DetId::HGCalHSi) &&
0169       (hitid.det() != DetId::HGCalHSc)) {
0170     IfLogTrace(debug_ > 0, "HGCalHitCalibration")
0171         << ">>> Wrong type of detid " << std::hex << hitid.rawId() << std::dec << " Det " << hitid.det() << " Subdet "
0172         << hitid.subdetId() << std::endl;
0173     return;
0174   }
0175 
0176   unsigned int layer = recHitTools_.getLayerWithOffset(hitid);
0177   assert(hitlayer == layer);
0178   Energy_layer_calib_fraction_[layer] += hitmap[hitid]->energy() * fraction;
0179   LayerOccupancy_->Fill(layer);
0180   if (seedEnergy < hitmap[hitid]->energy()) {
0181     seedEnergy = hitmap[hitid]->energy();
0182     seedDet = recHitTools_.getSiThickness(hitid);
0183     if (hitid.det() == DetId::HGCalHSc) {
0184       seedDet = scint_;
0185     }
0186   }
0187 }
0188 
0189 void HGCalHitCalibration::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0190   const edm::ESHandle<CaloGeometry>& geom = iSetup.getHandle(caloGeomToken_);
0191   recHitTools_.setGeometry(*geom);
0192 
0193   const edm::Handle<std::vector<CaloParticle> >& caloParticleHandle = iEvent.getHandle(caloParticles_);
0194   const std::vector<CaloParticle>& caloParticles = *caloParticleHandle;
0195 
0196   const edm::Handle<std::vector<reco::PFCluster> >& hgcalMultiClustersHandle = iEvent.getHandle(hgcalMultiClusters_);
0197 
0198   const edm::Handle<std::vector<reco::GsfElectron> >& PFElectronHandle = iEvent.getHandle(electrons_);
0199 
0200   const edm::Handle<std::vector<reco::Photon> >& PFPhotonHandle = iEvent.getHandle(photons_);
0201 
0202   // make a map detid-rechit
0203   std::map<DetId, const HGCRecHit*> hitmap;
0204   switch (algo_) {
0205     case 1: {
0206       const auto& rechitsEE = iEvent.get(recHitsEE_);
0207       const auto& rechitsFH = iEvent.get(recHitsFH_);
0208       const auto& rechitsBH = iEvent.get(recHitsBH_);
0209       for (unsigned int i = 0; i < rechitsEE.size(); ++i) {
0210         hitmap[rechitsEE[i].detid()] = &rechitsEE[i];
0211       }
0212       for (unsigned int i = 0; i < rechitsFH.size(); ++i) {
0213         hitmap[rechitsFH[i].detid()] = &rechitsFH[i];
0214       }
0215       for (unsigned int i = 0; i < rechitsBH.size(); ++i) {
0216         hitmap[rechitsBH[i].detid()] = &rechitsBH[i];
0217       }
0218       break;
0219     }
0220     case 2: {
0221       const HGCRecHitCollection& rechitsEE = iEvent.get(recHitsEE_);
0222       for (unsigned int i = 0; i < rechitsEE.size(); i++) {
0223         hitmap[rechitsEE[i].detid()] = &rechitsEE[i];
0224       }
0225       break;
0226     }
0227     case 3: {
0228       const auto& rechitsFH = iEvent.get(recHitsFH_);
0229       const auto& rechitsBH = iEvent.get(recHitsBH_);
0230       for (unsigned int i = 0; i < rechitsFH.size(); i++) {
0231         hitmap[rechitsFH[i].detid()] = &rechitsFH[i];
0232       }
0233       for (unsigned int i = 0; i < rechitsBH.size(); i++) {
0234         hitmap[rechitsBH[i].detid()] = &rechitsBH[i];
0235       }
0236       break;
0237     }
0238     default:
0239       assert(false);
0240       break;
0241   }
0242 
0243   // loop over caloParticles
0244   int seedDet = 0;
0245   float seedEnergy = 0.;
0246   IfLogTrace(debug_ > 0, "HGCalHitCalibration") << "Number of caloParticles: " << caloParticles.size() << std::endl;
0247   for (const auto& it_caloPart : caloParticles) {
0248     if (it_caloPart.g4Tracks()[0].eventId().event() != 0 or it_caloPart.g4Tracks()[0].eventId().bunchCrossing() != 0) {
0249       LogDebug("HGCalHitCalibration") << "Excluding CaloParticles from event: "
0250                                       << it_caloPart.g4Tracks()[0].eventId().event()
0251                                       << " with BX: " << it_caloPart.g4Tracks()[0].eventId().bunchCrossing()
0252                                       << std::endl;
0253       continue;
0254     }
0255 
0256     const SimClusterRefVector& simClusterRefVector = it_caloPart.simClusters();
0257     Energy_layer_calib_.fill(0.);
0258     Energy_layer_calib_fraction_.fill(0.);
0259 
0260     seedDet = 0;
0261     seedEnergy = 0.;
0262     for (const auto& it_sc : simClusterRefVector) {
0263       const SimCluster& simCluster = (*(it_sc));
0264       IfLogTrace(debug_ > 1, "HGCalHitCalibration")
0265           << ">>> SC.energy(): " << simCluster.energy() << " SC.simEnergy(): " << simCluster.simEnergy() << std::endl;
0266       const std::vector<std::pair<uint32_t, float> >& hits_and_fractions = simCluster.hits_and_fractions();
0267 
0268       // loop over hits
0269       for (const auto& it_haf : hits_and_fractions) {
0270         unsigned int hitlayer = recHitTools_.getLayerWithOffset(it_haf.first);
0271         DetId hitid = (it_haf.first);
0272         // dump raw RecHits and match
0273         if (rawRecHits_) {
0274           if ((hitid.det() == DetId::Forward) || (hitid.det() == DetId::HGCalEE) || (hitid.det() == DetId::HGCalHSi) ||
0275               (hitid.det() == DetId::HGCalHSc))
0276             fillWithRecHits(hitmap, hitid, hitlayer, it_haf.second, seedDet, seedEnergy);
0277         }
0278       }  // end simHit
0279     }    // end simCluster
0280 
0281     auto sumCalibRecHitCalib_fraction =
0282         std::accumulate(Energy_layer_calib_fraction_.begin(), Energy_layer_calib_fraction_.end(), 0.);
0283     IfLogTrace(debug_ > 0, "HGCalHitCalibration")
0284         << ">>> MC Energy: " << it_caloPart.energy() << " reco energy: " << sumCalibRecHitCalib_fraction << std::endl;
0285     if (h_EoP_CPene_calib_fraction_.count(seedDet))
0286       h_EoP_CPene_calib_fraction_[seedDet]->Fill(sumCalibRecHitCalib_fraction / it_caloPart.energy());
0287 
0288     // Loop on reconstructed SC.
0289     const auto& clusters = *hgcalMultiClustersHandle;
0290     float total_energy = 0.;
0291     float max_dR2 = 0.0025;
0292     auto closest =
0293         std::min_element(clusters.begin(), clusters.end(), [&](const reco::PFCluster& a, const reco::PFCluster& b) {
0294           auto dR2_a = reco::deltaR2(it_caloPart, a);
0295           auto dR2_b = reco::deltaR2(it_caloPart, b);
0296           auto ERatio_a = a.correctedEnergy() / it_caloPart.energy();
0297           auto ERatio_b = b.correctedEnergy() / it_caloPart.energy();
0298           // If both clusters are within 0.0025, mark as first (min) the
0299           // element with the highest ratio against the SimCluster
0300           if (dR2_a < max_dR2 && dR2_b < max_dR2)
0301             return ERatio_a > ERatio_b;
0302           return dR2_a < dR2_b;
0303         });
0304     if (closest != clusters.end() && reco::deltaR2(*closest, it_caloPart) < 0.01) {
0305       total_energy = closest->correctedEnergy();
0306       seedDet = recHitTools_.getSiThickness(closest->seed());
0307       if (closest->seed().det() == DetId::HGCalHSc) {
0308         seedDet = scint_;
0309       }
0310       if (hgcal_EoP_CPene_calib_fraction_.count(seedDet)) {
0311         hgcal_EoP_CPene_calib_fraction_[seedDet]->Fill(total_energy / it_caloPart.energy());
0312       }
0313     }
0314 
0315     auto closest_fcn = [&](auto const& a, auto const& b) {
0316       auto dR2_a = reco::deltaR2(it_caloPart, a);
0317       auto dR2_b = reco::deltaR2(it_caloPart, b);
0318       auto ERatio_a = a.energy() / it_caloPart.energy();
0319       auto ERatio_b = b.energy() / it_caloPart.energy();
0320       // If both clusters are within 0.0025, mark as first (min) the
0321       // element with the highest ratio against the SimCluster
0322       if (dR2_a < max_dR2 && dR2_b < max_dR2)
0323         return ERatio_a > ERatio_b;
0324       return dR2_a < dR2_b;
0325     };
0326     // ELECTRONS in HGCAL
0327     if (PFElectronHandle.isValid()) {
0328       auto const& ele = (*PFElectronHandle);
0329       auto closest = std::min_element(ele.begin(), ele.end(), closest_fcn);
0330       if (closest != ele.end() &&
0331           ((closest->superCluster()->seed()->seed().det() == DetId::Forward) ||
0332            (closest->superCluster()->seed()->seed().det() == DetId::HGCalEE) ||
0333            (closest->superCluster()->seed()->seed().det() == DetId::HGCalHSi)) &&
0334           reco::deltaR2(*closest, it_caloPart) < 0.01) {
0335         seedDet = recHitTools_.getSiThickness(closest->superCluster()->seed()->seed());
0336         if (closest->superCluster()->seed()->seed().det() == DetId::HGCalHSc) {
0337           seedDet = scint_;
0338         }
0339         if (hgcal_ele_EoP_CPene_calib_fraction_.count(seedDet)) {
0340           hgcal_ele_EoP_CPene_calib_fraction_[seedDet]->Fill(closest->energy() / it_caloPart.energy());
0341         }
0342       }
0343     }
0344 
0345     // PHOTONS in HGCAL
0346     if (PFPhotonHandle.isValid()) {
0347       auto const& photon = (*PFPhotonHandle);
0348       auto closest = std::min_element(photon.begin(), photon.end(), closest_fcn);
0349       if (closest != photon.end() &&
0350           ((closest->superCluster()->seed()->seed().det() == DetId::Forward) ||
0351            (closest->superCluster()->seed()->seed().det() == DetId::HGCalEE) ||
0352            (closest->superCluster()->seed()->seed().det() == DetId::HGCalHSi)) &&
0353           reco::deltaR2(*closest, it_caloPart) < 0.01) {
0354         seedDet = recHitTools_.getSiThickness(closest->superCluster()->seed()->seed());
0355         if (hgcal_photon_EoP_CPene_calib_fraction_.count(seedDet)) {
0356           hgcal_photon_EoP_CPene_calib_fraction_[seedDet]->Fill(closest->energy() / it_caloPart.energy());
0357         }
0358       }
0359     }
0360   }  // end caloparticle
0361 }
0362 
0363 // ------------ method fills 'descriptions' with the allowed parameters for the
0364 // module  ------------
0365 void HGCalHitCalibration::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0366   edm::ParameterSetDescription desc;
0367   desc.add<int>("debug", 0);
0368   desc.add<bool>("rawRecHits", true);
0369   desc.add<std::string>("detector", "all");
0370   desc.add<int>("depletionFine", 120);
0371   desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0372   desc.add<edm::InputTag>("recHitsEE", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0373   desc.add<edm::InputTag>("recHitsFH", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0374   desc.add<edm::InputTag>("recHitsBH", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0375   desc.add<edm::InputTag>("hgcalMultiClusters", edm::InputTag("particleFlowClusterHGCal"));
0376   desc.add<edm::InputTag>("electrons", edm::InputTag("ecalDrivenGsfElectronsHGC"));
0377   desc.add<edm::InputTag>("photons", edm::InputTag("photonsHGC"));
0378   descriptions.add("hgcalHitCalibrationDefault", desc);
0379 }
0380 
0381 // define this as a plug-in
0382 DEFINE_FWK_MODULE(HGCalHitCalibration);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!