DetLayers/src/MTDSectorForwardDoubleLayer.cc

0001 //#define EDM_ML_DEBUG
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),  // add back later
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   // Cache chamber pointers (the basic components_)
0029   // and find extension in R and Z
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   // mostly copied from ForwardDetLayer, except propagates to closest surface,
0072   // not to center
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   // take into account the thickness of the layer
0083   float deltaR = surface().bounds().thickness() / 2. * std::abs(tan(myState.localDirection().theta()));
0084
0085   // take into account the error on the predicted state
0086   constexpr float nSigma = 3.;
0087   if (myState.hasError()) {
0088     LocalError err = myState.localError().positionError();
0089     // ignore correlation for the moment...
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   // standard implementation of compatibleDets() for class which have
0113   // groupedCompatibleDets implemented.
0114   // This code should be moved in a common place intead of being
0115   // copied many times.
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   // this should be fixed either in RecoMTD/MeasurementDet/MTDDetLayerMeasurements or
0133   // RecoMTD/DetLayers/MTDSectorForwardDoubleLayer
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   // test that each front z is less than each back z
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 }