Jpp  17.2.1-pre0
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  cout.tie(&cerr);
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 
170  typedef map<JPMTIdentifier, TH1D*> map_type;
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 
302  // fit
303 
304  TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", x0 - 3 * sigma, x0 + 3 * sigma);
305 
306  if (result.Get() == NULL) {
307  FATAL("Invalid TFitResultPtr " << h1->GetName() << endl);
308  }
309 
310  if (debug >= notice_t || !result->IsValid() || f1.GetParameter(0) < Wmin) {
311  cout << "Histogram " << h1->GetName() << " fit " << (result->IsValid() ? "ok" : "failed") << endl;
312  cout << "\tw = " << FIXED(12,3) << f1.GetParameter(0) << endl;
313  cout << "\tx0 = " << FIXED(12,3) << x0 << endl;
314  cout << "\tt0 = " << FIXED(12,3) << f1.GetParameter(1) << endl;
315  }
316 
317 
318  // check for multiple peaks
319 
320  int number_of_peaks = 0;
321 
322  Double_t dx = 2.0 * f1.GetParameter(2);
323  Double_t W = 0.0;
324  Double_t Y = f1.GetParameter(3);
325 
326  if (dx < (xmax - xmin) / nx) {
327  dx = (xmax - xmin) / nx; // minimal width
328  }
329 
330  for (Int_t il = 1, ir = il; ir <= nx; ) {
331 
332  for ( ; ir <= nx && h1->GetBinCenter(ir) <= h1->GetBinCenter(il) + dx; ++ir) {
333  W += h1->GetBinContent(ir) - Y;
334  }
335 
336  if (W >= Wmin) {
337 
338  number_of_peaks += 1;
339 
340  W = 0.0; // reset
341  il = ir;
342  ir = il;
343 
344  } else {
345 
346  W -= h1->GetBinContent(il) - Y; // shift
347  il += 1;
348  }
349  }
350 
351  if (number_of_peaks != 1) {
352  WARNING("Number of peaks " << h1->GetName() << ' ' << number_of_peaks << " != 1" << endl);
353  }
354 
355  if (result->IsValid() && f1.GetParameter(0) >= Wmin) {
356  t0[pmt] = f1.GetParameter(1);
357  }
358  }
359 
360  out.Write();
361  out.Close();
362 
363 
364  if (counter != 0) {
365 
366  const double W = laserFrequency_Hz * counter * getFrameTime() * 1.0e-9;
367 
368  NOTICE("Expection values [npe]" << endl);
369 
370  for (map<JPMTIdentifier, int>::const_iterator i = counts.begin(); i != counts.end(); ++i) {
371  NOTICE(i->first << ' ' << FIXED(7,3) << i->second / W << endl);
372  }
373  }
374 
375 
376  if (overwriteDetector) {
377 
378  int errors = 0;
379 
380  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
381 
382  const JTDC_t::range_type range = TDC.equal_range(module->getID());
383 
384  if (range.first != range.second) {
385 
386  const JPMTIdentifier id(module->getID(), range.first->second); // select t0 of the first reference PMT in this optical module
387 
389 
390  if (p0 != t0.end()) {
391 
392  const double t1 = p0->second;
393 
394  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
395 
396  map<JPMTIdentifier, double>::const_iterator p1 = t0.find(JPMTIdentifier(module->getID(),pmt));
397 
398  if (p1 != t0.end())
399  module->getPMT(pmt).subT0(p1->second); // offset with t0 of this reference PMT
400  else
401  module->getPMT(pmt).subT0(t1); // offset with t0 of first reference PMT in same optical module
402  }
403 
404  } else {
405 
406  if (!module->getPMT(id.getPMTAddress()).has(PMT_DISABLE)) {
407  ++errors;
408  }
409 
410  ERROR("Module/PMT "
411  << setw(10) << module->getID() << "/" << FILL(2,'0') << id.getPMTAddress() << FILL() << ' '
412  << "missing or insufficient signal." << endl);
413  }
414  }
415  }
416 
417  if (errors == 0) {
418 
419  NOTICE("Store calibration data on file " << detectorFile << endl);
420 
421  detector.comment.add(JMeta(argc, argv));
422 
423  store(detectorFile, detector);
424 
425  } else {
426 
427  FATAL("Number of errors " << errors << " != 0" << endl);
428  }
429  }
430 }
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
#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)
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-1-gb6f251e 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
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62