File indexing completed on 2024-04-06 12:14:49
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <memory>
0020 #include <iostream>
0021 #include <fstream>
0022
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0032 #include "Geometry/Records/interface/HcalSimNumberingRecord.h"
0033 #include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
0034 #include "DataFormats/Math/interface/GeantUnits.h"
0035
0036 using namespace geant_units::operators;
0037
0038 class HcalScintillatorTester : public edm::one::EDAnalyzer<> {
0039 public:
0040 explicit HcalScintillatorTester(const edm::ParameterSet&);
0041 ~HcalScintillatorTester() override;
0042
0043 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0044
0045 private:
0046 edm::ESGetToken<HcalDDDSimConstants, HcalSimNumberingRecord> token_;
0047 };
0048
0049 HcalScintillatorTester::HcalScintillatorTester(const edm::ParameterSet&)
0050 : token_{esConsumes<HcalDDDSimConstants, HcalSimNumberingRecord>(edm::ESInputTag{})} {}
0051
0052 HcalScintillatorTester::~HcalScintillatorTester() {}
0053
0054
0055 void HcalScintillatorTester::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0056 const HcalDDDSimConstants hdc = iSetup.getData(token_);
0057
0058
0059 edm::LogVerbatim("HCalGeom") << "\n\nPrint out the parameters for HB"
0060 << "\n===============================\n";
0061 for (int zs = 0; zs <= 1; ++zs) {
0062 int zside = 2 * zs - 1;
0063 const auto etab = hdc.getiEtaRange(0);
0064 for (int eta = etab.first; eta <= etab.second; ++eta) {
0065 double etaL = hdc.parameter()->etaTable.at(eta - 1);
0066 double thetaL = 2. * atan(exp(-etaL));
0067 double etaH = hdc.parameter()->etaTable.at(eta);
0068 double thetaH = 2. * atan(exp(-etaH));
0069 for (int depth = hdc.getMinDepth(1, eta, 1, zside, false); depth <= hdc.getMaxDepth(1, eta, 1, zside, false);
0070 ++depth) {
0071 int layL = hdc.getLayerFront(1, eta, 1, zside, depth);
0072 int layH = hdc.getLayerBack(1, eta, 1, zside, depth);
0073 edm::LogVerbatim("HCalGeom") << "\nHB: zSide " << zside << " iEta " << eta << " Depth " << depth << " Layers "
0074 << layL << ":" << layH;
0075 for (int lay = (layL - 1); lay < layH; ++lay) {
0076 std::vector<double> area(2, 0), length(2, 0), width(2, 0);
0077 int kk(0);
0078 for (unsigned int k = 0; k < hdc.parameter()->layHB.size(); ++k) {
0079 if (lay == hdc.parameter()->layHB[k]) {
0080 double zmin = hdc.parameter()->rhoxHB[k] * std::cos(thetaL) / std::sin(thetaL);
0081 double zmax = hdc.parameter()->rhoxHB[k] * std::cos(thetaH) / std::sin(thetaH);
0082 double dz = (std::min(zmax, hdc.parameter()->dxHB[lay]) - zmin);
0083 if (dz > 0) {
0084 width[kk] = hdc.parameter()->dyHB[k];
0085 length[kk] = dz;
0086 area[kk] = dz * width[kk];
0087 ++kk;
0088 }
0089 }
0090 }
0091 if (kk > 1) {
0092 edm::LogVerbatim("HCalGeom") << "Layer " << (lay + 1) << " Length " << length[0] << ":" << length[1]
0093 << " width " << width[0] << ":" << width[1] << " Area " << area[0] << ":"
0094 << area[1];
0095 } else if (kk > 0) {
0096 edm::LogVerbatim("HCalGeom") << "Layer " << (lay + 1) << " Length " << length[0] << " width " << width[0]
0097 << " Area " << area[0];
0098 }
0099 }
0100 }
0101 }
0102 }
0103
0104
0105 edm::LogVerbatim("HCalGeom") << "\n\nPrint out the parameters for HE"
0106 << "\n===============================\n";
0107 for (int zs = 0; zs <= 1; ++zs) {
0108 int zside = 2 * zs - 1;
0109 const auto etae = hdc.getiEtaRange(1);
0110 for (int eta = etae.first; eta <= etae.second; ++eta) {
0111 double etaL = hdc.parameter()->etaTable.at(eta - 1);
0112 double thetaL = 2. * atan(exp(-etaL));
0113 double etaH = hdc.parameter()->etaTable.at(eta);
0114 double thetaH = 2. * atan(exp(-etaH));
0115 double phib = hdc.parameter()->phibin[eta - 1];
0116 int nphi = (phib > 6._deg) ? 1 : 2;
0117 for (int depth = hdc.getMinDepth(2, eta, 1, zside, false); depth <= hdc.getMaxDepth(2, eta, 1, zside, false);
0118 ++depth) {
0119 int layL = hdc.getLayerFront(2, eta, 1, zside, depth);
0120 int layH = hdc.getLayerBack(2, eta, 1, zside, depth);
0121 edm::LogVerbatim("HCalGeom") << "\nHE: zSide " << zside << " iEta " << eta << " Depth " << depth << " Layers "
0122 << layL << ":" << layH;
0123 for (int lay = (layL - 1); lay < layH; ++lay) {
0124 std::vector<double> area(4, 0), delr(4, 0), dely(4, 0);
0125 int kk(0);
0126 for (unsigned int k = 0; k < hdc.parameter()->layHE.size(); ++k) {
0127 if (lay == hdc.parameter()->layHE[k]) {
0128 double rmin = hdc.parameter()->zxHE[k] * std::tan(thetaH);
0129 double rmax = hdc.parameter()->zxHE[k] * std::tan(thetaL);
0130 if ((lay != 0 || eta == 18) &&
0131 (lay != 1 || (eta == 18 && (hdc.parameter()->rhoxHE[k] - hdc.parameter()->dyHE[k] > 1000)) ||
0132 (eta != 18 && (hdc.parameter()->rhoxHE[k] - hdc.parameter()->dyHE[k] < 1000))) &&
0133 ((rmin + 30) < (hdc.parameter()->rhoxHE[k] + hdc.parameter()->dyHE[k])) &&
0134 (rmax > (hdc.parameter()->rhoxHE[k] - hdc.parameter()->dyHE[k]))) {
0135 rmin = std::max(rmin, (hdc.parameter()->rhoxHE[k] - hdc.parameter()->dyHE[k]));
0136 rmax = std::min(rmax, (hdc.parameter()->rhoxHE[k] + hdc.parameter()->dyHE[k]));
0137 double dr = rmax - rmin;
0138 if (dr > 0) {
0139 double dx1 = rmin * std::tan(phib);
0140 double dx2 = rmax * std::tan(phib);
0141 if (nphi == 1) {
0142 double dy = 0.5 * (dx1 + dx2 - 4. * hdc.parameter()->dx1HE[k]);
0143 delr[kk] = dr;
0144 dely[kk] = dy;
0145 area[kk] = dr * dy;
0146 ++kk;
0147 } else {
0148 double dy1 = 0.5 * (dx1 + dx2 - 2. * hdc.parameter()->dx1HE[k]);
0149 delr[kk] = dr;
0150 dely[kk] = dy1;
0151 area[kk] = dr * dy1;
0152 double dy2 = 0.5 * ((rmax + rmin) * tan(10._deg) - 4 * hdc.parameter()->dx1HE[k]) - dy1;
0153 delr[kk + 2] = dr;
0154 dely[kk + 2] = dy2;
0155 area[kk + 2] = dr * dy2;
0156 ++kk;
0157 }
0158 }
0159 }
0160 }
0161 }
0162 if (area[0] > 0 && area[1] > 0) {
0163 if (nphi == 1) {
0164 edm::LogVerbatim("HCalGeom")
0165 << "Layer " << (lay + 1) << " Length " << delr[0] << ":" << delr[1] << " width " << dely[0] << ":"
0166 << dely[1] << " Area " << area[0] << ":" << area[1];
0167 } else {
0168 edm::LogVerbatim("HCalGeom")
0169 << "Layer " << (lay + 1) << " Length " << delr[0] << ":" << delr[1] << ":" << delr[2] << ":"
0170 << delr[3] << " width " << dely[0] << ":" << dely[1] << ":" << dely[2] << ":" << dely[3] << " Area "
0171 << area[0] << ":" << area[1] << ":" << area[2] << ":" << area[3];
0172 }
0173 }
0174 }
0175 }
0176 }
0177 }
0178 }
0179
0180
0181 DEFINE_FWK_MODULE(HcalScintillatorTester);