I am trying to perform calculations that scale as (n_obj choose 2 (or 3) )*N_events over a large data set where n_obj is on the order of 1000, and N_events is on the order of 10^6. Thus, it makes sense to multithread over each event.

My issue comes down to the following. I am trying to port some of my code over from C++ to Python and I am not as familiar with Python, and especially not the multiprocessing library. Currently I am doing in Python what I did in c++ which is to create an array/list with size N_events, where each element is an array corresponding to a binning in a histogram. I know that the calculated values in the process are correct, but when I try to access the variable in the main body of my program, I get all zeros. In c++, I just pass everything by reference, and that works ok, but it doesn't appear to work in Python. I was wondering how to modify the below code to preserve the values for further processing

ncore = 20 
def extractParticles(flatparts):
    particles=np.ones((int(len(flatparts)/4), 3))
    maxval=0
    if len(particles) == 0:
        return particles
    for i in range(len(particles)):
        partpt=flatparts[int(i*4)]
        parteta=flatparts[int(i*4+1)]
        partphi=flatparts[int(i*4+2)]
        if partpt==0:
            break
        maxval+=1
        particles[i]=[partpt, parteta, partphi]
    particles=particles[0:maxval]
    return particles
def calculateR2(particle1, particle2):
    deltaphi = abs(particle1[2]-particle2[2])
    if deltaphi > np.pi:
        deltaphi+=-2*np.pi
    return pow(deltaphi, 2)+pow(particle1[1]-particle2[1], 2)
    
def calculateEEC(particles, EECbins, output, njet):
    npart = len(particles)
    ncalc = npart*(npart-1)
    #if njet % ncore == 0: 
     #   print("Jet" + str(njet)+" has "+str(npart)+" particles, which means we need to do "+str(ncalc)+" calculations")
    for i in range(len(particles)-1):
       # if njet % ncore == 0:
      #      print("Doing calculations on particle "+str(i)+" [ "+str(ncalc)+" calculations remaining ]", end="\r", flush=True)
        for j in range(i, len(particles)):
            R2=calculateR2(particles[i], particles[j])
            e2=particles[i][0]*particles[j][0]
           # if njet == 20 and j ==i+1:
           #     print(e2)
            for k in range(len(EECbins)):
                if R2 > EECbins[k]:
                    continue
                else:
                    output[k]+=e2
                    break
            ncalc+=-1
    if njet == 20:
        print("Squares of energy are: "+str(output))
    return 
        
def getEEC(jets):
    minR=1e-2
    maxR=0.8
    minR2=pow(minR,2)
    maxR2=pow(maxR,2)
    minbin=np.log(minR2)
    maxbin=np.log(maxR2)
    #print(maxbin)
    nbins=abs((maxR2-minR2)/(minR2))
    print(nbins)
    bins=np.ones(int(nbins))
    for i in range(len(bins)):
        if i==0:
            bins[i]=minR2
        else:
            binval=minR2*(i+1)
            bins[i]=binval
    eec=np.zeros(int(nbins))
   
    jetsinverted=jets.iloc[:,8:].values
    alljetsoutput=np.zeros((len(jetsinverted), len(eec)))
    p=[]
    chunks=100*ncore #run 25 chunks before joining them and then moving onto the next group
    print("0% Complete [ 0 / "+str(len(jets))+" jets]", end='\r', flush=True) 
    for i1 in range(0, len(jetsinverted), chunks):
        for i in range(0, chunks):

            if i1+i >=len(jetsinverted):
                break
            jet=jetsinverted[i1+i]
            particles=extractParticles(jet)
            p1=Process(target=calculateEEC, args=(particles, bins, alljetsoutput[i1+i], i1+i))
            p.append(p1)
        lastmax=0
        for i in range(len(p)):
            p[i].start()
            if i % ncore == 0:
                for j in range(lastmax, i+1):
                    p[j].join()
                    lastmax=i+1
                    pctcomp=(i1+i)/len(jetsinverted)*100.
                    pctcompstr="{:.2f}".format(pctcomp)
                    nj="{:.2e}".format(i+i1)
                    nt="{:.2e}".format(len(jets))
                    print( pctcompstr+" % Complete [ " +nj+" / "+nt+" jets]", end='\r', flush=True)
                    if j == 15: 
                        print(alljetsoutput[j])
                    #print(f"{0:.2f}% Complete [ {1:.2e} / {2:.2e} jets]".format(pctcomp, i+i1, len(jetsinverted)), end='\r', flush=True)
        p.clear()

    print()
    for i in range(len(alljetsoutput)):    
        onejet=alljetsoutput[i]
        for j in range(len(onejet)):
            eec[j]+=onejet[j]#/pow(jets.iloc[i,1], 2) #per jet normalization
           # print("RL = "+str(bins[j])+" : weight = "+str(eec[j]))
    plt.hist(eec, bins=bins)
    return bins, eec