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