Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPulsar.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <map>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TF1.h"
10 #include "TFitResult.h"
11 
12 #include "JROOT/JMinimizer.hh"
13 
14 #include "JDAQ/JDAQTimesliceIO.hh"
15 
16 #include "JDetector/JDetector.hh"
20 
21 #include "JTrigger/JHit.hh"
22 #include "JTrigger/JFrame.hh"
23 #include "JTrigger/JTimeslice.hh"
24 #include "JTrigger/JHitR0.hh"
25 #include "JTrigger/JBuildL0.hh"
26 
27 #include "JCalibrate/JTDC_t.hh"
28 
30 #include "JSupport/JSupport.hh"
31 #include "JSupport/JMeta.hh"
32 
33 #include "JTools/JRange.hh"
36 
37 #include "Jeep/JProperties.hh"
38 #include "Jeep/JParser.hh"
39 #include "Jeep/JMessage.hh"
40 
41 
42 namespace {
43 
44  /**
45  * Get time within given period.
46  *
47  * The returned time is constraint to the interval <tt>[-T/2,+T/2]</tt>.
48  *
49  * \param t1 time
50  * \param T period
51  * \return time
52  */
53  inline double get_time(double t1, const double T)
54  {
55  for (int buffer[] = { 1000, 100, 10, 1, 0 }, *i = buffer; *i != 0; ++i) {
56 
57  const double xmin = -0.5 * (*i) * T;
58  const double xmax = +0.5 * (*i) * T;
59 
60  while (t1 > xmax) { t1 -= (*i) * T; }
61  while (t1 < xmin) { t1 += (*i) * T; }
62  }
63 
64  return t1;
65  }
66 }
67 
68 
69 /**
70  * \file
71  *
72  * Application for dark room time calibration using external laser.
73  *
74  * The time calibration is made for the specified reference PMTs (option <tt>-! "<module identifier> <PMT channel>"</tt>.\n
75  * In this, the wild card value <tt>-1</tt> can be used.
76  * A ROOT histogram corresponding to each reference PMT is created, filled, fitted and stored to the specified output file.\n
77  * The detector can be updated using option <tt>-A</tt>\n
78  * The laser should be synchronised with the same clock system as the detector and
79  * produce single pulses at a fixed frequency (option <tt>-l <laserFrequency_Hz></tt>.
80  *
81  * In addition, the transit-time distribution of the reference PMTs can be monitored (see JLegolas.cc).\n
82  * For this, the expectation values should be less than one and the time-over-threshold distribution nominal.
83  *
84  * \author mbouwhuis
85  */
86 int main(int argc, char **argv)
87 {
88  using namespace std;
89  using namespace JPP;
90  using namespace KM3NETDAQ;
91 
92  typedef JRange<double> JRange_t;
93 
95  JLimit_t& numberOfEvents = inputFile.getLimit();
96  string outputFile;
97  string detectorFile;
98  JTDC_t TDC;
99  double laserFrequency_Hz;
100  bool overwriteDetector;
101  JROOTClassSelector selector;
102  string option;
103  double Wmin = 1e3; // Minimal signal intensity
104  JRange_t T_ns;
105  JRange_t X;
106  int debug;
107 
108  try {
109 
110  JProperties properties;
111 
112  properties.insert(gmake_property(Wmin));
113 
114  JParser<> zap("Application for dark room time calibration.");
115 
116  zap['f'] = make_field(inputFile, "input file (time slice data from laser calibration).");
117  zap['o'] = make_field(outputFile, "output file.") = "pulsar.root";
118  zap['a'] = make_field(detectorFile, "detector file.");
119  zap['!'] = make_field(TDC,
120  "Set reference PMTs, e.g."
121  "\n-! \"808969848 0 808982077 23\" sets PMT 0 of module 808969848 and PMT 23 of module 808982077 as references.");
122  zap['n'] = make_field(numberOfEvents) = JLimit::max();
123  zap['l'] = make_field(laserFrequency_Hz, "laser frequency [Hz]") = 10000;
124  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
125  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "LS";
126  zap['@'] = make_field(properties) = JPARSER::initialised();
127  zap['C'] = make_field(selector, "timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
128  zap['T'] = make_field(T_ns, "time window for time-over-threshold monitor") = JRange_t(-10.0, +10.0);
129  zap['x'] = make_field(X, "time window for histogram") = JRange_t();
130  zap['d'] = make_field(debug, "debug.") = 1;
131 
132  zap(argc, argv);
133  }
134  catch(const exception& error) {
135  FATAL(error.what() << endl);
136  }
137 
138 
139  TDC.is_valid(true);
140 
141  if (laserFrequency_Hz <= 0.0) {
142  FATAL("Invalid laser frequency " << laserFrequency_Hz << endl);
143  }
144 
145  const double laserPeriod_ns = 1.0e9 / laserFrequency_Hz;
146 
147  if (option.find('R') == string::npos) { option += 'R'; }
148  if (option.find('S') == string::npos) { option += 'S'; }
149  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
150 
151 
153 
154  try {
155  load(detectorFile, detector);
156  }
157  catch(const JException& error) {
158  FATAL(error);
159  }
160 
161  const JModuleRouter moduleRouter(detector);
162 
163 
164  // ROOT I/O
165 
166  TFile out(outputFile.c_str(), "recreate");
167 
168  putObject(out, JMeta(argc, argv));
169 
171 
172  map_type zmap;
173 
174  const double xmin = (X == JRange_t() ? -0.5 * laserPeriod_ns : X.getLowerLimit());
175  const double xmax = (X == JRange_t() ? +0.5 * laserPeriod_ns : X.getUpperLimit());
176  const int nx = 2 * (int) (xmax - xmin);
177 
178  TH1D h0("h0", NULL, nx, xmin, xmax);
179  TH1D h1("h1", NULL, 256, -0.5, +255.5);
180 
182 
183  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
184 
185  const JTDC_t::range_type range = TDC.equal_range(module->getID());
186 
187  for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
188 
189  NOTICE("Reference PMT at"
190  << " string " << setw(3) << module->getString()
191  << " floor " << setw(2) << module->getFloor()
192  << " module " << setw(8) << module->getID()
193  << " channel " << setw(2) << i->second << endl);
194 
195  const JPMTIdentifier id(module->getID(),i->second);
196 
197  ostringstream os;
198 
199  os << getLabel(module->getLocation()) << ' ' << getLabel(id);
200 
201  zmap.insert(make_pair(id, new TH1D(os.str().c_str(), NULL, nx, xmin, xmax)));
202  }
203  }
204 
205 
206  typedef vector<JHitR0> JDataL0_t;
207 
208  const JBuildL0<JHitR0> buildL0;
209 
210 
212 
213  counter_type counter = 0;
214 
215  for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
216 
217  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
218 
219  const JDAQTimeslice* timeslice = in.next();
220 
221  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
222 
223  const JTDC_t::range_type range = TDC.equal_range(frame->getModuleID());
224 
225  if (range.first != range.second) {
226 
227  const double t0 = get_time(getTimeOfFrame(frame->getFrameIndex()), laserPeriod_ns);
228 
229  JDataL0_t dataL0;
230 
231  buildL0(*frame, moduleRouter.getModule(frame->getModuleID()), back_inserter(dataL0));
232 
233  for (JDataL0_t::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
234 
235  const JPMTIdentifier id(frame->getModuleID(), hit->getPMT());
236 
237  map_type::const_iterator p = zmap.find(id);
238 
239  if (p != zmap.end()) {
240 
241  const double t1 = get_time(t0 + hit->getT(), laserPeriod_ns);
242 
243  p->second->Fill(t1);
244 
245  h0.Fill(t1);
246 
247  if (T_ns(t1)) {
248 
249  h1.Fill(hit->getToT());
250 
251  counts[id] += 1;
252  }
253  }
254  }
255  }
256  }
257  }
258  STATUS(endl);
259 
260 
261  map<JPMTIdentifier, double> t0; // fitted time offsets
262 
263 
264  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2.0*TMath::Pi())*[2]) + [3]");
265 
266  for (map_type::iterator it = zmap.begin(); it != zmap.end(); ++it) {
267 
268  JPMTIdentifier pmt = it->first;
269  TH1D* h1 = it->second;
270 
271  if (h1->GetEntries() == 0) {
272  WARNING("Histogram " << h1->GetName() << " empty" << endl);
273  continue;
274  }
275 
276  STATUS("--- PMT = " << pmt << "; histogram " << h1->GetName() << endl);
277 
278  // start values
279 
280  Double_t ymax = 0.0;
281  Double_t x0 = 0.0; // [ns]
282  Double_t sigma = 2.0; // [ns]
283  Double_t ymin = 0.0;
284 
285  for (int i = 1; i != h1->GetNbinsX(); ++i) {
286 
287  const Double_t x = h1->GetBinCenter (i);
288  const Double_t y = h1->GetBinContent(i);
289 
290  if (y > ymax) {
291  ymax = y;
292  x0 = x;
293  }
294  }
295 
296  f1.SetParameter(0, ymax/sigma);
297  f1.SetParameter(1, x0);
298  f1.SetParameter(2, sigma);
299  f1.SetParameter(3, ymin);
300 
301  for (Int_t i = 0; i != f1.GetNpar(); ++i) {
302  f1.SetParError(i, 0.0);
303  }
304 
305 
306  // fit
307 
308  TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", x0 - 3 * sigma, x0 + 3 * sigma);
309 
310  if (result.Get() == NULL) {
311  FATAL("Invalid TFitResultPtr " << h1->GetName() << endl);
312  }
313 
314  if (debug >= notice_t || !result->IsValid() || f1.GetParameter(0) < Wmin) {
315  cout << "Histogram " << h1->GetName() << " fit " << (result->IsValid() ? "ok" : "failed") << endl;
316  cout << "\tw = " << FIXED(12,3) << f1.GetParameter(0) << endl;
317  cout << "\tx0 = " << FIXED(12,3) << x0 << endl;
318  cout << "\tt0 = " << FIXED(12,3) << f1.GetParameter(1) << endl;
319  }
320 
321 
322  // check for multiple peaks
323 
324  int number_of_peaks = 0;
325 
326  Double_t dx = 2.0 * f1.GetParameter(2);
327  Double_t W = 0.0;
328  Double_t Y = f1.GetParameter(3);
329 
330  if (dx < (xmax - xmin) / nx) {
331  dx = (xmax - xmin) / nx; // minimal width
332  }
333 
334  for (Int_t il = 1, ir = il; ir <= nx; ) {
335 
336  for ( ; ir <= nx && h1->GetBinCenter(ir) <= h1->GetBinCenter(il) + dx; ++ir) {
337  W += h1->GetBinContent(ir) - Y;
338  }
339 
340  if (W >= Wmin) {
341 
342  number_of_peaks += 1;
343 
344  W = 0.0; // reset
345  il = ir;
346  ir = il;
347 
348  } else {
349 
350  W -= h1->GetBinContent(il) - Y; // shift
351  il += 1;
352  }
353  }
354 
355  if (number_of_peaks != 1) {
356  WARNING("Number of peaks " << h1->GetName() << ' ' << number_of_peaks << " != 1" << endl);
357  }
358 
359  if (result->IsValid() && f1.GetParameter(0) >= Wmin) {
360  t0[pmt] = f1.GetParameter(1);
361  }
362  }
363 
364  out.Write();
365  out.Close();
366 
367 
368  if (counter != 0) {
369 
370  const double W = laserFrequency_Hz * counter * getFrameTime() * 1.0e-9;
371 
372  NOTICE("Expection values [npe]" << endl);
373 
374  for (map<JPMTIdentifier, int>::const_iterator i = counts.begin(); i != counts.end(); ++i) {
375  NOTICE(i->first << ' ' << FIXED(7,3) << i->second / W << endl);
376  }
377  }
378 
379 
380  if (overwriteDetector) {
381 
382  int errors = 0;
383 
384  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
385 
386  const JTDC_t::range_type range = TDC.equal_range(module->getID());
387 
388  if (range.first != range.second) {
389 
390  const JPMTIdentifier id(module->getID(), range.first->second); // select t0 of the first reference PMT in this optical module
391 
393 
394  if (p0 != t0.end()) {
395 
396  const double t1 = p0->second;
397 
398  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
399 
400  map<JPMTIdentifier, double>::const_iterator p1 = t0.find(JPMTIdentifier(module->getID(),pmt));
401 
402  if (p1 != t0.end())
403  module->getPMT(pmt).subT0(p1->second); // offset with t0 of this reference PMT
404  else
405  module->getPMT(pmt).subT0(t1); // offset with t0 of first reference PMT in same optical module
406  }
407 
408  } else {
409 
410  if (!module->getPMT(id.getPMTAddress()).has(PMT_DISABLE)) {
411  ++errors;
412  }
413 
414  ERROR("Module/PMT "
415  << setw(10) << module->getID() << "/" << FILL(2,'0') << id.getPMTAddress() << FILL() << ' '
416  << "missing or insufficient signal." << endl);
417  }
418  }
419  }
420 
421  if (errors == 0) {
422 
423  NOTICE("Store calibration data on file " << detectorFile << endl);
424 
425  detector.comment.add(JMeta(argc, argv));
426 
427  store(detectorFile, detector);
428 
429  } else {
430 
431  FATAL("Number of errors " << errors << " != 0" << endl);
432  }
433  }
434 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
#define WARNING(A)
Definition: JMessage.hh:65
debug
Definition: JMessage.hh:29
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
TPaveText * p1
Basic data structure for L0 hit.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
#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
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Auxiliary class to select ROOT class based on class name.
Router for direct addressing of module data in detector data structure.
Auxiliary class for TDC constraints.
Definition: JTDC_t.hh:37
notice
Definition: JMessage.hh:32
Utility class to parse parameter values.
Definition: JProperties.hh:497
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:84
Long64_t counter_type
Type definition for counter.
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
Auxiliary class for multiplexing object iterators.
Basic data structure for time and time over threshold information of hit.
string outputFile
Data structure for detector geometry and calibration.
Utility class to parse parameter values.
double getTimeOfFrame(const int frame_index)
Get start time of frame in ns since start of run for a given frame index.
Definition: JDAQClock.hh:185
const double sigma[]
Definition: JQuadrature.cc:74
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Type definition of range.
Definition: JHead.hh:41
Template specialisation of L0 builder for JHitR0 data type.
Definition: JBuildL0.hh:174
Detector file.
Definition: JHead.hh:226
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
do set_variable OUTPUT_DIRECTORY $WORKDIR T
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
static const int PMT_DISABLE
KM3NeT Data Definitions v3.3.0-2-g5cc95cf https://git.km3net.de/common/km3net-dataformat.
Definition: pmt_status.hh:12
ROOT I/O of application specific meta data.
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
Data time slice.
virtual const pointer_type & next() override
Get next element.
virtual bool hasNext() override
Check availability of next element.
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
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to module in detector data structure.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
const double xmin
Definition: JQuadrature.cc:23
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
General purpose class for object reading from a list of file names.
then fatal The output file must have the wildcard in the e g root fi 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:48
Utility class to parse command line options.
set_variable TDC
Definition: JPrintTDC.sh:20
no fit printf nominal n $STRING awk v X
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
std::map< int, range_type > map_type
do set_variable DETECTOR_TXT $WORKDIR detector
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62