88int main(
int argc,
char **argv)
95 JLimit_t& numberOfEvents = inputFile.getLimit();
98 bool overwriteDetector;
100 int numberOfTimeslices;
110 JParser<> zap(
"Example program to search for correlations between triggered events and timeslice data.");
112 zap[
'f'] =
make_field(inputFile,
"input file.");
114 zap[
'a'] =
make_field(detectorFile,
"detector file.");
115 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with module (PMT) status.");
117 zap[
'T'] =
make_field(TMaxLocal_ns,
"time window for local coincidences (L1),") = 20.0;
118 zap[
'N'] =
make_field(numberOfTimeslices,
"time slice difference between triggered event and L1.") = 100;
119 zap[
'W'] =
make_field(binWidth_ns,
"bin width for output histograms." ) = 10e3;
120 zap[
'D'] =
make_field(deadTime_us,
"L1 dead time (us)") = 200.0;
121 zap[
'p'] =
make_field(Pmin,
"minimal probability for background to be signal.") = 1.0e-7;
128 catch(
const exception& error) {
129 FATAL(error.what() << endl);
143 FATAL(
"Empty detector " << detectorFile << endl);
149 const double BOOST = 20.0;
150 const double deadTime_ns = deadTime_us * 1e3;
152 NOTICE(
"Time window " <<
FIXED(7,1) << TMax_ns <<
" [ns]" << endl);
158 vector <hit_type> data;
162 const Double_t xmin = -(numberOfTimeslices + 1) *
getFrameTime() - 0.5*binWidth_ns;
163 const Double_t xmax = +(numberOfTimeslices + 1) *
getFrameTime() + 0.5*binWidth_ns;
164 const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
166 JManager_t manager(
new TH1D(
"M_%", NULL, nx, xmin, xmax));
174 if (selector ==
"") {
182 FATAL(
"No timeslice data." << endl);
185 NOTICE(
"Selected class " << ps->getClassName() << endl);
191 ps->configure(inputFile);
194 ps->setLimit(inputFile.getLimit());
202 if (in.getCounter()%1000 == 0) {
203 STATUS(
"event: " << setw(10) << in.getCounter() <<
'\r');
DEBUG(endl);
214 if (router.
hasModule(hit->getModuleID())) {
223 buffer[
event->getFrameIndex()].push_back(t0);
228 while (ps->hasNext()) {
230 if (ps->getCounter()%1000 == 0) {
231 STATUS(
"event: " << setw(10) << ps->getCounter() <<
'\r');
DEBUG(endl);
236 map_type::const_iterator p = buffer.lower_bound(timeslice->
getFrameIndex() - numberOfTimeslices);
237 map_type::const_iterator q = buffer.upper_bound(timeslice->
getFrameIndex() + numberOfTimeslices);
239 int number_of_events = 0;
241 for (map_type::const_iterator i = p; i != q; ++i) {
242 number_of_events += i->second.size();
245 if (number_of_events == 0) {
249 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
251 if (router.
hasModule(frame->getModuleID())) {
255 buildL1(*frame, router.
getModule(frame->getModuleID()), back_inserter(data));
257 TH1D* h1 = manager[frame->getModuleID()];
259 double t1 = numeric_limits<double>::lowest();
261 for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
263 const double t2 = *hit + frame->getFrameIndex() *
getFrameTime();
265 if (t2 > t1 + deadTime_ns) {
267 for (map_type::const_iterator i = p; i != q; ++i) {
268 for (map_type::mapped_type::const_iterator t0 = i->second.begin(); t0 != i->second.end(); ++t0) {
284 TF1 f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
295 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
297 if (module->getFloor() != 0) {
299 status[module->getID()] = DEFAULT;
301 JManager_t::iterator p = manager.find(module->getID());
303 if (p == manager.end() || p->second->GetEntries() == 0) {
305 status[module->getID()] = NODATA;
310 TH1D* h1 = p->second;
316 Double_t sigma = 250.0;
319 for (
int i = 1; i != h1->GetNbinsX(); ++i) {
321 const Double_t x = h1->GetBinCenter (i);
322 const Double_t y = h1->GetBinContent(i);
332 ymin /= h1->GetNbinsX();
334 f1.SetParameter(0, ymax);
335 f1.SetParameter(1, x0);
336 if (binWidth_ns < sigma)
337 f1.SetParameter(2, sigma);
339 f1.FixParameter(2, binWidth_ns/sqrt(12.0));
340 f1.SetParameter(3, ymin);
342 for (Int_t i = 0; i != f1.GetNpar(); ++i) {
343 f1.SetParError(i, 0.0);
348 h1->Fit(&f1, option.c_str(),
"same", x0 - 5 * sigma, x0 + 5 * sigma);
350 if (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() &&
351 f1.GetParameter(0) >= f1.GetParameter(3)) {
353 status[module->getID()] = IN_SYNC;
356 double fm = fmod(f1.GetParameter(1), getFrameTime());
362 << setw(10) << module->getID() <<
' '
363 <<
FIXED(15,3) << f1.GetParameter(1) <<
' '
364 <<
FIXED(12,3) << f1.GetParameter(0) <<
' '
365 <<
FIXED(12,3) << f1.GetParameter(3) <<
' '
366 <<
FIXED(12,3) << fm <<
' '
367 << status[module->getID()] << endl);
374 for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
376 const Double_t x = h1->GetBinCenter (i);
377 const Double_t y = h1->GetBinContent(i);
379 while (x > (ns + 1) *
getFrameTime() - BOOST * TMax_ns) {
385 else if (fabs(x - ns *
getFrameTime()) < BOOST * TMax_ns)
391 const Double_t y = getBackground(i->second, bg[i->first]);
392 const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
394 if (
debug >=
debug_t || status[module->getID()] != IN_SYNC) {
395 cout <<
"Module/peak " << setw(10) << module->getID() <<
' '
396 << setw(4) << i->first <<
' '
403 <<
FIXED(7,1) << i->second.getTotal() <<
" / "
404 <<
FIXED(7,1) << y <<
' '
406 << (i->second.getTotal() > y && P < Pmin && i->first != 0 ?
"***" :
"") << endl;
409 if (i->second.getTotal() > y && P < Pmin)
415 if (!(sn.size() == 1 &&
416 sn.begin()->first == 0)) {
418 status[module->getID()] = (sn.size() == 1 ? OUT_SYNC :
ERROR);
420 ERROR(
"Module/error "
421 << setw(10) << module->getID() <<
' '
422 <<
"number of peaks " << sn.size() <<
' '
423 <<
"peak " << (sn.size() == 1 ?
MAKE_CSTRING(sn.begin()->first) :
"?") <<
' '
424 << status[module->getID()] << endl);
429 if (overwriteDetector) {
435 if (i->second != IN_SYNC &&
436 i->second != NODATA) {
438 NOTICE(
"Module " << setw(8) << i->first <<
" set out-of-sync." << endl);
442 module.getStatus().set(MODULE_OUT_OF_SYNC);
444 for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
461 if (i->second == IN_SYNC) {
464 if (i->second != IN_SYNC &&
465 i->second != NODATA) {
470 NOTICE(
"Number of modules out/in sync " << nout <<
'/' << nin << endl);
472 QAQC(nin <<
' ' << nout << endl);