62int main(
int argc,
char **argv)
69 JLimit_t& numberOfEvents = inputFile.getLimit();
74 int numberOfTimeslices;
85 JParser<> zap(
"Example program to search for correlations between triggered events and timeslice data.");
93 zap[
'N'] =
make_field(numberOfTimeslices) = 100;
102 zap.
read(argc, argv);
104 catch(
const exception& error) {
105 FATAL(error.what() << endl);
126 vector <hit_type> data;
130 const Double_t xmin = -(numberOfTimeslices + 1) *
getFrameTime() - 0.5*binWidth_ns;
131 const Double_t xmax = +(numberOfTimeslices + 1) *
getFrameTime() + 0.5*binWidth_ns;
132 const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
134 JManager_t manager(
new TH1D(
"M_%", NULL, nx, xmin, xmax));
142 if (selector ==
"") {
150 FATAL(
"No timeslice data." << endl);
153 NOTICE(
"Selected class " << ps->getClassName() << endl);
159 ps->configure(inputFile);
180 ps->setLimit(inputFile.getLimit());
187 map_type peaksPerDoms;
195 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
196 oosDomsId[module->getID()] =
false;
202 STATUS(
"event: " << setw(10) << in.getCounter() <<
'\r');
DEBUG(endl);
209 bool test_OOS =
false;
214 if(oosDomsId[hit->getModuleID()]) {
220 if(test_OOS)
continue;
225 buffer[
event->getFrameIndex()].push_back(t0);
232 while (ps->hasNext()) {
234 STATUS(
"event: " << setw(10) << ps->getCounter() <<
'\r');
DEBUG(endl);
238 map_type::const_iterator p = buffer.lower_bound(timeslice->
getFrameIndex() - numberOfTimeslices);
239 map_type::const_iterator q = buffer.upper_bound(timeslice->
getFrameIndex() + numberOfTimeslices);
241 int number_of_events = 0;
243 for (map_type::const_iterator i = p; i != q; ++i) {
244 number_of_events += i->second.size();
247 if (number_of_events == 0) {
251 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
255 buildL1(*frame, router.
getModule(frame->getModuleID()), back_inserter(data));
257 TH1D* h1 = manager[frame->getModuleID()];
259 for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
261 const double t1 = *hit + frame->getFrameIndex() *
getFrameTime();
263 for (map_type::const_iterator i = p; i != q; ++i) {
264 for (map_type::mapped_type::const_iterator
j = i->second.begin();
j != i->second.end(); ++
j) {
266 const double t0 = *
j;
279 Double_t most_OOS_peak = 0;
280 Int_t most_OOS_ID = 0;
282 TF1 f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
290 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
294 status2[module->getID()] = DEFAULT;
296 JManager_t::iterator p = manager.find(module->getID());
298 if (p == manager.end() || p->second->GetEntries() == 0) {
300 status2[module->getID()] = NODATA;
305 TH1D* h1 = p->second;
311 Double_t sigma = 250.0;
314 for (
int i = 1; i != h1->GetNbinsX(); ++i) {
316 const Double_t x = h1->GetBinCenter (i);
317 const Double_t y = h1->GetBinContent(i);
327 ymin /= h1->GetNbinsX();
329 f1.SetParameter(0, ymax);
330 f1.SetParameter(1, x0);
331 if (binWidth_ns < sigma)
332 f1.SetParameter(2, sigma);
334 f1.FixParameter(2, binWidth_ns/sqrt(12.0));
335 f1.SetParameter(3, ymin);
337 for (Int_t i = 0; i != f1.GetNpar(); ++i) {
338 f1.SetParError(i, 0.0);
343 h1->Fit(&f1, option.c_str(),
"same", x0 - 5 * sigma, x0 + 5 * sigma);
345 status = (fabs(f1.GetParameter(1)) <= 0.5*
getFrameTime() &&
346 f1.GetParameter(0) >= f1.GetParameter(3));
348 if (status) status2[module->getID()] = IN_SYNC;
355 for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
357 const Double_t x = h1->GetBinCenter (i);
358 const Double_t y = h1->GetBinContent(i);
360 while (x > (ns + 1) *
getFrameTime() - TMax_ns) { ++ns; }
368 DEBUG(
"Module " << setw(8) << module->getID() << endl);
372 const Double_t y = bg.getTotal() * i->second.getCount() / bg.getCount();
373 const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
375 DEBUG(
"Peak at " << setw(4) << i->first <<
" [frame time] "
377 <<
FIXED(7,1) << i->second.getTotal() <<
" / "
378 <<
FIXED(7,1) << y <<
' '
381 if (P < Pmin && i->second.getTotal() > y)
387 if (!(sn.size() == 1 &&
388 sn.begin()->first == 0)) {
390 WARNING(
"Module " << setw(8) << module->getID() <<
" number of peaks " << sn.size() << endl);
394 const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
396 WARNING(
"Peak at " << setw(4) << i->first <<
" [frame time] "
398 <<
FIXED(7,1) << i->second.getTotal() <<
" / "
399 <<
FIXED(7,1) << noise << endl);
401 peaksPerDoms[module->getID()].push_back(sn.size());
402 peaksPerDoms[module->getID()].push_back(i->first);
403 peaksPerDoms[module->getID()].push_back(i->second.getTotal());
404 peaksPerDoms[module->getID()].push_back(noise);
407 status2[module->getID()] = (sn.size() == 1 ? OUT_SYNC :
ERROR);
414 if ((f1.GetParameter(0) > most_OOS_peak) && !(oosDomsId[module->getID()])){
416 most_OOS_peak = f1.GetParameter(0);
417 most_OOS_ID =
module->getID();
422 if (most_OOS_ID != 0) {
423 oosDomsId[most_OOS_ID] =
true;
426 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
427 manager[module->getID()]->Reset(
"ICESM");
429 peaksPerDoms.clear();
437 NOTICE(
"Managing non-working and OOS DOMs." << endl);
443 parameters.
load(pmtFile.c_str());
449 if (i->second != IN_SYNC) {
451 NOTICE(
"Module " << setw(8) << i->first <<
" set QEs of all PMTs to zero." << endl);
459 ofstream out(pmtFile.c_str());
463 out << parameters << endl;
474 ofstream stream(peaks);
476 stream <<
"#DOM NB_PEAKS TIMESLICE_SHIFT SIGNAL NOISE\n";
478 for (map_type::const_iterator dom = peaksPerDoms.begin(); dom != peaksPerDoms.end(); ++dom) {
479 stream << dom->first <<
" ";
481 for (map_type::mapped_type::const_iterator data = dom->second.begin(); data != dom->second.end(); ++data) {
483 ((i%4 == 0) ? stream <<
"\n" : stream <<
" ");
495 if (i->second == IN_SYNC) {
498 if (i->second != IN_SYNC &&
499 i->second != NODATA) {
504 NOTICE(
"Number of modules out/in sync " << nout <<
'/' << nin << endl);
506 QAQC(nin <<
' ' << nout << endl);