From 7906c8dce6bfa7e6cf9774e3941d7dc3ff328138 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Fri, 31 May 2024 21:23:47 -0500 Subject: [PATCH 01/10] feat: add plugin usage info to `eicrecon -h` (#1479) ### Briefly, what does this PR introduce? This PR adds two extra lines to the usage info, explaining how to add extra plugins or disable default plugins. We now get: ``` $ eicrecon Usage: eicrecon [options] source1 source2 ... Description: Command-line interface for running JANA plugins. This can be used to read in events and process them. Command-line flags control configuration while additional arguments denote input files, which are to be loaded and processed by the appropriate EventSource plugin. Options: -h --help Display this message -v --version Display version information -j --janaversion Display JANA version information -c --configs Display configuration parameters -l --loadconfigs Load configuration parameters from file -d --dumpconfigs Dump configuration parameters to file -b --benchmark Run in benchmark mode -L --list-factories List all the factories without running -Pkey=value Specify a configuration parameter -Pplugin:param=value Specify a parameter value for a plugin --list-default-plugins List all the default plugins --list-available-plugins List plugins at $JANA_PLUGIN_PATH and $EICrecon_MY -Pplugins=values Comma-separated list of extra plugins to load -Pplugins_to_ignore=values Comma-separated list of plugins to ignore Example: eicrecon -Pplugins=plugin1,plugin2,plugin3 -Pnthreads=8 infile.root eicrecon -Ppodio:print_type_table=1 infile.root ``` ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [x] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [ ] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No. ### Does this PR change default behavior? No. --- src/utilities/eicrecon/eicrecon_cli.cpp | 935 ++++++++++++------------ 1 file changed, 474 insertions(+), 461 deletions(-) diff --git a/src/utilities/eicrecon/eicrecon_cli.cpp b/src/utilities/eicrecon/eicrecon_cli.cpp index bcbce93fdc..fd4f036222 100644 --- a/src/utilities/eicrecon/eicrecon_cli.cpp +++ b/src/utilities/eicrecon/eicrecon_cli.cpp @@ -8,7 +8,6 @@ #include #include #include -#include #include #include #include @@ -16,6 +15,7 @@ #include #include #include +#include #include #include @@ -29,508 +29,521 @@ #define STR(macro) QUOTE(macro) #ifndef EICRECON_APP_VERSION -# define EICRECON_APP_VERSION Error +#define EICRECON_APP_VERSION Error #endif #define EICRECON_APP_VERSION_STR STR(EICRECON_APP_VERSION) - namespace jana { - void PrintUsageOptions() { - std::cout << "Options:" << std::endl; - std::cout << " -h --help Display this message" << std::endl; - std::cout << " -v --version Display version information" << std::endl; - std::cout << " -j --janaversion Display JANA version information" << std::endl; - std::cout << " -c --configs Display configuration parameters" << std::endl; - std::cout << " -l --loadconfigs Load configuration parameters from file" << std::endl; - std::cout << " -d --dumpconfigs Dump configuration parameters to file" << std::endl; - std::cout << " -b --benchmark Run in benchmark mode" << std::endl; - std::cout << " -L --list-factories List all the factories without running" << std::endl; - std::cout << " -Pkey=value Specify a configuration parameter" << std::endl; - std::cout << " -Pplugin:param=value Specify a parameter value for a plugin" << std::endl; - std::cout << std::endl; - - std::cout << " --list-default-plugins List all the default plugins" << std::endl; - std::cout << " --list-available-plugins List plugins at $JANA_PLUGIN_PATH and $EICrecon_MY" << std::endl; - std::cout << std::endl << std::endl; - } - - void PrintUsageExample() { - - std::cout << "Example:" << std::endl; - std::cout << " eicrecon -Pplugins=plugin1,plugin2,plugin3 -Pnthreads=8 infile.root" << std::endl; - std::cout << " eicrecon -Ppodio:print_type_table=1 infile.root" << std::endl << std::endl; - std::cout << std::endl << std::endl; - } - - void PrintUsage() { - /// Prints jana.cc command-line options to stdout, for use by the CLI. - /// This does not include JANA parameters, which come from - /// JParameterManager::PrintParameters() instead. - - std::cout << std::endl; - std::cout << "Usage:" << std::endl; - std::cout << " eicrecon [options] source1 source2 ..." << std::endl; - std::cout << std::endl; - - std::cout << "Description:" << std::endl; - std::cout << " Command-line interface for running JANA plugins. This can be used to" << std::endl; - std::cout << " read in events and process them. Command-line flags control configuration" << std::endl; - std::cout << " while additional arguments denote input files, which are to be loaded and" << std::endl; - std::cout << " processed by the appropriate EventSource plugin." << std::endl; - std::cout << std::endl; - - PrintUsageOptions(); - PrintUsageExample(); - } - - void PrintVersion() { - std::cout << "EICrecon version: " << EICRECON_APP_VERSION_STR << std::endl; - } - - void PrintJANAVersion() { - std::cout << "JANA version: " << JVersion::GetVersion() << std::endl; - } - - void PrintDefaultPlugins(std::vector const& default_plugins) { - std::cout << "\n List default plugins:\n\n"; - printPluginNames(default_plugins); - std::cout << std::endl << std::endl; - } +void PrintUsageOptions() { + std::cout << "Options:" << std::endl; + std::cout << " -h --help Display this message" << std::endl; + std::cout << " -v --version Display version information" << std::endl; + std::cout << " -j --janaversion Display JANA version information" << std::endl; + std::cout << " -c --configs Display configuration parameters" << std::endl; + std::cout << " -l --loadconfigs Load configuration parameters from file" + << std::endl; + std::cout << " -d --dumpconfigs Dump configuration parameters to file" << std::endl; + std::cout << " -b --benchmark Run in benchmark mode" << std::endl; + std::cout << " -L --list-factories List all the factories without running" + << std::endl; + std::cout << " -Pkey=value Specify a configuration parameter" << std::endl; + std::cout << " -Pplugin:param=value Specify a parameter value for a plugin" + << std::endl; + std::cout << std::endl; + + std::cout << " --list-default-plugins List all the default plugins" << std::endl; + std::cout << " --list-available-plugins List plugins at $JANA_PLUGIN_PATH and $EICrecon_MY" + << std::endl; + std::cout << std::endl; + + std::cout << " -Pplugins=values Comma-separated list of extra plugins to load" + << std::endl; + std::cout << " -Pplugins_to_ignore=values Comma-separated list of plugins to ignore" + << std::endl; + + std::cout << std::endl << std::endl; +} - void GetPluginNamesInDir(std::set & plugin_names, std::string dir_str) { - // Edge case handler: taking care of invalid and empty dirs - if (std::filesystem::is_directory(dir_str) == false) - return; - if (std::filesystem::is_empty(dir_str)) - return; - - std::string full_path, filename; - for (const auto & entry : std::filesystem::directory_iterator(dir_str)) { - full_path = std::string(entry.path()); // Example: "/usr/local/plugins/Tutorial.so" - filename = full_path.substr(full_path.find_last_of("/") + 1); // Example: "Tutorial.so" - if (filename.substr(filename.size() - 3) == ".so") { - std::string s = filename.substr(0, filename.size() - 3); -// std::cout << filename << "==> " << s << std::endl; - plugin_names.insert(s); - } - } - } +void PrintUsageExample() { - /// Get the plugin names by searching for files named as *.so under $JANA_PLUGIN_PATH and $EICrecon_MY/plugins. - /// @note It does not guarantee any effectiveness of the plugins. - void GetPluginNamesFromEnvPath(std::set & plugin_names, const char* env_var) { - std::string dir_path, paths; - - const char* env_p = getenv(env_var); - if (env_p) { - if (strcmp(env_var, "EICrecon_MY") == 0) { - paths = std::string(env_p) + "/plugins"; - } - else { - paths = std::string(env_p); - } - - std::stringstream envvar_ss(paths); - while (getline(envvar_ss, dir_path, ':')) { - GetPluginNamesInDir(plugin_names, dir_path); - } - } - } + std::cout << "Example:" << std::endl; + std::cout << " eicrecon -Pplugins=plugin1,plugin2,plugin3 -Pnthreads=8 infile.root" + << std::endl; + std::cout << " eicrecon -Ppodio:print_type_table=1 infile.root" << std::endl << std::endl; + std::cout << std::endl << std::endl; +} - std::vector GetAvailablePluginNames(std::vector const& default_plugins) { - // Use set to remove duplicates. - /// @note The plugins will be override if you use the same plugin name at different paths. - std::set set_plugin_name; - for (std::string s : default_plugins) - set_plugin_name.insert(s); +void PrintUsage() { + /// Prints jana.cc command-line options to stdout, for use by the CLI. + /// This does not include JANA parameters, which come from + /// JParameterManager::PrintParameters() instead. + + std::cout << std::endl; + std::cout << "Usage:" << std::endl; + std::cout << " eicrecon [options] source1 source2 ..." << std::endl; + std::cout << std::endl; + + std::cout << "Description:" << std::endl; + std::cout << " Command-line interface for running JANA plugins. This can be used to" + << std::endl; + std::cout << " read in events and process them. Command-line flags control configuration" + << std::endl; + std::cout << " while additional arguments denote input files, which are to be loaded and" + << std::endl; + std::cout << " processed by the appropriate EventSource plugin." << std::endl; + std::cout << std::endl; + + PrintUsageOptions(); + PrintUsageExample(); +} - jana::GetPluginNamesFromEnvPath(set_plugin_name, "JANA_PLUGIN_PATH"); - jana::GetPluginNamesFromEnvPath(set_plugin_name, "EICrecon_MY"); +void PrintVersion() { std::cout << "EICrecon version: " << EICRECON_APP_VERSION_STR << std::endl; } - std::vector plugin_names(set_plugin_name.begin(), set_plugin_name.end()); - return plugin_names; - } +void PrintJANAVersion() { std::cout << "JANA version: " << JVersion::GetVersion() << std::endl; } - void PrintAvailablePlugins(std::vector const& default_plugins) { - std::cout << "\n List available plugins:\n\n"; - printPluginNames(GetAvailablePluginNames(default_plugins)); - std::cout << std::endl << std::endl; - } +void PrintDefaultPlugins(std::vector const& default_plugins) { + std::cout << "\n List default plugins:\n\n"; + printPluginNames(default_plugins); + std::cout << std::endl << std::endl; +} - bool HasPrintOnlyCliOptions(UserOptions& options, std::vector const& default_plugins) { - if (options.flags[jana::ShowUsage]) { - jana::PrintUsage(); - return true; - } - if (options.flags[jana::ShowVersion]) { - jana::PrintVersion(); - return true; - } - if (options.flags[jana::ShowJANAVersion]) { - jana::PrintJANAVersion(); - return true; - } - if (options.flags[jana::ShowDefaultPlugins]) { - jana::PrintDefaultPlugins(default_plugins); - return true; - } - if (options.flags[jana::ShowAvailablePlugins]) { - jana::PrintAvailablePlugins(default_plugins); - return true; - } - return false; +void GetPluginNamesInDir(std::set& plugin_names, std::string dir_str) { + // Edge case handler: taking care of invalid and empty dirs + if (std::filesystem::is_directory(dir_str) == false) + return; + if (std::filesystem::is_empty(dir_str)) + return; + + std::string full_path, filename; + for (const auto& entry : std::filesystem::directory_iterator(dir_str)) { + full_path = std::string(entry.path()); // Example: "/usr/local/plugins/Tutorial.so" + filename = full_path.substr(full_path.find_last_of("/") + 1); // Example: "Tutorial.so" + if (filename.substr(filename.size() - 3) == ".so") { + std::string s = filename.substr(0, filename.size() - 3); + // std::cout << filename << "==> " << s << std::endl; + plugin_names.insert(s); } + } +} - /// Detect whether the cli params @param options.params contain "-Pplugins_to_ignore=......". - /// If true, delete @param erase_str from the original cli string "-Pplugins_to_ignore". - bool HasExcludeDefaultPluginsInCliParams(UserOptions& options, const std::string erase_str) { - auto has_ignore_plugins = options.params.find("plugins_to_ignore"); - if (has_ignore_plugins == options.params.end()) - return false; - - // Has cli option "-Pplugins_to_ignore". Look for @param erase_str - size_t pos = has_ignore_plugins->second.find(erase_str); - if (pos == std::string::npos) // does not contain @param erase_str - return false; - - // Detect @param flag_str. Delete flag_str from the original cli option. - std::string ignore_str; - if (erase_str.length() + pos == has_ignore_plugins->second.length()) { // @param flag_str is at the end - ignore_str = has_ignore_plugins->second.erase(pos, erase_str.length()); - } else { // erase "," from "-Pplugins_to_ignore". - ignore_str = has_ignore_plugins->second.erase(pos, erase_str.length() + 1); - } - options.params["plugins_to_ignore"] = ignore_str; - return true; +/// Get the plugin names by searching for files named as *.so under $JANA_PLUGIN_PATH and +/// $EICrecon_MY/plugins. +/// @note It does not guarantee any effectiveness of the plugins. +void GetPluginNamesFromEnvPath(std::set& plugin_names, const char* env_var) { + std::string dir_path, paths; + + const char* env_p = getenv(env_var); + if (env_p) { + if (strcmp(env_var, "EICrecon_MY") == 0) { + paths = std::string(env_p) + "/plugins"; + } else { + paths = std::string(env_p); } - void AddAvailablePluginsToOptionParams(UserOptions& options, std::vector const& default_plugins) { - - std::set set_plugins; - // Add the plugins at $EICrecon_MY/plugins. -// jana::GetPluginNamesFromEnvPath(set_plugins, "EICrecon_MY"); // disabled as we do not want to automatically add these - - std::string plugins_str; // the complete plugins list - for (std::string s : set_plugins) - plugins_str += s + ","; - - // Add the default plugins into the plugin set if there is no - // "-Pplugins_to_ignore=default (exclude all default plugins)" option - if (HasExcludeDefaultPluginsInCliParams(options, "default") == false) - /// @note: The sequence of adding the default plugins matters. - /// Have to keep the original sequence to not causing troubles. - for (std::string s : default_plugins) { - plugins_str += s + ","; - } - - // Insert other plugins in cli option "-Pplugins=pl1,pl2,..." - auto has_cli_plugin_params = options.params.find("plugins"); - if (has_cli_plugin_params != options.params.end()) { - plugins_str += has_cli_plugin_params->second; - options.params["plugins"] = plugins_str; - } else { - options.params["plugins"] = plugins_str.substr(0, plugins_str.size() - 1); // exclude last "," - } + std::stringstream envvar_ss(paths); + while (getline(envvar_ss, dir_path, ':')) { + GetPluginNamesInDir(plugin_names, dir_path); } + } +} - JApplication* CreateJApplication(UserOptions& options) { - - auto *para_mgr = new JParameterManager(); // JApplication owns params_copy, does not own eventSources +std::vector GetAvailablePluginNames(std::vector const& default_plugins) { + // Use set to remove duplicates. + /// @note The plugins will be override if you use the same plugin name at different paths. + std::set set_plugin_name; + for (std::string s : default_plugins) + set_plugin_name.insert(s); - // Add the cli options based on the user inputs - for (auto pair : options.params) { - para_mgr->SetParameter(pair.first, pair.second); - } + jana::GetPluginNamesFromEnvPath(set_plugin_name, "JANA_PLUGIN_PATH"); + jana::GetPluginNamesFromEnvPath(set_plugin_name, "EICrecon_MY"); - // Shut down the [INFO] msg of adding plugins, printing cpu info - if (options.flags[ListFactories]) { - para_mgr->SetParameter( - "log:off", - "JPluginLoader,JArrowProcessingController,JArrow" - ); - } - - if (options.flags[LoadConfigs]) { - // If the user specified an external config file, we should definitely use that - try { - para_mgr->ReadConfigFile(options.load_config_file); - } - catch (JException &e) { - std::cout << "Problem loading config file '" << options.load_config_file << "'. Exiting." << std::endl - << std::endl; - exit(-1); - } - std::cout << "Loaded config file '" << options.load_config_file << "'." << std::endl << std::endl; - } + std::vector plugin_names(set_plugin_name.begin(), set_plugin_name.end()); + return plugin_names; +} - // If the user hasn't specified a timeout (on cmd line or in config file), set the timeout to something reasonably high - if (para_mgr->FindParameter("jana:timeout") == nullptr) { - para_mgr->SetParameter("jana:timeout", 180); // seconds - para_mgr->SetParameter("jana:warmup_timeout", 180); // seconds - } +void PrintAvailablePlugins(std::vector const& default_plugins) { + std::cout << "\n List available plugins:\n\n"; + printPluginNames(GetAvailablePluginNames(default_plugins)); + std::cout << std::endl << std::endl; +} - auto *app = new JApplication(para_mgr); +bool HasPrintOnlyCliOptions(UserOptions& options, std::vector const& default_plugins) { + if (options.flags[jana::ShowUsage]) { + jana::PrintUsage(); + return true; + } + if (options.flags[jana::ShowVersion]) { + jana::PrintVersion(); + return true; + } + if (options.flags[jana::ShowJANAVersion]) { + jana::PrintJANAVersion(); + return true; + } + if (options.flags[jana::ShowDefaultPlugins]) { + jana::PrintDefaultPlugins(default_plugins); + return true; + } + if (options.flags[jana::ShowAvailablePlugins]) { + jana::PrintAvailablePlugins(default_plugins); + return true; + } + return false; +} - const char* env_p = getenv("EICrecon_MY"); - if( env_p ){ - app->AddPluginPath( std::string(env_p) + "/plugins" ); - } +/// Detect whether the cli params @param options.params contain +/// "-Pplugins_to_ignore=......". If true, delete @param erase_str from the original cli +/// string "-Pplugins_to_ignore". +bool HasExcludeDefaultPluginsInCliParams(UserOptions& options, const std::string erase_str) { + auto has_ignore_plugins = options.params.find("plugins_to_ignore"); + if (has_ignore_plugins == options.params.end()) + return false; + + // Has cli option "-Pplugins_to_ignore". Look for @param erase_str + size_t pos = has_ignore_plugins->second.find(erase_str); + if (pos == std::string::npos) // does not contain @param erase_str + return false; + + // Detect @param flag_str. Delete flag_str from the original cli option. + std::string ignore_str; + if (erase_str.length() + pos == + has_ignore_plugins->second.length()) { // @param flag_str is at the end + ignore_str = has_ignore_plugins->second.erase(pos, erase_str.length()); + } else { // erase "," from "-Pplugins_to_ignore". + ignore_str = has_ignore_plugins->second.erase(pos, erase_str.length() + 1); + } + options.params["plugins_to_ignore"] = ignore_str; + return true; +} - for (auto event_src : options.eventSources) { - app->Add(event_src); - } - return app; +void AddAvailablePluginsToOptionParams(UserOptions& options, + std::vector const& default_plugins) { + + std::set set_plugins; + // Add the plugins at $EICrecon_MY/plugins. + // jana::GetPluginNamesFromEnvPath(set_plugins, "EICrecon_MY"); // disabled as we do not + // want to automatically add these + + std::string plugins_str; // the complete plugins list + for (std::string s : set_plugins) + plugins_str += s + ","; + + // Add the default plugins into the plugin set if there is no + // "-Pplugins_to_ignore=default (exclude all default plugins)" option + if (HasExcludeDefaultPluginsInCliParams(options, "default") == false) + /// @note: The sequence of adding the default plugins matters. + /// Have to keep the original sequence to not causing troubles. + for (std::string s : default_plugins) { + plugins_str += s + ","; } - void AddDefaultPluginsToJApplication(JApplication* app, std::vector const& default_plugins) { - for (std::string s : default_plugins) - app->AddPlugin(s); - } + // Insert other plugins in cli option "-Pplugins=pl1,pl2,..." + auto has_cli_plugin_params = options.params.find("plugins"); + if (has_cli_plugin_params != options.params.end()) { + plugins_str += has_cli_plugin_params->second; + options.params["plugins"] = plugins_str; + } else { + options.params["plugins"] = plugins_str.substr(0, plugins_str.size() - 1); // exclude last "," + } +} - void PrintFactories(JApplication* app) { - std::cout << std::endl << "List all the factories:" << std::endl << std::endl; - printFactoryTable(app->GetComponentSummary()); - std::cout << std::endl; +JApplication* CreateJApplication(UserOptions& options) { + + auto* para_mgr = + new JParameterManager(); // JApplication owns params_copy, does not own eventSources + + // Add the cli options based on the user inputs + for (auto pair : options.params) { + para_mgr->SetParameter(pair.first, pair.second); + } + + // Shut down the [INFO] msg of adding plugins, printing cpu info + if (options.flags[ListFactories]) { + para_mgr->SetParameter("log:off", "JPluginLoader,JArrowProcessingController,JArrow"); + } + + if (options.flags[LoadConfigs]) { + // If the user specified an external config file, we should definitely use that + try { + para_mgr->ReadConfigFile(options.load_config_file); + } catch (JException& e) { + std::cout << "Problem loading config file '" << options.load_config_file << "'. Exiting." + << std::endl + << std::endl; + exit(-1); } + std::cout << "Loaded config file '" << options.load_config_file << "'." << std::endl + << std::endl; + } + + // If the user hasn't specified a timeout (on cmd line or in config file), set the timeout to + // something reasonably high + if (para_mgr->FindParameter("jana:timeout") == nullptr) { + para_mgr->SetParameter("jana:timeout", 180); // seconds + para_mgr->SetParameter("jana:warmup_timeout", 180); // seconds + } + + auto* app = new JApplication(para_mgr); + + const char* env_p = getenv("EICrecon_MY"); + if (env_p) { + app->AddPluginPath(std::string(env_p) + "/plugins"); + } + + for (auto event_src : options.eventSources) { + app->Add(event_src); + } + return app; +} - void PrintPodioCollections(JApplication* app) { - if (app->GetJParameterManager()->Exists("PODIO:PRINT_TYPE_TABLE")) { - bool print_type_table = app->GetParameterValue("podio:print_type_table"); - - // cli criteria: Ppodio:print_type_table=1 - if (print_type_table) { - auto event_sources = app->GetService()->get_evt_srces(); - for (auto *event_source : event_sources) { -// std::cout << event_source->GetPluginName() << std::endl; // podio.so -// std::cout << event_source->GetResourceName() << std::endl; - if (event_source->GetPluginName().find("podio") != std::string::npos) - event_source->DoInitialize(); - } - } - - } - } +void AddDefaultPluginsToJApplication(JApplication* app, + std::vector const& default_plugins) { + for (std::string s : default_plugins) + app->AddPlugin(s); +} - void PrintConfigParameters(JApplication* app){ - /// Print a table of the currently defined configuration parameters. - /// n.b. this mostly duplicates a call to app->GetJParameterManager()->PrintParameters() - /// but avoids the issue it has of setting the values column to same - /// width for all parameters. (That leads to lots of whitespace being - /// printed due to the very long podio:output_include_collections param. - - // Determine column widths - auto params = app->GetJParameterManager()->GetAllParameters(); - size_t max_key_length = 0; - size_t max_val_length = 0; - size_t max_max_val_length = 32; // maximum width allowed for column. - for( auto &[key, p] : params ){ - if( key.length() > max_key_length ) max_key_length = key.length(); - if( p->GetValue().length() > max_val_length ){ - if( p->GetValue().length() <= max_max_val_length ) max_val_length = p->GetValue().length(); - } - } +void PrintFactories(JApplication* app) { + std::cout << std::endl << "List all the factories:" << std::endl << std::endl; + printFactoryTable(app->GetComponentSummary()); + std::cout << std::endl; +} - std::cout << "\nConfiguration Parameters:" << std::endl; - std::cout << "Name" + std::string(std::max(max_key_length, size_t(4)) - 4, ' ') << " : "; - std::cout << "Value" + std::string(std::max(max_val_length, size_t(5)) - 5, ' ') << " : "; - std::cout << "Description" << std::endl; - std::cout << std::string(max_key_length+max_val_length+20, '-') << std::endl; - for( auto &[key, p] : params ){ - std::stringstream ss; - int key_length_diff = max_key_length - key.length(); - if( key_length_diff>0 ) ss << std::string(key_length_diff, ' '); - ss << key; - ss << " | "; - - int val_length_diff = max_val_length - p->GetValue().length(); - if( val_length_diff>0 ) ss << std::string(val_length_diff, ' '); - ss << p->GetValue(); - ss << " | "; - ss << p->GetDescription(); - - std::cout << ss.str() << std::endl; - } - std::cout << std::string(max_key_length+max_val_length+20, '-') << std::endl; +void PrintPodioCollections(JApplication* app) { + if (app->GetJParameterManager()->Exists("PODIO:PRINT_TYPE_TABLE")) { + bool print_type_table = app->GetParameterValue("podio:print_type_table"); + + // cli criteria: Ppodio:print_type_table=1 + if (print_type_table) { + auto event_sources = app->GetService()->get_evt_srces(); + for (auto* event_source : event_sources) { + // std::cout << event_source->GetPluginName() << std::endl; // podio.so + // std::cout << event_source->GetResourceName() << std::endl; + if (event_source->GetPluginName().find("podio") != std::string::npos) + event_source->DoInitialize(); + } } + } +} - int Execute(JApplication* app, UserOptions &options) { +void PrintConfigParameters(JApplication* app) { + /// Print a table of the currently defined configuration parameters. + /// n.b. this mostly duplicates a call to app->GetJParameterManager()->PrintParameters() + /// but avoids the issue it has of setting the values column to same + /// width for all parameters. (That leads to lots of whitespace being + /// printed due to the very long podio:output_include_collections param. + + // Determine column widths + auto params = app->GetJParameterManager()->GetAllParameters(); + size_t max_key_length = 0; + size_t max_val_length = 0; + size_t max_max_val_length = 32; // maximum width allowed for column. + for (auto& [key, p] : params) { + if (key.length() > max_key_length) + max_key_length = key.length(); + if (p->GetValue().length() > max_val_length) { + if (p->GetValue().length() <= max_max_val_length) + max_val_length = p->GetValue().length(); + } + } + + std::cout << "\nConfiguration Parameters:" << std::endl; + std::cout << "Name" + std::string(std::max(max_key_length, size_t(4)) - 4, ' ') << " : "; + std::cout << "Value" + std::string(std::max(max_val_length, size_t(5)) - 5, ' ') << " : "; + std::cout << "Description" << std::endl; + std::cout << std::string(max_key_length + max_val_length + 20, '-') << std::endl; + for (auto& [key, p] : params) { + std::stringstream ss; + int key_length_diff = max_key_length - key.length(); + if (key_length_diff > 0) + ss << std::string(key_length_diff, ' '); + ss << key; + ss << " | "; + + int val_length_diff = max_val_length - p->GetValue().length(); + if (val_length_diff > 0) + ss << std::string(val_length_diff, ' '); + ss << p->GetValue(); + ss << " | "; + ss << p->GetDescription(); + + std::cout << ss.str() << std::endl; + } + std::cout << std::string(max_key_length + max_val_length + 20, '-') << std::endl; +} - std::cout << std::endl; +int Execute(JApplication* app, UserOptions& options) { - // std::cout << "JANA " << JVersion::GetVersion() << " [" << JVersion::GetRevision() << "]" << std::endl; + std::cout << std::endl; - if (options.flags[ShowConfigs]) { - // Load all plugins, collect all parameters, exit without running anything - app->Initialize(); - if (options.flags[Benchmark]) { - JBenchmarker benchmarker(app); // Show benchmarking configs only if benchmarking mode specified - } -// app->GetJParameterManager()->PrintParameters(true); - PrintConfigParameters(app); - } - else if (options.flags[DumpConfigs]) { - // Load all plugins, dump parameters to file, exit without running anything - app->Initialize(); - std::cout << std::endl << "Writing configuration options to file: " << options.dump_config_file - << std::endl; - app->GetJParameterManager()->WriteConfigFile(options.dump_config_file); - } - else if (options.flags[Benchmark]) { - JSignalHandler::register_handlers(app); - // Run JANA in benchmark mode - JBenchmarker benchmarker(app); // Benchmarking params override default params - benchmarker.RunUntilFinished(); // Benchmarker will control JApp Run/Stop - } - else if (options.flags[ListFactories]) { - app->Initialize(); - PrintFactories(app); + // std::cout << "JANA " << JVersion::GetVersion() << " [" << JVersion::GetRevision() << "]" << + // std::endl; - // TODO: more elegant processing here - PrintPodioCollections(app); - } - else { - // Run JANA in normal mode - try { - printJANAHeaderIMG(); - JSignalHandler::register_handlers(app); - app->Run(); - } - catch (JException& e) { - std::cout << "----------------------------------------------------------" << std::endl; - std::cout << e << std::endl; - app->SetExitCode(EXIT_FAILURE); - } - catch (std::runtime_error& e) { - std::cout << "----------------------------------------------------------" << std::endl; - std::cout << "Exception: " << e.what() << std::endl; - app->SetExitCode(EXIT_FAILURE); - } - } - return (int) app->GetExitCode(); + if (options.flags[ShowConfigs]) { + // Load all plugins, collect all parameters, exit without running anything + app->Initialize(); + if (options.flags[Benchmark]) { + JBenchmarker benchmarker( + app); // Show benchmarking configs only if benchmarking mode specified + } + // app->GetJParameterManager()->PrintParameters(true); + PrintConfigParameters(app); + } else if (options.flags[DumpConfigs]) { + // Load all plugins, dump parameters to file, exit without running anything + app->Initialize(); + std::cout << std::endl + << "Writing configuration options to file: " << options.dump_config_file << std::endl; + app->GetJParameterManager()->WriteConfigFile(options.dump_config_file); + } else if (options.flags[Benchmark]) { + JSignalHandler::register_handlers(app); + // Run JANA in benchmark mode + JBenchmarker benchmarker(app); // Benchmarking params override default params + benchmarker.RunUntilFinished(); // Benchmarker will control JApp Run/Stop + } else if (options.flags[ListFactories]) { + app->Initialize(); + PrintFactories(app); + + // TODO: more elegant processing here + PrintPodioCollections(app); + } else { + // Run JANA in normal mode + try { + printJANAHeaderIMG(); + JSignalHandler::register_handlers(app); + app->Run(); + } catch (JException& e) { + std::cout << "----------------------------------------------------------" << std::endl; + std::cout << e << std::endl; + app->SetExitCode(EXIT_FAILURE); + } catch (std::runtime_error& e) { + std::cout << "----------------------------------------------------------" << std::endl; + std::cout << "Exception: " << e.what() << std::endl; + app->SetExitCode(EXIT_FAILURE); } + } + return (int)app->GetExitCode(); +} +UserOptions GetCliOptions(int nargs, char* argv[], bool expect_extra) { + + UserOptions options; + + std::map tokenizer; + tokenizer["-h"] = ShowUsage; + tokenizer["--help"] = ShowUsage; + tokenizer["-v"] = ShowVersion; + tokenizer["--version"] = ShowVersion; + tokenizer["-j"] = ShowJANAVersion; + tokenizer["--janaversion"] = ShowJANAVersion; + tokenizer["-c"] = ShowConfigs; + tokenizer["--configs"] = ShowConfigs; + tokenizer["-l"] = LoadConfigs; + tokenizer["--loadconfigs"] = LoadConfigs; + tokenizer["-d"] = DumpConfigs; + tokenizer["--dumpconfigs"] = DumpConfigs; + tokenizer["-b"] = Benchmark; + tokenizer["--benchmark"] = Benchmark; + tokenizer["-L"] = ListFactories; + tokenizer["--list-factories"] = ListFactories; + tokenizer["--list-default-plugins"] = ShowDefaultPlugins; + tokenizer["--list-available-plugins"] = ShowAvailablePlugins; + + // `eicrecon` has the same effect with `eicrecon -h` + if (nargs == 1) { + options.flags[ShowUsage] = true; + } + + for (int i = 1; i < nargs; i++) { + + std::string arg = argv[i]; + // std::cout << "Found arg " << arg << std::endl; + + if (argv[i][0] != '-') { + options.eventSources.push_back(arg); + continue; + } - UserOptions GetCliOptions(int nargs, char *argv[], bool expect_extra) { - - UserOptions options; - - std::map tokenizer; - tokenizer["-h"] = ShowUsage; - tokenizer["--help"] = ShowUsage; - tokenizer["-v"] = ShowVersion; - tokenizer["--version"] = ShowVersion; - tokenizer["-j"] = ShowJANAVersion; - tokenizer["--janaversion"] = ShowJANAVersion; - tokenizer["-c"] = ShowConfigs; - tokenizer["--configs"] = ShowConfigs; - tokenizer["-l"] = LoadConfigs; - tokenizer["--loadconfigs"] = LoadConfigs; - tokenizer["-d"] = DumpConfigs; - tokenizer["--dumpconfigs"] = DumpConfigs; - tokenizer["-b"] = Benchmark; - tokenizer["--benchmark"] = Benchmark; - tokenizer["-L"] = ListFactories; - tokenizer["--list-factories"] = ListFactories; - tokenizer["--list-default-plugins"] = ShowDefaultPlugins; - tokenizer["--list-available-plugins"] = ShowAvailablePlugins; - - // `eicrecon` has the same effect with `eicrecon -h` - if (nargs == 1) { - options.flags[ShowUsage] = true; + switch (tokenizer[arg]) { + + case Benchmark: + options.flags[Benchmark] = true; + break; + + case ShowUsage: + options.flags[ShowUsage] = true; + break; + + case ShowVersion: + options.flags[ShowVersion] = true; + break; + + case ShowJANAVersion: + options.flags[ShowJANAVersion] = true; + break; + + case ShowConfigs: + options.flags[ShowConfigs] = true; + break; + + case LoadConfigs: + options.flags[LoadConfigs] = true; + if (i + 1 < nargs && argv[i + 1][0] != '-') { + options.load_config_file = argv[i + 1]; + i += 1; + } else { + options.load_config_file = "jana.config"; + } + break; + + case DumpConfigs: + options.flags[DumpConfigs] = true; + if (i + 1 < nargs && argv[i + 1][0] != '-') { + options.dump_config_file = argv[i + 1]; + i += 1; + } else { + options.dump_config_file = "jana.config"; + } + break; + + case ListFactories: + options.flags[ListFactories] = true; + break; + + case ShowDefaultPlugins: + options.flags[ShowDefaultPlugins] = true; + break; + + case ShowAvailablePlugins: + options.flags[ShowAvailablePlugins] = true; + break; + + // TODO: add exclude plugin options + case Unknown: + if (argv[i][0] == '-' && argv[i][1] == 'P') { + + size_t pos = arg.find("="); + if ((pos != std::string::npos) && (pos > 2)) { + std::string key = arg.substr(2, pos - 2); + std::string val = arg.substr(pos + 1); + if (options.params.find(key) != options.params.end()) { + std::cout << "Duplicate parameter '" << arg << "' ignored" << std::endl; + } else { + options.params.insert({key, val}); + } + } else { + std::cout << "Invalid JANA parameter '" << arg << "': Expected format -Pkey=value" + << std::endl; + options.flags[ShowConfigs] = true; } - - for (int i = 1; i < nargs; i++) { - - std::string arg = argv[i]; - // std::cout << "Found arg " << arg << std::endl; - - if (argv[i][0] != '-') { - options.eventSources.push_back(arg); - continue; - } - - switch (tokenizer[arg]) { - - case Benchmark: - options.flags[Benchmark] = true; - break; - - case ShowUsage: - options.flags[ShowUsage] = true; - break; - - case ShowVersion: - options.flags[ShowVersion] = true; - break; - - case ShowJANAVersion: - options.flags[ShowJANAVersion] = true; - break; - - case ShowConfigs: - options.flags[ShowConfigs] = true; - break; - - case LoadConfigs: - options.flags[LoadConfigs] = true; - if (i + 1 < nargs && argv[i + 1][0] != '-') { - options.load_config_file = argv[i + 1]; - i += 1; - } else { - options.load_config_file = "jana.config"; - } - break; - - case DumpConfigs: - options.flags[DumpConfigs] = true; - if (i + 1 < nargs && argv[i + 1][0] != '-') { - options.dump_config_file = argv[i + 1]; - i += 1; - } else { - options.dump_config_file = "jana.config"; - } - break; - - case ListFactories: - options.flags[ListFactories] = true; - break; - - case ShowDefaultPlugins: - options.flags[ShowDefaultPlugins] = true; - break; - - case ShowAvailablePlugins: - options.flags[ShowAvailablePlugins] = true; - break; - - // TODO: add exclude plugin options - case Unknown: - if (argv[i][0] == '-' && argv[i][1] == 'P') { - - size_t pos = arg.find("="); - if ((pos != std::string::npos) && (pos > 2)) { - std::string key = arg.substr(2, pos - 2); - std::string val = arg.substr(pos + 1); - if (options.params.find(key) != options.params.end()) { - std::cout << "Duplicate parameter '" << arg << "' ignored" << std::endl; - } else { - options.params.insert({key, val}); - } - } else { - std::cout << "Invalid JANA parameter '" << arg - << "': Expected format -Pkey=value" << std::endl; - options.flags[ShowConfigs] = true; - } - } else { - if (!expect_extra) { - std::cout << "Invalid command line flag '" << arg << "'" << std::endl; - options.flags[ShowUsage] = true; - } - } - } + } else { + if (!expect_extra) { + std::cout << "Invalid command line flag '" << arg << "'" << std::endl; + options.flags[ShowUsage] = true; } - return options; + } } + } + return options; } +} // namespace jana From cef88bf3e5ebf48a291db3f86ce1c263c99654d4 Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Sun, 2 Jun 2024 02:01:18 -0400 Subject: [PATCH 02/10] updated forward hcal insert recon to use hexplit and topo clustering (#1474) ### Briefly, what does this PR introduce? Use HEXPLIT and topo-clustering to the forward hadronic Calorimeter Insert (CALI). ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [X] New feature (issue #1480) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [X] Tests for the changes have been added - [X] Documentation has been added / updated - [X] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No ### Does this PR change default behavior? Yes. Changes which algorithms are used in Calorimeter Insert reconstruction --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Dmitry Kalinkin --- src/detectors/FHCAL/FHCAL.cc | 40 +++++++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/src/detectors/FHCAL/FHCAL.cc b/src/detectors/FHCAL/FHCAL.cc index 176fed0d7f..e58e146f5d 100644 --- a/src/detectors/FHCAL/FHCAL.cc +++ b/src/detectors/FHCAL/FHCAL.cc @@ -4,6 +4,7 @@ #include #include #include +#include #include #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h" @@ -14,6 +15,8 @@ #include "factories/calorimetry/CalorimeterHitsMerger_factory.h" #include "factories/calorimetry/CalorimeterIslandCluster_factory.h" #include "factories/calorimetry/CalorimeterTruthClustering_factory.h" +#include "factories/calorimetry/HEXPLIT_factory.h" +#include "factories/calorimetry/ImagingTopoCluster_factory.h" extern "C" { void InitPlugin(JApplication *app) { @@ -57,6 +60,7 @@ extern "C" { .sampFrac = "0.0098", .readout = "HcalEndcapPInsertHits", + .layerField="layer", }, app // TODO: Remove me once fixed )); @@ -73,20 +77,32 @@ extern "C" { "HcalEndcapPInsertTruthProtoClusters", {"HcalEndcapPInsertMergedHits", "HcalEndcapPInsertHits"}, {"HcalEndcapPInsertTruthProtoClusters"}, app // TODO: Remove me once fixed )); - app->Add(new JOmniFactoryGeneratorT( - "HcalEndcapPInsertIslandProtoClusters", {"HcalEndcapPInsertMergedHits"}, {"HcalEndcapPInsertIslandProtoClusters"}, + + app->Add(new JOmniFactoryGeneratorT( + "HcalEndcapPInsertSubcellHits", {"HcalEndcapPInsertRecHits"}, {"HcalEndcapPInsertSubcellHits"}, + { + .MIP = 800. * dd4hep::keV, + .Emin_in_MIPs=0.1, + .tmax=50 * dd4hep::ns, + }, + app // TODO: Remove me once fixed + )); + + double side_length=18.89 * dd4hep::mm; + app->Add(new JOmniFactoryGeneratorT( + "HcalEndcapPInsertImagingProtoClusters", {"HcalEndcapPInsertSubcellHits"}, {"HcalEndcapPInsertImagingProtoClusters"}, { - .sectorDist = 5.0 * dd4hep::cm, - .localDistXY = {15*dd4hep::mm, 15*dd4hep::mm}, - .dimScaledLocalDistXY = {15.0*dd4hep::mm, 15.0*dd4hep::mm}, - .splitCluster = true, - .minClusterHitEdep = 0.0 * dd4hep::MeV, - .minClusterCenterEdep = 30.0 * dd4hep::MeV, - .transverseEnergyProfileMetric = "globalDistEtaPhi", - .transverseEnergyProfileScale = 1., + .neighbourLayersRange = 1, + .localDistXY = {0.76*side_length, 0.76*side_length*sin(M_PI/3)}, + .layerDistEtaPhi = {17e-3, 18.1e-3}, + .sectorDist = 10.0 * dd4hep::cm, + .minClusterHitEdep = 100.0 * dd4hep::keV, + .minClusterCenterEdep = 11.0 * dd4hep::MeV, + .minClusterEdep = 11.0 * dd4hep::MeV, + .minClusterNhits = 10, }, app // TODO: Remove me once fixed - )); + )); app->Add( new JOmniFactoryGeneratorT( @@ -108,7 +124,7 @@ extern "C" { app->Add( new JOmniFactoryGeneratorT( "HcalEndcapPInsertClusters", - {"HcalEndcapPInsertIslandProtoClusters", // edm4eic::ProtoClusterCollection + {"HcalEndcapPInsertImagingProtoClusters", // edm4eic::ProtoClusterCollection "HcalEndcapPInsertHits"}, // edm4hep::SimCalorimeterHitCollection {"HcalEndcapPInsertClusters", // edm4eic::Cluster "HcalEndcapPInsertClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation From 018fa7700d9a3fd4f795833fcf3dc6d2decc8463 Mon Sep 17 00:00:00 2001 From: Simon Gardner Date: Sun, 2 Jun 2024 23:56:30 +0100 Subject: [PATCH 03/10] Add low Q2 NN inference (#1467) ### Briefly, what does this PR introduce? Adds the inference step, reconstructing the initial momentum of the low-Q2 electrons from their tracks. This relies on https://github.com/eic/epic-data/pull/17 to add the inference model to the expected default location (so far tested by manually copying the file there). This will be followed by the final PR which finally injects the low-Q2 electron tracks into the main tracking data flow so they get treated equally by the following algorithms. ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [x] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [ ] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No ### Does this PR change default behavior? No --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Dmitry Kalinkin Co-authored-by: Christopher Dilks Co-authored-by: Wouter Deconinck --- src/algorithms/fardetectors/CMakeLists.txt | 4 + .../FarDetectorMLReconstruction.cc | 130 ++++++++++++++++++ .../FarDetectorMLReconstruction.h | 67 +++++++++ .../FarDetectorMLReconstructionConfig.h | 15 ++ src/detectors/LOWQ2/LOWQ2.cc | 28 +++- .../FarDetectorMLReconstruction_factory.h | 55 ++++++++ src/services/io/podio/JEventProcessorPODIO.cc | 7 +- 7 files changed, 298 insertions(+), 8 deletions(-) create mode 100644 src/algorithms/fardetectors/FarDetectorMLReconstruction.cc create mode 100644 src/algorithms/fardetectors/FarDetectorMLReconstruction.h create mode 100644 src/algorithms/fardetectors/FarDetectorMLReconstructionConfig.h create mode 100644 src/factories/fardetectors/FarDetectorMLReconstruction_factory.h diff --git a/src/algorithms/fardetectors/CMakeLists.txt b/src/algorithms/fardetectors/CMakeLists.txt index 1596a83572..33ee97a05a 100644 --- a/src/algorithms/fardetectors/CMakeLists.txt +++ b/src/algorithms/fardetectors/CMakeLists.txt @@ -11,4 +11,8 @@ plugin_glob_all(${PLUGIN_NAME}) # Find dependencies plugin_add_dd4hep(${PLUGIN_NAME}) +plugin_add_eigen3(${PLUGIN_NAME}) plugin_add_event_model(${PLUGIN_NAME}) +plugin_add_cern_root(${PLUGIN_NAME}) + +plugin_link_libraries(${PLUGIN_NAME} ROOT::TMVA) diff --git a/src/algorithms/fardetectors/FarDetectorMLReconstruction.cc b/src/algorithms/fardetectors/FarDetectorMLReconstruction.cc new file mode 100644 index 0000000000..d92b382714 --- /dev/null +++ b/src/algorithms/fardetectors/FarDetectorMLReconstruction.cc @@ -0,0 +1,130 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Simon Gardner + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "FarDetectorMLReconstruction.h" +#include "algorithms/fardetectors/FarDetectorMLReconstructionConfig.h" + +namespace eicrecon { + + + void FarDetectorMLReconstruction::init() { + + m_reader = new TMVA::Reader( "!Color:!Silent" ); + // Create a set of variables and declare them to the reader + // - the variable names MUST corresponds in name and type to those given in the weight file(s) used + m_reader->AddVariable( "LowQ2Tracks[0].loc.a", &nnInput[FarDetectorMLNNIndexIn::PosY] ); + m_reader->AddVariable( "LowQ2Tracks[0].loc.b", &nnInput[FarDetectorMLNNIndexIn::PosZ] ); + m_reader->AddVariable( "sin(LowQ2Tracks[0].phi)*sin(LowQ2Tracks[0].theta)", &nnInput[FarDetectorMLNNIndexIn::DirX] ); + m_reader->AddVariable( "cos(LowQ2Tracks[0].phi)*sin(LowQ2Tracks[0].theta)", &nnInput[FarDetectorMLNNIndexIn::DirY] ); + + // Locate and load the weight file + // TODO - Add functionality to select passed by configuration + bool methodFound = false; + if(!m_cfg.modelPath.empty()){ + try{ + m_method = dynamic_cast(m_reader->BookMVA( m_cfg.methodName, m_cfg.modelPath )); + } + catch(std::exception &e){ + error(fmt::format("Failed to load method {} from file {}: {}", m_cfg.methodName, m_cfg.modelPath, e.what())); + } + + } else { + error("No model path provided for FarDetectorMLReconstruction"); + } + } + + + + void FarDetectorMLReconstruction::process( + const FarDetectorMLReconstruction::Input& input, + const FarDetectorMLReconstruction::Output& output) { + + const auto [inputTracks,beamElectrons] = input; + auto [outputFarDetectorMLTrajectories, outputFarDetectorMLTrackParameters, outputFarDetectorMLTracks] = output; + + //Set beam energy from first MCBeamElectron, using std::call_once + std::call_once(m_initBeamE,[&](){ + // Check if beam electrons are present + if(beamElectrons->size() == 0){ + error("No beam electrons found keeping default 10GeV beam energy."); + return; + } + m_beamE = beamElectrons->at(0).getEnergy(); + //Round beam energy to nearest GeV - Should be 5, 10 or 18GeV + m_beamE = round(m_beamE); + }); + + // Reconstructed particle members which don't change + std::int32_t type = 0; // Check? + float charge = -1; + + for(const auto& track: *inputTracks){ + + auto pos = track.getLoc(); + auto trackphi = track.getPhi(); + auto tracktheta = track.getTheta(); + + nnInput[FarDetectorMLNNIndexIn::PosY] = pos.a; + nnInput[FarDetectorMLNNIndexIn::PosZ] = pos.b; + nnInput[FarDetectorMLNNIndexIn::DirX] = sin(trackphi)*sin(tracktheta); + nnInput[FarDetectorMLNNIndexIn::DirY] = cos(trackphi)*sin(tracktheta); + + auto values = m_method->GetRegressionValues(); + + edm4hep::Vector3f momentum = {values[FarDetectorMLNNIndexOut::MomX],values[FarDetectorMLNNIndexOut::MomY],values[FarDetectorMLNNIndexOut::MomZ]}; + + // log out the momentum components and magnitude + trace("Prescaled Output Momentum: x {}, y {}, z {}",values[FarDetectorMLNNIndexOut::MomX],values[FarDetectorMLNNIndexOut::MomY],values[FarDetectorMLNNIndexOut::MomZ]); + trace("Prescaled Momentum: {}",edm4eic::magnitude(momentum)); + + // Scale momentum magnitude + momentum = momentum*m_beamE; + trace("Scaled Momentum: {}",edm4eic::magnitude(momentum)); + + // Track parameter variables + // TODO: Add time and momentum errors + // Plane Point + edm4hep::Vector2f loc(0,0); // Vertex estimate + uint64_t surface = 0; //Not used in this context + float theta = edm4eic::anglePolar(momentum); + float phi = edm4eic::angleAzimuthal(momentum); + float qOverP = charge/edm4eic::magnitude(momentum); + float time; + // PDG + int32_t pdg = 11; + // Point Error + edm4eic::Cov6f error; + + edm4eic::TrackParameters params = outputFarDetectorMLTrackParameters->create(type,surface,loc,theta,phi,qOverP,time,pdg,error); + + auto trajectory = outputFarDetectorMLTrajectories->create(); + trajectory.addToTrackParameters(params); + + int32_t trackType = 0; + edm4hep::Vector3f position = {0,0,0}; + float timeError; + float charge = -1; + float chi2 = 0; + uint32_t ndf = 0; + + auto outTrack = outputFarDetectorMLTracks->create(trackType,position,momentum,error,time,timeError,charge,chi2,ndf); + outTrack.setTrajectory(trajectory); + + } + + } + +} diff --git a/src/algorithms/fardetectors/FarDetectorMLReconstruction.h b/src/algorithms/fardetectors/FarDetectorMLReconstruction.h new file mode 100644 index 0000000000..a0f72141a7 --- /dev/null +++ b/src/algorithms/fardetectors/FarDetectorMLReconstruction.h @@ -0,0 +1,67 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2023, Simon Gardner + +#pragma once + +#include +#include +#include +#include +#include +#include +// Event Model related classes +#include +#include +#include +#include + +#include "FarDetectorMLReconstructionConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + + enum FarDetectorMLNNIndexIn{PosY,PosZ,DirX,DirY}; + enum FarDetectorMLNNIndexOut{MomX,MomY,MomZ}; + + using FarDetectorMLReconstructionAlgorithm = algorithms::Algorithm< + algorithms::Input< + edm4eic::TrackParametersCollection, + edm4hep::MCParticleCollection + >, + algorithms::Output< + edm4eic::TrajectoryCollection, + edm4eic::TrackParametersCollection, + edm4eic::TrackCollection + > + >; + + class FarDetectorMLReconstruction + : public FarDetectorMLReconstructionAlgorithm, + public WithPodConfig { + + public: + FarDetectorMLReconstruction(std::string_view name) + : FarDetectorMLReconstructionAlgorithm{name, + {"TrackParameters","BeamElectrons"}, + {"Trajectory","TrackParameters","Track"}, + "Reconstruct track parameters using ML method."} {} + + + /** One time initialization **/ + void init(); + + /** Event by event processing **/ + void process(const Input&, const Output&); + + //----- Define constants here ------ + + private: + TMVA::Reader* m_reader{nullptr}; + TMVA::MethodBase* m_method{nullptr}; + float m_beamE{10.0}; + std::once_flag m_initBeamE; + float nnInput[4] = {0.0,0.0,0.0,0.0}; + + }; + +} // eicrecon diff --git a/src/algorithms/fardetectors/FarDetectorMLReconstructionConfig.h b/src/algorithms/fardetectors/FarDetectorMLReconstructionConfig.h new file mode 100644 index 0000000000..29276a8bbe --- /dev/null +++ b/src/algorithms/fardetectors/FarDetectorMLReconstructionConfig.h @@ -0,0 +1,15 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2023, Simon Gardner + +#pragma once + +#include + +namespace eicrecon { + struct FarDetectorMLReconstructionConfig { + + std::string modelPath; + std::string methodName; + + }; +} diff --git a/src/detectors/LOWQ2/LOWQ2.cc b/src/detectors/LOWQ2/LOWQ2.cc index f512b19bc9..34e7a19676 100644 --- a/src/detectors/LOWQ2/LOWQ2.cc +++ b/src/detectors/LOWQ2/LOWQ2.cc @@ -6,11 +6,12 @@ #include #include #include +#include #include #include #include -#include #include +#include #include #include @@ -18,11 +19,12 @@ #include "algorithms/meta/SubDivideFunctors.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/digi/SiliconTrackerDigi_factory.h" +#include "factories/fardetectors/FarDetectorLinearProjection_factory.h" #include "factories/fardetectors/FarDetectorLinearTracking_factory.h" +#include "factories/fardetectors/FarDetectorMLReconstruction_factory.h" #include "factories/fardetectors/FarDetectorTrackerCluster_factory.h" -#include "factories/fardetectors/FarDetectorLinearProjection_factory.h" -#include "factories/meta/SubDivideCollection_factory.h" #include "factories/meta/CollectionCollector_factory.h" +#include "factories/meta/SubDivideCollection_factory.h" extern "C" { void InitPlugin(JApplication *app) { @@ -30,6 +32,8 @@ extern "C" { using namespace eicrecon; + std::string tracker_readout = "TaggerTrackerHits"; + // Digitization of silicon hits app->Add(new JOmniFactoryGeneratorT( "TaggerTrackerRawHits", @@ -119,9 +123,9 @@ extern "C" { // Combine the tracks from each module into one collection app->Add(new JOmniFactoryGeneratorT>( - "TaggerTrackerTracks", + "TaggerTrackerTrackSegments", outputTrackTags, - {"TaggerTrackerTracks"}, + {"TaggerTrackerTrackSegments"}, app ) ); @@ -129,7 +133,7 @@ extern "C" { // Project tracks onto a plane app->Add(new JOmniFactoryGeneratorT( "TaggerTrackerProjectedTracks", - {"TaggerTrackerTracks"}, + {"TaggerTrackerTrackSegments"}, {"TaggerTrackerProjectedTracks"}, { .plane_position = {0.0,0.0,0.0}, @@ -139,5 +143,17 @@ extern "C" { app )); + // Vector reconstruction at origin + app->Add(new JOmniFactoryGeneratorT( + "TaggerTrackerTrajectories", + {"TaggerTrackerProjectedTracks","MCBeamElectrons"}, + {"TaggerTrackerTrajectories","TaggerTrackerTrackParameters","TaggerTrackerTracks"}, + { + .modelPath = "calibrations/tmva/LowQ2_DNN_CPU.weights.xml", + .methodName = "DNN_CPU", + }, + app + )); + } } diff --git a/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h b/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h new file mode 100644 index 0000000000..13ace79e81 --- /dev/null +++ b/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h @@ -0,0 +1,55 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Simon Gardner + +#pragma once + +#include + +// Event Model related classes +#include +#include +#include + +#include "extensions/jana/JOmniFactory.h" +#include +#include +#include + +namespace eicrecon { + +class FarDetectorMLReconstruction_factory : + public JOmniFactory { + +public: + using AlgoT = eicrecon::FarDetectorMLReconstruction; +private: + std::unique_ptr m_algo; + + PodioInput m_trackparam_input {this}; + PodioInput m_beamelectrons_input {this}; + PodioOutput m_trajectory_output {this}; + PodioOutput m_trackparam_output {this}; + PodioOutput m_track_output {this}; + + + ParameterRef m_modelPath {this, "modelPath", config().modelPath }; + ParameterRef m_methodName {this, "methodName", config().methodName }; + + +public: + void Configure() { + m_algo = std::make_unique(GetPrefix()); + m_algo->level((algorithms::LogLevel)logger()->level()); + m_algo->applyConfig(config()); + m_algo->init(); + } + + void ChangeRun(int64_t run_number) { + } + + void Process(int64_t run_number, uint64_t event_number) { + m_algo->process({m_trackparam_input(),m_beamelectrons_input()}, {m_trajectory_output().get(), m_trackparam_output().get(), m_track_output().get()}); + } + }; + +} // eicrecon diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 0876fa8857..94e9642425 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -1,15 +1,15 @@ #include "JEventProcessorPODIO.h" -#include - #include #include #include #include +#include #include #include #include +#include #include #include @@ -143,6 +143,9 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "TaggerTrackerM1Tracks", "TaggerTrackerM2Tracks", "TaggerTrackerProjectedTracks", + "TaggerTrackerTracks", + "TaggerTrackerTrajectories", + "TaggerTrackerTrackParameters", // Forward & Far forward hits "B0TrackerRecHits", From d15e9ab6a476a37336b38efbd351638da813a67e Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Sun, 2 Jun 2024 19:09:31 -0400 Subject: [PATCH 04/10] Cluster axis direction (#1391) Added code to produce the direction of clusters in CalorimeterClusterRecoCoG ### Briefly, what does this PR introduce? adds the cluster direction (i.e., the eigenvector corresponding to the largest eigenvector of the matrix corresponding to the inertia matrix of the cluster) ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [X] New feature (issue #1462 ) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [X] Tests for the changes have been added - [X] Documentation has been added / updated - [X] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? no ### Does this PR change default behavior? Yes. While calculating the shape parameters, the hits are assigned weights, in a similar manner to how they are assigned weights when determining the cluster position. Previously, there were some inconsistencies in how these weights were determined when using certain options in the cluster fitting, which were fixed in order to make this feature (the cluster axis direction) work properly. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .../calorimetry/CalorimeterClusterRecoCoG.cc | 38 ++++++- src/tests/algorithms_test/CMakeLists.txt | 1 + .../calorimetry_CalorimeterClusterRecoCoG.cc | 107 ++++++++++++++++++ 3 files changed, 141 insertions(+), 5 deletions(-) create mode 100644 src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index 46a3a083ac..60716fc61e 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -26,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -246,12 +248,13 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co Eigen::Vector3f sum1_3D = Eigen::Vector3f::Zero(); Eigen::Vector2cf eigenValues_2D = Eigen::Vector2cf::Zero(); Eigen::Vector3cf eigenValues_3D = Eigen::Vector3cf::Zero(); - + //the axis is the direction of the eigenvalue corresponding to the largest eigenvalue. + double axis_x=0, axis_y=0, axis_z=0; if (cl.getNhits() > 1) { for (const auto& hit : pcl.getHits()) { - float w = weightFunc(hit.getEnergy(), cl.getEnergy(), m_cfg.logWeightBase, 0); + float w = weightFunc(hit.getEnergy(), totalE, logWeightBase, 0); // theta, phi Eigen::Vector2f pos2D( edm4hep::utils::anglePolar( hit.getPosition() ), edm4hep::utils::angleAzimuthal( hit.getPosition() ) ); @@ -273,8 +276,8 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co w_sum += w; } + radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); if( w_sum > 0 ) { - radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); dispersion = sqrt( dispersion / w_sum ); // normalize matrices @@ -289,11 +292,33 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co // Solve for eigenvalues. Corresponds to cluster's 2nd moments (widths) Eigen::EigenSolver es_2D(cov2, false); // set to true for eigenvector calculation - Eigen::EigenSolver es_3D(cov3, false); // set to true for eigenvector calculation + Eigen::EigenSolver es_3D(cov3, true); // set to true for eigenvector calculation // eigenvalues of symmetric real matrix are always real eigenValues_2D = es_2D.eigenvalues(); eigenValues_3D = es_3D.eigenvalues(); + //find the eigenvector corresponding to the largest eigenvalue + auto eigenvectors= es_3D.eigenvectors(); + auto max_eigenvalue_it = std::max_element( + eigenValues_3D.begin(), + eigenValues_3D.end(), + [](auto a, auto b) { + return std::real(a) < std::real(b); + } + ); + auto axis = eigenvectors.col(std::distance( + eigenValues_3D.begin(), + max_eigenvalue_it + )); + axis_x=axis(0,0).real(); + axis_y=axis(1,0).real(); + axis_z=axis(2,0).real(); + double norm=sqrt(axis_x*axis_x+axis_y*axis_y+axis_z*axis_z); + if (norm!=0){ + axis_x/=norm; + axis_y/=norm; + axis_z/=norm; + } } } @@ -304,7 +329,10 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co cl.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z cluster width 1 cl.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z cluster width 2 cl.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z cluster width 3 - + //last 3 shape parameters are the components of the axis direction + cl.addToShapeParameters( axis_x ); + cl.addToShapeParameters( axis_y ); + cl.addToShapeParameters( axis_z ); return std::move(cl); } diff --git a/src/tests/algorithms_test/CMakeLists.txt b/src/tests/algorithms_test/CMakeLists.txt index 3ca06a1d97..50ee83cdd0 100644 --- a/src/tests/algorithms_test/CMakeLists.txt +++ b/src/tests/algorithms_test/CMakeLists.txt @@ -8,6 +8,7 @@ add_executable( calorimetry_CalorimeterIslandCluster.cc tracking_SiliconSimpleCluster.cc calorimetry_CalorimeterHitDigi.cc + calorimetry_CalorimeterClusterRecoCoG.cc calorimetry_HEXPLIT.cc pid_MergeTracks.cc pid_MergeParticleID.cc diff --git a/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc new file mode 100644 index 0000000000..32f11942d3 --- /dev/null +++ b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc @@ -0,0 +1,107 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Sebouh Paul + + +#include // for MeV, mm, keV, ns +#include // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE +#include // for CalorimeterHitCollection, MutableCalorimeterHit, CalorimeterHitMutableCollectionIterator +#include +#include +#include +#include +#include // for Vector3f +#include +#include +#include // for level_enum +#include // for logger +#include // for default_logger +#include // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access +#include +#include + +#include "algorithms/calorimetry/CalorimeterClusterRecoCoG.h" // for CalorimeterClusterRecoCoG +#include "algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h" // for CalorimeterClusterRecoCoGConfig + +using eicrecon::CalorimeterClusterRecoCoG; +using eicrecon::CalorimeterClusterRecoCoGConfig; + +using edm4eic::CalorimeterHit; + +TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) { + CalorimeterClusterRecoCoG algo("CalorimeterClusterRecoCoG"); + + std::shared_ptr logger = spdlog::default_logger()->clone("CalorimeterClusterRecoCoG"); + logger->set_level(spdlog::level::trace); + + CalorimeterClusterRecoCoGConfig cfg; + cfg.energyWeight = "log"; + cfg.sampFrac = 0.0203, + cfg.logWeightBaseCoeffs={5.0,0.65,0.31}, + cfg.logWeightBase_Eref=50*dd4hep::GeV, + + + algo.applyConfig(cfg); + algo.init(); + + edm4eic::CalorimeterHitCollection hits_coll; + edm4eic::ProtoClusterCollection pclust_coll; + edm4hep::SimCalorimeterHitCollection simhits; + auto assoc = std::make_unique(); + auto clust_coll = std::make_unique(); + + //create a protocluster with 3 hits + auto pclust = pclust_coll.create(); + edm4hep::Vector3f position({0,0,0*dd4hep::mm}); + + auto hit1 = hits_coll.create(); + hit1.setCellID(0); + hit1.setEnergy(0.1*dd4hep::GeV); + hit1.setEnergyError(0); + hit1.setTime(0); + hit1.setTimeError(0); + hit1.setPosition(position); + hit1.setDimension({0,0,0}); + hit1.setLocal(position); + pclust.addToHits(hit1); + + position={0,0, 1*dd4hep::mm}; + auto hit2 = hits_coll.create(); + hit2.setCellID(0); + hit2.setEnergy(0.1*dd4hep::GeV); + hit2.setEnergyError(0); + hit2.setTime(0); + hit2.setTimeError(0); + hit2.setPosition(position); + hit2.setDimension({0,0,0}); + hit2.setLocal(position); + pclust.addToHits(hit2); + + position={0,0, 2*dd4hep::mm}; + auto hit3 = hits_coll.create();//0, 0.1*dd4hep::GeV, 0,0,0,position, {0,0,0}, 0,0, position); + hit3.setCellID(0); + hit3.setEnergy(0.1*dd4hep::GeV); + hit3.setEnergyError(0); + hit3.setTime(0); + hit3.setTimeError(0); + hit3.setPosition(position); + hit3.setDimension({0,0,0}); + hit3.setLocal(position); + pclust.addToHits(hit3); + pclust.addToWeights(1);pclust.addToWeights(1);pclust.addToWeights(1); + + // Constructing input and output as per the algorithm's expected signature + auto input = std::make_tuple(&pclust_coll, &simhits); + auto output = std::make_tuple(clust_coll.get(), assoc.get()); + + algo.process(input, output); + + + for (auto clust : *clust_coll){ + // require that this cluster's axis is 0,0,1 + REQUIRE(clust.getShapeParameters()[7] == 0); + REQUIRE(clust.getShapeParameters()[8] == 0); + REQUIRE(abs(clust.getShapeParameters()[9]) == 1); + } + + +} From 0546b276d5ae07c084db79abd319b96faca2d723 Mon Sep 17 00:00:00 2001 From: Tyler Kutz Date: Mon, 3 Jun 2024 12:13:49 -0400 Subject: [PATCH 05/10] Reconstructed inclusive kinematics algorithms now use centralized algorithms for scattered electron ID and hadronic final state calculation (#1453) ### Briefly, what does this PR introduce? All reconstructed inclusive kinematics algorithms now obtain: - Scattered electron from `ScatteredElectronsTruth` algorithm - Hadronic final state from `HadronicFinalState` algorithm Previously, these steps were repeated as-needed in each individual inclusive kinematics algorithm. Now, improvements/fixes made to `ScatteredElectronTruth` or `HadronicFinalState` will automatically be included in the inclusive kinematics algorithms. This also fixes issue #1318 , as the `HadronicFinalState` algorithm uses all reconstructed particles, not just reconstructed charged particles. ### What kind of change does this PR introduce? - [ x] Bug fix (issue #1318 ) - [ ] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [ ] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No, although the specific values of the reconstructed inclusive kinematics might be different (see below) the output collections are in the same format. ### Does this PR change default behavior? The electron reconstruction method should be identical, as the only difference is where the scattered electron is found. The reconstruction methods using any information from the hadronic final state (JB, DA, [e]Sigma) will be different, as the hadronic final state now includes all reconstructed particles, not just charged particles. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Dmitry Kalinkin --- src/algorithms/reco/HadronicFinalState.cc | 37 +++---- src/algorithms/reco/InclusiveKinematicsDA.cc | 101 ++++-------------- src/algorithms/reco/InclusiveKinematicsDA.h | 6 +- .../reco/InclusiveKinematicsElectron.cc | 47 ++------ .../reco/InclusiveKinematicsElectron.h | 6 +- src/algorithms/reco/InclusiveKinematicsJB.cc | 88 +++------------ src/algorithms/reco/InclusiveKinematicsJB.h | 6 +- .../reco/InclusiveKinematicsSigma.cc | 98 ++++------------- .../reco/InclusiveKinematicsSigma.h | 6 +- .../reco/InclusiveKinematicseSigma.cc | 100 ++++------------- .../reco/InclusiveKinematicseSigma.h | 6 +- ...InclusiveKinematicsReconstructed_factory.h | 6 +- src/global/reco/reco.cc | 46 ++++---- 13 files changed, 140 insertions(+), 413 deletions(-) diff --git a/src/algorithms/reco/HadronicFinalState.cc b/src/algorithms/reco/HadronicFinalState.cc index 96de76822e..9cc84121b9 100644 --- a/src/algorithms/reco/HadronicFinalState.cc +++ b/src/algorithms/reco/HadronicFinalState.cc @@ -16,7 +16,6 @@ #include #include #include -#include #include "Beam.h" #include "Boost.h" @@ -106,35 +105,27 @@ namespace eicrecon { auto hfs = hadronicfinalstate->create(0., 0., 0.); for (const auto& p: *rcparts) { - bool isHadron = true; // Check if it's the scattered electron - if (p.getObjectID().index == ef_rc_id) isHadron = false; - // Check for non-hadron PDG codes - if (p.getPDG() == 11) isHadron = false; - if (p.getPDG() == 22) isHadron = false; - if (p.getPDG() == 13) isHadron = false; - // If it's the scattered electron or not a hadron, skip - if(!isHadron) continue; - - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab); - - pxsum += hf_boosted.Px(); - pysum += hf_boosted.Py(); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - - hfs.addToHadrons(p); - + if (p.getObjectID().index != ef_rc_id) { + // Lorentz vector in lab frame + PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); + // Boost to colinear frame + PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab); + + pxsum += hf_boosted.Px(); + pysum += hf_boosted.Py(); + pzsum += hf_boosted.Pz(); + Esum += hf_boosted.E(); + + hfs.addToHadrons(p); + } } // Hadronic final state calculations auto sigma = Esum - pzsum; auto pT = sqrt(pxsum*pxsum + pysum*pysum); - auto gamma = (pT*pT - sigma*sigma)/(pT*pT + sigma*sigma); + auto gamma = acos((pT*pT - sigma*sigma)/(pT*pT + sigma*sigma)); hfs.setSigma(sigma); hfs.setPT(pT); diff --git a/src/algorithms/reco/InclusiveKinematicsDA.cc b/src/algorithms/reco/InclusiveKinematicsDA.cc index 07a411a90f..dd58f36a45 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.cc +++ b/src/algorithms/reco/InclusiveKinematicsDA.cc @@ -1,16 +1,17 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Wouter Deconinck +#include +#if EDM4EIC_VERSION_MAJOR >= 6 + #include #include #include #include -#include #include #include #include #include -#include #include #include @@ -35,7 +36,7 @@ namespace eicrecon { const InclusiveKinematicsDA::Input& input, const InclusiveKinematicsDA::Output& output) const { - const auto [mcparts, rcparts, rcassoc] = input; + const auto [mcparts, escat, hfs] = input; auto [kinematics] = output; // Get incoming electron beam @@ -66,100 +67,38 @@ namespace eicrecon { m_crossingAngle) ); - // Get first scattered electron - const auto ef_coll = find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); - return; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc->begin(), - // rcassoc->end(), - // [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); }); - auto ef_assoc = rcassoc->begin(); - for (; ef_assoc != rcassoc->end(); ++ef_assoc) { - if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) { - break; - } - } - if (!(ef_assoc != rcassoc->end())) { - m_log->debug("Truth scattered electron not in reconstructed particles"); - return; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to calculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - //Sums in colinear frame - double pxsum = 0; - double pysum = 0; - double pzsum = 0; - double Esum = 0; - double theta_e = 0; - // Get boost to colinear frame auto boost = determine_boost(ei, pi); - for (const auto& p: *rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - // Lorentz vector in lab frame - PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector e_boosted = apply_boost(boost, e_lab); - - theta_e = e_boosted.Theta(); - - // Sum over all particles other than scattered electron - } else { - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab); - - pxsum += hf_boosted.Px(); - pysum += hf_boosted.Py(); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } + // Get electron angle + auto kf = escat->at(0); + PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy()); + PxPyPzEVector e_boosted = apply_boost(boost, e_lab); + auto theta_e = e_boosted.Theta(); - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; + // Get hadronic final state variables + auto sigma_h = hfs->at(0).getSigma(); + auto ptsum = hfs->at(0).getPT(); + auto gamma_h = hfs->at(0).getGamma(); - // If no scattered hadron was found + // Sigma zero or negative if (sigma_h <= 0) { - m_log->debug("No scattered hadron found"); + m_log->debug("Sigma zero or negative"); return; } - auto ptsum = sqrt(pxsum*pxsum + pysum*pysum); - auto theta_h = 2.*atan(sigma_h/ptsum); - // Calculate kinematic variables - const auto y_da = tan(theta_h/2.) / ( tan(theta_e/2.) + tan(theta_h/2.) ); - const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(theta_h/2.)) ); + const auto y_da = tan(gamma_h/2.) / ( tan(theta_e/2.) + tan(gamma_h/2.) ); + const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(gamma_h/2.)) ); const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da); const auto nu_da = Q2_da / (2.*m_proton*x_da); const auto W_da = sqrt(m_proton*m_proton + 2*m_proton*nu_da - Q2_da); auto kin = kinematics->create(x_da, Q2_da, W_da, y_da, nu_da); - kin.setScat(ef_rc); + kin.setScat(kf); m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } -} // namespace Jug::Reco +} +#endif diff --git a/src/algorithms/reco/InclusiveKinematicsDA.h b/src/algorithms/reco/InclusiveKinematicsDA.h index af4cea4fb2..c1d3f72468 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.h +++ b/src/algorithms/reco/InclusiveKinematicsDA.h @@ -5,9 +5,9 @@ #include #include -#include #include #include +#include #include #include #include @@ -20,7 +20,7 @@ namespace eicrecon { algorithms::Input< edm4hep::MCParticleCollection, edm4eic::ReconstructedParticleCollection, - edm4eic::MCRecoParticleAssociationCollection + edm4eic::HadronicFinalStateCollection >, algorithms::Output< edm4eic::InclusiveKinematicsCollection @@ -33,7 +33,7 @@ namespace eicrecon { public: InclusiveKinematicsDA(std::string_view name) : InclusiveKinematicsDAAlgorithm{name, - {"MCParticles", "inputParticles", "inputAssociations"}, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, {"inclusiveKinematics"}, "Determine inclusive kinematics using double-angle method."} {} diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.cc b/src/algorithms/reco/InclusiveKinematicsElectron.cc index 4d27bb4c83..0a4781db91 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.cc +++ b/src/algorithms/reco/InclusiveKinematicsElectron.cc @@ -1,16 +1,17 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022, 2023 Wouter Deconinck, Tooba Ali +#include +#if EDM4EIC_VERSION_MAJOR >= 6 + #include #include #include #include -#include #include #include #include #include -#include #include #include #include @@ -35,7 +36,7 @@ namespace eicrecon { const InclusiveKinematicsElectron::Input& input, const InclusiveKinematicsElectron::Output& output) const { - const auto [mcparts, rcparts, rcassoc] = input; + const auto [mcparts, escat, hfs] = input; auto [kinematics] = output; // 1. find_if @@ -114,38 +115,11 @@ namespace eicrecon { m_crossingAngle) ); - // Get first scattered electron - const auto ef_coll = find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); - return; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc->begin(), - // rcassoc->end(), - // [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); }); - auto ef_assoc = rcassoc->begin(); - for (; ef_assoc != rcassoc->end(); ++ef_assoc) { - if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) { - break; - } - } - if (!(ef_assoc != rcassoc->end())) { - m_log->debug("Truth scattered electron not in reconstructed particles"); - return; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get outgoing scattered electron - // Use the true scattered electron from the MC information + // Get scattered electron std::vector electrons; - for (const auto& p: *rcparts) { - if (p.getObjectID().index == ef_rc_id) { - electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - break; - } + for (const auto& p: *escat) { + electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); + break; } // If no scattered electron was found @@ -164,10 +138,11 @@ namespace eicrecon { const auto x = Q2 / (2. * q_dot_pi); const auto W = sqrt(m_proton*m_proton + 2.*q_dot_pi - Q2); auto kin = kinematics->create(x, Q2, W, y, nu); - kin.setScat(ef_rc); + kin.setScat(escat->at(0)); m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } -} // namespace Jug::Reco +} +#endif diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.h b/src/algorithms/reco/InclusiveKinematicsElectron.h index 2425b8a79a..e95cc55138 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.h +++ b/src/algorithms/reco/InclusiveKinematicsElectron.h @@ -5,9 +5,9 @@ #include #include -#include #include #include +#include #include #include #include @@ -20,7 +20,7 @@ namespace eicrecon { algorithms::Input< edm4hep::MCParticleCollection, edm4eic::ReconstructedParticleCollection, - edm4eic::MCRecoParticleAssociationCollection + edm4eic::HadronicFinalStateCollection >, algorithms::Output< edm4eic::InclusiveKinematicsCollection @@ -33,7 +33,7 @@ namespace eicrecon { public: InclusiveKinematicsElectron(std::string_view name) : InclusiveKinematicsElectronAlgorithm{name, - {"MCParticles", "inputParticles", "inputAssociations"}, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, {"inclusiveKinematics"}, "Determine inclusive kinematics using electron method."} {} diff --git a/src/algorithms/reco/InclusiveKinematicsJB.cc b/src/algorithms/reco/InclusiveKinematicsJB.cc index cb8d19c2cf..324093b1aa 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.cc +++ b/src/algorithms/reco/InclusiveKinematicsJB.cc @@ -1,21 +1,19 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Wouter Deconinck +#include +#if EDM4EIC_VERSION_MAJOR >= 6 + #include #include #include +#include #include -#include -#include -#include -#include #include -#include #include #include #include "Beam.h" -#include "Boost.h" #include "InclusiveKinematicsJB.h" using ROOT::Math::PxPyPzEVector; @@ -35,7 +33,7 @@ namespace eicrecon { const InclusiveKinematicsJB::Input& input, const InclusiveKinematicsJB::Output& output) const { - const auto [mcparts, rcparts, rcassoc] = input; + const auto [mcparts, escat, hfs] = input; auto [kinematics] = output; // Get incoming electron beam @@ -66,73 +64,10 @@ namespace eicrecon { m_crossingAngle) ); - // Get first scattered electron - const auto ef_coll = find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); - return; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc->begin(), - // rcassoc->end(), - // [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); }); - auto ef_assoc = rcassoc->begin(); - for (; ef_assoc != rcassoc->end(); ++ef_assoc) { - if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) { - break; - } - } - if (!(ef_assoc != rcassoc->end())) { - m_log->debug("Truth scattered electron not in reconstructed particles"); - return; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles other than the scattered electron - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to calculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - // Sums in colinear frame - double pxsum = 0; - double pysum = 0; - double pzsum = 0; - double Esum = 0; - - // Get boost to colinear frame - auto boost = determine_boost(ei, pi); - - for (const auto& p: *rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - - // Sum over all particles other than scattered electron - } else { - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab); - - pxsum += hf_boosted.Px(); - pysum += hf_boosted.Py(); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } - - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; - auto ptsum = sqrt(pxsum*pxsum + pysum*pysum); + // Get hadronic final state variables + auto sigma_h = hfs->at(0).getSigma(); + auto ptsum = hfs->at(0).getPT(); + auto gamma_h = hfs->at(0).getGamma(); // Sigma zero or negative if (sigma_h <= 0) { @@ -147,10 +82,11 @@ namespace eicrecon { const auto nu_jb = Q2_jb / (2.*m_proton*x_jb); const auto W_jb = sqrt(m_proton*m_proton + 2*m_proton*nu_jb - Q2_jb); auto kin = kinematics->create(x_jb, Q2_jb, W_jb, y_jb, nu_jb); - kin.setScat(ef_rc); + kin.setScat(escat->at(0)); m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } -} // namespace Jug::Reco +} +#endif diff --git a/src/algorithms/reco/InclusiveKinematicsJB.h b/src/algorithms/reco/InclusiveKinematicsJB.h index 7eac50d22d..24526c04b2 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.h +++ b/src/algorithms/reco/InclusiveKinematicsJB.h @@ -5,9 +5,9 @@ #include #include -#include #include #include +#include #include #include #include @@ -20,7 +20,7 @@ namespace eicrecon { algorithms::Input< edm4hep::MCParticleCollection, edm4eic::ReconstructedParticleCollection, - edm4eic::MCRecoParticleAssociationCollection + edm4eic::HadronicFinalStateCollection >, algorithms::Output< edm4eic::InclusiveKinematicsCollection @@ -33,7 +33,7 @@ namespace eicrecon { public: InclusiveKinematicsJB(std::string_view name) : InclusiveKinematicsJBAlgorithm{name, - {"MCParticles", "inputParticles", "inputAssociations"}, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, {"inclusiveKinematics"}, "Determine inclusive kinematics using Jacquet-Blondel method."} {} diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.cc b/src/algorithms/reco/InclusiveKinematicsSigma.cc index d25b38555c..020d1dc5fb 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicsSigma.cc @@ -1,16 +1,16 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Wouter Deconinck, Barak Schmookler +#include +#if EDM4EIC_VERSION_MAJOR >= 6 + #include #include #include #include -#include -#include #include #include #include -#include #include #include @@ -35,7 +35,7 @@ namespace eicrecon { const InclusiveKinematicsSigma::Input& input, const InclusiveKinematicsSigma::Output& output) const { - const auto [mcparts, rcparts, rcassoc] = input; + const auto [mcparts, escat, hfs] = input; auto [kinematics] = output; // Get incoming electron beam @@ -66,89 +66,28 @@ namespace eicrecon { m_crossingAngle) ); - // Get first scattered electron - const auto ef_coll = find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); - return; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc->begin(), - // rcassoc->end(), - // [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); }); - auto ef_assoc = rcassoc->begin(); - for (; ef_assoc != rcassoc->end(); ++ef_assoc) { - if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) { - break; - } - } - if (!(ef_assoc != rcassoc->end())) { - m_log->debug("Truth scattered electron not in reconstructed particles"); - return; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to calculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - //Sums in colinear frame - double pxsum = 0; - double pysum = 0; - double pzsum = 0; - double Esum = 0; - - double pt_e = 0; - double sigma_e = 0; - // Get boost to colinear frame auto boost = determine_boost(ei, pi); - for (const auto& p: *rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - // Lorentz vector in lab frame - PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector e_boosted = apply_boost(boost, e_lab); - - pt_e = e_boosted.Pt(); - sigma_e = e_boosted.E() - e_boosted.Pz(); - - // Sum over all particles other than scattered electron - } else { - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab); - - pxsum += hf_boosted.Px(); - pysum += hf_boosted.Py(); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } + // Get electron variables + auto kf = escat->at(0); + PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy()); + PxPyPzEVector e_boosted = apply_boost(boost, e_lab); + auto pt_e = e_boosted.Pt(); + auto sigma_e = e_boosted.E() - e_boosted.Pz(); - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; - auto sigma_tot = sigma_e + sigma_h; + // Get hadronic final state variables + auto sigma_h = hfs->at(0).getSigma(); + auto ptsum = hfs->at(0).getPT(); + auto gamma_h = hfs->at(0).getGamma(); if (sigma_h <= 0) { m_log->debug("No scattered electron found or sigma zero or negative"); return; } + auto sigma_tot = sigma_e + sigma_h; + // Calculate kinematic variables const auto y_sig = sigma_h / sigma_tot; const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig); @@ -156,10 +95,11 @@ namespace eicrecon { const auto nu_sig = Q2_sig / (2.*m_proton*x_sig); const auto W_sig = sqrt(m_proton*m_proton + 2*m_proton*nu_sig - Q2_sig); auto kin = kinematics->create(x_sig, Q2_sig, W_sig, y_sig, nu_sig); - kin.setScat(ef_rc); + kin.setScat(kf); m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } -} // namespace Jug::Reco +} +#endif diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.h b/src/algorithms/reco/InclusiveKinematicsSigma.h index 8b2ca96be5..d3849806c3 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.h +++ b/src/algorithms/reco/InclusiveKinematicsSigma.h @@ -5,9 +5,9 @@ #include #include -#include #include #include +#include #include #include #include @@ -20,7 +20,7 @@ namespace eicrecon { algorithms::Input< edm4hep::MCParticleCollection, edm4eic::ReconstructedParticleCollection, - edm4eic::MCRecoParticleAssociationCollection + edm4eic::HadronicFinalStateCollection >, algorithms::Output< edm4eic::InclusiveKinematicsCollection @@ -33,7 +33,7 @@ namespace eicrecon { public: InclusiveKinematicsSigma(std::string_view name) : InclusiveKinematicsSigmaAlgorithm{name, - {"MCParticles", "inputParticles", "inputAssociations"}, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, {"inclusiveKinematics"}, "Determine inclusive kinematics using Sigma method."} {} diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.cc b/src/algorithms/reco/InclusiveKinematicseSigma.cc index 80880d5a90..b898e41da9 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicseSigma.cc @@ -1,16 +1,17 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Wouter Deconinck, Barak Schmookler +#include +#if EDM4EIC_VERSION_MAJOR >= 6 + #include #include #include #include -#include #include #include #include #include -#include #include #include @@ -35,7 +36,7 @@ namespace eicrecon { const InclusiveKinematicseSigma::Input& input, const InclusiveKinematicseSigma::Output& output) const { - const auto [mcparts, rcparts, rcassoc] = input; + const auto [mcparts, escat, hfs] = input; auto [kinematics] = output; // Get incoming electron beam @@ -66,90 +67,28 @@ namespace eicrecon { m_crossingAngle) ); - // Get first scattered electron - const auto ef_coll = find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); - return; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc->begin(), - // rcassoc->end(), - // [&ef_coll](const auto& a){ return a.getSim().getObjectID() == ef_coll[0].getObjectID(); }); - auto ef_assoc = rcassoc->begin(); - for (; ef_assoc != rcassoc->end(); ++ef_assoc) { - if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) { - break; - } - } - if (!(ef_assoc != rcassoc->end())) { - m_log->debug("Truth scattered electron not in reconstructed particles"); - return; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to calculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - //Sums in colinear frame - double pxsum = 0; - double pysum = 0; - double pzsum = 0; - double Esum = 0; - - double pt_e = 0; - double sigma_e = 0; - // Get boost to colinear frame auto boost = determine_boost(ei, pi); - for (const auto& p: *rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - // Lorentz vector in lab frame - PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector e_boosted = apply_boost(boost, e_lab); - - pt_e = e_boosted.Pt(); - sigma_e = e_boosted.E() - e_boosted.Pz(); - - // Sum over all particles other than scattered electron - } else { - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = apply_boost(boost, hf_lab); - - pxsum += hf_boosted.Px(); - pysum += hf_boosted.Py(); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } + // Get electron variables + auto kf = escat->at(0); + PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy()); + PxPyPzEVector e_boosted = apply_boost(boost, e_lab); + auto pt_e = e_boosted.Pt(); + auto sigma_e = e_boosted.E() - e_boosted.Pz(); - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; - auto sigma_tot = sigma_e + sigma_h; + // Get hadronic final state variables + auto sigma_h = hfs->at(0).getSigma(); + auto ptsum = hfs->at(0).getPT(); + auto gamma_h = hfs->at(0).getGamma(); - // If no scattered electron was found if (sigma_h <= 0) { - m_log->debug("No scattered electron found"); + m_log->debug("No scattered electron found or sigma zero or negative"); return; } + auto sigma_tot = sigma_e + sigma_h; + // Calculate kinematic variables const auto y_e = 1. - sigma_e / (2.*ei.energy()); const auto Q2_e = (pt_e*pt_e) / (1. - y_e); @@ -164,11 +103,12 @@ namespace eicrecon { const auto nu_esig = Q2_esig / (2.*m_proton*x_esig); const auto W_esig = sqrt(m_proton*m_proton + 2*m_proton*nu_esig - Q2_esig); auto kin = kinematics->create(x_esig, Q2_esig, W_esig, y_esig, nu_esig); - kin.setScat(ef_rc); + kin.setScat(kf); // Debugging output m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } -} // namespace Jug::Reco +} +#endif diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.h b/src/algorithms/reco/InclusiveKinematicseSigma.h index 9ed9fb7d0a..9a505f0eae 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.h +++ b/src/algorithms/reco/InclusiveKinematicseSigma.h @@ -5,9 +5,9 @@ #include #include -#include #include #include +#include #include #include #include @@ -20,7 +20,7 @@ namespace eicrecon { algorithms::Input< edm4hep::MCParticleCollection, edm4eic::ReconstructedParticleCollection, - edm4eic::MCRecoParticleAssociationCollection + edm4eic::HadronicFinalStateCollection >, algorithms::Output< edm4eic::InclusiveKinematicsCollection @@ -33,7 +33,7 @@ namespace eicrecon { public: InclusiveKinematicseSigma(std::string_view name) : InclusiveKinematicseSigmaAlgorithm{name, - {"MCParticles", "inputParticles", "inputAssociations"}, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, {"inclusiveKinematics"}, "Determine inclusive kinematics using e-Sigma method."} {} diff --git a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h index d86ee3d88e..7e99ba0253 100644 --- a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h +++ b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h @@ -25,8 +25,8 @@ class InclusiveKinematicsReconstructed_factory : std::unique_ptr m_algo; typename FactoryT::template PodioInput m_mc_particles_input {this}; - typename FactoryT::template PodioInput m_rc_particles_input {this}; - typename FactoryT::template PodioInput m_rc_particles_assoc_input {this}; + typename FactoryT::template PodioInput m_scattered_electron_input {this}; + typename FactoryT::template PodioInput m_hadronic_final_state_input {this}; typename FactoryT::template PodioOutput m_inclusive_kinematics_output {this}; public: @@ -40,7 +40,7 @@ class InclusiveKinematicsReconstructed_factory : } void Process(int64_t run_number, uint64_t event_number) { - m_algo->process({m_mc_particles_input(), m_rc_particles_input(), m_rc_particles_assoc_input()}, + m_algo->process({m_mc_particles_input(), m_scattered_electron_input(), m_hadronic_final_state_input()}, {m_inclusive_kinematics_output().get()}); } }; diff --git a/src/global/reco/reco.cc b/src/global/reco/reco.cc index 03e225d11b..1590a0091d 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -15,20 +15,23 @@ #include #include "algorithms/interfaces/WithPodConfig.h" + +#if EDM4EIC_VERSION_MAJOR >= 6 +#include "algorithms/reco/HadronicFinalState.h" #include "algorithms/reco/InclusiveKinematicsDA.h" #include "algorithms/reco/InclusiveKinematicsElectron.h" #include "algorithms/reco/InclusiveKinematicsJB.h" #include "algorithms/reco/InclusiveKinematicsSigma.h" #include "algorithms/reco/InclusiveKinematicseSigma.h" -#if EDM4EIC_VERSION_MAJOR >= 6 -#include "algorithms/reco/HadronicFinalState.h" #endif #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/meta/CollectionCollector_factory.h" #include "factories/meta/FilterMatching_factory.h" #include "factories/reco/FarForwardNeutronReconstruction_factory.h" #include "factories/reco/InclusiveKinematicsML_factory.h" +#if EDM4EIC_VERSION_MAJOR >= 6 #include "factories/reco/InclusiveKinematicsReconstructed_factory.h" +#endif #include "factories/reco/InclusiveKinematicsTruth_factory.h" #include "factories/reco/JetReconstruction_factory.h" #include "factories/reco/TransformBreitFrame_factory.h" @@ -102,26 +105,27 @@ void InitPlugin(JApplication *app) { )); - app->Add(new JOmniFactoryGeneratorT>( - "InclusiveKinematicsElectron", + app->Add(new JOmniFactoryGeneratorT( + "InclusiveKinematicsTruth", { - "MCParticles", - "ReconstructedChargedParticles", - "ReconstructedChargedParticleAssociations" + "MCParticles" }, { - "InclusiveKinematicsElectron" + "InclusiveKinematicsTruth" }, app )); - app->Add(new JOmniFactoryGeneratorT( - "InclusiveKinematicsTruth", +#if EDM4EIC_VERSION_MAJOR >= 6 + app->Add(new JOmniFactoryGeneratorT>( + "InclusiveKinematicsElectron", { - "MCParticles" + "MCParticles", + "ScatteredElectronsTruth", + "HadronicFinalState" }, { - "InclusiveKinematicsTruth" + "InclusiveKinematicsElectron" }, app )); @@ -130,8 +134,8 @@ void InitPlugin(JApplication *app) { "InclusiveKinematicsJB", { "MCParticles", - "ReconstructedChargedParticles", - "ReconstructedChargedParticleAssociations" + "ScatteredElectronsTruth", + "HadronicFinalState" }, { "InclusiveKinematicsJB" @@ -143,8 +147,8 @@ void InitPlugin(JApplication *app) { "InclusiveKinematicsDA", { "MCParticles", - "ReconstructedChargedParticles", - "ReconstructedChargedParticleAssociations" + "ScatteredElectronsTruth", + "HadronicFinalState" }, { "InclusiveKinematicsDA" @@ -156,8 +160,8 @@ void InitPlugin(JApplication *app) { "InclusiveKinematicseSigma", { "MCParticles", - "ReconstructedChargedParticles", - "ReconstructedChargedParticleAssociations" + "ScatteredElectronsTruth", + "HadronicFinalState" }, { "InclusiveKinematicseSigma" @@ -165,12 +169,13 @@ void InitPlugin(JApplication *app) { app )); + app->Add(new JOmniFactoryGeneratorT>( "InclusiveKinematicsSigma", { "MCParticles", - "ReconstructedChargedParticles", - "ReconstructedChargedParticleAssociations" + "ScatteredElectronsTruth", + "HadronicFinalState" }, { "InclusiveKinematicsSigma" @@ -189,6 +194,7 @@ void InitPlugin(JApplication *app) { }, app )); +#endif app->Add(new JOmniFactoryGeneratorT( "ReconstructedElectrons", From cec2ceba98f4c43a61aac66e0974a2d46bf9cee5 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Mon, 3 Jun 2024 15:55:19 -0500 Subject: [PATCH 06/10] fix: include vector_utils in ImagingClusterReco (#1482) ### Briefly, what does this PR introduce? ImagingClusterReco uses multiplication of Vector3f with float, which is defined in vector_utils. IWYU has trouble with operators. ### What kind of change does this PR introduce? - [x] Bug fix (issue: fails to compile when header is used separately, not in BEMC.cc) - [ ] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [ ] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No. ### Does this PR change default behavior? No. --- src/algorithms/calorimetry/ImagingClusterReco.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/algorithms/calorimetry/ImagingClusterReco.h b/src/algorithms/calorimetry/ImagingClusterReco.h index 45c13fce26..df9bce263b 100644 --- a/src/algorithms/calorimetry/ImagingClusterReco.h +++ b/src/algorithms/calorimetry/ImagingClusterReco.h @@ -23,6 +23,7 @@ // Event Model related classes #include #include +#include #include #include #include From 45cfa7ede1afc213b89333f406be5a986ac2eac2 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 3 Jun 2024 17:58:48 -0400 Subject: [PATCH 07/10] JEventProcessorPODIO: deprecate podio:output_include_collections parameter (#1466) This is a follow up on a resolution from #1323 ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [x] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [ ] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No ### Does this PR change default behavior? No --- .github/iwyu.imp | 1 + src/services/io/podio/JEventProcessorPODIO.cc | 37 ++++++++++++++++--- src/services/io/podio/JEventProcessorPODIO.h | 3 +- 3 files changed, 35 insertions(+), 6 deletions(-) diff --git a/.github/iwyu.imp b/.github/iwyu.imp index e7ef6d254c..e1ea3cd578 100644 --- a/.github/iwyu.imp +++ b/.github/iwyu.imp @@ -30,6 +30,7 @@ { include: ['', private, '', public] }, # the rest + { include: ['', private, '', public] }, { include: ['', private, '', public] }, { include: ['', private, '', public] }, # https://github.com/include-what-you-use/include-what-you-use/issues/166 diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 94e9642425..6ebcdd13b0 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -11,7 +11,9 @@ #include #include #include +#include #include +#include #include "services/log/Log_service.h" @@ -43,7 +45,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() { ); // Get the list of output collections to include/exclude - std::vector output_include_collections={ + std::vector output_collections={ // Header and other metadata "EventHeader", @@ -300,9 +302,20 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "DIRCSeededParticleIDs", }; std::vector output_exclude_collections; // need to get as vector, then convert to set + std::string output_include_collections = "DEPRECATED"; japp->SetDefaultParameter( "podio:output_include_collections", output_include_collections, + "DEPRECATED. Use podio:output_collections instead." + ); + if (output_include_collections != "DEPRECATED") { + output_collections.clear(); + JParameterManager::Parse(output_include_collections, output_collections); + m_output_include_collections_set = true; + } + japp->SetDefaultParameter( + "podio:output_collections", + output_collections, "Comma separated list of collection names to write out. If not set, all collections will be written (including ones from input file). Don't set this and use PODIO:OUTPUT_EXCLUDE_COLLECTIONS to write everything except a selection." ); japp->SetDefaultParameter( @@ -316,8 +329,8 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "Comma separated list of collection names to print to screen, e.g. for debugging." ); - m_output_include_collections = std::set(output_include_collections.begin(), - output_include_collections.end()); + m_output_collections = std::set(output_collections.begin(), + output_collections.end()); m_output_exclude_collections = std::set(output_exclude_collections.begin(), output_exclude_collections.end()); @@ -333,6 +346,13 @@ void JEventProcessorPODIO::Init() { // TODO: NWB: Verify that output file is writable NOW, rather than after event processing completes. // I definitely don't trust PODIO to do this for me. + if (m_output_include_collections_set) { + m_log->error("The podio:output_include_collections was provided, but is deprecated. Use podio:output_collections instead."); + // Adding a delay to ensure users notice the deprecation warning. + using namespace std::chrono_literals; + std::this_thread::sleep_for(10s); + } + } @@ -341,7 +361,7 @@ void JEventProcessorPODIO::FindCollectionsToWrite(const std::shared_ptr all_collections = event->GetAllCollectionNames(); - if (m_output_include_collections.empty()) { + if (m_output_collections.empty()) { // User has not specified an include list, so we include _all_ PODIO collections present in the first event. for (const std::string& col : all_collections) { if (m_output_exclude_collections.find(col) == m_output_exclude_collections.end()) { @@ -357,7 +377,7 @@ void JEventProcessorPODIO::FindCollectionsToWrite(const std::shared_ptr all_collections_set = std::set(all_collections.begin(), all_collections.end()); - for (const auto& col : m_output_include_collections) { + for (const auto& col : m_output_collections) { if (m_output_exclude_collections.find(col) == m_output_exclude_collections.end()) { // Included and not excluded if (all_collections_set.find(col) == all_collections_set.end()) { @@ -491,5 +511,12 @@ void JEventProcessorPODIO::Process(const std::shared_ptr &event) { } void JEventProcessorPODIO::Finish() { + if (m_output_include_collections_set) { + m_log->error("The podio:output_include_collections was provided, but is deprecated. Use podio:output_collections instead."); + // Adding a delay to ensure users notice the deprecation warning. + using namespace std::chrono_literals; + std::this_thread::sleep_for(10s); + } + m_writer->finish(); } diff --git a/src/services/io/podio/JEventProcessorPODIO.h b/src/services/io/podio/JEventProcessorPODIO.h index 4ab126a140..77d72679de 100644 --- a/src/services/io/podio/JEventProcessorPODIO.h +++ b/src/services/io/podio/JEventProcessorPODIO.h @@ -30,10 +30,11 @@ class JEventProcessorPODIO : public JEventProcessor { bool m_is_first_event = true; bool m_user_included_collections = false; std::shared_ptr m_log; + bool m_output_include_collections_set = false; std::string m_output_file = "podio_output.root"; std::string m_output_file_copy_dir = ""; - std::set m_output_include_collections; // config. parameter + std::set m_output_collections; // config. parameter std::set m_output_exclude_collections; // config. parameter std::vector m_collections_to_write; // derived from above config. parameters std::vector m_collections_to_print; From 204bfc0126c187b1950841499bc8e498151e67d5 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 3 Jun 2024 18:10:38 -0400 Subject: [PATCH 08/10] Update linux-eic-shell.yml: PIPELINE_NAME_CONTAINER -> PIPELINE_NAME (#1485) https://eicweb.phy.anl.gov/containers/eic_container/-/merge_requests/921 --- .github/workflows/linux-eic-shell.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/linux-eic-shell.yml b/.github/workflows/linux-eic-shell.yml index 16a247ca99..d7c6df253a 100644 --- a/.github/workflows/linux-eic-shell.yml +++ b/.github/workflows/linux-eic-shell.yml @@ -931,7 +931,7 @@ jobs: GITHUB_SHA=${{ github.event.pull_request.head.sha || github.sha }} GITHUB_PR=${{ github.event.pull_request.number }} EICRECON_VERSION="${{ github.event.pull_request.head.ref || github.ref_name }}" - PIPELINE_NAME_CONTAINER=${{ github.repository }}: ${{ github.event.pull_request.title || github.ref_name }} + PIPELINE_NAME=${{ github.repository }}: ${{ github.event.pull_request.title || github.ref_name }} - name: Post a GitHub CI status run: | gh api \ From 2e074ae55d393b572ee13143b9b2c8ef41607d3f Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Mon, 3 Jun 2024 18:05:57 -0500 Subject: [PATCH 09/10] fix: convert to algorithms logger, no m_log in each algorithm class (#1483) ### Briefly, what does this PR introduce? This PR moves a lot of algorithms and factories over to using the internal logger functionality that comes with the algorithms interface. It means that individual algorithms don't need to initialize their own `m_log`. It also means that these algorithms are now properly conformant to the algorithms interface, i.e. `init()` and `process(input, output) const`. ### What kind of change does this PR introduce? - [x] Bug fix (issue #__) - [ ] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [ ] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No. ### Does this PR change default behavior? No. --- src/algorithms/digi/SiliconTrackerDigi.cc | 29 +- src/algorithms/digi/SiliconTrackerDigi.h | 64 ++-- .../fardetectors/FarDetectorLinearTracking.cc | 14 +- .../fardetectors/FarDetectorLinearTracking.h | 69 ++--- .../fardetectors/FarDetectorTrackerCluster.cc | 277 +++++++++--------- .../fardetectors/FarDetectorTrackerCluster.h | 65 ++-- src/algorithms/meta/CollectionCollector.h | 7 +- src/algorithms/onnx/InclusiveKinematicsML.cc | 24 +- src/algorithms/onnx/InclusiveKinematicsML.h | 60 ++-- src/algorithms/pid/MergeTracks.cc | 45 ++- src/algorithms/pid/MergeTracks.h | 42 +-- src/algorithms/reco/HadronicFinalState.cc | 17 +- src/algorithms/reco/HadronicFinalState.h | 45 ++- src/algorithms/reco/InclusiveKinematicsDA.cc | 11 +- src/algorithms/reco/InclusiveKinematicsDA.h | 54 ++-- .../reco/InclusiveKinematicsElectron.cc | 13 +- .../reco/InclusiveKinematicsElectron.h | 54 ++-- src/algorithms/reco/InclusiveKinematicsJB.cc | 13 +- src/algorithms/reco/InclusiveKinematicsJB.h | 54 ++-- .../reco/InclusiveKinematicsSigma.cc | 13 +- .../reco/InclusiveKinematicsSigma.h | 54 ++-- .../reco/InclusiveKinematicsTruth.cc | 12 +- .../reco/InclusiveKinematicsTruth.h | 49 ++-- .../reco/InclusiveKinematicseSigma.cc | 13 +- .../reco/InclusiveKinematicseSigma.h | 54 ++-- src/algorithms/reco/JetReconstruction.cc | 253 ++++++++-------- src/algorithms/reco/JetReconstruction.h | 156 +++++----- src/algorithms/reco/TransformBreitFrame.cc | 26 +- src/algorithms/reco/TransformBreitFrame.h | 54 ++-- .../digi/SiliconTrackerDigi_factory.h | 2 +- .../FarDetectorLinearProjection_factory.h | 2 +- .../FarDetectorLinearTracking_factory.h | 2 +- .../FarDetectorMLReconstruction_factory.h | 2 +- .../FarDetectorTrackerCluster_factory.h | 2 +- .../meta/CollectionCollector_factory.h | 2 +- .../reco/HadronicFinalState_factory.h | 3 +- .../reco/InclusiveKinematicsML_factory.h | 2 +- ...InclusiveKinematicsReconstructed_factory.h | 2 +- .../reco/InclusiveKinematicsTruth_factory.h | 2 +- .../reco/JetReconstruction_factory.h | 2 +- .../reco/TransformBreitFrame_factory.h | 2 +- src/global/pid/MergeTrack_factory.h | 2 +- src/tests/algorithms_test/pid_MergeTracks.cc | 154 +++++----- .../tracking_SiliconSimpleCluster.cc | 163 +++++------ 44 files changed, 882 insertions(+), 1103 deletions(-) diff --git a/src/algorithms/digi/SiliconTrackerDigi.cc b/src/algorithms/digi/SiliconTrackerDigi.cc index 5833f5ae17..b5f21216eb 100644 --- a/src/algorithms/digi/SiliconTrackerDigi.cc +++ b/src/algorithms/digi/SiliconTrackerDigi.cc @@ -20,10 +20,7 @@ namespace eicrecon { -void SiliconTrackerDigi::init(std::shared_ptr& logger) { - // set logger - m_log = logger; - +void SiliconTrackerDigi::init() { // Create random gauss function m_gauss = [&](){ return m_random.Gaus(0, m_cfg.timeResolution); @@ -50,20 +47,20 @@ void SiliconTrackerDigi::process( double result_time = sim_hit.getTime() + time_smearing; auto hit_time_stamp = (std::int32_t) (result_time * 1e3); - m_log->debug("--------------------"); - m_log->debug("Hit cellID = {}", sim_hit.getCellID()); - m_log->debug(" position = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z); - m_log->debug(" xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y)); - m_log->debug(" momentum = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z); - m_log->debug(" edep = {:.2f}", sim_hit.getEDep()); - m_log->debug(" time = {:.4f}[ns]", sim_hit.getTime()); - m_log->debug(" particle time = {}[ns]", sim_hit.getMCParticle().getTime()); - m_log->debug(" time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time); - m_log->debug(" hit_time_stamp: {} [~ps]", hit_time_stamp); + debug("--------------------"); + debug("Hit cellID = {}", sim_hit.getCellID()); + debug(" position = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z); + debug(" xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y)); + debug(" momentum = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z); + debug(" edep = {:.2f}", sim_hit.getEDep()); + debug(" time = {:.4f}[ns]", sim_hit.getTime()); + debug(" particle time = {}[ns]", sim_hit.getMCParticle().getTime()); + debug(" time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time); + debug(" hit_time_stamp: {} [~ps]", hit_time_stamp); if (sim_hit.getEDep() < m_cfg.threshold) { - m_log->debug(" edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV); + debug(" edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV); continue; } @@ -77,7 +74,7 @@ void SiliconTrackerDigi::process( } else { // There is previous values in the cell auto& hit = cell_hit_map[sim_hit.getCellID()]; - m_log->debug(" Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp()); + debug(" Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp()); // keep earliest time for hit auto time_stamp = hit.getTimeStamp(); diff --git a/src/algorithms/digi/SiliconTrackerDigi.h b/src/algorithms/digi/SiliconTrackerDigi.h index 1280783a9f..52afcd44e1 100644 --- a/src/algorithms/digi/SiliconTrackerDigi.h +++ b/src/algorithms/digi/SiliconTrackerDigi.h @@ -5,12 +5,10 @@ #include #include -#include #include +#include #include -#include #include -#include #include #include @@ -19,45 +17,35 @@ namespace eicrecon { - using SiliconTrackerDigiAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::SimTrackerHitCollection - >, - algorithms::Output< - edm4eic::RawTrackerHitCollection, - edm4eic::MCRecoTrackerHitAssociationCollection - > - >; - - class SiliconTrackerDigi - : public SiliconTrackerDigiAlgorithm, - public WithPodConfig { - - public: - SiliconTrackerDigi(std::string_view name) - : SiliconTrackerDigiAlgorithm{name, - {"inputHitCollection"}, - {"outputRawHitCollection","outputHitAssociations"}, - "Apply threshold, digitize within ADC range, " - "convert time with smearing resolution."} {} +using SiliconTrackerDigiAlgorithm = + algorithms::Algorithm, + algorithms::Output>; - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; +class SiliconTrackerDigi : public SiliconTrackerDigiAlgorithm, + public WithPodConfig { - private: - /** algorithm logger */ - std::shared_ptr m_log; +public: + SiliconTrackerDigi(std::string_view name) + : SiliconTrackerDigiAlgorithm{name, + {"inputHitCollection"}, + {"outputRawHitCollection", "outputHitAssociations"}, + "Apply threshold, digitize within ADC range, " + "convert time with smearing resolution."} {} - /** Random number generation*/ - TRandomMixMax m_random; - std::function m_gauss; + void init() final; + void process(const Input&, const Output&) const final; - // FIXME replace with standard random engine - //std::default_random_engine generator; // TODO: need something more appropriate here - //std::normal_distribution m_normDist; // defaults to mean=0, sigma=1 +private: + /** Random number generation*/ + TRandomMixMax m_random; + std::function m_gauss; - //algorithms::Generator m_rng = algorithms::RandomSvc::instance().generator(); + // FIXME replace with standard random engine + // std::default_random_engine generator; // TODO: need something more appropriate here + // std::normal_distribution m_normDist; // defaults to mean=0, sigma=1 - }; + // algorithms::Generator m_rng = algorithms::RandomSvc::instance().generator(); +}; -} // eicrecon +} // namespace eicrecon diff --git a/src/algorithms/fardetectors/FarDetectorLinearTracking.cc b/src/algorithms/fardetectors/FarDetectorLinearTracking.cc index e323ff6db9..5d93d42bce 100644 --- a/src/algorithms/fardetectors/FarDetectorLinearTracking.cc +++ b/src/algorithms/fardetectors/FarDetectorLinearTracking.cc @@ -25,9 +25,7 @@ namespace eicrecon { - void FarDetectorLinearTracking::init(std::shared_ptr& logger) { - - m_log = logger; + void FarDetectorLinearTracking::init() { // For changing how strongly each layer hit is in contributing to the fit m_layerWeights = Eigen::VectorXd::Constant(m_cfg.n_layer,1); @@ -49,7 +47,7 @@ namespace eicrecon { // Check the number of input collections is correct int nCollections = inputhits.size(); if(nCollections!=m_cfg.n_layer){ - m_log->error("Wrong number of input collections passed to algorithm"); + error("Wrong number of input collections passed to algorithm"); return; } @@ -58,7 +56,7 @@ namespace eicrecon { // TODO - Implement more sensible solution for(const auto& layerHits: inputhits){ if((*layerHits).size()>m_cfg.layer_hits_max){ - m_log->info("Too many hits in layer"); + info("Too many hits in layer"); return; } } @@ -160,9 +158,9 @@ namespace eicrecon { double angle = std::acos(hitDiff.dot(m_optimumDirection)); - m_log->debug("Vector: x={}, y={}, z={}",hitDiff.x(),hitDiff.y(),hitDiff.z()); - m_log->debug("Optimum: x={}, y={}, z={}",m_optimumDirection.x(),m_optimumDirection.y(),m_optimumDirection.z()); - m_log->debug("Angle: {}, Tolerance {}",angle,m_cfg.step_angle_tolerance); + debug("Vector: x={}, y={}, z={}",hitDiff.x(),hitDiff.y(),hitDiff.z()); + debug("Optimum: x={}, y={}, z={}",m_optimumDirection.x(),m_optimumDirection.y(),m_optimumDirection.z()); + debug("Angle: {}, Tolerance {}",angle,m_cfg.step_angle_tolerance); if(angle>m_cfg.step_angle_tolerance) return false; diff --git a/src/algorithms/fardetectors/FarDetectorLinearTracking.h b/src/algorithms/fardetectors/FarDetectorLinearTracking.h index 5d4013373c..d9caea23f5 100644 --- a/src/algorithms/fardetectors/FarDetectorLinearTracking.h +++ b/src/algorithms/fardetectors/FarDetectorLinearTracking.h @@ -3,65 +3,54 @@ #pragma once +#include #include #include #include #include #include -#include #include #include #include -#include -#include + #include "FarDetectorLinearTrackingConfig.h" namespace eicrecon { - using FarDetectorLinearTrackingAlgorithm = algorithms::Algorithm< - algorithms::Input< - std::vector - >, - algorithms::Output< - edm4eic::TrackSegmentCollection - > - >; - - class FarDetectorLinearTracking - : public FarDetectorLinearTrackingAlgorithm, - public WithPodConfig { - - public: - FarDetectorLinearTracking(std::string_view name) - : FarDetectorLinearTrackingAlgorithm{name, - {"inputHitCollections"}, - {"outputTrackSegments"}, - "Fit track segments from hits in the tracker layers"} {} +using FarDetectorLinearTrackingAlgorithm = + algorithms::Algorithm>, + algorithms::Output>; - /** One time initialization **/ - void init(std::shared_ptr& logger); +class FarDetectorLinearTracking : public FarDetectorLinearTrackingAlgorithm, + public WithPodConfig { - /** Event by event processing **/ - void process(const Input&, const Output&) const final; +public: + FarDetectorLinearTracking(std::string_view name) + : FarDetectorLinearTrackingAlgorithm{name, + {"inputHitCollections"}, + {"outputTrackSegments"}, + "Fit track segments from hits in the tracker layers"} {} - private: - std::shared_ptr m_log; + /** One time initialization **/ + void init() final; - Eigen::VectorXd m_layerWeights; + /** Event by event processing **/ + void process(const Input&, const Output&) const final; - Eigen::Vector3d m_optimumDirection; +private: + Eigen::VectorXd m_layerWeights; - void buildMatrixRecursive(int level, - Eigen::MatrixXd* hitMatrix, - const std::vector>& hits, - gsl::not_null outputTracks) const; + Eigen::Vector3d m_optimumDirection; - void checkHitCombination(Eigen::MatrixXd* hitMatrix, - gsl::not_null outputTracks) const; + void + buildMatrixRecursive(int level, Eigen::MatrixXd* hitMatrix, + const std::vector>& hits, + gsl::not_null outputTracks) const; - bool checkHitPair(const Eigen::Vector3d& hit1, - const Eigen::Vector3d& hit2) const; + void checkHitCombination(Eigen::MatrixXd* hitMatrix, + gsl::not_null outputTracks) const; - }; + bool checkHitPair(const Eigen::Vector3d& hit1, const Eigen::Vector3d& hit2) const; +}; -} // eicrecon +} // namespace eicrecon diff --git a/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc b/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc index 287dea71e4..b12a61b738 100644 --- a/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc +++ b/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -15,192 +16,186 @@ #include #include #include -#include #include +#include #include "algorithms/fardetectors/FarDetectorTrackerCluster.h" #include "algorithms/fardetectors/FarDetectorTrackerClusterConfig.h" namespace eicrecon { - void FarDetectorTrackerCluster::init(std::shared_ptr& log) { +void FarDetectorTrackerCluster::init() { - m_log = log; - m_detector = algorithms::GeoSvc::instance().detector(); - m_cellid_converter = algorithms::GeoSvc::instance().cellIDPositionConverter(); + m_detector = algorithms::GeoSvc::instance().detector(); + m_cellid_converter = algorithms::GeoSvc::instance().cellIDPositionConverter(); - if (m_cfg.readout.empty()) { - throw JException("Readout is empty"); + if (m_cfg.readout.empty()) { + throw JException("Readout is empty"); + } + try { + m_seg = m_detector->readout(m_cfg.readout).segmentation(); + m_id_dec = m_detector->readout(m_cfg.readout).idSpec().decoder(); + if (!m_cfg.x_field.empty()) { + m_x_idx = m_id_dec->index(m_cfg.x_field); + debug("Find layer field {}, index = {}", m_cfg.x_field, m_x_idx); } - try { - m_seg = m_detector->readout(m_cfg.readout).segmentation(); - m_id_dec = m_detector->readout(m_cfg.readout).idSpec().decoder(); - if (!m_cfg.x_field.empty()) { - m_x_idx = m_id_dec->index(m_cfg.x_field); - m_log->debug("Find layer field {}, index = {}", m_cfg.x_field, m_x_idx); - } - if (!m_cfg.y_field.empty()) { - m_y_idx = m_id_dec->index(m_cfg.y_field); - m_log->debug("Find layer field {}, index = {}", m_cfg.y_field, m_y_idx); - } - } catch (...) { - m_log->error("Failed to load ID decoder for {}", m_cfg.readout); - throw JException("Failed to load ID decoder"); + if (!m_cfg.y_field.empty()) { + m_y_idx = m_id_dec->index(m_cfg.y_field); + debug("Find layer field {}, index = {}", m_cfg.y_field, m_y_idx); } - + } catch (...) { + error("Failed to load ID decoder for {}", m_cfg.readout); + throw JException("Failed to load ID decoder"); } +} - void FarDetectorTrackerCluster::process( - const FarDetectorTrackerCluster::Input& input, - const FarDetectorTrackerCluster::Output& output) const { - - const auto [inputHitsCollections] = input; - auto [outputClustersCollection] = output; +void FarDetectorTrackerCluster::process(const FarDetectorTrackerCluster::Input& input, + const FarDetectorTrackerCluster::Output& output) const { - //Loop over input and output collections - Any collection should only contain hits from a single surface - for(size_t i=0; isize() == 0) continue; - auto outputClusters = outputClustersCollection[i]; + const auto [inputHitsCollections] = input; + auto [outputClustersCollection] = output; - // Make clusters - auto clusters = ClusterHits(*inputHits); + // Loop over input and output collections - Any collection should only contain hits from a single + // surface + for (size_t i = 0; i < inputHitsCollections.size(); i++) { + auto inputHits = inputHitsCollections[i]; + if (inputHits->size() == 0) + continue; + auto outputClusters = outputClustersCollection[i]; - // Create TrackerHits from 2D cluster positions - ConvertClusters(clusters, *outputClusters); + // Make clusters + auto clusters = ClusterHits(*inputHits); - } + // Create TrackerHits from 2D cluster positions + ConvertClusters(clusters, *outputClusters); } +} - // Create vector of FDTrackerCluster from list of hits - std::vector FarDetectorTrackerCluster::ClusterHits(const edm4eic::RawTrackerHitCollection& inputHits) const { - - std::vector clusters; - - ROOT::VecOps::RVec id; - ROOT::VecOps::RVec x; - ROOT::VecOps::RVec y; - ROOT::VecOps::RVec e; - ROOT::VecOps::RVec t; - - // Gather detector id positions - for(const auto& hit: inputHits){ - auto cellID = hit.getCellID(); - id.push_back(cellID); - x.push_back (m_id_dec->get( cellID, m_x_idx )); - y.push_back (m_id_dec->get( cellID, m_y_idx )); - e.push_back (hit.getCharge()); - t.push_back (hit.getTimeStamp()); - } - - // Set up clustering variables - ROOT::VecOps::RVec available(id.size(), 1); - auto indices = Enumerate(id); - - // Loop while there are unclustered hits - while(ROOT::VecOps::Any(available)){ - - dd4hep::Position localPos = {0,0,0}; - float weightSum = 0; - - float esum = 0; - float t0 = 0; - float tError = 0; - auto maxIndex = ROOT::VecOps::ArgMax(e*available); - - available[maxIndex] = 0; - - ROOT::VecOps::RVec clusterList = {maxIndex}; - ROOT::VecOps::RVec clusterT; - std::vector clusterHits; +// Create vector of FDTrackerCluster from list of hits +std::vector +FarDetectorTrackerCluster::ClusterHits(const edm4eic::RawTrackerHitCollection& inputHits) const { + + std::vector clusters; + + ROOT::VecOps::RVec id; + ROOT::VecOps::RVec x; + ROOT::VecOps::RVec y; + ROOT::VecOps::RVec e; + ROOT::VecOps::RVec t; + + // Gather detector id positions + for (const auto& hit : inputHits) { + auto cellID = hit.getCellID(); + id.push_back(cellID); + x.push_back(m_id_dec->get(cellID, m_x_idx)); + y.push_back(m_id_dec->get(cellID, m_y_idx)); + e.push_back(hit.getCharge()); + t.push_back(hit.getTimeStamp()); + } - // Loop over hits, adding neighbouring hits as relevant - while(clusterList.size()){ + // Set up clustering variables + ROOT::VecOps::RVec available(id.size(), 1); + auto indices = Enumerate(id); - // Takes first remaining hit in cluster list - auto index = clusterList[0]; + // Loop while there are unclustered hits + while (ROOT::VecOps::Any(available)) { - // Finds neighbours of cluster within time limit - auto filter = available*(abs(x-x[index])<=1)*(abs(y-y[index])<=1)*(abs(t-t[index]) clusterList = {maxIndex}; + ROOT::VecOps::RVec clusterT; + std::vector clusterHits; - // Adds raw hit to TrackerHit contribution - clusterHits.push_back((inputHits)[index].getObjectID()); + // Loop over hits, adding neighbouring hits as relevant + while (clusterList.size()) { - // Energy - auto hitE = e[index]; - esum += hitE; - // TODO - See if now a single detector element is expected a better function is available. - auto pos = m_seg->position(id[index]); + // Takes first remaining hit in cluster list + auto index = clusterList[0]; - //Weighted position - float weight = hitE; // TODO - Calculate appropriate weighting based on sensor charge sharing - weightSum += weight; - localPos += pos*weight; + // Finds neighbours of cluster within time limit + auto filter = available * (abs(x - x[index]) <= 1) * (abs(y - y[index]) <= 1) * + (abs(t - t[index]) < m_cfg.hit_time_limit); - //Time - clusterT.push_back(t[index]); + // Adds the found hits to the cluster + clusterList = Concatenate(clusterList, indices[filter]); - } + // Removes the found hits from the list of still available hits + available = available * (!filter); - // Finalise position - localPos/=weightSum; + // Removes current hit from remaining found cluster hits + clusterList.erase(clusterList.begin()); - // Finalise time - t0 = Mean(clusterT); - tError = StdDev(clusterT); // TODO fold detector timing resolution into error + // Adds raw hit to TrackerHit contribution + clusterHits.push_back((inputHits)[index].getObjectID()); - // Create cluster - clusters.push_back(FDTrackerCluster{ - .cellID = id[maxIndex], - .x = localPos.x(), - .y = localPos.y(), - .energy = esum, - .time = t0, - .timeError = tError, - .rawHits = clusterHits - }); + // Energy + auto hitE = e[index]; + esum += hitE; + // TODO - See if now a single detector element is expected a better function is available. + auto pos = m_seg->position(id[index]); + // Weighted position + float weight = hitE; // TODO - Calculate appropriate weighting based on sensor charge sharing + weightSum += weight; + localPos += pos * weight; + // Time + clusterT.push_back(t[index]); } - return clusters; - + // Finalise position + localPos /= weightSum; + + // Finalise time + t0 = Mean(clusterT); + tError = StdDev(clusterT); // TODO fold detector timing resolution into error + + // Create cluster + clusters.push_back(FDTrackerCluster{.cellID = id[maxIndex], + .x = localPos.x(), + .y = localPos.y(), + .energy = esum, + .time = t0, + .timeError = tError, + .rawHits = clusterHits}); } - // Convert to global coordinates and create TrackerHits - void FarDetectorTrackerCluster::ConvertClusters(const std::vector& clusters, edm4hep::TrackerHitCollection& outputClusters) const { + return clusters; +} - // Get context of first hit - const dd4hep::VolumeManagerContext* context = m_cellid_converter->findContext(clusters[0].cellID); +// Convert to global coordinates and create TrackerHits +void FarDetectorTrackerCluster::ConvertClusters( + const std::vector& clusters, + edm4hep::TrackerHitCollection& outputClusters) const { - for(auto cluster: clusters){ - auto hitPos = outputClusters.create(); + // Get context of first hit + const dd4hep::VolumeManagerContext* context = m_cellid_converter->findContext(clusters[0].cellID); - auto globalPos = context->localToWorld({cluster.x,cluster.y,0}); + for (auto cluster : clusters) { + auto hitPos = outputClusters.create(); - // Set cluster members - hitPos.setCellID (cluster.cellID); - hitPos.setPosition(edm4hep::Vector3d(globalPos.x()/dd4hep::mm,globalPos.y()/dd4hep::mm,globalPos.z()/dd4hep::mm)); - hitPos.setEDep (cluster.energy); - hitPos.setTime (cluster.time); + auto globalPos = context->localToWorld({cluster.x, cluster.y, 0}); - // Add raw hits to cluster - for(auto hit: cluster.rawHits){ - hitPos.addToRawHits(hit); - } + // Set cluster members + hitPos.setCellID(cluster.cellID); + hitPos.setPosition(edm4hep::Vector3d(globalPos.x() / dd4hep::mm, globalPos.y() / dd4hep::mm, + globalPos.z() / dd4hep::mm)); + hitPos.setEDep(cluster.energy); + hitPos.setTime(cluster.time); + // Add raw hits to cluster + for (auto hit : cluster.rawHits) { + hitPos.addToRawHits(hit); } - } - - } + +} // namespace eicrecon diff --git a/src/algorithms/fardetectors/FarDetectorTrackerCluster.h b/src/algorithms/fardetectors/FarDetectorTrackerCluster.h index 381eb97886..6ac60e5067 100644 --- a/src/algorithms/fardetectors/FarDetectorTrackerCluster.h +++ b/src/algorithms/fardetectors/FarDetectorTrackerCluster.h @@ -11,8 +11,6 @@ #include #include #include -#include -#include #include #include #include @@ -32,48 +30,41 @@ struct FDTrackerCluster { }; namespace eicrecon { - using FarDetectorTrackerClusterAlgorithm = algorithms::Algorithm< - algorithms::Input< - std::vector - >, - algorithms::Output< - std::vector - > - >; +using FarDetectorTrackerClusterAlgorithm = + algorithms::Algorithm>, + algorithms::Output>>; - class FarDetectorTrackerCluster - : public FarDetectorTrackerClusterAlgorithm, - public WithPodConfig { +class FarDetectorTrackerCluster : public FarDetectorTrackerClusterAlgorithm, + public WithPodConfig { - public: - FarDetectorTrackerCluster(std::string_view name) +public: + FarDetectorTrackerCluster(std::string_view name) : FarDetectorTrackerClusterAlgorithm{name, - {"inputHitCollection"}, - {"outputClusterPositionCollection"}, - "Simple weighted clustering of hits by x-y component of single detector element segmentation"} {} + {"inputHitCollection"}, + {"outputClusterPositionCollection"}, + "Simple weighted clustering of hits by x-y component of " + "single detector element segmentation"} {} - /** One time initialization **/ - void init(std::shared_ptr& logger); + /** One time initialization **/ + void init() final; - /** Event by event processing **/ - void process(const Input&, const Output&) const final; + /** Event by event processing **/ + void process(const Input&, const Output&) const final; - /** Cluster hits **/ - std::vector ClusterHits(const edm4eic::RawTrackerHitCollection&) const; + /** Cluster hits **/ + std::vector ClusterHits(const edm4eic::RawTrackerHitCollection&) const; - /** Convert clusters to TrackerHits **/ - void ConvertClusters(const std::vector&, edm4hep::TrackerHitCollection&) const; + /** Convert clusters to TrackerHits **/ + void ConvertClusters(const std::vector&, edm4hep::TrackerHitCollection&) const; - private: - const dd4hep::Detector* m_detector{nullptr}; - const dd4hep::BitFieldCoder* m_id_dec{nullptr}; - std::shared_ptr m_log; - const dd4hep::rec::CellIDPositionConverter* m_cellid_converter{nullptr}; - dd4hep::Segmentation m_seg; +private: + const dd4hep::Detector* m_detector{nullptr}; + const dd4hep::BitFieldCoder* m_id_dec{nullptr}; + const dd4hep::rec::CellIDPositionConverter* m_cellid_converter{nullptr}; + dd4hep::Segmentation m_seg; - int m_x_idx{0}; - int m_y_idx{0}; - - }; + int m_x_idx{0}; + int m_y_idx{0}; +}; -} // eicrecon +} // namespace eicrecon diff --git a/src/algorithms/meta/CollectionCollector.h b/src/algorithms/meta/CollectionCollector.h index accf942f42..f6806366b8 100644 --- a/src/algorithms/meta/CollectionCollector.h +++ b/src/algorithms/meta/CollectionCollector.h @@ -28,9 +28,7 @@ namespace eicrecon { "Merge content of collections into one subset collection" }{} - void init(std::shared_ptr& logger){ // set logger - m_log = logger; - }; + void init() final { }; void process(const typename CollectionCollector::Input& input, const typename CollectionCollector::Output& output) const final{ @@ -46,8 +44,5 @@ namespace eicrecon { } } - private: - std::shared_ptr m_log; // logger - }; } // eicrecon diff --git a/src/algorithms/onnx/InclusiveKinematicsML.cc b/src/algorithms/onnx/InclusiveKinematicsML.cc index 6ef2e72a40..935efadb77 100644 --- a/src/algorithms/onnx/InclusiveKinematicsML.cc +++ b/src/algorithms/onnx/InclusiveKinematicsML.cc @@ -30,9 +30,7 @@ namespace eicrecon { return tensor; } - void InclusiveKinematicsML::init(std::shared_ptr& logger) { - m_log = logger; - + void InclusiveKinematicsML::init() { // onnxruntime setup Ort::Env env(ORT_LOGGING_LEVEL_WARNING, "inclusive-kinematics-ml"); Ort::SessionOptions session_options; @@ -41,19 +39,19 @@ namespace eicrecon { // print name/shape of inputs Ort::AllocatorWithDefaultOptions allocator; - m_log->debug("Input Node Name/Shape:"); + debug("Input Node Name/Shape:"); for (std::size_t i = 0; i < m_session.GetInputCount(); i++) { m_input_names.emplace_back(m_session.GetInputNameAllocated(i, allocator).get()); m_input_shapes.emplace_back(m_session.GetInputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); - m_log->debug("\t{} : {}", m_input_names.at(i), print_shape(m_input_shapes.at(i))); + debug("\t{} : {}", m_input_names.at(i), print_shape(m_input_shapes.at(i))); } // print name/shape of outputs - m_log->debug("Output Node Name/Shape:"); + debug("Output Node Name/Shape:"); for (std::size_t i = 0; i < m_session.GetOutputCount(); i++) { m_output_names.emplace_back(m_session.GetOutputNameAllocated(i, allocator).get()); m_output_shapes.emplace_back(m_session.GetOutputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); - m_log->debug("\t{} : {}", m_output_names.at(i), print_shape(m_output_shapes.at(i))); + debug("\t{} : {}", m_output_names.at(i), print_shape(m_output_shapes.at(i))); } // convert names to char* @@ -65,7 +63,7 @@ namespace eicrecon { [&](const std::string& str) { return str.c_str(); }); } catch(std::exception& e) { - m_log->error(e.what()); + error(e.what()); } } @@ -78,13 +76,13 @@ namespace eicrecon { // Require valid inputs if (electron->size() == 0 || da->size() == 0) { - m_log->debug("skipping because input collections have no entries"); + debug("skipping because input collections have no entries"); return; } // Assume model has 1 input nodes and 1 output node. if (m_input_names.size() != 1 || m_output_names.size() != 1) { - m_log->debug("skipping because model has incorrect input and output size"); + debug("skipping because model has incorrect input and output size"); return; } @@ -98,7 +96,7 @@ namespace eicrecon { // Double-check the dimensions of the input tensor if (! input_tensors[0].IsTensor() || input_tensors[0].GetTensorTypeAndShapeInfo().GetShape() != m_input_shapes.front()) { - m_log->debug("skipping because input tensor shape incorrect"); + debug("skipping because input tensor shape incorrect"); return; } @@ -109,7 +107,7 @@ namespace eicrecon { // Double-check the dimensions of the output tensors if (!output_tensors[0].IsTensor() || output_tensors.size() != m_output_names.size()) { - m_log->debug("skipping because output tensor shape incorrect"); + debug("skipping because output tensor shape incorrect"); return; } @@ -120,7 +118,7 @@ namespace eicrecon { kin.setX(x); } catch (const Ort::Exception& exception) { - m_log->error("error running model inference: {}", exception.what()); + error("error running model inference: {}", exception.what()); } } diff --git a/src/algorithms/onnx/InclusiveKinematicsML.h b/src/algorithms/onnx/InclusiveKinematicsML.h index e2d559dbf0..485fdb5d8b 100644 --- a/src/algorithms/onnx/InclusiveKinematicsML.h +++ b/src/algorithms/onnx/InclusiveKinematicsML.h @@ -4,11 +4,9 @@ #pragma once #include +#include #include #include -#include -#include -#include #include #include #include @@ -18,43 +16,35 @@ namespace eicrecon { - using InclusiveKinematicsMLAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4eic::InclusiveKinematicsCollection, - edm4eic::InclusiveKinematicsCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsML - : public InclusiveKinematicsMLAlgorithm, - public WithPodConfig { - - public: - InclusiveKinematicsML(std::string_view name) - : InclusiveKinematicsMLAlgorithm{name, - {"inclusiveKinematicsElectron", "inclusiveKinematicsDA"}, - {"inclusiveKinematicsML"}, - "Determine inclusive kinematics using combined ML method."} {} +using InclusiveKinematicsMLAlgorithm = + algorithms::Algorithm, + algorithms::Output>; - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; +class InclusiveKinematicsML : public InclusiveKinematicsMLAlgorithm, + public WithPodConfig { - private: - std::shared_ptr m_log; +public: + InclusiveKinematicsML(std::string_view name) + : InclusiveKinematicsMLAlgorithm{name, + {"inclusiveKinematicsElectron", "inclusiveKinematicsDA"}, + {"inclusiveKinematicsML"}, + "Determine inclusive kinematics using combined ML method."} { + } - mutable Ort::Session m_session{nullptr}; + void init() final; + void process(const Input&, const Output&) const final; - std::vector m_input_names; - std::vector m_input_names_char; - std::vector> m_input_shapes; +private: + mutable Ort::Session m_session{nullptr}; - std::vector m_output_names; - std::vector m_output_names_char; - std::vector> m_output_shapes; + std::vector m_input_names; + std::vector m_input_names_char; + std::vector> m_input_shapes; - }; + std::vector m_output_names; + std::vector m_output_names_char; + std::vector> m_output_shapes; +}; } // namespace eicrecon diff --git a/src/algorithms/pid/MergeTracks.cc b/src/algorithms/pid/MergeTracks.cc index 41087a91fb..847994332f 100644 --- a/src/algorithms/pid/MergeTracks.cc +++ b/src/algorithms/pid/MergeTracks.cc @@ -3,71 +3,64 @@ #include "MergeTracks.h" +#include +#include #include #include #include -#include -#include -#include +#include #include +#include #include #include namespace eicrecon { -void MergeTracks::init(std::shared_ptr& logger) -{ - m_log = logger; -} - -void MergeTracks::process( - const MergeTracks::Input& input, - const MergeTracks::Output& output) const { +void MergeTracks::process(const MergeTracks::Input& input, + const MergeTracks::Output& output) const { const auto [in_track_collections] = input; - auto [out_tracks] = output; + auto [out_tracks] = output; // logging - m_log->trace("{:=^70}"," call MergeTracks::AlgorithmProcess "); + trace("{:=^70}", " call MergeTracks::AlgorithmProcess "); // check that all input collections have the same size std::unordered_map in_track_collection_size_distribution; - for(const auto& in_track_collection : in_track_collections) { + for (const auto& in_track_collection : in_track_collections) { ++in_track_collection_size_distribution[in_track_collection->size()]; } if (in_track_collection_size_distribution.size() != 1) { std::vector in_track_collection_sizes; std::transform(in_track_collections.begin(), in_track_collections.end(), - std::back_inserter(in_track_collection_sizes), - [](const auto& in_track_collection) { return in_track_collection->size(); }); - m_log->error("cannot merge input track collections with different sizes {}", fmt::join(in_track_collection_sizes, ", ")); + std::back_inserter(in_track_collection_sizes), + [](const auto& in_track_collection) { return in_track_collection->size(); }); + error("cannot merge input track collections with different sizes {}", + fmt::join(in_track_collection_sizes, ", ")); return; } // loop over track collection elements std::size_t n_tracks = in_track_collection_size_distribution.begin()->first; - for(std::size_t i_track=0; i_trackcreate(); std::vector out_track_points; // loop over collections for this track, and add each track's points to `out_track_points` - for(const auto& in_track_collection : in_track_collections) { + for (const auto& in_track_collection : in_track_collections) { const auto& in_track = in_track_collection->at(i_track); - for(const auto& point : in_track.getPoints()) + for (const auto& point : in_track.getPoints()) out_track_points.push_back(point); } // sort all `out_track_points` by time - std::sort( - out_track_points.begin(), - out_track_points.end(), - [] (edm4eic::TrackPoint& a, edm4eic::TrackPoint& b) { return a.time < b.time; } - ); + std::sort(out_track_points.begin(), out_track_points.end(), + [](edm4eic::TrackPoint& a, edm4eic::TrackPoint& b) { return a.time < b.time; }); // add these sorted points to `out_track` - for(const auto& point : out_track_points) + for (const auto& point : out_track_points) out_track.addToPoints(point); /* FIXME: merge other members, such as `length` and `lengthError`; diff --git a/src/algorithms/pid/MergeTracks.h b/src/algorithms/pid/MergeTracks.h index abb5520e13..853ec1eb1c 100644 --- a/src/algorithms/pid/MergeTracks.h +++ b/src/algorithms/pid/MergeTracks.h @@ -12,36 +12,26 @@ // data model #include #include -#include -#include +#include +#include #include namespace eicrecon { - using MergeTracksAlgorithm = algorithms::Algorithm< - algorithms::Input< - std::vector - >, - algorithms::Output< - edm4eic::TrackSegmentCollection - > - >; +using MergeTracksAlgorithm = + algorithms::Algorithm>, + algorithms::Output>; - class MergeTracks - : public MergeTracksAlgorithm { +class MergeTracks : public MergeTracksAlgorithm { - public: - MergeTracks(std::string_view name) +public: + MergeTracks(std::string_view name) : MergeTracksAlgorithm{name, - {"inputTrackSegments"}, - {"outputTrackSegments"}, - "Effectively 'zip' the input track segments."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - - }; -} + {"inputTrackSegments"}, + {"outputTrackSegments"}, + "Effectively 'zip' the input track segments."} {} + + void init() final{}; + void process(const Input&, const Output&) const final; +}; +} // namespace eicrecon diff --git a/src/algorithms/reco/HadronicFinalState.cc b/src/algorithms/reco/HadronicFinalState.cc index 9cc84121b9..016770a85b 100644 --- a/src/algorithms/reco/HadronicFinalState.cc +++ b/src/algorithms/reco/HadronicFinalState.cc @@ -25,11 +25,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void HadronicFinalState::init(std::shared_ptr& logger) { - m_log = logger; + void HadronicFinalState::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -44,7 +43,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -58,7 +57,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -72,7 +71,7 @@ namespace eicrecon { // Get first scattered electron const auto ef_coll = find_first_scattered_electron(mcparts); if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); + debug("No truth scattered electron found"); return; } // Associate first scattered electron with reconstructed electrons @@ -87,7 +86,7 @@ namespace eicrecon { } } if (!(ef_assoc != rcassoc->end())) { - m_log->debug("Truth scattered electron not in reconstructed particles"); + debug("Truth scattered electron not in reconstructed particles"); return; } const auto ef_rc{ef_assoc->getRec()}; @@ -133,11 +132,11 @@ namespace eicrecon { // Sigma zero or negative if (sigma <= 0) { - m_log->debug("Sigma zero or negative"); + debug("Sigma zero or negative"); return; } - m_log->debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma()); + debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma()); } diff --git a/src/algorithms/reco/HadronicFinalState.h b/src/algorithms/reco/HadronicFinalState.h index e4de62ebc9..6b25f48d7f 100644 --- a/src/algorithms/reco/HadronicFinalState.h +++ b/src/algorithms/reco/HadronicFinalState.h @@ -8,41 +8,30 @@ #include #include #include -#include -#include #include #include - namespace eicrecon { - using HadronicFinalStateAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::MCRecoParticleAssociationCollection - >, - algorithms::Output< - edm4eic::HadronicFinalStateCollection - > - >; - - class HadronicFinalState - : public HadronicFinalStateAlgorithm { - - public: - HadronicFinalState(std::string_view name) +using HadronicFinalStateAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class HadronicFinalState : public HadronicFinalStateAlgorithm { + +public: + HadronicFinalState(std::string_view name) : HadronicFinalStateAlgorithm{name, - {"MCParticles", "inputParticles", "inputAssociations"}, - {"hadronicFinalState"}, - "Calculate summed quantities of the hadronic final state."} {} + {"MCParticles", "inputParticles", "inputAssociations"}, + {"hadronicFinalState"}, + "Calculate summed quantities of the hadronic final state."} {} - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; + void init() final; + void process(const Input&, const Output&) const final; - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsDA.cc b/src/algorithms/reco/InclusiveKinematicsDA.cc index dd58f36a45..d925f07f01 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.cc +++ b/src/algorithms/reco/InclusiveKinematicsDA.cc @@ -23,8 +23,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsDA::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsDA::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { // m_log->debug("Unable to locate Particle Service. " @@ -42,7 +41,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -56,7 +55,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -83,7 +82,7 @@ namespace eicrecon { // Sigma zero or negative if (sigma_h <= 0) { - m_log->debug("Sigma zero or negative"); + debug("Sigma zero or negative"); return; } @@ -96,7 +95,7 @@ namespace eicrecon { auto kin = kinematics->create(x_da, Q2_da, W_da, y_da, nu_da); kin.setScat(kf); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsDA.h b/src/algorithms/reco/InclusiveKinematicsDA.h index c1d3f72468..ab6b71f7de 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.h +++ b/src/algorithms/reco/InclusiveKinematicsDA.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsDAAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsDA - : public InclusiveKinematicsDAAlgorithm { - - public: - InclusiveKinematicsDA(std::string_view name) - : InclusiveKinematicsDAAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using double-angle method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsDAAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsDA : public InclusiveKinematicsDAAlgorithm { + +public: + InclusiveKinematicsDA(std::string_view name) + : InclusiveKinematicsDAAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using double-angle method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.cc b/src/algorithms/reco/InclusiveKinematicsElectron.cc index 0a4781db91..8ac014f2ff 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.cc +++ b/src/algorithms/reco/InclusiveKinematicsElectron.cc @@ -23,11 +23,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsElectron::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsElectron::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -90,7 +89,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -104,7 +103,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -124,7 +123,7 @@ namespace eicrecon { // If no scattered electron was found if (electrons.size() == 0) { - m_log->debug("No scattered electron found"); + debug("No scattered electron found"); return; } @@ -140,7 +139,7 @@ namespace eicrecon { auto kin = kinematics->create(x, Q2, W, y, nu); kin.setScat(escat->at(0)); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.h b/src/algorithms/reco/InclusiveKinematicsElectron.h index e95cc55138..8b4864d14a 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.h +++ b/src/algorithms/reco/InclusiveKinematicsElectron.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsElectronAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsElectron - : public InclusiveKinematicsElectronAlgorithm { - - public: - InclusiveKinematicsElectron(std::string_view name) - : InclusiveKinematicsElectronAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using electron method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsElectronAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsElectron : public InclusiveKinematicsElectronAlgorithm { + +public: + InclusiveKinematicsElectron(std::string_view name) + : InclusiveKinematicsElectronAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using electron method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsJB.cc b/src/algorithms/reco/InclusiveKinematicsJB.cc index 324093b1aa..48874d79bd 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.cc +++ b/src/algorithms/reco/InclusiveKinematicsJB.cc @@ -20,11 +20,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsJB::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsJB::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -39,7 +38,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -53,7 +52,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -71,7 +70,7 @@ namespace eicrecon { // Sigma zero or negative if (sigma_h <= 0) { - m_log->debug("Sigma zero or negative"); + debug("Sigma zero or negative"); return; } @@ -84,7 +83,7 @@ namespace eicrecon { auto kin = kinematics->create(x_jb, Q2_jb, W_jb, y_jb, nu_jb); kin.setScat(escat->at(0)); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsJB.h b/src/algorithms/reco/InclusiveKinematicsJB.h index 24526c04b2..bacc5a69f7 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.h +++ b/src/algorithms/reco/InclusiveKinematicsJB.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsJBAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsJB - : public InclusiveKinematicsJBAlgorithm { - - public: - InclusiveKinematicsJB(std::string_view name) - : InclusiveKinematicsJBAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using Jacquet-Blondel method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsJBAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsJB : public InclusiveKinematicsJBAlgorithm { + +public: + InclusiveKinematicsJB(std::string_view name) + : InclusiveKinematicsJBAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using Jacquet-Blondel method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.cc b/src/algorithms/reco/InclusiveKinematicsSigma.cc index 020d1dc5fb..c63b8f3981 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicsSigma.cc @@ -22,11 +22,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsSigma::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsSigma::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -41,7 +40,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -55,7 +54,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -82,7 +81,7 @@ namespace eicrecon { auto gamma_h = hfs->at(0).getGamma(); if (sigma_h <= 0) { - m_log->debug("No scattered electron found or sigma zero or negative"); + debug("No scattered electron found or sigma zero or negative"); return; } @@ -97,7 +96,7 @@ namespace eicrecon { auto kin = kinematics->create(x_sig, Q2_sig, W_sig, y_sig, nu_sig); kin.setScat(kf); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.h b/src/algorithms/reco/InclusiveKinematicsSigma.h index d3849806c3..6cf3d13cbe 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.h +++ b/src/algorithms/reco/InclusiveKinematicsSigma.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsSigmaAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsSigma - : public InclusiveKinematicsSigmaAlgorithm { - - public: - InclusiveKinematicsSigma(std::string_view name) - : InclusiveKinematicsSigmaAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using Sigma method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsSigmaAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsSigma : public InclusiveKinematicsSigmaAlgorithm { + +public: + InclusiveKinematicsSigma(std::string_view name) + : InclusiveKinematicsSigmaAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using Sigma method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.cc b/src/algorithms/reco/InclusiveKinematicsTruth.cc index 54bc232f88..e48bbc0d77 100644 --- a/src/algorithms/reco/InclusiveKinematicsTruth.cc +++ b/src/algorithms/reco/InclusiveKinematicsTruth.cc @@ -19,9 +19,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsTruth::init(std::shared_ptr& logger) { - m_log = logger; - + void InclusiveKinematicsTruth::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { // error() << "Unable to locate Particle Service. " @@ -47,7 +45,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const auto ei_p = ei_coll[0].getMomentum(); @@ -58,7 +56,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const auto pi_p = pi_coll[0].getMomentum(); @@ -73,7 +71,7 @@ namespace eicrecon { // the beam. const auto ef_coll = find_first_scattered_electron(mcparts); if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); + debug("No truth scattered electron found"); return; } const auto ef_p = ef_coll[0].getMomentum(); @@ -91,7 +89,7 @@ namespace eicrecon { const auto W = sqrt(pi_mass*pi_mass + 2.*q_dot_pi - Q2); auto kin = kinematics->create(x, Q2, W, y, nu); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.h b/src/algorithms/reco/InclusiveKinematicsTruth.h index 8ec7a58cf7..83a3f44576 100644 --- a/src/algorithms/reco/InclusiveKinematicsTruth.h +++ b/src/algorithms/reco/InclusiveKinematicsTruth.h @@ -6,39 +6,30 @@ #include #include #include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsTruthAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsTruth - : public InclusiveKinematicsTruthAlgorithm { - - public: - InclusiveKinematicsTruth(std::string_view name) - : InclusiveKinematicsTruthAlgorithm{name, - {"MCParticles"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics from truth information."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsTruthAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +class InclusiveKinematicsTruth : public InclusiveKinematicsTruthAlgorithm { + +public: + InclusiveKinematicsTruth(std::string_view name) + : InclusiveKinematicsTruthAlgorithm{ + name, + {"MCParticles"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics from truth information."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.cc b/src/algorithms/reco/InclusiveKinematicseSigma.cc index b898e41da9..56fb00eec1 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicseSigma.cc @@ -23,11 +23,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicseSigma::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicseSigma::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -42,7 +41,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -56,7 +55,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -83,7 +82,7 @@ namespace eicrecon { auto gamma_h = hfs->at(0).getGamma(); if (sigma_h <= 0) { - m_log->debug("No scattered electron found or sigma zero or negative"); + debug("No scattered electron found or sigma zero or negative"); return; } @@ -106,7 +105,7 @@ namespace eicrecon { kin.setScat(kf); // Debugging output - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.h b/src/algorithms/reco/InclusiveKinematicseSigma.h index 9a505f0eae..030a37453d 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.h +++ b/src/algorithms/reco/InclusiveKinematicseSigma.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicseSigmaAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicseSigma - : public InclusiveKinematicseSigmaAlgorithm { - - public: - InclusiveKinematicseSigma(std::string_view name) - : InclusiveKinematicseSigmaAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using e-Sigma method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicseSigmaAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicseSigma : public InclusiveKinematicseSigmaAlgorithm { + +public: + InclusiveKinematicseSigma(std::string_view name) + : InclusiveKinematicseSigmaAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using e-Sigma method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/JetReconstruction.cc b/src/algorithms/reco/JetReconstruction.cc index ebe60f55aa..09dad9e089 100644 --- a/src/algorithms/reco/JetReconstruction.cc +++ b/src/algorithms/reco/JetReconstruction.cc @@ -6,9 +6,10 @@ // for error handling #include -#include // IWYU pragma: keep +#include // IWYU pragma: keep #include #include +#include #include // for fastjet objects #include @@ -24,139 +25,141 @@ using namespace fastjet; namespace eicrecon { - template - void JetReconstruction::init(std::shared_ptr logger) { - - m_log = logger; - m_log->trace("Initialized"); - - // if specified algorithm, recomb. scheme, or area type - // are not defined, then issue error and throw exception - try { - m_mapJetAlgo.at(m_cfg.jetAlgo); - } catch (std::out_of_range &out) { - m_log->error(" Unknown jet algorithm \"{}\" specified!", m_cfg.jetAlgo); - throw JException(out.what()); - } - - try { - m_mapRecombScheme.at(m_cfg.recombScheme); - } catch (std::out_of_range &out) { - m_log->error(" Unknown recombination scheme \"{}\" specified!", m_cfg.recombScheme); - throw JException(out.what()); - } - - try { - m_mapAreaType.at(m_cfg.areaType); - } catch (std::out_of_range &out) { - m_log->error(" Unknown area type \"{}\" specified!", m_cfg.areaType); - throw JException(out.what()); +template void JetReconstruction::init() { + + this->trace("Initialized"); + + // if specified algorithm, recomb. scheme, or area type + // are not defined, then issue error and throw exception + try { + m_mapJetAlgo.at(m_cfg.jetAlgo); + } catch (std::out_of_range& out) { + this->error(" Unknown jet algorithm \"{}\" specified!", m_cfg.jetAlgo); + throw JException(out.what()); + } + + try { + m_mapRecombScheme.at(m_cfg.recombScheme); + } catch (std::out_of_range& out) { + this->error(" Unknown recombination scheme \"{}\" specified!", m_cfg.recombScheme); + throw JException(out.what()); + } + + try { + m_mapAreaType.at(m_cfg.areaType); + } catch (std::out_of_range& out) { + this->error(" Unknown area type \"{}\" specified!", m_cfg.areaType); + throw JException(out.what()); + } + + // Choose jet definition based on no. of parameters + switch (m_mapJetAlgo[m_cfg.jetAlgo]) { + + // contributed algorithms + case JetAlgorithm::plugin_algorithm: + + // expand to other algorithms as required + if (m_cfg.jetContribAlgo == "Centauro") { + m_jet_plugin = std::make_unique(m_cfg.rJet); + m_jet_def = std::make_unique(m_jet_plugin.get()); + } else { + this->error(" Unknown contributed FastJet algorithm \"{}\" specified!", m_cfg.jetContribAlgo); + throw JException("Invalid contributed FastJet algorithm"); } - - // Choose jet definition based on no. of parameters - switch (m_mapJetAlgo[m_cfg.jetAlgo]) { - - // contributed algorithms - case JetAlgorithm::plugin_algorithm: - - // expand to other algorithms as required - if(m_cfg.jetContribAlgo == "Centauro"){ - m_jet_plugin = std::make_unique(m_cfg.rJet); - m_jet_def = std::make_unique(m_jet_plugin.get()); - } - else { - m_log->error(" Unknown contributed FastJet algorithm \"{}\" specified!", m_cfg.jetContribAlgo); - throw JException("Invalid contributed FastJet algorithm"); - } - break; - - // 0 parameter algorithms - case JetAlgorithm::ee_kt_algorithm: - m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_mapRecombScheme[m_cfg.recombScheme]); - break; - - // 2 parameter algorithms - case JetAlgorithm::genkt_algorithm: - [[fallthrough]]; - - case JetAlgorithm::ee_genkt_algorithm: - m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, m_cfg.pJet, m_mapRecombScheme[m_cfg.recombScheme]); - break; - - // all others have only 1 parameter - default: - m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, m_mapRecombScheme[m_cfg.recombScheme]); - break; - - } // end switch (jet algorithm) - - // Define jet area - m_area_def = std::make_unique(m_mapAreaType[m_cfg.areaType], GhostedAreaSpec(m_cfg.ghostMaxRap, m_cfg.numGhostRepeat, m_cfg.ghostArea)); - - } // end 'init(std::shared_ptr)' - - template - void JetReconstruction::process( - const typename JetReconstructionAlgorithm::Input& input, - const typename JetReconstructionAlgorithm::Output& output - ) const { - // Grab input collections - const auto [input_collection] = input; - auto [jet_collection] = output; - - // extract input momenta and collect into pseudojets - std::vector particles; - for (unsigned iInput = 0; const auto& input : *input_collection) { - - // get 4-vector - const auto& momentum = input.getMomentum(); - const auto& energy = input.getEnergy(); - const auto pt = edm4hep::utils::magnitudeTransverse(momentum); - - // Only cluster particles within the given pt Range - if ((pt > m_cfg.minCstPt) && (pt < m_cfg.maxCstPt)) { - particles.emplace_back(momentum.x, momentum.y, momentum.z, energy); - particles.back().set_user_index(iInput); - } - ++iInput; + break; + + // 0 parameter algorithms + case JetAlgorithm::ee_kt_algorithm: + m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], + m_mapRecombScheme[m_cfg.recombScheme]); + break; + + // 2 parameter algorithms + case JetAlgorithm::genkt_algorithm: + [[fallthrough]]; + + case JetAlgorithm::ee_genkt_algorithm: + m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, m_cfg.pJet, + m_mapRecombScheme[m_cfg.recombScheme]); + break; + + // all others have only 1 parameter + default: + m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, + m_mapRecombScheme[m_cfg.recombScheme]); + break; + + } // end switch (jet algorithm) + + // Define jet area + m_area_def = std::make_unique( + m_mapAreaType[m_cfg.areaType], + GhostedAreaSpec(m_cfg.ghostMaxRap, m_cfg.numGhostRepeat, m_cfg.ghostArea)); + +} // end 'init()' + +template +void JetReconstruction::process( + const typename JetReconstructionAlgorithm::Input& input, + const typename JetReconstructionAlgorithm::Output& output) const { + // Grab input collections + const auto [input_collection] = input; + auto [jet_collection] = output; + + // extract input momenta and collect into pseudojets + std::vector particles; + for (unsigned iInput = 0; const auto& input : *input_collection) { + + // get 4-vector + const auto& momentum = input.getMomentum(); + const auto& energy = input.getEnergy(); + const auto pt = edm4hep::utils::magnitudeTransverse(momentum); + + // Only cluster particles within the given pt Range + if ((pt > m_cfg.minCstPt) && (pt < m_cfg.maxCstPt)) { + particles.emplace_back(momentum.x, momentum.y, momentum.z, energy); + particles.back().set_user_index(iInput); } + ++iInput; + } - // Skip empty - if (particles.empty()) { - m_log->trace(" Empty particle list."); - return; - } - m_log->trace(" Number of particles: {}", particles.size()); + // Skip empty + if (particles.empty()) { + this->trace(" Empty particle list."); + return; + } + this->trace(" Number of particles: {}", particles.size()); - // Run the clustering, extract the jets - fastjet::ClusterSequenceArea m_clus_seq(particles, *m_jet_def, *m_area_def); - std::vector jets = sorted_by_pt(m_clus_seq.inclusive_jets(m_cfg.minJetPt)); + // Run the clustering, extract the jets + fastjet::ClusterSequenceArea m_clus_seq(particles, *m_jet_def, *m_area_def); + std::vector jets = sorted_by_pt(m_clus_seq.inclusive_jets(m_cfg.minJetPt)); - // Print out some infos - m_log->trace(" Clustering with : {}", m_jet_def->description()); + // Print out some infos + this->trace(" Clustering with : {}", m_jet_def->description()); - // loop over jets - for (unsigned i = 0; i < jets.size(); i++) { + // loop over jets + for (unsigned i = 0; i < jets.size(); i++) { - m_log->trace(" jet {}: pt = {}, y = {}, phi = {}", i, jets[i].pt(), jets[i].rap(), jets[i].phi()); + this->trace(" jet {}: pt = {}, y = {}, phi = {}", i, jets[i].pt(), jets[i].rap(), + jets[i].phi()); - // create jet to store in output collection - edm4eic::MutableReconstructedParticle jet_output = jet_collection->create(); - jet_output.setMomentum(edm4hep::Vector3f(jets[i].px(), jets[i].py(), jets[i].pz())); - jet_output.setEnergy(jets[i].e()); - jet_output.setMass(jets[i].m()); + // create jet to store in output collection + edm4eic::MutableReconstructedParticle jet_output = jet_collection->create(); + jet_output.setMomentum(edm4hep::Vector3f(jets[i].px(), jets[i].py(), jets[i].pz())); + jet_output.setEnergy(jets[i].e()); + jet_output.setMass(jets[i].m()); - // link constituents to jet kinematic info - std::vector csts = jets[i].constituents(); - for (unsigned j = 0; j < csts.size(); j++) { - jet_output.addToParticles(input_collection->at(csts[j].user_index())); - } // for constituent j - } // for jet i + // link constituents to jet kinematic info + std::vector csts = jets[i].constituents(); + for (unsigned j = 0; j < csts.size(); j++) { + jet_output.addToParticles(input_collection->at(csts[j].user_index())); + } // for constituent j + } // for jet i - // return the jets - return; - } // end 'process(const T&)' + // return the jets + return; +} // end 'process(const T&)' - template class JetReconstruction; +template class JetReconstruction; -} // end namespace eicrecon +} // end namespace eicrecon diff --git a/src/algorithms/reco/JetReconstruction.h b/src/algorithms/reco/JetReconstruction.h index 460aa1d33d..9b2182479d 100644 --- a/src/algorithms/reco/JetReconstruction.h +++ b/src/algorithms/reco/JetReconstruction.h @@ -6,9 +6,7 @@ #include #include #include -#include #include -#include #include #include #include @@ -20,87 +18,73 @@ namespace eicrecon { - template - using JetReconstructionAlgorithm = algorithms::Algorithm< - algorithms::Input< - typename InputT::collection_type - >, - algorithms::Output< - edm4eic::ReconstructedParticleCollection - > - >; - - template - class JetReconstruction - : public JetReconstructionAlgorithm, - public WithPodConfig { - - public: - - JetReconstruction(std::string_view name) : - JetReconstructionAlgorithm { - name, - {"inputReconstructedParticles"}, - {"outputReconstructedParticles"}, - "Performs jet reconstruction using a FastJet algorithm." - } {} - - public: - - // algorithm initialization - void init(std::shared_ptr logger); - - // run algorithm - void process(const typename eicrecon::JetReconstructionAlgorithm::Input&, const typename eicrecon::JetReconstructionAlgorithm::Output&) const final; - - private: - - std::shared_ptr m_log; - - // fastjet components - std::unique_ptr m_jet_def; - std::unique_ptr m_area_def; - std::unique_ptr m_jet_plugin; - - // maps of user input onto fastjet options - std::map m_mapJetAlgo = { - {"kt_algorithm", fastjet::JetAlgorithm::kt_algorithm}, - {"cambridge_algorithm", fastjet::JetAlgorithm::cambridge_algorithm}, - {"antikt_algorithm", fastjet::JetAlgorithm::antikt_algorithm}, - {"genkt_algorithm", fastjet::JetAlgorithm::genkt_algorithm}, - {"cambridge_for_passive_algorithm", fastjet::JetAlgorithm::cambridge_for_passive_algorithm}, - {"genkt_for_passive_algorithm", fastjet::JetAlgorithm::genkt_for_passive_algorithm}, - {"ee_kt_algorithm", fastjet::JetAlgorithm::ee_kt_algorithm}, - {"ee_genkt_algorithm", fastjet::JetAlgorithm::ee_genkt_algorithm}, - {"plugin_algorithm", fastjet::JetAlgorithm::plugin_algorithm} - }; - std::map m_mapRecombScheme = { - {"E_scheme", fastjet::RecombinationScheme::E_scheme}, - {"pt_scheme", fastjet::RecombinationScheme::pt_scheme}, - {"pt2_scheme", fastjet::RecombinationScheme::pt2_scheme}, - {"Et_scheme", fastjet::RecombinationScheme::Et_scheme}, - {"Et2_scheme", fastjet::RecombinationScheme::Et2_scheme}, - {"BIpt_scheme", fastjet::RecombinationScheme::BIpt_scheme}, - {"BIpt2_scheme", fastjet::RecombinationScheme::BIpt2_scheme}, - {"WTA_pt_scheme", fastjet::RecombinationScheme::WTA_pt_scheme}, - {"WTA_modp_scheme", fastjet::RecombinationScheme::WTA_modp_scheme}, - {"external_scheme", fastjet::RecombinationScheme::external_scheme} - }; - std::map m_mapAreaType = { - {"active_area", fastjet::AreaType::active_area}, - {"active_area_explicit_ghosts", fastjet::AreaType::active_area_explicit_ghosts}, - {"one_ghost_passive_area", fastjet::AreaType::one_ghost_passive_area}, - {"passive_area", fastjet::AreaType::passive_area}, - {"voronoi_area", fastjet::AreaType::voronoi_area} - }; - - // default fastjet options - const struct defaults { - std::string jetAlgo; - std::string recombScheme; - std::string areaType; - } m_defaultFastjetOpts = {"antikt_algorithm", "E_scheme", "active_area"}; - - }; // end JetReconstruction definition - -} // end eicrecon namespace +template +using JetReconstructionAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +template +class JetReconstruction : public JetReconstructionAlgorithm, + public WithPodConfig { + +public: + JetReconstruction(std::string_view name) + : JetReconstructionAlgorithm{ + name, + {"inputReconstructedParticles"}, + {"outputReconstructedParticles"}, + "Performs jet reconstruction using a FastJet algorithm."} {} + +public: + // algorithm initialization + void init() final; + + // run algorithm + void process(const typename eicrecon::JetReconstructionAlgorithm::Input&, + const typename eicrecon::JetReconstructionAlgorithm::Output&) const final; + +private: + // fastjet components + std::unique_ptr m_jet_def; + std::unique_ptr m_area_def; + std::unique_ptr m_jet_plugin; + + // maps of user input onto fastjet options + std::map m_mapJetAlgo = { + {"kt_algorithm", fastjet::JetAlgorithm::kt_algorithm}, + {"cambridge_algorithm", fastjet::JetAlgorithm::cambridge_algorithm}, + {"antikt_algorithm", fastjet::JetAlgorithm::antikt_algorithm}, + {"genkt_algorithm", fastjet::JetAlgorithm::genkt_algorithm}, + {"cambridge_for_passive_algorithm", fastjet::JetAlgorithm::cambridge_for_passive_algorithm}, + {"genkt_for_passive_algorithm", fastjet::JetAlgorithm::genkt_for_passive_algorithm}, + {"ee_kt_algorithm", fastjet::JetAlgorithm::ee_kt_algorithm}, + {"ee_genkt_algorithm", fastjet::JetAlgorithm::ee_genkt_algorithm}, + {"plugin_algorithm", fastjet::JetAlgorithm::plugin_algorithm}}; + std::map m_mapRecombScheme = { + {"E_scheme", fastjet::RecombinationScheme::E_scheme}, + {"pt_scheme", fastjet::RecombinationScheme::pt_scheme}, + {"pt2_scheme", fastjet::RecombinationScheme::pt2_scheme}, + {"Et_scheme", fastjet::RecombinationScheme::Et_scheme}, + {"Et2_scheme", fastjet::RecombinationScheme::Et2_scheme}, + {"BIpt_scheme", fastjet::RecombinationScheme::BIpt_scheme}, + {"BIpt2_scheme", fastjet::RecombinationScheme::BIpt2_scheme}, + {"WTA_pt_scheme", fastjet::RecombinationScheme::WTA_pt_scheme}, + {"WTA_modp_scheme", fastjet::RecombinationScheme::WTA_modp_scheme}, + {"external_scheme", fastjet::RecombinationScheme::external_scheme}}; + std::map m_mapAreaType = { + {"active_area", fastjet::AreaType::active_area}, + {"active_area_explicit_ghosts", fastjet::AreaType::active_area_explicit_ghosts}, + {"one_ghost_passive_area", fastjet::AreaType::one_ghost_passive_area}, + {"passive_area", fastjet::AreaType::passive_area}, + {"voronoi_area", fastjet::AreaType::voronoi_area}}; + + // default fastjet options + const struct defaults { + std::string jetAlgo; + std::string recombScheme; + std::string areaType; + } m_defaultFastjetOpts = {"antikt_algorithm", "E_scheme", "active_area"}; + +}; // end JetReconstruction definition + +} // namespace eicrecon diff --git a/src/algorithms/reco/TransformBreitFrame.cc b/src/algorithms/reco/TransformBreitFrame.cc index c8eff8221d..5e168845d3 100644 --- a/src/algorithms/reco/TransformBreitFrame.cc +++ b/src/algorithms/reco/TransformBreitFrame.cc @@ -20,14 +20,6 @@ namespace eicrecon { - void TransformBreitFrame::init(std::shared_ptr logger) { - - m_log = logger; - m_log->trace("Initialized"); - - } // end 'init(std::shared_ptr)' - - void TransformBreitFrame::process( const TransformBreitFrame::Input& input, const TransformBreitFrame::Output& output @@ -42,7 +34,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcpart); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector e_initial( @@ -56,7 +48,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcpart); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector p_initial( @@ -67,11 +59,11 @@ namespace eicrecon { m_crossingAngle) ); - m_log->debug("electron energy, proton energy = {},{}",e_initial.E(),p_initial.E()); + debug("electron energy, proton energy = {},{}",e_initial.E(),p_initial.E()); // Get the event kinematics, set up transform if (kine->size() == 0) { - m_log->debug("No kinematics found"); + debug("No kinematics found"); return; } @@ -82,15 +74,15 @@ namespace eicrecon { // Use relation to get reconstructed scattered electron const PxPyPzEVector e_final = edm4hep::utils::detail::p4(evt_kin.getScat(),&edm4hep::utils::UseEnergy); - m_log->debug("scattered electron in lab frame px,py,pz,E = {},{},{},{}", + debug("scattered electron in lab frame px,py,pz,E = {},{},{},{}", e_final.Px(),e_final.Py(),e_final.Pz(),e_final.E()); // Set up the transformation const PxPyPzEVector virtual_photon = (e_initial - e_final); - m_log->debug("virtual photon in lab frame px,py,pz,E = {},{},{},{}", + debug("virtual photon in lab frame px,py,pz,E = {},{},{},{}", virtual_photon.Px(),virtual_photon.Py(),virtual_photon.Pz(),virtual_photon.E()); - m_log->debug("x, Q^2 = {},{}",meas_x,meas_Q2); + debug("x, Q^2 = {},{}",meas_x,meas_Q2); // Set up the transformation (boost) to the Breit frame const auto P3 = p_initial.Vect(); @@ -116,9 +108,9 @@ namespace eicrecon { e_final_breit = breitRot*e_final_breit; virtual_photon_breit = breitRot*virtual_photon_breit; - m_log->debug("incoming hadron in Breit frame px,py,pz,E = {},{},{},{}", + debug("incoming hadron in Breit frame px,py,pz,E = {},{},{},{}", p_initial_breit.Px(),p_initial_breit.Py(),p_initial_breit.Pz(),p_initial_breit.E()); - m_log->debug("virtual photon in Breit frame px,py,pz,E = {},{},{},{}", + debug("virtual photon in Breit frame px,py,pz,E = {},{},{},{}", virtual_photon_breit.Px(),virtual_photon_breit.Py(),virtual_photon_breit.Pz(),virtual_photon_breit.E()); // look over the input particles and transform diff --git a/src/algorithms/reco/TransformBreitFrame.h b/src/algorithms/reco/TransformBreitFrame.h index af654c2772..74a89dd69b 100644 --- a/src/algorithms/reco/TransformBreitFrame.h +++ b/src/algorithms/reco/TransformBreitFrame.h @@ -7,8 +7,6 @@ #include #include #include -#include -#include #include #include @@ -16,42 +14,30 @@ namespace eicrecon { - using TransformBreitFrameAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::InclusiveKinematicsCollection, - edm4eic::ReconstructedParticleCollection - >, - algorithms::Output< - edm4eic::ReconstructedParticleCollection - > - >; +using TransformBreitFrameAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; - class TransformBreitFrame - : public TransformBreitFrameAlgorithm, - public WithPodConfig { +class TransformBreitFrame : public TransformBreitFrameAlgorithm, public WithPodConfig { - public: +public: + TransformBreitFrame(std::string_view name) + : TransformBreitFrameAlgorithm{ + name, + {"inputMCParticles", "inputInclusiveKinematics", "inputReconstructedParticles"}, + {"outputReconstructedParticles"}, + "Transforms a set of particles from the lab frame to the Breit frame"} {} - TransformBreitFrame(std::string_view name) : - TransformBreitFrameAlgorithm { - name, - {"inputMCParticles", "inputInclusiveKinematics", "inputReconstructedParticles"}, - {"outputReconstructedParticles"}, - "Transforms a set of particles from the lab frame to the Breit frame" - } {} + // algorithm initialization + void init() final{}; - // algorithm initialization - void init(std::shared_ptr logger); + // run algorithm + void process(const Input&, const Output&) const final; - // run algorithm - void process(const Input&, const Output&) const final; +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - private: +}; // end TransformBreitFrame definition - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - - }; // end TransformBreitFrame definition - -} // end eicrecon namespace +} // namespace eicrecon diff --git a/src/factories/digi/SiliconTrackerDigi_factory.h b/src/factories/digi/SiliconTrackerDigi_factory.h index 3975dd34da..a708b3335c 100644 --- a/src/factories/digi/SiliconTrackerDigi_factory.h +++ b/src/factories/digi/SiliconTrackerDigi_factory.h @@ -30,7 +30,7 @@ class SiliconTrackerDigi_factory : public JOmniFactory(GetPrefix()); m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/fardetectors/FarDetectorLinearProjection_factory.h b/src/factories/fardetectors/FarDetectorLinearProjection_factory.h index 39510d37ca..ebaeab6539 100644 --- a/src/factories/fardetectors/FarDetectorLinearProjection_factory.h +++ b/src/factories/fardetectors/FarDetectorLinearProjection_factory.h @@ -29,7 +29,7 @@ public JOmniFactory(GetPrefix()); - m_algo->level((algorithms::LogLevel)logger()->level()); + m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); m_algo->init(); } diff --git a/src/factories/fardetectors/FarDetectorLinearTracking_factory.h b/src/factories/fardetectors/FarDetectorLinearTracking_factory.h index aa3dd32247..5611e1fc26 100644 --- a/src/factories/fardetectors/FarDetectorLinearTracking_factory.h +++ b/src/factories/fardetectors/FarDetectorLinearTracking_factory.h @@ -33,7 +33,7 @@ class FarDetectorLinearTracking_factory : m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h b/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h index 13ace79e81..2134c064b2 100644 --- a/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h +++ b/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h @@ -39,7 +39,7 @@ class FarDetectorMLReconstruction_factory : public: void Configure() { m_algo = std::make_unique(GetPrefix()); - m_algo->level((algorithms::LogLevel)logger()->level()); + m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); m_algo->init(); } diff --git a/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h b/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h index d45aa97e3e..fa4b400d63 100644 --- a/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h +++ b/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h @@ -36,7 +36,7 @@ public JOmniFactorylevel(static_cast(logger()->level())); // Setup algorithm m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } diff --git a/src/factories/meta/CollectionCollector_factory.h b/src/factories/meta/CollectionCollector_factory.h index 9edf056f58..df5f6da39c 100644 --- a/src/factories/meta/CollectionCollector_factory.h +++ b/src/factories/meta/CollectionCollector_factory.h @@ -22,7 +22,7 @@ namespace eicrecon { void Configure() { m_algo = std::make_unique(this->GetPrefix()); m_algo->level(static_cast(this->logger()->level())); - m_algo->init(this->logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/HadronicFinalState_factory.h b/src/factories/reco/HadronicFinalState_factory.h index 30b0dea850..00d03523e3 100644 --- a/src/factories/reco/HadronicFinalState_factory.h +++ b/src/factories/reco/HadronicFinalState_factory.h @@ -32,7 +32,8 @@ class HadronicFinalState_factory : public: void Configure() { m_algo = std::make_unique(this->GetPrefix()); - m_algo->init(this->logger()); + m_algo->level(static_cast(this->logger()->level())); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/InclusiveKinematicsML_factory.h b/src/factories/reco/InclusiveKinematicsML_factory.h index 9f4ccd10d2..78343cfbd9 100644 --- a/src/factories/reco/InclusiveKinematicsML_factory.h +++ b/src/factories/reco/InclusiveKinematicsML_factory.h @@ -34,7 +34,7 @@ class InclusiveKinematicsML_factory : m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h index 7e99ba0253..7c395f39dc 100644 --- a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h +++ b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h @@ -33,7 +33,7 @@ class InclusiveKinematicsReconstructed_factory : void Configure() { m_algo = std::make_unique(this->GetPrefix()); m_algo->level(static_cast(this->logger()->level())); - m_algo->init(this->logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/InclusiveKinematicsTruth_factory.h b/src/factories/reco/InclusiveKinematicsTruth_factory.h index 6c5f644fdc..b119c0541b 100644 --- a/src/factories/reco/InclusiveKinematicsTruth_factory.h +++ b/src/factories/reco/InclusiveKinematicsTruth_factory.h @@ -30,7 +30,7 @@ class InclusiveKinematicsTruth_factory : void Configure() { m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/JetReconstruction_factory.h b/src/factories/reco/JetReconstruction_factory.h index 98a665e826..6919ae5d08 100644 --- a/src/factories/reco/JetReconstruction_factory.h +++ b/src/factories/reco/JetReconstruction_factory.h @@ -46,7 +46,7 @@ namespace eicrecon { m_algo = std::make_unique(this->GetPrefix()); m_algo->level(static_cast(this->logger()->level())); m_algo->applyConfig(FactoryT::config()); - m_algo->init(FactoryT::logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/TransformBreitFrame_factory.h b/src/factories/reco/TransformBreitFrame_factory.h index 70836b075b..3a9852818b 100644 --- a/src/factories/reco/TransformBreitFrame_factory.h +++ b/src/factories/reco/TransformBreitFrame_factory.h @@ -38,7 +38,7 @@ namespace eicrecon { void Configure() { m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/global/pid/MergeTrack_factory.h b/src/global/pid/MergeTrack_factory.h index e3fa78b518..c59f4250ba 100644 --- a/src/global/pid/MergeTrack_factory.h +++ b/src/global/pid/MergeTrack_factory.h @@ -33,7 +33,7 @@ class MergeTrack_factory : public JOmniFactory { void Configure() { m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { } diff --git a/src/tests/algorithms_test/pid_MergeTracks.cc b/src/tests/algorithms_test/pid_MergeTracks.cc index 1bf10a8c09..ee2d3fd3e1 100644 --- a/src/tests/algorithms_test/pid_MergeTracks.cc +++ b/src/tests/algorithms_test/pid_MergeTracks.cc @@ -1,16 +1,15 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2023, Christopher Dilks +#include #include #include #include +#include #include #include #include -#include -#include -#include -#include +#include #include #include @@ -21,9 +20,8 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { // initialize algorithm //---------------------------------------------------------- - std::shared_ptr logger = spdlog::default_logger()->clone("MergeTracks"); - logger->set_level(spdlog::level::debug); - algo.init(logger); + algo.level(algorithms::LogLevel::kDebug); + algo.init(); // helper functions and objects //---------------------------------------------------------- @@ -32,12 +30,12 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { struct Point { decltype(edm4eic::TrackPoint::position) position; decltype(edm4eic::TrackPoint::momentum) momentum; - decltype(edm4eic::TrackPoint::time) time; + decltype(edm4eic::TrackPoint::time) time; }; - auto make_track = [] (auto& coll, std::vector points) { + auto make_track = [](auto& coll, std::vector points) { auto trk = coll->create(); - for(auto& point : points) { + for (auto& point : points) { edm4eic::TrackPoint trk_point; trk_point.position = point.position; trk_point.momentum = point.momentum; @@ -46,14 +44,11 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { } }; - auto track_length = [] (auto trk) { + auto track_length = [](auto trk) { auto a = trk.getPoints(0); - auto b = trk.getPoints(trk.points_size()-1); - return std::hypot( - a.position.x - b.position.x, - a.position.y - b.position.y, - a.position.z - b.position.z - ); + auto b = trk.getPoints(trk.points_size() - 1); + return std::hypot(a.position.x - b.position.x, a.position.y - b.position.y, + a.position.z - b.position.z); }; // inputs @@ -73,112 +68,101 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { auto collection2 = std::make_unique(); // track0 - make_track(collection0, { // horizontal - { {0, 0, 0}, {1, 0, 0}, 1 }, // { position, momentum, time } - { {1, 0, 0}, {1, 0, 0}, 2 } - }); - make_track(collection1, { // horizontal - { {2, 0, 0}, {1, 0, 0}, 3 }, - { {3, 0, 0}, {1, 0, 0}, 4 } - }); - make_track(collection2, { // horizontal - { {4, 0, 0}, {1, 0, 0}, 5 }, - { {5, 0, 0}, {1, 0, 0}, 6 } - }); + make_track(collection0, { // horizontal + {{0, 0, 0}, {1, 0, 0}, 1}, // { position, momentum, time } + {{1, 0, 0}, {1, 0, 0}, 2}}); + make_track(collection1, {// horizontal + {{2, 0, 0}, {1, 0, 0}, 3}, + {{3, 0, 0}, {1, 0, 0}, 4}}); + make_track(collection2, {// horizontal + {{4, 0, 0}, {1, 0, 0}, 5}, + {{5, 0, 0}, {1, 0, 0}, 6}}); // track1 - make_track(collection0, { // horizontal - { {0, 0, 0}, {1, 0, 0}, 1 }, - { {1, 0, 0}, {1, 0, 0}, 2 } - }); - make_track(collection1, { // empty - { } - }); - make_track(collection2, { // one point - { {2, 0, 0}, {1, 0, 0}, 3 } - }); + make_track(collection0, {// horizontal + {{0, 0, 0}, {1, 0, 0}, 1}, + {{1, 0, 0}, {1, 0, 0}, 2}}); + make_track(collection1, {// empty + {}}); + make_track(collection2, {// one point + {{2, 0, 0}, {1, 0, 0}, 3}}); // track2 - make_track(collection0, { // vertical - { {0, 0, 0}, {0, 1, 0}, 1 }, - { {0, 1, 0}, {0, 1, 0}, 2 } - }); - make_track(collection1, { // horizontal - { {0, 1, 0}, {1, 0, 0}, 3 }, - { {1, 1, 0}, {1, 0, 0}, 4 } - }); - make_track(collection2, { // vertical - { {1, 1, 0}, {0, 1, 0}, 5 }, - { {1, 2, 0}, {0, 1, 0}, 6 } - }); + make_track(collection0, {// vertical + {{0, 0, 0}, {0, 1, 0}, 1}, + {{0, 1, 0}, {0, 1, 0}, 2}}); + make_track(collection1, {// horizontal + {{0, 1, 0}, {1, 0, 0}, 3}, + {{1, 1, 0}, {1, 0, 0}, 4}}); + make_track(collection2, {// vertical + {{1, 1, 0}, {0, 1, 0}, 5}, + {{1, 2, 0}, {0, 1, 0}, 6}}); // track3 - make_track(collection0, { // horizontal - { {1, 0, 0}, {1, 0, 0}, 2 }, - { {0, 0, 0}, {1, 0, 0}, 1 } - }); - make_track(collection1, { // horizontal - { {3, 0, 0}, {1, 0, 0}, 4 }, - { {4, 0, 0}, {1, 0, 0}, 5 } - }); - make_track(collection2, { // horizontal - { {5, 0, 0}, {1, 0, 0}, 6 }, - { {2, 0, 0}, {1, 0, 0}, 3 } - }); + make_track(collection0, {// horizontal + {{1, 0, 0}, {1, 0, 0}, 2}, + {{0, 0, 0}, {1, 0, 0}, 1}}); + make_track(collection1, {// horizontal + {{3, 0, 0}, {1, 0, 0}, 4}, + {{4, 0, 0}, {1, 0, 0}, 5}}); + make_track(collection2, {// horizontal + {{5, 0, 0}, {1, 0, 0}, 6}, + {{2, 0, 0}, {1, 0, 0}, 3}}); // tests //---------------------------------------------------------- SECTION("merge tracks from 1 collection") { // run algorithm - std::vector> colls = { collection0.get() }; + std::vector> colls = {collection0.get()}; // create output collection auto trks = std::make_unique(); algo.process({colls}, {trks.get()}); // input collection(s) size == output collection size REQUIRE(trks->size() == colls.front()->size()); // track length: from endpoints - REQUIRE_THAT( track_length(trks->at(0)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(2)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(3)), Catch::Matchers::WithinAbs(1, EPSILON) ); + REQUIRE_THAT(track_length(trks->at(0)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(2)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(3)), Catch::Matchers::WithinAbs(1, EPSILON)); // track length: from algorithm // FIXME when implemented in `MergeTracks` - for(const auto& trk : *trks) - REQUIRE_THAT( trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON) ); + for (const auto& trk : *trks) + REQUIRE_THAT(trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON)); } SECTION("merge tracks from 2 collections") { // run algorithm - std::vector> colls = { collection0.get(), collection1.get() }; + std::vector> colls = {collection0.get(), + collection1.get()}; auto trks = std::make_unique(); algo.process({colls}, {trks.get()}); // input collection(s) size == output collection size REQUIRE(trks->size() == colls.front()->size()); // track length: from endpoints - REQUIRE_THAT( track_length(trks->at(0)), Catch::Matchers::WithinAbs(3, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1,1), EPSILON) ); - REQUIRE_THAT( track_length(trks->at(3)), Catch::Matchers::WithinAbs(4, EPSILON) ); + REQUIRE_THAT(track_length(trks->at(0)), Catch::Matchers::WithinAbs(3, EPSILON)); + REQUIRE_THAT(track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1, 1), EPSILON)); + REQUIRE_THAT(track_length(trks->at(3)), Catch::Matchers::WithinAbs(4, EPSILON)); // track length: from algorithm // FIXME when implemented in `MergeTracks` - for(const auto& trk : *trks) - REQUIRE_THAT( trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON) ); + for (const auto& trk : *trks) + REQUIRE_THAT(trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON)); } SECTION("merge tracks from 3 collections") { // run algorithm - std::vector> colls = { collection0.get(), collection1.get(), collection2.get() }; + std::vector> colls = { + collection0.get(), collection1.get(), collection2.get()}; auto trks = std::make_unique(); algo.process({colls}, {trks.get()}); // input collection(s) size == output collection size REQUIRE(trks->size() == colls.front()->size()); // track length: from endpoints - REQUIRE_THAT( track_length(trks->at(0)), Catch::Matchers::WithinAbs(5, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(1)), Catch::Matchers::WithinAbs(2, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1,2), EPSILON) ); - REQUIRE_THAT( track_length(trks->at(3)), Catch::Matchers::WithinAbs(5, EPSILON) ); + REQUIRE_THAT(track_length(trks->at(0)), Catch::Matchers::WithinAbs(5, EPSILON)); + REQUIRE_THAT(track_length(trks->at(1)), Catch::Matchers::WithinAbs(2, EPSILON)); + REQUIRE_THAT(track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1, 2), EPSILON)); + REQUIRE_THAT(track_length(trks->at(3)), Catch::Matchers::WithinAbs(5, EPSILON)); // track length: from algorithm // FIXME when implemented in `MergeTracks` - for(const auto& trk : *trks) - REQUIRE_THAT( trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON) ); + for (const auto& trk : *trks) + REQUIRE_THAT(trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON)); } - } diff --git a/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc b/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc index 64c612782c..fd66879028 100644 --- a/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc +++ b/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc @@ -5,160 +5,143 @@ #include #include #include +#include #include #include #include #include -#include -#include -#include -#include #include -#include +#include #include #include #include "algorithms/fardetectors/FarDetectorTrackerCluster.h" #include "algorithms/fardetectors/FarDetectorTrackerClusterConfig.h" -TEST_CASE( "the clustering algorithm runs", "[FarDetectorTrackerCluster]" ) { +TEST_CASE("the clustering algorithm runs", "[FarDetectorTrackerCluster]") { eicrecon::FarDetectorTrackerCluster algo("FarDetectorTrackerCluster"); - std::shared_ptr logger = spdlog::default_logger()->clone("FarDetectorTrackerCluster"); - logger->set_level(spdlog::level::trace); - eicrecon::FarDetectorTrackerClusterConfig cfg; cfg.hit_time_limit = 10.0 * edm4eic::unit::ns; - cfg.readout = "MockTrackerHits"; - cfg.x_field = "x"; - cfg.y_field = "y"; + cfg.readout = "MockTrackerHits"; + cfg.x_field = "x"; + cfg.y_field = "y"; auto detector = algorithms::GeoSvc::instance().detector(); - auto id_desc = detector->readout(cfg.readout).idSpec(); + auto id_desc = detector->readout(cfg.readout).idSpec(); algo.applyConfig(cfg); - algo.init(logger); + algo.level(algorithms::LogLevel::kTrace); + algo.init(); - SECTION( "on a single pixel" ) { + SECTION("on a single pixel") { edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 0.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 0.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 1 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[0].x == 0.0 ); - REQUIRE( clusterPositions[0].y == 0.0 ); - REQUIRE( clusterPositions[0].energy == 5.0 ); - REQUIRE( clusterPositions[0].time == 0.0 ); + REQUIRE(clusterPositions.size() == 1); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[0].x == 0.0); + REQUIRE(clusterPositions[0].y == 0.0); + REQUIRE(clusterPositions[0].energy == 5.0); + REQUIRE(clusterPositions[0].time == 0.0); } - SECTION( "on two separated pixels" ) { + SECTION("on two separated pixels") { edm4eic::RawTrackerHitCollection hits_coll; hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 10}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + id_desc.encode({{"system", 255}, {"x", 0}, {"y", 10}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 10}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + id_desc.encode({{"system", 255}, {"x", 10}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 2 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[1].rawHits.size() == 1 ); + REQUIRE(clusterPositions.size() == 2); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[1].rawHits.size() == 1); } - SECTION( "on two adjacent pixels" ) { + SECTION("on two adjacent pixels") { edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 1 ); - REQUIRE( clusterPositions[0].rawHits.size() == 2 ); - REQUIRE( clusterPositions[0].x == 0.5 ); - + REQUIRE(clusterPositions.size() == 1); + REQUIRE(clusterPositions[0].rawHits.size() == 2); + REQUIRE(clusterPositions[0].x == 0.5); } - SECTION( "on two adjacent pixels outwith the time separation" ) { + SECTION("on two adjacent pixels outwith the time separation") { edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 0.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 0.0 // int32 timeStamp ); - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 1.1 * cfg.hit_time_limit // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 1.1 * cfg.hit_time_limit // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 2 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[1].rawHits.size() == 1 ); + REQUIRE(clusterPositions.size() == 2); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[1].rawHits.size() == 1); } - - SECTION( "run on three adjacent pixels" ) { + SECTION("run on three adjacent pixels") { // Check I and L shape clusters - auto pixel3 = GENERATE(std::vector{2,0}, std::vector{1,1}); - auto pixelCharges = GENERATE(std::vector{5,10,5}, std::vector{10,5,5}); - float pixel2Time = GENERATE_COPY(0, 1.1 * cfg.hit_time_limit); + auto pixel3 = GENERATE(std::vector{2, 0}, std::vector{1, 1}); + auto pixelCharges = GENERATE(std::vector{5, 10, 5}, std::vector{10, 5, 5}); + float pixel2Time = GENERATE_COPY(0, 1.1 * cfg.hit_time_limit); edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - pixelCharges[0], // int32 charge, - 0.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + pixelCharges[0], // int32 charge, + 0.0 // int32 timeStamp ); - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, - pixelCharges[1], // int32 charge, - pixel2Time // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, + pixelCharges[1], // int32 charge, + pixel2Time // int32 timeStamp ); hits_coll.create( - id_desc.encode({{"system", 255}, {"x", pixel3[0]}, {"y", pixel3[1]}}), // std::uint64_t cellID, - pixelCharges[2], // int32 charge, - 0.0 // int32 timeStamp + id_desc.encode( + {{"system", 255}, {"x", pixel3[0]}, {"y", pixel3[1]}}), // std::uint64_t cellID, + pixelCharges[2], // int32 charge, + 0.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - if(pixel2Time < cfg.hit_time_limit) { - REQUIRE( clusterPositions.size() == 1 ); - REQUIRE( clusterPositions[0].rawHits.size() == 3 ); + if (pixel2Time < cfg.hit_time_limit) { + REQUIRE(clusterPositions.size() == 1); + REQUIRE(clusterPositions[0].rawHits.size() == 3); + } else if (pixel3[0] == 2) { + REQUIRE(clusterPositions.size() == 3); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[1].rawHits.size() == 1); + REQUIRE(clusterPositions[2].rawHits.size() == 1); + } else { + REQUIRE(clusterPositions.size() == 2); } - else if(pixel3[0] == 2){ - REQUIRE( clusterPositions.size() == 3 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[1].rawHits.size() == 1 ); - REQUIRE( clusterPositions[2].rawHits.size() == 1 ); - } - else { - REQUIRE( clusterPositions.size() == 2 ); - } - - } } From 1413111c28f5c4e95e3710f4bd7102197fc6e450 Mon Sep 17 00:00:00 2001 From: minjungkim12 Date: Tue, 4 Jun 2024 05:32:27 +0200 Subject: [PATCH 10/10] implementation of ambiguity resolution solver from ACTS (#1383) ### Briefly, what does this PR introduce? * implementation of ambiguity resolution solver from ACTS; removing duplicate tracks from realistic seeding * new factory called after CKFtracking; taking CKFtracking outputs as inputs; provides a set of output collections ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [x] New feature (issue #__) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [x] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? * No change required; ### Does this PR change default behavior? * adding 3 PodIO output collection --------- Co-authored-by: Minjung Kim Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Minjung Kim Co-authored-by: Wouter Deconinck Co-authored-by: Dmitry Kalinkin Co-authored-by: Minjung Kim Co-authored-by: Barak Schmookler --- src/algorithms/tracking/ActsToTracks.cc | 5 +- src/algorithms/tracking/AmbiguitySolver.cc | 121 ++++++++++++++++++ src/algorithms/tracking/AmbiguitySolver.h | 42 ++++++ .../tracking/AmbiguitySolverConfig.h | 17 +++ src/global/tracking/AmbiguitySolver_factory.h | 52 ++++++++ src/global/tracking/tracking.cc | 71 ++++++++-- src/services/io/podio/JEventProcessorPODIO.cc | 8 ++ 7 files changed, 304 insertions(+), 12 deletions(-) create mode 100644 src/algorithms/tracking/AmbiguitySolver.cc create mode 100644 src/algorithms/tracking/AmbiguitySolver.h create mode 100644 src/algorithms/tracking/AmbiguitySolverConfig.h create mode 100644 src/global/tracking/AmbiguitySolver_factory.h diff --git a/src/algorithms/tracking/ActsToTracks.cc b/src/algorithms/tracking/ActsToTracks.cc index 0bfeee2bd4..179d1541cf 100644 --- a/src/algorithms/tracking/ActsToTracks.cc +++ b/src/algorithms/tracking/ActsToTracks.cc @@ -3,7 +3,6 @@ #include #include -#include #include #include #include @@ -13,7 +12,6 @@ #include #include #include -#include #include #include #include @@ -58,8 +56,7 @@ void ActsToTracks::process(const Input& input, const Output& output) const { } // Loop over all trajectories in a multiTrajectory - // FIXME: we only retain the first trackTips entry - for (auto trackTip : decltype(trackTips){trackTips.front()}) { + for (auto trackTip : trackTips) { // Collect the trajectory summary info auto trajectoryState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip); diff --git a/src/algorithms/tracking/AmbiguitySolver.cc b/src/algorithms/tracking/AmbiguitySolver.cc new file mode 100644 index 0000000000..4be8faa411 --- /dev/null +++ b/src/algorithms/tracking/AmbiguitySolver.cc @@ -0,0 +1,121 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Minjung Kim, Barak Schmookler +#include "AmbiguitySolver.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Acts/Utilities/Logger.hpp" +#include "AmbiguitySolverConfig.h" +#include "extensions/spdlog/SpdlogFormatters.h" // IWYU pragma: keep +#include "extensions/spdlog/SpdlogToActs.h" + +namespace eicrecon { + +Acts::GreedyAmbiguityResolution::Config +transformConfig(const eicrecon::AmbiguitySolverConfig& cfg) { + Acts::GreedyAmbiguityResolution::Config result; + result.maximumSharedHits = cfg.maximum_shared_hits; + result.maximumIterations = cfg.maximum_iterations; + result.nMeasurementsMin = cfg.n_measurements_min; + return result; +} + +static std::size_t sourceLinkHash(const Acts::SourceLink& a) { + return static_cast(a.get().index()); +} + +static bool sourceLinkEquality(const Acts::SourceLink& a, const Acts::SourceLink& b) { + return a.get().index() == + b.get().index(); +} + + +AmbiguitySolver::AmbiguitySolver() {} + + +void AmbiguitySolver::init(std::shared_ptr log) { + + m_log = log; + m_acts_logger = eicrecon::getSpdlogLogger("AmbiguitySolver", m_log); + m_acts_cfg = transformConfig(m_cfg); + m_core = std::make_unique(m_acts_cfg, logger().clone()); +} + + +std::tuple, std::vector> +AmbiguitySolver::process(std::vector input_container, + const edm4eic::Measurement2DCollection& meas2Ds) { + + // Assuming ActsExamples::ConstTrackContainer is compatible with Acts::ConstVectorTrackContainer + // Create track container + std::vector output_trajectories; + std::vector output_tracks; + + auto& input_trks = input_container.front(); + Acts::GreedyAmbiguityResolution::State state; + m_core->computeInitialState(*input_trks, state, &sourceLinkHash, &sourceLinkEquality); + m_core->resolve(state); + + ActsExamples::TrackContainer solvedTracks{std::make_shared(), + std::make_shared()}; + solvedTracks.ensureDynamicColumns(*input_trks); + + for (auto iTrack : state.selectedTracks) { + + auto destProxy = solvedTracks.getTrack(solvedTracks.addTrack()); + auto srcProxy = input_trks->getTrack(state.trackTips.at(iTrack)); + destProxy.copyFrom(srcProxy, false); + destProxy.tipIndex() = srcProxy.tipIndex(); + + } + + output_tracks.push_back(new ActsExamples::ConstTrackContainer( + std::make_shared(std::move(solvedTracks.container())), + input_trks->trackStateContainerHolder())); + + //Make output trajectories + ActsExamples::Trajectories::IndexedParameters parameters; + std::vector tips; + + for (const auto& track : *(output_tracks.front())) { + + tips.clear(); + parameters.clear(); + + tips.push_back(track.tipIndex()); + parameters.emplace( + std::pair{track.tipIndex(), + ActsExamples::TrackParameters{track.referenceSurface().getSharedPtr(), + track.parameters(), track.covariance(), + track.particleHypothesis()}}); + + output_trajectories.push_back(new ActsExamples::Trajectories( + ((*output_tracks.front())).trackStateContainer(), + tips, parameters)); + + } + + return std::make_tuple(std::move(output_tracks), std::move(output_trajectories)); +} + +} // namespace eicrecon diff --git a/src/algorithms/tracking/AmbiguitySolver.h b/src/algorithms/tracking/AmbiguitySolver.h new file mode 100644 index 0000000000..6244626500 --- /dev/null +++ b/src/algorithms/tracking/AmbiguitySolver.h @@ -0,0 +1,42 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Minjung Kim, Barak Schmookler +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Acts/AmbiguityResolution/GreedyAmbiguityResolution.hpp" +#include "AmbiguitySolverConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + +/*Reco Track Filtering Based on Greedy ambiguity resolution solver adopted from ACTS*/ +class AmbiguitySolver : public WithPodConfig { +public: + AmbiguitySolver(); + + void init(std::shared_ptr log); + +std::tuple< + std::vector, + std::vector + > + process(std::vector input_container,const edm4eic::Measurement2DCollection& meas2Ds); + +private: + std::shared_ptr m_log; + Acts::GreedyAmbiguityResolution::Config m_acts_cfg; + std::unique_ptr m_core; + /// Private access to the logging instance + std::shared_ptr m_acts_logger{nullptr}; + const Acts::Logger& logger() const { return *m_acts_logger; } +}; + +} // namespace eicrecon diff --git a/src/algorithms/tracking/AmbiguitySolverConfig.h b/src/algorithms/tracking/AmbiguitySolverConfig.h new file mode 100644 index 0000000000..f16296f049 --- /dev/null +++ b/src/algorithms/tracking/AmbiguitySolverConfig.h @@ -0,0 +1,17 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Minjung Kim + +#pragma once + +#include + +namespace eicrecon { +struct AmbiguitySolverConfig { + /// Maximum amount of shared hits per track. + std::uint32_t maximum_shared_hits = 1; + /// Maximum number of iterations + std::uint32_t maximum_iterations = 100000; + /// Minimum number of measurement to form a track. + std::size_t n_measurements_min = 3; +}; +} // namespace eicrecon diff --git a/src/global/tracking/AmbiguitySolver_factory.h b/src/global/tracking/AmbiguitySolver_factory.h new file mode 100644 index 0000000000..a72dad6db0 --- /dev/null +++ b/src/global/tracking/AmbiguitySolver_factory.h @@ -0,0 +1,52 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Minjung Kim, Barak Schmookler +#pragma once + +#include "algorithms/tracking/AmbiguitySolver.h" +#include "algorithms/tracking/AmbiguitySolverConfig.h" +#include "extensions/jana/JOmniFactory.h" +#include "extensions/spdlog/SpdlogMixin.h" +#include +#include +#include +#include +#include +#include + +namespace eicrecon { + +class AmbiguitySolver_factory + : public JOmniFactory { + +private: + using AlgoT = eicrecon::AmbiguitySolver; + std::unique_ptr m_algo; + + Input m_acts_tracks_input {this}; + PodioInput m_measurements_input {this}; + Output m_acts_tracks_output {this}; + Output m_acts_trajectories_output {this}; + + ParameterRef m_maximumSharedHits{this, "maximumSharedHits", config().maximum_shared_hits, + "Maximum number of shared hits allowed"}; + ParameterRef m_maximumIterations{this, "maximumIterations", config().maximum_iterations, + "Maximum number of iterations"}; + ParameterRef m_nMeasurementsMin{ + this, "nMeasurementsMin", config().n_measurements_min, + "Number of measurements required for further reconstruction"}; + +public: + void Configure() { + m_algo = std::make_unique(); + m_algo->applyConfig(config()); + m_algo->init(logger()); + } + + void ChangeRun(int64_t run_number) {} + + void Process(int64_t run_number, uint64_t event_number) { + std::tie(m_acts_tracks_output(),m_acts_trajectories_output()) = m_algo->process(m_acts_tracks_input(),*m_measurements_input()); + } +} ; + +} // namespace eicrecon diff --git a/src/global/tracking/tracking.cc b/src/global/tracking/tracking.cc index fa37a317d8..06a259441c 100644 --- a/src/global/tracking/tracking.cc +++ b/src/global/tracking/tracking.cc @@ -14,6 +14,7 @@ #include "ActsToTracks.h" #include "ActsToTracks_factory.h" +#include "AmbiguitySolver_factory.h" #include "CKFTracking_factory.h" #include "IterativeVertexFinder_factory.h" #include "TrackParamTruthInit_factory.h" @@ -62,10 +63,10 @@ void InitPlugin(JApplication *app) { std::vector input_collections; auto readouts = app->GetService()->detector()->readouts(); for (const auto& [hit_collection, rec_collection] : possible_collections) { - if (readouts.find(hit_collection) != readouts.end()) { - // Add the collection to the list of input collections - input_collections.push_back(rec_collection); - } + if (readouts.find(hit_collection) != readouts.end()) { + // Add the collection to the list of input collections + input_collections.push_back(rec_collection); + } } // Tracker hits collector @@ -89,8 +90,35 @@ void InitPlugin(JApplication *app) { "CentralTrackerMeasurements" }, { - "CentralCKFActsTrajectories", - "CentralCKFActsTracks", + "CentralCKFActsTrajectoriesUnfiltered", + "CentralCKFActsTracksUnfiltered", + }, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "CentralCKFTracksUnfiltered", + { + "CentralTrackerMeasurements", + "CentralCKFActsTrajectoriesUnfiltered", + }, + { + "CentralCKFTrajectoriesUnfiltered", + "CentralCKFTrackParametersUnfiltered", + "CentralCKFTracksUnfiltered", + }, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "AmbiguityResolutionSolver", + { + "CentralCKFActsTracksUnfiltered", + "CentralTrackerMeasurements" + }, + { + "CentralCKFActsTracks", + "CentralCKFActsTrajectories", }, app )); @@ -124,8 +152,35 @@ void InitPlugin(JApplication *app) { "CentralTrackerMeasurements" }, { - "CentralCKFSeededActsTrajectories", - "CentralCKFSeededActsTracks", + "CentralCKFSeededActsTrajectoriesUnfiltered", + "CentralCKFSeededActsTracksUnfiltered", + }, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "CentralCKFSeededTracksUnfiltered", + { + "CentralTrackerMeasurements", + "CentralCKFSeededActsTrajectoriesUnfiltered", + }, + { + "CentralCKFSeededTrajectoriesUnfiltered", + "CentralCKFSeededTrackParametersUnfiltered", + "CentralCKFSeededTracksUnfiltered", + }, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "SeededAmbiguityResolutionSolver", + { + "CentralCKFSeededActsTracksUnfiltered", + "CentralTrackerMeasurements" + }, + { + "CentralCKFSeededActsTracks", + "CentralCKFSeededActsTrajectories", }, app )); diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 6ebcdd13b0..7b36d2f66f 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -186,6 +186,14 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "CentralCKFSeededTrajectories", "CentralCKFSeededTracks", "CentralCKFSeededTrackParameters", + //tracking properties - true seeding + "CentralCKFTrajectoriesUnfiltered", + "CentralCKFTracksUnfiltered", + "CentralCKFTrackParametersUnfiltered", + //tracking properties - realistic seeding + "CentralCKFSeededTrajectoriesUnfiltered", + "CentralCKFSeededTracksUnfiltered", + "CentralCKFSeededTrackParametersUnfiltered", "InclusiveKinematicsDA", "InclusiveKinematicsJB", "InclusiveKinematicsSigma",