MkFitCore/src/PropagationMPlexCommon.cc

0001 #include "FWCore/Utilities/interface/CMSUnrollLoop.h"
0002 #include "RecoTracker/MkFitCore/interface/PropagationConfig.h"
0003 #include "RecoTracker/MkFitCore/interface/Config.h"
0004 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0005
0006 #include "PropagationMPlex.h"
0007
0008 //#define DEBUG
0009 #include "Debug.h"
0010
0011 namespace mkfit {
0012
0013   // this version does not assume to know which elements are 0 or 1, so it does the full multiplication
0014   void MultHelixPropFull(const MPlexLL& A, const MPlexLS& B, MPlexLL& C) {
0015     for (int n = 0; n < NN; ++n) {
0016       for (int i = 0; i < 6; ++i) {
0017 // optimization reports indicate only the inner two loops are good
0018 // candidates for vectorization
0019 #pragma omp simd
0020         for (int j = 0; j < 6; ++j) {
0021           C(n, i, j) = 0.;
0022           for (int k = 0; k < 6; ++k)
0023             C(n, i, j) += A.constAt(n, i, k) * B.constAt(n, k, j);
0024         }
0025       }
0026     }
0027   }
0028
0029   // this version does not assume to know which elements are 0 or 1, so it does the full mupltiplication
0030   void MultHelixPropTranspFull(const MPlexLL& A, const MPlexLL& B, MPlexLS& C) {
0031     for (int n = 0; n < NN; ++n) {
0032       for (int i = 0; i < 6; ++i) {
0033 // optimization reports indicate only the inner two loops are good
0034 // candidates for vectorization
0035 #pragma omp simd
0036         for (int j = 0; j < 6; ++j) {
0037           C(n, i, j) = 0.;
0038           for (int k = 0; k < 6; ++k)
0039             C(n, i, j) += B.constAt(n, i, k) * A.constAt(n, j, k);
0040         }
0041       }
0042     }
0043   }
0044
0045 #ifdef UNUSED
0046   // this version does not assume to know which elements are 0 or 1, so it does the full multiplication
0047   void MultHelixPropFull(const MPlexLL& A, const MPlexLL& B, MPlexLL& C) {
0048 #pragma omp simd
0049     for (int n = 0; n < NN; ++n) {
0050       for (int i = 0; i < 6; ++i) {
0051         for (int j = 0; j < 6; ++j) {
0052           C(n, i, j) = 0.;
0053           for (int k = 0; k < 6; ++k)
0054             C(n, i, j) += A.constAt(n, i, k) * B.constAt(n, k, j);
0055         }
0056       }
0057     }
0058   }
0059
0060   // this version does not assume to know which elements are 0 or 1, so it does the full mupltiplication
0061   void MultHelixPropTranspFull(const MPlexLL& A, const MPlexLL& B, MPlexLL& C) {
0062 #pragma omp simd
0063     for (int n = 0; n < NN; ++n) {
0064       for (int i = 0; i < 6; ++i) {
0065         for (int j = 0; j < 6; ++j) {
0066           C(n, i, j) = 0.;
0067           for (int k = 0; k < 6; ++k)
0068             C(n, i, j) += B.constAt(n, i, k) * A.constAt(n, j, k);
0069         }
0070       }
0071     }
0072   }
0073 #endif
0074
0075   //==============================================================================
0076
0077   void applyMaterialEffects(const MPlexQF& hitsRl,
0078                             const MPlexQF& hitsXi,
0079                             const MPlexQF& propSign,
0080                             const MPlexHV& plNrm,
0081                             MPlexLS& outErr,
0082                             MPlexLV& outPar,
0083                             const int N_proc) {
0084 #pragma omp simd
0085     for (int n = 0; n < NN; ++n) {
0086       if (n >= N_proc)
0087         continue;
0088       float radL = hitsRl.constAt(n, 0, 0);
0089       if (radL < 1e-13f)
0090         continue;  //ugly, please fixme
0091       const float theta = outPar.constAt(n, 5, 0);
0092       // const float pt = 1.f / outPar.constAt(n, 3, 0);  //fixme, make sure it is positive?
0093       const float ipt = outPar.constAt(n, 3, 0);
0094       const float pt = 1.f / ipt;  //fixme, make sure it is positive?
0095       const float ipt2 = ipt * ipt;
0096       const float p = pt / std::sin(theta);
0097       const float pz = p * std::cos(theta);
0098       const float p2 = p * p;
0099       constexpr float mpi = 0.140;       // m=140 MeV, pion
0100       constexpr float mpi2 = mpi * mpi;  // m=140 MeV, pion
0101       const float beta2 = p2 / (p2 + mpi2);
0102       const float beta = std::sqrt(beta2);
0103       //radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum)
0104       const float invCos =
0105           p / std::abs(pt * std::cos(outPar.constAt(n, 4, 0)) * plNrm.constAt(n, 0, 0) +
0106                        pt * std::sin(outPar.constAt(n, 4, 0)) * plNrm.constAt(n, 1, 0) + pz * plNrm.constAt(n, 2, 0));
0107       radL = radL * invCos;  //fixme works only for barrel geom
0108       // multiple scattering
0109       //vary independently phi and theta by the rms of the planar multiple scattering angle
0110       // XXX-KMD radL < 0, see your fixme above! Repeating bailout
0111       if (radL < 1e-13f)
0112         continue;
0113       // const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15
0114       // const float thetaMSC2 = thetaMSC*thetaMSC;
0115       const float thetaMSC = 0.0136f * (1.f + 0.038f * std::log(radL)) / (beta * p);  // eq 32.15
0116       const float thetaMSC2 = thetaMSC * thetaMSC * radL;
0117       if /*constexpr*/ (Config::usePtMultScat) {
0118         outErr.At(n, 3, 3) += thetaMSC2 * pz * pz * ipt2 * ipt2;
0119         outErr.At(n, 3, 5) -= thetaMSC2 * pz * ipt2;
0120         outErr.At(n, 4, 4) += thetaMSC2 * p2 * ipt2;
0121         outErr.At(n, 5, 5) += thetaMSC2;
0122       } else {
0123         outErr.At(n, 4, 4) += thetaMSC2;
0124         outErr.At(n, 5, 5) += thetaMSC2;
0125       }
0126       //std::cout << "beta=" << beta << " p=" << p << std::endl;
0127       //std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl;
0128       //std::cout << "radL=" << hitsRl.constAt(n, 0, 0) << " beta=" << beta << " invCos=" << invCos << " radLCorr=" << radL << " thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << std::endl;
0129       // energy loss
0130       // XXX-KMD beta2 = 1 => 1 / sqrt(0)
0131       // const float gamma = 1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
0132       // const float gamma2 = gamma*gamma;
0133       const float gamma2 = (p2 + mpi2) / mpi2;
0134       const float gamma = std::sqrt(gamma2);  //1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
0135       constexpr float me = 0.0005;            // m=0.5 MeV, electron
0136       const float wmax = 2.f * me * beta2 * gamma2 / (1.f + 2.f * gamma * me / mpi + me * me / (mpi * mpi));
0137       constexpr float I = 16.0e-9 * 10.75;
0138       const float deltahalf = std::log(28.816e-9f * std::sqrt(2.33f * 0.498f) / I) + std::log(beta * gamma) - 0.5f;
0139       const float dEdx =
0140           beta < 1.f
0141               ? (2.f * (hitsXi.constAt(n, 0, 0) * invCos *
0142                         (0.5f * std::log(2.f * me * beta2 * gamma2 * wmax / (I * I)) - beta2 - deltahalf) / beta2))
0143               : 0.f;  //protect against infs and nans
0144       // dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg
0145       //std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.constAt(n,0,0) << std::endl;
0146       const float dP = propSign.constAt(n, 0, 0) * dEdx / beta;
0147       outPar.At(n, 3, 0) = p / (std::max(p - dP, 0.001f) * pt);  //stay above 1MeV
0148       //assume 100% uncertainty
0149       outErr.At(n, 3, 3) += dP * dP / (p2 * pt * pt);
0150     }
0151   }
0152
0153 }  // namespace mkfit