Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:22:51 

0001 #include "PixelForwardLayerPhase1.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 
0005 #include "DataFormats/GeometrySurface/interface/BoundingBox.h"
0006 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0007 
0008 #include "TrackingTools/DetLayers/interface/simple_stat.h"
0009 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing2Order.h"
0010 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0011 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0012 
0013 #include "LayerCrossingSide.h"
0014 #include "DetGroupMerger.h"
0015 #include "CompatibleDetToGroupAdder.h"
0016 
0017 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0018 
0019 #include <algorithm>
0020 #include <climits>
0021 
0022 using namespace std;
0023 
0024 typedef GeometricSearchDet::DetWithState DetWithState;
0025 
0026 PixelForwardLayerPhase1::PixelForwardLayerPhase1(vector<const Phase1PixelBlade*>& blades)
0027     : ForwardDetLayer(true), theComps(blades.begin(), blades.end()) {
0028   // Assumes blades are ordered inner first then outer; and within each by phi
0029   // where we go 0 -> pi, and then -pi -> 0
0030   // we also need the number of inner blades for the offset in theComps vector
0031   //
0032   // this->specificSurface() not yet available so need to calculate average R
0033   // we need some way to flag if FPIX is made of an inner and outer disk
0034   // or probably need to change the way this is done, e.g. a smarter binFinder
0035   float theRmin = (*(theComps.begin()))->surface().position().perp();
0036   float theRmax = theRmin;
0037   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0038     theRmin = std::min(theRmin, (*it)->surface().position().perp());
0039     theRmax = std::max(theRmax, (*it)->surface().position().perp());
0040   }
0041   //float split_inner_outer_radius = 0.5 * (theRmin + theRmax);
0042   // force the splitting rdius to be 10 cm to cope also with the FPIX disks with only the outer ring
0043   float split_inner_outer_radius = 10.;
0044   _num_innerpanels = 0;
0045   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0046     if ((**it).surface().position().perp() <= split_inner_outer_radius)
0047       ++_num_innerpanels;
0048   }
0049   _num_outerpanels = theComps.size() - _num_innerpanels;
0050   edm::LogInfo("TkDetLayers") << " Rmin, Rmax, R_average = " << theRmin << ", " << theRmax << ", "
0051                               << split_inner_outer_radius << std::endl
0052                               << " num inner, outer disks = " << _num_innerpanels << ", " << _num_outerpanels
0053                               << std::endl;
0054 
0055   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0056     theBasicComps.insert(theBasicComps.end(), (**it).basicComponents().begin(), (**it).basicComponents().end());
0057   }
0058 
0059   //They should be already phi-ordered. TO BE CHECKED!!
0060   //sort( theBlades.begin(), theBlades.end(), PhiLess());
0061   setSurface(computeSurface());
0062 
0063   theBinFinder_inner = BinFinderType(theComps.front()->surface().position().phi(), _num_innerpanels);
0064   theBinFinder_outer = BinFinderType(theComps[_num_innerpanels]->surface().position().phi(), _num_outerpanels);
0065 
0066   //--------- DEBUG INFO --------------
0067   LogDebug("TkDetLayers") << "DEBUG INFO for PixelForwardLayerPhase1"
0068                           << "\n"
0069                           << "Num of inner and outer panels = " << _num_innerpanels << ", " << _num_outerpanels << "\n"
0070                           << "Based on phi separation for inner: " << theComps.front()->surface().position().phi()
0071                           << " and on phi separation for outer: "
0072                           << theComps[_num_innerpanels]->surface().position().phi()
0073                           << "PixelForwardLayerPhase1.surfcace.phi(): " << std::endl
0074                           << this->surface().position().phi() << "\n"
0075                           << "PixelForwardLayerPhase1.surfcace.z(): " << this->surface().position().z() << "\n"
0076                           << "PixelForwardLayerPhase1.surfcace.innerR(): " << this->specificSurface().innerRadius()
0077                           << "\n"
0078                           << "PixelForwardLayerPhase1.surfcace.outerR(): " << this->specificSurface().outerRadius();
0079 
0080   for (vector<const GeometricSearchDet*>::const_iterator it = theComps.begin(); it != theComps.end(); it++) {
0081     LogDebug("TkDetLayers") << "blades phi,z,r: " << (*it)->surface().position().phi() << " , "
0082                             << (*it)->surface().position().z() << " , " << (*it)->surface().position().perp();
0083     //for(vector<const GeomDet*>::const_iterator iu=(**it).basicComponents().begin();
0084     //    iu!=(**it).basicComponents().end();iu++){
0085     //  std::cout << "   basic component rawId = " << hex << (**iu).geographicalId().rawId() << dec <<std::endl;
0086     //}
0087   }
0088   //-----------------------------------
0089 }
0090 
0091 PixelForwardLayerPhase1::~PixelForwardLayerPhase1() {
0092   vector<const GeometricSearchDet*>::const_iterator i;
0093   for (i = theComps.begin(); i != theComps.end(); i++) {
0094     delete *i;
0095   }
0096 }
0097 
0098 void PixelForwardLayerPhase1::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0099                                                      const Propagator& prop,
0100                                                      const MeasurementEstimator& est,
0101                                                      std::vector<DetGroup>& result) const {
0102   std::vector<DetGroup> closestResult_inner;
0103   std::vector<DetGroup> closestResult_outer;
0104   std::vector<DetGroup> nextResult_inner;
0105   std::vector<DetGroup> nextResult_outer;
0106   std::vector<DetGroup> result_inner;
0107   std::vector<DetGroup> result_outer;
0108   int frontindex_inner = 0;
0109   int frontindex_outer = 0;
0110   SubTurbineCrossings crossings_inner;
0111   SubTurbineCrossings crossings_outer;
0112 
0113   if (_num_innerpanels > 0)
0114     crossings_inner = computeCrossings(tsos, prop.propagationDirection(), true);
0115   crossings_outer = computeCrossings(tsos, prop.propagationDirection(), false);
0116   //  if (!crossings_inner.isValid){
0117   //    edm::LogInfo("TkDetLayers") << "inner computeCrossings returns invalid in PixelForwardLayerPhase1::groupedCompatibleDets:";
0118   //    return;
0119   //  }
0120   if (!crossings_outer.isValid) {
0121     edm::LogInfo("TkDetLayers")
0122         << "outer computeCrossings returns invalid in PixelForwardLayerPhase1::groupedCompatibleDets:";
0123     return;
0124   }
0125 
0126   typedef CompatibleDetToGroupAdder Adder;
0127   if (crossings_inner.isValid) {
0128     Adder::add(
0129         *theComps[theBinFinder_inner.binIndex(crossings_inner.closestIndex)], tsos, prop, est, closestResult_inner);
0130 
0131     if (closestResult_inner.empty()) {
0132       Adder::add(*theComps[theBinFinder_inner.binIndex(crossings_inner.nextIndex)], tsos, prop, est, result_inner);
0133       frontindex_inner = crossings_inner.nextIndex;
0134     } else {
0135       if (Adder::add(
0136               *theComps[theBinFinder_inner.binIndex(crossings_inner.nextIndex)], tsos, prop, est, nextResult_inner)) {
0137         int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0138         int theHelicity = computeHelicity(theComps[theBinFinder_inner.binIndex(crossings_inner.closestIndex)],
0139                                           theComps[theBinFinder_inner.binIndex(crossings_inner.nextIndex)]);
0140         std::vector<DetGroup> tmp99 = closestResult_inner;
0141         DetGroupMerger::orderAndMergeTwoLevels(
0142             std::move(tmp99), std::move(nextResult_inner), result_inner, theHelicity, crossingSide);
0143         if (theHelicity == crossingSide)
0144           frontindex_inner = crossings_inner.closestIndex;
0145         else
0146           frontindex_inner = crossings_inner.nextIndex;
0147       } else {
0148         result_inner.swap(closestResult_inner);
0149         frontindex_inner = crossings_inner.closestIndex;
0150       }
0151     }
0152   }
0153   if (!closestResult_inner.empty()) {
0154     DetGroupElement closestGel(closestResult_inner.front().front());
0155     float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0156     searchNeighbors(tsos, prop, est, crossings_inner, window, result_inner, true);
0157   }
0158 
0159   //DetGroupElement closestGel( closestResult.front().front());
0160   //float window = computeWindowSize( closestGel.det(), closestGel.trajectoryState(), est);
0161   //float detWidth = closestGel.det()->surface().bounds().width();
0162   //if (crossings.nextDistance < detWidth + window) {
0163 
0164   Adder::add(*theComps[(theBinFinder_outer.binIndex(crossings_outer.closestIndex)) + _num_innerpanels],
0165              tsos,
0166              prop,
0167              est,
0168              closestResult_outer);
0169 
0170   if (closestResult_outer.empty()) {
0171     Adder::add(*theComps[theBinFinder_outer.binIndex(crossings_outer.nextIndex) + _num_innerpanels],
0172                tsos,
0173                prop,
0174                est,
0175                result_outer);
0176     frontindex_outer = crossings_outer.nextIndex;
0177   } else {
0178     if (Adder::add(*theComps[theBinFinder_outer.binIndex(crossings_outer.nextIndex) + _num_innerpanels],
0179                    tsos,
0180                    prop,
0181                    est,
0182                    nextResult_outer)) {
0183       int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0184       int theHelicity =
0185           computeHelicity(theComps[theBinFinder_outer.binIndex(crossings_outer.closestIndex) + _num_innerpanels],
0186                           theComps[theBinFinder_outer.binIndex(crossings_outer.nextIndex) + _num_innerpanels]);
0187       std::vector<DetGroup> tmp99 = closestResult_outer;
0188       DetGroupMerger::orderAndMergeTwoLevels(
0189           std::move(tmp99), std::move(nextResult_outer), result_outer, theHelicity, crossingSide);
0190       if (theHelicity == crossingSide)
0191         frontindex_outer = crossings_outer.closestIndex;
0192       else
0193         frontindex_outer = crossings_outer.nextIndex;
0194     } else {
0195       result_outer.swap(closestResult_outer);
0196       frontindex_outer = crossings_outer.closestIndex;
0197     }
0198   }
0199   if (!closestResult_outer.empty()) {
0200     DetGroupElement closestGel(closestResult_outer.front().front());
0201     float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0202     searchNeighbors(tsos, prop, est, crossings_outer, window, result_outer, false);
0203   }
0204 
0205   if (result_inner.empty() && result_outer.empty())
0206     return;
0207   if (result_inner.empty())
0208     result.swap(result_outer);
0209   else if (result_outer.empty())
0210     result.swap(result_inner);
0211   else {
0212     int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0213     int theHelicity = computeHelicity(theComps[theBinFinder_inner.binIndex(frontindex_inner)],
0214                                       theComps[theBinFinder_outer.binIndex(frontindex_outer) + _num_innerpanels]);
0215     DetGroupMerger::orderAndMergeTwoLevels(
0216         std::move(result_inner), std::move(result_outer), result, theHelicity, crossingSide);
0217   }
0218 }
0219 
0220 void PixelForwardLayerPhase1::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0221                                               const Propagator& prop,
0222                                               const MeasurementEstimator& est,
0223                                               const SubTurbineCrossings& crossings,
0224                                               float window,
0225                                               vector<DetGroup>& result,
0226                                               bool innerDisk) const {
0227   typedef CompatibleDetToGroupAdder Adder;
0228   int crossingSide = LayerCrossingSide().endcapSide(tsos, prop);
0229   typedef DetGroupMerger Merger;
0230 
0231   int negStart = min(crossings.closestIndex, crossings.nextIndex) - 1;
0232   int posStart = max(crossings.closestIndex, crossings.nextIndex) + 1;
0233 
0234   int quarter = (innerDisk ? _num_innerpanels / 4 : _num_outerpanels / 4);
0235 
0236   for (int idet = negStart; idet >= negStart - quarter + 1; idet--) {
0237     std::vector<DetGroup> tmp1;
0238     std::vector<DetGroup> newResult;
0239     if (innerDisk) {
0240       const GeometricSearchDet* neighbor = theComps[theBinFinder_inner.binIndex(idet)];
0241       // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0242       // maybe also add shallow crossing angle test here???
0243       if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0244         break;
0245       int theHelicity =
0246           computeHelicity(theComps[theBinFinder_inner.binIndex(idet)], theComps[theBinFinder_inner.binIndex(idet + 1)]);
0247       std::vector<DetGroup> tmp2;
0248       tmp2.swap(result);
0249       Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
0250     } else {
0251       const GeometricSearchDet* neighbor = theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels];
0252       // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0253       // maybe also add shallow crossing angle test here???
0254       if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0255         break;
0256       int theHelicity = computeHelicity(theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels],
0257                                         theComps[(theBinFinder_outer.binIndex(idet + 1)) + _num_innerpanels]);
0258       std::vector<DetGroup> tmp2;
0259       tmp2.swap(result);
0260       Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
0261     }
0262     result.swap(newResult);
0263   }
0264   for (int idet = posStart; idet < posStart + quarter - 1; idet++) {
0265     std::vector<DetGroup> tmp1;
0266     std::vector<DetGroup> newResult;
0267     if (innerDisk) {
0268       const GeometricSearchDet* neighbor = theComps[theBinFinder_inner.binIndex(idet)];
0269       // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0270       // maybe also add shallow crossing angle test here???
0271       if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0272         break;
0273       int theHelicity =
0274           computeHelicity(theComps[theBinFinder_inner.binIndex(idet - 1)], theComps[theBinFinder_inner.binIndex(idet)]);
0275       std::vector<DetGroup> tmp2;
0276       tmp2.swap(result);
0277       Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
0278     } else {
0279       const GeometricSearchDet* neighbor = theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels];
0280       // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
0281       // maybe also add shallow crossing angle test here???
0282       if (!Adder::add(*neighbor, tsos, prop, est, tmp1))
0283         break;
0284       int theHelicity = computeHelicity(theComps[(theBinFinder_outer.binIndex(idet - 1)) + _num_innerpanels],
0285                                         theComps[(theBinFinder_outer.binIndex(idet)) + _num_innerpanels]);
0286       std::vector<DetGroup> tmp2;
0287       tmp2.swap(result);
0288       Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
0289     }
0290     result.swap(newResult);
0291   }
0292 }
0293 
0294 int PixelForwardLayerPhase1::computeHelicity(const GeometricSearchDet* firstBlade,
0295                                              const GeometricSearchDet* secondBlade) {
0296   return std::abs(firstBlade->position().z()) < std::abs(secondBlade->position().z()) ? 0 : 1;
0297 }
0298 
0299 PixelForwardLayerPhase1::SubTurbineCrossings PixelForwardLayerPhase1::computeCrossings(
0300     const TrajectoryStateOnSurface& startingState, PropagationDirection propDir, bool innerDisk) const {
0301   typedef MeasurementEstimator::Local2DVector Local2DVector;
0302 
0303   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0304   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0305   auto rho = startingState.transverseCurvature();
0306 
0307   HelixArbitraryPlaneCrossing turbineCrossing(startPos, startDir, rho, propDir);
0308 
0309   pair<bool, double> thePath = turbineCrossing.pathLength(specificSurface());
0310 
0311   if (!thePath.first) {
0312     return SubTurbineCrossings();
0313   }
0314 
0315   HelixPlaneCrossing::PositionType turbinePoint(turbineCrossing.position(thePath.second));
0316   HelixPlaneCrossing::DirectionType turbineDir(turbineCrossing.direction(thePath.second));
0317   int closestIndex = 0;
0318   int nextIndex = 0;
0319   if (innerDisk)
0320     closestIndex = theBinFinder_inner.binIndex(turbinePoint.barePhi());
0321   else
0322     closestIndex = theBinFinder_outer.binIndex(turbinePoint.barePhi());
0323 
0324   HelixArbitraryPlaneCrossing2Order theBladeCrossing(turbinePoint, turbineDir, rho);
0325 
0326   float closestDist = 0;
0327   if (innerDisk) {
0328     const BoundPlane& closestPlane(static_cast<const BoundPlane&>(theComps[closestIndex]->surface()));
0329 
0330     pair<bool, double> theClosestBladePath = theBladeCrossing.pathLength(closestPlane);
0331     LocalPoint closestPos = closestPlane.toLocal(GlobalPoint(theBladeCrossing.position(theClosestBladePath.second)));
0332     closestDist = closestPos.x();
0333     // use fact that local X perp to global Y
0334     nextIndex = Geom::phiLess(closestPlane.phi(), turbinePoint.barePhi()) ? closestIndex + 1 : closestIndex - 1;
0335   } else {
0336     const BoundPlane& closestPlane(
0337         static_cast<const BoundPlane&>(theComps[closestIndex + _num_innerpanels]->surface()));
0338 
0339     pair<bool, double> theClosestBladePath = theBladeCrossing.pathLength(closestPlane);
0340     LocalPoint closestPos = closestPlane.toLocal(GlobalPoint(theBladeCrossing.position(theClosestBladePath.second)));
0341     closestDist = closestPos.x();
0342     // use fact that local X perp to global Y
0343     nextIndex = Geom::phiLess(closestPlane.phi(), turbinePoint.barePhi()) ? closestIndex + 1 : closestIndex - 1;
0344   }
0345 
0346   float nextDist = 0;
0347   if (innerDisk) {
0348     const BoundPlane& nextPlane(
0349         static_cast<const BoundPlane&>(theComps[theBinFinder_inner.binIndex(nextIndex)]->surface()));
0350     pair<bool, double> theNextBladePath = theBladeCrossing.pathLength(nextPlane);
0351     LocalPoint nextPos = nextPlane.toLocal(GlobalPoint(theBladeCrossing.position(theNextBladePath.second)));
0352     nextDist = nextPos.x();
0353   } else {
0354     const BoundPlane& nextPlane(
0355         static_cast<const BoundPlane&>(theComps[theBinFinder_outer.binIndex(nextIndex) + _num_innerpanels]->surface()));
0356     pair<bool, double> theNextBladePath = theBladeCrossing.pathLength(nextPlane);
0357     LocalPoint nextPos = nextPlane.toLocal(GlobalPoint(theBladeCrossing.position(theNextBladePath.second)));
0358     nextDist = nextPos.x();
0359   }
0360   if (std::abs(closestDist) < std::abs(nextDist)) {
0361     return SubTurbineCrossings(closestIndex, nextIndex, nextDist);
0362   } else {
0363     return SubTurbineCrossings(nextIndex, closestIndex, closestDist);
0364   }
0365 }
0366 
0367 float PixelForwardLayerPhase1::computeWindowSize(const GeomDet* det,
0368                                                  const TrajectoryStateOnSurface& tsos,
0369                                                  const MeasurementEstimator& est) {
0370   return est.maximalLocalDisplacement(tsos, det->surface()).x();
0371 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!