TkHitPairs/src/InnerDeltaPhi.cc

0001 #include "RecoTracker/TkHitPairs/src/InnerDeltaPhi.h"
0002 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0003 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0004 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0005 #include "RecoTracker/TkMSParametrization/interface/PixelRecoPointRZ.h"
0006 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionBase.h"
0007 #include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h"
0008 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
0009 #include "DataFormats/Math/interface/ExtVec.h"
0010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0011
0012 #if !defined(__INTEL_COMPILER)
0013 #define USE_VECTORS_HERE
0014 #endif
0015
0016 using namespace std;
0017
0018 #ifdef VI_DEBUG
0019 namespace {
0020   struct Stat {
0021     float xmin = 1.1;
0022     float xmax = -1.1;
0023     int nt = 0;
0024     int nn = 0;
0025     int nl = 0;
0026
0027     ~Stat() { std::cout << "ASIN " << xmin << ',' << xmax << ',' << nt << ',' << nn << ',' << nl << std::endl; }
0028   };
0029
0030   Stat stat;
0031
0032 }  // namespace
0033
0034 namespace {
0035   template <class T>
0036   inline T cropped_asin(T x) {
0037     stat.nt++;
0038     if (x < 0.f)
0039       stat.nn++;
0040     if (x > 0.5f)
0041       stat.nl++;
0042     stat.xmin = std::min(x, stat.xmin);
0043     stat.xmax = std::max(x, stat.xmax);
0044
0045     return std::abs(x) <= 1 ? std::asin(x) : (x > 0 ? T(M_PI / 2) : -T(M_PI / 2));
0046   }
0047 }  // namespace
0048 #else  // for icc
0049 namespace {
0050   template <class T>
0051   inline T cropped_asin(T x) {
0052     return std::abs(x) <= 1 ? std::asin(x) : (x > 0 ? T(M_PI / 2) : -T(M_PI / 2));
0053   }
0054 }  // namespace
0055 #endif
0056
0057 namespace {
0058
0059   template <class T>
0060   inline T sqr(T t) {
0061     return t * t;
0062   }
0063
0064   // reasonable (5.e-4) only for |x|<0.7 then degrades..
0065   template <typename T>
0066   inline T f_asin07f(T x) {
0067     auto ret = 1.f + (x * x) * (0.157549798488616943359375f + (x * x) * 0.125192224979400634765625f);
0068
0069     return x * ret;
0070   }
0071
0072 }  // namespace
0073
0074 #include "DataFormats/Math/interface/approx_atan2.h"
0075
0076 namespace {
0077   inline float f_atan2f(float y, float x) { return unsafe_atan2f<7>(y, x); }
0078   template <typename V>
0079   inline float f_phi(V v) {
0080     return f_atan2f(v.y(), v.x());
0081   }
0082 }  // namespace
0083
0084 InnerDeltaPhi::InnerDeltaPhi(const DetLayer& outlayer,
0085                              const DetLayer& layer,
0086                              const TrackingRegion& region,
0087                              const MagneticField& field,
0088                              const MultipleScatteringParametrisationMaker& msmaker,
0089                              bool precise,
0090                              float extraTolerance)
0091     : innerIsBarrel(layer.isBarrel()),
0092       outerIsBarrel(outlayer.isBarrel()),
0093       thePrecise(precise),
0094       ol(outlayer.seqNum()),
0095       theROrigin(region.originRBound()),
0096       theRLayer(0),
0097       theThickness(0),
0098       theExtraTolerance(extraTolerance),
0099       theA(0),
0100       theB(0),
0101       theVtxZ(region.origin().z()),
0102       thePtMin(region.ptMin()),
0103       theVtx(region.origin().x(), region.origin().y()),
0104       sigma(msmaker.parametrisation(&layer)) {
0105   float zMinOrigin = theVtxZ - region.originZBound();
0106   float zMaxOrigin = theVtxZ + region.originZBound();
0107   theRCurvature = PixelRecoUtilities::bendingRadius(thePtMin, field);
0108
0109   if (innerIsBarrel)
0110     initBarrelLayer(layer);
0111   else
0112     initForwardLayer(layer, zMinOrigin, zMaxOrigin);
0113
0114   if (outerIsBarrel)
0115     initBarrelMS(outlayer);
0116   else
0117     initForwardMS(outlayer);
0118 }
0119
0120 void InnerDeltaPhi::initBarrelMS(const DetLayer& outLayer) {
0121   const BarrelDetLayer& bl = static_cast<const BarrelDetLayer&>(outLayer);
0122   float rLayer = bl.specificSurface().radius();
0123   auto zmax = 0.5f * outLayer.surface().bounds().length();
0124   PixelRecoPointRZ zero(0., 0.);
0125   PixelRecoPointRZ point1(rLayer, 0.);
0126   PixelRecoPointRZ point2(rLayer, zmax);
0127   auto scatt1 = 3.f * sigma(thePtMin, zero, point1, ol);
0128   auto scatt2 = 3.f * sigma(thePtMin, zero, point2, ol);
0129   theDeltaScatt = (scatt2 - scatt1) / zmax;
0130   theScatt0 = scatt1;
0131 }
0132
0133 void InnerDeltaPhi::initForwardMS(const DetLayer& outLayer) {
0134   const ForwardDetLayer& fl = static_cast<const ForwardDetLayer&>(outLayer);
0135   auto minR = fl.specificSurface().innerRadius();
0136   auto maxR = fl.specificSurface().outerRadius();
0137   auto layerZ = outLayer.position().z();
0138   // compute min and max multiple scattering correction
0139   PixelRecoPointRZ zero(0., theVtxZ);
0140   PixelRecoPointRZ point1(minR, layerZ);
0141   PixelRecoPointRZ point2(maxR, layerZ);
0142   auto scatt1 = 3.f * sigma(thePtMin, zero, point1, ol);
0143   auto scatt2 = 3.f * sigma(thePtMin, zero, point2, ol);
0144   theDeltaScatt = (scatt2 - scatt1) / (maxR - minR);
0145   theScatt0 = scatt1 - theDeltaScatt * minR;
0146 }
0147
0148 void InnerDeltaPhi::initBarrelLayer(const DetLayer& layer) {
0149   const BarrelDetLayer& bl = static_cast<const BarrelDetLayer&>(layer);
0150   float rLayer = bl.specificSurface().radius();
0151
0152   // the maximal delta phi will be for the innermost hits
0153   theThickness = layer.surface().bounds().thickness();
0154   theRLayer = rLayer - 0.5f * theThickness;
0155 }
0156
0157 void InnerDeltaPhi::initForwardLayer(const DetLayer& layer, float zMinOrigin, float zMaxOrigin) {
0158   const ForwardDetLayer& fl = static_cast<const ForwardDetLayer&>(layer);
0159   theRLayer = fl.specificSurface().innerRadius();
0160   float layerZ = layer.position().z();
0161   theThickness = layer.surface().bounds().thickness();
0162   float layerZmin = layerZ > 0 ? layerZ - 0.5f * theThickness : layerZ + 0.5f * theThickness;
0163   theB = layerZ > 0 ? zMaxOrigin : zMinOrigin;
0164   theA = layerZmin - theB;
0165 }
0166
0167 PixelRecoRange<float> InnerDeltaPhi::phiRange(const Point2D& hitXY, float hitZ, float errRPhi) const {
0168   float rLayer = theRLayer;
0169   Point2D crossing;
0170
0171   Point2D dHit = hitXY - theVtx;
0172   auto dHitmag = dHit.mag();
0173   float dLayer = 0.;
0174   float dL = 0.;
0175
0176   // track is crossing layer with angle such as:
0177   // this factor should be taken in computation of eror projection
0178   float cosCross = 0;
0179
0180   //
0181   // compute crossing of stright track with inner layer
0182   //
0183   if (!innerIsBarrel) {
0184     auto t = theA / (hitZ - theB);
0185     auto dt = std::abs(theThickness / (hitZ - theB));
0186     crossing = theVtx + t * dHit;
0187     rLayer = crossing.mag();
0188     dLayer = t * dHitmag;
0189     dL = dt * dHitmag;
0190     cosCross = std::abs(dHit.unit().dot(crossing.unit()));
0191   } else {
0192     //
0193     // compute crossing of track with layer
0194     // dHit - from VTX to outer hit
0195     // rLayer - layer radius
0196     // dLayer - distance from VTX to inner layer in direction of dHit
0197     // vect(rLayer) = vect(rVTX) + vect(dHit).unit * dLayer
0198     //     rLayer^2 = (vect(rVTX) + vect(dHit).unit * dLayer)^2 and we have square eqation for dLayer
0199     //
0200     // barrel case
0201     //
0202     auto vtxmag2 = theVtx.mag2();
0203     if (vtxmag2 < 1.e-10f) {
0204       dLayer = rLayer;
0205     } else {
0206       // there are cancellation here....
0207       double var_c = vtxmag2 - sqr(rLayer);
0208       double var_b = theVtx.dot(dHit.unit());
0209       double var_delta = sqr(var_b) - var_c;
0210       if (var_delta <= 0.)
0211         var_delta = 0;
0212       dLayer = -var_b + std::sqrt(var_delta);  //only the value along vector is OK.
0213     }
0214     crossing = theVtx + dHit.unit() * dLayer;
0215     cosCross = std::abs(dHit.unit().dot(crossing.unit()));
0216     dL = theThickness / cosCross;
0217   }
0218
0219 #ifdef USE_VECTORS_HERE
0220   cms_float32x4_t num{dHitmag, dLayer, theROrigin * (dHitmag - dLayer), 1.f};
0221   cms_float32x4_t den{2 * theRCurvature, 2 * theRCurvature, dHitmag * dLayer, 1.f};
0222   auto phis = f_asin07f(num / den);
0223   phis = phis * dLayer / (rLayer * cosCross);
0224   auto deltaPhi = std::abs(phis[0] - phis[1]);
0225   auto deltaPhiOrig = phis[2];
0226 #else
0227 #warning no vector!
0228   auto alphaHit = cropped_asin(dHitmag / (2 * theRCurvature));
0229   auto OdeltaPhi = std::abs(alphaHit - cropped_asin(dLayer / (2 * theRCurvature)));
0230   OdeltaPhi *= dLayer / (rLayer * cosCross);
0231   // compute additional delta phi due to origin radius
0232   auto OdeltaPhiOrig = cropped_asin(theROrigin * (dHitmag - dLayer) / (dHitmag * dLayer));
0233   OdeltaPhiOrig *= dLayer / (rLayer * cosCross);
0234   // std::cout << "dphi " << OdeltaPhi<<'/'<<OdeltaPhiOrig << ' ' << deltaPhi<<'/'<<deltaPhiOrig << std::endl;
0235
0236   auto deltaPhi = OdeltaPhi;
0237   auto deltaPhiOrig = OdeltaPhiOrig;
0238 #endif
0239
0240   // additinal angle due to not perpendicular stright line crossing  (for displaced beam)
0241   //  double dPhiCrossing = (cosCross > 0.9999) ? 0 : dL *  sqrt(1-sqr(cosCross))/ rLayer;
0242   Point2D crossing2 = theVtx + dHit.unit() * (dLayer + dL);
0243   auto phicross2 = f_phi(crossing2);
0244   auto phicross1 = f_phi(crossing);
0245   auto dphicross = phicross2 - phicross1;
0246   if (dphicross < -float(M_PI))
0247     dphicross += float(2 * M_PI);
0248   if (dphicross > float(M_PI))
0249     dphicross -= float(2 * M_PI);
0250   if (dphicross > float(M_PI / 2))
0251     dphicross = 0.;  // something wrong?
0252   phicross2 = phicross1 + dphicross;
0253
0254   // inner hit error taken as constant
0255   auto deltaPhiHit = theExtraTolerance / rLayer;
0256
0257   // outer hit error
0258   //   double deltaPhiHitOuter = errRPhi/rLayer;
0259   auto deltaPhiHitOuter = errRPhi / hitXY.mag();
0260
0261   auto margin = deltaPhi + deltaPhiOrig + deltaPhiHit + deltaPhiHitOuter;
0262
0263   if (thePrecise) {
0264     // add multiple scattering correction
0265
0266     /*
0267     PixelRecoPointRZ zero(0., theVtxZ);
0268     PixelRecoPointRZ point(hitXY.mag(), hitZ);
0269     auto scatt = 3.f*sigma(thePtMin,zero, point, ol);
0270     */
0271
0272     auto w = outerIsBarrel ? std::abs(hitZ) : hitXY.mag();
0273     auto nscatt = theScatt0 + theDeltaScatt * w;
0274
0275     // std::cout << "scatt " << (outerIsBarrel ? "B" : "F") << (innerIsBarrel ? "B " : "F ")
0276     //          << scatt << ' ' << nscatt << ' ' << nscatt/scatt << std::endl;
0277
0278     margin += nscatt / rLayer;
0279   }
0280
0281   return PixelRecoRange<float>(std::min(phicross1, phicross2) - margin, std::max(phicross1, phicross2) + margin);
0282 }