File indexing completed on 2024-04-06 12:15:28
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0032 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0033 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0034 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0035 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0036 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0037 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0038
0039 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0040 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0041 #include "DataFormats/GeometrySurface/interface/BoundSurface.h"
0042 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0043 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0044 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0045
0046
0047
0048
0049
0050
0051
0052 #define PRINT(X) std::cout << X << ':'
0053
0054 class TrackerDigiGeometryAnalyzer : public edm::one::EDAnalyzer<> {
0055 public:
0056 explicit TrackerDigiGeometryAnalyzer(const edm::ParameterSet&);
0057
0058 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0059
0060 private:
0061 void analyseTrapezoidal(const GeomDetUnit& det);
0062 void checkRotation(const GeomDetUnit& det);
0063 void checkTopology(const GeomDetUnit& det);
0064 std::ostream& cylindrical(std::ostream& os, const GlobalPoint& gp) const;
0065
0066 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> pDDToken_;
0067 };
0068
0069 TrackerDigiGeometryAnalyzer::TrackerDigiGeometryAnalyzer(const edm::ParameterSet& iConfig) : pDDToken_(esConsumes()) {}
0070
0071
0072 void TrackerDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0073 PRINT("TrackerDigiGeometryAnalyzer") << "Here I am";
0074
0075
0076
0077 auto const& pDD = iSetup.getData(pDDToken_);
0078 PRINT("TrackerDigiGeometryAnalyzer") << " Geometry node for TrackerGeom is " << &pDD << '\n';
0079 PRINT("TrackerDigiGeometryAnalyzer") << " I have " << pDD.detUnits().size() << " detectors" << '\n';
0080 PRINT("TrackerDigiGeometryAnalyzer") << " I have " << pDD.detTypes().size() << " types" << '\n';
0081
0082 for (auto const& it : pDD.detUnits()) {
0083 if (dynamic_cast<const PixelGeomDetUnit*>((it)) != nullptr) {
0084 const BoundPlane& p = (dynamic_cast<const PixelGeomDetUnit*>((it)))->specificSurface();
0085 PRINT("TrackerDigiGeometryAnalyzer") << it->geographicalId() << " RadLeng Pixel " << p.mediumProperties().radLen()
0086 << ' ' << " Xi Pixel " << p.mediumProperties().xi() << '\n';
0087 }
0088
0089 if (dynamic_cast<const StripGeomDetUnit*>((it)) != nullptr) {
0090 const BoundPlane& s = (dynamic_cast<const StripGeomDetUnit*>((it)))->specificSurface();
0091 PRINT("TrackerDigiGeometryAnalyzer") << it->geographicalId() << " RadLeng Strip " << s.mediumProperties().radLen()
0092 << " Xi Strip " << s.mediumProperties().xi() << '\n';
0093 }
0094
0095
0096 }
0097
0098 for (auto const& it : pDD.detTypes()) {
0099 if (dynamic_cast<const PixelGeomDetType*>((it)) != nullptr) {
0100 const PixelTopology& p = (dynamic_cast<const PixelGeomDetType*>((it)))->specificTopology();
0101 PRINT("TrackerDigiGeometryAnalyzer") << " PIXEL Det "
0102 << " Rows " << p.nrows() << " Columns " << p.ncolumns() << '\n';
0103 } else {
0104 const StripTopology& p = (dynamic_cast<const StripGeomDetType*>((it)))->specificTopology();
0105 PRINT("TrackerDigiGeometryAnalyzer") << " STRIP Det "
0106 << " Strips " << p.nstrips() << '\n';
0107 }
0108 }
0109 }
0110 void TrackerDigiGeometryAnalyzer::analyseTrapezoidal(const GeomDetUnit& det) {
0111
0112
0113 const double safety = 0.9999;
0114
0115 const Bounds& bounds = det.surface().bounds();
0116 const TrapezoidalPlaneBounds* tb = dynamic_cast<const TrapezoidalPlaneBounds*>(&bounds);
0117 if (tb == nullptr)
0118 return;
0119
0120 checkTopology(det);
0121
0122 const GlobalPoint& pos = det.position();
0123 double length = tb->length();
0124 double width = tb->width();
0125
0126 const std::array<const float, 4>& par = tb->parameters();
0127 double top = std::max(par[1], par[0]);
0128 double bot = std::min(par[1], par[0]);
0129
0130 std::cout << std::endl;
0131 std::cout << "Det at pos " << pos << " has length " << length << " width " << width << " pars ";
0132 for (int i = 0; i < 4; i++)
0133 std::cout << par[i] << ", ";
0134 std::cout << std::endl;
0135
0136 std::cout << "det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << std::endl;
0137
0138
0139
0140
0141 GlobalVector yShift = det.surface().toGlobal(LocalVector(0, safety * length / 2., 0));
0142 GlobalPoint outerMiddle = pos + yShift;
0143 GlobalPoint innerMiddle = pos + (-1. * yShift);
0144 if (outerMiddle.perp() < innerMiddle.perp())
0145 std::swap(outerMiddle, innerMiddle);
0146
0147 GlobalVector upperShift = det.surface().toGlobal(LocalVector(safety * top, 0, 0));
0148
0149 GlobalPoint ulc = outerMiddle + upperShift;
0150 GlobalPoint urc = outerMiddle + (-1. * upperShift);
0151 std::cout << "outerMiddle " << outerMiddle << " upperShift " << upperShift << " ulc "
0152 << "(" << ulc.perp() << "," << ulc.phi() << "," << ulc.z() << " urc "
0153 << "(" << urc.perp() << "," << urc.phi() << "," << urc.z() << std::endl;
0154
0155 std::cout << "upper corners inside bounds? " << tb->inside(det.surface().toLocal(ulc)) << " "
0156 << tb->inside(det.surface().toLocal(urc)) << std::endl;
0157
0158
0159
0160
0161 GlobalVector lowerShift = det.surface().toGlobal(LocalVector(safety * bot, 0, 0));
0162
0163 std::cout << "lower corners inside bounds? " << tb->inside(det.surface().toLocal(innerMiddle + lowerShift)) << " "
0164 << tb->inside(det.surface().toLocal(innerMiddle + (-1. * lowerShift))) << std::endl;
0165 }
0166
0167 void TrackerDigiGeometryAnalyzer::checkRotation(const GeomDetUnit& det) {
0168 const double eps = std::numeric_limits<float>::epsilon();
0169 static int first = 0;
0170 if (first == 0) {
0171 std::cout << "numeric_limits<float>::epsilon() " << std::numeric_limits<float>::epsilon() << std::endl;
0172 first = 1;
0173 }
0174
0175 const Surface::RotationType& rot(det.surface().rotation());
0176 GlobalVector a(rot.xx(), rot.xy(), rot.xz());
0177 GlobalVector b(rot.yx(), rot.yy(), rot.yz());
0178 GlobalVector c(rot.zx(), rot.zy(), rot.zz());
0179 GlobalVector cref = a.cross(b);
0180 GlobalVector aref = b.cross(c);
0181 GlobalVector bref = c.cross(a);
0182 if ((a - aref).mag() > eps || (b - bref).mag() > eps || (c - cref).mag() > eps) {
0183 std::cout << " Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag() << ", "
0184 << (c - cref).mag() << " for det at pos " << det.surface().position() << std::endl;
0185 }
0186 if (fabs(a.mag() - 1.) > eps || fabs(b.mag() - 1.) > eps || fabs(c.mag() - 1.) > eps) {
0187 std::cout << " Rotation not good by bector mag: " << (a).mag() << ", " << (b).mag() << ", " << (c).mag()
0188 << " for det at pos " << det.surface().position() << std::endl;
0189 }
0190 }
0191
0192 void TrackerDigiGeometryAnalyzer::checkTopology(const GeomDetUnit& det) {
0193 const StripTopology& topol = dynamic_cast<const StripTopology&>(det.topology());
0194 const int N = 5;
0195
0196 std::cout << std::endl << "Topology test along strip 0" << std::endl;
0197 LocalVector stripDir = topol.localPosition(MeasurementPoint(0, 0.5)) - topol.localPosition(MeasurementPoint(0, 0));
0198 std::cout << "StripDir " << stripDir << std::endl;
0199 for (int i = -N; i <= N; i++) {
0200 MeasurementPoint mp(0, float(i) / float(2 * N));
0201 LocalPoint lp = topol.localPosition(mp);
0202 double strlen = topol.localStripLength(lp);
0203 LocalError le = topol.localError(mp, MeasurementError(0.25, 0, 0.25));
0204 LocalError lep = le.rotate(stripDir.y(), stripDir.x());
0205 LocalError lem = le.rotate(stripDir.y(), -stripDir.x());
0206 GlobalPoint gp(det.surface().toGlobal(lp));
0207 std::cout << "gpos (r,phi) (" << gp.perp() << ", " << gp.phi() << ") lpos " << lp << " lerr.x (0,+,-) "
0208 << sqrt(le.xx()) << "," << sqrt(lep.xx()) << "," << sqrt(lem.xx()) << " strlen " << strlen << std::endl;
0209 }
0210
0211 std::cout << std::endl << "Topology test along middle strip" << std::endl;
0212 float midstrip = topol.strip(LocalPoint(0, 0));
0213 for (int i = -N; i <= N; i++) {
0214 MeasurementPoint mp(midstrip, float(i) / float(2 * N));
0215 LocalPoint lp = topol.localPosition(mp);
0216 double strlen = topol.localStripLength(lp);
0217 LocalError le = topol.localError(mp, MeasurementError(0.25, 0, 0.25));
0218 GlobalPoint gp(det.surface().toGlobal(lp));
0219 std::cout << "gpos (r,phi) (" << gp.perp() << ", " << gp.phi() << ") lpos " << lp << " lerr.x " << sqrt(le.xx())
0220 << " strlen " << strlen << std::endl;
0221 }
0222 }
0223
0224 std::ostream& TrackerDigiGeometryAnalyzer::cylindrical(std::ostream& os, const GlobalPoint& gp) const {
0225 os << "(" << gp.perp() << "," << gp.phi() << "," << gp.z();
0226 return os;
0227 }
0228
0229
0230 DEFINE_FWK_MODULE(TrackerDigiGeometryAnalyzer);