Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-07-28 03:10:15 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!