HGCalGeometry/test/HGCGeomAnalyzer.cc

0001 // user include files
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/StreamID.h"
0010
0011 // Geometry
0012 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0013 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0014 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0016 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0017
0018 #include <TH1F.h>
0019 #include <TH2F.h>
0020 #include <TProfile.h>
0021 #include <TProfile2D.h>
0022
0023 #include "DataFormats/Math/interface/Vector3D.h"
0024
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/PluginManager/interface/ModuleDef.h"
0027 #include "FWCore/ServiceRegistry/interface/Service.h"
0028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0029
0030 //STL headers
0031 #include <fstream>
0032 #include <iomanip>
0033 #include <map>
0034 #include <sstream>
0035 #include <string>
0036 #include <vector>
0037 #include "vdt/vdtMath.h"
0038
0039 //
0040 // class declaration
0041 //
0042 class HGCGeomAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0043 public:
0044   explicit HGCGeomAnalyzer(const edm::ParameterSet &);
0045   ~HGCGeomAnalyzer() override = default;
0046
0047 private:
0048   void beginJob() override {}
0049   void analyze(const edm::Event &, const edm::EventSetup &) override;
0050   void endJob() override {}
0051
0052   // ----------member data ---------------------------
0053   const std::string txtFileName_;
0054   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0055
0056   std::map<std::pair<DetId::Detector, int>, TProfile2D *> layerXYview_;
0057   std::map<std::pair<DetId::Detector, int>, TProfile *> layerThickR_;
0058   std::map<std::pair<DetId::Detector, int>, TProfile *> layerThickEta_;
0059 };
0060
0061 //
0062 HGCGeomAnalyzer::HGCGeomAnalyzer(const edm::ParameterSet &iConfig)
0063     : txtFileName_(iConfig.getParameter<std::string>("fileName")),
0064       geomToken_{esConsumes<CaloGeometry, CaloGeometryRecord>()} {
0065   usesResource("TFileService");
0066 }
0067
0068 //
0069 void HGCGeomAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &es) {
0070   //output text file
0071   std::ofstream boundaries;
0072   boundaries.open(txtFileName_);
0073
0074   //get geometry
0075   const CaloGeometry *caloGeom = &es.getData(geomToken_);
0076   edm::Service<TFileService> fs;
0077   fs->file().cd();
0078
0079   std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi, DetId::HGCalHSc};
0080   for (const auto &d : dets) {
0081     const HGCalGeometry *geom =
0082         static_cast<const HGCalGeometry *>(caloGeom->getSubdetectorGeometry(d, ForwardSubdetector::ForwardEmpty));
0083     const HGCalTopology *topo = &(geom->topology());
0084     const HGCalDDDConstants *ddd = &(topo->dddConstants());
0085
0086     //sub-detector boundaries
0087     unsigned int nlay = ddd->layers(true);
0088     unsigned int firstLay = ddd->firstLayer();
0089
0090     const std::vector<DetId> &detIdVec = geom->getValidDetIds();
0091     std::map<int, std::map<HGCSiliconDetId::waferType, std::map<double, double>>> typeRadMap;
0092
0093     //prepare header
0094     boundaries << "Subdetector: " << d << " has " << detIdVec.size() << " valid cells and " << nlay << " layers\n";
0095     boundaries << std::setw(5) << "layer" << std::setw(15) << "z" << std::setw(15) << "rIn" << std::setw(15) << "etaIn"
0096                << std::setw(15) << "rOut" << std::setw(15) << "etaOut" << std::setw(15);
0097
0098     if (d != DetId::HGCalHSc) {
0099       boundaries << "rInThin" << std::setw(15) << "etaInThin" << std::setw(15) << "rInThick" << std::setw(15)
0100                  << "etaInThick" << std::setw(15);
0101
0102       //book histograms
0103       TString baseName(Form("d%d_", d));
0104       TString title(d == DetId::HGCalEE ? "CEE" : "CEH_{Si}");
0105
0106       std::vector<int> sides = {-1, 1};
0107       for (const auto &zside : sides) {
0108         for (unsigned int ilay = firstLay; ilay < firstLay + nlay; ++ilay) {
0109           //layer and side histos
0110           int signedLayer = ilay * zside;
0111           std::pair<DetId::Detector, int> key(d, signedLayer);
0112           TString layerBaseName(Form("%slayer%d_", baseName.Data(), signedLayer));
0113           TString layerTitle(Form("%s %d", title.Data(), signedLayer));
0114           layerXYview_[key] = fs->make<TProfile2D>(
0115               layerBaseName + "xy_view", layerTitle + "; x [cm]; y [cm]; wafer type", 200, -200, 200, 200, -200, 200);
0116           layerThickR_[key] =
0117               fs->make<TProfile>(layerBaseName + "thickness_vs_r", layerTitle + "; r [cm]; wafer type", 200, 0, 200);
0118           layerThickEta_[key] = fs->make<TProfile>(
0119               layerBaseName + "thickness_vs_eta", layerTitle + "; abs(eta); wafer type", 200, 1.4, 3.3);
0120         }
0121       }
0122
0123       //for Si loop over the detIds to find the thickness transitions
0124       for (const auto &cellId : detIdVec) {
0125         HGCSiliconDetId id(cellId.rawId());
0126         GlobalPoint pt = geom->getPosition(id);
0127
0128         int layer = id.layer() * id.zside();
0129         double r(pt.perp());
0130         double eta(pt.eta());
0131         HGCSiliconDetId::waferType wt = (HGCSiliconDetId::waferType)id.type();
0132
0133         //fill histograms
0134         std::pair<DetId::Detector, int> key(d, layer);
0135         layerXYview_[key]->Fill(pt.x(), pt.y(), (int)wt + 1);
0136         layerThickR_[key]->Fill(r, (int)wt + 1);
0137         layerThickEta_[key]->Fill(fabs(eta), (int)wt + 1);
0138
0139         //fill map for boundary summary only for positive side
0140         if (pt.z() > 0)
0141           typeRadMap[layer][wt][r] = eta;
0142       }
0143     }
0144     boundaries << "\n";
0145
0146     //loop over map and print transitions
0147     for (unsigned int ilay = firstLay; ilay < firstLay + nlay; ++ilay) {
0148       double zz = ddd->waferZ(ilay, true);
0149       auto rr = ddd->rangeR(zz, true);
0150
0151       double rIn = rr.first;
0152       double rOut = rr.second;
0153
0154       math::XYZVector rInVec(0., rIn, zz);
0155       math::XYZVector rOutVec(0., rOut, zz);
0156       boundaries << std::setw(5) << ilay << std::setw(15) << rInVec.z() << std::setw(15) << rInVec.y() << std::setw(15)
0157                  << rInVec.eta() << std::setw(15) << rOutVec.y() << std::setw(15) << rOutVec.eta() << std::setw(15);
0158
0159       if (d != DetId::HGCalHSc) {
0160         boundaries << typeRadMap[ilay][HGCSiliconDetId::waferType::HGCalLD200].begin()->first << std::setw(15)
0161                    << typeRadMap[ilay][HGCSiliconDetId::waferType::HGCalLD200].begin()->second << std::setw(15)
0162                    << typeRadMap[ilay][HGCSiliconDetId::waferType::HGCalLD300].begin()->first << std::setw(15)
0163                    << typeRadMap[ilay][HGCSiliconDetId::waferType::HGCalLD300].begin()->second << std::setw(15);
0164       }
0165       boundaries << "\n";
0166     }
0167
0168     boundaries << "\n";
0169   }
0170 }
0171
0172 //define this as a plug-in
0173 DEFINE_FWK_MODULE(HGCGeomAnalyzer);