Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JAcousticsMonitor.cc
Go to the documentation of this file.
1#include <iostream>
2#include <iomanip>
3#include <limits>
4#include <map>
5#include <algorithm>
6
7#include "TROOT.h"
8#include "TFile.h"
9#include "TH1D.h"
10#include "TH2D.h"
11#include "TProfile2D.h"
12
14#include "JTools/JHashMap.hh"
15#include "JTools/JRange.hh"
16#include "JTools/JQuantile.hh"
17
20#include "JDetector/JTripod.hh"
22
25
26#include "JROOT/JRootToolkit.hh"
27#include "JROOT/JManager.hh"
28
29#include "JLang/JComparator.hh"
30
32#include "JAcoustics/JEvent.hh"
34
35#include "Jeep/JContainer.hh"
36#include "Jeep/JPrint.hh"
37#include "Jeep/JParser.hh"
38#include "Jeep/JMessage.hh"
39
40namespace {
41
44
45 /**
46 * Auxilisry data structure for module location and position.
47 */
48 struct JModule_t :
49 public JLocation,
50 public JPosition3D
51 {
52 /**
53 * Default constructor.
54 */
55 JModule_t()
56 {}
57
58 /**
59 * Constructor.
60 *
61 * \param location location
62 * \param position position
63 */
64 JModule_t(const JLocation& location,
65 const JPosition3D& position) :
66 JLocation (location),
67 JPosition3D(position)
68 {}
69 };
70}
71
72
73/**
74 * \file
75 *
76 * Example program to monitor acoustic events.
77 * \author mdejong
78 */
79int main(int argc, char **argv)
80{
81 using namespace std;
82 using namespace JPP;
83
86
87 JMultipleFileScanner<> inputFile;
88 JLimit_t& numberOfEvents = inputFile.getLimit();
89 string outputFile;
90 string detectorFile;
91 tripods_container tripods; // tripods
92 transmitters_container transmitters; // transmitters
93 double Q;
94 int debug;
95
96 try {
97
98 JParser<> zap("Example program to monitor acoustic events.");
99
100 zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
101 zap['n'] = make_field(numberOfEvents) = JLimit::max();
102 zap['o'] = make_field(outputFile) = "monitor.root";
103 zap['a'] = make_field(detectorFile);
104 zap['T'] = make_field(tripods, "tripod data") = JPARSER::initialised();
105 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
106 zap['Q'] = make_field(Q) = 0.9;
107 zap['d'] = make_field(debug) = 2;
108
109 zap(argc, argv);
110 }
111 catch(const exception &error) {
112 FATAL(error.what() << endl);
113 }
114
116
117 try {
118 load(detectorFile, detector);
119 }
120 catch(const JException& error) {
121 FATAL(error);
122 }
123
125
126 JHashMap<int, JModule_t> receivers;
128
129 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
130 receivers[i->getID()] = JModule_t(i->getLocation(), i->getPosition());
131 }
132
133 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
134 emitters[i->getID()] = (i->getUTMPosition() - detector.getUTMPosition());
135 }
136
137 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
138 try {
139 emitters[i->getID()] = (i->getPosition() + detector.getModule(i->getLocation()).getPosition());
140 }
141 catch(const exception&) {
142 continue; // if no module available, discard transmitter
143 }
144 }
145
147 const JRange <int> floor (make_array(detector.begin(), detector.end(), &JModule::getFloor));
148
149 JManager<int, TH1D> HA(new TH1D("H[%].size", NULL, detector.size() + 1, -0.5, detector.size() + 0.5));
150 JManager<int, TH1D> HB(new TH1D("H[%].overlays", NULL, 201, -0.5, 200.5));
151 JManager<int, TH1D> HC(new TH1D("H[%].time", NULL, 800, 0.0, 4.0));
152 JManager<int, TH1D> HD(new TH1D("H[%].rms", NULL, 500, 0.0, 1.0e-1));
153 JManager<int, TH1D> HE(new TH1D("H[%].quantile", NULL, 500, 0.0, 1.0e-1));
154 JManager<int, TH1D> HQ(new TH1D("H[%].quality", NULL, 200, 0.0, 8.0));
155 JManager<int, TH2D> HR(new TH2D("H[%].QR", NULL, 40, 1.5, 3.5, 40, 3.0, 7.0));
156 JManager<int, TH1D> H1(new TH1D("H[%].multiplicity", NULL, 101, -0.5, 100.5));
157 JManager<int, TH2D> H2(new TH2D("H[%].event", NULL,
158 string.size(), -0.5, string.size() - 0.5,
159 floor.getUpperLimit() + 1, - 0.5, floor.getUpperLimit() + 0.5));
160 JManager<int, TProfile2D> H4(new TProfile2D("H[%].log10(Q)", NULL,
161 string.size(), -0.5, string.size() - 0.5,
162 floor.getUpperLimit() + 1, - 0.5, floor.getUpperLimit() + 0.5));
163
164 for (Int_t i = 1; i <= H2->GetXaxis()->GetNbins(); ++i) {
165 H2->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
166 }
167
168 for (Int_t i = 1; i <= H2->GetYaxis()->GetNbins(); ++i) {
169 H2->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i-1));
170 }
171
172 for (Int_t i = 1; i <= H4->GetXaxis()->GetNbins(); ++i) {
173 H4->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
174 }
175
176 for (Int_t i = 1; i <= H4->GetYaxis()->GetNbins(); ++i) {
177 H4->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i-1));
178 }
179
180 JManager<int, TH2D> H3((TH2D*) H2->Clone("H[%].doubles"));
181
182 for (JHashMap<int, JPosition3D>::const_iterator i = emitters.begin(); i != emitters.end(); ++i) {
183 H2[i->first];
184 H3[i->first];
185 }
186
188
189 counter_type counter = 0;
190
191 for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
192
193 JTreeScanner_t in(*i);
194
196
197 for (JTreeScanner_t::iterator event = in.begin(); event != in.end() && counter != inputFile.getLimit(); ++event, ++counter) {
198
199 if (counter%1000 == 0) {
200 STATUS("event " << setw(8) << counter << '\r'); DEBUG(endl);
201 }
202
203 const JEvent evt = *event;
204
205 JQuantile Q1("", true);
206
207 for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
208 Q1.put(i->getToE());
209 }
210
211 HA[evt.getID()]->Fill((double) evt.size());
212 HB[evt.getID()]->Fill((double) evt.getOverlays());
213
214 if (toe.has(evt.getID())) {
215 HC[evt.getID()]->Fill(log10(evt.begin()->getToE() - toe[evt.getID()]));
216 }
217
218 HD[evt.getID()]->Fill(Q1.getSTDev());
219 HE[evt.getID()]->Fill(Q1.getQuantile(Q, JQuantile::symmetric_t));
220
221 toe[evt.getID()] = evt.begin()->getToE();
222
223 TH1D* hq = HQ[evt.getID()];
224 TH2D* hr = HR[evt.getID()];
225 TH1D* h1 = H1[evt.getID()];
226 TH2D* h2 = H2[evt.getID()];
227 TH2D* h3 = H3[evt.getID()];
228 TProfile2D* h4 = H4[evt.getID()];
229
230 for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
231 hq->Fill(log10(i->getQ()));
232 }
233
234 if (emitters.has(evt.getID())) {
235
236 const JPosition3D& p0 = emitters[evt.getID()];
237
238 for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
239 if (receivers.has(i->getID())) {
240 HR->Fill(log10(p0.getDistance(receivers[i->getID()])), log10(i->getQ()));
241 hr->Fill(log10(p0.getDistance(receivers[i->getID()])), log10(i->getQ()));
242 }
243 }
244 }
245
246 map<int, set<double> > buffer;
247
248 for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
249 buffer[i->getID()].insert(i->getQ());
250 }
251
252 for (map<int, set<double> >::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
253
254 h1->Fill((double) i->second.size());
255
256 if (receivers.has(i->first)) {
257
258 const JLocation& location = receivers[i->first];
259
260 const double x = string.getIndex(location.getString());
261 const double y = location.getFloor();
262
263 h2->Fill(x, y, 1.0);
264
265 if (i->second.size() >= 2u) {
266 h3->Fill(x, y, 1.0);
267 }
268
269 h4->Fill(x, y, log10(*(i->second.rbegin())));
270 }
271 }
272 }
273 }
274 STATUS(endl);
275
276
277 TFile out(outputFile.c_str(), "recreate");
278
279 out << HA << HB << HC << HD << HE << HQ << HR << *HR << H1 << H2 << H3 << H4;
280
281 out.Write();
282 out.Close();
283}
int main(int argc, char **argv)
Acoustic event.
ROOT TTree parameter settings.
Container I/O.
string outputFile
Data structure for detector geometry and calibration.
Acoustic emitter.
General purpose class for a hash collection of unique elements.
General purpose class for hash map of unique elements.
Dynamic ROOT object management.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
Auxiliary class to define a range between two values.
Data structure for transmitter.
Data structure for tripod.
Detector data structure.
Definition JDetector.hh:96
Logical location of module.
Definition JLocation.hh:40
const JLocation & getLocation() const
Get location.
Definition JLocation.hh:70
int getFloor() const
Get floor number.
Definition JLocation.hh:146
int getString() const
Get string number.
Definition JLocation.hh:135
Data structure for position in three dimensions.
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition JVector3D.hh:270
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition JManager.hh:304
Base class for JTreeScanner.
Template definition for direct access of elements in ROOT TChain.
General purpose class for hash collection of unique elements.
bool has(const T &value) const
Test whether given value is present.
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:79
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:81
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition JVectorize.hh:54
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Long64_t counter_type
Type definition for counter.
Detector file.
Definition JHead.hh:227
int getID() const
Get emitter identifier.
int getOverlays() const
Get number of overlayed events.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
General purpose class for hash map of unique keys.
Definition JHashMap.hh:75
container_type::const_iterator const_iterator
Definition JHashMap.hh:86
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
double getQuantile(const double Q, const Quantile_t option=forward_t) const
Get quantile.
Definition JQuantile.hh:349
@ symmetric_t
symmatric
Definition JQuantile.hh:337
double getSTDev() const
Get standard deviation.
Definition JQuantile.hh:281
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133