DetLayers/src/MTDSectorForwardLayer.cc

0001 //#define EDM_ML_DEBUG
0002
0003 #include <RecoMTD/DetLayers/interface/MTDSectorForwardLayer.h>
0004 #include <RecoMTD/DetLayers/interface/MTDDetSector.h>
0005 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0006 #include <DataFormats/GeometrySurface/interface/DiskSectorBounds.h>
0007 #include <TrackingTools/GeomPropagators/interface/Propagator.h>
0008 #include <TrackingTools/DetLayers/interface/MeasurementEstimator.h>
0009
0010 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0011
0012 #include <algorithm>
0013 #include <iostream>
0014 #include <vector>
0015
0016 using namespace std;
0017
0018 MTDSectorForwardLayer::MTDSectorForwardLayer(const vector<const MTDDetSector*>& sectors)
0019     : ForwardDetLayer(false), theSectors(sectors), theComponents(theSectors.begin(), theSectors.end()) {
0020   // Initial values for R, Z and Phi bounds
0021   float theRmin = sectors.front()->basicComponents().front()->position().perp();
0022   float theRmax = theRmin;
0023   float theZmin = sectors.front()->position().z();
0024   float theZmax = theZmin;
0025
0026   // Cache chamber pointers (the basic components_)
0027   // and find extension in R and Z
0028   for (const auto& isect : sectors) {
0029     vector<const GeomDet*> tmp2 = isect->basicComponents();
0030     theBasicComps.insert(theBasicComps.end(), tmp2.begin(), tmp2.end());
0031
0032     theRmin = min(theRmin, isect->specificSurface().innerRadius());
0033     theRmax = max(theRmax, isect->specificSurface().outerRadius());
0034     float halfThick = isect->surface().bounds().thickness() / 2.;
0035     float zCenter = isect->surface().position().z();
0036     theZmin = min(theZmin, zCenter - halfThick);
0037     theZmax = max(theZmax, zCenter + halfThick);
0038   }
0039
0040   // Build surface
0041
0042   float zPos = (theZmax + theZmin) / 2.;
0043   PositionType pos(0., 0., zPos);
0044   RotationType rot;
0045
0046   setSurface(new BoundDisk(pos, rot, new SimpleDiskBounds(theRmin, theRmax, theZmin - zPos, theZmax - zPos)));
0047
0048   LogTrace("MTDDetLayers") << "Constructing MTDSectorForwardLayer: " << std::fixed << std::setw(14)
0049                            << basicComponents().size() << " Dets, " << std::setw(14) << theSectors.size()
0050                            << " Sectors, "
0051                            << " Z: " << std::setw(14) << specificSurface().position().z() << " R1: " << std::setw(14)
0052                            << specificSurface().innerRadius() << " R2: " << std::setw(14)
0053                            << specificSurface().outerRadius();
0054 }
0055
0056 MTDSectorForwardLayer::~MTDSectorForwardLayer() {
0057   for (auto& i : theSectors) {
0058     delete i;
0059   }
0060 }
0061
0062 vector<GeometricSearchDet::DetWithState> MTDSectorForwardLayer::compatibleDets(
0063     const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0064   vector<DetWithState> result;
0065
0066   LogTrace("MTDDetLayers") << "MTDSectorForwardLayer::compatibleDets,"
0067                            << " R1 " << std::fixed << std::setw(14) << specificSurface().innerRadius()
0068                            << " R2: " << std::setw(14) << specificSurface().outerRadius()
0069                            << " FTS at R: " << std::setw(14) << startingState.globalPosition().perp();
0070
0071   pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0072
0073   if (!compat.first) {
0074     LogTrace("MTDDetLayers") << "     MTDSectorForwardLayer::compatibleDets: not compatible"
0075                              << " (should not have been selected!)";
0076     return result;
0077   }
0078
0079   TrajectoryStateOnSurface& tsos = compat.second;
0080
0081   // as there are either two or four sectors only, avoid complex logic and just loop on all of them
0082
0083   // Use state on layer surface. Note that local coordinates and errors
0084   // are the same on the layer and on all sectors surfaces, since
0085   // all BoundDisks are centered in 0,0 and have the same rotation.
0086   // CAVEAT: if the sectors are not at the same Z, the local position and error
0087   // will be "Z-projected" to the sectors. This is a fairly good approximation.
0088   // However in this case additional propagation will be done when calling
0089   // compatibleDets.
0090   GlobalPoint startPos = tsos.globalPosition();
0091
0092   for (unsigned int isect = 0; isect < theSectors.size(); isect++) {
0093     LocalPoint nextPos(theSectors[isect]->specificSurface().toLocal(startPos));
0094     LogDebug("MTDDetLayers") << "Global point = " << std::fixed << startPos << " local point = " << nextPos
0095                              << " global sector ref pos = " << theSectors[isect]->specificSurface().position();
0096     bool inside = false;
0097     if (tsos.hasError()) {
0098       inside = theSectors[isect]->specificSurface().bounds().inside(nextPos, tsos.localError().positionError(), 1.);
0099     } else {
0100       inside = theSectors[isect]->specificSurface().bounds().inside(nextPos);
0101     }
0102     if (inside) {
0103 #ifdef EDM_ML_DEBUG
0104       LogTrace("MTDDetLayers") << "     MTDSectorForwardLayer::fastCompatibleDets:NextSector " << std::fixed
0105                                << std::setw(14) << isect << "\n"
0106                                << (*theSectors[isect]) << "\n FTS at Z,R,phi: " << std::setw(14)
0107                                << tsos.globalPosition().z() << " , " << std::setw(14) << tsos.globalPosition().perp()
0108                                << "," << std::setw(14) << tsos.globalPosition().phi();
0109       if (tsos.hasError()) {
0110         LogTrace("MTDDetLayers") << " sR: " << sqrt(tsos.localError().positionError().yy())
0111                                  << " sX: " << sqrt(tsos.localError().positionError().xx());
0112       }
0113 #endif
0114       vector<DetWithState> nextRodDets(theSectors[isect]->compatibleDets(tsos, prop, est));
0115       if (!nextRodDets.empty()) {
0116         result.insert(result.end(), nextRodDets.begin(), nextRodDets.end());
0117       } else {
0118         break;
0119       }
0120     }
0121   }
0122
0123   LogTrace("MTDDetLayers") << "     MTDSectorForwardLayer::fastCompatibleDets: found: " << result.size();
0124
0125   return result;
0126 }
0127
0128 vector<DetGroup> MTDSectorForwardLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0129                                                               const Propagator& prop,
0130                                                               const MeasurementEstimator& est) const {
0131   // FIXME should return only 1 group
0132   edm::LogInfo("MTDDetLayers") << "dummy implementation of MTDSectorForwardLayer::groupedCompatibleDets()";
0133   return vector<DetGroup>();
0134 }
0135
0136 GeomDetEnumerators::SubDetector MTDSectorForwardLayer::subDetector() const {
0137   return theBasicComps.front()->subDetector();
0138 }
0139
0140 const vector<const GeometricSearchDet*>& MTDSectorForwardLayer::components() const { return theComponents; }