standalone/attic/seedtestMPlex.cc

0001 #include "seedtestMPlex.h"
0002 #include "oneapi/tbb/parallel_for.h"
0003
0004 // #define DEBUG
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;  // average third radius squared
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     // since we have two intersection points, arbitrate which one is closer to layer2 hit
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     // MIMI hack: Config::nlayers_per_seed = 4
0043     // const float seed_z2cut = (Config::nlayers_per_seed * Config::fRadialSpacing) / std::tan(2.0f*std::atan(std::exp(-1.0f*Config::dEtaSeedTrip)));
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     // 0 = first layer, 1 = second layer, 2 = third layer
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;  // pass by reference
0068             // MIMI lay0_hits.selectHitIndices(hit1_z/2.0f,hit1.phi(),Config::seed_z0cut/2.0f,Config::lay01angdiff,cand_hit0_indices,true,false);
0069             // loop over first layer hits
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               // center of negative curved track
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               // negative points of intersection with third layer
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               // center of positive curved track
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               // positive points of intersection with third layer
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               // MIMI lay2_hits.selectHitIndices((2.0f*hit1_z-hit0_z),(lay2_posphi+lay2_negphi)/2.0f,
0105               // MIMI seed_z2cut,(lay2_posphi-lay2_negphi)/2.0f,
0106               // MIMI cand_hit2_indices,true,false);
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           // loop over candidate third layer hits
0114           //temp_thr_seed_idcs.reserve(temp_thr_seed_idcs.size()+cand_hit2_indices.size());
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                 // filter by residual of second layer hit
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                 // now fit a circle, extract pT and d0 from center and radius
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                 // filter by d0 cut 5mm, pT cut 0.5 GeV (radius of 0.5 GeV track)
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               }  // end loop over third layer matches
0145             }    // end loop over first layer matches
0146           }      // end chunk of hits for parallel for
0147           seed_idcs.grow_by(temp_thr_seed_idcs.begin(), temp_thr_seed_idcs.end());
0148         });  // end parallel for loop over second layer hits
0149   }
0150
0151 }  // end namespace mkfit