Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:25:46

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