71 int main(
int argc,
char **argv)
75 using namespace KM3NETDAQ;
81 bool overwriteDetector;
83 int numberOfTimeslices;
92 JParser<> zap(
"Example program to search for correlations between triggered events and timeslice data.");
94 zap[
'f'] =
make_field(inputFile,
"input file.");
96 zap[
'a'] =
make_field(detectorFile,
"detector file.");
97 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with module (PMT) status.");
98 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
100 zap[
'N'] =
make_field(numberOfTimeslices,
"time slice difference between triggered event and L1.") = 100;
101 zap[
'W'] =
make_field(binWidth_ns,
"bin width for output histograms." ) = 10e3;
102 zap[
'p'] =
make_field(Pmin,
"minimal probability for background to be signal.") = 1.0e-7;
103 zap[
'C'] =
make_field(selector,
"data type selection (default is all).") =
"", getROOTClassSelection<JDAQTimesliceTypes_t>();
109 catch(
const exception& error) {
110 FATAL(error.what() << endl);
126 FATAL(
"Empty detector " << detectorFile << endl);
141 const Double_t xmin = -(numberOfTimeslices + 1) *
getFrameTime() - 0.5*binWidth_ns;
142 const Double_t xmax = +(numberOfTimeslices + 1) *
getFrameTime() + 0.5*binWidth_ns;
143 const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
145 JManager_t manager(
new TH1D(
"M_%", NULL, nx, xmin, xmax));
153 if (selector ==
"") {
161 FATAL(
"No timeslice data." << endl);
164 NOTICE(
"Selected class " << ps->getClassName() << endl);
170 ps->configure(inputFile);
173 ps->setLimit(inputFile.getLimit());
181 STATUS(
"event: " << setw(10) <<
in.getCounter() <<
'\r');
DEBUG(endl);
196 buffer[
event->getFrameIndex()].push_back(t0);
201 while (ps->hasNext()) {
203 STATUS(
"event: " << setw(10) << ps->getCounter() <<
'\r');
DEBUG(endl);
207 map_type::const_iterator p = buffer.lower_bound(timeslice->
getFrameIndex() - numberOfTimeslices);
208 map_type::const_iterator q = buffer.upper_bound(timeslice->
getFrameIndex() + numberOfTimeslices);
210 int number_of_events = 0;
212 for (map_type::const_iterator i = p; i != q; ++i) {
213 number_of_events += i->second.size();
216 if (number_of_events == 0) {
220 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
224 buildL1(*frame, router.
getModule(frame->getModuleID()), back_inserter(data));
226 TH1D*
h1 = manager[frame->getModuleID()];
230 const double t1 = *hit + frame->getFrameIndex() *
getFrameTime();
232 for (map_type::const_iterator i = p; i != q; ++i) {
233 for (map_type::mapped_type::const_iterator
j = i->second.begin();
j != i->second.end(); ++
j) {
235 const double t0 = *
j;
248 TF1 f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
259 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
261 if (module->getFloor() != 0) {
263 status[module->getID()] = DEFAULT;
265 JManager_t::iterator p = manager.find(module->getID());
267 if (p == manager.end() || p->second->GetEntries() == 0) {
269 status[module->getID()] = NODATA;
274 TH1D*
h1 = p->second;
280 Double_t sigma = 250.0;
283 for (
int i = 1; i != h1->GetNbinsX(); ++i) {
285 const Double_t
x = h1->GetBinCenter (i);
286 const Double_t y = h1->GetBinContent(i);
296 ymin /= h1->GetNbinsX();
298 f1.SetParameter(0, ymax);
299 f1.SetParameter(1, x0);
300 if (binWidth_ns < sigma)
301 f1.SetParameter(2, sigma);
303 f1.FixParameter(2, binWidth_ns/sqrt(12.0));
304 f1.SetParameter(3, ymin);
308 h1->Fit(&f1, option.c_str(),
"same", x0 - 5 * sigma, x0 + 5 * sigma);
312 f1.GetParameter(0) >= f1.GetParameter(3)) {
314 status[module->getID()] = IN_SYNC;
323 for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
325 const Double_t
x = h1->GetBinCenter (i);
326 const Double_t y = h1->GetBinContent(i);
328 while (x > (ns + 1) *
getFrameTime() - TMax_ns) { ++ns; }
332 else if (fabs(x - ns *
getFrameTime()) < 10.0 * TMax_ns)
336 DEBUG(
"Module " << setw(8) << module->getID() << endl);
340 const Double_t y = bg[i->first].getTotal() * i->second.getCount() / bg[i->first].getCount();
341 const Double_t
P = TMath::PoissonI(i->second.getTotal(), y);
343 DEBUG(
"Peak at " << setw(4) << i->first <<
" [frame time] "
345 <<
FIXED(7,1) << i->second.getTotal() <<
" / "
346 <<
FIXED(7,1) << y <<
' '
349 if (i->second.getTotal() > y && P < Pmin)
355 if (!(sn.size() == 1 &&
356 sn.begin()->first == 0)) {
358 ERROR(
"Module " << setw(8) << module->getID() <<
" number of peaks " << sn.size() << endl);
361 ERROR(
"Peak at " << setw(4) << i->first <<
" [frame time] "
363 <<
FIXED(7,1) << i->second.getTotal() <<
" / "
364 <<
FIXED(7,1) << bg[i->first].getTotal() * i->second.getCount() / bg[i->first].getCount() << endl);
367 status[module->getID()] = (sn.size() == 1 ? OUT_SYNC :
ERROR);
372 if (overwriteDetector) {
378 if (i->second != IN_SYNC) {
380 NOTICE(
"Module " << setw(8) << i->first <<
" set out-of-sync." << endl);
386 for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
403 if (i->second == IN_SYNC) {
406 if (i->second != IN_SYNC &&
407 i->second != NODATA) {
412 NOTICE(
"Number of modules out/in sync " << nout <<
'/' << nin << endl);
414 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 class to select ROOT class based on class name.
then for HISTOGRAM in h0 h1
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...
Auxiliary class for defining the range of iterations of objects.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
I/O formatting auxiliaries.
void set(const int bit)
Set PMT status.
#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.
const JStatus & getStatus() const
Get status.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
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.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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.
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
const JPMT & getPMT(const JDAQKeyHit &hit) const
Get PMT parameters.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Auxiliary data structure for floating point format specification.
int qaqc
QA/QC file descriptor.