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

added QECp and QECp_err branch

parent 2637b971
......@@ -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*, *graph-example-qec.py*
* Do TODO's in top of *datagetter.py*
# DONE
......
/*
* In this ROOT macro the ROOT Tree 'nuchart.root' generated by treemaker.py is expanded by
* adding Q-values for electron capture and subsequent proton emission to the tree via the formula
* QECp = QECP - SpD;
* QECP is the Q-value for the process,
* QECP is the QEC value of the Parent which is decaying,
* SpD is the proton separation energy of the Daughter which the Parent decays into.
*
* Note: The structure of this code is very similar to the structure of
* 'expand-nuchart-QEC.cxx'. Please refer to that for an overview.
*/
#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, SpD = 0.0;
Double_t QECP_err = 0.0, SpD_err = 0.0;
Double_t QECp = 0.0;
Double_t QECp_err = 0.0;
Bool_t QECp_est = 0, QECp_calc = 0;
Bool_t QECp_err_est = 0, QECp_err_calc = 0;
struct QECp_params { Double_t QECp; Bool_t QECp_est; Bool_t QECp_calc; };
QECp_params qecp = QECp_params();
struct QECp_err_params { Double_t QECp_err; Bool_t QECp_err_est; Bool_t QECp_err_calc; };
QECp_err_params qecp_err = QECp_err_params();
void QECp_NAN() {
QECp = NAN;
QECp_est = 0;
QECp_calc = 0;
QECp_err = NAN;
QECp_est = 0;
QECp_calc = 0;
qecp.QECp = QECp;
qecp.QECp_est = QECp_est;
qecp.QECp_calc = QECp_calc;
qecp_err.QECp_err = QECp_err;
qecp_err.QECp_err_est = QECp_err_est;
qecp_err.QECp_err_calc = QECp_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 Sp_params { Double_t Sp; Bool_t Sp_est; Bool_t Sp_calc; };
Sp_params Sp = Sp_params();
a->SetBranchAddress("Sp", &Sp);
struct Sp_err_params { Double_t Sp_err; Bool_t Sp_err_est; Bool_t Sp_err_calc; };
Sp_err_params Sp_err = Sp_err_params();
a->SetBranchAddress("Sp_err", &Sp_err);
TBranch *qecp_b = a->Branch("QECp", &qecp.QECp, "QECp/D:QECp_est/O:QECp_calc");
TBranch *qecp_err_b = a->Branch("QECp_err", &qecp_err.QECp_err, "QECp_err/D:QECp_err_est/O:QECp_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) {
QECp_calc = 1; QECp_err_calc = 1;
a->GetEntry(indP);
if (QEC.QEC_est == 1) {
QECp_est = 1; QECp_err_est = 1;
}
QECP = QEC.QEC;
QECP_err = QEC_err.QEC_err;
a->GetEntry(indD);
if (Sp.Sp_calc == 0) { // this is new compared to 'expand-nuchart-QEC.cxx': Sp might not be defined, so we need to check for this
QECp_NAN();
} else {
if (Sp.Sp_est == 1) {
QECp_est = 1; QECp_err_est = 1;
} else {
QECp_est = 0; QECp_err_est = 0;
}
SpD = Sp.Sp;
SpD_err = Sp_err.Sp_err;
QECp = QECP - SpD;
QECp_err = Sqrt(QECP_err*QECP_err + SpD_err*SpD_err);
qecp.QECp = QECp;
qecp.QECp_est = QECp_est;
qecp.QECp_calc = QECp_calc;
qecp_err.QECp_err = QECp_err;
qecp_err.QECp_err_est = QECp_err_est;
qecp_err.QECp_err_calc = QECp_err_calc;
}
} else if (indP != -1) {
QECp_NAN();
} else {
continue;
}
a->GetEntry(indP);
qecp_b->Fill();
qecp_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