CommonTopologies/test/ValidateRadial.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0003 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0009 #include "Geometry/TrackerGeometryBuilder/interface/ProxyStripTopology.h"
0010 #include "Geometry/CommonTopologies/interface/TkRadialStripTopology.h"
0011 #include "TFile.h"
0012 #include "TProfile.h"
0013
0014 class ValidateRadial : public edm::one::EDAnalyzer<> {
0015 public:
0016   ValidateRadial(const edm::ParameterSet&);
0017   ~ValidateRadial() override;
0018
0019   void beginJob() override {}
0020   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0021   void endJob() override {}
0022
0023 private:
0024   std::vector<const TkRadialStripTopology*> get_list_of_radial_topologies(const edm::Event&, const edm::EventSetup&);
0025   void test_topology(const TkRadialStripTopology*, unsigned);
0026   bool pass_frame_change_test(const TkRadialStripTopology* t, float strip, float stripErr2, bool);
0027   bool EQUAL(const double a, const double b) { return fabs(a - b) < epsilon_; }
0028   const double epsilon_;
0029   TFile* file_;
0030   const bool printOut_;
0031   const bool posOnly_;
0032   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tokTrackerGeometry_;
0033
0034   mutable float maxerrU = 0.;
0035   mutable float maxerrUV = 0.;
0036 };
0037
0038 ValidateRadial::ValidateRadial(const edm::ParameterSet& cfg)
0039     : epsilon_(cfg.getParameter<double>("Epsilon")),
0040       file_(new TFile(cfg.getParameter<std::string>("FileName").c_str(), "RECREATE")),
0041       printOut_(cfg.getParameter<bool>("PrintOut")),
0042       posOnly_(cfg.getParameter<bool>("PosOnly")),
0043       tokTrackerGeometry_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()) {
0044   std::cout << "I'm ALIVE" << std::endl;
0045 }
0046
0047 void ValidateRadial::analyze(const edm::Event& e, const edm::EventSetup& es) {
0048   std::vector<const TkRadialStripTopology*> topologies = get_list_of_radial_topologies(e, es);
0049   for (unsigned i = 0; i < topologies.size(); i++) {
0050     test_topology(topologies[i], i);
0051   }
0052   file_->Close();
0053 }
0054
0055 std::vector<const TkRadialStripTopology*> ValidateRadial::get_list_of_radial_topologies(const edm::Event& e,
0056                                                                                         const edm::EventSetup& es) {
0057   std::vector<const TkRadialStripTopology*> topos;
0058   auto theTrackerGeometry = es.getHandle(tokTrackerGeometry_);
0059   const uint32_t radial_detids[] = {402666125,   //TID r1
0060                                     402668833,   //TID r2
0061                                     402673476,   //TID r3
0062                                     470066725,   //TEC r1
0063                                     470390853,   //TEC r2
0064                                     470114664,   //TEC r3
0065                                     470131344,   //TEC r4
0066                                     470079661,   //TEC r5
0067                                     470049476,   //TEC r6
0068                                     470045428};  //TEC r7
0069   for (unsigned int radial_detid : radial_detids) {
0070     auto g = dynamic_cast<const StripGeomDetUnit*>(theTrackerGeometry->idToDet(radial_detid));
0071     if (!g)
0072       std::cout << "no geom for " << radial_detid << std::endl;
0073     auto const topol = &g->specificTopology();
0074     const TkRadialStripTopology* rt = nullptr;
0075     auto const proxyT = dynamic_cast<const ProxyStripTopology*>(topol);
0076     if (proxyT)
0077       rt = dynamic_cast<const TkRadialStripTopology*>(&(proxyT->specificTopology()));
0078     else
0079       rt = dynamic_cast<const TkRadialStripTopology*>(topol);
0080     if (!rt)
0081       std::cout << "no radial topology for " << radial_detid << std::endl;
0082     else
0083       topos.emplace_back(rt);
0084   }
0085   return topos;
0086 }
0087
0088 #include "Geometry/CommonTopologies/interface/CSCRadialStripTopology.h"
0089
0090 void compare(const TkRadialStripTopology& t,
0091              const CSCRadialStripTopology& ot,
0092              const float strip,
0093              const float stripErr2) {
0094   const LocalPoint lp = t.localPosition(strip);
0095   const LocalError le = t.localError(strip, stripErr2);
0096
0097   const MeasurementPoint mp = t.measurementPosition(lp);
0098   const MeasurementError me = t.measurementError(lp, le);
0099   const float newstrip = t.strip(lp);
0100
0101   const LocalPoint olp = ot.localPosition(strip);
0102   const LocalError ole = ot.localError(strip, stripErr2);
0103
0104   const MeasurementPoint omp = ot.measurementPosition(lp);
0105   const MeasurementError ome = ot.measurementError(lp, le);
0106   const float onewstrip = ot.strip(lp);
0107
0108   if (fabs(lp.x() - olp.x()) > 0.001)
0109     std::cout << "FAILED " << lp.x() << " " << olp.x() << std::endl;
0110   if (fabs(lp.y() - olp.y()) > 0.001)
0111     std::cout << "FAILED " << lp.y() << " " << olp.y() << std::endl;
0112   if (fabs(newstrip - onewstrip) > 0.001)
0113     std::cout << "FAILED " << newstrip << " " << onewstrip << std::endl;
0114   if (fabs(le.xx() - ole.xx()) > 0.001)
0115     std::cout << "FAILED " << le.xx() << " " << ole.xx() << std::endl;
0116   if (fabs(le.yy() - ole.yy()) > 0.001)
0117     std::cout << "FAILED " << le.yy() << " " << ole.xy() << std::endl;
0118   if (fabs(le.xy() - ole.xy()) > 0.001)
0119     std::cout << "FAILED " << le.xy() << " " << ole.yy() << std::endl;
0120   if (fabs(mp.x() - omp.x()) > 0.001)
0121     std::cout << "FAILED " << mp.x() << " " << omp.x() << std::endl;
0122   if (fabs(me.uu() - ome.uu()) > 0.001)
0123     std::cout << "FAILED " << me.uu() << " " << ome.uu() << std::endl;
0124 }
0125
0126 void ValidateRadial::test_topology(const TkRadialStripTopology* t, unsigned i) {
0127   CSCRadialStripTopology oldt(t->nstrips(),
0128                               t->angularWidth(),
0129                               t->detHeight(),
0130                               t->centreToIntersection(),
0131                               t->yAxisOrientation(),
0132                               t->yCentreOfStripPlane());
0133
0134   std::cout << "\nTK\n" << *t << std::endl;
0135   std::cout << "\nCSC\n" << oldt << std::endl;
0136
0137   TProfile prof(("se2limit1" + std::to_string(i)).c_str(),
0138                 "Precision Limit of recoverable strip error (1st order);strip;strip error",
0139                 t->nstrips() / 8,
0140                 0,
0141                 t->nstrips());
0142   TProfile prof2(("se2limit2" + std::to_string(i)).c_str(),
0143                  "Precision Limit of recoverable strip error (2nd order);strip;strip error",
0144                  t->nstrips() / 8,
0145                  0,
0146                  t->nstrips());
0147   for (float strip = 0; strip < t->nstrips(); strip += 0.5) {
0148     for (float stripErr2 = 0.03; stripErr2 > 1e-10; stripErr2 /= 1.1) {
0149       compare(*t, oldt, strip, stripErr2);
0150       if (!pass_frame_change_test(t, strip, stripErr2, false)) {
0151         prof.Fill(strip, sqrt(stripErr2));
0152         break;
0153       }
0154     }
0155     for (float stripErr2 = 0.03; stripErr2 > 1e-10; stripErr2 /= 1.1)
0156       if (!pass_frame_change_test(t, strip, stripErr2, true)) {
0157         prof2.Fill(strip, sqrt(stripErr2));
0158         break;
0159       }
0160   }
0161   prof.Write();
0162   prof2.Write();
0163 }
0164
0165 ValidateRadial::~ValidateRadial() {
0166   std::cout << "ValidateRadial max UU, max UV " << maxerrU << " " << maxerrUV << std::endl;
0167 }
0168
0169 bool ValidateRadial::pass_frame_change_test(const TkRadialStripTopology* t,
0170                                             const float strip,
0171                                             const float stripErr2,
0172                                             const bool secondOrder) {
0173   const LocalPoint lp = t->localPosition(strip);
0174   const LocalError le = t->localError(strip, stripErr2);
0175   const MeasurementPoint mp = t->measurementPosition(lp);
0176   const MeasurementError me = t->measurementError(lp, le);
0177   const float newstrip = t->strip(lp);
0178   const LocalPoint newlp = t->localPosition(mp);
0179   const LocalError newle = t->localError(mp, me);
0180   const MeasurementPoint newmp = t->measurementPosition(newlp);
0181   const MeasurementError newme = t->measurementError(newlp, newle);
0182
0183   maxerrU = std::max(maxerrU, std::abs(me.uu() - stripErr2) / stripErr2);
0184   maxerrUV = std::max(maxerrUV, std::abs(me.uv()));
0185
0186   const bool pass1p = fabs(strip - newstrip) < 0.001 && fabs(strip - mp.x()) < 0.001 && fabs(mp.x() - newstrip) < 0.001;
0187   const bool pass1e =
0188       (fabs(me.uu() - stripErr2) / stripErr2 < epsilon_) && (fabs(me.uv()) < epsilon_) && me.uu() > 0 && me.vv() > 0;
0189   const bool pass2p = fabs(strip - newmp.x()) < 0.001;
0190   const bool pass2e = (fabs(newme.uu() - stripErr2) / stripErr2 < epsilon_) && (fabs(newme.uv()) < epsilon_) &&
0191                       newme.uu() > 0 && newme.vv() > 0;
0192
0193   const bool passp = (secondOrder ? pass2p : pass1p);
0194   const bool passe = (secondOrder ? pass2e : pass1e);
0195
0196   if (printOut_ && ((!passp) || ((!posOnly_) && (!passe))))
0197     std::cout << "FAILED "
0198               << "(" << strip << ", " << newstrip << ", " << mp.x() << ")\t"
0199               << "(" << stripErr2 << ", " << me.uu() << ")\t\t" << (me.uv()) << std::endl;
0200   return passp & passe;
0201 }
0202
0203 DEFINE_FWK_MODULE(ValidateRadial);