Jpp  18.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPlatypus.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <set>
6 #include <map>
7 #include <deque>
8 #include <algorithm>
9 #include <limits>
10 
11 #include <type_traits>
12 #include <functional>
13 #include <future>
14 #include <mutex>
15 #include <thread>
16 #include <vector>
17 #include <queue>
18 
19 #include "TROOT.h"
20 #include "TFile.h"
21 
22 #include "JLang/JPredicate.hh"
23 #include "JLang/JComparator.hh"
24 #include "JLang/JComparison.hh"
25 #include "JLang/JFileStream.hh"
26 
27 #include "JDetector/JDetector.hh"
29 #include "JDetector/JTripod.hh"
31 #include "JDetector/JModule.hh"
32 #include "JDetector/JHydrophone.hh"
33 
34 #include "JTools/JHashMap.hh"
35 #include "JTools/JRange.hh"
36 #include "JTools/JQuantile.hh"
37 
39 
41 #include "JAcoustics/JEmitter.hh"
43 #include "JAcoustics/JHit.hh"
45 #include "JAcoustics/JKatoomba.hh"
46 #include "JAcoustics/JSuperEvt.hh"
48 #include "JAcoustics/JSupport.hh"
49 
50 #include "Jeep/JContainer.hh"
51 #include "Jeep/JParser.hh"
52 #include "Jeep/JMessage.hh"
53 
54 
55 namespace {
56 
57  using namespace std;
58  using namespace JPP;
59 
61  //typedef JKatoomba<JAbstractMinimiser> JKatoomba_t;
62 
63  /**
64  * Thread pool for global fits.
65  */
66  class JPlatypus {
67  public:
68  /**
69  * Constructor.
70  *
71  * \param katoomba global fit
72  * \param ns number of threads
73  */
74  JPlatypus(const JKatoomba_t& katoomba, const size_t ns) :
75  stop(false)
76  {
77  using namespace std;
78 
79  for (size_t i = 0; i < ns; ++i) {
80 
81  thread worker([this, katoomba]() {
82 
83  vector<JHit> data;
84 
85  for (JKatoomba_t platypus(katoomba); ; ) {
86 
87  {
88  unique_lock<mutex> lock(in);
89 
90  cv.wait(lock, [this]() { return stop || !input.empty(); });
91 
92  if (stop && input.empty()) {
93  return;
94  }
95 
96  const JSuperEvt& evt = input.front();
97 
98  data.clear();
99 
100  for (JSuperEvt::rx_type::const_iterator hit = evt.rx.begin(); hit != evt.rx.end(); ++hit) {
101 
102  data.push_back(JHit(emitters[hit->id],
103  hit->counter,
104  JLocation(hit->string, hit->floor),
105  hit->toa,
106  parameters.sigma_s,
107  hit->weight));
108  }
109 
110  platypus.value = getModel(evt);
111 
112  input.pop();
113  }
114 
115  const double chi2 = platypus (data.begin(), data.end()) / platypus.estimator->getRho(1.0);
116  const double ndf = getWeight(data.begin(), data.end()) - platypus.value.getN();
117 
118  {
119  unique_lock<mutex> lock(out);
120 
121  Q.put(chi2 / ndf);
122  }
123  }
124  });
125 
126  workers.emplace_back(move(worker));
127  }
128  }
129 
130 
131  /**
132  * Destructor.
133  */
134  ~JPlatypus()
135  {
136  using namespace std;
137 
138  {
139  unique_lock<mutex> lock(in);
140 
141  stop = true;
142  }
143 
144  cv.notify_all();
145 
146  for (auto& worker : workers) {
147  worker.join();
148  }
149  }
150 
151 
152  /**
153  * Get number of pending data.
154  *
155  * \return number of pending data
156  */
157  size_t backlog()
158  {
159  using namespace std;
160 
161  {
162  unique_lock<mutex> lock(in);
163 
164  return input.size();
165  }
166  }
167 
168 
169  /**
170  * Queue super event.
171  *
172  * \param evt super event
173  */
174  void enqueue(const JSuperEvt& evt)
175  {
176  using namespace std;
177 
178  {
179  unique_lock<mutex> lock(in);
180 
181  if (stop) {
182  throw runtime_error("The thread pool has been stopped.");
183  }
184 
185  input.emplace(evt);
186  }
187 
188  cv.notify_one();
189  }
190 
191  static std::string oid; //!< detector identifier
192  static JQuantile Q; //!< chi2/NDF
193  static JFitParameters parameters; //!< fit parameters
194  static::JHashMap<int, JEmitter> emitters; //!< emitters
195 
196 
197  private:
198  vector<thread> workers;
199  queue <JSuperEvt> input;
200  mutex in;
201  mutex out;
202  condition_variable cv;
203  bool stop;
204  };
205 
206  std::string JPlatypus::oid = "";
209  JHashMap<int, JEmitter> JPlatypus::emitters;
210 }
211 
212 
213 /**
214  * \file
215  *
216  * Application to make a global fit of the detector geometry to acoustic data.\n
217  * \author mdejong
218  */
219 int main(int argc, char **argv)
220 {
221  using namespace std;
222  using namespace JPP;
223 
224  typedef JContainer< vector<JTripod> > tripods_container;
225  typedef JContainer< vector<JTransmitter> > transmitters_container;
226  typedef JContainer< vector<JHydrophone> > hydrophones_container;
227 
229  string detectorFile;
230  JLimit_t& numberOfEvents = inputFile.getLimit();
231  JSoundVelocity V = getSoundVelocity; // default sound velocity
232  tripods_container tripods; // tripods
233  transmitters_container transmitters; // transmitters
234  hydrophones_container hydrophones; // hydrophones
235  size_t jobs; // number of parallel jobs
236  int debug;
237 
238  try {
239 
240  JParser<> zap("Application to fit position calibration model to acoustic data.");
241 
242  zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
243  zap['n'] = make_field(numberOfEvents) = JLimit::max();
244  zap['a'] = make_field(detectorFile);
246  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
247  zap['T'] = make_field(tripods, "tripod data");
248  zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
249  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
250  zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
251  zap['N'] = make_field(jobs, "number of parallel jobs") = 1;
252  zap['d'] = make_field(debug) = 1;
253 
254  zap(argc, argv);
255  }
256  catch(const exception &error) {
257  FATAL(error.what() << endl);
258  }
259 
260  ROOT::EnableThreadSafety();
261 
262  const int sleep_us = 100; // sleep time [us]
263 
265 
266  try {
267  load(detectorFile, detector);
268  }
269  catch(const JException& error) {
270  FATAL(error);
271  }
272 
273  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
274  JPlatypus::emitters[i->getID()] = JEmitter(i->getID(),
275  i->getUTMPosition() - detector.getUTMPosition());
276  }
277 
278  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
279  try {
280  JPlatypus::emitters[i->getID()] = JEmitter(i->getID(),
281  i->getPosition() + detector.getModule(i->getLocation()).getPosition());
282  }
283  catch(const exception&) {
284  continue; // if no module available, discard transmitter
285  }
286  }
287 
288  V.set(detector.getUTMZ()); // sound velocity at detector depth
289 
290  JGeometry geometry(detector, hydrophones);
291  JKatoomba_t katoomba(geometry, V, JPlatypus::parameters.option);
292 
295 
297 
298  try {
299 
300  JPlatypus platypus(katoomba, jobs);
301 
302  while (inputFile.hasNext()) {
303 
304  STATUS("event: " << setw(8) << inputFile.getCounter() << '\r'); DEBUG(endl);
305 
306  const JSuperEvt* evt = inputFile.next();
307 
308  while (platypus.backlog() > jobs) {
309  this_thread::sleep_for(chrono::microseconds(sleep_us));
310  }
311 
312  platypus.enqueue(*evt);
313  }
314  STATUS(endl);
315  }
316  catch(const exception& error) {
317  FATAL("main " << error.what());
318  }
319 
320  JFileOutputStream(3) << SCIENTIFIC(1,10) << JPlatypus::Q.getMean(numeric_limits<float>::max()) << endl;
321 }
Worker class for complete fit procedure of acoustic model.
Definition: JKatoomba_t.hh:29
Utility class to parse command line options.
Definition: JParser.hh:1514
Acoustic hit.
General exception.
Definition: JException.hh:23
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
Definition: Main.cc:15
JKatoomba< JEstimator > estimator
Definition: JKatoomba_t.hh:155
double getWeight(T __begin, T __end)
Get total weight of data points.
Definition: JKatoomba.hh:58
Sound velocity.
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition: JMechanics.hh:243
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
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
std::vector< size_t > ns
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
then fatal Number of tripods
Definition: JFootprint.sh:45
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
Acoustic emitter.
Data structure for detector geometry and calibration.
Acoustics hit.
Data structure for hydrophone.
then usage $script[port]< option > nPossible stop
Acoustic super event fit toolkit.
Acoustic fit parameters.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
static const JSoundVelocity getSoundVelocity(1541.0,-17.0e-3,-2000.0)
Function object for velocity of sound.
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
Definition: JKatoomba.hh:556
Detector file.
Definition: JHead.hh:226
Data structure for transmitter.
Acoustic emitter.
Definition: JEmitter.hh:27
Logical location of module.
Definition: JLocation.hh:37
Acoustics toolkit.
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:1989
Acoustic event fit.
Streaming of output.
Definition: JFileStream.hh:46
JPosition3D getPosition(const Vec &pos)
Get position.
then awk string
JMODEL::JString getModel(const JFit &fit)
Get model parameters of string.
General purpose messaging.
Implementation for depth dependend velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Acoustic super event fit.
Definition: JSuperEvt.hh:30
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
Fit functions of acoustic model.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
do set_variable DETECTOR_TXT $WORKDIR detector
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
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
Template definition of fit function of acoustic model.
Definition: JKatoomba.hh:105
Container I/O.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Data structure for optical module.