diff --git a/examples/ProcessHepMCTCSCombi.C b/examples/ProcessHepMCTCSCombi.C new file mode 100644 index 0000000..47bed4d --- /dev/null +++ b/examples/ProcessHepMCTCSCombi.C @@ -0,0 +1,125 @@ +#include "CommonDefines.h" +#include "HepMCElectro.h" +#include "KinematicsProcElectro.h" +#include "KineCalculation.h" +#include "Indicing.h" +#include "ElectronScatterKinematics.h" +#include "gammaN_2_Spin0Spin0SpinHalf.h" +#include "TCSKinematics.h" +#include + +/** + * @brief Example Script: TCS Analysis with Combinatorics and Missing Mass. + * Updated to use SetMesonParticles shortcut and CloneLinked. + */ +void ProcessHepMCTCSCombi(){ + + using namespace rad; + using namespace rad::consts::data_type; + + gBenchmark->Start("df"); + + // ================================================================================= + // 1. SETUP & INJECTION + // ================================================================================= + HepMCElectro hepmc{ + "hepmc3_tree", + "/w/work5/home/garyp/eic/Farm/data/EpIC_DDVCS_ee_18x275/rootfiles/18x275_ddvcs_edecay_hplus.root" + }; + hepmc.SetupMC(); + + // ================================================================================= + // 2. PARTICLE DEFINITIONS + // ================================================================================= + hepmc.SetBeamIonIndex(3); + hepmc.SetBeamElectronIndex(0); + //hepmc.SetScatElectronCandidates({1, 6}); + //hepmc.SetParticleCandidates("ele", {1, 6}); + hepmc.SetScatElectronIndex(1); + hepmc.SetParticleIndex("ele",6); + hepmc.SetParticleIndex("pos", 7); + hepmc.SetParticleIndex("pprime", 5); + + hepmc.MakeCombinations(); + + // ================================================================================= + // 3. KINEMATICS PROCESSOR (Standard Topology) + // ================================================================================= + KinematicsProcElectro kine{&hepmc, MC()}; + + // A. Define Particles FIRST (Topology Construction) + kine.Creator().Sum("gprime", {{"ele", "pos"}}); + + // Missing Mass: n_calc = (Beam + Target) - (ScatEle + Jpsi + pi+) + kine.Creator().Diff("pprime_calc",{ + {consts::BeamIon(), consts::BeamEle()}, + {"gprime", consts::ScatEle()} + }); + + // B. Define Groups NEXT (Lazy Definition) + // This uses the Processor's SetGroup wrapper, ensuring "Z" exists before the group is built. + kine.SetMesonParticles({"ele","pos"}); + kine.SetBaryonParticles({"pprime"}); + + // ================================================================================= + // 4. CALCULATIONS (Registered on Master) + // ================================================================================= + // These will be calculated for Master AND automatically copied to the Linked clone + + kine.Q2(); + kine.xbj(); + kine.y(); + kine.nu(); + kine.tau(); + kine.tauprime(); + + kine.RegisterCalc("GammaPol",rad::physics::ElS_PolVirtPhot); + kine.RegisterCalc("GammaPolCirc",rad::physics::ElS_CircPolVirtPhot); + + kine.RegisterCalc("CosThetaHel",rad::gn2s0s0s12::CosThetaHel); + kine.RegisterCalc("ThetaHel",rad::gn2s0s0s12::ThetaHel); + kine.RegisterCalc("PhiHel",rad::gn2s0s0s12::PhiHel); + + kine.CosThetaCM(); + kine.PhiCM(); + + kine.Mass("GMass", {"gprime"}); + kine.Mass2("Qp2", {"gprime"}); + kine.Mass2("MissMass2", + {consts::BeamIon(), consts::BeamEle()}, + {"gprime", consts::ScatEle(), "pprime"}); + + kine.RegisterCalc("t_top", rad::physics::TTop); + kine.RegisterCalc("t_bot", rad::physics::TBot); + kine.RegisterCalc("tp_top", rad::physics::TPrimeTop); + kine.RegisterCalc("tp_bot", rad::physics::TPrimeBot); + + kine.RegisterCalc("deltaT_top", rad::physics::DeltaTTop); + kine.RegisterCalc("deltaT_bot", rad::physics::DeltaTBot); + + // ================================================================================= + // 5. LINKED PROCESSOR (Cloned Hypothesis) + // ================================================================================= + + // 1. Clone: Copies all the calculations registered above (Q2, Phi, Mass, tp...) + // Creates a new processor for the "mc_" stream but with suffix "_ncalc" + // auto kine_miss = kine.CloneLinked("_ncalc"); + + // // 2. Customize: Override "Baryons" to be "n_calc" (missing neutron) for this hypothesis + // kine_miss->SetBaryonParticles({"n_calc"}); + + // ================================================================================= + // 6. INITIALIZATION & EXECUTION + // ================================================================================= + + // Init Master: Executes definitions and calculations for Standard Topology + kine.Init(); + + // Init Linked: Binds to 'kine' topology and executes copied calculations for Missing Topology + //kine_miss->InitLinked(kine); + + gBenchmark->Start("snapshot"); + hepmc.Snapshot("/w/work5/home/garyp/combirad_trees/HepMC_ddvcs_ee_18x275_hplus.root"); + gBenchmark->Stop("snapshot"); + gBenchmark->Print("snapshot"); +} diff --git a/include/AnalysisManager.h b/include/AnalysisManager.h index ffd1099..dbf05a4 100644 --- a/include/AnalysisManager.h +++ b/include/AnalysisManager.h @@ -59,7 +59,32 @@ namespace rad { * @param fileGlob Input file path/pattern. */ AnalysisManager(const std::string& name, const std::string& treeName, const std::string& fileGlob); - + + + /** + * @brief Constructor (multiple input files). + * @param name Name of the analysis (used for output filenames). + * @param treeName Input TTree name. + * @param files Vector of input file paths. + */ + AnalysisManager(const std::string& name, const std::string& treeName, const ROOT::RVec& files); + + + /** + * @brief Constructor (frome existing reaction). + * @param name Name of the analysis (used for output filenames). + * @param baseReaction Identifier of the base Reaction + */ + AnalysisManager(const std::string& name, const ReactionClass& baseReaction); + + + /** + * @brief Constructor clone method. + * @param newName Name of the new cloned analysis (used for output filenames). + * @param baseReaction Identifier of the base Reaction + */ + AnalysisManager Clone(const std::string& newName, bool copyStreams = false) const; + /** @brief Sets and creates the output directory. */ void SetOutputDir(const std::string& dir); @@ -78,13 +103,14 @@ namespace rad { */ void AddStream(const std::string& dataSource, const std::string& label = ""); + /** * @brief Shortcut to add standard default streams (label=""). * @details SetTypes("rec", "tru") creates streams named "rec" and "tru". * @param types Variadic list of stream prefixes (e.g. Rec(), Truth()). */ - template - void SetTypes(const std::string& name,Args... types); + // template + // void SetTypes(const std::string& name,Args... types); /** @return Reference to the underlying Reaction object. */ ReactionClass& Reaction(); @@ -145,7 +171,7 @@ namespace rad { * - Example: output/Y4260_rec_loose_Hist.root * @param suffix Suffix for the histogram file (default "Hist.root"). */ - void Run(const std::string& suffix = "Hist.root"); + void Run(const std::string& suffix = "_Hist.root"); /** * @brief Print comprehensive diagnostics for the entire analysis setup. @@ -181,13 +207,18 @@ namespace rad { fullName += lbl; outputSuffix = "_" + lbl; } + else{ + fullName.pop_back(); + //if no labels we dont want "rec_" + "_Tree.root" + //i.e. rec__Tree.root, rec__Hist.root + } // Ensure input prefix ends in "_" (e.g. "rec" -> "rec_") std::string inputPrefix = src; if(inputPrefix.back() != '_') inputPrefix += "_"; kine = std::make_unique(reaction, inputPrefix, outputSuffix); - sel = std::make_unique(*kine); + sel = std::make_unique(*kine); hist = std::make_unique(*kine, sel.get()); } }; @@ -211,17 +242,48 @@ namespace rad { // =========================================================================== template - inline AnalysisManager::AnalysisManager(const std::string& name, const std::string& treeName, const std::string& fileGlob) - : _reaction(treeName, fileGlob), _name(name) {} + inline AnalysisManager::AnalysisManager(const std::string& name, + const std::string& treeName, + const std::string& fileGlob) + : _reaction(treeName, fileGlob), _name(name) {} + + + template + inline AnalysisManager::AnalysisManager(const std::string& name, + const std::string& treeName, + const ROOT::RVec& files) + : _reaction(treeName, files), _name(name) {} template - inline void AnalysisManager::SetOutputDir(const std::string& dir) { - _outputDir = dir; - if (!_outputDir.empty() && !std::filesystem::exists(_outputDir)) { - std::filesystem::create_directories(_outputDir); - } + inline void AnalysisManager::SetOutputDir(const std::string& dir) { + _outputDir = dir; + if (!_outputDir.empty() && !std::filesystem::exists(_outputDir)) { + std::filesystem::create_directories(_outputDir); + } } - + + template + AnalysisManager::AnalysisManager(const std::string& name, + const R& baseReaction) + : _reaction(baseReaction), _name(name) + {}//if this works we probably want to have safety on being base i.e. nothing setup yet, and clearing triggers etc + + // --- Clone Management --- + + template + AnalysisManager AnalysisManager::Clone(const std::string& newName, + bool copyStreams) const + { + AnalysisManager out(newName, _reaction); + + if (copyStreams) + for (auto& [key, s] : _streams) + out.AddStream(s.source, s.label); + out.SetOutputDir(_outputDir); + return out; + } + + // --- Stream Management --- template @@ -241,12 +303,12 @@ namespace rad { if(_primaryStream.empty()) _primaryStream = key; } - template - template - inline void AnalysisManager::SetTypes(const std::string& name,Args... types) { - // Just add basic streams with empty labels - (AddStream(types, name), ...); - } + // template + // template + // inline void AnalysisManager::SetTypes(const std::string& name,Args... types) { + // // Just add basic streams with empty labels + // (AddStream(types, name), ...); + // } template inline R& AnalysisManager::Reaction() { return _reaction; } @@ -319,12 +381,17 @@ namespace rad { inline void AnalysisManager::Init() { if(_initialized) return; if(_primaryStream.empty()) throw std::runtime_error("[AnalysisManager] No streams defined! Call SetTypes() or AddStream()."); - + + /* // PASS 0: Ensure truth variable goes through prefix/suffix machinery */ + /* for(auto& [key, stream] : _streams) { */ + /* stream.kine->RegisterExternalVar(consts::TruthMatchedCombi()); */ + /* } */ + // PASS 1: Initialize Kinematics (Create Variables) for(auto& [key, stream] : _streams) { - stream.kine->Init(); - } - + stream.kine->Init(); + } + // PASS 2: Compile Selections (Create Masks) for(auto& [key, stream] : _streams) { if(stream.sel) stream.sel->Init(); @@ -381,7 +448,14 @@ namespace rad { // Always add the Event Counter cols.push_back("rdfentry_"); - + + //Add the isTruth flag for rec streams + if(stream.source.find(Rec()) == 0){ + cols.push_back(consts::TruthMatchedCombi()); + cols.push_back(Rec()+consts::TruthMatchedCombi()); + //cols.push_back(Rec()+consts::TruthMatchedCombi()+"_"+stream.label); + } + auto streamCols = CollectStreamColumns(*stream.kine); cols.insert(cols.end(), streamCols.begin(), streamCols.end()); @@ -410,7 +484,6 @@ namespace rad { std::cout << "[AnalysisManager] Snapshotting stream '" << stream.fullName << "' to " << finalPath << " (with Mask: " << mask << ")" << std::endl; - _reaction.BookSnapshotCombi(finalPath, "tree", cols, mask); } } diff --git a/include/BasicKinematics.h b/include/BasicKinematics.h index 1db62ac..498e62c 100644 --- a/include/BasicKinematics.h +++ b/include/BasicKinematics.h @@ -171,6 +171,20 @@ namespace rad { return psum.Pt(); } + template + inline ResultType_t FourVectorECalc(const RVecIndices &indices, const Tp &px, const Tp &py, const Tp &pz, const Tm &m) { + const auto& ipos = indices[0]; + const auto& ineg = indices[1]; + + // Note: Assuming error checking for invalid indices is performed in the wrapper or upstream. + + PxPyPzMVector psum(0,0,0,0); + SumFourVector(psum, ipos, px, py, pz, m); + SubtractFourVector(psum, ineg, px, py, pz, m); + + return psum.E(); + } + //---------------------- 3-Vector Components (RVec Ops) ---------------------- /** diff --git a/include/Combinatorics.h b/include/Combinatorics.h index 8c95f69..52ae48c 100644 --- a/include/Combinatorics.h +++ b/include/Combinatorics.h @@ -1,204 +1,204 @@ -/** - * @file Combinatorics.h - * @brief Combinatorial engines for generating particle pairs, triplets, etc. - */ +// /** +// * @file Combinatorics.h +// * @brief Combinatorial engines for generating particle pairs, triplets, etc. +// */ -#pragma once +// #pragma once -#include -#include -#include -#include -#include -#include +// #include +// #include +// #include +// #include +// #include +// #include -namespace rad { -namespace combinatorics { +// namespace rad { +// namespace combinatorics { - // Output Type: [ColumnIndex][CombinationIndex] - using Combis_t = ROOT::RVec>; +// // Output Type: [ColumnIndex][CombinationIndex] +// using Combis_t = ROOT::RVec>; - // ================================================================================= - // 1. STANDARD CARTESIAN PRODUCT (GenerateAllCombinations) - // ================================================================================= +// // ================================================================================= +// // 1. STANDARD CARTESIAN PRODUCT (GenerateAllCombinations) +// // ================================================================================= - namespace impl { - /** - * @brief Recursive helper to build the Cartesian product. - * @tparam I Current column index being processed. - * @tparam ColTypes Variadic types of the input index vectors. - */ - template - struct CartesianProduct { - static void Generate(const std::tuple& cols, - ROOT::RVec& current_indices, - Combis_t& result) { +// namespace impl { +// /** +// * @brief Recursive helper to build the Cartesian product. +// * @tparam I Current column index being processed. +// * @tparam ColTypes Variadic types of the input index vectors. +// */ +// template +// struct CartesianProduct { +// static void Generate(const std::tuple& cols, +// ROOT::RVec& current_indices, +// Combis_t& result) { - const auto& current_col = std::get(cols); +// const auto& current_col = std::get(cols); - // Loop over the candidates in the current column - for (const auto& idx : current_col) { - current_indices[I] = idx; +// // Loop over the candidates in the current column +// for (const auto& idx : current_col) { +// current_indices[I] = idx; - // Recurse to the next column - CartesianProduct::Generate(cols, current_indices, result); - } - } - }; - - // Base case: All columns processed, push to result - template - struct CartesianProduct { - static void Generate(const std::tuple&, - ROOT::RVec& current_indices, - Combis_t& result) { - for (std::size_t i = 0; i < sizeof...(ColTypes); ++i) { - result[i].push_back(current_indices[i]); - } - } - }; - } - - /** - * @brief Generates all possible combinations (Cartesian Product) of the input indices. - * @param args Variadic list of RVec columns (one per candidate). - * @return RVec> The combinations in SoA format. - */ - template - Combis_t GenerateAllCombinations(const Args&... args) { - // 1. Setup Output - const std::size_t N_cols = sizeof...(Args); - Combis_t result(N_cols); - - // 2. Prepare Recursion - auto cols_tuple = std::make_tuple(args...); - ROOT::RVec current_indices(N_cols); +// // Recurse to the next column +// CartesianProduct::Generate(cols, current_indices, result); +// } +// } +// }; + +// // Base case: All columns processed, push to result +// template +// struct CartesianProduct { +// static void Generate(const std::tuple&, +// ROOT::RVec& current_indices, +// Combis_t& result) { +// for (std::size_t i = 0; i < sizeof...(ColTypes); ++i) { +// result[i].push_back(current_indices[i]); +// } +// } +// }; +// } + +// /** +// * @brief Generates all possible combinations (Cartesian Product) of the input indices. +// * @param args Variadic list of RVec columns (one per candidate). +// * @return RVec> The combinations in SoA format. +// */ +// template +// Combis_t GenerateAllCombinations(const Args&... args) { +// // 1. Setup Output +// const std::size_t N_cols = sizeof...(Args); +// Combis_t result(N_cols); - impl::CartesianProduct<0, Args...>::Generate(cols_tuple, current_indices, result); +// // 2. Prepare Recursion +// auto cols_tuple = std::make_tuple(args...); +// ROOT::RVec current_indices(N_cols); - return result; - } +// impl::CartesianProduct<0, Args...>::Generate(cols_tuple, current_indices, result); +// cout<< " GenerateAllCombinations "<& indices, - const ROOT::RVec>& /*raw_groups_unused*/, - // Note: Logic below assumes groups are passed as vector> indices - const ROOT::RVec>& groups) { +// /** * @brief Checks if the current set of indices satisfies the symmetry constraints. +// * @details Ensures that for any defined symmetry group {A, B}, the values verify val(A) < val(B). +// */ +// inline bool CheckSymmetry(const ROOT::RVec& indices, +// const ROOT::RVec>& /*raw_groups_unused*/, +// // Note: Logic below assumes groups are passed as vector> indices +// const ROOT::RVec>& groups) { - for (const auto& group : groups) { - // Group contains column indices (e.g. {0, 1} for gamma1, gamma2) - // We strictly enforce indices[0] < indices[1] < ... to prevent swaps - for (size_t k = 0; k < group.size() - 1; ++k) { - int col_idx_a = group[k]; - int col_idx_b = group[k+1]; - if (indices[col_idx_a] >= indices[col_idx_b]) return false; - } - } - return true; - } +// for (const auto& group : groups) { +// // Group contains column indices (e.g. {0, 1} for gamma1, gamma2) +// // We strictly enforce indices[0] < indices[1] < ... to prevent swaps +// for (size_t k = 0; k < group.size() - 1; ++k) { +// int col_idx_a = group[k]; +// int col_idx_b = group[k+1]; +// if (indices[col_idx_a] >= indices[col_idx_b]) return false; +// } +// } +// return true; +// } - /** @brief Recursive generator with symmetry check hook. */ - template - struct SymmetricProduct { - static void Generate(const std::tuple& cols, - const ROOT::RVec>& groups, - ROOT::RVec& current_indices, - Combis_t& result) { +// /** @brief Recursive generator with symmetry check hook. */ +// template +// struct SymmetricProduct { +// static void Generate(const std::tuple& cols, +// const ROOT::RVec>& groups, +// ROOT::RVec& current_indices, +// Combis_t& result) { - const auto& current_col = std::get(cols); - for (const auto& idx : current_col) { - current_indices[I] = idx; - SymmetricProduct::Generate(cols, groups, current_indices, result); - } - } - }; - - template - struct SymmetricProduct { - static void Generate(const std::tuple&, - const ROOT::RVec>& groups, - ROOT::RVec& current_indices, - Combis_t& result) { - // Here is the "Symmetry Filter" - if (CheckSymmetry(current_indices, {}, groups)) { - for (std::size_t i = 0; i < sizeof...(ColTypes); ++i) { - result[i].push_back(current_indices[i]); - } - } - } - }; +// const auto& current_col = std::get(cols); +// for (const auto& idx : current_col) { +// current_indices[I] = idx; +// SymmetricProduct::Generate(cols, groups, current_indices, result); +// } +// } +// }; + +// template +// struct SymmetricProduct { +// static void Generate(const std::tuple&, +// const ROOT::RVec>& groups, +// ROOT::RVec& current_indices, +// Combis_t& result) { +// // Here is the "Symmetry Filter" +// if (CheckSymmetry(current_indices, {}, groups)) { +// for (std::size_t i = 0; i < sizeof...(ColTypes); ++i) { +// result[i].push_back(current_indices[i]); +// } +// } +// } +// }; - // --- Variadic Argument Peeler --- - // This machinery separates the parameter pack (columns) from the last argument (groups). +// // --- Variadic Argument Peeler --- +// // This machinery separates the parameter pack (columns) from the last argument (groups). - template struct SymHelper; +// template struct SymHelper; - // Recursive Step: Still have more than 2 args (Head + Tail...) - // We keep building the tuple of columns. - template - struct SymHelper { - static Combis_t Run(std::tuple head_tuple, const Tail&... tail) { - return SymHelper::RecurseAndRun(head_tuple, tail...); - } +// // Recursive Step: Still have more than 2 args (Head + Tail...) +// // We keep building the tuple of columns. +// template +// struct SymHelper { +// static Combis_t Run(std::tuple head_tuple, const Tail&... tail) { +// return SymHelper::RecurseAndRun(head_tuple, tail...); +// } - template - static Combis_t RecurseAndRun(std::tuple prev, const Head& current, const Tail&... rest) { - auto next_tuple = std::tuple_cat(prev, std::make_tuple(current)); - return SymHelper::RecurseAndRun(next_tuple, rest...); - } - }; - - // Base Case: Only 1 type left in the pack -> This MUST be the Groups vector - template - struct SymHelper { - template - static Combis_t RecurseAndRun(std::tuple cols, const Last& groups) { - // "Last" is expected to be ROOT::RVec> (or convertible) - // "cols" is the tuple of all RVec columns. +// template +// static Combis_t RecurseAndRun(std::tuple prev, const Head& current, const Tail&... rest) { +// auto next_tuple = std::tuple_cat(prev, std::make_tuple(current)); +// return SymHelper::RecurseAndRun(next_tuple, rest...); +// } +// }; + +// // Base Case: Only 1 type left in the pack -> This MUST be the Groups vector +// template +// struct SymHelper { +// template +// static Combis_t RecurseAndRun(std::tuple cols, const Last& groups) { +// // "Last" is expected to be ROOT::RVec> (or convertible) +// // "cols" is the tuple of all RVec columns. - const std::size_t N_cols = sizeof...(CollectedCols); - Combis_t result(N_cols); - ROOT::RVec current_indices(N_cols); +// const std::size_t N_cols = sizeof...(CollectedCols); +// Combis_t result(N_cols); +// ROOT::RVec current_indices(N_cols); - SymmetricProduct<0, CollectedCols...>::Generate(cols, groups, current_indices, result); - return result; - } - }; +// SymmetricProduct<0, CollectedCols...>::Generate(cols, groups, current_indices, result); +// return result; +// } +// }; - // Helper to extract the first argument and kickstart the peeling - // MOVED INSIDE IMPL NAMESPACE - template - Combis_t GenerateSymmetricCombinations_Launcher(const First& first, const Rest&... rest) { - return SymHelper::RecurseAndRun(std::make_tuple(first), rest...); - } +// // Helper to extract the first argument and kickstart the peeling +// // MOVED INSIDE IMPL NAMESPACE +// template +// Combis_t GenerateSymmetricCombinations_Launcher(const First& first, const Rest&... rest) { +// return SymHelper::RecurseAndRun(std::make_tuple(first), rest...); +// } - } // namespace impl - - /** - * @brief Generates combinations enforcing index ordering for symmetric groups. - * * usage: GenerateSymmetricCombinations(col1, col2, col3, {{0,1}, {2,3}}) - * @tparam Args... The columns followed by one vector>. - */ - template - Combis_t GenerateSymmetricCombinations(const Args&... args) { - // Start the recursive peeling process - // We assume at least 1 column + 1 group arg, so Args size >= 2 - static_assert(sizeof...(Args) >= 2, "GenerateSymmetricCombinations requires at least 1 column and 1 group vector."); +// } // namespace impl + +// /** +// * @brief Generates combinations enforcing index ordering for symmetric groups. +// * * usage: GenerateSymmetricCombinations(col1, col2, col3, {{0,1}, {2,3}}) +// * @tparam Args... The columns followed by one vector>. +// */ +// template +// Combis_t GenerateSymmetricCombinations(const Args&... args) { +// // Start the recursive peeling process +// // We assume at least 1 column + 1 group arg, so Args size >= 2 +// static_assert(sizeof...(Args) >= 2, "GenerateSymmetricCombinations requires at least 1 column and 1 group vector."); - // We need to unpack the first argument to start the tuple - return impl::GenerateSymmetricCombinations_Launcher(args...); - } +// // We need to unpack the first argument to start the tuple +// return impl::GenerateSymmetricCombinations_Launcher(args...); +// } -} // namespace combinatorics -} // namespace rad +// } // namespace combinatorics +// } // namespace rad // #pragma once @@ -334,3 +334,271 @@ namespace combinatorics { // } // namespace combinatorics // } // namespace rad +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include "DefineNames.h" +#include "Constants.h" +#include "CommonDefines.h" + +/** + * @file Combinatorics.h + * @brief Utilities for generating particle combinations. + * @details + * This file contains the logic to generate the "Cartesian Product" of all particle candidates. + * It converts a list of candidates into a set of unique combinatorial events. + */ + +namespace rad { + namespace combinatorics { + + using ROOT::RVecI; + using ROOT::RVec; + using std::string; + using rad::consts::InvalidEntry; + + //--------------------------------------------------------- + // Core Combinatorial Logic Helper + //--------------------------------------------------------- + + /** + * @brief Generates all valid, unique combinations of candidates. + * @details + * Output is a "Structure of Arrays" (SoA). + * `result[particle_role_index]` is a vector containing the candidate index for that role + * for every valid combination. + * * @param candidates_vec A vector of candidate lists. + * @return RVecIndices The structure of valid combinations. + */ + inline RVecIndices GenerateAllCombinations(const RVecIndices& candidates_vec) { + const size_t n_particles = candidates_vec.size(); + if (n_particles == 0) { + return RVecIndices(n_particles); + } + + // --- 1. Initial Setup and Pre-calculation --- + ROOT::RVec max_indices; + max_indices.reserve(n_particles); + + size_t total_combos = 1; + for (const auto& vec : candidates_vec) { + if (vec.empty()) return RVecIndices(n_particles); + max_indices.push_back(vec.size()); + total_combos *= vec.size(); + } + + // The N-dimensional counter (state of the current permutation) + ROOT::RVec current_indices(n_particles, 0); + + // --- 2. Initialize Transposed Output Structure --- + RVecIndices result_by_particle(n_particles); + for (size_t i = 0; i < n_particles; ++i) { + result_by_particle[i].reserve(total_combos); + } + + // --- 3. Iterative Combinatorial Loop --- + ROOT::RVec current_combination_buffer(n_particles); + + for (size_t i_combo = 0; i_combo < total_combos; ++i_combo) { + + bool is_unique = true; + + // --- A. Build Combination and Check Uniqueness --- + for (size_t i_particle = 0; i_particle < n_particles; ++i_particle) { + + auto particle_index = candidates_vec[i_particle][current_indices[i_particle]]; + + // Only check for duplicates if it's a valid physical track/cluster + if (!consts::IsInvalidEntry(particle_index )) { + auto start_it = current_combination_buffer.begin(); + auto end_it = start_it + i_particle; + + // Fast linear search over the tiny buffer + if (std::find(start_it, end_it, particle_index) != end_it) { + is_unique = false; + break; // Stop checking this combination immediately + } + } + + current_combination_buffer[i_particle] = particle_index; + } + + // --- B. Record Valid Combination --- + if (is_unique) { + for (size_t i_particle = 0; i_particle < n_particles; ++i_particle) { + result_by_particle[i_particle].push_back(current_combination_buffer[i_particle]); + } + } + + // --- C. Increment the N-dimensional counter --- + size_t k = 0; + while (k < n_particles) { + current_indices[k]++; + if (current_indices[k] < max_indices[k]) { + break; + } + current_indices[k] = 0; + k++; + } + } + return result_by_particle; + } + + } // namespace combinatorics +} // namespace rad +// #pragma once + +// #include +// #include +// #include +// #include +// #include "DefineNames.h" +// #include "Constants.h" +// #include "CommonDefines.h" + +// namespace rad { +// namespace combinatorics { + +// using Indices_t = ROOT::RVec; +// using RVecIndices = ROOT::RVec; +// using rad::consts::InvalidEntry; + +// // ================================================================================= +// // 1. CORE ITERATIVE ENGINE +// // ================================================================================= + +// /** +// * @brief Generates all valid, unique combinations, optionally applying symmetry constraints. +// * @param candidates_vec RVec of candidate lists for each role. +// * @param groups Optional RVec of groups enforcing indistinguishable permutations. +// * @return RVecIndices The structure of valid combinations (SoA format). +// */ +// inline RVecIndices GenerateCombinationsImpl( +// const RVecIndices& candidates_vec, +// const RVecIndices& groups = {}) +// { +// const size_t n_particles = candidates_vec.size(); +// if (n_particles == 0) return RVecIndices(0); + +// // --- A. Setup & Pre-calculation --- +// ROOT::RVec max_indices; +// max_indices.reserve(n_particles); +// size_t total_combos = 1; + +// for (const auto& vec : candidates_vec) { +// if (vec.empty()) return RVecIndices(n_particles); +// max_indices.push_back(vec.size()); +// total_combos *= vec.size(); +// } + +// ROOT::RVec current_indices(n_particles, 0); +// RVecIndices result_by_particle(n_particles); + +// // Reserving max possible size to prevent reallocation overhead +// for (size_t i = 0; i < n_particles; ++i) { +// result_by_particle[i].reserve(total_combos); +// } + +// Indices_t current_combination_buffer(n_particles); + +// // --- B. Iterative Odometer Loop --- +// for (size_t i_combo = 0; i_combo < total_combos; ++i_combo) { +// bool is_valid = true; + +// // 1. Build combination and check global uniqueness +// for (size_t i_particle = 0; i_particle < n_particles; ++i_particle) { +// int particle_index = candidates_vec[i_particle][current_indices[i_particle]]; + +// // Skip uniqueness check for dummies, but still record them +// if (!consts::IsInvalidEntry(particle_index)) { +// auto start_it = current_combination_buffer.begin(); +// auto end_it = start_it + i_particle; +// if (std::find(start_it, end_it, particle_index) != end_it) { +// is_valid = false; +// break; // Duplicate track used for multiple roles +// } +// } +// current_combination_buffer[i_particle] = particle_index; +// } + +// // 2. Symmetry Check (Prevents double counting e.g., gamma gamma) +// if (is_valid && !groups.empty()) { +// for (const auto& group : groups) { +// for (size_t k = 0; k < group.size() - 1; ++k) { +// if (current_combination_buffer[group[k]] >= current_combination_buffer[group[k+1]]) { +// is_valid = false; +// break; +// } +// } +// if (!is_valid) break; +// } +// } + +// // 3. Save Valid Combination +// if (is_valid) { +// for (size_t i_particle = 0; i_particle < n_particles; ++i_particle) { +// result_by_particle[i_particle].push_back(current_combination_buffer[i_particle]); +// } +// } + +// // 4. Increment the N-dimensional counter (Odometer) +// size_t k = 0; +// while (k < n_particles) { +// current_indices[k]++; +// if (current_indices[k] < max_indices[k]) break; +// current_indices[k] = 0; +// k++; +// } +// } +// return result_by_particle; +// } + +// // ================================================================================= +// // 2. RDATAFRAME INTERFACES +// // ================================================================================= + +// /** +// * @brief Interface for asymmetric combinations. +// * Connects to: Define(..., Form("GenerateAllCombinations({%s})", ...)) +// */ +// inline RVecIndices GenerateAllCombinations(const RVecIndices& candidates_vec) { +// return GenerateCombinationsImpl(candidates_vec); +// } + +// namespace impl { +// // Pack the first N-1 variadic arguments into an RVecIndices +// template +// RVecIndices PackArgsToRVec(const Tuple& t, std::index_sequence) { +// return RVecIndices{ std::get(t)... }; +// } +// } + +// /** +// * @brief Interface for symmetric combinations handling variadic inputs. +// * Connects to: Define(..., Form("GenerateSymmetricCombinations(%s, %s)", ...)) +// */ +// template +// RVecIndices GenerateSymmetricCombinations(const Args&... args) { +// static_assert(sizeof...(Args) >= 2, "Requires at least 1 column and 1 symmetry group."); + +// // Tie all arguments together +// auto all_args = std::tie(args...); +// constexpr std::size_t N = sizeof...(Args); + +// // The last argument is the symmetry groups RVec +// const auto& groups = std::get(all_args); + +// // Pack all preceding arguments (the columns) into RVecIndices +// RVecIndices candidates_vec = impl::PackArgsToRVec(all_args, std::make_index_sequence{}); + +// return GenerateCombinationsImpl(candidates_vec, groups); +// } + +// } // namespace combinatorics +// } // namespace rad diff --git a/include/ConfigReaction.h b/include/ConfigReaction.h index 0220032..5982598 100644 --- a/include/ConfigReaction.h +++ b/include/ConfigReaction.h @@ -3,6 +3,7 @@ * @brief Central configuration manager for RAD analysis. */ + #pragma once #include "RDFInterface.h" @@ -61,8 +62,8 @@ namespace rad { /** @brief Get the standard name for the matching column (e.g. "rec_match_id"). */ virtual std::string GetMatchName(const std::string& type) const { - // Standard convention: prefix + "match_id" - return type + "match_id"; + // Standard convention: prefix + consts::NameMatchId() + return type + consts::NameMatchId(); } /** @brief Define Signal Combination flag using standard naming. @@ -77,9 +78,16 @@ namespace rad { // define a default "Always True" signal flag so snapshots don't break. std::cout<<"[ConfigReaction] DefineTrueMatchCombi " < void SetSymmetryParticles(Args... args) { @@ -217,7 +225,7 @@ namespace rad { const std::string& name = match.first; int role = match.second; - std::string flagName = name + "_is_true";// + DoNotWriteTag(); + std::string flagName = name + "_is_true";// + DoNotWriteTag(); std::string colName = type + name; // e.g. "rec_ele" // Check: pIndices[i] is the candidate index for the i-th combination. @@ -238,17 +246,50 @@ namespace rad { }, {colName, matchIdCol}); - + + if (!logic.empty()) logic += " && "; logic += flagName; } if(logic.empty()) logic = "1"; Define(consts::TruthMatchedCombi(), logic); + Define(type + consts::TruthMatchedCombi(), consts::TruthMatchedCombi()); } + inline void ConfigReaction::DefineTruePID(const std::string& prefix) { + + std::string match_col = prefix + consts::NameMatchId(); + std::string truth_pid_col = consts::data_type::Truth() + consts::NamePid(); + + // Define the output using standard naming from DefineNames.h (if added), + // or just hardcode "true_pid" if you prefer. + std::string out_col = prefix + "true_pid"; + + // 1. Strict Robustness Checks + if (!ColumnExists(match_col)) { + throw std::runtime_error("[DefineTruePID] ERROR: Required matching column '" + match_col + "' does not exist. Did you forget to call SetupMatching()?"); + } + if (!ColumnExists(truth_pid_col)) { + throw std::runtime_error("[DefineTruePID] ERROR: Required truth PID column '" + truth_pid_col + "' does not exist. Verify SetupTruth() execution."); + } - // ... [Rest of implementation remains unchanged] ... - inline void ConfigReaction::SetParticleCandidatesExpr(const std::string& name, const std::string& type, const std::string& expression) { + // 2. Define the Column safely + Define(out_col, + [](const rad::Indices_t& match_id, const rad::Indices_t& tru_pdg) { + rad::Indices_t tpid(match_id.size(), 0); + for(size_t i = 0; i < match_id.size(); ++i) { + // Safe access: ensure match_id is valid and within bounds of the truth array + if(match_id[i] >= 0 && match_id[i] < (int)tru_pdg.size()) { + tpid[i] = tru_pdg[match_id[i]]; + } + } + return tpid; + }, + {match_col, truth_pid_col} + ); + } + + inline void ConfigReaction::SetParticleCandidatesExpr(const std::string& name, const std::string& type, const std::string& expression) { ValidateType(type); if (_typeCandidateExpressions[type].count(name) || _typeLambdaDependencies[type].count(name)) { throw std::invalid_argument("Candidate '" + name + "' already defined for type '" + type + "'."); @@ -321,7 +362,7 @@ namespace rad { if(colList.size() >= 2 && colList.front() == '{' && colList.back() == '}') colList = colList.substr(1, colList.size() - 2); if (_symmetryGroups.empty()) { - Define(comboColName, Form("rad::combinatorics::GenerateAllCombinations(%s)", colList.c_str())); + Define(comboColName, Form("rad::combinatorics::GenerateAllCombinations({%s})", colList.c_str())); } else { ROOT::RVec groupStrs; for(const auto& group : _symmetryGroups) { diff --git a/include/Constants.h b/include/Constants.h index 274e889..cbf479f 100644 --- a/include/Constants.h +++ b/include/Constants.h @@ -1,6 +1,9 @@ #pragma once #include +#include +#include + namespace rad{ namespace consts{ @@ -64,10 +67,16 @@ namespace rad{ return std::numeric_limits::max(); } - template + template inline bool IsInvalidEntry(const T& entry){ - return entry==InvalidEntry(); + if constexpr (std::is_floating_point_v) { + return std::isnan(entry); + } else { + return entry == InvalidEntry(); + } } + + // constexpr double InvalidEntry(){return std::numeric_limits::quiet_NaN();}; constexpr int InvalidIndex(){return -1;}; diff --git a/include/DefineNames.h b/include/DefineNames.h index b6aa250..35b7576 100644 --- a/include/DefineNames.h +++ b/include/DefineNames.h @@ -102,7 +102,9 @@ namespace rad{ const std::string NameM() {return "m";} const std::string NamePid() {return "pid";} - + // Standardized matching columns + const std::string NameMatchId() {return "match_id";} + const std::string NameTruePid() {return "true_pid";} /** * Types of data */ diff --git a/include/ElectroIonReaction.h b/include/ElectroIonReaction.h index 73ef53e..493403b 100644 --- a/include/ElectroIonReaction.h +++ b/include/ElectroIonReaction.h @@ -208,7 +208,7 @@ namespace rad { PxPyPzMVector _p4ion_beam; ///< Internal storage for Ion Beam P4 - bool _useBeamsFromMC = true; ///< Flag to determine if beams should be read from MC + bool _useBeamsFromMC = false; ///< Flag to determine if beams should be read from MC private: @@ -285,6 +285,7 @@ namespace rad { inline void ElectroIonReaction::SetMCBeamIndices(int eleIdx, int ionIdx) { _mcBeamEleIdx = eleIdx; _mcBeamIonIdx = ionIdx; + _useBeamsFromMC = true; } inline void ElectroIonReaction::SetBeamElectronColumns(const std::string& px, const std::string& py, diff --git a/include/ElectronScatterKinematics.h b/include/ElectronScatterKinematics.h index 7eb6bd6..e3960e0 100644 --- a/include/ElectronScatterKinematics.h +++ b/include/ElectronScatterKinematics.h @@ -23,6 +23,7 @@ namespace rad { * @brief Result structure for decay angles. */ struct DecayAngles_t { + double theta = 0; ///< Scattering angle in the proton rest frame (radians). double cosTheta = 0.; ///< Cosine of the decay angle (polar angle). double phi = 0.; ///< Azimuthal angle. }; @@ -43,7 +44,87 @@ namespace rad { auto phot = PhotoFourVector(react, px, py, pz, m); return -phot.M2(); } - + + /** + * @brief Calculates the struck quark momentum fraction, xbj = Q^2/(2p.q). + * @param react The fixed ReactionMap. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return ResultType_t The xbj value (always positive). + */ + inline ResultType_t ElS_xbj(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); + auto phot = PhotoFourVector(react, px, py, pz, m); + return -phot.M2() / (2*ion.Dot(phot)); + } + + /** + * @brief Calculates the inelasticity parameter, y = p.q/p.k + * @param react The fixed ReactionMap. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return ResultType_t The y value (always positive). + */ + inline ResultType_t ElS_y(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); + auto ebeam = FourVector(react[consts::OrderBeams()][consts::OrderBeamEle()], px, py, pz, m); + auto phot = PhotoFourVector(react, px, py, pz, m); + return ion.Dot(phot) / ion.Dot(ebeam); + } + + /** + * @brief Calculates the virtual photon energy (ion rest frame), nu = p.q / M + * @param react The fixed ReactionMap. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return ResultType_t The nu value (always positive). + */ + inline ResultType_t ElS_nu(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); + auto phot = PhotoFourVector(react, px, py, pz, m); + return ion.Dot(phot) / ion.M(); + } + + /** + * @brief Calculates the electroproduction variable, tau = Q^2/(4M^2). + * @param react The fixed ReactionMap. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return ResultType_t The tau value (always positive). + */ + inline ResultType_t ElS_tau(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); + auto phot = PhotoFourVector(react, px, py, pz, m); + //auto Q2 = -phot.M2(); //call ElS_Q2() here or just compute by hand? + return -phot.M2() / (4*ion.M2()); + } + + /** + * @brief Calculates the variable, tau' = Q'^2/(2p.q). + * @param react The fixed ReactionMap. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return ResultType_t The tau' value (always positive). + */ + inline ResultType_t ElS_tauprime(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + // Q^2 is defined as the negative invariant mass squared of the virtual photon. + auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); + auto phot = PhotoFourVector(react, px, py, pz, m); + auto mes = FourVector(react[consts::OrderMesons()], px, py, pz, m); // Assumes single particle or combined meson vector + return mes.M2() / (2*ion.Dot(phot)); + } + + /** * @brief Calculates the four-momentum of the initial state Center-of-Mass (CM) system. * @param react The fixed ReactionMap. @@ -104,7 +185,7 @@ namespace rad { return {theta, energy, Q2}; } - + /** * @brief Calculates the polarization parameter (epsilon) for the virtual photon. * @param react The fixed ReactionMap. @@ -128,6 +209,20 @@ namespace rad { return pol; } + + /** + * @brief Calculates the circular polarization parameter for the virtual photon sqrt(1-epsilon^2). + * @param react The fixed ReactionMap. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return ResultType_t The circular polarization parameter. + */ + inline ResultType_t ElS_CircPolVirtPhot(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + auto eps = ElS_PolVirtPhot(react,px,py,pz,m); + return sqrt(1-eps*eps); + } // --- Decay Frame Kinematics (CM and Proton Rest) --- @@ -229,9 +324,9 @@ namespace rad { XYZVector angles(prMes.Vect().Dot(xV), prMes.Vect().Dot(yV), prMes.Vect().Dot(zV)); DecayAngles_t result; + result.theta = (angles.Theta()); result.cosTheta = TMath::Cos(angles.Theta()); result.phi = angles.Phi(); - return result; } @@ -246,6 +341,17 @@ namespace rad { { return ElectroProtonRestDecay(react, px, py, pz, m).cosTheta; } + /** + * @brief Calculates the decay angle (Theta) in the Proton Rest frame. + * @param react, px, py, pz, m The inputs. + * @return ResultType_t polar decay angle. + */ + inline ResultType_t ElS_ThetaProtonRest(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + return ElectroProtonRestDecay(react, px, py, pz, m).theta; + } /** * @brief Calculates the decay angle Phi in the Proton Rest frame. diff --git a/include/Histogrammer.h b/include/Histogrammer.h index 873e78b..8f50608 100644 --- a/include/Histogrammer.h +++ b/include/Histogrammer.h @@ -58,7 +58,7 @@ #include namespace rad { -namespace histo { + namespace histo { /** * @struct HistoDef @@ -67,12 +67,25 @@ namespace histo { * during the Init phase (after Kinematics are defined). */ struct HistoDef { - std::string name; - std::string title; - int nbins; - double min; - double max; - std::string varBaseName; + std::string name; + std::string title; + int nbins; + double min; + double max; + std::string varBaseName; + }; + + struct HistoDef2D { + std::string name; + std::string title; + int nbinsX; + double xmin; + double xmax; + int nbinsY; + double ymin; + double ymax; + std::string XvarBaseName; + std::string YvarBaseName; }; /** @@ -80,11 +93,11 @@ namespace histo { * @brief Configuration for a single splitting axis (category). */ struct SplitDef { - std::string name; ///< Human-readable name (e.g., "Sector") - std::string baseCol; ///< Column name in Processor (e.g., "loc_sector") - int nbins; - double min; - double max; + std::string name; ///< Human-readable name (e.g., "Sector") + std::string baseCol; ///< Column name in Processor (e.g., "loc_sector") + int nbins; + double min; + double max; }; /** @@ -93,79 +106,92 @@ namespace histo { */ class Histogrammer { public: - /** - * @brief Constructor binding to a specific Processor stream. - * @param proc The processor (Rec, Sim, etc.) used to resolve variable names. - * @param sel Optional pointer to PhysicsSelection. If provided, histograms will - * automatically filter using the selection's compiled mask. - */ - Histogrammer(KinematicsProcessor& proc, PhysicsSelection* sel = nullptr); - - /** - * @brief Defines a split axis (category) for subsequent histograms. - * @details All histograms created after calling this will have this extra dimension. - * This allows you to easily categorize all plots by Sector, Helicity, etc. - * @param name Axis title (e.g. "Sector"). - * @param baseCol Variable name in RDF (e.g. "loc_sector"). Resolved via Processor. - * @param nbins Number of bins. - * @param min Axis min. - * @param max Axis max. - */ - void AddSplit(const std::string& name, const std::string& baseCol, int nbins, double min, double max); - - /** - * @brief queues a multidimensional histogram (THnSparse) for booking. - * @details - * Stores the configuration. Actual RDF Booking happens in Init(). - * @param name Output histogram base name. - * @param title Histogram title. - * @param nbins Number of bins for the main variable. - * @param min Minimum value for the main variable. - * @param max Maximum value for the main variable. - * @param varBaseName The RDF column name (e.g. "mm2"). Resolved via Processor. - */ - void Create(const std::string& name, const std::string& title, - int nbins, double min, double max, - const std::string& varBaseName); - - /** - * @brief Compiles definitions and Books actions in RDataFrame. - * @details MUST be called after PhysicsSelection::Compile(). - */ - void Init(); - - /** - * @brief Writes all histograms to a ROOT file. - * @details - * Performs "Unfolding": Recursively projects the N-D histograms into 1D TH1Ds - * for every bin combination of the split axes. - * @param filename Output filename. - * @param option ROOT File option ("RECREATE" or "UPDATE"). - */ - void File(const std::string& filename, const std::string& option = "RECREATE"); + /** + * @brief Constructor binding to a specific Processor stream. + * @param proc The processor (Rec, Sim, etc.) used to resolve variable names. + * @param sel Optional pointer to PhysicsSelection. If provided, histograms will + * automatically filter using the selection's compiled mask. + */ + Histogrammer(KinematicsProcessor& proc, PhysicsSelection* sel = nullptr); + + /** + * @brief Defines a split axis (category) for subsequent histograms. + * @details All histograms created after calling this will have this extra dimension. + * This allows you to easily categorize all plots by Sector, Helicity, etc. + * @param name Axis title (e.g. "Sector"). + * @param baseCol Variable name in RDF (e.g. "loc_sector"). Resolved via Processor. + * @param nbins Number of bins. + * @param min Axis min. + * @param max Axis max. + */ + void AddSplit(const std::string& name, const std::string& baseCol, int nbins, double min, double max); + + /** + * @brief queues a multidimensional histogram (THnSparse) for booking. + * @details + * Stores the configuration. Actual RDF Booking happens in Init(). + * @param name Output histogram base name. + * @param title Histogram title. + * @param nbins Number of bins for the main variable. + * @param min Minimum value for the main variable. + * @param max Maximum value for the main variable. + * @param varBaseName The RDF column name (e.g. "mm2"). Resolved via Processor. + */ + void Create(const std::string& name, const std::string& title, + int nbins, double min, double max, + const std::string& varBaseName); + + + void Create2D(const std::string& name, const std::string& title, + int nbinsX, double xmin, double xmax, + int nbinsY, double ymin, double ymax, + const std::string& XvarBaseName, + const std::string& YvarBaseName); + + /** + * @brief Compiles definitions and Books actions in RDataFrame. + * @details MUST be called after PhysicsSelection::Compile(). + */ + void Init(); + + /** + * @brief Writes all histograms to a ROOT file. + * @details + * Performs "Unfolding": Recursively projects the N-D histograms into 1D TH1Ds + * for every bin combination of the split axes. + * @param filename Output filename. + * @param option ROOT File option ("RECREATE" or "UPDATE"). + */ + void File(const std::string& filename, const std::string& option = "RECREATE"); private: - KinematicsProcessor& _proc; - ConfigReaction& _rad; // Underlying RDF Manager - PhysicsSelection* _sel = nullptr; // Store pointer to get mask later - std::string _maskCol; - bool _initialized = false; - - ROOT::RVec _splits; - ROOT::RVec _defs; // Queue for lazy init + KinematicsProcessor& _proc; + ConfigReaction& _rad; // Underlying RDF Manager + PhysicsSelection* _sel = nullptr; // Store pointer to get mask later + std::string _maskCol; + bool _initialized = false; - // Store RResultPtrs to ensure RDF executes the actions - ROOT::RVec> _results; - - /** - * @brief Internal helper to actually book the histogram in RDF. - */ - void BookInternal(const HistoDef& def); - - /** - * @brief Recursive helper to project N-D histogram to 1D slices. - */ - void UnfoldAndWrite(THnSparseD* hn, TDirectory* dir, std::string current_suffix = "", int axis_depth = 1); + ROOT::RVec _splits; + ROOT::RVec _defs; // Queue for lazy init + ROOT::RVec _defs2D; // Queue for lazy init + + // Store RResultPtrs to ensure RDF executes the actions + ROOT::RVec> _results; + + /** + * @brief Internal helper to actually book the histogram in RDF. + */ + void BookInternal(const HistoDef& def); + + /** + * @brief Internal helper to actually book the 2D histogram in RDF. + */ + void BookInternal2D(const HistoDef2D& def); + + /** + * @brief Recursive helper to project N-D histogram to 1D slices. + */ + void UnfoldAndWrite(THnSparseD* hn, TDirectory* dir, std::string current_suffix = "", int axis_depth = 1); }; // ================================================================================= @@ -173,159 +199,232 @@ namespace histo { // ================================================================================= inline Histogrammer::Histogrammer(KinematicsProcessor& proc, PhysicsSelection* sel) - : _proc(proc), _rad(*proc.Reaction()), _maskCol("") - { + : _proc(proc), _rad(*proc.Reaction()), _maskCol("") + { if (sel) { - // We just store the POINTER to selection, we don't ask for mask yet - // The mask column might not exist until PhysicsSelection::Compile() runs - _sel = sel; + // We just store the POINTER to selection, we don't ask for mask yet + // The mask column might not exist until PhysicsSelection::Compile() runs + _sel = sel; } - } + } inline void Histogrammer::AddSplit(const std::string& name, const std::string& baseCol, int nbins, double min, double max) { - _splits.push_back({name, baseCol, nbins, min, max}); + _splits.push_back({name, baseCol, nbins, min, max}); } - inline void Histogrammer::Create(const std::string& name, const std::string& title, - int nbins, double min, double max, - const std::string& varBaseName) + inline void Histogrammer::Create(const std::string& name, const std::string& title, int nbins, double min, double max, const std::string& varBaseName) { - // Lazy: Just store the request. Init() will execute it. - _defs.push_back({name, title, nbins, min, max, varBaseName}); + // Lazy: Just store the request. Init() will execute it. + _defs.push_back({name, title, nbins, min, max, varBaseName}); } - + + inline void Histogrammer::Create2D(const std::string& name, const std::string& title, int nbinsX, double xmin, double xmax, int nbinsY, double ymin, double ymax, const std::string& XvarBaseName, const std::string& YvarBaseName) + { + _defs2D.push_back({name, title, + nbinsX, xmin, xmax, + nbinsY, ymin, ymax, + XvarBaseName, YvarBaseName}); + } + inline void Histogrammer::Init() { - if(_initialized) return; - - // 1. Resolve Mask (Now it guarantees to exist if compiled) - if(_sel) _maskCol = _sel->GetMaskColumn(); - - // 2. Book Histograms - for(const auto& def : _defs) { - BookInternal(def); - } - _initialized = true; + if (_initialized) return; + + if (_sel) + _maskCol = _sel->GetMaskColumn(); + + // 1D histos + for (const auto& def : _defs) + BookInternal(def); + + // 2D histos + for (const auto& def2 : _defs2D) + BookInternal2D(def2); + + _initialized = true; } - + inline void Histogrammer::BookInternal(const HistoDef& def) { - // 1. Resolve Column Names using Processor (Safe because Kine::Init ran) - std::string fullVarName = _proc.FullName(def.varBaseName); - ROOT::RVec cols; - cols.push_back(fullVarName); - - // 2. Setup Dimensions (1 Main Var + N Splits) - int n_dims = 1 + _splits.size(); - ROOT::RVec bins_vec(n_dims); - ROOT::RVec xmin(n_dims), xmax(n_dims); + // 1. Resolve Column Names using Processor (Safe because Kine::Init ran) + std::string fullVarName = _proc.FullName(def.varBaseName); + ROOT::RVec cols; + cols.push_back(fullVarName); + + // 2. Setup Dimensions (1 Main Var + N Splits) + int n_dims = 1 + _splits.size(); + ROOT::RVec bins_vec(n_dims); + ROOT::RVec xmin(n_dims), xmax(n_dims); - // Axis 0: The Main Variable - bins_vec[0] = def.nbins; xmin[0] = def.min; xmax[0] = def.max; + // Axis 0: The Main Variable + bins_vec[0] = def.nbins; xmin[0] = def.min; xmax[0] = def.max; - // Axis 1..N: The Splits - for(size_t i=0; i<_splits.size(); ++i) { - std::string fullSplitCol = _proc.FullName(_splits[i].baseCol); - cols.push_back(fullSplitCol); + // Axis 1..N: The Splits + for(size_t i=0; i<_splits.size(); ++i) { + std::string fullSplitCol = _proc.FullName(_splits[i].baseCol); + cols.push_back(fullSplitCol); - bins_vec[i+1] = _splits[i].nbins; - xmin[i+1] = _splits[i].min; - xmax[i+1] = _splits[i].max; - } - - // 3. Construct the THnSparse Prototype - std::string hNameFull = _proc.GetPrefix() + def.name + _proc.GetSuffix(); - auto hist = std::make_shared(hNameFull.c_str(), def.title.c_str(), - n_dims, bins_vec.data(), xmin.data(), xmax.data()); + bins_vec[i+1] = _splits[i].nbins; + xmin[i+1] = _splits[i].min; + xmax[i+1] = _splits[i].max; + } + + // 3. Construct the THnSparse Prototype + std::string hNameFull = _proc.GetPrefix() + def.name + _proc.GetSuffix(); + auto hist = std::make_shared(hNameFull.c_str(), def.title.c_str(), n_dims, bins_vec.data(), xmin.data(), xmax.data()); - // Label Axes - hist->GetAxis(0)->SetName(def.name.c_str()); - for(size_t i=0; i<_splits.size(); ++i) { - hist->GetAxis(i+1)->SetName(_splits[i].name.c_str()); - } - - // 4. Book the Action (Masked vs Unmasked) - // This assumes THnCombi is defined elsewhere. - if(!_maskCol.empty()) { - cols.push_back(_maskCol); - // Use for Masked execution. No data types needed in template. - THnCombi action(hist); - _results.push_back(_rad.CurrFrame().Book(std::move(action), utils::as_stdvector(cols))); - } - else { - // Use for Standard execution. - THnCombi action(hist); - _results.push_back(_rad.CurrFrame().Book(std::move(action), utils::as_stdvector(cols))); - } + // Label Axes + hist->GetAxis(0)->SetName(def.name.c_str()); + for(size_t i=0; i<_splits.size(); ++i) { + hist->GetAxis(i+1)->SetName(_splits[i].name.c_str()); + } + + // 4. Book the Action (Masked vs Unmasked) + // This assumes THnCombi is defined elsewhere. + if(!_maskCol.empty()) { + cols.push_back(_maskCol); + // Use for Masked execution. No data types needed in template. + THnCombi action(hist); + _results.push_back(_rad.CurrFrame().Book(std::move(action), utils::as_stdvector(cols))); + } + else { + // Use for Standard execution. + THnCombi action(hist); + _results.push_back(_rad.CurrFrame().Book(std::move(action), utils::as_stdvector(cols))); + } } - - inline void Histogrammer::File(const std::string& filename, const std::string& option) { - auto file = std::unique_ptr(TFile::Open(filename.c_str(), option.c_str())); + + inline void Histogrammer::BookInternal2D(const HistoDef2D& def) + { + // Resolve variable names + std::string xName = _proc.FullName(def.XvarBaseName); + std::string yName = _proc.FullName(def.YvarBaseName); + + ROOT::RVec cols = {xName, yName}; + + // Dimensions: X,Y + Splits + int n_dims = 2 + _splits.size(); + ROOT::RVec bins_vec(n_dims); + ROOT::RVec xmin(n_dims), xmax(n_dims); + + // Main X axis + bins_vec[0] = def.nbinsX; xmin[0] = def.xmin; xmax[0] = def.xmax; + + // Main Y axis + bins_vec[1] = def.nbinsY; xmin[1] = def.ymin; xmax[1] = def.ymax; + + // Add split axes + for (size_t i = 0; i < _splits.size(); ++i) { + auto& s = _splits[i]; + std::string col = _proc.FullName(s.baseCol); + cols.push_back(col); + + bins_vec[i + 2] = s.nbins; + xmin[i + 2] = s.min; + xmax[i + 2] = s.max; + } + + // Construct THnSparse + std::string fullName = _proc.GetPrefix() + def.name + _proc.GetSuffix(); + auto hist = std::make_shared( + fullName.c_str(), def.title.c_str(), + n_dims, bins_vec.data(), xmin.data(), xmax.data()); + + hist->GetAxis(0)->SetName(def.XvarBaseName.c_str()); + hist->GetAxis(1)->SetName(def.YvarBaseName.c_str()); + + for (size_t i = 0; i < _splits.size(); ++i) + hist->GetAxis(i + 2)->SetName(_splits[i].name.c_str()); + + // Masked? + if (!_maskCol.empty()) { + cols.push_back(_maskCol); + THnCombi action(hist); + _results.push_back(_rad.CurrFrame().Book(std::move(action), + utils::as_stdvector(cols))); + } else { + THnCombi action(hist); + _results.push_back(_rad.CurrFrame().Book(std::move(action), + utils::as_stdvector(cols))); + } + } + + + inline void Histogrammer::File(const std::string& filename, const std::string& option) { + auto file = std::unique_ptr(TFile::Open(filename.c_str(), option.c_str())); - if (!file || file->IsZombie()) { - std::cerr << "Histogrammer Error: Could not open file " << filename << std::endl; - return; - } - - // Trigger the Event Loop here if it hasn't run yet! - for(auto& resultPtr : _results) { - // Accessing the result triggers the RDataFrame loop - auto hist_raw = resultPtr->Clone(); - auto hn = dynamic_cast(hist_raw); + if (!file || file->IsZombie()) { + std::cerr << "Histogrammer Error: Could not open file " << filename << std::endl; + return; + } + + // Trigger the Event Loop here if it hasn't run yet! + for(auto& resultPtr : _results) { + // Accessing the result triggers the RDataFrame loop + auto hist_raw = resultPtr->Clone(); + auto hn = dynamic_cast(hist_raw); - if (hn) { - // Create a directory for this variable (e.g. "rec_Mass_miss/") - TDirectory* subdir = file->mkdir(hn->GetName()); - if (!subdir) subdir = file->GetDirectory(hn->GetName()); - subdir->cd(); + if (hn) { + // Create a directory for this variable (e.g. "rec_Mass_miss/") + TDirectory* subdir = file->mkdir(hn->GetName()); + if (!subdir) subdir = file->GetDirectory(hn->GetName()); + subdir->cd(); - UnfoldAndWrite(hn, subdir); - } - delete hn; // Cleanup clone - } - file->Close(); + UnfoldAndWrite(hn, subdir); + } + delete hn; // Cleanup clone + } + file->Close(); } - inline void Histogrammer::UnfoldAndWrite(THnSparseD* hn, TDirectory* dir, std::string current_suffix, int axis_depth) { - // Base Case: We have constrained all split axes. Project the Main Variable (Axis 0). - if (axis_depth > (int)_splits.size()) { - dir->cd(); - - // Project Axis 0 (The Variable) - // "E" option ensures errors are computed correctly - TH1D* h1 = hn->Projection(0, "E"); - - std::string baseName = hn->GetName(); - std::string full_name = baseName + current_suffix; - - h1->SetName(full_name.c_str()); - h1->SetTitle((std::string(hn->GetTitle()) + current_suffix).c_str()); - - h1->Write(); - delete h1; - return; + inline void Histogrammer::UnfoldAndWrite(THnSparseD* hn, TDirectory* dir, + std::string current_suffix, + int axis_depth) + { + int ndim = hn->GetNdimensions(); + int nsplits = _splits.size(); + bool is2D = (ndim == nsplits + 2); + + // Base case: all split axes fixed ? now project + if (axis_depth > nsplits) { + dir->cd(); + + TH1* proj = nullptr; + + if (is2D) { + // Project 2D (axes 0 and 1) + proj = hn->Projection(0, 1, "E"); + } else { + // Project 1D (axis 0 only) + proj = hn->Projection(0, "E"); } - // Recursive Step: Iterate over all bins in the current Split Axis - auto axis = hn->GetAxis(axis_depth); - int nbins = axis->GetNbins(); - - // Save range to restore later (crucial for recursion) - int saved_min = axis->GetFirst(); - int saved_max = axis->GetLast(); + std::string full_name = std::string(hn->GetName()) + current_suffix; + proj->SetName(full_name.c_str()); + proj->SetTitle((std::string(hn->GetTitle()) + current_suffix).c_str()); - for(int i=1; i<=nbins; ++i) { - // Restrict this axis to bin 'i' - axis->SetRange(i, i); - - // Construct suffix: e.g. "_Sector1" - std::string suffix = current_suffix + "_" + axis->GetName() + std::to_string(i); - UnfoldAndWrite(hn, dir, suffix, axis_depth + 1); - } - - // Restore Range - axis->SetRange(saved_min, saved_max); + proj->Write(); + delete proj; + return; + } + + // Recursive case: iterate over bins of current split axis + auto axis = hn->GetAxis(axis_depth); + int nbins = axis->GetNbins(); + + int saved_min = axis->GetFirst(); + int saved_max = axis->GetLast(); + + for (int i = 1; i <= nbins; ++i) { + axis->SetRange(i, i); + std::string suffix = current_suffix + "_" + + axis->GetName() + std::to_string(i); + UnfoldAndWrite(hn, dir, suffix, axis_depth + 1); + } + + axis->SetRange(saved_min, saved_max); } -} // namespace histo + + } // namespace histo } // namespace rad // * @details Stores parameters passed to Create() so they can be re-applied diff --git a/include/Indicing.h b/include/Indicing.h index 3f58d67..e57babb 100644 --- a/include/Indicing.h +++ b/include/Indicing.h @@ -101,7 +101,24 @@ namespace rad { }; } - +/** + * @brief Creates a strategy to filter candidate indices by value AND a required flag. + * * **Usage:** `FilterIndicesWithFlag(11, 1)` -> Returns list of indices where `vec[i] == 11` AND `flags[i] == 1`. + * This is useful for finding specific particles that were detected in a specific subdetector. + * * @tparam T Type of value to compare (e.g., int for PDG). + * @tparam F Type of flag to compare (e.g., int for detector flag). + * @param val The value to match (e.g., 11 for electron). + * @param flagVal The required flag value (defaults to 1). + * @return A lambda with signature: `(const RVec& vec, const RVec& flags) -> RVecI` + */ + template + inline auto FilterIndicesWithFlag(const T val, const F flagVal = 1) { + return [val, flagVal](const ROOT::RVec& vec, const ROOT::RVec& flags) -> ROOT::RVecI { + // Element-wise logical AND across both vectors. + // Nonzero returns the indices where both conditions are true. + return ROOT::VecOps::Nonzero(vec == val && flags == flagVal); + }; + } // ========================================================================= // Index Utilities (Helper Functions) // ========================================================================= diff --git a/include/KinematicsProcElectro.h b/include/KinematicsProcElectro.h index 3cd20a8..c3a6bc0 100644 --- a/include/KinematicsProcElectro.h +++ b/include/KinematicsProcElectro.h @@ -37,6 +37,21 @@ namespace rad { /** @brief Calculates Q2 (Negative Squared 4-Momentum Transfer). */ void Q2(); + /** @brief Calculates xbj. */ + void xbj(); + + /** @brief Calculates y. */ + void y(); + + /** @brief Calculates nu. */ + void nu(); + + /** @brief Calculates tau. */ + void tau(); + + /** @brief Calculates tauprime. */ + void tauprime(); + /** @brief Calculates Cos(Theta) in the Center of Mass frame. */ void CosThetaCM(); @@ -78,7 +93,27 @@ namespace rad { inline void KinematicsProcElectro::Q2() { RegisterCalc("Q2", rad::physics::ElS_Q2); } - + + inline void KinematicsProcElectro::xbj() { + RegisterCalc("xbj", rad::physics::ElS_xbj); + } + + inline void KinematicsProcElectro::y() { + RegisterCalc("y", rad::physics::ElS_y); + } + + inline void KinematicsProcElectro::nu() { + RegisterCalc("nu", rad::physics::ElS_nu); + } + + inline void KinematicsProcElectro::tau() { + RegisterCalc("tau", rad::physics::ElS_tau); + } + + inline void KinematicsProcElectro::tauprime() { + RegisterCalc("tauprime", rad::physics::ElS_tauprime); + } + inline void KinematicsProcElectro::CosThetaCM() { RegisterCalc("CosThetaCM", rad::physics::ElS_CosThetaCM); } diff --git a/include/KinematicsProcessor.h b/include/KinematicsProcessor.h index dde3b8e..f2ba425 100644 --- a/include/KinematicsProcessor.h +++ b/include/KinematicsProcessor.h @@ -140,6 +140,7 @@ namespace rad { /** @brief Construct a full column name: prefix + base + suffix. */ std::string FullName(const std::string& baseName) const; + std::string CheckedFullName(const std::string& baseName) const; /** * @brief Get list of all variables defined via RegisterCalc/Mass/Pt etc. * @details This is used by `AnalysisManager::Snapshot` to auto-detect columns. @@ -175,13 +176,22 @@ namespace rad { void DefineKernel(const std::string& name, Lambda&& func); void DefineTruthFlag(); - // ================================================================================= + + /** + * @brief Passes raw detector arrays through the combinatorial engine to the flat output tree. + * @param pName The target particle (e.g., "ele"). + * @param rawArray The raw data array (e.g., "rec_charge"). + * @param compSuffix The name to append to the particle for the output (e.g., "_charge"). + */ + void PassThrough(const std::string& pName, const std::string& rawArray, const std::string& compSuffix); + // ================================================================================= // Physics Shortcuts // ================================================================================= void Mass(const std::string& name, const ParticleNames_t& particles_pos, const ParticleNames_t particles_neg={}); void Mass2(const std::string& name, const ParticleNames_t& particles_pos, const ParticleNames_t particles_neg={}); void Pt(const std::string& name, const ParticleNames_t& particles_pos, const ParticleNames_t particles_neg={}); + void Energy(const std::string& name, const ParticleNames_t& particles_pos, const ParticleNames_t particles_neg={}); void ParticleTheta(const ParticleNames_t& particles); void ParticlePhi(const ParticleNames_t& particles); @@ -225,7 +235,15 @@ namespace rad { }; ROOT::RVec _groupOverrides; - void ApplyGroupOverrides(); + void ApplyGroupOverrides(); + + struct PassThroughDef { + std::string pName; + std::string rawArray; + std::string compSuffix; + }; + ROOT::RVec _passThroughs; // Stores deferred passthrough requests + }; // ================================================================================= @@ -258,6 +276,18 @@ namespace rad { for(auto& calc : _calculations) { calc.Define(this); } + + // Execute deferred PassThroughs AFTER the clear! + for(const auto& pt : _passThroughs) { + std::string outNameBase = pt.pName + pt.compSuffix; + std::string colName = FullName(outNameBase); + std::string candCol = CheckedFullName(pt.pName); + + // Use FastTake to avoid the ROOT condition vector crash! + _reaction->Define(colName, "ROOT::VecOps::Take(" + pt.rawArray + ", " + candCol + ")"); + + _registered_vars.push_back(outNameBase); + } } inline void KinematicsProcessor::ApplyGroupOverrides() { @@ -313,7 +343,20 @@ namespace rad { inline std::string KinematicsProcessor::GetPrefix() const { return _prefix; } inline std::string KinematicsProcessor::FullName(const std::string& baseName) const { - return _prefix + baseName + _suffix; + return _prefix + baseName + _suffix; + } + inline std::string KinematicsProcessor::CheckedFullName(const std::string& baseName) const { + auto fullName = _prefix + baseName + _suffix; + if(_reaction->ColumnExists(fullName)==false){ + fullName = _prefix + baseName; + if(_reaction->ColumnExists(fullName)==false){ + fullName = baseName; + } + if(_reaction->ColumnExists(fullName)==false){ + throw std::runtime_error("KinematicsProcessor::FullName, Column '" + fullName + "' does not exist with any known perfic or suffix."); + } + } + return fullName; } // --- Core Operator --- @@ -331,7 +374,6 @@ namespace rad { const auto Ncombis = indices[0].size(); CombiOutputVec_t result(Ncombis, RVec(Ncomponents, RVecResultType(Nparticles))); - ROOT::RVecD temp_px(Nparticles, consts::InvalidEntry()); ROOT::RVecD temp_py(Nparticles, consts::InvalidEntry()); ROOT::RVecD temp_pz(Nparticles, consts::InvalidEntry()); @@ -341,7 +383,7 @@ namespace rad { AuxCacheI cache_pre_i(aux_pre_i.size(), ROOT::RVecI(Nparticles)); AuxCacheD cache_post_d(aux_post_d.size(), ROOT::RVecD(Nparticles)); AuxCacheI cache_post_i(aux_post_i.size(), ROOT::RVecI(Nparticles)); - + //cout<< "KinematicsProcessor::operator() "<> components = { - // {"_px", 0}, {"_py", 1}, {"_pz", 2}, {"_m", 3} - // }; - - // // Correctly use the Suffixed result name (e.g. rec_Components_loose) - // std::string resultColName = _prefix + consts::KineComponents() + _suffix; - - // for (const auto& pName : particle_names) { - // size_t idx = _creator.GetReactionIndex(pName); - // // Correctly construct suffixed output name (e.g. rec_ele_px_loose) - // std::string full_pName = FullName(pName); - - // for (const auto& comp : components) { - // std::string suffix = comp.first; - // auto compIdx = comp.second; - - // // IMPORTANT: Register these so AnalysisManager::CollectStreamColumns - // // knows they exist and adds them to the Snapshot list. - // // We only register the Base Name + Suffix (e.g., "ele_px") - // // The Manager adds the prefix later. - // _registered_vars.push_back(pName+suffix); - - // // Direct Component Copy - // _reaction->Define(full_pName + suffix, - // [idx,compIdx](const CombiOutputVec_t& res) { - // ROOT::RVec out(res.size()); - // for(size_t i=0; i>) - // // res[i][compIdx] = List of Particles (RVec) - // // res[i][compIdx][idx] = Value (double) - // out[i] = res[i][compIdx][idx]; - // } - // return out; - // }, - // {resultColName} - // ); - // } - // } - // } - + // --- Definitions --- /** @@ -488,8 +489,12 @@ namespace rad { // // and save it to the tree. // _registered_vars.push_back(baseName); // } - // } - + // } + + inline void KinematicsProcessor::PassThrough(const std::string& pName, const std::string& rawArray, const std::string& compSuffix) { + _passThroughs.push_back({pName, rawArray, compSuffix}); + } + inline void KinematicsProcessor::Define(const std::string& name, const std::string& func) { std::string colName = FullName(name); _reaction->Define(colName, util::createFunctionCallStringFromVec("rad::util::ApplyKinematics", @@ -521,6 +526,9 @@ namespace rad { inline void KinematicsProcessor::Pt(const std::string& name, const ParticleNames_t& particles_pos, const ParticleNames_t particles_neg) { RegisterCalc(name, rad::FourVectorPtCalc, {particles_pos, particles_neg}); } + inline void KinematicsProcessor::Energy(const std::string& name, const ParticleNames_t& particles_pos, const ParticleNames_t particles_neg) { + RegisterCalc(name, rad::FourVectorECalc, {particles_pos, particles_neg}); + } inline void KinematicsProcessor::ParticleTheta(const ParticleNames_t& particles) { //here we actually perform a loop over combies //so we need to call a dunction which returns @@ -541,12 +549,12 @@ namespace rad { RegisterCalc(p+"_pmag", rad::ThreeVectorMag, {{p}}); } } - inline void KinematicsProcessor::ParticleEta(const ParticleNames_t& particles) { + inline void KinematicsProcessor::ParticleEta(const ParticleNames_t& particles) { for(const auto& p: particles){ - RegisterCalc(p+"_eta", rad::ThreeVectorPhi, {{p}}); + RegisterCalc(p+"_eta", rad::ThreeVectorEta, {{p}}); } - } - + } + inline void KinematicsProcessor::PrintReactionMap() const { std::cout << "\n=== KinematicsProcessor [" << _prefix << "] " << _suffix << " Reaction Map ===" << std::endl; std::cout << std::left << std::setw(20) << "Particle Name" << "Index" << std::endl; diff --git a/include/ParticleCreator.h b/include/ParticleCreator.h index fc52a0a..b5c020f 100644 --- a/include/ParticleCreator.h +++ b/include/ParticleCreator.h @@ -313,9 +313,9 @@ namespace rad { inline int ParticleCreator::GetIndexSafe(const std::string& name) const { auto it = _nameIndex.find(name); if (it == _nameIndex.end()) { - std::cerr << "\n[ParticleCreator FATAL] Missing Particle Index: " << name << "\n"; - std::cerr << " Processor: " << _prefix << " (Suffix: '" << _suffix << "')\n"; - throw std::runtime_error("Particle '" + name + "' missing in map. Check inputs."); + std::cerr << "\n[ParticleCreator FATAL] Missing Particle Index: " << name << "\n"; + std::cerr << " Processor: " << _prefix << " (Suffix: '" << _suffix << "')\n"; + throw std::runtime_error("Particle '" + name + "' missing in map. Check inputs."); } return it->second; } diff --git a/include/ParticleModifier.h b/include/ParticleModifier.h index 4d668f5..6d28c6c 100644 --- a/include/ParticleModifier.h +++ b/include/ParticleModifier.h @@ -32,6 +32,16 @@ namespace rad { */ void ScaleMomentum(const std::string& name, double scale); + /** + * @brief Smear the momentum of a particle by a fixed width. + */ + void SmearMomentum(const std::string& name, double width); + + /** + * @brief Smears the energy of a particle by a fixed width, and assumes sqrt response. + */ + void SmearCalMomentum(const std::string& name, double width); + /** * @brief Forces the mass of a particle to a fixed value (modifies energy). */ @@ -43,7 +53,12 @@ namespace rad { */ void SetMomentumFrom(const std::string& name, const std::string& colName); - // --- Initialization & Execution --- + /** + * @brief Add user defined modifier inheriting from ModifierBase + */ + void AddModifier(const std::string& name,std::shared_ptr mod) ; + + // --- Initialization & Execution --- /** * @brief Resolves particle names to fixed array indices. @@ -96,19 +111,33 @@ namespace rad { // cout<< "ParticleModifier::ParticleModifier copy "< mod) { + _modi_configs.push_back({name, std::move(mod)}); + } + inline void ParticleModifier::ScaleMomentum(const std::string& name, double scale) { - auto mod = std::make_unique(scale); - _modi_configs.push_back({name, std::move(mod)}); + auto mod = std::make_shared(scale); + AddModifier(name,mod); + } + + inline void ParticleModifier::SmearMomentum(const std::string& name, double width) { + auto mod = std::make_shared(width); + AddModifier(name,mod); + } + + inline void ParticleModifier::SmearCalMomentum(const std::string& name, double width) { + auto mod = std::make_shared(width); + AddModifier(name,mod); } inline void ParticleModifier::FixMass(const std::string& name, double mass) { - auto mod = std::make_unique(mass); - _modi_configs.push_back({name, std::move(mod)}); + auto mod = std::make_shared(mass); + AddModifier(name,mod); } inline void ParticleModifier::SetMomentumFrom(const std::string& name, const std::string& colName) { Indice_t row = RegisterAuxDouble(colName); - auto mod = std::make_unique(row); + auto mod = std::make_shared(row); _modi_configs.push_back({name, std::move(mod)}); } diff --git a/include/ParticleModifierMethods.h b/include/ParticleModifierMethods.h index 772b6ba..8247a1d 100644 --- a/include/ParticleModifierMethods.h +++ b/include/ParticleModifierMethods.h @@ -1,5 +1,6 @@ #pragma once #include "CommonDefines.h" +#include "Random.h" #include #include #include @@ -52,6 +53,42 @@ namespace rad { double _scale; }; + /** + * @brief Scales the 3-momentum (px, py, pz) by a random fluctuation. + * scale factor is centered on 1 with given width + */ + class ModSmearMomentum : public ModifierBase { + public: + explicit ModSmearMomentum(double width); + ModSmearMomentum(const Indices_t&, const Indices_t&, double width); + + std::unique_ptr Clone() const override; + + void operator()(ROOT::RVecD& px, ROOT::RVecD& py, + ROOT::RVecD& pz, ROOT::RVecD& m, + const AuxCacheD&, const AuxCacheI&) const override; + private: + double _width; + }; + /** + * @brief Scales the 3-momentum (px, py, pz) by a random fluctuation. + * assuming calorimeter energy measurement + * scale factor is centered on 1 with given width + */ + class ModCalSmearMomentum : public ModifierBase { + public: + explicit ModCalSmearMomentum(double width); + ModCalSmearMomentum(const Indices_t&, const Indices_t&, double width); + + std::unique_ptr Clone() const override; + + void operator()(ROOT::RVecD& px, ROOT::RVecD& py, + ROOT::RVecD& pz, ROOT::RVecD& m, + const AuxCacheD&, const AuxCacheI&) const override; + private: + double _width; + }; + /** * @brief Constraints the mass of a particle to a fixed value. * Useful when the detector resolution suggests a specific particle (e.g. Proton) @@ -110,6 +147,51 @@ namespace rad { py[_target_idx] *= _scale; pz[_target_idx] *= _scale; } + // --- ModSmearMomentum --- + inline ModSmearMomentum::ModSmearMomentum(double width) : _width(width) {} + inline ModSmearMomentum::ModSmearMomentum(const Indices_t&, const Indices_t&, double width) : _width(width) {} + + inline std::unique_ptr ModSmearMomentum::Clone() const { + return std::make_unique(_width); + } + + inline void ModSmearMomentum::operator()(ROOT::RVecD& px, ROOT::RVecD& py, + ROOT::RVecD& pz, ROOT::RVecD& m, + const AuxCacheD&, const AuxCacheI&) const + { + if (_target_idx < 0 || _target_idx >= static_cast(px.size())) return; + auto scale = random::Generator().Gaus(1,_width); + px[_target_idx] *= scale; + py[_target_idx] *= scale; + pz[_target_idx] *= scale; + } + // --- ModCalSmearMomentum --- + inline ModCalSmearMomentum::ModCalSmearMomentum(double width) : _width(width) {} + inline ModCalSmearMomentum::ModCalSmearMomentum(const Indices_t&, const Indices_t&, double width) : _width(width) {} + + inline std::unique_ptr ModCalSmearMomentum::Clone() const { + return std::make_unique(_width); + } + + inline void ModCalSmearMomentum::operator()(ROOT::RVecD& px, ROOT::RVecD& py, + ROOT::RVecD& pz, ROOT::RVecD& m, + const AuxCacheD&, const AuxCacheI&) const + { + if (_target_idx < 0 || _target_idx >= static_cast(px.size())) return; + + auto kin_energy = TMath::Sqrt(px[_target_idx]*px[_target_idx]+py[_target_idx]*py[_target_idx]+pz[_target_idx]*pz[_target_idx] +m[_target_idx]*m[_target_idx])-m[_target_idx]; + + auto scale_K = random::Generator().Gaus(1,_width/TMath::Sqrt(kin_energy)); + kin_energy *=scale_K; + + auto tot_energy = kin_energy + m[_target_idx]; + auto momentum = TMath::Sqrt(tot_energy*tot_energy - m[_target_idx]*m[_target_idx]); + auto scale = momentum/TMath::Sqrt(px[_target_idx]*px[_target_idx]+py[_target_idx]*py[_target_idx]+pz[_target_idx]*pz[_target_idx]); + + px[_target_idx] *= scale; + py[_target_idx] *= scale; + pz[_target_idx] *= scale; + } // --- ModFixMass --- inline ModFixMass::ModFixMass(double mass) : _mass(mass) {} diff --git a/include/PhysicsSelection.h b/include/PhysicsSelection.h index 0535e91..1ac6984 100644 --- a/include/PhysicsSelection.h +++ b/include/PhysicsSelection.h @@ -24,12 +24,13 @@ namespace rad { * @brief Configuration storage for deferred definition. * @details Stores cut parameters so they can be compiled later. */ - struct CutDef { + struct CutDef { enum Type { Range, Min, Max, // Standard Comparison Equal, NotEqual, // Exact Match (IDs) - AbsRange, AbsMin, AbsMax // Magnitude Checks - }; + AbsRange, AbsMin, AbsMax, // Magnitude Checks + Bool + }; std::string name; ///< Name of the cut (e.g. "MassCut") std::string varBaseName; ///< Unresolved variable name (e.g. "MassJ") double min; @@ -91,6 +92,7 @@ namespace rad { * @param val Value to match. */ void AddCutEqual(const std::string& name, const std::string& var, double val); + void AddCutBool(const std::string& name, const std::string& var); /** * @brief Define an inequality cut: var != val. * @param name Unique name for this cut. @@ -118,8 +120,13 @@ namespace rad { * @param max Maximum absolute value. */ void AddCutAbsMax(const std::string& name, const std::string& var, double max); - + + // ===================================================================== + // // ===================================================================== + void AddTruthMatchedCut(const std::string& name = "truth_match_cut"); + + // ===================================================================== // Execution // ===================================================================== @@ -169,6 +176,10 @@ namespace rad { inline void PhysicsSelection::AddCutEqual(const std::string& name, const std::string& var, double val) { _config.push_back({name, var, val, 0.0, CutDef::Equal}); } + + inline void PhysicsSelection::AddCutBool(const std::string& name, const std::string& var) { + _config.push_back({name, var, 0.0, 0.0, CutDef::Bool}); + } inline void PhysicsSelection::AddCutNotEqual(const std::string& name, const std::string& var, double val) { _config.push_back({name, var, val, 0.0, CutDef::NotEqual}); @@ -181,64 +192,67 @@ namespace rad { inline void PhysicsSelection::AddCutAbsMax(const std::string& name, const std::string& var, double max) { _config.push_back({name, var, 0.0, max, CutDef::AbsMax}); } - + // --- Compilation (The Actual Work) --- inline void PhysicsSelection::Init() { _cutNames.clear(); _finalMask = _proc.GetPrefix() + "Analysis_Mask" + _proc.GetSuffix(); - std::cout << "[PhysicsSelection] Initializing Cuts for Stream: " + std::cout << "[PhysicsSelection] Initializing Cuts for Stream: " << _proc.GetPrefix() << "..." << _proc.GetSuffix() << " (" << _config.size() << " cuts configured)" << std::endl; // 1. Define individual cut columns for(const auto& def : _config) { - std::string col = _proc.FullName(def.varBaseName); - std::string cutName = _proc.GetPrefix() + def.name + _proc.GetSuffix(); - // DIAGNOSTIC: Print what we are defining - std::cout << " -> Defining Cut: " << cutName << " on Variable: " << col << std::endl; - - double min = def.min; - double max = def.max; - - switch(def.type) { - // Standard - case CutDef::Range: - _proc.Reaction()->Define(cutName, [min, max](const RVecResultType& val){ return val > min && val < max; }, {col}); - break; - case CutDef::Min: - _proc.Reaction()->Define(cutName, [min](const RVecResultType& val){ return val > min; }, {col}); - break; - case CutDef::Max: - _proc.Reaction()->Define(cutName, [max](const RVecResultType& val){ return val < max; }, {col}); - break; + std::string col = _proc.FullName(def.varBaseName); + std::string cutName = _proc.GetPrefix() + def.name + _proc.GetSuffix(); + // DIAGNOSTIC: Print what we are defining + std::cout << " -> Defining Cut: " << cutName << " on Variable: " << col << std::endl; + + double min = def.min; + double max = def.max; + + switch(def.type) { + // Standard + case CutDef::Range: + _proc.Reaction()->Define(cutName, [min, max](const RVecResultType& val){ return val > min && val < max; }, {col}); + break; + case CutDef::Min: + _proc.Reaction()->Define(cutName, [min](const RVecResultType& val){ return val > min; }, {col}); + break; + case CutDef::Max: + _proc.Reaction()->Define(cutName, [max](const RVecResultType& val){ return val < max; }, {col}); + break; - // Equality - case CutDef::Equal: - _proc.Reaction()->Define(cutName, [min](const RVecResultType& val){ return ROOT::VecOps::abs(val - min) < 1e-9; }, {col}); - break; - case CutDef::NotEqual: - _proc.Reaction()->Define(cutName, [min](const RVecResultType& val){ return ROOT::VecOps::abs(val - min) > 1e-9; }, {col}); - break; - - // Absolute - case CutDef::AbsRange: - _proc.Reaction()->Define(cutName, [min, max](const RVecResultType& val){ const RVecResultType& a=ROOT::VecOps::abs(val); return a > min && a < max; }, {col}); - break; - case CutDef::AbsMax: - _proc.Reaction()->Define(cutName, [max](const RVecResultType& val){ return ROOT::VecOps::abs(val) < max; }, {col}); - break; - default: break; - } - _cutNames.push_back(cutName); + // Equality + case CutDef::Equal: + _proc.Reaction()->Define(cutName, [min](const RVecResultType& val){ return ROOT::VecOps::abs(val - min) < 1e-9; }, {col}); + break; + case CutDef::Bool: + _proc.Reaction()->Define(cutName, [](const Indices_t& val){ return val; }, {col}); + break; + case CutDef::NotEqual: + _proc.Reaction()->Define(cutName, [min](const RVecResultType& val){ return ROOT::VecOps::abs(val - min) > 1e-9; }, {col}); + break; + + // Absolute + case CutDef::AbsRange: + _proc.Reaction()->Define(cutName, [min, max](const RVecResultType& val){ const RVecResultType& a=ROOT::VecOps::abs(val); return a > min && a < max; }, {col}); + break; + case CutDef::AbsMax: + _proc.Reaction()->Define(cutName, [max](const RVecResultType& val){ return ROOT::VecOps::abs(val) < max; }, {col}); + break; + default: break; + } + _cutNames.push_back(cutName); } // 2. Define Master Mask (AND logic) if (_cutNames.empty()) { // CASE: NO CUTS. // We cannot use "1" because it creates a scalar boolean. - // We need a VECTOR of 1s (all true) with the same length as the combinations. + // We need a VECTOR of increasing ints with the same length as the combinations. // We use the first particle's Px column to determine the size. auto particle_names = _proc.Creator().GetParticleNames(); @@ -248,10 +262,10 @@ namespace rad { // _proc.Reaction()->Define(_finalMask, "1"); } else { std::string refCol = _proc.FullName(particle_names[0] + "_px"); - _proc.Reaction()->Define(_finalMask, [](const ROOT::RVecD& v){ - // Return vector of same size, all true (1) - return ROOT::RVecI(v.size(), 1); - }, {refCol}); + _proc.Reaction()->Define(_finalMask, [](const ROOT::RVecD& v){ + // Return vector of same size, all true (1) + return rad::util::EnumerateIndicesFrom(0, v.size()); + }, {refCol}); } } else { @@ -268,15 +282,16 @@ namespace rad { _proc.Reaction()->Define(boolMaskName, ss.str()); // B. Convert Boolean -> Indices // Nonzero({1, 0, 1}) -> {0, 2} - _proc.Reaction()->Define(_finalMask, - [](const Indices_t& bools) { - return ROOT::VecOps::Nonzero(bools); - }, - {boolMaskName} - ); + _proc.Reaction()->Define(_finalMask, + [](const Indices_t& bools) { + //cout<<"PhysicsSelection "< void Filter(Lambda&& func, const ROOT::RDF::ColumnNames_t& columns = {}, std::string name = ""); - // --- Output & Snapshots --- + /** @brief Creates an alias for an existing column. */ + void SetBranchAlias(const string_view aliasName, const string_view columnName); + + // --- Output & Snapshots --- /** @brief Immediate Snapshot (Virtual, implementation specific). */ virtual void Snapshot(const string& filename) = 0; @@ -222,7 +225,9 @@ namespace rad { /** @brief Returns the C++ type string of a column (e.g. "double", "vector"). */ string ColObjTypeString(const string& name); - + /** @brief Extracts inner type string from "RVec". Needed for template metaprogramming. */ + std::string GetElementTypeName(const std::string& colName); + protected: ROOT::RDataFrame _orig_df; RDFstep _curr_df; @@ -326,6 +331,9 @@ namespace rad { SetCurrFrame(CurrFrame().Filter(func, columns, name)); } + inline void RDFInterface::SetBranchAlias(const string_view aliasName, const string_view columnName) { + SetCurrFrame(CurrFrame().Alias(columnName,aliasName));//Alias(New, Target) + } // --- Output & Snapshots --- inline void RDFInterface::RemoveSnapshotColumns(ROOT::RVec& cols) { @@ -396,5 +404,13 @@ namespace rad { inline string RDFInterface::ColObjTypeString(const string& name){ return CurrFrame().GetColumnType(name); } - + inline std::string RDFInterface::GetElementTypeName(const std::string& colName) { + std::string fullType = ColObjTypeString(colName); + auto start = fullType.find('<'); + auto end = fullType.rfind('>'); + if(start != std::string::npos && end != std::string::npos) { + return fullType.substr(start + 1, end - start - 1); + } + return fullType; + } } // namespace rad diff --git a/include/RDFUtils.h b/include/RDFUtils.h index b03a4b7..c0adede 100644 --- a/include/RDFUtils.h +++ b/include/RDFUtils.h @@ -6,6 +6,7 @@ #include #include +#include // Needed for thread-safe counting #include // For TString operations if preferred, or std::string namespace rad { @@ -86,7 +87,7 @@ namespace rad { // ========================================================================= - void PrintDefinedColumnNames(RNode df){ + inline void PrintDefinedColumnNames(RNode df){ std::cout<<"Print Column Names : "; auto cols = df.GetDefinedColumnNames(); for(auto& col:cols){ @@ -94,7 +95,7 @@ namespace rad { } cout<<"\n"; } - void PrintAllColumnNames(RNode df){ + inline void PrintAllColumnNames(RNode df){ std::cout<<"Print Column Names : "; auto cols = df.GetColumnNames(); for(auto& col:cols){ @@ -103,5 +104,38 @@ namespace rad { cout<<"\n"; } + template + inline void PrintColumnValues(T& radf, + const ROOT::RDF::ColumnNames_t& cols, + int max_events = -1, + const std::string& name = "column_print") { + + if (cols.empty() ) return; + + // 1. Inject a static, thread-safe atomic counter into the Cling environment. + // fetch_add(1) safely increments and returns the PREVIOUS value. + std::string code = "static std::atomic print_count{0}; "; + if(max_events!=-1){ + code += "if (print_count.fetch_add(1) < " + std::to_string(max_events) + ") { "; + } + else{ + code += "{"; + } + code += " std::cout << \"--- PrintColumns Event ---\" < void ResolutionFraction(T* const rdf, const std::string& var); - + + /** + * @brief Returns mass in GeV for detector-stable particles. + * @details Evaluated entirely at compile-time. Uses absolute value + * to automatically handle anti-particles perfectly. + */ + constexpr double PdgToMass(int pdg); + + /** + * @brief returns synched vector of masses given PDG codes + */ + ROOT::RVecD AssignMasses( const ROOT::RVecI &pid); + } // namespace util } // namespace rad @@ -61,7 +73,7 @@ namespace rad { * @param topo Vector of pairs: `{Required_Count, PID}`. */ template - void CountTopo(T* rdf, const std::string& type, const ROOT::RVec>& topo) { + inline void CountTopo(T* rdf, const std::string& type, const ROOT::RVec>& topo) { std::string col_name = type + "has_topo_debug"; @@ -103,7 +115,7 @@ namespace rad { * * **Lazy Execution:** Queues definitions and runs the event loop only ONCE at the end. */ template - void CheckCandidateStats(T* reaction, const std::string& prefix) { + inline void CheckCandidateStats(T* reaction, const std::string& prefix) { auto names = reaction->ParticleNames(); if(names.empty()) return; @@ -173,7 +185,7 @@ namespace rad { std::cout << "==========================================================\n" << std::endl; } template - void CountParticles(T* rdf, const std::string& type){ + inline void CountParticles(T* rdf, const std::string& type){ rdf->Define(type+"Ngamma", Form("rad::util::Count(%spid,22)", type.data()) ); rdf->Define(type+"Npip", Form("rad::util::Count(%spid,211)", type.data()) ); rdf->Define(type+"Npim", Form("rad::util::Count(%spid,-211)", type.data()) ); @@ -186,14 +198,50 @@ namespace rad { } template - void Resolution(T* const rdf, const std::string& var){ + inline void Resolution(T* const rdf, const std::string& var){ rdf->Define(string("res_")+var, Form("%s-%s", (Truth()+var).data(), (Rec()+var).data() )); } template - void ResolutionFraction(T* const rdf, const std::string& var){ + inline void ResolutionFraction(T* const rdf, const std::string& var){ rdf->Define(string("res_")+var, Form("(%s-%s)/%s", (Truth()+var).data(), (Rec()+var).data(), (Truth()+var).data() )); } + + + inline constexpr ResultType_t PdgToMass(int pdg) { + // Safe constexpr absolute value + int abs_pdg = (pdg < 0) ? -pdg : pdg; + + switch (abs_pdg) { + case 22: return 0.0; // Photon + case 11: return 0.00051099895; // Electron / Positron + case 13: return 0.10565837; // Muon / Anti-muon + case 211: return 0.13957039; // Charged Pion + case 321: return 0.493677; // Charged Kaon + case 2212: return 0.93827208; // Proton / Anti-proton + case 2112: return 0.93956541; // Neutron / Anti-neutron + + // Light Nuclei + case 45: return 1.875612; // Deuteron (CLAS12 old convention) + case 1000010020: return 1.875612; // Deuteron (Strict PDG convention) + case 1000020030: return 2.80839; // Helium-3 + case 1000020040: return 3.72738; // Alpha + + default: return 0.0; // Fallback for unmatched/undefined PID + } + } + + + /////////////////////////////////////////////////////// + inline RVecResultType AssignMasses( const ROOT::RVecI &pid){ + auto n = pid.size(); + RVecResultType masses(n); + for(size_t i=0;i(all_args); + if (!mask.empty()) { - for (auto idx : mask) { + for (const auto idx : mask) { fill_buffers(slot, idx, all_args, std::make_index_sequence{}); _trees[slot]->Fill(); _slotCounts[slot]++; diff --git a/include/gammaN_2_Spin0Spin0SpinHalf.h b/include/gammaN_2_Spin0Spin0SpinHalf.h index baa5a45..276ea14 100644 --- a/include/gammaN_2_Spin0Spin0SpinHalf.h +++ b/include/gammaN_2_Spin0Spin0SpinHalf.h @@ -1,130 +1,225 @@ #pragma once -//#include "Beams.h" -#include "BasicKinematics.h" -#include "ConfigReaction.h" +#include "ReactionKinematics.h" // FourVector, boost, PxPyPzMVector, etc. +#include "BasicKinematics.h" // If InitialFourVector / helpers live here +#include "ConfigReaction.h" // RVecIndexMap, names:: indices +#include "CommonDefines.h" // Ensures RVecResultType, ResultType_t, etc. +#include "TMath.h" +namespace rad { + namespace gn2s0s0s12 { + //namespace physics { -namespace rad{ - namespace gn2s0s0s12{ - - struct Angles_t{ - double CosTheta=0.; - double Phi=0.; + // --- Structures --- + + /** + * @brief Result structure for generic decay angles. + */ + struct DecayAngles_t { + double cosTheta = 0.; ///< Cosine of the polar decay angle. + double theta = 0.; ///< Polar decay angle. + double phi = 0.; ///< Azimuthal decay angle. }; - ///\brief functions to compute standard reaction kinematics + + // --- Center-of-Mass (CM) Kinematics --- + + /** + * @brief Calculates the initial CM four-momentum for photoproduction: + * P_CM = P_initial(bottom beam) + q (photon). + * + * @details + * Assumes `PhotoFourVector(...)` is defined by the reaction configuration + * (e.g. PhotoIonReaction.h or ElectroIonReaction.h). + * + * @param react The fixed reaction index map. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return PxPyPzMVector The CM four-momentum vector. + */ + inline PxPyPzMVector PhotoCMVector(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + // Assumes OrderBeams() [0] = ion, [1] = electron + // Assuming ion is at pos 0 in the beam group + auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); + auto phot = PhotoFourVector(react, px, py, pz, m); + + return ion + phot; + } /** - * calculate CM kinematics from beam - */ - template - PxPyPzMVector PhotoCMVector(const config::RVecIndexMap& react, const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - - //Note PhotoFourVector must be defined in the reaction config file - //e.g. in PhotoIonReaction.h or ElectroIonReaction.h - return beams::InitialFourVector(react[names::InitialBotIdx()][0],px,py,pz,m) + - PhotoFourVector(react,px,py,pz,m); - - } - - /** - * calculate CM decay angles - */ - template - Angles_t PhotoCMDecay(const config::RVecIndexMap& react,const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - - auto cm = PhotoCMVector(react,px,py,pz,m); - auto cmBoost = cm.BoostToCM(); - auto mes = FourVector(react[names::MesonsIdx()],px,py,pz,m); - PxPyPzMVector cm_mes=boost(mes,cmBoost); - - Angles_t result; - result.CosTheta=(TMath::Cos(cm_mes.Theta())); - result.Phi=cm_mes.Phi(); + * @brief Calculates decay angles in the overall CM frame of the initial state. + * + * @param react The fixed reaction index map. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return DecayAngles_t {cosTheta, phi} of the meson system in the CM frame. + */ + inline DecayAngles_t PhotoCMDecay(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + const auto cm = PhotoCMVector(react, px, py, pz, m); + const auto cmBoost = cm.BoostToCM(); + + auto mes = FourVector(react[consts::OrderMesons()], px, py, pz, m); // Assumes single particle or combined meson vector + const PxPyPzMVector cmMes = boost(mes, cmBoost); + + DecayAngles_t result; + result.cosTheta = TMath::Cos(cmMes.Theta()); + result.phi = cmMes.Phi(); + result.theta = cmMes.Theta(); return result; - } + } + + // --- Helicity Frame Decay Angles --- + /** - * calculate Helicity decay angles - * z-axis along -baryon in meson rest frame - */ - template - Angles_t PhotoHelicityDecay(const config::RVecIndexMap& react,const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - - auto baryon = reactkine::BaryonFourVector(react,px,py,pz,m); - auto meson = reactkine::MesonFourVector(react,px,py,pz,m); - - auto decBoost = meson.BoostToCM(); - //vectors in rest/decay frame of meson - auto decBar=boost(baryon,decBoost); - auto decGamma=boost(PhotoFourVector(react,px,py,pz,m),decBoost); + * @brief Calculates helicity decay angles in the meson rest frame. + * + * @details + * Helicity frame convention: + * - z-axis along **-baryon** direction in the meson rest frame. + * - y-axis along (baryon × photon). + * - x-axis completes the right-handed set: x = y × z. + * + * @param react The fixed reaction index map. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return DecayAngles_t {cosTheta, phi} for the first meson child (index 0). + */ + inline DecayAngles_t PhotoHelicityDecay(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + const auto baryon = FourVector(react[consts::OrderBaryons()], px, py, pz, m); + //const auto meson = MesonFourVector(react, px, py, pz, m); + const auto meson = FourVector(react[consts::OrderMesons()], px, py, pz, m); + const auto photon = PhotoFourVector(react, px, py, pz, m); - XYZVector zV=-decBar.Vect().Unit(); - XYZVector yV=decBar.Vect().Cross(decGamma.Vect()).Unit(); - XYZVector xV=yV.Cross(zV).Unit(); - - //four vector of first [0] decay product - auto child1 = FourVector(react[names::MesonsIdx()][0],px,py,pz,m); - auto decChild1=boost(child1,decBoost); - - //calculate decay angles - XYZVector angles(decChild1.Vect().Dot(xV),decChild1.Vect().Dot(yV),decChild1.Vect().Dot(zV)); - //store in angles struct - Angles_t result; - result.CosTheta=TMath::Cos(angles.Theta()); - result.Phi=angles.Phi(); + const auto decBoost = meson.BoostToCM(); + + // Four-vectors in the meson rest frame + const auto decBar = boost(baryon, decBoost); + const auto decGamma = boost(photon, decBoost); + + // Helicity axes + const XYZVector zV = -decBar.Vect().Unit(); + const XYZVector yV = decBar.Vect().Cross(decGamma.Vect()).Unit(); + const XYZVector xV = yV.Cross(zV).Unit(); + + // First decay product of the meson + const auto child1 = FourVector(react[consts::OrderMesons()][0], px, py, pz, m); + const auto decChild1 = boost(child1, decBoost); + + // Projections and angles + const XYZVector proj(decChild1.Vect().Dot(xV), + decChild1.Vect().Dot(yV), + decChild1.Vect().Dot(zV)); + + DecayAngles_t result; + result.cosTheta = TMath::Cos(proj.Theta()); + result.phi = proj.Phi(); + result.theta = proj.Theta(); return result; } - template - double CosThetaHel(const config::RVecIndexMap& react,const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - auto angles = PhotoHelicityDecay(react,px,py,pz,m); - return angles.CosTheta; + + /** + * @brief Convenience: returns cos(theta) in the helicity frame. + */ + inline ResultType_t CosThetaHel(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + return PhotoHelicityDecay(react, px, py, pz, m).cosTheta; } - template - double PhiHel(const config::RVecIndexMap& react,const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - auto angles = PhotoHelicityDecay(react,px,py,pz,m); - return angles.Phi; + + /** + * @brief Convenience: returns theta in the helicity frame. + */ + inline ResultType_t ThetaHel(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + return PhotoHelicityDecay(react, px, py, pz, m).theta; } - /** - * calculate GJ decay angles - * z-axis along gamma direction in meson rest frame - */ - template - Angles_t PhotoGJDecay(const config::RVecIndexMap& react,const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - - auto baryon = reactkine::BaryonFourVector(react,px,py,pz,m); - auto meson = reactkine::MesonFourVector(react,px,py,pz,m); - - auto decBoost = meson.BoostToCM(); - //vectors in rest/decay frame of meson - auto decBar=boost(baryon,decBoost); - auto decGamma=boost(PhotoFourVector(react,px,py,pz,m),decBoost); + + /** + * @brief Convenience: returns phi in the helicity frame. + */ + inline ResultType_t PhiHel(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + return PhotoHelicityDecay(react, px, py, pz, m).phi; + } + + // --- Gottfried–Jackson (GJ) Frame Decay Angles --- + + /** + * @brief Calculates Gottfried–Jackson (GJ) decay angles in the meson rest frame. + * + * @details + * GJ frame convention: + * - z-axis along **photon** direction in the meson rest frame. + * - y-axis along (baryon × photon). + * - x-axis completes the right-handed set: x = y × z. + * + * @param react The fixed reaction index map. + * @param px, py, pz, m The consolidated momentum component vectors. + * @return DecayAngles_t {cosTheta, phi} for the first meson child (index 0). + */ + inline DecayAngles_t PhotoGJDecay(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + const auto baryon = FourVector(react[consts::OrderBaryons()], px, py, pz, m); + const auto meson = FourVector(react[consts::OrderMesons()], px, py, pz, m); - XYZVector zV=decGamma.Vect().Unit(); - XYZVector yV=decBar.Vect().Cross(decGamma.Vect()).Unit(); - XYZVector xV=yV.Cross(zV).Unit(); - - //four vector of first [0] decay product - auto child1 = FourVector(react[names::MesonsIdx()][0],px,py,pz,m); - auto decChild1=boost(child1,decBoost); - - //calculate decay angles - XYZVector angles(decChild1.Vect().Dot(xV),decChild1.Vect().Dot(yV),decChild1.Vect().Dot(zV)); - //store in angles struct - Angles_t result; - result.CosTheta=TMath::Cos(angles.Theta()); - result.Phi=angles.Phi(); + const auto decBoost = meson.BoostToCM(); + + // Four-vectors in the meson rest frame + const auto decBar = boost(baryon, decBoost); + const auto decGamma = boost(PhotoFourVector(react, px, py, pz, m), decBoost); + + // GJ axes + const XYZVector zV = decGamma.Vect().Unit(); + const XYZVector yV = decBar.Vect().Cross(decGamma.Vect()).Unit(); + const XYZVector xV = yV.Cross(zV).Unit(); + + // First decay product of the meson + const auto child1 = FourVector(react[consts::OrderMesons()][0], px, py, pz, m); + const auto decChild1 = boost(child1, decBoost); + + // Projections and angles + const XYZVector proj(decChild1.Vect().Dot(xV), + decChild1.Vect().Dot(yV), + decChild1.Vect().Dot(zV)); + + DecayAngles_t result; + result.cosTheta = TMath::Cos(proj.Theta()); + result.phi = proj.Phi(); return result; } - template - double CosThetaGJ(const config::RVecIndexMap& react,const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - auto angles = PhotoGJDecay(react,px,py,pz,m); - return angles.CosTheta; + + /** + * @brief Convenience: returns cos(theta) in the GJ frame. + */ + inline ResultType_t CosThetaGJ(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + return PhotoGJDecay(react, px, py, pz, m).cosTheta; } - template - double PhiGJ(const config::RVecIndexMap& react,const RVec &px, const RVec &py, const RVec &pz, const RVec &m){ - auto angles = PhotoGJDecay(react,px,py,pz,m); - return angles.Phi; + + /** + * @brief Convenience: returns phi in the GJ frame. + */ + inline ResultType_t PhiGJ(const RVecIndexMap& react, + const RVecResultType& px, const RVecResultType& py, + const RVecResultType& pz, const RVecResultType& m) + { + return PhotoGJDecay(react, px, py, pz, m).phi; } - } -} + + } // namespace gn2s0s0s12 +} // namespace rad