Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-10-25 03:22:34 

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()
0137                                     << " type/RU/mod = " << std::setw(1) << modId.modType() << "/" << std::setw(1)
0138                                     << modId.runit() << "/" << std::setw(2) << modId.module() << std::setw(14)
0139                                     << " R = " << std::setprecision(4) << imod->position().perp() << std::setw(14)
0140                                     << " phi = " << imod->position().phi() << std::setw(14)
0141                                     << " Z = " << imod->position().z();
0142       }
0143     }
0144 
0145     const BoundCylinder& cyl = layer->specificSurface();
0146 
0147     double halfZ = cyl.bounds().length() / 2.;
0148 
0149     // Generate a random point on the cylinder
0150     double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi());
0151     double aZ = CLHEP::RandFlat::shoot(-halfZ, halfZ);
0152     GlobalPoint gp(GlobalPoint::Cylindrical(cyl.radius(), aPhi, aZ));
0153 
0154     // Momentum: 10 GeV, straight from the origin
0155     GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0156 
0157     //FIXME: only negative charge
0158     int charge = -1;
0159 
0160     GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0161     TrajectoryStateOnSurface tsos(gtp, cyl);
0162     LogVerbatim("MTDLayerDump") << "\ntestBTLLayers: at " << std::fixed << std::setw(14) << tsos.globalPosition()
0163                                 << " R=" << std::setw(14) << tsos.globalPosition().perp() << " phi=" << std::setw(14)
0164                                 << tsos.globalPosition().phi() << " Z=" << std::setw(14) << tsos.globalPosition().z()
0165                                 << " p = " << std::setw(14) << tsos.globalMomentum();
0166 
0167     SteppingHelixPropagator prop(field, anyDirection);
0168 
0169     pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, prop, *theEstimator);
0170     LogVerbatim("MTDLayerDump") << "is compatible: " << comp.first << " at: R=" << std::setw(14) << std::setprecision(4)
0171                                 << comp.second.globalPosition().perp() << " phi=" << std::setw(14)
0172                                 << comp.second.globalPosition().phi() << " Z=" << std::setw(14)
0173                                 << comp.second.globalPosition().z();
0174 
0175     vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0176     if (!compDets.empty()) {
0177       LogVerbatim("MTDLayerDump") << "compatibleDets: " << std::setw(14) << compDets.size() << "\n"
0178                                   << "  final state pos: " << std::setw(14) << compDets.front().second.globalPosition()
0179                                   << "\n"
0180                                   << "  det         pos: " << std::setw(14) << compDets.front().first->position()
0181                                   << " id: " << std::hex
0182                                   << BTLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec
0183                                   << "\n"
0184                                   << "  distance " << std::setw(14) << std::setprecision(4)
0185                                   << (tsos.globalPosition() - compDets.front().first->position()).mag();
0186     } else {
0187       LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found";
0188     }
0189   }
0190 }
0191 
0192 void MTDRecoGeometryAnalyzer::testETLLayers(const MTDDetLayerGeometry* geo, const MagneticField* field) const {
0193   const vector<const DetLayer*>& layers = geo->allETLLayers();
0194 
0195   for (const auto& ilay : layers) {
0196     const MTDRingForwardDoubleLayer* layer = static_cast<const MTDRingForwardDoubleLayer*>(ilay);
0197 
0198     LogVerbatim("MTDLayerDump") << std::fixed << "\nETL layer " << std::setw(4) << layer->subDetector()
0199                                 << " rings = " << std::setw(14) << layer->rings().size() << " dets = " << std::setw(14)
0200                                 << layer->basicComponents().size() << " front dets = " << std::setw(14)
0201                                 << layer->frontLayer()->basicComponents().size() << " back dets = " << std::setw(14)
0202                                 << layer->backLayer()->basicComponents().size();
0203 
0204     const BoundDisk& disk = layer->specificSurface();
0205 
0206     // Generate a random point on the disk
0207     double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi());
0208     double aR = CLHEP::RandFlat::shoot(disk.innerRadius(), disk.outerRadius());
0209     GlobalPoint gp(GlobalPoint::Cylindrical(aR, aPhi, disk.position().z()));
0210 
0211     // Momentum: 10 GeV, straight from the origin
0212     GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0213 
0214     //FIXME: only negative charge
0215     int charge = -1;
0216 
0217     GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0218     TrajectoryStateOnSurface tsos(gtp, disk);
0219     LogVerbatim("MTDLayerDump") << "\ntestETLLayers: at " << std::setw(14) << std::setprecision(4)
0220                                 << tsos.globalPosition() << " R=" << std::setw(14) << tsos.globalPosition().perp()
0221                                 << " phi=" << std::setw(14) << tsos.globalPosition().phi() << " Z=" << std::setw(14)
0222                                 << tsos.globalPosition().z() << " p = " << std::setw(14) << tsos.globalMomentum();
0223 
0224     SteppingHelixPropagator prop(field, anyDirection);
0225 
0226     pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, prop, *theEstimator);
0227     LogVerbatim("MTDLayerDump") << "is compatible: " << comp.first << " at: R=" << std::setw(14) << std::setprecision(4)
0228                                 << comp.second.globalPosition().perp() << " phi=" << std::setw(14)
0229                                 << comp.second.globalPosition().phi() << " Z=" << std::setw(14)
0230                                 << comp.second.globalPosition().z();
0231 
0232     vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0233     if (!compDets.empty()) {
0234       LogVerbatim("MTDLayerDump") << "compatibleDets: " << std::setw(14) << compDets.size() << "\n"
0235                                   << "  final state pos: " << std::setw(14) << std::setprecision(4)
0236                                   << compDets.front().second.globalPosition() << "\n"
0237                                   << "  det         pos: " << std::setw(14) << compDets.front().first->position()
0238                                   << " id: " << std::hex
0239                                   << ETLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec
0240                                   << "\n"
0241                                   << "  distance " << std::setw(14)
0242                                   << (tsos.globalPosition() - compDets.front().first->position()).mag();
0243     } else {
0244       if (layer->isCrack(gp)) {
0245         LogVerbatim("MTDLayerDump") << " MTD crack found ";
0246       } else {
0247         LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found in MTD"
0248                                     << " at: R=" << std::setw(14) << std::setprecision(4) << gp.perp()
0249                                     << " phi= " << std::setw(14) << gp.phi().degrees() << " Z= " << std::setw(14)
0250                                     << gp.z();
0251       }
0252     }
0253   }
0254 }
0255 
0256 void MTDRecoGeometryAnalyzer::testETLLayersNew(const MTDDetLayerGeometry* geo, const MagneticField* field) const {
0257   const vector<const DetLayer*>& layers = geo->allETLLayers();
0258 
0259   // dump of ETL layers structure
0260 
0261   for (const auto& ilay : layers) {
0262     const MTDSectorForwardDoubleLayer* layer = static_cast<const MTDSectorForwardDoubleLayer*>(ilay);
0263 
0264     LogVerbatim("MTDLayerDump") << std::fixed << "\nETL layer " << std::setw(4) << layer->subDetector()
0265                                 << " at z = " << std::setw(14) << std::setprecision(4)
0266                                 << layer->surface().position().z() << " sectors = " << std::setw(14)
0267                                 << layer->sectors().size() << " dets = " << std::setw(14)
0268                                 << layer->basicComponents().size() << " front dets = " << std::setw(14)
0269                                 << layer->frontLayer()->basicComponents().size() << " back dets = " << std::setw(14)
0270                                 << layer->backLayer()->basicComponents().size();
0271 
0272     unsigned int isectInd(0);
0273     for (const auto& isector : layer->sectors()) {
0274       isectInd++;
0275       LogVerbatim("MTDLayerDump") << std::fixed << "\nSector " << std::setw(4) << isectInd << "\n" << (*isector);
0276       for (const auto& imod : isector->basicComponents()) {
0277         ETLDetId modId(imod->geographicalId().rawId());
0278         LogVerbatim("MTDLayerDump") << std::fixed << "ETLDetId " << modId.rawId() << " side = " << std::setw(4)
0279                                     << modId.mtdSide() << " Disc/Side/Sector = " << std::setw(4) << modId.nDisc() << " "
0280                                     << std::setw(4) << modId.discSide() << " " << std::setw(4) << modId.sector()
0281                                     << " mod/type = " << std::setw(4) << modId.module() << " " << std::setw(4)
0282                                     << modId.modType() << " pos = " << std::setprecision(4) << imod->position();
0283       }
0284     }
0285   }
0286 
0287   // test propagation through layers
0288 
0289   for (const auto& ilay : layers) {
0290     const MTDSectorForwardDoubleLayer* layer = static_cast<const MTDSectorForwardDoubleLayer*>(ilay);
0291 
0292     const BoundDisk& disk = layer->specificSurface();
0293 
0294     // Generate a random point on the disk
0295     double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi());
0296     double aR = CLHEP::RandFlat::shoot(disk.innerRadius(), disk.outerRadius());
0297     GlobalPoint gp(GlobalPoint::Cylindrical(aR, aPhi, disk.position().z()));
0298 
0299     // Momentum: 10 GeV, straight from the origin
0300     GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0301 
0302     //FIXME: only negative charge
0303     int charge = -1;
0304 
0305     GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0306     TrajectoryStateOnSurface tsos(gtp, disk);
0307     LogVerbatim("MTDLayerDump") << "\ntestETLLayers: at " << std::setprecision(4) << std::fixed << tsos.globalPosition()
0308                                 << " R=" << std::setw(14) << tsos.globalPosition().perp() << " phi=" << std::setw(14)
0309                                 << tsos.globalPosition().phi() << " Z=" << std::setw(14) << tsos.globalPosition().z()
0310                                 << " p = " << tsos.globalMomentum();
0311 
0312     SteppingHelixPropagator prop(field, anyDirection);
0313 
0314     pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, prop, *theEstimator);
0315     LogVerbatim("MTDLayerDump") << std::fixed << "is compatible: " << comp.first << " at: R=" << std::setw(14)
0316                                 << std::setprecision(4) << comp.second.globalPosition().perp()
0317                                 << " phi=" << std::setw(14) << comp.second.globalPosition().phi()
0318                                 << " Z=" << std::setw(14) << comp.second.globalPosition().z();
0319 
0320     vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0321     if (!compDets.empty()) {
0322       LogVerbatim("MTDLayerDump")
0323           << std::fixed << "compatibleDets: " << std::setw(14) << compDets.size() << "\n"
0324           << "  final state pos: " << std::setprecision(4) << compDets.front().second.globalPosition() << "\n"
0325           << "  det         pos: " << compDets.front().first->position() << " id: " << std::hex
0326           << ETLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec << "\n"
0327           << "  distance " << std::setw(14)
0328           << (compDets.front().second.globalPosition() - compDets.front().first->position()).mag();
0329     } else {
0330       if (layer->isCrack(gp)) {
0331         LogVerbatim("MTDLayerDump") << " MTD crack found ";
0332       } else {
0333         LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found in MTD"
0334                                     << " at: R=" << std::setw(14) << std::setprecision(4) << gp.perp()
0335                                     << " phi= " << std::setw(14) << gp.phi().degrees() << " Z= " << std::setw(14)
0336                                     << gp.z();
0337       }
0338     }
0339   }
0340 }
0341 
0342 string MTDRecoGeometryAnalyzer::dumpLayer(const DetLayer* layer) const {
0343   stringstream output;
0344 
0345   const BoundSurface* sur = nullptr;
0346   const BoundCylinder* bc = nullptr;
0347   const BoundDisk* bd = nullptr;
0348 
0349   sur = &(layer->surface());
0350   if ((bc = dynamic_cast<const BoundCylinder*>(sur))) {
0351     output << std::fixed << " subdet = " << std::setw(4) << layer->subDetector() << " Barrel = " << layer->isBarrel()
0352            << " Forward = " << layer->isForward() << "  Cylinder of radius: " << std::setw(14) << std::setprecision(4)
0353            << bc->radius();
0354   } else if ((bd = dynamic_cast<const BoundDisk*>(sur))) {
0355     output << std::fixed << " subdet = " << std::setw(4) << layer->subDetector() << " Barrel = " << layer->isBarrel()
0356            << " Forward = " << layer->isForward() << "  Disk at: " << std::setw(14) << std::setprecision(4)
0357            << bd->position().z();
0358   }
0359   return output.str();
0360 }
0361 
0362 //define this as a plug-in
0363 #include <FWCore/Framework/interface/MakerMacros.h>
0364 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!