Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JMergeCalibrateK40.cc File Reference

Auxiliary program to merge K40 data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TKey.h"
#include "TString.h"
#include "km3net-dataformat/online/JDAQHeader.hh"
#include "JROOT/JRootFileReader.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JSupport/JMeta.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to merge K40 data.

Author
mkarel, mdejong

Definition in file JMergeCalibrateK40.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 40 of file JMergeCalibrateK40.cc.

41{
42 using namespace std;
43 using namespace JPP;
44 using namespace KM3NETDAQ;
45
46 vector<string> inputFile;
47 string outputFile;
48 int debug;
49
50 try {
51
52 JProperties properties;
53
54 properties.insert(gmake_property(LIVETIME_S));
55 properties.insert(gmake_property(ERROR));
56
57 JParser<> zap("Auxiliary program to merge K40 data.");
58
59 zap['@'] = make_field(properties) = JPARSER::initialised();
60 zap['f'] = make_field(inputFile, "input file (output from JCalibrateK40).");
61 zap['o'] = make_field(outputFile, "output file (input to JFitK40).") = "mergek40.root";
62 zap['d'] = make_field(debug, "debug.") = 1;
63
64 zap(argc, argv);
65 }
66 catch(const exception &error) {
67 FATAL(error.what() << endl);
68 }
69
70
71 gErrorIgnoreLevel = kError;
72
73 TH1* weights_hist = NULL;
74
75 typedef map<TString, TH1*> map_type;
76
77 map_type zmap;
78
79
80 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
81
82 DEBUG("Processing " << *i << endl) ;
83
84 TFile in(i->c_str(), "read");
85
86 // check if file contains weights vector
87
88 if (!in.GetListOfKeys()->Contains(weights_hist_t)) {
89 FATAL("Histogram does not contain weights histogram." << endl);
90 }
91
92 TH1* weights_hist_i = (TH1*) in.Get(weights_hist_t);
93
94 if (weights_hist == NULL) {
95
96 weights_hist = (TH1*) weights_hist_i->Clone();
97
98 } else {
99
100 // check normalization
101
102 TAxis* x_axis = weights_hist->GetXaxis();
103
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);
107 }
108
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);
112 }
113
114 // add normalization
115
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)));
119 }
120
121
122 TIter iter(in.GetListOfKeys());
123
124 for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
125
126 if (TString(key->GetName()).EndsWith(_2S) ||
127 TString(key->GetName()).EndsWith(_1B) ||
128 TString(key->GetName()).EndsWith(_1L)) {
129
130 TH1* h1 = dynamic_cast<TH1*>(key->ReadObj());
131
132 map_type::iterator p = zmap.find(h1->GetName());
133
134 if (p == zmap.end()) {
135
136 DEBUG("Clone " << h1->GetName() << endl);
137
138 p = zmap.insert(make_pair(h1->GetName(), (TH1*) h1->Clone())).first;
139
140 } else {
141
142 DEBUG("Add " << h1->GetName() << endl);
143
144 p->second->Add(h1);
145 }
146 }
147 }
148
149 weights_hist->SetDirectory(0);
150
151 for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
152 i->second->SetDirectory(0);
153 }
154
155 in.Close();
156 }
157
158
159 if (weights_hist == NULL) {
160 FATAL("Missing histogram " << weights_hist_t << endl);
161 }
162
163
164 // write file
165
166 TFile out(outputFile.c_str(), "recreate");
167
168 const Double_t WS = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WS_t)) ; // [ns]
169 const Double_t WB = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WB_t)) ; // [ns]
170
171 for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
172
173 if (i->first.EndsWith(_2S)) {
174
175 TH2D* h2s = (TH2D*) i->second;
176
177 if (h2s->GetEntries() != 0) {
178
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;
181
182 //----------------------------------------------------------
183 // normalisation and background estimation
184 //----------------------------------------------------------
185
186 for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
187
188 // normalised background estimation
189 // if livetime is compatible with zero then set bin contents
190 // and erros to zero and default value, respectively
191
192 if (h1L->GetBinContent(ix) <= LIVETIME_S) {
193
194 h1b->SetBinContent(ix, 0.0);
195 h1b->SetBinError (ix, 0.0);
196
197 for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
198 h2s->SetBinContent(ix, iy, 0.0);
199 h2s->SetBinError (ix, iy, ERROR);
200 }
201
202 } else {
203
204 const double Y = h1b->GetBinContent(ix) / h1L->GetBinContent(ix) / WB;
205 const double W = h1b->GetBinError (ix) / h1L->GetBinContent(ix) / WB;
206
207 h1b->SetBinContent(ix, Y);
208 h1b->SetBinError (ix, W);
209
210 for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
211
212 double y = h2s->GetBinContent(ix, iy) / h1L->GetBinContent(ix) / WS;
213 double w = h2s->GetBinError (ix, iy) / h1L->GetBinContent(ix) / WS;
214
215 if (w == 0.0) {
216
217 // for chi2 minimisation; accordingly set error of empty bins
218
219 w = 1.0 / h1L->GetBinContent(ix) / WS;
220 }
221
222 h2s->SetBinContent(ix, iy, y - Y);
223 h2s->SetBinError (ix, iy, sqrt(w*w + W*W));
224 }
225 }
226 }
227
228 // reset under and overflows.
229
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);
233 }
234
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);
238 }
239
240 NOTICE("Biased coincidence rate [Hz] " << h2s->GetName() << ' ' << h2s->GetSumOfWeights() * WS << endl);
241
242 h2s->SetName(TString(i->first).ReplaceAll(_2S, _2R));
243
244 h2s->Write();
245 }
246 }
247 }
248
249 putObject(out, JMeta(argc, argv));
250
251 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
252
253 JMeta::copy(i->c_str(), out);
254
255 for (JRootFileReader<JDAQHeader> in(i->c_str()); in.hasNext(); ) {
256 putObject(out, *in.next());
257 }
258 }
259
260 out.Write();
261 out.Close();
262}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition JProperties.hh:1293
JEEP::JProperties
Utility class to parse parameter values.
Definition JProperties.hh:501
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JCALIBRATE::_2S
static const char *const _2S
Name extension for 2D counts.
Definition JCalibrateK40.hh:33
JCALIBRATE::WB_t
static const char *const WB_t
Named bin for value of TMax_ns in JCalibrateK40.
Definition JCalibrateK40.hh:31
JCALIBRATE::weights_hist_t
static const char *const weights_hist_t
Histogram naming.
Definition JCalibrateK40.hh:27
JCALIBRATE::_1B
static const char *const _1B
Name extension for 1D background.
Definition JCalibrateK40.hh:34
JCALIBRATE::WS_t
static const char *const WS_t
Named bin for time residual bin width.
Definition JCalibrateK40.hh:30
JCALIBRATE::_2R
static const char *const _2R
Name extension for 2D rate measured.
Definition JCalibrateK40.hh:36
JCALIBRATE::W1_overall_t
static const char *const W1_overall_t
Named bin for duration of the run.
Definition JCalibrateK40.hh:29
JCALIBRATE::_1L
static const char *const _1L
Name extension for 1D live time.
Definition JCalibrateK40.hh:35
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JROOT::putObject
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Definition JRootFileWriter.hh:41
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
KM3NETDAQ::map_type
std::map< int, range_type > map_type
Definition JDQSimulator.cc:161
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
JSUPPORT::JMeta::copy
static void copy(const char *const file_name, TFile &out)
Copy meta data.
Definition JMeta.hh:421
Generated by doxygen 1.12.0