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