TkDetLayers/src/CompositeTECPetal.cc

0001 #include "CompositeTECPetal.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "ForwardDiskSectorBuilderFromWedges.h"
0006 #include "LayerCrossingSide.h"
0007 #include "DetGroupMerger.h"
0008 #include "CompatibleDetToGroupAdder.h"
0009
0010 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0011 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0012 #include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
0013
0014 #include "TkDetUtil.h"
0015 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0016
0017 #include <algorithm>
0018 #include <functional>
0019 #include <iterator>
0020 #include <numeric>
0021
0022 using namespace std;
0023
0024 typedef GeometricSearchDet::DetWithState DetWithState;
0025
0026 namespace {
0027   namespace details {
0028
0029     struct Mean {
0030       float operator()(const GeometricSearchDet* a, const GeometricSearchDet* b) const {
0031         return 0.5 * (b->position().perp() + a->position().perp());
0032       }
0033       float operator()(float a, float b) const { return 0.5 * (b + a); }
0034     };
0035
0036     void fillBoundaries(std::vector<const TECWedge*> const& dets, std::vector<float>& boundaries) {
0037       boundaries.resize(dets.size());
0038       std::transform(dets.begin(),
0039                      dets.end(),
0040                      boundaries.begin(),
0041                      std::bind(&GlobalPoint::perp, std::bind(&GeometricSearchDet::position, std::placeholders::_1)));
0042       std::adjacent_difference(boundaries.begin(), boundaries.end(), boundaries.begin(), Mean());
0043     }
0044
0045     inline int findBin(std::vector<float> const& boundaries, float r) {
0046       return std::lower_bound(boundaries.begin() + 1, boundaries.end(), r) - boundaries.begin() - 1;
0047     }
0048
0049     void fillPars(std::vector<const TECWedge*> const& dets, std::vector<CompositeTECPetal::WedgePar>& pars) {
0050       for (auto gsdet : dets) {
0051         const BoundDiskSector& wedgeSector = static_cast<const BoundDiskSector&>(gsdet->surface());
0052         float wedgeMinZ = std::abs(wedgeSector.position().z()) - 0.5 * wedgeSector.bounds().thickness();
0053         float wedgeMaxZ = std::abs(wedgeSector.position().z()) + 0.5 * wedgeSector.bounds().thickness();
0054         float thetaWedgeMin = wedgeSector.innerRadius() / wedgeMaxZ;
0055         float thetaWedgeMax = wedgeSector.outerRadius() / wedgeMinZ;
0056         CompositeTECPetal::WedgePar apar = {gsdet->position().perp(), thetaWedgeMin, thetaWedgeMax};
0057         pars.push_back(apar);
0058       }
0059     }
0060
0061     inline bool overlap(const GlobalPoint& gpos, const CompositeTECPetal::WedgePar& wpar, float ymax) {
0062       // this method is just a duplication of overlapInR
0063       // adapeted for groupedCompatibleDets() needs
0064
0065       // assume "fixed theta window", i.e. margin in local y = r is changing linearly with z
0066       auto tsRadius = gpos.perp();
0067       auto rmin = std::max(0.f, tsRadius - ymax);
0068       auto zmax = std::abs(gpos.z()) + 10.f;  // add 10 cm contingency
0069       auto rmax = (tsRadius + ymax);
0070       auto zmin = std::abs(gpos.z()) - 10.f;
0071
0072       // do the theta regions overlap ?
0073
0074       return !((rmin > zmax * wpar.thetaMax) | (zmin * wpar.thetaMin > rmax));
0075     }
0076
0077   }  // namespace details
0078 }  // namespace
0079
0080 CompositeTECPetal::CompositeTECPetal(vector<const TECWedge*>& innerWedges, vector<const TECWedge*>& outerWedges)
0081     : GeometricSearchDet(true), theFrontComps(innerWedges), theBackComps(outerWedges) {
0082   theComps.assign(theFrontComps.begin(), theFrontComps.end());
0083   theComps.insert(theComps.end(), theBackComps.begin(), theBackComps.end());
0084
0085   details::fillBoundaries(theFrontComps, theFrontBoundaries);
0086   details::fillBoundaries(theBackComps, theBackBoundaries);
0087   details::fillPars(theFrontComps, theFrontPars);
0088   details::fillPars(theBackComps, theBackPars);
0089
0090   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0091     theBasicComps.insert(theBasicComps.end(), (**it).basicComponents().begin(), (**it).basicComponents().end());
0092   }
0093
0094   //the Wedge are already R ordered
0095   //sort( theWedges.begin(), theWedges.end(), DetLessR());
0096   //sort( theFrontWedges.begin(), theFrontWedges.end(), DetLessR() );
0097   //sort( theBackWedges.begin(), theBackWedges.end(), DetLessR() );
0098   vector<const TECWedge*> allWedges;
0099   allWedges.assign(innerWedges.begin(), innerWedges.end());
0100   allWedges.insert(allWedges.end(), outerWedges.begin(), outerWedges.end());
0101
0102   theDiskSector = ForwardDiskSectorBuilderFromWedges()(allWedges);
0103   theFrontSector = ForwardDiskSectorBuilderFromWedges()(innerWedges);
0104   theBackSector = ForwardDiskSectorBuilderFromWedges()(outerWedges);
0105
0106   //--------- DEBUG INFO --------------
0107   LogDebug("TkDetLayers") << "DEBUG INFO for CompositeTECPetal";
0108
0109   for (auto it = theFrontComps.begin(); it != theFrontComps.end(); it++) {
0110     LogDebug("TkDetLayers") << "frontWedge phi,z,r: " << (*it)->surface().position().phi() << " , "
0111                             << (*it)->surface().position().z() << " , " << (*it)->surface().position().perp();
0112   }
0113
0114   for (auto it = theBackComps.begin(); it != theBackComps.end(); it++) {
0115     LogDebug("TkDetLayers") << "backWedge phi,z,r: " << (*it)->surface().position().phi() << " , "
0116                             << (*it)->surface().position().z() << " , " << (*it)->surface().position().perp();
0117   }
0118   //-----------------------------------
0119 }
0120
0121 CompositeTECPetal::~CompositeTECPetal() {
0122   vector<const GeometricSearchDet*>::const_iterator i;
0123   for (i = theComps.begin(); i != theComps.end(); i++) {
0124     delete *i;
0125   }
0126 }
0127
0128 pair<bool, TrajectoryStateOnSurface> CompositeTECPetal::compatible(const TrajectoryStateOnSurface& ts,
0129                                                                    const Propagator&,
0130                                                                    const MeasurementEstimator&) const {
0131   edm::LogError("TkDetLayers") << "temporary dummy implementation of CompositeTECPetal::compatible()!!";
0132   return pair<bool, TrajectoryStateOnSurface>();
0133 }
0134
0135 void CompositeTECPetal::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0136                                                const Propagator& prop,
0137                                                const MeasurementEstimator& est,
0138                                                std::vector<DetGroup>& result) const {
0139   vector<DetGroup> closestResult;
0140   SubLayerCrossings crossings;
0141   crossings = computeCrossings(tsos, prop.propagationDirection());
0142   if (!crossings.isValid())
0143     return;
0144
0145   addClosest(tsos, prop, est, crossings.closest(), closestResult);
0146   LogDebug("TkDetLayers") << "in TECPetal, closestResult.size(): " << closestResult.size();
0147
0148   if (closestResult.empty()) {
0149     vector<DetGroup> nextResult;
0150     addClosest(tsos, prop, est, crossings.other(), nextResult);
0151     LogDebug("TkDetLayers") << "in TECPetal, nextResult.size(): " << nextResult.size();
0152     if (nextResult.empty())
0153       return;
0154
0155     DetGroupElement nextGel(nextResult.front().front());
0156     int crossingSide = LayerCrossingSide::endcapSide(nextGel.trajectoryState(), prop);
0157     DetGroupMerger::orderAndMergeTwoLevels(
0158         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0159   } else {
0160     DetGroupElement closestGel(closestResult.front().front());
0161     float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0162
0163     searchNeighbors(tsos, prop, est, crossings.closest(), window, closestResult, false);
0164
0165     vector<DetGroup> nextResult;
0166     searchNeighbors(tsos, prop, est, crossings.other(), window, nextResult, true);
0167
0168     int crossingSide = LayerCrossingSide::endcapSide(closestGel.trajectoryState(), prop);
0169     DetGroupMerger::orderAndMergeTwoLevels(
0170         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0171   }
0172 }
0173
0174 SubLayerCrossings CompositeTECPetal::computeCrossings(const TrajectoryStateOnSurface& startingState,
0175                                                       PropagationDirection propDir) const {
0176   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0177   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0178
0179   auto rho = startingState.transverseCurvature();
0180
0181   HelixForwardPlaneCrossing crossing(startPos, startDir, rho, propDir);
0182
0183   pair<bool, double> frontPath = crossing.pathLength(*theFrontSector);
0184   if (!frontPath.first)
0185     return SubLayerCrossings();
0186
0187   pair<bool, double> backPath = crossing.pathLength(*theBackSector);
0188   if (!backPath.first)
0189     return SubLayerCrossings();
0190
0191   GlobalPoint gFrontPoint(crossing.position(frontPath.second));
0192   GlobalPoint gBackPoint(crossing.position(backPath.second));
0193
0194   LogDebug("TkDetLayers") << "in TECPetal,front crossing : r,z,phi: (" << gFrontPoint.perp() << "," << gFrontPoint.z()
0195                           << "," << gFrontPoint.phi() << ")";
0196
0197   LogDebug("TkDetLayers") << "in TECPetal,back crossing r,z,phi: (" << gBackPoint.perp() << "," << gBackPoint.z() << ","
0198                           << gBackPoint.phi() << ")";
0199
0200   int frontIndex = findBin(gFrontPoint.perp(), 0);
0201   SubLayerCrossing frontSLC(0, frontIndex, gFrontPoint);
0202
0203   int backIndex = findBin(gBackPoint.perp(), 1);
0204   SubLayerCrossing backSLC(1, backIndex, gBackPoint);
0205
0206   auto frontDist = std::abs(theFrontPars[frontIndex].theR - gFrontPoint.perp());
0207   auto backDist = std::abs(theBackPars[backIndex].theR - gBackPoint.perp());
0208
0209   // 0ss: frontDisk has index=0, backDisk has index=1
0210   if (frontDist < backDist) {
0211     return SubLayerCrossings(frontSLC, backSLC, 0);
0212   } else {
0213     return SubLayerCrossings(backSLC, frontSLC, 1);
0214   }
0215 }
0216
0217 bool CompositeTECPetal::addClosest(const TrajectoryStateOnSurface& tsos,
0218                                    const Propagator& prop,
0219                                    const MeasurementEstimator& est,
0220                                    const SubLayerCrossing& crossing,
0221                                    vector<DetGroup>& result) const {
0222   auto det = subLayer(crossing.subLayerIndex())[crossing.closestDetIndex()];
0223
0224   LogDebug("TkDetLayers") << "in TECPetal, adding Wedge at r,z,phi: (" << det->position().perp() << ","
0225                           << det->position().z() << "," << det->position().phi() << ")";
0226   LogDebug("TkDetLayers") << "wedge comps size: " << det->basicComponents().size();
0227
0228   return CompatibleDetToGroupAdder::add(*det, tsos, prop, est, result);
0229 }
0230
0231 void CompositeTECPetal::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0232                                         const Propagator& prop,
0233                                         const MeasurementEstimator& est,
0234                                         const SubLayerCrossing& crossing,
0235                                         float window,
0236                                         vector<DetGroup>& result,
0237                                         bool checkClosest) const {
0238   const GlobalPoint& gCrossingPos = crossing.position();
0239
0240   int closestIndex = crossing.closestDetIndex();
0241   int negStartIndex = closestIndex - 1;
0242   int posStartIndex = closestIndex + 1;
0243
0244   float detR = findPar(closestIndex, crossing.subLayerIndex()).theR;
0245
0246   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0247     if (gCrossingPos.perp2() < detR * detR) {
0248       posStartIndex = closestIndex;
0249     } else {
0250       negStartIndex = closestIndex;
0251     }
0252   }
0253
0254   const std::vector<const TECWedge*>& sLayer = subLayer(crossing.subLayerIndex());
0255
0256   //const BinFinderType& binFinder = (crossing.subLayerIndex()==0 ? theFrontBinFinder : theBackBinFinder);
0257   int theSize = crossing.subLayerIndex() == 0 ? theFrontComps.size() : theBackComps.size();
0258
0259   typedef CompatibleDetToGroupAdder Adder;
0260   for (int idet = negStartIndex; idet >= 0; idet--) {
0261     //if(idet<0 || idet>= theSize) {edm::LogInfo(TkDetLayers) << "===== error! gone out vector bounds.idet: " << idet ;exit;}
0262     const GeometricSearchDet& neighborWedge = *sLayer[idet];
0263     WedgePar const& wpar = findPar(idet, crossing.subLayerIndex());
0264     if (!details::overlap(gCrossingPos, wpar, window))
0265       break;  // --- to check
0266     if (!Adder::add(neighborWedge, tsos, prop, est, result))
0267       break;
0268     // maybe also add shallow crossing angle test here???
0269   }
0270   for (int idet = posStartIndex; idet < theSize; idet++) {
0271     //if(idet<0 || idet>= theSize) {edm::LogInfo(TkDetLayers) << "===== error! gone out vector bounds.idet: " << idet ;exit;}
0272     const GeometricSearchDet& neighborWedge = *sLayer[idet];
0273     WedgePar const& wpar = findPar(idet, crossing.subLayerIndex());
0274     if (!details::overlap(gCrossingPos, wpar, window))
0275       break;  // ---- to check
0276     if (!Adder::add(neighborWedge, tsos, prop, est, result))
0277       break;
0278     // maybe also add shallow crossing angle test here???
0279   }
0280 }
0281
0282 float CompositeTECPetal::computeWindowSize(const GeomDet* det,
0283                                            const TrajectoryStateOnSurface& tsos,
0284                                            const MeasurementEstimator& est) {
0285   return est.maximalLocalDisplacement(tsos, det->surface()).y();
0286 }
0287
0288 int CompositeTECPetal::findBin(float R, int diskSectorType) const {
0289   return details::findBin(diskSectorType == 0 ? theFrontBoundaries : theBackBoundaries, R);
0290 }