Skip to content

Commit

Permalink
Change attribute from no-inline to no-fast-math. Back out changes to …
Browse files Browse the repository at this point in the history 
…applyMaterialEffects since it does not vectorize with the conditional on radL.
  • Loading branch information
@gartung
gartung committed May 13, 2022
1 parent c19f9c8 commit ed58228
Showing 1 changed file with 49 additions and 163 deletions.
212 changes: 49 additions & 163 deletions RecoTracker/MkFitCore/src/PropagationMPlex.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -589,7 +589,7 @@ namespace mkfit {

//==============================================================================

void __attribute__((optimize("no-inline"))) propagateHelixToZMPlex(const MPlexLS& inErr,
void __attribute__((optimize("no-fast-math"))) propagateHelixToZMPlex(const MPlexLS& inErr,
const MPlexLV& inPar,
const MPlexQI& inChg,
const MPlexQF& msZ,
Expand Down Expand Up @@ -977,172 +977,58 @@ namespace mkfit {
MPlexLV& outPar,
const int N_proc,
const bool isBarrel) {
float radL[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
radL[n] = hitsRl.constAt(n, 0, 0);
}
float theta[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
theta[n] = outPar.constAt(n, 5, 0);
}
float pt[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
pt[n] = 1.f / outPar.constAt(n, 3, 0); //fixme, make sure it is positive?
}
float p[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
p[n] = pt[n] / std::sin(theta[n]);
}
float p2[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
p2[n] = p[n] * p[n];
}
constexpr float mpi = 0.140; // m=140 MeV, pion
constexpr float mpi2 = mpi * mpi; // m=140 MeV, pion
float beta2[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
beta2[n] = p2[n] / (p2[n] + mpi2);
}
float beta[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
beta[n] = std::sqrt(beta2[n]);
}
//radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum)
float invCos[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
invCos[n] = (isBarrel ? p[n] / pt[n] : 1.f / std::abs(std::cos(theta[n])));
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
radL[n] = radL[n] * invCos[n]; //fixme works only for barrel geom
}
// multiple scattering
//vary independently phi and theta by the rms of the planar multiple scattering angle
// XXX-KMD radL < 0, see your fixme above! Repeating bailout
// if (radL < 1e-13f)
// continue;
// const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15
// const float thetaMSC2 = thetaMSC*thetaMSC;
float thetaMSC[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
thetaMSC[n] = 0.0136f * (1.f + 0.038f * std::log(radL[n])) / (beta[n] * p[n]); // eq 32.15
}
float thetaMSC2[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
thetaMSC2[n] = thetaMSC[n] * thetaMSC[n] * radL[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
outErr.At(n, 4, 4) += thetaMSC2[n];
}
// outErr.At(n, 4, 5) += thetaMSC2;
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
outErr.At(n, 5, 5) += thetaMSC2[n];
}
//std::cout << "beta=" << beta << " p=" << p << std::endl;
//std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl;
// energy loss
// XXX-KMD beta2 = 1 => 1 / sqrt(0)
// const float gamma = 1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
// const float gamma2 = gamma*gamma;
float gamma2[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
gamma2[n] = (p2[n] + mpi2) / mpi2;
}
float gamma[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
gamma[n] = std::sqrt(gamma2[n]); //1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
}
constexpr float me = 0.0005; // m=0.5 MeV, electron
float wmax[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
wmax[n] = 2.f * me * beta2[n] * gamma2[n] / (1.f + 2.f * gamma[n] * me / mpi + me * me / (mpi * mpi));
}
constexpr float I = 16.0e-9 * 10.75;
float deltahalf[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
deltahalf[n] = std::log(28.816e-9f * std::sqrt(2.33f * 0.498f) / I) + std::log(beta[n] * gamma[n]) - 0.5f;
}
float dEdx[NN];
float dedxtmp[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
dedxtmp[n] =
2.f *
(hitsXi.constAt(n, 0, 0) * invCos[n] *
(0.5f * std::log(2.f * me * beta2[n] * gamma2[n] * wmax[n] / (I * I)) - beta2[n] - deltahalf[n]) / beta2[n]);
dEdx[n] = beta[n] < 1.f ? dedxtmp[n] : 0.f; //protect against infs and nans
}
// dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg
//std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.constAt(n,0,0) << std::endl;
float dP[NN];

#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
float radL = hitsRl.constAt(n, 0, 0);
if (radL < 1e-13f)
continue; //ugly, please fixme
dP[n] = propSign.constAt(n, 0, 0) * dEdx[n] / beta[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
outPar.At(n, 3, 0) = p[n] / (std::max(p[n] + dP[n], 0.001f) * pt[n]); //stay above 1MeV
const float theta = outPar.constAt(n, 5, 0);
const float pt = 1.f / outPar.constAt(n, 3, 0); //fixme, make sure it is positive?
const float p = pt / std::sin(theta);
const float p2 = p * p;
constexpr float mpi = 0.140; // m=140 MeV, pion
constexpr float mpi2 = mpi * mpi; // m=140 MeV, pion
const float beta2 = p2 / (p2 + mpi2);
const float beta = std::sqrt(beta2);
//radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum)
const float invCos = (isBarrel ? p / pt : 1.f / std::abs(std::cos(theta)));
radL = radL * invCos; //fixme works only for barrel geom
// multiple scattering
//vary independently phi and theta by the rms of the planar multiple scattering angle
// XXX-KMD radL < 0, see your fixme above! Repeating bailout
if (radL < 1e-13f)
continue;
// const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15
// const float thetaMSC2 = thetaMSC*thetaMSC;
const float thetaMSC = 0.0136f * (1.f + 0.038f * std::log(radL)) / (beta * p); // eq 32.15
const float thetaMSC2 = thetaMSC * thetaMSC * radL;
outErr.At(n, 4, 4) += thetaMSC2;
// outErr.At(n, 4, 5) += thetaMSC2;
outErr.At(n, 5, 5) += thetaMSC2;
//std::cout << "beta=" << beta << " p=" << p << std::endl;
//std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl;
// energy loss
// XXX-KMD beta2 = 1 => 1 / sqrt(0)
// const float gamma = 1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
// const float gamma2 = gamma*gamma;
const float gamma2 = (p2 + mpi2) / mpi2;
const float gamma = std::sqrt(gamma2); //1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
constexpr float me = 0.0005; // m=0.5 MeV, electron
const float wmax = 2.f * me * beta2 * gamma2 / (1.f + 2.f * gamma * me / mpi + me * me / (mpi * mpi));
constexpr float I = 16.0e-9 * 10.75;
const float deltahalf = std::log(28.816e-9f * std::sqrt(2.33f * 0.498f) / I) + std::log(beta * gamma) - 0.5f;
const float dEdx =
beta < 1.f
? (2.f * (hitsXi.constAt(n, 0, 0) * invCos *
(0.5f * std::log(2.f * me * beta2 * gamma2 * wmax / (I * I)) - beta2 - deltahalf) / beta2))
: 0.f; //protect against infs and nans
// dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg
//std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.constAt(n,0,0) << std::endl;
const float dP = propSign.constAt(n, 0, 0) * dEdx / beta;
outPar.At(n, 3, 0) = p / (std::max(p + dP, 0.001f) * pt); //stay above 1MeV
//assume 100% uncertainty
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
outErr.At(n, 3, 3) += dP[n] * dP[n] / (p2[n] * pt[n] * pt[n]);
outErr.At(n, 3, 3) += dP * dP / (p2 * pt * pt);
}
}

Expand Down

0 comments on commit ed58228

Please sign in to comment.