Commit 8a7e55e6 authored by Erik Asbjørn Mikkelsen Jensen's avatar Erik Asbjørn Mikkelsen Jensen
Browse files

added QECa and QECa_err branches

parent 0162b708
......@@ -11,6 +11,7 @@ dataprinter
expand-nuchart-QEC
expand-nuchart-QECp
expand-nuchart-QEC2p
expand-nuchart-QECa
# graphics output
*.pdf
......
......@@ -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
EXPANSION_EXS = expand-nuchart-QEC expand-nuchart-QECp expand-nuchart-QEC2p expand-nuchart-QECa
EXS = dataprinter
......@@ -25,6 +25,7 @@ else
python3 treemaker.py
endif
$(foreach exec, $(addprefix $(EXPANSION_DIR)/, $(EXPANSION_EXS)), ./$(exec); )
$(RM) -r __pycache__
example1: $(EXAMPLE_DIR)/graph-example-qec.py # PHONEY target avoids multiple executions with two real targets ('qec-all.pdf', 'qec-no-estimates.pdf'), but results in 'example1' always being executed if make is run on 'all' or 'example1'
python3 $<
......@@ -46,6 +47,9 @@ $(EXPANSION_DIR)/expand-nuchart-QECp: $(EXPANSION_DIR)/expand-nuchart-QECp.cxx
$(EXPANSION_DIR)/expand-nuchart-QEC2p: $(EXPANSION_DIR)/expand-nuchart-QEC2p.cxx
$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
$(EXPANSION_DIR)/expand-nuchart-QECa: $(EXPANSION_DIR)/expand-nuchart-QECa.cxx
$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
treemaker.py: datagetter.py
datagetter.py: $(DATATABLES)
......
......@@ -11,3 +11,4 @@
* Make generic script which writes desired values from root file to a text file ✔
* Combine subprojects into one big project ✔
* Make QEC example ✔
* Tidied root folder up ✔
\ No newline at end of file
/*
* QECa = QEC 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 QECP = 0.0, QaD = 0.0;
Double_t QECP_err = 0.0, QaD_err = 0.0;
Double_t QECa = 0.0;
Double_t QECa_err = 0.0;
Bool_t QECa_est = 0, QECa_calc = 0;
Bool_t QECa_err_est = 0, QECa_err_calc = 0;
struct QECa_params { Double_t QECa; Bool_t QECa_est; Bool_t QECa_calc; };
QECa_params qeca = QECa_params();
struct QECa_err_params { Double_t QECa_err; Bool_t QECa_err_est; Bool_t QECa_err_calc; };
QECa_err_params qeca_err = QECa_err_params();
void QECa_NAN() {
QECa = NAN;
QECa_est = 0;
QECa_calc = 0;
QECa_err = NAN;
QECa_est = 0;
QECa_calc = 0;
qeca.QECa = QECa;
qeca.QECa_est = QECa_est;
qeca.QECa_calc = QECa_calc;
qeca_err.QECa_err = QECa_err;
qeca_err.QECa_err_est = QECa_err_est;
qeca_err.QECa_err_calc = QECa_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 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 *qeca_b = a->Branch("QECa", &qeca.QECa, "QECa/D:QECa_est/O:QECa_calc");
TBranch *qeca_err_b = a->Branch("QECa_err", &qeca_err.QECa_err, "QECa_err/D:QECa_err_est/O:QECa_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) {
QECa_calc = 1; QECa_err_calc = 1;
a->GetEntry(indP);
if (QEC.QEC_calc == 0) {
QECa_NAN();
} else {
if (QEC.QEC_est == 1) {
QECa_est = 1; QECa_err_est = 1;
}
QECP = QEC.QEC;
QECP_err = QEC_err.QEC_err;
a->GetEntry(indD);
if (Qa.Qa_calc == 0 || Qa.Qa == 0.0) { // 4Li --(EC)--> 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.
QECa_NAN();
} else {
if (Qa.Qa_est == 1) {
QECa_est = 1; QECa_err_est = 1;
} else {
QECa_est = 0; QECa_err_est = 0;
}
QaD = Qa.Qa;
QaD_err = Qa_err.Qa_err;
QECa = QECP + QaD;
QECa_err = Sqrt(QECP_err*QECP_err + QaD_err*QaD_err);
qeca.QECa = QECa;
qeca.QECa_est = QECa_est;
qeca.QECa_calc = QECa_calc;
qeca_err.QECa_err = QECa_err;
qeca_err.QECa_err_est = QECa_err_est;
qeca_err.QECa_err_calc = QECa_err_calc;
}
}
} else if (indP != -1) {
QECa_NAN();
} else {
continue;
}
a->GetEntry(indP);
qeca_b->Fill();
qeca_err_b->Fill();
}
}
a->Write("", TObject::kOverwrite);
f.Close();
return EXIT_SUCCESS;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment