Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-12-05 02:48:10 

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   // If splitting jobs are requested then determine whether to process the event or not based on remainder
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();  // device side event calls are asynchronous: wait to measure time or print
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();  // device side event calls are asynchronous: wait to measure time or print
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();  // device side event calls are asynchronous: wait to measure time or print
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();  // device side event calls are asynchronous: wait to measure time or print
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();  // device side event calls are asynchronous: wait to measure time or print
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->pixelLineSegmentCleaning(no_pls_dupclean);
0225   event->wait();  // device side event calls are asynchronous: wait to measure time or print
0226   float pls_elapsed = my_timer.RealTime();
0227   if (ana.verbose >= 2)
0228     std::cout << "Reco Pixel Line Segment processing time: " << pls_elapsed << " secs" << std::endl;
0229 
0230   return pls_elapsed;
0231 }
0232 
0233 //___________________________________________________________________________________________________________________________________________________________________________________________
0234 float runPixelQuintuplet(LSTEvent *event) {
0235   TStopwatch my_timer;
0236   if (ana.verbose >= 2)
0237     std::cout << "Reco Pixel Quintuplet start" << std::endl;
0238   my_timer.Start();
0239   event->createPixelQuintuplets();
0240   event->wait();  // device side event calls are asynchronous: wait to measure time or print
0241   float pt5_elapsed = my_timer.RealTime();
0242   if (ana.verbose >= 2)
0243     std::cout << "Reco Pixel Quintuplet processing time: " << pt5_elapsed << " secs" << std::endl;
0244   if (ana.verbose >= 2)
0245     std::cout << "# of Pixel Quintuplets produced: " << event->getNumberOfPixelQuintuplets() << std::endl;
0246 
0247   return pt5_elapsed;
0248 }
0249 
0250 //___________________________________________________________________________________________________________________________________________________________________________________________
0251 float runTrackCandidate(LSTEvent *event, bool no_pls_dupclean, bool tc_pls_triplets) {
0252   TStopwatch my_timer;
0253   if (ana.verbose >= 2)
0254     std::cout << "Reco TrackCandidate start" << std::endl;
0255   my_timer.Start();
0256   event->createTrackCandidates(no_pls_dupclean, tc_pls_triplets);
0257   event->wait();  // device side event calls are asynchronous: wait to measure time or print
0258   float tc_elapsed = my_timer.RealTime();
0259   if (ana.verbose >= 2)
0260     std::cout << "Reco TrackCandidate processing time: " << tc_elapsed << " secs" << std::endl;
0261 
0262   if (ana.verbose >= 2)
0263     std::cout << "# of TrackCandidates produced: " << event->getNumberOfTrackCandidates() << std::endl;
0264   if (ana.verbose >= 2)
0265     std::cout << "# of Pixel TrackCandidates produced: " << event->getNumberOfPixelTrackCandidates() << std::endl;
0266   if (ana.verbose >= 2)
0267     std::cout << "    # of pT5 TrackCandidates produced: " << event->getNumberOfPT5TrackCandidates() << std::endl;
0268   if (ana.verbose >= 2)
0269     std::cout << "    # of pT3 TrackCandidates produced: " << event->getNumberOfPT3TrackCandidates() << std::endl;
0270   if (ana.verbose >= 2)
0271     std::cout << "    # of pLS TrackCandidates produced: " << event->getNumberOfPLSTrackCandidates() << std::endl;
0272   if (ana.verbose >= 2)
0273     std::cout << "# of T5 TrackCandidates produced: " << event->getNumberOfT5TrackCandidates() << std::endl;
0274 
0275   return tc_elapsed;
0276 }
0277 
0278 //  ---------------------------------- =========================================== ----------------------------------------------
0279 //  ---------------------------------- =========================================== ----------------------------------------------
0280 //  ---------------------------------- =========================================== ----------------------------------------------
0281 //  ---------------------------------- =========================================== ----------------------------------------------
0282 //  ---------------------------------- =========================================== ----------------------------------------------
0283 
0284 //___________________________________________________________________________________________________________________________________________________________________________________________
0285 std::vector<int> matchedSimTrkIdxs(std::vector<int> hitidxs, std::vector<int> hittypes, bool verbose) {
0286   std::vector<unsigned int> hitidxs_(std::begin(hitidxs), std::end(hitidxs));
0287   std::vector<unsigned int> hittypes_(std::begin(hittypes), std::end(hittypes));
0288   return matchedSimTrkIdxs(hitidxs_, hittypes_, verbose);
0289 }
0290 
0291 //___________________________________________________________________________________________________________________________________________________________________________________________
0292 std::vector<int> matchedSimTrkIdxs(std::vector<unsigned int> hitidxs,
0293                                    std::vector<unsigned int> hittypes,
0294                                    bool verbose,
0295                                    float *pmatched) {
0296   if (hitidxs.size() != hittypes.size()) {
0297     std::cout << "Error: matched_sim_trk_idxs()   hitidxs and hittypes have different lengths" << std::endl;
0298     std::cout << "hitidxs.size(): " << hitidxs.size() << std::endl;
0299     std::cout << "hittypes.size(): " << hittypes.size() << std::endl;
0300   }
0301 
0302   std::vector<std::pair<unsigned int, unsigned int>> to_check_duplicate;
0303   for (size_t i = 0; i < hitidxs.size(); ++i) {
0304     auto hitidx = hitidxs[i];
0305     auto hittype = hittypes[i];
0306     auto item = std::make_pair(hitidx, hittype);
0307     if (std::find(to_check_duplicate.begin(), to_check_duplicate.end(), item) == to_check_duplicate.end()) {
0308       to_check_duplicate.push_back(item);
0309     }
0310   }
0311 
0312   int nhits_input = to_check_duplicate.size();
0313 
0314   std::vector<std::vector<int>> simtrk_idxs;
0315   std::vector<int> unique_idxs;  // to aggregate which ones to count and test
0316 
0317   if (verbose) {
0318     std::cout << " '------------------------': "
0319               << "------------------------" << std::endl;
0320   }
0321 
0322   for (size_t ihit = 0; ihit < to_check_duplicate.size(); ++ihit) {
0323     auto ihitdata = to_check_duplicate[ihit];
0324     auto &&[hitidx, hittype] = ihitdata;
0325 
0326     if (verbose) {
0327       std::cout << " hitidx: " << hitidx << " hittype: " << hittype << std::endl;
0328     }
0329 
0330     std::vector<int> simtrk_idxs_per_hit;
0331 
0332     const std::vector<std::vector<int>> *simHitIdxs = hittype == 4 ? &trk.ph2_simHitIdx() : &trk.pix_simHitIdx();
0333 
0334     if (verbose) {
0335       std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl;
0336       std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl;
0337     }
0338 
0339     if (static_cast<const unsigned int>((*simHitIdxs).size()) <= hitidx) {
0340       std::cout << "ERROR" << std::endl;
0341       std::cout << " hittype: " << hittype << std::endl;
0342       std::cout << " trk.pix_simHitIdx().size(): " << trk.pix_simHitIdx().size() << std::endl;
0343       std::cout << " trk.ph2_simHitIdx().size(): " << trk.ph2_simHitIdx().size() << std::endl;
0344       std::cout << (*simHitIdxs).size() << " " << hittype << std::endl;
0345       std::cout << hitidx << " " << hittype << std::endl;
0346     }
0347 
0348     for (auto &simhit_idx : (*simHitIdxs).at(hitidx)) {
0349       if (static_cast<const int>(trk.simhit_simTrkIdx().size()) <= simhit_idx) {
0350         std::cout << (*simHitIdxs).size() << " " << hittype << std::endl;
0351         std::cout << hitidx << " " << hittype << std::endl;
0352         std::cout << trk.simhit_simTrkIdx().size() << " " << simhit_idx << std::endl;
0353       }
0354       int simtrk_idx = trk.simhit_simTrkIdx().at(simhit_idx);
0355       if (verbose) {
0356         std::cout << " hitidx: " << hitidx << " simhit_idx: " << simhit_idx << " simtrk_idx: " << simtrk_idx
0357                   << std::endl;
0358       }
0359       simtrk_idxs_per_hit.push_back(simtrk_idx);
0360       if (std::find(unique_idxs.begin(), unique_idxs.end(), simtrk_idx) == unique_idxs.end())
0361         unique_idxs.push_back(simtrk_idx);
0362     }
0363 
0364     if (simtrk_idxs_per_hit.size() == 0) {
0365       if (verbose) {
0366         std::cout << " hitidx: " << hitidx << " -1: " << -1 << std::endl;
0367       }
0368       simtrk_idxs_per_hit.push_back(-1);
0369       if (std::find(unique_idxs.begin(), unique_idxs.end(), -1) == unique_idxs.end())
0370         unique_idxs.push_back(-1);
0371     }
0372 
0373     simtrk_idxs.push_back(simtrk_idxs_per_hit);
0374   }
0375 
0376   if (verbose) {
0377     std::cout << " unique_idxs.size(): " << unique_idxs.size() << std::endl;
0378     for (auto &unique_idx : unique_idxs) {
0379       std::cout << " unique_idx: " << unique_idx << std::endl;
0380     }
0381   }
0382 
0383   // print
0384   if (verbose) {
0385     std::cout << "va print" << std::endl;
0386     for (auto &vec : simtrk_idxs) {
0387       for (auto &idx : vec) {
0388         std::cout << idx << " ";
0389       }
0390       std::cout << std::endl;
0391     }
0392     std::cout << "va print end" << std::endl;
0393   }
0394 
0395   // Compute all permutations
0396   std::function<void(std::vector<std::vector<int>> &, std::vector<int>, size_t, std::vector<std::vector<int>> &)> perm =
0397       [&](std::vector<std::vector<int>> &result,
0398           std::vector<int> intermediate,
0399           size_t n,
0400           std::vector<std::vector<int>> &va) {
0401         // std::cout <<  " 'called': " << "called" <<  std::endl;
0402         if (va.size() > n) {
0403           for (auto x : va[n]) {
0404             // std::cout <<  " n: " << n <<  std::endl;
0405             // std::cout <<  " intermediate.size(): " << intermediate.size() <<  std::endl;
0406             std::vector<int> copy_intermediate(intermediate);
0407             copy_intermediate.push_back(x);
0408             perm(result, copy_intermediate, n + 1, va);
0409           }
0410         } else {
0411           result.push_back(intermediate);
0412         }
0413       };
0414 
0415   std::vector<std::vector<int>> allperms;
0416   perm(allperms, std::vector<int>(), 0, simtrk_idxs);
0417 
0418   if (verbose) {
0419     std::cout << " allperms.size(): " << allperms.size() << std::endl;
0420     for (unsigned iperm = 0; iperm < allperms.size(); ++iperm) {
0421       std::cout << " allperms[iperm].size(): " << allperms[iperm].size() << std::endl;
0422       for (unsigned ielem = 0; ielem < allperms[iperm].size(); ++ielem) {
0423         std::cout << " allperms[iperm][ielem]: " << allperms[iperm][ielem] << std::endl;
0424       }
0425     }
0426   }
0427   int maxHitMatchCount = 0;  // ultimate maximum of the number of matched hits
0428   std::vector<int> matched_sim_trk_idxs;
0429   float max_percent_matched = 0.0f;
0430   for (auto &trkidx_perm : allperms) {
0431     std::vector<int> counts;
0432     for (auto &unique_idx : unique_idxs) {
0433       int cnt = std::count(trkidx_perm.begin(), trkidx_perm.end(), unique_idx);
0434       counts.push_back(cnt);
0435     }
0436     auto result = std::max_element(counts.begin(), counts.end());
0437     int rawidx = std::distance(counts.begin(), result);
0438     int trkidx = unique_idxs[rawidx];
0439     if (trkidx < 0)
0440       continue;
0441     float percent_matched = static_cast<float>(counts[rawidx]) / nhits_input;
0442     if (percent_matched > 0.75f)
0443       matched_sim_trk_idxs.push_back(trkidx);
0444     maxHitMatchCount = std::max(maxHitMatchCount, *std::max_element(counts.begin(), counts.end()));
0445     max_percent_matched = std::max(max_percent_matched, percent_matched);
0446   }
0447 
0448   // If pmatched is provided, set its value
0449   if (pmatched != nullptr) {
0450     *pmatched = max_percent_matched;
0451   }
0452 
0453   std::set<int> s;
0454   unsigned size = matched_sim_trk_idxs.size();
0455   for (unsigned i = 0; i < size; ++i)
0456     s.insert(matched_sim_trk_idxs[i]);
0457   matched_sim_trk_idxs.assign(s.begin(), s.end());
0458   return matched_sim_trk_idxs;
0459 }
0460 
0461 //__________________________________________________________________________________________
0462 int getDenomSimTrkType(int isimtrk) {
0463   if (isimtrk < 0)
0464     return 0;  // not a sim
0465   const int &q = trk.sim_q()[isimtrk];
0466   if (q == 0)
0467     return 1;  // sim
0468   const float &pt = trk.sim_pt()[isimtrk];
0469   const float &eta = trk.sim_eta()[isimtrk];
0470   if (pt < 1 or abs(eta) > 2.4)
0471     return 2;  // sim and charged
0472   const int &bunch = trk.sim_bunchCrossing()[isimtrk];
0473   const int &event = trk.sim_event()[isimtrk];
0474   const int &vtxIdx = trk.sim_parentVtxIdx()[isimtrk];
0475   const float &vtx_x = trk.simvtx_x()[vtxIdx];
0476   const float &vtx_y = trk.simvtx_y()[vtxIdx];
0477   const float &vtx_z = trk.simvtx_z()[vtxIdx];
0478   const float &vtx_perp = sqrt(vtx_x * vtx_x + vtx_y * vtx_y);
0479   if (vtx_perp > 2.5)
0480     return 3;  // pt > 1 and abs(eta) < 2.4
0481   if (abs(vtx_z) > 30)
0482     return 4;  // pt > 1 and abs(eta) < 2.4 and vtx < 2.5
0483   if (bunch != 0)
0484     return 5;  // pt > 1 and abs(eta) < 2.4 and vtx < 2.5 and vtx < 300
0485   if (event != 0)
0486     return 6;  // pt > 1 and abs(eta) < 2.4 and vtx 2.5/30 and bunch == 0
0487   return 7;    // pt > 1 and abs(eta) < 2.4 and vtx 2.5/30 and bunch == 0 and event == 0
0488 }
0489 
0490 //__________________________________________________________________________________________
0491 int getDenomSimTrkType(std::vector<int> simidxs) {
0492   int type = 0;
0493   for (auto &simidx : simidxs) {
0494     int this_type = getDenomSimTrkType(simidx);
0495     if (this_type > type) {
0496       type = this_type;
0497     }
0498   }
0499   return type;
0500 }
0501 
0502 //  ---------------------------------- =========================================== ----------------------------------------------
0503 //  ---------------------------------- =========================================== ----------------------------------------------
0504 //  ---------------------------------- =========================================== ----------------------------------------------
0505 //  ---------------------------------- =========================================== ----------------------------------------------
0506 //  ---------------------------------- =========================================== ----------------------------------------------
0507 
0508 //___________________________________________________________________________________________________________________________________________________________________________________________
0509 float drfracSimHitConsistentWithHelix(int isimtrk, int isimhitidx) {
0510   // Read track parameters
0511   float vx = trk.simvtx_x()[0];
0512   float vy = trk.simvtx_y()[0];
0513   float vz = trk.simvtx_z()[0];
0514   float pt = trk.sim_pt()[isimtrk];
0515   float eta = trk.sim_eta()[isimtrk];
0516   float phi = trk.sim_phi()[isimtrk];
0517   int charge = trk.sim_q()[isimtrk];
0518 
0519   // Construct helix object
0520   lst_math::Helix helix(pt, eta, phi, vx, vy, vz, charge);
0521 
0522   return drfracSimHitConsistentWithHelix(helix, isimhitidx);
0523 }
0524 
0525 //__________________________________________________________________________________________
0526 float drfracSimHitConsistentWithHelix(lst_math::Helix &helix, int isimhitidx) {
0527   // Sim hit vector
0528   std::vector<float> point = {trk.simhit_x()[isimhitidx], trk.simhit_y()[isimhitidx], trk.simhit_z()[isimhitidx]};
0529 
0530   // Inferring parameter t of helix parametric function via z position
0531   float t = helix.infer_t(point);
0532 
0533   // If the best fit is more than pi parameter away then it's a returning hit and should be ignored
0534   if (not(t <= M_PI))
0535     return 999;
0536 
0537   // Expected hit position with given z
0538   auto [x, y, z, r] = helix.get_helix_point(t);
0539 
0540   // ( expected_r - simhit_r ) / expected_r
0541   float drfrac = abs(helix.compare_radius(point)) / r;
0542 
0543   return drfrac;
0544 }
0545 
0546 //__________________________________________________________________________________________
0547 float distxySimHitConsistentWithHelix(int isimtrk, int isimhitidx) {
0548   // Read track parameters
0549   float vx = trk.simvtx_x()[0];
0550   float vy = trk.simvtx_y()[0];
0551   float vz = trk.simvtx_z()[0];
0552   float pt = trk.sim_pt()[isimtrk];
0553   float eta = trk.sim_eta()[isimtrk];
0554   float phi = trk.sim_phi()[isimtrk];
0555   int charge = trk.sim_q()[isimtrk];
0556 
0557   // Construct helix object
0558   lst_math::Helix helix(pt, eta, phi, vx, vy, vz, charge);
0559 
0560   return distxySimHitConsistentWithHelix(helix, isimhitidx);
0561 }
0562 
0563 //__________________________________________________________________________________________
0564 float distxySimHitConsistentWithHelix(lst_math::Helix &helix, int isimhitidx) {
0565   // Sim hit vector
0566   std::vector<float> point = {trk.simhit_x()[isimhitidx], trk.simhit_y()[isimhitidx], trk.simhit_z()[isimhitidx]};
0567 
0568   // Inferring parameter t of helix parametric function via z position
0569   float t = helix.infer_t(point);
0570 
0571   // If the best fit is more than pi parameter away then it's a returning hit and should be ignored
0572   if (not(t <= M_PI))
0573     return 999;
0574 
0575   // Expected hit position with given z
0576   //auto [x, y, z, r] = helix.get_helix_point(t);
0577 
0578   // ( expected_r - simhit_r ) / expected_r
0579   float distxy = helix.compare_xy(point);
0580 
0581   return distxy;
0582 }
0583 
0584 //__________________________________________________________________________________________
0585 TVector3 calculateR3FromPCA(const TVector3 &p3, const float dxy, const float dz) {
0586   const float pt = p3.Pt();
0587   const float p = p3.Mag();
0588   const float vz = dz * pt * pt / p / p;
0589 
0590   const float vx = -dxy * p3.y() / pt - p3.x() / p * p3.z() / p * dz;
0591   const float vy = dxy * p3.x() / pt - p3.y() / p * p3.z() / p * dz;
0592   return TVector3(vx, vy, vz);
0593 }
0594 
0595 //  ---------------------------------- =========================================== ----------------------------------------------
0596 //  ---------------------------------- =========================================== ----------------------------------------------
0597 //  ---------------------------------- =========================================== ----------------------------------------------
0598 //  ---------------------------------- =========================================== ----------------------------------------------
0599 //  ---------------------------------- =========================================== ----------------------------------------------
0600 
0601 //___________________________________________________________________________________________________________________________________________________________________________________________
0602 void addInputsToLineSegmentTrackingPreLoad(std::vector<std::vector<float>> &out_trkX,
0603                                            std::vector<std::vector<float>> &out_trkY,
0604                                            std::vector<std::vector<float>> &out_trkZ,
0605                                            std::vector<std::vector<unsigned int>> &out_hitId,
0606                                            std::vector<std::vector<unsigned int>> &out_hitIdxs,
0607                                            std::vector<std::vector<unsigned int>> &out_hitIndices_vec0,
0608                                            std::vector<std::vector<unsigned int>> &out_hitIndices_vec1,
0609                                            std::vector<std::vector<unsigned int>> &out_hitIndices_vec2,
0610                                            std::vector<std::vector<unsigned int>> &out_hitIndices_vec3,
0611                                            std::vector<std::vector<float>> &out_deltaPhi_vec,
0612                                            std::vector<std::vector<float>> &out_ptIn_vec,
0613                                            std::vector<std::vector<float>> &out_ptErr_vec,
0614                                            std::vector<std::vector<float>> &out_px_vec,
0615                                            std::vector<std::vector<float>> &out_py_vec,
0616                                            std::vector<std::vector<float>> &out_pz_vec,
0617                                            std::vector<std::vector<float>> &out_eta_vec,
0618                                            std::vector<std::vector<float>> &out_etaErr_vec,
0619                                            std::vector<std::vector<float>> &out_phi_vec,
0620                                            std::vector<std::vector<int>> &out_charge_vec,
0621                                            std::vector<std::vector<unsigned int>> &out_seedIdx_vec,
0622                                            std::vector<std::vector<int>> &out_superbin_vec,
0623                                            std::vector<std::vector<PixelType>> &out_pixelType_vec,
0624                                            std::vector<std::vector<char>> &out_isQuad_vec) {
0625   unsigned int count = 0;
0626   auto n_see = trk.see_stateTrajGlbPx().size();
0627   std::vector<float> px_vec;
0628   px_vec.reserve(n_see);
0629   std::vector<float> py_vec;
0630   py_vec.reserve(n_see);
0631   std::vector<float> pz_vec;
0632   pz_vec.reserve(n_see);
0633   std::vector<unsigned int> hitIndices_vec0;
0634   hitIndices_vec0.reserve(n_see);
0635   std::vector<unsigned int> hitIndices_vec1;
0636   hitIndices_vec1.reserve(n_see);
0637   std::vector<unsigned int> hitIndices_vec2;
0638   hitIndices_vec2.reserve(n_see);
0639   std::vector<unsigned int> hitIndices_vec3;
0640   hitIndices_vec3.reserve(n_see);
0641   std::vector<float> ptIn_vec;
0642   ptIn_vec.reserve(n_see);
0643   std::vector<float> ptErr_vec;
0644   ptErr_vec.reserve(n_see);
0645   std::vector<float> etaErr_vec;
0646   etaErr_vec.reserve(n_see);
0647   std::vector<float> eta_vec;
0648   eta_vec.reserve(n_see);
0649   std::vector<float> phi_vec;
0650   phi_vec.reserve(n_see);
0651   std::vector<int> charge_vec;
0652   charge_vec.reserve(n_see);
0653   std::vector<unsigned int> seedIdx_vec;
0654   seedIdx_vec.reserve(n_see);
0655   std::vector<float> deltaPhi_vec;
0656   deltaPhi_vec.reserve(n_see);
0657   std::vector<float> trkX = trk.ph2_x();
0658   std::vector<float> trkY = trk.ph2_y();
0659   std::vector<float> trkZ = trk.ph2_z();
0660   std::vector<unsigned int> hitId = trk.ph2_detId();
0661   std::vector<unsigned int> hitIdxs(trk.ph2_detId().size());
0662 
0663   std::vector<int> superbin_vec;
0664   std::vector<PixelType> pixelType_vec;
0665   std::vector<char> isQuad_vec;
0666   std::iota(hitIdxs.begin(), hitIdxs.end(), 0);
0667   const int hit_size = trkX.size();
0668 
0669   for (size_t iSeed = 0; iSeed < trk.see_stateTrajGlbPx().size(); ++iSeed) {
0670     //// track algorithm; partial copy from TrackBase.h
0671     // enum class TrackAlgorithm {
0672     //    undefAlgorithm = 0,
0673     //    ctf = 1,
0674     //    duplicateMerge = 2,
0675     //    cosmics = 3,
0676     //    initialStep = 4,
0677     //    lowPtTripletStep = 5,
0678     //    pixelPairStep = 6,
0679     //    detachedTripletStep = 7,
0680     //    mixedTripletStep = 8,
0681     //    pixelLessStep = 9,
0682     //    tobTecStep = 10,
0683     //    jetCoreRegionalStep = 11,
0684     //    conversionStep = 12,
0685     //    muonSeededStepInOut = 13,
0686     //    muonSeededStepOutIn = 14,
0687     //    outInEcalSeededConv = 15, inOutEcalSeededConv = 16,
0688     //    nuclInter = 17,
0689     //    standAloneMuon = 18, globalMuon = 19, cosmicStandAloneMuon = 20, cosmicGlobalMuon = 21,
0690     //    // Phase1
0691     //    highPtTripletStep = 22, lowPtQuadStep = 23, detachedQuadStep = 24,
0692     //    reservedForUpgrades1 = 25, reservedForUpgrades2 = 26,
0693     //    bTagGhostTracks = 27,
0694     //    beamhalo = 28,
0695     //    gsf = 29
0696     //};
0697     bool good_seed_type = false;
0698     if (trk.see_algo()[iSeed] == 4)
0699       good_seed_type = true;
0700     // if (trk.see_algo()[iSeed] == 5) good_seed_type = true;
0701     // if (trk.see_algo()[iSeed] == 7) good_seed_type = true;
0702     if (trk.see_algo()[iSeed] == 22)
0703       good_seed_type = true;
0704     // if (trk.see_algo()[iSeed] == 23) good_seed_type = true;
0705     // if (trk.see_algo()[iSeed] == 24) good_seed_type = true;
0706     if (not good_seed_type)
0707       continue;
0708 
0709     TVector3 p3LH(trk.see_stateTrajGlbPx()[iSeed], trk.see_stateTrajGlbPy()[iSeed], trk.see_stateTrajGlbPz()[iSeed]);
0710     float ptIn = p3LH.Pt();
0711     float eta = p3LH.Eta();
0712     float ptErr = trk.see_ptErr()[iSeed];
0713 
0714     if ((ptIn > ana.ptCut - 2 * ptErr)) {
0715       TVector3 r3LH(trk.see_stateTrajGlbX()[iSeed], trk.see_stateTrajGlbY()[iSeed], trk.see_stateTrajGlbZ()[iSeed]);
0716       TVector3 p3PCA(trk.see_px()[iSeed], trk.see_py()[iSeed], trk.see_pz()[iSeed]);
0717       TVector3 r3PCA(calculateR3FromPCA(p3PCA, trk.see_dxy()[iSeed], trk.see_dz()[iSeed]));
0718       TVector3 seedSD_mdRef_r3 = r3PCA;
0719       TVector3 seedSD_mdOut_r3 = r3LH;
0720       TVector3 seedSD_r3 = r3LH;
0721       TVector3 seedSD_p3 = p3LH;
0722 
0723       // The charge could be used directly in the line below
0724       float pixelSegmentDeltaPhiChange = r3LH.DeltaPhi(p3LH);
0725       float etaErr = trk.see_etaErr()[iSeed];
0726       float px = p3LH.X();
0727       float py = p3LH.Y();
0728       float pz = p3LH.Z();
0729       int charge = trk.see_q()[iSeed];
0730       unsigned int seedIdx = iSeed;
0731 
0732       PixelType pixtype = PixelType::kInvalid;
0733       if (ptIn >= 2.0) {
0734         pixtype = PixelType::kHighPt;
0735       } else if (ptIn >= (ana.ptCut - 2 * ptErr) and ptIn < 2.0) {
0736         if (pixelSegmentDeltaPhiChange >= 0) {
0737           pixtype = PixelType::kLowPtPosCurv;
0738         } else {
0739           pixtype = PixelType::kLowPtNegCurv;
0740         }
0741       } else {
0742         continue;
0743       }
0744 
0745       unsigned int hitIdx0 = hit_size + count;
0746       count++;
0747 
0748       unsigned int hitIdx1 = hit_size + count;
0749       count++;
0750 
0751       unsigned int hitIdx2 = hit_size + count;
0752       count++;
0753 
0754       unsigned int hitIdx3;
0755       if (trk.see_hitIdx()[iSeed].size() <= 3) {
0756         hitIdx3 = hitIdx2;
0757       } else {
0758         hitIdx3 = hit_size + count;
0759         count++;
0760       }
0761 
0762       trkX.push_back(r3PCA.X());
0763       trkY.push_back(r3PCA.Y());
0764       trkZ.push_back(r3PCA.Z());
0765       trkX.push_back(p3PCA.Pt());
0766       float p3PCA_Eta = p3PCA.Eta();
0767       trkY.push_back(p3PCA_Eta);
0768       float p3PCA_Phi = p3PCA.Phi();
0769       trkZ.push_back(p3PCA_Phi);
0770       trkX.push_back(r3LH.X());
0771       trkY.push_back(r3LH.Y());
0772       trkZ.push_back(r3LH.Z());
0773       hitId.push_back(1);
0774       hitId.push_back(1);
0775       hitId.push_back(1);
0776       if (trk.see_hitIdx()[iSeed].size() > 3) {
0777         trkX.push_back(r3LH.X());
0778         trkY.push_back(trk.see_dxy()[iSeed]);
0779         trkZ.push_back(trk.see_dz()[iSeed]);
0780         hitId.push_back(1);
0781       }
0782       px_vec.push_back(px);
0783       py_vec.push_back(py);
0784       pz_vec.push_back(pz);
0785 
0786       hitIndices_vec0.push_back(hitIdx0);
0787       hitIndices_vec1.push_back(hitIdx1);
0788       hitIndices_vec2.push_back(hitIdx2);
0789       hitIndices_vec3.push_back(hitIdx3);
0790       ptIn_vec.push_back(ptIn);
0791       ptErr_vec.push_back(ptErr);
0792       etaErr_vec.push_back(etaErr);
0793       eta_vec.push_back(eta);
0794       float phi = p3LH.Phi();
0795       phi_vec.push_back(phi);
0796       charge_vec.push_back(charge);
0797       seedIdx_vec.push_back(seedIdx);
0798       deltaPhi_vec.push_back(pixelSegmentDeltaPhiChange);
0799 
0800       // For matching with sim tracks
0801       hitIdxs.push_back(trk.see_hitIdx()[iSeed][0]);
0802       hitIdxs.push_back(trk.see_hitIdx()[iSeed][1]);
0803       hitIdxs.push_back(trk.see_hitIdx()[iSeed][2]);
0804       char isQuad = false;
0805       if (trk.see_hitIdx()[iSeed].size() > 3) {
0806         isQuad = true;
0807         hitIdxs.push_back(trk.see_hitIdx()[iSeed].size() > 3 ? trk.see_hitIdx()[iSeed][3] : trk.see_hitIdx()[iSeed][2]);
0808       }
0809       // if (pt < 0){ ptbin = 0;}
0810       float neta = 25.;
0811       float nphi = 72.;
0812       float nz = 25.;
0813       int etabin = (p3PCA_Eta + 2.6) / ((2 * 2.6) / neta);
0814       int phibin = (p3PCA_Phi + 3.14159265358979323846) / ((2. * 3.14159265358979323846) / nphi);
0815       int dzbin = (trk.see_dz()[iSeed] + 30) / (2 * 30 / nz);
0816       int isuperbin =
0817           /*(nz * nphi * neta) * ptbin + (removed since pt bin is determined by pixelType)*/ (nz * nphi) * etabin +
0818           (nz)*phibin + dzbin;
0819       // if(isuperbin<0 || isuperbin>=44900){printf("isuperbin %d %d %d %d %f\n",isuperbin,etabin,phibin,dzbin,p3PCA.Eta());}
0820       superbin_vec.push_back(isuperbin);
0821       pixelType_vec.push_back(pixtype);
0822       isQuad_vec.push_back(isQuad);
0823     }
0824   }
0825 
0826   out_trkX.push_back(trkX);
0827   out_trkY.push_back(trkY);
0828   out_trkZ.push_back(trkZ);
0829   out_hitId.push_back(hitId);
0830   out_hitIdxs.push_back(hitIdxs);
0831   out_hitIndices_vec0.push_back(hitIndices_vec0);
0832   out_hitIndices_vec1.push_back(hitIndices_vec1);
0833   out_hitIndices_vec2.push_back(hitIndices_vec2);
0834   out_hitIndices_vec3.push_back(hitIndices_vec3);
0835   out_deltaPhi_vec.push_back(deltaPhi_vec);
0836   out_ptIn_vec.push_back(ptIn_vec);
0837   out_ptErr_vec.push_back(ptErr_vec);
0838   out_px_vec.push_back(px_vec);
0839   out_py_vec.push_back(py_vec);
0840   out_pz_vec.push_back(pz_vec);
0841   out_eta_vec.push_back(eta_vec);
0842   out_etaErr_vec.push_back(etaErr_vec);
0843   out_phi_vec.push_back(phi_vec);
0844   out_charge_vec.push_back(charge_vec);
0845   out_seedIdx_vec.push_back(seedIdx_vec);
0846   out_superbin_vec.push_back(superbin_vec);
0847   out_pixelType_vec.push_back(pixelType_vec);
0848   out_isQuad_vec.push_back(isQuad_vec);
0849 
0850   //    float hit_loading_elapsed = my_timer.RealTime();
0851   //    if (ana.verbose >= 2) std::cout << "Loading inputs processing time: " << hit_loading_elapsed << " secs" << std::endl;
0852   //    return hit_loading_elapsed;
0853 }
0854 
0855 //___________________________________________________________________________________________________________________________________________________________________________________________
0856 float addInputsToEventPreLoad(LSTEvent *event,
0857                               bool useOMP,
0858                               std::vector<float> trkX,
0859                               std::vector<float> trkY,
0860                               std::vector<float> trkZ,
0861                               std::vector<unsigned int> hitId,
0862                               std::vector<unsigned int> hitIdxs,
0863                               std::vector<unsigned int> hitIndices_vec0,
0864                               std::vector<unsigned int> hitIndices_vec1,
0865                               std::vector<unsigned int> hitIndices_vec2,
0866                               std::vector<unsigned int> hitIndices_vec3,
0867                               std::vector<float> deltaPhi_vec,
0868                               std::vector<float> ptIn_vec,
0869                               std::vector<float> ptErr_vec,
0870                               std::vector<float> px_vec,
0871                               std::vector<float> py_vec,
0872                               std::vector<float> pz_vec,
0873                               std::vector<float> eta_vec,
0874                               std::vector<float> etaErr_vec,
0875                               std::vector<float> phi_vec,
0876                               std::vector<int> charge_vec,
0877                               std::vector<unsigned int> seedIdx_vec,
0878                               std::vector<int> superbin_vec,
0879                               std::vector<PixelType> pixelType_vec,
0880                               std::vector<char> isQuad_vec) {
0881   TStopwatch my_timer;
0882 
0883   if (ana.verbose >= 2)
0884     std::cout << "Loading Inputs (i.e. outer tracker hits, and pixel line segements) to the Line Segment Tracking.... "
0885               << std::endl;
0886 
0887   my_timer.Start();
0888 
0889   event->addHitToEvent(trkX, trkY, trkZ, hitId, hitIdxs);
0890 
0891   event->addPixelSegmentToEvent(hitIndices_vec0,
0892                                 hitIndices_vec1,
0893                                 hitIndices_vec2,
0894                                 hitIndices_vec3,
0895                                 deltaPhi_vec,
0896                                 ptIn_vec,
0897                                 ptErr_vec,
0898                                 px_vec,
0899                                 py_vec,
0900                                 pz_vec,
0901                                 eta_vec,
0902                                 etaErr_vec,
0903                                 phi_vec,
0904                                 charge_vec,
0905                                 seedIdx_vec,
0906                                 superbin_vec,
0907                                 pixelType_vec,
0908                                 isQuad_vec);
0909   event->wait();  // device side event calls are asynchronous: wait to measure time or print
0910   float hit_loading_elapsed = my_timer.RealTime();
0911 
0912   if (ana.verbose >= 2)
0913     std::cout << "Loading inputs processing time: " << hit_loading_elapsed << " secs" << std::endl;
0914 
0915   return hit_loading_elapsed;
0916 }
0917 
0918 //________________________________________________________________________________________________________________________________
0919 void printTimingInformation(std::vector<std::vector<float>> &timing_information, float fullTime, float fullavg) {
0920   if (ana.verbose == 0)
0921     return;
0922 
0923   std::cout << std::showpoint;
0924   std::cout << std::fixed;
0925   std::cout << std::setprecision(2);
0926   std::cout << std::right;
0927   std::cout << "Timing summary" << std::endl;
0928   std::cout << std::setw(6) << "Evt";
0929   std::cout << "   " << std::setw(6) << "Hits";
0930   std::cout << "   " << std::setw(6) << "MD";
0931   std::cout << "   " << std::setw(6) << "LS";
0932   std::cout << "   " << std::setw(6) << "T3";
0933   std::cout << "   " << std::setw(6) << "T5";
0934   std::cout << "   " << std::setw(6) << "pLS";
0935   std::cout << "   " << std::setw(6) << "pT5";
0936   std::cout << "   " << std::setw(6) << "pT3";
0937   std::cout << "   " << std::setw(6) << "TC";
0938   std::cout << "   " << std::setw(6) << "Reset";
0939   std::cout << "   " << std::setw(7) << "Total";
0940   std::cout << "   " << std::setw(7) << "Total(short)";
0941   std::cout << std::endl;
0942   std::vector<float> timing_sum_information(timing_information[0].size());
0943   std::vector<float> timing_shortlist;
0944   std::vector<float> timing_list;
0945   for (size_t ievt = 0; ievt < timing_information.size(); ++ievt) {
0946     auto timing = timing_information[ievt];
0947     float timing_total = 0.f;
0948     float timing_total_short = 0.f;
0949     timing_total += timing[0] * 1000;        // Hits
0950     timing_total += timing[1] * 1000;        // MD
0951     timing_total += timing[2] * 1000;        // LS
0952     timing_total += timing[3] * 1000;        // T3
0953     timing_total += timing[4] * 1000;        // T5
0954     timing_total += timing[5] * 1000;        // pLS
0955     timing_total += timing[6] * 1000;        // pT5
0956     timing_total += timing[7] * 1000;        // pT3
0957     timing_total += timing[8] * 1000;        // TC
0958     timing_total_short += timing[1] * 1000;  // MD
0959     timing_total_short += timing[2] * 1000;  // LS
0960     timing_total_short += timing[3] * 1000;  // T3
0961     timing_total_short += timing[4] * 1000;  // T5
0962     timing_total_short += timing[6] * 1000;  // pT5
0963     timing_total_short += timing[7] * 1000;  // pT3
0964     timing_total_short += timing[8] * 1000;  // TC
0965     timing_total_short += timing[9] * 1000;  // Reset
0966     std::cout << std::setw(6) << ievt;
0967     std::cout << "   " << std::setw(6) << timing[0] * 1000;    // Hits
0968     std::cout << "   " << std::setw(6) << timing[1] * 1000;    // MD
0969     std::cout << "   " << std::setw(6) << timing[2] * 1000;    // LS
0970     std::cout << "   " << std::setw(6) << timing[3] * 1000;    // T3
0971     std::cout << "   " << std::setw(6) << timing[4] * 1000;    // T5
0972     std::cout << "   " << std::setw(6) << timing[5] * 1000;    // pLS
0973     std::cout << "   " << std::setw(6) << timing[6] * 1000;    // pT5
0974     std::cout << "   " << std::setw(6) << timing[7] * 1000;    // pT3
0975     std::cout << "   " << std::setw(6) << timing[8] * 1000;    // TC
0976     std::cout << "   " << std::setw(6) << timing[9] * 1000;    // Reset
0977     std::cout << "   " << std::setw(7) << timing_total;        // Total time
0978     std::cout << "   " << std::setw(7) << timing_total_short;  // Total time
0979     std::cout << std::endl;
0980     timing_sum_information[0] += timing[0] * 1000;   // Hits
0981     timing_sum_information[1] += timing[1] * 1000;   // MD
0982     timing_sum_information[2] += timing[2] * 1000;   // LS
0983     timing_sum_information[3] += timing[3] * 1000;   // T3
0984     timing_sum_information[4] += timing[4] * 1000;   // T5
0985     timing_sum_information[5] += timing[5] * 1000;   // pLS
0986     timing_sum_information[6] += timing[6] * 1000;   // pT5
0987     timing_sum_information[7] += timing[7] * 1000;   // pT3
0988     timing_sum_information[8] += timing[8] * 1000;   // TC
0989     timing_sum_information[9] += timing[9] * 1000;   // Reset
0990     timing_shortlist.push_back(timing_total_short);  // short total
0991     timing_list.push_back(timing_total);             // short total
0992   }
0993   timing_sum_information[0] /= timing_information.size();  // Hits
0994   timing_sum_information[1] /= timing_information.size();  // MD
0995   timing_sum_information[2] /= timing_information.size();  // LS
0996   timing_sum_information[3] /= timing_information.size();  // T3
0997   timing_sum_information[4] /= timing_information.size();  // T5
0998   timing_sum_information[5] /= timing_information.size();  // pLS
0999   timing_sum_information[6] /= timing_information.size();  // pT5
1000   timing_sum_information[7] /= timing_information.size();  // pT3
1001   timing_sum_information[8] /= timing_information.size();  // TC
1002   timing_sum_information[9] /= timing_information.size();  // Reset
1003 
1004   float timing_total_avg = 0.0;
1005   timing_total_avg += timing_sum_information[0];  // Hits
1006   timing_total_avg += timing_sum_information[1];  // MD
1007   timing_total_avg += timing_sum_information[2];  // LS
1008   timing_total_avg += timing_sum_information[3];  // T3
1009   timing_total_avg += timing_sum_information[4];  // T5
1010   timing_total_avg += timing_sum_information[5];  // pLS
1011   timing_total_avg += timing_sum_information[6];  // pT5
1012   timing_total_avg += timing_sum_information[7];  // pT3
1013   timing_total_avg += timing_sum_information[8];  // TC
1014   timing_total_avg += timing_sum_information[9];  // Reset
1015   float timing_totalshort_avg = 0.0;
1016   timing_totalshort_avg += timing_sum_information[1];  // MD
1017   timing_totalshort_avg += timing_sum_information[2];  // LS
1018   timing_totalshort_avg += timing_sum_information[3];  // T3
1019   timing_totalshort_avg += timing_sum_information[4];  // T5
1020   timing_totalshort_avg += timing_sum_information[6];  // pT5
1021   timing_totalshort_avg += timing_sum_information[7];  // pT3
1022   timing_totalshort_avg += timing_sum_information[8];  // TC
1023   timing_totalshort_avg += timing_sum_information[9];  // Reset
1024 
1025   float standardDeviation = 0.0;
1026   for (auto shorttime : timing_shortlist) {
1027     standardDeviation += pow(shorttime - timing_totalshort_avg, 2);
1028   }
1029   float stdDev = sqrt(standardDeviation / timing_shortlist.size());
1030 
1031   std::cout << std::setprecision(1);
1032   std::cout << std::setw(6) << "Evt";
1033   std::cout << "   " << std::setw(6) << "Hits";
1034   std::cout << "   " << std::setw(6) << "MD";
1035   std::cout << "   " << std::setw(6) << "LS";
1036   std::cout << "   " << std::setw(6) << "T3";
1037   std::cout << "   " << std::setw(6) << "T5";
1038   std::cout << "   " << std::setw(6) << "pLS";
1039   std::cout << "   " << std::setw(6) << "pT5";
1040   std::cout << "   " << std::setw(6) << "pT3";
1041   std::cout << "   " << std::setw(6) << "TC";
1042   std::cout << "   " << std::setw(6) << "Reset";
1043   std::cout << "   " << std::setw(7) << "Total";
1044   std::cout << "   " << std::setw(7) << "Total(short)";
1045   std::cout << std::endl;
1046   std::cout << std::setw(6) << "avg";
1047   std::cout << "   " << std::setw(6) << timing_sum_information[0];  // Hits
1048   std::cout << "   " << std::setw(6) << timing_sum_information[1];  // MD
1049   std::cout << "   " << std::setw(6) << timing_sum_information[2];  // LS
1050   std::cout << "   " << std::setw(6) << timing_sum_information[3];  // T3
1051   std::cout << "   " << std::setw(6) << timing_sum_information[4];  // T5
1052   std::cout << "   " << std::setw(6) << timing_sum_information[5];  // pLS
1053   std::cout << "   " << std::setw(6) << timing_sum_information[6];  // pT5
1054   std::cout << "   " << std::setw(6) << timing_sum_information[7];  // pT3
1055   std::cout << "   " << std::setw(6) << timing_sum_information[8];  // TC
1056   std::cout << "   " << std::setw(6) << timing_sum_information[9];  // Reset
1057   std::cout << "   " << std::setw(7) << timing_total_avg;           // Average total time
1058   std::cout << "   " << std::setw(7) << timing_totalshort_avg;      // Average total time
1059   std::cout << "+/- " << std::setw(4) << stdDev;
1060   std::cout << "   " << std::setw(7) << fullavg;  // Average full time
1061   std::cout << "   " << ana.compilation_target;
1062   std::cout << "[s=" << ana.streams << "]";
1063   std::cout << std::endl;
1064 
1065   std::cout << std::left;
1066 }
1067 
1068 //  ---------------------------------- =========================================== ----------------------------------------------
1069 //  ---------------------------------- =========================================== ----------------------------------------------
1070 //  ---------------------------------- =========================================== ----------------------------------------------
1071 //  ---------------------------------- =========================================== ----------------------------------------------
1072 //  ---------------------------------- =========================================== ----------------------------------------------
1073 
1074 //__________________________________________________________________________________________
1075 TString get_absolute_path_after_check_file_exists(const std::string name) {
1076   std::filesystem::path fullpath = std::filesystem::absolute(name.c_str());
1077   // std::cout << "Checking file path = " << fullpath << std::endl;
1078   // std::cout <<  " fullpath.string().c_str(): " << fullpath.string().c_str() <<  std::endl;
1079   if (not std::filesystem::exists(fullpath)) {
1080     std::cout << "ERROR: Could not find the file = " << fullpath << std::endl;
1081     exit(2);
1082   }
1083   return TString(fullpath.string().c_str());
1084 }
1085 
1086 //_______________________________________________________________________________
1087 void writeMetaData() {
1088   // Write out metadata of the code to the output_tfile
1089   ana.output_tfile->cd();
1090   gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && echo '' && (cd - > /dev/null) ) > %s.gitversion.txt ",
1091                                 ana.output_tfile->GetName()));
1092   gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && git rev-parse HEAD && (cd - > /dev/null)) >> %s.gitversion.txt",
1093                                 ana.output_tfile->GetName()));
1094   gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && echo 'git status' && (cd - > /dev/null)) >> %s.gitversion.txt",
1095                                 ana.output_tfile->GetName()));
1096   gSystem->Exec(
1097       TString::Format("(cd $TRACKLOOPERDIR && git  --no-pager status && (cd - > /dev/null)) >> %s.gitversion.txt",
1098                       ana.output_tfile->GetName()));
1099   gSystem->Exec(TString::Format(
1100       "(cd $TRACKLOOPERDIR && echo 'git --no-pager log -n 100' && (cd - > /dev/null)) >> %s.gitversion.txt",
1101       ana.output_tfile->GetName()));
1102   gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && echo 'git diff' && (cd - > /dev/null)) >> %s.gitversion.txt",
1103                                 ana.output_tfile->GetName()));
1104   gSystem->Exec(
1105       TString::Format("(cd $TRACKLOOPERDIR && git --no-pager diff  && (cd - > /dev/null)) >> %s.gitversion.txt",
1106                       ana.output_tfile->GetName()));
1107 
1108   // Write gitversion info
1109   std::ifstream t(TString::Format("%s.gitversion.txt", ana.output_tfile->GetName()));
1110   std::string str((std::istreambuf_iterator<char>(t)), std::istreambuf_iterator<char>());
1111   TNamed code_tag_data("code_tag_data", str.c_str());
1112   ana.output_tfile->cd();
1113   code_tag_data.Write();
1114   gSystem->Exec(TString::Format("rm %s.gitversion.txt", ana.output_tfile->GetName()));
1115 
1116   // Write make script log
1117   TString make_log_path = TString::Format("%s/.make.log", ana.track_looper_dir_path.Data());
1118   std::ifstream makelog(make_log_path.Data());
1119   std::string makestr((std::istreambuf_iterator<char>(makelog)), std::istreambuf_iterator<char>());
1120   TNamed make_log("make_log", makestr.c_str());
1121   make_log.Write();
1122 
1123   // Write git diff output in a separate string to gauge the difference
1124   gSystem->Exec(TString::Format("(cd $TRACKLOOPERDIR && git --no-pager diff  && (cd - > /dev/null)) > %s.gitdiff.txt",
1125                                 ana.output_tfile->GetName()));
1126   std::ifstream gitdiff(TString::Format("%s.gitdiff.txt", ana.output_tfile->GetName()));
1127   std::string strgitdiff((std::istreambuf_iterator<char>(gitdiff)), std::istreambuf_iterator<char>());
1128   TNamed gitdifftnamed("gitdiff", strgitdiff.c_str());
1129   gitdifftnamed.Write();
1130   gSystem->Exec(TString::Format("rm %s.gitdiff.txt", ana.output_tfile->GetName()));
1131 
1132   // Write Parse from makestr the TARGET
1133   TString maketstr = makestr.c_str();
1134   TString rawstrdata = maketstr.ReplaceAll("MAKETARGET=", "%");
1135   TString targetrawdata = RooUtil::StringUtil::rsplit(rawstrdata, "%")[1];
1136   TString targetdata = RooUtil::StringUtil::split(targetrawdata)[0];
1137   ana.compilation_target = targetdata.Data();
1138 
1139   // Write out input sample or file name
1140   TNamed input("input", ana.input_raw_string.Data());
1141   input.Write();
1142 
1143   // Write the full command line used
1144   TNamed full_cmd_line("full_cmd_line", ana.full_cmd_line.Data());
1145   full_cmd_line.Write();
1146 
1147   // Write the TRACKLOOPERDIR
1148   TNamed tracklooper_path("tracklooper_path", ana.track_looper_dir_path.Data());
1149   tracklooper_path.Write();
1150 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!