Jpp  18.2.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JFitK40.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <cmath>
5 #include <set>
6 #include <vector>
7 #include <algorithm>
8 
9 #include "TROOT.h"
10 #include "TFile.h"
11 #include "TH1D.h"
12 #include "TH2D.h"
13 
16 
17 #include "JTools/JConstants.hh"
18 #include "JDetector/JDetector.hh"
22 
23 #include "JROOT/JRootFileWriter.hh"
24 #include "JROOT/JRootToolkit.hh"
25 
26 #include "JTools/JRange.hh"
27 
29 
30 #include "JSupport/JMeta.hh"
32 
34 #include "JCalibrate/JTDC_t.hh"
35 
36 #include "JFitK40.hh"
37 
38 #include "Jeep/JProperties.hh"
39 #include "Jeep/JPrint.hh"
40 #include "Jeep/JParser.hh"
41 #include "Jeep/JMessage.hh"
42 
43 
44 namespace {
45 
46  double MINIMAL_RATE_HZ = 1.0e-2; //!< minimal coincidence rate [Hz] (below, disable PMT)
47  double STDEV = 2.0; //!< number of standard deviations
48  double QE_MIN = 0.1; //!< minimal QE (below, disable PMT)
49 }
50 
51 
52 /**
53  * \file
54  *
55  * Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output.\n
56  * Note that the contebts of the 2D histograms are converted to standard deviations.
57  *
58  * \author mdejong
59  */
60 int main(int argc, char **argv)
61 {
62  using namespace std;
63  using namespace JPP;
64  using namespace KM3NETDAQ;
65 
66  typedef JRange<double> JRange_t;
67 
69 
70  string inputFile;
71  string outputFile;
72  string detectorFile;
73  string pmtFile;
74  JTDC_t TDC;
75  bool reverse;
76  bool overwriteDetector;
77  bool writeFits;
78  bool fitAngle;
79  bool fitModel;
80  JRange_t X;
81  JRange_t Y;
82  int debug;
83 
84  try {
85 
86  JProperties properties;
87 
88  properties.insert(gmake_property(MINIMAL_RATE_HZ));
89  properties.insert(gmake_property(STDEV));
90  properties.insert(gmake_property(QE_MIN));
91  properties.insert(gmake_property(K40.R));
92  properties.insert(gmake_property(K40.p1));
93  properties.insert(gmake_property(K40.p2));
94  properties.insert(gmake_property(K40.p3));
95  properties.insert(gmake_property(K40.p4));
96 
97  JParser<> zap("Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
98 
99  zap['@'] = make_field(properties) = JPARSER::initialised();
100  zap['f'] = make_field(inputFile, "input file (output from JMergeCalibrateK40).");
101  zap['o'] = make_field(outputFile, "output file.") = "fit.root";
102  zap['a'] = make_field(detectorFile, "detector file.");
103  zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
104  zap['!'] = make_field(TDC,
105  "Fix time offset(s) of PMT(s) of certain module(s), e.g."
106  "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
107  "\nSame PMT offset can be fixed for all optical modules, e.g."
108  "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") = JPARSER::initialised();
109  zap['r'] = make_field(reverse, "reverse TDC constraints due to option -! <TDC>.");
110  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with fitted time offsets.");
111  zap['w'] = make_field(writeFits, "write fit results.");
112  zap['D'] = make_field(fitAngle, "fit angular distribution; fix normalisation.");
113  zap['M'] = make_field(fitModel, "fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
114  zap['x'] = make_field(X, "fit range (PMT pairs).") = JRange_t();
115  zap['y'] = make_field(Y, "fit range (time residual [ns]).") = JRange_t();
116  zap['d'] = make_field(debug, "debug.") = 1;
117 
118  zap(argc, argv);
119  }
120  catch(const exception &error) {
121  FATAL(error.what() << endl);
122  }
123 
124 
125  const int option = (fitModel ? FIT_MODEL_t :
127  FIT_PMTS_t);
128 
129  if (reverse) {
130  TDC.reverse();
131  }
132 
133  for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
134  DEBUG("PMT " << setw(10) << i->first << ' ' << setw(2) << i->second << " constrain t0." << endl);
135  }
136 
137  try {
138  TDC.is_valid(true);
139  }
140  catch(const exception &error) {
141  FATAL(error.what() << endl);
142  }
143 
145 
146  try {
147  load(detectorFile, detector);
148  }
149  catch(const JException& error) {
150  FATAL(error);
151  }
152 
153 
155 
156  if (pmtFile != "") {
157  try {
158  parameters.load(pmtFile.c_str());
159  }
160  catch(const exception& error) {}
161  }
162 
163  parameters.comment.add(JMeta(argc, argv));
164 
165 
166  TFile* in = TFile::Open(inputFile.c_str(), "exist");
167 
168  if (in == NULL || !in->IsOpen()) {
169  FATAL("File: " << inputFile << " not opened." << endl);
170  }
171 
172 
173  TFile out(outputFile.c_str(), "recreate");
174 
175  TH1D h0("chi2", NULL, 500, 0.0, 5.0);
176  TH1D hn("hn", NULL, 501, -0.5, 500.0);
177  TH1D hr("rate", NULL, 500, 0.0, 25.0);
178  TH1D h1("p1", NULL, 500, -5.0, +5.0);
179  TH1D h2("p2", NULL, 500, -5.0, +5.0);
180  TH1D h3("p3", NULL, 500, -5.0, +5.0);
181  TH1D h4("p4", NULL, 500, -5.0, +5.0);
182  TH1D hb("bg", NULL, 500, -0.1, +0.1);
183  TH1D hc("cc", NULL, 500, -0.1, +0.1);
184 
187 
188  TH2D H2("detector", NULL,
189  string.size() + 0, -0.5, string.size() - 0.5,
190  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5);
191 
192  for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
193  H2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
194  }
195  for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
196  H2.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
197  }
198 
199 
200  JFit fit(0, debug);
201 
202  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
203 
204  if (module->getFloor() == 0) {
205  continue;
206  }
207 
208  const JTDC_t::range_type range = TDC.equal_range(module->getID());
209 
210  NOTICE("Module " << setw(10) << module->getID() << ' ' << getLabel(module->getLocation()) << " !" << distance(range.first, range.second) << endl);
211 
212  TH2D* h2d = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2R));
213 
214  if (h2d == NULL || h2d->GetEntries() == 0) {
215 
216  NOTICE("No data for module " << module->getID() << " -> set QEs to 0." << endl);
217 
218  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
219  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
220  }
221 
222  continue;
223  }
224 
225  JModel model(*module, K40, range, option);
226 
227  data_type data; // input data
228 
229  vector<size_t> count(NUMBER_OF_PMTS, 1); // counter for to low integrated rates
230 
231  for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
232 
233  const pair_type pair = model.getPair(ix - 1);
234 
235  auto& buffer = data[pair]; // storage
236 
237  double V = 0.0; // integrated value
238  double W = 0.0; // integrated error
239 
240  for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
241 
242  const double x = h2d->GetXaxis()->GetBinCenter(ix);
243  const double y = h2d->GetYaxis()->GetBinCenter(iy);
244 
245  if (X(x) && Y(y)) {
246 
247  double value = h2d->GetBinContent(ix,iy);
248  double error = h2d->GetBinError (ix,iy);
249 
250  buffer.push_back(rate_type(y, value, error));
251 
252  double width = h2d->GetYaxis()->GetBinWidth(iy);
253 
254  value *= width;
255  error *= width;
256 
257  V += value;
258  W += error * error;
259  }
260  }
261 
262  if (V <= MINIMAL_RATE_HZ + STDEV*sqrt(W)) {
263  count[pair.first] += 1;
264  count[pair.second] += 1;
265  }
266  }
267 
268  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
269 
270  if (count[pmt] == NUMBER_OF_PMTS) { // all paired rates to low
271 
272  WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << ' ' << "all rates to low -> disable" << endl);
273 
274  model.parameters[pmt].disable();
275  }
276  }
277 
278  DEBUG("Start value:" << endl << model << endl);
279 
280  try {
281 
282  fit.value = model; // start value
283 
284  auto result = fit(data);
285 
286  bool refit = false;
287 
288  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
289 
290  if (fit.value.parameters[pmt].QE() < QE_MIN && fit.value.parameters[pmt].status) {
291 
292  WARNING("PMT " << setw(10) << module->getID() << '.' << FILL(2,'0') << pmt << FILL() << ' ' << FIXED(5,3) << fit.value.parameters[pmt].QE() << " < " << FIXED(5,3) << QE_MIN << " -> disable" << (!refit ? " and refit" : "") << endl);
293 
294  fit.value.parameters[pmt].disable();
295 
296  refit = true;
297  }
298  }
299 
300  if (refit) {
301  result = fit(data);
302  }
303 
304  NOTICE("Fit result " << setw(10) << module->getID() << " chi2 / NDF " << FIXED(10,2) << result.chi2 << " / " << setw(5) << result.ndf << ' ' << setw(5) << fit.numberOfIterations << endl);
305 
306  DEBUG(fit.value);
307 
308  if (writeFits) {
309 
310  h0.Fill(result.chi2 / result.ndf);
311  hn.Fill(fit.numberOfIterations);
312  hr.Fill(fit.value.R );
313  h1.Fill(fit.value.p1);
314  h2.Fill(fit.value.p2);
315  h3.Fill(fit.value.p3);
316  h4.Fill(fit.value.p4);
317  hb.Fill(fit.value.bg);
318  hc.Fill(fit.value.cc);
319 
320  for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
321 
322  const pair_type pair = fit.value.getPair(ix - 1);
323 
324  for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
325 
326  const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
327 
328  h2d->SetBinContent(ix, iy, fit.value.getValue(pair, dt_ns));
329  h2d->SetBinError (ix, iy, 0.0);
330  }
331  }
332 
333  h2d->SetName(MAKE_CSTRING(module->getID() << _2F));
334  h2d->Write();
335 
336  H2.Fill((double) string.getIndex(module->getString()), (double) module->getFloor(), result.chi2 / result.ndf);
337  }
338 
339  const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
340 
341  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
342 
343  JPMTParameters& data = parameters[JPMTIdentifier(module->getID(), pmt)];
344 
345  const double P = getSurvivalProbability(data);
346 
347  if (P > 0.0)
348  data.QE = fit.value.parameters[pmt].QE / P;
349  else
350  data.QE = 0.0;
351 
352  module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
353  }
354  }
355  catch(const exception& error) {
356 
357  ERROR("Module " << setw(10) << module->getID() << ' ' << error.what() << " -> set QEs to 0." << endl);
358 
359  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
360  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
361  }
362  }
363  }
364 
365 
366  vector<JMeta> meta(1, JMeta(argc, argv));
367 
368  {
369  JSingleFileScanner<JMeta> reader(inputFile);
370  JSTDObjectWriter <JMeta> writer(meta);
371 
372  writer << reader;
373  }
374 
375  for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
376  parameters.comment.add(*i);
377  detector .comment.add(*i);
378  }
379 
380  if (overwriteDetector) {
381 
382  NOTICE("Store calibration data on file " << detectorFile << endl);
383 
384  store(detectorFile, detector);
385  }
386 
387  if (pmtFile != "") {
388  parameters.store(pmtFile.c_str());
389  }
390 
391 
392  for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
393  putObject(out,*i);
394  }
395 
396  for (JRootFileReader<JDAQHeader> in(inputFile.c_str()); in.hasNext(); ) {
397  putObject(out, *in.next());
398  }
399 
400  if (writeFits) {
401  out << h0 << hn << hr << h1 << h2 << h3 << h4 << hb << hc << H2;
402  }
403 
404  out.Close();
405 
406  return 0;
407 }
Data structure for measured coincidence rate of pair of PMTs.
Definition: JFitK40.hh:62
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:24
#define WARNING(A)
Definition: JMessage.hh:65
int main(int argc, char *argv[])
Definition: Main.cc:15
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:246
Detector data structure.
Definition: JDetector.hh:89
Auxiliary class for TDC constraints.
Definition: JTDC_t.hh:37
Utility class to parse parameter values.
Definition: JProperties.hh:497
Data structure for a pair of indices.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:136
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Data structure for detector geometry and calibration.
Utility class to parse parameter values.
Acoustic single fit.
Implementation of object output from STD container.
Scanning of objects from a single file according a format that follows from the extension of each fil...
Model for fit to acoustics data.
fit parameters of PMTs
Definition: JFitK40.hh:51
Type definition of range.
Definition: JHead.hh:41
static const char *const _2F
Name extension for 2F rate fitted.
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:226
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:597
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:75
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
ROOT I/O of application specific meta data.
fit parameters of PMTs and angular dependence of K40 rate
Definition: JFitK40.hh:53
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
then awk string
static const char *const _2R
Name extension for 2D rate measured.
Auxiliary class for map of PMT parameters.
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
General purpose messaging.
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:328
Direct access to string in detector data structure.
#define FATAL(A)
Definition: JMessage.hh:67
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
Constants.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
Router for mapping of string identifier to index.
Utility class to parse command line options.
fit parameters of K40 rate and TTSs of PMTs
Definition: JFitK40.hh:52
double getSurvivalProbability(const JPMTParameters &parameters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition: JFitK40.hh:98
set_variable TDC
Definition: JPrintTDC.sh:20
no fit printf nominal n $STRING awk v X
Object reading from a list of files.
Data structure for PMT parameters.
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
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
double QE
relative quantum efficiency
Livetime of noise data.
Definition: JHead.hh:1062
ROOT file reader.
int debug
debug level
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62