Jpp master_rocky-44-g75b7c4f75
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}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:69
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition JProperties.hh:1293
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition PhysicsEvent.hh:434
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JLocation
Logical location of module.
Definition JLocation.hh:40
JDETECTOR::JLocation::getLocation
const JLocation & getLocation() const
Get location.
Definition JLocation.hh:70
JDETECTOR::JLocation::getFloor
int getFloor() const
Get floor number.
Definition JLocation.hh:146
JDETECTOR::JLocation::getString
int getString() const
Get string number.
Definition JLocation.hh:135
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition JModuleRouter.hh:38
JEEP::JProperties
Utility class to parse parameter values.
Definition JProperties.hh:501
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition JPosition3D.hh:130
JGEOMETRY3D::JVector3D::getDistance
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition JVector3D.hh:270
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition JVector3D.hh:115
JLANG::JException
General exception.
Definition JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition JMultipleFileScanner.hh:306
JSUPPORT::JMultipleFileScanner::hasNext
virtual bool hasNext() override
Check availability of next element.
Definition JMultipleFileScanner.hh:409
JSUPPORT::JMultipleFileScanner::getCounter
counter_type getCounter() const
Get counter.
Definition JMultipleFileScanner.hh:382
JSUPPORT::JMultipleFileScanner::next
virtual const pointer_type & next() override
Get next element.
Definition JMultipleFileScanner.hh:457
JSUPPORT::JTreeScanner_t
Base class for JTreeScanner.
Definition JTreeScanner.hh:56
JSUPPORT::JTreeScanner
Template definition for direct access of elements in ROOT TChain.
Definition JTreeScanner.hh:698
JTOOLS::JHashCollection::has
bool has(const T &value) const
Test whether given value is present.
Definition JHashCollection.hh:309
JTOOLS::JRange::getLength
T getLength() const
Get length (difference between upper and lower limit).
Definition JRange.hh:289
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
JACOUSTICS::tripods_container
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:78
JACOUSTICS::transmitters_container
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:80
JACOUSTICS::hydrophones_container
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:79
JACOUSTICS::getSoundVelocity
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
JACOUSTICS::getMechanics
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition JMechanics.hh:242
JCOMPASS::getModel
JModel getModel(const JEvt &evt)
Get model.
Definition JCompass/JEvtToolkit.hh:51
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition JDetectorToolkit.hh:267
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:505
JDETECTOR::getNumberOfStrings
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Definition JStringCounter.hh:95
JEEP::debug_t
@ debug_t
debug
Definition JMessage.hh:29
JFIT::model
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition JGandalf.hh:56
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JTOOLS::result
return result
Definition JPolint.hh:853
JAANET::detector
Detector file.
Definition JHead.hh:227
JACOUSTICS::JEmitter
Acoustic emitter.
Definition JEmitter.hh:30
JACOUSTICS::JEvt
Acoustic event fit.
Definition JAcoustics/JEvt.hh:304
JACOUSTICS::JHead::UNIXTimeStop
double UNIXTimeStop
stop time
Definition JAcoustics/JEvt.hh:257
JACOUSTICS::JHead::UNIXTimeStart
double UNIXTimeStart
start time
Definition JAcoustics/JEvt.hh:256
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
JEEP::JContainer
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
JTOOLS::JHashMap
General purpose class for hash map of unique keys.
Definition JHashMap.hh:75
JTOOLS::JHashMap::const_iterator
container_type::const_iterator const_iterator
Definition JHashMap.hh:86
Generated by doxygen 1.12.0