Jpp  16.0.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JCanberra.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <vector>
4 #include <set>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TF1.h"
11 #include "TFitResult.h"
12 
13 #include "JDetector/JDetector.hh"
15 #include "JDetector/JTripod.hh"
16 #include "JDetector/JModule.hh"
17 #include "JDetector/JHydrophone.hh"
19 
21 #include "JSupport/JTreeScanner.hh"
22 
23 #include "JROOT/JRootToolkit.hh"
24 #include "JROOT/JManager.hh"
25 
27 #include "JAcoustics/JGeometry.hh"
28 #include "JAcoustics/JEmitter.hh"
30 #include "JAcoustics/JHit.hh"
31 #include "JAcoustics/JEvent.hh"
32 #include "JAcoustics/JEvt.hh"
34 #include "JAcoustics/JSupport.hh"
36 
37 #include "Jeep/JContainer.hh"
38 #include "Jeep/JPrint.hh"
39 #include "Jeep/JParser.hh"
40 #include "Jeep/JMessage.hh"
41 
42 
43 /**
44  * \file
45  *
46  * Example program to plot acoustic fit results.
47  * \author mdejong
48  */
49 int main(int argc, char **argv)
50 {
51  using namespace std;
52  using namespace JPP;
53 
54  typedef JContainer< vector<JTripod> > tripods_container;
55  typedef JContainer< vector<JHydrophone> > hydrophones_container;
56  typedef JContainer< set<JTransmission_t> > disable_container;
57 
59  JLimit_t& numberOfEvents = inputFile.getLimit();
60  string detectorFile;
61  string outputFile;
62  JSoundVelocity V = getSoundVelocity; // default sound velocity
63  tripods_container tripods; // tripods
64  hydrophones_container hydrophones; // hydrophones
65  int id; // emitter identifier
66  disable_container disable; // disable tansmissions
67  string option;
68  int debug;
69 
70  try {
71 
72  JParser<> zap("Example program to plot acoustic fit results.");
73 
74  zap['f'] = make_field(inputFile, "input file (output of JKatoomba)");
75  zap['n'] = make_field(numberOfEvents) = JLimit::max();
76  zap['a'] = make_field(detectorFile);
77  zap['o'] = make_field(outputFile) = "canberra.root";
78  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
79  zap['T'] = make_field(tripods, "tripod data");
80  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
81  zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
82  zap['E'] = make_field(id, "emitter identifier (-1 = all)") = -1;
83  zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
84  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "LS";
85  zap['d'] = make_field(debug) = 2;
86 
87  zap(argc, argv);
88  }
89  catch(const exception &error) {
90  FATAL(error.what() << endl);
91  }
92 
93 
94  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
95 
97 
98  try {
99  load(detectorFile, detector);
100  }
101  catch(const JException& error) {
102  FATAL(error);
103  }
104 
105  JHashMap<int, JLocation> receivers;
106  JHashMap<int, JEmitter> emitters;
107 
108  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
109  receivers[i->getID()] = i->getLocation();
110  }
111 
112  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
113  emitters[i->getID()] = JEmitter(i->getID(),
114  i->getUTMPosition() - detector.getUTMPosition());
115  }
116 
117  const JGeometry geometry(detector, hydrophones);
118 
119  V.set(detector.getUTMZ());
120 
121 
122  JManager<JTransmission_t, TH1D> H2(new TH1D("[%].toa", NULL, 1000, -1.0e-2, +1.0e-2));
123 
124 
126 
127  JTreeScanner_t in(inputFile);
128 
129  JTreeScanner_t::iterator p = in.begin();
130 
131  while (inputFile.hasNext()) {
132 
133  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
134 
135  const JEvt* evt = inputFile.next();
136  const JModel model = getModel(*evt);
137 
138  DEBUG("model" << endl << model << endl);
139 
140  for ( ; p != in.end() && p-> begin()->getToA() < evt->UNIXTimeStart - 0.5; ++p) {}
141 
142  JTreeScanner_t::iterator q = p;
143 
144  for ( ; q != in.end() && q->rbegin()->getToA() <= evt->UNIXTimeStop + 0.5; ++q) {}
145 
146  DEBUG("events " << distance(p, q) << endl);
147 
148  for (JTreeScanner_t::iterator i = p; i != q; ++i) {
149 
150  if (id == i->getID() || id == -1) {
151 
152  if (emitters.has(i->getID())) {
153 
154  const JEmitter& emitter = emitters[i->getID()];
155 
156  for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
157 
158  if (receivers.has(hit->getID())) {
159 
160  if (disable.count(JTransmission_t(i->getID(), hit->getID())) == 0 &&
161  disable.count(JTransmission_t(-1, hit->getID())) == 0) {
162 
163  const JLocation& location = receivers[hit->getID()];
164 
165  if (geometry .has(location.getString()) &&
166  model.string.has(location.getString())) {
167 
168  const JGEOMETRY::JString& string = geometry [location.getString()];
169  const JMODEL ::JString& parameters = model.string[location.getString()];
170  const JPosition3D position = string.getPosition(parameters, location.getFloor());
171 
172  const double D = emitter.getDistance(position);
173  const double Vi = V.getInverseVelocity(D, emitter.getZ(), position.getZ());
174  const double toa = hit->getToE() + D * Vi;
175 
176  H2[JTransmission_t(i->getID(), hit->getID())]->Fill(hit->getToA() - toa);
177  }
178  }
179  }
180  }
181  }
182  }
183  }
184 
185  p = q;
186  }
187  STATUS(endl);
188 
189 
191 
192  JManager<int, TH2D> HA(new TH2D("[%].mean", NULL,
194  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
195 
196  JManager<int, TH2D> HB(new TH2D("[%].sigma", NULL,
198  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
199 
200  for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
201  HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
202  HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
203  }
204 
205  const JModuleRouter router(detector);
206 
207  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))");
208 
209  const vector<Double_t> Q = { 0.33, 0.55, 0.66 };
210 
211  for (JManager<JTransmission_t, TH1D>::iterator i = H2.begin(); i != H2.end(); ++i) {
212 
213  TH1* h1 = i->second;
214 
215  vector<Double_t> X(Q.size(), 0.0);
216 
217  h1->GetQuantiles(Q.size(), X.data(), Q.data());
218 
219  f1.SetParameter(0, h1->GetMaximum());
220  f1.SetParameter(1, X[1]);
221  f1.SetParameter(2, 0.5*(X[2] - X[0]));
222 
223  TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", X[1] - 5.0 * (X[2] - X[0]), X[1] + 5.0 * (X[2] - X[0]));
224 
225  if (result.Get() == NULL) {
226 
227  ERROR("Invalid TFitResultPtr " << h1->GetName() << endl);
228 
229  } else {
230 
231  const JLocation location = router.getModule(i->first.rx).getLocation();
232 
233  HA[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), f1.GetParameter(1) * 1.0e3);
234  HB[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), f1.GetParameter(2) * 1.0e3);
235  }
236  }
237 
238 
239  TFile out(outputFile.c_str(), "recreate");
240 
241  out << H2 << HA << HB;
242 
243  out.Write();
244  out.Close();
245 }
Utility class to parse command line options.
Definition: JParser.hh:1500
Acoustic hit.
General exception.
Definition: JException.hh:23
Q(UTCMax_s-UTCMin_s)-livetime_s
JModel getModel(const JEvt &evt)
Get model.
debug
Definition: JMessage.hh:29
int main(int argc, char *argv[])
Definition: Main.cc:15
Sound velocity.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
const JModule & getModule(const JObjectID &id) const
Get module parameters.
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition: JMechanics.hh:243
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
#define STATUS(A)
Definition: JMessage.hh:63
ROOT TTree parameter settings.
Detector data structure.
Definition: JDetector.hh:89
Router for direct addressing of module data in detector data structure.
Acoustic geometries.
*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
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Acoustic event.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Dynamic ROOT object management.
int getIndex(const T &value) const
Get index of given value.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
Acoustic emitter.
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Template definition for direct access of elements in ROOT TChain.
Definition: JTreeScanner.hh:91
Data structure for detector geometry and calibration.
Data structure for hydrophone.
Model for fit to acoustics data.
Acoustic event fit.
double UNIXTimeStop
stop time
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.
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:224
Acoustic event fit.
Acoustic emitter.
Definition: JEmitter.hh:27
Logical location of module.
Definition: JLocation.hh:37
Acoustics toolkit.
const JLocation & getLocation() const
Get location.
Definition: JLocation.hh:69
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:1961
return result
Definition: JPolint.hh:743
#define ERROR(A)
Definition: JMessage.hh:66
then break fi done getCenter read X Y Z let X
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition: JManager.hh:295
int debug
debug level
Definition: JSirene.cc:63
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
General purpose messaging.
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
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...
Base class for JTreeScanner.
Definition: JTreeScanner.hh:52
Direct access to module in detector data structure.
JACOUSTICS::JModel::string_type string
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
int getString() const
Get string number.
Definition: JLocation.hh:134
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double UNIXTimeStart
start time
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
Acoustic transmission identifier.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for tripod.
then fatal Not enough tripods
Acoustic event fit.
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:42
Acoustic transmission identifier.
bool has(const T &value) const
Test whether given value is present.
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Container I/O.
Data structure for optical module.