feat: add dihadron analysis

macro/dihadrons/analysis.C

@@ -0,0 +1,44 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Christopher Dilks
+
+R__LOAD_LIBRARY(EpicAnalysis)
+
+void analysis(
+    TString source="delphes",
+    Double_t eleBeamEn=18, Double_t ionBeamEn=275
+    )
+{
+  Analysis *A;
+  TString configFile, outfilePrefix;
+  if(source=="delphes") {
+    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 {
+    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();
+};

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
 
@@ -107,6 +117,7 @@ Analysis::Analysis(
 
   weightInclusive = new WeightsUniform();
   weightTrack     = new WeightsUniform();
+  weightDihadron  = new WeightsUniform();
   weightJet       = new WeightsUniform();
 
   // miscellaneous
@@ -310,11 +321,11 @@ 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);
-  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
@@ -328,9 +339,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 +407,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 +546,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"]) {
@@ -562,6 +628,7 @@ void Analysis::Prepare() {
   // initialize total weights
   wInclusiveTotal = 0.;
   wTrackTotal     = 0.;
+  wDihadronTotal  = 0.;
   wJetTotal       = 0.;
 };
 
@@ -670,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
@@ -730,11 +799,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 +918,23 @@ void Analysis::FillHistos1h(Double_t wgt) {
   FillHistos(fill_payload);
 };
 
+// fill 2h (dihadron) histograms
+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
 void Analysis::FillHistosJets(Double_t wgt) {
 #ifndef EXCLUDE_DELPHES

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,18 +127,22 @@ 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 FillHistos2h(Double_t wgt);        // dihadron kinematics
     void FillHistosJets(Double_t wgt);      // jet kinematics
 
     // shared objects
     SimpleTree *ST;
     HFSTree *HFST;
     ParticleTree *PT;
     Kinematics *kin, *kinTrue;
+    Dihadron *dih, *dihTrue;
+    DihadronSet *dihSet;
     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 = "--------------------------------------------";
@@ -172,7 +181,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 +211,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`

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;

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,19 @@ 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"]) {
+      auto wDihadron = Q2weightFactor * weightDihadron->GetWeight(*kinTrue);
+      wDihadronTotal += wDihadron;
+      dihSet->CalculateKinematics(this, wDihadron); // calls FillHistos2h(wDihadron) and FillHistosInclusive(wDihadron);
+    }
 
     // jet loop - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     if(includeOutputSet["jets"]) {

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;

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)) continue;
 
       // calculate reconstructed hadron kinematics
       kin->vecHadron = part.vecPart;

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);