Back to home page

Project CMSSW displayed by LXR

 
 

    


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   // 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->addPixelSegmentToEventFinalize();
0225   event->pixelLineSegmentCleaning(no_pls_dupclean);
0226   event->wait();  // device side event calls are asynchronous: wait to measure time or print
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();  // device side event calls are asynchronous: wait to measure time or print
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();  // device side event calls are asynchronous: wait to measure time or print
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;  // to aggregate which ones to count and test
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   // print
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   // Compute all permutations
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         // std::cout <<  " 'called': " << "called" <<  std::endl;
0399         if (va.size() > n) {
0400           for (auto x : va[n]) {
0401             // std::cout <<  " n: " << n <<  std::endl;
0402             // std::cout <<  " intermediate.size(): " << intermediate.size() <<  std::endl;
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;  // ultimate maximum of the number of matched hits
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   // If pmatched is provided, set its value
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;  // not a sim
0471   const int& q = trk_sim_q[isimtrk];
0472   if (q == 0)
0473     return 1;  // sim
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;  // sim and charged
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;  // pt > 1 and abs(eta) < 2.4
0487   if (std::abs(vtx_z) > 30)
0488     return 4;  // pt > 1 and abs(eta) < 2.4 and vtx < 2.5
0489   if (bunch != 0)
0490     return 5;  // pt > 1 and abs(eta) < 2.4 and vtx < 2.5 and vtx < 300
0491   if (event != 0)
0492     return 6;  // pt > 1 and abs(eta) < 2.4 and vtx 2.5/30 and bunch == 0
0493   return 7;    // pt > 1 and abs(eta) < 2.4 and vtx 2.5/30 and bunch == 0 and event == 0
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   // Read track parameters
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   // Construct helix object
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   // Sim hit vector
0567   std::vector<float> point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]};
0568 
0569   // Inferring parameter t of helix parametric function via z position
0570   float t = helix.infer_t(point);
0571 
0572   // If the best fit is more than pi parameter away then it's a returning hit and should be ignored
0573   if (not(t <= M_PI))
0574     return 999;
0575 
0576   // Expected hit position with given z
0577   auto [x, y, z, r] = helix.get_helix_point(t);
0578 
0579   // ( expected_r - simhit_r ) / expected_r
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   // Read track parameters
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   // Construct helix object
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   // Sim hit vector
0620   std::vector<float> point = {trk_simhit_x[isimhitidx], trk_simhit_y[isimhitidx], trk_simhit_z[isimhitidx]};
0621 
0622   // Inferring parameter t of helix parametric function via z position
0623   float t = helix.infer_t(point);
0624 
0625   // If the best fit is more than pi parameter away then it's a returning hit and should be ignored
0626   if (not(t <= M_PI))
0627     return 999;
0628 
0629   // Expected hit position with given z
0630   //auto [x, y, z, r] = helix.get_helix_point(t);
0631 
0632   // ( expected_r - simhit_r ) / expected_r
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   // We can't use CopyToDevice because the device can be DevHost
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();  // device side event calls are asynchronous: wait to measure time or print
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;        // Hits
0717     timing_total += timing[1] * 1000;        // MD
0718     timing_total += timing[2] * 1000;        // LS
0719     timing_total += timing[3] * 1000;        // T3
0720     timing_total += timing[4] * 1000;        // T5
0721     timing_total += timing[5] * 1000;        // pLS
0722     timing_total += timing[6] * 1000;        // pT5
0723     timing_total += timing[7] * 1000;        // pT3
0724     timing_total += timing[8] * 1000;        // TC
0725     timing_total_short += timing[1] * 1000;  // MD
0726     timing_total_short += timing[2] * 1000;  // LS
0727     timing_total_short += timing[3] * 1000;  // T3
0728     timing_total_short += timing[4] * 1000;  // T5
0729     timing_total_short += timing[6] * 1000;  // pT5
0730     timing_total_short += timing[7] * 1000;  // pT3
0731     timing_total_short += timing[8] * 1000;  // TC
0732     timing_total_short += timing[9] * 1000;  // Reset
0733     std::cout << std::setw(6) << ievt;
0734     std::cout << "   " << std::setw(6) << timing[0] * 1000;    // Hits
0735     std::cout << "   " << std::setw(6) << timing[1] * 1000;    // MD
0736     std::cout << "   " << std::setw(6) << timing[2] * 1000;    // LS
0737     std::cout << "   " << std::setw(6) << timing[3] * 1000;    // T3
0738     std::cout << "   " << std::setw(6) << timing[4] * 1000;    // T5
0739     std::cout << "   " << std::setw(6) << timing[5] * 1000;    // pLS
0740     std::cout << "   " << std::setw(6) << timing[6] * 1000;    // pT5
0741     std::cout << "   " << std::setw(6) << timing[7] * 1000;    // pT3
0742     std::cout << "   " << std::setw(6) << timing[8] * 1000;    // TC
0743     std::cout << "   " << std::setw(6) << timing[9] * 1000;    // Reset
0744     std::cout << "   " << std::setw(7) << timing_total;        // Total time
0745     std::cout << "   " << std::setw(7) << timing_total_short;  // Total time
0746     std::cout << std::endl;
0747     timing_sum_information[0] += timing[0] * 1000;   // Hits
0748     timing_sum_information[1] += timing[1] * 1000;   // MD
0749     timing_sum_information[2] += timing[2] * 1000;   // LS
0750     timing_sum_information[3] += timing[3] * 1000;   // T3
0751     timing_sum_information[4] += timing[4] * 1000;   // T5
0752     timing_sum_information[5] += timing[5] * 1000;   // pLS
0753     timing_sum_information[6] += timing[6] * 1000;   // pT5
0754     timing_sum_information[7] += timing[7] * 1000;   // pT3
0755     timing_sum_information[8] += timing[8] * 1000;   // TC
0756     timing_sum_information[9] += timing[9] * 1000;   // Reset
0757     timing_shortlist.push_back(timing_total_short);  // short total
0758     timing_list.push_back(timing_total);             // short total
0759   }
0760   timing_sum_information[0] /= timing_information.size();  // Hits
0761   timing_sum_information[1] /= timing_information.size();  // MD
0762   timing_sum_information[2] /= timing_information.size();  // LS
0763   timing_sum_information[3] /= timing_information.size();  // T3
0764   timing_sum_information[4] /= timing_information.size();  // T5
0765   timing_sum_information[5] /= timing_information.size();  // pLS
0766   timing_sum_information[6] /= timing_information.size();  // pT5
0767   timing_sum_information[7] /= timing_information.size();  // pT3
0768   timing_sum_information[8] /= timing_information.size();  // TC
0769   timing_sum_information[9] /= timing_information.size();  // Reset
0770 
0771   float timing_total_avg = 0.0;
0772   timing_total_avg += timing_sum_information[0];  // Hits
0773   timing_total_avg += timing_sum_information[1];  // MD
0774   timing_total_avg += timing_sum_information[2];  // LS
0775   timing_total_avg += timing_sum_information[3];  // T3
0776   timing_total_avg += timing_sum_information[4];  // T5
0777   timing_total_avg += timing_sum_information[5];  // pLS
0778   timing_total_avg += timing_sum_information[6];  // pT5
0779   timing_total_avg += timing_sum_information[7];  // pT3
0780   timing_total_avg += timing_sum_information[8];  // TC
0781   timing_total_avg += timing_sum_information[9];  // Reset
0782   float timing_totalshort_avg = 0.0;
0783   timing_totalshort_avg += timing_sum_information[1];  // MD
0784   timing_totalshort_avg += timing_sum_information[2];  // LS
0785   timing_totalshort_avg += timing_sum_information[3];  // T3
0786   timing_totalshort_avg += timing_sum_information[4];  // T5
0787   timing_totalshort_avg += timing_sum_information[6];  // pT5
0788   timing_totalshort_avg += timing_sum_information[7];  // pT3
0789   timing_totalshort_avg += timing_sum_information[8];  // TC
0790   timing_totalshort_avg += timing_sum_information[9];  // Reset
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];  // Hits
0815   std::cout << "   " << std::setw(6) << timing_sum_information[1];  // MD
0816   std::cout << "   " << std::setw(6) << timing_sum_information[2];  // LS
0817   std::cout << "   " << std::setw(6) << timing_sum_information[3];  // T3
0818   std::cout << "   " << std::setw(6) << timing_sum_information[4];  // T5
0819   std::cout << "   " << std::setw(6) << timing_sum_information[5];  // pLS
0820   std::cout << "   " << std::setw(6) << timing_sum_information[6];  // pT5
0821   std::cout << "   " << std::setw(6) << timing_sum_information[7];  // pT3
0822   std::cout << "   " << std::setw(6) << timing_sum_information[8];  // TC
0823   std::cout << "   " << std::setw(6) << timing_sum_information[9];  // Reset
0824   std::cout << "   " << std::setw(7) << timing_total_avg;           // Average total time
0825   std::cout << "   " << std::setw(7) << timing_totalshort_avg;      // Average total time
0826   std::cout << "+/- " << std::setw(4) << stdDev;
0827   std::cout << "   " << std::setw(7) << fullavg;  // Average full time
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   // std::cout << "Checking file path = " << fullpath << std::endl;
0845   // std::cout <<  " fullpath.string().c_str(): " << fullpath.string().c_str() <<  std::endl;
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   // Write out metadata of the code to the output_tfile
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   // Write gitversion info
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   // Write make script log
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   // Write git diff output in a separate string to gauge the difference
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   // Write Parse from makestr the TARGET
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   // Write out input sample or file name
0907   TNamed input("input", ana.input_raw_string.Data());
0908   input.Write();
0909 
0910   // Write the full command line used
0911   TNamed full_cmd_line("full_cmd_line", ana.full_cmd_line.Data());
0912   full_cmd_line.Write();
0913 
0914   // Write the TRACKLOOPERDIR
0915   TNamed tracklooper_path("tracklooper_path", ana.track_looper_dir_path.Data());
0916   tracklooper_path.Write();
0917 }