diff --git a/CODEOWNERS b/CODEOWNERS index 565ee6e9..1c777213 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -9,10 +9,14 @@ ################################################################################# * @c-dilks +src/iguana/algorithms/clas12/EventBuilderFilter/* @c-dilks +src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar +src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3 src/iguana/algorithms/clas12/MomentumCorrection/* @RichCap @c-dilks -src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12 +src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3 src/iguana/algorithms/clas12/SectorFinder/* @rtysonCLAS12 -src/iguana/algorithms/clas12/FTEnergyCorrection/* @asligonulacar +src/iguana/algorithms/clas12/ZVertexFilter/* @rtysonCLAS12 +src/iguana/algorithms/physics/DihadronKinematics/* @c-dilks +src/iguana/algorithms/physics/InclusiveKinematics/* @c-dilks +src/iguana/algorithms/physics/SingleHadronKinematics/* @c-dilks src/iguana/services/YAMLReader.* @rtysonCLAS12 @c-dilks -src/iguana/algorithms/clas12/PhotonGBTFilter/* @Gregtom3 -src/iguana/algorithms/clas12/FiducialFilter/* @Gregtom3 \ No newline at end of file diff --git a/doc/gen/Doxyfile.in b/doc/gen/Doxyfile.in index fad3d8d9..c923bc87 100644 --- a/doc/gen/Doxyfile.in +++ b/doc/gen/Doxyfile.in @@ -302,6 +302,8 @@ ALIASES += end_doc_example="@}" # ignore things we don't want to document ALIASES += doxygen_off="@cond NO_DOC" ALIASES += doxygen_on="@endcond" +# latex +ALIASES += latex{1}="@f$\1@f$" # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources # only. Doxygen will then generate output that is more tailored for C. For diff --git a/examples/config/my_combined_config_file.yaml b/examples/config/my_combined_config_file.yaml index eb5f5252..8666c074 100644 --- a/examples/config/my_combined_config_file.yaml +++ b/examples/config/my_combined_config_file.yaml @@ -21,17 +21,15 @@ example::ExampleAlgorithm: clas12::ZVertexFilter: - # default cuts - - default: - electron_vz: [ -33.0, 11.0 ] + log: debug # NOTE: you may control algorithm log levels here, like so - #RG-A fall2018 inbending + electron: + - default: + vz: [ -33.0, 11.0 ] - runs: [ 4760, 5419 ] - electron_vz: [ -13.0, 12.0 ] - - #RG-A fall2018 outbending + vz: [ -13.0, 12.0 ] - runs: [ 5420, 5674 ] - electron_vz: [ -18.0, 10.0 ] + vz: [ -18.0, 10.0 ] another::Algorithm: #Cuts below are just examples: diff --git a/examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml b/examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml index 411b1103..ddb51358 100644 --- a/examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml +++ b/examples/config/my_config_directory/algorithms/clas12/ZVertexFilter/Config.yaml @@ -1,14 +1,17 @@ # Cut values for different run periods clas12::ZVertexFilter: - # default cuts - - default: - electron_vz: [ -15.0, 15.0 ] + # scattered electron cuts + electron: - # RG-A fall2018 inbending - - runs: [ 4760, 5419 ] - electron_vz: [ -5.0, 3.0 ] + # default cuts + - default: + vz: [ -15.0, 15.0 ] - # RG-A fall2018 outbending - - runs: [ 5420, 5674 ] - electron_vz: [ -8.0, 7.0 ] + # RG-A fall2018 inbending + - runs: [ 4760, 5419 ] + vz: [ -5.0, 3.0 ] + + # RG-A fall2018 outbending + - runs: [ 5420, 5674 ] + vz: [ -8.0, 7.0 ] diff --git a/examples/config/my_z_vertex_cuts.yaml b/examples/config/my_z_vertex_cuts.yaml index 8d5d8e37..8b959550 100644 --- a/examples/config/my_z_vertex_cuts.yaml +++ b/examples/config/my_z_vertex_cuts.yaml @@ -1,14 +1,22 @@ # Cut values for different run periods clas12::ZVertexFilter: - # default cuts - - default: - electron_vz: [ -1.5, 1.3 ] + ################################################################################### + # NOTE: for convenience, you can also set the log level from the configuration file + log: trace + ################################################################################### - # RG-A fall2018 inbending - - runs: [ 4760, 5419 ] - electron_vz: [ -0.5, 0.5 ] + # scattered electron cuts + electron: - # RG-A fall2018 outbending - - runs: [ 5420, 5674 ] - electron_vz: [ -0.8, 0.7 ] + # default cuts + - default: + vz: [ -1.5, 1.3 ] + + # RG-A fall2018 inbending + - runs: [ 4760, 5419 ] + vz: [ -0.5, 0.5 ] + + # RG-A fall2018 outbending + - runs: [ 5420, 5674 ] + vz: [ -0.8, 0.7 ] diff --git a/examples/iguana_ex_fortran_01_action_functions.f b/examples/iguana_ex_fortran_01_action_functions.f index 92a7766f..d3950c6b 100644 --- a/examples/iguana_ex_fortran_01_action_functions.f +++ b/examples/iguana_ex_fortran_01_action_functions.f @@ -51,6 +51,7 @@ program iguana_ex_fortran_01_action_functions logical(c_bool) accept(N_MAX) ! filter real(c_double) qx, qy, qz, qE ! q vector real(c_double) Q2, x, y, W, nu ! inclusive kinematics + real(c_double) beamPz, targetM ! beam and target integer(c_int) key_vz_filter ! key for Z-vertex filter integer(c_int) key_inc_kin ! key for inclusive kinematics @@ -273,7 +274,8 @@ program iguana_ex_fortran_01_action_functions & px(i_ele), py(i_ele), pz(i_ele), & key_inc_kin, & qx, qy, qz, qE, - & Q2, x, y, W, nu) + & Q2, x, y, W, nu, + & beamPz, targetM) print *, '===> inclusive kinematics:' print *, ' q = (', qx, qy, qz, qE, ')' print *, ' Q2 = ', Q2 diff --git a/meson/this_iguana.sh.in b/meson/this_iguana.sh.in index d9339652..35218115 100644 --- a/meson/this_iguana.sh.in +++ b/meson/this_iguana.sh.in @@ -22,6 +22,8 @@ for arg in "$@"; do unset arg_githubCI return 2 ;; + source*|.*) # handle callers which cause `$1` to be 'source this_iguana.sh' or '. this_iguana.sh' + ;; *) echo "ERROR: unknown option '$arg'" >&2 unset arg diff --git a/meson/ubsan.supp b/meson/ubsan.supp index fe8381cd..51719faa 100644 --- a/meson/ubsan.supp +++ b/meson/ubsan.supp @@ -3,3 +3,4 @@ alignment:hipo::structure::getFloatAt alignment:hipo::structure::getShortAt alignment:hipo::structure::getDoubleAt alignment:hipo::structure::putDoubleAt +alignment:hipo::structure::putIntAt diff --git a/src/iguana/algorithms/Algorithm.cc b/src/iguana/algorithms/Algorithm.cc index 0f18e2e9..343a40c9 100644 --- a/src/iguana/algorithms/Algorithm.cc +++ b/src/iguana/algorithms/Algorithm.cc @@ -16,17 +16,17 @@ namespace iguana { template OPTION_TYPE Algorithm::GetOptionScalar(std::string const& key, YAMLReader::node_path_t node_path) const { - try { - CompleteOptionNodePath(key, node_path); - auto opt = GetCachedOption(key); - auto val = opt ? opt.value() : m_yaml_config->GetScalar(node_path); - PrintOptionValue(key, val); - return val; + CompleteOptionNodePath(key, node_path); + auto opt = GetCachedOption(key); + if(!opt.has_value()) { + opt = m_yaml_config->GetScalar(node_path); } - catch(std::runtime_error const& ex) { - m_log->Error("Failed to `GetOptionScalar` for key '{}'", key); + if(!opt.has_value()) { + m_log->Error("Failed to `GetOptionScalar` for key {:?}", key); throw std::runtime_error("config file parsing issue"); } + PrintOptionValue(key, opt.value()); + return opt.value(); } template int Algorithm::GetOptionScalar(std::string const& key, YAMLReader::node_path_t node_path) const; template double Algorithm::GetOptionScalar(std::string const& key, YAMLReader::node_path_t node_path) const; @@ -37,17 +37,17 @@ namespace iguana { template std::vector Algorithm::GetOptionVector(std::string const& key, YAMLReader::node_path_t node_path) const { - try { - CompleteOptionNodePath(key, node_path); - auto opt = GetCachedOption>(key); - auto val = opt ? opt.value() : m_yaml_config->GetVector(node_path); - PrintOptionValue(key, val); - return val; + CompleteOptionNodePath(key, node_path); + auto opt = GetCachedOption>(key); + if(!opt.has_value()) { + opt = m_yaml_config->GetVector(node_path); } - catch(std::runtime_error const& ex) { - m_log->Error("Failed to `GetOptionVector` for key '{}'", key); + if(!opt.has_value()) { + m_log->Error("Failed to `GetOptionVector` for key {:?}", key); throw std::runtime_error("config file parsing issue"); } + PrintOptionValue(key, opt.value()); + return opt.value(); } template std::vector Algorithm::GetOptionVector(std::string const& key, YAMLReader::node_path_t node_path) const; template std::vector Algorithm::GetOptionVector(std::string const& key, YAMLReader::node_path_t node_path) const; @@ -108,19 +108,35 @@ namespace iguana { void Algorithm::ParseYAMLConfig() { + + // start YAMLReader instance, if not yet started if(!m_yaml_config) { + // set config files and directories specified by `::SetConfigFile`, `::SetConfigDirectory`, etc. o_user_config_file = GetCachedOption("config_file").value_or(""); o_user_config_dir = GetCachedOption("config_dir").value_or(""); m_log->Debug("Instantiating `YAMLReader`"); m_yaml_config = std::make_unique("config|" + m_name); - m_yaml_config->SetLogLevel(m_log->GetLevel()); + m_yaml_config->SetLogLevel(m_log->GetLevel()); // synchronize log levels m_yaml_config->AddDirectory(o_user_config_dir); m_yaml_config->AddFile(m_default_config_file); m_yaml_config->AddFile(o_user_config_file); } else m_log->Debug("`YAMLReader` already instantiated for this algorithm; using that"); + + // parse the files m_yaml_config->LoadFiles(); + + // if "log" was not set by `SetOption` (i.e., not in `m_option_cache`) + // - NB: not using `GetCachedOption` here, since `T` can be a few different types for key=='log' + if(m_option_cache.find("log") == m_option_cache.end()) { + // check if 'log' is set in the YAML node for this algorithm + auto log_level_from_yaml = m_yaml_config->GetScalar({m_class_name, "log"}); + if(log_level_from_yaml) { + m_log->SetLevel(log_level_from_yaml.value()); + m_yaml_config->SetLogLevel(log_level_from_yaml.value()); + } + } } /////////////////////////////////////////////////////////////////////////////// diff --git a/src/iguana/algorithms/Algorithm.h b/src/iguana/algorithms/Algorithm.h index 845d49d4..fdc21866 100644 --- a/src/iguana/algorithms/Algorithm.h +++ b/src/iguana/algorithms/Algorithm.h @@ -88,8 +88,7 @@ namespace iguana { else m_log->Error("Option '{}' must be a string or a Logger::Level", key); } - else - m_option_cache[key] = val; + m_option_cache[key] = val; return val; } @@ -136,17 +135,17 @@ namespace iguana { /// @param name the directory name void SetConfigDirectory(std::string const& name); - protected: // methods - - /// Parse YAML configuration files. Sets `m_yaml_config`. - void ParseYAMLConfig(); - /// Get the index of a bank in a `hipo::banklist`; throws an exception if the bank is not found /// @param banks the list of banks this algorithm will use /// @param bank_name the name of the bank /// returns the `hipo::banklist` index of the bank hipo::banklist::size_type GetBankIndex(hipo::banklist& banks, std::string const& bank_name) const noexcept(false); + protected: // methods + + /// Parse YAML configuration files. Sets `m_yaml_config`. + void ParseYAMLConfig(); + /// Get the reference to a bank from a `hipo::banklist`; optionally checks if the bank name matches the expectation /// @param banks the `hipo::banklist` from which to get the specified bank /// @param idx the index of `banks` of the specified bank diff --git a/src/iguana/algorithms/AlgorithmSequence.cc b/src/iguana/algorithms/AlgorithmSequence.cc index db75a9aa..a2da0395 100644 --- a/src/iguana/algorithms/AlgorithmSequence.cc +++ b/src/iguana/algorithms/AlgorithmSequence.cc @@ -67,4 +67,22 @@ namespace iguana { m_log->Print(level, " - {}", algo->GetName()); } + void AlgorithmSequence::SetConfigFileForEachAlgorithm(std::string const& name) + { + for(auto const& algo : m_sequence) + algo->SetConfigFile(name); + } + + void AlgorithmSequence::SetConfigDirectoryForEachAlgorithm(std::string const& name) + { + for(auto const& algo : m_sequence) + algo->SetConfigDirectory(name); + } + + void AlgorithmSequence::ForEachAlgorithm(std::function func) + { + for(auto& algo : m_sequence) + func(algo); + } + } diff --git a/src/iguana/algorithms/AlgorithmSequence.h b/src/iguana/algorithms/AlgorithmSequence.h index 5f620a75..04e9b61b 100644 --- a/src/iguana/algorithms/AlgorithmSequence.h +++ b/src/iguana/algorithms/AlgorithmSequence.h @@ -93,6 +93,29 @@ namespace iguana { /// @param level the log level of the printout void PrintSequence(Logger::Level level = Logger::info) const; + /// @brief Set a custom configuration file for each algorithm in the sequence + /// + /// Use this function if you have a single configuration file for all the + /// algorithms in your sequence + /// @param name the configuration file name + void SetConfigFileForEachAlgorithm(std::string const& name); + + /// @brief Set a custom configuration file directory for each algorithm in the sequence + /// + /// Use this function if you have a single configuration file directory for all the + /// algorithms in your sequence + /// @param name the directory name + void SetConfigDirectoryForEachAlgorithm(std::string const& name); + + /// @brief Call a function for each algorithm in the sequence + /// + /// Use as: + /// ```cpp + /// ForEachAlgorithm([](auto& algo){ algo->...; }); + /// ``` + /// @param func the function to call for each algorithm `algo` + void ForEachAlgorithm(std::function func); + private: /// The sequence of algorithms diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc b/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc index 16b872b3..f869c229 100644 --- a/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc +++ b/src/iguana/algorithms/clas12/ZVertexFilter/Algorithm.cc @@ -64,7 +64,7 @@ namespace iguana::clas12 { std::lock_guard const lock(m_mutex); // NOTE: be sure to lock successive `ConcurrentParam::Save` calls !!! m_log->Trace("-> calling Reload({}, {})", runnum, key); o_runnum->Save(runnum, key); - o_electron_vz_cuts->Save(GetOptionVector("electron_vz", {GetConfig()->InRange("runs", runnum), "electron_vz"}), key); + o_electron_vz_cuts->Save(GetOptionVector("electron_vz", {"electron", GetConfig()->InRange("runs", runnum), "vz"}), key); } bool ZVertexFilter::Filter(double const zvertex, int const pid, int const status, concurrent_key_t key) const diff --git a/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml b/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml index 557f9d33..8a7b1a2e 100644 --- a/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml +++ b/src/iguana/algorithms/clas12/ZVertexFilter/Config.yaml @@ -1,14 +1,16 @@ -# Cut values for different run periods clas12::ZVertexFilter: - # default cuts - - default: - electron_vz: [ -20.0, 20.0 ] + # scattered electron cuts + electron: - # RG-A fall2018 inbending - - runs: [ 4760, 5419 ] - electron_vz: [ -13.0, 12.0 ] + # default cuts + - default: + vz: [ -20.0, 20.0 ] - # RG-A fall2018 outbending - - runs: [ 5420, 5674 ] - electron_vz: [ -18.0, 10.0 ] + # RG-A fall2018 inbending + - runs: [ 4760, 5419 ] + vz: [ -13.0, 12.0 ] + + # RG-A fall2018 outbending + - runs: [ 5420, 5674 ] + vz: [ -18.0, 10.0 ] diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index a15f8838..7ac67bcb 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -74,6 +74,24 @@ algo_dict = [ 'algorithm_needs_ROOT': true, 'test_args': {'banks': [ 'REC::Particle', 'RUN::config' ]}, }, + { + 'name': 'physics::SingleHadronKinematics', + 'algorithm_needs_ROOT': true, + 'has_action_yaml': false, # FIXME: needs a vector action function + 'test_args': { + 'banks': [ 'REC::Particle', 'RUN::config' ], + 'prerequisites': [ 'physics::InclusiveKinematics' ], + }, + }, + { + 'name': 'physics::DihadronKinematics', + 'algorithm_needs_ROOT': true, + 'has_action_yaml': false, # FIXME: needs a vector action function + 'test_args': { + 'banks': [ 'REC::Particle', 'RUN::config' ], + 'prerequisites': [ 'physics::InclusiveKinematics' ], + }, + }, { 'name': 'clas12::PhotonGBTFilter', 'algorithm_needs_ROOT': true, @@ -83,8 +101,21 @@ algo_dict = [ ] # make lists of objects to build; inclusion depends on whether ROOT is needed or not, and if we have ROOT -algo_sources = [ 'Algorithm.cc', 'AlgorithmFactory.cc', 'AlgorithmSequence.cc' ] -algo_headers = [ 'Algorithm.h', 'AlgorithmBoilerplate.h', 'TypeDefs.h', 'AlgorithmSequence.h' ] +algo_sources = [ + 'Algorithm.cc', + 'AlgorithmFactory.cc', + 'AlgorithmSequence.cc', +] +algo_headers = [ + 'Algorithm.h', + 'AlgorithmBoilerplate.h', + 'TypeDefs.h', + 'AlgorithmSequence.h', +] +if ROOT_dep.found() + algo_sources += [ 'physics/Tools.cc' ] + algo_headers += [ 'physics/Tools.h' ] +endif vdor_sources = [ 'Validator.cc' ] vdor_headers = [ 'Validator.h' ] algo_configs = [] diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc new file mode 100644 index 00000000..0705d988 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc @@ -0,0 +1,245 @@ +#include "Algorithm.h" +#include "TypeDefs.h" + +#include + +namespace iguana::physics { + + REGISTER_IGUANA_ALGORITHM(DihadronKinematics, "physics::DihadronKinematics"); + + void DihadronKinematics::Start(hipo::banklist& banks) + { + b_particle = GetBankIndex(banks, "REC::Particle"); + b_inc_kin = GetBankIndex(banks, "physics::InclusiveKinematics"); + + // create the output bank + // FIXME: generalize the groupid and itemid + auto result_schema = CreateBank( + banks, + b_result, + GetClassName(), + { + "pindex_a/S", + "pindex_b/S", + "pdg_a/I", + "pdg_b/I", + "Mh/D", + "z/D", + "MX/D", + "xF/D", + "phiH/D", + "phiR/D", + "theta/D" + }, + 0xF000, + 5); + i_pindex_a = result_schema.getEntryOrder("pindex_a"); + i_pindex_b = result_schema.getEntryOrder("pindex_b"); + i_pdg_a = result_schema.getEntryOrder("pdg_a"); + i_pdg_b = result_schema.getEntryOrder("pdg_b"); + i_Mh = result_schema.getEntryOrder("Mh"); + i_z = result_schema.getEntryOrder("z"); + i_MX = result_schema.getEntryOrder("MX"); + i_xF = result_schema.getEntryOrder("xF"); + i_phiH = result_schema.getEntryOrder("phiH"); + i_phiR = result_schema.getEntryOrder("phiR"); + i_theta = result_schema.getEntryOrder("theta"); + + // parse config file + ParseYAMLConfig(); + o_hadron_a_pdgs = GetOptionSet("hadron_a_list"); + o_hadron_b_pdgs = GetOptionSet("hadron_b_list"); + o_phi_r_method = GetOptionScalar("phi_r_method"); + o_theta_method = GetOptionScalar("theta_method"); + + // check phiR method + if(o_phi_r_method == "RT_via_covariant_kT") + m_phi_r_method = e_RT_via_covariant_kT; + else + throw std::runtime_error(fmt::format("unknown phi_r_method: {:?}", o_phi_r_method)); + + // check theta method + if(o_theta_method == "hadron_a") + m_theta_method = e_hadron_a; + else + throw std::runtime_error(fmt::format("unknown theta_method: {:?}", o_theta_method)); + + } + + /////////////////////////////////////////////////////////////////////////////// + + void DihadronKinematics::Run(hipo::banklist& banks) const + { + auto& particle_bank = GetBank(banks, b_particle, "REC::Particle"); + auto& inc_kin_bank = GetBank(banks, b_inc_kin, "physics::InclusiveKinematics"); + auto& result_bank = GetBank(banks, b_result, GetClassName()); + ShowBank(particle_bank, Logger::Header("INPUT PARTICLES")); + + if(particle_bank.getRowList().empty() || inc_kin_bank.getRowList().empty()) { + m_log->Debug("skip this event, since not all required banks have entries"); + return; + } + + // get beam and target momenta + // FIXME: makes some assumptions about the beam; this should be generalized... + ROOT::Math::PxPyPzMVector p_beam( + 0.0, + 0.0, + inc_kin_bank.getDouble("beamPz", 0), + particle::mass.at(particle::electron)); + ROOT::Math::PxPyPzMVector p_target( + 0.0, + 0.0, + 0.0, + inc_kin_bank.getDouble("targetM", 0)); + + // get virtual photon momentum + ROOT::Math::PxPyPzEVector p_q( + inc_kin_bank.getDouble("qx", 0), + inc_kin_bank.getDouble("qy", 0), + inc_kin_bank.getDouble("qz", 0), + inc_kin_bank.getDouble("qE", 0)); + + // get additional inclusive variables + auto W = inc_kin_bank.getDouble("W", 0); + + // boosts + ROOT::Math::Boost boost__qp((p_q + p_target).BoostToCM()); // CoM frame of target and virtual photon + auto p_q__qp = boost__qp(p_q); + + // build list of dihadron rows (pindices) + auto dih_rows = PairHadrons(particle_bank); + + // loop over dihadrons + result_bank.setRows(dih_rows.size()); + int dih_row = 0; + for(const auto& [row_a, row_b] : dih_rows) { + + // get hadron momenta + Hadron had_a{.row = row_a}; + Hadron had_b{.row = row_b}; + for(auto& had : {&had_a, &had_b}) { + had->pdg = particle_bank.getInt("pid", had->row); + had->p = ROOT::Math::PxPyPzMVector( + particle_bank.getFloat("px", had->row), + particle_bank.getFloat("py", had->row), + particle_bank.getFloat("pz", had->row), + particle::mass.at(static_cast(had->pdg))); + } + + // calculate dihadron momenta and boosts + auto p_Ph = had_a.p + had_b.p; + auto p_Ph__qp = boost__qp(p_Ph); + ROOT::Math::Boost boost__dih(p_Ph.BoostToCM()); // CoM frame of dihadron + + // calculate z + double z = p_target.Dot(p_Ph) / p_target.Dot(p_q); + + // calculate Mh + double Mh = p_Ph.M(); + + // calculate MX + double MX = (p_target + p_q - p_Ph).M(); + + // calculate xF + double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R()); + + // calculate phiH + double phiH = tools::PlaneAngle( + p_q.Vect(), + p_beam.Vect(), + p_q.Vect(), + p_Ph.Vect()).value_or(tools::UNDEF); + + // calculate PhiR + double phiR = tools::UNDEF; + switch(m_phi_r_method) { + case e_RT_via_covariant_kT: + { + for(auto& had : {&had_a, &had_b}) { + had->z = p_target.Dot(had->p) / p_target.Dot(p_q); + had->p_perp = tools::RejectVector(had->p.Vect(), p_q.Vect()); + } + if(had_a.p_perp.has_value() && had_b.p_perp.has_value()) { + auto RT = (had_b.z * had_a.p_perp.value() - had_a.z * had_b.p_perp.value()) / (had_a.z + had_b.z); + phiR = tools::PlaneAngle( + p_q.Vect(), + p_beam.Vect(), + p_q.Vect(), + RT).value_or(tools::UNDEF); + } + break; + } + } + + // calculate theta + double theta = tools::UNDEF; + switch(m_theta_method) { + case e_hadron_a: + { + theta = tools::VectorAngle( + boost__dih(had_a.p).Vect(), + p_Ph.Vect()).value_or(tools::UNDEF); + break; + } + } + + result_bank.putShort(i_pindex_a, dih_row, had_a.row); + result_bank.putShort(i_pindex_b, dih_row, had_b.row); + result_bank.putInt(i_pdg_a, dih_row, had_a.pdg); + result_bank.putInt(i_pdg_b, dih_row, had_b.pdg); + result_bank.putDouble(i_Mh, dih_row, Mh); + result_bank.putDouble(i_z, dih_row, z); + result_bank.putDouble(i_MX, dih_row, MX); + result_bank.putDouble(i_xF, dih_row, xF); + result_bank.putDouble(i_phiH, dih_row, phiH); + result_bank.putDouble(i_phiR, dih_row, phiR); + result_bank.putDouble(i_theta, dih_row, theta); + + dih_row++; + } + + ShowBank(result_bank, Logger::Header("CREATED BANK")); + } + + /////////////////////////////////////////////////////////////////////////////// + + std::vector> DihadronKinematics::PairHadrons(hipo::bank const& particle_bank) const { + std::vector> result; + // loop over REC::Particle rows, for hadron A + for(auto const& row_a : particle_bank.getRowList()) { + // check PDG is in the hadron-A list + if(auto pdg_a{particle_bank.getInt("pid", row_a)}; o_hadron_a_pdgs.find(pdg_a) != o_hadron_a_pdgs.end()) { + // loop over REC::Particle rows, for hadron B + for(auto const& row_b : particle_bank.getRowList()) { + // don't pair a particle with itself + if(row_a == row_b) + continue; + // check PDG is in the hadron-B list + if(auto pdg_b{particle_bank.getInt("pid", row_b)}; o_hadron_b_pdgs.find(pdg_b) != o_hadron_b_pdgs.end()) { + // if the PDGs of hadrons A and B are the same, don't double count + if(pdg_a == pdg_b && row_b < row_a) + continue; + // we have a unique dihadron, add it to the list + result.push_back({row_a, row_b}); + } + } + } + } + // trace logging + if(m_log->GetLevel() <= Logger::Level::trace) { + if(result.empty()) + m_log->Trace("=> no dihadrons in this event"); + else + m_log->Trace("=> number of dihadrons found: {}", result.size()); + } + return result; + } + + /////////////////////////////////////////////////////////////////////////////// + + void DihadronKinematics::Stop() + { + } + +} diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h new file mode 100644 index 00000000..ca8dd510 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Algorithm.h @@ -0,0 +1,123 @@ +#pragma once + +#include "iguana/algorithms/Algorithm.h" +#include "iguana/algorithms/physics/Tools.h" +#include +#include + +namespace iguana::physics { + + /// Set of dihadron kinematics variables + struct DihadronKinematicsVars { + /// @brief `REC::Particle` row (`pindex`) of hadron A + int pindex_a; + /// @brief `REC::Particle` row (`pindex`) of hadron B + int pindex_b; + /// @brief PDG code of hadron A + int pdg_a; + /// @brief PDG code of hadron B + int pdg_b; + /// @brief @latex{M_h}: Invariant mass of the dihadron + double Mh; + /// @brief @latex{z}: Momentum fraction of the fragmenting parton carried by the dihadron + double z; + /// @brief @latex{M_X(ehhX)}: Missing mass of the dihadron + double MX; + /// @brief @latex{x_F}: Feynman-x of the dihadron + double xF; + /// @brief @latex{\phi_h}: @latex{q}-azimuthal angle between the lepton-scattering plane and the @latex{\vec{q}\times\vec{P}_h} plane; + /// if the value is `tools::UNDEF`, the calculation failed + double phiH; + /// @brief @latex{\phi_R}: @latex{q}-azimuthal angle between the lepton-scattering plane and dihadron plane; + /// if the value is `tools::UNDEF`, the calculation failed + double phiR; + /// @brief @latex{\theta}: The "decay" angle of hadron A in the dihadron rest frame, with respect; + /// to the dihadron momentum direction + double theta; + }; + + /// @brief_algo Calculate semi-inclusive dihadron kinematic quantities defined in `iguana::physics::DihadronKinematicsVars` + /// + /// @begin_doc_algo{physics::DihadronKinematics | Creator} + /// @input_banks{REC::Particle, %physics::InclusiveKinematics} + /// @output_banks{%physics::DihadronKinematics} + /// @end_doc + /// + /// @begin_doc_config + /// @config_param{hadron_a_list | list[int] | list of "hadron A" PDGs} + /// @config_param{hadron_b_list | list[int] | list of "hadron B" PDGs} + /// @config_param{phi_r_method | string | method used to calculate @latex{\phi_R} (see section "phiR calculation methods" below)} + /// @config_param{theta_method | string | method used to calculate @latex{\theta} (see section "theta calculation methods" below)} + /// @end_doc + /// + /// Dihadron PDGs will be formed from pairs from `hadron_a_list` and `hadron_b_list`. For example, + /// if you define: + /// ```yaml + /// hadron_a_list: [ 211 ] + /// hadron_b_list: [ -211, 2212 ] + /// ``` + /// then the algorithm will calculate kinematics for @latex{\pi^+\pi^-} and @latex{\pi^+p} dihadrons; hadron A + /// is the @latex{\pi^+} for both of these, whereas hadron B is the @latex{\pi^-} for the former and the proton + /// for the latter. + /// + /// @par phiR calculation methods + /// - `"RT_via_covariant_kT"`: use @latex{R_T} computed via covariant @latex{k_T} formula + /// + /// @par theta calculation methods + /// - `"hadron_a"`: use hadron A's "decay angle" in the dihadron rest frame + class DihadronKinematics : public Algorithm + { + + DEFINE_IGUANA_ALGORITHM(DihadronKinematics, physics::DihadronKinematics) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + /// @brief form dihadrons by pairing hadrons + /// @param particle_bank the particle bank (`REC::Particle`) + /// @returns a list of pairs of hadron rows + std::vector> PairHadrons(hipo::bank const& particle_bank) const; + + private: + + // banklist indices + hipo::banklist::size_type b_particle; + hipo::banklist::size_type b_inc_kin; + hipo::banklist::size_type b_result; + + // `b_result` bank item indices + int i_pindex_a; + int i_pindex_b; + int i_pdg_a; + int i_pdg_b; + int i_Mh; + int i_z; + int i_MX; + int i_xF; + int i_phiH; + int i_phiR; + int i_theta; + + // config options + std::set o_hadron_a_pdgs; + std::set o_hadron_b_pdgs; + std::string o_phi_r_method; + std::string o_theta_method; + enum {e_RT_via_covariant_kT} m_phi_r_method; + enum {e_hadron_a} m_theta_method; + + // storage for a single hadron + struct Hadron { + int row; + int pdg; + ROOT::Math::PxPyPzMVector p; + double z; + std::optional p_perp; + }; + + }; + +} diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml b/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml new file mode 100644 index 00000000..2469ee91 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Config.yaml @@ -0,0 +1,7 @@ +physics::DihadronKinematics: + + hadron_a_list: [ 211 ] + hadron_b_list: [ -211, 211, 2212 ] + + phi_r_method: RT_via_covariant_kT + theta_method: hadron_a diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc new file mode 100644 index 00000000..6e936c6d --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.cc @@ -0,0 +1,123 @@ +#include "Validator.h" +#include "TypeDefs.h" + +#include +#include + +namespace iguana::physics { + + REGISTER_IGUANA_VALIDATOR(DihadronKinematicsValidator); + + void DihadronKinematicsValidator::Start(hipo::banklist& banks) + { + // define the algorithm sequence + m_algo_seq = std::make_unique(); + m_algo_seq->Add("physics::InclusiveKinematics"); + m_algo_seq->Add("physics::DihadronKinematics"); + m_algo_seq->SetOption("physics::DihadronKinematics", "log", m_log->GetLevel()); + m_algo_seq->SetOption>("physics::DihadronKinematics", "hadron_a_list", {particle::pi_plus}); + m_algo_seq->SetOption>("physics::DihadronKinematics", "hadron_b_list", {particle::pi_minus}); + m_algo_seq->Start(banks); + + // get bank indices + b_result = GetBankIndex(banks, "physics::DihadronKinematics"); + + // set an output file + auto output_dir = GetOutputDirectory(); + if(output_dir) { + m_output_file_basename = output_dir.value() + "/dihadron_kinematics"; + m_output_file = new TFile(m_output_file_basename + ".root", "RECREATE"); + } + + // define plots + gStyle->SetOptStat(0); + const int n_bins = 100; + plot_list = { + { + new TH1D("Mh_dist", "invariant mass M_{h} [GeV]", n_bins, 0, 4), + [](auto const& b, auto const r) { return b.getDouble("Mh", r); } + }, + { + new TH1D("z_dist", "z", n_bins, 0, 1), + [](auto const& b, auto const r) { return b.getDouble("z", r); } + }, + { + new TH1D("MX_dist", "missing mass M_{X} [GeV];", n_bins, 0, 4), + [](auto const& b, auto const r) { return b.getDouble("MX", r); } + }, + { + new TH1D("xF_dist", "x_{F};", n_bins, -1, 1), + [](auto const& b, auto const r) { return b.getDouble("xF", r); } + }, + { + new TH1D("phiH_dist", "#phi_{h};", n_bins, -M_PI, M_PI), + [](auto const& b, auto const r) { return b.getDouble("phiH", r); } + }, + { + new TH1D("phiR_dist", "#phi_{R}", n_bins, -M_PI, M_PI), + [](auto const& b, auto const r) { return b.getDouble("phiR", r); } + }, + { + new TH1D("theta_dist", "#theta;", n_bins, 0, M_PI), + [](auto const& b, auto const r) { return b.getDouble("theta", r); } + } + }; + + // format plots + for(auto& plot : plot_list) { + plot.hist->SetLineColor(kRed); + plot.hist->SetFillColor(kRed); + plot.hist->SetTitle( + TString(particle::title.at(particle::pi_plus)) + + TString(particle::title.at(particle::pi_minus)) + + " " + plot.hist->GetTitle()); + } + } + + + void DihadronKinematicsValidator::Run(hipo::banklist& banks) const + { + // calculate kinematics + m_algo_seq->Run(banks); + auto& result_bank = GetBank(banks, b_result, "physics::DihadronKinematics"); + + // skip events with no dihadrons + if(result_bank.getRowList().size() == 0) { + m_log->Debug("skip this event, since it has no kinematics results"); + return; + } + + // lock mutex and fill the plots + std::scoped_lock lock(m_mutex); + for(auto const& row : result_bank.getRowList()) { + for(auto& plot : plot_list) + plot.hist->Fill(plot.get_val(result_bank, row)); + } + } + + + void DihadronKinematicsValidator::Stop() + { + if(GetOutputDirectory()) { + int const n_cols = 4; + int const n_rows = (plot_list.size() - 1) / n_cols + 1; + auto canv = new TCanvas("canv", "canv", n_cols * 800, n_rows * 600); + canv->Divide(n_cols, n_rows); + int pad_num = 0; + for(auto& plot : plot_list) { + auto pad = canv->GetPad(++pad_num); + pad->cd(); + pad->SetGrid(1, 1); + pad->SetLeftMargin(0.12); + pad->SetRightMargin(0.12); + pad->SetBottomMargin(0.12); + plot.hist->Draw(); + } + canv->SaveAs(m_output_file_basename + ".png"); + m_output_file->Write(); + m_log->Info("Wrote output file {}", m_output_file->GetName()); + m_output_file->Close(); + } + } + +} diff --git a/src/iguana/algorithms/physics/DihadronKinematics/Validator.h b/src/iguana/algorithms/physics/DihadronKinematics/Validator.h new file mode 100644 index 00000000..d7ca5ee6 --- /dev/null +++ b/src/iguana/algorithms/physics/DihadronKinematics/Validator.h @@ -0,0 +1,38 @@ +#pragma once + +#include "iguana/algorithms/Validator.h" + +#include +#include +#include +#include + +namespace iguana::physics { + + /// @brief `iguana::physics::DihadronKinematics` validator + class DihadronKinematicsValidator : public Validator + { + + DEFINE_IGUANA_VALIDATOR(DihadronKinematicsValidator, physics::DihadronKinematicsValidator) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + private: + + hipo::banklist::size_type b_result; + + struct Plot1D { + TH1D* hist; + std::function get_val; + }; + std::vector plot_list; + + TString m_output_file_basename; + TFile* m_output_file; + }; + +} diff --git a/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml b/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml index 81aeabb5..144d24d1 100644 --- a/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml +++ b/src/iguana/algorithms/physics/InclusiveKinematics/Action.yaml @@ -46,3 +46,7 @@ actions: type: double - name: nu type: double + - name: beamPz + type: double + - name: targetM + type: double diff --git a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc index ff5bde31..2dbc89f6 100644 --- a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc +++ b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc @@ -18,19 +18,21 @@ namespace iguana::physics { banks, b_result, GetClassName(), - {"pindex/S", "Q2/D", "x/D", "y/D", "W/D", "nu/D", "qx/D", "qy/D", "qz/D", "qE/D"}, + {"pindex/S", "Q2/D", "x/D", "y/D", "W/D", "nu/D", "qx/D", "qy/D", "qz/D", "qE/D", "beamPz/D", "targetM/D"}, 0xF000, 1); - i_pindex = result_schema.getEntryOrder("pindex"); - i_Q2 = result_schema.getEntryOrder("Q2"); - i_x = result_schema.getEntryOrder("x"); - i_y = result_schema.getEntryOrder("y"); - i_W = result_schema.getEntryOrder("W"); - i_nu = result_schema.getEntryOrder("nu"); - i_qx = result_schema.getEntryOrder("qx"); - i_qy = result_schema.getEntryOrder("qy"); - i_qz = result_schema.getEntryOrder("qz"); - i_qE = result_schema.getEntryOrder("qE"); + i_pindex = result_schema.getEntryOrder("pindex"); + i_Q2 = result_schema.getEntryOrder("Q2"); + i_x = result_schema.getEntryOrder("x"); + i_y = result_schema.getEntryOrder("y"); + i_W = result_schema.getEntryOrder("W"); + i_nu = result_schema.getEntryOrder("nu"); + i_qx = result_schema.getEntryOrder("qx"); + i_qy = result_schema.getEntryOrder("qy"); + i_qz = result_schema.getEntryOrder("qz"); + i_qE = result_schema.getEntryOrder("qE"); + i_beamPz = result_schema.getEntryOrder("beamPz"); + i_targetM = result_schema.getEntryOrder("targetM"); // instantiate RCDB reader m_rcdb = std::make_unique("RCDB|" + GetName(), m_log->GetLevel()); @@ -111,6 +113,8 @@ namespace iguana::physics { result_bank.putDouble(i_qy, 0, result_vars.qy); result_bank.putDouble(i_qz, 0, result_vars.qz); result_bank.putDouble(i_qE, 0, result_vars.qE); + result_bank.putDouble(i_beamPz, 0, result_vars.beamPz); + result_bank.putDouble(i_targetM, 0, result_vars.targetM); ShowBank(result_bank, Logger::Header("CREATED BANK")); } @@ -244,16 +248,18 @@ namespace iguana::physics { ROOT::Math::PxPyPzMVector vec_target(target[px], target[py], target[pz], target[m]); ROOT::Math::PxPyPzMVector vec_lepton(lepton_px, lepton_py, lepton_pz, beam[m]); - auto vec_q = vec_beam - vec_lepton; - result.qx = vec_q.Px(); - result.qy = vec_q.Py(); - result.qz = vec_q.Pz(); - result.qE = vec_q.E(); - result.Q2 = -1 * vec_q.M2(); - result.x = result.Q2 / (2 * vec_q.Dot(vec_target)); - result.y = vec_target.Dot(vec_q) / vec_target.Dot(vec_beam); - result.W = (vec_beam + vec_target - vec_lepton).M(); - result.nu = vec_target.Dot(vec_q) / target[m]; + auto vec_q = vec_beam - vec_lepton; + result.qx = vec_q.Px(); + result.qy = vec_q.Py(); + result.qz = vec_q.Pz(); + result.qE = vec_q.E(); + result.Q2 = -1 * vec_q.M2(); + result.x = result.Q2 / (2 * vec_q.Dot(vec_target)); + result.y = vec_target.Dot(vec_q) / vec_target.Dot(vec_beam); + result.W = (vec_beam + vec_target - vec_lepton).M(); + result.nu = vec_target.Dot(vec_q) / target[m]; + result.beamPz = beam[pz]; + result.targetM = target[m]; return result; } diff --git a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h index dda8da34..d91a52ac 100644 --- a/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h +++ b/src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.h @@ -9,24 +9,28 @@ namespace iguana::physics { /// Set of inclusive kinematics variables struct InclusiveKinematicsVars { - /// @f$x@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{x}-component of virtual photon momentum @latex{q} vector_element_t qx; - /// @f$y@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{y}-component of virtual photon momentum @latex{q} vector_element_t qy; - /// @f$z@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{z}-component of virtual photon momentum @latex{q} vector_element_t qz; - /// @f$E@f$-component of virtual photon momentum @f$q@f$ + /// @brief @latex{E}-component of virtual photon momentum @latex{q} vector_element_t qE; - /// @f$Q2@f$ (GeV@f$^2@f$) + /// @brief @latex{Q^2} (GeV@latex{^2}) double Q2; - /// @f$x_B@f$ + /// @brief @latex{x_B} double x; - /// @f$y@f$ + /// @brief @latex{y} double y; - /// @f$W@f$ (GeV) + /// @brief @latex{W} (GeV) double W; - /// @f$\nu@f$ + /// @brief @latex{\nu} double nu; + /// @brief beam momentum @latex{z}-component (GeV) + double beamPz; + /// @brief target mass (GeV) + double targetM; }; /// @brief_algo Calculate inclusive kinematics quantities defined in `iguana::physics::InclusiveKinematicsVars` @@ -63,9 +67,9 @@ namespace iguana::physics { concurrent_key_t PrepareEvent(int const runnum, double const beam_energy = -1) const; /// @action_function{scalar creator} compute kinematics from the scattered lepton. - /// @param lepton_px scattered lepton momentum component @f$p_x@f$ (GeV) - /// @param lepton_py scattered lepton momentum component @f$p_y@f$ (GeV) - /// @param lepton_pz scattered lepton momentum component @f$p_z@f$ (GeV) + /// @param lepton_px scattered lepton momentum component @latex{p_x} (GeV) + /// @param lepton_py scattered lepton momentum component @latex{p_y} (GeV) + /// @param lepton_pz scattered lepton momentum component @latex{p_z} (GeV) /// @param key the return value of `::PrepareEvent` /// @returns the reconstructed inclusive kinematics in a `iguana::physics::InclusiveKinematicsVars` instance InclusiveKinematicsVars ComputeFromLepton( @@ -104,6 +108,8 @@ namespace iguana::physics { int i_qy; int i_qz; int i_qE; + int i_beamPz; + int i_targetM; // config options mutable std::unique_ptr> o_runnum; diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc new file mode 100644 index 00000000..d0fb6f4c --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc @@ -0,0 +1,167 @@ +#include "Algorithm.h" +#include "TypeDefs.h" +#include "iguana/algorithms/physics/Tools.h" + +#include + +namespace iguana::physics { + + REGISTER_IGUANA_ALGORITHM(SingleHadronKinematics, "physics::SingleHadronKinematics"); + + void SingleHadronKinematics::Start(hipo::banklist& banks) + { + b_particle = GetBankIndex(banks, "REC::Particle"); + b_inc_kin = GetBankIndex(banks, "physics::InclusiveKinematics"); + + // create the output bank + // FIXME: generalize the groupid and itemid + auto result_schema = CreateBank( + banks, + b_result, + GetClassName(), + { + "pindex/S", + "pdg/I", + "z/D", + "MX/D", + "xF/D", + "phiH/D", + "xi/D" + }, + 0xF000, + 7); + i_pindex = result_schema.getEntryOrder("pindex"); + i_pdg = result_schema.getEntryOrder("pdg"); + i_z = result_schema.getEntryOrder("z"); + i_MX = result_schema.getEntryOrder("MX"); + i_xF = result_schema.getEntryOrder("xF"); + i_phiH = result_schema.getEntryOrder("phiH"); + i_xi = result_schema.getEntryOrder("xi"); + + // parse config file + ParseYAMLConfig(); + o_hadron_pdgs = GetOptionSet("hadron_list"); + + } + + /////////////////////////////////////////////////////////////////////////////// + + void SingleHadronKinematics::Run(hipo::banklist& banks) const + { + auto& particle_bank = GetBank(banks, b_particle, "REC::Particle"); + auto& inc_kin_bank = GetBank(banks, b_inc_kin, "physics::InclusiveKinematics"); + auto& result_bank = GetBank(banks, b_result, GetClassName()); + ShowBank(particle_bank, Logger::Header("INPUT PARTICLES")); + + if(particle_bank.getRowList().empty() || inc_kin_bank.getRowList().empty()) { + m_log->Debug("skip this event, since not all required banks have entries"); + return; + } + + // get beam and target momenta + // FIXME: makes some assumptions about the beam; this should be generalized... + ROOT::Math::PxPyPzMVector p_beam( + 0.0, + 0.0, + inc_kin_bank.getDouble("beamPz", 0), + particle::mass.at(particle::electron)); + ROOT::Math::PxPyPzMVector p_target( + 0.0, + 0.0, + 0.0, + inc_kin_bank.getDouble("targetM", 0)); + + // get virtual photon momentum + ROOT::Math::PxPyPzEVector p_q( + inc_kin_bank.getDouble("qx", 0), + inc_kin_bank.getDouble("qy", 0), + inc_kin_bank.getDouble("qz", 0), + inc_kin_bank.getDouble("qE", 0)); + + // get additional inclusive variables + auto W = inc_kin_bank.getDouble("W", 0); + + // boosts + ROOT::Math::Boost boost__qp((p_q + p_target).BoostToCM()); // CoM frame of target and virtual photon + auto p_q__qp = boost__qp(p_q); + + // banks' row lists + auto const& particle_bank_rowlist = particle_bank.getRowList(); + hipo::bank::rowlist::list_t result_bank_rowlist{}; + result_bank.setRows(particle_bank.getRows()); + + // loop over ALL rows of `particle_bank` + // - we will calculate kinematics for rows in `particle_bank_rowlist`, and zero out all the other rows + // - we want the `result_bank` to have the same number of rows as `particle_bank` and the same ordering, + // so that banks which reference `particle_bank` rows can be used to reference `result_bank` rows too + for(int row = 0; row < particle_bank.getRows(); row++) { + + // if the particle is in `o_hadron_pdgs` AND the row is in `particle_bank`'s filtered row list + if(auto pdg{particle_bank.getInt("pid", row)}; + o_hadron_pdgs.find(pdg) != o_hadron_pdgs.end() && + std::find(particle_bank_rowlist.begin(), particle_bank_rowlist.end(), row) != particle_bank_rowlist.end()) { + + // hadron momentum + auto p_Ph = ROOT::Math::PxPyPzMVector( + particle_bank.getFloat("px", row), + particle_bank.getFloat("py", row), + particle_bank.getFloat("pz", row), + particle::mass.at(static_cast(pdg))); + auto p_Ph__qp = boost__qp(p_Ph); + + // calculate z + double z = p_target.Dot(p_Ph) / p_target.Dot(p_q); + + // calculate xi + double xi = p_q.Dot(p_Ph) / p_target.Dot(p_q); + + // calculate MX + double MX = (p_target + p_q - p_Ph).M(); + + // calculate xF + double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R()); + + // calculate phiH + double phiH = tools::PlaneAngle( + p_q.Vect(), + p_beam.Vect(), + p_q.Vect(), + p_Ph.Vect()).value_or(tools::UNDEF); + + // put this particle in `result_bank`'s row list + result_bank_rowlist.push_back(row); + + // fill the bank + result_bank.putShort(i_pindex, row, row); + result_bank.putInt(i_pdg, row, pdg); + result_bank.putDouble(i_z, row, z); + result_bank.putDouble(i_MX, row, MX); + result_bank.putDouble(i_xF, row, xF); + result_bank.putDouble(i_phiH, row, phiH); + result_bank.putDouble(i_xi, row, xi); + } + else { + // zero the row + result_bank.putShort(i_pindex, row, row); + result_bank.putInt(i_pdg, row, pdg); + result_bank.putDouble(i_z, row, 0); + result_bank.putDouble(i_MX, row, 0); + result_bank.putDouble(i_xF, row, 0); + result_bank.putDouble(i_phiH, row, 0); + result_bank.putDouble(i_xi, row, 0); + } + } + + // apply the filtered rowlist to `result_bank` + result_bank.getMutableRowList().setList(result_bank_rowlist); + + ShowBank(result_bank, Logger::Header("CREATED BANK")); + } + + /////////////////////////////////////////////////////////////////////////////// + + void SingleHadronKinematics::Stop() + { + } + +} diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h new file mode 100644 index 00000000..f98d9adb --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.h @@ -0,0 +1,79 @@ +#pragma once + +#include "iguana/algorithms/Algorithm.h" +#include +#include + +namespace iguana::physics { + + /// Set of hadron kinematics variables + struct SingleHadronKinematicsVars { + /// @brief `REC::Particle` row (`pindex`) of the hadron + int pindex; + /// @brief PDG code of the hadron + int pdg; + /// @brief @latex{z}: Momentum fraction of the fragmenting parton carried by the hadron + double z; + /// @brief @latex{M_X(ehX)}: Missing mass of the hadron + double MX; + /// @brief @latex{x_F}: Feynman-x of the hadron + double xF; + /// @brief @latex{\phi_h}: @latex{q}-azimuthal angle between the lepton-scattering plane and the @latex{\vec{q}\times\vec{P}_h} plane; + /// if the value is `tools::UNDEF`, the calculation failed + double phiH; + /// @brief @latex{\xi_h}: Longitudinal momentum fraction of the nucleon carried by the hadron + double xi; + }; + + /// @brief_algo Calculate semi-inclusive hadron kinematic quantities defined in `iguana::physics::SingleHadronKinematicsVars` + /// + /// @begin_doc_algo{physics::SingleHadronKinematics | Creator} + /// @input_banks{REC::Particle, %physics::InclusiveKinematics} + /// @output_banks{%physics::SingleHadronKinematics} + /// @end_doc + /// + /// @begin_doc_config + /// @config_param{hadron_list | list[int] | list of hadron PDGs} + /// @end_doc + /// + /// The output bank `%physics::SingleHadronKinematics` will have the same number of rows as the input particle bank `REC::Particle` + /// - we want the output bank to have the same number of rows and ordering as the input + /// particle bank, so that banks which reference the input particle bank's rows (usually via `pindex`) can be used to + /// reference the output bank's rows too + /// - rows of the input particle bank which were filtered upstream will also be filtered out here, and all the values of the + /// corresponding row in the output bank will be zeroed, since no calculations are performed for + /// those particles + /// - particles which are not listed in the configuration parameter `hadron_list` will also be filtered out and zeroed + class SingleHadronKinematics : public Algorithm + { + + DEFINE_IGUANA_ALGORITHM(SingleHadronKinematics, physics::SingleHadronKinematics) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + private: + + // banklist indices + hipo::banklist::size_type b_particle; + hipo::banklist::size_type b_inc_kin; + hipo::banklist::size_type b_result; + + // `b_result` bank item indices + int i_pindex; + int i_pdg; + int i_z; + int i_MX; + int i_xF; + int i_phiH; + int i_xi; + + // config options + std::set o_hadron_pdgs; + + }; + +} diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Config.yaml b/src/iguana/algorithms/physics/SingleHadronKinematics/Config.yaml new file mode 100644 index 00000000..b2277af3 --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Config.yaml @@ -0,0 +1,3 @@ +physics::SingleHadronKinematics: + + hadron_list: [ 211, -211, 2212 ] diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc new file mode 100644 index 00000000..d84edeb9 --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.cc @@ -0,0 +1,113 @@ +#include "Validator.h" +#include "TypeDefs.h" + +#include +#include + +namespace iguana::physics { + + REGISTER_IGUANA_VALIDATOR(SingleHadronKinematicsValidator); + + void SingleHadronKinematicsValidator::Start(hipo::banklist& banks) + { + // define the algorithm sequence + m_algo_seq = std::make_unique(); + m_algo_seq->Add("physics::InclusiveKinematics"); + m_algo_seq->Add("physics::SingleHadronKinematics"); + m_algo_seq->SetOption("physics::SingleHadronKinematics", "log", m_log->GetLevel()); + m_algo_seq->SetOption>("physics::SingleHadronKinematics", "hadron_list", {particle::pi_plus}); + m_algo_seq->Start(banks); + + // get bank indices + b_result = GetBankIndex(banks, "physics::SingleHadronKinematics"); + + // set an output file + auto output_dir = GetOutputDirectory(); + if(output_dir) { + m_output_file_basename = output_dir.value() + "/single_hadron_kinematics"; + m_output_file = new TFile(m_output_file_basename + ".root", "RECREATE"); + } + + // define plots + gStyle->SetOptStat(0); + const int n_bins = 100; + plot_list = { + { + new TH1D("z_dist", "z", n_bins, 0, 1), + [](auto const& b, auto const r) { return b.getDouble("z", r); } + }, + { + new TH1D("MX_dist", "missing mass M_{X} [GeV];", n_bins, 0, 4), + [](auto const& b, auto const r) { return b.getDouble("MX", r); } + }, + { + new TH1D("xF_dist", "x_{F};", n_bins, -1, 1), + [](auto const& b, auto const r) { return b.getDouble("xF", r); } + }, + { + new TH1D("phiH_dist", "#phi_{h};", n_bins, -M_PI, M_PI), + [](auto const& b, auto const r) { return b.getDouble("phiH", r); } + }, + { + new TH1D("xi_dist", "#xi", n_bins, -1, 1), + [](auto const& b, auto const r) { return b.getDouble("xi", r); } + }, + }; + + // format plots + for(auto& plot : plot_list) { + plot.hist->SetLineColor(kGreen+1); + plot.hist->SetFillColor(kGreen+1); + plot.hist->SetTitle( + TString(particle::title.at(particle::pi_plus)) + + " " + plot.hist->GetTitle()); + } + } + + + void SingleHadronKinematicsValidator::Run(hipo::banklist& banks) const + { + // calculate kinematics + m_algo_seq->Run(banks); + auto& result_bank = GetBank(banks, b_result, "physics::SingleHadronKinematics"); + + // skip events with no hadrons + if(result_bank.getRowList().size() == 0) { + m_log->Debug("skip this event, since it has no kinematics results"); + return; + } + + // lock mutex and fill the plots + std::scoped_lock lock(m_mutex); + for(auto const& row : result_bank.getRowList()) { + for(auto& plot : plot_list) + plot.hist->Fill(plot.get_val(result_bank, row)); + } + } + + + void SingleHadronKinematicsValidator::Stop() + { + if(GetOutputDirectory()) { + int const n_cols = 4; + int const n_rows = (plot_list.size() - 1) / n_cols + 1; + auto canv = new TCanvas("canv", "canv", n_cols * 800, n_rows * 600); + canv->Divide(n_cols, n_rows); + int pad_num = 0; + for(auto& plot : plot_list) { + auto pad = canv->GetPad(++pad_num); + pad->cd(); + pad->SetGrid(1, 1); + pad->SetLeftMargin(0.12); + pad->SetRightMargin(0.12); + pad->SetBottomMargin(0.12); + plot.hist->Draw(); + } + canv->SaveAs(m_output_file_basename + ".png"); + m_output_file->Write(); + m_log->Info("Wrote output file {}", m_output_file->GetName()); + m_output_file->Close(); + } + } + +} diff --git a/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.h b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.h new file mode 100644 index 00000000..3017da57 --- /dev/null +++ b/src/iguana/algorithms/physics/SingleHadronKinematics/Validator.h @@ -0,0 +1,38 @@ +#pragma once + +#include "iguana/algorithms/Validator.h" + +#include +#include +#include +#include + +namespace iguana::physics { + + /// @brief `iguana::physics::SingleHadronKinematics` validator + class SingleHadronKinematicsValidator : public Validator + { + + DEFINE_IGUANA_VALIDATOR(SingleHadronKinematicsValidator, physics::SingleHadronKinematicsValidator) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + private: + + hipo::banklist::size_type b_result; + + struct Plot1D { + TH1D* hist; + std::function get_val; + }; + std::vector plot_list; + + TString m_output_file_basename; + TFile* m_output_file; + }; + +} diff --git a/src/iguana/algorithms/physics/Tools.cc b/src/iguana/algorithms/physics/Tools.cc new file mode 100644 index 00000000..e06c9748 --- /dev/null +++ b/src/iguana/algorithms/physics/Tools.cc @@ -0,0 +1,57 @@ +#include "Tools.h" + +namespace iguana::physics::tools { + + std::optional PlaneAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b, + ROOT::Math::XYZVector const v_c, + ROOT::Math::XYZVector const v_d) + { + auto cross_ab = v_a.Cross(v_b); // A x B + auto cross_cd = v_c.Cross(v_d); // C x D + + auto sgn = cross_ab.Dot(v_d); // (A x B) . D + if(!(std::abs(sgn) > 0)) + return std::nullopt; + sgn /= std::abs(sgn); // sign of (A x B) . D + + auto numer = cross_ab.Dot(cross_cd); // (A x B) . (C x D) + auto denom = cross_ab.R() * cross_cd.R(); // |A x B| * |C x D| + if(!(std::abs(denom) > 0)) + return std::nullopt; + return sgn * std::acos(numer / denom); + } + + std::optional ProjectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b) + { + auto denom = v_b.Dot(v_b); + if(!(std::abs(denom) > 0)) + return std::nullopt; + return v_b * ( v_a.Dot(v_b) / denom ); + } + + std::optional RejectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b) + { + auto v_c = ProjectVector(v_a, v_b); + if(v_c.has_value()) + return v_a - v_c.value(); + return std::nullopt; + } + + std::optional VectorAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b) + { + double m = v_a.R() * v_b.R(); + if(m > 0) + return std::acos(v_a.Dot(v_b) / m); + return std::nullopt; + } + +} + diff --git a/src/iguana/algorithms/physics/Tools.h b/src/iguana/algorithms/physics/Tools.h new file mode 100644 index 00000000..a25686f8 --- /dev/null +++ b/src/iguana/algorithms/physics/Tools.h @@ -0,0 +1,49 @@ +/// @file Tools.h + +#include +#include + +namespace iguana::physics::tools { + + /// a value used when some calculation fails + double const UNDEF{-10000}; + + /// @brief calculate the angle between two planes + /// + /// The two planes are transverse to @latex{\vec{v}_a\times\vec{v}_b} and @latex{\vec{v}_c\times\vec{v}_d} + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @param v_c vector @latex{\vec{v}_c} + /// @param v_d vector @latex{\vec{v}_d} + /// @returns the angle between the planes, in radians, if the calculation is successful + std::optional PlaneAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b, + ROOT::Math::XYZVector const v_c, + ROOT::Math::XYZVector const v_d); + + /// @brief projection of one vector onto another + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @returns the vector @latex{\vec{v}_a} projected onto vector @latex{\vec{v}_b}, if the calculation is successful + std::optional ProjectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b); + + /// @brief projection of one vector onto the plane transverse to another vector + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @returns the vector @latex{\vec{v}_a} projected onto the plane transverse to @latex{\vec{v}_b}, if the calculation is successful + std::optional RejectVector( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b); + + /// @brief calculate the angle between two vectors + /// @param v_a vector @latex{\vec{v}_a} + /// @param v_b vector @latex{\vec{v}_b} + /// @returns the angle between @latex{\vec{v}_a} and @latex{\vec{v}_b}, if the calculation is successful + std::optional VectorAngle( + ROOT::Math::XYZVector const v_a, + ROOT::Math::XYZVector const v_b); + +} diff --git a/src/iguana/services/YAMLReader.cc b/src/iguana/services/YAMLReader.cc index 9c232e1e..2f417e76 100644 --- a/src/iguana/services/YAMLReader.cc +++ b/src/iguana/services/YAMLReader.cc @@ -22,7 +22,7 @@ namespace iguana { /////////////////////////////////////////////////////////////////////////////// template - SCALAR YAMLReader::GetScalar(YAML::Node node) + std::optional YAMLReader::GetScalar(YAML::Node node) { if(node.IsDefined() && !node.IsNull()) { try { @@ -35,32 +35,32 @@ namespace iguana { m_log->Error("YAML Misc. Exception: {}", e.what()); } } - throw std::runtime_error("Failed `GetScalar`"); + return std::nullopt; } - template int YAMLReader::GetScalar(YAML::Node node); - template double YAMLReader::GetScalar(YAML::Node node); - template std::string YAMLReader::GetScalar(YAML::Node node); + template std::optional YAMLReader::GetScalar(YAML::Node node); + template std::optional YAMLReader::GetScalar(YAML::Node node); + template std::optional YAMLReader::GetScalar(YAML::Node node); /////////////////////////////////////////////////////////////////////////////// template - SCALAR YAMLReader::GetScalar(node_path_t node_path) + std::optional YAMLReader::GetScalar(node_path_t node_path) { for(auto const& [config, filename] : m_configs) { auto node = FindNode(config, node_path); if(node.IsDefined() && !node.IsNull()) return GetScalar(node); } - throw std::runtime_error("Failed `GetScalar`"); + return std::nullopt; } - template int YAMLReader::GetScalar(node_path_t node_path); - template double YAMLReader::GetScalar(node_path_t node_path); - template std::string YAMLReader::GetScalar(node_path_t node_path); + template std::optional YAMLReader::GetScalar(node_path_t node_path); + template std::optional YAMLReader::GetScalar(node_path_t node_path); + template std::optional YAMLReader::GetScalar(node_path_t node_path); /////////////////////////////////////////////////////////////////////////////// template - std::vector YAMLReader::GetVector(YAML::Node node) + std::optional> YAMLReader::GetVector(YAML::Node node) { if(node.IsDefined() && !node.IsNull() && node.IsSequence()) { try { @@ -76,27 +76,27 @@ namespace iguana { m_log->Error("YAML Misc. Exception: {}", e.what()); } } - throw std::runtime_error("Failed `GetVector`"); + return std::nullopt; } - template std::vector YAMLReader::GetVector(YAML::Node node); - template std::vector YAMLReader::GetVector(YAML::Node node); - template std::vector YAMLReader::GetVector(YAML::Node node); + template std::optional> YAMLReader::GetVector(YAML::Node node); + template std::optional> YAMLReader::GetVector(YAML::Node node); + template std::optional> YAMLReader::GetVector(YAML::Node node); /////////////////////////////////////////////////////////////////////////////// template - std::vector YAMLReader::GetVector(node_path_t node_path) + std::optional> YAMLReader::GetVector(node_path_t node_path) { for(auto const& [config, filename] : m_configs) { auto node = FindNode(config, node_path); if(node.IsDefined() && !node.IsNull()) return GetVector(node); } - throw std::runtime_error("Failed `GetVector`"); + return std::nullopt; } - template std::vector YAMLReader::GetVector(node_path_t node_path); - template std::vector YAMLReader::GetVector(node_path_t node_path); - template std::vector YAMLReader::GetVector(node_path_t node_path); + template std::optional> YAMLReader::GetVector(node_path_t node_path); + template std::optional> YAMLReader::GetVector(node_path_t node_path); + template std::optional> YAMLReader::GetVector(node_path_t node_path); /////////////////////////////////////////////////////////////////////////////// @@ -114,7 +114,7 @@ namespace iguana { auto bounds_node = sub_node[key]; if(bounds_node.IsDefined()) { auto bounds = GetVector(bounds_node); - if(bounds.size() == 2 && bounds[0] <= val && bounds[1] >= val) + if(bounds.value().size() == 2 && bounds.value()[0] <= val && bounds.value()[1] >= val) return sub_node; } } diff --git a/src/iguana/services/YAMLReader.h b/src/iguana/services/YAMLReader.h index 3748a359..9aa73e26 100644 --- a/src/iguana/services/YAMLReader.h +++ b/src/iguana/services/YAMLReader.h @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -37,27 +38,27 @@ namespace iguana { /// Read a scalar value from a `YAML::Node` /// @param node the `YAML::Node` to read - /// @return the scalar + /// @return the scalar, if found template - SCALAR GetScalar(YAML::Node node); + std::optional GetScalar(YAML::Node node); /// Read a scalar value from a `YAML::Node` path; searches all currently loaded config files. /// @param node_path the `YAML::Node` path - /// @return the scalar + /// @return the scalar, if found template - SCALAR GetScalar(node_path_t node_path); + std::optional GetScalar(node_path_t node_path); /// Read a vector value from a `YAML::Node` /// @param node the `YAML::Node` to read - /// @return the vector + /// @return the vector, if found template - std::vector GetVector(YAML::Node node); + std::optional> GetVector(YAML::Node node); /// Read a vector value from a `YAML::Node` path; searches all currently loaded config files. /// @param node_path the `YAML::Node` path - /// @return the vector + /// @return the vector, if found template - std::vector GetVector(node_path_t node_path); + std::optional> GetVector(node_path_t node_path); /// Create a function to search a `YAML::Node` for a sub-`YAML::Node` such that /// the scalar `val` is within a range specified by `key` diff --git a/src/iguana/tests/iguana_test.cc b/src/iguana/tests/iguana_test.cc index cad9a6c4..605dbfe3 100644 --- a/src/iguana/tests/iguana_test.cc +++ b/src/iguana/tests/iguana_test.cc @@ -26,7 +26,7 @@ int main(int argc, char** argv) auto exe = std::string(argv[0]); auto UsageCommands = [&](int exit_code) { - fmt::print(stderr, "\nUSAGE: {} [COMMAND] [OPTIONS]...\n", exe); + fmt::print("\nUSAGE: {} [COMMAND] [OPTIONS]...\n", exe); fmt::print("\n COMMANDS:\n\n"); fmt::print(" {:<20} {}\n", "algorithm", "call `Run` on an algorithm"); fmt::print(" {:<20} {}\n", "multithreading", "call `Run` on an algorithm, multithreaded"); @@ -44,12 +44,13 @@ int main(int argc, char** argv) return UsageCommands(2); command = std::string(argv[1]); if(command == "--help" || command == "-h") - return UsageCommands(2); + return UsageCommands(0); // omit the command, for getopt argv++; argc--; // usage options + // clang-format off auto UsageOptions = [&](int exit_code) { std::unordered_map> print_option = { @@ -115,39 +116,36 @@ int main(int argc, char** argv) { fmt::print(" {:<20} {}\n", "-v", "increase verbosity"); }}}; - std::vector available_options; - if(command == "algorithm" || command == "unit") { - available_options = {"f", "n", "a-algo", "b", "p"}; - } - if(command == "multithreading") { - available_options = {"f", "n", "a-algo", "b", "p", "j", "m", "V"}; - } - else if(command == "validator") { - available_options = {"f", "n", "a-vdor", "b", "o"}; - } - else if(command == "config") { - available_options = {"t"}; - } - else if(command == "logger") { - available_options = {}; + std::map> available_options = { + {"algorithm", {"f", "n", "a-algo", "b", "p"}}, + {"unit", {"f", "n", "a-algo", "b", "p"}}, + {"multithreading", {"f", "n", "a-algo", "b", "p", "j", "m", "V"}}, + {"validator", {"f", "n", "a-vdor", "b", "o"}}, + {"config", {"t"}}, + {"logger", {}} + }; + for(auto& it : available_options) + it.second.push_back("v"); + if(auto it{available_options.find(command)}; it != available_options.end()) { + fmt::print("\nUSAGE: {} {} [OPTIONS]...\n", exe, command); + fmt::print("\n OPTIONS:\n\n"); + for(auto available_opt : it->second) { + print_option.at(available_opt)(); + fmt::print("\n"); + } + return exit_code; } else { fmt::print(stderr, "ERROR: unknown command '{}'\n", command); return 1; } - available_options.push_back("v"); - fmt::print(stderr, "\nUSAGE: {} {} [OPTIONS]...\n", exe, command); - fmt::print("\n OPTIONS:\n\n"); - for(auto available_opt : available_options) { - print_option.at(available_opt)(); - fmt::print("\n"); - } - return exit_code; }; - if(argc <= 2 && command != "logger") - return UsageOptions(2); - auto first_option = std::string(argv[2]); + // clang-format on + + auto first_option = argc >= 2 ? std::string(argv[1]) : ""; if(first_option == "--help" || first_option == "-h") + return UsageOptions(0); + if(argc <= 2 && command != "logger") return UsageOptions(2); // parse option arguments @@ -155,7 +153,7 @@ int main(int argc, char** argv) while((opt = getopt(argc, argv, "hf:n:a:b:p:t:j:m:Vo:v|")) != -1) { switch(opt) { case 'h': - return UsageOptions(2); + return UsageOptions(0); case 'f': data_file = std::string(optarg); break; @@ -192,7 +190,7 @@ int main(int argc, char** argv) verbose = true; break; default: - return UsageOptions(1); + return UsageOptions(2); } } diff --git a/src/iguana/tests/include/TestAlgorithm.h b/src/iguana/tests/include/TestAlgorithm.h index 1ccd7f06..be446b3d 100644 --- a/src/iguana/tests/include/TestAlgorithm.h +++ b/src/iguana/tests/include/TestAlgorithm.h @@ -62,15 +62,15 @@ inline int TestAlgorithm( return 1; } // print the banks, before and after - if(verbose) { - for(decltype(bank_names)::size_type it_bank = 0; it_bank < bank_names.size(); it_bank++) { - fmt::print("{:=^70}\n", fmt::format(" BEFORE: {} ", bank_names.at(it_bank))); - banks_before.at(it_bank).show(); - fmt::print("{:=^70}\n", fmt::format(" AFTER: {} ", bank_names.at(it_bank))); - banks_after.at(it_bank).show(); - fmt::print("\n"); - } - } + // if(verbose) { + // for(decltype(bank_names)::size_type it_bank = 0; it_bank < bank_names.size(); it_bank++) { + // fmt::print("{:=^70}\n", fmt::format(" BEFORE: {} ", bank_names.at(it_bank))); + // banks_before.at(it_bank).show(); + // fmt::print("{:=^70}\n", fmt::format(" AFTER: {} ", bank_names.at(it_bank))); + // banks_after.at(it_bank).show(); + // fmt::print("\n"); + // } + // } } // stop the algorithm