Jpp  15.0.5
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 
10 
11 #include "JAAnet/JAAnetToolkit.hh"
12 #include "JAAnet/JEvtEvaluator.hh"
13 
15 #include "JDetector/JDetector.hh"
21 
23 #include "JDAQ/JDAQTimesliceIO.hh"
26 
28 #include "JSupport/JMeta.hh"
31 #include "JSupport/JSupport.hh"
32 #include "JSupport/JTreeScanner.hh"
34 
35 #include "Jeep/JParser.hh"
36 #include "Jeep/JMessage.hh"
37 
38 
39 /**
40  * \file
41  * Auxiliary program to convert multiple Monte Carlo events to DAQ timeslices.
42  * Events are timed according to a given rate (if > 0), otherwise they are timed according to Evt::mc_t.
43  * \author mdejong, mlincett
44  */
45 int main(int argc, char **argv)
46 {
47  using namespace std;
48  using namespace JPP;
49  using namespace KM3NETDAQ;
50 
51 
54  JLimit_t& numberOfEvents = inputFile.getLimit();
55  string detectorFile;
56  int run;
57  JPMTParametersMap pmtParameters;
58  JK40Rates rates_Hz;
59  double eventRate_Hz;
60  UInt_t seed;
61  int debug;
62 
63  try {
64 
65  JParser<> zap("Auxiliary program to convert multiple Monte Carlo events to time slices.");
66 
67  zap['f'] = make_field(inputFile);
68  zap['o'] = make_field(outputFile) = "timeslice.root";
69  zap['n'] = make_field(numberOfEvents) = JLimit::max();
70  zap['a'] = make_field(detectorFile);
71  zap['R'] = make_field(run) = 1;
72  zap['P'] = make_field(pmtParameters) = JPARSER::initialised();
73  zap['B'] = make_field(rates_Hz) = JPARSER::initialised();
74  zap['r'] = make_field(eventRate_Hz);
75  zap['S'] = make_field(seed) = 0;
76  zap['d'] = make_field(debug) = 0;
77 
78  zap(argc, argv);
79  }
80  catch(const exception &error) {
81  FATAL(error.what() << endl);
82  }
83 
84  gRandom->SetSeed(seed);
85 
86  cout.tie(&cerr);
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:1500
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:66
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.
Definition: JTreeScanner.hh:91
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:196
static const JEvtEvaluator getEvtValue
Function object for evaluation of DAQ objects.
JAANET::livetime livetime
Definition: JHead.hh:1412
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:1961
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.
int debug
debug level
Definition: JSirene.cc:63
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:1113
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:67
#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.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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:829
Timeslice with random data.
Timeslice with Monte Carlo event.
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:19