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

added QEC_err branch

parent 3fd73b20
......@@ -2,7 +2,7 @@
* Make example with N:Z:E colour plot (see Uffe Bergmann) using matshow or pcolor or something else in matplotlib
* Make QBa:A figure
* Do TODO's in top of *datagetter.py*, *expand-nuchart-QBxn-QBxp-QBa.cxx*, *expand-nuchart-QEC.cxx*, *graph-example-qec.py*
* Do TODO's in top of *datagetter.py*, *expand-nuchart-QBxn-QBxp-QBa.cxx*, *graph-example-qec.py*
# DONE
......
// TODO: Add an uncertainty branch, QEC_err, perhaps simply by the error propagation law
/*
* In this ROOT macro the ROOT Tree 'nuchart.root' generated by treemaker.py is expanded by
* adding Q-values for electron capture to the tree via the formula
......@@ -15,8 +13,10 @@
#include <TFile.h>
#include <TTree.h>
#include <TMath.h>
using namespace std;
using namespace TMath;
int main() {
// Open tree file and retrieve the tree named 'a'
......@@ -28,17 +28,28 @@ int main() {
struct ME_params { Double_t ME; Bool_t ME_est; };
ME_params ME = ME_params();
a->SetBranchAddress("ME", &ME);
// Same for errors on mass excesses
struct ME_err_params {Double_t ME_err; Bool_t ME_err_est; };
ME_err_params ME_err = ME_err_params();
a->SetBranchAddress("ME_err", &ME_err);
// Add a new branch to the tree which is to contain QEC values and the above-mentioned boolean values
struct QEC_params { Double_t QEC; Bool_t QEC_est; Bool_t QEC_calc; };
QEC_params qec = QEC_params();
TBranch *newBranch = a->Branch("QEC", &qec.QEC, "QEC/D:QEC_est/O:QEC_calc");
// Same for errors on QEC values
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();
TBranch *newErrorBranch = a->Branch("QEC_err", &qec_err.QEC_err, "QEC_err/D:QEC_err_est/O:QEC_err_calc");
a->BuildIndex("A", "Z"); // prerequisite for searching in the tree based on A and Z values in GetEntryNumberWithIndex() below
Double_t MEP = 0.0, MED = 0.0;
Double_t MEP_err = 0.0, MED_err = 0.0;
Double_t QEC = 0.0;
Double_t QEC_err = 0.0;
Bool_t QEC_est = 0, QEC_calc = 0;
Bool_t QEC_err_est = 0, QEC_err_calc = 0;
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++) {
......@@ -48,36 +59,53 @@ int main() {
if (indP != -1 && indD != -1) { // we can get the mass excesses of both the Parent and the Daughter nuclei
QEC_calc = 1;
QEC_err_calc = 1;
a->GetEntry(indP);
if (ME.ME_est == 1) {
QEC_est = 1;
QEC_err_est = 1;
}
MEP = ME.ME;
MEP_err = ME_err.ME_err;
a->GetEntry(indD);
if (ME.ME_est == 1) {
QEC_est = 1;
QEC_err_est = 1;
} else { // if neither of the ME's were estimates (ME_est = 0) then QEC is not an estimate (QEC_est = 0)
QEC_est = 0;
QEC_err_est = 0;
}
MED = ME.ME;
MED_err = ME_err.ME_err;
QEC = MEP - MED;
QEC_err = Sqrt(MEP_err*MEP_err + MED_err*MED_err); // error propagation law; probably not the best, but hopefully good enough
qec.QEC = QEC; // values might not be recorded correctly, unless the values are assigned in this order
qec.QEC_est = QEC_est;
qec.QEC_calc = QEC_calc;
qec_err.QEC_err = QEC_err;
qec_err.QEC_err_est = QEC_err_est;
qec_err.QEC_err_calc = QEC_err_calc;
} else if (indP != -1) { // only the mass excess of the Parent is known; record the fact that we cannot compute QEC
QEC = NAN; // mimicking assignment order and appearance from just above
QEC_est = 0;
QEC_calc = 0;
QEC_err = NAN;
QEC_err_est = 0;
QEC_err_calc = 0;
qec.QEC = QEC;
qec.QEC_est = QEC_est;
qec.QEC_calc = QEC_calc;
qec_err.QEC_err = QEC_err;
qec_err.QEC_err_est = QEC_err_est;
qec_err.QEC_err_calc = QEC_err_calc;
} else { // no entry in the tree corresponds to this case
continue;
}
a->GetEntry(indP); // navigate (back) to the index of the Parent nucleus ...
newBranch->Fill(); // ... and record the QEC value and the boolean values there
newErrorBranch->Fill(); // Same for error on QEC value
}
}
a->Write("", TObject::kOverwrite); // save only the new version of the tree
......
Supports Markdown
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