Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JTurbot.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <limits>
5#include <map>
6#include <vector>
7
8#include "TROOT.h"
9#include "TFile.h"
10#include "TH1D.h"
11#include "TF1.h"
12#include "TMath.h"
13
14#include "JROOT/JMinimizer.hh"
15
19
20#include "JDAQ/JDAQEventIO.hh"
22#include "JDAQ/JDAQEvaluator.hh"
23
27#include "JTrigger/JBuildL0.hh"
28#include "JTrigger/JBuildL1.hh"
30#include "JROOT/JManager.hh"
31#include "JTools/JWeight.hh"
37#include "JSupport/JSupport.hh"
38#include "JSupport/JMeta.hh"
39
40#include "Jeep/JPrint.hh"
41#include "Jeep/JParser.hh"
42#include "Jeep/JMessage.hh"
43
44namespace {
45
46 using JTOOLS::JWeight;
47
48 // status of module
49
50 enum JStatus_t {
51 DEFAULT = 0, //!< module exists in detector
52 NODATA = -2, //!< module not present in data or no entries
53 OUT_SYNC = -3, //!< module is out-of-sync
54 ERROR = -4, //!< module is out-of-sync, possibly multiple time offsets
55 IN_SYNC = +1 //!< module is in-sync
56 };
57
58 /**
59 * Get estimated number of background events.
60 *
61 * \param sn weight of signal
62 * \param bg weight of background
63 * \return background
64 */
65 double getBackground(const JWeight& sn, const JWeight& bg)
66 {
67 return std::max(bg.getTotal(), 1.0) * (double) sn.getCount() / (double) bg.getCount();
68 }
69}
70
71/**
72 * \file
73 *
74 * Example program to search for correlations between triggered events and time slice data.\n
75 *
76 * This application can be used to identify optical modules that are out-of-sync with respect to the central White Rabbit clock.\n
77 * To this end, time correlations are searched for between triggered events and local coincidences, referred to as L1 hits.
78 * The local coincidences are independently obtained from time slice data and can thereby cross the time slice border of a triggered event.\n
79 * The MODULE_OUT_OF_SYNC of status the modules that are found to be out-of-sync are set.
80 *
81 * The number of time slices to look for L1 hits away from the triggered event can be set via command line option <tt>-N <number of time slices></tt>.\n
82 * The larger this value, the bigger time offsets can be covered.
83 * To identify the modules that are out-of-sync as fast as possible (without determining the time offset),
84 * this value can be set to 0 with the caveat that a out-of-sync occurring during the data taking run may then be missed.
85 *
86 * \author mdejong
87 */
88int main(int argc, char **argv)
89{
90 using namespace std;
91 using namespace JPP;
92 using namespace KM3NETDAQ;
93
94 JSingleFileScanner<> inputFile;
95 JLimit_t& numberOfEvents = inputFile.getLimit();
96 string outputFile;
97 string detectorFile;
98 bool overwriteDetector;
99 double TMaxLocal_ns;
100 int numberOfTimeslices;
101 double binWidth_ns;
102 double deadTime_us;
103 double Pmin;
104 JROOTClassSelector selector;
105 int qaqc;
106 int debug;
107
108 try {
109
110 JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
111
112 zap['f'] = make_field(inputFile, "input file.");
113 zap['o'] = make_field(outputFile, "optional output file containing ROOT histograms.") = "";
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();
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;
122 zap['C'] = make_field(selector, "data type selection (default is all).") = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
123 zap['Q'] = make_field(qaqc) = 0;
124 zap['d'] = make_field(debug) = 2;
125
126 zap(argc, argv);
127 }
128 catch(const exception& error) {
129 FATAL(error.what() << endl);
130 }
131
132
134
135 try {
136 load(detectorFile, detector);
137 }
138 catch(const JException& error) {
139 FATAL(error);
140 }
141
142 if (detector.empty()) {
143 FATAL("Empty detector " << detectorFile << endl);
144 }
145
146 const JDAQHitRouter router(detector);
147
148 const double TMax_ns = max(binWidth_ns, getMaximalTime(detector)); // signal time window
149 const double BOOST = 20.0; // scaling of time window for background
150 const double deadTime_ns = deadTime_us * 1e3;
151
152 NOTICE("Time window " << FIXED(7,1) << TMax_ns << " [ns]" << endl);
153
154 typedef double hit_type;
155
156 JBuildL0<hit_type> buildL0;
157 JBuildL1<hit_type> buildL1(JBuildL1Parameters(TMaxLocal_ns, true));
158 vector <hit_type> data;
159
160 typedef JManager<int, TH1D> JManager_t;
161
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);
165
166 JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));
167
169
171
173
174 if (selector == "") {
175
176 if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
177 (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
178 (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
179 (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
180 (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
181 } else {
182 FATAL("No timeslice data." << endl);
183 }
184
185 NOTICE("Selected class " << ps->getClassName() << endl);
186
187 } else {
188
189 ps = zmap[selector];
190
191 ps->configure(inputFile);
192 }
193
194 ps->setLimit(inputFile.getLimit());
195
196 typedef map<int, vector<double> > map_type; // frame index -> event times
197
198 map_type buffer;
199
200 for (JTreeScanner<JDAQEvent> in(inputFile.getFilename()); in.hasNext(); ) {
201
202 if (in.getCounter()%1000 == 0) {
203 STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
204 }
205
206 JDAQEvent* event = in.next();
207
208 // Average time of triggered hits
209
210 double t0 = 0.0;
211 size_t n0 = 0;
212
213 for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
214 if (router.hasModule(hit->getModuleID())) {
215 t0 += getTime(*hit, router.getPMT(*hit));
216 n0 += 1;
217 }
218 }
219
220 t0 /= n0;
221 t0 += event->getFrameIndex() * getFrameTime();
222
223 buffer[event->getFrameIndex()].push_back(t0);
224 }
225 STATUS(endl);
226
227
228 while (ps->hasNext()) {
229
230 if (ps->getCounter()%1000 == 0) {
231 STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
232 }
233
234 JDAQTimeslice* timeslice = ps->next();
235
236 map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
237 map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
238
239 int number_of_events = 0;
240
241 for (map_type::const_iterator i = p; i != q; ++i) {
242 number_of_events += i->second.size();
243 }
244
245 if (number_of_events == 0) {
246 continue;
247 }
248
249 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
250
251 if (router.hasModule(frame->getModuleID())) {
252
253 data.clear();
254
255 buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
256
257 TH1D* h1 = manager[frame->getModuleID()];
258
259 double t1 = numeric_limits<double>::lowest();
260
261 for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
262
263 const double t2 = *hit + frame->getFrameIndex() * getFrameTime();
264
265 if (t2 > t1 + deadTime_ns) {
266
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) {
269 h1->Fill(t2 - *t0);
270 }
271 }
272
273 t1 = t2;
274 }
275 }
276 }
277 }
278 }
279 STATUS(endl);
280
281
282 // fit function
283
284 TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
285
286 string option = "L";
287
288 if (debug < debug_t) {
289 option += "Q";
290 }
291
292
293 map<int, JStatus_t> status;
294
295 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
296
297 if (module->getFloor() != 0) {
298
299 status[module->getID()] = DEFAULT;
300
301 JManager_t::iterator p = manager.find(module->getID());
302
303 if (p == manager.end() || p->second->GetEntries() == 0) {
304
305 status[module->getID()] = NODATA;
306
307 continue;
308 }
309
310 TH1D* h1 = p->second;
311
312 // start values
313
314 Double_t ymax = 0.0;
315 Double_t x0 = 0.0; // [ns]
316 Double_t sigma = 250.0; // [ns]
317 Double_t ymin = 0.0;
318
319 for (int i = 1; i != h1->GetNbinsX(); ++i) {
320
321 const Double_t x = h1->GetBinCenter (i);
322 const Double_t y = h1->GetBinContent(i);
323
324 if (y > ymax) {
325 ymax = y;
326 x0 = x;
327 }
328
329 ymin += y;
330 }
331
332 ymin /= h1->GetNbinsX();
333
334 f1.SetParameter(0, ymax);
335 f1.SetParameter(1, x0);
336 if (binWidth_ns < sigma)
337 f1.SetParameter(2, sigma);
338 else
339 f1.FixParameter(2, binWidth_ns/sqrt(12.0));
340 f1.SetParameter(3, ymin);
341
342 for (Int_t i = 0; i != f1.GetNpar(); ++i) {
343 f1.SetParError(i, 0.0);
344 }
345
346 // fit
347
348 h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
349
350 if (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() && // position
351 f1.GetParameter(0) >= f1.GetParameter(3)) { // S/N
352
353 status[module->getID()] = IN_SYNC;
354 }
355
356 double fm = fmod(f1.GetParameter(1), getFrameTime());
357
358 if (fm < -0.5 * getFrameTime()) { fm += getFrameTime(); }
359 if (fm > +0.5 * getFrameTime()) { fm -= getFrameTime(); }
360
361 STATUS("Module/fit "
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);
368
369 // look for peaks at time offsets equal to a multiple of frame time
370
371 map<int, JWeight> bg; // background
372 map<int, JWeight> sn; // signal(s)
373
374 for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
375
376 const Double_t x = h1->GetBinCenter (i);
377 const Double_t y = h1->GetBinContent(i);
378
379 while (x > (ns + 1) * getFrameTime() - BOOST * TMax_ns) {
380 ++ns;
381 }
382
383 if (fabs(x - ns * getFrameTime()) < TMax_ns)
384 sn[ns].put(y);
385 else if (fabs(x - ns * getFrameTime()) < BOOST * TMax_ns)
386 bg[ns].put(y);
387 }
388
389 for (map<int, JWeight>::iterator i = sn.begin(); i != sn.end(); ) {
390
391 const Double_t y = getBackground(i->second, bg[i->first]);
392 const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
393
394 if (debug >= debug_t || status[module->getID()] != IN_SYNC) {
395 cout << "Module/peak " << setw(10) << module->getID() << ' '
396 << setw(4) << i->first << ' '
397 << "["
398 << FIXED(15,1) << i->first * getFrameTime() - TMax_ns
399 << ","
400 << FIXED(15,1) << i->first * getFrameTime() + TMax_ns
401 << "]"
402 << " S/N = "
403 << FIXED(7,1) << i->second.getTotal() << " / "
404 << FIXED(7,1) << y << ' '
405 << SCIENTIFIC(12,3) << P << ' '
406 << (i->second.getTotal() > y && P < Pmin && i->first != 0 ? "***" : "") << endl;
407 }
408
409 if (i->second.getTotal() > y && P < Pmin)
410 ++i; // keep
411 else
412 sn.erase(i++); // remove
413 }
414
415 if (!(sn.size() == 1 &&
416 sn.begin()->first == 0)) {
417
418 status[module->getID()] = (sn.size() == 1 ? OUT_SYNC : ERROR);
419
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);
425 }
426 }
427 }
428
429 if (overwriteDetector) {
430
431 detector.comment.add(JMeta(argc, argv));
432
433 for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
434
435 if (i->second != IN_SYNC &&
436 i->second != NODATA) {
437
438 NOTICE("Module " << setw(8) << i->first << " set out-of-sync." << endl);
439
440 JModule& module = detector.getModule(router.getAddress(i->first));
441
442 module.getStatus().set(MODULE_OUT_OF_SYNC);
443
444 for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
445 pmt->getStatus().set(OUT_OF_SYNC);
446 }
447 }
448 }
449
450 store(detectorFile.c_str(), detector);
451 }
452
453 if (outputFile != "") {
454 manager.Write(outputFile.c_str());
455 }
456
457 int nin = 0;
458 int nout = 0;
459
460 for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
461 if (i->second == IN_SYNC) {
462 ++nin;
463 }
464 if (i->second != IN_SYNC &&
465 i->second != NODATA) {
466 ++nout;
467 }
468 }
469
470 NOTICE("Number of modules out/in sync " << nout << '/' << nin << endl);
471
472 QAQC(nin << ' ' << nout << endl);
473
474 return 0;
475}
Direct access to PMT data in detector data structure for DAQ hits.
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Dynamic ROOT object management.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define QAQC(A)
QA/QC output macro.
Definition JMessage.hh:100
int qaqc
QA/QC file descriptor.
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
ROOT I/O of application specific meta data.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
Scanning of objects from a single file according a format that follows from the extension of each fil...
ROOT TTree parameter settings of various packages.
int main(int argc, char **argv)
Definition JTurbot.cc:88
Simple wrapper around JModuleRouter class for direct addressing of PMT data in detector data structur...
const JModule & getModule(const JDAQKeyHit &hit) const
Get module parameters.
const JPMT & getPMT(const JDAQKeyHit &hit) const
Get PMT parameters.
Detector data structure.
Definition JDetector.hh:96
bool hasModule(const JObjectID &id) const
Has module.
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
Data structure for a composite optical module.
Definition JModule.hh:75
Auxialiary class to assert type conversion.
General exception.
Definition JException.hh:24
The template JSinglePointer class can be used to hold a pointer to an object.
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
Template definition for direct access of elements in ROOT TChain.
Weight calculator.
Definition JWeight.hh:25
Template L0 hit builder.
Definition JBuildL0.hh:38
Template L1 hit builder.
Definition JBuildL1.hh:90
int getFrameIndex() const
Get frame index.
Template const_iterator.
Definition JDAQEvent.hh:68
JTriggerCounter_t next()
Increment trigger counter.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
@ debug_t
debug
Definition JMessage.hh:29
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::set< JROOTClassSelector > getROOTClassSelection(const bool option=false)
Get ROOT class selection.
const char * getTime()
Get current local time conform ISO-8601 standard.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
double getFrameTime()
Get frame time duration.
Definition JDAQClock.hh:162
double getMaximalTime(const double R_Hz)
Get maximal time for given rate.
std::map< int, range_type > map_type
static const int OUT_OF_SYNC
Enable (disable) synchronous signal from this PMT if this status bit is 0 (1);.
Definition pmt_status.hh:17
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Detector file.
Definition JHead.hh:227
Auxiliary data structure to unify weights of acoustics data according to the number of pings per emit...
Acoustic hit.
Definition JBillabong.cc:70
Auxiliary class for a type holder.
Definition JType.hh:19
Auxiliary class to select ROOT class based on class name.
Auxiliary class to select JTreeScanner based on ROOT class name.
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
Auxiliary data structure for L1 build parameters.
Definition JBuildL1.hh:38
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488