|
| 1 | +// MIT Licenced, Copyright (c) 2023-2025 CERN |
| 2 | +// |
| 3 | +// Code to run and time the jet finding of against various |
| 4 | +// HepMC3 input files |
| 5 | + |
| 6 | +// Original version of this code Philippe Gras, IRFU |
| 7 | +// Modified by Graeme A Stewart, CERN |
| 8 | + |
| 9 | +#include <cassert> |
| 10 | +#include <chrono> |
| 11 | +#include <cmath> |
| 12 | +#include <cstddef> |
| 13 | +#include <cstdio> // needed for io |
| 14 | +#include <iostream> // needed for io |
| 15 | +#include <stdlib.h> |
| 16 | +#include <string> |
| 17 | +#include <unistd.h> |
| 18 | +#include <vector> |
| 19 | + |
| 20 | +// Program Options Parser Library (https://github.com/badaix/popl) |
| 21 | +#include "popl.hpp" |
| 22 | + |
| 23 | +#include "JetReconstruction.h" |
| 24 | + |
| 25 | +#include "HepMC3/GenEvent.h" |
| 26 | +#include "HepMC3/GenParticle.h" |
| 27 | +#include "HepMC3/ReaderAscii.h" |
| 28 | + |
| 29 | +#include "cjetreconstruction-utils.hh" |
| 30 | +extern "C" { |
| 31 | +#include "julia_init.h" |
| 32 | +} |
| 33 | +using namespace std; |
| 34 | +using namespace popl; |
| 35 | +using Time = std::chrono::high_resolution_clock; |
| 36 | +using us = std::chrono::microseconds; |
| 37 | + |
| 38 | +// sorted into decreasing transverse momentum as in fastjet::sorted_by_pt |
| 39 | +void sorted_by_pt(jetreconstruction_JetsResult &jets) { |
| 40 | + std::sort(jets.data, jets.data + jets.length, |
| 41 | + [](const auto &left, const auto &right) { |
| 42 | + return left._pt2 > right._pt2; |
| 43 | + }); |
| 44 | +} |
| 45 | + |
| 46 | +jetreconstruction_ClusterSequence |
| 47 | +run_clustering(std::vector<jetreconstruction_PseudoJet> input_particles, |
| 48 | + jetreconstruction_RecoStrategy strategy, |
| 49 | + jetreconstruction_JetAlgorithm algorithm, double R, double p) { |
| 50 | + |
| 51 | + auto clust_seq = jetreconstruction_ClusterSequence{}; |
| 52 | + auto retv = jetreconstruction_jet_reconstruct( |
| 53 | + input_particles.data(), input_particles.size(), algorithm, R, strategy, |
| 54 | + &clust_seq); |
| 55 | + assert(retv == jetreconstruction_StatusCode::JETRECONSTRUCTION_STATUSCODE_OK); |
| 56 | + return clust_seq; |
| 57 | +} |
| 58 | + |
| 59 | +int main(int argc, char *argv[]) { |
| 60 | + init_julia(0, nullptr); |
| 61 | + // Default values |
| 62 | + int maxevents = -1; |
| 63 | + int skip_events = 0; |
| 64 | + int trials = 1; |
| 65 | + string mystrategy = "Best"; |
| 66 | + double power = -1.0; |
| 67 | + string alg = ""; |
| 68 | + double R = 0.4; |
| 69 | + string dump_file = ""; |
| 70 | + |
| 71 | + OptionParser opts("Allowed options"); |
| 72 | + auto help_option = opts.add<Switch>("h", "help", "produce help message"); |
| 73 | + auto max_events_option = opts.add<Value<int>>("m", "maxevents", "Maximum events in file to process (-1 = all events)", maxevents, &maxevents); |
| 74 | + auto skip_events_option = opts.add<Value<int>>("", "skipevents", "Number of events to skip over (0 = none)", skip_events, &skip_events); |
| 75 | + auto trials_option = opts.add<Value<int>>("n", "trials", "Number of repeated trials", trials, &trials); |
| 76 | + auto strategy_option = opts.add<Value<string>>("s", "strategy", "Valid values are 'Best' (default), 'N2Plain', 'N2Tiled'", mystrategy, &mystrategy); |
| 77 | + auto power_option = opts.add<Value<double>>("p", "power", "Algorithm p value: -1=antikt, 0=cambridge_aachen, 1=inclusive kt; otherwise generalised Kt", power, &power); |
| 78 | + auto alg_option = opts.add<Value<string>>("A", "algorithm", "Algorithm: AntiKt CA Kt GenKt EEKt Durham (overrides power)", alg, &alg); |
| 79 | + auto radius_option = opts.add<Value<double>>("R", "radius", "Algorithm R parameter", R, &R); |
| 80 | + auto ptmin_option = opts.add<Value<double>>("", "ptmin", "pt cut for inclusive jets"); |
| 81 | + auto dijmax_option = opts.add<Value<double>>("", "dijmax", "dijmax value for exclusive jets"); |
| 82 | + auto njets_option = opts.add<Value<int>>("", "njets", "njets value for exclusive jets"); |
| 83 | + auto dump_option = opts.add<Value<string>>("d", "dump", "Filename to dump jets to"); |
| 84 | + auto debug_clusterseq_option = opts.add<Switch>("c", "debug-clusterseq", "Dump cluster sequence jet and history content"); |
| 85 | + |
| 86 | + |
| 87 | + opts.parse(argc, argv); |
| 88 | + |
| 89 | + if (help_option->count() == 1) { |
| 90 | + cout << argv[0] << " [options] HEPMC3_INPUT_FILE" << endl; |
| 91 | + cout << endl; |
| 92 | + cout << opts << "\n"; |
| 93 | + cout << "Note the only one of ptmin, dijmax or njets can be specified!\n" << endl; |
| 94 | + exit(EXIT_SUCCESS); |
| 95 | + } |
| 96 | + |
| 97 | + const auto extra_args = opts.non_option_args(); |
| 98 | + std::string input_file{}; |
| 99 | + if (extra_args.size() == 1) { |
| 100 | + input_file = extra_args[0]; |
| 101 | + } else if (extra_args.size() == 0) { |
| 102 | + std::cerr << "No <HepMC3_input_file> argument after options" << std::endl; |
| 103 | + } else { |
| 104 | + std::cerr << "Only one <HepMC3_input_file> supported" << std::endl; |
| 105 | + } |
| 106 | + |
| 107 | + // Check we only have 1 option for final jet selection |
| 108 | + auto sum = int(njets_option->is_set()) + int(dijmax_option->is_set()) + int(ptmin_option->is_set()); |
| 109 | + if (sum != 1) { |
| 110 | + cerr << "One, and only one, of ptmin, dijmax or njets needs to be specified (currently " << |
| 111 | + sum << ")" << endl; |
| 112 | + exit(EXIT_FAILURE); |
| 113 | + } |
| 114 | + |
| 115 | + // read in input events |
| 116 | + //---------------------------------------------------------- |
| 117 | + auto events = read_input_events(input_file.c_str(), maxevents); |
| 118 | + |
| 119 | + // Set strategy |
| 120 | + auto strategy = |
| 121 | + jetreconstruction_RecoStrategy::JETRECONSTRUCTION_RECOSTRATEGY_BEST; |
| 122 | + if (mystrategy == string("N2Plain")) { |
| 123 | + strategy = |
| 124 | + jetreconstruction_RecoStrategy::JETRECONSTRUCTION_RECOSTRATEGY_N2PLAIN; |
| 125 | + } else if (mystrategy == string("N2Tiled")) { |
| 126 | + strategy = |
| 127 | + jetreconstruction_RecoStrategy::JETRECONSTRUCTION_RECOSTRATEGY_N2TILTED; |
| 128 | + } |
| 129 | + |
| 130 | + auto algorithm = |
| 131 | + jetreconstruction_JetAlgorithm::JETRECONSTRUCTION_JETALGORITHM_ANTIKT; |
| 132 | + if (alg != "") { |
| 133 | + if (alg == "AntiKt") { |
| 134 | + algorithm = |
| 135 | + jetreconstruction_JetAlgorithm::JETRECONSTRUCTION_JETALGORITHM_ANTIKT; |
| 136 | + power = -1.0; |
| 137 | + } else if (alg == "CA") { |
| 138 | + algorithm = |
| 139 | + jetreconstruction_JetAlgorithm::JETRECONSTRUCTION_JETALGORITHM_CA; |
| 140 | + power = 0.0; |
| 141 | + } else if (alg == "Kt") { |
| 142 | + algorithm = |
| 143 | + jetreconstruction_JetAlgorithm::JETRECONSTRUCTION_JETALGORITHM_KT; |
| 144 | + power = 1.0; |
| 145 | + //} else if (alg == "GenKt") { |
| 146 | + // algorithm = fastjet::genkt_algorithm; |
| 147 | + //} else if (alg == "Durham") { |
| 148 | + // algorithm = fastjet::ee_kt_algorithm; |
| 149 | + // power = 1.0; |
| 150 | + //} else if (alg == "EEKt") { |
| 151 | + // algorithm = fastjet::ee_genkt_algorithm; |
| 152 | + } else { |
| 153 | + std::cout << "Unknown algorithm type: " << alg << std::endl; |
| 154 | + exit(1); |
| 155 | + } |
| 156 | + } else { |
| 157 | + if (power == -1.0) { |
| 158 | + algorithm = |
| 159 | + jetreconstruction_JetAlgorithm::JETRECONSTRUCTION_JETALGORITHM_ANTIKT; |
| 160 | + } else if (power == 0.0) { |
| 161 | + algorithm = |
| 162 | + jetreconstruction_JetAlgorithm::JETRECONSTRUCTION_JETALGORITHM_CA; |
| 163 | + } else if (power == 1.0) { |
| 164 | + algorithm = |
| 165 | + jetreconstruction_JetAlgorithm::JETRECONSTRUCTION_JETALGORITHM_KT; |
| 166 | + //} else { |
| 167 | + // algorithm = fastjet::genkt_algorithm; |
| 168 | + } |
| 169 | + } |
| 170 | + |
| 171 | + std::cout << "Strategy: " << mystrategy << "; Power: " << power << "; Algorithm " << algorithm << std::endl; |
| 172 | + |
| 173 | + auto dump_fh = stdout; |
| 174 | + if (dump_option->is_set()) { |
| 175 | + if (dump_option->value() != "-") { |
| 176 | + dump_fh = fopen(dump_option->value().c_str(), "w"); |
| 177 | + } |
| 178 | + } |
| 179 | + |
| 180 | + double time_total = 0.0; |
| 181 | + double time_total2 = 0.0; |
| 182 | + double sigma = 0.0; |
| 183 | + double time_lowest = 1.0e20; |
| 184 | + for (long trial = 0; trial < trials; ++trial) { |
| 185 | + std::cout << "Trial " << trial << " "; |
| 186 | + auto start_t = std::chrono::steady_clock::now(); |
| 187 | + for (size_t ievt = skip_events_option->value(); ievt < events.size(); ++ievt) { |
| 188 | + auto cluster_sequence = |
| 189 | + run_clustering(events[ievt], strategy, algorithm, R, power); |
| 190 | + auto final_jets = jetreconstruction_JetsResult{nullptr,0}; |
| 191 | + if (ptmin_option->is_set()) { |
| 192 | + auto retv = jetreconstruction_inclusive_jets( |
| 193 | + &cluster_sequence, ptmin_option->value(), &final_jets); |
| 194 | + assert(retv == |
| 195 | + jetreconstruction_StatusCode::JETRECONSTRUCTION_STATUSCODE_OK); |
| 196 | + } else if (dijmax_option->is_set()) { |
| 197 | + auto retv = jetreconstruction_exclusive_jets_dcut( |
| 198 | + &cluster_sequence, dijmax_option->value(), &final_jets); |
| 199 | + assert(retv == |
| 200 | + jetreconstruction_StatusCode::JETRECONSTRUCTION_STATUSCODE_OK); |
| 201 | + } else if (njets_option->is_set()) { |
| 202 | + auto retv = jetreconstruction_exclusive_jets_njets( |
| 203 | + &cluster_sequence, njets_option->value(), &final_jets); |
| 204 | + assert(retv == |
| 205 | + jetreconstruction_StatusCode::JETRECONSTRUCTION_STATUSCODE_OK); |
| 206 | + } |
| 207 | + sorted_by_pt(final_jets); |
| 208 | + |
| 209 | + if (dump_option->is_set() && trial == 0) { |
| 210 | + fprintf(dump_fh, "Jets in processed event %zu\n", ievt + 1); |
| 211 | + |
| 212 | + // print out the details for each jet |
| 213 | + for (unsigned int i = 0; i < final_jets.length; i++) { |
| 214 | + const auto &jet = final_jets.data[i]; |
| 215 | + auto perp = std::sqrt(jet.px * jet.px + jet.py * jet.py); |
| 216 | + fprintf(dump_fh, "%5u %15.10f %15.10f %15.10f\n", i, jet._rap, |
| 217 | + jet._phi, perp); |
| 218 | + } |
| 219 | + |
| 220 | + // Dump the cluster sequence history content as well? |
| 221 | + if (debug_clusterseq_option->is_set()) { |
| 222 | + // dump_clusterseq(cluster_sequence); |
| 223 | + } |
| 224 | + jetreconstruction_ClusterSequence_free_members(&cluster_sequence); |
| 225 | + jetreconstruction_JetsResult_free_members(&final_jets); |
| 226 | + } |
| 227 | + } |
| 228 | + auto stop_t = std::chrono::steady_clock::now(); |
| 229 | + auto elapsed = stop_t - start_t; |
| 230 | + auto us_elapsed = double(chrono::duration_cast<chrono::microseconds>(elapsed).count()); |
| 231 | + std::cout << us_elapsed << " us" << endl; |
| 232 | + time_total += us_elapsed; |
| 233 | + time_total2 += us_elapsed*us_elapsed; |
| 234 | + if (us_elapsed < time_lowest) time_lowest = us_elapsed; |
| 235 | + } |
| 236 | + time_total /= trials; |
| 237 | + time_total2 /= trials; |
| 238 | + if (trials > 1) { |
| 239 | + sigma = std::sqrt(double(trials)/(trials-1) * (time_total2 - time_total*time_total)); |
| 240 | + } else { |
| 241 | + sigma = 0.0; |
| 242 | + } |
| 243 | + double mean_per_event = time_total / events.size(); |
| 244 | + double sigma_per_event = sigma / events.size(); |
| 245 | + time_lowest /= events.size(); |
| 246 | + std::cout << "Processed " << events.size() << " events, " << trials << " times" << endl; |
| 247 | + std::cout << "Total time " << time_total << " us" << endl; |
| 248 | + std::cout << "Time per event " << mean_per_event << " +- " << sigma_per_event << " us" << endl; |
| 249 | + std::cout << "Lowest time per event " << time_lowest << " us" << endl; |
| 250 | + shutdown_julia(0); |
| 251 | + return 0; |
| 252 | +} |
0 commit comments