HGCalGeometry/test/HGCalWaferTypeTester.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    HGCalWaferTypeTester
0004 // Class:      HGCalWaferTypeTester
0005 //
0006 /**\class HGCalWaferTypeTester HGCalWaferTypeTester.cc
0007  test/HGCalWaferTypeTester.cc
0008
0009  Description: <one line class summary>
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author:  Sunanda Banerjee
0016 //         Created:  Mon 2020/06/08
0017 //
0018 //
0019
0020 // system include files
0021 #include <fstream>
0022 #include <iostream>
0023 #include <memory>
0024 #include <string>
0025 #include <vector>
0026
0027 // user include files
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/EventSetup.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0036
0037 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0038 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0039 #include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
0040 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0041 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0042
0043 class HGCalWaferTypeTester : public edm::one::EDAnalyzer<> {
0044 public:
0045   explicit HGCalWaferTypeTester(const edm::ParameterSet&);
0046
0047   static void fillDescriptions(edm::ConfigurationDescriptions&);
0048
0049   void beginJob() override {}
0050   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0051   void endJob() override {}
0052
0053 private:
0054   const std::string nameSense_, nameDetector_;
0055   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0056 };
0057
0058 HGCalWaferTypeTester::HGCalWaferTypeTester(const edm::ParameterSet& iC)
0059     : nameSense_(iC.getParameter<std::string>("NameSense")),
0060       nameDetector_(iC.getParameter<std::string>("NameDevice")),
0061       geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameSense_})) {
0062   edm::LogVerbatim("HGCalGeomX") << "Test wafer types for " << nameDetector_ << " using constants of " << nameSense_
0063                                  << " for  RecoFlag true";
0064 }
0065
0066 void HGCalWaferTypeTester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0067   edm::ParameterSetDescription desc;
0068   desc.add<std::string>("NameSense", "HGCalEESensitive");
0069   desc.add<std::string>("NameDevice", "HGCal EE");
0070   descriptions.add("hgcalEEWaferTypeTester", desc);
0071 }
0072
0073 // ------------ method called to produce the data  ------------
0074 void HGCalWaferTypeTester::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0075   const auto& geomR = iSetup.getData(geomToken_);
0076   const HGCalGeometry* geom = &geomR;
0077   const HGCalDDDConstants& hgdc = geom->topology().dddConstants();
0078   HGCalGeometryMode::GeometryMode mode = hgdc.geomMode();
0079   edm::LogVerbatim("HGCalGeomX") << nameDetector_ << "\n Mode = " << mode;
0080   if (hgdc.waferHexagon8()) {
0081     double r = hgdc.waferParameters(true).first;
0082     double R = hgdc.waferParameters(true).second;
0083     edm::LogVerbatim("HGCalGeomX") << "Wafer Parameters " << r << ":" << R << std::endl;
0084     // Determine if the 24 points on the perphery of the wafer is within range
0085     // These are the 6 corners; the middle of the edges and the positions
0086     // which determines the positions of choptwoMinus and semiMinus
0087     // Offset of these points wrt the center of the wafer is given in dx, dy
0088     static const unsigned int nc = 24;
0089     double dx[nc] = {0.0, 0.25 * r,  0.50 * r,  0.75 * r,  r,  r,  r,  r,  r,  0.75 * r,  0.50 * r,  0.25 * r,
0090                      0.0, -0.25 * r, -0.50 * r, -0.75 * r, -r, -r, -r, -r, -r, -0.75 * r, -0.50 * r, -0.25 * r};
0091     double dy[nc] = {-R,       -0.875 * R, -0.75 * R, -0.625 * R, -0.50 * R, -0.25 * R,  0.0,       0.25 * R,
0092                      0.50 * R, 0.625 * R,  0.75 * R,  0.875 * R,  R,         0.875 * R,  0.75 * R,  0.625 * R,
0093                      0.50 * R, 0.25 * R,   0.0,       -0.25 * R,  -0.50 * R, -0.625 * R, -0.75 * R, -0.875 * R};
0094     // There are 43 valid patterns corresponding to 7 types of partial wafers
0095     // in 6 orientation and one full wafer. These are are minmum requirerement
0096     // for these patterns
0097     static const unsigned int np = 43;
0098     unsigned int pat[np] = {0xFFFF01, 0xFFF01F, 0xFF01FF, 0xF01FFF, 0x01FFFF, 0x1FFFF0, 0xFFFC07, 0xFFC07F, 0xFC07FF,
0099                             0xC07FFF, 0x07FFFC, 0x7FFFC0, 0xFFF803, 0xFF803F, 0xF803FF, 0x803FFF, 0x03FFF8, 0x3FFF80,
0100                             0xFFF001, 0xFF001F, 0xF001FF, 0x001FFF, 0x01FFF0, 0x1FFF00, 0xFFC007, 0xFC007F, 0xC007FF,
0101                             0x007FFC, 0x07FFC0, 0x7FFC00, 0xFF8003, 0xF8003F, 0x8003FF, 0x003FF8, 0x03FF80, 0x3FF800,
0102                             0xFF0001, 0xF0001F, 0x0001FF, 0x001FF0, 0x01FF00, 0x1FF000, 0xFFFFFF};
0103     const std::vector<int> partTypeNew = {HGCalTypes::WaferLDBottom,
0104                                           HGCalTypes::WaferLDLeft,
0105                                           HGCalTypes::WaferLDRight,
0106                                           HGCalTypes::WaferLDFive,
0107                                           HGCalTypes::WaferLDThree,
0108                                           HGCalTypes::WaferHDTop,
0109                                           HGCalTypes::WaferHDBottom,
0110                                           HGCalTypes::WaferHDLeft,
0111                                           HGCalTypes::WaferHDRight,
0112                                           HGCalTypes::WaferHDFive};
0113     const std::vector<int> partTypeOld = {HGCalTypes::WaferHalf,
0114                                           HGCalTypes::WaferHalf,
0115                                           HGCalTypes::WaferSemi,
0116                                           HGCalTypes::WaferSemi,
0117                                           HGCalTypes::WaferFive,
0118                                           HGCalTypes::WaferThree,
0119                                           HGCalTypes::WaferHalf2,
0120                                           HGCalTypes::WaferChopTwoM,
0121                                           HGCalTypes::WaferSemi2,
0122                                           HGCalTypes::WaferSemi2,
0123                                           HGCalTypes::WaferFive2};
0124     const std::vector<DetId>& ids = geom->getValidGeomDetIds();
0125     int all(0), total(0), good(0), bad(0);
0126     for (auto id : ids) {
0127       HGCSiliconDetId hid(id);
0128       auto type = hgdc.waferTypeRotation(hid.layer(), hid.waferU(), hid.waferV(), false, false);
0129       if (hid.zside() > 0)
0130         ++all;
0131       // Not a full wafer
0132       int part = type.first;
0133       if (part > 10) {
0134         auto itr = std::find(partTypeNew.begin(), partTypeNew.end(), part);
0135         if (itr != partTypeNew.end()) {
0136           unsigned int indx = static_cast<unsigned int>(itr - partTypeNew.begin());
0137           part = partTypeOld[indx];
0138         }
0139       }
0140       if (part > 0 && part < 10 && hid.zside() > 0) {
0141         ++total;
0142         int wtype = hgdc.waferType(hid.layer(), hid.waferU(), hid.waferV(), false);
0143         int indx = (part - 1) * 6 + type.second;
0144         GlobalPoint xyz = geom->getWaferPosition(id);
0145         auto range = hgdc.rangeRLayer(hid.layer(), true);
0146         unsigned int ipat(0), ii(1);
0147         for (unsigned int i = 0; i < nc; ++i) {
0148           double rp = std::sqrt((xyz.x() + dx[i]) * (xyz.x() + dx[i]) + (xyz.y() + dy[i]) * (xyz.y() + dy[i]));
0149           if ((rp >= range.first) && (rp <= range.second))
0150             ipat += ii;
0151           ii *= 2;
0152         }
0153         bool match = (ipat == pat[indx]);
0154         if (!match) {
0155           // Make sure the minimum requirement is satisfied
0156           ii = 1;
0157           match = true;
0158           for (unsigned int i = 0; i < nc; ++i) {
0159             if ((((pat[indx] / ii) & 1) != 0) && (((ipat / ii) & 1) == 0)) {
0160               edm::LogVerbatim("HGCalGeomX") << "Fail at " << i << ":" << ii << " Expect " << ((pat[indx] / ii) & 1)
0161                                              << " Found " << ((ipat / ii) & 1);
0162               match = false;
0163               break;
0164             }
0165             ii *= 2;
0166           }
0167           if (match) {
0168             // and it doe not satify the higher ups
0169             if (wtype == 0) {
0170               match = (static_cast<unsigned int>(std::find(pat, pat + np, ipat) - pat) >= np);
0171             } else {
0172               // for coarse wafers the "minus" types are not allowed
0173               for (unsigned int i = 0; i < np; ++i) {
0174                 if (i < 12 || (i >= 18 && i < 30) || (i >= 36)) {
0175                   if (ipat == pat[i]) {
0176                     match = false;
0177                     break;
0178                   }
0179                 }
0180               }
0181             }
0182           }
0183         }
0184         std::string cherr = (!match) ? " ***** ERROR *****" : "";
0185         edm::LogVerbatim("HGCalGeomX") << "Wafer[" << wtype << ", " << hid.layer() << ", " << hid.waferU() << ", "
0186                                        << hid.waferV() << "]  with type: rotation " << type.first << "(" << part
0187                                        << "):" << type.second << " Pattern " << std::hex << pat[indx] << ":" << ipat
0188                                        << std::dec << cherr;
0189         if (!match) {
0190           ++bad;
0191           // Need debug information here
0192           hgdc.waferTypeRotation(hid.layer(), hid.waferU(), hid.waferV(), true, false);
0193           HGCalWaferMask::getTypeMode(xyz.x(), xyz.y(), r, R, range.first, range.second, wtype, 0, true);
0194           for (unsigned int i = 0; i < 24; ++i) {
0195             double rp = std::sqrt((xyz.x() + dx[i]) * (xyz.x() + dx[i]) + (xyz.y() + dy[i]) * (xyz.y() + dy[i]));
0196             edm::LogVerbatim("HGCalGeomX")
0197                 << "Corner[" << i << "] (" << xyz.x() << ":" << xyz.y() << ") (" << (xyz.x() + dx[i]) << ":"
0198                 << (xyz.y() + dy[i]) << " ) R " << rp << " Limit " << range.first << ":" << range.second;
0199           }
0200         } else {
0201           ++good;
0202         }
0203         ++total;
0204       }
0205     }
0206     edm::LogVerbatim("HGCalGeomX") << "\n\nExamined " << ids.size() << ":" << all << " wafers " << total
0207                                    << " partial wafers of which " << good << " are good and " << bad << " are bad\n";
0208   }
0209 }
0210
0211 // define this as a plug-in
0212 DEFINE_FWK_MODULE(HGCalWaferTypeTester);