Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-17 22:53:47

0001 #include <iostream>
0002 #include <sstream>
0003 #include <string>
0004 #include <vector>
0005 
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 
0014 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0015 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0016 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0017 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0018 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0019 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0020 #include "CoralBase/Exception.h"
0021 
0022 class HGCalGeometryTester : public edm::one::EDAnalyzer<> {
0023 public:
0024   explicit HGCalGeometryTester(const edm::ParameterSet&);
0025   ~HGCalGeometryTester() override = default;
0026 
0027   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0028 
0029   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0030 
0031 private:
0032   void doTest(const HGCalGeometry* geom, ForwardSubdetector subdet);
0033   void doTestWafer(const HGCalGeometry* geom, DetId::Detector det);
0034   void doTestScint(const HGCalGeometry* geom, DetId::Detector det);
0035 
0036   const std::string name;
0037   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0038 };
0039 
0040 HGCalGeometryTester::HGCalGeometryTester(const edm::ParameterSet& iC)
0041     : name{iC.getParameter<std::string>("Detector")},
0042       geomToken_{esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", name})} {}
0043 
0044 void HGCalGeometryTester::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0045   const HGCalGeometry* geom = &(iSetup.getData(geomToken_));
0046   if (geom->topology().waferHexagon6()) {
0047     ForwardSubdetector subdet;
0048     if (name == "HGCalHESiliconSensitive")
0049       subdet = HGCHEF;
0050     else if (name == "HGCalHEScintillatorSensitive")
0051       subdet = HGCHEB;
0052     else
0053       subdet = HGCEE;
0054     edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name << " Detector:Subdetector " << DetId::Forward << ":"
0055                                    << subdet << " Mode " << geom->topology().dddConstants().geomMode();
0056     doTest(geom, subdet);
0057   } else {
0058     DetId::Detector det;
0059     if (name == "HGCalHESiliconSensitive")
0060       det = DetId::HGCalHSi;
0061     else if (name == "HGCalHEScintillatorSensitive")
0062       det = DetId::HGCalHSc;
0063     else
0064       det = DetId::HGCalEE;
0065     edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name << " Detector " << det << " Mode "
0066                                    << geom->topology().dddConstants().geomMode();
0067     if (name == "HGCalHEScintillatorSensitive") {
0068       doTestScint(geom, det);
0069     } else {
0070       doTestWafer(geom, det);
0071     }
0072   }
0073 }
0074 
0075 void HGCalGeometryTester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0076   edm::ParameterSetDescription desc;
0077   desc.add<std::string>("Detector", "HGCalEESensitive");
0078   descriptions.add("hgcalGeometryTesterEE", desc);
0079 }
0080 
0081 void HGCalGeometryTester::doTest(const HGCalGeometry* geom, ForwardSubdetector subdet) {
0082   const std::vector<DetId>& ids = geom->getValidDetIds();
0083   edm::LogVerbatim("HGCalGeomX") << "doTest: " << ids.size() << " valid ids for " << geom->cellElement();
0084 
0085   int layers[] = {1, 5, 10};
0086   int zsides[] = {1, -1};
0087   int cells[] = {1, 51, 101};
0088   int wafers[] = {1, 101, 201, 301, 401};
0089   const int ismax(5);
0090   for (int zside : zsides) {
0091     for (int is = 0; is < ismax; ++is) {
0092       int sector = wafers[is];
0093       int type = geom->topology().dddConstants().waferTypeT(sector);
0094       if (type != 1)
0095         type = 0;
0096       for (int layer : layers) {
0097         for (int cell : cells) {
0098           DetId id1;
0099           id1 = static_cast<DetId>(HGCalDetId(subdet, zside, layer, type, sector, cell));
0100           if (geom->topology().valid(id1)) {
0101             auto icell1 = geom->getGeometry(id1);
0102             GlobalPoint global1 = geom->getPosition(id1);
0103             DetId idc1 = geom->getClosestCell(global1);
0104             GlobalPoint global2 = geom->getPosition(idc1);
0105             std::string cherr = (id1.rawId() != idc1.rawId()) ? " ***** ERROR *****" : "";
0106             edm::LogVerbatim("HGCalGeomX")
0107                 << "DetId (" << subdet << ":" << zside << ":" << layer << ":" << sector << ":0:" << cell << ") Geom "
0108                 << icell1 << " position (" << global1.x() << ", " << global1.y() << ", " << global1.z() << ") ids "
0109                 << std::hex << id1.rawId() << ":" << idc1.rawId() << std::dec << ":" << HGCalDetId(id1) << ":"
0110                 << HGCalDetId(idc1) << " new position (" << global2.x() << ", " << global2.y() << ", " << global2.z()
0111                 << ") parameter[3] = " << icell1->param()[2] << ":" << icell1->param()[2] << cherr;
0112             std::vector<GlobalPoint> corners = geom->getCorners(idc1);
0113             std::ostringstream st1;
0114             st1 << corners.size() << " corners";
0115             for (auto const& cor : corners)
0116               st1 << " [" << cor.x() << "," << cor.y() << "," << cor.z() << "]";
0117             edm::LogVerbatim("HGCalGeomX") << st1.str();
0118           }
0119         }
0120       }
0121     }
0122   }
0123 }
0124 
0125 void HGCalGeometryTester::doTestWafer(const HGCalGeometry* geom, DetId::Detector det) {
0126   const std::vector<DetId>& ids = geom->getValidDetIds();
0127   edm::LogVerbatim("HGCalGeomX") << "doTestWafer:: " << ids.size() << " valid ids for " << geom->cellElement();
0128   int layers[] = {1, 5, 10};
0129   int zsides[] = {1, -1};
0130   int cells[] = {1, 4, 7};
0131   int wafers[] = {7, 5, 3, -3, -5, -7};
0132   for (int zside : zsides) {
0133     for (int layer : layers) {
0134       for (int waferU : wafers) {
0135         for (int waferV : wafers) {
0136           int type = geom->topology().dddConstants().getTypeHex(layer, waferU, waferV);
0137           edm::LogVerbatim("HGCalGeomX") << "zside " << zside << " layer " << layer << " wafer " << waferU << ":"
0138                                          << waferV << " type " << type;
0139           for (int cellU : cells) {
0140             for (int cellV : cells) {
0141               edm::LogVerbatim("HGCalGeomX") << "det " << det << " cell " << cellU << ":" << cellV;
0142               DetId id1 = static_cast<DetId>(HGCSiliconDetId(det, zside, type, layer, waferU, waferV, cellU, cellV));
0143               edm::LogVerbatim("HGCalGeomX") << HGCSiliconDetId(id1);
0144               if (geom->topology().valid(id1)) {
0145                 auto icell1 = geom->getGeometry(id1);
0146                 GlobalPoint global1 = geom->getPosition(id1);
0147                 DetId idc1 = geom->getClosestCell(global1);
0148                 GlobalPoint global2 = geom->getPosition(idc1);
0149                 std::string cherr = (id1.rawId() != idc1.rawId()) ? " ***** ERROR *****" : "";
0150                 edm::LogVerbatim("HGCalGeomX")
0151                     << "DetId (" << det << ":" << zside << ":" << type << ":" << layer << ":" << waferU << ":" << waferV
0152                     << ":" << cellU << ":" << cellV << ") Geom " << icell1 << " position (" << global1.x() << ", "
0153                     << global1.y() << ", " << global1.z() << ") ids " << std::hex << id1.rawId() << ":" << idc1.rawId()
0154                     << std::dec << ":" << HGCSiliconDetId(id1) << ":" << HGCSiliconDetId(idc1) << " new position ("
0155                     << global2.x() << ", " << global2.y() << ", " << global2.z()
0156                     << ") parameter[3] = " << icell1->param()[2] << ":" << icell1->param()[2] << cherr;
0157                 std::vector<GlobalPoint> corners = geom->getCorners(idc1);
0158                 std::ostringstream st1;
0159                 st1 << corners.size() << " corners";
0160                 for (auto const& cor : corners)
0161                   st1 << " [" << cor.x() << "," << cor.y() << "," << cor.z() << "]";
0162                 edm::LogVerbatim("HGCalGeomX") << st1.str();
0163               }
0164             }
0165           }
0166         }
0167       }
0168     }
0169   }
0170 }
0171 
0172 void HGCalGeometryTester::doTestScint(const HGCalGeometry* geom, DetId::Detector det) {
0173   const std::vector<DetId>& ids = geom->getValidDetIds();
0174   edm::LogVerbatim("HGCalGeomX") << "doTestScint: " << ids.size() << " valid ids for " << geom->cellElement();
0175   int layers[] = {9, 14, 21};
0176   int zsides[] = {1, -1};
0177   int iphis[] = {1, 51, 101, 151, 201};
0178   int ietas[] = {11, 20, 29};
0179   for (int zside : zsides) {
0180     for (int layer : layers) {
0181       int type = geom->topology().dddConstants().getTypeTrap(layer);
0182       for (int ieta : ietas) {
0183         std::pair<int, int> typm = geom->topology().dddConstants().tileType(layer, ieta, 0);
0184         for (int iphi : iphis) {
0185           HGCScintillatorDetId detId(type, layer, zside * ieta, iphi);
0186           if (typm.first >= 0) {
0187             detId.setType(typm.first);
0188             detId.setSiPM(typm.second);
0189           }
0190           DetId id1 = static_cast<DetId>(detId);
0191           if (geom->topology().valid(id1)) {
0192             auto icell1 = geom->getGeometry(id1);
0193             GlobalPoint global1 = geom->getPosition(id1);
0194             DetId idc1 = geom->getClosestCell(global1);
0195             GlobalPoint global2 = geom->getPosition(idc1);
0196             std::string cherr = (id1.rawId() != idc1.rawId()) ? " ***** ERROR *****" : "";
0197             edm::LogVerbatim("HGCalGeomX")
0198                 << "DetId (" << det << ":" << zside << ":" << type << ":" << layer << ":" << ieta << ":" << iphi
0199                 << ") Geom " << icell1 << " position (" << global1.x() << ", " << global1.y() << ", " << global1.z()
0200                 << ":" << global1.perp() << ") ids " << std::hex << id1.rawId() << ":" << idc1.rawId() << std::dec
0201                 << ":" << HGCScintillatorDetId(id1) << ":" << HGCScintillatorDetId(idc1) << " new position ("
0202                 << global2.x() << ", " << global2.y() << ", " << global2.z() << ":" << global2.perp()
0203                 << ") parameter[11] = " << icell1->param()[10] << ":" << icell1->param()[10] << cherr;
0204             std::vector<GlobalPoint> corners = geom->getCorners(idc1);
0205             std::ostringstream st1;
0206             st1 << corners.size() << " corners";
0207             for (auto const& cor : corners)
0208               st1 << " [" << cor.x() << "," << cor.y() << "," << cor.z() << "]";
0209             edm::LogVerbatim("HGCalGeomX") << st1.str();
0210           }
0211         }
0212       }
0213     }
0214   }
0215   for (int layer = geom->topology().dddConstants().firstLayer();
0216        layer <= geom->topology().dddConstants().lastLayer(true);
0217        ++layer) {
0218     HGCScintillatorDetId id(geom->topology().dddConstants().getTypeTrap(layer), layer, 20, 1);
0219     GlobalPoint global1 = geom->getPosition(id);
0220     edm::LogVerbatim("HGCalGeomX") << "Layer " << layer << " DetId " << id << " position (" << global1.x() << ", "
0221                                    << global1.y() << ", " << global1.z() << ", " << global1.perp() << ")";
0222   }
0223 }
0224 
0225 //define this as a plug-in
0226 DEFINE_FWK_MODULE(HGCalGeometryTester);