Empty bins after calibration

Newbie
1

Greetings,

I am having difficulties calibrating an energy spectrum.

I have a vector ENERGY which I fill an histogram with. The result looks fine (uncal.pdf).

However, after I apply the calibration
ENERGY_CAL[i] = channel2Energy(ENERGY[i]),
where channel2Energy is a linear function, I fill another histogram with ENERGY_CAL and get a weird looking result (cal.pdf).

I already tried using a different binning for the second histogram but I cannot get it to look nice.

Thanks for your time!
Christoph

uncal.pdf (77.4 KB)
cal.pdf (92.8 KB)

2

Hi Christoph, it is because you pass from integer values (channel) to continuous values (energy). To smooth your spectrum, in your function channel2Energy, you should add a random number (between 0 and 1) to the channel value.
Doing so smooths your calibrated spectrum because when you have channel with value 10, the “exact” value is between 10 and 11. Adding a random value solves this “issue”.

3

Maybe it will be sufficient to try:

Double_t dE = hSpectrumCalibrated->GetBinWidth(1) / 2.0; // half of the "fixed" bin width
// ...
hSpectrumCalibrated->Fill(ENERGY_CAL[i] + dE);

However, in principle, you should create a histogram with variable bin sizes.

See also:

1 Like
4

Actually ENERGY does not have integer values. The function channel2Energy is named very badly, I know. It was used to convert channels once, that’s why it is still called like that.
Some examples:

i=0, ENERGY=752.125345, ENERGY_CAL=173.638841
i=10000, ENERGY=11236.404687, ENERGY_CAL=2612.239957
i=20000, ENERGY=9045.618917, ENERGY_CAL=2102.608610
i=30000, ENERGY=281.856741, ENERGY_CAL=64.083220
i=40000, ENERGY=9664.323836, ENERGY_CAL=2246.589351
i=50000, ENERGY=2040.249410, ENERGY_CAL=473.230432
i=60000, ENERGY=3721.239336, ENERGY_CAL=864.234674
i=70000, ENERGY=2446.173003, ENERGY_CAL=567.666911
i=80000, ENERGY=3925.829752, ENERGY_CAL=911.685516
i=90000, ENERGY=11248.942738, ENERGY_CAL=2615.031183
i=100000, ENERGY=3133.295953, ENERGY_CAL=727.464600
5

Ok, but this effect comes clearly from the fact you started with integer values that you converted to continuous values. So, if you do not have the channel values anymore, maybe you can try what @Wile_E_Coyote has proposed.

6

The easy solution with dE did not work for me.

I looked at your link and tried out the scaling for histograms with fixed size. Explicitly I did

TH1F* hSpectrumCalibrated = (TH1F*) hSpectrumUncalibrated->Clone("hSpectrumCalibrated");
auto a =  hSpectrumCalibrated->GetXaxis();
a->Set( a->GetNbins(), 0, channel2Energy(fCalibration, a->GetXmax()) );

which produces exactly the spectrum I want! very nice.

Sadly this does break something further down in my code, what could I have changed?
Why should I create a histogram with variable bin size? At least something seems to have been fixed with this approach.

7

In your current approach, you should use (instead of “0”):
channel2Energy(fCalibration, a->GetXmin())
and in the end you need to add:
hSpectrumCalibrated->ResetStats();

If at some point your “channel2Energy” function is no longer linear then the energy bins’ widths will not be constant / equal (so you will need a histogram with variable bin sizes).

BTW. If “+ dE” didn’t work well then you might try “- dE” (or e.g. “+ dE / 2.0” or “- dE / 2.0”).

8

Thank you very much!
The problem has already been solved by your previous answer. The crash was my fault.

However, now the number of entries has been reduced even though I did ResetStats() and everything looks fine now.

Something else also still bothers me. I don’t understand why it should make a difference if I calculate ENERGY_CAL and fill a histogram with it, or if I simply scale an axis if I use the same calibration/scaling function in both cases. I have constructed an example creating a histogram hSpectrumCalibrated2 which displays the difference in both cases (see attached pdfs):

hSpectrumUncalibrated = new TH1F("hSpectrumUncalibrated", "Th-228 -- uncalibrated", nEnergyBins, 0, nEnergyBins-1);
// ...
TH1F *  hSpectrumCalibrated = new TH1F("hSpectrumCalibrated", "Th-228 -- calibrated", nEnergyBins, channel2Energy(fCalibration, 0), channel2Energy(fCalibration, nEnergyBins-1));
  for(Int_t i=0; i<nEVENTS; i++) {
    ENERGY_CAL[i] = channel2Energy(fCalibration, ENERGY[i]);
    hSpectrumCalibrated->Fill(ENERGY_CAL[i]);
  }

TH1F* hSpectrumCalibrated2 = (TH1F*) hSpectrumUncalibrated->Clone("hSpectrumCalibrated2");
hSpectrumCalibrated2->SetTitle("Th-228 calibrated2");
hSpectrumCalibrated2->SetXTitle("energy (keV)");
hSpectrumCalibrated2->SetYTitle("counts");
auto a =  hSpectrumCalibrated2->GetXaxis();
a->Set( a->GetNbins(), channel2Energy(fCalibration, a->GetXmin()), channel2Energy(fCalibration, a->GetXmax()) );
hSpectrumCalibrated2->ResetStats();

hSpectrumCalibrated2->Add(hSpectrumCalibrated, -1);

It seems like the difference is highest at the peaks. Can someone elaborate on that and explain which of the two approaches is the correct one to use?

NoE shows that the number of entries has changed.

calibrated2.pdf (90.2 KB) calibrated2_zoom.pdf (16.2 KB)NoE.pdf (77.2 KB)

9

TH1::ResetStats

10

Well, nothing unexpected there.

But out of curiosity I just did ResetStats() on hSpectrumUncalibrated and now it also shows me the same number as for hSpectrumCalibrated.

Therfore after filling hSpectrumUncalibrated has somehow stored a wrong number of entries. How come?

11

Try to print the underflow:
hSpectrumUncalibrated->GetBinContent(0)
and the overflow:
hSpectrumUncalibrated->GetBinContent(nEnergyBins+1)

BTW. I think you want the upper edge “nEnergyBins” and NOT “nEnergyBins - 1”.

12

The under- and overflow don’t suffice to explain the difference:

root [1] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetEntries()
(double) 8665377.000000
root [2] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->ResetStats()
root [3] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetEntries()
(double) 5543289.000000
root [4] 8665377.000000-5543289.000000
(double) 3122088.000000
root [5] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetBinContent(0)
(double) 0.000000
root [6] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetBinContent(12000+1)
(double) 42154.000000
// ...
root [15] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetBinContent(12000)
(double) 10.000000

BWT: Yeah, you are right.

13

Make sure that you did not use: TH1::SetAxisRange
Try to print: TH1::GetEffectiveEntries
Try to print: TH1::Integral

14

There is no string containing “SetRange” or “SetAxis” in my code above the relevant lines.

Here are the print-outs:

root [1] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetEntries()
(double) 8665377.000000
root [2] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetEffectiveEntries()
(double) 5543289.000000
root [3] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->Integral()
(double) 5543289.000000
root [4] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->ResetStats()
root [5] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetEntries()
(double) 5543289.000000
root [6] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->GetEffectiveEntries()
(double) 5543289.000000
root [7] ((TH1F*)E_cal->Get("hSpectrumUncalibrated"))->Integral()
(double) 5543289.000000

Can you infer anything from that?

15

I guess you fill the “hSpectrumUncalibrated” with weights not equal to 1.

16

I simply to the following for every event:
hSpectrumUncalibrated->Fill(energy);

17

Maybe @couet has some another idea.

In the meantime, print (just to make sure that it is less than about 1.6777e7):
hSpectrumUncalibrated->GetBinContent(hSpectrumUncalibrated->GetMaximumBin())

18

Do you have some code (we could run) reproducing the problem ?

19

The one generating the histograms I guess?
I will try to produce an example, but this will take a while.

20

Yes if you can that might be useful.
Or the histogram saved in a ROOT file ?