Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JProcessKexing2D.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <map>
5#include <cmath>
6
7#include "TROOT.h"
8#include "TFile.h"
9#include "TH1D.h"
10#include "TH2D.h"
11#include "TF2.h"
12#include "TMath.h"
13#include "TFitResult.h"
14
17
18#include "JROOT/JRootToolkit.hh"
19#include "JTools/JRange.hh"
20#include "JSupport/JMeta.hh"
21
23#include "JCalibrate/JFitK40.hh"
24
25#include "JROOT/JManager.hh"
26
27#include "Jeep/JPrint.hh"
28#include "Jeep/JParser.hh"
29#include "Jeep/JMessage.hh"
30
31#include <string>
32
34
35#include "JROOT/JRootToolkit.hh"
36#include "JMath/JMathToolkit.hh"
37
38/**
39 * \file
40 *
41 * Auxiliary program to build supernova background from JKexing2D output
42 * \author mlincett
43 */
44int main(int argc, char **argv) {
45 using namespace std;
46 using namespace JPP;
47
48 string inputFile;
49 string outputFile;
50 string summaryFile;
51 int windowSize;
52 double fractionThreshold;
53 int debug;
54 JRange<int> multiplicity;
55
56 try {
57
58 JParser<> zap("Auxiliary program to build supernova background from JKexing2D output");
59
60 zap['f'] = make_field(inputFile, "input file (JKexing2D).");
61 zap['o'] = make_field(outputFile, "output file.") = "sn_background.root";
62 zap['w'] = make_field(windowSize, "size of the sliding window to test") = 5;
63 zap['M'] = make_field(multiplicity, "final multiplicity range") = JRange<int>(7,11);
64 zap['F'] = make_field(fractionThreshold, "minimum fraction of active channels to compute distribution") = 0.99;
65 zap['d'] = make_field(debug) = 1;
66
67 zap(argc, argv);
68 }
69 catch(const exception &error) {
70 FATAL(error.what() << endl);
71 }
72
73 typedef JManager<int, TH1D> JManager_i1D_t;
74 typedef JManager<string, TH1D> JManager_s1D_t;
75 typedef JManager<string, TH2D> JManager2D_t;
76
77 // input
78
79 TFile in(inputFile.c_str(), "exist");
80
81 TParameter<int>* runNumber;
82
83 in.GetObject("RUNNR", runNumber);
84
85 JManager2D_t MT = JManager2D_t::Read(in, mul_p , '%');
86 JManager_s1D_t ST = JManager_s1D_t::Read(in, status_p, '%');
87
88 in.Close();
89
90 // output
91
92 const int factoryLimit_peak = 250;
93 const int factoryLimit_runs = 10000;
94
95 // distributions of individual multiplicities
96
97 vector<string> filters = { "RAW", "TAV", "TA1" };
98
100
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));
104 }
105
106 const int nx = 1 + NUMBER_OF_PMTS;
107 const double xmin = -0.5;
108 const double xmax = nx - 0.5;
109
110 // distributions of the trigger level
111
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));
113
114 double epsilon = 1e-4;
115
116 // bioluminescence benchmarking
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));
119
120 // history (value vs run number)
121 JManager_s1D_t H(new TH1D("H_%", NULL, factoryLimit_runs, 0, factoryLimit_runs));
122
123 // livetime vs fraction of active channels
124
125 TH1D* LT = new TH1D("LIVETIME_ACF", NULL, 100, epsilon, 1 + epsilon);
126
127 // fraction of active PMT
128
129 TGraph* activeChannelFraction = histogramToGraph(*ST["PMT"]);
130
131 const int nb = activeChannelFraction->GetN();
132 const Double_t* arr_acf_y = activeChannelFraction->GetY();
133
134 const vector<Double_t> vec_acf_y(arr_acf_y, arr_acf_y + nb);
135
136 LT->FillN(nb, arr_acf_y, NULL);
137
138 for (int M = 0; M <= NUMBER_OF_PMTS; M++) {
139
140 map<string, TGraph*> count_vs_frame;
141
142 // individual multiplicities
143
144 for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
145
146 // number of events with multiplicity M as a function of the frame index
147
148 const TH1* px = projectHistogram(*MT[*f], M, M, 'x');
149
150 count_vs_frame[*f] = histogramToGraph(*px);
151
152 delete px;
153
154 // distribution of the number of events with multiplicity M
155 // selection criterion based on the minimum fraction of active PMTs
156
157 vector<double> select;
158 vector<double> threshold(nb, fractionThreshold);
159
160 // stores in 'select' the result of the element-wise comparison between 'vec_acf_y' and 'threshold'
161
162 transform(vec_acf_y.begin(),
163 vec_acf_y.end(),
164 threshold.begin(),
165 back_inserter(select),
166 greater_equal<double> {});
167
168 mmap[*f][M]->FillN(nb, count_vs_frame[*f]->GetY(), &select[0]);
169
170 // multiplicity spectrum as a function of number of active PMTs
171
172 const vector<double> nb_M(nb, M);
173
174 MUL_EFF[*f]->FillN(nb,
175 arr_acf_y,
176 &nb_M[0],
177 count_vs_frame[*f]->GetY());
178
179 }
180 }
181
182 int RUNNR = runNumber->GetVal();
183
184 for (vector<string>::const_iterator f = filters.begin(); f != filters.end(); ++f) {
185
186 // distribution of number of events per frame in the selected multiplicity range (A)
187
188 const TH1* px = projectHistogram(*MT[*f], multiplicity.getLowerLimit(), multiplicity.getUpperLimit(), 'x');
189
190 TGraph* events_per_frame = histogramToGraph(*px);
191
192 delete px;
193
194 TD[*f]->FillN(events_per_frame->GetN(), events_per_frame->GetY(), NULL);
195
196 // distribution of number of events per frame in the selected multiplicity range (A)
197 // sliding window applied
198
199 vector<double> vec(events_per_frame->GetY(), events_per_frame->GetY() + nb);
200 vector<double> ker(windowSize, 1.0);
201
202 vector<double> events_per_frame_sliding = convolve(vec, ker);
203
204 TD[(*f) + "_nTS"]->FillN(events_per_frame_sliding.size(),
205 &events_per_frame_sliding[0],
206 NULL);
207
208 // distribution of number of events per frame in the selected multiplicity range (A)
209 // as a function of the fraction of active channels
210 BL[*f]->FillN(nb,
211 arr_acf_y,
212 events_per_frame->GetY(),
213 NULL, 1);
214
215 // peak as a function of the run number
216
217 px = projectHistogram(*MT[*f], multiplicity.getLowerLimit(), multiplicity.getUpperLimit(), 'x');
218
219 int peak = px->GetMaximum();
220
221 delete px;
222
223 H["PK_" + (*f)]->Fill(RUNNR, peak);
224
225 }
226
227 // history of bioluminescence
228
229 double BI = (ST["PMT"]->GetSumOfWeights() / ST["PMT"]->GetEntries());
230
231 H["BIOLUM"]->Fill(RUNNR, BI);
232
233 // write output file
234
235 TFile out(outputFile.c_str(), "recreate");
236
237 TD.Write(out);
238
239 TDirectory* bm = out.mkdir("BIOLUM");
240 TDirectory* hs = out.mkdir("HISTORY");
241
242 H.Write(*hs);
243 BL.Write(*bm);
244
245 MUL_EFF.Write(*bm);
246
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);
251 }
252
253 LT->Write();
254
255 out.Close();
256
257}
258
259
string outputFile
static const char * status_p
Definition JKexing2D.hh:12
static const char * mul_p
Definition JKexing2D.hh:11
Dynamic ROOT object management.
Auxiliary methods for geometrical methods.
General purpose messaging.
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
ROOT I/O of application specific meta data.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
I/O formatting auxiliaries.
int main(int argc, char **argv)
Auxiliary class to define a range between two values.
Utility class to parse command line options.
Definition JParser.hh:1698
JKey_t first
Definition JPair.hh:128
JValue_t second
Definition JPair.hh:129
Range of values.
Definition JRange.hh:42
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).