Home | News | Documentation | Download

Consistency between a histogram and a flux (/cm^2/MeV/s) plot?

Hello everyone:

I am comparing the cosmic ray flux generated by CRY (https://nuclear.llnl.gov/simulation/main.html) and would like to compared that with some existing known flux plots e.g. \Phi /cm^2/s and/or e.g. \Phi /cm^2/MeV/s. This is, therefore, more an exercise to reproduce the plots, rather than a discovery. After I generate, say muons, I filled a histogram with uniform bins (between 10^{-5} MeV to 10^{5} MeV) embedded in a canvas set with log scale in x and y directions:

h=ROOT.TH1F(particle+'_'+altitude+' m_histogm',particle+'_'+altitude+' m',10000000,0.0,100.0)

However, the plots do not match with the expected flux as seen from EXPACS page: https://phits.jaea.go.jp/expacs/. Nice tables of energy (MeV) vs flux (/cm^2/MeV/s) or
energy (MeV) vs flux (/cm^2/s) are available and one can easily get a TGraph to see the consistency with the histogram-view.

Later, I tried do change to non-uniform bins:

        nbins = 100;
        xmin  = 1e-5;
        xmax  = 1e+5;
        logxmin = ROOT.TMath.log10(xmin);
        logxmax = ROOT.TMath.log10(xmax);
        binwidth = (logxmax-logxmin)/nbins;
        xbins = [None]*(nbins+1)
        xbins[0] = xmin

        for i in range(nbins+1):
                xbins[i] = xmin + ROOT.TMath.Power(10,logxmin+i*binwidth)
                print(str(xbins[i])+' = '+str(xmin)+' + 10^{'+str(logxmin)+'+'+str(i)+'*'+str(binwidth)+'}, '+str(ROOT.TMath.Power(10.0,logxmin+i*binwidth))+', type '+str(type(xbins[i])))

        #print(str(xbins))

        h=ROOT.TH1F(particle+'_'+altitude+' m_histogm',particle+'_'+altitude+' m',nbins,array('f',xbins))

This one seems to be more close to EXPACS plots, but I am still not convinced, as not all the plots for all the particles show similar behavior. In this case also, I cannot understand if the bin-widths are chosen wisely. Should they all be 1 MeV?

There seems to be a couple of issues going on here:

(a) The first part is: I cannot understand how to compare the frequency of the TH1F with the plot for which the well defined units are available.

(b) There are other normalization related issues, that is specific to CRY simulation - but these are minor to me, at this moment.

Can anyone please throw light on how to check this consistency? I understand that this must be a simple problem - just I cannot understand how to do it.

Any help will be appreciated,
Regards,
Kolahal

It seems to me that your questions are not really related to ROOT. I think, you should ask these questions to your physics supervisor.

mmm…yes…I agree…

Perhaps I posted this question here, since this is a subform on general discussion in ROOT and I thought someone might show where in my code I was going wrong. But I agree, this is more a confusion about the concept and I apologize. This post can be removed, if necessary.

Regards,
Kolahal

Hello @kolahalb,

You may indeed ask a colleague/supervisor, but in theory you should be able to use all kinds of bin widths. However, if you use non-equidistant bins, just remember that you have to decide whether you want to divide by the bin width/volume, i.e. plot a density, or you don’t want that.
For everything else, we are probably not the experts.

Hello @StephanH and @Wile_E_Coyote: yes. I guess my issue is happening in normalizing w.r.t. energy and I will figure that out. Thanks.

Regards,
Kolahal