Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-06-28 02:36:36

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/global/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/ESTransientHandle.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/StreamID.h"
0009 
0010 #include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h"
0011 #include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h"
0012 
0013 #include "MagneticField/Engine/interface/MagneticField.h"
0014 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0015 
0016 #include "Geometry/Records/interface/MTDTopologyRcd.h"
0017 #include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
0018 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0019 
0020 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0021 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0022 
0023 #include "RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h"
0024 #include "RecoMTD/DetLayers/interface/MTDDetTray.h"
0025 #include "RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h"
0026 #include "RecoMTD/DetLayers/interface/MTDDetSector.h"
0027 
0028 #include <DataFormats/MuonDetId/interface/CSCDetId.h>
0029 
0030 #include <DataFormats/ForwardDetId/interface/BTLDetId.h>
0031 #include <DataFormats/ForwardDetId/interface/ETLDetId.h>
0032 
0033 #include <memory>
0034 #include <sstream>
0035 
0036 #include "CLHEP/Random/RandFlat.h"
0037 #include "DataFormats/Math/interface/angle.h"
0038 #include "DataFormats/Math/interface/Rounding.h"
0039 
0040 using namespace std;
0041 using namespace edm;
0042 using namespace angle_units;
0043 using namespace cms_rounding;
0044 
0045 class MTDRecoGeometryAnalyzer : public global::EDAnalyzer<> {
0046 public:
0047   MTDRecoGeometryAnalyzer(const ParameterSet& pset);
0048 
0049   void analyze(edm::StreamID, edm::Event const&, edm::EventSetup const&) const override;
0050 
0051   void testBTLLayers(const MTDDetLayerGeometry*, const MagneticField* field) const;
0052   void testETLLayersNew(const MTDDetLayerGeometry*, const MagneticField* field) const;
0053 
0054   string dumpLayer(const DetLayer* layer) const;
0055 
0056 private:
0057   std::unique_ptr<MeasurementEstimator> theEstimator;
0058 
0059   inline std::string fround(const double in, const size_t prec) const {
0060     std::stringstream ss;
0061     ss << std::setprecision(prec) << std::fixed << std::setw(14) << roundIfNear0(in);
0062     return ss.str();
0063   }
0064 
0065   inline std::string fvecround(const GlobalPoint vecin, const size_t prec) const {
0066     std::stringstream ss;
0067     ss << std::setprecision(prec) << std::fixed << std::setw(14) << roundVecIfNear0(vecin);
0068     return ss.str();
0069   }
0070 
0071   inline std::string fvecround(const GlobalVector vecin, const size_t prec) const {
0072     std::stringstream ss;
0073     ss << std::setprecision(prec) << std::fixed << std::setw(14) << roundVecIfNear0(vecin);
0074     return ss.str();
0075   }
0076 
0077   const edm::ESInputTag tag_;
0078   edm::ESGetToken<MTDDetLayerGeometry, MTDRecoGeometryRecord> geomToken_;
0079   edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0080   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magfieldToken_;
0081 };
0082 
0083 MTDRecoGeometryAnalyzer::MTDRecoGeometryAnalyzer(const ParameterSet& iConfig) : tag_(edm::ESInputTag{"", ""}) {
0084   geomToken_ = esConsumes<MTDDetLayerGeometry, MTDRecoGeometryRecord>(tag_);
0085   mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>(tag_);
0086   magfieldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>(tag_);
0087 
0088   float theMaxChi2 = 25.;
0089   float theNSigma = 3.;
0090   theEstimator = std::make_unique<Chi2MeasurementEstimator>(theMaxChi2, theNSigma);
0091 }
0092 
0093 void MTDRecoGeometryAnalyzer::analyze(edm::StreamID, edm::Event const&, edm::EventSetup const& es) const {
0094   auto geo = es.getTransientHandle(geomToken_);
0095   auto mtdtopo = es.getTransientHandle(mtdtopoToken_);
0096   auto magfield = es.getTransientHandle(magfieldToken_);
0097 
0098   // Some printouts
0099 
0100   LogVerbatim("MTDLayerDump") << "\n*** allBTLLayers(): " << std::fixed << std::setw(14) << geo->allBTLLayers().size();
0101   for (auto dl = geo->allBTLLayers().begin(); dl != geo->allBTLLayers().end(); ++dl) {
0102     LogVerbatim("MTDLayerDumpFull") << "  " << static_cast<int>(dl - geo->allBTLLayers().begin()) << " "
0103                                     << dumpLayer(*dl);
0104     LogVerbatim("MTDLayerDump") << "  " << static_cast<int>(dl - geo->allBTLLayers().begin()) << " " << dumpLayer(*dl);
0105   }
0106 
0107   LogVerbatim("MTDLayerDump") << "\n*** allETLLayers(): " << std::fixed << std::setw(14) << geo->allETLLayers().size();
0108   for (auto dl = geo->allETLLayers().begin(); dl != geo->allETLLayers().end(); ++dl) {
0109     LogVerbatim("MTDLayerDumpFull") << "  " << static_cast<int>(dl - geo->allETLLayers().begin()) << " "
0110                                     << dumpLayer(*dl);
0111     LogVerbatim("MTDLayerDump") << "  " << static_cast<int>(dl - geo->allETLLayers().begin()) << " " << dumpLayer(*dl);
0112   }
0113 
0114   LogVerbatim("MTDLayerDumpFull") << "\n*** allForwardLayers(): " << std::fixed << std::setw(14)
0115                                   << geo->allForwardLayers().size();
0116   LogVerbatim("MTDLayerDump") << "\n*** allForwardLayers(): " << std::fixed << std::setw(14)
0117                               << geo->allForwardLayers().size();
0118   for (auto dl = geo->allForwardLayers().begin(); dl != geo->allForwardLayers().end(); ++dl) {
0119     LogVerbatim("MTDLayerDumpFull") << "  " << static_cast<int>(dl - geo->allForwardLayers().begin()) << " "
0120                                     << dumpLayer(*dl);
0121     LogVerbatim("MTDLayerDump") << "  " << static_cast<int>(dl - geo->allForwardLayers().begin()) << " "
0122                                 << dumpLayer(*dl);
0123   }
0124 
0125   LogVerbatim("MTDLayerDumpFull") << "\n*** allBackwardLayers(): " << std::fixed << std::setw(14)
0126                                   << geo->allBackwardLayers().size();
0127   LogVerbatim("MTDLayerDump") << "\n*** allBackwardLayers(): " << std::fixed << std::setw(14)
0128                               << geo->allBackwardLayers().size();
0129   for (auto dl = geo->allBackwardLayers().begin(); dl != geo->allBackwardLayers().end(); ++dl) {
0130     LogVerbatim("MTDLayerDumpFull") << "  " << static_cast<int>(dl - geo->allBackwardLayers().begin()) << " "
0131                                     << dumpLayer(*dl);
0132     LogVerbatim("MTDLayerDump") << "  " << static_cast<int>(dl - geo->allBackwardLayers().begin()) << " "
0133                                 << dumpLayer(*dl);
0134   }
0135 
0136   LogVerbatim("MTDLayerDumpFull") << "\n*** allLayers(): " << std::fixed << std::setw(14) << geo->allLayers().size();
0137   LogVerbatim("MTDLayerDump") << "\n*** allLayers(): " << std::fixed << std::setw(14) << geo->allLayers().size();
0138   for (auto dl = geo->allLayers().begin(); dl != geo->allLayers().end(); ++dl) {
0139     LogVerbatim("MTDLayerDumpFull") << "  " << static_cast<int>(dl - geo->allLayers().begin()) << " " << dumpLayer(*dl);
0140     LogVerbatim("MTDLayerDump") << "  " << static_cast<int>(dl - geo->allLayers().begin()) << " " << dumpLayer(*dl);
0141   }
0142 
0143   testBTLLayers(geo.product(), magfield.product());
0144   testETLLayersNew(geo.product(), magfield.product());
0145 }
0146 
0147 void MTDRecoGeometryAnalyzer::testBTLLayers(const MTDDetLayerGeometry* geo, const MagneticField* field) const {
0148   const vector<const DetLayer*>& layers = geo->allBTLLayers();
0149 
0150   for (const auto& ilay : layers) {
0151     const MTDTrayBarrelLayer* layer = static_cast<const MTDTrayBarrelLayer*>(ilay);
0152 
0153     LogVerbatim("MTDLayerDumpFull") << std::fixed << "\nBTL layer " << std::setw(4) << layer->subDetector()
0154                                     << " rods = " << std::setw(14) << layer->rods().size()
0155                                     << " dets = " << std::setw(14) << layer->basicComponents().size();
0156     LogVerbatim("MTDLayerDump") << std::fixed << "\nBTL layer " << std::setw(4) << layer->subDetector()
0157                                 << " rods = " << std::setw(14) << layer->rods().size() << " dets = " << std::setw(14)
0158                                 << layer->basicComponents().size();
0159 
0160     unsigned int irodInd(0);
0161     for (const auto& irod : layer->rods()) {
0162       irodInd++;
0163       LogVerbatim("MTDLayerDumpFull") << std::fixed << "\nRod " << irodInd
0164                                       << " dets = " << irod->basicComponents().size() << "\n";
0165       LogVerbatim("MTDLayerDump") << std::fixed << "\nRod " << irodInd << " dets = " << irod->basicComponents().size()
0166                                   << "\n";
0167       for (const auto& imod : irod->basicComponents()) {
0168         BTLDetId modId(imod->geographicalId().rawId());
0169         LogVerbatim("MTDLayerDumpFull") << std::fixed << "BTLDetId " << modId.rawId() << " side = " << std::setw(4)
0170                                         << modId.mtdSide() << " rod = " << modId.mtdRR()
0171                                         << " type/RU/mod = " << std::setw(1) << modId.modType() << "/" << std::setw(1)
0172                                         << modId.runit() << "/" << std::setw(2) << modId.module()
0173                                         << " R = " << fround(imod->position().perp(), 4)
0174                                         << " phi = " << fround(imod->position().phi(), 4)
0175                                         << " Z = " << fround(imod->position().z(), 4);
0176         LogVerbatim("MTDLayerDump") << std::fixed << "BTLDetId " << modId.rawId() << " side = " << std::setw(4)
0177                                     << modId.mtdSide() << " rod = " << modId.mtdRR()
0178                                     << " type/RU/mod = " << std::setw(1) << modId.modType() << "/" << std::setw(1)
0179                                     << modId.runit() << "/" << std::setw(2) << modId.module()
0180                                     << " R = " << fround(imod->position().perp(), 2)
0181                                     << " phi = " << fround(imod->position().phi(), 2)
0182                                     << " Z = " << fround(imod->position().z(), 2);
0183       }
0184     }
0185 
0186     const BoundCylinder& cyl = layer->specificSurface();
0187 
0188     double halfZ = cyl.bounds().length() / 2.;
0189 
0190     // Generate a random point on the cylinder
0191     double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi());
0192     double aZ = CLHEP::RandFlat::shoot(-halfZ, halfZ);
0193     GlobalPoint gp(GlobalPoint::Cylindrical(cyl.radius(), aPhi, aZ));
0194 
0195     // Momentum: 10 GeV, straight from the origin
0196     GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0197 
0198     //FIXME: only negative charge
0199     int charge = -1;
0200 
0201     GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0202     TrajectoryStateOnSurface tsos(gtp, cyl);
0203     LogVerbatim("MTDLayerDumpFull") << "\ntestBTLLayers: at " << fvecround(tsos.globalPosition(), 4)
0204                                     << " R=" << fround(tsos.globalPosition().perp(), 4)
0205                                     << " phi=" << fround(tsos.globalPosition().phi(), 4)
0206                                     << " Z=" << fround(tsos.globalPosition().z(), 4)
0207                                     << " p = " << fvecround(tsos.globalMomentum(), 4);
0208     LogVerbatim("MTDLayerDump") << "\ntestBTLLayers: at " << fvecround(tsos.globalPosition(), 2)
0209                                 << " R=" << fround(tsos.globalPosition().perp(), 2)
0210                                 << " phi=" << fround(tsos.globalPosition().phi(), 2)
0211                                 << " Z=" << fround(tsos.globalPosition().z(), 2)
0212                                 << " p = " << fvecround(tsos.globalMomentum(), 2);
0213 
0214     SteppingHelixPropagator prop(field, anyDirection);
0215 
0216     pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, prop, *theEstimator);
0217     LogVerbatim("MTDLayerDumpFull") << "is compatible: " << comp.first
0218                                     << " at: R=" << fround(comp.second.globalPosition().perp(), 4)
0219                                     << " phi=" << fround(comp.second.globalPosition().phi(), 4)
0220                                     << " Z=" << fround(comp.second.globalPosition().z(), 4);
0221     LogVerbatim("MTDLayerDump") << "is compatible: " << comp.first
0222                                 << " at: R=" << fround(comp.second.globalPosition().perp(), 2)
0223                                 << " phi=" << fround(comp.second.globalPosition().phi(), 2)
0224                                 << " Z=" << fround(comp.second.globalPosition().z(), 2);
0225 
0226     vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0227     if (!compDets.empty()) {
0228       LogVerbatim("MTDLayerDumpFull") << "compatibleDets: " << std::setw(14) << compDets.size() << "\n"
0229                                       << "  final state pos: " << fvecround(compDets.front().second.globalPosition(), 4)
0230                                       << "\n"
0231                                       << "  det         pos: " << fvecround(compDets.front().first->position(), 4)
0232                                       << " id: " << std::hex
0233                                       << BTLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec
0234                                       << "\n"
0235                                       << "  distance "
0236                                       << fround((tsos.globalPosition() - compDets.front().first->position()).mag(), 2);
0237       LogVerbatim("MTDLayerDump") << "compatibleDets: " << std::setw(14) << compDets.size() << "\n"
0238                                   << "  final state pos: " << fvecround(compDets.front().second.globalPosition(), 2)
0239                                   << "\n"
0240                                   << "  det         pos: " << fvecround(compDets.front().first->position(), 2)
0241                                   << " id: " << std::hex
0242                                   << BTLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec
0243                                   << "\n"
0244                                   << "  distance "
0245                                   << fround((tsos.globalPosition() - compDets.front().first->position()).mag(), 2);
0246     } else {
0247       LogVerbatim("MTDLayerDumpFull") << " ERROR : no compatible det found";
0248       LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found";
0249     }
0250 
0251     // scan in phi at the given z
0252     LogVerbatim("MTDLayerDumpFull") << "\nBTL phi scan at Z = " << fround(aZ, 4) << "\n";
0253     LogVerbatim("MTDLayerDump") << "\nBTL phi scan at Z = " << fround(aZ, 2) << "\n";
0254     aPhi = (int)(-piRadians * 1000) / 1000.;
0255     double dPhi = 0.005;
0256     uint32_t nTot(0);
0257     uint32_t nComp(0);
0258     while (aPhi <= piRadians) {
0259       nTot++;
0260       GlobalPoint gp(GlobalPoint::Cylindrical(cyl.radius(), aPhi, aZ));
0261       GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0262       GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0263       TrajectoryStateOnSurface tsos(gtp, cyl);
0264       SteppingHelixPropagator prop(field, anyDirection);
0265       vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0266       std::stringstream ss;
0267       if (!compDets.empty()) {
0268         nComp++;
0269         for (const auto& dets : compDets) {
0270           ss << " " << BTLDetId(dets.first->geographicalId().rawId()).rawId();
0271         }
0272       }
0273       LogVerbatim("MTDLayerDumpFull") << "BTL scan at phi = " << std::fixed << std::setw(5) << aPhi
0274                                       << " compatible dets = " << std::setw(14) << compDets.size() << ss.str();
0275       LogVerbatim("MTDLayerDump") << "BTL scan at phi = " << std::fixed << std::setw(5) << aPhi
0276                                   << " compatible dets = " << std::setw(14) << compDets.size() << ss.str();
0277       aPhi += dPhi;
0278     }
0279     LogVerbatim("MTDLayerDumpFull") << "\nBTL scan total points = " << nTot << " compatible = " << nComp
0280                                     << " fraction = " << double(nComp) / double(nTot);
0281     LogVerbatim("MTDLayerDump") << "\nBTL scan total points = " << nTot << " compatible = " << nComp
0282                                 << " fraction = " << double(nComp) / double(nTot);
0283   }
0284 }
0285 
0286 void MTDRecoGeometryAnalyzer::testETLLayersNew(const MTDDetLayerGeometry* geo, const MagneticField* field) const {
0287   const vector<const DetLayer*>& layers = geo->allETLLayers();
0288 
0289   // dump of ETL layers structure
0290 
0291   for (const auto& ilay : layers) {
0292     const MTDSectorForwardDoubleLayer* layer = static_cast<const MTDSectorForwardDoubleLayer*>(ilay);
0293 
0294     LogVerbatim("MTDLayerDumpFull") << std::fixed << "\nETL layer " << std::setw(4) << layer->subDetector()
0295                                     << " at z = " << fround(layer->surface().position().z(), 4)
0296                                     << " sectors = " << std::setw(14) << layer->sectors().size()
0297                                     << " dets = " << std::setw(14) << layer->basicComponents().size()
0298                                     << " front dets = " << std::setw(14)
0299                                     << layer->frontLayer()->basicComponents().size() << " back dets = " << std::setw(14)
0300                                     << layer->backLayer()->basicComponents().size();
0301     LogVerbatim("MTDLayerDump") << std::fixed << "\nETL layer " << std::setw(4) << layer->subDetector()
0302                                 << " at z = " << fround(layer->surface().position().z(), 2)
0303                                 << " sectors = " << std::setw(14) << layer->sectors().size()
0304                                 << " dets = " << std::setw(14) << layer->basicComponents().size()
0305                                 << " front dets = " << std::setw(14) << layer->frontLayer()->basicComponents().size()
0306                                 << " back dets = " << std::setw(14) << layer->backLayer()->basicComponents().size();
0307 
0308     unsigned int isectInd(0);
0309     for (const auto& isector : layer->sectors()) {
0310       isectInd++;
0311       LogVerbatim("MTDLayerDumpFull") << std::fixed << "\nSector " << std::setw(4) << isectInd << "\n" << (*isector);
0312       LogVerbatim("MTDLayerDump") << std::fixed << "\nSector " << std::setw(4) << isectInd << "\n" << (*isector);
0313       for (const auto& imod : isector->basicComponents()) {
0314         ETLDetId modId(imod->geographicalId().rawId());
0315         LogVerbatim("MTDLayerDumpFull") << std::fixed << "ETLDetId " << modId.rawId() << " side = " << std::setw(4)
0316                                         << modId.mtdSide() << " Disc/Side/Sector = " << std::setw(4) << modId.nDisc()
0317                                         << " " << std::setw(4) << modId.discSide() << " " << std::setw(4)
0318                                         << modId.sector() << " mod/type = " << std::setw(4) << modId.module() << " "
0319                                         << std::setw(4) << modId.modType()
0320                                         << " pos = " << fvecround(imod->position(), 4);
0321         LogVerbatim("MTDLayerDump") << std::fixed << "ETLDetId " << modId.rawId() << " side = " << std::setw(4)
0322                                     << modId.mtdSide() << " Disc/Side/Sector = " << std::setw(4) << modId.nDisc() << " "
0323                                     << std::setw(4) << modId.discSide() << " " << std::setw(4) << modId.sector()
0324                                     << " mod/type = " << std::setw(4) << modId.module() << " " << std::setw(4)
0325                                     << modId.modType() << " pos = " << fvecround(imod->position(), 2);
0326       }
0327     }
0328   }
0329 
0330   // test propagation through layers
0331 
0332   for (const auto& ilay : layers) {
0333     const MTDSectorForwardDoubleLayer* layer = static_cast<const MTDSectorForwardDoubleLayer*>(ilay);
0334 
0335     const BoundDisk& disk = layer->specificSurface();
0336 
0337     // Generate a random point on the disk
0338     double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi());
0339     double aR = CLHEP::RandFlat::shoot(disk.innerRadius(), disk.outerRadius());
0340     GlobalPoint gp(GlobalPoint::Cylindrical(aR, aPhi, disk.position().z()));
0341 
0342     // Momentum: 10 GeV, straight from the origin
0343     GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0344 
0345     //FIXME: only negative charge
0346     int charge = -1;
0347 
0348     GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0349     TrajectoryStateOnSurface tsos(gtp, disk);
0350     LogVerbatim("MTDLayerDumpFull") << "\ntestETLLayers: at " << fvecround(tsos.globalPosition(), 4)
0351                                     << " R=" << fround(tsos.globalPosition().perp(), 4)
0352                                     << " phi=" << fround(tsos.globalPosition().phi(), 4)
0353                                     << " Z=" << fround(tsos.globalPosition().z(), 4)
0354                                     << " p = " << fvecround(tsos.globalMomentum(), 4);
0355     LogVerbatim("MTDLayerDump") << "\ntestETLLayers: at " << fvecround(tsos.globalPosition(), 2)
0356                                 << " R=" << fround(tsos.globalPosition().perp(), 2)
0357                                 << " phi=" << fround(tsos.globalPosition().phi(), 2)
0358                                 << " Z=" << fround(tsos.globalPosition().z(), 2)
0359                                 << " p = " << fvecround(tsos.globalMomentum(), 2);
0360 
0361     SteppingHelixPropagator prop(field, anyDirection);
0362 
0363     pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, prop, *theEstimator);
0364     LogVerbatim("MTDLayerDumpFull") << std::fixed << "is compatible: " << comp.first
0365                                     << " at: R=" << fround(comp.second.globalPosition().perp(), 4)
0366                                     << " phi=" << fround(comp.second.globalPosition().phi(), 4)
0367                                     << " Z=" << fround(comp.second.globalPosition().z(), 4);
0368     LogVerbatim("MTDLayerDump") << std::fixed << "is compatible: " << comp.first
0369                                 << " at: R=" << fround(comp.second.globalPosition().perp(), 2)
0370                                 << " phi=" << fround(comp.second.globalPosition().phi(), 2)
0371                                 << " Z=" << fround(comp.second.globalPosition().z(), 2);
0372 
0373     vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0374     if (!compDets.empty()) {
0375       LogVerbatim("MTDLayerDumpFull")
0376           << std::fixed << "compatibleDets: " << std::setw(14) << compDets.size() << "\n"
0377           << "  final state pos: " << fvecround(compDets.front().second.globalPosition(), 4) << "\n"
0378           << "  det         pos: " << fvecround(compDets.front().first->position(), 4) << " id: " << std::hex
0379           << ETLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec << "\n"
0380           << "  distance "
0381           << fround((compDets.front().second.globalPosition() - compDets.front().first->position()).mag(), 4);
0382       LogVerbatim("MTDLayerDump")
0383           << std::fixed << "compatibleDets: " << std::setw(14) << compDets.size() << "\n"
0384           << "  final state pos: " << fvecround(compDets.front().second.globalPosition(), 2) << "\n"
0385           << "  det         pos: " << fvecround(compDets.front().first->position(), 2) << " id: " << std::hex
0386           << ETLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec << "\n"
0387           << "  distance "
0388           << fround((compDets.front().second.globalPosition() - compDets.front().first->position()).mag(), 2);
0389     } else {
0390       if (layer->isCrack(gp)) {
0391         LogVerbatim("MTDLayerDumpFull") << " MTD crack found ";
0392         LogVerbatim("MTDLayerDump") << " MTD crack found ";
0393       } else {
0394         LogVerbatim("MTDLayerDumpFull") << " ERROR : no compatible det found in MTD"
0395                                         << " at: R=" << fround(gp.perp(), 4)
0396                                         << " phi= " << fround(gp.phi().degrees(), 4) << " Z= " << fround(gp.z(), 4);
0397         LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found in MTD"
0398                                     << " at: R=" << fround(gp.perp(), 4) << " phi= " << fround(gp.phi().degrees(), 2)
0399                                     << " Z= " << fround(gp.z(), 2);
0400       }
0401     }
0402   }
0403 }
0404 
0405 string MTDRecoGeometryAnalyzer::dumpLayer(const DetLayer* layer) const {
0406   stringstream output;
0407 
0408   const BoundSurface* sur = nullptr;
0409   const BoundCylinder* bc = nullptr;
0410   const BoundDisk* bd = nullptr;
0411 
0412   sur = &(layer->surface());
0413   if ((bc = dynamic_cast<const BoundCylinder*>(sur))) {
0414     output << std::fixed << " subdet = " << std::setw(4) << layer->subDetector() << " Barrel = " << layer->isBarrel()
0415            << " Forward = " << layer->isForward() << "  Cylinder of radius: " << fround(bc->radius(), 4);
0416   } else if ((bd = dynamic_cast<const BoundDisk*>(sur))) {
0417     output << std::fixed << " subdet = " << std::setw(4) << layer->subDetector() << " Barrel = " << layer->isBarrel()
0418            << " Forward = " << layer->isForward() << "  Disk at: " << fround(bd->position().z(), 4);
0419   }
0420   return output.str();
0421 }
0422 
0423 //define this as a plug-in
0424 #include <FWCore/Framework/interface/MakerMacros.h>
0425 DEFINE_FWK_MODULE(MTDRecoGeometryAnalyzer);