86int main(
int argc,
char **argv)
95 JLimit_t& numberOfEvents = inputFile.getLimit();
99 double laserFrequency_Hz;
100 bool overwriteDetector;
114 JParser<> zap(
"Application for dark room time calibration.");
116 zap[
'f'] =
make_field(inputFile,
"input file (time slice data from laser calibration).");
118 zap[
'a'] =
make_field(detectorFile,
"detector file.");
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.");
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";
128 zap[
'T'] =
make_field(T_ns,
"time window for time-over-threshold monitor") =
JRange_t(-10.0, +10.0);
134 catch(
const exception& error) {
135 FATAL(error.what() << endl);
141 if (laserFrequency_Hz <= 0.0) {
142 FATAL(
"Invalid laser frequency " << laserFrequency_Hz << endl);
145 const double laserPeriod_ns = 1.0e9 / laserFrequency_Hz;
147 if (option.find(
'R') == string::npos) { option +=
'R'; }
148 if (option.find(
'S') == string::npos) { option +=
'S'; }
176 const int nx = 2 * (int) (xmax - xmin);
178 TH1D h0(
"h0", NULL, nx, xmin, xmax);
179 TH1D h1(
"h1", NULL, 256, -0.5, +255.5);
183 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
187 for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
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);
201 zmap.insert(make_pair(
id,
new TH1D(os.str().c_str(), NULL, nx, xmin, xmax)));
215 for ( ; in.
hasNext() && counter != inputFile.getLimit(); ++counter) {
217 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
221 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
225 if (range.first != range.second) {
227 const double t0 = get_time(
getTimeOfFrame(frame->getFrameIndex()), laserPeriod_ns);
231 buildL0(*frame, moduleRouter.
getModule(frame->getModuleID()), back_inserter(dataL0));
233 for (JDataL0_t::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
237 map_type::const_iterator p = zmap.find(
id);
239 if (p != zmap.end()) {
241 const double t1 = get_time(t0 + hit->getT(), laserPeriod_ns);
249 h1.Fill(hit->getToT());
264 TF1 f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2.0*TMath::Pi())*[2]) + [3]");
266 for (map_type::iterator it = zmap.begin(); it != zmap.end(); ++it) {
269 TH1D* h1 = it->second;
271 if (h1->GetEntries() == 0) {
272 WARNING(
"Histogram " << h1->GetName() <<
" empty" << endl);
276 STATUS(
"--- PMT = " << pmt <<
"; histogram " << h1->GetName() << endl);
282 Double_t sigma = 2.0;
285 for (
int i = 1; i != h1->GetNbinsX(); ++i) {
287 const Double_t x = h1->GetBinCenter (i);
288 const Double_t y = h1->GetBinContent(i);
296 f1.SetParameter(0, ymax/sigma);
297 f1.SetParameter(1, x0);
298 f1.SetParameter(2, sigma);
299 f1.SetParameter(3, ymin);
301 for (Int_t i = 0; i != f1.GetNpar(); ++i) {
302 f1.SetParError(i, 0.0);
308 TFitResultPtr
result = h1->Fit(&f1, option.c_str(),
"same", x0 - 3 * sigma, x0 + 3 * sigma);
310 if (
result.Get() == NULL) {
311 FATAL(
"Invalid TFitResultPtr " << h1->GetName() << endl);
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;
324 int number_of_peaks = 0;
326 Double_t dx = 2.0 * f1.GetParameter(2);
328 Double_t Y = f1.GetParameter(3);
330 if (dx < (xmax - xmin) / nx) {
331 dx = (xmax - xmin) / nx;
334 for (Int_t il = 1, ir = il; ir <= nx; ) {
336 for ( ; ir <= nx && h1->GetBinCenter(ir) <= h1->GetBinCenter(il) + dx; ++ir) {
337 W += h1->GetBinContent(ir) - Y;
342 number_of_peaks += 1;
350 W -= h1->GetBinContent(il) - Y;
355 if (number_of_peaks != 1) {
356 WARNING(
"Number of peaks " << h1->GetName() <<
' ' << number_of_peaks <<
" != 1" << endl);
359 if (
result->IsValid() && f1.GetParameter(0) >= Wmin) {
360 t0[pmt] = f1.GetParameter(1);
370 const double W = laserFrequency_Hz * counter *
getFrameTime() * 1.0e-9;
372 NOTICE(
"Expection values [npe]" << endl);
375 NOTICE(i->first <<
' ' <<
FIXED(7,3) << i->second / W << endl);
380 if (overwriteDetector) {
384 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
388 if (range.first != range.second) {
394 if (p0 != t0.end()) {
396 const double t1 = p0->second;
403 module->getPMT(pmt).subT0(p1->second);
405 module->getPMT(pmt).subT0(t1);
410 if (!module->getPMT(
id.getPMTAddress()).has(
PMT_DISABLE)) {
415 << setw(10) << module->getID() <<
"/" <<
FILL(2,
'0') <<
id.getPMTAddress() <<
FILL() <<
' '
416 <<
"missing or insufficient signal." << endl);
423 NOTICE(
"Store calibration data on file " << detectorFile << endl);
431 FATAL(
"Number of errors " << errors <<
" != 0" << endl);