Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:26:50

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