Combi dev

examples/ProcessHepMCTCSCombi.C

@@ -0,0 +1,125 @@
+#include "CommonDefines.h"
+#include "HepMCElectro.h"
+#include "KinematicsProcElectro.h"
+#include "KineCalculation.h"
+#include "Indicing.h"
+#include "ElectronScatterKinematics.h"
+#include "gammaN_2_Spin0Spin0SpinHalf.h"
+#include "TCSKinematics.h"
+#include <TBenchmark.h>
+
+/**
+ * @brief Example Script: TCS Analysis with Combinatorics and Missing Mass.
+ * Updated to use SetMesonParticles shortcut and CloneLinked.
+ */
+void ProcessHepMCTCSCombi(){
+
+  using namespace rad;
+  using namespace rad::consts::data_type; 
+
+  gBenchmark->Start("df");
+
+  // =================================================================================
+  // 1. SETUP & INJECTION
+  // =================================================================================
+  HepMCElectro hepmc{
+      "hepmc3_tree", 
+      "/w/work5/home/garyp/eic/Farm/data/EpIC_DDVCS_ee_18x275/rootfiles/18x275_ddvcs_edecay_hplus.root"
+  };
+  hepmc.SetupMC();
+
+  // =================================================================================
+  // 2. PARTICLE DEFINITIONS
+  // =================================================================================
+  hepmc.SetBeamIonIndex(3);
+  hepmc.SetBeamElectronIndex(0);
+  //hepmc.SetScatElectronCandidates({1, 6}); 
+  //hepmc.SetParticleCandidates("ele", {1, 6}); 
+  hepmc.SetScatElectronIndex(1);
+  hepmc.SetParticleIndex("ele",6);
+  hepmc.SetParticleIndex("pos", 7); 
+  hepmc.SetParticleIndex("pprime", 5);   
+
+  hepmc.MakeCombinations();
+
+  // =================================================================================
+  // 3. KINEMATICS PROCESSOR (Standard Topology)
+  // =================================================================================
+  KinematicsProcElectro kine{&hepmc, MC()}; 
+
+  // A. Define Particles FIRST (Topology Construction)
+  kine.Creator().Sum("gprime", {{"ele", "pos"}});       
+
+  // Missing Mass: n_calc = (Beam + Target) - (ScatEle + Jpsi + pi+)
+  kine.Creator().Diff("pprime_calc",{
+        {consts::BeamIon(), consts::BeamEle()},  
+	{"gprime", consts::ScatEle()}
+    });
+
+  // B. Define Groups NEXT (Lazy Definition)
+  // This uses the Processor's SetGroup wrapper, ensuring "Z" exists before the group is built.
+  kine.SetMesonParticles({"ele","pos"});
+  kine.SetBaryonParticles({"pprime"});
+
+  // =================================================================================
+  // 4. CALCULATIONS (Registered on Master)
+  // =================================================================================
+  // These will be calculated for Master AND automatically copied to the Linked clone
+
+  kine.Q2();         
+  kine.xbj();
+  kine.y();         
+  kine.nu();
+  kine.tau();
+  kine.tauprime();
+
+  kine.RegisterCalc("GammaPol",rad::physics::ElS_PolVirtPhot);
+  kine.RegisterCalc("GammaPolCirc",rad::physics::ElS_CircPolVirtPhot);
+
+  kine.RegisterCalc("CosThetaHel",rad::gn2s0s0s12::CosThetaHel);
+  kine.RegisterCalc("ThetaHel",rad::gn2s0s0s12::ThetaHel);
+  kine.RegisterCalc("PhiHel",rad::gn2s0s0s12::PhiHel);
+
+  kine.CosThetaCM(); 
+  kine.PhiCM();       
+
+  kine.Mass("GMass", {"gprime"});
+  kine.Mass2("Qp2", {"gprime"});
+  kine.Mass2("MissMass2", 
+	     {consts::BeamIon(), consts::BeamEle()}, 
+	     {"gprime", consts::ScatEle(), "pprime"});
+
+  kine.RegisterCalc("t_top", rad::physics::TTop);
+  kine.RegisterCalc("t_bot", rad::physics::TBot);
+  kine.RegisterCalc("tp_top", rad::physics::TPrimeTop);
+  kine.RegisterCalc("tp_bot", rad::physics::TPrimeBot);
+
+  kine.RegisterCalc("deltaT_top", rad::physics::DeltaTTop);
+  kine.RegisterCalc("deltaT_bot", rad::physics::DeltaTBot);
+
+  // =================================================================================
+  // 5. LINKED PROCESSOR (Cloned Hypothesis)
+  // =================================================================================
+
+  // 1. Clone: Copies all the calculations registered above (Q2, Phi, Mass, tp...)
+  //    Creates a new processor for the "mc_" stream but with suffix "_ncalc"
+  // auto kine_miss = kine.CloneLinked("_ncalc");
+
+  // // 2. Customize: Override "Baryons" to be "n_calc" (missing neutron) for this hypothesis
+  // kine_miss->SetBaryonParticles({"n_calc"});
+
+  // =================================================================================
+  // 6. INITIALIZATION & EXECUTION
+  // =================================================================================
+
+  // Init Master: Executes definitions and calculations for Standard Topology
+  kine.Init();
+
+  // Init Linked: Binds to 'kine' topology and executes copied calculations for Missing Topology
+  //kine_miss->InitLinked(kine);
+
+  gBenchmark->Start("snapshot");
+  hepmc.Snapshot("/w/work5/home/garyp/combirad_trees/HepMC_ddvcs_ee_18x275_hplus.root");
+  gBenchmark->Stop("snapshot");
+  gBenchmark->Print("snapshot");
+}

include/AnalysisManager.h

@@ -59,7 +59,19 @@ namespace rad {
      * @param fileGlob Input file path/pattern.
      */
     AnalysisManager(const std::string& name, const std::string& treeName, const std::string& fileGlob);
-
+
+
+    /**
+     * @brief Constructor (multiple input files).
+     * @param name Name of the analysis (used for output filenames).
+     * @param treeName Input TTree name.
+     * @param files Vector of input file paths.
+     */
+    AnalysisManager(const std::string& name,
+		    const std::string& treeName,
+		    const std::vector<std::string>& files);
+
+
     /** @brief Sets and creates the output directory. */
     void SetOutputDir(const std::string& dir);
 
@@ -212,8 +224,17 @@ namespace rad {
   // ===========================================================================
 
   template <typename R, typename P>
-  inline AnalysisManager<R,P>::AnalysisManager(const std::string& name, const std::string& treeName, const std::string& fileGlob) 
-      : _reaction(treeName, fileGlob), _name(name) {}
+    inline AnalysisManager<R,P>::AnalysisManager(const std::string& name,
+						 const std::string& treeName,
+						 const std::string& fileGlob) 
+    : _reaction(treeName, fileGlob), _name(name) {}
+
+
+  template <typename R, typename P>
+    inline AnalysisManager<R,P>::AnalysisManager(const std::string& name,
+						 const std::string& treeName,
+						 const std::vector<std::string>& files)
+    : _reaction(treeName, files), _name(name) {}
 
   template <typename R, typename P>
   inline void AnalysisManager<R,P>::SetOutputDir(const std::string& dir) {

include/ElectroIonReaction.h

@@ -201,13 +201,7 @@ namespace rad {
 
         /** @brief Helper to set the internal index for the Beam Ion (Always 0). */
         void SetBeamIonIndex(const int idx, const std::string& type = "");
-
-       /** @brief Helper to set the internal index for the Beam electron via JIT function expression */
-      void SetBeamElectronExpr(const std::string& type, const std::string& expression);
-
-       /** @brief Helper to set the internal index for the Beam Ion via JIT function expression */
-      void SetBeamIonExpr(const std::string& type, const std::string& expression);
-
+
     protected:
 
       PxPyPzMVector _p4el_beam;   ///< Internal storage for Electron Beam P4
@@ -292,7 +286,7 @@ namespace rad {
         _mcBeamEleIdx = eleIdx;
         _mcBeamIonIdx = ionIdx;
 	_useBeamsFromMC = true;
-   }
+    }
 
     inline void ElectroIonReaction::SetBeamElectronColumns(const std::string& px, const std::string& py, 
                                                            const std::string& pz, const std::string& m, 
@@ -389,24 +383,14 @@ namespace rad {
     inline PxPyPzMVector ElectroIonReaction::P4BeamIon() const { return _p4ion_beam; }
     inline PxPyPzMVector ElectroIonReaction::P4BeamEle() const { return _p4el_beam; }
 
-
-  // --- Beam Electron via function---
-  inline void ElectroIonReaction::SetBeamElectronExpr(const std::string& type, const std::string& expression){
-    SetParticleCandidatesExpr(consts::BeamEle(),type,expression);
-  }
-  // --- Beam Ion via function---
-  inline void ElectroIonReaction::SetBeamIonExpr(const std::string& type, const std::string& expression){
-    SetParticleCandidatesExpr(consts::BeamIon(),type,expression);
-  }
+    inline void ElectroIonReaction::SetBeamElectronIndex(const int idx, const std::string& type) {
+        SetParticleIndex(consts::BeamEle().data(), type.empty() ? GetDefaultType() : type, idx);
+    }
 
-  inline void ElectroIonReaction::SetBeamElectronIndex(const int idx, const std::string& type) {
-    SetParticleIndex(consts::BeamEle().data(), type.empty() ? GetDefaultType() : type, idx);
-  }
+    inline void ElectroIonReaction::SetBeamIonIndex(const int idx, const std::string& type) {
+        SetParticleIndex(consts::BeamIon().data(), type.empty() ? GetDefaultType() : type, idx);
+    }
 
-  inline void ElectroIonReaction::SetBeamIonIndex(const int idx, const std::string& type) {
-    SetParticleIndex(consts::BeamIon().data(), type.empty() ? GetDefaultType() : type, idx);
-  }
-
     // --- Physics Helpers ---
 
     namespace electroion {

include/ElectronScatterKinematics.h

@@ -43,7 +43,87 @@ namespace rad {
       auto phot = PhotoFourVector(react, px, py, pz, m);
       return -phot.M2();
     }
-
+
+    /**
+     * @brief Calculates the struck quark momentum fraction, xbj = Q^2/(2p.q).
+     * @param react The fixed ReactionMap.
+     * @param px, py, pz, m The consolidated momentum component vectors.
+     * @return ResultType_t The xbj value (always positive).
+     */
+    inline ResultType_t ElS_xbj(const RVecIndexMap& react, 
+			       const RVecResultType& px, const RVecResultType& py, 
+			       const RVecResultType& pz, const RVecResultType& m) 
+    {
+      auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); 
+      auto phot = PhotoFourVector(react, px, py, pz, m);
+      return -phot.M2() / (2*ion.Dot(phot));
+    }
+
+    /**
+     * @brief Calculates the inelasticity parameter, y = p.q/p.k
+     * @param react The fixed ReactionMap.
+     * @param px, py, pz, m The consolidated momentum component vectors.
+     * @return ResultType_t The y value (always positive).
+     */
+    inline ResultType_t ElS_y(const RVecIndexMap& react, 
+			       const RVecResultType& px, const RVecResultType& py, 
+			       const RVecResultType& pz, const RVecResultType& m) 
+    {
+      auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); 
+      auto ebeam = FourVector(react[consts::OrderBeams()][consts::OrderBeamEle()], px, py, pz, m);
+      auto phot = PhotoFourVector(react, px, py, pz, m);
+      return ion.Dot(phot) / ion.Dot(ebeam);
+    }
+
+    /**
+     * @brief Calculates the virtual photon energy (ion rest frame), nu = p.q / M
+     * @param react The fixed ReactionMap.
+     * @param px, py, pz, m The consolidated momentum component vectors.
+     * @return ResultType_t The nu value (always positive).
+     */
+    inline ResultType_t ElS_nu(const RVecIndexMap& react, 
+			       const RVecResultType& px, const RVecResultType& py, 
+			       const RVecResultType& pz, const RVecResultType& m) 
+    {
+      auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); 
+      auto phot = PhotoFourVector(react, px, py, pz, m);
+      return ion.Dot(phot) / ion.M();
+    }
+
+    /**
+     * @brief Calculates the electroproduction variable, tau  = Q^2/(4M^2).
+     * @param react The fixed ReactionMap.
+     * @param px, py, pz, m The consolidated momentum component vectors.
+     * @return ResultType_t The tau value (always positive).
+     */
+    inline ResultType_t ElS_tau(const RVecIndexMap& react, 
+			       const RVecResultType& px, const RVecResultType& py, 
+			       const RVecResultType& pz, const RVecResultType& m) 
+    {
+      auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); 
+      auto phot = PhotoFourVector(react, px, py, pz, m);
+      //auto Q2 = -phot.M2(); //call ElS_Q2() here or just compute by hand?
+      return -phot.M2() / (4*ion.M2());
+    }
+
+    /**
+     * @brief Calculates the  variable, tau' = Q'^2/(2p.q).
+     * @param react The fixed ReactionMap.
+     * @param px, py, pz, m The consolidated momentum component vectors.
+     * @return ResultType_t The tau' value (always positive).
+     */
+    inline ResultType_t ElS_tauprime(const RVecIndexMap& react, 
+			       const RVecResultType& px, const RVecResultType& py, 
+			       const RVecResultType& pz, const RVecResultType& m) 
+    {
+      // Q^2 is defined as the negative invariant mass squared of the virtual photon.
+      auto ion = FourVector(react[consts::OrderBeams()][consts::OrderBeamIon()], px, py, pz, m); 
+      auto phot = PhotoFourVector(react, px, py, pz, m);
+      auto mes = FourVector(react[consts::OrderMesons()], px, py, pz, m); // Assumes single particle or combined meson vector
+      return mes.M2() / (2*ion.Dot(phot));
+    }
+
+
     /**
      * @brief Calculates the four-momentum of the initial state Center-of-Mass (CM) system.
      * @param react The fixed ReactionMap.
@@ -104,7 +184,7 @@ namespace rad {
 
       return {theta, energy, Q2};
     }
-
+    
     /**
      * @brief Calculates the polarization parameter (epsilon) for the virtual photon.
      * @param react The fixed ReactionMap.
@@ -128,6 +208,20 @@ namespace rad {
 
       return pol;
     }
+
+    /**
+     * @brief Calculates the circular polarization parameter for the virtual photon sqrt(1-epsilon^2).
+     * @param react The fixed ReactionMap.
+     * @param px, py, pz, m The consolidated momentum component vectors.
+     * @return ResultType_t The circular polarization parameter.
+     */
+    inline ResultType_t ElS_CircPolVirtPhot(const RVecIndexMap& react, 
+					const RVecResultType& px, const RVecResultType& py, 
+					const RVecResultType& pz, const RVecResultType& m) 
+    {
+      auto eps = ElS_PolVirtPhot(react,px,py,pz,m);
+      return sqrt(1-eps*eps);
+    }
 
     // --- Decay Frame Kinematics (CM and Proton Rest) ---
 

include/KinematicsProcElectro.h

@@ -37,6 +37,21 @@ namespace rad {
     /** @brief Calculates Q2 (Negative Squared 4-Momentum Transfer). */
     void Q2();
 
+    /** @brief Calculates xbj. */
+    void xbj();
+
+    /** @brief Calculates y. */
+    void y();
+
+    /** @brief Calculates nu. */
+    void nu();
+
+    /** @brief Calculates tau. */
+    void tau();
+
+    /** @brief Calculates tauprime. */
+    void tauprime();
+
     /** @brief Calculates Cos(Theta) in the Center of Mass frame. */
     void CosThetaCM();
 
@@ -78,7 +93,27 @@ namespace rad {
   inline void KinematicsProcElectro::Q2() {
     RegisterCalc("Q2", rad::physics::ElS_Q2);
   }
-
+
+  inline void KinematicsProcElectro::xbj() {
+    RegisterCalc("xbj", rad::physics::ElS_xbj);
+  }
+
+  inline void KinematicsProcElectro::y() {
+    RegisterCalc("y", rad::physics::ElS_y);
+  }
+
+  inline void KinematicsProcElectro::nu() {
+    RegisterCalc("nu", rad::physics::ElS_nu);
+  }
+
+  inline void KinematicsProcElectro::tau() {
+    RegisterCalc("tau", rad::physics::ElS_tau);
+  }
+
+  inline void KinematicsProcElectro::tauprime() {
+    RegisterCalc("tauprime", rad::physics::ElS_tauprime);
+  }
+
   inline void KinematicsProcElectro::CosThetaCM() {
     RegisterCalc("CosThetaCM", rad::physics::ElS_CosThetaCM);
   }