Jpp  18.5.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 
89 
90  if (pmtParameters.getQE() != 1.0) {
91 
92  WARNING("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
93 
94  rates_Hz.correct(pmtParameters.getQE());
95  }
96 
97  // parameters
98 
99  DEBUG("PMT paramaters: " << endl << pmtParameters << endl);
100  DEBUG("K40 rates: " << endl << rates_Hz << endl);
101 
102  if (eventRate_Hz < 0.0) {
103  FATAL("Invalid event rate " << eventRate_Hz << "; consider using JRandomTimesliceWriter." << endl);
104  }
105 
106  // detector
107 
109 
110  try {
111  load(detectorFile, detector);
112  }
113  catch(const JException& error) {
114  FATAL(error);
115  }
116 
117  // header
118 
119  Head header;
120 
121  try {
122  header = getHeader(inputFile);
123  }
124  catch (const exception& error) {
125  FATAL("Monte Carlo header is invalid");
126  }
127 
128  double liveTimeMC = JHead(header).livetime.numberOfSeconds;
129 
130  if (liveTimeMC < 0) {
131  FATAL("Monte Carlo live time is negative; input file may be corrupted.");
132  }
133 
134  // background simulator
135 
136  JPMTParametersMap::Throw(false);
137 
139 
140  simbad.reset(new JPMTDefaultSimulator(pmtParameters, detector));
141  simbad.reset(new JK40DefaultSimulator(rates_Hz));
142  simbad.reset(new JCLBDefaultSimulator());
143 
144  // output file
145 
146  outputFile.open();
147 
148  if (!outputFile.is_open()) {
149  FATAL("Error opening file " << outputFile << endl);
150  }
151 
152  outputFile.put(*gRandom);
153  outputFile.put(JMeta(argc, argv));
154  outputFile.put(header);
155 
156  // input stream;
157 
158  bool absTime = false;
159 
161 
162  if (eventRate_Hz == 0.0 && liveTimeMC == 0.0) {
163 
164  NOTICE("Event will be timed according to absolute MC time." << endl);
165 
166  absTime = true;
167 
168  // sorted access;
169  scan = new JTreeScanner<Evt, JEvtEvaluator>(inputFile);
171 
172  } else {
173 
174  // unsorted access;
175  scan = new JTreeScanner<Evt>(inputFile);
177 
178  if (eventRate_Hz == 0.0 && liveTimeMC > 0.0) {
179  int nEntries = scan->getEntries();
180  eventRate_Hz = nEntries / liveTimeMC;
181  DEBUG(nEntries << " events to be written." << endl);
182  NOTICE("MC live time is " << liveTimeMC << " s" << endl);
183  NOTICE("Event rate set to " << eventRate_Hz << " Hz from MC header." << endl);
184  }
185  }
186 
187  // begin timeslice generation;
188 
189  // state variables;
190  bool pendingEvt = false;
191  bool pendingTimeslice = false;
192 
193  Evt* event = NULL; // pending event;
194 
195  JTimeRange timeRange;
196  double tDAQ = 0; // DAQ event time;
197  double tOff = 0; // hit production time offset;
198 
199  int frame_index = 0;
200  JDAQChronometer chronometer;
201  JRandomTimeslice timeslice;
202 
203  counter_type evtCount = 0;
204 
205  while (pendingEvt || scan->hasNext()) {
206 
207  if (!pendingTimeslice) {
208  // update background timeslice;
209  frame_index++;
210  chronometer = JDAQChronometer(detector.getID(), run, frame_index, JDAQUTCExtended(getTimeOfFrame(frame_index)));
211  timeslice = JRandomTimeslice(chronometer, simbad);
212  DEBUG("evt count: " << setw(10) << evtCount << endl);
213  STATUS("frame index: " << setw(10) << frame_index << " | evt count: " << setw(10) << evtCount << "\r"); DEBUG(endl);
214  pendingTimeslice = true;
215  }
216 
217  if (!pendingEvt) {
218  // get new event;
219  event = scan->next();
220  timeRange = getTimeRange(*event);
221  tOff = timeRange.is_valid() ? timeRange.getLowerLimit() : 0;
222  // update DAQ time according to absolute or random;
223  if (absTime)
224  tDAQ = getEvtValue(*event) + tOff;
225  else
226  tDAQ += gRandom->Exp(1.0e9 / eventRate_Hz);
227 
228  DEBUG("event time [s] " << setprecision(5) << tDAQ * 1.0e-9 << endl);
229  pendingEvt = true;
230  }
231 
232  // DEBUG("tDAQ : " << tDAQ << " - frame index: " << frame_index << " - DAQ frame index " << getFrameIndex(tDAQ) << endl);
233  if (getFrameIndex(tDAQ) == frame_index) {
234 
235  // process event;
236  if (timeRange.is_valid()) {
237  DEBUG(*event << endl);
238  // real start time of the hits will be always tDAQ;
239  event->mc_t = tDAQ - tOff;
240  timeslice.add(JEventTimeslice(chronometer, simbad, *event));
241  evtCount++;
242  }
243  // event processed;
244  outputFile.put(*event);
245  pendingEvt = false;
246  } else {
247  // event time is past end of timeslice;
248  DEBUG(timeslice << endl);
249  outputFile.put(timeslice);
250  pendingTimeslice = false;
251  }
252 
253  }
254 
255  if (pendingTimeslice) {
256  DEBUG(timeslice << endl);
257  outputFile.put(timeslice);
258  pendingTimeslice = false;
259  }
260 
261  STATUS(endl);
262 
263  outputFile.put(*gRandom);
264  outputFile.close();
265 
266  NOTICE(evtCount << " events written over " << frame_index << " timeslices. " << endl);
267 }
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:1514
General exception.
Definition: JException.hh:24
#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.
Auxiliary interface for direct access of elements in ROOT TChain.
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.
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
Monte Carlo run header.
Definition: JHead.hh:1234
static const JEvtEvaluator getEvtValue
Function object for evaluation of DAQ objects.
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:1989
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.
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:84
do set_variable DETECTOR_TXT $WORKDIR detector
double numberOfSeconds
Live time [s].
Definition: JHead.hh:896
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
JAANET::livetime livetime
Definition: JHead.hh:1601