From 31a95916d2e00b3b5056105ffec9bdaa2e39e45d Mon Sep 17 00:00:00 2001 From: garypenman Date: Sat, 21 Feb 2026 10:00:04 +0000 Subject: [PATCH 1/8] Added extra DIS variables to electron scattering: xbj, nu, tau, tau prime, circular polarisation of virt photon. Added shortcuts to KinematicsProcElectro --- include/ElectronScatterKinematics.h | 98 ++++++++++++++++++++++++++++- include/KinematicsProcElectro.h | 37 ++++++++++- 2 files changed, 132 insertions(+), 3 deletions(-) diff --git a/include/ElectronScatterKinematics.h b/include/ElectronScatterKinematics.h index 7eb6bd6..adc9dcb 100644 --- a/include/ElectronScatterKinematics.h +++ b/include/ElectronScatterKinematics.h @@ -43,7 +43,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 +184,7 @@ namespace rad { return {theta, energy, Q2}; } - + /** * @brief Calculates the polarization parameter (epsilon) for the virtual photon. * @param react The fixed ReactionMap. @@ -128,6 +208,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) --- diff --git a/include/KinematicsProcElectro.h b/include/KinematicsProcElectro.h index c41a952..89c0c78 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(); @@ -74,7 +89,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); } From 8c7b13a9cdb3815c65a4c602c8a0bbad6307aa0f Mon Sep 17 00:00:00 2001 From: garypenman Date: Sat, 21 Feb 2026 10:09:27 +0000 Subject: [PATCH 2/8] updated photo helicity decay header with inline functions and doxygen --- include/gammaN_2_Spin0Spin0SpinHalf.h | 313 +++++++++++++++++--------- 1 file changed, 203 insertions(+), 110 deletions(-) diff --git a/include/gammaN_2_Spin0Spin0SpinHalf.h b/include/gammaN_2_Spin0Spin0SpinHalf.h index baa5a45..131bad6 100644 --- a/include/gammaN_2_Spin0Spin0SpinHalf.h +++ b/include/gammaN_2_Spin0Spin0SpinHalf.h @@ -1,130 +1,223 @@ #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 --- /** - * 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 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; + } + + /** + * @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(); 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(); 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; + } + + /** + * @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; } - 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 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; } - /** - * 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); + + // --- 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 From d9c51492a1398abf31d4454ef3adb3a8e3120927 Mon Sep 17 00:00:00 2001 From: garypenman Date: Sat, 21 Feb 2026 10:10:16 +0000 Subject: [PATCH 3/8] Added TCS HepMC example --- examples/ProcessHepMCTCSCombi.C | 122 ++++++++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) create mode 100644 examples/ProcessHepMCTCSCombi.C diff --git a/examples/ProcessHepMCTCSCombi.C b/examples/ProcessHepMCTCSCombi.C new file mode 100644 index 0000000..f26cfd6 --- /dev/null +++ b/examples/ProcessHepMCTCSCombi.C @@ -0,0 +1,122 @@ +#include "CommonDefines.h" +#include "HepMCElectro.h" +#include "KinematicsProcElectro.h" +#include "KineCalculation.h" +#include "Indicing.h" +#include "ElectronScatterKinematics.h" +#include "gammaN_2_Spin0Spin0SpinHalf.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("MissMass2", + {consts::BeamIon(), consts::BeamEle()}, + {"gprime", consts::ScatEle(), "pprime"}); + + + // Custom Calc (t_prime) + 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); + + // ================================================================================= + // 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"); +} From f45d55cab436fa7a0be12f2ccd4b726902e6512c Mon Sep 17 00:00:00 2001 From: garypenman Date: Tue, 24 Feb 2026 12:37:51 +0000 Subject: [PATCH 4/8] indicing tests TCS --- examples/ProcessHepMCTCSCombi.C | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/ProcessHepMCTCSCombi.C b/examples/ProcessHepMCTCSCombi.C index f26cfd6..47bed4d 100644 --- a/examples/ProcessHepMCTCSCombi.C +++ b/examples/ProcessHepMCTCSCombi.C @@ -5,6 +5,7 @@ #include "Indicing.h" #include "ElectronScatterKinematics.h" #include "gammaN_2_Spin0Spin0SpinHalf.h" +#include "TCSKinematics.h" #include /** @@ -83,17 +84,19 @@ void ProcessHepMCTCSCombi(){ kine.PhiCM(); kine.Mass("GMass", {"gprime"}); + kine.Mass2("Qp2", {"gprime"}); kine.Mass2("MissMass2", {consts::BeamIon(), consts::BeamEle()}, {"gprime", consts::ScatEle(), "pprime"}); - - // Custom Calc (t_prime) 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) // ================================================================================= From b5985a7d4a6ae493a4794dfaa1f7a1ccab70ecb1 Mon Sep 17 00:00:00 2001 From: garypenman Date: Mon, 2 Mar 2026 13:16:20 +0000 Subject: [PATCH 5/8] added overlad to analysis manager constructor to read in vector of strings for file lists --- include/AnalysisManager.h | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/include/AnalysisManager.h b/include/AnalysisManager.h index d98df18..7c496a3 100644 --- a/include/AnalysisManager.h +++ b/include/AnalysisManager.h @@ -59,7 +59,19 @@ 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 std::vector& files); + + /** @brief Sets and creates the output directory. */ void SetOutputDir(const std::string& dir); @@ -211,8 +223,17 @@ 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 std::vector& files) + : _reaction(treeName, files), _name(name) {} template inline void AnalysisManager::SetOutputDir(const std::string& dir) { From 0e5eb64d72735333be771f911c43584d00a85d62 Mon Sep 17 00:00:00 2001 From: garypenman Date: Mon, 2 Mar 2026 13:20:20 +0000 Subject: [PATCH 6/8] added theta to DecayAngles_t struct, for brufit and other purposes --- include/gammaN_2_Spin0Spin0SpinHalf.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/gammaN_2_Spin0Spin0SpinHalf.h b/include/gammaN_2_Spin0Spin0SpinHalf.h index 131bad6..276ea14 100644 --- a/include/gammaN_2_Spin0Spin0SpinHalf.h +++ b/include/gammaN_2_Spin0Spin0SpinHalf.h @@ -68,6 +68,7 @@ namespace rad { DecayAngles_t result; result.cosTheta = TMath::Cos(cmMes.Theta()); result.phi = cmMes.Phi(); + result.theta = cmMes.Theta(); return result; } @@ -118,6 +119,7 @@ namespace rad { DecayAngles_t result; result.cosTheta = TMath::Cos(proj.Theta()); result.phi = proj.Phi(); + result.theta = proj.Theta(); return result; } From 1970b2c5fd05019ab67db4f2beb7f359e6705b54 Mon Sep 17 00:00:00 2001 From: garypenman Date: Mon, 2 Mar 2026 13:21:10 +0000 Subject: [PATCH 7/8] default behaviour for useBeamsFromMC should be false --- include/ElectroIonReaction.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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, From b7679ad3089c4bd3ff7aa03c16043c42ea2fcea7 Mon Sep 17 00:00:00 2001 From: garypenman Date: Tue, 3 Mar 2026 09:39:21 +0000 Subject: [PATCH 8/8] fixed erronous pickup of file from old branch in last commit --- include/ElectroIonReaction.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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,