File indexing completed on 2023-03-17 11:20:16
0001
0002
0003 #include <RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h>
0004 #include <RecoMTD/DetLayers/interface/MTDDetSector.h>
0005 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0006 #include <DataFormats/GeometrySurface/interface/SimpleDiskBounds.h>
0007 #include <TrackingTools/GeomPropagators/interface/Propagator.h>
0008 #include <TrackingTools/DetLayers/interface/MeasurementEstimator.h>
0009 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0010
0011 #include <algorithm>
0012 #include <iostream>
0013 #include <vector>
0014
0015 using namespace std;
0016
0017 MTDSectorForwardDoubleLayer::MTDSectorForwardDoubleLayer(const vector<const MTDDetSector*>& frontSectors,
0018 const vector<const MTDDetSector*>& backSectors)
0019 : ForwardDetLayer(true),
0020 theFrontLayer(frontSectors),
0021 theBackLayer(backSectors),
0022 theSectors(frontSectors),
0023 theComponents(),
0024 theBasicComponents() {
0025 theSectors.insert(theSectors.end(), backSectors.begin(), backSectors.end());
0026 theComponents = std::vector<const GeometricSearchDet*>(theSectors.begin(), theSectors.end());
0027
0028
0029
0030 for (const auto& isect : theSectors) {
0031 vector<const GeomDet*> tmp2 = isect->basicComponents();
0032 theBasicComponents.insert(theBasicComponents.end(), tmp2.begin(), tmp2.end());
0033 }
0034
0035 setSurface(computeSurface());
0036
0037 LogTrace("MTDDetLayers") << "Constructing MTDSectorForwardDoubleLayer: " << std::fixed << std::setw(14)
0038 << basicComponents().size() << " Dets " << std::setw(14) << theSectors.size() << " Sectors "
0039 << " Z: " << std::setw(14) << specificSurface().position().z() << " R1: " << std::setw(14)
0040 << specificSurface().innerRadius() << " R2: " << std::setw(14)
0041 << specificSurface().outerRadius();
0042
0043 selfTest();
0044 }
0045
0046 BoundDisk* MTDSectorForwardDoubleLayer::computeSurface() {
0047 const BoundDisk& frontDisk = theFrontLayer.specificSurface();
0048 const BoundDisk& backDisk = theBackLayer.specificSurface();
0049
0050 float rmin = min(frontDisk.innerRadius(), backDisk.innerRadius());
0051 float rmax = max(frontDisk.outerRadius(), backDisk.outerRadius());
0052 float zmin = frontDisk.position().z();
0053 float halfThickness = frontDisk.bounds().thickness() / 2.;
0054 zmin = (zmin > 0) ? zmin - halfThickness : zmin + halfThickness;
0055 float zmax = backDisk.position().z();
0056 halfThickness = backDisk.bounds().thickness() / 2.;
0057 zmax = (zmax > 0) ? zmax + halfThickness : zmax - halfThickness;
0058 float zPos = (zmax + zmin) / 2.;
0059 PositionType pos(0., 0., zPos);
0060 RotationType rot;
0061
0062 return new BoundDisk(pos, rot, new SimpleDiskBounds(rmin, rmax, zmin - zPos, zmax - zPos));
0063 }
0064
0065 bool MTDSectorForwardDoubleLayer::isInsideOut(const TrajectoryStateOnSurface& tsos) const {
0066 return tsos.globalPosition().basicVector().dot(tsos.globalMomentum().basicVector()) > 0;
0067 }
0068
0069 std::pair<bool, TrajectoryStateOnSurface> MTDSectorForwardDoubleLayer::compatible(
0070 const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0071
0072
0073
0074 bool insideOut = isInsideOut(startingState);
0075 const MTDSectorForwardLayer& closerLayer = (insideOut) ? theFrontLayer : theBackLayer;
0076 LogTrace("MTDDetLayers") << "MTDSectorForwardDoubleLayer::compatible is assuming inside-out direction: " << insideOut;
0077
0078 TrajectoryStateOnSurface myState = prop.propagate(startingState, closerLayer.specificSurface());
0079 if (!myState.isValid())
0080 return make_pair(false, myState);
0081
0082
0083 float deltaR = surface().bounds().thickness() / 2. * std::abs(tan(myState.localDirection().theta()));
0084
0085
0086 constexpr float nSigma = 3.;
0087 if (myState.hasError()) {
0088 LocalError err = myState.localError().positionError();
0089
0090 deltaR += nSigma * sqrt(err.xx() + err.yy());
0091 }
0092
0093 float zPos = (zmax() + zmin()) / 2.;
0094 SimpleDiskBounds tmp(rmin() - deltaR, rmax() + deltaR, zmin() - zPos, zmax() - zPos);
0095
0096 return make_pair(tmp.inside(myState.localPosition()), myState);
0097 }
0098
0099 vector<GeometricSearchDet::DetWithState> MTDSectorForwardDoubleLayer::compatibleDets(
0100 const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0101 vector<DetWithState> result;
0102 pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0103
0104 if (!compat.first) {
0105 LogTrace("MTDDetLayers") << " MTDSectorForwardDoubleLayer::compatibleDets: not compatible"
0106 << " (should not have been selected!)";
0107 return result;
0108 }
0109
0110 TrajectoryStateOnSurface& tsos = compat.second;
0111
0112
0113
0114
0115
0116 vector<DetGroup> vectorGroups(groupedCompatibleDets(tsos, prop, est));
0117 for (const auto& thisDG : vectorGroups) {
0118 for (const auto& thisDGE : thisDG) {
0119 result.emplace_back(DetWithState(thisDGE.det(), thisDGE.trajectoryState()));
0120 }
0121 }
0122
0123 return result;
0124 }
0125
0126 vector<DetGroup> MTDSectorForwardDoubleLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0127 const Propagator& prop,
0128 const MeasurementEstimator& est) const {
0129 vector<GeometricSearchDet::DetWithState> detWithStates1, detWithStates2;
0130
0131 LogTrace("MTDDetLayers") << "groupedCompatibleDets are currently given always in inside-out order";
0132
0133
0134
0135 detWithStates1 = theFrontLayer.compatibleDets(startingState, prop, est);
0136 detWithStates2 = theBackLayer.compatibleDets(startingState, prop, est);
0137
0138 vector<DetGroup> result;
0139 if (!detWithStates1.empty())
0140 result.push_back(DetGroup(detWithStates1));
0141 if (!detWithStates2.empty())
0142 result.push_back(DetGroup(detWithStates2));
0143 LogTrace("MTDDetLayers") << "DoubleLayer Compatible dets: " << result.size();
0144 return result;
0145 }
0146
0147 bool MTDSectorForwardDoubleLayer::isCrack(const GlobalPoint& gp) const {
0148 bool result = false;
0149 LogTrace("MTDDetLayers")
0150 << "MTDSectorForwardDoubleLayer::isCrack kept only for backward compatibility, no real implementation";
0151 return result;
0152 }
0153
0154 void MTDSectorForwardDoubleLayer::selfTest() const {
0155 const std::vector<const GeomDet*>& frontDets = theFrontLayer.basicComponents();
0156 const std::vector<const GeomDet*>& backDets = theBackLayer.basicComponents();
0157
0158
0159 float frontz(0.);
0160 float backz(1e3f);
0161 for (const auto& thisFront : frontDets) {
0162 float tmpz(std::abs(thisFront->surface().position().z()));
0163 if (tmpz > frontz) {
0164 frontz = tmpz;
0165 }
0166 }
0167 for (const auto& thisBack : backDets) {
0168 float tmpz(std::abs(thisBack->surface().position().z()));
0169 if (tmpz < backz) {
0170 backz = tmpz;
0171 }
0172 }
0173 assert(frontz < backz);
0174 }