Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
JFeedback.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <set>
5
6#include "TROOT.h"
7#include "TFile.h"
8
9#include "km3net-dataformat/offline/Head.hh"
10#include "km3net-dataformat/offline/Evt.hh"
11#include "km3net-dataformat/offline/Evt.hh"
12#include "km3net-dataformat/definitions/weightlist.hh"
13#include "km3net-dataformat/online/JDAQClock.hh"
14
15#include "flux/Flux.hh"
16
17#include "JAAnet/JHead.hh"
18#include "JLang/JManip.hh"
19
20#include "JSupport/JSingleFileScanner.hh"
21#include "JSupport/JSupport.hh"
22#include "JSupport/JMonteCarloFileSupportkit.hh"
23
24#include "Jeep/JParser.hh"
25#include "Jeep/JMessage.hh"
26
27namespace {
28
29 using JLANG::JObjectIterator;
30
31 const Flux_Atmospheric conventional;
32
33
34 /**
35 * Auxiliary data structure for run statistics.
36 */
37 struct stats_t {
38
39 std::string name;
40 double events = 0; //!< number of events
41 double range_s = 0.0; //!< time difference between first and last event
42 double livetime_s = 0.0; //!< livetime of events
43 double daq_s = 0.0; //!< livetime of data taking run
44 double gap_s = 0.0; //!< maximal time between two consecutive events
45 long int tin_s = 0.0; //!< time interval
46
47 friend inline std::ostream& operator<<(std::ostream& out, const stats_t& object)
48 {
49 using namespace std;
50
51 out << object.name << endl;
52 out << "events: " << FIXED(9,1) << object.events << endl;
53 out << "range: " << FIXED(7,1) << object.range_s << " [s]" << endl;
54 out << "livetime: " << FIXED(7,1) << object.livetime_s << " [s]" << endl;
55 out << "DAQ: " << FIXED(7,1) << object.daq_s << " [s]" << endl;
56 out << "gap: " << FIXED(7,1) << object.gap_s << " [s]" << endl;
57 out << "interval: " << setw(7) << object.tin_s << " [s]" << endl;
58
59 return out;
60 }
61 };
62
63
64 /**
65 * Get run statistics for real data.
66 *
67 * \param head header
68 * \param in event input
69 * \return statistics
70 */
71 inline stats_t getStatsFromData(const JHead& head, JObjectIterator<Evt>& in)
72 {
73 using namespace std;
74 using namespace JPP;
75 using namespace KM3NETDAQ;
76
77 stats_t object;
78
79 set<int> buffer;
80 double events = 0;
81
82 while (in.hasNext()) {
83
84 const Evt* evt = in.next();
85
86 buffer.insert(evt->frame_index);
87
88 events += 1;
89 }
90
91 int gap = -1;
92
93 if (buffer.size() >= 2) {
94 for (set<int>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
95 if (*q - *p > gap) {
96 gap = *q - *p;
97 }
98 }
99 }
100
101 object.name = "data";
102 object.events = events;
103 object.range_s = (*buffer.rbegin() - *buffer.begin()) * getFrameTime() * 1.0e-9;
104 object.livetime_s = buffer.size() * getFrameTime() * 1.0e-9;
105 object.daq_s = head.DAQ.livetime_s;
106 object.gap_s = gap * getFrameTime() * 1.0e-9;
107 object.tin_s = 0;
108
109 return object;
110 }
111
112
113 /**
114 * Get run statistics for MUPAGE.
115 *
116 * \param head header
117 * \param in event input
118 * \return statistics
119 */
120 inline stats_t getStatsFromMUPAGE(const JHead& head, JObjectIterator<Evt>& in)
121 {
122 using namespace std;
123 using namespace JPP;
124
125 stats_t object;
126
127 set<double> buffer;
128 double events = 0;
129
130 while (in.hasNext()) {
131
132 const Evt* evt = in.next();
133
134 buffer.insert(evt->mc_event_time.AsDouble());
135
136 events += 1;
137 }
138
139 double gap = -1.0;
140
141 if (buffer.size() >= 2) {
142 for (set<double>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
143 if (*q - *p > gap) {
144 gap = *q - *p;
145 }
146 }
147 }
148
149 object.name = "MUPAGE";
150 object.events = events * head.DAQ.livetime_s / head.livetime.numberOfSeconds;
151 object.range_s = (*buffer.rbegin() - *buffer.begin());
152 object.livetime_s = head.livetime.numberOfSeconds;
153 object.daq_s = head.DAQ.livetime_s;
154 object.gap_s = gap;
155 object.tin_s = head.time_interval.t2 - head.time_interval.t1;
156
157 return object;
158 }
159
160
161 /**
162 * Get run statistics for gSeaGen.
163 *
164 * \param head header
165 * \param in event input
166 * \return statistics
167 */
168 inline stats_t getStatsFromGSEAGEN(const JHead& head, JObjectIterator<Evt>& in)
169 {
170 using namespace std;
171 using namespace JPP;
172
173 stats_t object;
174
175 set<double> buffer;
176 double events = 0;
177
178 while (in.hasNext()) {
179
180 const Evt* evt = in.next();
181
182 buffer.insert(evt->mc_event_time.AsDouble());
183
184 const Trk trk = get_neutrino(*evt);
185
186 events += evt->w[WEIGHTLIST_DIFFERENTIAL_EVENT_RATE] * evt->w[WEIGHTLIST_NORMALISATION] * conventional.dNdEdOmega(trk.type, log10(trk.E), trk.dir.z);
187 }
188
189 double gap = -1.0;
190
191 if (buffer.size() >= 2) {
192 for (set<double>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
193 if (*q - *p > gap) {
194 gap = *q - *p;
195 }
196 }
197 }
198
199 object.name = "gSeaGen";
200 object.events = events * head.DAQ.livetime_s;
201 object.range_s = (*buffer.rbegin() - *buffer.begin());
202 object.livetime_s = 0.0;
203 object.daq_s = head.DAQ.livetime_s;
204 object.gap_s = gap;
205 object.tin_s = head.time_interval.t2 - head.time_interval.t1;
206
207 return object;
208 }
209
210
211 /**
212 * Get run statistics for pure-noise.
213 *
214 * \param head header
215 * \param in event input
216 * \return statistics
217 */
218 inline stats_t getStatsFromK40(const JHead& head, JObjectIterator<Evt>& in)
219 {
220 using namespace std;
221 using namespace JPP;
222
223 stats_t object;
224
225 set<int> buffer;
226 double events = 0;
227
228 while (in.hasNext()) {
229
230 const Evt* evt = in.next();
231
232 buffer.insert(evt->frame_index);
233
234 events += 1;
235 }
236
237 int gap = -1;
238
239 if (buffer.size() >= 2) {
240 for (set<int>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {
241 if (*q - *p > gap) {
242 gap = *q - *p;
243 }
244 }
245 }
246
247 object.name = "noise";
248 object.events = events * head.DAQ.livetime_s / head.K40.livetime_s;
249 object.range_s = (*buffer.rbegin() - *buffer.begin()) * getFrameTime() * 1.0e-9;
250 object.livetime_s = head.K40.livetime_s;
251 object.daq_s = head.DAQ.livetime_s;
252 object.gap_s = gap * getFrameTime() * 1.0e-9;
253 object.tin_s = 0;
254
255 return object;
256 }
257
258
259 /**
260 * Get run statistics.
261 *
262 * \param head header
263 * \param in event input
264 * \return statistics
265 */
266 inline stats_t getStats(const JHead& head, JObjectIterator<Evt>& in)
267 {
268 using namespace std;
269 using namespace JPP;
270
271 const bool mupage = is_mupage (head);
272 const bool gseagen = is_gseagen(head);
273 const bool k40 = head.is_valid(&JHead::K40);
274
275 if (mupage)
276 return getStatsFromMUPAGE (head, in);
277 else if (gseagen)
278 return getStatsFromGSEAGEN(head, in);
279 else if (k40)
280 return getStatsFromK40 (head, in);
281 else
282 return getStatsFromData (head, in);
283 }
284}
285
286
287/**
288 * \file
289 *
290 * Example program to analyse track fit results from Evt formatted data.
291 * \author mdejong
292 */
293int main(int argc, char **argv)
294{
295 using namespace std;
296 using namespace JPP;
297
298 JSingleFileScanner<Evt> inputFile;
299 JLimit_t& numberOfEvents = inputFile.getLimit();
300 int debug;
301
302 try {
303
304 JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
305
306 zap['f'] = make_field(inputFile);
307 zap['n'] = make_field(numberOfEvents) = JLimit::max();
308 zap['d'] = make_field(debug) = 2;
309
310 zap(argc, argv);
311 }
312 catch(const exception& error) {
313 FATAL(error.what() << endl);
314 }
315
316
317 const JHead head = getHeader(inputFile);
318 const stats_t stats = getStats(head, inputFile);
319
320 DEBUG(inputFile << endl << head << endl << stats << endl);
321
322 cout << setw(132) << left << inputFile << right << ' '
323 << FIXED(9,1) << stats.events << ' '
324 << FIXED(7,1) << stats.range_s << ' '
325 << FIXED(7,1) << stats.livetime_s << ' '
326 << FIXED(7,1) << stats.daq_s << ' '
327 << FIXED(7,1) << stats.gap_s << ' '
328 << setw(7) << stats.tin_s << endl;
329}
main
int main(int argc, char **argv)
Definition JFeedback.cc:293
JManip.hh
I/O manipulators.
JMessage.hh
General purpose messaging.
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JSingleFileScanner.hh
Scanning of objects from a single file according a format that follows from the extension of each fil...
JSupport.hh
ROOT TTree parameter settings of various packages.
JAANET::JHead
Monte Carlo run header.
Definition JHead.hh:1236
JAANET::JHead::K40
JAANET::K40 K40
Definition JHead.hh:1626
JLANG::JObjectIterator
Interface of object iteration for a single data type.
Definition JObjectIterator.hh:45
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JTOOLS::JMultiFunction::insert
void insert(const JMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Insert multidimensional input.
Definition JMultiFunction.hh:121
JAANET::is_gseagen
bool is_gseagen(const JHead &header)
Check for generator.
Definition JHeadToolkit.hh:61
JAANET::is_mupage
bool is_mupage(const JHead &header)
Check for generator.
Definition JHeadToolkit.hh:85
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition JAAnetToolkit.hh:542
JASTRONOMY::operator<<
std::ostream & operator<<(std::ostream &out, const morphology_type &object)
Write morphology to output stream.
Definition JMorphology.hh:570
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition JMonteCarloFileSupportkit.hh:547
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
KM3NETDAQ::getFrameTime
double getFrameTime()
Get frame time duration.
Definition JDAQClock.hh:162
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
Evt::w
std::vector< double > w
MC: Weights w[0]=w1, w[1]=w2, w[2]=w3 (see e.g. Tag list or km3net-dataformat/definitions)
Definition Evt.hh:42
Evt::frame_index
int frame_index
from the raw data
Definition Evt.hh:29
Evt::mc_event_time
TTimeStamp mc_event_time
MC: true generation time (UTC) of the event, (default: 01 Jan 1970 00:00:00)
Definition Evt.hh:46
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
WEIGHTLIST_DIFFERENTIAL_EVENT_RATE
static const int WEIGHTLIST_DIFFERENTIAL_EVENT_RATE
Event rate per unit of flux (c.f. taglist document) [GeV m2 sr].
Definition weightlist.hh:15
WEIGHTLIST_NORMALISATION
static const int WEIGHTLIST_NORMALISATION
Event rate normalisation.
Definition weightlist.hh:17
Generated by doxygen 1.10.0