Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-06 01:54:03

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/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/ForwardDetId/interface/HGCalDetId.h"
0015 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0016 
0017 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0018 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0019 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0020 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0021 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0022 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0023 
0024 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0025 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0026 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0027 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0028 
0029 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0032 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0033 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0034 
0035 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0036 
0037 #include <map>
0038 #include <array>
0039 #include <string>
0040 #include <numeric>
0041 
0042 class HGCalShowerSeparation : public DQMEDAnalyzer {
0043 public:
0044   explicit HGCalShowerSeparation(const edm::ParameterSet&);
0045   ~HGCalShowerSeparation() override;
0046 
0047   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0048 
0049 private:
0050   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0051   void analyze(const edm::Event&, const edm::EventSetup&) override;
0052 
0053   void fillWithRecHits(std::unordered_map<DetId, const HGCRecHit*>&, DetId, unsigned int, float, int&, float&);
0054 
0055   edm::EDGetTokenT<std::unordered_map<DetId, const HGCRecHit*>> hitMap_;
0056   edm::EDGetTokenT<std::vector<CaloParticle>> caloParticles_;
0057   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0058 
0059   int debug_;
0060   bool filterOnEnergyAndCaloP_;
0061   hgcal::RecHitTools recHitTools_;
0062 
0063   MonitorElement* eta1_;
0064   MonitorElement* eta2_;
0065   MonitorElement* energy1_;
0066   MonitorElement* energy2_;
0067   MonitorElement* energytot_;
0068   MonitorElement* scEnergy_;
0069   MonitorElement* showerProfile_;
0070   MonitorElement* layerEnergy_;
0071   MonitorElement* layerDistance_;
0072   MonitorElement* etaPhi_;
0073   MonitorElement* deltaEtaPhi_;
0074   std::vector<MonitorElement*> profileOnLayer_;
0075   std::vector<MonitorElement*> globalProfileOnLayer_;
0076   std::vector<MonitorElement*> distanceOnLayer_;
0077   std::vector<MonitorElement*> idealDistanceOnLayer_;
0078   std::vector<MonitorElement*> idealDeltaXY_;
0079   std::vector<MonitorElement*> centers_;
0080 
0081   static constexpr int layers_ = 52;
0082 };
0083 
0084 HGCalShowerSeparation::HGCalShowerSeparation(const edm::ParameterSet& iConfig)
0085     : tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0086       debug_(iConfig.getParameter<int>("debug")),
0087       filterOnEnergyAndCaloP_(iConfig.getParameter<bool>("filterOnEnergyAndCaloP")) {
0088   auto hitMapInputTag = iConfig.getParameter<edm::InputTag>("hitMapTag");
0089   auto caloParticles = iConfig.getParameter<edm::InputTag>("caloParticles");
0090   hitMap_ = consumes<std::unordered_map<DetId, const HGCRecHit*>>(hitMapInputTag);
0091   caloParticles_ = consumes<std::vector<CaloParticle>>(caloParticles);
0092 }
0093 
0094 HGCalShowerSeparation::~HGCalShowerSeparation() {
0095   // do anything here that needs to be done at desctruction time
0096   // (e.g. close files, deallocate resources etc.)
0097 }
0098 
0099 void HGCalShowerSeparation::bookHistograms(DQMStore::IBooker& ibooker,
0100                                            edm::Run const& iRun,
0101                                            edm::EventSetup const& /* iSetup */) {
0102   ibooker.cd();
0103   ibooker.setCurrentFolder("HGCalShowerSeparation");
0104   scEnergy_ = ibooker.book1D("SCEnergy", "SCEnergy", 240, 0., 120.);
0105   eta1_ = ibooker.book1D("eta1", "eta1", 80, 0., 4.);
0106   eta2_ = ibooker.book1D("eta2", "eta2", 80, 0., 4.);
0107   energy1_ = ibooker.book1D("energy1", "energy1", 240, 0., 120.);
0108   energy2_ = ibooker.book1D("energy2", "energy2", 240, 0., 120.);
0109   energytot_ = ibooker.book1D("energytot", "energytot", 200, 100., 200.);
0110   showerProfile_ = ibooker.book2D("ShowerProfile", "ShowerProfile", 800, -400., 400., layers_, 0., (float)layers_);
0111   layerEnergy_ = ibooker.book2D("LayerEnergy", "LayerEnergy", 60, 0., 60., 50, 0., 0.1);
0112   layerDistance_ = ibooker.book2D("LayerDistance", "LayerDistance", 60, 0., 60., 400, -400., 400.);
0113   etaPhi_ = ibooker.book2D("EtaPhi", "EtaPhi", 800, -4., 4., 800, -4., 4.);
0114   deltaEtaPhi_ = ibooker.book2D("DeltaEtaPhi", "DeltaEtaPhi", 100, -0.5, 0.5, 100, -0.5, 0.5);
0115   for (int i = 0; i < layers_; ++i) {
0116     profileOnLayer_.push_back(ibooker.book2D(std::string("ProfileOnLayer_") + std::to_string(i),
0117                                              std::string("ProfileOnLayer_") + std::to_string(i),
0118                                              120,
0119                                              -600.,
0120                                              600.,
0121                                              120,
0122                                              -600.,
0123                                              600.));
0124     globalProfileOnLayer_.push_back(ibooker.book2D(std::string("GlobalProfileOnLayer_") + std::to_string(i),
0125                                                    std::string("GlobalProfileOnLayer_") + std::to_string(i),
0126                                                    320,
0127                                                    -160.,
0128                                                    160.,
0129                                                    320,
0130                                                    -160.,
0131                                                    160.));
0132     distanceOnLayer_.push_back(ibooker.book1D(std::string("DistanceOnLayer_") + std::to_string(i),
0133                                               std::string("DistanceOnLayer_") + std::to_string(i),
0134                                               120,
0135                                               -600.,
0136                                               600.));
0137     idealDistanceOnLayer_.push_back(ibooker.book1D(std::string("IdealDistanceOnLayer_") + std::to_string(i),
0138                                                    std::string("IdealDistanceOnLayer_") + std::to_string(i),
0139                                                    120,
0140                                                    -600.,
0141                                                    600.));
0142     idealDeltaXY_.push_back(ibooker.book2D(std::string("IdealDeltaXY_") + std::to_string(i),
0143                                            std::string("IdealDeltaXY_") + std::to_string(i),
0144                                            800,
0145                                            -400.,
0146                                            400.,
0147                                            800,
0148                                            -400.,
0149                                            400.));
0150     centers_.push_back(ibooker.book2D(std::string("Centers_") + std::to_string(i),
0151                                       std::string("Centers_") + std::to_string(i),
0152                                       320,
0153                                       -1600.,
0154                                       1600.,
0155                                       320,
0156                                       -1600.,
0157                                       1600.));
0158   }
0159 }
0160 
0161 void HGCalShowerSeparation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0162   using namespace edm;
0163 
0164   recHitTools_.setGeometry(iSetup.getData(tok_geom_));
0165 
0166   Handle<std::vector<CaloParticle>> caloParticleHandle;
0167   iEvent.getByToken(caloParticles_, caloParticleHandle);
0168   const std::vector<CaloParticle>& caloParticles = *caloParticleHandle;
0169 
0170   Handle<std::unordered_map<DetId, const HGCRecHit*>> hitMapHandle;
0171   iEvent.getByToken(hitMap_, hitMapHandle);
0172   const auto hitmap = *hitMapHandle;
0173 
0174   // loop over caloParticles
0175   IfLogTrace(debug_ > 0, "HGCalShowerSeparation") << "Number of caloParticles: " << caloParticles.size() << std::endl;
0176   if (caloParticles.size() == 2) {
0177     auto eta1 = caloParticles[0].eta();
0178     auto phi1 = caloParticles[0].phi();
0179     auto theta1 = 2. * atan(exp(-eta1));
0180     auto eta2 = caloParticles[1].eta();
0181     auto phi2 = caloParticles[1].phi();
0182     auto theta2 = 2. * atan(exp(-eta2));
0183     eta1_->Fill(eta1);
0184     eta2_->Fill(eta2);
0185 
0186     // Select event only if the sum of the energy of its recHits
0187     // is close enough to the gen energy
0188     int count = 0;
0189     int size = 0;
0190     float energy = 0.;
0191     float energy_tmp = 0.;
0192     for (const auto& it_caloPart : caloParticles) {
0193       count++;
0194       const SimClusterRefVector& simClusterRefVector = it_caloPart.simClusters();
0195       size += simClusterRefVector.size();
0196       for (const auto& it_sc : simClusterRefVector) {
0197         const SimCluster& simCluster = (*(it_sc));
0198         const std::vector<std::pair<uint32_t, float>>& hits_and_fractions = simCluster.hits_and_fractions();
0199         for (const auto& it_haf : hits_and_fractions) {
0200           if (hitmap.count(it_haf.first))
0201             energy += hitmap.at(it_haf.first)->energy() * it_haf.second;
0202         }  //hits and fractions
0203       }    // simcluster
0204       if (count == 1) {
0205         energy1_->Fill(energy);
0206         energy_tmp = energy;
0207       } else {
0208         energy2_->Fill(energy - energy_tmp);
0209       }
0210     }  // caloParticle
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             profileOnLayer_[hitlayer]->Fill(10. * (globalx - half_point_x),
0256                                             10. * (globaly - half_point_y),
0257                                             hitmap.at(it_haf.first)->energy() * it_haf.second);
0258             profileOnLayer_[55]->Fill(10. * (globalx - half_point_x),
0259                                       10. * (globaly - half_point_y),
0260                                       hitmap.at(it_haf.first)->energy() * it_haf.second);
0261             globalProfileOnLayer_[hitlayer]->Fill(globalx, globaly, hitmap.at(it_haf.first)->energy() * it_haf.second);
0262             globalProfileOnLayer_[55]->Fill(globalx, globaly, hitmap.at(it_haf.first)->energy() * it_haf.second);
0263             layerEnergy_->Fill(hitlayer, hitmap.at(it_haf.first)->energy());
0264             layerDistance_->Fill(hitlayer, std::abs(10. * distance), hitmap.at(it_haf.first)->energy() * it_haf.second);
0265             etaPhi_->Fill(global.eta(), global.phi());
0266             distanceOnLayer_[hitlayer]->Fill(10. * distance);                  //,
0267             idealDistanceOnLayer_[hitlayer]->Fill(10. * idealDistance);        //,
0268             idealDeltaXY_[hitlayer]->Fill(10. * (x1 - x2), 10. * (y1 - y2));   //,
0269             centers_[hitlayer]->Fill(10. * half_point_x, 10. * half_point_y);  //,
0270             IfLogTrace(debug_ > 0, "HGCalShowerSeparation")
0271                 << ">>> " << distance << " " << hitlayer << " " << hitmap.at(it_haf.first)->energy() * it_haf.second
0272                 << std::endl;
0273             showerProfile_->Fill(10. * distance, hitlayer, hitmap.at(it_haf.first)->energy() * it_haf.second);
0274           }
0275         }  // end simHit
0276       }    // end simCluster
0277     }      // end caloparticle
0278   }
0279 }
0280 
0281 // ------------ method fills 'descriptions' with the allowed parameters for the
0282 // module  ------------
0283 void HGCalShowerSeparation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0284   edm::ParameterSetDescription desc;
0285   desc.add<int>("debug", 1);
0286   desc.add<bool>("filterOnEnergyAndCaloP", false);
0287   desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0288   desc.add<edm::InputTag>("hitMapTag", edm::InputTag("hgcalRecHitMapProducer"));
0289   descriptions.add("hgcalShowerSeparationDefault", desc);
0290 }
0291 
0292 // define this as a plug-in
0293 DEFINE_FWK_MODULE(HGCalShowerSeparation);