Jpp  19.1.0-rc.1
the software that should make you happy
JRate.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 
5 #include "TROOT.h"
6 #include "TFile.h"
7 #include "TH1D.h"
8 #include "TProfile.h"
9 
12 #include "JAAnet/JHead.hh"
13 #include "JAAnet/JAAnetToolkit.hh"
14 #include "JDAQ/JDAQEventIO.hh"
15 
16 #include "JAstronomy/JAstronomy.hh"
18 
21 #include "JSupport/JSupport.hh"
22 
23 #include "Jeep/JPrint.hh"
24 #include "Jeep/JParser.hh"
25 #include "Jeep/JMessage.hh"
26 
27 
28 namespace {
29 
30  /**
31  * Neutrino flux for RXJ1713.
32  *
33  * \param E neutrino energy [GeV]
34  * \return flux [GeV^-1 s^-1 m^-2]
35  */
36  inline double rxj1713(const double E)
37  {
38  static const double k = 16.80e-15; // [GeV^-1 s^-1 cm^-2]
39  static const double a = 1.72;
40  static const double e = 2.10; // [TeV]
41 
42  return 1e4 * k * pow(E*1e-3, -a) * exp(-sqrt(E*1e-3/e));
43  }
44 }
45 
46 
47 /**
48  * \file
49  *
50  * Example program to plot event rates using JASTRONOMY::JStarTrek.
51  * \author mdejong
52  */
53 int main(int argc, char **argv)
54 {
55  using namespace std;
56  using namespace JPP;
57  using namespace KM3NETDAQ;
58 
59  JTriggeredFileScanner<> inputFile;
60  JLimit_t& numberOfEvents = inputFile.getLimit();
61  string outputFile;
62  double ct;
63  double numberOfBlocks;
64  bool join;
65  int debug;
66 
67  try {
68 
69  JParser<> zap("Example program to plot event rates using JASTRONOMY::JStarTrek.");
70 
71  zap['f'] = make_field(inputFile);
72  zap['o'] = make_field(outputFile) = "volume.root";
73  zap['n'] = make_field(numberOfEvents) = JLimit::max();
74  zap['Z'] = make_field(ct) = -0.25;
75  zap['N'] = make_field(numberOfBlocks) = 1.0;
76  zap['J'] = make_field(join); // join neutrino and anti-neutrino
77  zap['d'] = make_field(debug) = 1;
78 
79  zap(argc, argv);
80  }
81  catch(const exception &error) {
82  FATAL(error.what() << endl);
83  }
84  using namespace std;
85 
86 
87  const double radius = 1.0; // size of source [deg]
88  const double omega = 2*PI * (1.0 - cos(radius*PI/180)); // solid angle of source
89 
90  const JStarTrek startrek(Sicily, RXJ1713);
91 
92 
93  Head head;
94 
95  try {
96  head = getHeader(inputFile);
97  }
98  catch(const JException& error) {
99  FATAL(error);
100  }
101 
102 
103  double Wall = numberOfBlocks;
104 
105  {
106  JHead buffer(head);
107 
108  if (!buffer.is_valid(&JHead::genvol) || buffer.genvol.numberOfEvents == 0) {
109  FATAL("No generated events." << endl);
110  }
111 
112  Wall /= buffer.genvol.numberOfEvents;
113 
114  STATUS("Generation: " << buffer.genvol.numberOfEvents << endl);
115 
116  if (buffer.is_valid(&JHead::cut_nu)) {
117  STATUS("Solid angle: " << 2 * buffer.cut_nu.cosT.getLength() << "pi" << endl);
118  }
119  }
120 
121  if (join) {
122  Wall *= 0.5; // join neutrino and anti-neutrino
123  }
124 
125 
126  // output
127 
128  TFile out(outputFile.c_str(), "recreate");
129 
130  TH1D hc("ct [live time]", NULL, 1000, -1.0, +1.0);
131 
132  for (int i = 1; i != hc.GetNbinsX(); ++i) {
133 
134  const double x = hc.GetBinCenter(i);
135  const double y = startrek(x);
136 
137  hc.SetBinContent(i, y);
138  }
139 
140  TH1D h0("h0 [background]", NULL, 20, +1.0, +6.0);
141  TH1D h1("h1 [signal]", NULL, 20, +1.0, +6.0);
142 
143  TH1D* H[] = { &h0, &h1 };
144  double W[] = { 0.0, 0.0 };
145 
146  for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
147  H[i]->Sumw2();
148  }
149 
150 
151  while (inputFile.hasNext()) {
152 
153  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
154 
156 
157  //const JDAQEvent* tev = ps;
158  const Evt* event = ps;
159 
160  if (has_neutrino(*event)) {
161 
162  const Trk& neutrino = get_neutrino(*event);
163 
164  const double E = neutrino.E; // GeV
165  const double dz = neutrino.dir.z; // cos(zenith angle)
166  const double P = 1.0; // P(Earth)?
167  const double W2 = event->w[1]; // generation weight GeV x m^2 x sr x (sec/year)
168  const double W3 = event->w[2]; // number of atmospheric neutrino events / year
169 
170  // number of events / year
171 
172  const double x = log10(E);
173  const double y[] = { P * W3 * omega * startrek(dz),
174  P * W2 * rxj1713(E) * startrek(dz) };
175 
176  if (dz >= ct) {
177  for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
178  H[i]->Fill(x, y[i] * Wall);
179  W[i] += y[i] * Wall;
180  }
181  }
182  }
183  }
184  STATUS(endl);
185 
186  for (int i = 0; i != sizeof(W)/sizeof(W[0]); ++i) {
187  NOTICE("W[" << i << "] = " << FIXED(9,2) << W[i] << endl);
188  }
189 
190  out.Write();
191  out.Close();
192 }
193 
194 
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Interface methods for slalib and auxiliary classes and methods for astronomy.
string outputFile
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define STATUS(A)
Definition: JMessage.hh:63
#define NOTICE(A)
Definition: JMessage.hh:64
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
I/O formatting auxiliaries.
int main(int argc, char **argv)
Definition: JRate.cc:53
ROOT TTree parameter settings of various packages.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
Monte Carlo run header.
Definition: JHead.hh:1236
JAANET::cut_nu cut_nu
Definition: JHead.hh:1596
JAANET::genvol genvol
Definition: JHead.hh:1600
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
Definition: JHead.hh:1319
General exception.
Definition: JException.hh:24
Utility class to parse command line options.
Definition: JParser.hh:1714
counter_type getCounter() const
Get counter.
T getLength() const
Get length (difference between upper and lower limit).
Definition: JRange.hh:289
const double a
Definition: JQuadrature.cc:42
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
bool has_neutrino(const Evt &evt)
Test whether given event has an 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 .
Definition: JMath.hh:97
static const double PI
Mathematical constants.
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.
JRange< T, JComparator_t > join(const JRange< T, JComparator_t > &first, const JRange< T, JComparator_t > &second)
Join ranges.
Definition: JRange.hh:659
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
Definition: JSTDTypes.hh:14
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:21
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
JRange_t cosT
Cosine zenith angle range
Definition: JHead.hh:420
double numberOfEvents
Number of events.
Definition: JHead.hh:721
Auxiliary class for source tracking.
General purpose class for multiple pointers.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:45
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.
Definition: Trk.hh:15
Vec dir
track direction
Definition: Trk.hh:18
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
double z
Definition: Vec.hh:14