40int main(
int argc,
char **argv)
57 JParser<> zap(
"Auxiliary program to merge K40 data.");
60 zap[
'f'] =
make_field(inputFile,
"input file (output from JCalibrateK40).");
66 catch(
const exception &error) {
67 FATAL(error.what() << endl);
71 gErrorIgnoreLevel = kError;
73 TH1* weights_hist = NULL;
80 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
82 DEBUG(
"Processing " << *i << endl) ;
84 TFile in(i->c_str(),
"read");
89 FATAL(
"Histogram does not contain weights histogram." << endl);
94 if (weights_hist == NULL) {
96 weights_hist = (TH1*) weights_hist_i->Clone();
102 TAxis* x_axis = weights_hist->GetXaxis();
104 if (weights_hist_i->GetBinContent(x_axis->FindBin(
WS_t)) !=
105 weights_hist ->GetBinContent(x_axis->FindBin(
WS_t))) {
106 FATAL(
"Number of y-bins does not match for this file " << *i << endl);
109 if (weights_hist_i->GetBinContent(x_axis->FindBin(
WB_t)) !=
110 weights_hist ->GetBinContent(x_axis->FindBin(
WB_t))) {
111 FATAL(
"TMax_ns is not the same for this file " << *i << endl);
116 weights_hist->SetBinContent(x_axis->FindBin(
W1_overall_t),
117 weights_hist ->GetBinContent(x_axis->FindBin(
W1_overall_t)) +
118 weights_hist_i->GetBinContent(x_axis->FindBin(
W1_overall_t)));
122 TIter iter(in.GetListOfKeys());
124 for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
126 if (TString(key->GetName()).EndsWith(
_2S) ||
127 TString(key->GetName()).EndsWith(
_1B) ||
128 TString(key->GetName()).EndsWith(
_1L)) {
130 TH1* h1 =
dynamic_cast<TH1*
>(key->ReadObj());
132 map_type::iterator p = zmap.find(h1->GetName());
134 if (p == zmap.end()) {
136 DEBUG(
"Clone " << h1->GetName() << endl);
138 p = zmap.insert(make_pair(h1->GetName(), (TH1*) h1->Clone())).first;
142 DEBUG(
"Add " << h1->GetName() << endl);
149 weights_hist->SetDirectory(0);
151 for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
152 i->second->SetDirectory(0);
159 if (weights_hist == NULL) {
168 const Double_t WS = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(
WS_t)) ;
169 const Double_t WB = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(
WB_t)) ;
171 for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
173 if (i->first.EndsWith(
_2S)) {
175 TH2D* h2s = (TH2D*) i->second;
177 if (h2s->GetEntries() != 0) {
179 TH1D* h1b = (TH1D*) zmap.find(TString(i->first).ReplaceAll(
_2S,
_1B))->second;
180 TH1D* h1L = (TH1D*) zmap.find(TString(i->first).ReplaceAll(
_2S,
_1L))->second;
186 for (
int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
192 if (h1L->GetBinContent(ix) <= LIVETIME_S) {
194 h1b->SetBinContent(ix, 0.0);
195 h1b->SetBinError (ix, 0.0);
197 for (
int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
198 h2s->SetBinContent(ix, iy, 0.0);
199 h2s->SetBinError (ix, iy,
ERROR);
204 const double Y = h1b->GetBinContent(ix) / h1L->GetBinContent(ix) / WB;
205 const double W = h1b->GetBinError (ix) / h1L->GetBinContent(ix) / WB;
207 h1b->SetBinContent(ix, Y);
208 h1b->SetBinError (ix, W);
210 for (
int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
212 double y = h2s->GetBinContent(ix, iy) / h1L->GetBinContent(ix) / WS;
213 double w = h2s->GetBinError (ix, iy) / h1L->GetBinContent(ix) / WS;
219 w = 1.0 / h1L->GetBinContent(ix) / WS;
222 h2s->SetBinContent(ix, iy, y - Y);
223 h2s->SetBinError (ix, iy, sqrt(w*w + W*W));
230 for (
int ix = 0; ix <= h2s->GetXaxis()->GetNbins() + 1; ++ix) {
231 h2s->SetBinContent(ix, 0, 0.0); h2s->SetBinContent(ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
232 h2s->SetBinError (ix, 0, 0.0); h2s->SetBinError (ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
235 for (
int iy = 0; iy <= h2s->GetYaxis()->GetNbins() + 1; ++iy) {
236 h2s->SetBinContent(0, iy, 0.0); h2s->SetBinContent(h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
237 h2s->SetBinError (0, iy, 0.0); h2s->SetBinError (h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
240 NOTICE(
"Biased coincidence rate [Hz] " << h2s->GetName() <<
' ' << h2s->GetSumOfWeights() * WS << endl);
242 h2s->SetName(TString(i->first).ReplaceAll(
_2S,
_2R));
251 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {