File indexing completed on 2022-07-12 02:58:21
0001 #include <iostream>
0002 #include <fstream>
0003 #include <string>
0004 #include <utility>
0005 #include <vector>
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESTransientHandle.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015
0016 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0017
0018 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0019 #include "DetectorDescription/Core/interface/DDSolid.h"
0020 #include "DetectorDescription/Core/interface/DDSolidShapes.h"
0021
0022 #include "Geometry/MTDCommonData/interface/MTDBaseNumber.h"
0023 #include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h"
0024 #include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h"
0025
0026 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0027 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0028
0029 #include "DataFormats/Math/interface/angle_units.h"
0030 #include "DataFormats/Math/interface/Rounding.h"
0031
0032
0033
0034 class TestMTDIdealGeometry : public edm::one::EDAnalyzer<> {
0035 public:
0036 explicit TestMTDIdealGeometry(const edm::ParameterSet&);
0037 ~TestMTDIdealGeometry() override;
0038
0039 void beginJob() override {}
0040 void analyze(edm::Event const&, edm::EventSetup const&) override;
0041 void endJob() override {}
0042
0043 void theBaseNumber(const DDGeoHistory& gh);
0044
0045 private:
0046 int nNodes_;
0047 std::string ddTopNodeName_;
0048
0049 MTDBaseNumber thisN_;
0050 BTLNumberingScheme btlNS_;
0051 ETLNumberingScheme etlNS_;
0052
0053 edm::ESGetToken<DDCompactView, IdealGeometryRecord> cpvToken_;
0054 };
0055
0056 TestMTDIdealGeometry::TestMTDIdealGeometry(const edm::ParameterSet& iConfig)
0057 : ddTopNodeName_(iConfig.getUntrackedParameter<std::string>("ddTopNodeName", "BarrelTimingLayer")),
0058 thisN_(),
0059 btlNS_(),
0060 etlNS_() {
0061 cpvToken_ = esConsumes<DDCompactView, IdealGeometryRecord>();
0062 }
0063
0064 TestMTDIdealGeometry::~TestMTDIdealGeometry() {}
0065
0066 using angle_units::operators::convertRadToDeg;
0067 using cms_rounding::roundIfNear0;
0068
0069 void TestMTDIdealGeometry::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0070 using namespace edm;
0071
0072 if (ddTopNodeName_ != "BarrelTimingLayer" && ddTopNodeName_ != "EndcapTimingLayer") {
0073 edm::LogWarning("TestMTDIdealGeometry") << ddTopNodeName_ << "Not valid top MTD volume";
0074 return;
0075 }
0076
0077 auto pDD = iSetup.getTransientHandle(cpvToken_);
0078
0079 if (!pDD.isValid()) {
0080 edm::LogError("TestMTDIdealGeometry") << "ESTransientHandle<DDCompactView> pDD is not valid!";
0081 return;
0082 }
0083 if (pDD.description()) {
0084 edm::LogInfo("TestMTDIdealGeometry") << pDD.description()->type_ << " label: " << pDD.description()->label_;
0085 } else {
0086 edm::LogWarning("TestMTDIdealGeometry") << "NO label found pDD.description() returned false.";
0087 }
0088
0089 DDPassAllFilter filter;
0090 DDFilteredView fv(*pDD, filter);
0091
0092 edm::LogInfo("TestMTDIdealGeometry") << "Top Most LogicalPart = " << fv.logicalPart();
0093
0094 using nav_type = DDFilteredView::nav_type;
0095 using id_type = std::map<nav_type, int>;
0096 id_type idMap;
0097 int id = 0;
0098
0099 bool write = false;
0100 bool isBarrel = true;
0101 size_t limit = 0;
0102
0103 do {
0104 nav_type pos = fv.navPos();
0105 idMap[pos] = id;
0106
0107 size_t num = fv.geoHistory().size();
0108
0109 if (num <= limit) {
0110 write = false;
0111 }
0112 if (fv.geoHistory()[num - 1].logicalPart().name() == "btl:BarrelTimingLayer") {
0113 isBarrel = true;
0114 limit = num;
0115 write = true;
0116 #ifdef EDM_ML_DEBUG
0117 edm::LogInfo("TestMTDIdealGeometry") << "isBarrel = " << isBarrel;
0118 #endif
0119 } else if (fv.geoHistory()[num - 1].logicalPart().name() == "etl:EndcapTimingLayer") {
0120 isBarrel = false;
0121 limit = num;
0122 write = true;
0123 #ifdef EDM_ML_DEBUG
0124 edm::LogInfo("TestMTDIdealGeometry") << "isBarrel = " << isBarrel;
0125 #endif
0126 }
0127
0128
0129
0130 std::stringstream ss;
0131 auto print_path = [&]() {
0132 ss << " - OCMS[0]/";
0133 for (uint i = 1; i < fv.geoHistory().size(); i++) {
0134 ss << fv.geoHistory()[i].logicalPart().name().fullname();
0135 ss << "[";
0136 ss << std::to_string(fv.geoHistory()[i].copyno());
0137 ss << "]/";
0138 }
0139 };
0140
0141 if (write && fv.geoHistory()[limit - 1].logicalPart().name().name() == ddTopNodeName_) {
0142 print_path();
0143 edm::LogInfo("TestMTDPath") << ss.str();
0144
0145 bool isSens = false;
0146
0147 if (!fv.geoHistory()[num - 1].logicalPart().specifics().empty()) {
0148 for (auto vec : fv.geoHistory()[num - 1].logicalPart().specifics()) {
0149 for (const auto& elem : *vec) {
0150 if (elem.second.name() == "SensitiveDetector") {
0151 isSens = true;
0152 break;
0153 }
0154 }
0155 }
0156 }
0157
0158 if (isSens) {
0159 theBaseNumber(fv.geoHistory());
0160
0161
0162
0163
0164
0165 std::stringstream sunitt;
0166 std::stringstream snum;
0167
0168 if (isBarrel) {
0169 BTLDetId theId(btlNS_.getUnitID(thisN_));
0170 sunitt << theId.rawId();
0171 snum << theId;
0172 snum << "\n";
0173 } else {
0174 ETLDetId theId(etlNS_.getUnitID(thisN_));
0175 sunitt << theId.rawId();
0176 snum << theId;
0177 #ifdef EDM_ML_DEBUG
0178 edm::LogInfo("TestMTDNumbering")
0179 << " ETLDetId = " << theId << "\n geographicalId = " << theId.geographicalId();
0180 #endif
0181 }
0182 edm::LogInfo("TestMTDNumbering") << snum.str();
0183
0184
0185
0186
0187
0188 std::stringstream spos;
0189
0190 auto fround = [&](double in) {
0191 std::stringstream ss;
0192 ss << std::fixed << std::setw(14) << roundIfNear0(in);
0193 return ss.str();
0194 };
0195
0196 DDBox mySens = fv.geoHistory()[num - 1].logicalPart().solid();
0197 spos << "Solid shape name: " << DDSolidShapesName::name(mySens.shape()) << "\n";
0198 if (static_cast<int>(mySens.shape()) != 1) {
0199 throw cms::Exception("TestMTDPosition") << "MTD sensitive element not a DDBox";
0200 break;
0201 }
0202 spos << "Box dimensions: " << fround(mySens.halfX()) << " " << fround(mySens.halfY()) << " "
0203 << fround(mySens.halfZ()) << "\n";
0204
0205 DD3Vector x, y, z;
0206 fv.rotation().GetComponents(x, y, z);
0207 spos << "Translation vector components: " << fround(fv.translation().x()) << " " << fround(fv.translation().y())
0208 << " " << fround(fv.translation().z()) << " "
0209 << "\n";
0210 spos << "Rotation matrix components: " << fround(x.X()) << " " << fround(x.Y()) << " " << fround(x.Z()) << " "
0211 << fround(y.X()) << " " << fround(y.Y()) << " " << fround(y.Z()) << " " << fround(z.X()) << " "
0212 << fround(z.Y()) << " " << fround(z.Z()) << "\n";
0213
0214 DD3Vector zeroLocal(0., 0., 0.);
0215 DD3Vector cn1Local(mySens.halfX(), mySens.halfY(), mySens.halfZ());
0216 double distLocal = cn1Local.R();
0217 DD3Vector zeroGlobal = (fv.rotation())(zeroLocal) + fv.translation();
0218 DD3Vector cn1Global = (fv.rotation())(cn1Local) + fv.translation();
0219 ;
0220 double distGlobal =
0221 std::sqrt(std::pow(zeroGlobal.X() - cn1Global.X(), 2) + std::pow(zeroGlobal.Y() - cn1Global.Y(), 2) +
0222 std::pow(zeroGlobal.Z() - cn1Global.Z(), 2));
0223
0224 spos << "Center global = " << fround(zeroGlobal.X()) << fround(zeroGlobal.Y()) << fround(zeroGlobal.Z())
0225 << " r = " << fround(zeroGlobal.Rho()) << " phi = " << fround(convertRadToDeg(zeroGlobal.Phi())) << "\n";
0226
0227 spos << "Corner 1 local = " << fround(cn1Local.X()) << fround(cn1Local.Y()) << fround(cn1Local.Z())
0228 << " DeltaR = " << fround(distLocal) << "\n";
0229
0230 spos << "Corner 1 global = " << fround(cn1Global.X()) << fround(cn1Global.Y()) << fround(cn1Global.Z())
0231 << " DeltaR = " << fround(distGlobal) << "\n";
0232
0233 spos << "\n";
0234 if (std::fabs(distGlobal - distLocal) > 1.e-6) {
0235 spos << "DIFFERENCE IN DISTANCE \n";
0236 }
0237 sunitt << fround(zeroGlobal.X()) << fround(zeroGlobal.Y()) << fround(zeroGlobal.Z()) << fround(cn1Global.X())
0238 << fround(cn1Global.Y()) << fround(cn1Global.Z());
0239 edm::LogInfo("TestMTDPosition") << spos.str();
0240
0241 edm::LogVerbatim("MTDUnitTest") << sunitt.str();
0242 }
0243 }
0244 ++id;
0245 } while (fv.next());
0246 }
0247
0248 void TestMTDIdealGeometry::theBaseNumber(const DDGeoHistory& gh) {
0249 thisN_.reset();
0250 thisN_.setSize(gh.size());
0251
0252 for (uint i = gh.size(); i-- > 0;) {
0253 thisN_.addLevel(gh[i].logicalPart().name().name(), gh[i].copyno());
0254 #ifdef EDM_ML_DEBUG
0255 edm::LogInfo("TestMTDIdealGeometry") << gh[i].logicalPart().name().name() << " " << gh[i].copyno();
0256 #endif
0257 }
0258 }
0259
0260 DEFINE_FWK_MODULE(TestMTDIdealGeometry);