File indexing completed on 2023-03-17 11:22:37
0001 #include "seedtestMPlex.h"
0002 #include "oneapi/tbb/parallel_for.h"
0003
0004
0005 #include "Debug.h"
0006
0007 namespace mkfit {
0008
0009 inline void intersectThirdLayer(
0010 const float a, const float b, const float hit1_x, const float hit1_y, float& lay2_x, float& lay2_y) {
0011 const float a2 = a * a;
0012 const float b2 = b * b;
0013 const float a2b2 = a2 + b2;
0014 const float lay2rad2 = (Config::fRadialSpacing * Config::fRadialSpacing) * 9.0f;
0015 const float maxCurvR2 = Config::maxCurvR * Config::maxCurvR;
0016
0017 const float quad =
0018 std::sqrt(2.0f * maxCurvR2 * (a2b2 + lay2rad2) - (a2b2 - lay2rad2) * (a2b2 - lay2rad2) - maxCurvR2 * maxCurvR2);
0019 const float pos[2] = {(a2 * a + a * (b2 + lay2rad2 - maxCurvR2) - b * quad) / a2b2,
0020 (b2 * b + b * (a2 + lay2rad2 - maxCurvR2) + a * quad) / a2b2};
0021 const float neg[2] = {(a2 * a + a * (b2 + lay2rad2 - maxCurvR2) + b * quad) / a2b2,
0022 (b2 * b + b * (a2 + lay2rad2 - maxCurvR2) - a * quad) / a2b2};
0023
0024
0025 if (getHypot(pos[0] - hit1_x, pos[1] - hit1_y) < getHypot(neg[0] - hit1_x, neg[1] - hit1_y)) {
0026 lay2_x = pos[0];
0027 lay2_y = pos[1];
0028 } else {
0029 lay2_x = neg[0];
0030 lay2_y = neg[1];
0031 }
0032 }
0033
0034 void findSeedsByRoadSearch(TripletIdxConVec& seed_idcs,
0035 std::vector<LayerOfHits>& evt_lay_hits,
0036 int lay1_size,
0037 Event*& ev) {
0038 #ifdef DEBUG
0039 bool debug(false);
0040 #endif
0041
0042
0043
0044 #ifdef DEBUG
0045 const float seed_z2cut =
0046 (4 * Config::fRadialSpacing) / std::tan(2.0f * std::atan(std::exp(-1.0f * Config::dEtaSeedTrip)));
0047 #endif
0048
0049
0050 const LayerOfHits& lay1_hits = evt_lay_hits[1];
0051 LayerOfHits& lay0_hits = evt_lay_hits[0];
0052 LayerOfHits& lay2_hits = evt_lay_hits[2];
0053
0054 tbb::parallel_for(
0055 tbb::blocked_range<int>(0, lay1_size, std::max(1, Config::numHitsPerTask)),
0056 [&](const tbb::blocked_range<int>& i) {
0057 TripletIdxVec temp_thr_seed_idcs;
0058 for (int ihit1 = i.begin(); ihit1 < i.end(); ++ihit1) {
0059 const Hit& hit1 = lay1_hits.refHit(ihit1);
0060 const float hit1_z = hit1.z();
0061
0062 dprint("ihit1: " << ihit1 << " mcTrackID: " << hit1.mcTrackID(ev->simHitsInfo_) << " phi: " << hit1.phi()
0063 << " z: " << hit1.z());
0064 dprint(" predphi: " << hit1.phi() << "+/-" << Config::lay01angdiff << " predz: " << hit1.z() / 2.0f << "+/-"
0065 << Config::seed_z0cut / 2.0f << std::endl);
0066
0067 std::vector<int> cand_hit0_indices;
0068
0069
0070 for (auto&& ihit0 : cand_hit0_indices) {
0071 const Hit& hit0 = lay0_hits.refHit(ihit0);
0072 const float hit0_z = hit0.z();
0073 const float hit0_x = hit0.x();
0074 const float hit0_y = hit0.y();
0075 const float hit1_x = hit1.x();
0076 const float hit1_y = hit1.y();
0077 const float hit01_r2 = getRad2(hit0_x - hit1_x, hit0_y - hit1_y);
0078
0079 const float quad = std::sqrt((4.0f * Config::maxCurvR * Config::maxCurvR - hit01_r2) / hit01_r2);
0080
0081
0082 const float aneg = 0.5f * ((hit0_x + hit1_x) - (hit0_y - hit1_y) * quad);
0083 const float bneg = 0.5f * ((hit0_y + hit1_y) + (hit0_x - hit1_x) * quad);
0084
0085
0086 float lay2_negx = 0.0f, lay2_negy = 0.0f;
0087 intersectThirdLayer(aneg, bneg, hit1_x, hit1_y, lay2_negx, lay2_negy);
0088 #ifdef DEBUG
0089 const float lay2_negphi = getPhi(lay2_negx, lay2_negy);
0090 #endif
0091
0092
0093 const float apos = 0.5f * ((hit0_x + hit1_x) + (hit0_y - hit1_y) * quad);
0094 const float bpos = 0.5f * ((hit0_y + hit1_y) - (hit0_x - hit1_x) * quad);
0095
0096
0097 float lay2_posx = 0.0f, lay2_posy = 0.0f;
0098 intersectThirdLayer(apos, bpos, hit1_x, hit1_y, lay2_posx, lay2_posy);
0099 #ifdef DEBUG
0100 const float lay2_posphi = getPhi(lay2_posx, lay2_posy);
0101 #endif
0102
0103 std::vector<int> cand_hit2_indices;
0104
0105
0106
0107
0108 dprint(" ihit0: " << ihit0 << " mcTrackID: " << hit0.mcTrackID(ev->simHitsInfo_) << " phi: " << hit0.phi()
0109 << " z: " << hit0.z());
0110 dprint(" predphi: " << (lay2_posphi + lay2_negphi) / 2.0f << "+/-" << (lay2_posphi - lay2_negphi) / 2.0f
0111 << " predz: " << 2.0f * hit1_z - hit0_z << "+/-" << seed_z2cut << std::endl);
0112
0113
0114
0115 #pragma omp simd
0116 for (size_t idx = 0; idx < cand_hit2_indices.size(); ++idx) {
0117 const int ihit2 = cand_hit2_indices[idx];
0118 const Hit& hit2 = lay2_hits.refHit(ihit2);
0119
0120 const float lay1_predz = (hit0_z + hit2.z()) / 2.0f;
0121
0122 if (std::abs(lay1_predz - hit1_z) > Config::seed_z1cut)
0123 continue;
0124
0125 const float hit2_x = hit2.x();
0126 const float hit2_y = hit2.y();
0127
0128
0129 const float mr = (hit1_y - hit0_y) / (hit1_x - hit0_x);
0130 const float mt = (hit2_y - hit1_y) / (hit2_x - hit1_x);
0131 const float a = (mr * mt * (hit2_y - hit0_y) + mr * (hit1_x + hit2_x) - mt * (hit0_x + hit1_x)) /
0132 (2.0f * (mr - mt));
0133 const float b = -1.0f * (a - (hit0_x + hit1_x) / 2.0f) / mr + (hit0_y + hit1_y) / 2.0f;
0134 const float r = getHypot(hit0_x - a, hit0_y - b);
0135
0136
0137 if ((r < Config::maxCurvR) || (std::abs(getHypot(a, b) - r) > Config::seed_d0cut))
0138 continue;
0139
0140 dprint(" ihit2: " << ihit2 << " mcTrackID: " << hit2.mcTrackID(ev->simHitsInfo_)
0141 << " phi: " << hit2.phi() << " z: " << hit2.z());
0142
0143 temp_thr_seed_idcs.emplace_back(TripletIdx{{ihit0, ihit1, ihit2}});
0144 }
0145 }
0146 }
0147 seed_idcs.grow_by(temp_thr_seed_idcs.begin(), temp_thr_seed_idcs.end());
0148 });
0149 }
0150
0151 }