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