44int main(
int argc,
char **argv) {
52 double fractionThreshold;
58 JParser<> zap(
"Auxiliary program to build supernova background from JKexing2D output");
60 zap[
'f'] =
make_field(inputFile,
"input file (JKexing2D).");
62 zap[
'w'] =
make_field(windowSize,
"size of the sliding window to test") = 5;
64 zap[
'F'] =
make_field(fractionThreshold,
"minimum fraction of active channels to compute distribution") = 0.99;
69 catch(
const exception &error) {
70 FATAL(error.what() << endl);
73 typedef JManager<int, TH1D> JManager_i1D_t;
74 typedef JManager<string, TH1D> JManager_s1D_t;
75 typedef JManager<string, TH2D> JManager2D_t;
79 TFile in(inputFile.c_str(),
"exist");
81 TParameter<int>* runNumber;
83 in.GetObject(
"RUNNR", runNumber);
85 JManager2D_t MT = JManager2D_t::Read(in,
mul_p ,
'%');
86 JManager_s1D_t ST = JManager_s1D_t::Read(in,
status_p,
'%');
92 const int factoryLimit_peak = 250;
93 const int factoryLimit_runs = 10000;
101 for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
102 string title =
"MD_" + (*f) +
"_%";
103 mmap[*f] = JManager_i1D_t(
new TH1D(title.c_str(), NULL, factoryLimit_peak, 0, factoryLimit_peak));
106 const int nx = 1 + NUMBER_OF_PMTS;
107 const double xmin = -0.5;
108 const double xmax = nx - 0.5;
112 JManager_s1D_t TD(
new TH1D(Form(
"SNT_[%d,%d]_", multiplicity.
first, multiplicity.
second) + TString(
"%"), NULL, factoryLimit_peak, -0.5, -0.5 + factoryLimit_peak));
114 double epsilon = 1e-4;
117 JManager2D_t BL(
new TH2D(
"BL_%", NULL, 100, epsilon, 1 + epsilon, factoryLimit_peak, -0.5, -0.5 + factoryLimit_peak));
118 JManager2D_t MUL_EFF(
new TH2D(
"MUL_EFF_%", NULL, 100, epsilon, 1 + epsilon, nx, xmin, xmax));
121 JManager_s1D_t H(
new TH1D(
"H_%", NULL, factoryLimit_runs, 0, factoryLimit_runs));
125 TH1D* LT =
new TH1D(
"LIVETIME_ACF", NULL, 100, epsilon, 1 + epsilon);
129 TGraph* activeChannelFraction = histogramToGraph(*ST[
"PMT"]);
131 const int nb = activeChannelFraction->GetN();
132 const Double_t* arr_acf_y = activeChannelFraction->GetY();
136 LT->FillN(nb, arr_acf_y, NULL);
138 for (
int M = 0; M <= NUMBER_OF_PMTS; M++) {
144 for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
148 const TH1* px = projectHistogram(*MT[*f], M, M,
'x');
150 count_vs_frame[*f] = histogramToGraph(*px);
162 transform(vec_acf_y.begin(),
165 back_inserter(select),
166 greater_equal<double> {});
168 mmap[*f][M]->FillN(nb, count_vs_frame[*f]->GetY(), &select[0]);
174 MUL_EFF[*f]->FillN(nb,
177 count_vs_frame[*f]->GetY());
182 int RUNNR = runNumber->GetVal();
184 for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
190 TGraph* events_per_frame = histogramToGraph(*px);
194 TD[*f]->FillN(events_per_frame->GetN(), events_per_frame->GetY(), NULL);
199 vector<double> vec(events_per_frame->GetY(), events_per_frame->GetY() + nb);
204 TD[(*f) +
"_nTS"]->FillN(events_per_frame_sliding.size(),
205 &events_per_frame_sliding[0],
212 events_per_frame->GetY(),
219 int peak = px->GetMaximum();
223 H[
"PK_" + (*f)]->Fill(RUNNR, peak);
229 double BI = (ST[
"PMT"]->GetSumOfWeights() / ST[
"PMT"]->GetEntries());
231 H[
"BIOLUM"]->Fill(RUNNR, BI);
239 TDirectory* bm = out.mkdir(
"BIOLUM");
240 TDirectory* hs = out.mkdir(
"HISTORY");
247 for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
248 string dir_name =
"MUL_" + (*f);
249 TDirectory* dir = out.mkdir(dir_name.c_str());
250 mmap[*f].Write(*dir);