diff --git a/Configuration/Eras/python/Era_Phase2C11_etlV4_cff.py b/Configuration/Eras/python/Era_Phase2C11_etlV4_cff.py new file mode 100644 index 0000000000000..2e7ba054f8db8 --- /dev/null +++ b/Configuration/Eras/python/Era_Phase2C11_etlV4_cff.py @@ -0,0 +1,6 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.Eras.Era_Phase2C11_cff import Phase2C11 +from Configuration.Eras.Modifier_phase2_etlV4_cff import phase2_etlV4 + +Phase2C11_etlV4 = cms.ModifierChain(Phase2C11, phase2_etlV4) diff --git a/Configuration/Eras/python/Modifier_phase2_etlV4_cff.py b/Configuration/Eras/python/Modifier_phase2_etlV4_cff.py new file mode 100644 index 0000000000000..8e9e091355803 --- /dev/null +++ b/Configuration/Eras/python/Modifier_phase2_etlV4_cff.py @@ -0,0 +1,3 @@ +import FWCore.ParameterSet.Config as cms + +phase2_etlV4 = cms.Modifier() diff --git a/Configuration/Geometry/python/dict2026Geometry.py b/Configuration/Geometry/python/dict2026Geometry.py index 6ea0aeebe76d4..fb4b2b9ff9b7a 100644 --- a/Configuration/Geometry/python/dict2026Geometry.py +++ b/Configuration/Geometry/python/dict2026Geometry.py @@ -1186,7 +1186,7 @@ 'from Geometry.MTDGeometryBuilder.idealForDigiMTDGeometry_cff import *', 'mtdGeometry.applyAlignment = cms.bool(False)' ], - "era" : "phase2_timing, phase2_timing_layer", + "era" : "phase2_timing, phase2_timing_layer, phase2_etlV4", }, "I13" : { 1 : [ @@ -1214,7 +1214,7 @@ 'from Geometry.MTDGeometryBuilder.idealForDigiMTDGeometry_cff import *', 'mtdGeometry.applyAlignment = cms.bool(False)' ], - "era" : "phase2_timing, phase2_timing_layer", + "era" : "phase2_timing, phase2_timing_layer, phase2_etlV4", }, } diff --git a/Configuration/PyReleaseValidation/python/relval_2026.py b/Configuration/PyReleaseValidation/python/relval_2026.py index 1b012ca73b0f5..5d18d5f19e907 100644 --- a/Configuration/PyReleaseValidation/python/relval_2026.py +++ b/Configuration/PyReleaseValidation/python/relval_2026.py @@ -29,6 +29,8 @@ numWFIB.extend([31834.0]) #2026D69 numWFIB.extend([32234.0]) #2026D70 numWFIB.extend([32634.0]) #2026D71 +numWFIB.extend([33034.0]) #2026D72 +numWFIB.extend([33434.0]) #2026D73 numWFIB.extend([33834.0]) #2026D74 for numWF in numWFIB: diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index 63c8e53253e7a..3aacb69bb51cf 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -1020,14 +1020,14 @@ def condition(self, fragment, stepList, key, hasHarvest): 'Geom' : 'Extended2026D72', 'HLTmenu': '@fake2', 'GT' : 'auto:phase2_realistic_T21', - 'Era' : 'Phase2C11', + 'Era' : 'Phase2C11_etlV4', 'ScenToRun' : ['GenSimHLBeamSpot','DigiTrigger','RecoGlobal', 'HARVESTGlobal'], }, '2026D73' : { 'Geom' : 'Extended2026D73', 'HLTmenu': '@fake2', 'GT' : 'auto:phase2_realistic_T21', - 'Era' : 'Phase2C11', + 'Era' : 'Phase2C11_etlV4', 'ScenToRun' : ['GenSimHLBeamSpot','DigiTrigger','RecoGlobal', 'HARVESTGlobal'], }, '2026D74' : { diff --git a/Configuration/StandardSequences/python/Eras.py b/Configuration/StandardSequences/python/Eras.py index 5f5524eb6176a..39a50fd5a08a5 100644 --- a/Configuration/StandardSequences/python/Eras.py +++ b/Configuration/StandardSequences/python/Eras.py @@ -42,6 +42,7 @@ def __init__(self): 'Phase2C9_dd4hep', 'Phase2C10_dd4hep', 'Phase2C11_dd4hep', + 'Phase2C11_etlV4', 'Phase2C12_dd4hep', 'Phase2C11M9', ] diff --git a/RecoTracker/TkDetLayers/src/BoundDiskSector.h b/DataFormats/GeometrySurface/interface/BoundDiskSector.h similarity index 79% rename from RecoTracker/TkDetLayers/src/BoundDiskSector.h rename to DataFormats/GeometrySurface/interface/BoundDiskSector.h index 4a5108b4e4765..8f15d188d21d2 100644 --- a/RecoTracker/TkDetLayers/src/BoundDiskSector.h +++ b/DataFormats/GeometrySurface/interface/BoundDiskSector.h @@ -1,10 +1,9 @@ -#ifndef RecoTracker_TkDetLayers_BoundDiskSector_h -#define RecoTracker_TkDetLayers_BoundDiskSector_h +#ifndef DataFormats_GeometrySurface_BoundDiskSector_h +#define DataFormats_GeometrySurface_BoundDiskSector_h #include "DataFormats/GeometrySurface/interface/Plane.h" #include "DiskSectorBounds.h" -#pragma GCC visibility push(hidden) class BoundDiskSector final : public Plane { public: ~BoundDiskSector() override {} @@ -19,5 +18,4 @@ class BoundDiskSector final : public Plane { DiskSectorBounds const& bounds() const { return static_cast(Plane::bounds()); } }; -#pragma GCC visibility pop #endif diff --git a/RecoTracker/TkDetLayers/src/DiskSectorBounds.h b/DataFormats/GeometrySurface/interface/DiskSectorBounds.h similarity index 89% rename from RecoTracker/TkDetLayers/src/DiskSectorBounds.h rename to DataFormats/GeometrySurface/interface/DiskSectorBounds.h index 74ff6f399ca7f..35ce7e6912b4e 100644 --- a/RecoTracker/TkDetLayers/src/DiskSectorBounds.h +++ b/DataFormats/GeometrySurface/interface/DiskSectorBounds.h @@ -1,5 +1,5 @@ -#ifndef RecoTracker_TkDetLayers_DiskSectorBounds_h -#define RecoTracker_TkDetLayers_DiskSectorBounds_h +#ifndef DataFormats_GeometrySurface_DiskSectorBounds_h +#define DataFormats_GeometrySurface_DiskSectorBounds_h #include "DataFormats/GeometryVector/interface/LocalPoint.h" #include "DataFormats/GeometrySurface/interface/LocalError.h" @@ -8,7 +8,6 @@ #include #include -#pragma GCC visibility push(hidden) class DiskSectorBounds final : public Bounds { public: DiskSectorBounds(float rmin, float rmax, float zmin, float zmax, float phiExt) @@ -44,5 +43,4 @@ class DiskSectorBounds final : public Bounds { float theOffset; }; -#pragma GCC visibility pop #endif diff --git a/RecoTracker/TkDetLayers/src/DiskSectorBounds.cc b/DataFormats/GeometrySurface/src/DiskSectorBounds.cc similarity index 95% rename from RecoTracker/TkDetLayers/src/DiskSectorBounds.cc rename to DataFormats/GeometrySurface/src/DiskSectorBounds.cc index 323f4bf1a5243..3ce960d797450 100644 --- a/RecoTracker/TkDetLayers/src/DiskSectorBounds.cc +++ b/DataFormats/GeometrySurface/src/DiskSectorBounds.cc @@ -1,4 +1,4 @@ -#include "DiskSectorBounds.h" +#include "DataFormats/GeometrySurface/interface/DiskSectorBounds.h" using namespace std; diff --git a/Geometry/MTDCommonData/data/mtdParameters/v2/mtdParameters.xml b/Geometry/MTDCommonData/data/mtdParameters/v2/mtdParameters.xml index 1b433a94bc794..0a2b4fc3a560a 100644 --- a/Geometry/MTDCommonData/data/mtdParameters/v2/mtdParameters.xml +++ b/Geometry/MTDCommonData/data/mtdParameters/v2/mtdParameters.xml @@ -9,7 +9,7 @@ 22, 24, 16, 10, 0x1, 0x3, 0x3F, 0x3F, 1, 16, 3, 1 - 22, 24, 16, 7, 0x1, 0x3, 0x3F, 0xFF, 24, 4, 2, 8 + 22, 24, 16, 7, 0x1, 0x3, 0x3F, 0xFF, 16, 16, 1, 2 diff --git a/Geometry/MTDGeometryBuilder/src/MTDTopologyBuilder.cc b/Geometry/MTDGeometryBuilder/src/MTDTopologyBuilder.cc index f85ddd6bc8245..ee0206551dcdb 100644 --- a/Geometry/MTDGeometryBuilder/src/MTDTopologyBuilder.cc +++ b/Geometry/MTDGeometryBuilder/src/MTDTopologyBuilder.cc @@ -1,8 +1,12 @@ +//#define EDM_ML_DEBUG + // Make the change for "big" pixels. 3/06 d.k. #include "Geometry/MTDGeometryBuilder/interface/MTDTopologyBuilder.h" #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h" #include "DataFormats/GeometrySurface/interface/Bounds.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + MTDTopologyBuilder::MTDTopologyBuilder(void) {} PixelTopology* MTDTopologyBuilder::build(const Bounds* bs, @@ -26,6 +30,19 @@ PixelTopology* MTDTopologyBuilder::build(const Bounds* bs, // 2 big pixels per ROC float pitchY = length / (ncols + pixelROCsInY * BIG_PIX_PER_ROC_Y); +#ifdef EDM_ML_DEBUG + edm::LogInfo("MTDTopologyBuilder") << std::fixed << "Building topology for module of width(X) = " << std::setw(10) + << width << " length(Y) = " << std::setw(10) << length + << "\n Rows per ROC = " << std::setw(10) << pixelROCRows + << " Cols per ROC = " << std::setw(10) << pixelROCCols + << "\n ROCs in X = " << std::setw(10) << pixelROCsInX + << " ROCs in Y = " << std::setw(10) << pixelROCsInY + << "\n # pixel rows X = " << std::setw(10) << nrows + << " # pixel cols Y = " << std::setw(10) << ncols + << "\n pitch in X = " << std::setw(10) << pitchX + << " # pitch in Y = " << std::setw(10) << pitchY; +#endif + return (new RectangularMTDTopology(nrows, ncols, pitchX, diff --git a/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc b/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc index 4fcd999f27f89..3c0e8068c9262 100644 --- a/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc +++ b/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc @@ -59,6 +59,8 @@ void MTDDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::Event << "Geometry node for MTDGeom is " << &(*pDD) << "\n" << " # detectors = " << pDD->detUnits().size() << "\n" << " # types = " << pDD->detTypes().size() << "\n" + << " # BTL dets = " << pDD->detsBTL().size() << "\n" + << " # ETL dets = " << pDD->detsETL().size() << "\n" << " # layers " << pDD->geomDetSubDetector(1) << " = " << pDD->numberOfLayers(1) << "\n" << " # layers " << pDD->geomDetSubDetector(2) << " = " << pDD->numberOfLayers(2) << "\n"; sunitt << std::fixed << std::setw(7) << pDD->detUnits().size() << std::setw(7) << pDD->detTypes().size() << "\n"; diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc index e00e2d506e028..f76c4a5510134 100644 --- a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc @@ -211,7 +211,7 @@ void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeomet if (mtdId.mtdSubDetector() == MTDDetId::BTL) { subDet = GeomDetEnumerators::barrel; BTLDetId id = itr->id(); - DetId geoId = id.geographicalId((BTLDetId::CrysLayout)topo->getMTDTopologyMode()); + DetId geoId = id.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode())); const auto& det = geom->idToDet(geoId); const ProxyMTDTopology& topoproxy = static_cast(det->topology()); const RectangularMTDTopology& topol = static_cast(topoproxy.specificTopology()); diff --git a/RecoLocalFastTime/FTLRecProducers/python/mtdRecHits_cfi.py b/RecoLocalFastTime/FTLRecProducers/python/mtdRecHits_cfi.py index d5b098b510ca2..6c86ef07a3515 100644 --- a/RecoLocalFastTime/FTLRecProducers/python/mtdRecHits_cfi.py +++ b/RecoLocalFastTime/FTLRecProducers/python/mtdRecHits_cfi.py @@ -13,6 +13,8 @@ calibrationConstant = cms.double(0.085), # MeV/MIP ) +from Configuration.Eras.Modifier_phase2_etlV4_cff import phase2_etlV4 +phase2_etlV4.toModify(_endcapAlgo, thresholdToKeep = 0.005, calibrationConstant = 0.001 ) mtdRecHits = cms.EDProducer( "MTDRecHitProducer", diff --git a/RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h b/RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h index 3a012d84b8bc4..0fc82b62fae1a 100644 --- a/RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h +++ b/RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h @@ -1,5 +1,5 @@ -#ifndef DetLayers_MTDDetLayerGeometry_h -#define DetLayers_MTDDetLayerGeometry_h +#ifndef RecoMTD_DetLayers_MTDDetLayerGeometry_h +#define RecoMTD_DetLayers_MTDDetLayerGeometry_h /** \class MTDDetLayerGeometry * diff --git a/RecoMTD/DetLayers/interface/MTDDetRing.h b/RecoMTD/DetLayers/interface/MTDDetRing.h index 28d57bab3e6d2..47bac389e52f5 100644 --- a/RecoMTD/DetLayers/interface/MTDDetRing.h +++ b/RecoMTD/DetLayers/interface/MTDDetRing.h @@ -1,5 +1,5 @@ -#ifndef DetLayers_MTDDetRing_H -#define DetLayers_MTDDetRing_H +#ifndef RecoMTD_DetLayers_MTDDetRing_H +#define RecoMTD_DetLayers_MTDDetRing_H /** \class MTDDetRing * A ring of periodic, possibly overlapping vertical detectors. diff --git a/RecoMTD/DetLayers/interface/MTDDetSector.h b/RecoMTD/DetLayers/interface/MTDDetSector.h new file mode 100644 index 0000000000000..7dac4a2e2dc96 --- /dev/null +++ b/RecoMTD/DetLayers/interface/MTDDetSector.h @@ -0,0 +1,77 @@ +#ifndef RecoMTD_DetLayers_MTDDetSector_H +#define RecoMTD_DetLayers_MTDDetSector_H + +#include "TrackingTools/DetLayers/interface/GeometricSearchDet.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" + +#include + +class GeomDet; + +class MTDDetSector : public GeometricSearchDet { +public: + using GeometricSearchDet::GeometricSearchDet; + + /// Construct from iterators on GeomDet* + MTDDetSector(std::vector::const_iterator first, std::vector::const_iterator last); + + /// Construct from a vector of GeomDet* + MTDDetSector(const std::vector& dets); + + ~MTDDetSector() override{}; + + // GeometricSearchDet structure + + const std::vector& basicComponents() const override { return theDets; } + + const BoundSurface& surface() const final { return *theDiskS; } + + const std::vector& components() const override; + + std::pair compatible(const TrajectoryStateOnSurface& ts, + const Propagator& prop, + const MeasurementEstimator& est) const override; + + std::vector compatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const override; + + void compatibleDetsV(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est, + std::vector& result) const override; + + std::vector groupedCompatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const override; + + // GeometricSearchDet extension + + const BoundDiskSector& specificSurface() const { return *theDiskS; } + +protected: + void setDisk(BoundDiskSector* diskS) { theDiskS = diskS; } + + bool add(size_t idet, + std::vector& result, + const TrajectoryStateOnSurface& tsos, + const Propagator& prop, + const MeasurementEstimator& est) const; + +private: + ReferenceCountingPointer theDiskS; + std::vector theDets; + + // Window of detid ordered modules around that closest to the track extrapolation on the sector surface + // needed to limit the size of the vector of distances to sort + // value 50 based on the possible mismatch of module number between adjacent + // modules, due to left-right type imparity + + static constexpr size_t detsRange = 50; + + void init(); +}; + +std::ostream& operator<<(std::ostream&, const MTDDetSector&); + +#endif diff --git a/RecoMTD/DetLayers/interface/MTDDetTray.h b/RecoMTD/DetLayers/interface/MTDDetTray.h index 749fe2d93e80b..f79168322a973 100644 --- a/RecoMTD/DetLayers/interface/MTDDetTray.h +++ b/RecoMTD/DetLayers/interface/MTDDetTray.h @@ -1,5 +1,5 @@ -#ifndef DetLayers_MTDDetTray_H -#define DetLayers_MTDDetTray_H +#ifndef RecoMTD_DetLayers_MTDDetTray_H +#define RecoMTD_DetLayers_MTDDetTray_H /** \class MTDDetTray * A tray of aligned equal-sized non-overlapping detectors. diff --git a/RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h b/RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h index 41ac03461f38e..4ca5dff51a2d8 100644 --- a/RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h +++ b/RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h @@ -1,5 +1,5 @@ -#ifndef DetLayers_MTDRingForwardDoubleLayer_H -#define DetLayers_MTDRingForwardDoubleLayer_H +#ifndef RecoMTD_DetLayers_MTDRingForwardDoubleLayer_H +#define RecoMTD_DetLayers_MTDRingForwardDoubleLayer_H /** \class MTDRingForwardDoubleLayer * A plane composed two layers of disks. The Endcap Timing Layer. diff --git a/RecoMTD/DetLayers/interface/MTDRingForwardLayer.h b/RecoMTD/DetLayers/interface/MTDRingForwardLayer.h index a5a2e4f63f086..bd78ae363e0ec 100644 --- a/RecoMTD/DetLayers/interface/MTDRingForwardLayer.h +++ b/RecoMTD/DetLayers/interface/MTDRingForwardLayer.h @@ -1,5 +1,5 @@ -#ifndef DetLayers_MTDRingForwardLayer_H -#define DetLayers_MTDRingForwardLayer_H +#ifndef RecoMTD_DetLayers_MTDRingForwardLayer_H +#define RecoMTD_DetLayers_MTDRingForwardLayer_H /** \class MTDRingForwardLayer * A plane composed of disks (MTDRingForwardDisk). Represents ETL. diff --git a/RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h b/RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h new file mode 100644 index 0000000000000..e6b1efd26d986 --- /dev/null +++ b/RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h @@ -0,0 +1,65 @@ +#ifndef RecoMTD_DetLayers_MTDSectorForwardDoubleLayer_H +#define RecoMTD_DetLayers_MTDSectorForwardDoubleLayer_H + +#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h" +#include "Utilities/BinningTools/interface/BaseBinFinder.h" +#include "RecoMTD/DetLayers/interface/MTDSectorForwardLayer.h" + +class MTDDetSector; +class GeomDet; + +class MTDSectorForwardDoubleLayer : public ForwardDetLayer { +public: + /// Constructor, takes ownership of pointers + MTDSectorForwardDoubleLayer(const std::vector& frontSectors, + const std::vector& backSectors); + + ~MTDSectorForwardDoubleLayer() override {} + + // GeometricSearchDet interface + + const std::vector& basicComponents() const override { return theBasicComponents; } + + const std::vector& components() const override { return theComponents; } + + bool isInsideOut(const TrajectoryStateOnSurface& tsos) const; + + // tries closest layer first + std::pair compatible(const TrajectoryStateOnSurface&, + const Propagator&, + const MeasurementEstimator&) const override; + + std::vector compatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const override; + + std::vector groupedCompatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const override; + + // DetLayer interface + SubDetector subDetector() const override { return theBackLayer.subDetector(); } + + // Extension of the interface + + /// Return the vector of sectors. + virtual const std::vector& sectors() const { return theSectors; } + + bool isCrack(const GlobalPoint& gp) const; + + const MTDSectorForwardLayer* frontLayer() const { return &theFrontLayer; } + const MTDSectorForwardLayer* backLayer() const { return &theBackLayer; } + + void selfTest() const; + +protected: + BoundDisk* computeSurface() override; + +private: + MTDSectorForwardLayer theFrontLayer; + MTDSectorForwardLayer theBackLayer; + std::vector theSectors; + std::vector theComponents; // duplication of the above + std::vector theBasicComponents; // All chambers +}; +#endif diff --git a/RecoMTD/DetLayers/interface/MTDSectorForwardLayer.h b/RecoMTD/DetLayers/interface/MTDSectorForwardLayer.h new file mode 100644 index 0000000000000..0242b560e2e96 --- /dev/null +++ b/RecoMTD/DetLayers/interface/MTDSectorForwardLayer.h @@ -0,0 +1,43 @@ +#ifndef RecoMTD_DetLayers_MTDSectorForwardLayer_H +#define RecoMTD_DetLayers_MTDSectorForwardLayer_H + +#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h" + +class MTDDetSector; +class GeomDet; + +class MTDSectorForwardLayer : public ForwardDetLayer { +public: + /// Constructor, takes ownership of pointers + MTDSectorForwardLayer(const std::vector& sectors); + + ~MTDSectorForwardLayer() override; + + // GeometricSearchDet interface + + const std::vector& basicComponents() const override { return theBasicComps; } + + const std::vector& components() const override; + + std::vector compatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const override; + + std::vector groupedCompatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const override; + + // DetLayer interface + SubDetector subDetector() const override; + + // Extension of the interface + + /// Return the vector of sectors + virtual const std::vector& sectors() const { return theSectors; } + +private: + std::vector theSectors; + std::vector theComponents; // duplication of the above + std::vector theBasicComps; // All chambers +}; +#endif diff --git a/RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h b/RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h index dbe3e62ff9fea..613bc3fb55bbc 100644 --- a/RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h +++ b/RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h @@ -1,5 +1,5 @@ -#ifndef DetLayers_MTDTrayBarrelLayer_H -#define DetLayers_MTDTrayBarrelLayer_H +#ifndef RecoMTD_DetLayers_MTDTrayBarrelLayer_H +#define RecoMTD_DetLayers_MTDTrayBarrelLayer_H /** \class MTDTrayBarrelLayer * A cylinder composed of half-trays. Represents Barrel Timing Layer. diff --git a/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.cc b/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.cc index 540bd50e2b79d..810adc6fc948f 100644 --- a/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.cc +++ b/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.cc @@ -4,8 +4,11 @@ #include #include +#include +#include #include #include +#include #include #include @@ -15,19 +18,52 @@ using namespace std; -pair, vector > ETLDetLayerGeometryBuilder::buildLayers(const MTDGeometry& geo) { +pair, vector > ETLDetLayerGeometryBuilder::buildLayers(const MTDGeometry& geo, + const int mtdTopologyMode) { vector result[2]; // one for each endcap - for (unsigned endcap = 0; endcap < 2; ++endcap) { - // there is only one layer for ETL right now, maybe more later - for (unsigned layer = 0; layer <= 0; ++layer) { - vector rings; - for (unsigned ring = 1; ring <= 12; ++ring) { - rings.push_back(ring); + if (mtdTopologyMode <= static_cast(MTDTopologyMode::Mode::barphiflat)) { + for (unsigned endcap = 0; endcap < 2; ++endcap) { + // there is only one layer for ETL right now, maybe more later + for (unsigned layer = 0; layer < ETLDetId::kETLv1nDisc; ++layer) { + vector rings; + rings.reserve(ETLDetId::kETLv1maxRing + 1); + for (unsigned ring = 1; ring <= ETLDetId::kETLv1maxRing; ++ring) { + rings.push_back(ring); + } + MTDRingForwardDoubleLayer* thelayer = buildLayer(endcap, layer, rings, geo); + if (thelayer) + result[endcap].push_back(thelayer); + } + } + } else { + // number of layers is identical for post TDR scenarios, pick v4 + // loop on number of sectors per face, two faces per disc (i.e. layer) taken into account in layer building (front/back) + unsigned int nSector(1); + switch (mtdTopologyMode) { + case static_cast(MTDTopologyMode::Mode::btlv1etlv4): + nSector *= ETLDetId::kETLv4maxSector; + break; + case static_cast(MTDTopologyMode::Mode::btlv1etlv5): + nSector *= ETLDetId::kETLv5maxSector; + break; + default: + throw cms::Exception("MTDDetLayers") << "Not implemented scenario " << mtdTopologyMode; + break; + } + + for (unsigned endcap = 0; endcap < 2; ++endcap) { + // number of layers is two, identical for post TDR scenarios, pick v4 + for (unsigned layer = 1; layer <= ETLDetId::kETLv4nDisc; ++layer) { + vector sectors; + sectors.reserve(nSector + 1); + for (unsigned sector = 1; sector <= nSector; ++sector) { + sectors.push_back(sector); + } + MTDSectorForwardDoubleLayer* thelayer = buildLayerNew(endcap, layer, sectors, geo); + if (thelayer) + result[endcap].push_back(thelayer); } - MTDRingForwardDoubleLayer* thelayer = buildLayer(endcap, layer, rings, geo); - if (thelayer) - result[endcap].push_back(thelayer); } } pair, vector > res_pair(result[0], result[1]); @@ -71,7 +107,7 @@ MTDRingForwardDoubleLayer* ETLDetLayerGeometryBuilder::buildLayer(int endcap, assert(!backGeomDets.empty()); float frontz = frontRings[0]->position().z(); float backz = backRings[0]->position().z(); - assert(fabs(frontz) < fabs(backz)); + assert(std::abs(frontz) < std::abs(backz)); } } @@ -96,3 +132,85 @@ MTDDetRing* ETLDetLayerGeometryBuilder::makeDetRing(vector& geom << " R2: " << result->specificSurface().outerRadius(); return result; } + +MTDSectorForwardDoubleLayer* ETLDetLayerGeometryBuilder::buildLayerNew(int endcap, + int layer, + vector& sectors, + const MTDGeometry& geo) { + MTDSectorForwardDoubleLayer* result = nullptr; + + std::vector frontSectors, backSectors; + + LogDebug("MTDDetLayers") << "ETL dets array size = " << geo.detsETL().size(); + + for (unsigned sector : sectors) { + std::vector frontGeomDets, backGeomDets; + LogDebug("MTDDetLayers") << "endcap = " << endcap << " layer = " << layer << " sector = " << sector; +#ifdef EDM_ML_DEBUG + unsigned int nfront(0), nback(0); +#endif + for (auto det : geo.detsETL()) { + ETLDetId theMod(det->geographicalId().rawId()); + if (theMod.mtdSide() == endcap && theMod.nDisc() == layer && theMod.sector() == static_cast(sector)) { + LogTrace("MTDLayerDump") << std::fixed << "ETLDetId " << theMod.rawId() << " side = " << std::setw(4) + << theMod.mtdSide() << " Disc/Side/Sector = " << std::setw(4) << theMod.nDisc() << " " + << std::setw(4) << theMod.discSide() << " " << std::setw(4) << theMod.sector() + << " mod/type = " << std::setw(4) << theMod.module() << " " << std::setw(4) + << theMod.modType() << " pos = " << det->position(); + // front layer face + if (theMod.discSide() == 0) { +#ifdef EDM_ML_DEBUG + nfront++; + LogTrace("MTDDetLayers") << "Front " << theMod.discSide() << " " << nfront; +#endif + frontGeomDets.emplace_back(det); + // back layer face + } else if (theMod.discSide() == 1) { +#ifdef EDM_ML_DEBUG + nback++; + LogTrace("MTDDetLayers") << "Back " << theMod.discSide() << " " << nback; +#endif + backGeomDets.emplace_back(det); + } + } + } + + if (!backGeomDets.empty()) { + std::sort(backGeomDets.begin(), backGeomDets.end(), orderGeomDets); + LogDebug("MTDDetLayers") << "backGeomDets size = " << backGeomDets.size(); + backSectors.emplace_back(makeDetSector(backGeomDets)); + } + + if (!frontGeomDets.empty()) { + std::sort(frontGeomDets.begin(), frontGeomDets.end(), orderGeomDets); + LogDebug("MTDDetLayers") << "frontGeomDets size = " << frontGeomDets.size(); + frontSectors.emplace_back(makeDetSector(frontGeomDets)); + assert(!backGeomDets.empty()); + float frontz = frontSectors.back()->position().z(); + float backz = backSectors.back()->position().z(); + assert(std::abs(frontz) < std::abs(backz)); + } + } + + result = new MTDSectorForwardDoubleLayer(frontSectors, backSectors); + LogTrace("MTDDetLayers") << "New MTDSectorForwardDoubleLayer with " << std::fixed << std::setw(14) + << frontSectors.size() << " and " << std::setw(14) << backSectors.size() << " rings, at Z " + << std::setw(14) << result->specificSurface().position().z() << " R1: " << std::setw(14) + << result->specificSurface().innerRadius() << " R2: " << std::setw(14) + << result->specificSurface().outerRadius(); + + return result; +} + +MTDDetSector* ETLDetLayerGeometryBuilder::makeDetSector(vector& geomDets) { + MTDDetSector* result = new MTDDetSector(geomDets); + LogTrace("MTDDetLayers") << "ETLDetLayerGeometryBuilder::makeDetSector new MTDDetSector with " << std::fixed + << std::setw(14) << geomDets.size() << " modules \n" + << (*result); + + return result; +} + +bool ETLDetLayerGeometryBuilder::orderGeomDets(const GeomDet*& gd1, const GeomDet*& gd2) { + return gd1->geographicalId().rawId() < gd2->geographicalId().rawId(); +} diff --git a/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.h b/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.h index ab760a666205e..d218ded5ed597 100644 --- a/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.h +++ b/RecoMTD/DetLayers/plugins/ETLDetLayerGeometryBuilder.h @@ -14,12 +14,15 @@ class DetLayer; class MTDRingForwardDoubleLayer; class MTDDetRing; +class MTDSectorForwardDoubleLayer; +class MTDDetSector; class ETLDetLayerGeometryBuilder { public: /// return.first=forward (+Z), return.second=backward (-Z) /// both vectors are sorted inside-out - static std::pair, std::vector > buildLayers(const MTDGeometry& geo); + static std::pair, std::vector > buildLayers(const MTDGeometry& geo, + const int mtdTopologyMode); private: // Disable constructor - only static access is allowed. @@ -30,7 +33,14 @@ class ETLDetLayerGeometryBuilder { std::vector& rings, const MTDGeometry& geo); + static MTDSectorForwardDoubleLayer* buildLayerNew(int endcap, + int layer, + std::vector& sectors, + const MTDGeometry& geo); + static MTDDetRing* makeDetRing(std::vector& geomDets); static bool isFront(int layer, int ring, int module); + static MTDDetSector* makeDetSector(std::vector& geomDets); + static bool orderGeomDets(const GeomDet*&, const GeomDet*&); }; #endif diff --git a/RecoMTD/DetLayers/plugins/MTDDetLayerGeometryESProducer.cc b/RecoMTD/DetLayers/plugins/MTDDetLayerGeometryESProducer.cc index d3f2c7a1e5bdc..93a83d0d95bee 100644 --- a/RecoMTD/DetLayers/plugins/MTDDetLayerGeometryESProducer.cc +++ b/RecoMTD/DetLayers/plugins/MTDDetLayerGeometryESProducer.cc @@ -15,6 +15,8 @@ #include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h" #include "Geometry/Records/interface/MTDDigiGeometryRecord.h" #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" +#include "Geometry/Records/interface/MTDTopologyRcd.h" +#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h" #include "ETLDetLayerGeometryBuilder.h" #include "BTLDetLayerGeometryBuilder.h" @@ -36,13 +38,17 @@ class MTDDetLayerGeometryESProducer : public edm::ESProducer { static void fillDescriptions(edm::ConfigurationDescriptions&); private: - const edm::ESGetToken geomToken_; + edm::ESGetToken geomToken_; + edm::ESGetToken mtdtopoToken_; }; using namespace edm; -MTDDetLayerGeometryESProducer::MTDDetLayerGeometryESProducer(const edm::ParameterSet& p) - : geomToken_(setWhatProduced(this).consumes()) {} +MTDDetLayerGeometryESProducer::MTDDetLayerGeometryESProducer(const edm::ParameterSet& p) { + auto cc = setWhatProduced(this); + geomToken_ = cc.consumes(); + mtdtopoToken_ = cc.consumes(); +} std::unique_ptr MTDDetLayerGeometryESProducer::produce(const MTDRecoGeometryRecord& record) { auto mtdDetLayerGeometry = std::make_unique(); @@ -50,8 +56,12 @@ std::unique_ptr MTDDetLayerGeometryESProducer::produce(cons if (auto mtd = record.getHandle(geomToken_)) { // Build BTL layers mtdDetLayerGeometry->addBTLLayers(BTLDetLayerGeometryBuilder::buildLayers(*mtd)); - // Build ETL layers - mtdDetLayerGeometry->addETLLayers(ETLDetLayerGeometryBuilder::buildLayers(*mtd)); + // Build ETL layers, depends on the scenario + if (auto mtdtopo = record.getHandle(mtdtopoToken_)) { + mtdDetLayerGeometry->addETLLayers(ETLDetLayerGeometryBuilder::buildLayers(*mtd, mtdtopo->getMTDTopologyMode())); + } else { + LogWarning("MTDDetLayers") << "No MTD topology is available."; + } } else { LogWarning("MTDDetLayers") << "No MTD geometry is available."; } diff --git a/RecoMTD/DetLayers/src/MTDDetRing.cc b/RecoMTD/DetLayers/src/MTDDetRing.cc index c3496cb3717d7..47599d9b87542 100644 --- a/RecoMTD/DetLayers/src/MTDDetRing.cc +++ b/RecoMTD/DetLayers/src/MTDDetRing.cc @@ -53,10 +53,7 @@ pair MTDDetRing::compatible(const TrajectoryStat } #endif - if (ms.isValid()) - return make_pair(est.estimate(ms, specificSurface()) != 0, ms); - else - return make_pair(false, ms); + return make_pair(ms.isValid() and est.estimate(ms, specificSurface()) != 0, ms); } vector MTDDetRing::compatibleDets(const TrajectoryStateOnSurface& startingState, diff --git a/RecoMTD/DetLayers/src/MTDDetSector.cc b/RecoMTD/DetLayers/src/MTDDetSector.cc new file mode 100644 index 0000000000000..286844aee4ab3 --- /dev/null +++ b/RecoMTD/DetLayers/src/MTDDetSector.cc @@ -0,0 +1,174 @@ +//#define EDM_ML_DEBUG + +#include "RecoMTD/DetLayers/interface/MTDDetSector.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/DetLayers/interface/MeasurementEstimator.h" +#include "MTDDiskSectorBuilderFromDet.h" +#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include +#include +#include + +using namespace std; + +MTDDetSector::MTDDetSector(vector::const_iterator first, vector::const_iterator last) + : GeometricSearchDet(false), theDets(first, last) { + init(); +} + +MTDDetSector::MTDDetSector(const vector& vdets) : GeometricSearchDet(false), theDets(vdets) { init(); } + +void MTDDetSector::init() { + // Add here the sector build based on a collection of GeomDets, mimic what done in ForwardDetRingOneZ + // using the code from tracker BladeShapeBuilderFromDet + // simple initial version, no sorting for the time being + setDisk(MTDDiskSectorBuilderFromDet()(theDets)); +} + +const vector& MTDDetSector::components() const { + // FIXME dummy impl. + edm::LogError("MTDDetLayers") << "temporary dummy implementation of MTDDetSector::components()!!"; + static const vector result; + return result; +} + +pair MTDDetSector::compatible(const TrajectoryStateOnSurface& ts, + const Propagator& prop, + const MeasurementEstimator& est) const { + TrajectoryStateOnSurface ms = prop.propagate(ts, specificSurface()); + +#ifdef EDM_ML_DEBUG + LogTrace("MTDDetLayers") << "MTDDetSector::compatible, sector: \n" + << (*this) << "\n TS at Z,R,phi: " << std::fixed << std::setw(14) << ts.globalPosition().z() + << " , " << std::setw(14) << ts.globalPosition().perp() << " , " << std::setw(14) + << ts.globalPosition().phi(); + if (ms.isValid()) { + LogTrace("MTDDetLayers") << " DEST at Z,R,phi: " << std::fixed << std::setw(14) << ms.globalPosition().z() << " , " + << std::setw(14) << ms.globalPosition().perp() << " , " << std::setw(14) + << ms.globalPosition().phi() << " local Z: " << std::setw(14) << ms.localPosition().z(); + } else { + LogTrace("MTDDetLayers") << " DEST: not valid"; + } +#endif + + return make_pair(ms.isValid() and est.estimate(ms, specificSurface()) != 0, ms); +} + +vector MTDDetSector::compatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const { + LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, sector: \n" + << (*this) << "\n TS at Z,R,phi: " << std::fixed << std::setw(14) + << startingState.globalPosition().z() << " , " << std::setw(14) + << startingState.globalPosition().perp() << " , " << std::setw(14) + << startingState.globalPosition().phi(); + + vector result; + + // Propagate and check that the result is within bounds + pair compat = compatible(startingState, prop, est); + if (!compat.first) { + LogTrace("MTDDetLayers") << " MTDDetSector::compatibleDets: not compatible" + << " (should not have been selected!)"; + return result; + } + + TrajectoryStateOnSurface& tsos = compat.second; + GlobalPoint startPos = tsos.globalPosition(); + + LogTrace("MTDDetLayers") << "Starting position: " << startPos; + + // determine distance of det center from extrapolation on the surface, sort dets accordingly + + size_t idetMin = basicComponents().size(); + double dist2Min = std::numeric_limits::max(); + std::vector > tmpDets; + tmpDets.reserve(basicComponents().size()); + + for (size_t idet = 0; idet < basicComponents().size(); idet++) { + double dist2 = (startPos - theDets[idet]->position()).mag2(); + tmpDets.emplace_back(dist2, idet); + if (dist2 < dist2Min) { + dist2Min = dist2; + idetMin = idet; + } + LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets " << std::fixed << std::setw(8) << idet << " " + << theDets[idet]->geographicalId().rawId() << " dist = " << std::setw(10) + << std::sqrt(dist2) << " Min idet/dist = " << std::setw(8) << idetMin << " " + << std::setw(10) << std::sqrt(dist2Min) << " " << theDets[idet]->position(); + } + + // loop on an interval od ordered detIds around the minimum + // set a range of GeomDets around the minimum compatible with the geometry of ETL + + size_t iniPos(idetMin > detsRange ? idetMin - detsRange : static_cast(0)); + size_t endPos(std::min(idetMin + detsRange, basicComponents().size() - 1)); + tmpDets.erase(tmpDets.begin() + endPos, tmpDets.end()); + tmpDets.erase(tmpDets.begin(), tmpDets.begin() + iniPos); + std::sort(tmpDets.begin(), tmpDets.end()); + + for (const auto& thisDet : tmpDets) { + if (add(thisDet.second, result, tsos, prop, est)) { + LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets found compatible det " << thisDet.second + << " detId = " << theDets[thisDet.second]->geographicalId().rawId() << " at " + << theDets[thisDet.second]->position() << " dist = " << std::sqrt(thisDet.first); + } else { + break; + } + } +#ifdef EDM_ML_DEBUG + if (result.empty()) { + LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, closest not compatible!"; + } else { + LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, found " << result.size() << " compatible dets"; + } +#endif + + return result; +} + +void MTDDetSector::compatibleDetsV(const TrajectoryStateOnSurface&, + const Propagator&, + const MeasurementEstimator&, + std::vector&) const { + edm::LogError("MTDDetLayers") << "At the moment not a real implementation"; +} + +vector MTDDetSector::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const { + // FIXME should be implemented to allow returning overlapping chambers + // as separate groups! + edm::LogInfo("MTDDetLayers") << "dummy implementation of MTDDetSector::groupedCompatibleDets()"; + vector result; + return result; +} + +bool MTDDetSector::add(size_t idet, + vector& result, + const TrajectoryStateOnSurface& tsos, + const Propagator& prop, + const MeasurementEstimator& est) const { + pair compat = theCompatibilityChecker.isCompatible(theDets[idet], tsos, prop, est); + + if (compat.first) { + result.push_back(DetWithState(theDets[idet], compat.second)); + } + + return compat.first; +} + +std::ostream& operator<<(std::ostream& os, const MTDDetSector& id) { + os << " MTDDetSector at " << std::fixed << id.specificSurface().position() << std::endl + << " L/W/T : " << std::setw(14) << id.specificSurface().bounds().length() << " / " << std::setw(14) + << id.specificSurface().bounds().width() << " / " << std::setw(14) << id.specificSurface().bounds().thickness() + << std::endl + << " rmin : " << std::setw(14) << id.specificSurface().innerRadius() << std::endl + << " rmax : " << std::setw(14) << id.specificSurface().outerRadius() << std::endl + << " phi ref : " << std::setw(14) << id.specificSurface().position().phi() << std::endl + << " phi w/2 : " << std::setw(14) << id.specificSurface().phiHalfExtension() << std::endl; + return os; +} diff --git a/RecoMTD/DetLayers/src/MTDDetTray.cc b/RecoMTD/DetLayers/src/MTDDetTray.cc index 306bad38dd087..be983b653a865 100644 --- a/RecoMTD/DetLayers/src/MTDDetTray.cc +++ b/RecoMTD/DetLayers/src/MTDDetTray.cc @@ -37,10 +37,7 @@ pair MTDDetTray::compatible(const TrajectoryStat const Propagator& prop, const MeasurementEstimator& est) const { TrajectoryStateOnSurface ms = prop.propagate(ts, specificSurface()); - if (ms.isValid()) - return make_pair(est.estimate(ms, specificSurface()) != 0, ms); - else - return make_pair(false, ms); + return make_pair(ms.isValid() and est.estimate(ms, specificSurface()) != 0, ms); } vector MTDDetTray::compatibleDets(const TrajectoryStateOnSurface& startingState, diff --git a/RecoMTD/DetLayers/src/MTDDiskSectorBuilderFromDet.cc b/RecoMTD/DetLayers/src/MTDDiskSectorBuilderFromDet.cc new file mode 100644 index 0000000000000..da8bc0caadb5f --- /dev/null +++ b/RecoMTD/DetLayers/src/MTDDiskSectorBuilderFromDet.cc @@ -0,0 +1,99 @@ +//#define EDM_ML_DEBUG + +#include "MTDDiskSectorBuilderFromDet.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DataFormats/GeometryVector/interface/VectorUtil.h" +#include "DataFormats/GeometrySurface/interface/BoundingBox.h" + +#include + +using namespace std; + +namespace { + + pair computeBounds(const vector& dets) { + // go over all corners and compute maximum deviations + float rmin(dets.front()->surface().position().perp()); + float rmax(rmin); + float zmin(dets.front()->surface().position().z()); + float zmax(zmin); + float phimin(dets.front()->surface().position().phi()); + float phimax(phimin); + + for (auto const& idet : dets) { + vector corners = BoundingBox().corners(idet->specificSurface()); + for (auto const& i : corners) { + float r = i.perp(); + float z = i.z(); + float phi = i.phi(); + rmin = min(rmin, r); + rmax = max(rmax, r); + zmin = min(zmin, z); + zmax = max(zmax, z); + if (Geom::phiLess(phi, phimin)) + phimin = phi; + if (Geom::phiLess(phimax, phi)) + phimax = phi; + } + } + + if (!Geom::phiLess(phimin, phimax)) + edm::LogError("MTDDetLayers") << " MTDDiskSectorBuilderFromDet : " + << "Something went wrong with Phi Sorting !"; + float zPos = (zmax + zmin) / 2.; + float phiWin = phimax - phimin; + float phiPos = (phimax + phimin) / 2.; + float rmed = (rmin + rmax) / 2.; + if (phiWin < 0.) { + if ((phimin < Geom::pi() / 2.) || (phimax > -Geom::pi() / 2.)) { + edm::LogError("MTDDetLayers") << " something strange going on, please check " << phimin << " " << phimax << " " + << phiWin; + } + phiWin += 2. * Geom::pi(); + phiPos += Geom::pi(); + } + + GlobalVector pos(rmed * cos(phiPos), rmed * sin(phiPos), zPos); + + LogTrace("MTDDetLayers") << "MTDDiskSectorBuilderFromDet::computeBounds sector at: " << std::fixed << pos << "\n" + << "zmin : " << std::setw(14) << zmin << "\n" + << "zmax : " << std::setw(14) << zmax << "\n" + << "rmin : " << std::setw(14) << rmin << "\n" + << "rmax : " << std::setw(14) << rmax << "\n" + << "phi ref : " << std::setw(14) << phiPos << "\n" + << "phi win : " << std::setw(14) << phiWin; + + return make_pair(new DiskSectorBounds(rmin, rmax, zmin - zPos, zmax - zPos, phiWin), pos); + } + + Surface::RotationType computeRotation(const vector& dets, const Surface::PositionType pos) { + GlobalVector yAxis = (GlobalVector(pos.x(), pos.y(), 0.)).unit(); + + GlobalVector zAxis(0., 0., 1.); + GlobalVector xAxis = yAxis.cross(zAxis); + + return Surface::RotationType(xAxis, yAxis); + } + +} // namespace + +BoundDiskSector* MTDDiskSectorBuilderFromDet::operator()(const vector& dets) const { + // check that the dets are all at about the same z + float zcheck = dets.front()->surface().position().z(); + constexpr double tol(0.5); // minimal safety check on z position of modules within a sector, width ~ 10 mm + for (auto const& idet : dets) { + float zdiff = zcheck - (*idet).surface().position().z(); + if (std::abs(zdiff) > tol) { + edm::LogError("MTDDetLayers") + << " MTDDiskSectorBuilderFromDet: Trying to build sector from Dets at different z positions !! Delta_z = " + << zdiff; + } + } + + auto bo = computeBounds(dets); + + Surface::PositionType pos(bo.second.x(), bo.second.y(), bo.second.z()); + Surface::RotationType rot = computeRotation(dets, pos); + return new BoundDiskSector(pos, rot, bo.first); +} diff --git a/RecoMTD/DetLayers/src/MTDDiskSectorBuilderFromDet.h b/RecoMTD/DetLayers/src/MTDDiskSectorBuilderFromDet.h new file mode 100644 index 0000000000000..504321c927029 --- /dev/null +++ b/RecoMTD/DetLayers/src/MTDDiskSectorBuilderFromDet.h @@ -0,0 +1,20 @@ +#ifndef RecoMTD_DetLayers_MTDDiskSectorBuilderFromDet_H +#define RecoMTD_DetLayers_MTDDiskSectorBuilderFromDet_H + +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/DiskSectorBounds.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include +#include +#include + +/** The trapezoid has the minimal size fully containing all Dets. + */ + +class MTDDiskSectorBuilderFromDet { +public: + BoundDiskSector* operator()(const std::vector& dets) const; +}; + +#endif diff --git a/RecoMTD/DetLayers/src/MTDRingForwardDoubleLayer.cc b/RecoMTD/DetLayers/src/MTDRingForwardDoubleLayer.cc index a7eebd3128f44..f64c225bf25e1 100644 --- a/RecoMTD/DetLayers/src/MTDRingForwardDoubleLayer.cc +++ b/RecoMTD/DetLayers/src/MTDRingForwardDoubleLayer.cc @@ -13,6 +13,8 @@ #include #include +#include + #include #include #include @@ -118,9 +120,9 @@ vector MTDRingForwardDoubleLayer::compatibleDe // This code should be moved in a common place intead of being // copied many times. vector vectorGroups = groupedCompatibleDets(tsos, prop, est); - for (vector::const_iterator itDG = vectorGroups.begin(); itDG != vectorGroups.end(); itDG++) { - for (vector::const_iterator itDGE = itDG->begin(); itDGE != itDG->end(); itDGE++) { - result.push_back(DetWithState(itDGE->det(), itDGE->trajectoryState())); + for (const auto& thisDG : vectorGroups) { + for (const auto& thisDGE : thisDG) { + result.emplace_back(DetWithState(thisDGE.det(), thisDGE.trajectoryState())); } } return result; @@ -170,15 +172,25 @@ void MTDRingForwardDoubleLayer::selfTest() const { const std::vector& frontDets = theFrontLayer.basicComponents(); const std::vector& backDets = theBackLayer.basicComponents(); - std::vector::const_iterator frontItr = frontDets.begin(), lastFront = frontDets.end(), - backItr = backDets.begin(), lastBack = backDets.end(); - - // test that each front z is less than each back z - for (; frontItr != lastFront; ++frontItr) { - float frontz = fabs((**frontItr).surface().position().z()); - for (; backItr != lastBack; ++backItr) { - float backz = fabs((**backItr).surface().position().z()); - assert(frontz < backz); + for (int iring = 0; iring <= ETLDetId::kETLv1maxRing; iring++) { + float frontz(0.); + float backz(1e3f); + for (const auto& thisFront : frontDets) { + if (static_cast(thisFront->geographicalId().rawId()).mtdRR() == iring) { + float tmpz(std::abs(thisFront->surface().position().z())); + if (tmpz > frontz) { + frontz = tmpz; + } + } + } + for (const auto& thisBack : backDets) { + if (static_cast(thisBack->geographicalId().rawId()).mtdRR() == iring) { + float tmpz(std::abs(thisBack->surface().position().z())); + if (tmpz < backz) { + backz = tmpz; + } + } } + assert(frontz < backz); } } diff --git a/RecoMTD/DetLayers/src/MTDRingForwardLayer.cc b/RecoMTD/DetLayers/src/MTDRingForwardLayer.cc index 70ca27484cb55..10420c544f341 100644 --- a/RecoMTD/DetLayers/src/MTDRingForwardLayer.cc +++ b/RecoMTD/DetLayers/src/MTDRingForwardLayer.cc @@ -36,14 +36,14 @@ MTDRingForwardLayer::MTDRingForwardLayer(const vector& ri // Cache chamber pointers (the basic components_) // and find extension in R and Z - for (vector::const_iterator it = rings.begin(); it != rings.end(); it++) { - vector tmp2 = (*it)->basicComponents(); + for (const auto& it : rings) { + vector tmp2 = it->basicComponents(); theBasicComps.insert(theBasicComps.end(), tmp2.begin(), tmp2.end()); - theRmin = min(theRmin, (*it)->specificSurface().innerRadius()); - theRmax = max(theRmax, (*it)->specificSurface().outerRadius()); - float halfThick = (*it)->surface().bounds().thickness() / 2.; - float zCenter = (*it)->surface().position().z(); + theRmin = min(theRmin, it->specificSurface().innerRadius()); + theRmax = max(theRmax, it->specificSurface().outerRadius()); + float halfThick = it->surface().bounds().thickness() / 2.; + float zCenter = it->surface().position().z(); theZmin = min(theZmin, zCenter - halfThick); theZmax = max(theZmax, zCenter + halfThick); } @@ -69,8 +69,8 @@ MTDRingForwardLayer::MTDRingForwardLayer(const vector& ri MTDRingForwardLayer::~MTDRingForwardLayer() { delete theBinFinder; - for (vector::iterator i = theRings.begin(); i < theRings.end(); i++) { - delete *i; + for (auto& i : theRings) { + delete i; } } diff --git a/RecoMTD/DetLayers/src/MTDSectorForwardDoubleLayer.cc b/RecoMTD/DetLayers/src/MTDSectorForwardDoubleLayer.cc new file mode 100644 index 0000000000000..3045885412cc7 --- /dev/null +++ b/RecoMTD/DetLayers/src/MTDSectorForwardDoubleLayer.cc @@ -0,0 +1,174 @@ +//#define EDM_ML_DEBUG + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace std; + +MTDSectorForwardDoubleLayer::MTDSectorForwardDoubleLayer(const vector& frontSectors, + const vector& backSectors) + : ForwardDetLayer(true), + theFrontLayer(frontSectors), + theBackLayer(backSectors), + theSectors(frontSectors), // add back later + theComponents(), + theBasicComponents() { + theSectors.insert(theSectors.end(), backSectors.begin(), backSectors.end()); + theComponents = std::vector(theSectors.begin(), theSectors.end()); + + // Cache chamber pointers (the basic components_) + // and find extension in R and Z + for (const auto& isect : theSectors) { + vector tmp2 = isect->basicComponents(); + theBasicComponents.insert(theBasicComponents.end(), tmp2.begin(), tmp2.end()); + } + + setSurface(computeSurface()); + + LogTrace("MTDDetLayers") << "Constructing MTDSectorForwardDoubleLayer: " << std::fixed << std::setw(14) + << basicComponents().size() << " Dets " << std::setw(14) << theSectors.size() << " Sectors " + << " Z: " << std::setw(14) << specificSurface().position().z() << " R1: " << std::setw(14) + << specificSurface().innerRadius() << " R2: " << std::setw(14) + << specificSurface().outerRadius(); + + selfTest(); +} + +BoundDisk* MTDSectorForwardDoubleLayer::computeSurface() { + const BoundDisk& frontDisk = theFrontLayer.specificSurface(); + const BoundDisk& backDisk = theBackLayer.specificSurface(); + + float rmin = min(frontDisk.innerRadius(), backDisk.innerRadius()); + float rmax = max(frontDisk.outerRadius(), backDisk.outerRadius()); + float zmin = frontDisk.position().z(); + float halfThickness = frontDisk.bounds().thickness() / 2.; + zmin = (zmin > 0) ? zmin - halfThickness : zmin + halfThickness; + float zmax = backDisk.position().z(); + halfThickness = backDisk.bounds().thickness() / 2.; + zmax = (zmax > 0) ? zmax + halfThickness : zmax - halfThickness; + float zPos = (zmax + zmin) / 2.; + PositionType pos(0., 0., zPos); + RotationType rot; + + return new BoundDisk(pos, rot, new SimpleDiskBounds(rmin, rmax, zmin - zPos, zmax - zPos)); +} + +bool MTDSectorForwardDoubleLayer::isInsideOut(const TrajectoryStateOnSurface& tsos) const { + return tsos.globalPosition().basicVector().dot(tsos.globalMomentum().basicVector()) > 0; +} + +std::pair MTDSectorForwardDoubleLayer::compatible( + const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const { + // mostly copied from ForwardDetLayer, except propagates to closest surface, + // not to center + + bool insideOut = isInsideOut(startingState); + const MTDSectorForwardLayer& closerLayer = (insideOut) ? theFrontLayer : theBackLayer; + LogTrace("MTDDetLayers") << "MTDSectorForwardDoubleLayer::compatible is assuming inside-out direction: " << insideOut; + + TrajectoryStateOnSurface myState = prop.propagate(startingState, closerLayer.specificSurface()); + if (!myState.isValid()) + return make_pair(false, myState); + + // take into account the thickness of the layer + float deltaR = surface().bounds().thickness() / 2. * std::abs(tan(myState.localDirection().theta())); + + // take into account the error on the predicted state + constexpr float nSigma = 3.; + if (myState.hasError()) { + LocalError err = myState.localError().positionError(); + // ignore correlation for the moment... + deltaR += nSigma * sqrt(err.xx() + err.yy()); + } + + float zPos = (zmax() + zmin()) / 2.; + SimpleDiskBounds tmp(rmin() - deltaR, rmax() + deltaR, zmin() - zPos, zmax() - zPos); + + return make_pair(tmp.inside(myState.localPosition()), myState); +} + +vector MTDSectorForwardDoubleLayer::compatibleDets( + const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const { + vector result; + pair compat = compatible(startingState, prop, est); + + if (!compat.first) { + LogTrace("MTDDetLayers") << " MTDSectorForwardDoubleLayer::compatibleDets: not compatible" + << " (should not have been selected!)"; + return result; + } + + TrajectoryStateOnSurface& tsos = compat.second; + + // standard implementation of compatibleDets() for class which have + // groupedCompatibleDets implemented. + // This code should be moved in a common place intead of being + // copied many times. + vector vectorGroups(groupedCompatibleDets(tsos, prop, est)); + for (const auto& thisDG : vectorGroups) { + for (const auto& thisDGE : thisDG) { + result.emplace_back(DetWithState(thisDGE.det(), thisDGE.trajectoryState())); + } + } + + return result; +} + +vector MTDSectorForwardDoubleLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const { + vector detWithStates1, detWithStates2; + + LogTrace("MTDDetLayers") << "groupedCompatibleDets are currently given always in inside-out order"; + // this should be fixed either in RecoMTD/MeasurementDet/MTDDetLayerMeasurements or + // RecoMTD/DetLayers/MTDSectorForwardDoubleLayer + + detWithStates1 = theFrontLayer.compatibleDets(startingState, prop, est); + detWithStates2 = theBackLayer.compatibleDets(startingState, prop, est); + + vector result; + if (!detWithStates1.empty()) + result.push_back(DetGroup(detWithStates1)); + if (!detWithStates2.empty()) + result.push_back(DetGroup(detWithStates2)); + LogTrace("MTDDetLayers") << "DoubleLayer Compatible dets: " << result.size(); + return result; +} + +bool MTDSectorForwardDoubleLayer::isCrack(const GlobalPoint& gp) const { + bool result = false; + LogTrace("MTDDetLayers") + << "MTDSectorForwardDoubleLayer::isCrack kept only for backward compatibility, no real implementation"; + return result; +} + +void MTDSectorForwardDoubleLayer::selfTest() const { + const std::vector& frontDets = theFrontLayer.basicComponents(); + const std::vector& backDets = theBackLayer.basicComponents(); + + // test that each front z is less than each back z + float frontz(0.); + float backz(1e3f); + for (const auto& thisFront : frontDets) { + float tmpz(std::abs(thisFront->surface().position().z())); + if (tmpz > frontz) { + frontz = tmpz; + } + } + for (const auto& thisBack : backDets) { + float tmpz(std::abs(thisBack->surface().position().z())); + if (tmpz < backz) { + backz = tmpz; + } + } + assert(frontz < backz); +} diff --git a/RecoMTD/DetLayers/src/MTDSectorForwardLayer.cc b/RecoMTD/DetLayers/src/MTDSectorForwardLayer.cc new file mode 100644 index 0000000000000..f2ba187d7a434 --- /dev/null +++ b/RecoMTD/DetLayers/src/MTDSectorForwardLayer.cc @@ -0,0 +1,140 @@ +//#define EDM_ML_DEBUG + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +using namespace std; + +MTDSectorForwardLayer::MTDSectorForwardLayer(const vector& sectors) + : ForwardDetLayer(false), theSectors(sectors), theComponents(theSectors.begin(), theSectors.end()) { + // Initial values for R, Z and Phi bounds + float theRmin = sectors.front()->basicComponents().front()->position().perp(); + float theRmax = theRmin; + float theZmin = sectors.front()->position().z(); + float theZmax = theZmin; + + // Cache chamber pointers (the basic components_) + // and find extension in R and Z + for (const auto& isect : sectors) { + vector tmp2 = isect->basicComponents(); + theBasicComps.insert(theBasicComps.end(), tmp2.begin(), tmp2.end()); + + theRmin = min(theRmin, isect->specificSurface().innerRadius()); + theRmax = max(theRmax, isect->specificSurface().outerRadius()); + float halfThick = isect->surface().bounds().thickness() / 2.; + float zCenter = isect->surface().position().z(); + theZmin = min(theZmin, zCenter - halfThick); + theZmax = max(theZmax, zCenter + halfThick); + } + + // Build surface + + float zPos = (theZmax + theZmin) / 2.; + PositionType pos(0., 0., zPos); + RotationType rot; + + setSurface(new BoundDisk(pos, rot, new SimpleDiskBounds(theRmin, theRmax, theZmin - zPos, theZmax - zPos))); + + LogTrace("MTDDetLayers") << "Constructing MTDSectorForwardLayer: " << std::fixed << std::setw(14) + << basicComponents().size() << " Dets, " << std::setw(14) << theSectors.size() + << " Sectors, " + << " Z: " << std::setw(14) << specificSurface().position().z() << " R1: " << std::setw(14) + << specificSurface().innerRadius() << " R2: " << std::setw(14) + << specificSurface().outerRadius(); +} + +MTDSectorForwardLayer::~MTDSectorForwardLayer() { + for (auto& i : theSectors) { + delete i; + } +} + +vector MTDSectorForwardLayer::compatibleDets( + const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const { + vector result; + + LogTrace("MTDDetLayers") << "MTDSectorForwardLayer::compatibleDets," + << " R1 " << std::fixed << std::setw(14) << specificSurface().innerRadius() + << " R2: " << std::setw(14) << specificSurface().outerRadius() + << " FTS at R: " << std::setw(14) << startingState.globalPosition().perp(); + + pair compat = compatible(startingState, prop, est); + + if (!compat.first) { + LogTrace("MTDDetLayers") << " MTDSectorForwardLayer::compatibleDets: not compatible" + << " (should not have been selected!)"; + return result; + } + + TrajectoryStateOnSurface& tsos = compat.second; + + // as there are either two or four sectors only, avoid complex logic and just loop on all of them + + // Use state on layer surface. Note that local coordinates and errors + // are the same on the layer and on all sectors surfaces, since + // all BoundDisks are centered in 0,0 and have the same rotation. + // CAVEAT: if the sectors are not at the same Z, the local position and error + // will be "Z-projected" to the sectors. This is a fairly good approximation. + // However in this case additional propagation will be done when calling + // compatibleDets. + GlobalPoint startPos = tsos.globalPosition(); + + for (unsigned int isect = 0; isect < theSectors.size(); isect++) { + LocalPoint nextPos(theSectors[isect]->specificSurface().toLocal(startPos)); + LogDebug("MTDDetLayers") << "Global point = " << std::fixed << startPos << " local point = " << nextPos + << " global sector ref pos = " << theSectors[isect]->specificSurface().position(); + bool inside = false; + if (tsos.hasError()) { + inside = theSectors[isect]->specificSurface().bounds().inside(nextPos, tsos.localError().positionError(), 1.); + } else { + inside = theSectors[isect]->specificSurface().bounds().inside(nextPos); + } + if (inside) { +#ifdef EDM_ML_DEBUG + LogTrace("MTDDetLayers") << " MTDSectorForwardLayer::fastCompatibleDets:NextSector " << std::fixed + << std::setw(14) << isect << "\n" + << (*theSectors[isect]) << "\n FTS at Z,R,phi: " << std::setw(14) + << tsos.globalPosition().z() << " , " << std::setw(14) << tsos.globalPosition().perp() + << "," << std::setw(14) << tsos.globalPosition().phi(); + if (tsos.hasError()) { + LogTrace("MTDDetLayers") << " sR: " << sqrt(tsos.localError().positionError().yy()) + << " sX: " << sqrt(tsos.localError().positionError().xx()); + } +#endif + vector nextRodDets(theSectors[isect]->compatibleDets(tsos, prop, est)); + if (!nextRodDets.empty()) { + result.insert(result.end(), nextRodDets.begin(), nextRodDets.end()); + } else { + break; + } + } + } + + LogTrace("MTDDetLayers") << " MTDSectorForwardLayer::fastCompatibleDets: found: " << result.size(); + + return result; +} + +vector MTDSectorForwardLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) const { + // FIXME should return only 1 group + edm::LogInfo("MTDDetLayers") << "dummy implementation of MTDSectorForwardLayer::groupedCompatibleDets()"; + return vector(); +} + +GeomDetEnumerators::SubDetector MTDSectorForwardLayer::subDetector() const { + return theBasicComps.front()->subDetector(); +} + +const vector& MTDSectorForwardLayer::components() const { return theComponents; } diff --git a/RecoMTD/DetLayers/src/MTDTrayBarrelLayer.cc b/RecoMTD/DetLayers/src/MTDTrayBarrelLayer.cc index 43bea1d27b1ab..a845805e77f2e 100644 --- a/RecoMTD/DetLayers/src/MTDTrayBarrelLayer.cc +++ b/RecoMTD/DetLayers/src/MTDTrayBarrelLayer.cc @@ -33,8 +33,8 @@ MTDTrayBarrelLayer::MTDTrayBarrelLayer(vector& rods) std::copy(theRods.begin(), theRods.end(), back_inserter(theComponents)); // Cache chamber pointers (the basic components_) - for (vector::const_iterator it = rods.begin(); it != rods.end(); it++) { - vector tmp2 = (*it)->basicComponents(); + for (const auto& it : rods) { + vector tmp2 = it->basicComponents(); theBasicComps.insert(theBasicComps.end(), tmp2.begin(), tmp2.end()); } @@ -59,8 +59,8 @@ MTDTrayBarrelLayer::MTDTrayBarrelLayer(vector& rods) MTDTrayBarrelLayer::~MTDTrayBarrelLayer() { delete theBinFinder; - for (vector::iterator i = theRods.begin(); i < theRods.end(); i++) { - delete *i; + for (auto& i : theRods) { + delete i; } } diff --git a/RecoMTD/DetLayers/test/BuildFile.xml b/RecoMTD/DetLayers/test/BuildFile.xml index 480d6ff8e48f7..2d8892ffe2fd8 100644 --- a/RecoMTD/DetLayers/test/BuildFile.xml +++ b/RecoMTD/DetLayers/test/BuildFile.xml @@ -8,4 +8,6 @@ + + diff --git a/RecoMTD/DetLayers/test/MTDRecoGeometryAnalyzer.cc b/RecoMTD/DetLayers/test/MTDRecoGeometryAnalyzer.cc index 4ff45152ed8e2..f7f50a5454b41 100644 --- a/RecoMTD/DetLayers/test/MTDRecoGeometryAnalyzer.cc +++ b/RecoMTD/DetLayers/test/MTDRecoGeometryAnalyzer.cc @@ -2,7 +2,7 @@ #include "FWCore/Framework/interface/EDAnalyzer.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/EventSetup.h" -#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/ESTransientHandle.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -12,6 +12,10 @@ #include "MagneticField/Engine/interface/MagneticField.h" #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "Geometry/Records/interface/MTDTopologyRcd.h" +#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h" +#include "Geometry/MTDCommonData/interface/MTDTopologyMode.h" + #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h" #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" @@ -19,6 +23,8 @@ #include "RecoMTD/DetLayers/interface/MTDDetTray.h" #include "RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h" #include "RecoMTD/DetLayers/interface/MTDDetRing.h" +#include "RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h" +#include "RecoMTD/DetLayers/interface/MTDDetSector.h" #include @@ -40,51 +46,68 @@ class MTDRecoGeometryAnalyzer : public EDAnalyzer { void testBTLLayers(const MTDDetLayerGeometry*, const MagneticField* field); void testETLLayers(const MTDDetLayerGeometry*, const MagneticField* field); + void testETLLayersNew(const MTDDetLayerGeometry*, const MagneticField* field); string dumpLayer(const DetLayer* layer) const; private: MeasurementEstimator* theEstimator; + + const edm::ESInputTag tag_; + edm::ESGetToken geomToken_; + edm::ESGetToken mtdtopoToken_; + edm::ESGetToken magfieldToken_; }; -MTDRecoGeometryAnalyzer::MTDRecoGeometryAnalyzer(const ParameterSet& iConfig) { +MTDRecoGeometryAnalyzer::MTDRecoGeometryAnalyzer(const ParameterSet& iConfig) : tag_(edm::ESInputTag{"", ""}) { + geomToken_ = esConsumes(tag_); + mtdtopoToken_ = esConsumes(tag_); + magfieldToken_ = esConsumes(tag_); + float theMaxChi2 = 25.; float theNSigma = 3.; theEstimator = new Chi2MeasurementEstimator(theMaxChi2, theNSigma); } void MTDRecoGeometryAnalyzer::analyze(const Event& ev, const EventSetup& es) { - ESHandle geo; - es.get().get(geo); + auto geo = es.getTransientHandle(geomToken_); + auto mtdtopo = es.getTransientHandle(mtdtopoToken_); + auto magfield = es.getTransientHandle(magfieldToken_); - ESHandle magfield; - es.get().get(magfield); // Some printouts - LogVerbatim("MTDLayerDump") << "\n*** allBTLLayers(): " << geo->allBTLLayers().size(); + LogVerbatim("MTDLayerDump") << "\n*** allBTLLayers(): " << std::fixed << std::setw(14) << geo->allBTLLayers().size(); for (auto dl = geo->allBTLLayers().begin(); dl != geo->allBTLLayers().end(); ++dl) { - LogVerbatim("MTDLayerDump") << " " << (int)(dl - geo->allBTLLayers().begin()) << " " << dumpLayer(*dl); + LogVerbatim("MTDLayerDump") << " " << static_cast(dl - geo->allBTLLayers().begin()) << " " << dumpLayer(*dl); } - LogVerbatim("MTDLayerDump") << "\n*** allETLLayers(): " << geo->allETLLayers().size(); + LogVerbatim("MTDLayerDump") << "\n*** allETLLayers(): " << std::fixed << std::setw(14) << geo->allETLLayers().size(); for (auto dl = geo->allETLLayers().begin(); dl != geo->allETLLayers().end(); ++dl) { - LogVerbatim("MTDLayerDump") << " " << (int)(dl - geo->allETLLayers().begin()) << " " << dumpLayer(*dl); + LogVerbatim("MTDLayerDump") << " " << static_cast(dl - geo->allETLLayers().begin()) << " " << dumpLayer(*dl); } - LogVerbatim("MTDLayerDump") << "\n*** allLayers(): " << geo->allLayers().size(); + LogVerbatim("MTDLayerDump") << "\n*** allLayers(): " << std::fixed << std::setw(14) << geo->allLayers().size(); for (auto dl = geo->allLayers().begin(); dl != geo->allLayers().end(); ++dl) { - LogVerbatim("MTDLayerDump") << " " << (int)(dl - geo->allLayers().begin()) << " " << dumpLayer(*dl); + LogVerbatim("MTDLayerDump") << " " << static_cast(dl - geo->allLayers().begin()) << " " << dumpLayer(*dl); } testBTLLayers(geo.product(), magfield.product()); - testETLLayers(geo.product(), magfield.product()); + if (mtdtopo->getMTDTopologyMode() <= static_cast(MTDTopologyMode::Mode::barphiflat)) { + testETLLayers(geo.product(), magfield.product()); + } else { + testETLLayersNew(geo.product(), magfield.product()); + } } void MTDRecoGeometryAnalyzer::testBTLLayers(const MTDDetLayerGeometry* geo, const MagneticField* field) { const vector& layers = geo->allBTLLayers(); - for (auto ilay = layers.begin(); ilay != layers.end(); ++ilay) { - const MTDTrayBarrelLayer* layer = (const MTDTrayBarrelLayer*)(*ilay); + for (const auto& ilay : layers) { + const MTDTrayBarrelLayer* layer = static_cast(ilay); + + LogVerbatim("MTDLayerDump") << std::fixed << "\nBTL layer " << std::setw(4) << layer->subDetector() + << " rods = " << std::setw(14) << layer->rods().size() << " dets = " << std::setw(14) + << layer->basicComponents().size(); const BoundCylinder& cyl = layer->specificSurface(); @@ -103,25 +126,29 @@ void MTDRecoGeometryAnalyzer::testBTLLayers(const MTDDetLayerGeometry* geo, cons GlobalTrajectoryParameters gtp(gp, gv, charge, field); TrajectoryStateOnSurface tsos(gtp, cyl); - LogVerbatim("MTDLayerDump") << "\ntestBTLLayers: at " << tsos.globalPosition() - << " R=" << tsos.globalPosition().perp() << " phi=" << tsos.globalPosition().phi() - << " Z=" << tsos.globalPosition().z() << " p = " << tsos.globalMomentum(); + LogVerbatim("MTDLayerDump") << "\ntestBTLLayers: at " << std::fixed << std::setw(14) << tsos.globalPosition() + << " R=" << std::setw(14) << tsos.globalPosition().perp() << " phi=" << std::setw(14) + << tsos.globalPosition().phi() << " Z=" << std::setw(14) << tsos.globalPosition().z() + << " p = " << std::setw(14) << tsos.globalMomentum(); SteppingHelixPropagator prop(field, anyDirection); pair comp = layer->compatible(tsos, prop, *theEstimator); - LogVerbatim("MTDLayerDump") << "is compatible: " << comp.first << " at: R=" << comp.second.globalPosition().perp() - << " phi=" << comp.second.globalPosition().phi() - << " Z=" << comp.second.globalPosition().z(); + LogVerbatim("MTDLayerDump") << "is compatible: " << comp.first << " at: R=" << std::setw(14) + << comp.second.globalPosition().perp() << " phi=" << std::setw(14) + << comp.second.globalPosition().phi() << " Z=" << std::setw(14) + << comp.second.globalPosition().z(); vector compDets = layer->compatibleDets(tsos, prop, *theEstimator); if (compDets.size()) { - LogVerbatim("MTDLayerDump") << "compatibleDets: " << compDets.size() << "\n" - << " final state pos: " << compDets.front().second.globalPosition() << "\n" - << " det pos: " << compDets.front().first->position() << " id: " << std::hex + LogVerbatim("MTDLayerDump") << "compatibleDets: " << std::setw(14) << compDets.size() << "\n" + << " final state pos: " << std::setw(14) << compDets.front().second.globalPosition() + << "\n" + << " det pos: " << std::setw(14) << compDets.front().first->position() + << " id: " << std::hex << BTLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec << "\n" - << " distance " + << " distance " << std::setw(14) << (tsos.globalPosition() - compDets.front().first->position()).mag(); } else { LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found"; @@ -132,8 +159,14 @@ void MTDRecoGeometryAnalyzer::testBTLLayers(const MTDDetLayerGeometry* geo, cons void MTDRecoGeometryAnalyzer::testETLLayers(const MTDDetLayerGeometry* geo, const MagneticField* field) { const vector& layers = geo->allETLLayers(); - for (auto ilay = layers.begin(); ilay != layers.end(); ++ilay) { - const MTDRingForwardDoubleLayer* layer = (const MTDRingForwardDoubleLayer*)(*ilay); + for (const auto& ilay : layers) { + const MTDRingForwardDoubleLayer* layer = static_cast(ilay); + + LogVerbatim("MTDLayerDump") << std::fixed << "\nETL layer " << std::setw(4) << layer->subDetector() + << " rings = " << std::setw(14) << layer->rings().size() << " dets = " << std::setw(14) + << layer->basicComponents().size() << " front dets = " << std::setw(14) + << layer->frontLayer()->basicComponents().size() << " back dets = " << std::setw(14) + << layer->backLayer()->basicComponents().size(); const BoundDisk& disk = layer->specificSurface(); @@ -150,32 +183,121 @@ void MTDRecoGeometryAnalyzer::testETLLayers(const MTDDetLayerGeometry* geo, cons GlobalTrajectoryParameters gtp(gp, gv, charge, field); TrajectoryStateOnSurface tsos(gtp, disk); - LogVerbatim("MTDLayerDump") << "\ntestETLLayers: at " << tsos.globalPosition() - << " R=" << tsos.globalPosition().perp() << " phi=" << tsos.globalPosition().phi() - << " Z=" << tsos.globalPosition().z() << " p = " << tsos.globalMomentum(); + LogVerbatim("MTDLayerDump") << "\ntestETLLayers: at " << std::setw(14) << tsos.globalPosition() + << " R=" << std::setw(14) << tsos.globalPosition().perp() << " phi=" << std::setw(14) + << tsos.globalPosition().phi() << " Z=" << std::setw(14) << tsos.globalPosition().z() + << " p = " << std::setw(14) << tsos.globalMomentum(); SteppingHelixPropagator prop(field, anyDirection); pair comp = layer->compatible(tsos, prop, *theEstimator); - LogVerbatim("MTDLayerDump") << "is compatible: " << comp.first << " at: R=" << comp.second.globalPosition().perp() - << " phi=" << comp.second.globalPosition().phi() - << " Z=" << comp.second.globalPosition().z(); + LogVerbatim("MTDLayerDump") << "is compatible: " << comp.first << " at: R=" << std::setw(14) + << comp.second.globalPosition().perp() << " phi=" << std::setw(14) + << comp.second.globalPosition().phi() << " Z=" << std::setw(14) + << comp.second.globalPosition().z(); vector compDets = layer->compatibleDets(tsos, prop, *theEstimator); if (compDets.size()) { - LogVerbatim("MTDLayerDump") << "compatibleDets: " << compDets.size() << "\n" - << " final state pos: " << compDets.front().second.globalPosition() << "\n" - << " det pos: " << compDets.front().first->position() << " id: " << std::hex + LogVerbatim("MTDLayerDump") << "compatibleDets: " << std::setw(14) << compDets.size() << "\n" + << " final state pos: " << std::setw(14) << compDets.front().second.globalPosition() + << "\n" + << " det pos: " << std::setw(14) << compDets.front().first->position() + << " id: " << std::hex << ETLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec << "\n" - << " distance " + << " distance " << std::setw(14) << (tsos.globalPosition() - compDets.front().first->position()).mag(); } else { if (layer->isCrack(gp)) { LogVerbatim("MTDLayerDump") << " MTD crack found "; } else { LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found in MTD" - << " at: R=" << gp.perp() << " phi= " << gp.phi().degrees() << " Z= " << gp.z(); + << " at: R=" << std::setw(14) << gp.perp() << " phi= " << std::setw(14) + << gp.phi().degrees() << " Z= " << std::setw(14) << gp.z(); + } + } + } +} + +void MTDRecoGeometryAnalyzer::testETLLayersNew(const MTDDetLayerGeometry* geo, const MagneticField* field) { + const vector& layers = geo->allETLLayers(); + + // dump of ETL layers structure + + for (const auto& ilay : layers) { + const MTDSectorForwardDoubleLayer* layer = static_cast(ilay); + + LogVerbatim("MTDLayerDump") << std::fixed << "\nETL layer " << std::setw(4) << layer->subDetector() + << " at z = " << std::setw(14) << layer->surface().position().z() + << " sectors = " << std::setw(14) << layer->sectors().size() + << " dets = " << std::setw(14) << layer->basicComponents().size() + << " front dets = " << std::setw(14) << layer->frontLayer()->basicComponents().size() + << " back dets = " << std::setw(14) << layer->backLayer()->basicComponents().size(); + + unsigned int isectInd(0); + for (const auto& isector : layer->sectors()) { + isectInd++; + LogVerbatim("MTDLayerDump") << std::fixed << "\nSector " << std::setw(4) << isectInd << "\n" << (*isector); + for (const auto& imod : isector->basicComponents()) { + ETLDetId modId(imod->geographicalId().rawId()); + LogVerbatim("MTDLayerDump") << std::fixed << "ETLDetId " << modId.rawId() << " side = " << std::setw(4) + << modId.mtdSide() << " Disc/Side/Sector = " << std::setw(4) << modId.nDisc() << " " + << std::setw(4) << modId.discSide() << " " << std::setw(4) << modId.sector() + << " mod/type = " << std::setw(4) << modId.module() << " " << std::setw(4) + << modId.modType() << " pos = " << imod->position(); + } + } + } + + // test propagation through layers + + for (const auto& ilay : layers) { + const MTDSectorForwardDoubleLayer* layer = static_cast(ilay); + + const BoundDisk& disk = layer->specificSurface(); + + // Generate a random point on the disk + double aPhi = CLHEP::RandFlat::shoot(-Geom::pi(), Geom::pi()); + double aR = CLHEP::RandFlat::shoot(disk.innerRadius(), disk.outerRadius()); + GlobalPoint gp(GlobalPoint::Cylindrical(aR, aPhi, disk.position().z())); + + // Momentum: 10 GeV, straight from the origin + GlobalVector gv(GlobalVector::Spherical(gp.theta(), aPhi, 10.)); + + //FIXME: only negative charge + int charge = -1; + + GlobalTrajectoryParameters gtp(gp, gv, charge, field); + TrajectoryStateOnSurface tsos(gtp, disk); + LogVerbatim("MTDLayerDump") << "\ntestETLLayers: at " << std::fixed << tsos.globalPosition() + << " R=" << std::setw(14) << tsos.globalPosition().perp() << " phi=" << std::setw(14) + << tsos.globalPosition().phi() << " Z=" << std::setw(14) << tsos.globalPosition().z() + << " p = " << tsos.globalMomentum(); + + SteppingHelixPropagator prop(field, anyDirection); + + pair comp = layer->compatible(tsos, prop, *theEstimator); + LogVerbatim("MTDLayerDump") << std::fixed << "is compatible: " << comp.first << " at: R=" << std::setw(14) + << comp.second.globalPosition().perp() << " phi=" << std::setw(14) + << comp.second.globalPosition().phi() << " Z=" << std::setw(14) + << comp.second.globalPosition().z(); + + vector compDets = layer->compatibleDets(tsos, prop, *theEstimator); + if (compDets.size()) { + LogVerbatim("MTDLayerDump") + << std::fixed << "compatibleDets: " << std::setw(14) << compDets.size() << "\n" + << " final state pos: " << compDets.front().second.globalPosition() << "\n" + << " det pos: " << compDets.front().first->position() << " id: " << std::hex + << ETLDetId(compDets.front().first->geographicalId().rawId()).rawId() << std::dec << "\n" + << " distance " << std::setw(14) + << (compDets.front().second.globalPosition() - compDets.front().first->position()).mag(); + } else { + if (layer->isCrack(gp)) { + LogVerbatim("MTDLayerDump") << " MTD crack found "; + } else { + LogVerbatim("MTDLayerDump") << " ERROR : no compatible det found in MTD" + << " at: R=" << std::setw(14) << gp.perp() << " phi= " << std::setw(14) + << gp.phi().degrees() << " Z= " << std::setw(14) << gp.z(); } } } @@ -190,9 +312,11 @@ string MTDRecoGeometryAnalyzer::dumpLayer(const DetLayer* layer) const { sur = &(layer->surface()); if ((bc = dynamic_cast(sur))) { - output << " Cylinder of radius: " << bc->radius(); + output << std::fixed << " subdet = " << std::setw(4) << layer->subDetector() << " Barrel = " << layer->isBarrel() + << " Forward = " << layer->isForward() << " Cylinder of radius: " << std::setw(14) << bc->radius(); } else if ((bd = dynamic_cast(sur))) { - output << " Disk at: " << bd->position().z(); + output << std::fixed << " subdet = " << std::setw(4) << layer->subDetector() << " Barrel = " << layer->isBarrel() + << " Forward = " << layer->isForward() << " Disk at: " << std::setw(14) << bd->position().z(); } return output.str(); } diff --git a/RecoMTD/Records/interface/MTDRecoGeometryRecord.h b/RecoMTD/Records/interface/MTDRecoGeometryRecord.h index ba76abff88e18..b2fb975eaade6 100644 --- a/RecoMTD/Records/interface/MTDRecoGeometryRecord.h +++ b/RecoMTD/Records/interface/MTDRecoGeometryRecord.h @@ -11,11 +11,13 @@ #include "FWCore/Framework/interface/EventSetupRecordImplementation.h" #include "FWCore/Framework/interface/DependentRecordImplementation.h" #include "Geometry/Records/interface/MTDDigiGeometryRecord.h" +#include "Geometry/Records/interface/MTDTopologyRcd.h" #include "FWCore/Utilities/interface/mplVector.h" class MTDRecoGeometryRecord : public edm::eventsetup::DependentRecordImplementation > {}; + edm::mpl::Vector > { +}; #endif diff --git a/RecoTracker/TkDetLayers/src/BladeShapeBuilderFromDet.h b/RecoTracker/TkDetLayers/src/BladeShapeBuilderFromDet.h index 49689238188c3..107b4bedeb290 100644 --- a/RecoTracker/TkDetLayers/src/BladeShapeBuilderFromDet.h +++ b/RecoTracker/TkDetLayers/src/BladeShapeBuilderFromDet.h @@ -1,8 +1,8 @@ #ifndef RecoTracker_TkDetLayers_BladeShapeBuilderFromDet_h #define RecoTracker_TkDetLayers_BladeShapeBuilderFromDet_h -#include "BoundDiskSector.h" -#include "DiskSectorBounds.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/DiskSectorBounds.h" #include "Geometry/CommonDetUnit/interface/GeomDet.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include diff --git a/RecoTracker/TkDetLayers/src/BoundDiskSector.cc b/RecoTracker/TkDetLayers/src/BoundDiskSector.cc deleted file mode 100644 index de50f4945bc59..0000000000000 --- a/RecoTracker/TkDetLayers/src/BoundDiskSector.cc +++ /dev/null @@ -1 +0,0 @@ -#include "BoundDiskSector.h" diff --git a/RecoTracker/TkDetLayers/src/CompositeTECPetal.h b/RecoTracker/TkDetLayers/src/CompositeTECPetal.h index de86ffa075be3..b89bb2247a5f4 100644 --- a/RecoTracker/TkDetLayers/src/CompositeTECPetal.h +++ b/RecoTracker/TkDetLayers/src/CompositeTECPetal.h @@ -5,7 +5,7 @@ #include "TECWedge.h" -#include "BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" #include "SubLayerCrossings.h" diff --git a/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromDet.h b/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromDet.h index 0b85f3fd49e37..f95e4b30fa900 100644 --- a/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromDet.h +++ b/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromDet.h @@ -1,8 +1,8 @@ #ifndef RecoTracker_TkDetLayers_ForwardDiskSectorBuilderFromDet_h #define RecoTracker_TkDetLayers_ForwardDiskSectorBuilderFromDet_h -#include "BoundDiskSector.h" -#include "DiskSectorBounds.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/DiskSectorBounds.h" #include "Geometry/CommonDetUnit/interface/GeomDet.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include diff --git a/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromWedges.h b/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromWedges.h index b6e778ca4b530..b1741743461f7 100644 --- a/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromWedges.h +++ b/RecoTracker/TkDetLayers/src/ForwardDiskSectorBuilderFromWedges.h @@ -1,8 +1,8 @@ #ifndef RecoTracker_TkDetLayers_ForwardDiskSectorBuilderFromWedges_h #define RecoTracker_TkDetLayers_ForwardDiskSectorBuilderFromWedges_h -#include "BoundDiskSector.h" -#include "DiskSectorBounds.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/DiskSectorBounds.h" #include "TECWedge.h" #include #include diff --git a/RecoTracker/TkDetLayers/src/Phase1PixelBlade.h b/RecoTracker/TkDetLayers/src/Phase1PixelBlade.h index 0b4935d3172a4..012bda3a843ba 100644 --- a/RecoTracker/TkDetLayers/src/Phase1PixelBlade.h +++ b/RecoTracker/TkDetLayers/src/Phase1PixelBlade.h @@ -1,7 +1,7 @@ #ifndef TkDetLayers_Phase1PixelBlade_h #define TkDetLayers_Phase1PixelBlade_h -#include "BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" #include "TrackingTools/DetLayers/interface/GeometricSearchDet.h" #include "TrackingTools/DetLayers/interface/PeriodicBinFinderInZ.h" diff --git a/RecoTracker/TkDetLayers/src/PixelBlade.h b/RecoTracker/TkDetLayers/src/PixelBlade.h index 2d91ae373cbe0..b78b6382a8b65 100644 --- a/RecoTracker/TkDetLayers/src/PixelBlade.h +++ b/RecoTracker/TkDetLayers/src/PixelBlade.h @@ -1,7 +1,7 @@ #ifndef TkDetLayers_PixelBlade_h #define TkDetLayers_PixelBlade_h -#include "BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" #include "TrackingTools/DetLayers/interface/GeometricSearchDet.h" #include "TrackingTools/DetLayers/interface/PeriodicBinFinderInZ.h" diff --git a/RecoTracker/TkDetLayers/src/TECWedge.h b/RecoTracker/TkDetLayers/src/TECWedge.h index 02dc71cb3fc77..b37dd943bacf5 100644 --- a/RecoTracker/TkDetLayers/src/TECWedge.h +++ b/RecoTracker/TkDetLayers/src/TECWedge.h @@ -2,7 +2,7 @@ #define TkDetLayers_TECWedge_h #include "TrackingTools/DetLayers/interface/GeometricSearchDet.h" -#include "BoundDiskSector.h" +#include "DataFormats/GeometrySurface/interface/BoundDiskSector.h" /** A concrete implementation for TEC layer * built out of TECPetals diff --git a/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py b/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py index e59c556deba78..c3b26ccc8d283 100644 --- a/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py +++ b/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py @@ -34,8 +34,8 @@ SigmaClock = cms.double(0.015), # [ns] CorrelationCoefficient = cms.double(1.), SmearTimeForOOTtails = cms.bool(True), - Npe_to_pC = cms.double(0.016), # [pC] - Npe_to_V = cms.double(0.0064),# [V] + Npe_to_pC = cms.double(0.016), # [pC] + Npe_to_V = cms.double(0.0064),# [V] # n bits for the ADC adcNbits = cms.uint32(10), @@ -84,6 +84,9 @@ ) ) +from Configuration.Eras.Modifier_phase2_etlV4_cff import phase2_etlV4 +phase2_etlV4.toModify(_endcap_MTDDigitizer.DeviceSimulation, meVPerMIP = 0.001 ) + from Configuration.ProcessModifiers.premix_stage1_cff import premix_stage1 for _m in [_barrel_MTDDigitizer, _endcap_MTDDigitizer]: premix_stage1.toModify(_m, premixStage1 = True) @@ -97,4 +100,3 @@ barrelDigitizer = _barrel_MTDDigitizer, endcapDigitizer = _endcap_MTDDigitizer ) -