Jpp  17.2.1-pre0
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 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 
12 #include "JTools/JHashMap.hh"
13 #include "JTools/JRange.hh"
14 #include "JTools/JQuantile.hh"
15 
16 #include "JDetector/JDetector.hh"
18 #include "JDetector/JTripod.hh"
20 
22 #include "JSupport/JTreeScanner.hh"
23 
24 #include "JROOT/JRootToolkit.hh"
25 #include "JROOT/JManager.hh"
26 
27 #include "JAcoustics/JEmitter.hh"
28 #include "JAcoustics/JEvent.hh"
29 #include "JAcoustics/JSupport.hh"
30 
31 #include "Jeep/JContainer.hh"
32 #include "Jeep/JPrint.hh"
33 #include "Jeep/JParser.hh"
34 #include "Jeep/JMessage.hh"
35 
36 namespace {
37 
40 
41  /**
42  * Auxilisry data structure for module location and position.
43  */
44  struct JModule_t :
45  public JLocation,
46  public JPosition3D
47  {
48  /**
49  * Default constructor.
50  */
51  JModule_t()
52  {}
53 
54  /**
55  * Constructor.
56  *
57  * \param location location
58  * \param position position
59  */
60  JModule_t(const JLocation& location,
61  const JPosition3D& position) :
62  JLocation (location),
63  JPosition3D(position)
64  {}
65  };
66 }
67 
68 
69 /**
70  * \file
71  *
72  * Example program to monitor acoustic events.
73  * \author mdejong
74  */
75 int main(int argc, char **argv)
76 {
77  using namespace std;
78  using namespace JPP;
79 
80  typedef JContainer< vector<JTripod> > tripods_container;
81  typedef JContainer< vector<JTransmitter> > transmitters_container;
82 
83  JMultipleFileScanner<> inputFile;
84  JLimit_t& numberOfEvents = inputFile.getLimit();
85  string outputFile;
86  string detectorFile;
87  tripods_container tripods; // tripods
88  transmitters_container transmitters; // transmitters
89  double Q;
90  int debug;
91 
92  try {
93 
94  JParser<> zap("Example program to monitor acoustic events.");
95 
96  zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
97  zap['n'] = make_field(numberOfEvents) = JLimit::max();
98  zap['o'] = make_field(outputFile) = "monitor.root";
99  zap['a'] = make_field(detectorFile);
100  zap['T'] = make_field(tripods, "tripod data") = JPARSER::initialised();
101  zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
102  zap['Q'] = make_field(Q) = 0.9;
103  zap['d'] = make_field(debug) = 2;
104 
105  zap(argc, argv);
106  }
107  catch(const exception &error) {
108  FATAL(error.what() << endl);
109  }
110 
112 
113  try {
114  load(detectorFile, detector);
115  }
116  catch(const JException& error) {
117  FATAL(error);
118  }
119 
120 
121  JHashMap<int, JModule_t> receivers;
123 
124  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
125  receivers[i->getID()] = JModule_t(i->getLocation(), i->getPosition());
126  }
127 
128  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
129  emitters[i->getID()] = (i->getUTMPosition() - detector.getUTMPosition());
130  }
131 
132  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
133  try {
134  emitters[i->getID()] = (i->getPosition() + detector.getModule(i->getLocation()).getPosition());
135  }
136  catch(const exception&) {
137  continue; // if no module available, discard transmitter
138  }
139  }
140 
141  const JHashCollection<int> string(make_array(detector.begin(), detector.end(), &JModule::getString));
142  const JRange <int> floor (make_array(detector.begin(), detector.end(), &JModule::getFloor));
143 
144  JManager<int, TH1D> HA(new TH1D("H[%].size", NULL, detector.size() + 1, -0.5, detector.size() + 0.5));
145  JManager<int, TH1D> HB(new TH1D("H[%].overlays", NULL, 101, -0.5, 100.5));
146  JManager<int, TH1D> HC(new TH1D("H[%].time", NULL, 800, 0.0, 4.0));
147  JManager<int, TH1D> HD(new TH1D("H[%].rms", NULL, 500, 0.0, 1.0e-1));
148  JManager<int, TH1D> HE(new TH1D("H[%].quantile", NULL, 500, 0.0, 1.0e-1));
149  JManager<int, TH1D> HQ(new TH1D("H[%].quality", NULL, 200, 0.0, 8.0));
150  JManager<int, TH2D> HR(new TH2D("H[%].QR", NULL, 40, 1.5, 3.5, 40, 3.0, 7.0));
151  JManager<int, TH1D> H1(new TH1D("H[%].multiplicity", NULL, 101, -0.5, 100.5));
152  JManager<int, TH2D> H2(new TH2D("H[%].event", NULL,
153  string.size(), -0.5, string.size() - 0.5,
154  floor.getUpperLimit() + 1, - 0.5, floor.getUpperLimit() + 0.5));
155 
156  for (Int_t i = 1; i <= H2->GetXaxis()->GetNbins(); ++i) {
157  H2->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
158  }
159 
160  for (Int_t i = 1; i <= H2->GetYaxis()->GetNbins(); ++i) {
161  H2->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i-1));
162  }
163 
164  JManager<int, TH2D> H3((TH2D*) H2->Clone("H[%].doubles"));
165 
166 
168 
169  counter_type counter = 0;
170 
171  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
172 
173  JTreeScanner_t in(*i);
174 
176 
177  for (JTreeScanner_t::iterator event = in.begin(); event != in.end() && counter != inputFile.getLimit(); ++event, ++counter) {
178 
179  if (counter%1000 == 0) {
180  STATUS("event " << setw(8) << counter << '\r'); DEBUG(endl);
181  }
182 
183  const JEvent evt = *event;
184 
185  JQuantile Q1("", true);
186 
187  for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
188  Q1.put(i->getToE());
189  }
190 
191  HA[evt.getID()]->Fill((double) evt.size());
192  HB[evt.getID()]->Fill((double) evt.getOverlays());
193 
194  if (toe.has(evt.getID())) {
195  HC[evt.getID()]->Fill(log10(evt.begin()->getToE() - toe[evt.getID()]));
196  }
197 
198  HD[evt.getID()]->Fill(Q1.getSTDev());
199  HE[evt.getID()]->Fill(Q1.getQuantile(Q, JQuantile::symmetric_t));
200 
201  toe[evt.getID()] = evt.begin()->getToE();
202 
203  TH1D* hq = HQ[evt.getID()];
204  TH2D* hr = HR[evt.getID()];
205  TH1D* h1 = H1[evt.getID()];
206  TH2D* h2 = H2[evt.getID()];
207  TH2D* h3 = H3[evt.getID()];
208 
209  for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
210  hq->Fill(log10(i->getQ()));
211  }
212 
213  if (emitters.has(evt.getID())) {
214 
215  const JPosition3D& p0 = emitters[evt.getID()];
216 
217  for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
218  if (receivers.has(i->getID())) {
219  HR->Fill(log10(p0.getDistance(receivers[i->getID()])), log10(i->getQ()));
220  hr->Fill(log10(p0.getDistance(receivers[i->getID()])), log10(i->getQ()));
221  }
222  }
223  }
224 
225  map<int, set<double> > buffer;
226 
227  for (JEvent::const_iterator i = evt.begin(); i != evt.end(); ++i) {
228  buffer[i->getID()].insert(i->getQ());
229  }
230 
231  for (map<int, set<double> >::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
232 
233  h1->Fill((double) i->second.size());
234 
235  if (receivers.has(i->first)) {
236 
237  const JLocation& location = receivers[i->first];
238 
239  h2->Fill((double) string.getIndex(location.getString()), (double) location.getFloor(), 1.0);
240 
241  if (i->second.size() >= 2u) {
242  h3->Fill((double) string.getIndex(location.getString()), (double) location.getFloor(), 1.0);
243  }
244  }
245  }
246  }
247  }
248  STATUS(endl);
249 
250 
251  TFile out(outputFile.c_str(), "recreate");
252 
253  out << HA << HB << HC << HD << HE << HQ << HR << *HR << H1 << H2 << H3;
254 
255  out.Write();
256  out.Close();
257 }
Utility class to parse command line options.
Definition: JParser.hh:1517
General exception.
Definition: JException.hh:23
int getOverlays() const
Get overlays.
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
Definition: Main.cc:15
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:266
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:83
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.
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:1993
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:319
then awk string
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
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
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.
Utility class to parse command line options.
Acoustic event.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
do set_variable DETECTOR_TXT $WORKDIR detector
int getID() const
Get identifier.
double u[N+1]
Definition: JPolint.hh:776
Data structure for tripod.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] 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:46
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