Skip to content

Commit

Permalink
code format
Browse files Browse the repository at this point in the history
  • Loading branch information
@gartung
gartung committed May 13, 2022
1 parent 373ea0d commit c19f9c8
Show file tree
Hide file tree
Showing 2 changed files with 308 additions and 255 deletions.
157 changes: 89 additions & 68 deletions RecoTracker/MkFitCore/src/PropagationMPlex.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -588,16 +588,16 @@ namespace mkfit {
}

//==============================================================================
void __attribute__ ((optimize("no-inline"))) propagateHelixToZMPlex(const MPlexLS& inErr,
const MPlexLV& inPar,
const MPlexQI& inChg,
const MPlexQF& msZ,
MPlexLS& outErr,
MPlexLV& outPar,
const int N_proc,
const PropagationFlags pflags,
const MPlexQI* noMatEffPtr) {

void __attribute__((optimize("no-inline"))) propagateHelixToZMPlex(const MPlexLS& inErr,
const MPlexLV& inPar,
const MPlexQI& inChg,
const MPlexQF& msZ,
MPlexLS& outErr,
MPlexLV& outPar,
const int N_proc,
const PropagationFlags pflags,
const MPlexQI* noMatEffPtr) {
// debug = true;

outErr = inErr;
Expand Down Expand Up @@ -708,7 +708,7 @@ namespace mkfit {
errorProp(n, 3, 3) = 1.f;
errorProp(n, 4, 4) = 1.f;
errorProp(n, 5, 5) = 1.f;
}
}
float zout[NN];
float zin[NN];
float ipt[NN];
Expand All @@ -728,14 +728,13 @@ namespace mkfit {
if (pflags.use_param_b_field) {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol *
Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
k[n] = inChg.constAt(n, 0, 0) * 100.f /
(-Const::sol * Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
}
} else {
} else {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol *
Config::Bfield);
k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
}
}

Expand Down Expand Up @@ -763,8 +762,8 @@ namespace mkfit {
for (int n = 0; n < NN; ++n) {
pt[n] = 1.f / ipt[n];
}
//no trig approx here, phi can be large

//no trig approx here, phi can be large
float cosP[NN];
float sinP[NN];
#pragma omp simd
Expand Down Expand Up @@ -802,7 +801,6 @@ namespace mkfit {
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {

//fixme, make this printout useful for propagation to z
dprint_np(n,
std::endl
Expand Down Expand Up @@ -847,7 +845,7 @@ namespace mkfit {
cosa[n] = 1.f - 2.f * sinah[n] * sinah[n];
sina[n] = 2.f * sinah[n] * cosah[n];
}
//update parameters
//update parameters
#pragma omp simd
for (int n = 0; n < NN; ++n) {
outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f * k[n] * sinah[n] * (pxin[n] * cosah[n] - pyin[n] * sinah[n]);
Expand All @@ -873,7 +871,9 @@ namespace mkfit {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 0, 2) = -tanT[n] * ipt[n] * pxcaMpysa[n];
errorProp(n, 0, 3) = k[n] * pt[n] * pt[n] * (cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
errorProp(n, 0, 3) =
k[n] * pt[n] * pt[n] *
(cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
errorProp(n, 0, 4) = -2.f * k[n] * pt[n] * sinah[n] * (sinP[n] * cosah[n] + cosP[n] * sinah[n]);
errorProp(n, 0, 5) = deltaZ[n] * ipt[n] * pxcaMpysa[n] * icos2T[n];
}
Expand All @@ -887,7 +887,9 @@ namespace mkfit {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 1, 2) = -tanT[n] * ipt[n] * pycaPpxsa[n];
errorProp(n, 1, 3) = k[n] * pt[n] * pt[n] * (sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
errorProp(n, 1, 3) =
k[n] * pt[n] * pt[n] *
(sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
errorProp(n, 1, 4) = 2.f * k[n] * pt[n] * sinah[n] * (cosP[n] * cosah[n] - sinP[n] * sinah[n]);
errorProp(n, 1, 5) = deltaZ[n] * ipt[n] * pycaPpxsa[n] * icos2T[n];
}
Expand Down Expand Up @@ -975,7 +977,6 @@ namespace mkfit {
MPlexLV& outPar,
const int N_proc,
const bool isBarrel) {

float radL[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
Expand All @@ -984,143 +985,163 @@ namespace mkfit {
float theta[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f) continue; //ugly, please fixme
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
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
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
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
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
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)
//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
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
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;
// 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
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
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
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
outErr.At(n, 4, 4) += thetaMSC2[n];
}
// outErr.At(n, 4, 5) += thetaMSC2;
// 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
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;
//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
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));
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
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
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
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
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;
// 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) continue; //ugly, please fixme
dP[n]= propSign.constAt(n, 0, 0) * dEdx[n] / beta[n];
if (radL[n] < 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
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
//assume 100% uncertainty
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (radL[n] < 1e-13f) continue; //ugly, please fixme
if (radL[n] < 1e-13f)
continue; //ugly, please fixme
outErr.At(n, 3, 3) += dP[n] * dP[n] / (p2[n] * pt[n] * pt[n]);
}
}
Expand Down
Loading

0 comments on commit c19f9c8

Please sign in to comment.