Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
10 #include "evt/Head.hh"
11 #include "evt/Evt.hh"
12 #include "JAAnet/JHead.hh"
13 #include "JAAnet/JHeadToolkit.hh"
14 #include "JAAnet/JAAnetToolkit.hh"
15 #include "JDAQ/JDAQEvent.hh"
16 
17 #include "JAstronomy/JAstronomy.hh"
19 
22 #include "JSupport/JSupport.hh"
23 
24 #include "Jeep/JPrint.hh"
25 #include "Jeep/JParser.hh"
26 #include "Jeep/JMessage.hh"
27 
28 
29 namespace {
30 
31  /**
32  * Neutrino flux for RXJ1713.
33  *
34  * \param E neutrino energy [GeV]
35  * \return flux [GeV^-1 s^-1 m^-2]
36  */
37  inline double rxj1713(const double E)
38  {
39  static const double k = 16.80e-15; // [GeV^-1 s^-1 cm^-2]
40  static const double a = 1.72;
41  static const double e = 2.10; // [TeV]
42 
43  return 1e4 * k * pow(E*1e-3, -a) * exp(-sqrt(E*1e-3/e));
44  }
45 }
46 
47 
48 /**
49  * \file
50  *
51  * Example program to plot event rates using JASTRONOMY::JStarTrek.
52  * \author mdejong
53  */
54 int main(int argc, char **argv)
55 {
56  using namespace std;
57  using namespace JPP;
58  using namespace KM3NETDAQ;
59 
60  JTriggeredFileScanner<> inputFile;
61  JLimit_t& numberOfEvents = inputFile.getLimit();
62  string outputFile;
63  double ct;
64  double numberOfBlocks;
65  bool join;
66  int debug;
67 
68  try {
69 
70  JParser<> zap("Example program to plot event rates using JASTRONOMY::JStarTrek.");
71 
72  zap['f'] = make_field(inputFile);
73  zap['o'] = make_field(outputFile) = "volume.root";
74  zap['n'] = make_field(numberOfEvents) = JLimit::max();
75  zap['Z'] = make_field(ct) = -0.25;
76  zap['N'] = make_field(numberOfBlocks) = 1.0;
77  zap['J'] = make_field(join); // join neutrino and anti-neutrino
78  zap['d'] = make_field(debug) = 1;
79 
80  zap(argc, argv);
81  }
82  catch(const exception &error) {
83  FATAL(error.what() << endl);
84  }
85  using namespace std;
86 
87 
88  const double radius = 1.0; // size of source [deg]
89  const double omega = 2*PI * (1.0 - cos(radius*PI/180)); // solid angle of source
90 
91  const JStarTrek startrek(Sicily, RXJ1713);
92 
93 
94  Head head;
95 
96  try {
97  head = getHeader(inputFile);
98  }
99  catch(const JException& error) {
100  FATAL(error);
101  }
102 
103 
104  double Wall = numberOfBlocks;
105 
106  {
107  JHead buffer(head);
108 
109  if (!is_valid(buffer.genvol) || buffer.genvol.numberOfEvents == 0) {
110  FATAL("No generated events." << endl);
111  }
112 
113  Wall /= buffer.genvol.numberOfEvents;
114 
115  STATUS("Generation: " << buffer.genvol.numberOfEvents << endl);
116 
117  if (is_valid(buffer.cut_nu)) {
118  STATUS("Solid angle: " << 2 * (buffer.cut_nu.cosTmax - buffer.cut_nu.cosTmin) << "pi" << endl);
119  }
120  }
121 
122  if (join) {
123  Wall *= 0.5; // join neutrino and anti-neutrino
124  }
125 
126 
127  // output
128 
129  TFile out(outputFile.c_str(), "recreate");
130 
131  TH1D hc("ct [live time]", NULL, 1000, -1.0, +1.0);
132 
133  for (int i = 1; i != hc.GetNbinsX(); ++i) {
134 
135  const double x = hc.GetBinCenter(i);
136  const double y = startrek(x);
137 
138  hc.SetBinContent(i, y);
139  }
140 
141  TH1D h0("h0 [background]", NULL, 20, +1.0, +6.0);
142  TH1D h1("h1 [signal]", NULL, 20, +1.0, +6.0);
143 
144  TH1D* H[] = { &h0, &h1 };
145  double W[] = { 0.0, 0.0 };
146 
147  for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
148  H[i]->Sumw2();
149  }
150 
151 
152  while (inputFile.hasNext()) {
153 
154  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
155 
156  JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
157 
158  //const JDAQEvent* tev = ps;
159  const Evt* event = ps;
160 
161  if (has_neutrino(*event)) {
162 
163  const Trk& neutrino = get_neutrino(*event);
164 
165  const double E = neutrino.E; // GeV
166  const double dz = neutrino.dir.z; // cos(zenith angle)
167  const double P = 1.0; // P(Earth)?
168  const double W2 = event->w[1]; // generation weight GeV x m^2 x sr x (sec/year)
169  const double W3 = event->w[2]; // number of atmospheric neutrino events / year
170 
171  // number of events / year
172 
173  const double x = log10(E);
174  const double y[] = { P * W3 * omega * startrek(dz),
175  P * W2 * rxj1713(E) * startrek(dz) };
176 
177  if (dz >= ct) {
178  for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
179  H[i]->Fill(x, y[i] * Wall);
180  W[i] += y[i] * Wall;
181  }
182  }
183  }
184  }
185  STATUS(endl);
186 
187  for (int i = 0; i != sizeof(W)/sizeof(W[0]); ++i) {
188  NOTICE("W[" << i << "] = " << FIXED(9,2) << W[i] << endl);
189  }
190 
191  out.Write();
192  out.Close();
193 }
194 
195 
static const JGeographicalLocation Sicily(36, 16, 16, 06)
static const double H
Planck constant [eV s].
Definition: JConstants.hh:25
Utility class to parse command line options.
Definition: JParser.hh:1410
JAANET::genvol genvol
Definition: JHead.hh:864
JRange< T, JComparator_t > join(const JRange< T, JComparator_t > &first, const JRange< T, JComparator_t > &second)
Join ranges.
Definition: JRange.hh:573
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
#define STATUS(A)
Definition: JMessage.hh:61
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
double numberOfEvents
Number of events.
Definition: JHead.hh:423
static const double PI
Constants.
Definition: JConstants.hh:20
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:461
string outputFile
static const JSourceLocation RXJ1713(getRadians(-39,-46, 0.0), getHourAngle(17, 13, 7))
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:62
int debug
debug level
Definition: JSirene.cc:59
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:727
#define FATAL(A)
Definition: JMessage.hh:65
Interface methods for slalib and auxiliary classes and methods for astronomy.
Utility class to parse command line options.
ROOT TTree parameter settings.
double cosTmax
Maximal cosine zenith angle.
Definition: JHead.hh:209
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
JAANET::cut_nu cut_nu
Definition: JHead.hh:861
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
double cosTmin
Minimal cosine zenith angle.
Definition: JHead.hh:208
bool is_valid(const T &value)
Check validity of given value.
Definition: JHead.hh:711
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
int main(int argc, char *argv[])
Definition: Main.cpp:15