Jpp 20.0.0-rc.9-29-gccc23c492-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JCanberra.cc File Reference

Example program to plot acoustic fit results. More...

#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TFitResult.h"
#include "JROOT/JMinimizer.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.hh"
#include "JDetector/JHydrophone.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JModuleRouter.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JManager.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JGeometry.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_t.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to plot acoustic fit results.

Author
mdejong

Definition in file JCanberra.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 53 of file JCanberra.cc.

54{
55 using namespace std;
56 using namespace JPP;
57
58 typedef JContainer< vector<JTripod> > tripods_container;
59 typedef JContainer< vector<JTransmitter> > transmitters_container;
60 typedef JContainer< vector<JHydrophone> > hydrophones_container;
61 typedef JContainer< set<JTransmission_t> > disable_container;
62
63 JMultipleFileScanner<JEvt> inputFile;
64 JLimit_t& numberOfEvents = inputFile.getLimit();
65 string detectorFile;
66 string outputFile;
67 JSoundVelocity V = getSoundVelocity; // default sound velocity
68 tripods_container tripods; // tripods
69 transmitters_container transmitters; // transmitters
70 hydrophones_container hydrophones; // hydrophones
71 int id; // emitter identifier
72 disable_container disable; // disable tansmissions
73 bool revert;
74 string option;
75 int numberOfEntries = 21;
76 int debug;
77
78 try {
79
80 JParser<> zap("Example program to plot acoustic fit results.");
81
82 JProperties properties;
83
84 properties.insert(gmake_property(numberOfEntries));
85
86 zap['f'] = make_field(inputFile, "input file (output of JKatoomba[.sh])");
87 zap['n'] = make_field(numberOfEvents) = JLimit::max();
88 zap['a'] = make_field(detectorFile);
89 zap['o'] = make_field(outputFile) = "canberra.root";
90 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
91 zap['T'] = make_field(tripods, "tripod data");
92 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
93 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
94 zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
95 zap['E'] = make_field(id, "emitter identifier (-1 = all)") = -1;
96 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
97 zap['r'] = make_field(revert, "revert disabling of transmissions");
98 zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "LS";
99 zap['@'] = make_field(properties) = JPARSER::initialised();
100 zap['d'] = make_field(debug) = 2;
101
102 zap(argc, argv);
103 }
104 catch(const exception &error) {
105 FATAL(error.what() << endl);
106 }
107
108
109 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
110
111 JDetector detector;
112
113 try {
114 load(detectorFile, detector);
115 }
116 catch(const JException& error) {
117 FATAL(error);
118 }
119
120 JHashMap<int, JLocation> receivers;
121 JHashMap<int, JEmitter> emitters;
122
123 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
124 receivers[i->getID()] = i->getLocation();
125 }
126
127 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
128 emitters[i->getID()] = JEmitter(i->getID(),
129 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()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition());
135 }
136 catch(const exception&) {} // if no module available, discard transmitter
137 }
138
139 const JGeometry geometry(detector, hydrophones);
140
141 V.set(detector.getUTMZ());
142
143
144 JManager<JTransmission_t, TH1D> H2(new TH1D("[%].toa", NULL, 1000, -1.0e-2, +1.0e-2));
145
146 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
147 for (JHashMap<int, JEmitter>::const_iterator i = emitters.begin(); i != emitters.end(); ++i) {
148 if (id == i->first || id == -1) {
149 H2[JTransmission_t(i->first, module->getID())];
150 }
151 }
152 }
153
154 typedef JTreeScanner<JEvent, JEvent::JEvaluator> JTreeScanner_t;
155
156 JTreeScanner_t in(inputFile);
157
158 JTreeScanner_t::iterator p = in.begin();
159
160 while (inputFile.hasNext()) {
161
162 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
163
164 const JEvt* evt = inputFile.next();
165 const JModel model = getModel(*evt);
166
167 DEBUG("model" << endl << model << endl);
168
169 for ( ; p != in.end() && p-> begin()->getToA() < evt->UNIXTimeStart - 0.5; ++p) {}
170
171 JTreeScanner_t::iterator q = p;
172
173 for ( ; q != in.end() && q->rbegin()->getToA() <= evt->UNIXTimeStop + 0.5; ++q) {}
174
175 DEBUG("events " << distance(p, q) << endl);
176
177 for (JTreeScanner_t::iterator i = p; i != q; ++i) {
178
179 if (id == i->getID() || id == -1) {
180
181 if (emitters.has(i->getID())) {
182
183 const JEmitter& emitter = emitters[i->getID()];
184
185 for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
186
187 if (receivers.has(hit->getID()) && geometry.hasLocation(receivers[hit->getID()])) {
188
189 if ((disable.count(JTransmission_t(i->getID(), hit->getID())) == 0 &&
190 disable.count(JTransmission_t(-1, hit->getID())) == 0) == !revert) {
191
192 const JLocation& location = receivers[hit->getID()];
193
194 if (geometry .has(location.getString()) &&
195 model.string.has(location.getString())) {
196
197 const JGEOMETRY::JString& string = geometry [location.getString()];
198 const JMODEL ::JString& parameters = model.string[location.getString()];
199 const JPosition3D position = string.getPosition(parameters, location.getFloor());
200
201 const double D = emitter.getDistance(position);
202 const double Vi = V.getInverseVelocity(D, emitter.getZ(), position.getZ());
203 const double toa = hit->getToE() + D * Vi;
204
205 H2[JTransmission_t(i->getID(), hit->getID())]->Fill(hit->getToA() - toa);
206 }
207 }
208 }
209 }
210 }
211 }
212 }
213
214 p = q;
215 }
216 STATUS(endl);
217
218
219 const floor_range range = getRangeOfFloors(detector);
220
221 JManager<int, TH2D> HA(new TH2D("[%].mean", NULL,
222 getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
223 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
224
225 JManager<int, TH2D> HB(new TH2D("[%].sigma", NULL,
226 getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
227 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
228
229 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
230 HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
231 HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
232 }
233
234 const JModuleRouter router(detector);
235
236 TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))");
237
238 const vector<Double_t> Q = { 0.33, 0.55, 0.66 };
239
240 for (JManager<JTransmission_t, TH1D>::iterator i = H2.begin(); i != H2.end(); ++i) {
241
242 TH1* h1 = i->second;
243
244 const JLocation location = router.getModule(i->first.rx).getLocation();
245
246 if (h1->GetEntries() >= numberOfEntries) {
247
248 vector<Double_t> X(Q.size(), 0.0);
249
250 h1->GetQuantiles(Q.size(), X.data(), Q.data());
251
252 f1.SetParameter(0, h1->GetMaximum());
253 f1.SetParameter(1, X[1]);
254 f1.SetParameter(2, 0.5*(X[2] - X[0]));
255
256 f1.SetParError(0, 0.1);
257 f1.SetParError(1, 1.0e-6);
258 f1.SetParError(2, 1.0e-6);
259
260 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]));
261
262 if (result.Get() == NULL) {
263
264 ERROR("Invalid TFitResultPtr " << h1->GetName() << endl);
265
266 } else {
267
268 HA[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), f1.GetParameter(1) * 1.0e3);
269 HB[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), f1.GetParameter(2) * 1.0e3);
270 }
271
272 } else {
273
274 HA[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), numeric_limits<double>::lowest());
275 }
276 }
277
278
279 TFile out(outputFile.c_str(), "recreate");
280
281 out << H2 << HA << HB;
282
283 out.Write();
284 out.Close();
285}
JACOUSTICS::tripods_container
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:79
JACOUSTICS::transmitters_container
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:81
JACOUSTICS::hydrophones_container
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:80
JFIT::model
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition JGandalf.hh:57
JACOUSTICS::JEmitter
Acoustic emitter.
Definition JEmitter.hh:30
JACOUSTICS::JEvt
Acoustic event fit.
Definition JEvt.hh:310
JACOUSTICS::JHead::UNIXTimeStop
double UNIXTimeStop
stop time
Definition JEvt.hh:262
JACOUSTICS::JHead::UNIXTimeStart
double UNIXTimeStart
start time
Definition JEvt.hh:261
JACOUSTICS::JModel
Model for fit to acoustics data.
Definition JAcoustics/JModel.hh:480
JACOUSTICS::JSoundVelocity
Implementation for depth dependend velocity of sound.
Definition JSoundVelocity.hh:37
JACOUSTICS::JSoundVelocity::set
JSoundVelocity & set(const double z0)
Set depth.
Definition JSoundVelocity.hh:62
JACOUSTICS::JSoundVelocity::getInverseVelocity
virtual double getInverseVelocity(const double D_m, const double z1, const double z2) const override
Get inverse velocity of sound.
Definition JSoundVelocity.hh:175
JACOUSTICS::JTransmission_t
Acoustic transmission identifier.
Definition JTransmission_t.hh:28
JACOUSTICS::getModel
Auxiliary data structure to convert event to model.
Definition JEvtToolkit.hh:126
Generated by doxygen 1.10.0