diff --git a/math/physics/CMakeLists.txt b/math/physics/CMakeLists.txt index 2b55d54742819..b33afcc316aa7 100644 --- a/math/physics/CMakeLists.txt +++ b/math/physics/CMakeLists.txt @@ -10,6 +10,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Physics HEADERS + TFANG.h TFeldmanCousins.h TGenPhaseSpace.h TLorentzRotation.h @@ -20,7 +21,8 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Physics TRotation.h TVector2.h TVector3.h - SOURCES + SOURCES + src/TFANG.cxx src/TFeldmanCousins.cxx src/TGenPhaseSpace.cxx src/TLorentzRotation.cxx @@ -32,9 +34,13 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Physics src/TVector2.cxx src/TVector3.cxx DEPENDENCIES + Core Matrix MathCore GenVector DICTIONARY_OPTIONS -writeEmptyRootPCM ) +if(testing) + add_subdirectory(test) +endif() diff --git a/math/physics/inc/TFANG.h b/math/physics/inc/TFANG.h new file mode 100644 index 0000000000000..6c5e561e8d777 --- /dev/null +++ b/math/physics/inc/TFANG.h @@ -0,0 +1,488 @@ +// @(#)root/physics:$Id$ +// Author: Arik Kreisel, Itay Horin + +/************************************************************************* + * Copyright (C) 1995-2025, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +//////////////////////////////////////////////////////////////////////////////// +/// \file TFANG.h +/// +/// \brief TFANG - Focused Angular N-body event Generator +/// \authors Arik Kreisel, Itay Horin +/// +/// TFANG is a Monte Carlo tool for efficient event generation in restricted +/// (or full) Lorentz-Invariant Phase Space (LIPS). Unlike conventional approaches +/// that always sample the full 4pi solid angle, TFANG can also directly generate +/// events in which selected final-state particles are constrained to fixed +/// directions or finite angular regions in the laboratory frame. +/// +/// This file contains: +/// - namespace FANG: Core generator functions and utilities +/// - class TFANG: User-friendly interface class +/// +/// Reference: Horin, I., Kreisel, A. & Alon, O. Focused angular N -body event generator (FANG). +/// J. High Energ. Phys. 2025, 137 (2025). +/// https://doi.org/10.1007/JHEP12(2025)137 +/// https://arxiv.org/abs/2509.11105 +//////////////////////////////////////////////////////////////////////////////// + +#ifndef ROOT_TFANG +#define ROOT_TFANG + +#include "Rtypes.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" + +#include + +class TRandom3; + +//////////////////////////////////////////////////////////////////////////////// +// FANG namespace - Core generator functions +//////////////////////////////////////////////////////////////////////////////// +namespace FANG { + +//////////////////////////////////////////////////////////////////////////////// +/// Mathematical constants +//////////////////////////////////////////////////////////////////////////////// +constexpr Double_t kPi = 3.14159265358979323846; +constexpr Double_t kTwoPi = 2.0 * kPi; +constexpr Double_t kFourPi = 4.0 * kPi; + +//////////////////////////////////////////////////////////////////////////////// +/// Physics constants +//////////////////////////////////////////////////////////////////////////////// +constexpr Double_t kDipoleMassSq = 0.71; ///< GeV^2, for form factor +constexpr Double_t kProtonMagneticMoment = 2.793; ///< mu_p in nuclear magnetons + +//////////////////////////////////////////////////////////////////////////////// +/// Numerical tolerances +//////////////////////////////////////////////////////////////////////////////// +constexpr Double_t kPositionTolerance = 1E-6; +constexpr Double_t kMomentumTolerance = 1E-12; + +//////////////////////////////////////////////////////////////////////////////// +/// Generation mode constants +/// +/// These correspond to the Shape parameter values: +/// - POINT (2): Fixed direction for all events +/// - RING (<0): Fixed polar angle, uniform azimuthal +/// - CIRCLE (0): Uniform within a cone +/// - STRIP (0 0.0 && shape <= 1.0; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Check if shape parameter indicates ring generation mode +/// \param[in] shape Shape parameter value +/// \return kTRUE if ring mode +//////////////////////////////////////////////////////////////////////////////// +inline Bool_t IsRing(Double_t shape) +{ + return shape < 0.0; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \struct Node_t +/// \brief Binary tree node for tracking multiple kinematic solutions +/// +/// When a constrained particle can reach a detector via two different +/// momenta, both solutions are stored in a binary tree structure. +//////////////////////////////////////////////////////////////////////////////// +struct Node_t { + ROOT::Math::PxPyPzMVector fPV; ///< Virtual system 4-momentum + ROOT::Math::PxPyPzMVector fPDet; ///< Detected particle 4-momentum + Double_t fWeight; ///< Weight (Jacobian contribution) + Node_t *fLeft; ///< Left child (second solution) + Node_t *fRight; ///< Right child (first solution) + Node_t *fParent; ///< Parent node + + /// Constructor + Node_t(const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight, Node_t *parent); + + /// Prevent copying to avoid ownership issues + Node_t(const Node_t &) = delete; + Node_t &operator=(const Node_t &) = delete; +}; + +//////////////////////////////////////////////////////////////////////////////// +// Tree Management Functions +//////////////////////////////////////////////////////////////////////////////// + +/// \brief Recursively delete a tree and free all memory +/// \param[in] node Root of the tree to delete +void DeleteTree(Node_t *node); + +/// \brief Create the first (root) node of the tree +Node_t *CreateFirst(Node_t *node, + const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight); + +/// \brief Create a right child node +Node_t *CreateRight(Node_t *node, Node_t *tmp, + const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight); + +/// \brief Create a left child node +Node_t *CreateLeft(Node_t *node, Node_t *tmp, + const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight); + +/// \brief Collect all root-to-leaf paths for 4-momenta +void CollectPaths(Int_t nBody, Node_t *node, + std::vector &path, + std::vector> &paths); + +/// \brief Collect all root-to-leaf paths for weights +void CollectPathsWeights(Int_t nBody, Node_t *node, + std::vector &path, + std::vector> &paths); + +//////////////////////////////////////////////////////////////////////////////// +// Utility Functions +//////////////////////////////////////////////////////////////////////////////// + +/// \brief Phase space kinematic function F(x,y) = sqrt((1-x-y)^2 - 4xy) +/// \param[in] x First mass ratio squared (m1^2/M^2) +/// \param[in] y Second mass ratio squared (m2^2/M^2) +/// \return Kinematic function value +Double_t CalcKMFactor(Double_t x, Double_t y); + +//////////////////////////////////////////////////////////////////////////////// +// Core Physics Functions +//////////////////////////////////////////////////////////////////////////////// + +/// \brief Generate isotropic two-body decay +/// +/// Performs a two-body decay isotropically in the rest frame of S, +/// then boosts results back to the lab frame. +/// +/// \param[in] S Total 4-momentum of decaying system +/// \param[in] m1 Mass of first decay product +/// \param[in] m2 Mass of second decay product +/// \param[out] p1 4-momentum of first decay product (lab frame) +/// \param[out] p2 4-momentum of second decay product (lab frame) +/// \param[in] rng Pointer to TRandom3 random number generator (thread-safe) +void TwoBody(const ROOT::Math::PxPyPzMVector &S, + Double_t m1, Double_t m2, + ROOT::Math::PxPyPzMVector &p1, + ROOT::Math::PxPyPzMVector &p2, + TRandom3 *rng); + +/// \brief Calculate 4-momentum for particle constrained to a lab-frame direction +/// +/// Given a two-body system S1 decaying to masses m1 and m2, with m1 constrained +/// to travel in direction vDet, calculate the possible 4-momenta. +/// +/// \param[in] S1 Total 4-momentum of the decaying system +/// \param[in] m1 Mass of constrained particle +/// \param[in] m2 Mass of other particle +/// \param[in] vDet Unit vector specifying lab-frame direction for m1 +/// \param[out] solutions Number of physical solutions (0, 1, or 2) +/// \param[out] jackPDF Array of Jacobian * PDF values for each solution +/// \param[out] pDet Array of 4-momenta for constrained particle +/// \param[out] pD2 Array of 4-momenta for other particle +/// \return kTRUE if at least one physical solution exists +Bool_t TGenPointSpace(const ROOT::Math::PxPyPzMVector &S1, + Double_t m1, Double_t m2, + ROOT::Math::XYZVector vDet, + Int_t &solutions, + Double_t *jackPDF, + ROOT::Math::PxPyPzMVector *pDet, + ROOT::Math::PxPyPzMVector *pD2); + +/// \brief Generate random direction vector within specified solid angle +/// +/// \param[in] Omega Solid angle size [steradians] +/// \param[in] Ratio Shape parameter determining generation mode +/// \param[in] Vcenter Central direction vector +/// \param[out] vPoint Generated direction vector +/// \param[in] rng Pointer to TRandom3 random number generator (thread-safe) +void TGenVec(Double_t Omega, Double_t Ratio, + ROOT::Math::XYZVector Vcenter, + ROOT::Math::XYZVector &vPoint, + TRandom3 *rng); + +//////////////////////////////////////////////////////////////////////////////// +// Main Generator Function +//////////////////////////////////////////////////////////////////////////////// + +/// \brief Generate phase-space events with angular constraints +/// +/// Main FANG generator function. Generates n-body phase space events +/// where selected particles are constrained to specified detector directions. +/// +/// \param[in] nBody Number of outgoing particles +/// \param[in] S Total 4-momentum of the system +/// \param[in] masses Array of outgoing particle masses [GeV], length nBody +/// \param[in] Om Array of solid angles for constrained detectors [sr] +/// \param[in] Ratio Array of shape parameters for each detector: +/// - = 2: Point generation (fixed direction) +/// - = 0: Circle generation (uniform in cone) +/// - 0 < Ratio[] <= 1: Strip generation (rectangular region) +/// Dphi = Ratio[] * TwoPi; +/// Dcos = Omega / Dphi; +/// - < 0: Ring generation (fixed theta, uniform phi) +/// \param[in] V3Det Vector of direction vectors for constrained detectors +/// \param[out] VecVecP Output: vector of 4-momenta vectors for each solution +/// \param[out] vecWi Output: weight for each solution +/// \param[in] rng Pointer to TRandom3 random number generator (thread-safe) +/// \return 1 on success, 0 if no physical solution exists +Int_t GenFANG(Int_t nBody, + const ROOT::Math::PxPyPzMVector &S, + const Double_t *masses, + const Double_t *Om, + const Double_t *Ratio, + std::vector V3Det, + std::vector> &VecVecP, + std::vector &vecWi, + TRandom3 *rng); + +} // namespace FANG + +//////////////////////////////////////////////////////////////////////////////// +/// \class TFANG +/// \brief User interface class for the FANG phase space generator +/// +/// TFANG wraps the FANG namespace functions to provide an object-oriented +/// interface for phase space event generation with optional angular constraints. +/// +/// Example usage (unconstrained): +/// ~~~{.cpp} +/// // Create generator for 3-body decay +/// TFANG gen; +/// Double_t masses[3] = {0.139, 0.139, 0.139}; // pion masses +/// ROOT::Math::PxPyPzMVector S(0, 0, 0, 1.0); // parent at rest +/// gen.SetDecay(S, 3, masses); +/// +/// // Generate single event +/// gen.Generate(); +/// auto p0 = gen.GetDecay(0); // 4-momentum of first particle +/// +/// // Calculate phase space integral +/// Double_t ps, err; +/// gen.GetPhaseSpace(1000000, ps, err); +/// ~~~ +/// +/// Example usage (constrained): +/// ~~~{.cpp} +/// TFANG gen; +/// Double_t masses[3] = {0.139, 0.139, 0.139}; +/// ROOT::Math::PxPyPzMVector S(0, 0, 0, 1.0); +/// gen.SetDecay(S, 3, masses); +/// +/// // Add detector constraint +/// ROOT::Math::XYZVector detDir(0, 0, 1); +/// Double_t omega = 0.01; // solid angle +/// Double_t ratio = 0.0; // circle mode +/// gen.AddConstraint(detDir, omega, ratio); +/// +/// // Generate and loop over solutions +/// if (gen.Generate() > 0) { +/// for (Int_t i = 0; i < gen.GetNSolutions(); ++i) { +/// Double_t w = gen.GetWeight(i); +/// auto p0 = gen.GetDecay(i, 0); +/// } +/// } +/// ~~~ +//////////////////////////////////////////////////////////////////////////////// +class TFANG { + +private: + Int_t fNBody; ///< Number of outgoing particles + ROOT::Math::PxPyPzMVector fS; ///< Total 4-momentum of system + std::vector fMasses; ///< Particle masses + std::vector fOmega; ///< Solid angles for constraints + std::vector fRatio; ///< Shape parameters for constraints + std::vector fV3Det; ///< Direction vectors for constraints + + std::vector> fVecVecP; ///< Output 4-momenta for all solutions + std::vector fVecWi; ///< Weights for each solution + + TRandom3 *fRng; ///< Random number generator + Bool_t fOwnRng; ///< Whether we own the RNG + +public: + /// Default constructor + TFANG(); + + /// Constructor with external random number generator + /// \param[in] rng Pointer to external TRandom3 (ownership not transferred) + explicit TFANG(TRandom3 *rng); + + /// Destructor + ~TFANG(); + + /// Prevent copying + TFANG(const TFANG &) = delete; + TFANG &operator=(const TFANG &) = delete; + + //////////////////////////////////////////////////////////////////////////// + // Configuration methods + //////////////////////////////////////////////////////////////////////////// + + /// \brief Set decay configuration + /// \param[in] S Total 4-momentum of the decaying system + /// \param[in] nBody Number of outgoing particles + /// \param[in] masses Array of particle masses, length nBody + /// \return kTRUE if configuration is valid + Bool_t SetDecay(const ROOT::Math::PxPyPzMVector &S, Int_t nBody, const Double_t *masses); + + /// \brief Add angular constraint for a particle + /// \param[in] direction Direction vector for the constraint + /// \param[in] omega Solid angle in steradians + /// \param[in] ratio Shape parameter (2=point, 0=circle, 0(fVecVecP.size()); } + + /// \brief Get weight for a specific solution (unconstrained mode) + /// \return Weight for the single solution + Double_t GetWeight() const; + + /// \brief Get weight for a specific solution (constrained mode) + /// \param[in] iSolution Solution index + /// \return Weight for the specified solution + Double_t GetWeight(Int_t iSolution) const; + + /// \brief Get 4-momentum of particle (unconstrained mode, single solution) + /// \param[in] iParticle Particle index (0 to nBody-1) + /// \return 4-momentum of the particle + ROOT::Math::PxPyPzMVector GetDecay(Int_t iParticle) const; + + /// \brief Get 4-momentum of particle for specific solution + /// \param[in] iSolution Solution index + /// \param[in] iParticle Particle index (0 to nBody-1) + /// \return 4-momentum of the particle + ROOT::Math::PxPyPzMVector GetDecay(Int_t iSolution, Int_t iParticle) const; + + /// \brief Get all 4-momenta for a solution + /// \param[in] iSolution Solution index + /// \return Vector of 4-momenta + const std::vector &GetDecays(Int_t iSolution) const; + + //////////////////////////////////////////////////////////////////////////// + // Phase space integration methods + //////////////////////////////////////////////////////////////////////////// + + /// \brief Calculate full (unconstrained) phase space integral with uncertainty + /// + /// Runs Monte Carlo integration to estimate the full phase space volume. + /// This method ignores any constraints and calculates the Lorentz-invariant + /// phase space for the decay. + /// + /// \param[in] nEvents Number of events for Monte Carlo integration (default 1E6) + /// \param[out] phaseSpace Estimated phase space value + /// \param[out] error Statistical uncertainty on phase space + /// \return kTRUE if calculation succeeded + Bool_t GetPhaseSpace(Long64_t nEvents, Double_t &phaseSpace, Double_t &error) const; + + /// \brief Calculate full phase space with default 1E6 events + /// \param[out] phaseSpace Estimated phase space value + /// \param[out] error Statistical uncertainty on phase space + /// \return kTRUE if calculation succeeded + Bool_t GetPhaseSpace(Double_t &phaseSpace, Double_t &error) const; + + /// \brief Calculate partial (constrained) phase space integral with uncertainty + /// + /// Runs Monte Carlo integration to estimate the partial phase space volume + /// when angular constraints are applied. The result is the phase space + /// restricted to the detector solid angles, multiplied by the product of + /// all omega values (total solid angle factor). + /// + /// Requires constraints to be set via AddConstraint() before calling. + /// + /// \param[in] nEvents Number of events for Monte Carlo integration (default 1E6) + /// \param[out] phaseSpace Estimated partial phase space value + /// \param[out] error Statistical uncertainty on phase space + /// \return kTRUE if calculation succeeded, kFALSE if no constraints set + Bool_t GetPartialPhaseSpace(Long64_t nEvents, Double_t &phaseSpace, Double_t &error) const; + + /// \brief Calculate partial phase space with default 1E6 events + /// \param[out] phaseSpace Estimated partial phase space value + /// \param[out] error Statistical uncertainty on phase space + /// \return kTRUE if calculation succeeded, kFALSE if no constraints set + Bool_t GetPartialPhaseSpace(Double_t &phaseSpace, Double_t &error) const; + + //////////////////////////////////////////////////////////////////////////// + // Utility methods + //////////////////////////////////////////////////////////////////////////// + + /// \brief Check if generator is in constrained mode + /// \return kTRUE if constraints are set + Bool_t IsConstrained() const { return !fV3Det.empty(); } + + /// \brief Get number of constraints + /// \return Number of detector constraints + Int_t GetNConstraints() const { return static_cast(fV3Det.size()); } + + /// \brief Get number of particles in decay + /// \return Number of outgoing particles + Int_t GetNBody() const { return fNBody; } +}; + +#endif // ROOT_TFANG diff --git a/math/physics/src/TFANG.cxx b/math/physics/src/TFANG.cxx new file mode 100644 index 0000000000000..915749c2e4e26 --- /dev/null +++ b/math/physics/src/TFANG.cxx @@ -0,0 +1,1165 @@ +// @(#)root/physics:$Id$ +// Author: Arik Kreisel, Itay Horin + +/************************************************************************* + * Copyright (C) 1995-2025, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +//////////////////////////////////////////////////////////////////////////////// +/// \file TFANG.cxx +/// +/// \brief Implementation of TFANG (Focused Angular N-body event Generator) +/// \authors Arik Kreisel, Itay Horin +/// +/// TFANG is a Monte Carlo tool for efficient event generation in restricted +/// (or full) Lorentz-Invariant Phase Space (LIPS). Unlike conventional approaches +/// that always sample the full 4pi solid angle, TFANG can also directly generate +/// events in which selected final-state particles are constrained to fixed +/// directions or finite angular regions in the laboratory frame. +/// +/// Reference: Horin, I., Kreisel, A. & Alon, O. Focused angular N -body event generator (FANG). +/// J. High Energ. Phys. 2025, 137 (2025). +/// https://doi.org/10.1007/JHEP12(2025)137 +/// https://arxiv.org/abs/2509.11105 +//////////////////////////////////////////////////////////////////////////////// + + +#include "TFANG.h" +#include "TRandom3.h" +#include "TMath.h" +#include "TError.h" +#include "Math/GenVector/Rotation3D.h" +#include "Math/GenVector/Boost.h" +#include "Math/GenVector/Polar3D.h" + +#include +#include + +//////////////////////////////////////////////////////////////////////////////// +// FANG namespace - implementation +//////////////////////////////////////////////////////////////////////////////// +namespace FANG { + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Node_t constructor implementation +//////////////////////////////////////////////////////////////////////////////// +Node_t::Node_t(const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight, Node_t *parent) + : fPV(p2) + , fPDet(p1) + , fWeight(weight) + , fLeft(nullptr) + , fRight(nullptr) + , fParent(parent) +{ +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Recursively delete a tree and free all memory +//////////////////////////////////////////////////////////////////////////////// +void DeleteTree(Node_t *node) +{ + if (node == nullptr) + return; + DeleteTree(node->fLeft); + DeleteTree(node->fRight); + delete node; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Create the first (root) node of the tree +//////////////////////////////////////////////////////////////////////////////// +Node_t *CreateFirst(Node_t *node, + const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight) +{ + if (node == nullptr) { + return new Node_t(p1, p2, weight, nullptr); + } + return node; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Create a right child node +//////////////////////////////////////////////////////////////////////////////// +Node_t *CreateRight(Node_t *node, Node_t *tmp, + const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight) +{ + if (node == nullptr) { + return new Node_t(p1, p2, weight, tmp); + } + node->fRight = CreateRight(node->fRight, node, p1, p2, weight); + return node; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Create a left child node +//////////////////////////////////////////////////////////////////////////////// +Node_t *CreateLeft(Node_t *node, Node_t *tmp, + const ROOT::Math::PxPyPzMVector &p1, + const ROOT::Math::PxPyPzMVector &p2, + Double_t weight) +{ + if (node == nullptr) { + return new Node_t(p1, p2, weight, tmp); + } + node->fLeft = CreateLeft(node->fLeft, node, p1, p2, weight); + return node; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Collect all root-to-leaf paths for 4-momenta +//////////////////////////////////////////////////////////////////////////////// +void CollectPaths(Int_t nBody, Node_t *node, + std::vector &path, + std::vector> &paths) +{ + if (node == nullptr) + return; + + path.push_back(node->fPDet); + + if (node->fLeft == nullptr && node->fRight == nullptr && + path.size() == static_cast(nBody + 1)) { + paths.push_back(path); + } else { + CollectPaths(nBody, node->fLeft, path, paths); + CollectPaths(nBody, node->fRight, path, paths); + } + + path.pop_back(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Collect all root-to-leaf paths for weights +//////////////////////////////////////////////////////////////////////////////// +void CollectPathsWeights(Int_t nBody, Node_t *node, + std::vector &path, + std::vector> &paths) +{ + if (node == nullptr) + return; + + path.push_back(node->fWeight); + + if (node->fLeft == nullptr && node->fRight == nullptr && + path.size() == static_cast(nBody + 1)) { + paths.push_back(path); + } else { + CollectPathsWeights(nBody, node->fLeft, path, paths); + CollectPathsWeights(nBody, node->fRight, path, paths); + } + + path.pop_back(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Phase space kinematic function F(x,y) = sqrt((1-x-y)^2 - 4xy) +//////////////////////////////////////////////////////////////////////////////// +Double_t CalcKMFactor(Double_t x, Double_t y) +{ + Double_t arg = (1.0 - x - y) * (1.0 - x - y) - 4.0 * x * y; + if (arg < 0) { + ::Warning("TFANG::CalcKMFactor", "Received negative sqrt argument: %g", arg); + return 0.0; + } + return std::sqrt(arg); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Generate isotropic two-body decay +//////////////////////////////////////////////////////////////////////////////// +void TwoBody(const ROOT::Math::PxPyPzMVector &S, + Double_t m1, Double_t m2, + ROOT::Math::PxPyPzMVector &p1, + ROOT::Math::PxPyPzMVector &p2, + TRandom3 *rng) +{ + Double_t cst = rng->Uniform(-1.0, 1.0); + Double_t snt = std::sqrt(1.0 - cst * cst); + Double_t phi = rng->Uniform(0.0, kTwoPi); + + Double_t E1 = (S.M2() - m2 * m2 + m1 * m1) / (2.0 * S.M()); + + if ((E1 * E1 - m1 * m1) < 0) { + ::Error("TFANG::TwoBody", "E1^2 - m1^2 < 0, E1=%g, m1=%g", E1, m1); + return; + } + + Double_t sp = std::sqrt(E1 * E1 - m1 * m1); + + ROOT::Math::PxPyPzMVector p1CM(sp * snt * std::cos(phi), + sp * snt * std::sin(phi), + sp * cst, m1); + ROOT::Math::PxPyPzMVector p2CM(-sp * snt * std::cos(phi), + -sp * snt * std::sin(phi), + -sp * cst, m2); + + ROOT::Math::XYZVector betaVS = S.BoostToCM(); + ROOT::Math::Boost bstCM; + bstCM.SetComponents(betaVS); + ROOT::Math::Boost bstLAB = bstCM.Inverse(); + + p1 = bstLAB(p1CM); + p1.SetM(m1); + p2 = bstLAB(p2CM); + p2.SetM(m2); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate 4-momentum for particle constrained to a lab-frame direction +//////////////////////////////////////////////////////////////////////////////// +Bool_t TGenPointSpace(const ROOT::Math::PxPyPzMVector &S1, + Double_t m1, Double_t m2, + ROOT::Math::XYZVector vDet, + Int_t &solutions, + Double_t *jackPDF, + ROOT::Math::PxPyPzMVector *pDet, + ROOT::Math::PxPyPzMVector *pD2) +{ + ROOT::Math::XYZVector VSu(S1.Px() / S1.P(), + S1.Py() / S1.P(), + S1.Pz() / S1.P()); + VSu = VSu.Unit(); + + ROOT::Math::XYZVector betaVS(-S1.Beta() * VSu.X(), + -S1.Beta() * VSu.Y(), + -S1.Beta() * VSu.Z()); + ROOT::Math::Boost bstCM; + bstCM.SetComponents(betaVS); + ROOT::Math::Boost bstLAB = bstCM.Inverse(); + + vDet = vDet.Unit(); + LongDouble_t cosLAB = VSu.Dot(vDet); + LongDouble_t sinLAB = std::sqrt(1.0 - cosLAB * cosLAB); + + LongDouble_t mCM = S1.M(); + LongDouble_t ECM = S1.E(); + LongDouble_t pCM = S1.P(); + LongDouble_t gamma1 = S1.Gamma(); + + LongDouble_t CME3 = (mCM * mCM - m2 * m2 + m1 * m1) / (2.0 * mCM); + + if ((CME3 * CME3 - m1 * m1) < 0) { + ::Error("TFANG::TGenPointSpace", "CMp3 negative sqrt"); + ::Error("TFANG::TGenPointSpace", "S1.M()=%g S1.P()=%g", (Double_t)S1.M(), (Double_t)S1.P()); + ::Error("TFANG::TGenPointSpace", "m1=%g m2=%g", m1, m2); + solutions = 0; + return kFALSE; + } + + LongDouble_t CMp3 = std::sqrt(CME3 * CME3 - m1 * m1); + + LongDouble_t aa = pCM * pCM * cosLAB * cosLAB - ECM * ECM; + LongDouble_t bb = 2.0 * pCM * cosLAB * CME3 * mCM; + LongDouble_t cc = CME3 * mCM * CME3 * mCM - m1 * m1 * ECM * ECM; + + LongDouble_t discriminant = bb * bb - 4.0 * aa * cc; + + jackPDF[0] = 0.0; + jackPDF[1] = 0.0; + solutions = 1; + + if (discriminant < 0) { + solutions = 0; + return kFALSE; + } + + LongDouble_t p3LAB[2]; + LongDouble_t sqrtDisc = std::sqrt(discriminant); + + p3LAB[0] = (-bb + sqrtDisc) / (2.0 * aa); + + if (p3LAB[0] <= 0) { + p3LAB[0] = (-bb - sqrtDisc) / (2.0 * aa); + if (p3LAB[0] <= 0) { + solutions = 0; + return kFALSE; + } + } else { + p3LAB[1] = (-bb - sqrtDisc) / (2.0 * aa); + if (p3LAB[1] > 0) { + solutions = 2; + } + } + + LongDouble_t pdfCM = 1.0 / kFourPi; + + for (Int_t l = 0; l < solutions; l++) { + pDet[l].SetCoordinates(p3LAB[l] * vDet.X(), + p3LAB[l] * vDet.Y(), + p3LAB[l] * vDet.Z(), m1); + + ROOT::Math::PxPyPzMVector p3CM = bstCM(pDet[l]); + p3CM.SetM(m1); + + ROOT::Math::PxPyPzMVector p4CM(-p3CM.Px(), -p3CM.Py(), -p3CM.Pz(), m2); + pD2[l] = bstLAB(p4CM); + pD2[l].SetM(m2); + + if (std::abs(pD2[l].M() - m2) > kMomentumTolerance) { + ::Warning("TFANG::TGenPointSpace", "Mass mismatch: %g != %g", + pD2[l].M(), m2); + } + + LongDouble_t cosCM = p3CM.Vect().Dot(VSu) / CMp3; + LongDouble_t qqq = pCM * p3CM.E() / (ECM * p3CM.P()); + + LongDouble_t Jack; + + if (std::abs(cosLAB) > 0.99 && std::abs(cosCM) > 0.99) { + Jack = gamma1 * gamma1 * cosCM * cosCM * + (1.0 + qqq * cosCM) * (1.0 + qqq * cosCM) / + (cosLAB * cosLAB); + } else { + Jack = ((1.0 - cosCM) * (1.0 + cosCM)) / + ((1.0 - cosLAB) * (1.0 + cosLAB)) * + std::sqrt(((1.0 - cosCM) * (1.0 + cosCM)) / + ((1.0 - cosLAB) * (1.0 + cosLAB))) / + gamma1 / (1.0 + qqq * cosCM); + } + + jackPDF[l] = std::abs(pdfCM * Jack); + } + + return kTRUE; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Generate random direction vector within specified solid angle +//////////////////////////////////////////////////////////////////////////////// +void TGenVec(Double_t Omega, Double_t Ratio, + ROOT::Math::XYZVector Vcenter, + ROOT::Math::XYZVector &vPoint, + TRandom3 *rng) +{ + ROOT::Math::XYZVector newZ, newX, newY, Vz; + ROOT::Math::Polar3DVector Vgen; + Double_t cst, phi, Dphi, Dcos, phi0, cst0; + + if (Omega > kFourPi || Omega < 0) { + Omega = kFourPi; + ::Warning("TFANG::TGenVec", "Omega out of range, set to 4pi"); + } + + if (Ratio > 1.0) { + Ratio = 0.0; + ::Warning("TFANG::TGenVec", "Ratio out of range, set to 0"); + } + + if (IsCircle(Ratio)) { + cst = rng->Uniform(1.0 - Omega / kTwoPi, 1.0); + phi = rng->Uniform(0.0, kTwoPi); + + if (std::abs(Vcenter.X()) < kPositionTolerance && + std::abs(Vcenter.Y()) < kPositionTolerance) { + if (Vcenter.Z() > 0) { + Vgen.SetCoordinates(1.0, std::acos(cst), phi); + } else { + Vgen.SetCoordinates(1.0, std::acos(-cst), phi); + } + vPoint = Vgen; + } else { + Vz.SetXYZ(0, 0, 1); + newZ = Vcenter.Unit(); + newY = newZ.Cross(Vz).Unit(); + newX = newY.Cross(newZ).Unit(); + ROOT::Math::Rotation3D m(newX.X(), newY.X(), newZ.X(), + newX.Y(), newY.Y(), newZ.Y(), + newX.Z(), newY.Z(), newZ.Z()); + Vgen.SetCoordinates(1.0, std::acos(cst), phi); + vPoint = m * Vgen; + } + } else if (IsStrip(Ratio)) { + Dphi = Ratio * kTwoPi; + Dcos = Omega / Dphi; + phi0 = Vcenter.Phi(); + cst0 = std::cos(Vcenter.Theta()); + + if (cst0 > (1.0 - Dcos / 2.0)) { + cst0 = 1.0 - Dcos / 2.0; + ::Warning("TFANG::TGenVec", "Center moved to agree with Omega (near +1)"); + } + if (cst0 < (-1.0 + Dcos / 2.0)) { + cst0 = -1.0 + Dcos / 2.0; + ::Warning("TFANG::TGenVec", "Center moved to agree with Omega (near -1)"); + } + + cst = rng->Uniform(cst0 - Dcos / 2.0, cst0 + Dcos / 2.0); + phi = rng->Uniform(phi0 - Dphi / 2.0, phi0 + Dphi / 2.0); + Vgen.SetCoordinates(1.0, std::acos(cst), phi); + vPoint = Vgen; + } else if (IsRing(Ratio)) { + cst = 1.0 - Omega / kTwoPi; + phi = rng->Uniform(0.0, kTwoPi); + + if (std::abs(Vcenter.X()) < kPositionTolerance && + std::abs(Vcenter.Y()) < kPositionTolerance) { + if (Vcenter.Z() > 0) { + Vgen.SetCoordinates(1.0, std::acos(cst), phi); + } else { + Vgen.SetCoordinates(1.0, std::acos(-cst), phi); + } + vPoint = Vgen; + } else { + Vz.SetXYZ(0, 0, 1); + newZ = Vcenter.Unit(); + newY = newZ.Cross(Vz).Unit(); + newX = newY.Cross(newZ).Unit(); + ROOT::Math::Rotation3D m(newX.X(), newY.X(), newZ.X(), + newX.Y(), newY.Y(), newZ.Y(), + newX.Z(), newY.Z(), newZ.Z()); + Vgen.SetCoordinates(1.0, std::acos(cst), phi); + vPoint = m * Vgen; + } + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Main FANG generator function +//////////////////////////////////////////////////////////////////////////////// +Int_t GenFANG(Int_t nBody, + const ROOT::Math::PxPyPzMVector &S, + const Double_t *masses, + const Double_t *Om, + const Double_t *Ratio, + std::vector V3Det, + std::vector> &VecVecP, + std::vector &vecWi, + TRandom3 *rng) +{ + Int_t nDet = static_cast(V3Det.size()); + Double_t mS = S.M(); + Double_t wh = 1.0; + Double_t mB, mA, mall, whPS; + + mall = 0.0; + for (Int_t l = 0; l < nBody; l++) { + mall += masses[l]; + } + + if (mall >= mS) { + ::Error("TFANG::GenFANG", "Sum of decay masses (%g) >= parent mass (%g)", + mall, mS); + return 0; + } + + std::vector> pathsP; + std::vector> pathsJ; + std::vector vecP; + std::vector vecJ; + std::vector branch; + + vecP.clear(); + vecJ.clear(); + vecWi.clear(); + pathsJ.clear(); + + Bool_t Hit; + ROOT::Math::XYZVector V3; + ROOT::Math::PxPyPzMVector p1; + ROOT::Math::PxPyPzMVector p2; + ROOT::Math::PxPyPzMVector pV; + Int_t solutions; + Double_t jackPDF[2]; + ROOT::Math::PxPyPzMVector pDet[2]; + ROOT::Math::PxPyPzMVector pD2[2]; + + std::vector mV(nBody - 2); + std::vector rrr(nBody - 2); + + //========================================================================== + // Two-body decay case + //========================================================================== + if (nBody == 2) { + whPS = CalcKMFactor(masses[0] * masses[0] / S.M2(), + masses[1] * masses[1] / S.M2()) * kPi / 2.0; + + if (nDet == 1) { + if (IsPoint(Ratio[0])) { + V3 = V3Det[0].Unit(); + } else { + TGenVec(Om[0], Ratio[0], V3Det[0].Unit(), V3, rng); + } + + Hit = TGenPointSpace(S, masses[0], masses[1], V3, solutions, + jackPDF, pDet, pD2); + if (!Hit) + return 0; + + for (Int_t l = 0; l < solutions; l++) { + vecP.clear(); + vecP.push_back(pDet[l]); + vecP.push_back(pD2[l]); + VecVecP.push_back(vecP); + vecWi.push_back(jackPDF[l] * whPS); + } + } else { + TwoBody(S, masses[0], masses[1], p1, p2, rng); + vecP.push_back(p1); + vecP.push_back(p2); + wh = 1.0; + vecWi.push_back(wh * whPS); + VecVecP.push_back(vecP); + } + + return 1; + } + + //========================================================================== + // N-body decay case (nBody > 2) + //========================================================================== + + rng->RndmArray(nBody - 2, rrr.data()); + std::sort(rrr.begin(), rrr.end()); + + for (Int_t i = 0; i < nBody - 2; i++) { + mB = 0.0; + mA = 0.0; + for (Int_t l = 0; l < i + 1; l++) { + mB += masses[l]; + } + for (Int_t l = i + 1; l < nBody; l++) { + mA += masses[l]; + } + mV[i] = rrr[nBody - 2 - i - 1] * (mS - mall) + mA; + + if (i > 0 && (mV[i - 1] - mV[i]) < masses[i]) { + ::Error("TFANG::GenFANG", + "Virtual mass constraint violated at i=%d, mV[i-1]=%g, mV[i]=%g, masses[i]=%g", + i, mV[i - 1], mV[i], masses[i]); + } + } + + whPS = mV[0] * CalcKMFactor(masses[0] * masses[0] / S.M2(), + mV[0] * mV[0] / S.M2()); + + for (Int_t i = 0; i < nBody - 3; i++) { + whPS *= mV[i + 1] * CalcKMFactor(masses[i + 1] * masses[i + 1] / (mV[i] * mV[i]), + mV[i + 1] * mV[i + 1] / (mV[i] * mV[i])); + } + + if (nBody > 2) { + whPS *= CalcKMFactor(masses[nBody - 2] * masses[nBody - 2] / (mV[nBody - 3] * mV[nBody - 3]), + masses[nBody - 1] * masses[nBody - 1] / (mV[nBody - 3] * mV[nBody - 3])); + } + + whPS *= std::pow(kPi, nBody - 1) / 2.0 * + std::pow(mS - mall, nBody - 2) / TMath::Factorial(nBody - 2); + + //========================================================================== + // No detector constraints + //========================================================================== + if (nDet == 0) { + TwoBody(S, masses[0], mV[0], p1, p2, rng); + vecP.push_back(p1); + pV = p2; + + for (Int_t i = 0; i < nBody - 3; i++) { + TwoBody(pV, masses[i + 1], mV[i + 1], p1, p2, rng); + vecP.push_back(p1); + pV = p2; + } + + TwoBody(pV, masses[nBody - 2], masses[nBody - 1], p1, p2, rng); + vecP.push_back(p1); + vecP.push_back(p2); + + wh = 1.0; + vecWi.push_back(wh * whPS); + VecVecP.push_back(vecP); + + return 1; + } + + //========================================================================== + // With detector constraints - use tree to track solutions + //========================================================================== + Node_t *root = CreateFirst(nullptr, S, S, 1.0); + Node_t *cur = root; + Int_t level = 0; + branch.clear(); + + while (level < (nBody - 1)) { + + // Case 1: Constrained particle, not the last two-body decay + if (level < nDet && level < (nBody - 2)) { + pV = cur->fPV; + + if (IsPoint(Ratio[level])) { + V3 = V3Det[level].Unit(); + } else { + TGenVec(Om[level], Ratio[level], V3Det[level].Unit(), V3, rng); + } + + Hit = TGenPointSpace(pV, masses[level], mV[level], V3, + solutions, jackPDF, pDet, pD2); + + if (solutions == 0 && branch.empty()) { + level = nBody - 1; + continue; + } + if (solutions == 0 && !branch.empty()) { + while (level > branch.back()) { + level--; + cur = cur->fParent; + } + cur = cur->fLeft; + level++; + branch.pop_back(); + continue; + } + if (solutions == 1) { + CreateRight(cur, nullptr, pDet[0], pD2[0], jackPDF[0]); + cur = cur->fRight; + level++; + continue; + } else if (solutions == 2) { + branch.push_back(level); + CreateLeft(cur, nullptr, pDet[1], pD2[1], jackPDF[1]); + CreateRight(cur, nullptr, pDet[0], pD2[0], jackPDF[0]); + cur = cur->fRight; + level++; + continue; + } + } + + // Case 2: Constrained particle, last two-body decay + if (level < nDet && level == (nBody - 2)) { + pV = cur->fPV; + + if (IsPoint(Ratio[level])) { + V3 = V3Det[level].Unit(); + } else { + TGenVec(Om[level], Ratio[level], V3Det[level].Unit(), V3, rng); + } + + Hit = TGenPointSpace(pV, masses[level], masses[level + 1], V3, + solutions, jackPDF, pDet, pD2); + + if (solutions == 0 && branch.empty()) { + level = nBody - 1; + continue; + } + if (solutions == 0 && !branch.empty()) { + while (level > branch.back()) { + level--; + cur = cur->fParent; + } + cur = cur->fLeft; + level++; + branch.pop_back(); + continue; + } + if (solutions == 1) { + CreateRight(cur, nullptr, pDet[0], pD2[0], jackPDF[0]); + cur = cur->fRight; + CreateRight(cur, nullptr, pD2[0], S, 1.0); + cur = cur->fRight; + level++; + continue; + } else if (solutions == 2) { + CreateRight(cur, nullptr, pDet[0], pD2[0], jackPDF[0]); + cur = cur->fRight; + CreateRight(cur, nullptr, pD2[0], S, 1.0); + cur = cur->fParent; + CreateLeft(cur, nullptr, pDet[1], pD2[1], jackPDF[1]); + cur = cur->fRight; + CreateRight(cur, nullptr, pD2[1], S, 1.0); + cur = cur->fRight; + level++; + + if (level == (nBody - 1) && !branch.empty()) { + while (level > branch.back()) { + level--; + cur = cur->fParent; + } + cur = cur->fParent; + cur = cur->fLeft; + level++; + branch.pop_back(); + continue; + } + continue; + } + } + + // Case 3: Unconstrained particle, not the last two-body decay + if (level >= nDet && level < nBody - 2) { + pV = cur->fPV; + TwoBody(pV, masses[level], mV[level], p1, p2, rng); + CreateRight(cur, nullptr, p1, p2, 1.0); + cur = cur->fRight; + level++; + continue; + } + + // Case 4: Unconstrained particle, last two-body decay + if (level >= nDet && level == nBody - 2) { + pV = cur->fPV; + TwoBody(pV, masses[level], masses[level + 1], p1, p2, rng); + CreateRight(cur, nullptr, p1, p2, 1.0); + cur = cur->fRight; + CreateRight(cur, nullptr, p2, S, 1.0); + cur = cur->fRight; + level++; + + if (level == (nBody - 1) && !branch.empty()) { + while (level > branch.back()) { + level--; + cur = cur->fParent; + } + cur = cur->fParent; + cur = cur->fLeft; + level++; + branch.pop_back(); + continue; + } + continue; + } + + // Backtrack if needed + if (level == (nBody - 1) && !branch.empty()) { + while (level > branch.back()) { + level--; + cur = cur->fParent; + } + cur = cur->fParent; + cur = cur->fLeft; + level++; + branch.pop_back(); + continue; + } + + } // end while + + CollectPathsWeights(nBody, root, vecJ, pathsJ); + CollectPaths(nBody, root, vecP, pathsP); + + DeleteTree(root); + root = nullptr; + + if (pathsP.empty()) { + return 0; + } + + for (size_t i = 0; i < pathsJ.size(); i++) { + vecJ = pathsJ[i]; + vecP = pathsP[i]; + vecP.erase(vecP.begin()); + VecVecP.push_back(vecP); + + wh = 1.0; + for (Int_t j = 1; j < nDet + 1; j++) { + wh *= vecJ[j]; + } + vecWi.push_back(wh * whPS); + } + + return 1; +} + +} // namespace FANG + +//////////////////////////////////////////////////////////////////////////////// +// TFANG class implementation +//////////////////////////////////////////////////////////////////////////////// + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Default constructor +//////////////////////////////////////////////////////////////////////////////// +TFANG::TFANG() + : fNBody(0) + , fS() + , fMasses() + , fOmega() + , fRatio() + , fV3Det() + , fVecVecP() + , fVecWi() + , fRng(new TRandom3(0)) + , fOwnRng(kTRUE) +{ +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Constructor with external random number generator +//////////////////////////////////////////////////////////////////////////////// +TFANG::TFANG(TRandom3 *rng) + : fNBody(0) + , fS() + , fMasses() + , fOmega() + , fRatio() + , fV3Det() + , fVecVecP() + , fVecWi() + , fRng(rng) + , fOwnRng(kFALSE) +{ + if (fRng == nullptr) { + fRng = new TRandom3(0); + fOwnRng = kTRUE; + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Destructor +//////////////////////////////////////////////////////////////////////////////// +TFANG::~TFANG() +{ + if (fOwnRng && fRng != nullptr) { + delete fRng; + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Set decay configuration +//////////////////////////////////////////////////////////////////////////////// +Bool_t TFANG::SetDecay(const ROOT::Math::PxPyPzMVector &S, Int_t nBody, const Double_t *masses) +{ + if (nBody < 2) { + ::Error("TFANG::SetDecay", "nBody must be >= 2, got %d", nBody); + return kFALSE; + } + + if (masses == nullptr) { + ::Error("TFANG::SetDecay", "masses array is null"); + return kFALSE; + } + + Double_t totalMass = 0.0; + for (Int_t i = 0; i < nBody; ++i) { + if (masses[i] < 0.0) { + ::Error("TFANG::SetDecay", "Negative mass at index %d: %g", i, masses[i]); + return kFALSE; + } + totalMass += masses[i]; + } + + if (totalMass >= S.M()) { + ::Error("TFANG::SetDecay", "Sum of decay masses (%g) >= parent mass (%g)", + totalMass, S.M()); + return kFALSE; + } + + fS = S; + fNBody = nBody; + fMasses.assign(masses, masses + nBody); + + fVecVecP.clear(); + fVecWi.clear(); + + return kTRUE; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Add angular constraint for a particle +//////////////////////////////////////////////////////////////////////////////// +void TFANG::AddConstraint(const ROOT::Math::XYZVector &direction, Double_t omega, Double_t ratio) +{ + if (fNBody == 0) { + ::Warning("TFANG::AddConstraint", "SetDecay must be called before AddConstraint"); + } + + fV3Det.push_back(direction); + fOmega.push_back(omega); + fRatio.push_back(ratio); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Clear all constraints +//////////////////////////////////////////////////////////////////////////////// +void TFANG::ClearConstraints() +{ + fV3Det.clear(); + fOmega.clear(); + fRatio.clear(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Set random number generator seed +//////////////////////////////////////////////////////////////////////////////// +void TFANG::SetSeed(UInt_t seed) +{ + if (fRng != nullptr) { + fRng->SetSeed(seed); + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Generate a single phase space event +//////////////////////////////////////////////////////////////////////////////// +Int_t TFANG::Generate() +{ + if (fNBody < 2) { + ::Error("TFANG::Generate", "SetDecay must be called before Generate"); + return 0; + } + + fVecVecP.clear(); + fVecWi.clear(); + + const Double_t *omPtr = fOmega.empty() ? nullptr : fOmega.data(); + const Double_t *ratioPtr = fRatio.empty() ? nullptr : fRatio.data(); + + Int_t result = FANG::GenFANG(fNBody, fS, fMasses.data(), + omPtr, ratioPtr, fV3Det, + fVecVecP, fVecWi, fRng); + + if (result == 0) { + return 0; + } + + return static_cast(fVecVecP.size()); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Get weight (unconstrained mode) +//////////////////////////////////////////////////////////////////////////////// +Double_t TFANG::GetWeight() const +{ + if (fVecWi.empty()) { + ::Warning("TFANG::GetWeight", "No event generated yet"); + return 0.0; + } + return fVecWi[0]; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Get weight for specific solution +//////////////////////////////////////////////////////////////////////////////// +Double_t TFANG::GetWeight(Int_t iSolution) const +{ + if (iSolution < 0 || iSolution >= static_cast(fVecWi.size())) { + ::Error("TFANG::GetWeight", "Solution index %d out of range [0, %d)", + iSolution, static_cast(fVecWi.size())); + return 0.0; + } + return fVecWi[iSolution]; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Get 4-momentum of particle (single solution) +//////////////////////////////////////////////////////////////////////////////// +ROOT::Math::PxPyPzMVector TFANG::GetDecay(Int_t iParticle) const +{ + return GetDecay(0, iParticle); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Get 4-momentum of particle for specific solution +//////////////////////////////////////////////////////////////////////////////// +ROOT::Math::PxPyPzMVector TFANG::GetDecay(Int_t iSolution, Int_t iParticle) const +{ + if (iSolution < 0 || iSolution >= static_cast(fVecVecP.size())) { + ::Error("TFANG::GetDecay", "Solution index %d out of range [0, %d)", + iSolution, static_cast(fVecVecP.size())); + return ROOT::Math::PxPyPzMVector(); + } + + const auto &vecP = fVecVecP[iSolution]; + if (iParticle < 0 || iParticle >= static_cast(vecP.size())) { + ::Error("TFANG::GetDecay", "Particle index %d out of range [0, %d)", + iParticle, static_cast(vecP.size())); + return ROOT::Math::PxPyPzMVector(); + } + + return vecP[iParticle]; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Get all 4-momenta for a solution +//////////////////////////////////////////////////////////////////////////////// +const std::vector &TFANG::GetDecays(Int_t iSolution) const +{ + static const std::vector empty; + + if (iSolution < 0 || iSolution >= static_cast(fVecVecP.size())) { + ::Error("TFANG::GetDecays", "Solution index %d out of range [0, %d)", + iSolution, static_cast(fVecVecP.size())); + return empty; + } + + return fVecVecP[iSolution]; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate phase space integral with uncertainty +//////////////////////////////////////////////////////////////////////////////// +Bool_t TFANG::GetPhaseSpace(Long64_t nEvents, Double_t &phaseSpace, Double_t &error) const +{ + if (fNBody < 2) { + ::Error("TFANG::GetPhaseSpace", "SetDecay must be called before GetPhaseSpace"); + phaseSpace = 0.0; + error = 0.0; + return kFALSE; + } + + if (nEvents <= 0) { + ::Error("TFANG::GetPhaseSpace", "nEvents must be positive, got %lld", nEvents); + phaseSpace = 0.0; + error = 0.0; + return kFALSE; + } + + TRandom3 localRng(fRng->GetSeed() + 1); + + Double_t sumW = 0.0; + Double_t sumW2 = 0.0; + Long64_t nSuccess = 0; + + std::vector> vecVecP; + std::vector vecWi; + + // Always run unconstrained - pass empty vectors for constraints + std::vector emptyV3Det; + + for (Long64_t i = 0; i < nEvents; ++i) { + vecVecP.clear(); + vecWi.clear(); + + Int_t result = FANG::GenFANG(fNBody, fS, fMasses.data(), + nullptr, nullptr, emptyV3Det, + vecVecP, vecWi, &localRng); + + if (result > 0) { + for (size_t j = 0; j < vecWi.size(); ++j) { + Double_t w = vecWi[j]; + sumW += w; + sumW2 += w * w; + ++nSuccess; + } + } + } + + if (nSuccess == 0) { + ::Warning("TFANG::GetPhaseSpace", "No successful events generated"); + phaseSpace = 0.0; + error = 0.0; + return kFALSE; + } + + Double_t mean = sumW / static_cast(nSuccess); + Double_t variance = (sumW2 / static_cast(nSuccess) - mean * mean); + + phaseSpace = mean; + error = std::sqrt(variance / static_cast(nSuccess)); + + return kTRUE; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate phase space with default 1E6 events +//////////////////////////////////////////////////////////////////////////////// +Bool_t TFANG::GetPhaseSpace(Double_t &phaseSpace, Double_t &error) const +{ + return GetPhaseSpace(1000000, phaseSpace, error); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate partial (constrained) phase space integral with uncertainty +/// +/// Runs Monte Carlo integration to estimate the partial phase space volume +/// when angular constraints are applied. The result is multiplied by the +/// product of all omega values (total solid angle factor). +//////////////////////////////////////////////////////////////////////////////// +Bool_t TFANG::GetPartialPhaseSpace(Long64_t nEvents, Double_t &phaseSpace, Double_t &error) const +{ + if (fNBody < 2) { + ::Error("TFANG::GetPartialPhaseSpace", "SetDecay must be called before GetPartialPhaseSpace"); + phaseSpace = 0.0; + error = 0.0; + return kFALSE; + } + + if (fV3Det.empty()) { + ::Error("TFANG::GetPartialPhaseSpace", "No constraints set. Use AddConstraint() first or use GetPhaseSpace() for unconstrained calculation"); + phaseSpace = 0.0; + error = 0.0; + return kFALSE; + } + + if (nEvents <= 0) { + ::Error("TFANG::GetPartialPhaseSpace", "nEvents must be positive, got %lld", nEvents); + phaseSpace = 0.0; + error = 0.0; + return kFALSE; + } + + // Calculate total solid angle factor based on constraint type + Double_t totalOmega = 1.0; + for (size_t i = 0; i < fOmega.size(); ++i) { + if (FANG::IsPoint(fRatio[i])) { + // Point constraint: factor of 1 (no change) + } else if (FANG::IsRing(fRatio[i])) { + // Ring constraint: factor of 2*pi + totalOmega *= TMath::TwoPi(); + } else { + // Cone constraint: use actual omega value + totalOmega *= fOmega[i]; + } + } + + TRandom3 localRng(fRng->GetSeed() + 1); + + Double_t sumW = 0.0; + Double_t sumW2 = 0.0; + Long64_t nSuccess = 0; + + std::vector> vecVecP; + std::vector vecWi; + + const Double_t *omPtr = fOmega.data(); + const Double_t *ratioPtr = fRatio.data(); + + for (Long64_t i = 0; i < nEvents; ++i) { + vecVecP.clear(); + vecWi.clear(); + + Int_t result = FANG::GenFANG(fNBody, fS, fMasses.data(), + omPtr, ratioPtr, fV3Det, + vecVecP, vecWi, &localRng); + + if (result > 0) { + for (size_t j = 0; j < vecWi.size(); ++j) { + Double_t w = vecWi[j]; + sumW += w; + sumW2 += w * w; + ++nSuccess; + } + } + } + + if (nSuccess == 0) { + ::Warning("TFANG::GetPartialPhaseSpace", "No successful events generated"); + phaseSpace = 0.0; + error = 0.0; + return kFALSE; + } + + Double_t mean = sumW / static_cast(nSuccess); + Double_t variance = (sumW2 / static_cast(nSuccess) - mean * mean); + + // Multiply by total solid angle factor + phaseSpace = totalOmega * mean; + error = totalOmega * std::sqrt(variance / static_cast(nSuccess)); + + return kTRUE; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate partial phase space with default 1E6 events +//////////////////////////////////////////////////////////////////////////////// +Bool_t TFANG::GetPartialPhaseSpace(Double_t &phaseSpace, Double_t &error) const +{ + return GetPartialPhaseSpace(1000000, phaseSpace, error); +} diff --git a/math/physics/src/TGenPhaseSpace.cxx b/math/physics/src/TGenPhaseSpace.cxx index 6413974920072..b7c2699cde632 100644 --- a/math/physics/src/TGenPhaseSpace.cxx +++ b/math/physics/src/TGenPhaseSpace.cxx @@ -2,7 +2,7 @@ // Author: Rene Brun , Valerio Filippini 06/09/2000 /** \class TGenPhaseSpace - \legacy{TGenPhaseSpace, No alternative classes are currently available.} + \legacy{TGenPhaseSpace, Consider using FANG for angular constraints and exact LIPS weights.} \ingroup Physics Utility class to generate n-body event, diff --git a/math/physics/test/CMakeLists.txt b/math/physics/test/CMakeLists.txt new file mode 100644 index 0000000000000..305e92c3c7885 --- /dev/null +++ b/math/physics/test/CMakeLists.txt @@ -0,0 +1,7 @@ +ROOT_ADD_GTEST(testTFANG + testTFANG.cxx + LIBRARIES + Core + MathCore + Physics +) diff --git a/math/physics/test/testTFANG.cxx b/math/physics/test/testTFANG.cxx new file mode 100644 index 0000000000000..97fc45f45a067 --- /dev/null +++ b/math/physics/test/testTFANG.cxx @@ -0,0 +1,622 @@ +// @(#)root/physics:$Id$ +// Author: Arik Kreisel, Itay Horin + +/************************************************************************* + * Copyright (C) 1995-2025, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +//////////////////////////////////////////////////////////////////////////////// +/// \file testTFANG.cxx +/// \ingroup Physics +/// \brief Unit tests for TFANG (Focused Angular N-body event Generator) +/// \authors Arik Kreisel, Itay Horin +/// +/// This file contains gtest unit tests for the TFANG class interface: +/// 1. Full phase space calculation validation using GetPhaseSpace() +/// 2. Partial phase space with detector constraints using GetPartialPhaseSpace() +/// 3. Elastic ep scattering differential cross section vs Rosenbluth formula +/// +/// Reference: Horin, I., Kreisel, A. & Alon, O. Focused angular N -body event generator (FANG). +/// J. High Energ. Phys. 2025, 137 (2025). +/// https://doi.org/10.1007/JHEP12(2025)137 +/// https://arxiv.org/abs/2509.11105 +//////////////////////////////////////////////////////////////////////////////// + +#include "TFANG.h" + +#include "gtest/gtest.h" +#include "TError.h" +#include "TRandom3.h" +#include "TMath.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" + +#include +#include + +using namespace FANG; + +//////////////////////////////////////////////////////////////////////////////// +/// Test fixture for TFANG tests +//////////////////////////////////////////////////////////////////////////////// +class TFANGTest : public ::testing::Test { +protected: + TRandom3 rng; + + void SetUp() override + { + // Set random seed for reproducibility in tests + rng.SetSeed(12345); + } +}; + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test CalcKMFactor with known values +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, CalcKMFactor_KnownValues) +{ + // F(0,0) = sqrt((1-0-0)^2 - 4*0*0) = 1 + EXPECT_DOUBLE_EQ(CalcKMFactor(0.0, 0.0), 1.0); + + // F(0.25, 0.25) = sqrt((1-0.5)^2 - 4*0.0625) = sqrt(0.25 - 0.25) = 0 + EXPECT_NEAR(CalcKMFactor(0.25, 0.25), 0.0, 1e-10); + + // F(0.1, 0.1) = sqrt((0.8)^2 - 0.04) = sqrt(0.64 - 0.04) = sqrt(0.6) + EXPECT_NEAR(CalcKMFactor(0.1, 0.1), std::sqrt(0.6), 1e-10); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TwoBody decay conserves 4-momentum +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, TwoBody_MomentumConservation) +{ + ROOT::Math::PxPyPzMVector S(1.0, 2.0, 3.0, 5.0); + Double_t m1 = 1.0; + Double_t m2 = 1.5; + ROOT::Math::PxPyPzMVector p1, p2; + + TwoBody(S, m1, m2, p1, p2, &rng); + + // Check 4-momentum conservation + ROOT::Math::PxPyPzMVector pSum = p1 + p2; + EXPECT_NEAR(pSum.Px(), S.Px(), 1e-10); + EXPECT_NEAR(pSum.Py(), S.Py(), 1e-10); + EXPECT_NEAR(pSum.Pz(), S.Pz(), 1e-10); + EXPECT_NEAR(pSum.E(), S.E(), 1e-10); + + // Check masses are correct + EXPECT_NEAR(p1.M(), m1, 1e-10); + EXPECT_NEAR(p2.M(), m2, 1e-10); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TwoBody decay produces physical momenta +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, TwoBody_PhysicalMomenta) +{ + ROOT::Math::PxPyPzMVector S(0.0, 0.0, 5.0, 10.0); + Double_t m1 = 2.0; + Double_t m2 = 3.0; + ROOT::Math::PxPyPzMVector p1, p2; + + TwoBody(S, m1, m2, p1, p2, &rng); + + // Check energies are positive and >= mass + EXPECT_GE(p1.E(), m1); + EXPECT_GE(p2.E(), m2); + + // Check 3-momentum magnitudes are positive + EXPECT_GE(p1.P(), 0.0); + EXPECT_GE(p2.P(), 0.0); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG full phase space with known reference value using GetPhaseSpace +/// +/// Uses P(0,0,5,M=12) decaying to 5 particles of mass 1 each. +/// Reference value from FANG paper Table I: 26628.1 ± 3.0 +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, FullPhaseSpace_ReferenceValue) +{ + const Int_t kNBody = 5; + Double_t masses[kNBody] = {1.0, 1.0, 1.0, 1.0, 1.0}; + ROOT::Math::PxPyPzMVector pTotal(0, 0, 5, 12); // Note: E=13 as in paper + + TFANG gen(&rng); + Bool_t valid = gen.SetDecay(pTotal, kNBody, masses); + EXPECT_TRUE(valid) << "SetDecay should succeed for valid configuration"; + + Double_t phaseSpace, error; + const Long64_t nEvents = 1000000; + + Bool_t success = gen.GetPhaseSpace(nEvents, phaseSpace, error); + EXPECT_TRUE(success) << "GetPhaseSpace should succeed"; + + // Reference value from paper: 26628.1 ± 3.0 + // Allow 0.5% tolerance for Monte Carlo fluctuations + Double_t expectedValue = 26628.1; + Double_t tolerance = 0.005 * expectedValue; + + EXPECT_NEAR(phaseSpace, expectedValue, tolerance) + << "Phase space = " << phaseSpace << " +/- " << error + << ", expected = " << expectedValue; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG partial phase space with detector constraints using GetPartialPhaseSpace +/// +/// Tests that constrained particles are within specified solid angles. +/// Reference value 4.764 +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, PartialPhaseSpace_Constraints) +{ + const Int_t kNBody = 5; + Double_t masses[kNBody] = {1.0, 1.0, 1.0, 1.0, 1.0}; + ROOT::Math::PxPyPzMVector pTotal(0, 0, 5, 12); + + TFANG gen(&rng); + Bool_t valid = gen.SetDecay(pTotal, kNBody, masses); + EXPECT_TRUE(valid) << "SetDecay should succeed for valid configuration"; + + // Detector 1: Circle at (0, 0, 0.5), radius 0.2 + ROOT::Math::XYZVector v3_1(0.0, 0.0, 0.5); + Double_t radius1 = TMath::Sqrt(v3_1.Mag2() + 0.2 * 0.2); + Double_t omega1 = kTwoPi * radius1 * (radius1 - v3_1.R()); + gen.AddConstraint(v3_1, omega1, 0.0); // Circle mode + + // Detector 2: Circle at (0.5, 0, 0), radius 0.3 + ROOT::Math::XYZVector v3_2(0.5, 0.0, 0.0); + Double_t radius2 = TMath::Sqrt(v3_2.Mag2() + 0.3 * 0.3); + Double_t omega2 = kTwoPi * radius2 * (radius2 - v3_2.R()); + gen.AddConstraint(v3_2, omega2, 0.0); // Circle mode + + // Detector 3: Strip at (0, 0.5, 0) + ROOT::Math::XYZVector v3_3(0.0, 0.5, 0.0); + Double_t omega3 = 1.2 * kPi; + gen.AddConstraint(v3_3, omega3, 0.4); // Strip mode + + Double_t partialPhaseSpace, error; + const Long64_t nEvents = 1000000; + + Bool_t success = gen.GetPartialPhaseSpace(nEvents, partialPhaseSpace, error); + EXPECT_TRUE(success) << "GetPartialPhaseSpace should succeed"; + + // Reference value: 4.764 + // Allow 5% tolerance for Monte Carlo fluctuations + Double_t expectedValue = 4.764; + Double_t tolerance = 0.05 * expectedValue; + + EXPECT_NEAR(partialPhaseSpace, expectedValue, tolerance) + << "Partial Phase space = " << partialPhaseSpace << " +/- " << error + << ", expected = " << expectedValue; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG Generate verifies momentum conservation +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, Generate_MomentumConservation) +{ + const Int_t kNBody = 5; + Double_t masses[kNBody] = {1.0, 1.0, 1.0, 1.0, 1.0}; + ROOT::Math::PxPyPzMVector pTotal(0, 0, 5, 12); + + TFANG gen(&rng); + gen.SetDecay(pTotal, kNBody, masses); + + // Add constraints + ROOT::Math::XYZVector v3_1(0.0, 0.0, 0.5); + Double_t radius1 = TMath::Sqrt(v3_1.Mag2() + 0.2 * 0.2); + Double_t omega1 = kTwoPi * radius1 * (radius1 - v3_1.R()); + gen.AddConstraint(v3_1, omega1, 0.0); + + ROOT::Math::XYZVector v3_2(0.5, 0.0, 0.0); + Double_t radius2 = TMath::Sqrt(v3_2.Mag2() + 0.3 * 0.3); + Double_t omega2 = kTwoPi * radius2 * (radius2 - v3_2.R()); + gen.AddConstraint(v3_2, omega2, 0.0); + + ROOT::Math::XYZVector v3_3(0.0, 0.5, 0.0); + gen.AddConstraint(v3_3, 1.2 * kPi, 0.4); + + const Int_t nLoop = 1000; + Int_t nSuccess = 0; + + for (Int_t k = 0; k < nLoop; k++) { + Int_t nSolutions = gen.Generate(); + if (nSolutions == 0) continue; + + nSuccess++; + for (Int_t i = 0; i < nSolutions; i++) { + // Verify momentum conservation + ROOT::Math::PxPyPzMVector pSum; + for (Int_t j = 0; j < kNBody; j++) { + pSum = pSum + gen.GetDecay(i, j); + } + EXPECT_NEAR(pSum.Px(), pTotal.Px(), 1e-8); + EXPECT_NEAR(pSum.Py(), pTotal.Py(), 1e-8); + EXPECT_NEAR(pSum.Pz(), pTotal.Pz(), 1e-8); + EXPECT_NEAR(pSum.E(), pTotal.E(), 1e-6); + } + } + + EXPECT_GT(nSuccess, 0) << "Should have at least some successful generations"; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG two-body constrained decay +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, TwoBody_Constrained) +{ + const Int_t kNBody = 2; + Double_t masses[kNBody] = {1.0, 2.0}; + ROOT::Math::PxPyPzMVector pTotal(0, 0, 3, 6); + + TFANG gen(&rng); + Bool_t valid = gen.SetDecay(pTotal, kNBody, masses); + EXPECT_TRUE(valid) << "SetDecay should succeed"; + + // Constrain first particle to z-axis (point mode) + ROOT::Math::XYZVector v3(0.0, 0.0, 1.0); + gen.AddConstraint(v3, 0.0, kModePoint); + + Int_t nSolutions = gen.Generate(); + EXPECT_GE(nSolutions, 1) << "Should have at least one solution"; + + for (Int_t i = 0; i < nSolutions; i++) { + // Check masses + EXPECT_NEAR(gen.GetDecay(i, 0).M(), masses[0], 1e-10); + EXPECT_NEAR(gen.GetDecay(i, 1).M(), masses[1], 1e-10); + + // Check momentum conservation + ROOT::Math::PxPyPzMVector pSum = gen.GetDecay(i, 0) + gen.GetDecay(i, 1); + EXPECT_NEAR(pSum.Px(), pTotal.Px(), 1e-10); + EXPECT_NEAR(pSum.Py(), pTotal.Py(), 1e-10); + EXPECT_NEAR(pSum.Pz(), pTotal.Pz(), 1e-10); + EXPECT_NEAR(pSum.E(), pTotal.E(), 1e-8); + + // First particle should be along z-axis (within numerical precision) + ROOT::Math::PxPyPzMVector p0 = gen.GetDecay(i, 0); + if (p0.P() > 1e-10) { + Double_t cosTheta = p0.Pz() / p0.P(); + EXPECT_NEAR(std::abs(cosTheta), 1.0, 1e-10) + << "Constrained particle should be along z-axis"; + } + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Rosenbluth Cross Section Helper Functions +//////////////////////////////////////////////////////////////////////////////// + +namespace { + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate Rosenbluth cross section for elastic ep scattering +/// \param[in] cosTheta cos(theta) in lab frame +/// \param[in] kineticE electron kinetic energy [GeV] +/// \return Differential cross section dsigma/dOmega [GeV^-2] +//////////////////////////////////////////////////////////////////////////////// +Double_t RosenbluthCrossSection(Double_t cosTheta, Double_t kineticE) +{ + Double_t sigma = 0.0; + const Double_t alpha = 1.0 / 137.0; + const Double_t massProton = 0.938272029; + const Double_t massElectron = 0.000511; + + ROOT::Math::XYZVector vDir(TMath::Sqrt(1.0 - cosTheta * cosTheta), 0.0, cosTheta); + + ROOT::Math::PxPyPzMVector pProton(0.0, 0.0, 0.0, massProton); + Double_t gamma = kineticE / massElectron + 1.0; + Double_t beta = TMath::Sqrt(1.0 - 1.0 / (gamma * gamma)); + ROOT::Math::PxPyPzMVector pElectron(0.0, 0.0, gamma * beta * massElectron, massElectron); + + ROOT::Math::PxPyPzMVector pTotal = pProton + pElectron; + + LongDouble_t massCM = pTotal.M(); + LongDouble_t energyCM = pTotal.E(); + LongDouble_t momCM = pTotal.P(); + LongDouble_t energyCM3 = (massCM * massCM - massProton * massProton + + massElectron * massElectron) / (2.0 * massCM); + + LongDouble_t aa = momCM * momCM * cosTheta * cosTheta - energyCM * energyCM; + LongDouble_t bb = 2.0 * momCM * cosTheta * energyCM3 * massCM; + LongDouble_t cc = energyCM3 * massCM * energyCM3 * massCM - + massElectron * massElectron * energyCM * energyCM; + + if (bb * bb - 4.0 * aa * cc < 0.0) { + return 0.0; + } + + LongDouble_t momLAB = (-bb + TMath::Sqrt(bb * bb - 4.0 * aa * cc)) / (2.0 * aa); + if (momLAB > 0.0) { + ROOT::Math::PxPyPzMVector pElectronOut(momLAB * vDir.X(), momLAB * vDir.Y(), + momLAB * vDir.Z(), massElectron); + ROOT::Math::PxPyPzMVector pMomentumTransfer = pElectronOut - pElectron; + Double_t qSquared = -pMomentumTransfer.M2(); + Double_t formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * + (1.0 + qSquared / kDipoleMassSq)); + Double_t formGM = kProtonMagneticMoment * formGE; + Double_t tau = qSquared / (4.0 * massProton * massProton); + Double_t mottXS = alpha * alpha / + (pElectron.E() * pElectron.E() * (1.0 - cosTheta) * (1.0 - cosTheta)) * + pElectronOut.E() / pElectron.E() * (1.0 + cosTheta) / 2.0; + sigma = mottXS * ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) + + 2.0 * tau * formGM * formGM * (1.0 - cosTheta) / (1.0 + cosTheta)); + } + + momLAB = (-bb - TMath::Sqrt(bb * bb - 4.0 * aa * cc)) / (2.0 * aa); + if (momLAB > 0.0) { + ROOT::Math::PxPyPzMVector pElectronOut(momLAB * vDir.X(), momLAB * vDir.Y(), + momLAB * vDir.Z(), massElectron); + ROOT::Math::PxPyPzMVector pMomentumTransfer = pElectronOut - pElectron; + Double_t qSquared = -pMomentumTransfer.M2(); + Double_t formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * + (1.0 + qSquared / kDipoleMassSq)); + Double_t formGM = kProtonMagneticMoment * formGE; + Double_t tau = qSquared / (4.0 * massProton * massProton); + Double_t mottXS = alpha * alpha / + (pElectron.E() * pElectron.E() * (1.0 - cosTheta) * (1.0 - cosTheta)) * + pElectronOut.E() / pElectron.E() * (1.0 + cosTheta) / 2.0; + sigma += mottXS * ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) + + 2.0 * tau * formGM * formGM * (1.0 - cosTheta) / (1.0 + cosTheta)); + } + + return sigma; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate FANG cross section at a specific angle using TFANG class +/// \param[in] cosTheta cos(theta) in lab frame +/// \param[in] kineticE electron kinetic energy [GeV] +/// \param[in] nLoop number of Monte Carlo iterations +/// \param[out] error statistical error estimate +/// \param[in] rng pointer to random number generator +/// \return Differential cross section dsigma/dOmega [GeV^-2] +//////////////////////////////////////////////////////////////////////////////// +Double_t TFANGCrossSection(Double_t cosTheta, Double_t kineticE, Int_t nLoop, Double_t &error, TRandom3 *rng) +{ + const Int_t kNBody = 2; + const Double_t massElectron = 0.000511; + const Double_t massProton = 0.938272029; + const Double_t alphaQED = 1.0 / 137.0; + + Double_t masses[kNBody] = {massElectron, massProton}; + + // Setup kinematics + ROOT::Math::PxPyPzMVector pTarget(0.0, 0.0, 0.0, massProton); + Double_t gamma = kineticE / massElectron + 1.0; + Double_t beta = TMath::Sqrt(1.0 - 1.0 / (gamma * gamma)); + ROOT::Math::PxPyPzMVector pBeam(0.0, 0.0, gamma * beta * massElectron, massElectron); + ROOT::Math::PxPyPzMVector pTotal = pBeam + pTarget; + + Double_t flux = 1.0 / (16.0 * kPi * kPi * + TMath::Sqrt(pBeam.Dot(pTarget) * pBeam.Dot(pTarget) - + massElectron * massElectron * massProton * massProton)); + + // Setup TFANG generator + TFANG gen(rng); + gen.SetDecay(pTotal, kNBody, masses); + + // Add point constraint for electron direction + ROOT::Math::XYZVector v3(TMath::Sqrt(1.0 - cosTheta * cosTheta), 0.0, cosTheta); + gen.AddConstraint(v3, 0.0, kModePoint); + + Double_t sumW = 0.0; + Double_t sumW2 = 0.0; + Int_t nEvents = 0; + + for (Int_t k = 0; k < nLoop; k++) { + Int_t nSolutions = gen.Generate(); + if (nSolutions == 0) continue; + + for (Int_t i = 0; i < nSolutions; i++) { + Double_t weight = gen.GetWeight(i); + + ROOT::Math::PxPyPzMVector pElectronOut = gen.GetDecay(i, 0); + ROOT::Math::PxPyPzMVector pMomTransfer = pBeam - pElectronOut; + ROOT::Math::PxPyPzMVector pU = pTarget - pElectronOut; + Double_t qSquared = -pMomTransfer.M2(); + + Double_t formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * + (1.0 + qSquared / kDipoleMassSq)); + Double_t formGM = kProtonMagneticMoment * formGE; + Double_t tau = qSquared / (4.0 * massProton * massProton); + Double_t lambda = (pTotal.M2() - pU.M2()) / (4.0 * massProton * massProton); + + Double_t ampSquared = 16.0 * kPi * kPi * alphaQED * alphaQED / (tau * tau) * + ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) * + (lambda * lambda - tau * tau - tau) + + 2.0 * tau * tau * formGM * formGM); + + weight *= ampSquared; + nEvents++; + sumW += weight; + sumW2 += weight * weight; + } + } + + error = flux * TMath::Sqrt(sumW2) / nEvents; + return flux * sumW / nEvents; +} + +} // anonymous namespace + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG differential cross section against Rosenbluth formula +/// +/// Tests elastic ep scattering at 3 GeV for various angles. +/// Skips cos(theta) = ±1 where numerical issues may occur. +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, Rosenbluth_ElasticEP) +{ + const Double_t kineticE = 3.0; // GeV + const Int_t nLoop = 50000; + + // Test angles: cos(theta) from -0.8 to 0.8 (skip ±1) + std::vector testAngles = {-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}; + + for (Double_t cosTheta : testAngles) { + Double_t rosenbluth = RosenbluthCrossSection(cosTheta, kineticE); + Double_t fangError; + Double_t fang = TFANGCrossSection(cosTheta, kineticE, nLoop, fangError, &rng); + + // Calculate ratio + Double_t ratio = fang / rosenbluth; + + // Allow 10% tolerance for Monte Carlo fluctuations at this statistics + EXPECT_NEAR(ratio, 1.0, 0.10) + << "cos(theta) = " << cosTheta + << ": TFANG = " << fang << " +/- " << fangError + << ", Rosenbluth = " << rosenbluth + << ", ratio = " << ratio; + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG cross section precision at a single angle +/// +/// Uses higher statistics to verify agreement with Rosenbluth within 5%. +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, Rosenbluth_HighPrecision) +{ + const Double_t kineticE = 3.0; // GeV + const Double_t cosTheta = 0.0; // 90 degree scattering + const Int_t nLoop = 100000; + + Double_t rosenbluth = RosenbluthCrossSection(cosTheta, kineticE); + Double_t fangError; + Double_t fang = TFANGCrossSection(cosTheta, kineticE, nLoop, fangError, &rng); + + Double_t ratio = fang / rosenbluth; + + // At high statistics, expect agreement within 5% + EXPECT_NEAR(ratio, 1.0, 0.05) + << "High precision test at cos(theta) = 0" + << ": TFANG = " << fang << " +/- " << fangError + << ", Rosenbluth = " << rosenbluth + << ", ratio = " << ratio; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test that TFANG returns 0 solutions for unphysical configurations +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, UnphysicalConfiguration) +{ + // Temporarily suppress error messages since we expect an error condition + Int_t oldLevel = gErrorIgnoreLevel; + gErrorIgnoreLevel = kFatal; // Only show Fatal messages + + const Int_t kNBody = 3; + Double_t masses[kNBody] = {5.0, 5.0, 5.0}; // Total mass = 15 + ROOT::Math::PxPyPzMVector pTotal(0, 0, 0, 10); // M = 10 < 15, unphysical + + TFANG gen(&rng); + Bool_t valid = gen.SetDecay(pTotal, kNBody, masses); + + gErrorIgnoreLevel = oldLevel; // Restore previous error level + EXPECT_FALSE(valid) << "SetDecay should fail for unphysical mass configuration"; +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test mode detection functions +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, ModeDetection) +{ + EXPECT_TRUE(IsPoint(kModePoint)); + EXPECT_FALSE(IsPoint(0.0)); + EXPECT_FALSE(IsPoint(0.5)); + EXPECT_FALSE(IsPoint(-1.0)); + + EXPECT_TRUE(IsCircle(kModeCircle)); + EXPECT_FALSE(IsCircle(kModePoint)); + EXPECT_FALSE(IsCircle(0.5)); + EXPECT_FALSE(IsCircle(-1.0)); + + EXPECT_TRUE(IsStrip(0.5)); + EXPECT_TRUE(IsStrip(0.1)); + EXPECT_TRUE(IsStrip(1.0)); + EXPECT_FALSE(IsStrip(0.0)); + EXPECT_FALSE(IsStrip(kModePoint)); + EXPECT_FALSE(IsStrip(-0.5)); + + EXPECT_TRUE(IsRing(-0.5)); + EXPECT_TRUE(IsRing(-1.0)); + EXPECT_FALSE(IsRing(0.0)); + EXPECT_FALSE(IsRing(0.5)); + EXPECT_FALSE(IsRing(kModePoint)); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG utility methods +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, UtilityMethods) +{ + const Int_t kNBody = 3; + Double_t masses[kNBody] = {1.0, 1.0, 1.0}; + ROOT::Math::PxPyPzMVector pTotal(0, 0, 0, 5); + + TFANG gen(&rng); + gen.SetDecay(pTotal, kNBody, masses); + + // Test unconstrained state + EXPECT_FALSE(gen.IsConstrained()); + EXPECT_EQ(gen.GetNConstraints(), 0); + EXPECT_EQ(gen.GetNBody(), kNBody); + + // Add constraint and verify + ROOT::Math::XYZVector v3(0.0, 0.0, 1.0); + gen.AddConstraint(v3, 0.1, 0.0); + EXPECT_TRUE(gen.IsConstrained()); + EXPECT_EQ(gen.GetNConstraints(), 1); + + // Clear constraints + gen.ClearConstraints(); + EXPECT_FALSE(gen.IsConstrained()); + EXPECT_EQ(gen.GetNConstraints(), 0); +} + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Test TFANG unconstrained generation +//////////////////////////////////////////////////////////////////////////////// +TEST_F(TFANGTest, UnconstrainedGeneration) +{ + const Int_t kNBody = 3; + Double_t masses[kNBody] = {0.139, 0.139, 0.139}; // pion masses + ROOT::Math::PxPyPzMVector pTotal(0, 0, 0, 1.0); // 1 GeV at rest + + TFANG gen(&rng); + gen.SetDecay(pTotal, kNBody, masses); + + // Generate events and verify properties + Int_t nSuccess = 0; + const Int_t nLoop = 100; + + for (Int_t k = 0; k < nLoop; k++) { + Int_t nSolutions = gen.Generate(); + if (nSolutions == 0) continue; + + nSuccess++; + EXPECT_EQ(nSolutions, 1) << "Unconstrained should have 1 solution"; + + // Check weight is positive + Double_t weight = gen.GetWeight(); + EXPECT_GT(weight, 0.0) << "Weight should be positive"; + + // Verify momentum conservation + ROOT::Math::PxPyPzMVector pSum; + for (Int_t j = 0; j < kNBody; j++) { + ROOT::Math::PxPyPzMVector p = gen.GetDecay(j); + pSum = pSum + p; + EXPECT_NEAR(p.M(), masses[j], 1e-10) << "Mass should match"; + } + EXPECT_NEAR(pSum.Px(), pTotal.Px(), 1e-10); + EXPECT_NEAR(pSum.Py(), pTotal.Py(), 1e-10); + EXPECT_NEAR(pSum.Pz(), pTotal.Pz(), 1e-10); + EXPECT_NEAR(pSum.E(), pTotal.E(), 1e-8); + } + + EXPECT_EQ(nSuccess, nLoop) << "All unconstrained generations should succeed"; +} diff --git a/montecarlo/CMakeLists.txt b/montecarlo/CMakeLists.txt index feaa74ad7cf5f..d4906974ba6f3 100644 --- a/montecarlo/CMakeLists.txt +++ b/montecarlo/CMakeLists.txt @@ -3,7 +3,6 @@ # # For the licensing terms see $ROOTSYS/LICENSE. # For the list of contributors see $ROOTSYS/README/CREDITS. - add_subdirectory(eg) if(pythia8) add_subdirectory(pythia8) diff --git a/tutorials/evegen/runTFANG.C b/tutorials/evegen/runTFANG.C new file mode 100644 index 0000000000000..469efc6ee160e --- /dev/null +++ b/tutorials/evegen/runTFANG.C @@ -0,0 +1,839 @@ +// @(#)root/fang:$Id$ +// Author: Arik Kreisel + +//////////////////////////////////////////////////////////////////////////////// +/// \file runTFANG.C +/// \ingroup Physics +/// \brief Demonstration and validation of FANG using the TFANG class interface +/// \author Arik Kreisel +/// +/// TFANG is a Monte Carlo tool for efficient event generation in restricted +/// (or full) Lorentz-Invariant Phase Space (LIPS). Unlike conventional approaches +/// that always sample the full 4pi solid angle, TFANG can also directly generate +/// events in which selected final-state particles are constrained to fixed +/// directions or finite angular regions in the laboratory frame. +/// +/// Reference: Horin, I., Kreisel, A. & Alon, O. Focused angular N -body event generator (FANG). +/// J. High Energ. Phys. 2025, 137 (2025). +/// https://doi.org/10.1007/JHEP12(2025)137 +/// https://arxiv.org/abs/2509.11105 +/// +/// This file contains: +/// 1. Rosenbluth cross section function for elastic ep scattering +/// 2. runTFANG() - main demonstration function that validates TFANG against: +/// - Full phase space calculation +/// - Partial phase space with detector constraints (vs TFANG unconstrained with cuts) +/// - Partial phase space with detector constraints (vs TGenPhaseSpace) - optional +/// - Elastic ep differential cross section (vs Rosenbluth formula) +//////////////////////////////////////////////////////////////////////////////// + +#include "TFANG.h" +#include "TStyle.h" +#include "TCanvas.h" +#include "TH1D.h" +#include "TH1F.h" +#include "TF1.h" +#include "TLegend.h" +#include "TGraphErrors.h" +#include "TGenPhaseSpace.h" +#include "TLorentzVector.h" +#include "TVector3.h" +#include "TRandom3.h" + +//////////////////////////////////////////////////////////////////////////////// +// Configuration: Set to false to skip TGenPhaseSpace comparison +//////////////////////////////////////////////////////////////////////////////// +const Bool_t kRunTGenPhaseSpace = true; + +//////////////////////////////////////////////////////////////////////////////// +// Rosenbluth Cross Section for Elastic ep Scattering +//////////////////////////////////////////////////////////////////////////////// + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Calculate Rosenbluth cross section for elastic ep scattering +/// +/// This is a ROOT TF1-compatible function that calculates the differential +/// cross section dsigma/dOmega for elastic electron-proton scattering. +/// +/// \param[in] x Array where x[0] = cos(theta_lab) +/// \param[in] par Array where par[0] = electron kinetic energy W [GeV] +/// \return Differential cross section in GeV^-2 +//////////////////////////////////////////////////////////////////////////////// +Double_t fElastic(Double_t *x, Double_t *par) +{ + using namespace FANG; + + Double_t sigma = 0.0; + Double_t alpha = 1.0 / 137.0; + + // Direction unit vector for scattering angle + ROOT::Math::XYZVector vDir(TMath::Sqrt(1.0 - x[0] * x[0]), 0.0, x[0]); + + // Particle masses + Double_t massProton = 0.938272029; // proton mass in GeV + Double_t massElectron = 0.000511; // electron mass in GeV + + // Setup kinematics + ROOT::Math::PxPyPzMVector pProton(0.0, 0.0, 0.0, massProton); // proton at rest + Double_t kineticE = par[0]; // electron kinetic energy + Double_t gamma = kineticE / massElectron + 1.0; + Double_t beta = TMath::Sqrt(1.0 - 1.0 / (gamma * gamma)); + ROOT::Math::PxPyPzMVector pElectron(0.0, 0.0, gamma * beta * massElectron, massElectron); + + ROOT::Math::PxPyPzMVector pElectronOut, pMomentumTransfer; + ROOT::Math::PxPyPzMVector pTotal = pProton + pElectron; // total 4-momentum + + Double_t mottXS, tau, formGE, formGM, qSquared; + + // CM frame quantities + LongDouble_t massCM = pTotal.M(); + LongDouble_t energyCM = pTotal.E(); + LongDouble_t momCM = pTotal.P(); + LongDouble_t energyCM3 = (massCM * massCM - massProton * massProton + + massElectron * massElectron) / (2.0 * massCM); + + // Quadratic equation coefficients + LongDouble_t aa = momCM * momCM * x[0] * x[0] - energyCM * energyCM; + LongDouble_t bb = 2.0 * momCM * x[0] * energyCM3 * massCM; + LongDouble_t cc = energyCM3 * massCM * energyCM3 * massCM - + massElectron * massElectron * energyCM * energyCM; + + // Check for physical solutions + if (bb * bb - 4.0 * aa * cc < 0.0) { + return 0.0; + } + + // First solution + LongDouble_t momLAB = (-bb + TMath::Sqrt(bb * bb - 4.0 * aa * cc)) / (2.0 * aa); + if (momLAB > 0.0) { + pElectronOut.SetCoordinates(momLAB * vDir.X(), momLAB * vDir.Y(), + momLAB * vDir.Z(), massElectron); + pMomentumTransfer = pElectronOut - pElectron; + qSquared = -pMomentumTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + mottXS = alpha * alpha / (pElectron.E() * pElectron.E() * (1.0 - x[0]) * (1.0 - x[0])) * + pElectronOut.E() / pElectron.E() * (1.0 + x[0]) / 2.0; + sigma = mottXS * ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) + + 2.0 * tau * formGM * formGM * (1.0 - x[0]) / (1.0 + x[0])); + } + + // Second solution + momLAB = (-bb - TMath::Sqrt(bb * bb - 4.0 * aa * cc)) / (2.0 * aa); + if (momLAB > 0.0) { + pElectronOut.SetCoordinates(momLAB * vDir.X(), momLAB * vDir.Y(), + momLAB * vDir.Z(), massElectron); + pMomentumTransfer = pElectronOut - pElectron; + qSquared = -pMomentumTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + mottXS = alpha * alpha / (pElectron.E() * pElectron.E() * (1.0 - x[0]) * (1.0 - x[0])) * + pElectronOut.E() / pElectron.E() * (1.0 + x[0]) / 2.0; + sigma += mottXS * ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) + + 2.0 * tau * formGM * formGM * (1.0 - x[0]) / (1.0 + x[0])); + } + + return sigma; +} + +//////////////////////////////////////////////////////////////////////////////// +// Main Demonstration Function +//////////////////////////////////////////////////////////////////////////////// + +//////////////////////////////////////////////////////////////////////////////// +/// \brief Main demonstration and validation function for TFANG +/// +/// Performs three validation tests: +/// 1. Full phase space calculation for 5-body decay +/// 2. Partial phase space with 3 detector constraints, compared to: +/// - TFANG unconstrained with geometric cuts +/// - TGenPhaseSpace with cuts (if kRunTGenPhaseSpace is true) +/// 3. Elastic ep scattering differential cross section vs Rosenbluth formula +//////////////////////////////////////////////////////////////////////////////// +void runTFANG() +{ + using namespace FANG; + + gStyle->SetOptStat(0); + + Int_t nEvents = 0; + + //========================================================================== + // Setup for 5-body decay test + //========================================================================== + const Int_t kNBody = 5; + Double_t masses[kNBody] = {1.0, 1.0, 1.0, 1.0, 1.0}; // outgoing masses + ROOT::Math::PxPyPzMVector pTotal(0, 0, 5, 12); // total 4-momentum + + Double_t weight = 0.0; + Double_t sumW = 0.0; + Double_t sumW2 = 0.0; + + //========================================================================== + // Test 1: TFANG Full Phase Space Calculation + //========================================================================== + std::cout << "========================================" << std::endl; + std::cout << "Test 1: Full Phase Space Calculation (TFANG)" << std::endl; + std::cout << "========================================" << std::endl; + + Int_t nLoop = 1E6; + + // Create TFANG generator for unconstrained phase space + TFANG genFull; + genFull.SetDecay(pTotal, kNBody, masses); + + for (Int_t k = 0; k < nLoop; k++) { + if (genFull.Generate() == 0) continue; + + for (Int_t i = 0; i < genFull.GetNSolutions(); i++) { + weight = genFull.GetWeight(i); + nEvents++; + sumW += weight; + sumW2 += weight * weight; + } + } + Double_t mean = sumW / nEvents; + Double_t variance = sumW2 / nEvents - mean * mean; + + // Also get phase space using GetPhaseSpace for comparison + Double_t phaseSpace, phaseSpaceErr; + genFull.GetPhaseSpace(static_cast(nLoop), phaseSpace, phaseSpaceErr); + + std::cout << "nEvents = " << nEvents << std::endl; + std::cout << "Total Phase Space from user loop= " << mean + << " +/- " << TMath::Sqrt(variance / nEvents) << std::endl; + std::cout << "Total Phase Space from GetPhaseSpace = " << phaseSpace + << " +/- " << phaseSpaceErr << std::endl; + + //========================================================================== + // Test 2: Partial Phase Space with Detector Constraints + //========================================================================== + std::cout << "\n========================================" << std::endl; + std::cout << "Test 2: Partial Phase Space (TFANG)" << std::endl; + std::cout << " - TFANG constrained vs TFANG unconstrained with cuts" << std::endl; + if (kRunTGenPhaseSpace) { + std::cout << " - TFANG constrained vs TGenPhaseSpace with cuts" << std::endl; + } + std::cout << "========================================" << std::endl; + + const Int_t kNDet = 3; + Double_t omega[kNDet]; + Double_t shape[kNDet]; + + // Detector positions and radii + Double_t detPosX[kNDet - 1] = {0.0, 0.5}; + Double_t detPosY[kNDet - 1] = {0.0, 0.0}; + Double_t detPosZ[kNDet - 1] = {0.5, 0.0}; + Double_t detRadius[kNDet - 1] = {0.2, 0.3}; + + ROOT::Math::XYZVector v3; + std::vector v3Det; + Double_t radius; + Double_t totalOmega = 1.0; + + // Setup first two detectors (circular) + for (Int_t i = 0; i < kNDet - 1; i++) { + v3.SetXYZ(detPosX[i], detPosY[i], detPosZ[i]); + v3Det.push_back(v3); + radius = TMath::Sqrt(v3.Mag2() + detRadius[i] * detRadius[i]); + omega[i] = kTwoPi * radius * (radius - v3.R()); + shape[i] = 0.0; // Circle generation + totalOmega *= omega[i]; + } + + // Setup third detector (strip) + v3.SetXYZ(0, 0.5, 0); + v3Det.push_back(v3); + omega[2] = 1.2 * kPi; + shape[2] = 0.4; // Strip generation + totalOmega *= omega[2]; + + std::cout << "Detector configurations:" << std::endl; + std::cout << " Det 1: Circle, Omega = " << omega[0] << " sr" << std::endl; + std::cout << " Det 2: Circle, Omega = " << omega[1] << " sr" << std::endl; + std::cout << " Det 3: Strip, Omega = " << omega[2] << " sr" << std::endl; + std::cout << " Total solid angle factor = " << totalOmega << std::endl; + + // Calculate total available mass for kinetic energy + Double_t totalMass = 0.0; + for (Int_t l = 0; l < kNBody; l++) { + totalMass += masses[l]; + } + + // Create histograms for TFANG constrained results + TH1D *hFangE[kNBody]; + TH1D *hFangCos[kNBody]; + TH1D *hFangPhi[kNBody]; + TH1D *hFullE[kNBody]; + TH1D *hFullCos[kNBody]; + + for (Int_t i = 0; i < kNBody; i++) { + hFangE[i] = new TH1D(Form("hFangE_%d", i), "", 100, 0, pTotal.E() - totalMass); + hFangE[i]->SetMarkerStyle(20); + hFangE[i]->SetLineColor(6); + hFangE[i]->SetMinimum(0); + hFangE[i]->GetXaxis()->SetTitle(Form("p_{%d} Energy", i + 1)); + hFangE[i]->GetXaxis()->SetTitleSize(0.07); + hFangE[i]->GetXaxis()->SetLabelSize(0.06); + hFangE[i]->GetYaxis()->SetLabelSize(0.05); + hFangE[i]->GetYaxis()->SetTitle("Events"); + hFangE[i]->GetYaxis()->SetTitleSize(0.07); + hFangE[i]->GetYaxis()->SetTitleOffset(0.5); + hFangE[i]->GetXaxis()->SetTitleOffset(0.9); + + hFangCos[i] = new TH1D(Form("hFangCos_%d", i), "", 50, -1, 1); + hFangCos[i]->SetMarkerStyle(20); + hFangCos[i]->SetLineColor(6); + hFangCos[i]->SetMinimum(0); + hFangCos[i]->GetXaxis()->SetTitle(Form("p_{%d} cos(#theta)", i + 1)); + hFangCos[i]->GetXaxis()->SetTitleSize(0.07); + hFangCos[i]->SetTitleOffset(0.7); + hFangCos[i]->GetYaxis()->SetTitle("Events"); + hFangCos[i]->GetYaxis()->SetTitleSize(0.07); + hFangCos[i]->GetYaxis()->SetTitleOffset(0.5); + hFangCos[i]->GetXaxis()->SetLabelSize(0.06); + hFangCos[i]->GetYaxis()->SetLabelSize(0.05); + hFangCos[i]->GetXaxis()->SetTitleOffset(0.9); + + hFangPhi[i] = new TH1D(Form("hFangPhi_%d", i), "", 50, -kPi, kPi); + hFangPhi[i]->SetMarkerStyle(20); + hFangPhi[i]->SetLineColor(6); + hFangPhi[i]->SetMinimum(0); + hFangPhi[i]->GetXaxis()->SetTitle(Form("p_{%d} #varphi", i + 1)); + hFangPhi[i]->GetXaxis()->SetTitleSize(0.07); + hFangPhi[i]->SetTitleOffset(0.7); + hFangPhi[i]->GetYaxis()->SetTitle("Events"); + hFangPhi[i]->GetYaxis()->SetTitleSize(0.07); + hFangPhi[i]->GetYaxis()->SetTitleOffset(0.5); + hFangPhi[i]->GetXaxis()->SetLabelSize(0.06); + hFangPhi[i]->GetYaxis()->SetLabelSize(0.05); + hFangPhi[i]->GetXaxis()->SetTitleOffset(0.9); + + hFullE[i] = new TH1D(Form("hFullE_%d", i), "hFullE", 100, 0, pTotal.E() - totalMass); + hFullCos[i] = new TH1D(Form("hFullCos_%d", i), "hFullCos", 50, -1, 1); + hFullE[i]->SetMarkerStyle(20); + hFullCos[i]->SetMarkerStyle(20); + } + + // Create histograms for TFANG unconstrained with cuts comparison + TH1D *hFangCutsE[kNBody]; + TH1D *hFangCutsCos[kNBody]; + TH1D *hFangCutsPhi[kNBody]; + + for (Int_t i = 0; i < kNBody; i++) { + hFangCutsE[i] = new TH1D(Form("hFangCutsE_%d", i), "hFangCutsE", 100, 0, pTotal.E() - totalMass); + hFangCutsCos[i] = new TH1D(Form("hFangCutsCos_%d", i), "hFangCutsCos", 50, -1, 1); + hFangCutsPhi[i] = new TH1D(Form("hFangCutsPhi_%d", i), "hFangCutsPhi", 50, -kPi, kPi); + hFangCutsE[i]->SetMarkerStyle(21); + hFangCutsE[i]->SetMarkerColor(kBlue); + hFangCutsCos[i]->SetMarkerStyle(21); + hFangCutsCos[i]->SetMarkerColor(kBlue); + hFangCutsPhi[i]->SetMarkerStyle(21); + hFangCutsPhi[i]->SetMarkerColor(kBlue); + } + + // Create histograms for TGenPhaseSpace (CERN/GENBOD) comparison + TH1D *hGenbodE[kNBody]; + TH1D *hGenbodCos[kNBody]; + TH1D *hGenbodPhi[kNBody]; + + for (Int_t i = 0; i < kNBody; i++) { + hGenbodE[i] = new TH1D(Form("hGenbodE_%d", i), "hGenbodE", 100, 0, pTotal.E() - totalMass); + hGenbodCos[i] = new TH1D(Form("hGenbodCos_%d", i), "hGenbodCos", 50, -1, 1); + hGenbodPhi[i] = new TH1D(Form("hGenbodPhi_%d", i), "hGenbodPhi", 50, -kPi, kPi); + hGenbodE[i]->SetMarkerStyle(20); + hGenbodCos[i]->SetMarkerStyle(20); + hGenbodPhi[i]->SetMarkerStyle(20); + } + + // Create TFANG generator with detector constraints + TFANG genConstrained; + genConstrained.SetDecay(pTotal, kNBody, masses); + for (Int_t i = 0; i < kNDet; i++) { + genConstrained.AddConstraint(v3Det[i], omega[i], shape[i]); + } + + // Run TFANG with detector constraints + weight = 0.0; + sumW = 0.0; + sumW2 = 0.0; + nEvents = 0; + nLoop = 1E6; + + TH1D *hWeight = new TH1D("hWeight", "hWeight", 100, 0, 10); + + for (Int_t k = 0; k < nLoop; k++) { + if (genConstrained.Generate() == 0) continue; + + for (Int_t i = 0; i < genConstrained.GetNSolutions(); i++) { + weight = genConstrained.GetWeight(i); + nEvents++; + sumW += weight; + sumW2 += weight * weight; + + for (Int_t j = 0; j < kNBody; j++) { + ROOT::Math::PxPyPzMVector p = genConstrained.GetDecay(i, j); + hFangE[j]->Fill(p.E() - masses[j], weight * totalOmega); + hFangCos[j]->Fill(TMath::Cos(p.Theta()), weight * totalOmega); + hFangPhi[j]->Fill(p.Phi(), weight * totalOmega); + } + } + } + + // get partial phaseSpace with GetPartialPhaseSpace for comparison + genConstrained.GetPartialPhaseSpace(static_cast(nLoop), phaseSpace, phaseSpaceErr); + + std::cout << "\nTFANG Constrained Results:" << std::endl; + std::cout << " nEvents = " << nEvents << std::endl; + mean = sumW / nEvents; + variance = sumW2 / nEvents - mean * mean; + std::cout << " Partial Phase Space from user loop= " << totalOmega * mean + << " +/- " << totalOmega * TMath::Sqrt(variance / nEvents) << std::endl; + std::cout << " Partial Phase Space from GetPartialPhaseSpace = " <Integral() = " << hFangE[0]->Integral() << std::endl; + + // Draw TFANG results + TCanvas *c1 = new TCanvas("c1", "c1 En", 10, 10, 1800, 1500); + c1->Divide(2, static_cast(TMath::Floor(kNBody / 2.0 + 0.6))); + for (Int_t i = 0; i < kNBody; i++) { + c1->cd(i + 1); + gPad->SetBottomMargin(0.15); + hFangE[i]->Draw("hist"); + } + + TCanvas *c2 = new TCanvas("c2", "c2 cos", 10, 10, 1800, 1500); + c2->Divide(2, static_cast(TMath::Floor(kNBody / 2.0 + 0.6))); + for (Int_t i = 0; i < kNBody; i++) { + c2->cd(i + 1); + gPad->SetBottomMargin(0.15); + hFangCos[i]->Draw("hist"); + } + + TCanvas *c3 = new TCanvas("c3", "c3 phi", 10, 10, 1800, 1500); + c3->Divide(2, static_cast(TMath::Floor(kNBody / 2.0 + 0.6))); + for (Int_t i = 0; i < kNBody; i++) { + c3->cd(i + 1); + gPad->SetBottomMargin(0.15); + hFangPhi[i]->Draw("hist"); + } + + //========================================================================== + // TFANG Unconstrained with Cuts Comparison + //========================================================================== + std::cout << "\n--- TFANG Unconstrained with Cuts ---" << std::endl; + + // Direction vectors for cut comparison + TVector3 tv3[kNDet]; + for (Int_t i = 0; i < kNDet; i++) { + tv3[i].SetXYZ(v3Det[i].X(), v3Det[i].Y(), v3Det[i].Z()); + tv3[i] = tv3[i].Unit(); + } + + Double_t scaleFactor = 10.0; // Need more events since most will be rejected by cuts + Int_t outsideCut = 0; + Int_t nPassedCuts = 0; + Int_t nTotalGenerated = 0; + + // Create unconstrained TFANG generator + TFANG genUnconstrained; + genUnconstrained.SetDecay(pTotal, kNBody, masses); + + for (Int_t k = 0; k < nLoop * scaleFactor; k++) { + // Generate unconstrained events + if (genUnconstrained.Generate() == 0) continue; + + nTotalGenerated++; + + for (Int_t i = 0; i < genUnconstrained.GetNSolutions(); i++) { + weight = genUnconstrained.GetWeight(i)/scaleFactor; + outsideCut = 0; + + // Apply geometric cuts (same as TGenPhaseSpace comparison) + for (Int_t j = 0; j < kNDet; j++) { + ROOT::Math::PxPyPzMVector p = genUnconstrained.GetDecay(i, j); + TVector3 pVec(p.Px(), p.Py(), p.Pz()); + + if (shape[j] == 0.0 && + (1.0 - TMath::Cos(tv3[j].Angle(pVec))) > omega[j] / kTwoPi) { + outsideCut = 1; + } + if (shape[j] > 0.0 && + (TMath::Abs(tv3[j].Phi() - p.Phi()) > kPi * shape[j] || + TMath::Abs(TMath::Cos(tv3[j].Theta()) - TMath::Cos(p.Theta())) > + omega[j] / (4.0 * kPi * shape[j]))) { + outsideCut = 1; + } + } + + if (outsideCut == 1) continue; + + nPassedCuts++; + + for (Int_t j = 0; j < kNBody; j++) { + ROOT::Math::PxPyPzMVector p = genUnconstrained.GetDecay(i, j); + hFangCutsE[j]->Fill(p.E() - masses[j], weight); + hFangCutsCos[j]->Fill(TMath::Cos(p.Theta()), weight); + hFangCutsPhi[j]->Fill(p.Phi(), weight); + } + } + } + + std::cout << " Total events generated: " << nTotalGenerated << std::endl; + std::cout << " Events passing cuts: " << nPassedCuts << std::endl; + if (nTotalGenerated > 0) { + std::cout << " Cut efficiency: " << 100.0 * nPassedCuts / nTotalGenerated << "%" << std::endl; + } + std::cout << " hFangCutsE[0]->Integral() = " << hFangCutsE[0]->Integral() << std::endl; + + //========================================================================== + // TGenPhaseSpace comparison (GENBOD with cuts) - Optional + //========================================================================== + if (kRunTGenPhaseSpace) { + std::cout << "\n--- TGenPhaseSpace (GENBOD) with Cuts ---" << std::endl; + + TLorentzVector pTotalCern; + pTotalCern.SetPxPyPzE(0, 0, 5, 13); + TGenPhaseSpace genPhaseSpace; + Double_t genWeight; + genPhaseSpace.SetDecay(pTotalCern, kNBody, masses); + + Double_t normFactor = 2050032.6; // Normalizing factor + + for (Int_t k = 0; k < nLoop * scaleFactor; k++) { + genWeight = genPhaseSpace.Generate() / scaleFactor * normFactor; + outsideCut = 0; + + // Apply geometric cuts + for (Int_t i = 0; i < kNDet; i++) { + if (shape[i] == 0.0 && + (1.0 - TMath::Cos(tv3[i].Angle(genPhaseSpace.GetDecay(i)->Vect()))) > omega[i] / kTwoPi) { + outsideCut = 1; + } + if (shape[i] > 0.0 && + (TMath::Abs(tv3[i].Phi() - genPhaseSpace.GetDecay(i)->Phi()) > kPi * shape[i] || + TMath::Abs(TMath::Cos(tv3[i].Theta()) - TMath::Cos(genPhaseSpace.GetDecay(i)->Theta())) > + omega[i] / (4.0 * kPi * shape[i]))) { + outsideCut = 1; + } + } + + if (outsideCut == 1) continue; + + for (Int_t i = 0; i < kNBody; i++) { + hGenbodE[i]->Fill(genPhaseSpace.GetDecay(i)->E() - masses[i], genWeight); + hGenbodCos[i]->Fill(TMath::Cos(genPhaseSpace.GetDecay(i)->Theta()), genWeight); + hGenbodPhi[i]->Fill(genPhaseSpace.GetDecay(i)->Phi(), genWeight); + } + } + + std::cout << " hGenbodE[0]->Integral() = " << hGenbodE[0]->Integral() << std::endl; + } + + // Setup legends + TLegend *leg[3 * kNBody]; + for (Int_t i = 0; i < kNBody * 3; i++) { + leg[i] = new TLegend(0.52, 0.62, 0.85, 0.88); + } + + // Adjust legend positions for some plots + leg[10] = new TLegend(0.12, 0.12, 0.45, 0.38); + leg[11] = new TLegend(0.56, 0.62, 0.89, 0.88); + leg[12] = new TLegend(0.12, 0.62, 0.45, 0.88); + for (Int_t i = 5; i <= 9; i++) { + leg[i] = new TLegend(0.12, 0.62, 0.45, 0.88); + } + + for (Int_t i = 0; i < kNBody; i++) { + leg[i]->AddEntry(hFangE[i], "TFANG constrained", "l"); + leg[i]->AddEntry(hFangCutsE[i], "TFANG unconstrained with cuts", "p"); + if (kRunTGenPhaseSpace) { + leg[i]->AddEntry(hGenbodE[i], "GENBOD with cuts", "p"); + } + + leg[i + kNBody]->AddEntry(hFangCos[i], "TFANG constrained", "l"); + leg[i + kNBody]->AddEntry(hFangCutsCos[i], "TFANG unconstrained with cuts", "p"); + if (kRunTGenPhaseSpace) { + leg[i + kNBody]->AddEntry(hGenbodCos[i], "GENBOD with cuts", "p"); + } + + leg[i + 2 * kNBody]->AddEntry(hFangPhi[i], "TFANG constrained", "l"); + leg[i + 2 * kNBody]->AddEntry(hFangCutsPhi[i], "TFANG unconstrained with cuts", "p"); + if (kRunTGenPhaseSpace) { + leg[i + 2 * kNBody]->AddEntry(hGenbodPhi[i], "GENBOD with cuts", "p"); + } + } + + // Overlay comparison results + for (Int_t i = 0; i < kNBody; i++) { + c1->cd(i + 1); + hFangCutsE[i]->DrawCopy("ep same"); + if (kRunTGenPhaseSpace) { + hGenbodE[i]->DrawCopy("ep same"); + } + leg[i]->Draw(); + + c2->cd(i + 1); + hFangCutsCos[i]->DrawCopy("ep same"); + if (kRunTGenPhaseSpace) { + hGenbodCos[i]->DrawCopy("ep same"); + } + leg[i + kNBody]->Draw(); + + c3->cd(i + 1); + hFangCutsPhi[i]->DrawCopy("ep same"); + if (kRunTGenPhaseSpace) { + hGenbodPhi[i]->DrawCopy("ep same"); + } + leg[i + 2 * kNBody]->Draw(); + } + + //========================================================================== + // Test 3: Elastic ep Scattering Cross Section + //========================================================================== + std::cout << "\n========================================" << std::endl; + std::cout << "Test 3: Elastic ep Differential Cross Section (TFANG)" << std::endl; + std::cout << "========================================" << std::endl; + + const Int_t kNBody2 = 2; + nLoop = 1E5; + nEvents = 0; + + Double_t massElectron = 0.000511; // GeV + Double_t massProton = 0.938272029; // proton mass in GeV + + // Setup kinematics + ROOT::Math::PxPyPzMVector pTarget(0.0, 0.0, 0.0, massProton); + Double_t kineticE = 3.0; // GeV electron kinetic energy + Double_t gamma = kineticE / massElectron + 1.0; + Double_t beta = TMath::Sqrt(1.0 - 1.0 / (gamma * gamma)); + ROOT::Math::PxPyPzMVector pBeam(0.0, 0.0, gamma * beta * massElectron, massElectron); + ROOT::Math::PxPyPzMVector pTotal2 = pBeam + pTarget; + + Double_t masses2[kNBody2] = {massElectron, massProton}; + + Double_t alphaQED = 1.0 / 137.0; + Double_t ampSquared = 0.0; + weight = 0.0; + sumW = 0.0; + sumW2 = 0.0; + + ROOT::Math::PxPyPzMVector pProtonIn, pElectronIn, pProtonOut, pElectronOut, pMomTransfer; + Double_t lambda, tau, formGE, formGM, qSquared; + + pElectronIn = pBeam; + pProtonIn = pTarget; + Double_t flux = 1.0 / (16.0 * kPi * kPi * + TMath::Sqrt(pElectronIn.Dot(pProtonIn) * pElectronIn.Dot(pProtonIn) - + massElectron * massElectron * massProton * massProton)); + + // Setup Rosenbluth function for comparison + TF1 *fRosenbluth = new TF1("fRosenbluth", fElastic, -1, 0.9992, 1); + Double_t parElastic[1] = {kineticE}; + fRosenbluth->SetParameters(parElastic); + + //========================================================================== + // TFANG Point Generation: Differential Cross Section at Specific Angles + //========================================================================== + Double_t sigmaArr[11]; + Double_t sigmaErrArr[11]; + Double_t cosThetaArr[11]; + Double_t cosThetaErrArr[11]; + Double_t cosTheta; + + std::cout << "\nCalculating differential cross section at specific angles:" << std::endl; + + for (Int_t l = 0; l < 11; l++) { + ampSquared = 0.0; + weight = 0.0; + sumW = 0.0; + sumW2 = 0.0; + nEvents = 0; + + cosTheta = -0.99 + l * 0.2; + if (l == 10) cosTheta = 0.95; + cosThetaArr[l] = cosTheta; + cosThetaErrArr[l] = 0.0; + + v3.SetXYZ(TMath::Sqrt(1.0 - cosTheta * cosTheta), 0.0, cosTheta); + + // Create TFANG generator with point constraint + TFANG genPoint; + genPoint.SetDecay(pTotal2, kNBody2, masses2); + genPoint.AddConstraint(v3, 1.0, kModePoint); + + for (Int_t k = 0; k < nLoop; k++) { + if (genPoint.Generate() == 0) continue; + + for (Int_t i = 0; i < genPoint.GetNSolutions(); i++) { + weight = genPoint.GetWeight(i); + pElectronOut = genPoint.GetDecay(i, 0); + pProtonOut = genPoint.GetDecay(i, 1); + pMomTransfer = pElectronIn - pElectronOut; + ROOT::Math::PxPyPzMVector pU = pTarget - pElectronOut; + qSquared = -pMomTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * + (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + lambda = (pTotal2.M2() - pU.M2()) / (4.0 * massProton * massProton); + + // Calculate squared amplitude + ampSquared = 16.0 * kPi * kPi * alphaQED * alphaQED / (tau * tau) * + ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) * + (lambda * lambda - tau * tau - tau) + + 2.0 * tau * tau * formGM * formGM); + + weight *= ampSquared; + nEvents++; + sumW += weight; + sumW2 += weight * weight; + } + } + + sigmaArr[l] = flux * sumW / nEvents; + sigmaErrArr[l] = flux * TMath::Sqrt(sumW2) / nEvents; + + std::cout << " cos(theta) = " << cosTheta + << ": dsigma/dOmega = " << sigmaArr[l] << " +/- " << sigmaErrArr[l] + << " (TFANG/Rosenbluth = " << sigmaArr[l] / fRosenbluth->Eval(cosTheta) << ")" + << std::endl; + } + + TGraphErrors *grElastic = new TGraphErrors(11, cosThetaArr, sigmaArr, cosThetaErrArr, sigmaErrArr); + grElastic->SetMarkerStyle(20); + grElastic->SetMarkerSize(1.3); + + //========================================================================== + // TFANG Event Generation: Full Angular Distribution + //========================================================================== + std::cout << "\nGenerating full angular distribution..." << std::endl; + + Double_t sinTheta, phi, r1; + + TH1D *hXsec = new TH1D("hXsec", "hXsec", 440, -1.1, 1.1); + TH1D *hNorm = new TH1D("hNorm", "hNorm", 440, -1.1, 1.1); + TH1D *hCount = new TH1D("hCount", "hCount", 440, -1.1, 1.1); + TH1D *hError = new TH1D("hError", "hError", 440, -1.1, 1.1); + hXsec->SetMinimum(1E-17); + hNorm->SetMinimum(1E-17); + hError->SetMinimum(0.999); + + // Generate events in multiple r1 ranges to cover full angular range + // Using 1/r^2 distribution to importance-sample forward angles + + struct Range_t { + Double_t fC1; + Double_t fC2; + }; + Range_t ranges[4] = {{1.0, 2.0}, {0.4, 1.0}, {0.12, 0.4}, {0.01, 0.12}}; + + // Create TFANG generator once outside the loop + TFANG genAngular; + genAngular.SetDecay(pTotal2, kNBody2, masses2); + + for (Int_t rangeIdx = 0; rangeIdx < 4; rangeIdx++) { + Double_t c1 = ranges[rangeIdx].fC1; + Double_t c2 = ranges[rangeIdx].fC2; + + std::cout << " Range " << rangeIdx + 1 << ": cos(theta) in [" + << 1.0 - c2 << ", " << 1.0 - c1 << "]" << std::endl; + + for (Int_t k = 0; k < 1000000; k++) { + // Generate r1 with 1/r^2 distribution + r1 = c1 * c1 * c2 / (c2 * c1 - gRandom->Uniform(0, 1) * c1 * (c2 - c1)); + cosTheta = 1.0 - r1; + sinTheta = TMath::Sqrt(1.0 - cosTheta * cosTheta); + phi = gRandom->Uniform(0, kTwoPi); + + v3.SetXYZ(sinTheta * TMath::Cos(phi), sinTheta * TMath::Sin(phi), cosTheta); + + // Update constraint for new direction + genAngular.ClearConstraints(); + genAngular.AddConstraint(v3, 1.0, kModePoint); + + if (genAngular.Generate() == 0) continue; + + for (Int_t i = 0; i < genAngular.GetNSolutions(); i++) { + weight = genAngular.GetWeight(i); + pElectronOut = genAngular.GetDecay(i, 0); + pProtonOut = genAngular.GetDecay(i, 1); + pMomTransfer = pElectronIn - pElectronOut; + ROOT::Math::PxPyPzMVector pU = pTarget - pElectronOut; + qSquared = -pMomTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * + (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + lambda = (pTotal2.M2() - pU.M2()) / (4.0 * massProton * massProton); + + ampSquared = 16.0 * kPi * kPi * alphaQED * alphaQED / (tau * tau) * + ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) * + (lambda * lambda - tau * tau - tau) + + 2.0 * tau * tau * formGM * formGM); + + weight *= ampSquared; + + // Reweight from 1/r^2 to flat distribution + Double_t reweight = r1 * r1 * (c2 - c1) / c1 / c2; + hXsec->Fill(TMath::Cos(pElectronOut.Theta()), reweight * weight); + hNorm->Fill(TMath::Cos(pElectronOut.Theta()), reweight); + hCount->Fill(TMath::Cos(pElectronOut.Theta()), 1.0); + } + } + } + + // Scale and compute errors + Double_t scaleXsec = flux * 2.0 / hXsec->GetBinWidth(2) / hNorm->Integral(); + + for (Int_t l = 1; l <= hXsec->GetNbinsX(); l++) { + Double_t signal = hXsec->GetBinContent(l); + Double_t error = hXsec->GetBinError(l); + Double_t entries = hCount->GetBinContent(l); + + if (entries > 0) { + hError->SetBinContent(l, signal / error / TMath::Sqrt(entries)); + } + hXsec->SetBinContent(l, signal * scaleXsec); + hXsec->SetBinError(l, error * scaleXsec); + } + + // Configure histogram appearance + hXsec->SetYTitle("#frac{d#sigma}{d#Omega}(ep -> ep) [GeV^{-2}]"); + hXsec->SetXTitle("cos(#theta_{LAB})"); + hXsec->SetTitle("Electron Energy E=3 GeV"); + + // Create final comparison plot + TLegend *legFinal = new TLegend(0.12, 0.68, 0.42, 0.88); + + TCanvas *cFinal = new TCanvas("cFinal", "cFinal his", 10, 10, 1800, 1500); + gPad->SetLogy(); + TH1F *vFrame = gPad->DrawFrame(-1.2, 5E-10, 1, 5E-3); + hXsec->Draw("hist same"); + grElastic->Draw("P"); + fRosenbluth->Draw("same"); + + legFinal->AddEntry(hXsec, "TFANG event generation", "l"); + legFinal->AddEntry(grElastic, "TFANG point calculation", "p"); + legFinal->AddEntry(fRosenbluth, "Rosenbluth cross section", "l"); + legFinal->Draw(); + + vFrame->SetYTitle("#frac{d#sigma}{d#Omega}(ep -> ep) [GeV^{-2}]"); + vFrame->SetXTitle("cos(#theta_{LAB})"); + vFrame->SetTitle("Electron Energy E=3 GeV"); + + // Additional diagnostic plots + TCanvas *cDiag = new TCanvas("cDiag", "cDiag Wi error", 10, 10, 1800, 1500); + cDiag->Divide(2, 1); + cDiag->cd(1); + hNorm->Draw("hist"); + cDiag->cd(2); + hCount->Draw("hist"); + + std::cout << "\n========================================" << std::endl; + std::cout << "runTFANG() completed successfully" << std::endl; + std::cout << "J. High Energ. Phys. 2025, 137 (2025)" << std::endl; + std::cout << "https://doi.org/10.1007/JHEP12(2025)137" << std::endl; + std::cout << "========================================" << std::endl; +} diff --git a/tutorials/evegen/runTFANG_parallel.C b/tutorials/evegen/runTFANG_parallel.C new file mode 100644 index 0000000000000..3e777c9241414 --- /dev/null +++ b/tutorials/evegen/runTFANG_parallel.C @@ -0,0 +1,1170 @@ +// @(#)root/fang:$Id$ +// Author: Arik Kreisel +// Parallelized version using std::thread (no OpenMP required) + +//////////////////////////////////////////////////////////////////////////////// +/// \file runTFANG_parallel.C +/// \ingroup Physics +/// \brief Parallelized demonstration and validation of FANG using the TFANG class +/// \author Arik Kreisel +/// +/// This is the multi-threaded version of runTFANG.C using std::thread. +/// All event generation loops are parallelized for improved performance. +/// +/// TFANG is a Monte Carlo tool for efficient event generation in restricted +/// (or full) Lorentz-Invariant Phase Space (LIPS). +/// +/// Reference: Horin, I., Kreisel, A. & Alon, O. Focused angular N -body event generator (FANG). +/// J. High Energ. Phys. 2025, 137 (2025). +/// https://doi.org/10.1007/JHEP12(2025)137 +/// https://arxiv.org/abs/2509.11105 +/// +/// Features: +/// - Comparison of TFANG constrained vs TFANG unconstrained with cuts +//////////////////////////////////////////////////////////////////////////////// + +#include "TFANG.h" +#include "TStyle.h" +#include "TCanvas.h" +#include "TH1D.h" +#include "TH1F.h" +#include "TF1.h" +#include "TLegend.h" +#include "TGraphErrors.h" +#include "TLorentzVector.h" +#include "TVector3.h" +#include "TRandom3.h" +#include "TROOT.h" + +// Threading includes +#include +#include +#include +#include +#include +#include + +//////////////////////////////////////////////////////////////////////////////// +// Thread-safe Work Queue +//////////////////////////////////////////////////////////////////////////////// + +class WorkQueue { +public: + void Push(Int_t item) { + std::lock_guard lock(fMutex); + fQueue.push(item); + fCV.notify_one(); + } + + bool Pop(Int_t& item) { + std::lock_guard lock(fMutex); + if (fQueue.empty()) return false; + item = fQueue.front(); + fQueue.pop(); + return true; + } + + bool Empty() { + std::lock_guard lock(fMutex); + return fQueue.empty(); + } + +private: + std::queue fQueue; + std::mutex fMutex; + std::condition_variable fCV; +}; + +//////////////////////////////////////////////////////////////////////////////// +// Structure to hold accumulated results from workers +//////////////////////////////////////////////////////////////////////////////// + +struct AccumulatorResult { + Double_t fSumW; + Double_t fSumW2; + Int_t fNEvents; + + AccumulatorResult() : fSumW(0.0), fSumW2(0.0), fNEvents(0) {} +}; + +//////////////////////////////////////////////////////////////////////////////// +// Structure to hold TFANG unconstrained with cuts results +//////////////////////////////////////////////////////////////////////////////// + +struct TFangCutsResult { + Int_t fNTotalGenerated; + Int_t fNPassedCuts; + + TFangCutsResult() : fNTotalGenerated(0), fNPassedCuts(0) {} +}; + +//////////////////////////////////////////////////////////////////////////////// +// Structure to hold point calculation results +//////////////////////////////////////////////////////////////////////////////// + +struct PointResult { + Double_t fCosTheta; + Double_t fSigma; + Double_t fSigmaErr; +}; + +//////////////////////////////////////////////////////////////////////////////// +// Rosenbluth Cross Section for Elastic ep Scattering +//////////////////////////////////////////////////////////////////////////////// + +Double_t fElastic(Double_t *x, Double_t *par) +{ + using namespace FANG; + + Double_t sigma = 0.0; + Double_t alpha = 1.0 / 137.0; + + ROOT::Math::XYZVector vDir(TMath::Sqrt(1.0 - x[0] * x[0]), 0.0, x[0]); + + Double_t massProton = 0.938272029; + Double_t massElectron = 0.000511; + + ROOT::Math::PxPyPzMVector pProton(0.0, 0.0, 0.0, massProton); + Double_t kineticE = par[0]; + Double_t gamma = kineticE / massElectron + 1.0; + Double_t beta = TMath::Sqrt(1.0 - 1.0 / (gamma * gamma)); + ROOT::Math::PxPyPzMVector pElectron(0.0, 0.0, gamma * beta * massElectron, massElectron); + + ROOT::Math::PxPyPzMVector pElectronOut, pMomentumTransfer; + ROOT::Math::PxPyPzMVector pTotal = pProton + pElectron; + + Double_t mottXS, tau, formGE, formGM, qSquared; + + LongDouble_t massCM = pTotal.M(); + LongDouble_t energyCM = pTotal.E(); + LongDouble_t momCM = pTotal.P(); + LongDouble_t energyCM3 = (massCM * massCM - massProton * massProton + + massElectron * massElectron) / (2.0 * massCM); + + LongDouble_t aa = momCM * momCM * x[0] * x[0] - energyCM * energyCM; + LongDouble_t bb = 2.0 * momCM * x[0] * energyCM3 * massCM; + LongDouble_t cc = energyCM3 * massCM * energyCM3 * massCM - + massElectron * massElectron * energyCM * energyCM; + + if (bb * bb - 4.0 * aa * cc < 0.0) { + return 0.0; + } + + LongDouble_t momLAB = (-bb + TMath::Sqrt(bb * bb - 4.0 * aa * cc)) / (2.0 * aa); + if (momLAB > 0.0) { + pElectronOut.SetCoordinates(momLAB * vDir.X(), momLAB * vDir.Y(), + momLAB * vDir.Z(), massElectron); + pMomentumTransfer = pElectronOut - pElectron; + qSquared = -pMomentumTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + mottXS = alpha * alpha / (pElectron.E() * pElectron.E() * (1.0 - x[0]) * (1.0 - x[0])) * + pElectronOut.E() / pElectron.E() * (1.0 + x[0]) / 2.0; + sigma = mottXS * ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) + + 2.0 * tau * formGM * formGM * (1.0 - x[0]) / (1.0 + x[0])); + } + + momLAB = (-bb - TMath::Sqrt(bb * bb - 4.0 * aa * cc)) / (2.0 * aa); + if (momLAB > 0.0) { + pElectronOut.SetCoordinates(momLAB * vDir.X(), momLAB * vDir.Y(), + momLAB * vDir.Z(), massElectron); + pMomentumTransfer = pElectronOut - pElectron; + qSquared = -pMomentumTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + mottXS = alpha * alpha / (pElectron.E() * pElectron.E() * (1.0 - x[0]) * (1.0 - x[0])) * + pElectronOut.E() / pElectron.E() * (1.0 + x[0]) / 2.0; + sigma += mottXS * ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) + + 2.0 * tau * formGM * formGM * (1.0 - x[0]) / (1.0 + x[0])); + } + + return sigma; +} + +//////////////////////////////////////////////////////////////////////////////// +// Worker Functions for Parallel Loops +//////////////////////////////////////////////////////////////////////////////// + +//------------------------------------------------------------------------------ +// Worker for Test 1: Full Phase Space Calculation +//------------------------------------------------------------------------------ +void WorkerTest1( + Int_t threadId, + WorkQueue& workQueue, + const Int_t kNBody, + const ROOT::Math::PxPyPzMVector& pTotal, + const Double_t* masses, + std::vector& results, + std::mutex& resultsMutex +) { + using namespace FANG; + + AccumulatorResult localResult; + + // Thread-local RNG with unique seed + TRandom3 rng(threadId + 100); + + // Thread-local TFANG instance with thread-local RNG + TFANG gen(&rng); + gen.SetDecay(pTotal, kNBody, masses); + + Double_t weight; + Int_t workItem; + + while (workQueue.Pop(workItem)) { + if (gen.Generate() == 0) continue; + + for (Int_t i = 0; i < gen.GetNSolutions(); i++) { + weight = gen.GetWeight(i); + localResult.fNEvents++; + localResult.fSumW += weight; + localResult.fSumW2 += weight * weight; + } + } + + // Store results + { + std::lock_guard lock(resultsMutex); + results.push_back(localResult); + } +} + +//------------------------------------------------------------------------------ +// Worker for Test 2: TFANG with Detector Constraints +//------------------------------------------------------------------------------ +void WorkerTest2TFANG( + Int_t threadId, + WorkQueue& workQueue, + const Int_t kNBody, + const Int_t kNDet, + const ROOT::Math::PxPyPzMVector& pTotal, + const Double_t* masses, + const Double_t* omega, + const Double_t* shape, + const std::vector& v3DetConst, + Double_t totalOmega, + std::vector& histsE, + std::vector& histsCos, + std::vector& histsPhi, + std::vector& results, + std::mutex& resultsMutex +) { + using namespace FANG; + + AccumulatorResult localResult; + + // Thread-local RNG with unique seed + TRandom3 rng(threadId + 200); + + // Thread-local TFANG instance + TFANG gen(&rng); + gen.SetDecay(pTotal, kNBody, masses); + for (Int_t i = 0; i < kNDet; i++) { + gen.AddConstraint(v3DetConst[i], omega[i], shape[i]); + } + + Double_t weight; + Int_t workItem; + + while (workQueue.Pop(workItem)) { + if (gen.Generate() == 0) continue; + + for (Int_t i = 0; i < gen.GetNSolutions(); i++) { + weight = gen.GetWeight(i); + localResult.fNEvents++; + localResult.fSumW += weight; + localResult.fSumW2 += weight * weight; + + for (Int_t j = 0; j < kNBody; j++) { + ROOT::Math::PxPyPzMVector p = gen.GetDecay(i, j); + histsE[threadId * kNBody + j]->Fill(p.E() - masses[j], weight * totalOmega); + histsCos[threadId * kNBody + j]->Fill(TMath::Cos(p.Theta()), weight * totalOmega); + histsPhi[threadId * kNBody + j]->Fill(p.Phi(), weight * totalOmega); + } + } + } + + { + std::lock_guard lock(resultsMutex); + results.push_back(localResult); + } +} + +//------------------------------------------------------------------------------ +// Worker for Test 2: TFANG Unconstrained with Cuts +//------------------------------------------------------------------------------ +void WorkerTest2TFANGCuts( + Int_t threadId, + WorkQueue& workQueue, + const Int_t kNBody, + const Int_t kNDet, + const ROOT::Math::PxPyPzMVector& pTotal, + const Double_t* masses, + const Double_t* omega, + const Double_t* shape, + const std::vector& tv3, + Double_t scaleFactor, + std::vector& histsE, + std::vector& histsCos, + std::vector& histsPhi, + std::vector& results, + std::mutex& resultsMutex +) { + using namespace FANG; + + TFangCutsResult localResult; + + // Thread-local RNG with unique seed + TRandom3 rng(threadId + 300); + + // Thread-local unconstrained TFANG instance + TFANG gen(&rng); + gen.SetDecay(pTotal, kNBody, masses); + // No constraints added for unconstrained generation + + Double_t weight; + Int_t outsideCut; + Int_t workItem; + + while (workQueue.Pop(workItem)) { + if (gen.Generate() == 0) continue; + + localResult.fNTotalGenerated++; + + for (Int_t i = 0; i < gen.GetNSolutions(); i++) { + weight = gen.GetWeight(i); + outsideCut = 0; + + // Apply geometric cuts + for (Int_t j = 0; j < kNDet; j++) { + ROOT::Math::PxPyPzMVector p = gen.GetDecay(i, j); + TVector3 pVec(p.Px(), p.Py(), p.Pz()); + + if (shape[j] == 0.0 && + (1.0 - TMath::Cos(tv3[j].Angle(pVec))) > omega[j] / kTwoPi) { + outsideCut = 1; + } + if (shape[j] > 0.0 && + (TMath::Abs(tv3[j].Phi() - p.Phi()) > kPi * shape[j] || + TMath::Abs(TMath::Cos(tv3[j].Theta()) - TMath::Cos(p.Theta())) > + omega[j] / (4.0 * kPi * shape[j]))) { + outsideCut = 1; + } + } + + if (outsideCut == 1) continue; + + localResult.fNPassedCuts++; + + for (Int_t j = 0; j < kNBody; j++) { + ROOT::Math::PxPyPzMVector p = gen.GetDecay(i, j); + histsE[threadId * kNBody + j]->Fill(p.E() - masses[j], weight / scaleFactor); + histsCos[threadId * kNBody + j]->Fill(TMath::Cos(p.Theta()), weight / scaleFactor); + histsPhi[threadId * kNBody + j]->Fill(p.Phi(), weight / scaleFactor); + } + } + } + + { + std::lock_guard lock(resultsMutex); + results.push_back(localResult); + } +} + +//------------------------------------------------------------------------------ +// Worker for Test 3: Point generation at specific angles +//------------------------------------------------------------------------------ +void WorkerTest3Point( + Int_t threadId, + WorkQueue& workQueue, + Int_t nLoop, + Double_t massElectron, + Double_t massProton, + const ROOT::Math::PxPyPzMVector& pTotal2, + const ROOT::Math::PxPyPzMVector& pTarget, + const ROOT::Math::PxPyPzMVector& pElectronIn, + Double_t flux, + std::vector& pointResults, + std::mutex& resultsMutex +) { + using namespace FANG; + + const Int_t kNBody2 = 2; + Double_t masses2[kNBody2] = {massElectron, massProton}; + Double_t alphaQED = 1.0 / 137.0; + + // Thread-local RNG + TRandom3 rng(threadId + 500); + + ROOT::Math::XYZVector v3; + ROOT::Math::PxPyPzMVector pElectronOut, pProtonOut, pMomTransfer; + Double_t qSquared, formGE, formGM, tau, lambda, ampSquared, weight; + + Int_t angleIdx; + while (workQueue.Pop(angleIdx)) { + Double_t sumW = 0.0; + Double_t sumW2 = 0.0; + Int_t nEvents = 0; + + Double_t cosTheta = -0.99 + angleIdx * 0.2; + if (angleIdx == 10) cosTheta = 0.95; + + v3.SetXYZ(TMath::Sqrt(1.0 - cosTheta * cosTheta), 0.0, cosTheta); + + // Thread-local TFANG instance with point constraint + TFANG gen(&rng); + gen.SetDecay(pTotal2, kNBody2, masses2); + gen.AddConstraint(v3, 1.0, kModePoint); + + for (Int_t k = 0; k < nLoop; k++) { + if (gen.Generate() == 0) continue; + + for (Int_t i = 0; i < gen.GetNSolutions(); i++) { + weight = gen.GetWeight(i); + pElectronOut = gen.GetDecay(i, 0); + pProtonOut = gen.GetDecay(i, 1); + pMomTransfer = pElectronIn - pElectronOut; + ROOT::Math::PxPyPzMVector pU = pTarget - pElectronOut; + qSquared = -pMomTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * + (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + lambda = (pTotal2.M2() - pU.M2()) / (4.0 * massProton * massProton); + + ampSquared = 16.0 * kPi * kPi * alphaQED * alphaQED / (tau * tau) * + ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) * + (lambda * lambda - tau * tau - tau) + + 2.0 * tau * tau * formGM * formGM); + + weight *= ampSquared; + nEvents++; + sumW += weight; + sumW2 += weight * weight; + } + } + + PointResult result; + result.fCosTheta = cosTheta; + result.fSigma = (nEvents > 0) ? flux * sumW / nEvents : 0.0; + result.fSigmaErr = (nEvents > 0) ? flux * TMath::Sqrt(sumW2) / nEvents : 0.0; + + { + std::lock_guard lock(resultsMutex); + pointResults.push_back(result); + } + } +} + +//------------------------------------------------------------------------------ +// Worker for Test 3: Full angular distribution +//------------------------------------------------------------------------------ +void WorkerTest3Angular( + Int_t threadId, + WorkQueue& workQueue, + Int_t nEventsPerRange, + Double_t massElectron, + Double_t massProton, + const ROOT::Math::PxPyPzMVector& pTotal2, + const ROOT::Math::PxPyPzMVector& pTarget, + const ROOT::Math::PxPyPzMVector& pElectronIn, + Double_t c1, + Double_t c2, + TH1D* hXsec, + TH1D* hNorm, + TH1D* hCount +) { + using namespace FANG; + + // Thread-local random generator for direction sampling + TRandom3 rng(threadId + 600); + + const Int_t kNBody2 = 2; + Double_t masses2[kNBody2] = {massElectron, massProton}; + Double_t alphaQED = 1.0 / 137.0; + + // Thread-local TFANG instance + TFANG gen(&rng); + gen.SetDecay(pTotal2, kNBody2, masses2); + + ROOT::Math::XYZVector v3; + ROOT::Math::PxPyPzMVector pElectronOut, pProtonOut, pMomTransfer; + Double_t qSquared, formGE, formGM, tau, lambda, ampSquared, weight; + Double_t r1, cosTheta, sinTheta, phi; + + Int_t workItem; + while (workQueue.Pop(workItem)) { + // Generate r1 with 1/r^2 distribution + r1 = c1 * c1 * c2 / (c2 * c1 - rng.Uniform(0, 1) * c1 * (c2 - c1)); + cosTheta = 1.0 - r1; + sinTheta = TMath::Sqrt(1.0 - cosTheta * cosTheta); + phi = rng.Uniform(0, kTwoPi); + + v3.SetXYZ(sinTheta * TMath::Cos(phi), sinTheta * TMath::Sin(phi), cosTheta); + + // Update constraint for new direction + gen.ClearConstraints(); + gen.AddConstraint(v3, 1.0, kModePoint); + + if (gen.Generate() == 0) continue; + + for (Int_t i = 0; i < gen.GetNSolutions(); i++) { + weight = gen.GetWeight(i); + pElectronOut = gen.GetDecay(i, 0); + pProtonOut = gen.GetDecay(i, 1); + pMomTransfer = pElectronIn - pElectronOut; + ROOT::Math::PxPyPzMVector pU = pTarget - pElectronOut; + qSquared = -pMomTransfer.M2(); + formGE = 1.0 / ((1.0 + qSquared / kDipoleMassSq) * + (1.0 + qSquared / kDipoleMassSq)); + formGM = kProtonMagneticMoment * formGE; + tau = qSquared / (4.0 * massProton * massProton); + lambda = (pTotal2.M2() - pU.M2()) / (4.0 * massProton * massProton); + + ampSquared = 16.0 * kPi * kPi * alphaQED * alphaQED / (tau * tau) * + ((formGE * formGE + tau * formGM * formGM) / (1.0 + tau) * + (lambda * lambda - tau * tau - tau) + + 2.0 * tau * tau * formGM * formGM); + + weight *= ampSquared; + + Double_t reweight = r1 * r1 * (c2 - c1) / c1 / c2; + hXsec->Fill(TMath::Cos(pElectronOut.Theta()), reweight * weight); + hNorm->Fill(TMath::Cos(pElectronOut.Theta()), reweight); + hCount->Fill(TMath::Cos(pElectronOut.Theta()), 1.0); + } + } +} + +//////////////////////////////////////////////////////////////////////////////// +// Main Demonstration Function - Parallelized +//////////////////////////////////////////////////////////////////////////////// + +void runTFANG() +{ + using namespace FANG; + + // Enable ROOT thread safety - CRITICAL + ROOT::EnableThreadSafety(); + + gStyle->SetOptStat(0); + + Int_t nThreads = std::thread::hardware_concurrency(); + std::cout << "Using std::thread with " << nThreads << " threads" << std::endl; + + Int_t nEvents = 0; + + //========================================================================== + // Setup for 5-body decay test + //========================================================================== + const Int_t kNBody = 5; + Double_t masses[kNBody] = {1.0, 1.0, 1.0, 1.0, 1.0}; + ROOT::Math::PxPyPzMVector pTotal(0, 0, 5, 12); + + Double_t weight = 0.0; + Double_t sumW = 0.0; + Double_t sumW2 = 0.0; + + //========================================================================== + // Test 1: TFANG Full Phase Space Calculation (Parallelized) + //========================================================================== + std::cout << "========================================" << std::endl; + std::cout << "Test 1: Full Phase Space Calculation (TFANG Parallel)" << std::endl; + std::cout << "========================================" << std::endl; + + Int_t nLoop = 1000000; + { + WorkQueue workQueue; + for (Int_t k = 0; k < nLoop; k++) { + workQueue.Push(k); + } + + std::vector results; + std::mutex resultsMutex; + std::vector threads; + + for (Int_t t = 0; t < nThreads; t++) { + threads.emplace_back(WorkerTest1, + t, std::ref(workQueue), kNBody, std::cref(pTotal), masses, + std::ref(results), std::ref(resultsMutex)); + } + + for (auto& t : threads) { + t.join(); + } + + // Aggregate results + sumW = 0.0; + sumW2 = 0.0; + nEvents = 0; + for (const auto& r : results) { + sumW += r.fSumW; + sumW2 += r.fSumW2; + nEvents += r.fNEvents; + } + } + Double_t mean = sumW / nEvents; + Double_t variance = sumW2 / nEvents - mean * mean; + + // Also get phase space using GetPhaseSpace for comparison + Double_t phaseSpace, phaseSpaceErr; + { + TFANG genFull; + genFull.SetDecay(pTotal, kNBody, masses); + genFull.GetPhaseSpace(static_cast(nLoop), phaseSpace, phaseSpaceErr); + } + + std::cout << "nEvents = " << nEvents << std::endl; + std::cout << "Total Phase Space from parallel loop = " << sumW / nEvents + << " +/- " << TMath::Sqrt(variance / nEvents) << std::endl; + std::cout << "Total Phase Space from GetPhaseSpace = " << phaseSpace + << " +/- " << phaseSpaceErr << std::endl; + + //========================================================================== + // Test 2: Partial Phase Space with Detector Constraints (Parallelized) + //========================================================================== + std::cout << "\n========================================" << std::endl; + std::cout << "Test 2: Partial Phase Space (TFANG Parallel)" << std::endl; + std::cout << " - TFANG constrained vs TFANG unconstrained with cuts" << std::endl; + std::cout << "========================================" << std::endl; + + const Int_t kNDet = 3; + Double_t omega[kNDet]; + Double_t shape[kNDet]; + + Double_t detPosX[kNDet - 1] = {0.0, 0.5}; + Double_t detPosY[kNDet - 1] = {0.0, 0.0}; + Double_t detPosZ[kNDet - 1] = {0.5, 0.0}; + Double_t detRadius[kNDet - 1] = {0.2, 0.3}; + + std::vector v3Det; + ROOT::Math::XYZVector v3; + Double_t radius; + Double_t totalOmega = 1.0; + + for (Int_t i = 0; i < kNDet - 1; i++) { + v3.SetXYZ(detPosX[i], detPosY[i], detPosZ[i]); + v3Det.push_back(v3); + radius = TMath::Sqrt(v3.Mag2() + detRadius[i] * detRadius[i]); + omega[i] = kTwoPi * radius * (radius - v3.R()); + shape[i] = 0.0; + totalOmega *= omega[i]; + } + + v3.SetXYZ(0, 0.5, 0); + v3Det.push_back(v3); + omega[2] = 1.2 * kPi; + shape[2] = 0.4; + totalOmega *= omega[2]; + + std::cout << "Detector configurations:" << std::endl; + std::cout << " Det 1: Circle, Omega = " << omega[0] << " sr" << std::endl; + std::cout << " Det 2: Circle, Omega = " << omega[1] << " sr" << std::endl; + std::cout << " Det 3: Strip, Omega = " << omega[2] << " sr" << std::endl; + std::cout << " Total solid angle factor = " << totalOmega << std::endl; + + Double_t totalMass = 0.0; + for (Int_t l = 0; l < kNBody; l++) { + totalMass += masses[l]; + } + + // Create per-thread histograms for TFANG constrained + std::vector hFangE_vec(nThreads * kNBody); + std::vector hFangCos_vec(nThreads * kNBody); + std::vector hFangPhi_vec(nThreads * kNBody); + + for (Int_t t = 0; t < nThreads; t++) { + for (Int_t i = 0; i < kNBody; i++) { + hFangE_vec[t * kNBody + i] = new TH1D(Form("hFangE_%d_%d", t, i), "", + 100, 0, pTotal.E() - totalMass); + hFangCos_vec[t * kNBody + i] = new TH1D(Form("hFangCos_%d_%d", t, i), "", 50, -1, 1); + hFangPhi_vec[t * kNBody + i] = new TH1D(Form("hFangPhi_%d_%d", t, i), "", 50, -kPi, kPi); + hFangE_vec[t * kNBody + i]->SetDirectory(0); + hFangCos_vec[t * kNBody + i]->SetDirectory(0); + hFangPhi_vec[t * kNBody + i]->SetDirectory(0); + } + } + + // Run TFANG constrained (parallel) + nLoop = 1000000; + { + WorkQueue workQueue; + for (Int_t k = 0; k < nLoop; k++) { + workQueue.Push(k); + } + + std::vector results; + std::mutex resultsMutex; + std::vector threads; + + for (Int_t t = 0; t < nThreads; t++) { + threads.emplace_back(WorkerTest2TFANG, + t, std::ref(workQueue), kNBody, kNDet, + std::cref(pTotal), masses, omega, shape, std::cref(v3Det), + totalOmega, + std::ref(hFangE_vec), std::ref(hFangCos_vec), std::ref(hFangPhi_vec), + std::ref(results), std::ref(resultsMutex)); + } + + for (auto& t : threads) { + t.join(); + } + + // Aggregate results + sumW = 0.0; + sumW2 = 0.0; + nEvents = 0; + for (const auto& r : results) { + sumW += r.fSumW; + sumW2 += r.fSumW2; + nEvents += r.fNEvents; + } + } + + // Merge per-thread histograms for TFANG constrained + TH1D *hFangE[kNBody]; + TH1D *hFangCos[kNBody]; + TH1D *hFangPhi[kNBody]; + + for (Int_t i = 0; i < kNBody; i++) { + hFangE[i] = new TH1D(Form("hFangE_%d", i), "", 100, 0, pTotal.E() - totalMass); + hFangE[i]->SetMarkerStyle(20); + hFangE[i]->SetLineColor(6); + hFangE[i]->SetMinimum(0); + hFangE[i]->GetXaxis()->SetTitle(Form("p_{%d} Energy", i + 1)); + hFangE[i]->GetXaxis()->SetTitleSize(0.07); + hFangE[i]->GetXaxis()->SetLabelSize(0.06); + hFangE[i]->GetYaxis()->SetLabelSize(0.05); + hFangE[i]->GetYaxis()->SetTitle("Events"); + hFangE[i]->GetYaxis()->SetTitleSize(0.07); + hFangE[i]->GetYaxis()->SetTitleOffset(0.5); + hFangE[i]->GetXaxis()->SetTitleOffset(0.9); + + hFangCos[i] = new TH1D(Form("hFangCos_%d", i), "", 50, -1, 1); + hFangCos[i]->SetMarkerStyle(20); + hFangCos[i]->SetLineColor(6); + hFangCos[i]->SetMinimum(0); + hFangCos[i]->GetXaxis()->SetTitle(Form("p_{%d} cos(#theta)", i + 1)); + hFangCos[i]->GetXaxis()->SetTitleSize(0.07); + hFangCos[i]->SetTitleOffset(0.7); + hFangCos[i]->GetYaxis()->SetTitle("Events"); + hFangCos[i]->GetYaxis()->SetTitleSize(0.07); + hFangCos[i]->GetYaxis()->SetTitleOffset(0.5); + hFangCos[i]->GetXaxis()->SetLabelSize(0.06); + hFangCos[i]->GetYaxis()->SetLabelSize(0.05); + hFangCos[i]->GetXaxis()->SetTitleOffset(0.9); + + hFangPhi[i] = new TH1D(Form("hFangPhi_%d", i), "", 50, -kPi, kPi); + hFangPhi[i]->SetMarkerStyle(20); + hFangPhi[i]->SetLineColor(6); + hFangPhi[i]->SetMinimum(0); + hFangPhi[i]->GetXaxis()->SetTitle(Form("p_{%d} #varphi", i + 1)); + hFangPhi[i]->GetXaxis()->SetTitleSize(0.07); + hFangPhi[i]->SetTitleOffset(0.7); + hFangPhi[i]->GetYaxis()->SetTitle("Events"); + hFangPhi[i]->GetYaxis()->SetTitleSize(0.07); + hFangPhi[i]->GetYaxis()->SetTitleOffset(0.5); + hFangPhi[i]->GetXaxis()->SetLabelSize(0.06); + hFangPhi[i]->GetYaxis()->SetLabelSize(0.05); + hFangPhi[i]->GetXaxis()->SetTitleOffset(0.9); + + for (Int_t t = 0; t < nThreads; t++) { + hFangE[i]->Add(hFangE_vec[t * kNBody + i]); + hFangCos[i]->Add(hFangCos_vec[t * kNBody + i]); + hFangPhi[i]->Add(hFangPhi_vec[t * kNBody + i]); + delete hFangE_vec[t * kNBody + i]; + delete hFangCos_vec[t * kNBody + i]; + delete hFangPhi_vec[t * kNBody + i]; + } + } + + // Get partial phase space using GetPartialPhaseSpace for comparison + { + TFANG genConstrained; + genConstrained.SetDecay(pTotal, kNBody, masses); + for (Int_t i = 0; i < kNDet; i++) { + genConstrained.AddConstraint(v3Det[i], omega[i], shape[i]); + } + genConstrained.GetPartialPhaseSpace(static_cast(nLoop), phaseSpace, phaseSpaceErr); + } + + mean = sumW / nEvents; + variance = sumW2 / nEvents - mean * mean; + std::cout << "\nTFANG Constrained Results:" << std::endl; + std::cout << " nEvents = " << nEvents << std::endl; + std::cout << " Partial Phase Space from parallel loop = " << totalOmega * sumW / nEvents + << " +/- " << totalOmega * TMath::Sqrt(variance / nEvents) << std::endl; + std::cout << " Partial Phase Space from GetPartialPhaseSpace = " << phaseSpace + << " +/- " << phaseSpaceErr << std::endl; + std::cout << " hFangE[0]->Integral() = " << hFangE[0]->Integral() << std::endl; + + // Draw TFANG results + TCanvas *c1 = new TCanvas("c1", "c1 En", 10, 10, 1800, 1500); + c1->Divide(2, static_cast(TMath::Floor(kNBody / 2.0 + 0.6))); + for (Int_t i = 0; i < kNBody; i++) { + c1->cd(i + 1); + gPad->SetBottomMargin(0.15); + hFangE[i]->Draw("hist"); + } + + TCanvas *c2 = new TCanvas("c2", "c2 cos", 10, 10, 1800, 1500); + c2->Divide(2, static_cast(TMath::Floor(kNBody / 2.0 + 0.6))); + for (Int_t i = 0; i < kNBody; i++) { + c2->cd(i + 1); + gPad->SetBottomMargin(0.15); + hFangCos[i]->Draw("hist"); + } + + TCanvas *c3 = new TCanvas("c3", "c3 phi", 10, 10, 1800, 1500); + c3->Divide(2, static_cast(TMath::Floor(kNBody / 2.0 + 0.6))); + for (Int_t i = 0; i < kNBody; i++) { + c3->cd(i + 1); + gPad->SetBottomMargin(0.15); + hFangPhi[i]->Draw("hist"); + } + + //========================================================================== + // TFANG Unconstrained with Cuts Comparison (Parallelized) + //========================================================================== + std::cout << "\n--- TFANG Unconstrained with Cuts (Parallel) ---" << std::endl; + + // Direction vectors for cut comparison + std::vector tv3(kNDet); + for (Int_t i = 0; i < kNDet; i++) { + tv3[i].SetXYZ(v3Det[i].X(), v3Det[i].Y(), v3Det[i].Z()); + tv3[i] = tv3[i].Unit(); + } + + Double_t scaleFactor = 100.0; + + // Create per-thread histograms for TFANG unconstrained with cuts + std::vector hFangCutsE_vec(nThreads * kNBody); + std::vector hFangCutsCos_vec(nThreads * kNBody); + std::vector hFangCutsPhi_vec(nThreads * kNBody); + + for (Int_t t = 0; t < nThreads; t++) { + for (Int_t i = 0; i < kNBody; i++) { + hFangCutsE_vec[t * kNBody + i] = new TH1D(Form("hFangCutsE_%d_%d", t, i), "", + 100, 0, pTotal.E() - totalMass); + hFangCutsCos_vec[t * kNBody + i] = new TH1D(Form("hFangCutsCos_%d_%d", t, i), "", 50, -1, 1); + hFangCutsPhi_vec[t * kNBody + i] = new TH1D(Form("hFangCutsPhi_%d_%d", t, i), "", 50, -kPi, kPi); + hFangCutsE_vec[t * kNBody + i]->SetDirectory(0); + hFangCutsCos_vec[t * kNBody + i]->SetDirectory(0); + hFangCutsPhi_vec[t * kNBody + i]->SetDirectory(0); + } + } + + // Run TFANG unconstrained with cuts (parallel) + { + WorkQueue workQueue; + for (Int_t k = 0; k < static_cast(nLoop * scaleFactor); k++) { + workQueue.Push(k); + } + + std::vector results; + std::mutex resultsMutex; + std::vector threads; + + for (Int_t t = 0; t < nThreads; t++) { + threads.emplace_back(WorkerTest2TFANGCuts, + t, std::ref(workQueue), kNBody, kNDet, + std::cref(pTotal), masses, omega, shape, std::cref(tv3), + scaleFactor, + std::ref(hFangCutsE_vec), std::ref(hFangCutsCos_vec), std::ref(hFangCutsPhi_vec), + std::ref(results), std::ref(resultsMutex)); + } + + for (auto& t : threads) { + t.join(); + } + + // Aggregate results + Int_t nTotalGenerated = 0; + Int_t nPassedCuts = 0; + for (const auto& r : results) { + nTotalGenerated += r.fNTotalGenerated; + nPassedCuts += r.fNPassedCuts; + } + + std::cout << " Total events generated: " << nTotalGenerated << std::endl; + std::cout << " Events passing cuts: " << nPassedCuts << std::endl; + if (nTotalGenerated > 0) { + std::cout << " Cut efficiency: " << 100.0 * nPassedCuts / nTotalGenerated << "%" << std::endl; + } + } + + // Merge per-thread histograms for TFANG unconstrained with cuts + TH1D *hFangCutsE[kNBody]; + TH1D *hFangCutsCos[kNBody]; + TH1D *hFangCutsPhi[kNBody]; + + for (Int_t i = 0; i < kNBody; i++) { + hFangCutsE[i] = new TH1D(Form("hFangCutsE_%d", i), "", 100, 0, pTotal.E() - totalMass); + hFangCutsCos[i] = new TH1D(Form("hFangCutsCos_%d", i), "", 50, -1, 1); + hFangCutsPhi[i] = new TH1D(Form("hFangCutsPhi_%d", i), "", 50, -kPi, kPi); + hFangCutsE[i]->SetMarkerStyle(21); + hFangCutsE[i]->SetMarkerColor(kBlue); + hFangCutsCos[i]->SetMarkerStyle(21); + hFangCutsCos[i]->SetMarkerColor(kBlue); + hFangCutsPhi[i]->SetMarkerStyle(21); + hFangCutsPhi[i]->SetMarkerColor(kBlue); + + for (Int_t t = 0; t < nThreads; t++) { + hFangCutsE[i]->Add(hFangCutsE_vec[t * kNBody + i]); + hFangCutsCos[i]->Add(hFangCutsCos_vec[t * kNBody + i]); + hFangCutsPhi[i]->Add(hFangCutsPhi_vec[t * kNBody + i]); + delete hFangCutsE_vec[t * kNBody + i]; + delete hFangCutsCos_vec[t * kNBody + i]; + delete hFangCutsPhi_vec[t * kNBody + i]; + } + } + + std::cout << " hFangCutsE[0]->Integral() = " << hFangCutsE[0]->Integral() << std::endl; + + // Setup legends + TLegend *leg[3 * kNBody]; + for (Int_t i = 0; i < kNBody * 3; i++) { + leg[i] = new TLegend(0.52, 0.62, 0.85, 0.88); + } + + // Adjust legend positions for some plots + leg[10] = new TLegend(0.12, 0.12, 0.45, 0.38); + leg[11] = new TLegend(0.56, 0.62, 0.89, 0.88); + leg[12] = new TLegend(0.12, 0.62, 0.45, 0.88); + for (Int_t i = 5; i <= 9; i++) { + leg[i] = new TLegend(0.12, 0.62, 0.45, 0.88); + } + + for (Int_t i = 0; i < kNBody; i++) { + leg[i]->AddEntry(hFangE[i], "TFANG constrained", "l"); + leg[i]->AddEntry(hFangCutsE[i], "TFANG unconstrained with cuts", "p"); + + leg[i + kNBody]->AddEntry(hFangCos[i], "TFANG constrained", "l"); + leg[i + kNBody]->AddEntry(hFangCutsCos[i], "TFANG unconstrained with cuts", "p"); + + leg[i + 2 * kNBody]->AddEntry(hFangPhi[i], "TFANG constrained", "l"); + leg[i + 2 * kNBody]->AddEntry(hFangCutsPhi[i], "TFANG unconstrained with cuts", "p"); + } + + // Overlay comparison results + for (Int_t i = 0; i < kNBody; i++) { + c1->cd(i + 1); + hFangCutsE[i]->DrawCopy("ep same"); + leg[i]->Draw(); + + c2->cd(i + 1); + hFangCutsCos[i]->DrawCopy("ep same"); + leg[i + kNBody]->Draw(); + + c3->cd(i + 1); + hFangCutsPhi[i]->DrawCopy("ep same"); + leg[i + 2 * kNBody]->Draw(); + } + + //========================================================================== + // Test 3: Elastic ep Scattering Cross Section (Parallelized) + //========================================================================== + std::cout << "\n========================================" << std::endl; + std::cout << "Test 3: Elastic ep Differential Cross Section (TFANG Parallel)" << std::endl; + std::cout << "========================================" << std::endl; + + const Int_t kNBody2 = 2; + nLoop = 100000; + + Double_t massElectron = 0.000511; + Double_t massProton = 0.938272029; + + ROOT::Math::PxPyPzMVector pTarget(0.0, 0.0, 0.0, massProton); + Double_t kineticE = 3.0; + Double_t gamma = kineticE / massElectron + 1.0; + Double_t beta = TMath::Sqrt(1.0 - 1.0 / (gamma * gamma)); + ROOT::Math::PxPyPzMVector pBeam(0.0, 0.0, gamma * beta * massElectron, massElectron); + ROOT::Math::PxPyPzMVector pTotal2 = pBeam + pTarget; + + ROOT::Math::PxPyPzMVector pElectronIn = pBeam; + Double_t flux = 1.0 / (16.0 * kPi * kPi * + TMath::Sqrt(pElectronIn.Dot(pTarget) * pElectronIn.Dot(pTarget) - + massElectron * massElectron * massProton * massProton)); + + TF1 *fRosenbluth = new TF1("fRosenbluth", fElastic, -1, 0.9992, 1); + Double_t parElastic[1] = {kineticE}; + fRosenbluth->SetParameters(parElastic); + + //========================================================================== + // TFANG Point Generation: Differential Cross Section at Specific Angles + //========================================================================== + Double_t sigmaArr[11]; + Double_t sigmaErrArr[11]; + Double_t cosThetaArr[11]; + Double_t cosThetaErrArr[11]; + + std::cout << "\nCalculating differential cross section at specific angles:" << std::endl; + + { + WorkQueue workQueue; + for (Int_t l = 0; l < 11; l++) { + workQueue.Push(l); + } + + std::vector pointResults; + std::mutex resultsMutex; + std::vector threads; + + for (Int_t t = 0; t < nThreads; t++) { + threads.emplace_back(WorkerTest3Point, + t, std::ref(workQueue), nLoop, + massElectron, massProton, + std::cref(pTotal2), std::cref(pTarget), std::cref(pElectronIn), + flux, std::ref(pointResults), std::ref(resultsMutex)); + } + + for (auto& t : threads) { + t.join(); + } + + // Sort results by cosTheta and extract + std::sort(pointResults.begin(), pointResults.end(), + [](const PointResult& a, const PointResult& b) { + return a.fCosTheta < b.fCosTheta; + }); + + for (Int_t l = 0; l < 11; l++) { + cosThetaArr[l] = pointResults[l].fCosTheta; + cosThetaErrArr[l] = 0.0; + sigmaArr[l] = pointResults[l].fSigma; + sigmaErrArr[l] = pointResults[l].fSigmaErr; + + std::cout << " cos(theta) = " << cosThetaArr[l] + << ": dsigma/dOmega = " << sigmaArr[l] << " +/- " << sigmaErrArr[l] + << " (TFANG/Rosenbluth = " << sigmaArr[l] / fRosenbluth->Eval(cosThetaArr[l]) << ")" + << std::endl; + } + } + + TGraphErrors *grElastic = new TGraphErrors(11, cosThetaArr, sigmaArr, cosThetaErrArr, sigmaErrArr); + grElastic->SetMarkerStyle(20); + grElastic->SetMarkerSize(1.3); + + //========================================================================== + // TFANG Event Generation: Full Angular Distribution (Parallelized) + //========================================================================== + std::cout << "\nGenerating full angular distribution..." << std::endl; + + struct Range_t { + Double_t fC1; + Double_t fC2; + }; + Range_t ranges[4] = {{1.0, 2.0}, {0.4, 1.0}, {0.12, 0.4}, {0.01, 0.12}}; + + // Create per-thread histograms for each range + std::vector hXsec_vec(nThreads); + std::vector hNorm_vec(nThreads); + std::vector hCount_vec(nThreads); + + for (Int_t t = 0; t < nThreads; t++) { + hXsec_vec[t] = new TH1D(Form("hXsec_%d", t), "", 440, -1.1, 1.1); + hNorm_vec[t] = new TH1D(Form("hNorm_%d", t), "", 440, -1.1, 1.1); + hCount_vec[t] = new TH1D(Form("hCount_%d", t), "", 440, -1.1, 1.1); + hXsec_vec[t]->SetDirectory(0); + hNorm_vec[t]->SetDirectory(0); + hCount_vec[t]->SetDirectory(0); + } + + for (Int_t rangeIdx = 0; rangeIdx < 4; rangeIdx++) { + Double_t c1 = ranges[rangeIdx].fC1; + Double_t c2 = ranges[rangeIdx].fC2; + + std::cout << " Range " << rangeIdx + 1 << ": cos(theta) in [" + << 1.0 - c2 << ", " << 1.0 - c1 << "]" << std::endl; + + WorkQueue workQueue; + Int_t nEventsPerRange = 1000000; + for (Int_t k = 0; k < nEventsPerRange; k++) { + workQueue.Push(k); + } + + std::vector threads; + for (Int_t t = 0; t < nThreads; t++) { + threads.emplace_back(WorkerTest3Angular, + t, std::ref(workQueue), nEventsPerRange, + massElectron, massProton, + std::cref(pTotal2), std::cref(pTarget), std::cref(pElectronIn), + c1, c2, + hXsec_vec[t], hNorm_vec[t], hCount_vec[t]); + } + + for (auto& t : threads) { + t.join(); + } + } + + // Merge per-thread histograms + TH1D *hXsec = new TH1D("hXsec", "hXsec", 440, -1.1, 1.1); + TH1D *hNorm = new TH1D("hNorm", "hNorm", 440, -1.1, 1.1); + TH1D *hCount = new TH1D("hCount", "hCount", 440, -1.1, 1.1); + TH1D *hError = new TH1D("hError", "hError", 440, -1.1, 1.1); + hXsec->SetMinimum(1E-17); + hNorm->SetMinimum(1E-17); + hError->SetMinimum(0.999); + + for (Int_t t = 0; t < nThreads; t++) { + hXsec->Add(hXsec_vec[t]); + hNorm->Add(hNorm_vec[t]); + hCount->Add(hCount_vec[t]); + delete hXsec_vec[t]; + delete hNorm_vec[t]; + delete hCount_vec[t]; + } + + // Scale and compute errors + Double_t scaleXsec = flux * 2.0 / hXsec->GetBinWidth(2) / hNorm->Integral(); + + for (Int_t l = 1; l <= hXsec->GetNbinsX(); l++) { + Double_t signal = hXsec->GetBinContent(l); + Double_t error = hXsec->GetBinError(l); + Double_t entries = hCount->GetBinContent(l); + + if (entries > 0) { + hError->SetBinContent(l, signal / error / TMath::Sqrt(entries)); + } + hXsec->SetBinContent(l, signal * scaleXsec); + hXsec->SetBinError(l, error * scaleXsec); + } + + hXsec->SetYTitle("#frac{d#sigma}{d#Omega}(ep -> ep) [GeV^{-2}]"); + hXsec->SetXTitle("cos(#theta_{LAB})"); + hXsec->SetTitle("Electron Energy E=3 GeV"); + + TLegend *legFinal = new TLegend(0.12, 0.68, 0.42, 0.88); + + TCanvas *cFinal = new TCanvas("cFinal", "cFinal his", 10, 10, 1800, 1500); + gPad->SetLogy(); + TH1F *vFrame = gPad->DrawFrame(-1.2, 5E-10, 1, 5E-3); + hXsec->Draw("hist same"); + grElastic->Draw("P"); + fRosenbluth->Draw("same"); + + legFinal->AddEntry(hXsec, "TFANG event generation", "l"); + legFinal->AddEntry(grElastic, "TFANG point calculation", "p"); + legFinal->AddEntry(fRosenbluth, "Rosenbluth cross section", "l"); + legFinal->Draw(); + + vFrame->SetYTitle("#frac{d#sigma}{d#Omega}(ep -> ep) [GeV^{-2}]"); + vFrame->SetXTitle("cos(#theta_{LAB})"); + vFrame->SetTitle("Electron Energy E=3 GeV"); + + TCanvas *cDiag = new TCanvas("cDiag", "cDiag Wi error", 10, 10, 1800, 1500); + cDiag->Divide(2, 1); + cDiag->cd(1); + hNorm->Draw("hist"); + cDiag->cd(2); + hCount->Draw("hist"); + + std::cout << "\n========================================" << std::endl; + std::cout << "runTFANG() completed successfully (parallel version)" << std::endl; + std::cout << "J. High Energ. Phys. 2025, 137 (2025)" << std::endl; + std::cout << "https://doi.org/10.1007/JHEP12(2025)137" << std::endl; + std::cout << "========================================" << std::endl; +}