Hi,
Apparently I am having some issues computing integral and quantiles of simple functions (a quarter of a circle). The bug is related to a storage area which is overwritten. It appears from time to time and shows the following :
Fatal in <operator delete>: storage area overwritten
aborting
#0 0x00007fe6525cf4ca in __GI___waitpid (pid=9486, stat_loc=stat_loc
entry=0x7fffbf55f930, options=options
entry=0) at ../sysdeps/unix/sysv/linux/waitpid.c:29
#1 0x00007fe652548fbb in do_system (line=<optimised out>) at ../sysdeps/posix/system.c:148
#2 0x00007fe6587259bd in TUnixSystem::StackTrace() () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#3 0x00007fe65863ff6a in DefaultErrorHandler(int, bool, char const*, char const*) () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#4 0x00007fe65863f902 in ErrorHandler () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#5 0x00007fe65863fdff in Fatal(char const*, char const*, ...) () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#6 0x00007fe6553db180 in operator delete(void*) () from /opt/cern/root/root_v6.10.02/lib/libNew.so
#7 0x00007fe65799410e in TF1::IntegralOneDim(double, double, double, double, double&) () from /opt/cern/root/root_v6.10.02/lib/libHist.so
#8 0x00007fe657987479 in TF1::Integral(double, double, double) () from /opt/cern/root/root_v6.10.02/lib/libHist.so
#9 0x0000000000402546 in main (argc=<optimised out>, argv=<optimised out>) at main.cpp:105
Aborted (core dumped)
Whereas when the output is correct, it gives the lines below where something looks wrong in the TF1::GetQuantiles function :
Curve length is : 9.42478 units
Probas:
0 0.125
1 0.25
2 0.375
3 0.5
4 0.625
5 0.75
6 0.875
7 1
Error in <TArrayD::operator[]>: index 10000 out of bounds (size: 10000, this: 0x7fffc91c0330)
Error in <TArrayD::operator[]>: index 10001 out of bounds (size: 10001, this: 0x7fffc91c0310)
Quantiles:
0 0.125 2.11529
1 0.25 2.45672
2 0.375 3.01118
3 0.5 3.75736
4 0.625 4.66658
5 0.75 5.7039
6 0.875 6.82946
7 1 8
ROfZ at quantiles:
0 2.11529 34.8295
1 2.45672 33.7039
2 3.01118 32.6666
3 3.75736 31.7574
4 4.66658 31.0112
5 5.7039 30.4567
6 6.82946 30.1153
7 8 30
Interval length:
1 2.45672 1.1781
2 3.01118 1.1781
3 3.75736 1.1781
4 4.66658 1.1781
5 5.7039 1.1781
6 6.82946 1.1781
7 8 1.1781
The simple and short program, compiled on Ubuntu 16.04 LTS using root 6.10/02 Built for linuxx8664gcc From tag v6-10-02, 6 July 2017 is the following :
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
#include <algorithm>
#include <cmath>
#include <ctype.h>
#include <float.h> // required for LDBL_EPSILON, DBL_MAX
#include <fstream>
#include <iomanip>
#include <iostream>
#include <math.h> // required for fabs(), fabsl(), sqrtl(), and M_PI_2
#include <sstream>
#include <stdexcept>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <utility> // std::pair, std::make_pair
#include <vector>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_ellint.h>
// ROOT (http://root.cern.ch) specific files
#include "TROOT.h"
#include "Riostream.h"
#include "TArrayD.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TF1.h"
#include "TGraph.h"
#include "TH1D.h"
#include "TNamed.h"
#include "TObject.h"
#include "TRint.h"
#include "TString.h"
#include "TMath.h"
#include "Math/IFunction.h"
#include "Math/SpecFuncMathMore.h"
#include "Math/WrappedTF1.h"
#include "Math/GSLIntegrator.h"
// Use functions from namespace std
using namespace std;
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
TROOT root("Rint", "The ROOT Interactive Interface");
int main(int argc, char **argv)
{
TRint theApp("Electrobem", NULL, NULL);
gROOT->Reset();
Int_t icell = 0;
Double_t zmin = 2;
Double_t zmax = 8;
const Int_t nz = 8;
Double_t dz[nz+1], PointAlongZ[nz+1];
TCanvas *cnv = new TCanvas("cnv", "cnv", 700., 700.);
TF1 * ROfZ = new TF1("ROfZ","36-sqrt(36-(x-8)*(x-8))", zmin, zmax);
ROfZ->SetNpx(2000);
ROfZ->SetLineColor(kRed);
ROfZ->Draw();
TF1 * SOfZ = new TF1("SOfZ","sqrt(1+(x-8)*(x-8)/(36-(x-8)*(x-8)))/[0]", zmin, zmax);
SOfZ->SetParameter(0, 1.);
SOfZ->SetNpx(10000);// useful for integral !!
Double_t CurveL = SOfZ->Integral(zmin, zmax, 1.e-07);
cout << "Curve length is : " << CurveL << " units" << endl;
SOfZ->SetParameter(0, CurveL); // PDF normalization
SOfZ->SetNpx(10000);
// Number of desired equidistant point along the curve
Double_t q[nz], pr[nz], foq[nz];
// Array of probabilities in which to divide Y axis before taking
// reciproqual of the distribution called quantiles
cout << "Probas:" << endl;
for (Int_t ipr = 0; ipr < nz; ++ipr){
pr[ipr] = (ipr + 1) / ((double)nz);
cout << ipr << " " << pr[ipr] << endl;
}
// Quantiles determination
SOfZ->GetQuantiles(nz, q, pr); // does not like the quantile for xq=1
q[nz-1] = zmax;
cout << "Quantiles:" << endl;
for (Int_t ipr = 0; ipr<nz; ipr++){
cout << ipr << " " << pr[ipr] << " " << q[ipr] << endl;
}
cout << "ROfZ at quantiles:" << endl;
for (Int_t ipr = 0; ipr < nz; ++ipr){
foq[ipr] = ROfZ->Eval(q[ipr]);
cout << ipr << " " << q[ipr] << " " << foq[ipr] << endl;
}
cout << "Interval length:" << endl;
SOfZ->SetParameter(0, 1.);
for (Int_t ipr = 1; ipr < nz; ++ipr){
cout << ipr << " " << q[ipr] << " " << SOfZ->Integral(q[ipr-1], q[ipr]) << endl;
}
TGraph *gr = new TGraph(nz, q, foq);
gr->SetLineColor(1);
gr->SetLineWidth(1);
gr->SetMarkerColor(1);
gr->SetMarkerSize(1);
gr->SetMarkerStyle(20);
gr->Draw("PLsame");
theApp.Run();
return 0;
}
I have not been able to find a bug on my side, i.e. in the above lines, therefore I am wondering if something is wrong in the functions TF1::Integral() and TF1::GetQuantiles().
Any hint would be helpful.
Thanks a lot,
Gilles