86 transmitters_container transmitters;
87 hydrophones_container hydrophones;
92 disable_container disable;
93 waveform_container waveforms;
98 numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
99 waveforms [WILDCARD] = EMITTER_WAVEFORM;
103 JParser<> zap(
"Example application to test fit of model to simulated acoustic data.");
114 zap[
'F'] =
make_field(fit,
"fit type") = linear_t, simplex_t, gandalf_t;
115 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
118 zap[
'B'] =
make_field(background,
"background probability") = 0.0;
124 catch(
const exception &error) {
125 FATAL(error.what() << endl);
129 gRandom->SetSeed(seed);
131 if (numberOfEvents <= 0) {
FATAL(
"Invalid number of events " << numberOfEvents << endl); }
144 for (tripods_container::const_iterator i =
tripods.begin(); i !=
tripods.end(); ++i) {
145 emitters.push_back(
JEmitter(i->getID(),
146 i->getUTMPosition() -
detector.getUTMPosition()));
149 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
151 emitters.push_back(
JEmitter(i->getID(),
154 catch(
const exception&) {
159 if (
detector.empty()) {
FATAL(
"No modules in detector." << endl); }
160 if (emitters.empty()) {
FATAL(
"No emitters in system." << endl); }
165 const double D_m = -
detector.getUTMZ();
180 simplex.debug =
debug;
181 gandalf.debug =
debug;
186 TH1D h0(
"cpu", NULL, 100, 1.0, 7.0);
187 TH1D h1(
"chi2/NDF", NULL, 100, 0.0, 5.0);
195 for (
int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
197 STATUS(
"event: " << setw(10) << number_of_events <<
'\r');
DEBUG(endl);
201 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
206 gRandom->Uniform(-2.0e-2, +2.0e-2));
210 DEBUG(
"Model" << endl << model << endl);
214 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
215 if (geometry.hasLocation(module->getLocation())) {
216 module->set(estimator.geometry[module->getString()].getPosition(model.
string[module->getString()], module->getFloor()));
223 int minimum_number_of_pings = numeric_limits<int>::max();
226 minimum_number_of_pings = min(minimum_number_of_pings,
getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
233 const int number_of_pings =
getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
235 const double weight = (unify ? (double) minimum_number_of_pings / (
double) number_of_pings : 1.0);
237 for (
int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
241 const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
243 model.
emission[emitter->getID()][count] = toe_s;
245 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
247 if (disable.count(
JTransmission_t(emitter->getID(), module->getID())) == 0) {
249 if (geometry.hasLocation(module->getLocation())) {
251 const JWaveform& waveform =
getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
253 const double d_m =
getDistance(module->getPosition(), emitter->getPosition());
254 const double toa_s = toe_s +
V.getTime(d_m, emitter->getZ(), module->getZ());
255 const double Q = waveform.
getQ(D_m, d_m);
261 if (gRandom->Rndm() >= background)
262 t1_s = gRandom->Gaus(toa_s,
parameters.sigma_s);
264 t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
265 toa_s + T_s.getUpperLimit());
267 const JHit hit(*emitter,
269 module->getLocation(),
274 DEBUG(
"hit: " << hit <<
' ' <<
FIXED(7,1) << Q << endl);
289 double chi2 = numeric_limits<double>::max();
295 chi2 = evaluator(result, data.begin(), data.end());
302 chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
303 result = simplex.value;
310 chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
311 result = gandalf.value;
323 W += hit->getWeight();
326 const int ndf = data.size() - result.
getN();
328 DEBUG(
"Final values" << endl
329 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
333 h0.Fill(
log10((
double) timer.usec_wall));
334 h1.Fill(chi2 / (
double) (W - result.
getN()));
337 H2[i->first]->Fill((i->second.tx - result.
string [i->first].tx) * 1.0e3,
338 (i->second.ty - result.
string [i->first].ty) * 1.0e3);
342 H1[i->first.getID()]->Fill(i->second.t1 - result.
emission[i->first].t1);
Utility class to parse command line options.
size_t getN() const
Get number of fit parameters.
Q(UTCMax_s-UTCMin_s)-livetime_s
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Template specialisation of fit function of acoustic model based on linear approximation.
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
*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
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary data structure for floating point format specification.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Model for fit to acoustics data.
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...
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})"
Auxiliary class for CPU timing and usage.
JPosition3D getPosition(const Vec &pos)
Get position.
JACOUSTICS::JModel::emission_type emission
Implementation for depth dependend velocity of sound.
JACOUSTICS::JModel::string_type string
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
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...
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do echo n Creating graphics for string $STRING for((FLOOR=$FIRST_FLOOR;$FLOOR<=$LAST_FLOOR;FLOOR+=1))
do set_variable DETECTOR_TXT $WORKDIR detector
Acoustic transmission identifier.
#define DEBUG(A)
Message macros.