added QBa and QBa_err branches

.gitignore

@@ -13,6 +13,7 @@ expand-nuchart-QECp
 expand-nuchart-QEC2p
 expand-nuchart-QECa
 expand-nuchart-QB2n
+expand-nuchart-QBa
 
 # 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 expand-nuchart-QB2n
+EXPANSION_EXS = expand-nuchart-QEC expand-nuchart-QECp expand-nuchart-QEC2p expand-nuchart-QECa expand-nuchart-QB2n expand-nuchart-QBa
 EXS = dataprinter
 
 
@@ -53,6 +53,9 @@ $(EXPANSION_DIR)/expand-nuchart-QECa: $(EXPANSION_DIR)/expand-nuchart-QECa.cxx
 $(EXPANSION_DIR)/expand-nuchart-QB2n: $(EXPANSION_DIR)/expand-nuchart-QB2n.cxx
 	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
 
+$(EXPANSION_DIR)/expand-nuchart-QBa: $(EXPANSION_DIR)/expand-nuchart-QBa.cxx
+	$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
+
 treemaker.py: datagetter.py
 
 datagetter.py: $(DATATABLES)

expansion/expand-nuchart-QBa.cxx

+/*
+ * QBa = QB of parent + Qa 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>
+#include <iostream>
+
+using namespace std;
+using namespace TMath;
+
+Double_t QBP = 0.0, QaD = 0.0;
+Double_t QBP_err = 0.0, QaD_err = 0.0;
+Double_t QBa = 0.0;
+Double_t QBa_err = 0.0;
+Bool_t QBa_est = 0, QBa_calc = 0;
+Bool_t QBa_err_est = 0, QBa_err_calc = 0;
+struct QBa_params { Double_t QBa; Bool_t QBa_est; Bool_t QBa_calc; };
+QBa_params qba = QBa_params();
+struct QBa_err_params { Double_t QBa_err; Bool_t QBa_err_est; Bool_t QBa_err_calc; };
+QBa_err_params qba_err = QBa_err_params();
+
+void QBa_NAN() {
+  QBa = NAN;
+  QBa_est = 0;
+  QBa_calc = 0;
+  QBa_err = NAN;
+  QBa_est = 0;
+  QBa_calc = 0;
+
+  qba.QBa = QBa;
+  qba.QBa_est = QBa_est;
+  qba.QBa_calc = QBa_calc;
+  qba_err.QBa_err = QBa_err;
+  qba_err.QBa_err_est = QBa_err_est;
+  qba_err.QBa_err_calc = QBa_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 Qa_params { Double_t Qa; Bool_t Qa_est; Bool_t Qa_calc; };
+  Qa_params Qa = Qa_params();
+  a->SetBranchAddress("Qa", &Qa);
+  struct Qa_err_params { Double_t Qa_err; Bool_t Qa_err_est; Bool_t Qa_err_calc; };
+  Qa_err_params Qa_err = Qa_err_params();
+  a->SetBranchAddress("Qa_err", &Qa_err);
+
+  TBranch *qba_b = a->Branch("QBa", &qba.QBa, "QBa/D:QBa_est/O:QBa_calc");
+  TBranch *qba_err_b = a->Branch("QBa_err", &qba_err.QBa_err, "QBa_err/D:QBa_err_est/O:QBa_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) {
+        QBa_calc = 1; QBa_err_calc = 1;
+        a->GetEntry(indP);
+        if (QB.QB_calc == 0) {
+          QBa_NAN();
+        } else {
+          if (QB.QB_est == 1) {
+            QBa_est = 1; QBa_err_est = 1;
+          }
+          QBP = QB.QB;
+          QBP_err = QB_err.QB_err;
+          a->GetEntry(indD);
+          if (Qa.Qa_calc == 0 || Qa.Qa == 0.0) { // 4H --(beta)--> 4He. Qa defined as zero (and calculable) in nuchart.root for 4He; gives nonsense result here. In principle this also happens e.g. when finding QEC for the neutron in 'expand-nuchart-QEC.cxx', but AME's convention is to define the "pure decay" Q-values even in such a case. Skik følge, land fly.
+            QBa_NAN();
+          } else {
+            if (Qa.Qa_est == 1) {
+              QBa_est = 1; QBa_err_est = 1;
+            } else {
+              QBa_est = 0; QBa_err_est = 0;
+            }
+            QaD = Qa.Qa;
+            QaD_err = Qa_err.Qa_err;
+            QBa = QBP + QaD;
+            QBa_err = Sqrt(QBP_err*QBP_err + QaD_err*QaD_err);
+
+            qba.QBa = QBa;
+            qba.QBa_est = QBa_est;
+            qba.QBa_calc = QBa_calc;
+            qba_err.QBa_err = QBa_err;
+            qba_err.QBa_err_est = QBa_err_est;
+            qba_err.QBa_err_calc = QBa_err_calc;
+          }
+        }
+      } else if (indP != -1) {
+        QBa_NAN();
+      } else {
+        continue;
+      }
+      a->GetEntry(indP);
+      qba_b->Fill();
+      qba_err_b->Fill();
+    }
+  }
+  a->Write("", TObject::kOverwrite);
+  f.Close();
+  
+  return EXIT_SUCCESS;
+}