Go to the documentation of this file.
10 #include "evt/Head.hh"
36 inline double rxj1713(
const double E)
38 static const double k = 16.80e-15;
39 static const double a = 1.72;
40 static const double e = 2.10;
42 return 1e4 * k * pow(E*1e-3, -a) * exp(-sqrt(E*1e-3/e));
53 int main(
int argc,
char **argv)
59 JTriggeredFileScanner<> inputFile;
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));
98 catch(
const JException& error) {
103 double Wall = numberOfBlocks;
109 FATAL(
"No generated events." << 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);
155 JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
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);
Auxiliary data structure for floating point format specification.
static const JGeographicalLocation Sicily(36, 16, 16, 06)
bool is_valid(const T &value)
Check validity of given value.
JLimit JLimit_t
Type definition of limit.
Utility class to parse command line options.
double numberOfEvents
Number of events.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
const JLimit & getLimit() const
Get limit.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
double cosTmin
Minimal cosine zenith angle.
double cosTmax
Maximal cosine zenith angle.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
#define DEBUG(A)
Message macros.
KM3NeT DAQ data structures and auxiliaries.
static const JSourceLocation RXJ1713(getRadians(-39, -46, 0.0), getHourAngle(17, 13, 7))
int main(int argc, char **argv)