Jpp  17.2.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JEventTimesliceWriter.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <fstream>
4 
5 #include "TRandom3.h"
6 
11 
12 #include "JAAnet/JAAnetToolkit.hh"
13 #include "JAAnet/JEvtEvaluator.hh"
14 
16 #include "JDetector/JDetector.hh"
22 
24 #include "JDAQ/JDAQTimesliceIO.hh"
27 
29 #include "JSupport/JMeta.hh"
32 #include "JSupport/JSupport.hh"
33 #include "JSupport/JTreeScanner.hh"
35 
36 #include "Jeep/JParser.hh"
37 #include "Jeep/JMessage.hh"
38 
39 
40 /**
41  * \file
42  * Auxiliary program to convert multiple Monte Carlo events to DAQ timeslices.
43  * Events are timed according to a given rate (if > 0), otherwise they are timed according to Evt::mc_t.
44  * \author mdejong, mlincett
45  */
46 int main(int argc, char **argv)
47 {
48  using namespace std;
49  using namespace JPP;
50  using namespace KM3NETDAQ;
51 
52 
55  JLimit_t& numberOfEvents = inputFile.getLimit();
56  string detectorFile;
57  int run;
58  JPMTParametersMap pmtParameters;
59  JK40Rates rates_Hz;
60  double eventRate_Hz;
61  UInt_t seed;
62  int debug;
63 
64  try {
65 
66  JParser<> zap("Auxiliary program to convert multiple Monte Carlo events to time slices.");
67 
68  zap['f'] = make_field(inputFile);
69  zap['o'] = make_field(outputFile) = "timeslice.root";
70  zap['n'] = make_field(numberOfEvents) = JLimit::max();
71  zap['a'] = make_field(detectorFile);
72  zap['R'] = make_field(run) = 1;
73  zap['P'] = make_field(pmtParameters) = JPARSER::initialised();
74  zap['B'] = make_field(rates_Hz) = JPARSER::initialised();
75  zap['r'] = make_field(eventRate_Hz);
76  zap['S'] = make_field(seed) = 0;
77  zap['d'] = make_field(debug) = 0;
78 
79  zap(argc, argv);
80  }
81  catch(const exception &error) {
82  FATAL(error.what() << endl);
83  }
84 
85  gRandom->SetSeed(seed);
86 
87  cout.tie(&cerr);
88 
90 
91  if (pmtParameters.getQE() != 1.0) {
92 
93  WARNING("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
94 
95  rates_Hz.correct(pmtParameters.getQE());
96  }
97 
98  // parameters
99 
100  DEBUG("PMT paramaters: " << endl << pmtParameters << endl);
101  DEBUG("K40 rates: " << endl << rates_Hz << endl);
102 
103  if (eventRate_Hz < 0.0) {
104  FATAL("Invalid event rate " << eventRate_Hz << "; consider using JRandomTimesliceWriter." << endl);
105  }
106 
107  // detector
108 
110 
111  try {
112  load(detectorFile, detector);
113  }
114  catch(const JException& error) {
115  FATAL(error);
116  }
117 
118  // header
119 
120  Head header;
121 
122  try {
123  header = getHeader(inputFile);
124  }
125  catch (const exception& error) {
126  FATAL("Monte Carlo header is invalid");
127  }
128 
129  double liveTimeMC = JHead(header).livetime.numberOfSeconds;
130 
131  if (liveTimeMC < 0) {
132  FATAL("Monte Carlo live time is negative; input file may be corrupted.");
133  }
134 
135  // background simulator
136 
137  JPMTParametersMap::Throw(false);
138 
140 
141  simbad.reset(new JPMTDefaultSimulator(pmtParameters, detector));
142  simbad.reset(new JK40DefaultSimulator(rates_Hz));
143  simbad.reset(new JCLBDefaultSimulator());
144 
145  // output file
146 
147  outputFile.open();
148 
149  if (!outputFile.is_open()) {
150  FATAL("Error opening file " << outputFile << endl);
151  }
152 
153  outputFile.put(*gRandom);
154  outputFile.put(JMeta(argc, argv));
155  outputFile.put(header);
156 
157  // input stream;
158 
159  bool absTime = false;
160 
162 
163  if (eventRate_Hz == 0.0 && liveTimeMC == 0.0) {
164 
165  NOTICE("Event will be timed according to absolute MC time." << endl);
166 
167  absTime = true;
168 
169  // sorted access;
170  scan = new JTreeScanner<Evt, JEvtEvaluator>(inputFile);
172 
173  } else {
174 
175  // unsorted access;
176  scan = new JTreeScanner<Evt>(inputFile);
178 
179  if (eventRate_Hz == 0.0 && liveTimeMC > 0.0) {
180  int nEntries = scan->getEntries();
181  eventRate_Hz = nEntries / liveTimeMC;
182  DEBUG(nEntries << " events to be written." << endl);
183  NOTICE("MC live time is " << liveTimeMC << " s" << endl);
184  NOTICE("Event rate set to " << eventRate_Hz << " Hz from MC header." << endl);
185  }
186  }
187 
188  // begin timeslice generation;
189 
190  // state variables;
191  bool pendingEvt = false;
192  bool pendingTimeslice = false;
193 
194  Evt* event = NULL; // pending event;
195 
196  JTimeRange timeRange;
197  double tDAQ = 0; // DAQ event time;
198  double tOff = 0; // hit production time offset;
199 
200  int frame_index = 0;
201  JDAQChronometer chronometer;
202  JRandomTimeslice timeslice;
203 
204  counter_type evtCount = 0;
205 
206  while (pendingEvt || scan->hasNext()) {
207 
208  if (!pendingTimeslice) {
209  // update background timeslice;
210  frame_index++;
211  chronometer = JDAQChronometer(detector.getID(), run, frame_index, JDAQUTCExtended(getTimeOfFrame(frame_index)));
212  timeslice = JRandomTimeslice(chronometer, simbad);
213  DEBUG("evt count: " << setw(10) << evtCount << endl);
214  STATUS("frame index: " << setw(10) << frame_index << " | evt count: " << setw(10) << evtCount << "\r"); DEBUG(endl);
215  pendingTimeslice = true;
216  }
217 
218  if (!pendingEvt) {
219  // get new event;
220  event = scan->next();
221  timeRange = getTimeRange(*event);
222  tOff = timeRange.is_valid() ? timeRange.getLowerLimit() : 0;
223  // update DAQ time according to absolute or random;
224  if (absTime)
225  tDAQ = getEvtValue(*event) + tOff;
226  else
227  tDAQ += gRandom->Exp(1.0e9 / eventRate_Hz);
228 
229  DEBUG("event time [s] " << setprecision(5) << tDAQ * 1.0e-9 << endl);
230  pendingEvt = true;
231  }
232 
233  // DEBUG("tDAQ : " << tDAQ << " - frame index: " << frame_index << " - DAQ frame index " << getFrameIndex(tDAQ) << endl);
234  if (getFrameIndex(tDAQ) == frame_index) {
235 
236  // process event;
237  if (timeRange.is_valid()) {
238  DEBUG(*event << endl);
239  // real start time of the hits will be always tDAQ;
240  event->mc_t = tDAQ - tOff;
241  timeslice.add(JEventTimeslice(chronometer, simbad, *event));
242  evtCount++;
243  }
244  // event processed;
245  outputFile.put(*event);
246  pendingEvt = false;
247  } else {
248  // event time is past end of timeslice;
249  DEBUG(timeslice << endl);
250  outputFile.put(timeslice);
251  pendingTimeslice = false;
252  }
253 
254  }
255 
256  if (pendingTimeslice) {
257  DEBUG(timeslice << endl);
258  outputFile.put(timeslice);
259  pendingTimeslice = false;
260  }
261 
262  STATUS(endl);
263 
264  outputFile.put(*gRandom);
265  outputFile.close();
266 
267  NOTICE(evtCount << " events written over " << frame_index << " timeslices. " << endl);
268 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1517
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
debug
Definition: JMessage.hh:29
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Default implementation of the simulation of K40 background.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Recording of objects on file according a format that follows from the file name extension.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
Long64_t counter_type
Type definition for counter.
string outputFile
Data structure for UTC time.
Template definition for direct access of elements in ROOT TChain.
Data structure for detector geometry and calibration.
Auxiliary interface for direct access of elements in ROOT TChain.
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
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e. time between earliest and latest hit) of Monte Carlo event.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
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
Detector file.
Definition: JHead.hh:226
static const JEvtEvaluator getEvtValue
Function object for evaluation of DAQ objects.
JAANET::livetime livetime
Definition: JHead.hh:1529
int getFrameIndex(const double t_ns)
Get frame index for a given time in ns.
Definition: JDAQClock.hh:251
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
ROOT I/O of application specific meta data.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:64
Auxiliary class for map of PMT parameters.
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:1167
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
void reset(JK40Simulator *k40Simulator)
Reset K40 simulator.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
do set_variable DETECTOR_TXT $WORKDIR detector
double numberOfSeconds
Live time [s].
Definition: JHead.hh:883
Timeslice with random data.
Timeslice with Monte Carlo event.
int debug
debug level
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62