Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JTriggerEfficiency.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <vector>
6 #include <set>
7 
8 #include "TH1D.h"
9 #include "TRandom3.h"
10 
11 #include "evt/Head.hh"
12 #include "evt/Evt.hh"
13 #include "evt/Hit.hh"
14 #include "JDAQ/JDAQTimeslice.hh"
15 #include "JDAQ/JDAQEvent.hh"
16 #include "JDAQ/JDAQSummaryslice.hh"
18 
19 #include "JAAnet/JHead.hh"
20 #include "JAAnet/JHeadToolkit.hh"
21 
22 #include "JDetector/JDetector.hh"
26 #include "JDetector/JPMTRouter.hh"
27 #include "JDetector/JTimeRange.hh"
28 #include "JDetector/JK40Rates.hh"
35 
36 #include "JTrigger/JHit.hh"
37 #include "JTrigger/JHitToolkit.hh"
38 #include "JTrigger/JTimeslice.hh"
41 #include "JTrigger/JHitL0.hh"
42 #include "JTrigger/JHitL1.hh"
43 #include "JTrigger/JBuildL1.hh"
44 #include "JTrigger/JBuildL2.hh"
48 #include "JTrigger/JTriggerBits.hh"
52 #include "JTrigger/JTimesliceL1.hh"
56 #include "JTrigger/JRunByRun.hh"
58 
63 #include "JSupport/JSupport.hh"
66 #include "JSupport/JMeta.hh"
67 
68 #include "Jeep/JPrint.hh"
69 #include "Jeep/JParser.hh"
70 #include "Jeep/JMessage.hh"
71 
72 
73 /**
74  * \file
75  * Auxiliary program to trigger Monte Carlo events.
76  *
77  * Note that option
78  * - <tt>-a</tt> refers to the detector configuration that is used to convert Monte Carlo true information to raw data; and
79  * - <tt>-b</tt> refers to the detector configuration that is used to convert raw data to calibrated data.
80  *
81  * Note that the event counter of the triggered events, as returned by JDAQEvent::getCounter(),
82  * are set to the corresponding index of the Monte Carlo event in the output.
83  * \author mdejong
84  */
85 int main(int argc, char **argv)
86 {
87  using namespace std;
88  using namespace JPP;
89  using namespace KM3NETDAQ;
90 
92 
93  JMultipleFileScanner<> inputFile;
94  JFileRecorder<typelist> outputFile;
95  JLimit_t& numberOfEvents = inputFile.getLimit();
96  string detectorFileA;
97  string detectorFileB;
98  int run;
99  JTriggerParameters parameters;
100  bool triggeredEventsOnly;
101  JPMTParametersMap pmtParameters;
102  JK40Rates rates_Hz;
103  JRunByRun runbyrun;
104  UInt_t seed;
105  int debug;
106 
107  try {
108 
109  JParser<> zap("Auxiliary program to trigger Monte Carlo events.");
110 
111  zap['f'] = make_field(inputFile, "input file (output of detector simulation)");
112  zap['o'] = make_field(outputFile, "output file") = "trigger_efficieny.root";
113  zap['n'] = make_field(numberOfEvents) = JLimit::max();
114  zap['a'] = make_field(detectorFileA, "detector used for converion from Monte Carlo truth to raw data.");
115  zap['b'] = make_field(detectorFileB, "detector used for conversion of raw data to calibrated data.") = "";
116  zap['R'] = make_field(run, "run number") = -1;
117  zap['r'] = make_field(runbyrun, "option for run-by-run mode") = JPARSER::initialised();
118  zap['@'] = make_field(parameters, "Trigger parameters (or corresponding file name)") = JPARSER::initialised();
119  zap['O'] = make_field(triggeredEventsOnly, "optionally write only triggered events.");
120  zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised();
121  zap['B'] = make_field(rates_Hz, "background rates [Hz]") = JPARSER::initialised();
122  zap['S'] = make_field(seed, "seed") = 0;
123  zap['d'] = make_field(debug, "debug") = 0;
124 
125  zap(argc, argv);
126  }
127  catch(const exception &error) {
128  FATAL(error.what() << endl);
129  }
130 
131  gRandom->SetSeed(seed);
132 
133  cout.tie(&cerr);
134 
136 
137  if (detectorFileB == "") {
138  detectorFileB = detectorFileA;
139  }
140 
141 
142  JDetector detectorA;
143  JDetector detectorB;
144 
145  try {
146  load(detectorFileA, detectorA);
147  load(detectorFileB, detectorB);
148  }
149  catch(const JException& error) {
150  FATAL(error);
151  }
152 
153  DEBUG("Detector " << detectorFileA << ' ' << detectorA.size() << endl);
154  DEBUG("Detector " << detectorFileB << ' ' << detectorB.size() << endl);
155 
156  JPMTParametersMap::Throw(true);
157 
158  if (!pmtParameters.is_valid()) {
159  FATAL("Invalid PMT parameters " << pmtParameters << endl);
160  }
161 
162  if (pmtParameters.getQE() != 1.0) {
163 
164  WARNING("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
165 
166  rates_Hz.correct(pmtParameters.getQE());
167  }
168 
169  const JModuleRouter moduleRouter(detectorB);
170  JDetectorSimulator simbad (detectorA);
171  JSummaryRouter summaryRouter;
172 
173  if (runbyrun.is_valid()) {
174 
175  NOTICE("Using run-by-run:" << endl << runbyrun << endl);
176 
177  if (!runbyrun.hasNext()) {
178  FATAL("Run-by-run simulation yields no input." << endl);
179  }
180 
181  if (rates_Hz.getSinglesRate() != 0.0) {
182  WARNING("Run-by-run simulation discards singles rate [Hz] " << rates_Hz.getSinglesRate() << endl);
183  }
184 
185  try {
186  simbad.reset(new JK40RunByRunSimulator(summaryRouter, rates_Hz));
187  simbad.reset(new JPMTRunByRunSimulator(summaryRouter, pmtParameters, detectorA, runbyrun.range_Hz));
188  simbad.reset(new JCLBDefaultSimulator());
189  }
190  catch(const JException& error) {
191  FATAL(error.what() << endl);
192  }
193 
194  try {
195 
196  parameters = getTriggerParameters(runbyrun.get<JMultipleFileScanner<>,true>());
197 
198  NOTICE("Set trigger parameters from run-by-run input." << endl);
199  }
200  catch(const JException& error) {
201  WARNING("No trigger parameters from run-by-run input;\nrun with default/user input." << endl);
202  }
203 
204  } else {
205 
206  NOTICE("Using fixed rates [Hz]: " << rates_Hz << endl);
207 
208  try {
209  simbad.reset(new JK40DefaultSimulator(rates_Hz));
210  simbad.reset(new JPMTDefaultSimulator(pmtParameters, detectorA));
211  simbad.reset(new JCLBDefaultSimulator());
212  }
213  catch(const JException& error) {
214  FATAL(error.what() << endl);
215  }
216  }
217 
218  parameters.set(getMaximalDistance(detectorB));
219 
220  DEBUG("Trigger:" << endl << parameters << endl);
221  DEBUG("PMT parameters:" << endl << pmtParameters << endl);
222 
223  const double Tmax = getMaximalTime(detectorB);
224 
225  const JTimeRange period(-(Tmax + parameters.TMaxLocal_ns),
226  +(Tmax + parameters.TMaxLocal_ns));
227 
228  typedef double hit_type;
229 
230  typedef JSuperFrame1D<hit_type> JSuperFrame1D_t;
231  typedef JSuperFrame2D<hit_type> JSuperFrame2D_t;
232  typedef JTimeslice <hit_type> JTimeslice_t;
233  typedef JBuildL1 <hit_type> JBuildL1_t;
234  typedef JBuildL2 <hit_type> JBuildL2_t;
235 
236  const JBuildL1_t buildL1(parameters);
237  const JBuildL2_t buildL2(parameters.L2);
238  const JBuildL2_t buildSN(parameters.SN);
239 
240  JTimesliceRouter timesliceRouter(parameters.numberOfBins);
241 
242  const JTrigger3DMuon trigger3DMuon (parameters);
243  const JTrigger3DShower trigger3DShower(parameters);
244  const JTriggerMXShower triggerMXShower(parameters, detectorB);
245 
246 
247  Head header;
248 
249  try {
250  header = getHeader(inputFile);
251  }
252  catch(const JException& error) {
253  FATAL(error);
254  }
255 
256  const JPosition3D center = get<JPosition3D>(header);
257 
258  NOTICE("Apply detector offset " << center << endl);
259 
260  detectorA -= center;
261  detectorB -= center;
262 
263  TH1D h1("Trigger bits", NULL, NUMBER_OF_TRIGGER_BITS, -0.5, NUMBER_OF_TRIGGER_BITS - 0.5);
264 
265 
266  outputFile.open();
267 
268  if (!outputFile.is_open()) {
269  FATAL("Error opening file " << outputFile << endl);
270  }
271 
272  outputFile.put(*gRandom);
273  outputFile.put(JMeta(argc, argv));
274  outputFile.put(header);
275  outputFile.put(parameters);
276 
277  JLimit_t limit = inputFile.getLimit();
278  counter_type number_of_events = 0;
279  JTriggerCounter_t trigger_counter = 0;
280 
281  for (JMultipleFileScanner<>::const_iterator file = inputFile.begin(); file != inputFile.end(); ++file) {
282 
283  int mc_run_id = 0;
284 
285  try {
286 
287  const JHead head(getHeader(*file));
288 
289  mc_run_id = head.start_run.run_id;
290  }
291  catch(const JException& error) {
292  FATAL(error);
293  }
294 
295  JMultipleFileScanner<Evt> in(*file);
296 
297  limit.setLowerLimit(limit.getLowerLimit() - in.skip(limit.getLowerLimit()));
298 
299  for ( ; in.hasNext() && number_of_events != limit; ++number_of_events) {
300 
301  STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
302 
303  Evt* event = in.next();
304 
305  event->mc_run_id = mc_run_id;
306 
307  DEBUG(*event << endl);
308 
309  bool trigger = false;
310 
311  if (!event->mc_hits.empty()) {
312 
313  int frame_index = (int) in.getCounter() + 1; // 1 event / slice
314 
315  if (runbyrun.is_valid() && runbyrun.hasNext()) {
316 
317  summaryRouter.update(runbyrun.next());
318 
319  summaryRouter.correct(dynamic_cast<const JPMTDefaultSimulatorInterface&>(simbad.getPMTSimulator()));
320 
321  frame_index = summaryRouter.getFrameIndex();
322  run = summaryRouter.getRunNumber();
323  }
324 
325  // set the event time!
326 
327  JTimeRange timeRange = getTimeRange(*event, period);
328 
329  if (!timeRange.is_valid()) {
330  timeRange.setRange(0.0,0.0);
331  }
332 
333  const double t0 = 0.5 * (timeRange.getLowerLimit() + timeRange.getUpperLimit());
334  const double t1 = getTimeOfFrame(frame_index) + gRandom->Rndm() * getFrameTime();
335 
336  event->mc_t = t1 - t0; // time since start of data taking run
337 
338  timeRange.add(event->mc_t);
339  timeRange.add(period);
340 
341  const JDAQChronometer chronometer(detectorB.getID(), (run != -1 ? run : mc_run_id) , frame_index, JDAQUTCExtended(getTimeOfFrame(frame_index)));
342 
343  const JEventTimeslice timeslice(chronometer, simbad, *event, period);
344 
345  DEBUG(timeslice << endl);
346 
347 
348  timesliceRouter.configure(timeslice);
349 
350  JTimeslice_t timesliceL0(timeslice.getDAQChronometer());
351  JTimeslice_t timesliceL1(timeslice.getDAQChronometer());
352  JTimeslice_t timesliceL2(timeslice.getDAQChronometer());
353  JTimeslice_t timesliceSN(timeslice.getDAQChronometer());
354 
355  for (JDAQTimeslice::const_iterator super_frame = timeslice.begin(); super_frame != timeslice.end(); ++super_frame) {
356 
357  if (moduleRouter.hasModule(super_frame->getModuleID())) {
358 
359  // calibration
360 
361  const JModule& module = moduleRouter.getModule(super_frame->getModuleID());
362  const JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, module);
363 
364  // L0
365 
366  timesliceL0.push_back(JSuperFrame1D_t(buffer));
367 
368  // L1
369 
370  timesliceL1.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
371  super_frame->getModuleIdentifier(),
372  module.getPosition()));
373 
374  buildL1(*timesliceL0.rbegin(), back_inserter(*timesliceL1.rbegin()));
375 
376  // L2
377 
378  timesliceL2.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
379  super_frame->getModuleIdentifier(),
380  module.getPosition()));
381 
382  buildL2(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceL2.rbegin()));
383 
384  // SN
385 
386  timesliceSN.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
387  super_frame->getModuleIdentifier(),
388  module.getPosition()));
389 
390  buildSN(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceSN.rbegin()));
391 
392  DEBUG("L0 " << setw(8) << timesliceL0.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL0.rbegin()->size() << endl);
393  DEBUG("L1 " << setw(8) << timesliceL1.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL1.rbegin()->size() << endl);
394  DEBUG("L2 " << setw(8) << timesliceL2.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL2.rbegin()->size() << endl);
395  DEBUG("SN " << setw(8) << timesliceSN.rbegin()->getModuleID() << ' ' << setw(8) << timesliceSN.rbegin()->size() << endl);
396  }
397  }
398 
399 
400  // Trigger
401 
402  JTriggerInput trigger_input(timesliceL2);
403  JTriggerOutput trigger_output;
404 
405  trigger3DMuon (trigger_input, back_inserter(trigger_output));
406  trigger3DShower(trigger_input, back_inserter(trigger_output));
407  triggerMXShower(trigger_input, timesliceL0, back_inserter(trigger_output));
408 
409  trigger_output.merge(JEventOverlap(parameters.TMaxEvent_ns));
410 
411  for (JTriggerOutput::const_iterator to = trigger_output.begin(); to != trigger_output.end(); ++to) {
412 
413  for (int i = 0; i != h1.GetNbinsX(); ++i) {
414  if (to->hasTriggerBit(i)) {
415  h1.Fill((double) i);
416  }
417  }
418 
419  JTimeRange eventTime = getTimeRange(*to).add(getTimeOfRTS(*to));
420 
421  DEBUG("Event time: "
422  << to->getFrameIndex() << ' '
423  << eventTime << ' '
424  << timeRange << endl);
425 
426  if (timeRange.overlap(eventTime)) {
427 
428  JTriggeredEvent tev(*to,
429  timesliceRouter,
430  moduleRouter,
431  parameters.TMaxLocal_ns,
432  getTimeRange(parameters));
433 
434  // Set the event counter to the index of the Monte Carlo event in the output.
435 
436  tev.setCounter(trigger_counter);
437 
438  outputFile.put(tev);
439 
440  trigger = true;
441  }
442  }
443 
444 
445  if (!triggeredEventsOnly || trigger) {
446 
447  if (parameters.writeL0()) {
448  outputFile.put(timeslice);
449  }
450 
451  if (parameters.writeL1()) {
452  outputFile.put(JTimesliceL1<JDAQTimesliceL1>(timesliceL1, timesliceRouter, moduleRouter, parameters.TMaxLocal_ns));
453  }
454 
455  if (parameters.writeL2()) {
456  outputFile.put(JTimesliceL1<JDAQTimesliceL2>(timesliceL2, timesliceRouter, moduleRouter, parameters.L2.TMaxLocal_ns));
457  }
458 
459  if (parameters.writeSN()) {
460  outputFile.put(JTimesliceL1<JDAQTimesliceSN>(timesliceSN, timesliceRouter, moduleRouter, parameters.SN.TMaxLocal_ns));
461  }
462 
463  if (parameters.writeSummary()) {
464  outputFile.put(JSummaryslice(chronometer,simbad));
465  }
466  }
467  }
468 
469  if (!triggeredEventsOnly || trigger) {
470 
471  outputFile.put(*event);
472 
473  ++trigger_counter;
474  }
475  }
476  }
477  STATUS(endl);
478 
479 
480  JMultipleFileScanner<JMetaTypes_t> io(inputFile);
481 
482  io >> outputFile;
483 
484  outputFile.put(h1);
485  outputFile.put(*gRandom);
486  outputFile.close();
487 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
Utility class to parse command line options.
Definition: JParser.hh:1410
#define WARNING(A)
Definition: JMessage.hh:63
Data structure for all trigger parameters.
debug
Definition: JMessage.hh:27
static const unsigned int NUMBER_OF_TRIGGER_BITS
Number of trigger bits.
JBuildL2< hit_type > JBuildL2_t
Definition: JDataFilter.cc:96
Timeslice with Monte Carlo event.
#define STATUS(A)
Definition: JMessage.hh:61
Recording of objects on file according a format that follows from the file name extension.
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.
Definition: JParser.hh:64
Long64_t counter_type
Type definition for counter.
Definition: JCounter.hh:24
JSuperFrame2D< hit_type > JSuperFrame2D_t
Definition: JDataFilter.cc:93
double getTimeOfRTS(const JDAQChronometer &chronometer)
Get time of last RTS in ns since start of run for a given chronometer.
string outputFile
Data structure for UTC time.
Data structure for detector geometry and calibration.
Tools for handling different hit types.
unsigned long long int JTriggerCounter_t
Type definition of trigger counter.
double getTimeOfFrame(const int frame_index)
Get start time of frame in ns since start of run for a given frame index.
Definition: JDAQClock.hh:185
JRange< double > JTimeRange
Type definition for time range.
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Basic data structure for L0 hit.
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
JBuildL1< hit_type > JBuildL1_t
Definition: JDataFilter.cc:95
Basic data structure for time and time over threshold information of hit.
Type list.
Definition: JTypeList.hh:22
void setDAQLongprint(const bool option)
Set DAQ print option.
Definition: JDAQPrint.hh:28
Template L2 builder.
Definition: JBuildL2.hh:47
double getMaximalDistance(const JDetector &detector)
Get maximal distance between modules in detector.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
JAANET::start_run start_run
Definition: JHead.hh:851
ROOT I/O of application specific meta data.
#define NOTICE(A)
Definition: JMessage.hh:62
const JDAQChronometer & getDAQChronometer() const
Get DAQ chronometer.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Direct access to PMT in detector data structure.
int debug
debug level
Definition: JSirene.cc:59
JTimeslice< hit_type > JTimeslice_t
Definition: JDataFilter.cc:94
double getMaximalTime(const JDetector &detector)
Get maximal time between modules in detector following causality.
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:727
int run_id
MC run number.
Definition: JHead.hh:89
Template L1 hit builder.
Definition: JBuildL1.hh:76
#define FATAL(A)
Definition: JMessage.hh:65
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Auxiliaries for creation of summary data.
Utility class to parse command line options.
Map of associated modules in detector.
ROOT TTree parameter settings.
2-dimensional frame with time calibrated data from one optical module.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
Setting of trigger bits.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
Basic data structure for L1 hit.
double hit_type
Definition: JDataFilter.cc:89
Time slice with calibrated data.
Definition: JTimeslice.hh:26
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JSuperFrame1D< hit_type > JSuperFrame1D_t
Definition: JDataFilter.cc:92
int main(int argc, char *argv[])
Definition: Main.cpp:15