diff --git a/Analysis.py b/Analysis.py index 291ea55..de9ebf8 100755 --- a/Analysis.py +++ b/Analysis.py @@ -7,14 +7,18 @@ from Modules.Spectrum import spectrum from Modules.Angle import radial, angle_vs_energy from Modules.Timing import timing -def error(error_type): - print error_type - exit() +def InProgress(fileId,i,Nmax): + print "Work on ", fileId, ": ", (i*100)/Nmax, "% done" -if shape(argv)[0] < 2: - error("not enough arguments (at least 1)") +if shape(argv)[0] < 2: + print "not enough arguments (at least 1)" + exit() +#================================================================# +# Parameters +#================================================================# folder = "Results/" +PowerSpectrum=[1,1.5,2,2.5] for fileId in argv[1:]: # read files @@ -23,6 +27,17 @@ for fileId in argv[1:]: weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5]) nbPhotonsEmitted=ReadProfile(fileId,[3]) + lim =100000 + time=time[:lim] + energy=energy[:lim] + weightini=weightini[:lim] + generation = generation[:lim] + theta_arrival = theta_arrival[:lim] + Esource=Esource[:lim] + + cond = (generation==0) + test=energy[cond]/Esource[cond] + Gen_contrib = [] Gen_cont = zeros((int(max(generation))+1)) Source = [] @@ -31,9 +46,15 @@ for fileId in argv[1:]: Angle_Energy = [] Timing = [] - for powerlaw_index in [1,1.5,2,2.5]: + for powerlaw_index in PowerSpectrum: # apply source spectrum - weight = weightini*(Esource)**(1-powerlaw_index) + weight_source = Esource**(1-powerlaw_index) + weight = weightini* weight_source + + print "%e" %(max(weight_source)/min(weight_source)) + + #cond = (energy>95) & (energy<105) & (generation==0) + #print weightini[cond], weight[cond] #=============================================================================# # GENERATION CONTRIBUTION @@ -52,8 +73,9 @@ for fileId in argv[1:]: nbBins = 100 # draw source spectrum Es=array(list(set(Esource))) + print shape(Es) Ws= (Es)**(1-powerlaw_index) - Es,Fs = spectrum(Es,Ws,nbBins) + Es,Fs = spectrum(Es,Ws,nbBins=nbBins) Es=Es[:,newaxis] Fs=Fs[:,newaxis] if Source==[]: @@ -62,11 +84,8 @@ for fileId in argv[1:]: else: Source = append(Source,Fs,axis=1) - # primary gamma-rays contribution - cond = (generation==0) - ener,flux_0 = spectrum(energy[cond],weight[cond],nbBins) - # full spectrum - ener,flux = spectrum(energy,weight,nbBins) + # primary gamma-rays contribution and full spectrum + ener,flux,flux_0 = spectrum(energy,weight,generation,nbBins) ener=ener[:,newaxis] flux=flux[:,newaxis] flux_0=flux_0[:,newaxis] @@ -113,6 +132,8 @@ for fileId in argv[1:]: else: Timing = append(Timing,dNdt,axis=1) + InProgress(fileId,PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0]) + savetxt(folder+fileId+"/Generation.txt",Gen_contrib) savetxt(folder+fileId+"/Spectrum.txt",Spectrum) savetxt(folder+fileId+"/Source_spectrum.txt",Source) diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc index 6a0d4b6..0adbeef 100644 Binary files a/Modules/Angle.pyc and b/Modules/Angle.pyc differ diff --git a/Modules/Generation.py b/Modules/Generation.py index 6ff8b4f..92c2f68 100644 --- a/Modules/Generation.py +++ b/Modules/Generation.py @@ -15,12 +15,11 @@ def drawGeneration(files): i=1 for powerlaw_index in [1,1.5,2,2.5]: generation, ratio = ReadGeneration(fileId,[0,i]) - ax.plot(generation-0.5,ratio,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index) + ax.plot(generation,ratio,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index) xtick=arange(9) ax.legend(loc='best') - ax.set_xticks(xtick+0.5) - ax.set_xticklabels(xtick) + ax.set_yscale('log') ax.set_xlabel("Generation") ax.set_ylabel("Ratio of events [%]") diff --git a/Modules/Generation.pyc b/Modules/Generation.pyc index 0bbe0d1..b71b40b 100644 Binary files a/Modules/Generation.pyc and b/Modules/Generation.pyc differ diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py index f44c12d..48e6d88 100644 --- a/Modules/Spectrum.py +++ b/Modules/Spectrum.py @@ -5,20 +5,32 @@ from Read import ReadSpectrum, ReadSourceSpectrum from Constants import degre from Analytic import Ethreshold_gg -def spectrum(energy,weight,nbBins=100): +def spectrum(energy,weight,generation=[],nbBins=100): ''' Compute flux versus energy Input: events energy and weight (optional: number of bins) Output: energy, flux (E^2 dN/dE) ''' + Emin=min(energy) + Emax=max(energy) energy =log10(energy) - flux,ener=histogram(energy,nbBins,weights=weight) + flux,ener=histogram(energy,nbBins,range=[log10(Emin),log10(Emax)],weights=weight) ener = 10**ener enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 binSize=ener[1:nbBins+1]-ener[0:nbBins] flux = enercenter**2 * flux/binSize - flux /= max(flux) #nbPhotonsEmitted - return enercenter,flux + + if generation == []: + return enercenter,flux + else: + cond=(generation==0) + flux_0,ener_0=histogram(energy[cond],nbBins,range=[log10(Emin),log10(Emax)],weights=weight[cond]) + ener_0 = 10**ener_0 + enercenter_0=(ener_0[1:nbBins+1]+ener_0[0:nbBins])/2 + binSize=ener_0[1:nbBins+1]-ener_0[0:nbBins] + flux_0 = enercenter_0**2 * flux_0/binSize + return enercenter,flux,flux_0 + def drawSpectrum(files,PlotAnalytic=False): ''' @@ -34,19 +46,19 @@ def drawSpectrum(files,PlotAnalytic=False): ax2 = fig.add_subplot(gs[1],sharex=ax1) for fileId in files: - i=1 - j=1 + i=1#3 + j=1#5 for powerlaw_index in [1,1.5,2,2.5]: # read files Es,Fs = ReadSourceSpectrum(fileId,[0,i]) - p = ax1.plot(Es,Fs,drawstyle='--') + p = ax1.plot(Es,Fs,drawstyle='.') # draw full spectrum ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1]) - ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='-',label="p=%.1f"%powerlaw_index) + ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index) # primary gamma-rays contribution - ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='-.') + ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid') contrib=flux_0/flux - ax2.plot(ener,contrib,'-',color=p[0].get_color()) + ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color()) i+=1 j+=2 diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc index 4ea2bcc..33b6326 100644 Binary files a/Modules/Spectrum.pyc and b/Modules/Spectrum.pyc differ diff --git a/Modules/Timing.py b/Modules/Timing.py index 385a1ef..eff938a 100644 --- a/Modules/Timing.py +++ b/Modules/Timing.py @@ -56,9 +56,9 @@ def drawDelay_vs_angle(fileId): ax = fig.add_subplot(111) nbBins = 100 - theta = ReadExtraFile(fileId,[4]) - energy = ReadEnergy(fileId) - delay = ReadTime(fileId) + theta = ReadExtraFile(fileId,[4]) + energy = ReadEnergy(fileId) + delay = ReadTime(fileId) for n in arange(0,8,2): Emin = 10**(3-n) #GeV @@ -74,14 +74,14 @@ def drawDelay_vs_angle(fileId): dtheta=10**dtheta thetacenter=(dtheta[1:nbBins+1]+dtheta[0:nbBins])/2 dt=dt/dN - ax.plot(thetacenter,dt,color='k',linestyle="steps-mid",linewidth=2) + ax.plot(thetacenter,dt,color='m',linestyle="steps-mid",linewidth=2) thmin = 1.e-8 thmax = pi/2. nth = 1000 th = thmin*(thmax/thmin)**(arange(nth)/(nth-1.)) delay_fit = Analytic_delay_vs_theta(th,fileId) - ax.plot(th,delay_fit,'--m',linewidth=2) + ax.plot(th,delay_fit,'--k',linewidth=2) ax.set_xscale('log') ax.set_yscale('log') diff --git a/Modules/Timing.pyc b/Modules/Timing.pyc index f0b9fb5..adf71e1 100644 Binary files a/Modules/Timing.pyc and b/Modules/Timing.pyc differ -- libgit2 0.21.2