import numpy as np
import math

#massmax = 1.2e9
massmin = 1.0e6
bins =31  # 300
binwidth = 0.1 #math.log10(massmax/massmin)/bins
massmax = 10**(math.log10(massmin)+bins*binwidth)
#eV_erg = 6.24e11

f=open('PopII_haloes.txt','r')

fo = open('PopII_histogram.txt','w')

halolist = []
PopII = []
PopIImass = []

line = f.readline()
while line != '':
   halomass = float(line.split()[0])
   halolist.append(halomass)
   PopII.append(int(line.split()[1]))
   PopIImass.append(float(line.split()[2])*float(line.split()[1]))
   line = f.readline()

halocount = np.zeros(bins+1)
totalPopII = np.zeros(bins+1)
totalmass = np.zeros(bins+1)
totalhalo = np.zeros(bins+1)

for halomass,star,mass in zip(halolist,PopII,PopIImass):
   index = int((math.log10(halomass/massmin)+0.5*binwidth)/binwidth)
   if (star) > 0:
     halocount[index] = halocount[index] + 1
     if star >0:
         totalPopII[index] = totalPopII[index] + star
         totalmass[index] = totalmass[index] +mass
   totalhalo[index] = totalhalo[index] + 1

cumulationPopII = 0.0
cumulationmass = 0.0
bin = math.log10(massmin) + 0.5*binwidth
for count,star,mass,halo in zip(halocount,totalPopII,totalmass,totalhalo):
    cumulationPopII = cumulationPopII + star
    cumulationmass = cumulationmass + mass
    fo.write(str(10**bin)+" "+str(count)+" "+str(halo)+" "+str(star)+" "+str(cumulationPopII)+\
     " "+str(mass)+" "+str(cumulationmass)+" "+str(((star)/count if count > 0 else 0))+"\n")
    bin =  bin + binwidth

f.close()
fo.close()


