64 tripods_container tripods;
65 hydrophones_container hydrophones;
69 double background = 1.0e-4;
70 double Tmax_s = 300.0;
88 JParser<> zap(
"Application to fit position calibration model to acoustic data.");
90 zap[
'f'] =
make_field(inputFile,
"output of JAcousticEventBuilder[.sh]");
92 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
100 zap[
'u'] =
make_field(unique,
"one ping per cycle");
105 catch(
const exception &error) {
106 FATAL(error.what() << endl);
125 for (JDetector::const_iterator i =
detector.begin(); i !=
detector.end(); ++i) {
126 receivers[i->getID()] = i->getLocation();
129 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
130 emitters[i->getID()] =
JEmitter(i->getID(),
131 i->getUTMPosition() -
detector.getUTMPosition());
138 i->second.mechanics = mechanics[i->first];
152 simplex.debug =
debug;
153 gandalf.debug =
debug;
159 TH1D h0(
"chi2/NDF", NULL, 5000, 0.0, 100.0);
169 for (Int_t i = 1; i <= ha.GetXaxis()->GetNbins(); ++i) {
170 ha.GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
171 hb.GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
180 while (inputFile.hasNext()) {
182 STATUS(
"input " << setw(6) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
184 const JEvent* evt = inputFile.next();
190 if (oid != evt->
getOID()) {
191 FATAL(
"Invalid detector identifier " << evt->
getOID() <<
" != " << oid << endl);
194 zbuf.push_back(*evt);
198 sort(zbuf.begin(), zbuf.end());
210 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + Tmax_s; ) {}
217 numberOfPings[i->getID()] += 1;
220 int weight = numeric_limits<int>::max();
223 DEBUG(
"Number of pings " << setw(2) << i->first <<
' ' << setw(3) << i->second << endl);
227 if (i->second < weight) {
237 const JEmitter& emitter = emitters[i->getID()];
239 const double signal = (unique ? (double) weight / (
double) numberOfPings[i->getID()] : 1.0);
245 for (
JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
246 buffer[hit->getID()].push_back(*hit);
251 if (receivers.has(ps->first) && geometry.hasLocation(receivers[ps->first])) {
253 sort(ps->second.begin(), ps->second.end(),
make_comparator(&JTransmission::getToA));
255 data.push_back(hit_type(emitter,
257 receivers[ps->first],
265 DEBUG(
"hit: " << *hit << endl);
272 DEBUG(
"prefit:" << endl
277 ha.Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
293 const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;
303 DEBUG(
"Remove outlier " << out->getLocation() <<
' ' << xmax << endl);
307 result = estimator(data.begin(), data.end());
317 hb.Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
321 double chi2 = numeric_limits<double>::max();
327 result = estimator(data.begin(), data.end());
328 chi2 = evaluator(result, data.begin(), data.end());
335 chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
336 result = simplex.value;
343 chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
344 result = gandalf.value;
354 W += hit->getWeight();
357 const int ndf = data.size() - result.
getN();
359 DEBUG(
"result:" << endl
360 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
363 h0.Fill(chi2 / (W - result.
getN()));
367 const JEvt evt =
getEvt(
JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, chi2), result);
375 const double toe = result.
emitter[
JEKey(i->getID(), i->getCounter())].t1;
377 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
Utility class to parse command line options.
size_t getN() const
Get number of fit parameters.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Template specialisation of fit function of acoustic model based on linear approximation.
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
Utility class to parse parameter values.
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
#define MAKE_CSTRING(A)
Make C-string.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary data structure for floating point format specification.
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Model for fit to acoustics data.
Auxiliary class for defining the range of iterations of objects.
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.
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
const std::string & getOID() const
Get detector identifier.
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt &&-f $JPP_LIB/${KEY}_ ${(l:8::0::0:) ID}.txt]]
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Implementation for velocity of sound.
JACOUSTICS::JModel::emitter_type emitter
then $DIR JKatoomba a $DETECTOR o $WORKDIR katoomba root T $TRIPOD n sigma_s
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++...
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Custom probability density function of time-of-arrival.
const JLimit & getLimit() const
Get limit.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector
do for((RUN=${RANGE%%-*};$RUN<=${RANGE##*-};RUN+=1))
std::vector< double > weight
#define DEBUG(A)
Message macros.