Jpp  17.2.1-pre0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
15 #include "JPhysics/JConstants.hh"
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 
33 #include "JSupernova/JKexing2D.hh"
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  */
44 int 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 
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1517
int main(int argc, char *argv[])
Definition: Main.cc:15
Auxiliary methods for geometrical methods.
static const char * mul_p
Definition: JKexing2D.hh:11
static const char * status_p
Definition: JKexing2D.hh:12
o $QUALITY_ROOT d $DEBUG!JPlot1D f
Definition: JDataQuality.sh:66
static const double H
Planck constant [eV s].
Dynamic ROOT object management.
string outputFile
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
ROOT I/O of application specific meta data.
Physics constants.
General purpose messaging.
TGraph * histogramToGraph(const TH1 &h1)
Helper method to convert a 1D histogram to a graph.
#define FATAL(A)
Definition: JMessage.hh:67
const double xmin
Definition: JQuadrature.cc:23
Auxiliary class to define a range between two values.
Utility class to parse command line options.
std::vector< T > convolve(const std::vector< T > &input, const std::vector< T > &kernel)
Convolute data with given kernel.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
const double epsilon
Definition: JQuadrature.cc:21
int debug
debug level
TH1 * projectHistogram(const TH2 &h2, const Double_t xmin, const Double_t xmax, const char projection)
Helper method for ROOT histogram projections.