DetLayers/src/MTDRingForwardDoubleLayer.cc

0001 //#define EDM_ML_DEBUG
0002
0003 /** \file
0004  *
0005  *  \author L. Gray
0006  */
0007
0008 #include <RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h>
0009 #include <RecoMTD/DetLayers/interface/MTDDetRing.h>
0010 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0011 #include <DataFormats/GeometrySurface/interface/SimpleDiskBounds.h>
0012 #include <TrackingTools/GeomPropagators/interface/Propagator.h>
0013 #include <TrackingTools/DetLayers/interface/MeasurementEstimator.h>
0014 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0015
0016 #include <DataFormats/ForwardDetId/interface/ETLDetId.h>
0017
0018 #include <algorithm>
0019 #include <iostream>
0020 #include <vector>
0021
0022 using namespace std;
0023
0024 MTDRingForwardDoubleLayer::MTDRingForwardDoubleLayer(const vector<const ForwardDetRing*>& frontRings,
0025                                                      const vector<const ForwardDetRing*>& backRings)
0026     : RingedForwardLayer(true),
0027       theFrontLayer(frontRings),
0028       theBackLayer(backRings),
0029       theRings(frontRings),  // add back later
0030       theComponents(),
0031       theBasicComponents() {
0032   theRings.insert(theRings.end(), backRings.begin(), backRings.end());
0033   theComponents = std::vector<const GeometricSearchDet*>(theRings.begin(), theRings.end());
0034
0035   // Cache chamber pointers (the basic components_)
0036   // and find extension in R and Z
0037   for (vector<const ForwardDetRing*>::const_iterator it = theRings.begin(); it != theRings.end(); it++) {
0038     vector<const GeomDet*> tmp2 = (*it)->basicComponents();
0039     theBasicComponents.insert(theBasicComponents.end(), tmp2.begin(), tmp2.end());
0040   }
0041
0042   setSurface(computeSurface());
0043
0044   LogTrace("MTDDetLayers") << "Constructing MTDRingForwardDoubleLayer: " << basicComponents().size() << " Dets "
0045                            << theRings.size() << " Rings "
0046                            << " Z: " << specificSurface().position().z() << " R1: " << specificSurface().innerRadius()
0047                            << " R2: " << specificSurface().outerRadius();
0048
0049   selfTest();
0050 }
0051
0052 BoundDisk* MTDRingForwardDoubleLayer::computeSurface() {
0053   const BoundDisk& frontDisk = theFrontLayer.specificSurface();
0054   const BoundDisk& backDisk = theBackLayer.specificSurface();
0055
0056   float rmin = min(frontDisk.innerRadius(), backDisk.innerRadius());
0057   float rmax = max(frontDisk.outerRadius(), backDisk.outerRadius());
0058   float zmin = frontDisk.position().z();
0059   float halfThickness = frontDisk.bounds().thickness() / 2.;
0060   zmin = (zmin > 0) ? zmin - halfThickness : zmin + halfThickness;
0061   float zmax = backDisk.position().z();
0062   halfThickness = backDisk.bounds().thickness() / 2.;
0063   zmax = (zmax > 0) ? zmax + halfThickness : zmax - halfThickness;
0064   float zPos = (zmax + zmin) / 2.;
0065   PositionType pos(0., 0., zPos);
0066   RotationType rot;
0067
0068   return new BoundDisk(pos, rot, new SimpleDiskBounds(rmin, rmax, zmin - zPos, zmax - zPos));
0069 }
0070
0071 bool MTDRingForwardDoubleLayer::isInsideOut(const TrajectoryStateOnSurface& tsos) const {
0072   return tsos.globalPosition().basicVector().dot(tsos.globalMomentum().basicVector()) > 0;
0073 }
0074
0075 std::pair<bool, TrajectoryStateOnSurface> MTDRingForwardDoubleLayer::compatible(
0076     const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0077   // mostly copied from ForwardDetLayer, except propagates to closest surface,
0078   // not to center
0079
0080   bool insideOut = isInsideOut(startingState);
0081   const MTDRingForwardLayer& closerLayer = (insideOut) ? theFrontLayer : theBackLayer;
0082   LogTrace("MTDDetLayers") << "MTDRingForwardDoubleLayer::compatible is assuming inside-out direction: " << insideOut;
0083
0084   TrajectoryStateOnSurface myState = prop.propagate(startingState, closerLayer.specificSurface());
0085   if (!myState.isValid())
0086     return make_pair(false, myState);
0087
0088   // take into account the thickness of the layer
0089   float deltaR = surface().bounds().thickness() / 2. * fabs(tan(myState.localDirection().theta()));
0090
0091   // take into account the error on the predicted state
0092   const float nSigma = 3.;
0093   if (myState.hasError()) {
0094     LocalError err = myState.localError().positionError();
0095     // ignore correlation for the moment...
0096     deltaR += nSigma * sqrt(err.xx() + err.yy());
0097   }
0098
0099   float zPos = (zmax() + zmin()) / 2.;
0100   SimpleDiskBounds tmp(rmin() - deltaR, rmax() + deltaR, zmin() - zPos, zmax() - zPos);
0101
0102   return make_pair(tmp.inside(myState.localPosition()), myState);
0103 }
0104
0105 vector<GeometricSearchDet::DetWithState> MTDRingForwardDoubleLayer::compatibleDets(
0106     const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0107   vector<DetWithState> result;
0108   pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0109
0110   if (!compat.first) {
0111     LogTrace("MTDDetLayers") << "     MTDRingForwardDoubleLayer::compatibleDets: not compatible"
0112                              << " (should not have been selected!)";
0113     return result;
0114   }
0115
0116   TrajectoryStateOnSurface& tsos = compat.second;
0117
0118   // standard implementation of compatibleDets() for class which have
0119   // groupedCompatibleDets implemented.
0120   // This code should be moved in a common place intead of being
0121   // copied many times.
0122   vector<DetGroup> vectorGroups = groupedCompatibleDets(tsos, prop, est);
0123   for (const auto& thisDG : vectorGroups) {
0124     for (const auto& thisDGE : thisDG) {
0125       result.emplace_back(DetWithState(thisDGE.det(), thisDGE.trajectoryState()));
0126     }
0127   }
0128   return result;
0129 }
0130
0131 vector<DetGroup> MTDRingForwardDoubleLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0132                                                                   const Propagator& prop,
0133                                                                   const MeasurementEstimator& est) const {
0134   vector<GeometricSearchDet::DetWithState> detWithStates1, detWithStates2;
0135
0136   LogTrace("MTDDetLayers") << "groupedCompatibleDets are currently given always in inside-out order";
0137   // this should be fixed either in RecoMTD/MeasurementDet/MTDDetLayerMeasurements or
0138   // RecoMTD/DetLayers/MTDRingForwardDoubleLayer
0139
0140   detWithStates1 = theFrontLayer.compatibleDets(startingState, prop, est);
0141   detWithStates2 = theBackLayer.compatibleDets(startingState, prop, est);
0142
0143   vector<DetGroup> result;
0144   if (!detWithStates1.empty())
0145     result.push_back(DetGroup(detWithStates1));
0146   if (!detWithStates2.empty())
0147     result.push_back(DetGroup(detWithStates2));
0148   LogTrace("MTDDetLayers") << "DoubleLayer Compatible dets: " << result.size();
0149   return result;
0150 }
0151
0152 bool MTDRingForwardDoubleLayer::isCrack(const GlobalPoint& gp) const {
0153   // approximate
0154   bool result = false;
0155   double r = gp.perp();
0156   const std::vector<const ForwardDetRing*>& backRings = theBackLayer.rings();
0157   if (backRings.size() > 1) {
0158     const MTDDetRing* innerRing = dynamic_cast<const MTDDetRing*>(backRings[0]);
0159     const MTDDetRing* outerRing = dynamic_cast<const MTDDetRing*>(backRings[1]);
0160     assert(innerRing && outerRing);
0161     float crackInner = innerRing->specificSurface().outerRadius();
0162     float crackOuter = outerRing->specificSurface().innerRadius();
0163     LogTrace("MTDDetLayers") << "In a crack:" << crackInner << " " << r << " " << crackOuter;
0164     if (r > crackInner && r < crackOuter)
0165       return true;
0166   }
0167   // non-overlapping rings
0168   return result;
0169 }
0170
0171 void MTDRingForwardDoubleLayer::selfTest() const {
0172   const std::vector<const GeomDet*>& frontDets = theFrontLayer.basicComponents();
0173   const std::vector<const GeomDet*>& backDets = theBackLayer.basicComponents();
0174
0175   for (int iring = 0; iring <= ETLDetId::kETLv1maxRing; iring++) {
0176     float frontz(0.);
0177     float backz(1e3f);
0178     for (const auto& thisFront : frontDets) {
0179       if (static_cast<ETLDetId>(thisFront->geographicalId().rawId()).mtdRR() == iring) {
0180         float tmpz(std::abs(thisFront->surface().position().z()));
0181         if (tmpz > frontz) {
0182           frontz = tmpz;
0183         }
0184       }
0185     }
0186     for (const auto& thisBack : backDets) {
0187       if (static_cast<ETLDetId>(thisBack->geographicalId().rawId()).mtdRR() == iring) {
0188         float tmpz(std::abs(thisBack->surface().position().z()));
0189         if (tmpz < backz) {
0190           backz = tmpz;
0191         }
0192       }
0193     }
0194     assert(frontz < backz);
0195   }
0196 }