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
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
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
0155 GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0156
0157
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
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
0212 GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0213
0214
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
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
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
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
0300 GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.));
0301
0302
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
0363 #include <FWCore/Framework/interface/MakerMacros.h>
0364 DEFINE_FWK_MODULE(MTDRecoGeometryAnalyzer);