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
Functions
JFitK40.cc File Reference

Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <set>
#include <vector>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQHeader.hh"
#include "JTools/JConstants.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JStringRouter.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JLang/JSTDObjectWriter.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JTDC_t.hh"
#include "JFitK40.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.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 fit PMT parameters from JMergeCalibrateK40.cc output.


Note that the contebts of the 2D histograms are converted to standard deviations.

Author
mdejong

Definition in file JFitK40.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 60 of file JFitK40.cc.

61 {
62  using namespace std;
63  using namespace JPP;
64  using namespace KM3NETDAQ;
65 
66  typedef JRange<double> JRange_t;
67 
68  JK40Parameters K40 = JK40Parameters::getInstance();
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 :
126  fitAngle ? FIT_PMTS_AND_ANGULAR_DEPENDENCE_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 
144  JDetector detector;
145 
146  try {
147  load(detectorFile, detector);
148  }
149  catch(const JException& error) {
150  FATAL(error);
151  }
152 
153 
154  JPMTParametersMap parameters;
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 
185  const floor_range range = getRangeOfFloors(detector);
186  const JStringRouter string(detector);
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 }
JCALIBRATE::rate_type
Data structure for measured coincidence rate of pair of PMTs.
Definition: JFitK40.hh:62
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1514
JLANG::JException
General exception.
Definition: JException.hh:24
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1292
JDETECTOR::getLabel
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:246
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:89
JCALIBRATE::JTDC_t
Auxiliary class for TDC constraints.
Definition: JTDC_t.hh:37
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:497
JTOOLS::JCombinatorics::pair_type
Data structure for a pair of indices.
Definition: JCombinatorics.hh:27
parameters
*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
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:136
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
Y
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
Definition: detector-XY:fit1d.sh:33
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
V
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition: JDetectorToolkit.hh:443
JACOUSTICS::JFit
Acoustic single fit.
Definition: JAcoustics/JEvt.hh:35
JLANG::JSTDObjectWriter
Implementation of object output from STD container.
Definition: JSTDObjectWriter.hh:24
JACOUSTICS::JModel
Model for fit to acoustics data.
Definition: JAcoustics/JModel.hh:476
JCALIBRATE::FIT_PMTS_t
fit parameters of PMTs
Definition: JFitK40.hh:51
JAANET::JRange_t
Type definition of range.
Definition: JHead.hh:41
JCALIBRATE::_2F
static const char *const _2F
Name extension for 2F rate fitted.
Definition: JCalibrateK40.hh:37
JAANET::detector
Detector file.
Definition: JHead.hh:226
JCALIBRATE::JK40Parameters
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:597
JLANG::getInstance
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:75
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:578
JCALIBRATE::FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
fit parameters of PMTs and angular dependence of K40 rate
Definition: JFitK40.hh:53
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
string
then awk string
Definition: pre-calibration.sh:45
JCALIBRATE::_2R
static const char *const _2R
Name extension for 2D rate measured.
Definition: JCalibrateK40.hh:36
JDETECTOR::JPMTParametersMap
Auxiliary class for map of PMT parameters.
Definition: JPMTParametersMap.hh:89
JLANG::getIndex
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
JROOT::putObject
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Definition: JRootFileWriter.hh:41
FILL
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:328
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
range
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
Definition: module-Z:fit.sh:118
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:485
JDETECTOR::JStringRouter
Router for mapping of string identifier to index.
Definition: JStringRouter.hh:21
JCALIBRATE::FIT_MODEL_t
fit parameters of K40 rate and TTSs of PMTs
Definition: JFitK40.hh:52
JDETECTOR::getSurvivalProbability
double getSurvivalProbability(const JPMTParameters &parameters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
Definition: JPMTParametersToolkit.hh:32
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition: JFitK40.hh:98
TDC
set_variable TDC
Definition: JPrintTDC.sh:20
X
no fit printf nominal n $STRING awk v X
Definition: detector-XY:fit.sh:102
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition: JSingleFileScanner.hh:75
JDETECTOR::JPMTParameters
Data structure for PMT parameters.
Definition: JPMTParameters.hh:29
detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:38
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
in
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
JDETECTOR::JPMTParameters::QE
double QE
relative quantum efficiency
Definition: JPMTParameters.hh:233
JAANET::K40
Livetime of noise data.
Definition: JHead.hh:1062
JROOT::JRootFileReader
ROOT file reader.
Definition: JRootFileReader.hh:54
debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92
P
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5