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   bool v17OrLess = hgdc.v17OrLess();
0080   edm::LogVerbatim("HGCalGeomX") << nameDetector_ << "\n Mode = " << mode;
0081   if (hgdc.waferHexagon8()) {
0082     double r = hgdc.waferParameters(true).first;
0083     double R = hgdc.waferParameters(true).second;
0084     edm::LogVerbatim("HGCalGeomX") << "Wafer Parameters " << r << ":" << R << std::endl;
0085     // Determine if the 24 points on the perphery of the wafer is within range
0086     // These are the 6 corners; the middle of the edges and the positions
0087     // which determines the positions of choptwoMinus and semiMinus
0088     // Offset of these points wrt the center of the wafer is given in dx, dy
0089     static const unsigned int nc = 24;
0090     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,
0091                      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};
0092     double dy[nc] = {-R,       -0.875 * R, -0.75 * R, -0.625 * R, -0.50 * R, -0.25 * R,  0.0,       0.25 * R,
0093                      0.50 * R, 0.625 * R,  0.75 * R,  0.875 * R,  R,         0.875 * R,  0.75 * R,  0.625 * R,
0094                      0.50 * R, 0.25 * R,   0.0,       -0.25 * R,  -0.50 * R, -0.625 * R, -0.75 * R, -0.875 * R};
0095     // There are 43 valid patterns corresponding to 7 types of partial wafers
0096     // in 6 orientation and one full wafer. These are are minmum requirerement
0097     // for these patterns
0098     static const unsigned int np = 43;
0099     unsigned int pat[np] = {0xFFFF01, 0xFFF01F, 0xFF01FF, 0xF01FFF, 0x01FFFF, 0x1FFFF0, 0xFFFC07, 0xFFC07F, 0xFC07FF,
0100                             0xC07FFF, 0x07FFFC, 0x7FFFC0, 0xFFF803, 0xFF803F, 0xF803FF, 0x803FFF, 0x03FFF8, 0x3FFF80,
0101                             0xFFF001, 0xFF001F, 0xF001FF, 0x001FFF, 0x01FFF0, 0x1FFF00, 0xFFC007, 0xFC007F, 0xC007FF,
0102                             0x007FFC, 0x07FFC0, 0x7FFC00, 0xFF8003, 0xF8003F, 0x8003FF, 0x003FF8, 0x03FF80, 0x3FF800,
0103                             0xFF0001, 0xF0001F, 0x0001FF, 0x001FF0, 0x01FF00, 0x1FF000, 0xFFFFFF};
0104     const std::vector<int> partTypeNew = {HGCalTypes::WaferLDBottom,
0105                                           HGCalTypes::WaferLDLeft,
0106                                           HGCalTypes::WaferLDRight,
0107                                           HGCalTypes::WaferLDFive,
0108                                           HGCalTypes::WaferLDThree,
0109                                           HGCalTypes::WaferHDTop,
0110                                           HGCalTypes::WaferHDBottom,
0111                                           HGCalTypes::WaferHDLeft,
0112                                           HGCalTypes::WaferHDRight,
0113                                           HGCalTypes::WaferHDFive};
0114     const std::vector<int> partTypeOld = {HGCalTypes::WaferHalf,
0115                                           HGCalTypes::WaferHalf,
0116                                           HGCalTypes::WaferSemi,
0117                                           HGCalTypes::WaferSemi,
0118                                           HGCalTypes::WaferFive,
0119                                           HGCalTypes::WaferThree,
0120                                           HGCalTypes::WaferHalf2,
0121                                           HGCalTypes::WaferChopTwoM,
0122                                           HGCalTypes::WaferSemi2,
0123                                           HGCalTypes::WaferSemi2,
0124                                           HGCalTypes::WaferFive2};
0125     const std::vector<DetId>& ids = geom->getValidGeomDetIds();
0126     int all(0), total(0), good(0), bad(0);
0127     for (auto id : ids) {
0128       HGCSiliconDetId hid(id);
0129       auto type = hgdc.waferTypeRotation(hid.layer(), hid.waferU(), hid.waferV(), false, false);
0130       if (hid.zside() > 0)
0131         ++all;
0132       // Not a full wafer
0133       int part = type.first;
0134       if (part > 10) {
0135         auto itr = std::find(partTypeNew.begin(), partTypeNew.end(), part);
0136         if (itr != partTypeNew.end()) {
0137           unsigned int indx = static_cast<unsigned int>(itr - partTypeNew.begin());
0138           part = partTypeOld[indx];
0139         }
0140       }
0141       if (part > 0 && part < 10 && hid.zside() > 0) {
0142         ++total;
0143         int wtype = hgdc.waferType(hid.layer(), hid.waferU(), hid.waferV(), false);
0144         int indx = (part - 1) * 6 + type.second;
0145         GlobalPoint xyz = geom->getWaferPosition(id);
0146         auto range = hgdc.rangeRLayer(hid.layer(), true);
0147         unsigned int ipat(0), ii(1);
0148         for (unsigned int i = 0; i < nc; ++i) {
0149           double rp = std::sqrt((xyz.x() + dx[i]) * (xyz.x() + dx[i]) + (xyz.y() + dy[i]) * (xyz.y() + dy[i]));
0150           if ((rp >= range.first) && (rp <= range.second))
0151             ipat += ii;
0152           ii *= 2;
0153         }
0154         bool match = (ipat == pat[indx]);
0155         if (!match) {
0156           // Make sure the minimum requirement is satisfied
0157           ii = 1;
0158           match = true;
0159           for (unsigned int i = 0; i < nc; ++i) {
0160             if ((((pat[indx] / ii) & 1) != 0) && (((ipat / ii) & 1) == 0)) {
0161               edm::LogVerbatim("HGCalGeomX") << "Fail at " << i << ":" << ii << " Expect " << ((pat[indx] / ii) & 1)
0162                                              << " Found " << ((ipat / ii) & 1);
0163               match = false;
0164               break;
0165             }
0166             ii *= 2;
0167           }
0168           if (match) {
0169             // and it doe not satify the higher ups
0170             if (wtype == 0) {
0171               match = (static_cast<unsigned int>(std::find(pat, pat + np, ipat) - pat) >= np);
0172             } else {
0173               // for coarse wafers the "minus" types are not allowed
0174               for (unsigned int i = 0; i < np; ++i) {
0175                 if (i < 12 || (i >= 18 && i < 30) || (i >= 36)) {
0176                   if (ipat == pat[i]) {
0177                     match = false;
0178                     break;
0179                   }
0180                 }
0181               }
0182             }
0183           }
0184         }
0185         std::string cherr = (!match) ? " ***** ERROR *****" : "";
0186         edm::LogVerbatim("HGCalGeomX") << "Wafer[" << wtype << ", " << hid.layer() << ", " << hid.waferU() << ", "
0187                                        << hid.waferV() << "]  with type: rotation " << type.first << "(" << part
0188                                        << "):" << type.second << " Pattern " << std::hex << pat[indx] << ":" << ipat
0189                                        << std::dec << cherr;
0190         if (!match) {
0191           ++bad;
0192           // Need debug information here
0193           hgdc.waferTypeRotation(hid.layer(), hid.waferU(), hid.waferV(), true, false);
0194           HGCalWaferMask::getTypeMode(xyz.x(), xyz.y(), r, R, range.first, range.second, wtype, 0, v17OrLess, true);
0195           for (unsigned int i = 0; i < 24; ++i) {
0196             double rp = std::sqrt((xyz.x() + dx[i]) * (xyz.x() + dx[i]) + (xyz.y() + dy[i]) * (xyz.y() + dy[i]));
0197             edm::LogVerbatim("HGCalGeomX")
0198                 << "Corner[" << i << "] (" << xyz.x() << ":" << xyz.y() << ") (" << (xyz.x() + dx[i]) << ":"
0199                 << (xyz.y() + dy[i]) << " ) R " << rp << " Limit " << range.first << ":" << range.second;
0200           }
0201         } else {
0202           ++good;
0203         }
0204         ++total;
0205       }
0206     }
0207     edm::LogVerbatim("HGCalGeomX") << "\n\nExamined " << ids.size() << ":" << all << " wafers " << total
0208                                    << " partial wafers of which " << good << " are good and " << bad << " are bad\n";
0209   }
0210 }
0211
0212 // define this as a plug-in
0213 DEFINE_FWK_MODULE(HGCalWaferTypeTester);