expand-nuchart-QECp.cxx 4.14 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
/*
 * 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);
83
        if (QEC.QEC_calc == 0) { // this is new compared to 'expand-nuchart-QEC.cxx': QEC might not be defined, so we need to check for this
84
85
          QECp_NAN();
        } else {
86
          if (QEC.QEC_est == 1) {
87
88
            QECp_est = 1; QECp_err_est = 1;
          }
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
          QECP = QEC.QEC;
          QECP_err = QEC_err.QEC_err;
          a->GetEntry(indD);
          if (Sp.Sp_calc == 0) { // same goes for Sp
            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);
104

105
106
107
108
109
110
111
            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;
          }
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
        }
      } 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;
}