Total number of particles from ParticleBinning

[mlyto]

Hello everyone,

Could you help me with a simple question I'm still unclear on?
One of the diagnostics I use looks like this

DiagParticleBinning(                        # 11                                                                      
    deposited_quantity = "weight",
    every = [every_in, every_step],
    species = ["electron"],
    axes = [ ["ekin",   0.48,    10.,   100],
             ["pz",     -12.,    12.,     2],
             ["z",  z_foc - 30 * l0,  z_foc + 30 * l0,  3] ]    
)

Let's say I want to count the number of electrons at a given time (it = 13) by integrating over all energies, and only for the region pz > 0, z > z_foc + 10 * l0.
Here's how I do it.

S = happi.Open("results", show=True)
V_r = N_r * L_r**3   # it's 3D simulations
Diag_ekin_integ = S.ParticleBinning(11)
ekin_axis = Diag_ekin_integ.getAxis('ekin')
pz_axis = Diag_ekin_integ.getAxis('pz')
z_axis = Diag_ekin_integ.getAxis('z')
dE = ekin_axis[1] - ekin_axis[0]
dpz = pz_axis[1] - pz_axis[0]
dz = z_axis[1] - z_axis[0]
tot_ekin_spec_13 = Diag_ekin_integ.getData()[13]
ekin_forward_13 = tot_ekin_spec_13[:,1,2]   # pz>0 , z > z_foc + 10*l0
ekin_all_everywhere_13 = np.sum(np.sum(tot_ekin_spec_13, axis=1),axis=1)
N_forward = np.sum(ekin_forward_13) * dE * dpz * dz * V_r
N_all = np.sum(ekin_all_everywhere_13) * dE *dpz * dz * V_r
print("N_forward = ", N_forward)
print("N_all = ", N_all)

and I get
N_forward = 54148.884201168294
N_all = 70354.86803573875

Am I right in this approach? (Particle count below expected)

Thank you

[mccoys]

Please do not make an issue to ask a question. Use discussions

I don't understand why you sum energies when you just want a number. Just count the size of the array.

[mlyto]

Sorry for asking my question in the wrong place.
Thank you for your reply.
"I don't understand why you sum energies when you just want a number."
After binning, I have a distribution (over particle weights) of the type w(E,pz,z), then I need number of particles N = \int w(E, pz, z) dEdpzdz (over all energies E, pz > 0, i.e., only the second half, z>z_foc+10*l0 is the third part). Then I multiply the result by the hypervolume.
Am I understanding something wrong?
"Just count the size of the array." - the size of array is 100 (at fixed pz and z) - but this is not what I need.

[mccoys]

Ah ok you want a number of real particles. Sorry I misunderstood

If you are in 2D or 1D, then you have no way to get a number of particles in a well defined way. I suggest to use happi 's average option to sum over axes then check the units you obtain. You should only need to multiply by a volume. Happi takes into account the summing, the bin size, even when non uniform. It just does not work if you have custom axes but that is not your case

[mlyto]

Since this might be useful to someone, I leave my answer here.
In the code above, I've missed the multiplication by the transverse propagation lengths Lx and Ly (like dE, dpz and dz, these two quantities are expressed in Smiley units), so the correct formula for the number of particles would be:

N_forward = np.sum(ekin_forward_13) * dE * dpz * dz * V_r * Lx * Ly
N_all = np.sum(ekin_all_everywhere_13) * dE *dpz * dz * V_r * Lx * Ly

This applies when you use happi to read a ParticleBinning file. Thus, you can undo happi's actions when it reading the data.
If you're reading a ParticleBinning file without happi's intervention, as is typically done for hdf5/h5 files, simply multiplying the value by V_r will suffice.