48 static const int NUMBER_OF_PINGS = 11;
49 static const JWaveform EMITTER_WAVEFORM(2.5e6,
69 int main(
int argc,
char **argv)
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) {
146 i->getUTMPosition() -
detector.getUTMPosition()));
149 for (transmitters_container::const_iterator
i = transmitters.begin();
i != transmitters.end(); ++
i) {
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();
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) {
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) {
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);
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);
static int debug
debug level (default is off).
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.
int main(int argc, char *argv[])
JContainer< std::vector< JTransmitter > > transmitters_container
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)...
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
Data structure for detector geometry and calibration.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Data structure for hydrophone.
static const char WILDCARD
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...
JContainer< std::vector< JHydrophone > > hydrophones_container
Data structure for transmitter.
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})"
static int debug
debug level
Auxiliary class for CPU timing and usage.
JPosition3D getPosition(const Vec &pos)
Get position.
JACOUSTICS::JModel::emission_type emission
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
JContainer< std::vector< JTripod > > tripods_container
General purpose messaging.
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++...
Utility class to parse command line options.
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Acoustic transmission identifier.
Fit functions of acoustic model.
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
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
JFit_t
Enumeration for fit algorithms.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for tripod.
Acoustic transmission identifier.
#define DEBUG(A)
Message macros.
Data structure for optical module.