57 using namespace KM3NETDAQ;
63 double numberOfBlocks;
69 JParser<> zap(
"Example program to plot event rates using JASTRONOMY::JStarTrek.");
73 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
81 catch(
const exception &error) {
82 FATAL(error.what() << endl);
87 const double radius = 1.0;
88 const double omega = 2*
PI * (1.0 - cos(radius*
PI/180));
103 double Wall = numberOfBlocks;
108 if (!buffer.is_valid(&JHead::genvol) || buffer.genvol.numberOfEvents == 0) {
109 FATAL(
"No generated events." << endl);
112 Wall /= buffer.genvol.numberOfEvents;
114 STATUS(
"Generation: " << buffer.genvol.numberOfEvents << endl);
116 if (buffer.is_valid(&JHead::cut_nu)) {
117 STATUS(
"Solid angle: " << 2 * (buffer.cut_nu.cosTmax - buffer.cut_nu.cosTmin) <<
"pi" << endl);
130 TH1D hc(
"ct [live time]", NULL, 1000, -1.0, +1.0);
132 for (
int i = 1; i != hc.GetNbinsX(); ++i) {
134 const double x = hc.GetBinCenter(i);
135 const double y = startrek(x);
137 hc.SetBinContent(i, y);
140 TH1D h0(
"h0 [background]", NULL, 20, +1.0, +6.0);
141 TH1D
h1(
"h1 [signal]", NULL, 20, +1.0, +6.0);
143 TH1D*
H[] = { &h0, &
h1 };
144 double W[] = { 0.0, 0.0 };
146 for (
int i = 0; i !=
sizeof(
H)/
sizeof(
H[0]); ++i) {
151 while (inputFile.hasNext()) {
153 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
158 const Evt*
event = ps;
164 const double E = neutrino.
E;
165 const double dz = neutrino.
dir.
z;
166 const double P = 1.0;
167 const double W2 =
event->w[1];
168 const double W3 =
event->w[2];
172 const double x = log10(E);
173 const double y[] = { P * W3 * omega * startrek(dz),
174 P * W2 * rxj1713(E) * startrek(dz) };
177 for (
int i = 0; i !=
sizeof(
H)/
sizeof(
H[0]); ++i) {
178 H[i]->Fill(x, y[i] * Wall);
186 for (
int i = 0; i !=
sizeof(W)/
sizeof(W[0]); ++i) {
187 NOTICE(
"W[" << i <<
"] = " <<
FIXED(9,2) << W[i] << endl);
static const JGeographicalLocation Sicily(36, 16, 16, 06)
Auxiliary class for source tracking.
Utility class to parse command line options.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
then for HISTOGRAM in h0 h1
Auxiliary data structure for floating point format specification.
static const JSourceLocation RXJ1713(getRadians(-39,-46, 0.0), getHourAngle(17, 13, 7))
double E
Energy [GeV] (either MC truth or reconstructed)
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for defining the range of iterations of objects.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
const JLimit & getLimit() const
Get limit.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
General purpose class for multiple pointers.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
then usage $script[input file[working directory[option]]] nWhere option can be E
#define DEBUG(A)
Message macros.