added QEC2p and QEC2p_err branches

.gitignore

@@ -8,8 +8,9 @@ __pycache__
 
 # exectuable output
 dataprinter
-expand-nuchart-QBxn-QBxp-QBa
 expand-nuchart-QEC
+expand-nuchart-QECp
+expand-nuchart-QEC2p
 
 # graphics output
 *.pdf

Makefile

@@ -5,7 +5,7 @@ LDLIBS = $$(root-config --glibs)
 DATATABLES = ame16/mass16.txt ame16/rct1-16.txt ame16/rct2-16.txt nubase16/nubase2016.txt
 EXAMPLE_FIGURES = qec-all.pdf qec-no-estimates.pdf
 EXAMPLE_DATAFILES = qec.dat qbn.dat
-EXECUTABLES = dataprinter expand-nuchart-QEC expand-nuchart-QBxn-QBxp-QBa
+EXECUTABLES = dataprinter expand-nuchart-QEC expand-nuchart-QECp expand-nuchart-QEC2p
 
 
 
@@ -15,25 +15,26 @@ default: nuchart.root
 
 all: default $(EXAMPLE_FIGURES) qbn.dat
 
-nuchart.root: treemaker.py expand-nuchart-QEC expand-nuchart-QBxn-QBxp-QBa
+nuchart.root: treemaker.py expand-nuchart-QEC expand-nuchart-QECp expand-nuchart-QEC2p
 ifeq ($(ROOTSYS), ) # source thisroot.sh if ROOTSYS is not defined
 	source $$(root-config --bindir)/thisroot.sh && python3 treemaker.py
 else
 	python3 treemaker.py
 endif
 	./expand-nuchart-QEC
-	./expand-nuchart-QBxn-QBxp-QBa
+	./expand-nuchart-QECp
+	./expand-nuchart-QEC2p
 
 qec-all.pdf qec-no-estimates.pdf: graph-example-qec.py
 	python3 $<
 
 graph-example-qec.py: qec.dat
 
-qec.dat: dataprinter Makefile
+qec.dat: dataprinter nuchart.root Makefile
 	./$< nuchart.root a A Z QEC QEC_est "A <= 70" "QEC >= 0" > $@
 
-qbn.dat: dataprinter Makefile
-	./$< nuchart.root a A Z QB QEC Sp QBp "A <= 70" > $@
+qbn.dat: dataprinter nuchart.root Makefile
+	./$< nuchart.root a A Z QB QEC Sp QECp "A <= 70" > $@
 
 dataprinter: dataprinter.cxx
 	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
@@ -41,7 +42,10 @@ dataprinter: dataprinter.cxx
 expand-nuchart-QEC: expand-nuchart-QEC.cxx
 	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
 
-expand-nuchart-QBxn-QBxp-QBa: expand-nuchart-QBxn-QBxp-QBa.cxx
+expand-nuchart-QECp: expand-nuchart-QECp.cxx
+	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
+
+expand-nuchart-QEC2p: expand-nuchart-QEC2p.cxx
 	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
 
 treemaker.py: datagetter.py

expand-nuchart-QEC2p.cxx

+/*
+ * QEC2p = QEC of parent - S2p of daughter
+ *
+ * Please refer to 'expand-nuchart-QEC.cxx' and 'expand-nuchartQECp.cxx' for
+ * an explanation.
+ */
+
+#include <TFile.h>
+#include <TTree.h>
+#include <TMath.h>
+
+using namespace std;
+using namespace TMath;
+
+/*
+ * As a new thing, compared to 'expand-nuchart-QEC.cxx', we define our parameters and structs
+ * here to be able to access them globally.
+ */
+Double_t QECP = 0.0, S2pD = 0.0;
+Double_t QECP_err = 0.0, S2pD_err = 0.0;
+Double_t QEC2p = 0.0;
+Double_t QEC2p_err = 0.0;
+Bool_t QEC2p_est = 0, QEC2p_calc = 0;
+Bool_t QEC2p_err_est = 0, QEC2p_err_calc = 0;
+struct QEC2p_params { Double_t QEC2p; Bool_t QEC2p_est; Bool_t QEC2p_calc; };
+QEC2p_params qec2p = QEC2p_params();
+struct QEC2p_err_params { Double_t QEC2p_err; Bool_t QEC2p_err_est; Bool_t QEC2p_err_calc; };
+QEC2p_err_params qec2p_err = QEC2p_err_params();
+
+void QEC2p_NAN() {
+  QEC2p = NAN;
+  QEC2p_est = 0;
+  QEC2p_calc = 0;
+  QEC2p_err = NAN;
+  QEC2p_est = 0;
+  QEC2p_calc = 0;
+
+  qec2p.QEC2p = QEC2p;
+  qec2p.QEC2p_est = QEC2p_est;
+  qec2p.QEC2p_calc = QEC2p_calc;
+  qec2p_err.QEC2p_err = QEC2p_err;
+  qec2p_err.QEC2p_err_est = QEC2p_err_est;
+  qec2p_err.QEC2p_err_calc = QEC2p_err_calc;
+}
+
+int main() {
+  TFile f("nuchart.root", "update");
+  TTree *a;
+  f.GetObject("a", a);
+
+  struct QEC_params { Double_t QEC; Bool_t QEC_est; Bool_t QEC_calc; };
+  QEC_params QEC = QEC_params();
+  a->SetBranchAddress("QEC", &QEC);
+  struct QEC_err_params { Double_t QEC_err; Bool_t QEC_err_est; Bool_t QEC_err_calc; };
+  QEC_err_params QEC_err = QEC_err_params();
+  a->SetBranchAddress("QEC_err", &QEC_err);
+  struct S2p_params { Double_t S2p; Bool_t S2p_est; Bool_t S2p_calc; };
+  S2p_params S2p = S2p_params();
+  a->SetBranchAddress("S2p", &S2p);
+  struct S2p_err_params { Double_t S2p_err; Bool_t S2p_err_est; Bool_t S2p_err_calc; };
+  S2p_err_params S2p_err = S2p_err_params();
+  a->SetBranchAddress("S2p_err", &S2p_err);
+
+  TBranch *qec2p_b = a->Branch("QEC2p", &qec2p.QEC2p, "QEC2p/D:QEC2p_est/O:QEC2p_calc");
+  TBranch *qec2p_err_b = a->Branch("QEC2p_err", &qec2p_err.QEC2p_err, "QEC2p_err/D:QEC2p_err_est/O:QEC2p_err_calc");
+
+  a->BuildIndex("A", "Z");
+  Long64_t indP = -1, indD = -1;
+  Long64_t minA = a->GetMinimum("A"), maxA = a->GetMaximum("A"), minZ = a->GetMinimum("Z"), maxZ = a->GetMaximum("Z");
+  for (Long64_t A = minA; A < maxA; A++) {
+    for (Long64_t Z = minZ; Z < maxZ; Z++) {
+      indP = a->GetEntryNumberWithIndex(A, Z);
+      indD = a->GetEntryNumberWithIndex(A, Z - 1);
+      
+      if (indP != -1 && indD != -1) {
+        QEC2p_calc = 1; QEC2p_err_calc = 1;
+        a->GetEntry(indP);
+        if (QEC.QEC_est == 1) {
+          QEC2p_est = 1; QEC2p_err_est = 1;
+        }
+        QECP = QEC.QEC;
+        QECP_err = QEC_err.QEC_err;
+        a->GetEntry(indD);
+        if (S2p.S2p_calc == 0) { // this is new compared to 'expand-nuchart-QEC.cxx': Sp might not be defined, so we need to check for this
+          QEC2p_NAN();
+        } else {
+          if (S2p.S2p_est == 1) {
+            QEC2p_est = 1; QEC2p_err_est = 1;
+          } else {
+            QEC2p_est = 0; QEC2p_err_est = 0;
+          }
+          S2pD = S2p.S2p;
+          S2pD_err = S2p_err.S2p_err;
+          QEC2p = QECP - S2pD;
+          QEC2p_err = Sqrt(QECP_err*QECP_err + S2pD_err*S2pD_err);
+
+          qec2p.QEC2p = QEC2p;
+          qec2p.QEC2p_est = QEC2p_est;
+          qec2p.QEC2p_calc = QEC2p_calc;
+          qec2p_err.QEC2p_err = QEC2p_err;
+          qec2p_err.QEC2p_err_est = QEC2p_err_est;
+          qec2p_err.QEC2p_err_calc = QEC2p_err_calc;
+        }
+      } else if (indP != -1) {
+        QEC2p_NAN();
+      } else {
+        continue;
+      }
+      a->GetEntry(indP);
+      qec2p_b->Fill();
+      qec2p_err_b->Fill();
+    }
+  }
+  a->Write("", TObject::kOverwrite);
+  f.Close();
+  
+  return EXIT_SUCCESS;
+}