From 374f3d52ed4a3587b3d9f2cb7c21bbd0af4787a1 Mon Sep 17 00:00:00 2001 From: Thomas Fitoussi Date: Tue, 28 Apr 2015 15:10:50 +0200 Subject: [PATCH] add Weight distribution and versus energy spectrum can be used with generation selection --- Analysis.py | 7 +++++-- EBL-spectrums.py | 2 +- Modules/Analytic.py | 4 ++-- Modules/Analytic.pyc | Bin 3862 -> 0 bytes Modules/Read.py | 4 ++-- Modules/Read.pyc | Bin 4603 -> 0 bytes Modules/Spectrum.py | 94 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------- Modules/Spectrum.pyc | Bin 2824 -> 0 bytes Modules/Weight.py | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 157 insertions(+), 20 deletions(-) create mode 100644 Modules/Weight.py diff --git a/Analysis.py b/Analysis.py index adc6c4a..259e824 100755 --- a/Analysis.py +++ b/Analysis.py @@ -2,7 +2,7 @@ from sys import argv from numpy import shape -from Modules.Spectrum import drawSpectrum +from Modules.Spectrum import drawSpectrum, drawSpectrum_normaLepton from Modules.Angle import drawAngle, drawRadial from Modules.Timing import drawTiming from Modules.Generation import drawGeneration @@ -25,9 +25,12 @@ def error(error_type): if shape(argv)[0] < 3: error("not enough argument") -if argv[1] == "spectrum": +elif argv[1] == "spectrum": drawSpectrum(argv[2:]) +elif argv[1] == "spectrumLept": + drawSpectrum_normaLepton(argv[2:]) + elif argv[1] == "angle": drawAngle(argv[2:]) diff --git a/EBL-spectrums.py b/EBL-spectrums.py index 9e675df..a1d8b4a 100755 --- a/EBL-spectrums.py +++ b/EBL-spectrums.py @@ -3,7 +3,7 @@ import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as plt -from Constantes import * +from Modules.Constants import * fig = plt.figure() ax = fig.add_subplot(111) diff --git a/Modules/Analytic.py b/Modules/Analytic.py index 2927817..af5dfa4 100644 --- a/Modules/Analytic.py +++ b/Modules/Analytic.py @@ -22,8 +22,8 @@ def Analytic_theta(E_ic,fileId): return abs(arcsin(delta_to_theta*sin(delta))*degre) # Compton accumulation -def ECompton_threshold(Compton_threshold = 0.05): - return Compton_threshold/(4/3*2.7*Tcmbp) *me*1e-6 #GeV +def ECompton_threshold(Compton_threshold = 0.005): + return Compton_threshold/(4/3*Ecmb/me*1e-3) *me*1e-6 #GeV # Compton scattering def Dic(Ee=Ee_default): # Ee (GeV) diff --git a/Modules/Analytic.pyc b/Modules/Analytic.pyc index ac15ee1..aebef08 100644 Binary files a/Modules/Analytic.pyc and b/Modules/Analytic.pyc differ diff --git a/Modules/Read.py b/Modules/Read.py index 3ffcf9d..0066acf 100644 --- a/Modules/Read.py +++ b/Modules/Read.py @@ -40,8 +40,8 @@ def ReadGeneration(fileName): def ReadElmag(fileName): Elmag=resultDirectory+fileName+ElmagFile energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2]) - energy *= 10**(-9) - flux=fluxGamma+fluxLepton + energy *= 10**(-9) #GeV + flux=(fluxGamma+fluxLepton)*1e-9 #GeV normalized to 1 initial photon return energy,flux def ReadE_source(fileName): # Emax (GeV), redshift, distance (Mpc), Nb photons diff --git a/Modules/Read.pyc b/Modules/Read.pyc index 55a26a4..e80e7d5 100644 Binary files a/Modules/Read.pyc and b/Modules/Read.pyc differ diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py index 5a442c6..fdaad0c 100644 --- a/Modules/Spectrum.py +++ b/Modules/Spectrum.py @@ -1,14 +1,23 @@ -from numpy import histogram, log10, shape, logspace +from numpy import histogram, log10, shape, logspace, sqrt from matplotlib.pyplot import figure, show, ylim, setp from matplotlib import gridspec from os.path import isfile from Read import ReadEnergy, ReadWeight, ReadElmag, ReadNbLeptonsProd, resultDirectory, ElmagFile -from Analytic import Analytic_flux +from Read import ReadNphot_source, ReadGeneration +from Analytic import Analytic_flux, Ethreshold_gg + +def spectrum(fileId,gen=-1): + energy = ReadEnergy(fileId) + weight = ReadWeight(fileId) + + if gen != -1: + generation = ReadGeneration(fileId) + energy=energy[generation==gen] + weight=weight[generation==gen] + print "file", fileId, "gen=", gen,"->", shape(energy)[0], "events" + else: + print "file", fileId, "->", shape(energy)[0], "events" -def spectrum(fileName): - energy = ReadEnergy(fileName) - weight = ReadWeight(fileName) - print "file", fileName, "->", shape(energy)[0], "events" nbBins=100 energy =log10(energy) flux,ener=histogram(energy,nbBins,weights=weight) @@ -16,9 +25,46 @@ def spectrum(fileName): enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 binSize=ener[1:nbBins+1]-ener[0:nbBins] flux = enercenter**2 * flux/binSize + maxflux = ReadNphot_source(fileId) + flux /= maxflux return enercenter,flux -def drawSpectrum(files): +########################################################################################### +def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False): + fig = figure() + ax = fig.add_subplot(111) + + for fileId in files: + energy,flux = spectrum(fileId) + if PlotElmag: + maxflux = max(flux)*(.8) + p=ax.plot(energy,flux,drawstyle='steps-mid',label=fileId) + + if PlotElmag & isfile(resultDirectory+fileId+ElmagFile): + energy, flux = ReadElmag(fileId) + maxflux = max(flux) + flux /= maxflux + ax.plot(energy,flux,"--",color=p[0].get_color(),label="Elmag - "+fileId) + + if PlotAnalytic: + alpha=5e3 + ax.plot(energy,alpha*sqrt(energy),'--m',linewidth=2) + ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2) + ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2) + + ax.set_xscale('log') + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + if PlotElmag: + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [normalized]") + else: + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV / photon]") + ax.set_xlabel("energy [GeV]") + + show() + +def drawSpectrum_normaLepton(files): fig = figure() gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) fig.subplots_adjust(hspace=0.001) @@ -37,16 +83,11 @@ def drawSpectrum(files): error=flux/yfit-1 ax2.plot(energy,error,'+',color=p[0].get_color()) - Elmag=resultDirectory+fileId+ElmagFile - if isfile(Elmag): - energy, flux = ReadElmag(fileId) - ax1.plot(energy,flux,"-",color=p[0].get_color(),label="Elmag - "+fileId) - ax1.set_xscale('log') ax1.set_yscale('log') ax1.grid(b=True,which='major') ax1.legend(loc="best") - ax1.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV]") + ax1.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV] / lepton") ax2.set_xscale('log') ax2.grid(b=True,which='major') ax2.set_xlabel("energy [GeV]") @@ -57,3 +98,30 @@ def drawSpectrum(files): setp(xticklabels, visible=False) show() + +######################################################################################## +def drawSpectrumGen(fileId): + fig = figure() + ax = fig.add_subplot(111) + + for gen in list(set(ReadGeneration(fileId))): + energy,flux = spectrum(fileId,gen) + p=ax.plot(energy,flux,drawstyle='steps-mid',label="gen="+str(int(gen))) + + energy,flux = spectrum(fileId) + p=ax.plot(energy,flux,drawstyle='steps-mid',color='k',label="gen = all") + alpha=2e3 + + ax.plot(energy,sqrt(energy)*alpha,'--m',linewidth=2) + ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2) + ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2) + + ax.set_xscale('log') + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + ax.set_title(fileId) + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV / photon]") + ax.set_xlabel("energy [GeV]") + + show() diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc index bc0563e..597e17c 100644 Binary files a/Modules/Spectrum.pyc and b/Modules/Spectrum.pyc differ diff --git a/Modules/Weight.py b/Modules/Weight.py new file mode 100644 index 0000000..1a94bf1 --- /dev/null +++ b/Modules/Weight.py @@ -0,0 +1,66 @@ +from numpy import shape, arange, log10, histogram +from matplotlib.pyplot import figure, show, hist +from Read import ReadGeneration, ReadWeight, ReadEnergy + + +def drawWeightHisto(fileId): + fig = figure() + ax = fig.add_subplot(111) + nbBins =100 + + weight = ReadWeight(fileId) + generation = ReadGeneration(fileId) + print "file", fileId, "->", shape(generation)[0], "events" + + for gen in list(set(generation)): + hist(log10(weight[generation==gen]),nbBins,log=1,range=[0,7],alpha=.5, + label="gen="+str(int(gen))) + + ax.legend(loc='best') + ax.set_xlabel("Weight [log]") + ax.set_ylabel("Number of events") + + show() + +####################################################################################### +def weightVsEnergy(fileId,gen=-1): + energy = ReadEnergy(fileId) + weight = ReadWeight(fileId) + nbBins = 100 + + if gen != -1: + generation = ReadGeneration(fileId) + energy=energy[generation==gen] + weight=weight[generation==gen] + nbBins =40 + + weight,ener=histogram(log10(energy),nbBins,weights=log10(weight)) + nb,ener=histogram(log10(energy),nbBins) + ener=10**ener + enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 + weight /= nb + + return enercenter, weight + +def drawWeightVsEnergy(fileId): + fig = figure() + ax = fig.add_subplot(111) + + generation = ReadGeneration(fileId) + print "file", fileId, "->", shape(generation)[0], "events" + + for gen in list(set(generation)): + energy,weight = weightVsEnergy(fileId,gen) + ax.plot(energy,weight,".",label="gen="+str(int(gen))) + + energy,weight = weightVsEnergy(fileId) + ax.plot(energy,weight,"k",label="gen=all") + + ax.set_xscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + ax.set_title(fileId) + ax.set_xlabel("energy [GeV]") + ax.set_ylabel("mean weight [log]") + + show() -- libgit2 0.21.2