60 inline double getWeight(
T __begin,
T __end)
64 for (
T i = __begin; i != __end; ++i) {
79 int main(
int argc,
char **argv)
95 transmitters_container transmitters;
96 hydrophones_container hydrophones;
100 disable_container disable;
105 JParser<> zap(
"Application to fit position calibration model to acoustic data.");
107 zap[
'f'] =
make_field(inputFile,
"output of JAcousticEventBuilder[.sh]");
109 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
118 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
124 catch(
const exception &error) {
125 FATAL(error.what() << endl);
145 for (JDetector::const_iterator i =
detector.begin(); i !=
detector.end(); ++i) {
146 receivers[i->getID()] = i->getLocation();
149 for (tripods_container::const_iterator i =
tripods.begin(); i !=
tripods.end(); ++i) {
150 emitters[i->getID()] =
JEmitter(i->getID(),
151 i->getUTMPosition() -
detector.getUTMPosition());
154 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
156 emitters[i->getID()] =
JEmitter(i->getID(),
159 catch(
const exception&) {
187 TH1D h0(
"chi2/NDF", NULL, 50000, 0.0, 1000.0);
188 TH1D
h1(
"h1", NULL, 51, -0.5, 50.5);
189 TH1D hn(
"hn", NULL, 100, 0.0, 6.0);
193 range.getLength() + 1,
range.getLowerLimit() - 0.5,
range.getUpperLimit() + 0.5));
197 range.getLength() + 1,
range.getLowerLimit() - 0.5,
range.getUpperLimit() + 0.5));
199 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
200 HA->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
201 HB->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
209 for (
const string& file_name : inputFile) {
219 else if (oid != evt->
getOID())
220 FATAL(
"Invalid detector identifier " << evt->
getOID() <<
" != " << oid << endl);
223 zmap[evt->begin()->getToE()] = file_name;
232 inputFile.push_back(i->second);
242 int counter[] = { 0, 0 };
244 typedef deque<JEvent> buffer_type;
246 for (buffer_type zbuf; inputFile.hasNext(); ) {
248 STATUS(inputFile.getFilename() <<
'\r');
DEBUG(endl);
252 for (
const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
254 const JEvent* evt = inputFile.next();
256 zbuf.push_back(*evt);
259 sort(zbuf.begin(), zbuf.end());
261 for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
263 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() +
parameters.Tmax_s; ) {}
265 if (q == zbuf.end()) {
267 if (inputFile.hasNext()) {
269 zbuf.erase(zbuf.begin(), p);
281 for (buffer_type::const_iterator i = p; i != q; ++i) {
282 numberOfPings[i->getID()] += 1;
286 DEBUG(
"Number of pings " << setw(2) << i->first <<
' ' << setw(3) << i->second << endl);
289 int minimum_number_of_pings = numeric_limits<int>::max();
292 minimum_number_of_pings = min(minimum_number_of_pings, i->second);
297 for (buffer_type::iterator i = p; i != q; ++i) {
299 if (emitters.has(i->getID())) {
301 const JEmitter& emitter = emitters[i->getID()];
303 const double signal = (unify ? (double) minimum_number_of_pings / (
double) numberOfPings[i->getID()] : 1.0);
309 for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
312 buffer[hit->getID()].push_back(*hit);
318 if (receivers.has(ps->first) && geometry.
hasLocation(receivers[ps->first])) {
322 if (ps->second.begin()->getQ() >=
parameters.Qmin) {
326 receivers[ps->first],
335 DEBUG(
"hit: " << *hit << endl);
341 double chi2 = evaluator(result, data.begin(), data.end());
342 double ndf = getWeight(data.begin(), data.end()) - result.
getN();
344 DEBUG(
"prefit:" << endl
345 <<
FIXED(12,3) << chi2 <<
'/' <<
FIXED(6,1) << ndf << endl
350 HA[hit->getID()]->Fill(geometry.
getIndex(hit->getString()), hit->getFloor(), 1.0);
358 double chi2_old = evaluator(result, data.begin(), data.end());
368 const double x = fabs(hit->getValue() - estimator.
getToA(result, *hit)) / hit->sigma;
378 DEBUG(
"Remove outlier " << out->getEKey() <<
' ' << out->getLocation() <<
' ' <<
FIXED(12,3) << xmax <<
"?" << flush);
382 iter_swap(out, --__end);
384 result = estimator(data.begin(), __end);
386 double chi2_new = evaluator(result, data.begin(), __end);
390 DEBUG(
" yes " <<
FIXED(12,3) << (chi2_old - chi2_new) << endl);
400 DEBUG(
" no " << endl);
402 result = estimator(data.begin(), ++__end);
411 DEBUG(
"hit: " << *hit << endl);
414 DEBUG(
"prefit:" << endl
418 HB[hit->getID()]->Fill(geometry.
getIndex(hit->getString()), hit->getFloor(), 1.0);
422 chi2 = numeric_limits<double>::max();
428 result = estimator(data.begin(), data.end());
429 chi2 = evaluator(result, data.begin(), data.end());
438 chi2 = simplex(data.begin(), data.end()) / simplex.
estimator->getRho(1.0);
439 result = simplex.
value;
448 chi2 = gandalf(data.begin(), data.end()) / gandalf.
estimator->getRho(1.0);
449 result = gandalf.
value;
458 ndf = getWeight(data.begin(), data.end()) - result.
getN();
460 DEBUG(
"result:" << endl
461 <<
FIXED(12,3) << chi2 <<
'/' <<
FIXED(6,1) << ndf << endl
472 data.rbegin()->getValue(),
481 for (buffer_type::iterator i = p; i != q; ++i) {
487 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
498 WARNING(endl <<
"Event not written: " <<
FIXED(12,3) << chi2 <<
'/' <<
FIXED(6,1) << ndf << endl);
507 STATUS(
"Number of events written / rejected: " << counter[0] <<
" / " << counter[1] << endl);
static int debug
debug level (default is off).
Utility class to parse command line options.
size_t getN() const
Get number of fit parameters.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
int main(int argc, char *argv[])
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
static JDetectorMechanics getMechanics
Function object to get string mechanics.
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 elements.
ROOT TTree parameter settings.
Recording of objects on file according a format that follows from the file name extension.
Template specialisation of fit function of acoustic model based on linear approximation.
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
*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
#define MAKE_CSTRING(A)
Make C-string.
then for HISTOGRAM in h0 h1
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Data structure for detector geometry and calibration.
Data structure for hydrophone.
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.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Data structure for transmitter.
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
set_variable E_E log10(E_{fit}/E_{#mu})"
const std::string & getOID() const
Get detector identifier.
static int debug
debug level
do set_variable OUTPUT_DIRECTORY $WORKDIR T
JPosition3D getPosition(const Vec &pos)
Get position.
double getQ(const double D_m, const double f_kHz, const double d_m)
Get relative quality for given frequency at given distance.
General purpose messaging.
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Implementation for depth dependend velocity of sound.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
JACOUSTICS::JModel::emitter_type emitter
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.
Utility class to parse command line options.
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Acoustic transmission identifier.
double getToA(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-arrival for given hit.
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
Custom probability density function of time-of-arrival.
const JLimit & getLimit() const
Get limit.
Fit functions of acoustic model.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for tripod.
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 JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Acoustic transmission identifier.
Template definition of fit function of acoustic model.
#define DEBUG(A)
Message macros.
Data structure for optical module.