TkMSParametrization/src/MSLayersKeeperX0AtEta.cc

0001 #include "MSLayersKeeperX0AtEta.h"
0002 #include "MultipleScatteringX0Data.h"
0003 #include "MultipleScatteringGeometry.h"
0004 #include <algorithm>
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   // float unsafe_asinhf(float x) { return std::abs(x)>10.f ? std::copysign(3.f,x) : unsafe_logf<3>(x+std::sqrt(1.f+x*x)); }
0015   inline float unsafe_asinhf(float x) { return unsafe_logf<3>(x + std::sqrt(1.f + x * x)); }
0016 }  // namespace
0017
0018 //------------------------------------------------------------------------------
0019 const MSLayersAtAngle& MSLayersKeeperX0AtEta::layers(float cotTheta) const {
0020   float eta = unsafe_asinhf(cotTheta);
0021   return theLayersData[idxBin(eta)];
0022 }
0023
0024 //------------------------------------------------------------------------------
0025 float MSLayersKeeperX0AtEta::eta(int idxBin) const { return (idxBin + 0.5 - theHalfNBins) * theDeltaEta; }
0026
0027 //------------------------------------------------------------------------------
0028 int MSLayersKeeperX0AtEta::idxBin(float eta) const {
0029   float ieta = eta / theDeltaEta;
0030   if (std::abs(ieta) >= theHalfNBins - 1.e-3)
0031     return (eta > 0) ? max(2 * theHalfNBins - 1, 0) : 0;
0032   else
0033     return int(ieta + theHalfNBins);
0034 }
0035
0036 //------------------------------------------------------------------------------
0037 MSLayersKeeperX0AtEta::MSLayersKeeperX0AtEta(const GeometricSearchTracker& tracker, const MagneticField& bfield) {
0038   const float BIG = 99999.;
0039
0040   // set size from data file
0041   MultipleScatteringX0Data msX0data;
0042   theHalfNBins = msX0data.nBinsEta();
0043   float etaMax = msX0data.maxEta();
0044   theDeltaEta = (theHalfNBins != 0) ? etaMax / theHalfNBins : BIG;
0045
0046   theLayersData = vector<MSLayersAtAngle>(max(2 * theHalfNBins, 1));
0047   MultipleScatteringGeometry layout(tracker);
0048   for (int idxbin = 0; idxbin < 2 * theHalfNBins; idxbin++) {
0049     float etaValue = eta(idxbin);
0050     float cotTheta = sinh(etaValue);
0051
0052     vector<MSLayer> layers = layout.detLayers(etaValue, 0, bfield);
0053     vector<MSLayer> tmplay = layout.otherLayers(etaValue);
0054     layers.insert(layers.end(), tmplay.begin(), tmplay.end());
0055     sort(layers.begin(), layers.end());
0056     setX0(layers, etaValue, msX0data);
0057     theLayersData[idxbin] = MSLayersAtAngle(layers);
0058     PixelRecoPointRZ zero(0., 0.);
0059     PixelRecoLineRZ line(zero, cotTheta);
0060     vector<MSLayer>::iterator it;
0061     for (it = layers.begin(); it != layers.end(); it++) {
0062       float x0 = getDataX0(*it).x0;
0063       float sumX0D = theLayersData[idxbin].sumX0D(zero, it->crossing(line).first);
0064       setDataX0(*it, DataX0(x0, sumX0D, cotTheta));
0065       theLayersData[idxbin].update(*it);
0066     }
0067   }
0068
0069   // add layers not seen from nominal vertex but crossed if
0070   // vertex seperated from nominal by less than 3 sigma
0071   for (int idxbin = 0; idxbin < 2 * theHalfNBins; idxbin++) {
0072     float etaValue = eta(idxbin);
0073     for (int isign = -1; isign <= 1; isign += 2) {
0074       float z = isign * 15.9;  //3 sigma from zero
0075       const MSLayersAtAngle& layersAtAngle = theLayersData[idxbin];
0076       vector<MSLayer> candidates = layout.detLayers(etaValue, z, bfield);
0077       vector<MSLayer>::iterator it;
0078       for (it = candidates.begin(); it != candidates.end(); it++) {
0079         if (layersAtAngle.findLayer(*it))
0080           continue;
0081         const MSLayer* found = nullptr;
0082         int bin = idxbin;
0083         while (!found) {
0084           bin--;
0085           if (bin < 0)
0086             break;
0087           found = theLayersData[bin].findLayer(*it);
0088         }
0089         bin = idxbin;
0090         while (!found) {
0091           bin++;
0092           if (bin > 2 * theHalfNBins - 1)
0093             break;
0094           found = theLayersData[bin].findLayer(*it);
0095         }
0096         if (found)
0097           theLayersData[idxbin].update(*found);
0098       }
0099     }
0100   }
0101
0102   // cout << "LAYERS, size=: "<<theLayersData.size()<< endl;
0103   /*
0104   for (int idxbin = 0; idxbin <= theHalfNBins; idxbin+=25) {
0105     float etaValue = eta(idxbin);
0106     const MSLayersAtAngle & layers= theLayersData[idxbin];
0107     cout << "ETA: "<< etaValue <<" (bin:"<<idxbin<<") #layers:"
0108          <<layers.size()<<endl;
0109     layers.print();
0110   }
0111   for (int idxbin = 2*theHalfNBins-1; idxbin > theHalfNBins; idxbin-=25) {
0112     float etaValue = eta(idxbin);
0113     const MSLayersAtAngle & layers= theLayersData[idxbin];
0114     cout << "ETA: "<< etaValue <<" (bin:"<<idxbin<<") #layers:"
0115          <<layers.size()<<endl;
0116     layers.print();
0117   }
0118 */
0119 }
0120
0121 //------------------------------------------------------------------------------
0122 void MSLayersKeeperX0AtEta::setX0(vector<MSLayer>& layers, float eta, const SumX0AtEtaDataProvider& sumX0) {
0123   const float BIG = 99999.;
0124   float cotTheta = sinh(eta);
0125   float sinTheta = 1 / sqrt(1 + sqr(cotTheta));
0126   float cosTheta = cotTheta * sinTheta;
0127
0128   float sumX0atAngleLast = 0.;
0129   vector<MSLayer>::iterator il;
0130   for (il = layers.begin(); il != layers.end(); il++) {
0131     PixelRecoLineRZ line(PixelRecoPointRZ(0., 0.), cotTheta);
0132     float rL = (*il).crossing(line).first.r();
0133     float rN = (il + 1 != layers.end()) ? (il + 1)->crossing(line).first.r() : BIG;
0134     float rBound = (rL + rN) / 2.;
0135     float sumX0atAngle = sumX0.sumX0atEta(eta, rBound);
0136
0137     float dX0 = (il->face() == GeomDetEnumerators::barrel) ? (sumX0atAngle - sumX0atAngleLast) * sinTheta
0138                                                            : (sumX0atAngle - sumX0atAngleLast) * fabs(cosTheta);
0139
0140     setDataX0(*il, DataX0(dX0, 0, cotTheta));
0141     sumX0atAngleLast = sumX0atAngle;
0142   }
0143 }
0144
0145 //------------------------------------------------------------------------------
0146 MSLayersKeeperX0AtEta::~MSLayersKeeperX0AtEta() = default;