11 #include "evt/Head.hh"
85 int main(
int argc,
char **argv)
95 JLimit_t& numberOfEvents = inputFile.getLimit();
100 bool triggeredEventsOnly;
101 JPMTParametersMap pmtParameters;
109 JParser<> zap(
"Auxiliary program to trigger Monte Carlo events.");
111 zap[
'f'] =
make_field(inputFile,
"input file (output of detector simulation)");
113 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
114 zap[
'a'] =
make_field(detectorFileA,
"detector used for conversion 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;
119 zap[
'O'] =
make_field(triggeredEventsOnly,
"optionally write only triggered events.");
127 catch(
const exception &error) {
128 FATAL(error.what() << endl);
131 gRandom->SetSeed(seed);
137 if (detectorFileB ==
"") {
138 detectorFileB = detectorFileA;
146 load(detectorFileA, detectorA);
147 load(detectorFileB, detectorB);
149 catch(
const JException& error) {
153 JPMTParametersMap::Throw(
true);
155 if (!pmtParameters.is_valid()) {
156 FATAL(
"Invalid PMT parameters " << pmtParameters << endl);
159 if (pmtParameters.getQE() != 1.0) {
161 WARNING(
"Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
163 rates_Hz.correct(pmtParameters.getQE());
166 const JModuleRouter moduleRouter(detectorB);
167 JDetectorSimulator simbad (detectorA);
168 JSummaryRouter summaryRouter;
170 if (runbyrun.is_valid()) {
172 NOTICE(
"Using run-by-run:" << endl << runbyrun << endl);
174 if (!runbyrun.hasNext()) {
175 FATAL(
"Run-by-run simulation yields no input." << endl);
178 if (rates_Hz.getSinglesRate() != 0.0) {
179 WARNING(
"Run-by-run simulation discards singles rate [Hz] " << rates_Hz.getSinglesRate() << endl);
183 simbad.reset(
new JK40RunByRunSimulator(summaryRouter, rates_Hz));
184 simbad.reset(
new JPMTRunByRunSimulator(summaryRouter, pmtParameters, detectorA, runbyrun.range_Hz));
185 simbad.reset(
new JCLBDefaultSimulator());
187 catch(
const JException& error) {
188 FATAL(error.what() << endl);
195 NOTICE(
"Set trigger parameters from run-by-run input." << endl);
197 catch(
const JException& error) {
198 WARNING(
"No trigger parameters from run-by-run input;\nrun with default/user input." << endl);
203 NOTICE(
"Using fixed rates [Hz]: " << rates_Hz << endl);
206 simbad.reset(
new JK40DefaultSimulator(rates_Hz));
207 simbad.reset(
new JPMTDefaultSimulator(pmtParameters, detectorA));
208 simbad.reset(
new JCLBDefaultSimulator());
210 catch(
const JException& error) {
211 FATAL(error.what() << endl);
217 DEBUG(
"Trigger:" << endl << parameters << endl);
218 DEBUG(
"PMT parameters:" << endl << pmtParameters << endl);
225 typedef double hit_type;
227 typedef JSuperFrame1D<hit_type> JSuperFrame1D_t;
228 typedef JSuperFrame2D<hit_type> JSuperFrame2D_t;
229 typedef JTimeslice <hit_type> JTimeslice_t;
230 typedef JBuildL1 <hit_type> JBuildL1_t;
231 typedef JBuildL2 <hit_type> JBuildL2_t;
233 const JBuildL1_t buildL1(parameters);
234 const JBuildL2_t buildL2(parameters.
L2);
235 const JBuildL2_t buildSN(parameters.
SN);
237 JTimesliceRouter timesliceRouter(parameters.
numberOfBins);
239 const JTrigger3DMuon trigger3DMuon (parameters);
240 const JTrigger3DShower trigger3DShower(parameters);
241 const JTriggerMXShower triggerMXShower(parameters, detectorB);
249 catch(
const JException& error) {
253 const JPosition3D center = get<JPosition3D>(header);
255 NOTICE(
"Apply detector offset from Monte Carlo run header (" << center <<
")" << endl);
274 JLimit_t limit = inputFile.getLimit();
278 for (JMultipleFileScanner<>::const_iterator file = inputFile.begin(); file != inputFile.end(); ++file) {
288 catch(
const JException& error) {
292 JMultipleFileScanner<Evt> in(*file);
294 limit.setLowerLimit(limit.getLowerLimit() - in.skip(limit.getLowerLimit()));
296 for ( ; in.hasNext() && number_of_events != limit; ++number_of_events) {
298 STATUS(
"event: " << setw(10) << number_of_events <<
'\r');
DEBUG(endl);
300 Evt*
event = in.next();
302 event->mc_run_id = mc_run_id;
304 DEBUG(*event << endl);
306 bool trigger =
false;
308 if (!event->mc_hits.empty()) {
310 int frame_index = (int) in.getCounter() + 1;
312 if (runbyrun.is_valid() && runbyrun.hasNext()) {
314 summaryRouter.update(runbyrun.next());
316 summaryRouter.correct(dynamic_cast<const JPMTDefaultSimulatorInterface&>(simbad.getPMTSimulator()));
318 frame_index = summaryRouter.getFrameIndex();
319 run = summaryRouter.getRunNumber();
326 if (!timeRange.is_valid()) {
327 timeRange.setRange(0.0,0.0);
330 const double t0 = 0.5 * (timeRange.getLowerLimit() + timeRange.getUpperLimit());
333 event->mc_t = t1 - t0;
335 timeRange.add(event->mc_t);
336 timeRange.add(period);
342 DEBUG(timeslice << endl);
345 timesliceRouter.configure(timeslice);
352 for (JDAQTimeslice::const_iterator super_frame = timeslice.begin(); super_frame != timeslice.end(); ++super_frame) {
354 if (moduleRouter.hasModule(super_frame->getModuleID())) {
358 const JModule& module = moduleRouter.getModule(super_frame->getModuleID());
359 const JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, module);
363 timesliceL0.push_back(JSuperFrame1D_t(buffer));
367 timesliceL1.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
368 super_frame->getModuleIdentifier(),
369 module.getPosition()));
371 buildL1(*timesliceL0.rbegin(), back_inserter(*timesliceL1.rbegin()));
375 timesliceL2.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
376 super_frame->getModuleIdentifier(),
377 module.getPosition()));
379 buildL2(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceL2.rbegin()));
383 timesliceSN.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
384 super_frame->getModuleIdentifier(),
385 module.getPosition()));
387 buildSN(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceSN.rbegin()));
389 DEBUG(
"L0 " << setw(8) << timesliceL0.rbegin()->getModuleID() <<
' ' << setw(8) << timesliceL0.rbegin()->size() << endl);
390 DEBUG(
"L1 " << setw(8) << timesliceL1.rbegin()->getModuleID() <<
' ' << setw(8) << timesliceL1.rbegin()->size() << endl);
391 DEBUG(
"L2 " << setw(8) << timesliceL2.rbegin()->getModuleID() <<
' ' << setw(8) << timesliceL2.rbegin()->size() << endl);
392 DEBUG(
"SN " << setw(8) << timesliceSN.rbegin()->getModuleID() <<
' ' << setw(8) << timesliceSN.rbegin()->size() << endl);
399 JTriggerInput trigger_input(timesliceL2);
400 JTriggerOutput trigger_output;
402 trigger3DMuon (trigger_input, back_inserter(trigger_output));
403 trigger3DShower(trigger_input, back_inserter(trigger_output));
404 triggerMXShower(trigger_input, timesliceL0, back_inserter(trigger_output));
406 trigger_output.merge(JEventOverlap(parameters.
TMaxEvent_ns));
408 for (JTriggerOutput::const_iterator to = trigger_output.begin(); to != trigger_output.end(); ++to) {
410 for (
int i = 0; i != h1.GetNbinsX(); ++i) {
411 if (to->hasTriggerBit(i)) {
419 << to->getFrameIndex() <<
' '
421 << timeRange << endl);
423 if (timeRange.overlap(eventTime)) {
425 JTriggeredEvent tev(*to,
433 tev.setCounter(trigger_counter);
442 if (!triggeredEventsOnly || trigger) {
449 outputFile.put(JTimesliceL1<JDAQTimesliceL1>(timesliceL1, timesliceRouter, moduleRouter, parameters.
TMaxLocal_ns));
453 outputFile.put(JTimesliceL1<JDAQTimesliceL2>(timesliceL2, timesliceRouter, moduleRouter, parameters.
L2.
TMaxLocal_ns));
457 outputFile.put(JTimesliceL1<JDAQTimesliceSN>(timesliceSN, timesliceRouter, moduleRouter, parameters.
SN.
TMaxLocal_ns));
466 if (!triggeredEventsOnly || trigger) {
477 JMultipleFileScanner<JMetaTypes_t> io(inputFile);