TH1F variable bin width

[sara.c]

I am trying to create a histogram with variable bin width, but I believe that the constructor for TH1F states that the widths must be binwidth+1.

Here are my lower bin bounds:

0.0, 0.4, 0.58, 0.90, 1.08, 1.50, 2.20.

Is it possible for me to plot these variable bins?

My input so far is:

root [0] Float_t Lower[7];
root [1] Lower[0]=0.0;
root [2] Lower[1]=0.40001;
root [3] Lower[2]=0.58;
root [4] Lower[3]=0.900001;
root [5] Lower[4]=1.08000001;
root [6] Lower[5]=1.5000001;
root [7] Lower[6]=2.2000001;
root [8] TH1F* h=new TH1F(“h”,“Gaussian PDF”,7,Lower);

I get the error “Error in TAxis::TAxis::Set: bins must be in increasing order,” which I believe is due to the fact that the widths must be binwidth+1 apart.

I am working with version 5.32.

Thanks in advance for any help!

[Wile_E_Coyote]

See the TH1:TH1 constructor description (in particular see what the expected size of the “xbins” array is).
Search for “variable bin width” in the ROOT User’s Guide -> Histograms chapter.

[sara.c]

Ah, yes I see the problem. This fixes it. Thank you!

[spyhunter0066]

Hi guys,

     I need to fill my histogram with a data where I don't have fix width. The data is basically consist of bin number, energy, counts as in the attachment. I tried to write two codes. One is using histogram, and another one is using graph. I couldn't managed to figure out how I should write the histogram part even if I checked so many things on web. There is one option by using the statement below.

TH1* h = new TH1D(
/* name / “h1”,
/ title / “Hist with variable bin width”,
/ number of bins / NBINS,
/ edge array */ edges
);

    The problem with this part is there is no way to  know my bin edges, the data gives me only the bin center which corresponds to my energy, and bin width is varying. 

However, it seems the graph works. I can fit Gaussian and 1st degree polynomial, but it doesn’t give me the right integral. I compared the real result with the MAESTRO or SCINTIVISION programs which ORTEC company gives with detectors. They don’t match at all.

       Could you help me to sort this issue by giving me little script about finding integrals for this kind of case where the binning changes?


     Sincerely yours

EnergyData.txt (92 Bytes)
can.cxx (14.7 KB)
sample.pdf (18.8 KB)
60Codata.txt (17.3 KB)

[couet]

I repeat what I told you by direct email:

[quote] don’t know my bin edges, the data gives me only the bin center which corresponds to my energy, and bin width is varying.
[/quote]
so if you know how the bin width is varying and you know the bin center, you can compute the bin edges and use variable bin size histogram ?

2. the example (sample.C) you posted is the C macro you obtained saving a TCanvas. This does not show the actual Fit methods you used to perform the fit. Can you send a small real example doing the Fit command ?

[Pepe_Le_Pew]

{
  const char *filename = "60Codata.txt";
  TGraph *g0 = new TGraph(filename, "%lg %lg %*lg"); // X axis = channels
  g0->SetTitle("g0;Channel number;Energy"); // energy calibration
  TGraph *g1 = new TGraph(filename, "%lg %*lg %lg"); // X axis = channels
  g1->SetTitle("g1;Channel number;Counts");
  TGraph *g2 = new TGraph(filename, "%*lg %lg %lg"); // X axis = energies
  g2->SetTitle("g2;Energy;Counts");
  if ((g1->GetN() < 2) || (g2->GetN() < 2)) return; // just a precaution
  
  TF1 *fg = new TF1("fg", "gaus");
  
  Int_t ch_min = 360, ch_max = 460; // "channel numbers"
  TF1 *f1 = new TF1("f1", "gaus(0) + pol1(3)", ch_min, ch_max);
  f1->SetParameters(30000, 430, 10, 5000, 0);
  g1->Fit(f1, "WEMR");
  fg->SetParameters(f1->GetParameters());
  std::cout << "gaus peak integral (mean -+ 5 sigma) = "
            << fg->Integral(f1->GetParameter(1) - 5.0 * f1->GetParameter(2),
                            f1->GetParameter(1) + 5.0 * f1->GetParameter(2))
            << std::endl;
  
  TF1 *f2 = new TF1("f2", "gaus(0) + pol1(3)",
                    g2->GetX()[ch_min], g2->GetX()[ch_max]);
  f2->SetParameters(30000, 1180, 25, 5000, 0);
  g2->Fit(f2, "WEMR");
  Double_t s, sa;
  fg->SetParameters(f2->GetParameters());
  s = fg->Integral(f2->GetParameter(1) - 5.0 * f2->GetParameter(2),
                   f2->GetParameter(1) + 5.0 * f2->GetParameter(2));
  // note: sa = s / "average X bin width"
  sa = s / ((g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min));
  std::cout << "gaus peak integral (mean -+ 5 sigma) = " << s
            << " ( adjusted = " << sa << " )"
            << std::endl;
  
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(1, 2);
  c->cd(1);
  g1->Draw("AP");
  c->cd(2);
  g2->Draw("AP");
  c->cd(0);
}

[spyhunter0066]

Hi,

   Can you remind me what GetParameters() is for ?  I checked, it returns to 

“virtual Double_t* TF1::GetParameters() const” . I mean why we use it or need it. I also check your version of integral function which is "Double_t TF1::Integral ( Double_t a,Double_t b,Double_t epsrel =1.e-12 ) " in ROOT website as a usage. I have no idea about what epsrel is , but it was saying somewhere “specified relative accuracy”. Is that correct? If yes, what does it mean?

  I only used the version of Integral(...) and IntegralError(....) in my attachment which is called "can.cxx". Sorry, but  I am not familiar the way you used it.   I need an explanation of the line where integral function appears in your code. Could you also put how you find integral error? 

When  I ran your code , I got an error saying what is below.

Error:

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/home/ic00025/Desktop/NumuneSpectrum/SourceData/./cernsample.cpp:48:167: error: cannot initialize a parameter of type ‘Double_t’ (aka ‘double’) with an rvalue
of type ‘Double_t *’ (aka ‘double *’)
…= "<< fg->Integral(f1->GetParameter(1) - 5.0 * f1->GetParameter(2), f1->GetParameter(1) + 5.0 * f1->GetParameter(2), f1->GetParameters() ) << endl;
^~~~~~~~~~~~~~~~~~~
/home/ic00025/root_install/include/TF1.h:395:63: note: passing argument to parameter ‘epsrel’ here
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12);
^
In file included from input_line_12:9:
/home/ic00025/Desktop/NumuneSpectrum/SourceData/./cernsample.cpp:55:118: error: cannot initialize a parameter of type ‘Double_t’ (aka ‘double’) with an rvalue
of type ‘Double_t *’ (aka ‘double *’)
s = fg->Integral(f2->GetParameter(1) - 5.0 * f2->GetParameter(2), f2->GetParameter(1) + 5.0 * f2->GetParameter(2), f2->GetParameters());
^~~~~~~~~~~~~~~~~~~
/home/ic00025/root_install/include/TF1.h:395:63: note: passing argument to parameter ‘epsrel’ here
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12);
^
Error in : Dictionary generation failed!
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
cernsample.cpp (2.15 KB)

[Pepe_Le_Pew]

It seems that TF1::Integral(Double_t a, Double_t b, const Double_t* params = 0, Double_t epsilon = 1e-12) is no longer available in ROOT 6, so I modified my macro and it doesn’t use it now.

[spyhunter0066]

Hi again,

        Thanks for the quick response. I am appreciated about that.  Could  you please add the "integral error part" into your code since I also need that part? 

          I have one question about what you did for integral correction in the line below.

////////////
sa = s / ((g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min));
cout << “gaus peak integral (mean -+ 5 sigma) = " << s << " ( adjusted = " << sa << " )”<< endl;
///////////////////////

         As we know, the energy calibration is just shifting our histogram to the right in this case. The real deal is about counts we see. And, the count doesn't change before and after calibration. Thus, what's the logic behind 

your integral correction here. Why you think they are not the same with the first integral in your first pad you drew
by using channel versus count.

        By the way. I was getting 2.16351e+06  for the count with my code before we seemed it's fixed in some way. However, programs like MAESTRO or SCINTIVISION was giving 716850+/-1588.  This depends surely how you decide you Gaussian limits, so I'll say the results are pretty close.

       Thanks a lot. I'll wait the last bit of your answer.

[Pepe_Le_Pew]

@spyhunter0066

Could you please add the “integral error part” into your code since I also need that part?

I don’t want to take all the fun away from you

@spyhunter0066

Thus, what’s the logic behind your integral correction here.

You’ll need to talk to a mathematician about what a “definite integral” means

BTW. Get familiar with the TSpectrum class and “Spectrum tutorials”.

[spyhunter0066]

Ok. I liked this figure, https://en.wikipedia.org/wiki/Integral#/media/File:Riemann_Integration_and_Darboux_Lower_Sums.gif , on your link. Btw, surely, it’s a definite integral. I’ll check TSpectrum. Physics is fun. That’s why we all are doing it .

[Pepe_Le_Pew]

Using the "60Codata.txt" file (attached in another post here), try this:

{
  const Bool_t fit_with_fixed_background = kFALSE; // kFALSE or kTRUE
  
  const char *filename = "60Codata.txt";
  TGraph *g0 = new TGraph(filename, "%lg %lg %*lg"); // X axis = channels
  g0->SetTitle("g0;Channel number;Energy"); // energy calibration
  TGraph *g1 = new TGraph(filename, "%lg %*lg %lg"); // X axis = channels
  g1->SetTitle("g1;Channel number;Counts");
  TGraph *g2 = new TGraph(filename, "%*lg %lg %lg"); // X axis = energies
  g2->SetTitle("g2;Energy;Counts");
  if ((g1->GetN() < 2) || (g2->GetN() < 2)) return; // just a precaution
  
  Int_t ch_min = 360, ch_max = 460; // "channel numbers"
  TF1 *f1 = new TF1("f1", "gausn(0) + pol1(3)", ch_min, ch_max);
  f1->SetParameters((30000. * 10. * 2.5), 430., 10., 5000., 0.);
  g1->Fit(f1, "WEMR");
  std::cout << "gaus peak integral = " << f1->GetParameter(0)
            << " +/- " << f1->GetParError(0)
            << std::endl;
  if (fit_with_fixed_background) {
    f1->FixParameter(3, f1->GetParameter(3));
    f1->FixParameter(4, f1->GetParameter(4));
    g1->Fit(f1, "WEMR");
    std::cout << "note: background parameters are fixed" << std::endl;
    std::cout << "gaus peak integral = " << f1->GetParameter(0)
              << " +/- " << f1->GetParError(0)
              << std::endl;
  }
  
  TF1 *f2 = new TF1("f2", "gausn(0) + pol1(3)",
                    g2->GetX()[ch_min], g2->GetX()[ch_max]);
  f2->SetParameters((30000. * 28. * 2.5), 1180., 28., 5000., 0.);
  g2->Fit(f2, "WEMR");
  std::cout << "gaus peak integral = " << f2->GetParameter(0)
            << " +/- " << f2->GetParError(0)
            << std::endl;
  Double_t dx; // "average X bin width"
  dx = (g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min);
  std::cout << "gaus peak integral adjusted = " << f2->GetParameter(0) / dx
            << " +/- " << f2->GetParError(0) / dx
            << std::endl;
  if (fit_with_fixed_background) {
    f2->FixParameter(3, f2->GetParameter(3));
    f2->FixParameter(4, f2->GetParameter(4));
    g2->Fit(f2, "WEMR");
    std::cout << "note: background parameters are fixed" << std::endl;
    std::cout << "gaus peak integral = " << f2->GetParameter(0)
              << " +/- " << f2->GetParError(0)
              << std::endl;
    std::cout << "gaus peak integral adjusted = " << f2->GetParameter(0) / dx
              << " +/- " << f2->GetParError(0) / dx
              << std::endl;
  }
  
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(1, 2);
  c->cd(1);
  g1->Draw("AP");
  c->cd(2);
  g2->Draw("AP");
  c->cd(0);
}

and / or that:

{
  const Bool_t use_partial_covariance = kFALSE; // kFALSE or kTRUE
  
  const char *filename = "60Codata.txt";
  TGraph *g0 = new TGraph(filename, "%lg %lg %*lg"); // X axis = channels
  g0->SetTitle("g0;Channel number;Energy"); // energy calibration
  TGraph *g1 = new TGraph(filename, "%lg %*lg %lg"); // X axis = channels
  g1->SetTitle("g1;Channel number;Counts");
  TGraph *g2 = new TGraph(filename, "%*lg %lg %lg"); // X axis = energies
  g2->SetTitle("g2;Energy;Counts");
  if ((g1->GetN() < 2) || (g2->GetN() < 2)) return; // just a precaution
  
  TF1 *fg = new TF1("fg", "gaus");
  Double_t sb = 5.0; // sigma band for the fg integral (e.g. 3 or 5 or 7)
  
  Int_t ch_min = 360, ch_max = 460; // "channel numbers"
  TF1 *f1 = new TF1("f1", "gaus(0) + pol1(3)", ch_min, ch_max);
  f1->SetParameters(30000., 430., 10., 5000., 0.);
  TFitResultPtr r1 = g1->Fit(f1, "WEMRS");
  TMatrixD c1 = r1->GetCovarianceMatrix();
  TMatrixD c1g = c1.GetSub(0, 2, 0, 2);
  Double_t s1, s1e;
  fg->SetParameters(f1->GetParameters());
  // fg->SetParErrors(f1->GetParErrors());
  s1 = fg->Integral(fg->GetParameter(1) - sb * fg->GetParameter(2),
                    fg->GetParameter(1) + sb * fg->GetParameter(2));
  s1e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                          fg->GetParameter(1) + sb * fg->GetParameter(2),
                          fg->GetParameters(),
                          c1g.GetMatrixArray());
  std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
            << s1 << " +/- " << s1e << std::endl;
  if (use_partial_covariance) {
    c1g -=
      c1.GetSub(0, 2, 3, 4) *
      c1.GetSub(3, 4, 3, 4).InvertFast() *
      c1.GetSub(3, 4, 0, 2);
    s1e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                            fg->GetParameter(1) + sb * fg->GetParameter(2),
                            fg->GetParameters(),
                            c1g.GetMatrixArray());
    std::cout << "note: using partial covariance" << std::endl;
    std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
              << s1 << " +/- " << s1e << std::endl;
  }
  
  TF1 *f2 = new TF1("f2", "gaus(0) + pol1(3)",
                    g2->GetX()[ch_min], g2->GetX()[ch_max]);
  f2->SetParameters(30000., 1180., 28., 5000., 0.);
  TFitResultPtr r2 = g2->Fit(f2, "WEMRS");
  TMatrixD c2 = r2->GetCovarianceMatrix();
  TMatrixD c2g = c2.GetSub(0, 2, 0, 2);
  Double_t s2, s2e;
  fg->SetParameters(f2->GetParameters());
  // fg->SetParErrors(f2->GetParErrors());
  s2 = fg->Integral(fg->GetParameter(1) - sb * fg->GetParameter(2),
                    fg->GetParameter(1) + sb * fg->GetParameter(2));
  s2e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                          fg->GetParameter(1) + sb * fg->GetParameter(2),
                          fg->GetParameters(),
                          c2g.GetMatrixArray());
  std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
            << s2 << " +/- " << s2e << std::endl;
  Double_t dx; // "average X bin width"
  dx = (g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min);
  std::cout << "gaus peak integral adjusted (mean +/- " << sb << " sigma) = "
            << s2 / dx << " +/- " << s2e / dx << std::endl;
  if (use_partial_covariance) {
    c2g -=
      c2.GetSub(0, 2, 3, 4) *
      c2.GetSub(3, 4, 3, 4).InvertFast() *
      c2.GetSub(3, 4, 0, 2);
    s2e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
                            fg->GetParameter(1) + sb * fg->GetParameter(2),
                            fg->GetParameters(),
                            c2g.GetMatrixArray());
    std::cout << "note: using partial covariance" << std::endl;
    std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
              << s2 << " +/- " << s2e << std::endl;
    std::cout << "gaus peak integral adjusted (mean +/- " << sb << " sigma) = "
              << s2 / dx << " +/- " << s2e / dx << std::endl;
  }
  
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(1, 2);
  c->cd(1);
  g1->Draw("AP");
  c->cd(2);
  g2->Draw("AP");
  c->cd(0);
}

See also: Wikipedia → Covariance matrix → Partial covariance matrix