56 int main(
int argc,
char **argv)
71 hydrophones_container hydrophones;
75 disable_container disable;
80 JParser<> zap(
"Application to fit position calibration model to acoustic data.");
82 zap[
'f'] =
make_field(inputFile,
"output of JAcousticEventBuilder[.sh]");
84 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
92 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
98 catch(
const exception &error) {
99 FATAL(error.what() << endl);
119 for (JDetector::const_iterator i =
detector.begin(); i !=
detector.end(); ++i) {
120 receivers[i->getID()] = i->getLocation();
123 for (tripods_container::const_iterator i =
tripods.begin(); i !=
tripods.end(); ++i) {
124 emitters[i->getID()] =
JEmitter(i->getID(),
125 i->getUTMPosition() -
detector.getUTMPosition());
148 JKatoomba_t::setOption(
false);
156 TH1D h0(
"chi2/NDF", NULL, 50000, 0.0, 1000.0);
157 TH1D
h1(
"h1", NULL, 51, -0.5, 50.5);
158 TH1D hn(
"hn", NULL, 100, 0.0, 6.0);
162 range.getLength() + 1,
range.getLowerLimit() - 0.5,
range.getUpperLimit() + 0.5));
166 range.getLength() + 1,
range.getLowerLimit() - 0.5,
range.getUpperLimit() + 0.5));
168 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
169 HA->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
170 HB->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
178 for (
const string& file_name : inputFile) {
188 else if (oid != evt->
getOID())
189 FATAL(
"Invalid detector identifier " << evt->
getOID() <<
" != " << oid << endl);
192 zmap[evt->begin()->getToE()] = file_name;
201 inputFile.push_back(i->second);
210 int counter[] = { 0, 0 };
212 typedef deque<JEvent> buffer_type;
214 for (buffer_type zbuf; inputFile.hasNext(); ) {
216 STATUS(inputFile.getFilename() <<
'\r');
DEBUG(endl);
220 for (
const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
222 const JEvent* evt = inputFile.next();
224 zbuf.push_back(*evt);
227 sort(zbuf.begin(), zbuf.end());
229 for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
231 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() +
parameters.Tmax_s; ) {}
233 if (q == zbuf.end()) {
235 if (inputFile.hasNext()) {
237 zbuf.erase(zbuf.begin(), p);
249 for (buffer_type::const_iterator i = p; i != q; ++i) {
250 numberOfPings[i->getID()] += 1;
254 DEBUG(
"Number of pings " << setw(2) << i->first <<
' ' << setw(3) << i->second << endl);
257 int minimum_number_of_pings = numeric_limits<int>::max();
260 minimum_number_of_pings = min(minimum_number_of_pings, i->second);
265 for (buffer_type::iterator i = p; i != q; ++i) {
267 const JEmitter& emitter = emitters[i->getID()];
269 const double signal = (unify ? (double) minimum_number_of_pings / (
double) numberOfPings[i->getID()] : 1.0);
275 for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
278 buffer[hit->getID()].push_back(*hit);
284 if (receivers.has(ps->first) && geometry.
hasLocation(receivers[ps->first])) {
288 if (ps->second.begin()->getQ() >=
parameters.Qmin) {
292 receivers[ps->first],
300 DEBUG(
"hit: " << *hit << endl);
306 DEBUG(
"prefit:" << endl
311 HA[hit->getID()]->Fill(geometry.
getIndex(hit->getString()), hit->getFloor(), 1.0);
319 double chi2_old = evaluator(result, data.begin(), data.end());
329 const double x = fabs(hit->getValue() - estimator.
getToA(result, *hit)) / hit->sigma;
341 iter_swap(out, --__end);
343 result = estimator(data.begin(), __end);
345 double chi2_new = evaluator(result, data.begin(), __end);
349 DEBUG(
"Remove outlier " << __end->getEKey() <<
' ' << __end->getLocation() <<
' ' << xmax << endl);
359 result = estimator(data.begin(), ++__end);
368 DEBUG(
"hit: " << *hit << endl);
371 DEBUG(
"prefit:" << endl
375 HB[hit->getID()]->Fill(geometry.
getIndex(hit->getString()), hit->getFloor(), 1.0);
379 double chi2 = numeric_limits<double>::max();
385 result = estimator(data.begin(), data.end());
386 chi2 = evaluator(result, data.begin(), data.end());
393 chi2 = simplex(data.begin(), data.end()) / simplex.
estimator->getRho(1.0);
394 result = simplex.
value;
403 chi2 = gandalf(data.begin(), data.end()) / gandalf.
estimator->getRho(1.0);
404 result = gandalf.
value;
416 W += hit->getWeight();
419 const int ndf = data.size() - result.
getN();
421 DEBUG(
"result:" << endl
422 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
425 h0.Fill(chi2 / (W - result.
getN()));
431 const JEvt evt =
getEvt(
JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, (W - result.
getN()), chi2),
result);
435 for (buffer_type::iterator i = p; i != q; ++i) {
441 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
452 WARNING(endl <<
"Event not written: " << chi2 <<
'/' << ndf << endl);
461 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.
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...
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
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.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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.
Data structure for tripod.
then fatal Not enough tripods
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Acoustic transmission identifier.
Template definition of fit function of acoustic model.
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt]]
Data structure for optical module.