67 return std::max(bg.getTotal(), 1.0) * (double) sn.getCount() / (double) bg.getCount();
88 int main(
int argc,
char **argv)
92 using namespace KM3NETDAQ;
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.");
116 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
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;
122 zap[
'C'] =
make_field(selector,
"data type selection (default is all).") =
"", getROOTClassSelection<JDAQTimesliceTypes_t>();
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);
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 STATUS(
"event: " << setw(10) <<
in.getCounter() <<
'\r');
DEBUG(endl);
217 buffer[
event->getFrameIndex()].push_back(t0);
222 while (ps->hasNext()) {
224 STATUS(
"event: " << setw(10) << ps->getCounter() <<
'\r');
DEBUG(endl);
228 map_type::const_iterator p = buffer.lower_bound(timeslice->
getFrameIndex() - numberOfTimeslices);
229 map_type::const_iterator q = buffer.upper_bound(timeslice->
getFrameIndex() + numberOfTimeslices);
231 int number_of_events = 0;
233 for (map_type::const_iterator
i = p;
i != q; ++
i) {
234 number_of_events +=
i->second.size();
237 if (number_of_events == 0) {
241 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
245 buildL1(*frame, router.
getModule(frame->getModuleID()), back_inserter(
data));
247 TH1D* h1 = manager[frame->getModuleID()];
249 double t1 = numeric_limits<double>::lowest();
253 const double t2 = *hit + frame->getFrameIndex() *
getFrameTime();
255 if (t2 > t1 + deadTime_ns) {
257 for (map_type::const_iterator
i = p;
i != q; ++
i) {
258 for (map_type::mapped_type::const_iterator t0 =
i->second.begin(); t0 !=
i->second.end(); ++t0) {
273 TF1
f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
284 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
286 if (module->getFloor() != 0) {
288 status[module->getID()] = DEFAULT;
290 JManager_t::iterator p = manager.find(module->getID());
292 if (p == manager.end() || p->second->GetEntries() == 0) {
294 status[module->getID()] = NODATA;
299 TH1D* h1 = p->second;
305 Double_t
sigma = 250.0;
308 for (
int i = 1;
i != h1->GetNbinsX(); ++
i) {
310 const Double_t
x = h1->GetBinCenter (
i);
311 const Double_t
y = h1->GetBinContent(
i);
321 ymin /= h1->GetNbinsX();
323 f1.SetParameter(0, ymax);
324 f1.SetParameter(1, x0);
325 if (binWidth_ns < sigma)
326 f1.SetParameter(2, sigma);
328 f1.FixParameter(2, binWidth_ns/sqrt(12.0));
329 f1.SetParameter(3, ymin);
331 for (Int_t
i = 0;
i != f1.GetNpar(); ++
i) {
332 f1.SetParError(
i, 0.0);
337 h1->Fit(&f1, option.c_str(),
"same", x0 - 5 *
sigma, x0 + 5 *
sigma);
341 f1.GetParameter(0) >= f1.GetParameter(3)) {
343 status[module->getID()] = IN_SYNC;
352 for (Int_t
i = 1,
ns = -(numberOfTimeslices + 1);
i <= h1->GetNbinsX(); ++
i) {
354 const Double_t
x = h1->GetBinCenter (
i);
355 const Double_t
y = h1->GetBinContent(
i);
369 const Double_t
y = getBackground(
i->second, bg[
i->first]);
370 const Double_t
P = TMath::PoissonI(
i->second.getTotal(),
y);
372 STATUS(
"Module " << setw(10) << module->getID()
374 << setw(4) <<
i->first <<
" [frame time] "
376 <<
FIXED(7,1) <<
i->second.getTotal() <<
" / "
377 <<
FIXED(7,1) << y <<
' '
380 if (
i->second.getTotal() > y && P < Pmin)
386 if (!(sn.size() == 1 &&
387 sn.begin()->first == 0)) {
389 ERROR(
"Module " << setw(8) << module->getID() <<
" number of peaks " << sn.size() << endl);
393 ERROR(
"Peak at " << setw(4) <<
i->first <<
" [frame time] "
395 <<
FIXED(7,1) <<
i->second.getTotal() <<
" / "
396 <<
FIXED(7,1) << getBackground(
i->second, bg[
i->first]) << endl);
399 status[module->getID()] = (sn.size() == 1 ? OUT_SYNC :
ERROR);
404 if (overwriteDetector) {
410 if (
i->second != IN_SYNC &&
411 i->second != NODATA) {
413 NOTICE(
"Module " << setw(8) <<
i->first <<
" set out-of-sync." << endl);
419 for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
436 if (
i->second == IN_SYNC) {
439 if (
i->second != IN_SYNC &&
440 i->second != NODATA) {
445 NOTICE(
"Number of modules out/in sync " << nout <<
'/' << nin << endl);
447 QAQC(nin <<
' ' << nout << endl);
Utility class to parse command line options.
Direct access to PMT data in detector data structure for DAQ hits.
int main(int argc, char *argv[])
ROOT TTree parameter settings of various packages.
Data structure for a composite optical module.
Auxiliary data structure to unify weights of acoustics data according to the number of pings per emit...
Auxiliary class to select ROOT class based on class name.
Auxialiary class to assert type conversion.
const JModule & getModule(const JDAQKeyHit &hit) const
Get module parameters.
Dynamic ROOT object management.
Auxiliary class for a type holder.
Auxiliary data structure for floating point format specification.
double getTime(const Hit &hit)
Get true time of hit.
Template definition for direct access of elements in ROOT TChain.
Data structure for detector geometry and calibration.
Simple wrapper around JModuleRouter class for direct addressing of PMT data in detector data structur...
int getFrameIndex() const
Get frame index.
The template JSinglePointer class can be used to hold a pointer to an object.
Scanning of objects from a single file according a format that follows from the extension of each fil...
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Auxiliary class for defining the range of iterations of objects.
int getStatus() const
Get status.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getFrameTime()
Get frame time duration.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
double getMaximalTime(const double R_Hz)
Get maximal time for given rate.
General purpose messaging.
Auxiliary class to select JTreeScanner based on ROOT class name.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Utility class to parse command line options.
#define QAQC(A)
QA/QC output macro.
static const int OUT_OF_SYNC
Enable (disable) synchronous signal from this PMT if this status bit is 0 (1);.
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
static const int MODULE_OUT_OF_SYNC
Enable (disable) synchronous signal from this module if this status bit is 0 (1);.
Auxiliary data structure for L1 build parameters.
const JLimit & getLimit() const
Get limit.
std::map< int, range_type > map_type
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
const JPMT & getPMT(const JDAQKeyHit &hit) const
Get PMT parameters.
Auxiliary data structure for floating point format specification.
int qaqc
QA/QC file descriptor.
#define DEBUG(A)
Message macros.