File indexing completed on 2024-07-02 00:53:46
0001
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
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
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
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
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
0071 std::ofstream boundaries;
0072 boundaries.open(txtFileName_);
0073
0074
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
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
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
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
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
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
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
0140 if (pt.z() > 0)
0141 typeRadMap[layer][wt][r] = eta;
0142 }
0143 }
0144 boundaries << "\n";
0145
0146
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
0173 DEFINE_FWK_MODULE(HGCGeomAnalyzer);