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