Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JTurbotIter.cc
Go to the documentation of this file.
1
2
3#include <string>
4#include <iostream>
5#include <iomanip>
6#include <limits>
7#include <map>
8#include <vector>
9
10#include "TROOT.h"
11#include "TFile.h"
12#include "TH1D.h"
13#include "TF1.h"
14#include "TMath.h"
15
16#include "JROOT/JMinimizer.hh"
17
19#include "JDAQ/JDAQEventIO.hh"
21#include "JDAQ/JDAQEvaluator.hh"
26#include "JTrigger/JBuildL0.hh"
27#include "JTrigger/JBuildL1.hh"
29#include "JROOT/JManager.hh"
30#include "JTools/JWeight.hh"
36#include "JSupport/JSupport.hh"
37#include "JSupport/JMeta.hh"
38
39#include "Jeep/JPrint.hh"
40#include "Jeep/JParser.hh"
41#include "Jeep/JMessage.hh"
42
43namespace {
44
45 // status of module
46 enum JStatus_t {
47 DEFAULT = 0, //!< module exists in detector
48 NODATA = -2, //!< module not present in data or no entries
49 OUT_SYNC = -3, //!< module is out-of-sync
50 ERROR = -4, //!< module is out-of-sync, possibly multiple time offsets
51 IN_SYNC = +1 //!< module is in-sync
52 };
53}
54
55
56/**
57 * \file
58 *
59 * Example program to search for correlations between triggered events and timeslice data.
60 * \author mdejong
61 */
62int main(int argc, char **argv)
63{
64 using namespace std;
65 using namespace JPP;
66 using namespace KM3NETDAQ;
67
68 JSingleFileScanner<> inputFile;
69 JLimit_t& numberOfEvents = inputFile.getLimit();
70 string outputFile;
71 string detectorFile;
72 double TMaxLocal_ns;
73 string pmtFile;
74 int numberOfTimeslices;
75 double binWidth_ns;
76 double Pmin;
77 JROOTClassSelector selector;
78 int qaqc;
79 int debug;
80
81 string peaks;
82
83 try {
84
85 JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
86
87 zap['f'] = make_field(inputFile);
88 zap['o'] = make_field(outputFile) = "";
89 zap['a'] = make_field(detectorFile);
90 zap['n'] = make_field(numberOfEvents) = JLimit::max();
91 zap['T'] = make_field(TMaxLocal_ns) = 20.0;
92 zap['P'] = make_field(pmtFile) = "";
93 zap['N'] = make_field(numberOfTimeslices) = 100;
94 zap['W'] = make_field(binWidth_ns) = 10e3;
95 zap['p'] = make_field(Pmin) = 1.0e-6;
97 zap['Q'] = make_field(qaqc) = 0;
98 zap['d'] = make_field(debug) = 2;
99
100 zap['O'] = make_field(peaks) = ""; //Make a small summary file, with the list of the peaks
101
102 zap.read(argc, argv);
103 }
104 catch(const exception& error) {
105 FATAL(error.what() << endl);
106 }
107
108
110
111 try {
112 load(detectorFile, detector);
113 }
114 catch(const JException& error) {
115 FATAL(error);
116 }
117
118 const JDAQHitRouter router(detector);
119
120 const double TMax_ns = max(binWidth_ns, getMaximalTime(detector));
121
122 typedef double hit_type;
123
124 JBuildL0<hit_type> buildL0;
125 JBuildL1<hit_type> buildL1(JBuildL1Parameters(TMaxLocal_ns, true));
126 vector <hit_type> data;
127
128 typedef JManager<int, TH1D> JManager_t;
129
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);
133
134 JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));
135
137
139
141
142 if (selector == "") {
143
144 if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
145 (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
146 (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
147 (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
148 (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
149 } else {
150 FATAL("No timeslice data." << endl);
151 }
152
153 NOTICE("Selected class " << ps->getClassName() << endl);
154
155 } else {
156
157 ps = zmap[selector];
158
159 ps->configure(inputFile);
160 }
161
162 /* Updated method :
163 * OOS DOMs can have an impact on the histogram of other DOMs
164 * especialy OOS DOMs on adjacent In Sync or OOS DOMs (presence of ghost peaks)
165 * To remove these ghost peaks, the idea is to iterate the method, removing
166 * at each loop the most OOS DOM, and recomputing the histograms without
167 * them (removing events where they are involved)
168 * Once there is no more OOS DOMs in the remaining DOMs, end of precedure.
169 * Steps of procedure (from J. Hofestaedt, code updated by R. Le Breton)
170 * 1) Check if any DOM shows a peak away from 0 (OOS)
171 * If no OOS -> finish procedure
172 * 2) If at least one OOS, then identify the DOM that is most OOS,
173 * ie. highest peak (other than delta time = 0)
174 * 3) Remove the most OOS DOM and calculate the histograms again but do not
175 * consider the OOS DOM (in later iterations more than one DOM is OOS and not
176 * considered for histogram filling).
177 * 4) start with 1) again
178 */
179
180 ps->setLimit(inputFile.getLimit());
181
182 typedef map<int, vector<double> > map_type; // frame index -> event times
183
184 bool OOS = true; //Boolean: to continue or not the procedure
185 map_type buffer;
186
187 map_type peaksPerDoms;
188
189 map<int, JStatus_t> status2;
190
191 //OOS DOMs IDs : if at least one OOS DOM, one ID is added into oosDomsId at each iteration
192 map<int, bool> oosDomsId;
193
194 //Default value for all modules : false
195 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
196 oosDomsId[module->getID()] = false;
197 }
198
199 do{
200 for (JTreeScanner<JDAQEvent> in(inputFile.getFilename()); in.hasNext(); ) {
201
202 STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
203
204 JDAQEvent* event = in.next();
205
206 // Average time of triggered hits
207
208 double t0 = 0.0;
209 bool test_OOS = false; //If an OOS DOM is involved in an event, remove this event
210
211 for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
212 t0 += getTime(*hit, router.getPMT(*hit));
213
214 if(oosDomsId[hit->getModuleID()]) {
215 test_OOS = true;
216 break;
217 }
218 }
219
220 if(test_OOS) continue; //If true, event not taken into account
221
222 t0 /= event->size<JDAQTriggeredHit>();
223 t0 += event->getFrameIndex() * getFrameTime();
224
225 buffer[event->getFrameIndex()].push_back(t0);
226 }
227 STATUS(endl);
228
229 //Events : start from the beginning
230 ps->rewind();
231
232 while (ps->hasNext()) {
233
234 STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
235
236 JDAQTimeslice* timeslice = ps->next();
237
238 map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
239 map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
240
241 int number_of_events = 0;
242
243 for (map_type::const_iterator i = p; i != q; ++i) {
244 number_of_events += i->second.size();
245 }
246
247 if (number_of_events == 0) {
248 continue;
249 }
250
251 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
252
253 data.clear();
254
255 buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
256
257 TH1D* h1 = manager[frame->getModuleID()];
258
259 for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
260
261 const double t1 = *hit + frame->getFrameIndex() * getFrameTime();
262
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) {
265
266 const double t0 = *j;
267
268 h1->Fill(t1 - t0);
269 }
270 }
271 }
272 }
273 }
274
275 buffer.clear(); //Clear the buffer for the next loop, if any
276
277 STATUS(endl);
278
279 Double_t most_OOS_peak = 0;
280 Int_t most_OOS_ID = 0;
281
282 TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
283
284 string option = "L";
285
286 if (debug < debug_t) {
287 option += "Q";
288 }
289
290 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
291
292 bool status = false;
293
294 status2[module->getID()] = DEFAULT;
295
296 JManager_t::iterator p = manager.find(module->getID());
297
298 if (p == manager.end() || p->second->GetEntries() == 0) {
299
300 status2[module->getID()] = NODATA;
301
302 continue;
303 }
304
305 TH1D* h1 = p->second;
306
307 // start values
308
309 Double_t ymax = 0.0;
310 Double_t x0 = 0.0; // [ns]
311 Double_t sigma = 250.0; // [ns]
312 Double_t ymin = 0.0;
313
314 for (int i = 1; i != h1->GetNbinsX(); ++i) {
315
316 const Double_t x = h1->GetBinCenter (i);
317 const Double_t y = h1->GetBinContent(i);
318
319 if (y > ymax) {
320 ymax = y;
321 x0 = x;
322 }
323
324 ymin += y;
325 }
326
327 ymin /= h1->GetNbinsX();
328
329 f1.SetParameter(0, ymax);
330 f1.SetParameter(1, x0);
331 if (binWidth_ns < sigma)
332 f1.SetParameter(2, sigma);
333 else
334 f1.FixParameter(2, binWidth_ns/sqrt(12.0));
335 f1.SetParameter(3, ymin);
336
337 for (Int_t i = 0; i != f1.GetNpar(); ++i) {
338 f1.SetParError(i, 0.0);
339 }
340
341 // fit
342
343 h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
344
345 status = (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() && // position
346 f1.GetParameter(0) >= f1.GetParameter(3)); // S/N
347
348 if (status) status2[module->getID()] = IN_SYNC;
349
350
351 // look for peaks at offsets equal to a multiple of frame time
352 JWeight bg; // background
353 map<int, JWeight> sn; // signal(s)
354
355 for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
356
357 const Double_t x = h1->GetBinCenter (i);
358 const Double_t y = h1->GetBinContent(i);
359
360 while (x > (ns + 1) * getFrameTime() - TMax_ns) { ++ns; }
361
362 if (fabs(x - ns * getFrameTime()) < TMax_ns)
363 sn[ns].put(y);
364 else
365 bg .put(y);
366 }
367
368 DEBUG("Module " << setw(8) << module->getID() << endl);
369
370 for (map<int, JWeight>::iterator i = sn.begin(); i != sn.end(); ) {
371
372 const Double_t y = bg.getTotal() * i->second.getCount() / bg.getCount();
373 const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
374
375 DEBUG("Peak at " << setw(4) << i->first << " [frame time] "
376 << "S/N = "
377 << FIXED(7,1) << i->second.getTotal() << " / "
378 << FIXED(7,1) << y << ' '
379 << SCIENTIFIC(12,3) << P << endl);
380
381 if (P < Pmin && i->second.getTotal() > y) //Avoid edge effect (second bolean)
382 ++i; // keep
383 else
384 sn.erase(i++); // remove
385 }
386
387 if (!(sn.size() == 1 &&
388 sn.begin()->first == 0)) {
389
390 WARNING("Module " << setw(8) << module->getID() << " number of peaks " << sn.size() << endl);
391
392 for (map<int, JWeight>::const_iterator i = sn.begin(); i != sn.end(); ++i) {
393
394 const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
395
396 WARNING("Peak at " << setw(4) << i->first << " [frame time] "
397 << "S/N = "
398 << FIXED(7,1) << i->second.getTotal() << " / "
399 << FIXED(7,1) << noise << endl);
400
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);
405 }
406
407 status2[module->getID()] = (sn.size() == 1 ? OUT_SYNC : ERROR);
408 }
409
410
411 //if DOM is OOS
412 if (!status) {
413 //if OOS peak greater than the last and not in OOS DOMs vector
414 if ((f1.GetParameter(0) > most_OOS_peak) && !(oosDomsId[module->getID()])){
415
416 most_OOS_peak = f1.GetParameter(0);
417 most_OOS_ID = module->getID();
418 }
419 }
420 }
421
422 if (most_OOS_ID != 0) {
423 oosDomsId[most_OOS_ID] = true;//exclude the DOM to compute the histograms
424
425 //Reseting all histograms
426 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
427 manager[module->getID()]->Reset("ICESM");
428 }
429 peaksPerDoms.clear();
430 }
431 else{
432 OOS = false;
433 }
434 }while(OOS);
435
436 //To manage non working and OOS DOMs
437 NOTICE("Managing non-working and OOS DOMs." << endl);
438 if (pmtFile != "") {
439
440 JPMTParametersMap parameters;
441
442 try {
443 parameters.load(pmtFile.c_str());
444 }
445 catch(const JException& error) {}
446
447 for (map<int, JStatus_t>::const_iterator i = status2.begin(); i != status2.end(); ++i) {
448
449 if (i->second != IN_SYNC) {
450
451 NOTICE("Module " << setw(8) << i->first << " set QEs of all PMTs to zero." << endl);
452
453 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
454 parameters[JPMTIdentifier(i->first, pmt)].QE = 0.0;
455 }
456 }
457 }
458
459 ofstream out(pmtFile.c_str());
460
461 parameters.comment.add(JMeta(argc, argv));
462
463 out << parameters << endl;
464
465 out.close();
466 }
467
468 if (outputFile != "") {
469 manager.Write(outputFile.c_str());
470 }
471
472 //print small summary file of the peaks
473 if (peaks != ""){
474 ofstream stream(peaks);
475
476 stream << "#DOM NB_PEAKS TIMESLICE_SHIFT SIGNAL NOISE\n";
477
478 for (map_type::const_iterator dom = peaksPerDoms.begin(); dom != peaksPerDoms.end(); ++dom) {
479 stream << dom->first << " ";
480 int i = 1;
481 for (map_type::mapped_type::const_iterator data = dom->second.begin(); data != dom->second.end(); ++data) {
482 stream << *data;
483 ((i%4 == 0) ? stream << "\n" : stream << " ");
484 i++;
485 }
486 }
487 stream.close();
488 }
489
490 //QAQC
491 int nin = 0;
492 int nout = 0;
493
494 for (map<int, JStatus_t>::const_iterator i = status2.begin(); i != status2.end(); ++i) {
495 if (i->second == IN_SYNC) {
496 ++nin;
497 }
498 if (i->second != IN_SYNC &&
499 i->second != NODATA) {
500 ++nout;
501 }
502 }
503
504 NOTICE("Number of modules out/in sync " << nout << '/' << nin << endl);
505
506 QAQC(nin << ' ' << nout << endl);
507
508 return 0;
509
510}
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
#define WARNING(A)
Definition JMessage.hh:65
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.
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)
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
Auxiliary class for map of PMT parameters.
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
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition JParser.hh:1992
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.
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.
@ 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.
int j
Definition JPolint.hh:801
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.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition JDAQ.hh:26
std::map< int, range_type > map_type
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
JComment & add(const std::string &comment)
Add comment.
Definition JComment.hh:100
void load(const char *file_name)
Load from input file.
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