Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
18 #include "JDetector/JDetector.hh"
20 #include "JDetector/JTripod.hh"
22 
24 #include "JSupport/JTreeScanner.hh"
25 
26 #include "JROOT/JRootToolkit.hh"
27 #include "JROOT/JManager.hh"
28 
29 #include "JLang/JComparator.hh"
30 
31 #include "JAcoustics/JEmitter.hh"
32 #include "JAcoustics/JEvent.hh"
33 #include "JAcoustics/JSupport.hh"
34 
35 #include "Jeep/JContainer.hh"
36 #include "Jeep/JPrint.hh"
37 #include "Jeep/JParser.hh"
38 #include "Jeep/JMessage.hh"
39 
40 namespace {
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  */
79 int 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 
124  sort(detector.begin(), detector.end(), make_comparator(&JModule::getLocation));
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 
146  const JHashCollection<int> string(make_array(detector.begin(), detector.end(), &JModule::getString));
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 }
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
int getOverlays() const
Get number of overlayed events.
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
Definition: Main.cc:15
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
JContainer< std::vector< JTransmitter > > transmitters_container
Definition: JSydney.cc:80
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:281
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
General purpose class for hash map of unique elements.
#define STATUS(A)
Definition: JMessage.hh:63
ROOT TTree parameter settings.
Detector data structure.
Definition: JDetector.hh:89
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:136
Acoustic event.
then fatal Number of tripods
Definition: JFootprint.sh:45
General purpose class for a hash collection of unique elements.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:84
Long64_t counter_type
Type definition for counter.
Dynamic ROOT object management.
string outputFile
Acoustic emitter.
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
Template definition for direct access of elements in ROOT TChain.
Data structure for detector geometry and calibration.
JMODEL::JString getString(const JFit &fit)
Get model parameters of string.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:226
Data structure for transmitter.
Logical location of module.
Definition: JLocation.hh:37
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition: JContainer.hh:39
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
set_variable E_E log10(E_{fit}/E_{#mu})"
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
JPosition3D getPosition(const Vec &pos)
Get position.
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
double getQuantile(const double Q, const Quantile_t option=forward_t) const
Get quantile.
Definition: JQuantile.hh:349
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:133
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition: JManager.hh:295
JContainer< std::vector< JTripod > > tripods_container
Definition: JSydney.cc:78
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Base class for JTreeScanner.
Definition: JTreeScanner.hh:54
int getString() const
Get string number.
Definition: JLocation.hh:134
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:48
Utility class to parse command line options.
Acoustic event.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86
do set_variable DETECTOR_TXT $WORKDIR detector
int getID() const
Get emitter identifier.
double u[N+1]
Definition: JPolint.hh:865
Data structure for tripod.
bool has(const T &value) const
Test whether given value is present.
Container I/O.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62