Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JTriggerEfficiency.cc File Reference

Auxiliary program to trigger Monte Carlo events. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include "TH1D.h"
#include "TRandom3.h"
#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "evt/Hit.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQSummaryslice.hh"
#include "JTimeslice/JEventTimeslice.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorSimulator.hh"
#include "JDetector/JModuleMapper.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JTimeRange.hh"
#include "JDetector/JK40Rates.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JK40DefaultSimulator.hh"
#include "JDetector/JPMTDefaultSimulator.hh"
#include "JDetector/JCLBDefaultSimulator.hh"
#include "JDetector/JK40RunByRunSimulator.hh"
#include "JDetector/JPMTRunByRunSimulator.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JBuildL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JTrigger3DShower.hh"
#include "JTrigger/JTriggerMXShower.hh"
#include "JTrigger/JTrigger3DMuon.hh"
#include "JTrigger/JTriggerBits.hh"
#include "JTrigger/JEventOverlap.hh"
#include "JTrigger/JTimesliceRouter.hh"
#include "JTrigger/JTriggeredEvent.hh"
#include "JTrigger/JTimesliceL1.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JEventToolkit.hh"
#include "JTrigger/JSummaryRouter.hh"
#include "JTrigger/JRunByRun.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JRandomSampler.hh"
#include "JSupport/JSummarysliceSupportkit.hh"
#include "JSupport/JMeta.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to trigger Monte Carlo events.

Note that option

Note that the event counter of the triggered events, as returned by JDAQEvent::getCounter(), are set to the corresponding index of the Monte Carlo event in the output.

Author
mdejong

Definition in file JTriggerEfficiency.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 85 of file JTriggerEfficiency.cc.

86 {
87  using namespace std;
88  using namespace JPP;
89  using namespace KM3NETDAQ;
90 
91  typedef JTYPELIST<JAllTypes_t, JRootTypes_t>::typelist typelist;
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 
135  setDAQLongprint(debug >= JEEP::debug_t);
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 }
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1410
WARNING
#define WARNING(A)
Definition: JMessage.hh:63
JTRIGGER::JTriggerParameters
Data structure for all trigger parameters.
Definition: JTriggerParameters.hh:115
JEEP::debug_t
debug
Definition: JMessage.hh:27
KM3NETDAQ::NUMBER_OF_TRIGGER_BITS
static const unsigned int NUMBER_OF_TRIGGER_BITS
Number of trigger bits.
Definition: JDAQTriggerMask.hh:31
JBuildL2_t
JBuildL2< hit_type > JBuildL2_t
Definition: JDataFilter.cc:96
KM3NETDAQ::JEventTimeslice
Timeslice with Monte Carlo event.
Definition: JEventTimeslice.hh:34
STATUS
#define STATUS(A)
Definition: JMessage.hh:61
event
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
Definition: JVoltageOptimizer_utils.hh:26
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:64
JROOT::counter_type
Long64_t counter_type
Type definition for counter.
Definition: JCounter.hh:24
JSuperFrame2D_t
JSuperFrame2D< hit_type > JSuperFrame2D_t
Definition: JDataFilter.cc:93
KM3NETDAQ::getTimeOfRTS
double getTimeOfRTS(const JDAQChronometer &chronometer)
Get time of last RTS in ns since start of run for a given chronometer.
Definition: JDAQChronometer.hh:273
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
KM3NETDAQ::JDAQUTCExtended
Data structure for UTC time.
Definition: JDAQUTCExtended.hh:26
KM3NETDAQ::JTriggerCounter_t
unsigned long long int JTriggerCounter_t
Type definition of trigger counter.
Definition: JDAQTriggerCounter.hh:21
KM3NETDAQ::getTimeOfFrame
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
JTOOLS::JTimeRange
JRange< double > JTimeRange
Type definition for time range.
Definition: JTools/JTimeRange.hh:19
JAANET::getTimeRange
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e.
Definition: JAAnetToolkit.hh:134
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:417
JSUPPORT::JLimit_t
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
JBuildL1_t
JBuildL1< hit_type > JBuildL1_t
Definition: JDataFilter.cc:95
JLANG::JTypeList
Type list.
Definition: JTypeList.hh:22
KM3NETDAQ::setDAQLongprint
void setDAQLongprint(const bool option)
Set DAQ print option.
Definition: JDAQPrint.hh:28
JTRIGGER::JBuildL2
Template L2 builder.
Definition: JBuildL2.hh:47
JDETECTOR::getMaximalDistance
double getMaximalDistance(const JDetector &detector)
Get maximal distance between modules in detector.
Definition: JDetectorToolkit.hh:73
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
KM3NETDAQ::getFrameTime
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:62
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:444
debug
int debug
debug level
Definition: JSirene.cc:59
JTimeslice_t
JTimeslice< hit_type > JTimeslice_t
Definition: JDataFilter.cc:94
JDETECTOR::getMaximalTime
double getMaximalTime(const JDetector &detector)
Get maximal time between modules in detector following causality.
Definition: JDetectorToolkit.hh:274
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:727
JTRIGGER::JBuildL1
Template L1 hit builder.
Definition: JBuildL1.hh:76
FATAL
#define FATAL(A)
Definition: JMessage.hh:65
JTRIGGER::JSuperFrame2D
2-dimensional frame with time calibrated data from one optical module.
Definition: JSuperFrame2D.hh:40
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
JSUPPORT::getTriggerParameters
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
Definition: JTriggerParametersSupportkit.hh:148
hit_type
double hit_type
Definition: JDataFilter.cc:89
JTRIGGER::JTimeslice
Time slice with calibrated data.
Definition: JTimeslice.hh:26
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JSuperFrame1D_t
JSuperFrame1D< hit_type > JSuperFrame1D_t
Definition: JDataFilter.cc:92
Generated by   doxygen 1.8.5