Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:14:49

0001 // -*- C++ -*-
0002 //
0003 // Package:    HcalScintillatorTester
0004 // Class:      HcalScintillatorTester
0005 //
0006 /**\class HcalScintillatorTester HcalScintillatorTester.cc test/HcalScintillatorTester.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Sunanda Banerjee
0015 //         Created:  Mon 2020/01/24
0016 //
0017 
0018 // system include files
0019 #include <memory>
0020 #include <iostream>
0021 #include <fstream>
0022 // user include files
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 // ------------ method called to produce the data  ------------
0055 void HcalScintillatorTester::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0056   const HcalDDDSimConstants hdc = iSetup.getData(token_);
0057 
0058   // Layers for HB
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   //Layers for HE
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 //define this as a plug-in
0181 DEFINE_FWK_MODULE(HcalScintillatorTester);