TkMSParametrization/src/MSLayersAtAngle.cc

0001
0002 #include "MSLayersAtAngle.h"
0003
0004 #include "RecoTracker/TkMSParametrization/interface/PixelRecoLineRZ.h"
0005 #include "DataFormats/Math/interface/approx_log.h"
0006
0007 using namespace std;
0008
0009 namespace {
0010   template <class T>
0011   inline T sqr(T t) {
0012     return t * t;
0013   }
0014 }  // namespace
0015
0016 void MSLayersAtAngle::init() {
0017   std::sort(theLayers.begin(), theLayers.end());
0018   int i = -1;
0019   indeces.clear();
0020   for (auto const& l : theLayers) {
0021     ++i;
0022     int sq = l.seqNum();
0023     if (sq < 0)
0024       continue;
0025     if (sq >= int(indeces.size()))
0026       indeces.resize(sq + 1, -1);
0027     indeces[sq] = i;
0028   }
0029 }
0030
0031 //------------------------------------------------------------------------------
0032 MSLayersAtAngle::MSLayersAtAngle(const vector<MSLayer>& layers) : theLayers(layers) { init(); }
0033 //------------------------------------------------------------------------------
0034 const MSLayer* MSLayersAtAngle::findLayer(const MSLayer& layer) const {
0035   vector<MSLayer>::const_iterator it = find(theLayers.begin(), theLayers.end(), layer);
0036   return it == theLayers.end() ? nullptr : &(*it);
0037 }
0038
0039 //------------------------------------------------------------------------------
0040 void MSLayersAtAngle::update(const MSLayer& layer) {
0041   vector<MSLayer>::iterator it = find(theLayers.begin(), theLayers.end(), layer);
0042   if (it == theLayers.end()) {
0043     theLayers.push_back(layer);
0044     init();
0045   } else {
0046     *it = layer;
0047   }
0048 }
0049
0050 //------------------------------------------------------------------------------
0051 float MSLayersAtAngle::sumX0D(const PixelRecoPointRZ& pointI, const PixelRecoPointRZ& pointO) const {
0052   LayerItr iO = findLayer(pointO, theLayers.begin(), theLayers.end());
0053   //  cout << "outer Layer: "<<*iO<<endl;
0054   LayerItr iI = findLayer(pointI, theLayers.begin(), iO);
0055   //  cout << "inner Layer: "<<*iI<<endl;
0056
0057   return sqrt(sum2RmRn(iI, iO, pointO.r(), SimpleLineRZ(pointI, pointO)));
0058 }
0059
0060 float MSLayersAtAngle::sumX0D(int il, int ol, const PixelRecoPointRZ& pointI, const PixelRecoPointRZ& pointO) const {
0061   // if layer not at this angle (WHY???) revert to slow comp
0062   if (il >= int(indeces.size()) || ol >= int(indeces.size()) || indeces[il] < 0 || indeces[ol] < 0)
0063     return sumX0D(pointI, pointO);
0064
0065   LayerItr iI = theLayers.begin() + indeces[il];
0066   LayerItr iO = theLayers.begin() + indeces[ol];
0067
0068   return sqrt(sum2RmRn(iI, iO, pointO.r(), SimpleLineRZ(pointI, pointO)));
0069 }
0070
0071 //------------------------------------------------------------------------------
0072
0073 static const bool doPrint = false;
0074
0075 float MSLayersAtAngle::sumX0D(const PixelRecoPointRZ& pointI,
0076                               const PixelRecoPointRZ& pointM,
0077                               const PixelRecoPointRZ& pointO) const {
0078   LayerItr iO = findLayer(pointO, theLayers.begin(), theLayers.end());
0079   LayerItr iI = findLayer(pointI, theLayers.begin(), iO);
0080   LayerItr iM = findLayer(pointM, iI, iO);
0081
0082   float drOI = pointO.r() - pointI.r();
0083   float drMO = pointO.r() - pointM.r();
0084   float drMI = pointM.r() - pointI.r();
0085
0086   SimpleLineRZ line(pointI, pointO);
0087   float sum2I = sum2RmRn(iI + 1, iM, pointI.r(), line);
0088   float sum2O = sum2RmRn(iM, iO, pointO.r(), line);
0089
0090   float sum = std::sqrt(sum2I * sqr(drMO) + sum2O * sqr(drMI)) / drOI;
0091
0092   // if (iI!=theLayers.begin() )
0093   // doPrint = ((*iM).seqNum()<0 || (*iO).seqNum()<0);
0094   if (doPrint)
0095     std::cout << "old " << (*iI).seqNum() << " " << iI - theLayers.begin() << ", " << (*iM).seqNum() << " "
0096               << iM - theLayers.begin() << ", " << (*iO).seqNum() << " " << iO - theLayers.begin() << "  " << pointI.r()
0097               << " " << pointI.z() << " " << sum << std::endl;
0098
0099   return sum;
0100 }
0101
0102 float MSLayersAtAngle::sumX0D(
0103     float zV, int il, int ol, const PixelRecoPointRZ& pointI, const PixelRecoPointRZ& pointO) const {
0104   PixelRecoPointRZ pointV(0.f, zV);
0105
0106   // if layer not at this angle (WHY???) revert to slow comp
0107   if (il >= int(indeces.size()) || ol >= int(indeces.size()) || indeces[il] < 0 || indeces[ol] < 0)
0108     return sumX0D(pointV, pointI, pointO);
0109
0110   LayerItr iI = theLayers.begin() + indeces[il];
0111   LayerItr iO = theLayers.begin() + indeces[ol];
0112
0113   float drOI = pointO.r();
0114   float drMO = pointO.r() - pointI.r();
0115   float drMI = pointI.r();
0116
0117   SimpleLineRZ line(pointV, pointO);
0118   float sum2I = sum2RmRn(theLayers.begin() + 1, iI, pointV.r(), line);
0119   float sum2O = sum2RmRn(iI, iO, pointO.r(), line);
0120
0121   float sum = std::sqrt(sum2I * sqr(drMO) + sum2O * sqr(drMI)) / drOI;
0122
0123   if (doPrint)
0124     std::cout << "new " << il << " " << (*iI).seqNum() << " " << iI - theLayers.begin() << ", " << ol << " "
0125               << (*iO).seqNum() << " " << iO - theLayers.begin() << "  " << zV << " " << sum << std::endl;
0126
0127   return sum;
0128 }
0129
0130 //------------------------------------------------------------------------------
0131 float MSLayersAtAngle::sum2RmRn(MSLayersAtAngle::LayerItr i1,
0132                                 MSLayersAtAngle::LayerItr i2,
0133                                 float rTarget,
0134                                 const SimpleLineRZ& line) const {
0135   float sum2 = 0.f;
0136   float cotTh = line.cotLine();
0137   for (LayerItr it = i1; it < i2; it++) {
0138     std::pair<PixelRecoPointRZ, bool> cross = it->crossing(line);
0139     if (cross.second) {
0140       float x0 = it->x0(cotTh);
0141       float dr = rTarget - cross.first.r();
0142       if (x0 > 1.e-5f)
0143         dr *= 1.f + 0.038f * unsafe_logf<2>(x0);
0144       sum2 += x0 * dr * dr;
0145     }
0146     //  cout << *it << " crossing: "<<cross.second<<endl;
0147   }
0148   return sum2;
0149 }
0150 //------------------------------------------------------------------------------
0151 MSLayersAtAngle::LayerItr MSLayersAtAngle::findLayer(const PixelRecoPointRZ& point,
0152                                                      MSLayersAtAngle::LayerItr ibeg,
0153                                                      MSLayersAtAngle::LayerItr iend) const {
0154   const float BIG = 99999.f;
0155   const float EPSILON = 1.e-8f;
0156   LayerItr theIt = ibeg;
0157   float dist = BIG;
0158   for (LayerItr it = ibeg; it < iend; it++) {
0159     float d = it->distance2(point);
0160     if (d < dist) {
0161       if (d < EPSILON)
0162         return it;
0163       dist = d;
0164       theIt = it;
0165     }
0166   }
0167   return theIt;
0168 }
0169
0170 //------------------------------------------------------------------------------
0171 void MSLayersAtAngle::print() const {
0172   for (LayerItr it = theLayers.begin(); it != theLayers.end(); it++)
0173     cout << *it << endl;
0174 }