Using an iterative algorithm to generate non-linear energy bins

Dear Forum,

I start with a linear distribion of energy bins from 0.1 to 15.4 MeV in 100 keV bins. My goal is to fit experimental data and determine the chi-square. My next step is to perturb an energy bin say [...., 10.1, 10.2, 10.35, 10.4, ...] where I made the bin width between 10.2 and 10.3 to be 10.35 or I changed this bin width to be 150 keV. Once I make a change to the energy distribution I will then fit my data and find a new chi-square. The idea is to either accept or reject the new distribution. I have written a test code to do this but I think something is incorrect and I cannot locate my error. If someone has an idea on what my issue that would be very helpful. I tried to cover all the possible scenarios but I must be missing something. Here is my code:

function test_rnd_xbin

dNdE_energy = [1:0.1:15.4]';
dNdE = ones(length(dNdE_energy),1);
energy_binning = [-.200 -.150 -.100 -.050 .050 .100 .150 .200];

for i=1:200
new_val_energy_bin = randperm(length(energy_binning),1);
new_val_energy_adj = energy_binning(new_val_energy_bin)
new_val_energy = randperm(length(dNdE_energy),1);
energy_bin_region = dNdE_energy(new_val_energy);

if new_val_energy_adj > 0 & new_val_energy < length(dNdE_energy)
    if (new_val_energy_adj + energy_bin_region) > dNdE_energy(new_val_energy+1)
        energy_plus = find(dNdE_energy > round((new_val_energy_adj + energy_bin_region),3))
        dNdE_energy = [dNdE_energy(1:new_val_energy)' dNdE_energy(new_val_energy)+new_val_energy_adj' dNdE_energy(energy_plus)']'     
        
    elseif (new_val_energy_adj + energy_bin_region) == dNdE_energy(new_val_energy+1)
        
    else (new_val_energy_adj + energy_bin_region) < dNdE_energy(new_val_energy+1)
        dNdE_energy = [dNdE_energy(1:new_val_energy)' dNdE_energy(new_val_energy)+new_val_energy_adj' dNdE_energy(new_val_energy+1:end)']'  
    end
else new_val_energy_adj < 0 & new_val_energy > 1
    if (new_val_energy_adj + energy_bin_region) < dNdE_energy(new_val_energy-1)
        energy_minus = find(dNdE_energy < round((new_val_energy_adj + energy_bin_region),3))
        dNdE_energy = [dNdE_energy(energy_minus)' dNdE_energy(new_val_energy)+new_val_energy_adj' dNdE_energy(new_val_energy:end)']'   
        
    elseif (new_val_energy_adj + energy_bin_region) == dNdE_energy(new_val_energy-1)
        
    else (new_val_energy_adj + energy_bin_region) < dNdE_energy(new_val_energy-1)
        dNdE_energy = [dNdE_energy(1:new_val_energy-1)' dNdE_energy(new_val_energy)+new_val_energy_adj' dNdE_energy(new_val_energy:end)']'  
    end
end
dNdE = ones(length(dNdE_energy),1);
bar(dNdE_energy,dNdE)
pause(.05)
end

[David Goodmanson]

HI Chad,

Some of the values in energy_binning are considerably larger than the inital 0.1 bin width. Is the idea to iteratively vary the width of the bins and keep the total number of bins the same, and have the total domain still between 1 and 15.4? Or can the number of bins and/or the domain vary?

[Chad]

Dear David,

The number of bins can vary and in fact I have shown that a predefined non-linear binning works well. For example, I have this in the past

dNdE_energy = 1 + (delta_E*num_bins)^(3/2); Where delta_E was 0.05 MeV and the numbers of bins was ~ 200

I think the problem is that my logic is somwhat correct but I must missing a scenario that is causing two bins to have the same energy.