Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-12 00:01:21

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackerDigiGeometryAnalyzer
0004 // Class:      TrackerDigiGeometryAnalyzer
0005 //
0006 /**\class TrackerDigiGeometryAnalyzer TrackerDigiGeometryAnalyzer.cc 
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Filippo Ambroglini
0015 //         Created:  Tue Jul 26 08:47:57 CEST 2005
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
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 // class decleration
0049 //
0050 
0051 // #define PRINT(X) edm::LogInfo(X)
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 // ------------ method called to produce the data  ------------
0072 void TrackerDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0073   PRINT("TrackerDigiGeometryAnalyzer") << "Here I am";
0074   //
0075   // get the TrackerGeom
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     //analyseTrapezoidal(*it);
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 "  // << it->geographicalId()
0102                                            << "   isBricked    " << p.isBricked() << "    Rows    " << p.nrows()
0103                                            << "    Columns " << p.ncolumns() << '\n';
0104     } else {
0105       const StripTopology& p = (dynamic_cast<const StripGeomDetType*>((it)))->specificTopology();
0106       PRINT("TrackerDigiGeometryAnalyzer") << " STRIP Det "  // << it->geographicalId()
0107                                            << "    Strips    " << p.nstrips() << '\n';
0108     }
0109   }
0110 }
0111 void TrackerDigiGeometryAnalyzer::analyseTrapezoidal(const GeomDetUnit& det) {
0112   // checkRotation( det);
0113 
0114   const double safety = 0.9999;
0115 
0116   const Bounds& bounds = det.surface().bounds();
0117   const TrapezoidalPlaneBounds* tb = dynamic_cast<const TrapezoidalPlaneBounds*>(&bounds);
0118   if (tb == nullptr)
0119     return;  // not trapezoidal
0120 
0121   checkTopology(det);
0122 
0123   const GlobalPoint& pos = det.position();
0124   double length = tb->length();
0125   double width = tb->width();
0126 
0127   const std::array<const float, 4>& par = tb->parameters();
0128   double top = std::max(par[1], par[0]);
0129   double bot = std::min(par[1], par[0]);
0130 
0131   std::cout << std::endl;
0132   std::cout << "Det at pos " << pos << " has length " << length << " width " << width << " pars ";
0133   for (int i = 0; i < 4; i++)
0134     std::cout << par[i] << ", ";
0135   std::cout << std::endl;
0136 
0137   std::cout << "det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << std::endl;
0138 
0139   //   double outerScale = (pos.perp()+safety*length/2.) / pos.perp();
0140   //   GlobalPoint outerMiddle = GlobalPoint( outerScale*pos.x(), outerScale*pos.y(), pos.z());
0141 
0142   GlobalVector yShift = det.surface().toGlobal(LocalVector(0, safety * length / 2., 0));
0143   GlobalPoint outerMiddle = pos + yShift;
0144   GlobalPoint innerMiddle = pos + (-1. * yShift);
0145   if (outerMiddle.perp() < innerMiddle.perp())
0146     std::swap(outerMiddle, innerMiddle);
0147 
0148   GlobalVector upperShift = det.surface().toGlobal(LocalVector(safety * top, 0, 0));
0149 
0150   GlobalPoint ulc = outerMiddle + upperShift;
0151   GlobalPoint urc = outerMiddle + (-1. * upperShift);
0152   std::cout << "outerMiddle " << outerMiddle << " upperShift " << upperShift << " ulc "
0153             << "(" << ulc.perp() << "," << ulc.phi() << "," << ulc.z() << " urc "
0154             << "(" << urc.perp() << "," << urc.phi() << "," << urc.z() << std::endl;
0155 
0156   std::cout << "upper corners inside bounds? " << tb->inside(det.surface().toLocal(ulc)) << " "
0157             << tb->inside(det.surface().toLocal(urc)) << std::endl;
0158 
0159   //   double innerScale = (pos.perp()-safety*length/2.) / pos.perp();
0160   //   GlobalPoint innerMiddle = GlobalPoint( innerScale*pos.x(), innerScale*pos.y(), pos.z());
0161 
0162   GlobalVector lowerShift = det.surface().toGlobal(LocalVector(safety * bot, 0, 0));
0163 
0164   std::cout << "lower corners inside bounds? " << tb->inside(det.surface().toLocal(innerMiddle + lowerShift)) << " "
0165             << tb->inside(det.surface().toLocal(innerMiddle + (-1. * lowerShift))) << std::endl;
0166 }
0167 
0168 void TrackerDigiGeometryAnalyzer::checkRotation(const GeomDetUnit& det) {
0169   const double eps = std::numeric_limits<float>::epsilon();
0170   static int first = 0;
0171   if (first == 0) {
0172     std::cout << "numeric_limits<float>::epsilon() " << std::numeric_limits<float>::epsilon() << std::endl;
0173     first = 1;
0174   }
0175 
0176   const Surface::RotationType& rot(det.surface().rotation());
0177   GlobalVector a(rot.xx(), rot.xy(), rot.xz());
0178   GlobalVector b(rot.yx(), rot.yy(), rot.yz());
0179   GlobalVector c(rot.zx(), rot.zy(), rot.zz());
0180   GlobalVector cref = a.cross(b);
0181   GlobalVector aref = b.cross(c);
0182   GlobalVector bref = c.cross(a);
0183   if ((a - aref).mag() > eps || (b - bref).mag() > eps || (c - cref).mag() > eps) {
0184     std::cout << " Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag() << ", "
0185               << (c - cref).mag() << " for det at pos " << det.surface().position() << std::endl;
0186   }
0187   if (fabs(a.mag() - 1.) > eps || fabs(b.mag() - 1.) > eps || fabs(c.mag() - 1.) > eps) {
0188     std::cout << " Rotation not good by bector mag: " << (a).mag() << ", " << (b).mag() << ", " << (c).mag()
0189               << " for det at pos " << det.surface().position() << std::endl;
0190   }
0191 }
0192 
0193 void TrackerDigiGeometryAnalyzer::checkTopology(const GeomDetUnit& det) {
0194   const StripTopology& topol = dynamic_cast<const StripTopology&>(det.topology());
0195   const int N = 5;
0196 
0197   std::cout << std::endl << "Topology test along strip 0" << std::endl;
0198   LocalVector stripDir = topol.localPosition(MeasurementPoint(0, 0.5)) - topol.localPosition(MeasurementPoint(0, 0));
0199   std::cout << "StripDir " << stripDir << std::endl;
0200   for (int i = -N; i <= N; i++) {
0201     MeasurementPoint mp(0, float(i) / float(2 * N));
0202     LocalPoint lp = topol.localPosition(mp);
0203     double strlen = topol.localStripLength(lp);
0204     LocalError le = topol.localError(mp, MeasurementError(0.25, 0, 0.25));
0205     LocalError lep = le.rotate(stripDir.y(), stripDir.x());
0206     LocalError lem = le.rotate(stripDir.y(), -stripDir.x());
0207     GlobalPoint gp(det.surface().toGlobal(lp));
0208     std::cout << "gpos (r,phi) (" << gp.perp() << ", " << gp.phi() << ") lpos " << lp << " lerr.x (0,+,-) "
0209               << sqrt(le.xx()) << "," << sqrt(lep.xx()) << "," << sqrt(lem.xx()) << " strlen " << strlen << std::endl;
0210   }
0211 
0212   std::cout << std::endl << "Topology test along middle strip" << std::endl;
0213   float midstrip = topol.strip(LocalPoint(0, 0));
0214   for (int i = -N; i <= N; i++) {
0215     MeasurementPoint mp(midstrip, float(i) / float(2 * N));
0216     LocalPoint lp = topol.localPosition(mp);
0217     double strlen = topol.localStripLength(lp);
0218     LocalError le = topol.localError(mp, MeasurementError(0.25, 0, 0.25));
0219     GlobalPoint gp(det.surface().toGlobal(lp));
0220     std::cout << "gpos (r,phi) (" << gp.perp() << ", " << gp.phi() << ") lpos " << lp << " lerr.x " << sqrt(le.xx())
0221               << " strlen " << strlen << std::endl;
0222   }
0223 }
0224 
0225 std::ostream& TrackerDigiGeometryAnalyzer::cylindrical(std::ostream& os, const GlobalPoint& gp) const {
0226   os << "(" << gp.perp() << "," << gp.phi() << "," << gp.z();
0227   return os;
0228 }
0229 
0230 //define this as a plug-in
0231 DEFINE_FWK_MODULE(TrackerDigiGeometryAnalyzer);