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));
53int main(
int argc,
char **argv)
60 JLimit_t& numberOfEvents = inputFile.getLimit();
63 double numberOfBlocks;
69 JParser<> zap(
"Example program to plot event rates using JASTRONOMY::JStarTrek.");
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;
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) {
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);
Interface methods for slalib and auxiliary classes and methods for astronomy.
General purpose messaging.
#define DEBUG(A)
Message macros.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
I/O formatting auxiliaries.
int main(int argc, char **argv)
ROOT TTree parameter settings of various packages.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
Utility class to parse command line options.
counter_type getCounter() const
Get counter.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
static const JGeographicalLocation Sicily(36, 16, 16, 06)
static const JSourceLocation RXJ1713(getRadians(-39, -46, 0.0), getHourAngle(17, 13, 7))
T pow(const T &x, const double y)
Power .
static const double H
Planck constant [eV s].
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
KM3NeT DAQ data structures and auxiliaries.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Auxiliary data structure for floating point format specification.
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
JRange_t cosT
Cosine zenith angle range
double numberOfEvents
Number of events.
Auxiliary class for source tracking.
General purpose class for multiple pointers.
Auxiliary class for defining the range of iterations of objects.
static counter_type max()
Get maximum counter value.
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
virtual bool hasNext() override
Check availability of next element.
virtual const multi_pointer_type & next() override
Get next element.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
double E
Energy [GeV] (either MC truth or reconstructed)