27 using TMath::PoissonI;
29 static const int NUMBER_OF_DIMENSIONS = 3;
35 typedef JArray<NUMBER_OF_DIMENSIONS, double> point_type;
42 public JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>
48 JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>()
58 template<
class ...Args>
59 model_type(argument_type value,
const Args& ...args) :
60 JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>(value, args...)
70 inline double operator()(
const point_type& point)
const
72 double value = (*this)[0];
74 for (
int i = 0; i != NUMBER_OF_DIMENSIONS;++i) {
75 value *=
gauss(point[i] - (*
this)[2*i + 1], (*
this)[2*i + 2]);
89 static inline double gauss(
const double x,
const double sigma)
91 const double u = x /
sigma;
114 hit_type(
const point_type& point,
const double value,
const double sigma) :
132 inline double chi2(
const model_type& model,
const hit_type& hit)
134 const double u = (model(hit) - hit.value) / hit.sigma;
147 inline double likelihood(
const model_type& model,
const hit_type& hit)
149 return -
log(PoissonI(hit.value, model(hit)));
160 int main(
int argc,
char **argv)
165 typedef JArray<2, double> array_type;
175 JParser<> zap(
"Program to test JSimplex algorithm.");
185 catch(
const exception& error) {
186 FATAL(error.what() << endl);
190 gRandom->SetSeed(seed);
200 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
201 model[i*2 + 1] = gRandom->Uniform(-sigma, +sigma);
202 model[i*2 + 2] = gRandom->Gaus(sigma, 0.2*sigma);
205 const JGrid<double> grid(11, -3.0*sigma, +3.0*sigma);
208 JSimplex<model_type> simplex;
211 JSimplex<model_type>::MAXIMUM_ITERATIONS = 10000;
213 double (*fit)(
const model_type&,
const hit_type&) = likelihood;
218 TH1D h0(
"chi2/NDF", NULL, 200, 0.0, 10.0);
220 H1.push_back(
new TH1D(
MAKE_CSTRING(
"p" << 0), NULL, 101, -20.0, +20.0));
222 for (
size_t i =1; i != model.size(); ++i) {
223 H1.push_back(
new TH1D(
MAKE_CSTRING(
"p" << i), NULL, 101, -0.1, +0.1));
227 for (
int counter = 0; counter != numberOfEvents; ++counter) {
229 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
236 for (JArrayIterator<NUMBER_OF_DIMENSIONS, double>
g1(grid);
g1; ++
g1) {
238 const JArray<NUMBER_OF_DIMENSIONS, double>
p1 = *
g1;
240 const double value = gRandom->Poisson(model(p1));
241 const double sigma = (value > 1.0 ? sqrt(value) : 0.7);
243 data.push_back(hit_type(p1, value, sigma));
249 JQuantile
Q[NUMBER_OF_DIMENSIONS];
253 for (
const auto& hit : data) {
255 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
256 Q[i].put(hit[i], hit.value);
259 if (hit.value > value) {
261 simplex.value[0] = hit.value;
263 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
264 simplex.value[2*i + 1] = hit[i];
271 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
272 simplex.value[2*i + 2] = 1.0 * Q[i].getSTDev();
278 simplex.step.resize(2*NUMBER_OF_DIMENSIONS + 1);
282 simplex.step[0] = model_type();
283 simplex.step[0][0] = 1.0;
285 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
287 simplex.step[2*i + 1] = model_type();
288 simplex.step[2*i + 1][2*i + 1] = 0.1;
290 simplex.step[2*i + 2] = model_type();
291 simplex.step[2*i + 2][2*i + 2] = 0.1;
294 DEBUG(
"start values " << simplex.value << endl);
296 for (
size_t i = 0; i != simplex.step.size(); ++i) {
297 DEBUG(
"step: " << setw(2) << i <<
' ' << simplex.step[i] << endl);
303 const double chi2 = simplex(fit, data.begin(), data.end());
304 const int ndf = data.size() - simplex.step.size();
309 for (
size_t i = 0; i != model.size(); ++i) {
310 H1[i]->Fill(model[i] - simplex.value[i]);
313 DEBUG(
"chi2 / NDF" <<
FIXED(12,3) << chi2 <<
" / " << ndf << endl);
315 DEBUG(
"final values " << simplex.value << endl);
316 DEBUG(
"model values " << model << endl);
320 for (
const auto& hit : data) {
324 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
325 cout <<
FIXED(7,3) << hit[i] <<
' ';
328 const double value = simplex.value(hit);
330 cout <<
FIXED(7,3) << hit.value <<
' '
331 <<
FIXED(7,3) << value <<
' '
332 <<
FIXED(7,3) << (hit.value - value) / hit.sigma << endl;
346 for (
size_t i = 0; i !=
sizeof(H1)/
sizeof(H1[0]); ++i) {
Utility class to parse command line options.
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
*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.
Auxiliary data structure for floating point format specification.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
General purpose messaging.
Utility class to parse command line options.
then set_variable NUMBER_OF_TESTS else set_variable NUMBER_OF_TESTS fi function gauss()
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
then cat $TRIPOD_INITIAL<< EOF1 256877.5 4743716.7-2438.42 256815.5 4743395.0-2435.53 257096.2 4743636.0-2439.5EOFfiif[[!-f $DETECTOR]];then JEditDetector-a $DETECTOR_INITIAL-s"-1 addz -6.9"-o $DETECTOR--!eval`JPrintDetector-a $DETECTOR-O SUMMARY`for STRING in ${STRINGS[*]};do set_variable MODULE`getModule-a $DETECTOR-L"$STRING 0"`JEditDetector-a $DETECTOR-M"$MODULE setz -2.9"-o $DETECTOR--!donefiif[[!-f $TRIPOD]];then cp-p $TRIPOD_INITIAL $TRIPODfiJAcoustics.sh $DETECTOR_IDcat > acoustics_trigger_parameters txt<< EOFQ=0.0;TMax_s=0.020;numberOfHits=90;EOFJAcousticsEventBuilder.sh $DETECTOR $RUNS[*]INPUT_FILES=(`ls KM3NeT_ ${(l:8::0::0:) DETECTOR_ID}_0 *${^RUNS}_event.root`) cd $WORKDIRif[!$HOMEDIR-ef $WORKDIR];then cp-p $HOMEDIR/$DETECTOR $WORKDIR cp-p $HOMEDIR/${^ACOUSTICS_KEYS}.txt $WORKDIR cp-p $HOMEDIR/${^INPUT_FILES}$WORKDIRfisource $JPP_DIR/examples/JAcoustics/acoustics-fit-toolkit.shtimer_startinitialise stage_b 1 0 100.0e-6 10.0 0.002 0.1 0 > &stage log
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
#define DEBUG(A)
Message macros.
Double_t g1(const Double_t x)
Function.