TkMSParametrization/src/MSLayer.cc

0001 #include "RecoTracker/TkMSParametrization/interface/MSLayer.h"
0002 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0003 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0004 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0005 #include "MSLayersKeeper.h"
0006 #include "FWCore/Utilities/interface/Likely.h"
0007
0008 #include <iostream>
0009
0010 using namespace GeomDetEnumerators;
0011 using namespace std;
0012 template <class T>
0013 T sqr(T t) {
0014   return t * t;
0015 }
0016
0017 //----------------------------------------------------------------------
0018 ostream& operator<<(ostream& s, const MSLayer& l) {
0019   s << " face: " << l.face() << " pos:" << l.position() << ", "
0020     << " range:" << l.range() << ", " << l.theX0Data;
0021   return s;
0022 }
0023 //----------------------------------------------------------------------
0024 ostream& operator<<(ostream& s, const MSLayer::DataX0& d) {
0025   if (d.hasX0)
0026     s << "x0=" << d.x0 << " sumX0D=" << d.sumX0D;
0027   else if (d.allLayers)
0028     s << "x0 by MSLayersKeeper";
0029   else
0030     s << "empty DataX0";
0031   return s;
0032 }
0033 //----------------------------------------------------------------------
0034 //MP
0035 MSLayer::MSLayer(const DetLayer* layer, const DataX0& dataX0)
0036     : theFace(layer->location()), theSeqNum(layer->seqNum()), theX0Data(dataX0) {
0037   const BarrelDetLayer* bl;
0038   const ForwardDetLayer* fl;
0039   theHalfThickness = layer->surface().bounds().thickness() / 2;
0040
0041   switch (theFace) {
0042     case barrel:
0043       bl = static_cast<const BarrelDetLayer*>(layer);
0044       thePosition = bl->specificSurface().radius();
0045       theRange = Range(-bl->surface().bounds().length() / 2, bl->surface().bounds().length() / 2);
0046       break;
0047     case endcap:
0048       fl = static_cast<const ForwardDetLayer*>(layer);
0049       thePosition = fl->position().z();
0050       theRange = Range(fl->specificSurface().innerRadius(), fl->specificSurface().outerRadius());
0051       break;
0052     default:
0053       // should throw or simimal
0054       cout << " ** MSLayer ** unknown part - will not work!" << endl;
0055       break;
0056   }
0057 }
0058 //----------------------------------------------------------------------
0059 MSLayer::MSLayer(Location part, float position, Range range, float halfThickness, const DataX0& dataX0)
0060     : theFace(part),
0061       thePosition(position),
0062       theRange(range),
0063       theHalfThickness(halfThickness),
0064       theSeqNum(-1),
0065       theX0Data(dataX0) {}
0066
0067 //----------------------------------------------------------------------
0068 bool MSLayer::operator==(const MSLayer& o) const {
0069   return theFace == o.theFace && std::abs(thePosition - o.thePosition) < 1.e-3f;
0070 }
0071 //----------------------------------------------------------------------
0072 bool MSLayer::operator<(const MSLayer& o) const {
0073   if (theFace == barrel && o.theFace == barrel)
0074     return thePosition < o.thePosition;
0075   else if (theFace == barrel && o.theFace == endcap)
0076     return thePosition < o.range().max();
0077   else if (theFace == endcap && o.theFace == endcap)
0078     return std::abs(thePosition) < std::abs(o.thePosition);
0079   else
0080     return range().max() < o.thePosition;
0081 }
0082
0083 //----------------------------------------------------------------------
0084 pair<PixelRecoPointRZ, bool> MSLayer::crossing(const PixelRecoLineRZ& line) const {
0085   const float eps = 1.e-5;
0086   bool inLayer = true;
0087   float value = (theFace == barrel) ? line.zAtR(thePosition) : line.rAtZ(thePosition);
0088   if (value > theRange.max()) {
0089     value = theRange.max() - eps;
0090     inLayer = false;
0091   } else if (value < theRange.min()) {
0092     value = theRange.min() + eps;
0093     inLayer = false;
0094   }
0095   float z = thePosition;
0096   if (theFace == barrel)
0097     std::swap(z, value);  // if barrel value is z
0098   return make_pair(PixelRecoPointRZ(value, z), inLayer);
0099 }
0100 pair<PixelRecoPointRZ, bool> MSLayer::crossing(const SimpleLineRZ& line) const {
0101   const float eps = 1.e-5;
0102   bool inLayer = true;
0103   float value = (theFace == barrel) ? line.zAtR(thePosition) : line.rAtZ(thePosition);
0104   if (value > theRange.max()) {
0105     value = theRange.max() - eps;
0106     inLayer = false;
0107   } else if (value < theRange.min()) {
0108     value = theRange.min() + eps;
0109     inLayer = false;
0110   }
0111   float z = thePosition;
0112   if (theFace == barrel)
0113     std::swap(z, value);  // if barrel value is z
0114   return make_pair(PixelRecoPointRZ(value, z), inLayer);
0115 }
0116
0117 //----------------------------------------------------------------------
0118 float MSLayer::distance2(const PixelRecoPointRZ& point) const {
0119   float u = (theFace == barrel) ? point.r() : point.z();
0120   float v = (theFace == barrel) ? point.z() : point.r();
0121
0122   float du = std::abs(u - thePosition);
0123   if (theRange.inside(v))
0124     return (du < theHalfThickness) ? 0.f : du * du;
0125
0126   float dv = (v > theRange.max()) ? v - theRange.max() : theRange.min() - v;
0127   return sqr(du) + sqr(dv);
0128 }
0129
0130 //----------------------------------------------------------------------
0131 float MSLayer::x0(float cotTheta) const {
0132   if LIKELY (theX0Data.hasX0) {
0133     float OverSinTheta = std::sqrt(1.f + cotTheta * cotTheta);
0134     return (theFace == barrel) ? theX0Data.x0 * OverSinTheta : theX0Data.x0 * OverSinTheta / std::abs(cotTheta);
0135   } else if (theX0Data.allLayers) {
0136     const MSLayer* dataLayer = theX0Data.allLayers->layers(cotTheta).findLayer(*this);
0137     if (dataLayer)
0138       return dataLayer->x0(cotTheta);
0139   }
0140   return 0.;
0141 }
0142
0143 //----------------------------------------------------------------------
0144 float MSLayer::sumX0D(float cotTheta) const {
0145   if LIKELY (theX0Data.hasX0) {
0146     switch (theFace) {
0147       case barrel:
0148         return theX0Data.sumX0D *
0149                std::sqrt(std::sqrt((1.f + cotTheta * cotTheta) / (1.f + theX0Data.cotTheta * theX0Data.cotTheta)));
0150       case endcap:
0151         return (theX0Data.hasFSlope)
0152                    ? theX0Data.sumX0D + theX0Data.slopeSumX0D * (1.f / cotTheta - 1.f / theX0Data.cotTheta)
0153                    : theX0Data.sumX0D;
0154       case invalidLoc:
0155         break;  // make gcc happy
0156     }
0157   } else if (theX0Data.allLayers) {
0158     const MSLayer* dataLayer = theX0Data.allLayers->layers(cotTheta).findLayer(*this);
0159     if (dataLayer)
0160       return dataLayer->sumX0D(cotTheta);
0161   }
0162   return 0.;
0163 }