Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JPulsar.cc File Reference

Application for dark room time calibration using external laser. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "TFitResult.h"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTIdentifier.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JCalibrate/JTDC_t.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JTools/JRange.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "Jeep/JProperties.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

Application for dark room time calibration using external laser.

The time calibration is made for the specified reference PMTs (option -! "\<module identifier\> \<PMT channel\>".
In this, the wild card value -1 can be used. A ROOT histogram corresponding to each reference PMT is created, filled, fitted and stored to the specified output file.
The detector can be updated using option -A
The laser should be synchronised with the same clock system as the detector and produce single pulses at a fixed frequency (option -l <laserFrequency_Hz>.

In addition, the transit-time distribution of the reference PMTs can be monitored (see JLegolas.cc).
For this, the expectation values should be less than one and the time-over-threshold distribution nominal.

Author
mbouwhuis

Definition in file JPulsar.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

< Minimal signal intensity

Definition at line 84 of file JPulsar.cc.

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  int debug;
104 
105  try {
106 
107  JProperties properties;
108 
109  properties.insert(gmake_property(Wmin));
110 
111  JParser<> zap("Application for dark room time calibration.");
112 
113  zap['f'] = make_field(inputFile, "input file (time slice data from laser calibration).");
114  zap['o'] = make_field(outputFile, "output file.") = "pulsar.root";
115  zap['a'] = make_field(detectorFile, "detector file.");
116  zap['!'] = make_field(TDC,
117  "Set reference PMTs, e.g."
118  "\n-! \"808969848 0 808982077 23\" sets PMT 0 of module 808969848 and PMT 23 of module 808982077 as references.");
119  zap['n'] = make_field(numberOfEvents) = JLimit::max();
120  zap['l'] = make_field(laserFrequency_Hz, "laser frequency [Hz]") = 10000;
121  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
122  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "LS";
123  zap['@'] = make_field(properties) = JPARSER::initialised();
124  zap['C'] = make_field(selector, "timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
125  zap['T'] = make_field(T_ns) = JRange_t(-10.0, +10.0);
126  zap['d'] = make_field(debug, "debug.") = 1;
127 
128  zap(argc, argv);
129  }
130  catch(const exception& error) {
131  FATAL(error.what() << endl);
132  }
133 
134  cout.tie(&cerr);
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 int nx = 2 * (int) laserPeriod_ns;
173  const double xmin = -0.5 * laserPeriod_ns;
174  const double xmax = 0.5 * laserPeriod_ns;
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 
305  if (debug >= notice_t || !result->IsValid() || f1.GetParameter(0) < Wmin) {
306  cout << "Histogram " << h1->GetName() << " fit " << (result->IsValid() ? "ok" : "failed") << endl;
307  cout << "\tw = " << FIXED(12,3) << f1.GetParameter(0) << endl;
308  cout << "\tx0 = " << FIXED(12,3) << x0 << endl;
309  cout << "\tt0 = " << FIXED(12,3) << f1.GetParameter(1) << endl;
310  }
311 
312 
313  // check for multiple peaks
314 
315  int number_of_peaks = 0;
316 
317  Double_t dx = 2.0 * f1.GetParameter(2);
318  Double_t W = 0.0;
319  Double_t Y = f1.GetParameter(3);
320 
321  if (dx < (xmax - xmin) / nx) {
322  dx = (xmax - xmin) / nx; // minimal width
323  }
324 
325  for (Int_t il = 1, ir = il; ir <= nx; ) {
326 
327  for ( ; ir <= nx && h1->GetBinCenter(ir) <= h1->GetBinCenter(il) + dx; ++ir) {
328  W += h1->GetBinContent(ir) - Y;
329  }
330 
331  if (W >= Wmin) {
332 
333  number_of_peaks += 1;
334 
335  W = 0.0; // reset
336  il = ir;
337  ir = il;
338 
339  } else {
340 
341  W -= h1->GetBinContent(il) - Y; // shift
342  il += 1;
343  }
344  }
345 
346  if (number_of_peaks != 1) {
347  WARNING("Number of peaks " << h1->GetName() << ' ' << number_of_peaks << " != 1" << endl);
348  }
349 
350  if (result->IsValid() && f1.GetParameter(0) >= Wmin) {
351  t0[pmt] = f1.GetParameter(1);
352  }
353  }
354 
355  out.Write();
356  out.Close();
357 
358 
359  if (counter != 0) {
360 
361  const double W = laserFrequency_Hz * counter * getFrameTime() * 1.0e-9;
362 
363  NOTICE("Expection values [npe]" << endl);
364 
365  for (map<JPMTIdentifier, int>::const_iterator i = counts.begin(); i != counts.end(); ++i) {
366  NOTICE(i->first << ' ' << FIXED(7,3) << i->second / W << endl);
367  }
368  }
369 
370 
371  if (overwriteDetector) {
372 
373  int errors = 0;
374 
375  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
376 
377  const JTDC_t::range_type range = TDC.equal_range(module->getID());
378 
379  if (range.first != range.second) {
380 
381  const JPMTIdentifier id(module->getID(), range.first->second); // select t0 of the first reference PMT in this optical module
382 
384 
385  if (p0 != t0.end()) {
386 
387  const double t1 = p0->second;
388 
389  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
390 
391  map<JPMTIdentifier, double>::const_iterator p1 = t0.find(JPMTIdentifier(module->getID(),pmt));
392 
393  if (p1 != t0.end())
394  module->getPMT(pmt).subT0(p1->second); // offset with t0 of this reference PMT
395  else
396  module->getPMT(pmt).subT0(t1); // offset with t0 of first reference PMT in same optical module
397  }
398 
399  } else {
400 
401  ++errors;
402 
403  ERROR("Module/PMT "
404  << setw(10) << module->getID() << "/" << FILL(2,'0') << id.getPMTAddress() << FILL() << ' '
405  << "missing or insufficient signal." << endl);
406  }
407  }
408  }
409 
410  if (errors == 0) {
411 
412  NOTICE("Store calibration data on file " << detectorFile << endl);
413 
414  detector.comment.add(JMeta(argc, argv));
415 
416  store(detectorFile, detector);
417 
418  } else {
419 
420  FATAL("Number of errors " << errors << " != 0" << endl);
421  }
422  }
423 }
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:1500
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
debug
Definition: JMessage.hh:29
TPaveText * p1
#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:80
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:34
notice
Definition: JMessage.hh:32
Utility class to parse parameter values.
Definition: JProperties.hh:496
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
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:445
Auxiliary class for multiplexing object iterators.
string outputFile
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
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Template specialisation of L0 builder for JHitR0 data type.
Definition: JBuildL0.hh:174
Detector file.
Definition: JHead.hh:196
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
return result
Definition: JPolint.hh:727
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
Data time slice.
int debug
debug level
Definition: JSirene.cc:63
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:327
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi mv $WORKDIR/fit.root $MODULE_ROOT typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
Definition: module-Z:fit.sh:84
#define FATAL(A)
Definition: JMessage.hh:67
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
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 source JAcoustics sh $DETECTOR_ID typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:36
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
bool putObject(TDirectory *dir, const T &object)
Write object to ROOT directory.