File indexing completed on 2025-06-12 23:30:08
0001 #include "trkCore.h"
0002
0003
0004 bool goodEvent() {
0005 if (ana.specific_event_index >= 0) {
0006 if ((int)ana.looper.getCurrentEventIndex() != ana.specific_event_index)
0007 return false;
0008 }
0009
0010
0011 if (ana.nsplit_jobs >= 0 and ana.job_index >= 0) {
0012 if (ana.looper.getNEventsProcessed() % ana.nsplit_jobs != (unsigned int)ana.job_index)
0013 return false;
0014 }
0015
0016 if (ana.verbose >= 2)
0017 std::cout << " ana.looper.getCurrentEventIndex(): " << ana.looper.getCurrentEventIndex() << std::endl;
0018
0019 return true;
0020 }
0021
0022
0023 float runMiniDoublet(LSTEvent* event, int evt) {
0024 TStopwatch my_timer;
0025 if (ana.verbose >= 2)
0026 std::cout << "Reco Mini-Doublet start " << evt << std::endl;
0027 my_timer.Start();
0028 event->createMiniDoublets();
0029 event->wait();
0030 float md_elapsed = my_timer.RealTime();
0031
0032 if (ana.verbose >= 2)
0033 std::cout << evt << " Reco Mini-doublet processing time: " << md_elapsed << " secs" << std::endl;
0034
0035 if (ana.verbose >= 2)
0036 std::cout << "# of Mini-doublets produced: " << event->getNumberOfMiniDoublets() << std::endl;
0037 if (ana.verbose >= 2)
0038 std::cout << "# of Mini-doublets produced barrel layer 1: " << event->getNumberOfMiniDoubletsByLayerBarrel(0)
0039 << std::endl;
0040 if (ana.verbose >= 2)
0041 std::cout << "# of Mini-doublets produced barrel layer 2: " << event->getNumberOfMiniDoubletsByLayerBarrel(1)
0042 << std::endl;
0043 if (ana.verbose >= 2)
0044 std::cout << "# of Mini-doublets produced barrel layer 3: " << event->getNumberOfMiniDoubletsByLayerBarrel(2)
0045 << std::endl;
0046 if (ana.verbose >= 2)
0047 std::cout << "# of Mini-doublets produced barrel layer 4: " << event->getNumberOfMiniDoubletsByLayerBarrel(3)
0048 << std::endl;
0049 if (ana.verbose >= 2)
0050 std::cout << "# of Mini-doublets produced barrel layer 5: " << event->getNumberOfMiniDoubletsByLayerBarrel(4)
0051 << std::endl;
0052 if (ana.verbose >= 2)
0053 std::cout << "# of Mini-doublets produced barrel layer 6: " << event->getNumberOfMiniDoubletsByLayerBarrel(5)
0054 << std::endl;
0055
0056 if (ana.verbose >= 2)
0057 std::cout << "# of Mini-doublets produced endcap layer 1: " << event->getNumberOfMiniDoubletsByLayerEndcap(0)
0058 << std::endl;
0059 if (ana.verbose >= 2)
0060 std::cout << "# of Mini-doublets produced endcap layer 2: " << event->getNumberOfMiniDoubletsByLayerEndcap(1)
0061 << std::endl;
0062 if (ana.verbose >= 2)
0063 std::cout << "# of Mini-doublets produced endcap layer 3: " << event->getNumberOfMiniDoubletsByLayerEndcap(2)
0064 << std::endl;
0065 if (ana.verbose >= 2)
0066 std::cout << "# of Mini-doublets produced endcap layer 4: " << event->getNumberOfMiniDoubletsByLayerEndcap(3)
0067 << std::endl;
0068 if (ana.verbose >= 2)
0069 std::cout << "# of Mini-doublets produced endcap layer 5: " << event->getNumberOfMiniDoubletsByLayerEndcap(4)
0070 << std::endl;
0071
0072 return md_elapsed;
0073 }
0074
0075
0076 float runSegment(LSTEvent* event) {
0077 TStopwatch my_timer;
0078 if (ana.verbose >= 2)
0079 std::cout << "Reco Segment start" << std::endl;
0080 my_timer.Start();
0081 event->createSegmentsWithModuleMap();
0082 event->wait();
0083 float sg_elapsed = my_timer.RealTime();
0084 if (ana.verbose >= 2)
0085 std::cout << "Reco Segment processing time: " << sg_elapsed << " secs" << std::endl;
0086
0087 if (ana.verbose >= 2)
0088 std::cout << "# of Segments produced: " << event->getNumberOfSegments() << std::endl;
0089 if (ana.verbose >= 2)
0090 std::cout << "# of Segments produced layer 1-2: " << event->getNumberOfSegmentsByLayerBarrel(0) << std::endl;
0091 if (ana.verbose >= 2)
0092 std::cout << "# of Segments produced layer 2-3: " << event->getNumberOfSegmentsByLayerBarrel(1) << std::endl;
0093 if (ana.verbose >= 2)
0094 std::cout << "# of Segments produced layer 3-4: " << event->getNumberOfSegmentsByLayerBarrel(2) << std::endl;
0095 if (ana.verbose >= 2)
0096 std::cout << "# of Segments produced layer 4-5: " << event->getNumberOfSegmentsByLayerBarrel(3) << std::endl;
0097 if (ana.verbose >= 2)
0098 std::cout << "# of Segments produced layer 5-6: " << event->getNumberOfSegmentsByLayerBarrel(4) << std::endl;
0099 if (ana.verbose >= 2)
0100 std::cout << "# of Segments produced endcap layer 1: " << event->getNumberOfSegmentsByLayerEndcap(0) << std::endl;
0101 if (ana.verbose >= 2)
0102 std::cout << "# of Segments produced endcap layer 2: " << event->getNumberOfSegmentsByLayerEndcap(1) << std::endl;
0103 if (ana.verbose >= 2)
0104 std::cout << "# of Segments produced endcap layer 3: " << event->getNumberOfSegmentsByLayerEndcap(2) << std::endl;
0105 if (ana.verbose >= 2)
0106 std::cout << "# of Segments produced endcap layer 4: " << event->getNumberOfSegmentsByLayerEndcap(3) << std::endl;
0107 if (ana.verbose >= 2)
0108 std::cout << "# of Segments produced endcap layer 5: " << event->getNumberOfSegmentsByLayerEndcap(4) << std::endl;
0109
0110 return sg_elapsed;
0111 }
0112
0113
0114 float runT3(LSTEvent* event) {
0115 TStopwatch my_timer;
0116 if (ana.verbose >= 2)
0117 std::cout << "Reco T3 start" << std::endl;
0118 my_timer.Start();
0119 event->createTriplets();
0120 event->wait();
0121 float t3_elapsed = my_timer.RealTime();
0122 if (ana.verbose >= 2)
0123 std::cout << "Reco T3 processing time: " << t3_elapsed << " secs" << std::endl;
0124
0125 if (ana.verbose >= 2)
0126 std::cout << "# of T3s produced: " << event->getNumberOfTriplets() << std::endl;
0127 if (ana.verbose >= 2)
0128 std::cout << "# of T3s produced layer 1-2-3: " << event->getNumberOfTripletsByLayerBarrel(0) << std::endl;
0129 if (ana.verbose >= 2)
0130 std::cout << "# of T3s produced layer 2-3-4: " << event->getNumberOfTripletsByLayerBarrel(1) << std::endl;
0131 if (ana.verbose >= 2)
0132 std::cout << "# of T3s produced layer 3-4-5: " << event->getNumberOfTripletsByLayerBarrel(2) << std::endl;
0133 if (ana.verbose >= 2)
0134 std::cout << "# of T3s produced layer 4-5-6: " << event->getNumberOfTripletsByLayerBarrel(3) << std::endl;
0135 if (ana.verbose >= 2)
0136 std::cout << "# of T3s produced endcap layer 1-2-3: " << event->getNumberOfTripletsByLayerEndcap(0) << std::endl;
0137 if (ana.verbose >= 2)
0138 std::cout << "# of T3s produced endcap layer 2-3-4: " << event->getNumberOfTripletsByLayerEndcap(1) << std::endl;
0139 if (ana.verbose >= 2)
0140 std::cout << "# of T3s produced endcap layer 3-4-5: " << event->getNumberOfTripletsByLayerEndcap(2) << std::endl;
0141 if (ana.verbose >= 2)
0142 std::cout << "# of T3s produced endcap layer 1: " << event->getNumberOfTripletsByLayerEndcap(0) << std::endl;
0143 if (ana.verbose >= 2)
0144 std::cout << "# of T3s produced endcap layer 2: " << event->getNumberOfTripletsByLayerEndcap(1) << std::endl;
0145 if (ana.verbose >= 2)
0146 std::cout << "# of T3s produced endcap layer 3: " << event->getNumberOfTripletsByLayerEndcap(2) << std::endl;
0147 if (ana.verbose >= 2)
0148 std::cout << "# of T3s produced endcap layer 4: " << event->getNumberOfTripletsByLayerEndcap(3) << std::endl;
0149 if (ana.verbose >= 2)
0150 std::cout << "# of T3s produced endcap layer 5: " << event->getNumberOfTripletsByLayerEndcap(4) << std::endl;
0151
0152 return t3_elapsed;
0153 }
0154
0155
0156 float runpT3(LSTEvent* event) {
0157 TStopwatch my_timer;
0158 if (ana.verbose >= 2)
0159 std::cout << "Reco Pixel Triplet pT3 start" << std::endl;
0160 my_timer.Start();
0161 event->createPixelTriplets();
0162 event->wait();
0163 float pt3_elapsed = my_timer.RealTime();
0164 if (ana.verbose >= 2)
0165 std::cout << "Reco pT3 processing time: " << pt3_elapsed << " secs" << std::endl;
0166 if (ana.verbose >= 2)
0167 std::cout << "# of Pixel T3s produced: " << event->getNumberOfPixelTriplets() << std::endl;
0168
0169 return pt3_elapsed;
0170 }
0171
0172
0173 float runQuintuplet(LSTEvent* event) {
0174 TStopwatch my_timer;
0175 if (ana.verbose >= 2)
0176 std::cout << "Reco Quintuplet start" << std::endl;
0177 my_timer.Start();
0178 event->createQuintuplets();
0179 event->wait();
0180 float t5_elapsed = my_timer.RealTime();
0181 if (ana.verbose >= 2)
0182 std::cout << "Reco Quintuplet processing time: " << t5_elapsed << " secs" << std::endl;
0183
0184 if (ana.verbose >= 2)
0185 std::cout << "# of Quintuplets produced: " << event->getNumberOfQuintuplets() << std::endl;
0186 if (ana.verbose >= 2)
0187 std::cout << "# of Quintuplets produced layer 1-2-3-4-5-6: " << event->getNumberOfQuintupletsByLayerBarrel(0)
0188 << std::endl;
0189 if (ana.verbose >= 2)
0190 std::cout << "# of Quintuplets produced layer 2: " << event->getNumberOfQuintupletsByLayerBarrel(1) << std::endl;
0191 if (ana.verbose >= 2)
0192 std::cout << "# of Quintuplets produced layer 3: " << event->getNumberOfQuintupletsByLayerBarrel(2) << std::endl;
0193 if (ana.verbose >= 2)
0194 std::cout << "# of Quintuplets produced layer 4: " << event->getNumberOfQuintupletsByLayerBarrel(3) << std::endl;
0195 if (ana.verbose >= 2)
0196 std::cout << "# of Quintuplets produced layer 5: " << event->getNumberOfQuintupletsByLayerBarrel(4) << std::endl;
0197 if (ana.verbose >= 2)
0198 std::cout << "# of Quintuplets produced layer 6: " << event->getNumberOfQuintupletsByLayerBarrel(5) << std::endl;
0199 if (ana.verbose >= 2)
0200 std::cout << "# of Quintuplets produced endcap layer 1: " << event->getNumberOfQuintupletsByLayerEndcap(0)
0201 << std::endl;
0202 if (ana.verbose >= 2)
0203 std::cout << "# of Quintuplets produced endcap layer 2: " << event->getNumberOfQuintupletsByLayerEndcap(1)
0204 << std::endl;
0205 if (ana.verbose >= 2)
0206 std::cout << "# of Quintuplets produced endcap layer 3: " << event->getNumberOfQuintupletsByLayerEndcap(2)
0207 << std::endl;
0208 if (ana.verbose >= 2)
0209 std::cout << "# of Quintuplets produced endcap layer 4: " << event->getNumberOfQuintupletsByLayerEndcap(3)
0210 << std::endl;
0211 if (ana.verbose >= 2)
0212 std::cout << "# of Quintuplets produced endcap layer 5: " << event->getNumberOfQuintupletsByLayerEndcap(4)
0213 << std::endl;
0214
0215 return t5_elapsed;
0216 }
0217
0218
0219 float runPixelLineSegment(LSTEvent* event, bool no_pls_dupclean) {
0220 TStopwatch my_timer;
0221 if (ana.verbose >= 2)
0222 std::cout << "Reco Pixel Line Segment start" << std::endl;
0223 my_timer.Start();
0224 event->addPixelSegmentToEventFinalize();
0225 event->pixelLineSegmentCleaning(no_pls_dupclean);
0226 event->wait();
0227 float pls_elapsed = my_timer.RealTime();
0228 if (ana.verbose >= 2)
0229 std::cout << "Reco Pixel Line Segment processing time: " << pls_elapsed << " secs" << std::endl;
0230
0231 return pls_elapsed;
0232 }
0233
0234
0235 float runPixelQuintuplet(LSTEvent* event) {
0236 TStopwatch my_timer;
0237 if (ana.verbose >= 2)
0238 std::cout << "Reco Pixel Quintuplet start" << std::endl;
0239 my_timer.Start();
0240 event->createPixelQuintuplets();
0241 event->wait();
0242 float pt5_elapsed = my_timer.RealTime();
0243 if (ana.verbose >= 2)
0244 std::cout << "Reco Pixel Quintuplet processing time: " << pt5_elapsed << " secs" << std::endl;
0245 if (ana.verbose >= 2)
0246 std::cout << "# of Pixel Quintuplets produced: " << event->getNumberOfPixelQuintuplets() << std::endl;
0247
0248 return pt5_elapsed;
0249 }
0250
0251
0252 float runTrackCandidate(LSTEvent* event, bool no_pls_dupclean, bool tc_pls_triplets) {
0253 TStopwatch my_timer;
0254 if (ana.verbose >= 2)
0255 std::cout << "Reco TrackCandidate start" << std::endl;
0256 my_timer.Start();
0257 event->createTrackCandidates(no_pls_dupclean, tc_pls_triplets);
0258 event->wait();
0259 float tc_elapsed = my_timer.RealTime();
0260 if (ana.verbose >= 2)
0261 std::cout << "Reco TrackCandidate processing time: " << tc_elapsed << " secs" << std::endl;
0262
0263 if (ana.verbose >= 2)
0264 std::cout << "# of TrackCandidates produced: " << event->getNumberOfTrackCandidates() << std::endl;
0265 if (ana.verbose >= 2)
0266 std::cout << "# of Pixel TrackCandidates produced: " << event->getNumberOfPixelTrackCandidates() << std::endl;
0267 if (ana.verbose >= 2)
0268 std::cout << " # of pT5 TrackCandidates produced: " << event->getNumberOfPT5TrackCandidates() << std::endl;
0269 if (ana.verbose >= 2)
0270 std::cout << " # of pT3 TrackCandidates produced: " << event->getNumberOfPT3TrackCandidates() << std::endl;
0271 if (ana.verbose >= 2)
0272 std::cout << " # of pLS TrackCandidates produced: " << event->getNumberOfPLSTrackCandidates() << std::endl;
0273 if (ana.verbose >= 2)
0274 std::cout << "# of T5 TrackCandidates produced: " << event->getNumberOfT5TrackCandidates() << std::endl;
0275
0276 return tc_elapsed;
0277 }
0278
0279
0280
0281
0282
0283
0284
0285
0286 std::vector<int> matchedSimTrkIdxs(std::vector<unsigned int> hitidxs,
0287 std::vector<unsigned int> hittypes,
0288 std::vector<int> const& trk_simhit_simTrkIdx,
0289 std::vector<std::vector<int>> const& trk_ph2_simHitIdx,
0290 std::vector<std::vector<int>> const& trk_pix_simHitIdx,
0291 bool verbose,
0292 float* pmatched) {
0293 if (hitidxs.size() != hittypes.size()) {
0294 std::cout << "Error: matched_sim_trk_idxs() hitidxs and hittypes have different lengths" << std::endl;
0295 std::cout << "hitidxs.size(): " << hitidxs.size() << std::endl;
0296 std::cout << "hittypes.size(): " << hittypes.size() << std::endl;
0297 }
0298
0299 std::vector<std::pair<unsigned int, unsigned int>> to_check_duplicate;
0300 for (size_t i = 0; i < hitidxs.size(); ++i) {
0301 auto hitidx = hitidxs[i];
0302 auto hittype = hittypes[i];
0303 auto item = std::make_pair(hitidx, hittype);
0304 if (std::find(to_check_duplicate.begin(), to_check_duplicate.end(), item) == to_check_duplicate.end()) {
0305 to_check_duplicate.push_back(item);
0306 }
0307 }
0308
0309 int nhits_input = to_check_duplicate.size();
0310
0311 std::vector<std::vector<int>> simtrk_idxs;
0312 std::vector<int> unique_idxs;
0313
0314 if (verbose) {
0315 std::cout << " '------------------------': "
0316 << "------------------------" << std::endl;
0317 }
0318
0319 for (size_t ihit = 0; ihit < to_check_duplicate.size(); ++ihit) {
0320 auto ihitdata = to_check_duplicate[ihit];
0321 auto&& [hitidx, hittype] = ihitdata;
0322
0323 if (verbose) {
0324 std::cout << " hitidx: " << hitidx << " hittype: " << hittype << std::endl;
0325 }
0326
0327 std::vector<int> simtrk_idxs_per_hit;
0328
0329 const std::vector<std::vector<int>>* simHitIdxs = hittype == 4 ? &trk_ph2_simHitIdx : &trk_pix_simHitIdx;
0330
0331 if (verbose) {
0332 std::cout << " trk_ph2_simHitIdx.size(): " << trk_ph2_simHitIdx.size() << std::endl;
0333 std::cout << " trk_pix_simHitIdx.size(): " << trk_pix_simHitIdx.size() << std::endl;
0334 }
0335
0336 if (static_cast<const unsigned int>((*simHitIdxs).size()) <= hitidx) {
0337 std::cout << "ERROR" << std::endl;
0338 std::cout << " hittype: " << hittype << std::endl;
0339 std::cout << " trk_pix_simHitIdx.size(): " << trk_pix_simHitIdx.size() << std::endl;
0340 std::cout << " trk_ph2_simHitIdx.size(): " << trk_ph2_simHitIdx.size() << std::endl;
0341 std::cout << (*simHitIdxs).size() << " " << hittype << std::endl;
0342 std::cout << hitidx << " " << hittype << std::endl;
0343 }
0344
0345 for (auto& simhit_idx : (*simHitIdxs).at(hitidx)) {
0346 if (static_cast<const int>(trk_simhit_simTrkIdx.size()) <= simhit_idx) {
0347 std::cout << (*simHitIdxs).size() << " " << hittype << std::endl;
0348 std::cout << hitidx << " " << hittype << std::endl;
0349 std::cout << trk_simhit_simTrkIdx.size() << " " << simhit_idx << std::endl;
0350 }
0351 int simtrk_idx = trk_simhit_simTrkIdx[simhit_idx];
0352 if (verbose) {
0353 std::cout << " hitidx: " << hitidx << " simhit_idx: " << simhit_idx << " simtrk_idx: " << simtrk_idx
0354 << std::endl;
0355 }
0356 simtrk_idxs_per_hit.push_back(simtrk_idx);
0357 if (std::find(unique_idxs.begin(), unique_idxs.end(), simtrk_idx) == unique_idxs.end())
0358 unique_idxs.push_back(simtrk_idx);
0359 }
0360
0361 if (simtrk_idxs_per_hit.size() == 0) {
0362 if (verbose) {
0363 std::cout << " hitidx: " << hitidx << " -1: " << -1 << std::endl;
0364 }
0365 simtrk_idxs_per_hit.push_back(-1);
0366 if (std::find(unique_idxs.begin(), unique_idxs.end(), -1) == unique_idxs.end())
0367 unique_idxs.push_back(-1);
0368 }
0369
0370 simtrk_idxs.push_back(simtrk_idxs_per_hit);
0371 }
0372
0373 if (verbose) {
0374 std::cout << " unique_idxs.size(): " << unique_idxs.size() << std::endl;
0375 for (auto& unique_idx : unique_idxs) {
0376 std::cout << " unique_idx: " << unique_idx << std::endl;
0377 }
0378 }
0379
0380
0381 if (verbose) {
0382 std::cout << "va print" << std::endl;
0383 for (auto& vec : simtrk_idxs) {
0384 for (auto& idx : vec) {
0385 std::cout << idx << " ";
0386 }
0387 std::cout << std::endl;
0388 }
0389 std::cout << "va print end" << std::endl;
0390 }
0391
0392
0393 std::function<void(std::vector<std::vector<int>>&, std::vector<int>, size_t, std::vector<std::vector<int>>&)> perm =
0394 [&](std::vector<std::vector<int>>& result,
0395 std::vector<int> intermediate,
0396 size_t n,
0397 std::vector<std::vector<int>>& va) {
0398
0399 if (va.size() > n) {
0400 for (auto x : va[n]) {
0401
0402
0403 std::vector<int> copy_intermediate(intermediate);
0404 copy_intermediate.push_back(x);
0405 perm(result, copy_intermediate, n + 1, va);
0406 }
0407 } else {
0408 result.push_back(intermediate);
0409 }
0410 };
0411
0412 std::vector<std::vector<int>> allperms;
0413 perm(allperms, std::vector<int>(), 0, simtrk_idxs);
0414
0415 if (verbose) {
0416 std::cout << " allperms.size(): " << allperms.size() << std::endl;
0417 for (unsigned iperm = 0; iperm < allperms.size(); ++iperm) {
0418 std::cout << " allperms[iperm].size(): " << allperms[iperm].size() << std::endl;
0419 for (unsigned ielem = 0; ielem < allperms[iperm].size(); ++ielem) {
0420 std::cout << " allperms[iperm][ielem]: " << allperms[iperm][ielem] << std::endl;
0421 }
0422 }
0423 }
0424 int maxHitMatchCount = 0;
0425 std::vector<int> matched_sim_trk_idxs;
0426 float max_percent_matched = 0.0f;
0427 for (auto& trkidx_perm : allperms) {
0428 std::vector<int> counts;
0429 for (auto& unique_idx : unique_idxs) {
0430 int cnt = std::count(trkidx_perm.begin(), trkidx_perm.end(), unique_idx);
0431 counts.push_back(cnt);
0432 }
0433 auto result = std::max_element(counts.begin(), counts.end());
0434 int rawidx = std::distance(counts.begin(), result);
0435 int trkidx = unique_idxs[rawidx];
0436 if (trkidx < 0)
0437 continue;
0438 float percent_matched = static_cast<float>(counts[rawidx]) / nhits_input;
0439 if (percent_matched > 0.75f)
0440 matched_sim_trk_idxs.push_back(trkidx);
0441 maxHitMatchCount = std::max(maxHitMatchCount, *std::max_element(counts.begin(), counts.end()));
0442 max_percent_matched = std::max(max_percent_matched, percent_matched);
0443 }
0444
0445
0446 if (pmatched != nullptr) {
0447 *pmatched = max_percent_matched;
0448 }
0449
0450 std::set<int> s;
0451 unsigned size = matched_sim_trk_idxs.size();
0452 for (unsigned i = 0; i < size; ++i)
0453 s.insert(matched_sim_trk_idxs[i]);
0454 matched_sim_trk_idxs.assign(s.begin(), s.end());
0455 return matched_sim_trk_idxs;
0456 }
0457
0458
0459 int getDenomSimTrkType(int isimtrk,
0460 std::vector<int> const& trk_sim_q,
0461 std::vector<float> const& trk_sim_pt,
0462 std::vector<float> const& trk_sim_eta,
0463 std::vector<int> const& trk_sim_bunchCrossing,
0464 std::vector<int> const& trk_sim_event,
0465 std::vector<int> const& trk_sim_parentVtxIdx,
0466 std::vector<float> const& trk_simvtx_x,
0467 std::vector<float> const& trk_simvtx_y,
0468 std::vector<float> const& trk_simvtx_z) {
0469 if (isimtrk < 0)
0470 return 0;
0471 const int& q = trk_sim_q[isimtrk];
0472 if (q == 0)
0473 return 1;
0474 const float& pt = trk_sim_pt[isimtrk];
0475 const float& eta = trk_sim_eta[isimtrk];
0476 if (pt < 1 or std::abs(eta) > 2.4)
0477 return 2;
0478 const int& bunch = trk_sim_bunchCrossing[isimtrk];
0479 const int& event = trk_sim_event[isimtrk];
0480 const int& vtxIdx = trk_sim_parentVtxIdx[isimtrk];
0481 const float& vtx_x = trk_simvtx_x[isimtrk];
0482 const float& vtx_y = trk_simvtx_y[isimtrk];
0483 const float& vtx_z = trk_simvtx_z[isimtrk];
0484 const float& vtx_perp = sqrt(vtx_x * vtx_x + vtx_y * vtx_y);
0485 if (vtx_perp > 2.5)
0486 return 3;
0487 if (std::abs(vtx_z) > 30)
0488 return 4;
0489 if (bunch != 0)
0490 return 5;
0491 if (event != 0)
0492 return 6;
0493 return 7;
0494 }
0495
0496
0497 int getDenomSimTrkType(std::vector<int> simidxs,
0498 std::vector<int> const& trk_sim_q,
0499 std::vector<float> const& trk_sim_pt,
0500 std::vector<float> const& trk_sim_eta,
0501 std::vector<int> const& trk_sim_bunchCrossing,
0502 std::vector<int> const& trk_sim_event,
0503 std::vector<int> const& trk_sim_parentVtxIdx,
0504 std::vector<float> const& trk_simvtx_x,
0505 std::vector<float> const& trk_simvtx_y,
0506 std::vector<float> const& trk_simvtx_z) {
0507 int type = 0;
0508 for (auto& simidx : simidxs) {
0509 int this_type = getDenomSimTrkType(simidx,
0510 trk_sim_q,
0511 trk_sim_pt,
0512 trk_sim_eta,
0513 trk_sim_bunchCrossing,
0514 trk_sim_event,
0515 trk_sim_parentVtxIdx,
0516 trk_simvtx_x,
0517 trk_simvtx_y,
0518 trk_simvtx_z);
0519 if (this_type > type) {
0520 type = this_type;
0521 }
0522 }
0523 return type;
0524 }
0525
0526
0527
0528
0529
0530
0531
0532
0533 float drfracSimHitConsistentWithHelix(int isimtrk,
0534 int isimhitidx,
0535 std::vector<float> const& trk_simvtx_x,
0536 std::vector<float> const& trk_simvtx_y,
0537 std::vector<float> const& trk_simvtx_z,
0538 std::vector<float> const& trk_sim_pt,
0539 std::vector<float> const& trk_sim_eta,
0540 std::vector<float> const& trk_sim_phi,
0541 std::vector<int> const& trk_sim_q,
0542 std::vector<float> const& trk_simhit_x,
0543 std::vector<float> const& trk_simhit_y,
0544 std::vector<float> const& trk_simhit_z) {
0545
0546 float vx = trk_simvtx_x[0];
0547 float vy = trk_simvtx_y[0];
0548 float vz = trk_simvtx_z[0];
0549 float pt = trk_sim_pt[isimtrk];
0550 float eta = trk_sim_eta[isimtrk];
0551 float phi = trk_sim_phi[isimtrk];
0552 int charge = trk_sim_q[isimtrk];
0553
0554
0555 lst_math::Helix helix(pt, eta, phi, vx, vy, vz, charge);
0556
0557 return drfracSimHitConsistentWithHelix(helix, isimhitidx, trk_simhit_x, trk_simhit_y, trk_simhit_z);
0558 }
0559
0560
0561 float drfracSimHitConsistentWithHelix(lst_math::Helix& helix,
0562 int isimhitidx,
0563 std::vector<float> const& trk_simhit_x,
0564 std::vector<float> const& trk_simhit_y,
0565 std::vector<float> const& trk_simhit_z) {
0566
0567 std::vector<float> point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]};
0568
0569
0570 float t = helix.infer_t(point);
0571
0572
0573 if (not(t <= M_PI))
0574 return 999;
0575
0576
0577 auto [x, y, z, r] = helix.get_helix_point(t);
0578
0579
0580 float drfrac = std::abs(helix.compare_radius(point)) / r;
0581
0582 return drfrac;
0583 }
0584
0585
0586 float distxySimHitConsistentWithHelix(int isimtrk,
0587 int isimhitidx,
0588 std::vector<float> const& trk_simvtx_x,
0589 std::vector<float> const& trk_simvtx_y,
0590 std::vector<float> const& trk_simvtx_z,
0591 std::vector<float> const& trk_sim_pt,
0592 std::vector<float> const& trk_sim_eta,
0593 std::vector<float> const& trk_sim_phi,
0594 std::vector<int> const& trk_sim_q,
0595 std::vector<float> const& trk_simhit_x,
0596 std::vector<float> const& trk_simhit_y,
0597 std::vector<float> const& trk_simhit_z) {
0598
0599 float vx = trk_simvtx_x[0];
0600 float vy = trk_simvtx_y[0];
0601 float vz = trk_simvtx_z[0];
0602 float pt = trk_sim_pt[isimtrk];
0603 float eta = trk_sim_eta[isimtrk];
0604 float phi = trk_sim_phi[isimtrk];
0605 int charge = trk_sim_q[isimtrk];
0606
0607
0608 lst_math::Helix helix(pt, eta, phi, vx, vy, vz, charge);
0609
0610 return distxySimHitConsistentWithHelix(helix, isimhitidx, trk_simhit_x, trk_simhit_y, trk_simhit_z);
0611 }
0612
0613
0614 float distxySimHitConsistentWithHelix(lst_math::Helix& helix,
0615 int isimhitidx,
0616 std::vector<float> const& trk_simhit_x,
0617 std::vector<float> const& trk_simhit_y,
0618 std::vector<float> const& trk_simhit_z) {
0619
0620 std::vector<float> point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]};
0621
0622
0623 float t = helix.infer_t(point);
0624
0625
0626 if (not(t <= M_PI))
0627 return 999;
0628
0629
0630
0631
0632
0633 float distxy = helix.compare_xy(point);
0634
0635 return distxy;
0636 }
0637
0638
0639 TVector3 calculateR3FromPCA(const TVector3& p3, const float dxy, const float dz) {
0640 const float pt = p3.Pt();
0641 const float p = p3.Mag();
0642 const float vz = dz * pt * pt / p / p;
0643
0644 const float vx = -dxy * p3.y() / pt - p3.x() / p * p3.z() / p * dz;
0645 const float vy = dxy * p3.x() / pt - p3.y() / p * p3.z() / p * dz;
0646 return TVector3(vx, vy, vz);
0647 }
0648
0649
0650
0651
0652
0653
0654
0655
0656 float addInputsToEventPreLoad(LSTEvent* event,
0657 lst::LSTInputHostCollection* lstInputHC,
0658 LSTInputDeviceCollection* lstInputDC,
0659 ALPAKA_ACCELERATOR_NAMESPACE::Queue& queue) {
0660 TStopwatch my_timer;
0661
0662 if (ana.verbose >= 2)
0663 std::cout << "Loading Inputs (i.e. outer tracker hits, and pixel line segements) to the Line Segment Tracking.... "
0664 << std::endl;
0665
0666 my_timer.Start();
0667
0668
0669 alpaka::memcpy(queue, lstInputDC->buffer(), lstInputHC->buffer());
0670 alpaka::wait(queue);
0671
0672 event->addInputToEvent(lstInputDC);
0673 event->addHitToEvent();
0674
0675 event->addPixelSegmentToEventStart();
0676 event->wait();
0677 float hit_loading_elapsed = my_timer.RealTime();
0678
0679 if (ana.verbose >= 2)
0680 std::cout << "Loading inputs processing time: " << hit_loading_elapsed << " secs" << std::endl;
0681
0682 return hit_loading_elapsed;
0683 }
0684
0685
0686 void printTimingInformation(std::vector<std::vector<float>>& timing_information, float fullTime, float fullavg) {
0687 if (ana.verbose == 0)
0688 return;
0689
0690 std::cout << std::showpoint;
0691 std::cout << std::fixed;
0692 std::cout << std::setprecision(2);
0693 std::cout << std::right;
0694 std::cout << "Timing summary" << std::endl;
0695 std::cout << std::setw(6) << "Evt";
0696 std::cout << " " << std::setw(6) << "Hits";
0697 std::cout << " " << std::setw(6) << "MD";
0698 std::cout << " " << std::setw(6) << "LS";
0699 std::cout << " " << std::setw(6) << "T3";
0700 std::cout << " " << std::setw(6) << "T5";
0701 std::cout << " " << std::setw(6) << "pLS";
0702 std::cout << " " << std::setw(6) << "pT5";
0703 std::cout << " " << std::setw(6) << "pT3";
0704 std::cout << " " << std::setw(6) << "TC";
0705 std::cout << " " << std::setw(6) << "Reset";
0706 std::cout << " " << std::setw(7) << "Total";
0707 std::cout << " " << std::setw(7) << "Total(short)";
0708 std::cout << std::endl;
0709 std::vector<float> timing_sum_information(timing_information[0].size());
0710 std::vector<float> timing_shortlist;
0711 std::vector<float> timing_list;
0712 for (size_t ievt = 0; ievt < timing_information.size(); ++ievt) {
0713 auto timing = timing_information[ievt];
0714 float timing_total = 0.f;
0715 float timing_total_short = 0.f;
0716 timing_total += timing[0] * 1000;
0717 timing_total += timing[1] * 1000;
0718 timing_total += timing[2] * 1000;
0719 timing_total += timing[3] * 1000;
0720 timing_total += timing[4] * 1000;
0721 timing_total += timing[5] * 1000;
0722 timing_total += timing[6] * 1000;
0723 timing_total += timing[7] * 1000;
0724 timing_total += timing[8] * 1000;
0725 timing_total_short += timing[1] * 1000;
0726 timing_total_short += timing[2] * 1000;
0727 timing_total_short += timing[3] * 1000;
0728 timing_total_short += timing[4] * 1000;
0729 timing_total_short += timing[6] * 1000;
0730 timing_total_short += timing[7] * 1000;
0731 timing_total_short += timing[8] * 1000;
0732 timing_total_short += timing[9] * 1000;
0733 std::cout << std::setw(6) << ievt;
0734 std::cout << " " << std::setw(6) << timing[0] * 1000;
0735 std::cout << " " << std::setw(6) << timing[1] * 1000;
0736 std::cout << " " << std::setw(6) << timing[2] * 1000;
0737 std::cout << " " << std::setw(6) << timing[3] * 1000;
0738 std::cout << " " << std::setw(6) << timing[4] * 1000;
0739 std::cout << " " << std::setw(6) << timing[5] * 1000;
0740 std::cout << " " << std::setw(6) << timing[6] * 1000;
0741 std::cout << " " << std::setw(6) << timing[7] * 1000;
0742 std::cout << " " << std::setw(6) << timing[8] * 1000;
0743 std::cout << " " << std::setw(6) << timing[9] * 1000;
0744 std::cout << " " << std::setw(7) << timing_total;
0745 std::cout << " " << std::setw(7) << timing_total_short;
0746 std::cout << std::endl;
0747 timing_sum_information[0] += timing[0] * 1000;
0748 timing_sum_information[1] += timing[1] * 1000;
0749 timing_sum_information[2] += timing[2] * 1000;
0750 timing_sum_information[3] += timing[3] * 1000;
0751 timing_sum_information[4] += timing[4] * 1000;
0752 timing_sum_information[5] += timing[5] * 1000;
0753 timing_sum_information[6] += timing[6] * 1000;
0754 timing_sum_information[7] += timing[7] * 1000;
0755 timing_sum_information[8] += timing[8] * 1000;
0756 timing_sum_information[9] += timing[9] * 1000;
0757 timing_shortlist.push_back(timing_total_short);
0758 timing_list.push_back(timing_total);
0759 }
0760 timing_sum_information[0] /= timing_information.size();
0761 timing_sum_information[1] /= timing_information.size();
0762 timing_sum_information[2] /= timing_information.size();
0763 timing_sum_information[3] /= timing_information.size();
0764 timing_sum_information[4] /= timing_information.size();
0765 timing_sum_information[5] /= timing_information.size();
0766 timing_sum_information[6] /= timing_information.size();
0767 timing_sum_information[7] /= timing_information.size();
0768 timing_sum_information[8] /= timing_information.size();
0769 timing_sum_information[9] /= timing_information.size();
0770
0771 float timing_total_avg = 0.0;
0772 timing_total_avg += timing_sum_information[0];
0773 timing_total_avg += timing_sum_information[1];
0774 timing_total_avg += timing_sum_information[2];
0775 timing_total_avg += timing_sum_information[3];
0776 timing_total_avg += timing_sum_information[4];
0777 timing_total_avg += timing_sum_information[5];
0778 timing_total_avg += timing_sum_information[6];
0779 timing_total_avg += timing_sum_information[7];
0780 timing_total_avg += timing_sum_information[8];
0781 timing_total_avg += timing_sum_information[9];
0782 float timing_totalshort_avg = 0.0;
0783 timing_totalshort_avg += timing_sum_information[1];
0784 timing_totalshort_avg += timing_sum_information[2];
0785 timing_totalshort_avg += timing_sum_information[3];
0786 timing_totalshort_avg += timing_sum_information[4];
0787 timing_totalshort_avg += timing_sum_information[6];
0788 timing_totalshort_avg += timing_sum_information[7];
0789 timing_totalshort_avg += timing_sum_information[8];
0790 timing_totalshort_avg += timing_sum_information[9];
0791
0792 float standardDeviation = 0.0;
0793 for (auto shorttime : timing_shortlist) {
0794 standardDeviation += pow(shorttime - timing_totalshort_avg, 2);
0795 }
0796 float stdDev = sqrt(standardDeviation / timing_shortlist.size());
0797
0798 std::cout << std::setprecision(1);
0799 std::cout << std::setw(6) << "Evt";
0800 std::cout << " " << std::setw(6) << "Hits";
0801 std::cout << " " << std::setw(6) << "MD";
0802 std::cout << " " << std::setw(6) << "LS";
0803 std::cout << " " << std::setw(6) << "T3";
0804 std::cout << " " << std::setw(6) << "T5";
0805 std::cout << " " << std::setw(6) << "pLS";
0806 std::cout << " " << std::setw(6) << "pT5";
0807 std::cout << " " << std::setw(6) << "pT3";
0808 std::cout << " " << std::setw(6) << "TC";
0809 std::cout << " " << std::setw(6) << "Reset";
0810 std::cout << " " << std::setw(7) << "Total";
0811 std::cout << " " << std::setw(7) << "Total(short)";
0812 std::cout << std::endl;
0813 std::cout << std::setw(6) << "avg";
0814 std::cout << " " << std::setw(6) << timing_sum_information[0];
0815 std::cout << " " << std::setw(6) << timing_sum_information[1];
0816 std::cout << " " << std::setw(6) << timing_sum_information[2];
0817 std::cout << " " << std::setw(6) << timing_sum_information[3];
0818 std::cout << " " << std::setw(6) << timing_sum_information[4];
0819 std::cout << " " << std::setw(6) << timing_sum_information[5];
0820 std::cout << " " << std::setw(6) << timing_sum_information[6];
0821 std::cout << " " << std::setw(6) << timing_sum_information[7];
0822 std::cout << " " << std::setw(6) << timing_sum_information[8];
0823 std::cout << " " << std::setw(6) << timing_sum_information[9];
0824 std::cout << " " << std::setw(7) << timing_total_avg;
0825 std::cout << " " << std::setw(7) << timing_totalshort_avg;
0826 std::cout << "+/- " << std::setw(4) << stdDev;
0827 std::cout << " " << std::setw(7) << fullavg;
0828 std::cout << " " << ana.compilation_target;
0829 std::cout << "[s=" << ana.streams << "]";
0830 std::cout << std::endl;
0831
0832 std::cout << std::left;
0833 }
0834
0835
0836
0837
0838
0839
0840
0841
0842 TString get_absolute_path_after_check_file_exists(const std::string name) {
0843 std::filesystem::path fullpath = std::filesystem::absolute(name.c_str());
0844
0845
0846 if (not std::filesystem::exists(fullpath)) {
0847 std::cout << "ERROR: Could not find the file = " << fullpath << std::endl;
0848 exit(2);
0849 }
0850 return TString(fullpath.string().c_str());
0851 }
0852
0853
0854 void writeMetaData() {
0855
0856 ana.output_tfile->cd();
0857 gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && echo '' && (cd - > /dev/null) ) > %s.gitversion.txt ",
0858 ana.output_tfile->GetName()));
0859 gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && git rev-parse HEAD && (cd - > /dev/null)) >> %s.gitversion.txt",
0860 ana.output_tfile->GetName()));
0861 gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && echo 'git status' && (cd - > /dev/null)) >> %s.gitversion.txt",
0862 ana.output_tfile->GetName()));
0863 gSystem->Exec(
0864 TString::Format("(cd $TRACKLOOPERDIR && git --no-pager status && (cd - > /dev/null)) >> %s.gitversion.txt",
0865 ana.output_tfile->GetName()));
0866 gSystem->Exec(TString::Format(
0867 "(cd $TRACKLOOPERDIR && echo 'git --no-pager log -n 100' && (cd - > /dev/null)) >> %s.gitversion.txt",
0868 ana.output_tfile->GetName()));
0869 gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && echo 'git diff' && (cd - > /dev/null)) >> %s.gitversion.txt",
0870 ana.output_tfile->GetName()));
0871 gSystem->Exec(
0872 TString::Format("(cd $TRACKLOOPERDIR && git --no-pager diff && (cd - > /dev/null)) >> %s.gitversion.txt",
0873 ana.output_tfile->GetName()));
0874
0875
0876 std::ifstream t(TString::Format("%s.gitversion.txt", ana.output_tfile->GetName()));
0877 std::string str((std::istreambuf_iterator<char>(t)), std::istreambuf_iterator<char>());
0878 TNamed code_tag_data("code_tag_data", str.c_str());
0879 ana.output_tfile->cd();
0880 code_tag_data.Write();
0881 gSystem->Exec(TString::Format("rm %s.gitversion.txt", ana.output_tfile->GetName()));
0882
0883
0884 TString make_log_path = TString::Format("%s/.make.log", ana.track_looper_dir_path.Data());
0885 std::ifstream makelog(make_log_path.Data());
0886 std::string makestr((std::istreambuf_iterator<char>(makelog)), std::istreambuf_iterator<char>());
0887 TNamed make_log("make_log", makestr.c_str());
0888 make_log.Write();
0889
0890
0891 gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && git --no-pager diff && (cd - > /dev/null)) > %s.gitdiff.txt",
0892 ana.output_tfile->GetName()));
0893 std::ifstream gitdiff(TString::Format("%s.gitdiff.txt", ana.output_tfile->GetName()));
0894 std::string strgitdiff((std::istreambuf_iterator<char>(gitdiff)), std::istreambuf_iterator<char>());
0895 TNamed gitdifftnamed("gitdiff", strgitdiff.c_str());
0896 gitdifftnamed.Write();
0897 gSystem->Exec(TString::Format("rm %s.gitdiff.txt", ana.output_tfile->GetName()));
0898
0899
0900 TString maketstr = makestr.c_str();
0901 TString rawstrdata = maketstr.ReplaceAll("MAKETARGET=", "%");
0902 TString targetrawdata = RooUtil::StringUtil::rsplit(rawstrdata, "%")[1];
0903 TString targetdata = RooUtil::StringUtil::split(targetrawdata)[0];
0904 ana.compilation_target = targetdata.Data();
0905
0906
0907 TNamed input("input", ana.input_raw_string.Data());
0908 input.Write();
0909
0910
0911 TNamed full_cmd_line("full_cmd_line", ana.full_cmd_line.Data());
0912 full_cmd_line.Write();
0913
0914
0915 TNamed tracklooper_path("tracklooper_path", ana.track_looper_dir_path.Data());
0916 tracklooper_path.Write();
0917 }