DetLayers/test/MuonRecoGeometryAnalyzer.cc

0001 /** \file
0002  *
0003  */
0004
0005 #include <FWCore/Framework/interface/Frameworkfwd.h>
0006 #include <FWCore/Framework/interface/one/EDAnalyzer.h>
0007 #include <FWCore/Framework/interface/Event.h>
0008 #include <FWCore/Framework/interface/EventSetup.h>
0009 #include <FWCore/Framework/interface/ESHandle.h>
0010 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0011
0012 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
0013 #include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h"
0014
0015 #include "MagneticField/Engine/interface/MagneticField.h"
0016 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0017
0018 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0019 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0020
0021 #include "RecoMuon/DetLayers/interface/MuRodBarrelLayer.h"
0022 #include "RecoMuon/DetLayers/interface/MuDetRod.h"
0023 #include "RecoMuon/DetLayers/interface/MuRingForwardDoubleLayer.h"
0024 #include "RecoMuon/DetLayers/interface/MuDetRing.h"
0025
0026 #include <DataFormats/MuonDetId/interface/DTWireId.h>
0027 #include <DataFormats/MuonDetId/interface/CSCDetId.h>
0028
0029 #include <sstream>
0030
0031 #include "CLHEP/Random/RandFlat.h"
0032
0033 using namespace std;
0034 using namespace edm;
0035
0036 class MuonRecoGeometryAnalyzer : public one::EDAnalyzer<> {
0037 public:
0038   MuonRecoGeometryAnalyzer(const ParameterSet& pset);
0039
0040   void analyze(const Event& ev, const EventSetup& es) override;
0041
0042   void testDTLayers(const MuonDetLayerGeometry*, const MagneticField* field);
0043   void testCSCLayers(const MuonDetLayerGeometry*, const MagneticField* field);
0044
0045   string dumpLayer(const DetLayer* layer) const;
0046
0047 private:
0048   MeasurementEstimator* theEstimator;
0049   ESGetToken<MuonDetLayerGeometry, MuonRecoGeometryRecord> theGeomToken;
0050   ESGetToken<MagneticField, IdealMagneticFieldRecord> theMagfieldToken;
0051 };
0052
0053 MuonRecoGeometryAnalyzer::MuonRecoGeometryAnalyzer(const ParameterSet& iConfig) {
0054   float theMaxChi2 = 25.;
0055   float theNSigma = 3.;
0056   theEstimator = new Chi2MeasurementEstimator(theMaxChi2, theNSigma);
0057   theGeomToken = esConsumes();
0058   theMagfieldToken = esConsumes();
0059 }
0060
0061 void MuonRecoGeometryAnalyzer::analyze(const Event& ev, const EventSetup& es) {
0062   ESHandle<MuonDetLayerGeometry> geo = es.getHandle(theGeomToken);
0063   ESHandle<MagneticField> magfield = es.getHandle(theMagfieldToken);
0064   // Some printouts
0065
0066   cout << "*** allDTLayers(): " << geo->allDTLayers().size() << endl;
0067   for (auto dl = geo->allDTLayers().begin(); dl != geo->allDTLayers().end(); ++dl) {
0068     cout << "  " << (int)(dl - geo->allDTLayers().begin()) << " " << dumpLayer(*dl);
0069   }
0070   cout << endl << endl;
0071
0072   cout << "*** allCSCLayers(): " << geo->allCSCLayers().size() << endl;
0073   for (auto dl = geo->allCSCLayers().begin(); dl != geo->allCSCLayers().end(); ++dl) {
0074     cout << "  " << (int)(dl - geo->allCSCLayers().begin()) << " " << dumpLayer(*dl);
0075   }
0076   cout << endl << endl;
0077
0078   cout << "*** forwardCSCLayers(): " << geo->forwardCSCLayers().size() << endl;
0079   for (auto dl = geo->forwardCSCLayers().begin(); dl != geo->forwardCSCLayers().end(); ++dl) {
0080     cout << "  " << (int)(dl - geo->forwardCSCLayers().begin()) << " " << dumpLayer(*dl);
0081   }
0082   cout << endl << endl;
0083
0084   cout << "*** backwardCSCLayers(): " << geo->backwardCSCLayers().size() << endl;
0085   for (auto dl = geo->backwardCSCLayers().begin(); dl != geo->backwardCSCLayers().end(); ++dl) {
0086     cout << "  " << (int)(dl - geo->backwardCSCLayers().begin()) << " " << dumpLayer(*dl);
0087   }
0088   cout << endl << endl;
0089
0090   cout << "*** allRPCLayers(): " << geo->allRPCLayers().size() << endl;
0091   for (auto dl = geo->allRPCLayers().begin(); dl != geo->allRPCLayers().end(); ++dl) {
0092     cout << "  " << (int)(dl - geo->allRPCLayers().begin()) << " " << dumpLayer(*dl);
0093   }
0094   cout << endl << endl;
0095
0096   cout << "*** endcapRPCLayers(): " << geo->endcapRPCLayers().size() << endl;
0097   for (auto dl = geo->endcapRPCLayers().begin(); dl != geo->endcapRPCLayers().end(); ++dl) {
0098     cout << "  " << (int)(dl - geo->endcapRPCLayers().begin()) << " " << dumpLayer(*dl);
0099   }
0100   cout << endl << endl;
0101
0102   cout << "*** barrelRPCLayers(): " << geo->barrelRPCLayers().size() << endl;
0103   for (auto dl = geo->barrelRPCLayers().begin(); dl != geo->barrelRPCLayers().end(); ++dl) {
0104     cout << "  " << (int)(dl - geo->barrelRPCLayers().begin()) << " " << dumpLayer(*dl);
0105   }
0106   cout << endl << endl;
0107
0108   cout << "*** forwardRPCLayers(): " << geo->forwardRPCLayers().size() << endl;
0109   for (auto dl = geo->forwardRPCLayers().begin(); dl != geo->forwardRPCLayers().end(); ++dl) {
0110     cout << "  " << (int)(dl - geo->forwardRPCLayers().begin()) << " " << dumpLayer(*dl);
0111   }
0112   cout << endl << endl;
0113
0114   cout << "*** backwardRPCLayers(): " << geo->backwardRPCLayers().size() << endl;
0115   for (auto dl = geo->backwardRPCLayers().begin(); dl != geo->backwardRPCLayers().end(); ++dl) {
0116     cout << "  " << (int)(dl - geo->backwardRPCLayers().begin()) << " " << dumpLayer(*dl);
0117   }
0118   cout << endl << endl;
0119
0120   cout << "*** allBarrelLayers(): " << geo->allBarrelLayers().size() << endl;
0121   for (auto dl = geo->allBarrelLayers().begin(); dl != geo->allBarrelLayers().end(); ++dl) {
0122     cout << "  " << (int)(dl - geo->allBarrelLayers().begin()) << " " << dumpLayer(*dl);
0123   }
0124   cout << endl << endl;
0125
0126   cout << "*** allEndcapLayers(): " << geo->allEndcapLayers().size() << endl;
0127   for (auto dl = geo->allEndcapLayers().begin(); dl != geo->allEndcapLayers().end(); ++dl) {
0128     cout << "  " << (int)(dl - geo->allEndcapLayers().begin()) << " " << dumpLayer(*dl);
0129   }
0130   cout << endl << endl;
0131
0132   cout << "*** allForwardLayers(): " << geo->allForwardLayers().size() << endl;
0133   for (auto dl = geo->allForwardLayers().begin(); dl != geo->allForwardLayers().end(); ++dl) {
0134     cout << "  " << (int)(dl - geo->allForwardLayers().begin()) << " " << dumpLayer(*dl);
0135   }
0136   cout << endl << endl;
0137
0138   cout << "*** allBackwardLayers(): " << geo->allBackwardLayers().size() << endl;
0139   for (auto dl = geo->allBackwardLayers().begin(); dl != geo->allBackwardLayers().end(); ++dl) {
0140     cout << "  " << (int)(dl - geo->allBackwardLayers().begin()) << " " << dumpLayer(*dl);
0141   }
0142   cout << endl << endl;
0143
0144   cout << "*** allLayers(): " << geo->allLayers().size() << endl;
0145   for (auto dl = geo->allLayers().begin(); dl != geo->allLayers().end(); ++dl) {
0146     cout << "  " << (int)(dl - geo->allLayers().begin()) << " " << dumpLayer(*dl);
0147   }
0148   cout << endl << endl;
0149
0150   testDTLayers(geo.product(), magfield.product());
0151   testCSCLayers(geo.product(), magfield.product());
0152 }
0153
0154 void MuonRecoGeometryAnalyzer::testDTLayers(const MuonDetLayerGeometry* geo, const MagneticField* field) {
0155   const vector<const DetLayer*>& layers = geo->allDTLayers();
0156
0157   for (auto ilay = layers.begin(); ilay != layers.end(); ++ilay) {
0158     const MuRodBarrelLayer* layer = (const MuRodBarrelLayer*)(*ilay);
0159
0160     const BoundCylinder& cyl = layer->specificSurface();
0161
0162     double halfZ = cyl.bounds().length() / 2.;
0163
0164     // Generate a random point on the cylinder
0165     double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi());
0166     double aZ = CLHEP::RandFlat::shoot(-halfZ, halfZ);
0167     GlobalPoint gp(GlobalPoint::Cylindrical(cyl.radius(), aPhi, aZ));
0168
0169     // Momentum: 10 GeV, straight from the origin
0170     GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0171
0172     //FIXME: only negative charge
0173     int charge = -1;
0174
0175     GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0176     TrajectoryStateOnSurface tsos(gtp, cyl);
0177     cout << "testDTLayers: at " << tsos.globalPosition() << " R=" << tsos.globalPosition().perp()
0178          << " phi=" << tsos.globalPosition().phi() << " Z=" << tsos.globalPosition().z()
0179          << " p = " << tsos.globalMomentum() << endl;
0180
0181     SteppingHelixPropagator prop(field, anyDirection);
0182
0183     pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, prop, *theEstimator);
0184     cout << "is compatible: " << comp.first << " at: R=" << comp.second.globalPosition().perp()
0185          << " phi=" << comp.second.globalPosition().phi() << " Z=" << comp.second.globalPosition().z() << endl;
0186
0187     vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0188     if (compDets.size()) {
0189       cout << "compatibleDets: " << compDets.size() << endl
0190
0191            << "  final state pos: " << compDets.front().second.globalPosition() << endl
0192            << "  det         pos: " << compDets.front().first->position()
0193            << " id: " << DTWireId(compDets.front().first->geographicalId().rawId()) << endl
0194            << "  distance " << (tsos.globalPosition() - compDets.front().first->position()).mag()
0195
0196            << endl
0197            << endl;
0198     } else {
0199       cout << " ERROR : no compatible det found" << endl;
0200     }
0201   }
0202 }
0203
0204 void MuonRecoGeometryAnalyzer::testCSCLayers(const MuonDetLayerGeometry* geo, const MagneticField* field) {
0205   const vector<const DetLayer*>& layers = geo->allCSCLayers();
0206
0207   for (auto ilay = layers.begin(); ilay != layers.end(); ++ilay) {
0208     const MuRingForwardDoubleLayer* layer = (const MuRingForwardDoubleLayer*)(*ilay);
0209
0210     const BoundDisk& disk = layer->specificSurface();
0211
0212     // Generate a random point on the disk
0213     double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi());
0214     double aR = CLHEP::RandFlat::shoot(disk.innerRadius(), disk.outerRadius());
0215     GlobalPoint gp(GlobalPoint::Cylindrical(aR, aPhi, disk.position().z()));
0216
0217     // Momentum: 10 GeV, straight from the origin
0218     GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0219
0220     //FIXME: only negative charge
0221     int charge = -1;
0222
0223     GlobalTrajectoryParameters gtp(gp, gv, charge, field);
0224     TrajectoryStateOnSurface tsos(gtp, disk);
0225     cout << "testCSCLayers: at " << tsos.globalPosition() << " R=" << tsos.globalPosition().perp()
0226          << " phi=" << tsos.globalPosition().phi() << " Z=" << tsos.globalPosition().z()
0227          << " p = " << tsos.globalMomentum() << endl;
0228
0229     SteppingHelixPropagator prop(field, anyDirection);
0230
0231     pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, prop, *theEstimator);
0232     cout << "is compatible: " << comp.first << " at: R=" << comp.second.globalPosition().perp()
0233          << " phi=" << comp.second.globalPosition().phi() << " Z=" << comp.second.globalPosition().z() << endl;
0234
0235     vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, prop, *theEstimator);
0236     if (compDets.size()) {
0237       cout << "compatibleDets: " << compDets.size() << endl
0238
0239            << "  final state pos: " << compDets.front().second.globalPosition() << endl
0240            << "  det         pos: " << compDets.front().first->position()
0241            << " id: " << CSCDetId(compDets.front().first->geographicalId().rawId()) << endl
0242            << "  distance " << (tsos.globalPosition() - compDets.front().first->position()).mag()
0243
0244            << endl
0245            << endl;
0246     } else {
0247       if (layer->isCrack(gp)) {
0248         cout << " CSC crack found ";
0249       } else {
0250         cout << " ERROR : no compatible det found in CSC"
0251              << " at: R=" << gp.perp() << " phi= " << gp.phi().degrees() << " Z= " << gp.z();
0252       }
0253     }
0254   }
0255 }
0256
0257 string MuonRecoGeometryAnalyzer::dumpLayer(const DetLayer* layer) const {
0258   stringstream output;
0259
0260   const BoundSurface* sur = 0;
0261   const BoundCylinder* bc = 0;
0262   const BoundDisk* bd = 0;
0263
0264   sur = &(layer->surface());
0265   if ((bc = dynamic_cast<const BoundCylinder*>(sur))) {
0266     output << "  Cylinder of radius: " << bc->radius() << endl;
0267   } else if ((bd = dynamic_cast<const BoundDisk*>(sur))) {
0268     output << "  Disk at: " << bd->position().z() << endl;
0269   }
0270   return output.str();
0271 }
0272
0273 //define this as a plug-in
0274 #include <FWCore/Framework/interface/MakerMacros.h>
0275 DEFINE_FWK_MODULE(MuonRecoGeometryAnalyzer);