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