added QB2n and QB2n_err branches

.gitignore

@@ -12,6 +12,7 @@ expand-nuchart-QEC
 expand-nuchart-QECp
 expand-nuchart-QEC2p
 expand-nuchart-QECa
+expand-nuchart-QB2n
 
 # graphics output
 *.pdf

Makefile

@@ -7,7 +7,7 @@ EXAMPLE_DIR = examples
 EXAMPLE_FIGURES = qec-all.pdf qec-no-estimates.pdf
 EXAMPLE_DATAFILES = qec.dat
 EXPANSION_DIR = expansion
-EXPANSION_EXS = expand-nuchart-QEC expand-nuchart-QECp expand-nuchart-QEC2p expand-nuchart-QECa
+EXPANSION_EXS = expand-nuchart-QEC expand-nuchart-QECp expand-nuchart-QEC2p expand-nuchart-QECa expand-nuchart-QB2n
 EXS = dataprinter
 
 
@@ -50,6 +50,9 @@ $(EXPANSION_DIR)/expand-nuchart-QEC2p: $(EXPANSION_DIR)/expand-nuchart-QEC2p.cxx
 $(EXPANSION_DIR)/expand-nuchart-QECa: $(EXPANSION_DIR)/expand-nuchart-QECa.cxx
 	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
 
+$(EXPANSION_DIR)/expand-nuchart-QB2n: $(EXPANSION_DIR)/expand-nuchart-QB2n.cxx
+	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
+
 treemaker.py: datagetter.py
 
 datagetter.py: $(DATATABLES)

expansion/expand-nuchart-QB2n.cxx

+/*
+ * QB2n = QB of parent - S2n 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;
+
+Double_t QBP = 0.0, S2nD = 0.0;
+Double_t QBP_err = 0.0, S2nD_err = 0.0;
+Double_t QB2n = 0.0;
+Double_t QB2n_err = 0.0;
+Bool_t QB2n_est = 0, QB2n_calc = 0;
+Bool_t QB2n_err_est = 0, QB2n_err_calc = 0;
+struct QB2n_params { Double_t QB2n; Bool_t QB2n_est; Bool_t QB2n_calc; };
+QB2n_params qb2n = QB2n_params();
+struct QB2n_err_params { Double_t QB2n_err; Bool_t QB2n_err_est; Bool_t QB2n_err_calc; };
+QB2n_err_params qb2n_err = QB2n_err_params();
+
+void QB2n_NAN() {
+  QB2n = NAN;
+  QB2n_est = 0;
+  QB2n_calc = 0;
+  QB2n_err = NAN;
+  QB2n_est = 0;
+  QB2n_calc = 0;
+
+  qb2n.QB2n = QB2n;
+  qb2n.QB2n_est = QB2n_est;
+  qb2n.QB2n_calc = QB2n_calc;
+  qb2n_err.QB2n_err = QB2n_err;
+  qb2n_err.QB2n_err_est = QB2n_err_est;
+  qb2n_err.QB2n_err_calc = QB2n_err_calc;
+}
+
+int main() {
+  TFile f("nuchart.root", "update");
+  TTree *a;
+  f.GetObject("a", a);
+
+  struct QB_params { Double_t QB; Bool_t QB_est; Bool_t QB_calc; };
+  QB_params QB = QB_params();
+  a->SetBranchAddress("QB", &QB);
+  struct QB_err_params { Double_t QB_err; Bool_t QB_err_est; Bool_t QB_err_calc; };
+  QB_err_params QB_err = QB_err_params();
+  a->SetBranchAddress("QB_err", &QB_err);
+  struct S2n_params { Double_t S2n; Bool_t S2n_est; Bool_t S2n_calc; };
+  S2n_params S2n = S2n_params();
+  a->SetBranchAddress("S2n", &S2n);
+  struct S2n_err_params { Double_t S2n_err; Bool_t S2n_err_est; Bool_t S2n_err_calc; };
+  S2n_err_params S2n_err = S2n_err_params();
+  a->SetBranchAddress("S2n_err", &S2n_err);
+
+  TBranch *qb2n_b = a->Branch("QB2n", &qb2n.QB2n, "QB2n/D:QB2n_est/O:QB2n_calc");
+  TBranch *qb2n_err_b = a->Branch("QB2n_err", &qb2n_err.QB2n_err, "QB2n_err/D:QB2n_err_est/O:QB2n_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) {
+        QB2n_calc = 1; QB2n_err_calc = 1;
+        a->GetEntry(indP);
+        if (QB.QB_calc == 0) {
+          QB2n_NAN();
+        } else {
+          if (QB.QB_est == 1) {
+            QB2n_est = 1; QB2n_err_est = 1;
+          }
+          QBP = QB.QB;
+          QBP_err = QB_err.QB_err;
+          a->GetEntry(indD);
+          if (S2n.S2n_calc == 0) {
+            QB2n_NAN();
+          } else {
+            if (S2n.S2n_est == 1) {
+              QB2n_est = 1; QB2n_err_est = 1;
+            } else {
+              QB2n_est = 0; QB2n_err_est = 0;
+            }
+            S2nD = S2n.S2n;
+            S2nD_err = S2n_err.S2n_err;
+            QB2n = QBP - S2nD;
+            QB2n_err = Sqrt(QBP_err*QBP_err + S2nD_err*S2nD_err);
+
+            qb2n.QB2n = QB2n;
+            qb2n.QB2n_est = QB2n_est;
+            qb2n.QB2n_calc = QB2n_calc;
+            qb2n_err.QB2n_err = QB2n_err;
+            qb2n_err.QB2n_err_est = QB2n_err_est;
+            qb2n_err.QB2n_err_calc = QB2n_err_calc;
+          }
+        }
+      } else if (indP != -1) {
+        QB2n_NAN();
+      } else {
+        continue;
+      }
+      a->GetEntry(indP);
+      qb2n_b->Fill();
+      qb2n_err_b->Fill();
+    }
+  }
+  a->Write("", TObject::kOverwrite);
+  f.Close();
+  
+  return EXIT_SUCCESS;
+}