From 89a01a1c8d08cea7fa6d6d3842fe76e23121a01c Mon Sep 17 00:00:00 2001
From: Christopher Dilks <christopher.j.dilks@gmail.com>
Date: Fri, 14 Oct 2022 18:46:23 -0400
Subject: [PATCH 1/3] feat: add dihadron analysis

---
 macro/dihadrons/analysis.C   |  42 +++++++++++++
 src/Analysis.cxx             |  96 +++++++++++++++++++++++++++-
 src/Analysis.h               |  13 +++-
 src/AnalysisAthena.cxx       |   2 +-
 src/AnalysisDelphes.cxx      |   5 +-
 src/AnalysisEcce.cxx         |   2 +-
 src/AnalysisEpic.cxx         |   2 +-
 src/Dihadrons.cxx            | 119 +++++++++++++++++++++++++++++++++++
 src/Dihadrons.h              |  36 +++++++++++
 src/Kinematics.h             |   6 ++
 src/LinkDef.h                |   2 +
 tutorial/analysis_template.C |   1 +
 12 files changed, 318 insertions(+), 8 deletions(-)
 create mode 100644 macro/dihadrons/analysis.C
 create mode 100644 src/Dihadrons.cxx
 create mode 100644 src/Dihadrons.h

diff --git a/macro/dihadrons/analysis.C b/macro/dihadrons/analysis.C
new file mode 100644
index 00000000..41b6b897
--- /dev/null
+++ b/macro/dihadrons/analysis.C
@@ -0,0 +1,42 @@
+R__LOAD_LIBRARY(Sidis-eic)
+
+void analysis(
+    TString source="delphes",
+    Double_t eleBeamEn=5,  Double_t ionBeamEn=41
+    // Double_t eleBeamEn=18, Double_t ionBeamEn=275
+    )
+{
+  Analysis *A;
+  TString configFile, outfilePrefix;
+  if(source=="delphes") {
+    configFile = Form("datarec/delphes/%dx%d/delphes.config",(int)eleBeamEn,(int)ionBeamEn);
+    outfilePrefix = Form("dihadrons.delphes.%dx%d",(int)eleBeamEn,(int)ionBeamEn);
+    A = new AnalysisDelphes(configFile, outfilePrefix);
+  } else {
+    fmt::print(stderr,"ERROR: source '{}' not implemented\n",source);
+    return;
+  }
+
+  A->maxEvents = 30000; // limiter
+  A->writeSimpleTree = true;
+  A->SetReconMethod("Ele");
+  A->includeOutputSet["1h"] = true; // optionally output single-hadron plots
+
+  // dihadron final state ==================================
+  A->includeOutputSet["2h"] = true; // include the output set
+  A->AddFinalState("pipTrack");     // call `AddFinalState` exactly 2 times: once for each hadron
+  A->AddFinalState("pimTrack");
+
+  // cuts ==================================================
+  // - inclusive cuts
+  A->AddBinScheme("w");     A->BinScheme("w")->BuildBin("Min",3.0);         // W > 3 GeV
+  A->AddBinScheme("y");     A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95
+  // - single-hadron cuts:
+  /* don't use these for the dihadron final state, since they will only apply to the second hadron
+  A->AddBinScheme("z");     A->BinScheme("z")->BuildBin("Range",0.2,0.9);   // 0.2 < z < 0.9
+  A->AddBinScheme("xF");    A->BinScheme("xF")->BuildBin("Min",0.0);        // xF > 0
+  A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1);     // pT_lab > 0.1 GeV (tracking limit)
+  */
+
+  A->Execute();
+};
diff --git a/src/Analysis.cxx b/src/Analysis.cxx
index 502e9804..430699b5 100644
--- a/src/Analysis.cxx
+++ b/src/Analysis.cxx
@@ -53,6 +53,15 @@ Analysis::Analysis(
   availableBinSchemes.insert({ "JetEta", "jet eta" });
   availableBinSchemes.insert({ "JetE", "jet energy" });
 #endif
+  /* dihadrons */
+  availableBinSchemes.insert({ "DihMh",     "M_{h}"    });
+  availableBinSchemes.insert({ "DihMX",     "M_{X}"    });
+  availableBinSchemes.insert({ "DihZ",      "Z"        });
+  availableBinSchemes.insert({ "DihPhPerp", "P_{h,T}"  });
+  availableBinSchemes.insert({ "DihTheta",  "#theta"   });
+  availableBinSchemes.insert({ "DihPhiH",   "#phi_{h}" });
+  availableBinSchemes.insert({ "DihPhiR",   "#phi_{R}" });
+  availableBinSchemes.insert({ "DihPhiS",   "#phi_{S}" });
 
   // available final states
   // - specify which final states you want to include using `AddFinalState(TString name)`
@@ -87,6 +96,7 @@ Analysis::Analysis(
   // - the default settings are set here; override them at the macro level
   includeOutputSet.insert({ "inclusive",      true  }); // inclusive kinematics
   includeOutputSet.insert({ "1h",             true  }); // single hadron kinematics
+  includeOutputSet.insert({ "2h",             false }); // dihadron kinematics
   includeOutputSet.insert({ "jets",           false }); // jet kinematics
   includeOutputSet.insert({ "depolarization", false }); // depolarization factors & ratios
 
@@ -310,7 +320,7 @@ void Analysis::Prepare() {
   outFile = new TFile(outfileName,"RECREATE");
 
   // instantiate shared objects
-  kin = new Kinematics(eleBeamEn,ionBeamEn,crossingAngle);
+  kin     = new Kinematics(eleBeamEn, ionBeamEn, crossingAngle);
   kinTrue = new Kinematics(eleBeamEn, ionBeamEn, crossingAngle);
   ST   = new SimpleTree("tree",kin,kinTrue);
   HFST = new HFSTree("hfstree",kin,kinTrue);
@@ -328,9 +338,36 @@ void Analysis::Prepare() {
   includeOutputSet.insert({ "inclusive_only",
       includeOutputSet["inclusive"]
       && !includeOutputSet["1h"]
+      && !includeOutputSet["2h"]
       && !includeOutputSet["jets"]
       });
 
+  // if including dihadrons, define a dihadron final state
+  if(includeOutputSet["2h"]) {
+    if(activeFinalStates.size()!=2) {
+      fmt::print(stderr,"ERROR: cannot include dihadron outputSet, since there should only be 2 final states defined\n");
+      includeOutputSet["2h"] = false;
+    } else {
+      // add to dihSet
+      dihSet = new DihadronSet();
+      // set finalStateID and title
+      TString dihadronFinalState = "";
+      TString dihadronTitle      = "";
+      for(auto state : activeFinalStates) {
+        dihadronFinalState += state + "_";
+        dihadronTitle      += finalStateToTitle.at(state);
+        dihadronTitle(TRegexp(" .*")) = "";
+        dihSet->IncludeHadron(state);
+      }
+      dihadronFinalState(TRegexp("_$")) = "";
+      dihadronTitle += " dihadrons";
+      // add the new dihadron final state
+      finalStateToTitle.insert({dihadronFinalState,dihadronTitle});
+      AddFinalState(dihadronFinalState);
+      dihSet->SetFinalStateID(dihadronFinalState);
+    }
+  }
+
   // if there are no final states defined, default to definitions here:
   if(BinScheme("finalState")->GetNumBins()==0) {
     std::cout << "NOTE: adding pi+ tracks for final state, since you specified none" << std::endl;
@@ -369,6 +406,15 @@ void Analysis::Prepare() {
   HD->SetBinSchemeValue("phiS",  [this](){ return kin->phiS;                    });
   HD->SetBinSchemeValue("tSpin", [this](){ return (Double_t)kin->tSpin;         });
   HD->SetBinSchemeValue("lSpin", [this](){ return (Double_t)kin->lSpin;         });
+  /* dihadron */
+  HD->SetBinSchemeValue("DihMh",     [this](){ return dih->Mh;     });
+  HD->SetBinSchemeValue("DihMX",     [this](){ return dih->MX;     });
+  HD->SetBinSchemeValue("DihZ",      [this](){ return dih->Z;      });
+  HD->SetBinSchemeValue("DihPhPerp", [this](){ return dih->PhPerp; });
+  HD->SetBinSchemeValue("DihTheta",  [this](){ return dih->Theta;  });
+  HD->SetBinSchemeValue("DihPhiH",   [this](){ return dih->PhiH;   });
+  HD->SetBinSchemeValue("DihPhiR",   [this](){ return dih->PhiR;   });
+  HD->SetBinSchemeValue("DihPhiS",   [this](){ return dih->PhiS;   });
   /* jets */
 #ifndef EXCLUDE_DELPHES
   HD->SetBinSchemeValue("JetPT", [this](){ return kin->pTjet;                   });
@@ -499,7 +545,26 @@ void Analysis::Prepare() {
           NBINS,-TMath::Pi(),TMath::Pi()
           );
     }
-
+    // -- dihadron kinematics
+    if(includeOutputSet["2h"]) {
+      HS->DefineHist1D("DihMh",     "M_{h}",    "GeV", 2*NBINS, 0,            5);
+      HS->DefineHist1D("DihMX",     "M_{X}",    "GeV", NBINS,   0,            40);
+      HS->DefineHist1D("DihZ",      "Z",        "",    NBINS,   0,            1);
+      HS->DefineHist1D("DihPhPerp", "P_{h,T}",  "GeV", NBINS,   1e-2,         3, true);
+      HS->DefineHist1D("DihTheta",  "#theta",   "",    NBINS,   0,            TMath::Pi());
+      HS->DefineHist1D("DihPhiH",   "#phi_{h}", "",    NBINS,   -TMath::Pi(), TMath::Pi());
+      HS->DefineHist1D("DihPhiR",   "#phi_{R}", "",    NBINS,   -TMath::Pi(), TMath::Pi());
+      HS->DefineHist1D("DihPhiS",   "#phi_{S}", "",    NBINS,   -TMath::Pi(), TMath::Pi());
+      HS->DefineHist2D("DihPhiHvsPhiR", "#phi_{R}", "#phi_{h}", "", "",
+          NBINS, -TMath::Pi(), TMath::Pi(),
+          NBINS, -TMath::Pi(), TMath::Pi()
+          );
+      HS->DefineHist2D("DihThetaVsPh", "P_{h}", "#theta", "GeV", "",
+          NBINS, 1e-2, 40,
+          NBINS, 0, TMath::Pi(),
+          true
+          );
+    }
     // -- jet kinematics
 #ifndef EXCLUDE_DELPHES
     if(includeOutputSet["jets"]) {
@@ -730,11 +795,18 @@ void Analysis::AddFinalState(TString finalStateN) {
     return;
   };
   BinScheme("finalState")->BuildExternalBin(finalStateN,finalStateT);
-  activeFinalStates.insert(finalStateN);
+  activeFinalStates.push_back(finalStateN);
   fmt::print("AddFinalState: name='{}'\n               title='{}'\n",finalStateN,finalStateT);
 };
 
 
+// check if this final state bin has been added
+//----------------------------------------------
+Bool_t Analysis::IsFinalState(TString finalState) {
+  return std::find( activeFinalStates.begin(), activeFinalStates.end(), finalState) != activeFinalStates.end();
+};
+
+
 // FillHistos methods: check bins and fill associated histograms
 // - checks which bins the track/jet/etc. falls in
 // - fills the histograms in the associated Histos objects
@@ -842,6 +914,24 @@ void Analysis::FillHistos1h(Double_t wgt) {
   FillHistos(fill_payload);
 };
 
+// fill 2h (dihadron) histograms
+void Analysis::FillHistosDihadrons() {
+  HD->CheckBins();
+  HD->Payload([this](Histos *H){
+    H->FillHist1D("DihMh",     dih->Mh,     wTrack);
+    H->FillHist1D("DihMX",     dih->MX,     wTrack);
+    H->FillHist1D("DihZ",      dih->Z,      wTrack);
+    H->FillHist1D("DihPhPerp", dih->PhPerp, wTrack);
+    H->FillHist1D("DihTheta",  dih->Theta,  wTrack);
+    H->FillHist1D("DihPhiH",   dih->PhiH,   wTrack);
+    H->FillHist1D("DihPhiR",   dih->PhiR,   wTrack);
+    H->FillHist1D("DihPhiS",   dih->PhiS,   wTrack);
+    H->FillHist2D("DihPhiHvsPhiR", dih->PhiR, dih->PhiH,  wTrack);
+    H->FillHist2D("DihThetaVsPh",  dih->Ph,   dih->Theta, wTrack);
+    });
+  HD->ExecuteOps(true);
+}
+
 // fill jet histograms
 void Analysis::FillHistosJets(Double_t wgt) {
 #ifndef EXCLUDE_DELPHES
diff --git a/src/Analysis.h b/src/Analysis.h
index e98b1832..f6fc635f 100644
--- a/src/Analysis.h
+++ b/src/Analysis.h
@@ -32,12 +32,16 @@
 #include "Histos.h"
 #include "HistosDAG.h"
 #include "Kinematics.h"
+#include "Dihadrons.h"
 #include "SimpleTree.h"
 #include "HFSTree.h"
 #include "ParticleTree.h"
 #include "Weights.h"
 #include "CommonConstants.h"
 
+class Dihadron;
+class DihadronSet;
+
 class Analysis : public TNamed
 {
   public:
@@ -58,6 +62,7 @@ class Analysis : public TNamed
 
     // add a new final state bin
     void AddFinalState(TString finalStateN);
+    Bool_t IsFinalState(TString finalState); // check if it's been added
 
     // common settings
     Bool_t verbose; // if true, print a lot more information
@@ -122,6 +127,7 @@ class Analysis : public TNamed
     // `FillHistos(weight)` methods: fill histograms
     void FillHistosInclusive(Double_t wgt); // inclusive kinematics
     void FillHistos1h(Double_t wgt);        // single-hadron kinematics
+    void FillHistosDihadrons();
     void FillHistosJets(Double_t wgt);      // jet kinematics
 
     // shared objects
@@ -129,6 +135,8 @@ class Analysis : public TNamed
     HFSTree *HFST;
     ParticleTree *PT;
     Kinematics *kin, *kinTrue;
+    Dihadron *dih, *dihTrue;
+    DihadronSet *dihSet;
     HistosDAG *HD;
     Weights const* weightInclusive;
     Weights const* weightTrack;
@@ -172,7 +180,7 @@ class Analysis : public TNamed
     std::map<TString,TString> reconMethodToTitle;
     std::map<TString,TString> finalStateToTitle;
     std::map<int,TString> PIDtoFinalState;
-    std::set<TString> activeFinalStates;
+    std::vector<TString> activeFinalStates;
 
     // check if Q2 `val` is between `min` and `max`; if `max==0`, only `val>=min` is checked
     template<class T> bool InQ2Range(T val, T min, T max, bool ignoreZero=false) {
@@ -202,6 +210,9 @@ class Analysis : public TNamed
       fmt::print("}}\n");
     }
 
+    // expose protected members to `DihadronSet`, such as `kin` and `kinTrue`
+    friend class DihadronSet;
+
   private:
 
     // fill histograms, according to `fill_payload`
diff --git a/src/AnalysisAthena.cxx b/src/AnalysisAthena.cxx
index c7c1c909..f92de2ec 100644
--- a/src/AnalysisAthena.cxx
+++ b/src/AnalysisAthena.cxx
@@ -233,7 +233,7 @@ void AnalysisAthena::Execute()
       //   histograms; if not, proceed to next track
       auto kv = PIDtoFinalState.find(pid_);
       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
-      if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
+      if(!IsFinalState(finalStateID)) continue;
 
       // calculate reconstructed hadron kinematics
       kin->vecHadron = part.vecPart;
diff --git a/src/AnalysisDelphes.cxx b/src/AnalysisDelphes.cxx
index cfff9534..5e35c7fa 100644
--- a/src/AnalysisDelphes.cxx
+++ b/src/AnalysisDelphes.cxx
@@ -183,7 +183,7 @@ void AnalysisDelphes::Execute() {
           );
       auto kv = PIDtoFinalState.find(pid);
       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
-      if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
+      if(!IsFinalState(finalStateID)) continue;
 
       // get parent particle, to check if pion is from vector meson
       GenParticle *trkParticle = (GenParticle*)trk->Particle.GetObject();
@@ -228,12 +228,15 @@ void AnalysisDelphes::Execute() {
         // - `IsActiveEvent()` is only true if at least one bin gets filled for this track
         if( writeSimpleTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
       }
+      if(includeOutputSet["2h"]) dihSet->AddHadron(this);
 
       // tests
       //kin->ValidateHeadOnFrame();
 
     }; // end track loop
 
+    // dihadrons - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+    if(includeOutputSet["2h"]) dihSet->CalculateKinematics(this);
 
     // jet loop - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     if(includeOutputSet["jets"]) {
diff --git a/src/AnalysisEcce.cxx b/src/AnalysisEcce.cxx
index ae78b325..efff1068 100644
--- a/src/AnalysisEcce.cxx
+++ b/src/AnalysisEcce.cxx
@@ -448,7 +448,7 @@ void AnalysisEcce::Execute()
       //   histograms; if not, proceed to next track
       auto kv = PIDtoFinalState.find(pid_);
       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
-      if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
+      if(!IsFinalState(finalStateID)) continue;
 
       // calculate reconstructed hadron kinematics
       kin->vecHadron = part.vecPart;
diff --git a/src/AnalysisEpic.cxx b/src/AnalysisEpic.cxx
index 47d8af5c..1cb1cbcf 100644
--- a/src/AnalysisEpic.cxx
+++ b/src/AnalysisEpic.cxx
@@ -322,7 +322,7 @@ void AnalysisEpic::Execute()
       //   histograms; if not, proceed to next track
       auto kv = PIDtoFinalState.find(pid_);
       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
-      if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
+      if(!IsFinalState(finalStateID)) return;
 
       // calculate reconstructed hadron kinematics
       kin->vecHadron = part.vecPart;
diff --git a/src/Dihadrons.cxx b/src/Dihadrons.cxx
new file mode 100644
index 00000000..5b84c7df
--- /dev/null
+++ b/src/Dihadrons.cxx
@@ -0,0 +1,119 @@
+#include "Dihadrons.h"
+
+ClassImp(DihadronSet)
+ClassImp(Dihadron)
+
+Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K)
+{
+  const TLorentzVector &vecE   = K->vecEleBeam;
+  const TLorentzVector &vecP   = K->vecIonBeam;
+  const TLorentzVector &vecL   = K->vecElectron;
+  const TLorentzVector &vecQ   = K->vecQ;
+  const TLorentzVector &vecW   = K->vecW;
+  const TLorentzVector &vecPh  = vecHad1 + vecHad2;
+  const TLorentzVector &vecR   = 0.5 * (vecHad1-vecHad2);
+  const TLorentzVector vecH[2] = {vecHad1, vecHad2};
+
+  // boost to ion rest frame
+  TLorentzVector IvecH[2];
+  TLorentzVector IvecL, IvecQ, IvecPh;
+  K->BoostToIonFrame(vecL,IvecL);
+  K->BoostToIonFrame(vecQ,IvecQ);
+  K->BoostToIonFrame(vecPh,IvecPh);
+  for(int h=0; h<2; h++)
+    K->BoostToIonFrame(vecH[h],IvecH[h]);
+
+  // invariant mass, missing mass, Z, PhPerp, Zeta
+  Mh = vecPh.M();
+  MX = TMath::Abs((vecW-vecPh).M());
+  Z  = vecP.Dot(vecPh) / vecP.Dot(vecQ);
+  for(int h=0; h<2; h++)
+    hadZ[h] = vecP.Dot(vecH[h]) / vecP.Dot(vecQ);
+  Ph     = vecPh.Vect().Mag();
+  PhPerp = Kinematics::Reject(IvecPh.Vect(),IvecQ.Vect()).Mag();
+  Zeta   = 2 * vecR.Dot(vecP) / vecPh.Dot(vecP);
+
+  // PhiH
+  PhiH = Kinematics::PlaneAngle(
+      IvecQ.Vect(), IvecL.Vect(),
+      IvecQ.Vect(), IvecPh.Vect()
+      );
+
+  // PhiR
+  Double_t coeff =
+    K->x * (Zeta*Mh*Mh - (vecH[0].M2() - vecH[1].M2()) ) /
+    ( K->Q2 * Z );
+  TLorentzVector vecRT = vecR - (Zeta/2)*vecPh + coeff*vecP;
+  TLorentzVector IvecRT;
+  K->BoostToIonFrame(vecRT,IvecRT);
+  PhiR = Kinematics::PlaneAngle(
+      IvecQ.Vect(), IvecL.Vect(),
+      IvecQ.Vect(), IvecRT.Vect()
+      );
+
+  // PhiS
+  PhiS = K->phiS;
+
+  // Theta
+  Double_t MRterm[2];
+  for(int h=0; h<2; h++)
+    MRterm[h] = TMath::Sqrt( vecH[h].M2() + vecR.Vect().Mag2() );
+  Theta = TMath::ACos(
+      ( MRterm[0] - MRterm[1] - Mh*Zeta ) / 
+      ( 2 * vecR.Vect().Mag() )  
+      );
+}
+
+////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////
+
+void DihadronSet::IncludeHadron(TString hadName) {
+  includedHadrons.push_back(hadName);
+}
+
+void DihadronSet::AddHadron(Analysis *A) {
+  if(hadSetRec.find(A->finalStateID)==hadSetRec.end()) {
+    hadSetRec.insert({A->finalStateID,{}});
+    hadSetGen.insert({A->finalStateID,{}});
+  }
+  if(debug) fmt::print("DihadronSet: AddHadron '{}'\n",A->finalStateID);
+  hadSetRec[A->finalStateID].push_back( A->kin->vecHadron     );
+  hadSetGen[A->finalStateID].push_back( A->kinTrue->vecHadron );
+}
+
+void DihadronSet::CalculateKinematics(Analysis *A) {
+
+  // TODO: this does not yet handle generalized dihadrons
+  if(includedHadrons.size()!=2) fmt::print("ERROR: more or less than 2 final states defined for DihadronSet\n");
+
+  // hadron pairing
+  auto PairHadrons = [this] (auto hadSet, auto &dihList, Kinematics *K) {
+    TString hadNames[2] = { includedHadrons[0], includedHadrons[1] };
+    for(const auto &had0vec : hadSet[hadNames[0]]) {
+      for(const auto &had1vec : hadSet[hadNames[1]]) {
+        dihList.push_back(Dihadron(had0vec, had1vec, K));
+        if(debug) fmt::print(" pair: {}, {}\n",hadNames[0],hadNames[1]);
+      }
+    }
+  };
+  if(debug) fmt::print("DihadronSet::CalculateKinematics reconstructed\n");
+  PairHadrons(hadSetRec,dihListRec,A->kin);
+  if(debug) fmt::print("DihadronSet::CalculateKinematics generated\n");
+  PairHadrons(hadSetGen,dihListGen,A->kinTrue);
+
+  // fill output data structures
+  for(std::size_t i=0; i<dihListRec.size(); i++) {
+    A->dih     = &(dihListRec[i]);
+    A->dihTrue = &(dihListGen[i]);
+    auto finalStateID_tmp = A->finalStateID; // temporarily change Analysis::finalStateID
+    A->finalStateID = dihadronFinalStateID;  // to that of this DihadronSet
+    A->FillHistosDihadrons();
+    A->finalStateID = finalStateID_tmp; // revert Analysis::finalStateID
+  }
+
+  // reset internal storage
+  hadSetRec.clear();
+  hadSetGen.clear();
+  dihListRec.clear();
+  dihListGen.clear();
+}
diff --git a/src/Dihadrons.h b/src/Dihadrons.h
new file mode 100644
index 00000000..c9ca6822
--- /dev/null
+++ b/src/Dihadrons.h
@@ -0,0 +1,36 @@
+#pragma once
+
+#include <vector>
+#include "Analysis.h"
+#include "Kinematics.h"
+
+class Analysis;
+
+class Dihadron {
+  public:
+    Dihadron() {};
+    Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K);
+    ~Dihadron() {};
+    Double_t Mh, MX, Z, Ph, PhPerp, Theta, PhiH, PhiR, PhiS, Zeta;
+    Double_t hadZ[2];
+  private:
+  ClassDef(Dihadron,1);
+};
+
+class DihadronSet {
+  public:
+    DihadronSet() : debug(false), dihadronFinalStateID("") {};
+    ~DihadronSet() {};
+    void IncludeHadron(TString hadName);
+    void AddHadron(Analysis *A);
+    void SetFinalStateID(TString state) { dihadronFinalStateID = state; }
+    void CalculateKinematics(Analysis *A);
+  private:
+    std::vector<TString> includedHadrons;
+    std::map<TString,std::vector<TLorentzVector>> hadSetRec, hadSetGen;
+    std::vector<Dihadron> dihListRec, dihListGen;
+    TString dihadronFinalStateID;
+    bool debug;
+  ClassDef(DihadronSet,1);
+};
+
diff --git a/src/Kinematics.h b/src/Kinematics.h
index aee5e953..a5b80d6f 100644
--- a/src/Kinematics.h
+++ b/src/Kinematics.h
@@ -204,6 +204,12 @@ class Kinematics : public TObject
     // - convert energy,mass to momentum
     static Double_t EMtoP(Double_t energy, Double_t mass) {
       return TMath::Sqrt( TMath::Power(energy,2) - TMath::Power(mass,2) );
+    }; 
+    // - get angle between two vectors
+    static Double_t AngleSubtend(TVector3 vA, TVector3 vB) {
+      Double_t m = vA.Mag() * vB.Mag();
+      if(m>0) return TMath::ACos( vA.Dot(vB) / m );
+      return -10000;
     };
     // - vector projection: returns vA projected onto vB
     static TVector3 Project(TVector3 vA, TVector3 vB) {
diff --git a/src/LinkDef.h b/src/LinkDef.h
index e9c023ff..49d23908 100644
--- a/src/LinkDef.h
+++ b/src/LinkDef.h
@@ -16,6 +16,8 @@
 
 // analysis objects
 #pragma link C++ class Kinematics+;
+#pragma link C++ class Dihadron+;
+#pragma link C++ class DihadronSet+;
 #pragma link C++ class SimpleTree+;
 #pragma link C++ class HFSTree+;
 #pragma link C++ class ParticleTree+;
diff --git a/tutorial/analysis_template.C b/tutorial/analysis_template.C
index 35f82185..bcb9e740 100644
--- a/tutorial/analysis_template.C
+++ b/tutorial/analysis_template.C
@@ -32,6 +32,7 @@ void analysis_template(
   A->includeOutputSet["jets"] = true;
   // - additional example settings; see `src/Analysis.cxx` for more
   //A->includeOutputSet["1h"] = false;
+  //A->includeOutputSet["2h"] = true;
   //A->includeOutputSet["inclusive"] = false;
   //A->includeOutputSet["depolarization"] = true;
 

From 17deebed2693533457598626c0382fbfc1e3fc77 Mon Sep 17 00:00:00 2001
From: Christopher Dilks <christopher.j.dilks@gmail.com>
Date: Mon, 9 Jan 2023 11:54:44 -0500
Subject: [PATCH 2/3] fix: updates w.r.t. `main`

---
 macro/dihadrons/analysis.C |  5 ++++-
 src/Analysis.cxx           | 35 +++++++++++++++++------------------
 src/Analysis.h             |  2 +-
 src/AnalysisDelphes.cxx    |  3 ++-
 src/AnalysisEpicPodio.cxx  |  2 +-
 src/Dihadrons.cxx          |  7 +++++--
 src/Dihadrons.h            |  5 ++++-
 7 files changed, 34 insertions(+), 25 deletions(-)

diff --git a/macro/dihadrons/analysis.C b/macro/dihadrons/analysis.C
index 41b6b897..e514cfdc 100644
--- a/macro/dihadrons/analysis.C
+++ b/macro/dihadrons/analysis.C
@@ -1,4 +1,7 @@
-R__LOAD_LIBRARY(Sidis-eic)
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Christopher Dilks
+
+R__LOAD_LIBRARY(EpicAnalysis)
 
 void analysis(
     TString source="delphes",
diff --git a/src/Analysis.cxx b/src/Analysis.cxx
index 430699b5..4fecc8e1 100644
--- a/src/Analysis.cxx
+++ b/src/Analysis.cxx
@@ -322,9 +322,9 @@ void Analysis::Prepare() {
   // instantiate shared objects
   kin     = new Kinematics(eleBeamEn, ionBeamEn, crossingAngle);
   kinTrue = new Kinematics(eleBeamEn, ionBeamEn, crossingAngle);
-  ST   = new SimpleTree("tree",kin,kinTrue);
-  HFST = new HFSTree("hfstree",kin,kinTrue);
-  PT   = new ParticleTree("ptree");
+  ST      = new SimpleTree("tree",kin,kinTrue);
+  HFST    = new HFSTree("hfstree",kin,kinTrue);
+  PT      = new ParticleTree("ptree");
 
   // if including jets, define a `jet` final state
 #ifndef EXCLUDE_DELPHES
@@ -915,21 +915,20 @@ void Analysis::FillHistos1h(Double_t wgt) {
 };
 
 // fill 2h (dihadron) histograms
-void Analysis::FillHistosDihadrons() {
-  HD->CheckBins();
-  HD->Payload([this](Histos *H){
-    H->FillHist1D("DihMh",     dih->Mh,     wTrack);
-    H->FillHist1D("DihMX",     dih->MX,     wTrack);
-    H->FillHist1D("DihZ",      dih->Z,      wTrack);
-    H->FillHist1D("DihPhPerp", dih->PhPerp, wTrack);
-    H->FillHist1D("DihTheta",  dih->Theta,  wTrack);
-    H->FillHist1D("DihPhiH",   dih->PhiH,   wTrack);
-    H->FillHist1D("DihPhiR",   dih->PhiR,   wTrack);
-    H->FillHist1D("DihPhiS",   dih->PhiS,   wTrack);
-    H->FillHist2D("DihPhiHvsPhiR", dih->PhiR, dih->PhiH,  wTrack);
-    H->FillHist2D("DihThetaVsPh",  dih->Ph,   dih->Theta, wTrack);
-    });
-  HD->ExecuteOps(true);
+void Analysis::FillHistos2h(Double_t wgt) {
+  auto fill_payload = [this,wgt] (Histos *H) {
+    H->FillHist1D("DihMh",     dih->Mh,     wgt);
+    H->FillHist1D("DihMX",     dih->MX,     wgt);
+    H->FillHist1D("DihZ",      dih->Z,      wgt);
+    H->FillHist1D("DihPhPerp", dih->PhPerp, wgt);
+    H->FillHist1D("DihTheta",  dih->Theta,  wgt);
+    H->FillHist1D("DihPhiH",   dih->PhiH,   wgt);
+    H->FillHist1D("DihPhiR",   dih->PhiR,   wgt);
+    H->FillHist1D("DihPhiS",   dih->PhiS,   wgt);
+    H->FillHist2D("DihPhiHvsPhiR", dih->PhiR, dih->PhiH,  wgt);
+    H->FillHist2D("DihThetaVsPh",  dih->Ph,   dih->Theta, wgt);
+  };
+  FillHistos(fill_payload);
 }
 
 // fill jet histograms
diff --git a/src/Analysis.h b/src/Analysis.h
index f6fc635f..2696af49 100644
--- a/src/Analysis.h
+++ b/src/Analysis.h
@@ -127,7 +127,7 @@ class Analysis : public TNamed
     // `FillHistos(weight)` methods: fill histograms
     void FillHistosInclusive(Double_t wgt); // inclusive kinematics
     void FillHistos1h(Double_t wgt);        // single-hadron kinematics
-    void FillHistosDihadrons();
+    void FillHistos2h(Double_t wgt);        // dihadron kinematics
     void FillHistosJets(Double_t wgt);      // jet kinematics
 
     // shared objects
diff --git a/src/AnalysisDelphes.cxx b/src/AnalysisDelphes.cxx
index 5e35c7fa..fdbc3f67 100644
--- a/src/AnalysisDelphes.cxx
+++ b/src/AnalysisDelphes.cxx
@@ -236,7 +236,8 @@ void AnalysisDelphes::Execute() {
     }; // end track loop
 
     // dihadrons - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-    if(includeOutputSet["2h"]) dihSet->CalculateKinematics(this);
+    if(includeOutputSet["2h"])
+      dihSet->CalculateKinematics(this, Q2weightFactor * weightTrack->GetWeight(*kinTrue)); // FIXME: need separate `Weights` object?
 
     // jet loop - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     if(includeOutputSet["jets"]) {
diff --git a/src/AnalysisEpicPodio.cxx b/src/AnalysisEpicPodio.cxx
index 245fb185..1bcbdd0a 100644
--- a/src/AnalysisEpicPodio.cxx
+++ b/src/AnalysisEpicPodio.cxx
@@ -241,7 +241,7 @@ void AnalysisEpicPodio::Execute()
       // - check PID, to see if it's a final state we're interested in
       auto kv = PIDtoFinalState.find(recPDG);
       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else return;
-      if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) return;
+      if(!IsFinalState(finalStateID)) return;
 
       // set SIDIS particle 4-momenta, and calculate their kinematics
       kinTrue->vecHadron = GetP4(simPart);
diff --git a/src/Dihadrons.cxx b/src/Dihadrons.cxx
index 5b84c7df..c839ddb8 100644
--- a/src/Dihadrons.cxx
+++ b/src/Dihadrons.cxx
@@ -1,3 +1,6 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Christopher Dilks
+
 #include "Dihadrons.h"
 
 ClassImp(DihadronSet)
@@ -81,7 +84,7 @@ void DihadronSet::AddHadron(Analysis *A) {
   hadSetGen[A->finalStateID].push_back( A->kinTrue->vecHadron );
 }
 
-void DihadronSet::CalculateKinematics(Analysis *A) {
+void DihadronSet::CalculateKinematics(Analysis *A, Double_t wgt) {
 
   // TODO: this does not yet handle generalized dihadrons
   if(includedHadrons.size()!=2) fmt::print("ERROR: more or less than 2 final states defined for DihadronSet\n");
@@ -107,7 +110,7 @@ void DihadronSet::CalculateKinematics(Analysis *A) {
     A->dihTrue = &(dihListGen[i]);
     auto finalStateID_tmp = A->finalStateID; // temporarily change Analysis::finalStateID
     A->finalStateID = dihadronFinalStateID;  // to that of this DihadronSet
-    A->FillHistosDihadrons();
+    A->FillHistos2h(wgt);
     A->finalStateID = finalStateID_tmp; // revert Analysis::finalStateID
   }
 
diff --git a/src/Dihadrons.h b/src/Dihadrons.h
index c9ca6822..fb483975 100644
--- a/src/Dihadrons.h
+++ b/src/Dihadrons.h
@@ -1,3 +1,6 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Christopher Dilks
+
 #pragma once
 
 #include <vector>
@@ -24,7 +27,7 @@ class DihadronSet {
     void IncludeHadron(TString hadName);
     void AddHadron(Analysis *A);
     void SetFinalStateID(TString state) { dihadronFinalStateID = state; }
-    void CalculateKinematics(Analysis *A);
+    void CalculateKinematics(Analysis *A, Double_t wgt=1.0);
   private:
     std::vector<TString> includedHadrons;
     std::map<TString,std::vector<TLorentzVector>> hadSetRec, hadSetGen;

From 2a4babde40bd5433cff02dbecdf392ecbf9c3e5e Mon Sep 17 00:00:00 2001
From: Christopher Dilks <christopher.j.dilks@gmail.com>
Date: Mon, 9 Jan 2023 14:51:08 -0500
Subject: [PATCH 3/3] fix: separate weight for dihadrons

---
 macro/dihadrons/analysis.C | 7 +++----
 src/Analysis.cxx           | 4 ++++
 src/Analysis.h             | 3 ++-
 src/AnalysisDelphes.cxx    | 7 +++++--
 src/AnalysisEpic.cxx       | 2 +-
 src/Dihadrons.cxx          | 1 +
 6 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/macro/dihadrons/analysis.C b/macro/dihadrons/analysis.C
index e514cfdc..d1d16139 100644
--- a/macro/dihadrons/analysis.C
+++ b/macro/dihadrons/analysis.C
@@ -5,14 +5,13 @@ R__LOAD_LIBRARY(EpicAnalysis)
 
 void analysis(
     TString source="delphes",
-    Double_t eleBeamEn=5,  Double_t ionBeamEn=41
-    // Double_t eleBeamEn=18, Double_t ionBeamEn=275
+    Double_t eleBeamEn=18, Double_t ionBeamEn=275
     )
 {
   Analysis *A;
   TString configFile, outfilePrefix;
   if(source=="delphes") {
-    configFile = Form("datarec/delphes/%dx%d/delphes.config",(int)eleBeamEn,(int)ionBeamEn);
+    configFile = Form("datarec/hepmc.pythia6/radcor/%dx%d/files.config",(int)eleBeamEn,(int)ionBeamEn);
     outfilePrefix = Form("dihadrons.delphes.%dx%d",(int)eleBeamEn,(int)ionBeamEn);
     A = new AnalysisDelphes(configFile, outfilePrefix);
   } else {
@@ -20,7 +19,7 @@ void analysis(
     return;
   }
 
-  A->maxEvents = 30000; // limiter
+  // A->maxEvents = 30000; // limiter
   A->writeSimpleTree = true;
   A->SetReconMethod("Ele");
   A->includeOutputSet["1h"] = true; // optionally output single-hadron plots
diff --git a/src/Analysis.cxx b/src/Analysis.cxx
index 4fecc8e1..091f42e4 100644
--- a/src/Analysis.cxx
+++ b/src/Analysis.cxx
@@ -117,6 +117,7 @@ Analysis::Analysis(
 
   weightInclusive = new WeightsUniform();
   weightTrack     = new WeightsUniform();
+  weightDihadron  = new WeightsUniform();
   weightJet       = new WeightsUniform();
 
   // miscellaneous
@@ -627,6 +628,7 @@ void Analysis::Prepare() {
   // initialize total weights
   wInclusiveTotal = 0.;
   wTrackTotal     = 0.;
+  wDihadronTotal  = 0.;
   wJetTotal       = 0.;
 };
 
@@ -735,10 +737,12 @@ void Analysis::Finish() {
   HD->Payload([this](Histos *H){ H->Write(); }); HD->ExecuteAndClearOps();
   std::vector<Double_t> vec_wInclusiveTotal { wInclusiveTotal };
   std::vector<Double_t> vec_wTrackTotal     { wTrackTotal     };
+  std::vector<Double_t> vec_wDihadronTotal  { wDihadronTotal  };
   std::vector<Double_t> vec_wJetTotal       { wJetTotal       };
   outFile->WriteObject(&Q2xsecs, "XsTotal");
   outFile->WriteObject(&vec_wInclusiveTotal, "WeightInclusiveTotal");
   outFile->WriteObject(&vec_wTrackTotal,     "WeightTrackTotal");
+  outFile->WriteObject(&vec_wDihadronTotal,  "WeightDihadronTotal");
   outFile->WriteObject(&vec_wJetTotal,       "WeightJetTotal");
 
   // write binning schemes
diff --git a/src/Analysis.h b/src/Analysis.h
index 2696af49..4ab6aec4 100644
--- a/src/Analysis.h
+++ b/src/Analysis.h
@@ -140,8 +140,9 @@ class Analysis : public TNamed
     HistosDAG *HD;
     Weights const* weightInclusive;
     Weights const* weightTrack;
+    Weights const* weightDihadron;
     Weights const* weightJet;
-    Double_t wInclusiveTotal, wTrackTotal, wJetTotal;
+    Double_t wInclusiveTotal, wTrackTotal, wDihadronTotal, wJetTotal;
     Long64_t entriesTot;
     Long64_t errorCnt;
     const TString sep = "--------------------------------------------";
diff --git a/src/AnalysisDelphes.cxx b/src/AnalysisDelphes.cxx
index fdbc3f67..afae7017 100644
--- a/src/AnalysisDelphes.cxx
+++ b/src/AnalysisDelphes.cxx
@@ -236,8 +236,11 @@ void AnalysisDelphes::Execute() {
     }; // end track loop
 
     // dihadrons - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-    if(includeOutputSet["2h"])
-      dihSet->CalculateKinematics(this, Q2weightFactor * weightTrack->GetWeight(*kinTrue)); // FIXME: need separate `Weights` object?
+    if(includeOutputSet["2h"]) {
+      auto wDihadron = Q2weightFactor * weightDihadron->GetWeight(*kinTrue);
+      wDihadronTotal += wDihadron;
+      dihSet->CalculateKinematics(this, wDihadron); // calls FillHistos2h(wDihadron) and FillHistosInclusive(wDihadron);
+    }
 
     // jet loop - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     if(includeOutputSet["jets"]) {
diff --git a/src/AnalysisEpic.cxx b/src/AnalysisEpic.cxx
index 1cb1cbcf..8341214d 100644
--- a/src/AnalysisEpic.cxx
+++ b/src/AnalysisEpic.cxx
@@ -322,7 +322,7 @@ void AnalysisEpic::Execute()
       //   histograms; if not, proceed to next track
       auto kv = PIDtoFinalState.find(pid_);
       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
-      if(!IsFinalState(finalStateID)) return;
+      if(!IsFinalState(finalStateID)) continue;
 
       // calculate reconstructed hadron kinematics
       kin->vecHadron = part.vecPart;
diff --git a/src/Dihadrons.cxx b/src/Dihadrons.cxx
index c839ddb8..3417759b 100644
--- a/src/Dihadrons.cxx
+++ b/src/Dihadrons.cxx
@@ -111,6 +111,7 @@ void DihadronSet::CalculateKinematics(Analysis *A, Double_t wgt) {
     auto finalStateID_tmp = A->finalStateID; // temporarily change Analysis::finalStateID
     A->finalStateID = dihadronFinalStateID;  // to that of this DihadronSet
     A->FillHistos2h(wgt);
+    A->FillHistosInclusive(wgt);
     A->finalStateID = finalStateID_tmp; // revert Analysis::finalStateID
   }