Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-29 07:44:26

0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004 
0005 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0006 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0008 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0010 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0011 
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0020 
0021 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0022 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0023 #include "MagneticField/Engine/interface/MagneticField.h"
0024 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0025 #include "CoralBase/Exception.h"
0026 
0027 class HGCalTestNeighbor : public edm::one::EDAnalyzer<> {
0028 public:
0029   explicit HGCalTestNeighbor(const edm::ParameterSet&);
0030   ~HGCalTestNeighbor() override = default;
0031 
0032   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0033   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0034 
0035 private:
0036   void doTest(const HGCalGeometry* geom, ForwardSubdetector subdet, const MagneticField* bField);
0037   void doTestWafer(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField);
0038   void doTestScint(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField);
0039 
0040   const std::string name_;
0041   const std::vector<double> px_, py_, pz_;
0042   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> fieldToken_;
0043   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0044 };
0045 
0046 HGCalTestNeighbor::HGCalTestNeighbor(const edm::ParameterSet& iC)
0047     : name_(iC.getParameter<std::string>("detector")),
0048       px_(iC.getParameter<std::vector<double> >("pX")),
0049       py_(iC.getParameter<std::vector<double> >("pY")),
0050       pz_(iC.getParameter<std::vector<double> >("pZ")),
0051       fieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>(edm::ESInputTag{})),
0052       geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", name_})) {}
0053 
0054 void HGCalTestNeighbor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0055   edm::ParameterSetDescription desc;
0056   desc.add<std::string>("detector", "HGCalEESensitive");
0057   std::vector<double> pxy = {2.0, 5.0, 10.0, 100.0};
0058   std::vector<double> pz = {10.0, 5.0, 20.0, 50.0};
0059   desc.add<std::vector<double> >("pX", pxy);
0060   desc.add<std::vector<double> >("pY", pxy);
0061   desc.add<std::vector<double> >("pZ", pz);
0062   descriptions.add("hgcalEETestNeighbor", desc);
0063 }
0064 
0065 void HGCalTestNeighbor::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0066   const auto& bFieldR = iSetup.getData(fieldToken_);
0067   const MagneticField* bField = &bFieldR;
0068 
0069   const auto& geomR = iSetup.getData(geomToken_);
0070   const HGCalGeometry* geom = &geomR;
0071   if (geom->topology().waferHexagon6()) {
0072     ForwardSubdetector subdet;
0073     if (name_ == "HGCalHESiliconSensitive")
0074       subdet = HGCHEF;
0075     else if (name_ == "HGCalHEScintillatorSensitive")
0076       subdet = HGCHEB;
0077     else
0078       subdet = HGCEE;
0079     edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name_ << " Detector:Subdetector " << DetId::Forward << ":"
0080                                    << subdet;
0081     doTest(geom, subdet, bField);
0082   } else {
0083     DetId::Detector det;
0084     if (name_ == "HGCalHESiliconSensitive")
0085       det = DetId::HGCalHSi;
0086     else if (name_ == "HGCalHEScintillatorSensitive")
0087       det = DetId::HGCalHSc;
0088     else
0089       det = DetId::HGCalEE;
0090     edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name_ << " Detector " << det;
0091     if (name_ == "HGCalHEScintillatorSensitive") {
0092       doTestScint(geom, det, bField);
0093     } else {
0094       doTestWafer(geom, det, bField);
0095     }
0096   }
0097 }
0098 
0099 void HGCalTestNeighbor::doTest(const HGCalGeometry* geom, ForwardSubdetector subdet, const MagneticField* bField) {
0100   const std::vector<DetId>& ids = geom->getValidDetIds();
0101   edm::LogVerbatim("HGCalGeomX") << "doTest: " << ids.size() << " valid ids for " << geom->cellElement();
0102 
0103   int layers[] = {1, 5, 10};
0104   int zsides[] = {1, -1};
0105   int cells[] = {1, 51, 101};
0106   int wafers[] = {1, 101, 201, 301, 401};
0107   const int ismax(5);
0108   for (int zside : zsides) {
0109     for (int is = 0; is < ismax; ++is) {
0110       int sector = wafers[is];
0111       int type = geom->topology().dddConstants().waferTypeT(sector);
0112       if (type != 1)
0113         type = 0;
0114       for (int layer : layers) {
0115         for (int cell : cells) {
0116           DetId id1;
0117           id1 = (DetId)(HGCalDetId(subdet, zside, layer, type, sector, cell));
0118           if (geom->topology().valid(id1)) {
0119             auto icell1 = geom->getGeometry(id1);
0120             GlobalPoint global1 = geom->getPosition(id1);
0121             for (unsigned int k = 0; k < px_.size(); ++k) {
0122               GlobalVector p(px_[k], py_[k], zside * pz_[k]);
0123               DetId id2 = geom->neighborZ(id1, p);
0124               DetId id3 = geom->neighborZ(id1, bField, 1, p);
0125               edm::LogVerbatim("HGCalGeomX")
0126                   << "DetId" << HGCalDetId(id1) << " :" << global1 << " p" << p << " ID2" << HGCalDetId(id2) << " :"
0127                   << geom->getPosition(id2) << " ID3" << HGCalDetId(id3) << " :" << geom->getPosition(id3);
0128             }
0129           }
0130         }
0131       }
0132     }
0133   }
0134 }
0135 
0136 void HGCalTestNeighbor::doTestWafer(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField) {
0137   const std::vector<DetId>& ids = geom->getValidDetIds();
0138   edm::LogVerbatim("HGCalGeomX") << "doTestWafer:: " << ids.size() << " valid ids for " << geom->cellElement();
0139   int layers[] = {1, 5, 10};
0140   int zsides[] = {1, -1};
0141   int cells[] = {1, 4, 7};
0142   int wafers[] = {7, 5, 3, -3, -5, -7};
0143   for (int zside : zsides) {
0144     for (int layer : layers) {
0145       for (int waferU : wafers) {
0146         for (int waferV : wafers) {
0147           int type = geom->topology().dddConstants().getTypeHex(layer, waferU, waferV);
0148           for (int cellU : cells) {
0149             for (int cellV : cells) {
0150               DetId id1 = (DetId)(HGCSiliconDetId(det, zside, type, layer, waferU, waferV, cellU, cellV));
0151               if (geom->topology().valid(id1)) {
0152                 auto icell1 = geom->getGeometry(id1);
0153                 GlobalPoint global1 = geom->getPosition(id1);
0154                 for (unsigned int k = 0; k < px_.size(); ++k) {
0155                   GlobalVector p(px_[k], py_[k], zside * pz_[k]);
0156                   DetId id2 = geom->neighborZ(id1, p);
0157                   DetId id3 = geom->neighborZ(id1, bField, 1, p);
0158                   edm::LogVerbatim("HGCalGeomX") << "DetId" << HGCSiliconDetId(id1) << " :" << global1 << " p" << p
0159                                                  << " ID2" << HGCSiliconDetId(id2) << " :" << geom->getPosition(id2)
0160                                                  << " ID3" << HGCSiliconDetId(id3) << " :" << geom->getPosition(id3);
0161                 }
0162               }
0163             }
0164           }
0165         }
0166       }
0167     }
0168   }
0169 }
0170 
0171 void HGCalTestNeighbor::doTestScint(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField) {
0172   const std::vector<DetId>& ids = geom->getValidDetIds();
0173   edm::LogVerbatim("HGCalGeomX") << "doTestScint: " << ids.size() << " valid ids for " << geom->cellElement();
0174   int layers[] = {9, 15, 22};
0175   int zsides[] = {1, -1};
0176   int iphis[] = {1, 51, 101, 151, 201};
0177   int ietas[] = {11, 20, 29};
0178   for (int zside : zsides) {
0179     for (int layer : layers) {
0180       int type = geom->topology().dddConstants().getTypeTrap(layer);
0181       for (int ieta : ietas) {
0182         for (int iphi : iphis) {
0183           DetId id1 = (DetId)(HGCScintillatorDetId(type, layer, zside * ieta, iphi));
0184           if (geom->topology().valid(id1)) {
0185             auto icell1 = geom->getGeometry(id1);
0186             GlobalPoint global1 = geom->getPosition(id1);
0187             for (unsigned int k = 0; k < px_.size(); ++k) {
0188               GlobalVector p(px_[k], py_[k], zside * pz_[k]);
0189               DetId id2 = geom->neighborZ(id1, p);
0190               DetId id3 = geom->neighborZ(id1, bField, 1, p);
0191               edm::LogVerbatim("HGCalGeomX") << "DetId" << HGCScintillatorDetId(id1) << " :" << global1 << " p" << p
0192                                              << " ID2" << HGCScintillatorDetId(id2) << " :" << geom->getPosition(id2)
0193                                              << " ID3" << HGCScintillatorDetId(id3) << " :" << geom->getPosition(id3);
0194             }
0195           }
0196         }
0197       }
0198     }
0199   }
0200 }
0201 
0202 //define this as a plug-in
0203 DEFINE_FWK_MODULE(HGCalTestNeighbor);