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 "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JTimeslice/JEventTimeslice.hh"
#include "JSummaryslice/JSummaryslice.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.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 "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 "JTrigger/JK40RunByRunSimulator.hh"
#include "JTrigger/JPMTRunByRunSimulator.hh"
#include "JTrigger/JCLBRunByRunSimulator.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/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 KM3NETDAQ::JDAQEvent::getCounter(), are set to the corresponding index of the Monte Carlo event in the output.

In run-by-run mode, the trigger parameters as well as the singles rates and status of the PMTs are read from the given raw data file. The two-fold (and higher) coincidence rates should still be provided on the command line.

Author
mdejong

Definition in file JTriggerEfficiency.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 92 of file JTriggerEfficiency.cc.

93 {
94  using namespace std;
95  using namespace JPP;
96  using namespace KM3NETDAQ;
97 
99 
100  JMultipleFileScanner<> inputFile;
102  JLimit_t& numberOfEvents = inputFile.getLimit();
103  string detectorFileA;
104  string detectorFileB;
105  int run;
107  bool triggeredEventsOnly;
108  JPMTParametersMap pmtParameters;
109  JK40Rates rates_Hz;
110  JRunByRun runbyrun;
111  UInt_t seed;
112  int debug;
113 
114  try {
115 
116  JParser<> zap("Auxiliary program to trigger Monte Carlo events.");
117 
118  zap['f'] = make_field(inputFile, "input file (output of detector simulation)");
119  zap['o'] = make_field(outputFile, "output file") = "trigger_efficieny.root";
120  zap['n'] = make_field(numberOfEvents) = JLimit::max();
121  zap['a'] = make_field(detectorFileA, "detector used for conversion from Monte Carlo truth to raw data.");
122  zap['b'] = make_field(detectorFileB, "detector used for conversion of raw data to calibrated data.") = "";
123  zap['R'] = make_field(run, "run number") = -1;
124  zap['r'] = make_field(runbyrun, "option for run-by-run mode") = JPARSER::initialised();
125  zap['@'] = make_field(parameters, "Trigger parameters (or corresponding file name)") = JPARSER::initialised();
126  zap['O'] = make_field(triggeredEventsOnly, "optionally write only triggered events.");
127  zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised();
128  zap['B'] = make_field(rates_Hz, "background rates [Hz]") = JPARSER::initialised();
129  zap['S'] = make_field(seed, "seed") = 0;
130  zap['d'] = make_field(debug, "debug") = 0;
131 
132  zap(argc, argv);
133  }
134  catch(const exception &error) {
135  FATAL(error.what() << endl);
136  }
137 
138  gRandom->SetSeed(seed);
139 
140  cout.tie(&cerr);
141 
143 
144  if (detectorFileB == "") {
145  detectorFileB = detectorFileA;
146  }
147 
148 
149  JDetector detectorA;
150  JDetector detectorB;
151 
152  try {
153  load(detectorFileA, detectorA);
154  load(detectorFileB, detectorB);
155  }
156  catch(const JException& error) {
157  FATAL(error);
158  }
159 
160  JPMTParametersMap::Throw(true);
161 
162  if (!pmtParameters.is_valid()) {
163  FATAL("Invalid PMT parameters " << pmtParameters << endl);
164  }
165 
166  if (pmtParameters.getQE() != 1.0) {
167 
168  WARNING("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
169 
170  rates_Hz.correct(pmtParameters.getQE());
171  }
172 
173  const JModuleRouter moduleRouter(detectorB);
174  JDetectorSimulator simbad (detectorA);
175  JSummaryRouter summaryRouter;
176 
177  if (runbyrun.is_valid()) {
178 
179  NOTICE("Using run-by-run:" << endl << runbyrun << endl);
180 
181  if (!runbyrun.hasNext()) {
182  FATAL("Run-by-run simulation yields no input." << endl);
183  }
184 
185  if (rates_Hz.getSinglesRate() != 0.0) {
186  WARNING("Run-by-run simulation discards singles rate [Hz] " << rates_Hz.getSinglesRate() << endl);
187  }
188 
189  try {
190  simbad.reset(new JK40RunByRunSimulator(summaryRouter, rates_Hz));
191  simbad.reset(new JPMTRunByRunSimulator(summaryRouter, pmtParameters, detectorA));
192  simbad.reset(new JCLBRunByRunSimulator(summaryRouter));
193  }
194  catch(const JException& error) {
195  FATAL(error.what() << endl);
196  }
197 
198  try {
199 
201 
202  NOTICE("Set trigger parameters from run-by-run input." << endl);
203  }
204  catch(const JException& error) {
205  WARNING("No trigger parameters from run-by-run input;\nrun with default/user input." << endl);
206  }
207 
208  } else {
209 
210  NOTICE("Using fixed rates [Hz]: " << rates_Hz << endl);
211 
212  try {
213  simbad.reset(new JK40DefaultSimulator(rates_Hz));
214  simbad.reset(new JPMTDefaultSimulator(pmtParameters, detectorA));
215  simbad.reset(new JCLBDefaultSimulator());
216  }
217  catch(const JException& error) {
218  FATAL(error.what() << endl);
219  }
220  }
221 
222  parameters.set(getMaximalDistance(detectorB));
223 
224  DEBUG("Trigger:" << endl << parameters << endl);
225  DEBUG("PMT parameters:" << endl << pmtParameters << endl);
226 
227  const double Tmax = max(getMaximalTime(detectorA),
228  getMaximalTime(detectorB));
229 
230  const JTimeRange period(-(Tmax + parameters.TMaxLocal_ns),
231  +(Tmax + parameters.TMaxLocal_ns));
232 
233  typedef double hit_type;
234 
235  typedef JSuperFrame1D<hit_type> JSuperFrame1D_t;
236  typedef JSuperFrame2D<hit_type> JSuperFrame2D_t;
237  typedef JTimeslice <hit_type> JTimeslice_t;
238  typedef JBuildL1 <hit_type> JBuildL1_t;
239  typedef JBuildL2 <hit_type> JBuildL2_t;
240 
241  const JBuildL1_t buildL1(parameters);
242  const JBuildL2_t buildL2(parameters.L2);
243  const JBuildL2_t buildSN(parameters.SN);
244 
245  JTimesliceRouter timesliceRouter(parameters.numberOfBins);
246 
247  const JTrigger3DMuon trigger3DMuon (parameters);
248  const JTrigger3DShower trigger3DShower(parameters);
249  const JTriggerMXShower triggerMXShower(parameters, detectorB);
250 
251 
252  Head header;
253 
254  try {
255  header = getHeader(inputFile);
256  }
257  catch(const JException& error) {
258  FATAL(error);
259  }
260 
261  const JPosition3D center = get<JPosition3D>(header);
262 
263  NOTICE("Apply detector offset from Monte Carlo run header (" << center << ")" << endl);
264 
265  detectorA -= center;
266  detectorB -= center;
267 
268  TH1D h1("Trigger bits", NULL, NUMBER_OF_TRIGGER_BITS, -0.5, NUMBER_OF_TRIGGER_BITS - 0.5);
269 
270 
271  outputFile.open();
272 
273  if (!outputFile.is_open()) {
274  FATAL("Error opening file " << outputFile << endl);
275  }
276 
277  outputFile.put(*gRandom);
278  outputFile.put(JMeta(argc, argv));
279  outputFile.put(header);
280  outputFile.put(parameters);
281 
282  JLimit_t limit = inputFile.getLimit();
283  counter_type number_of_events = 0;
284  JTriggerCounter_t trigger_counter = 0;
285 
286  for (JMultipleFileScanner<>::const_iterator file = inputFile.begin(); file != inputFile.end(); ++file) {
287 
288  int mc_run_id = 0;
289 
290  try {
291 
292  const JHead head(getHeader(*file));
293 
294  mc_run_id = head.start_run.run_id;
295  }
296  catch(const JException& error) {
297  FATAL(error);
298  }
299 
301 
302  limit.setLowerLimit(limit.getLowerLimit() - in.skip(limit.getLowerLimit()));
303 
304  for ( ; in.hasNext() && number_of_events != limit; ++number_of_events) {
305 
306  STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
307 
308  Evt* event = in.next();
309 
310  event->mc_run_id = mc_run_id;
311 
312  DEBUG(*event << endl);
313 
314  bool trigger = false;
315 
316  if (!event->mc_hits.empty()) {
317 
318  int frame_index = (int) in.getCounter() + 1; // 1 event / slice
319 
320  if (runbyrun.is_valid() && runbyrun.hasNext()) {
321 
322  summaryRouter.update(runbyrun.next());
323 
324  summaryRouter.correct(dynamic_cast<const JPMTDefaultSimulatorInterface&>(simbad.getPMTSimulator()));
325 
326  frame_index = summaryRouter.getFrameIndex();
327  run = summaryRouter.getRunNumber();
328  }
329 
330  // set the event time!
331 
332  JTimeRange timeRange = getTimeRange(*event, period);
333 
334  if (!timeRange.is_valid()) {
335  timeRange.setRange(0.0,0.0);
336  }
337 
338  const double t0 = 0.5 * (timeRange.getLowerLimit() + timeRange.getUpperLimit());
339  const double t1 = getTimeOfFrame(frame_index) + gRandom->Rndm() * getFrameTime();
340 
341  event->mc_t = t1 - t0; // time since start of data taking run
342 
343  timeRange.add(event->mc_t);
344  timeRange.add(period);
345 
346  const JDAQChronometer chronometer(detectorB.getID(), (run != -1 ? run : mc_run_id) , frame_index, JDAQUTCExtended(getTimeOfFrame(frame_index)));
347 
348  const JEventTimeslice timeslice(chronometer, simbad, *event, period);
349 
350  DEBUG(timeslice << endl);
351 
352 
353  timesliceRouter.configure(timeslice);
354 
355  JTimeslice_t timesliceL0(timeslice.getDAQChronometer());
356  JTimeslice_t timesliceL1(timeslice.getDAQChronometer());
357  JTimeslice_t timesliceL2(timeslice.getDAQChronometer());
358  JTimeslice_t timesliceSN(timeslice.getDAQChronometer());
359 
360  for (JDAQTimeslice::const_iterator super_frame = timeslice.begin(); super_frame != timeslice.end(); ++super_frame) {
361 
362  if (moduleRouter.hasModule(super_frame->getModuleID())) {
363 
364  // calibration
365 
366  const JModule& module = moduleRouter.getModule(super_frame->getModuleID());
367  const JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, module);
368 
369  // L0
370 
371  timesliceL0.push_back(JSuperFrame1D_t(buffer));
372 
373  // L1
374 
375  timesliceL1.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
376  super_frame->getModuleIdentifier(),
377  module.getPosition()));
378 
379  buildL1(*timesliceL0.rbegin(), back_inserter(*timesliceL1.rbegin()));
380 
381  // L2
382 
383  timesliceL2.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
384  super_frame->getModuleIdentifier(),
385  module.getPosition()));
386 
387  buildL2(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceL2.rbegin()));
388 
389  // SN
390 
391  timesliceSN.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
392  super_frame->getModuleIdentifier(),
393  module.getPosition()));
394 
395  buildSN(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceSN.rbegin()));
396 
397  DEBUG("L0 " << setw(8) << timesliceL0.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL0.rbegin()->size() << endl);
398  DEBUG("L1 " << setw(8) << timesliceL1.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL1.rbegin()->size() << endl);
399  DEBUG("L2 " << setw(8) << timesliceL2.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL2.rbegin()->size() << endl);
400  DEBUG("SN " << setw(8) << timesliceSN.rbegin()->getModuleID() << ' ' << setw(8) << timesliceSN.rbegin()->size() << endl);
401  }
402  }
403 
404 
405  // Trigger
406 
407  JTriggerInput trigger_input(timesliceL2);
408  JTriggerOutput trigger_output;
409 
410  trigger3DMuon (trigger_input, back_inserter(trigger_output));
411  trigger3DShower(trigger_input, back_inserter(trigger_output));
412  triggerMXShower(trigger_input, timesliceL0, back_inserter(trigger_output));
413 
414  trigger_output.merge(JEventOverlap(parameters.TMaxEvent_ns));
415 
416  for (JTriggerOutput::const_iterator to = trigger_output.begin(); to != trigger_output.end(); ++to) {
417 
418  for (int i = 0; i != h1.GetNbinsX(); ++i) {
419  if (to->hasTriggerBit(i)) {
420  h1.Fill((double) i);
421  }
422  }
423 
424  JTimeRange eventTime = getTimeRange(*to).add(getTimeOfRTS(*to));
425 
426  DEBUG("Event time: "
427  << to->getFrameIndex() << ' '
428  << eventTime << ' '
429  << timeRange << endl);
430 
431  if (timeRange.overlap(eventTime)) {
432 
433  JTriggeredEvent tev(*to,
434  timesliceRouter,
435  moduleRouter,
436  parameters.TMaxLocal_ns,
438 
439  // Set the event counter to the index of the Monte Carlo event in the output.
440 
441  tev.setCounter(trigger_counter);
442 
443  outputFile.put(tev);
444 
445  trigger = true;
446  }
447  }
448 
449 
450  if (!triggeredEventsOnly || trigger) {
451 
452  if (parameters.writeL0()) {
453  outputFile.put(timeslice);
454  }
455 
456  if (parameters.writeL1()) {
457  outputFile.put(JTimesliceL1<JDAQTimesliceL1>(timesliceL1, timesliceRouter, moduleRouter, parameters.TMaxLocal_ns));
458  }
459 
460  if (parameters.writeL2()) {
461  outputFile.put(JTimesliceL1<JDAQTimesliceL2>(timesliceL2, timesliceRouter, moduleRouter, parameters.L2.TMaxLocal_ns));
462  }
463 
464  if (parameters.writeSN()) {
465  outputFile.put(JTimesliceL1<JDAQTimesliceSN>(timesliceSN, timesliceRouter, moduleRouter, parameters.SN.TMaxLocal_ns));
466  }
467 
468  if (parameters.writeSummary()) {
469  outputFile.put(JSummaryslice(chronometer,simbad));
470  }
471  }
472  }
473 
474  if (!triggeredEventsOnly || trigger) {
475 
476  outputFile.put(*event);
477 
478  ++trigger_counter;
479  }
480  }
481  }
482  STATUS(endl);
483 
484 
486 
487  io >> outputFile;
488 
489  outputFile.put(h1);
490  outputFile.put(*gRandom);
491  outputFile.close();
492 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1493
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
debug
Definition: JMessage.hh:29
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:215
static const unsigned int NUMBER_OF_TRIGGER_BITS
Number of trigger bits.
Default implementation of the simulation of K40 background.
Data structure for a composite optical module.
Definition: JModule.hh:50
Auxiliary class to handle run by run options.
Definition: JRunByRun.hh:47
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
Router for direct addressing of module data in detector data structure.
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Definition: JSirene.sh:45
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:63
Long64_t counter_type
Type definition for counter.
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.
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
1-dimensional frame with time calibrated data from one optical module.
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e. time between earliest and latest hit) of Monte Carlo event.
K40 simulation based on run-by-run information.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Type list.
Definition: JTypeList.hh:22
void setDAQLongprint(const bool option)
Set DAQ print option.
Definition: JDAQPrint.hh:28
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Template L2 builder.
Definition: JBuildL2.hh:45
double getMaximalDistance(const JDetector &detector)
Get maximal distance between modules in detector.
void merge(const JMatch_t &match)
Merge events.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
CLB simulation based on run-by-run information.
Router for fast addressing of hits in KM3NETDAQ::JDAQTimeslice data structure as a function of the op...
range_type & add(argument_type x)
Add offset.
Definition: JRange.hh:448
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
void setRange(const range_type &range)
Set range.
Definition: JRange.hh:139
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:226
#define NOTICE(A)
Definition: JMessage.hh:64
Auxiliary class to create summary data.
Auxiliary class to build JDAQEvent for a triggered event.
Auxiliary class for map of PMT parameters.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:61
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
double getMaximalTime(const double R_Hz)
Get maximal time for given rate.
Router for fast addressing of summary data in JDAQSummaryslice data structure as a function of the op...
Monte Carlo run header.
Definition: JHead.hh:836
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:66
Template L1 hit builder.
Definition: JBuildL1.hh:85
#define FATAL(A)
Definition: JMessage.hh:67
PMT simulation based on run-by-run information.
bool is_valid() const
Check validity of range.
Definition: JRange.hh:324
Auxiliary class to build JDAQTimeslice for L1 timeslice.
Definition: JTimesliceL1.hh:36
2-dimensional frame with time calibrated data from one optical module.
Data structure for input to trigger algorithm.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
virtual const char * what() const
Get error message.
Definition: JException.hh:48
Match of two events considering overlap in time.
bool overlap(const range_type &range) const
Test overlap with given range.
Definition: JRange.hh:384
Timeslice with Monte Carlo event.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
Time slice with calibrated data.
Definition: JTimeslice.hh:26
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41