ROOT::Interpolation: Busy flag cleared

[Deb]

Hello everyone,

I’m trying to follow the example of Interpolation function that you can find here: https://root.cern.ch/function-interpolation.
My code is:

#include <stdint.h>
#include <stdio.h>
#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
#include <vector>
#include <iostream>
#include "Math/Polynomial.h"
#include "Math/Interpolator.h"
#include "TCanvas.h"
#include "TGraph.h"

using namespace std;

void Interpolation(){
  
  TH1D *h;
  int xmin=0, xmax=200;
  double X_max=0, Y_max=0;
  TFile *f=new TFile("run_4/result.root" ,"READ");
  
  const Int_t ni = xmax-xmin/6.;
  Double_t xi[ni], yi[ni];
  
  const Int_t nb = ni*10;
  Double_t x[nb], y[nb], iy[nb];
  
  f->cd();
  h=(TH1D*)gDirectory->Get("h1x");
  int n=0;

  for(int i_bin=1; i_bin<200; i_bin++){
    if(((TH1D*)gDirectory->Get("h1x"))->GetBinContent(i_bin)!=0 && n==0){
      xmin=i_bin;
      n++;
    }
    if(((TH1D*)gDirectory->Get("h1x"))->GetBinContent(i_bin)==0 && n>0){
      xmax=i_bin;
      i_bin=200;
    }
  }
  
  ROOT::Math::Polynomial polyf(2);
  double p[3];
  X_max=(((TH1D*)gDirectory->Get("h1x"))->GetMaximumBin());
  Y_max=((TH1D*)gDirectory->Get("h1x"))->GetBinContent(((TH1D*)gDirectory->Get("h1x"))->GetMaximumBin());
  
  p[1]=Y_max;
  p[2]=X_max*p[1]*2;
  p[0]=p[1]-p[1]*pow(X_max,2)-p[2]*X_max;
  
  polyf.SetParameters(p);
  ROOT::Math::IBaseFunctionOneDim & f1 = polyf;
  ROOT::Math::Interpolator inter(ni,  ROOT::Math::Interpolation::kPOLYNOMIAL);
  
  for ( Int_t i = 0; i < ni; ++i )
   {   
      xi[i] = (Double_t) i*(xmax-xmin)/(ni-1) + xmin;
      yi[i] = f1(xi[i]);
   }
   inter.SetData(ni, xi, yi);
 
   for ( Int_t i = 0; i < nb; ++i )
   {
      x[i]  = (Double_t) i*(xmax-xmin)/(nb-1) + xmin;
      y[i]  = f1(x[i]);
      iy[i] = inter.Eval(x[i]);
   }
  
   new TCanvas("c1", "Two Graphs", xmin, xmax); 
   
   TGraph* gf = new TGraph(ni, xi, yi);
   gf->SetLineColor(14);
   gf->SetLineWidth(3);
   gf->SetTitle("Function: sin(x)");
   gf->SetMarkerStyle(22);
   gf->Draw("AP");
   
   TGraph* gi = new TGraph(nb, x, iy);
   gi->SetLineColor(28);
   gi->SetLineWidth(1);
   gi->SetTitle("Integral");
   gi->Draw("SAME L"); 

   xi.clear();
   yi.clear();
   x.clear();
   y.clear();
   iy.clear();
   
}

The problem is that it breaks immediately, giving me this:

root [0] .L Interpolation.C
root [1] Interpolation()

 *** Break *** segmentation violation



===========================================================
There was a crash.
This is the entire stack trace of all threads:
===========================================================
#0  0x00007faf913e307a in __GI___waitpid (pid=6305, stat_loc=stat_loc
entry=0x7fff94e7d500, options=options
entry=0) at ../sysdeps/unix/sysv/linux/waitpid.c:29
#1  0x00007faf9135bfbb in do_system (line=<optimized out>) at ../sysdeps/posix/system.c:148
#2  0x00007faf920bba04 in TUnixSystem::Exec (shellcmd=<optimized out>, this=0xaa5520) at /home/debora/Documents/root/core/unix/src/TUnixSystem.cxx:2172
#3  TUnixSystem::StackTrace (this=0xaa5520) at /home/debora/Documents/root/core/unix/src/TUnixSystem.cxx:2419
#4  0x00007faf920bdcac in TUnixSystem::DispatchSignals (this=0xaa5520, sig=kSigSegmentationViolation) at /home/debora/Documents/root/core/unix/src/TUnixSystem.cxx:1294
#5  <signal handler called>
#6  __memmove_avx_unaligned () at ../sysdeps/x86_64/multiarch/memcpy-avx-unaligned.S:126
#7  0x00007faf8f22df03 in ROOT::Math::IBaseFunctionOneDim::operator() (x=23.010050251256281, this=<optimized out>) at /home/debora/Documents/root-build/include/Math/IFunction.h:154
#8  G__G__MathCore_222_0_3 (result7=result7
entry=0xad25f0, funcname=funcname
entry=0xffffffffffffffff <error: Cannot access memory at address 0xffffffffffffffff>, libp=libp
entry=0x7fff94e80760, hash=<optimized out>) at /home/debora/Documents/root-build/math/mathcore/G__MathCore.cxx:11666
#9  0x00007faf90acb932 in Cint::G__ExceptionWrapper (funcp=funcp
entry=0x7faf8f22dec0 <G__G__MathCore_222_0_3(G__value*, char const*, G__param*, int)>, result7=result7
entry=0xad25f0, funcname=funcname
entry=0xffffffffffffffff <error: Cannot access memory at address 0xffffffffffffffff>, libp=libp
entry=0x7fff94e80760, hash=<optimized out>) at /home/debora/Documents/root/cint/cint/src/Api.cxx:393
#10 0x00007faf90ab997c in G__exec_asm (start=<optimized out>, stack=stack
entry=0, presult=presult
entry=0x7fff94e8b290, localmem=localmem
entry=0) at /home/debora/Documents/root/cint/cint/src/bc_exec_asm.h:916
#11 0x00007faf90a868ed in G__exec_loop (condition=0xae68f0 "i<ni", foraction=..., forinit=0x0) at /home/debora/Documents/root/cint/cint/src/parse.cxx:3318
#12 0x00007faf90a8305b in G__exec_for () at /home/debora/Documents/root/cint/cint/src/parse.cxx:2764
#13 G__exec_statement (mparen=mparen
entry=0x7fff94e8c35c) at /home/debora/Documents/root/cint/cint/src/parse.cxx:7262
#14 0x00007faf90ae81c0 in G__interpret_func (result7=result7
entry=0x7fff94e8cf50, funcname=<optimized out>, libp=libp
entry=0x7fff94e8cf90, hash=<optimized out>, p_ifunc=0x14a7e10, funcmatch=funcmatch
entry=1, memfunc_flag=0) at /home/debora/Documents/root/cint/cint/src/ifunc.cxx:6710
#15 0x00007faf90a06b90 in G__getfunction (item=item
entry=0xabaee0 "Interpolation()", known3=known3
entry=0x7fff94e97a44, memfunc_flag=<optimized out>, memfunc_flag
entry=0) at /home/debora/Documents/root/cint/cint/src/func.cxx:2901
#16 0x00007faf90a64e76 in G__getitem (item=0xabaee0 "Interpolation()") at /home/debora/Documents/root/cint/cint/src/expr.cxx:1918
#17 0x00007faf90a6b2b0 in G__getexpr (expression=<optimized out>) at /home/debora/Documents/root/cint/cint/src/expr.cxx:1488
#18 0x00007faf90a7f2ec in G__exec_function (presult=0x7fff94e9a4d0, plargestep=0x7fff94e9a34c, piout=0x7fff94e9a348, pc=<synthetic pointer>, statement=...) at /home/debora/Documents/root/cint/cint/src/parse.cxx:645
#19 G__exec_statement (mparen=mparen
entry=0x7fff94e9b05c) at /home/debora/Documents/root/cint/cint/src/parse.cxx:7371
#20 0x00007faf909f7914 in G__exec_tempfile_core (file=file
entry=0x0, fp=<optimized out>) at /home/debora/Documents/root/cint/cint/src/debug.cxx:266
#21 0x00007faf909f8e42 in G__exec_tempfile_fp (fp=<optimized out>) at /home/debora/Documents/root/cint/cint/src/debug.cxx:807
#22 0x00007faf90a91ab4 in G__process_cmd (line=line
entry=0x1049eb0 "Interpolation()", prompt=prompt
entry=0xaabc68 "", more=more
entry=0xaabc60, err=err
entry=0x7fff94e9c26c, rslt=rslt
entry=0x7fff94e9c270) at /home/debora/Documents/root/cint/cint/src/pause.cxx:3341
#23 0x00007faf920864c0 in TCint::ProcessLine (this=0xaabc20, line=0x1049eb0 "Interpolation()", error=0x0) at /home/debora/Documents/root/core/meta/src/TCint.cxx:549
#24 0x00007faf921370ef in TApplication::ProcessLine (this=0xc37e60, line=<optimized out>, sync=<optimized out>, err=0x0) at /home/debora/Documents/root/core/base/src/TApplication.cxx:978
#25 0x00007faf925d758b in TRint::HandleTermInput (this=0xc37e60) at /home/debora/Documents/root/core/rint/src/TRint.cxx:583
#26 0x00007faf920bd37c in TUnixSystem::CheckDescriptors (this=this
entry=0xaa5520) at /home/debora/Documents/root/core/unix/src/TUnixSystem.cxx:1396
#27 0x00007faf920bdfda in TUnixSystem::DispatchOneEvent (this=0xaa5520, pendingOnly=<optimized out>) at /home/debora/Documents/root/core/unix/src/TUnixSystem.cxx:1103
#28 0x00007faf92106034 in TSystem::InnerLoop (this=0xaa5520) at /home/debora/Documents/root/core/base/src/TSystem.cxx:410
#29 0x00007faf92104acf in TSystem::Run (this=0xaa5520) at /home/debora/Documents/root/core/base/src/TSystem.cxx:360
#30 0x00007faf92134a8f in TApplication::Run (this=this
entry=0xc37e60, retrn=retrn
entry=false) at /home/debora/Documents/root/core/base/src/TApplication.cxx:1126
#31 0x00007faf925d86cb in TRint::Run (this=this
entry=0xc37e60, retrn=retrn
entry=false) at /home/debora/Documents/root/core/rint/src/TRint.cxx:454
#32 0x000000000040101c in main (argc=1, argv=0x7fff94e9e828) at /home/debora/Documents/root/main/src/rmain.cxx:29
===========================================================


The lines below might hint at the cause of the crash.
If they do not help you then please submit a bug report at
http://root.cern.ch/bugs. Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#6  __memmove_avx_unaligned () at ../sysdeps/x86_64/multiarch/memcpy-avx-unaligned.S:126
#7  0x00007faf8f22df03 in ROOT::Math::IBaseFunctionOneDim::operator() (x=23.010050251256281, this=<optimized out>) at /home/debora/Documents/root-build/include/Math/IFunction.h:154
#8  G__G__MathCore_222_0_3 (result7=result7
entry=0xad25f0, funcname=funcname
entry=0xffffffffffffffff <error: Cannot access memory at address 0xffffffffffffffff>, libp=libp
entry=0x7fff94e80760, hash=<optimized out>) at /home/debora/Documents/root-build/math/mathcore/G__MathCore.cxx:11666
===========================================================


Root > Function Interpolation() busy flag cleared

Any ideas, hints? Thank you very much in advance for all the help,

Ciao,

Debora

[ksmith]

Just FYI:

Getting the same histograms after cuts

Remark on the way to post code snippets: When you post pieces of code on the forum enclose it with ``` . That makes your post much more readable (I systematically edit your posts that way) . So do: ``` // this is a C++ piece of code. ``` And it will display as: // this is a C++ piece of code.

[ksmith]

You at least need some error checking:

...
  TFile *f=new TFile("run_4/result.root" ,"READ"); //The file might not have opened correctly.
  if (!f->IsOpen()) {
    cerr<<"ERROR: The file didn't open!\n";
    return -1;
  }
  
  const Int_t ni = xmax-xmin/6.;
  Double_t xi[ni], yi[ni];
  
  const Int_t nb = ni*10;
  Double_t x[nb], y[nb], iy[nb];
  
  f->cd();
  h=(TH1D*)gDirectory->Get("h1x"); //This could be NULL
  if (!h) {
    cerr<<"ERROR: The hist wasn't found!\n";
    return -1;
  }
  int n=0;

  for(int i_bin=1; i_bin<200; i_bin++){
    if(h->GetBinContent(i_bin)!=0 && n==0){ //Use h since you error checked the pointer above.
      xmin=i_bin;
      n++;
    }
...

Also, it is common practice to only include what you need. This reduces code bloat, compile times, etc. At least one of the following are not being used:

#include <stdint.h>
#include <stdio.h>
#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
#include <vector>
#include <iostream>
#include "Math/Polynomial.h"
#include "Math/Interpolator.h"
#include "TCanvas.h"
#include "TGraph.h"

[Deb]

Thanks very much for your help, It breaks anyway, but I guess it’s more complete now, always good to know how to better code. I will look into the error hint,thanks!

[Deb]

I found the line where I have a problem, and I corrected some more errors, just in case someone is interested in helping me, the line is
yi[i] = f1(xi[i]);
And the corrected code is:

#include <stdio.h>
#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
#include <iostream>
#include "Math/Polynomial.h"
#include "Math/Interpolator.h"
#include "TCanvas.h"
#include "TGraph.h"

using namespace std;

void Interpolation(){
  
  TH1D *h;
  int binmin=0, binmax=200;
  double X_max=0, Y_max=0;
  double xmax=200,xmin=-200;
  
  TFile *f=new TFile("run_4/result.root" ,"READ");
  if (!f->IsOpen()) {
    cerr<<"ERROR: The file didn't open!\n";
    return -1;
  }
  
  f->cd();
  h=(TH1D*)gDirectory->Get("h1x");
  int n=0;
  if (!h) {
    cerr<<"ERROR: The hist wasn't found!\n";
    return -1;
  }
  
  for(int i_bin=1; i_bin<200; i_bin++){
    if(h->GetBinContent(i_bin)!=0 && n==0){
      binmin=i_bin;
      n++;
    }
    if(h->GetBinContent(i_bin)==0 && h->GetBinContent(i_bin+5)==0 && n>0){
      binmax=i_bin;
      i_bin=200;
    }
  }
  
  xmin=(h->GetXaxis()->GetBinCenter(binmin))+200;
  xmax=(h->GetXaxis()->GetBinCenter(binmax))+200;
  cout<<"binmin and xmin: "<<xmin<<" "<<binmin<<" binmax and xmax: "<<xmax<<" "<<binmax<<endl;
  X_max=(h->GetXaxis()->GetBinCenter(h->GetMaximumBin()))+200;
  Y_max=h->GetBinContent(h->GetMaximumBin());
  cout<<"Xmax: "<<X_max<<" Ymax: "<<Y_max<<endl;

  const Int_t ni = xmax-xmin;
  Double_t xi[ni], yi[ni];
  
  const Int_t nb = ni*10;
  Double_t x[nb], y[nb], iy[nb];

  double p[3];
  p[0]=Y_max;
  p[1]=X_max*p[1]*2;
  p[2]=p[1]-p[1]*pow(X_max,2)-p[2]*X_max;
  
  ROOT::Math::Polynomial polyf(3);
  polyf.SetParameters(p);
  ROOT::Math::IBaseFunctionOneDim & f1 = polyf;
  ROOT::Math::Interpolator inter(ni,  ROOT::Math::Interpolation::kPOLYNOMIAL);

  for ( Int_t i = 0; i < ni; ++i )
   {   
     xi[i] = (Double_t) i*(xmax-xmin)/(ni-1) + xmin;
     yi[i] = f1(xi[i]);
   }
  inter.SetData(ni, xi, yi);
 
   for ( Int_t i = 0; i < nb; ++i )
   {
     x[i]  = (Double_t) i*(xmax-xmin)/(nb-1) + xmin;
     y[i]  = f1(x[i]);
     iy[i] = inter.Eval(x[i]);
   }
  
   new TCanvas("c1", "A Graph", xmin, xmax); 
   
   TGraph* gf = new TGraph(ni, xi, yi);
   gf->SetLineColor(14);
   gf->SetLineWidth(3);
   gf->SetTitle("Function polynomial ^2");
   gf->SetMarkerStyle(22);
   gf->Draw("AP");
   
   TGraph* gi = new TGraph(nb, x, iy);
   gi->SetLineColor(28);
   gi->SetLineWidth(1);
   gi->SetTitle("Integral");
   gi->Draw("SAME L"); 

   xi.clear();
   yi.clear();
   x.clear();
   y.clear();
   iy.clear();
   
}```

[ksmith]

@Deb

for(int i_bin=1; i_bin<200; i_bin++){
    if(h->GetBinContent(i_bin)!=0 && n==0){
      binmin=i_bin;
      n++;
    }
    if(h->GetBinContent(i_bin)==0 && h->GetBinContent(i_bin+5)==0 && n>0){
      binmax=i_bin;
      i_bin=200;
    }
  }

What about this piece of code? You don’t check that the histogram actually has that number of bins, and the i_bin+5 is bound to be problematic. Maybe this for loop works better for you?

for(int i_bin=1; i_bin<200 && i_bin<=h->GetNbinsX()-5; i_bin++){

[ksmith]

@ksmith

You don’t check that the histogram actually has that number of bins, and the i_bin+5 is bound to be problematic. Maybe this for loop works better for you?

for(int i_bin=1; i_bin<200 && i_bin<=h->GetNbinsX()-5; i_bin++){

After the above correction and removing the following block I was able to run your script without error.

   //The following is not valid for arrays:
/*
   xi.clear();
   yi.clear();
   x.clear();
   y.clear();
   iy.clear();
*/

I also replaced your histogram with a random gaussian as I don’t have your input:

   TH1F *h = new TH1F("h","",200,-2,2);
   h->FillRandom("gaus",5000);

It generated this image:

[Deb]

Thank you very much, I knew the number of bins in the histogram but I didn’t think about the -5, thanks.
Except it doesn’t work yet, at this point I think it’s something in my rootfile, something that has to do with my histos I guess. Anyway thanks again, have a nice day.

Debora