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.
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 JAcoustics.sh -- typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
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 usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
Double_t g1(const Double_t x)
Function.