52 using namespace KM3NETDAQ;
60 int numberOfTimeslices;
64 JROOTClassSelector selector;
71 JParser<> zap(
"Example program to search for correlations between triggered events and timeslice data.");
76 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
79 zap[
'N'] =
make_field(numberOfTimeslices) = 100;
83 zap[
'C'] =
make_field(selector) =
"", getROOTClassSelection<JDAQTimesliceTypes_t>();
88 if (zap.read(argc, argv) != 0)
91 catch(
const exception& error) {
92 FATAL(error.what() << endl);
100 const int WR = 0x80000000;
107 JBit(i->second).
set(MASK[i->first]);
111 DEBUG(
"Mask " << setw(8) << i->first <<
" 0x" << hex << i->second << dec << endl);
117 load(detectorFile, detector);
119 catch(
const JException& error) {
123 const JDAQHitRouter router(detector);
125 const double TMax_ns = max(binWidth_ns,
getMaximalTime(detector));
129 JBuildL0<hit_type> buildL0;
130 JBuildL1<hit_type> buildL1(TMaxLocal_ns,
true);
135 const Double_t xmin = -(numberOfTimeslices + 1) *
getFrameTime() - 0.5*binWidth_ns;
136 const Double_t xmax = +(numberOfTimeslices + 1) *
getFrameTime() + 0.5*binWidth_ns;
137 const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
139 JManager_t manager(
new TH1D(
"M_%", NULL, nx, xmin, xmax));
143 JSinglePointer< JTreeScannerInterface<JDAQTimeslice> > ps;
147 if (selector ==
"") {
149 if ((ps =
new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
150 (ps =
new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
151 (ps =
new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
152 (ps =
new JTreeScanner< JAssertConversion<JDAQTimeslice, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
153 (ps =
new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
155 FATAL(
"No timeslice data." << endl);
158 NOTICE(
"Selected class " << ps->getClassName() << endl);
164 ps->configure(inputFile);
191 JPMTParametersMap parameters;
198 parameters.load(pmtFile.c_str());
200 catch(
const JException& error) {}
204 for (JTreeScanner<JDAQEvent> in(inputFile.getFilename(), inputFile.getLimit()); in.hasNext(); ) {
206 STATUS(
"event: " << setw(10) << in.getCounter() <<
'\r');
DEBUG(endl);
213 bool test_OOS =
false;
216 t0 +=
getTime(*hit, router.getPMT(*hit));
218 if(std::find(oosDomsId.begin(), oosDomsId.end(), hit->getModuleIdentifier()) != oosDomsId.end()) {
224 if(test_OOS)
continue;
229 buffer[
event->getFrameIndex()].push_back(t0);
236 while (ps->hasNext()) {
238 STATUS(
"event: " << setw(10) << ps->getCounter() <<
'\r');
DEBUG(endl);
242 map_type::const_iterator p = buffer.lower_bound(timeslice->
getFrameIndex() - numberOfTimeslices);
243 map_type::const_iterator q = buffer.upper_bound(timeslice->
getFrameIndex() + numberOfTimeslices);
245 int number_of_events = 0;
247 for (map_type::const_iterator i = p; i != q; ++i) {
248 number_of_events += i->second.size();
251 if (number_of_events == 0) {
255 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
259 if ((frame->getStatus() & ~MASK[frame->getModuleID()] & ~WR) == 0) {
263 buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
265 TH1D* h1 = manager[frame->getModuleID()];
269 const double t1 = *hit + frame->getFrameIndex() *
getFrameTime();
271 for (map_type::const_iterator i = p; i != q; ++i) {
272 for (map_type::mapped_type::const_iterator j = i->second.begin(); j != i->second.end(); ++j) {
274 const double t0 = *j;
290 Double_t most_OOS_peak = 0;
291 Int_t most_OOS_ID = 0;
293 TF1 f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
301 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
305 JManager_t::iterator p = manager.find(module->getID());
307 if (p != manager.end()) {
309 TH1D* h1 = p->second;
311 if (h1->GetEntries() != 0) {
317 Double_t sigma = 250.0;
320 for (
int i = 1; i != h1->GetNbinsX(); ++i) {
322 const Double_t x = h1->GetBinCenter (i);
323 const Double_t y = h1->GetBinContent(i);
333 ymin /= h1->GetNbinsX();
335 f1.SetParameter(0, ymax);
336 f1.SetParameter(1, x0);
337 if (binWidth_ns < sigma)
338 f1.SetParameter(2, sigma);
340 f1.FixParameter(2, binWidth_ns/sqrt(12.0));
341 f1.SetParameter(3, ymin);
345 h1->Fit(&f1, option.c_str(),
"same", x0 - 5 * sigma, x0 + 5 * sigma);
347 status = (fabs(f1.GetParameter(1)) <= 0.5*
getFrameTime() &&
348 f1.GetParameter(0) >= f1.GetParameter(3));
351 if (h1->GetEntries() != 0) {
357 for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
359 const Double_t x = h1->GetBinCenter (i);
360 const Double_t y = h1->GetBinContent(i);
362 while (x > (ns + 1) *
getFrameTime() - TMax_ns) { ++ns; }
372 const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
373 const Double_t P = TMath::PoissonI(i->second.getTotal(), noise);
375 if (P < Pmin && i->second.getTotal() > noise)
381 if (ps.size() != 1) {
382 ERROR(
"Module " << setw(8) << module->getID() <<
" number of peaks " << ps.size() <<
" != 1" << endl);
388 const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
390 if (ps.size() != 1) {
391 WARNING(
"Peak at " << setw(4) << i->first <<
" [frame time] S/N = "
392 <<
FIXED(7,1) << i->second.getTotal() <<
" / "
393 <<
FIXED(7,1) << noise << endl);
396 subPeaksPerDoms.push_back(module->getID());
397 subPeaksPerDoms.push_back(ps.size());
398 subPeaksPerDoms.push_back(i->first);
399 subPeaksPerDoms.push_back(i->second.getTotal());
400 subPeaksPerDoms.push_back(noise);
401 peaksPerDoms.push_back(subPeaksPerDoms);
408 if ((f1.GetParameter(0) > most_OOS_peak) && std::find(oosDomsId.begin(), oosDomsId.end(), module->getID()) == oosDomsId.end()){
410 most_OOS_peak = f1.GetParameter(0);
411 most_OOS_ID = module->getID();
417 if (most_OOS_ID != 0) {
418 oosDomsId.push_back(most_OOS_ID);
420 NOTICE(
"Module " << most_OOS_ID <<
" set QEs of all PMTs to zero." << endl);
423 parameters[JPMTIdentifier(most_OOS_ID, pmt)].QE = 0.0;
426 peaksPerDoms.clear();
438 if (multiPeaks == 1) {
440 NOTICE(
"Multi-peaks analysis activated" << endl);
445 for(
int i = 0; i != int(peaksPerDoms.size()); i++) {
446 if((std::find(workingDOMs.begin(), workingDOMs.end(), peaksPerDoms[i][0]) == workingDOMs.end()) &&
447 (std::find(oosDomsId.begin(), oosDomsId.end(), peaksPerDoms[i][0]) == oosDomsId.end())) {
448 workingDOMs.push_back(peaksPerDoms[i][0]);
452 for(
int i = 0; i != int(workingDOMs.size()); i++) {
455 double sumAllPeaks = 0.;
457 for(
int j = 0; j != int(peaksPerDoms.size()); j++) {
458 if (peaksPerDoms[j][0] == workingDOMs[i]){
459 if (
int(peaksPerDoms[j][2]) != 0) {
460 ratio += peaksPerDoms[j][3] - peaksPerDoms[j][4];
462 sumAllPeaks += peaksPerDoms[j][3] - peaksPerDoms[j][4];
466 ratio /= sumAllPeaks;
469 NOTICE(
"Module " << workingDOMs[i] <<
" set QEs of all PMTs to zero." << endl);
470 NOTICE(
"----> Ratio = " << 100*ratio <<
"%" << endl);
472 parameters[JPMTIdentifier(workingDOMs[i], pmt)].QE = 0.0;
480 NOTICE(
"Managing non-working DOMs." << endl);
481 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
482 JManager_t::iterator p = manager.find(module->getID());
483 if (p == manager.end()){
484 NOTICE(
"Module " << module->getID() <<
" set QEs of all PMTs to zero." << endl);
486 parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
491 ofstream out(pmtFile.c_str());
493 parameters.comment.add(
JMeta(argc, argv));
495 out << parameters << endl;
505 ofstream stream(peaks);
507 stream <<
"#DOM NB_PEAKS TIMESLICE_SHIFT SIGNAL NOISE\n";
509 for(
int i = 0; i != int(peaksPerDoms.size()); i++) {
510 for(
int j = 0; j != int(peaksPerDoms[i].size()); j++){
511 if (j == 0 || j == 1 || j == 2 || j == 3) stream << int(peaksPerDoms[i][j]) <<
" " ;
512 else stream << peaksPerDoms[i][j] <<
" " ;
Utility class to parse command line options.
void set(int &mask) const
Set bit in given bit mask.
Auxiliary class to manage set of histograms.
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
Empty structure for specification of parser element that is initialised (i.e.
Auxiliary class for a type holder.
Auxiliary data structure for floating point format specification.
double getTime(const Hit &hit)
Get true time of hit.
int getFrameIndex() const
Get frame index.
JLimit JLimit_t
Type definition of limit.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getFrameTime()
Get frame time duration.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
double getMaximalTime(const JDetector &detector)
Get maximal time between modules in detector following causality.
Auxiliary data structure for single bit.
const JLimit & getLimit() const
Get limit.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
JTriggerCounter_t next()
Increment trigger counter.
#define DEBUG(A)
Message macros.