add separate finalStateID for dihadrons

Author: c-dilks

macro/dihadrons/analysis.C

@@ -0,0 +1,37 @@
+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;
+  A->includeOutputSet["2h"] = true;
+
+  A->AddFinalState("pipTrack");
+  A->AddFinalState("pimTrack");
+
+  // single-hadron 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
+  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

@@ -67,6 +67,9 @@ Analysis::Analysis(
   PIDtoFinalState.insert({ -321,  "KmTrack"  });
   PIDtoFinalState.insert({  2212, "pTrack"   });
 
+  // dihadrons
+  availableBinSchemes.insert({ "dihMh", "M_{h}" });
+
   // kinematics reconstruction methods
   // - choose one of these methods using `SetReconMethod(TString name)`
   // - if you specify none, a default method will be chosen
@@ -101,7 +104,6 @@ Analysis::Analysis(
   infiles.clear();
   entriesTot = 0;
   errorCnt = 0;
-  dihSet = new DihadronSet();
 };
 
 
@@ -320,6 +322,38 @@ void Analysis::Prepare() {
     SetReconMethod("Ele");
   };
 
+  // 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("_$")) = "";
+      // aesthetic quick fix: re-order the name and title for pi+pi-
+      if(dihadronFinalState=="pimTrack_pipTrack") {
+        dihadronFinalState = "pipTrack_pimTrack";
+        dihadronTitle      = "#pi^{+}#pi^{-}";
+      }
+      dihadronTitle += " dihadrons";
+      // add the new dihadron final state
+      finalStateToTitle.insert({dihadronFinalState,dihadronTitle});
+      fmt::print("DEFINE DIHADRON finalStateID='{}'  title='{}'\n",dihadronFinalState,dihadronTitle);
+      AddFinalState(dihadronFinalState);
+      dihSet->SetFinalStateID(dihadronFinalState);
+    }
+  }
+
   // build HistosDAG with specified binning
   HD = new HistosDAG();
   HD->Build(binSchemes);
@@ -347,7 +381,7 @@ void Analysis::Prepare() {
   HD->SetBinSchemeValue("tSpin", [this](){ return (Double_t)kin->tSpin;         });
   HD->SetBinSchemeValue("lSpin", [this](){ return (Double_t)kin->lSpin;         });
   /* dihadron */
-  HD->SetBinSchemeValue("Mh",    [this](){ return dih->Mh;                      });
+  HD->SetBinSchemeValue("dihMh", [this](){ return dih->Mh;                      });
   /* jets */
 #ifndef EXCLUDE_DELPHES
   HD->SetBinSchemeValue("JetPT", [this](){ return kin->pTjet;                   });
@@ -478,7 +512,7 @@ void Analysis::Prepare() {
     }
     // -- dihadron kinematics
     if(includeOutputSet["2h"]) {
-      HS->DefineHist1D("Mh","M_{h}","GeV",NBINS,0,5);
+      HS->DefineHist1D("dihMh", "M_{h}", "GeV", 2*NBINS, 0, 5);
     }
     // -- jet kinematics
 #ifndef EXCLUDE_DELPHES
@@ -690,7 +724,6 @@ void Analysis::AddFinalState(TString finalStateN) {
   };
   BinScheme("finalState")->BuildExternalBin(finalStateN,finalStateT);
   activeFinalStates.insert(finalStateN);
-  if(includeOutputSet["2h"]) dihSet->IncludeHadron(finalStateN);
   fmt::print("AddFinalState: name='{}'\n               title='{}'\n",finalStateN,finalStateT);
 };
 
@@ -795,7 +828,7 @@ void Analysis::FillHistosTracks() {
 void Analysis::FillHistosDihadrons() {
   HD->CheckBins();
   HD->Payload([this](Histos *H){
-    H->FillHist1D("Mh", dih->Mh, wTrack);
+    H->FillHist1D("dihMh", dih->Mh, wTrack);
     });
   HD->ExecuteOps(true);
 }

src/Dihadrons.cxx

@@ -49,10 +49,10 @@ void DihadronSet::CalculateKinematics(Analysis *A) {
   for(std::size_t i=0; i<dihListRec.size(); i++) {
     A->dih     = &(dihListRec[i]);
     A->dihTrue = &(dihListGen[i]);
-    auto finalStateID_tmp = A->finalStateID; // FIXME: this is a bad workaround to fill a single final state
-    A->finalStateID = includedHadrons[0];
+    auto finalStateID_tmp = A->finalStateID; // temporarily change Analysis::finalStateID
+    A->finalStateID = dihadronFinalStateID;  // to that of this DihadronSet
     A->FillHistosDihadrons();
-    A->finalStateID = finalStateID_tmp; // revert bad workaround
+    A->finalStateID = finalStateID_tmp; // revert Analysis::finalStateID
   }
 
   // reset internal storage

src/Dihadrons.h

@@ -17,15 +17,17 @@ class Dihadron {
 
 class DihadronSet {
   public:
-    DihadronSet() : debug(false) {};
+    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);
 };