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/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 #include "DQMServices/Core/interface/DQMStore.h"
0013
0014 #include "DataFormats/DetId/interface/DetId.h"
0015 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0016 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0017 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0018 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0019 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0020
0021 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0022 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0023 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0024 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0025
0026 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0028
0029 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0030
0031 #include <map>
0032 #include <array>
0033 #include <string>
0034 #include <numeric>
0035
0036 class HGCalShowerSeparation : public DQMEDAnalyzer {
0037 public:
0038 explicit HGCalShowerSeparation(const edm::ParameterSet&);
0039 ~HGCalShowerSeparation() override = default;
0040
0041 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0042
0043 private:
0044 void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0045 void analyze(const edm::Event&, const edm::EventSetup&) override;
0046
0047 void fillWithRecHits(std::unordered_map<DetId, const unsigned int>&, DetId, unsigned int, float, int&, float&);
0048
0049 edm::EDGetTokenT<std::unordered_map<DetId, const unsigned int>> hitMap_;
0050 edm::EDGetTokenT<std::vector<CaloParticle>> caloParticles_;
0051 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0052 std::vector<edm::InputTag> hits_label;
0053 std::vector<edm::EDGetTokenT<HGCRecHitCollection>> hits_token_;
0054
0055 int debug_;
0056 bool filterOnEnergyAndCaloP_;
0057 hgcal::RecHitTools recHitTools_;
0058 std::vector<HGCRecHit> hits_;
0059
0060 MonitorElement* eta1_;
0061 MonitorElement* eta2_;
0062 MonitorElement* energy1_;
0063 MonitorElement* energy2_;
0064 MonitorElement* energytot_;
0065 MonitorElement* scEnergy_;
0066 MonitorElement* showerProfile_;
0067 MonitorElement* layerEnergy_;
0068 MonitorElement* layerDistance_;
0069 MonitorElement* etaPhi_;
0070 MonitorElement* deltaEtaPhi_;
0071 std::vector<MonitorElement*> profileOnLayer_;
0072 std::vector<MonitorElement*> globalProfileOnLayer_;
0073 std::vector<MonitorElement*> distanceOnLayer_;
0074 std::vector<MonitorElement*> idealDistanceOnLayer_;
0075 std::vector<MonitorElement*> idealDeltaXY_;
0076 std::vector<MonitorElement*> centers_;
0077
0078 static constexpr int layers_ = 52;
0079 };
0080
0081 HGCalShowerSeparation::HGCalShowerSeparation(const edm::ParameterSet& iConfig)
0082 : tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0083 debug_(iConfig.getParameter<int>("debug")),
0084 filterOnEnergyAndCaloP_(iConfig.getParameter<bool>("filterOnEnergyAndCaloP")) {
0085 auto hitMapInputTag = iConfig.getParameter<edm::InputTag>("hitMapTag");
0086 auto caloParticles = iConfig.getParameter<edm::InputTag>("caloParticles");
0087 hitMap_ = consumes<std::unordered_map<DetId, const unsigned int>>(hitMapInputTag);
0088 caloParticles_ = consumes<std::vector<CaloParticle>>(caloParticles);
0089 auto hits_label_ = iConfig.getParameter<std::vector<edm::InputTag>>("hits");
0090
0091 for (auto& label : hits_label_) {
0092 hits_token_.push_back(consumes<HGCRecHitCollection>(label));
0093 }
0094 }
0095
0096 void HGCalShowerSeparation::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const&) {
0097 ibooker.cd();
0098 ibooker.setCurrentFolder("HGCalShowerSeparation");
0099 scEnergy_ = ibooker.book1D("SCEnergy", "SCEnergy", 240, 0., 120.);
0100 eta1_ = ibooker.book1D("eta1", "eta1", 80, 0., 4.);
0101 eta2_ = ibooker.book1D("eta2", "eta2", 80, 0., 4.);
0102 energy1_ = ibooker.book1D("energy1", "energy1", 240, 0., 120.);
0103 energy2_ = ibooker.book1D("energy2", "energy2", 240, 0., 120.);
0104 energytot_ = ibooker.book1D("energytot", "energytot", 200, 100., 200.);
0105 showerProfile_ = ibooker.book2D("ShowerProfile", "ShowerProfile", 800, -400., 400., layers_, 0., (float)layers_);
0106 layerEnergy_ = ibooker.book2D("LayerEnergy", "LayerEnergy", 60, 0., 60., 50, 0., 0.1);
0107 layerDistance_ = ibooker.book2D("LayerDistance", "LayerDistance", 60, 0., 60., 400, -400., 400.);
0108 etaPhi_ = ibooker.book2D("EtaPhi", "EtaPhi", 800, -4., 4., 800, -4., 4.);
0109 deltaEtaPhi_ = ibooker.book2D("DeltaEtaPhi", "DeltaEtaPhi", 100, -0.5, 0.5, 100, -0.5, 0.5);
0110 for (int i = 0; i < layers_; ++i) {
0111 profileOnLayer_.push_back(ibooker.book2D(std::string("ProfileOnLayer_") + std::to_string(i),
0112 std::string("ProfileOnLayer_") + std::to_string(i),
0113 120,
0114 -600.,
0115 600.,
0116 120,
0117 -600.,
0118 600.));
0119 globalProfileOnLayer_.push_back(ibooker.book2D(std::string("GlobalProfileOnLayer_") + std::to_string(i),
0120 std::string("GlobalProfileOnLayer_") + std::to_string(i),
0121 320,
0122 -160.,
0123 160.,
0124 320,
0125 -160.,
0126 160.));
0127 distanceOnLayer_.push_back(ibooker.book1D(std::string("DistanceOnLayer_") + std::to_string(i),
0128 std::string("DistanceOnLayer_") + std::to_string(i),
0129 120,
0130 -600.,
0131 600.));
0132 idealDistanceOnLayer_.push_back(ibooker.book1D(std::string("IdealDistanceOnLayer_") + std::to_string(i),
0133 std::string("IdealDistanceOnLayer_") + std::to_string(i),
0134 120,
0135 -600.,
0136 600.));
0137 idealDeltaXY_.push_back(ibooker.book2D(std::string("IdealDeltaXY_") + std::to_string(i),
0138 std::string("IdealDeltaXY_") + std::to_string(i),
0139 800,
0140 -400.,
0141 400.,
0142 800,
0143 -400.,
0144 400.));
0145 centers_.push_back(ibooker.book2D(std::string("Centers_") + std::to_string(i),
0146 std::string("Centers_") + std::to_string(i),
0147 320,
0148 -1600.,
0149 1600.,
0150 320,
0151 -1600.,
0152 1600.));
0153 }
0154 }
0155
0156 void HGCalShowerSeparation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0157 recHitTools_.setGeometry(iSetup.getData(tok_geom_));
0158
0159 for (auto& token : hits_token_) {
0160 edm::Handle<HGCRecHitCollection> hits_handle;
0161 iEvent.getByToken(token, hits_handle);
0162 hits_.insert(hits_.end(), (*hits_handle).begin(), (*hits_handle).end());
0163 }
0164
0165 const edm::Handle<std::vector<CaloParticle>>& caloParticleHandle = iEvent.getHandle(caloParticles_);
0166 const std::vector<CaloParticle>& caloParticles = *(caloParticleHandle.product());
0167
0168 edm::Handle<std::unordered_map<DetId, const unsigned int>> hitMapHandle;
0169 iEvent.getByToken(hitMap_, hitMapHandle);
0170 const auto hitmap = *hitMapHandle;
0171
0172
0173 IfLogTrace(debug_ > 0, "HGCalShowerSeparation") << "Number of caloParticles: " << caloParticles.size() << std::endl;
0174 if (caloParticles.size() == 2) {
0175 auto eta1 = caloParticles[0].eta();
0176 auto phi1 = caloParticles[0].phi();
0177 auto theta1 = 2. * std::atan(exp(-eta1));
0178 auto eta2 = caloParticles[1].eta();
0179 auto phi2 = caloParticles[1].phi();
0180 auto theta2 = 2. * std::atan(exp(-eta2));
0181 eta1_->Fill(eta1);
0182 eta2_->Fill(eta2);
0183
0184
0185
0186 int count = 0;
0187 int size = 0;
0188 float energy = 0.;
0189 float energy_tmp = 0.;
0190 for (const auto& it_caloPart : caloParticles) {
0191 count++;
0192 const SimClusterRefVector& simClusterRefVector = it_caloPart.simClusters();
0193 size += simClusterRefVector.size();
0194 for (const auto& it_sc : simClusterRefVector) {
0195 const SimCluster& simCluster = (*(it_sc));
0196 const std::vector<std::pair<uint32_t, float>>& hits_and_fractions = simCluster.hits_and_fractions();
0197 for (const auto& it_haf : hits_and_fractions) {
0198 if (hitmap.count(it_haf.first)) {
0199 const HGCRecHit* hit = &(hits_[hitmap.at(it_haf.first)]);
0200 energy += hit->energy() * it_haf.second;
0201 }
0202 }
0203 }
0204 if (count == 1) {
0205 energy1_->Fill(energy);
0206 energy_tmp = energy;
0207 } else {
0208 energy2_->Fill(energy - energy_tmp);
0209 }
0210 }
0211 energytot_->Fill(energy);
0212 if (filterOnEnergyAndCaloP_ && (energy < 2. * 0.8 * 80 or size != 2))
0213 return;
0214
0215 deltaEtaPhi_->Fill(eta1 - eta2, phi1 - phi2);
0216
0217 for (const auto& it_caloPart : caloParticles) {
0218 const SimClusterRefVector& simClusterRefVector = it_caloPart.simClusters();
0219 IfLogTrace(debug_ > 0, "HGCalShowerSeparation") << ">>> " << simClusterRefVector.size() << std::endl;
0220 for (const auto& it_sc : simClusterRefVector) {
0221 const SimCluster& simCluster = (*(it_sc));
0222 if (simCluster.energy() < 80 * 0.8)
0223 continue;
0224 scEnergy_->Fill(simCluster.energy());
0225 IfLogTrace(debug_ > 1, "HGCalShowerSeparation")
0226 << ">>> SC.energy(): " << simCluster.energy() << " SC.simEnergy(): " << simCluster.simEnergy() << std::endl;
0227 const std::vector<std::pair<uint32_t, float>>& hits_and_fractions = simCluster.hits_and_fractions();
0228
0229 for (const auto& it_haf : hits_and_fractions) {
0230 if (!hitmap.count(it_haf.first))
0231 continue;
0232 unsigned int hitlayer = recHitTools_.getLayerWithOffset(it_haf.first);
0233 auto global = recHitTools_.getPosition(it_haf.first);
0234 float globalx = global.x();
0235 float globaly = global.y();
0236 float globalz = global.z();
0237 if (globalz == 0)
0238 continue;
0239 auto rho1 = globalz * tan(theta1);
0240 auto rho2 = globalz * tan(theta2);
0241 auto x1 = rho1 * cos(phi1);
0242 auto y1 = rho1 * sin(phi1);
0243 auto x2 = rho2 * cos(phi2);
0244 auto y2 = rho2 * sin(phi2);
0245 auto half_point_x = (x1 + x2) / 2.;
0246 auto half_point_y = (y1 + y2) / 2.;
0247 auto half_point = sqrt((x1 - half_point_x) * (x1 - half_point_x) + (y1 - half_point_y) * (y1 - half_point_y));
0248 auto d_len = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
0249 auto dn_x = (x2 - x1) / d_len;
0250 auto dn_y = (y2 - y1) / d_len;
0251 auto distance = (globalx - x1) * dn_x + (globaly - y1) * dn_y;
0252 distance -= half_point;
0253 auto idealDistance = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
0254 if (hitmap.count(it_haf.first)) {
0255 const HGCRecHit* hit = &(hits_[hitmap.at(it_haf.first)]);
0256 profileOnLayer_[hitlayer]->Fill(
0257 10. * (globalx - half_point_x), 10. * (globaly - half_point_y), hit->energy() * it_haf.second);
0258 profileOnLayer_[55]->Fill(
0259 10. * (globalx - half_point_x), 10. * (globaly - half_point_y), hit->energy() * it_haf.second);
0260 globalProfileOnLayer_[hitlayer]->Fill(globalx, globaly, hit->energy() * it_haf.second);
0261 globalProfileOnLayer_[55]->Fill(globalx, globaly, hit->energy() * it_haf.second);
0262 layerEnergy_->Fill(hitlayer, hit->energy());
0263 layerDistance_->Fill(hitlayer, std::abs(10. * distance), hit->energy() * it_haf.second);
0264 etaPhi_->Fill(global.eta(), global.phi());
0265 distanceOnLayer_[hitlayer]->Fill(10. * distance);
0266 idealDistanceOnLayer_[hitlayer]->Fill(10. * idealDistance);
0267 idealDeltaXY_[hitlayer]->Fill(10. * (x1 - x2), 10. * (y1 - y2));
0268 centers_[hitlayer]->Fill(10. * half_point_x, 10. * half_point_y);
0269 IfLogTrace(debug_ > 0, "HGCalShowerSeparation")
0270 << ">>> " << distance << " " << hitlayer << " " << hit->energy() * it_haf.second << std::endl;
0271 showerProfile_->Fill(10. * distance, hitlayer, hit->energy() * it_haf.second);
0272 }
0273 }
0274 }
0275 }
0276 }
0277 }
0278
0279
0280
0281 void HGCalShowerSeparation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0282 edm::ParameterSetDescription desc;
0283 desc.add<int>("debug", 1);
0284 desc.add<bool>("filterOnEnergyAndCaloP", false);
0285 desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0286 desc.add<edm::InputTag>("hitMapTag", edm::InputTag("recHitMapProducer", "hgcalRecHitMap"));
0287 desc.add<std::vector<edm::InputTag>>("hits",
0288 {edm::InputTag("HGCalRecHit", "HGCEERecHits"),
0289 edm::InputTag("HGCalRecHit", "HGCHEFRecHits"),
0290 edm::InputTag("HGCalRecHit", "HGCHEBRecHits")});
0291 descriptions.add("hgcalShowerSeparationDefault", desc);
0292 }
0293
0294
0295 DEFINE_FWK_MODULE(HGCalShowerSeparation);
