Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:06:45 

0001 // system include files
0002 #include <cmath>
0003 #include <iostream>
0004 #include <fstream>
0005 #include <vector>
0006 #include <map>
0007 #include <string>
0008 
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010 
0011 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0012 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0013 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0014 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0015 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0016 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0017 
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "FWCore/Utilities/interface/transform.h"
0027 #include "FWCore/Utilities/interface/InputTag.h"
0028 
0029 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0030 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0031 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0032 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0033 
0034 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0035 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0036 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0037 
0038 #include "CLHEP/Geometry/Point3D.h"
0039 #include "CLHEP/Geometry/Transform3D.h"
0040 #include "CLHEP/Geometry/Vector3D.h"
0041 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0042 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0043 
0044 #include "TH2D.h"
0045 
0046 class HGCalWaferStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0047 public:
0048   explicit HGCalWaferStudy(const edm::ParameterSet&);
0049   ~HGCalWaferStudy() override = default;
0050   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0051 
0052 protected:
0053   void beginJob() override {}
0054   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0055   void analyze(edm::Event const&, edm::EventSetup const&) override;
0056   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0057 
0058 private:
0059   void analyzeHits(int, const std::string&, const std::vector<PCaloHit>&);
0060 
0061   // ----------member data ---------------------------
0062   const std::vector<std::string> nameDetectors_, caloHitSources_;
0063   const std::vector<edm::InputTag> digiSources_;
0064   const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> ddTokens_;
0065   const std::vector<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>> geomTokens_;
0066   const int verbosity_, nBinHit_, nBinDig_;
0067   const double xyMinHit_, xyMaxHit_;
0068   const double xyMinDig_, xyMaxDig_;
0069   const bool ifNose_;
0070   std::vector<int> layerMnSim_, layerMxSim_;
0071   std::vector<int> layerMnDig_, layerMxDig_;
0072   std::vector<int> layerSim_, layerDig_, layerFront_;
0073   std::vector<const HGCalDDDConstants*> hgcons_;
0074   std::vector<const HGCalGeometry*> hgeoms_;
0075   std::vector<edm::EDGetTokenT<edm::PCaloHitContainer>> tok_hits_;
0076   std::vector<edm::EDGetTokenT<HGCalDigiCollection>> tok_digi_;
0077 
0078   //histogram related stuff
0079   static const int nType = 2;
0080   std::vector<TH2D*> h_XYsi1_[nType], h_XYsi2_[nType];
0081   std::vector<TH2D*> h_XYdi1_[nType], h_XYdi2_[nType];
0082 };
0083 
0084 HGCalWaferStudy::HGCalWaferStudy(const edm::ParameterSet& iConfig)
0085     : nameDetectors_(iConfig.getParameter<std::vector<std::string>>("detectorNames")),
0086       caloHitSources_(iConfig.getParameter<std::vector<std::string>>("caloHitSources")),
0087       digiSources_(iConfig.getParameter<std::vector<edm::InputTag>>("digiSources")),
0088       ddTokens_{
0089           edm::vector_transform(nameDetectors_,
0090                                 [this](const std::string& name) {
0091                                   return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0092                                       edm::ESInputTag{"", name});
0093                                 })},
0094       geomTokens_{edm::vector_transform(
0095           nameDetectors_,
0096           [this](const std::string& name) {
0097             return esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0098           })},
0099       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0100       nBinHit_(iConfig.getUntrackedParameter<int>("nBinHit", 600)),
0101       nBinDig_(iConfig.getUntrackedParameter<int>("nBinDig", 600)),
0102       xyMinHit_(iConfig.getUntrackedParameter<double>("xyMinHit", -3000.0)),
0103       xyMaxHit_(iConfig.getUntrackedParameter<double>("xyMaxHit", 3000.0)),
0104       xyMinDig_(iConfig.getUntrackedParameter<double>("xyMinDig", -300.0)),
0105       xyMaxDig_(iConfig.getUntrackedParameter<double>("xyMaxDig", 300.0)),
0106       ifNose_(iConfig.getUntrackedParameter<bool>("ifNose", false)),
0107       layerMnSim_(iConfig.getUntrackedParameter<std::vector<int>>("layerMinSim")),
0108       layerMxSim_(iConfig.getUntrackedParameter<std::vector<int>>("layerMaxSim")),
0109       layerMnDig_(iConfig.getUntrackedParameter<std::vector<int>>("layerMinDig")),
0110       layerMxDig_(iConfig.getUntrackedParameter<std::vector<int>>("layerMaxDig")) {
0111   usesResource(TFileService::kSharedResource);
0112 
0113   for (auto const& source : caloHitSources_) {
0114     tok_hits_.emplace_back(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", source)));
0115     edm::LogVerbatim("HGCalValidation") << "SimHitSource: " << source;
0116   }
0117   for (auto const& source : digiSources_) {
0118     tok_digi_.emplace_back(consumes<HGCalDigiCollection>(source));
0119     edm::LogVerbatim("HGCalValidation") << "DigiSource: " << source;
0120   }
0121 }
0122 
0123 void HGCalWaferStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0124   edm::ParameterSetDescription desc;
0125   std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive"};
0126   std::vector<std::string> simSources = {"HGCHitsEE", "HGCHitsHEfront"};
0127   std::vector<edm::InputTag> digSources = {edm::InputTag("simHGCalUnsuppressedDigis", "EE"),
0128                                            edm::InputTag("simHGCalUnsuppressedDigis", "HEfront")};
0129   std::vector<int> layers = {28, 24};
0130   std::vector<int> layerMin = {1, 1};
0131   desc.add<std::vector<std::string>>("detectorNames", names);
0132   desc.add<std::vector<std::string>>("caloHitSources", simSources);
0133   desc.add<std::vector<edm::InputTag>>("digiSources", digSources);
0134   desc.addUntracked<int>("verbosity", 0);
0135   desc.addUntracked<int>("nBinHit", 600);
0136   desc.addUntracked<int>("nBinDig", 600);
0137   desc.addUntracked<double>("xyMinHit", -3000.0);
0138   desc.addUntracked<double>("xyMaxHit", 3000.0);
0139   desc.addUntracked<double>("xyMinDig", -300.0);
0140   desc.addUntracked<double>("xyMaxDig", 300.0);
0141   desc.addUntracked<bool>("ifNose", false);
0142   desc.addUntracked<std::vector<int>>("layerMaxSim", layers);
0143   desc.addUntracked<std::vector<int>>("layerMaxDig", layers);
0144   desc.addUntracked<std::vector<int>>("layerMinSim", layerMin);
0145   desc.addUntracked<std::vector<int>>("layerMinDig", layerMin);
0146   descriptions.add("hgcalWaferStudy", desc);
0147 }
0148 
0149 void HGCalWaferStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0150   //First the hits
0151   for (unsigned int k = 0; k < tok_hits_.size(); ++k) {
0152     const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hits_[k]);
0153     if (theCaloHitContainers.isValid()) {
0154       if (verbosity_ > 0)
0155         edm::LogVerbatim("HGCalValidation")
0156             << " PcalohitItr = " << theCaloHitContainers->size() << " Geometry Pointer " << hgcons_[k] << " # of Hists "
0157             << h_XYsi1_[k].size() << ":" << h_XYsi2_[k].size();
0158       int kount(0);
0159       for (auto const& hit : *(theCaloHitContainers)) {
0160         unsigned int id = hit.id();
0161         std::pair<float, float> xy;
0162         int layer(0), zside(1);
0163         bool wvtype(true);
0164         if (ifNose_) {
0165           HFNoseDetId detId = HFNoseDetId(id);
0166           layer = detId.layer();
0167           zside = detId.zside();
0168           wvtype = hgcons_[k]->waferVirtual(layer, detId.waferU(), detId.waferV());
0169           xy = hgcons_[k]->locateCell(
0170               zside, layer, detId.waferU(), detId.waferV(), detId.cellU(), detId.cellV(), false, true, false, false);
0171         } else if (hgcons_[k]->waferHexagon8()) {
0172           HGCSiliconDetId detId = HGCSiliconDetId(id);
0173           layer = detId.layer();
0174           zside = detId.zside();
0175           wvtype = hgcons_[k]->waferVirtual(layer, detId.waferU(), detId.waferV());
0176           xy = hgcons_[k]->locateCell(
0177               zside, layer, detId.waferU(), detId.waferV(), detId.cellU(), detId.cellV(), false, true, false, false);
0178         } else {
0179           int subdet, sector, type, cell;
0180           HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, type, cell);
0181           xy = hgcons_[k]->locateCell(cell, layer, sector, false);
0182           wvtype = hgcons_[k]->waferVirtual(layer, sector, 0);
0183         }
0184         double xp = (zside < 0) ? -xy.first : xy.first;
0185         double yp = xy.second;
0186         int ll = layer - layerMnSim_[k];
0187         if (verbosity_ > 1)
0188           edm::LogVerbatim("HGCalValidation")
0189               << "Hit [" << kount << "] Layer " << layer << ":" << ll << " x|y " << xp << ":" << yp << " Type "
0190               << wvtype << " Hist pointers " << h_XYsi1_[k][ll] << ":" << h_XYsi2_[k][ll];
0191         if (layer >= layerMnSim_[k] && layer <= layerMxSim_[k]) {
0192           if (wvtype)
0193             h_XYsi2_[k][ll]->Fill(xp, yp);
0194           else
0195             h_XYsi1_[k][ll]->Fill(xp, yp);
0196         }
0197         ++kount;
0198       }
0199     } else if (verbosity_ > 0) {
0200       edm::LogVerbatim("HGCalValidation") << "PCaloHitContainer does not "
0201                                           << "exist for " << nameDetectors_[k];
0202     }
0203   }
0204 
0205   //Then the digis
0206   for (unsigned int k = 0; k < tok_digi_.size(); ++k) {
0207     const edm::Handle<HGCalDigiCollection>& theHGCDigiContainer = iEvent.getHandle(tok_digi_[k]);
0208     if (theHGCDigiContainer.isValid()) {
0209       if (verbosity_ > 0)
0210         edm::LogVerbatim("HGCalValidation")
0211             << nameDetectors_[k] << " with " << theHGCDigiContainer->size() << " element(s) Geometry Pointer "
0212             << hgeoms_[k] << " # of Hists " << h_XYdi1_[k].size() << ":" << h_XYdi2_[k].size();
0213       int kount(0);
0214       for (const auto& it : *(theHGCDigiContainer.product())) {
0215         DetId id = it.id();
0216         int layer(0);
0217         bool wvtype(true);
0218         if (ifNose_) {
0219           HFNoseDetId detId = HFNoseDetId(id);
0220           layer = detId.layer();
0221           wvtype = hgcons_[k]->waferVirtual(layer, detId.waferU(), detId.waferV());
0222         } else if (hgcons_[k]->waferHexagon8()) {
0223           HGCSiliconDetId detId = HGCSiliconDetId(id);
0224           layer = detId.layer();
0225           wvtype = hgcons_[k]->waferVirtual(layer, detId.waferU(), detId.waferV());
0226         } else {
0227           HGCalDetId detId = HGCalDetId(id);
0228           layer = detId.layer();
0229           int wafer = detId.wafer();
0230           wvtype = hgcons_[k]->waferVirtual(layer, wafer, 0);
0231         }
0232         int ll = layer - layerMnDig_[k];
0233         const GlobalPoint& gcoord = hgeoms_[k]->getPosition(id);
0234         double xp = gcoord.x();
0235         double yp = gcoord.y();
0236         if (verbosity_ > 1)
0237           edm::LogVerbatim("HGCalValidation")
0238               << "Digi [" << kount << "] Layer " << layer << ":" << ll << " x|y " << xp << ":" << yp << " Type "
0239               << wvtype << " Hist pointers " << h_XYdi1_[k][ll] << ":" << h_XYdi2_[k][ll];
0240         if (layer >= layerMnDig_[k] && layer <= layerMxDig_[k]) {
0241           if (wvtype)
0242             h_XYdi2_[k][ll]->Fill(xp, yp);
0243           else
0244             h_XYdi1_[k][ll]->Fill(xp, yp);
0245         }
0246         ++kount;
0247       }
0248     } else {
0249       edm::LogVerbatim("HGCalValidation")
0250           << "DigiCollection handle " << digiSources_[k] << " does not exist for " << nameDetectors_[k] << " !!!";
0251     }
0252   }
0253 }
0254 
0255 // ------------ method called when starting to processes a run  ------------
0256 void HGCalWaferStudy::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0257   for (unsigned int k = 0; k < nameDetectors_.size(); ++k) {
0258     const auto& pHGDC = iSetup.getData(ddTokens_[k]);
0259     hgcons_.emplace_back(&pHGDC);
0260     layerSim_.emplace_back(hgcons_.back()->layers(false));
0261     layerDig_.emplace_back(hgcons_.back()->layers(true));
0262     layerFront_.emplace_back(hgcons_.back()->firstLayer());
0263     layerMnSim_[k] = std::max(layerMnSim_[k], layerFront_[k]);
0264     layerMnDig_[k] = std::max(layerMnDig_[k], layerFront_[k]);
0265     layerMxSim_[k] = std::min((layerFront_[k] + layerSim_[k] - 1), layerMxSim_[k]);
0266     layerMxDig_[k] = std::min((layerFront_[k] + layerDig_[k] - 1), layerMxDig_[k]);
0267 
0268     const auto& geom = iSetup.getData(geomTokens_[k]);
0269     hgeoms_.emplace_back(&geom);
0270     if (verbosity_ > 0)
0271       edm::LogVerbatim("HGCalValidation")
0272           << nameDetectors_[k] << " defined with " << layerFront_[k] << ":" << layerSim_[k] << ":" << layerDig_[k]
0273           << " Layers which gives limits " << layerMnSim_[k] << ":" << layerMxSim_[k] << " for Sim " << layerMnDig_[k]
0274           << ":" << layerMxDig_[k] << " for Digi "
0275           << "\nLimits for plots " << nBinHit_ << ":" << xyMinHit_ << ":" << xyMaxHit_ << " (Sim) " << nBinDig_ << ":"
0276           << xyMinDig_ << ":" << xyMaxDig_ << " (Digi) with Pointers " << hgcons_.back() << ":" << hgeoms_.back();
0277   }
0278 
0279   // Now define the histograms
0280   edm::Service<TFileService> fs;
0281   std::ostringstream name, title;
0282   for (unsigned int ih = 0; ih <= nameDetectors_.size(); ++ih) {
0283     for (int i = layerMnSim_[ih]; i <= layerMxSim_[ih]; ++i) {
0284       name.str("");
0285       title.str("");
0286       name << "XY_" << nameDetectors_[ih] << "L" << i << "SimR";
0287       title << "y vs x (Layer " << i << ") real wafers for " << nameDetectors_[ih] << " at SimHit level";
0288       if (verbosity_ > 0)
0289         edm::LogVerbatim("HGCalValidation") << i << " book " << name.str() << ":" << title.str();
0290       h_XYsi1_[ih].emplace_back(fs->make<TH2D>(
0291           name.str().c_str(), title.str().c_str(), nBinHit_, xyMinHit_, xyMaxHit_, nBinHit_, xyMinHit_, xyMaxHit_));
0292       name.str("");
0293       title.str("");
0294       name << "XY_" << nameDetectors_[ih] << "L" << i << "SimV";
0295       title << "y vs x (Layer " << i << ") virtual wafers for " << nameDetectors_[ih] << " at SimHit level";
0296       if (verbosity_ > 0)
0297         edm::LogVerbatim("HGCalValidation") << i << " book " << name.str() << ":" << title.str();
0298       h_XYsi2_[ih].emplace_back(fs->make<TH2D>(
0299           name.str().c_str(), title.str().c_str(), nBinHit_, xyMinHit_, xyMaxHit_, nBinHit_, xyMinHit_, xyMaxHit_));
0300     }
0301     edm::LogVerbatim("HGCalValidation") << "Complete booking of Sim Plots for " << nameDetectors_[ih];
0302 
0303     for (int i = layerMnDig_[ih]; i <= layerMxDig_[ih]; ++i) {
0304       name.str("");
0305       title.str("");
0306       name << "XY_" << nameDetectors_[ih] << "L" << i << "DigR";
0307       title << "y vs x (Layer " << i << ") real wafers for " << nameDetectors_[ih] << " at Digi level";
0308       if (verbosity_ > 0)
0309         edm::LogVerbatim("HGCalValidation") << i << " book " << name.str() << ":" << title.str();
0310       h_XYdi1_[ih].emplace_back(fs->make<TH2D>(
0311           name.str().c_str(), title.str().c_str(), nBinDig_, xyMinDig_, xyMaxDig_, nBinDig_, xyMinDig_, xyMaxDig_));
0312       name.str("");
0313       title.str("");
0314       name << "XY_" << nameDetectors_[ih] << "L" << i << "DigV";
0315       title << "y vs x (Layer " << i << ") virtual wafers for " << nameDetectors_[ih] << " at Digi level";
0316       if (verbosity_ > 0)
0317         edm::LogVerbatim("HGCalValidation") << i << " book " << name.str() << ":" << title.str();
0318       h_XYdi2_[ih].emplace_back(fs->make<TH2D>(
0319           name.str().c_str(), title.str().c_str(), nBinDig_, xyMinDig_, xyMaxDig_, nBinDig_, xyMinDig_, xyMaxDig_));
0320     }
0321     edm::LogVerbatim("HGCalValidation") << "Complete booking of Digi Plots for " << nameDetectors_[ih];
0322   }
0323 }
0324 
0325 #include "FWCore/Framework/interface/MakerMacros.h"
0326 //define this as a plug-in
0327 DEFINE_FWK_MODULE(HGCalWaferStudy);
[ 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!