From 074c1729389f8fae1c62b0dbb1f1245a04272783 Mon Sep 17 00:00:00 2001 From: Thomas Fitoussi Date: Fri, 17 Jul 2015 11:53:35 +0200 Subject: [PATCH] Analysis.py used to pre-compute spectrum, timing and radial distribution computing done for power source spectrum (p= 1, 1.5, 2, 2.5) speed-up graph drawing, no need to recompute each time everything --- Analysis.py | 105 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------- Modules/Angle.py | 77 ++++++++++++++++++++++++++++------------------------------------------------- Modules/Angle.pyc | Bin 4135 -> 0 bytes Modules/Generation.py | 17 ++++++++++++++--- Modules/Generation.pyc | Bin 1358 -> 0 bytes Modules/Read.py | 15 +++++++++++++++ Modules/Read.pyc | Bin 3526 -> 0 bytes Modules/Spectrum.py | 33 ++++++++++++++++----------------- Modules/Spectrum.pyc | Bin 3186 -> 0 bytes Modules/Timing.py | 19 +++++++------------ Modules/Timing.pyc | Bin 3857 -> 0 bytes 11 files changed, 175 insertions(+), 91 deletions(-) diff --git a/Analysis.py b/Analysis.py index 110a4ec..ebe4354 100755 --- a/Analysis.py +++ b/Analysis.py @@ -1,22 +1,107 @@ #!/bin/python from sys import argv -from numpy import shape +from numpy import append,savetxt, shape, array, newaxis, empty +from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile 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() if shape(argv)[0] < 2: - error("not enough arguments (at least 2)") + error("not enough arguments (at least 1)") -# choose what to print -print "What will be printed?" -print " 1. Spectrum" -print " 2. Image" -print " 3. Radial distribution VS energy" -print " 4. Timing distribution VS arrival angle" -selection = raw_input("selection (number): ") +folder = "Results/" -if (selection == "1"): # draw spectrum +for fileId in argv[1:]: + # read files + time = ReadTime(fileId) + energy = ReadEnergy(fileId) + weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5]) + nbPhotonsEmitted=ReadProfile(fileId,[3]) + + Source = [] + Spectrum = [] + Radial = [] + Angle_Energy = [] + Timing = [] + + for powerlaw_index in [1,1.5,2,2.5]: + # apply source spectrum + weight = weightini*(Esource)**(1-powerlaw_index) + + #=============================================================================# + # SPECTRUM (SOURCE AND MEASURED) + #=============================================================================# + nbBins = 100 + # draw source spectrum + Es=array(list(set(Esource))) + Ws= (Es)**(1-powerlaw_index) + Es,Fs = spectrum(Es,Ws,nbBins) + Es=Es[:,newaxis] + Fs=Fs[:,newaxis] + if Source==[]: + Source = Es + Source = append(Source,Fs,axis=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) + ener=ener[:,newaxis] + flux=flux[:,newaxis] + flux_0=flux_0[:,newaxis] + if Spectrum==[]: + Spectrum = ener + Spectrum = append(Spectrum,flux,axis=1) + Spectrum = append(Spectrum,flux_0,axis=1) + else: + Spectrum = append(Spectrum,flux,axis=1) + Spectrum = append(Spectrum,flux_0,axis=1) + + #=============================================================================# + # RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY + #=============================================================================# + nbBins = 100 + ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted) + ener=ener[:,newaxis] + angle=angle[:,newaxis] + if Angle_Energy==[]: + Angle_Energy = ener + Angle_Energy = append(Angle_Energy,angle,axis=1) + else: + Angle_Energy = append(Angle_Energy,angle,axis=1) + + theta,dndtheta2 = radial(theta_arrival,weight,nbPhotonsEmitted) + theta=theta[:,newaxis] + dndtheta2=dndtheta2[:,newaxis] + if Radial==[]: + Radial = theta + Radial = append(Radial,dndtheta2,axis=1) + else: + Radial = append(Radial,dndtheta2,axis=1) + + #=============================================================================# + # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE + #=============================================================================# + nbBins = 100 + delta_t,dNdt = timing(time,weight,nbBins) + delta_t=delta_t[:,newaxis] + dNdt=dNdt[:,newaxis] + if Timing==[]: + Timing = delta_t + Timing = append(Timing,dNdt,axis=1) + else: + Timing = append(Timing,dNdt,axis=1) + + savetxt(folder+fileId+"/Spectrum.txt",Spectrum) + savetxt(folder+fileId+"/Source_spectrum.txt",Source) + savetxt(folder+fileId+"/Angle_vs_Energy.txt",Angle_Energy) + savetxt(folder+fileId+"/Radial_distribution.txt",Radial) + savetxt(folder+fileId+"/Timing.txt",Timing) diff --git a/Modules/Angle.py b/Modules/Angle.py index 42960dd..a73edf4 100644 --- a/Modules/Angle.py +++ b/Modules/Angle.py @@ -1,7 +1,7 @@ from numpy import log10, histogram from matplotlib.pyplot import figure, show from matplotlib import gridspec -from Read import ReadEnergy, ReadExtraFile, ReadProfile +from Read import ReadRadial, ReadAngleVsEnergy from Analytic import Analytic_theta from Constants import degre @@ -23,6 +23,24 @@ def angle_vs_energy(theta,energy,weight,nbPhotonsEmitted=1,nbBins=40): return enercenter, angle +def radial(theta,weight,nbPhotonsEmitted=1,nbBins=50,thetamax=90): + ''' + Compute flux versus arrival angle + Input: directory name + Input (optionnal): maximal angle desired (by default full sky) + Output: arrival angle, flux (dn/dtheta) + ''' + + theta=log10(theta) + dn,angle = histogram(theta,nbBins,range=[log10(1e-6),log10(thetamax)],weights=weight) + angle=10**angle + thetacenter = (angle[1:nbBins+1]+angle[0:nbBins])/2 + dtheta = angle[1:nbBins+1]-angle[0:nbBins] + dndtheta = dn/dtheta*thetacenter + dndtheta /= nbPhotonsEmitted + + return thetacenter, dndtheta + def drawAngle_vs_energy(files): ''' Plot arrival angle versus energy + analytic expression @@ -33,27 +51,11 @@ def drawAngle_vs_energy(files): ax1 = fig.add_subplot(111) for fileId in files: - energy = ReadEnergy(fileId) - weightini, theta_arrival, Esource = ReadExtraFile(fileId,[2,4,5]) - nbPhotonsEmitted=ReadProfile(fileId,[3]) - - # apply source spectrum - powerlaw_index = 1.5 - weight = weightini*(Esource/min(Esource))**(1-powerlaw_index) - ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted) - p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index) - - # apply source spectrum - powerlaw_index = 2 - weight = weightini*(Esource/min(Esource))**(1-powerlaw_index) - ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted) - p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index) - - # apply source spectrum - powerlaw_index = 2.5 - weight = weightini*(Esource/min(Esource))**(1-powerlaw_index) - ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted) - p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index) + i=1 + for powerlaw_index in [1,1.5,2,2.5]: + # apply source spectrum + ener,angle = ReadAngleVsEnergy(fileId,[0,i]) + p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index) yfit = Analytic_theta(ener,fileId) ax1.plot(ener,yfit,'--',color="m",linewidth=2,label="analytic") @@ -68,42 +70,19 @@ def drawAngle_vs_energy(files): show() -def radial(theta,weight,nbPhotonsEmitted=1,nbBins=50,thetamax=90): - ''' - Compute flux versus arrival angle - Input: directory name - Input (optionnal): maximal angle desired (by default full sky) - Output: arrival angle, flux (dn/dtheta) - ''' - - theta=log10(theta) - dn,angle = histogram(theta,nbBins,range=[log10(1e-6),log10(thetamax)],weights=weight) - angle=10**angle - thetacenter = (angle[1:nbBins+1]+angle[0:nbBins])/2 - dtheta = angle[1:nbBins+1]-angle[0:nbBins] - dndtheta = dn/dtheta*thetacenter - dndtheta /= nbPhotonsEmitted - - return thetacenter, dndtheta - -def drawRadial(files,thetamax=90): +def drawRadial(files): ''' Plot flux versus arrival angle Input: list of directories - Input (optionnal): maximal angle desired (by default full sky) Output: graph of flux versus angle ''' fig = figure() ax = fig.add_subplot(111) for fileId in files: - weightini, theta_arrival, Esource = ReadExtraFile(fileId,[2,4,5]) - nbPhotonsEmitted=ReadProfile(fileId,[3]) - - for powerlaw_index in [1.5,2,2.5]: - # apply source spectrum - weight = weightini*(Esource/min(Esource))**(1-powerlaw_index) - theta,dndtheta2=radial(theta_arrival,weight,nbPhotonsEmitted,thetamax) + i=1 + for powerlaw_index in [1,1.5,2,2.5]: + theta,dndtheta2 = ReadRadial(fileId,[0,i]) ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index) ax.legend(loc="best") diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc index cfe5530..6a0d4b6 100644 Binary files a/Modules/Angle.pyc and b/Modules/Angle.pyc differ diff --git a/Modules/Generation.py b/Modules/Generation.py index cea5271..2c303f0 100644 --- a/Modules/Generation.py +++ b/Modules/Generation.py @@ -1,6 +1,6 @@ from numpy import shape, arange from matplotlib.pyplot import figure, show, hist -from Read import ReadGeneration, ReadWeight +from Read import ReadExtraFile def drawGeneration(files): @@ -13,8 +13,7 @@ def drawGeneration(files): ax = fig.add_subplot(111) nbBins = 7 for fileId in files: - weight = ReadWeight(fileId) - generation = ReadGeneration(fileId) + weight, generation = ReadExtraFile(fileId,[2,3]) print "file", fileId, "->", shape(generation)[0], "events" hist(generation,nbBins,range=[2,9],log=1,alpha=.5, weights=weight,label=fileId) @@ -26,3 +25,15 @@ def drawGeneration(files): ax.set_ylabel("Number of events") show() + +def GenContribution(fileId): + ''' + compute the contribution to the final events of each + generation + ''' + weight, generation = ReadExtraFile(fileId,[2,3]) + + NbTotEvents = sum(weight) + + for gen in list(set(generation)): + print "Generation", gen, " -> ", sum(weight[generation==gen])/NbTotEvents *100, "\percent" diff --git a/Modules/Generation.pyc b/Modules/Generation.pyc index 257b80e..39fb30c 100644 Binary files a/Modules/Generation.pyc and b/Modules/Generation.pyc differ diff --git a/Modules/Read.py b/Modules/Read.py index 46b0501..85566d7 100644 --- a/Modules/Read.py +++ b/Modules/Read.py @@ -59,3 +59,18 @@ def ReadElmag(fileName): energy *= 10**(-9) #GeV flux=(fluxGamma+fluxLepton)*1e-9 #GeV normalized to 1 initial photon return energy,flux + +def ReadSpectrum(fileName,cols=[0,1,2,3,4,5,6,7,8]): + return loadtxt(resultDirectory+fileName+"/Spectrum.txt",unpack=True,usecols=cols) + +def ReadSourceSpectrum(fileName,cols=[0,1,2,3,4]): + return loadtxt(resultDirectory+fileName+"/Source_spectrum.txt",unpack=True,usecols=cols) + +def ReadAngleVsEnergy(fileName,cols=[0,1,2,3,4]): + return loadtxt(resultDirectory+fileName+"/Angle_vs_Energy.txt",unpack=True,usecols=cols) + +def ReadRadial(fileName,cols=[0,1,2,3,4]): + return loadtxt(resultDirectory+fileName+"/Radial_distribution.txt",unpack=True,usecols=cols) + +def ReadTiming(fileName,cols=[0,1,2,3,4]): + return loadtxt(resultDirectory+fileName+"/Timing.txt",unpack=True,usecols=cols) diff --git a/Modules/Read.pyc b/Modules/Read.pyc index 106b9c5..34b3354 100644 Binary files a/Modules/Read.pyc and b/Modules/Read.pyc differ diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py index fdec444..f44c12d 100644 --- a/Modules/Spectrum.py +++ b/Modules/Spectrum.py @@ -1,11 +1,11 @@ from numpy import histogram, log10, sqrt, array from matplotlib.pyplot import figure, show, setp from matplotlib import gridspec -from Read import ReadEnergy, ReadExtraFile, ReadProfile +from Read import ReadSpectrum, ReadSourceSpectrum from Constants import degre from Analytic import Ethreshold_gg -def spectrum(energy,weight,nbPhotonsEmitted=1,nbBins=100): +def spectrum(energy,weight,nbBins=100): ''' Compute flux versus energy Input: events energy and weight (optional: number of bins) @@ -17,7 +17,7 @@ def spectrum(energy,weight,nbPhotonsEmitted=1,nbBins=100): enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 binSize=ener[1:nbBins+1]-ener[0:nbBins] flux = enercenter**2 * flux/binSize - flux /= nbPhotonsEmitted + flux /= max(flux) #nbPhotonsEmitted return enercenter,flux def drawSpectrum(files,PlotAnalytic=False): @@ -34,23 +34,22 @@ def drawSpectrum(files,PlotAnalytic=False): ax2 = fig.add_subplot(gs[1],sharex=ax1) for fileId in files: - # read files - energy = ReadEnergy(fileId) - weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5]) - nbPhotonsEmitted=ReadProfile(fileId,[3]) - - for powerlaw_index in [1.5,2,2.5]: - # apply source spectrum - weight = weightini*(Esource/min(Esource))**(1-powerlaw_index) + i=1 + j=1 + for powerlaw_index in [1,1.5,2,2.5]: + # read files + Es,Fs = ReadSourceSpectrum(fileId,[0,i]) + p = ax1.plot(Es,Fs,drawstyle='--') # draw full spectrum - ener,flux = spectrum(energy,weight,nbPhotonsEmitted) - p = ax1.plot(ener,flux,drawstyle='-',label="p=%.1f"%powerlaw_index) + 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) # primary gamma-rays contribution - cond = (generation==0) - ener,flux_cond = spectrum(energy[cond],weight[cond],nbPhotonsEmitted) - ax1.plot(ener,flux_cond,color=p[0].get_color(),linestyle='-.') - contrib=flux_cond/flux + ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='-.') + contrib=flux_0/flux ax2.plot(ener,contrib,'-',color=p[0].get_color()) + + i+=1 + j+=2 # add spectrum for Fermi/LAT and CTA (FoV?) diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc index 2ce9e11..4ea2bcc 100644 Binary files a/Modules/Spectrum.pyc and b/Modules/Spectrum.pyc differ diff --git a/Modules/Timing.py b/Modules/Timing.py index 8386b55..1eb4a5d 100644 --- a/Modules/Timing.py +++ b/Modules/Timing.py @@ -1,10 +1,10 @@ from numpy import histogram, log10, shape, arange, pi from matplotlib.pyplot import figure, show -from Read import ReadTime, ReadExtraFile, ReadProfile, ReadEnergy +from Read import ReadTime, ReadTiming, ReadExtraFile, ReadProfile, ReadEnergy from Constants import yr, degre from Analytic import Analytic_delay_vs_theta -def timing(time,weight,NbPhotonsEmitted=1,nbBins=100): +def timing(time,weight,nbBins=100): ''' Compute flux versus time delay Input: directory name @@ -23,7 +23,7 @@ def timing(time,weight,NbPhotonsEmitted=1,nbBins=100): dN,dt=histogram(time,nbBins,weights=weight) timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2 binSize=dt[1:nbBins+1]-dt[0:nbBins] - dNdt=dN/(NbPhotonsEmitted*binSize) + dNdt=dN/(max(dN)*binSize) return 10**timecenter, dNdt @@ -36,14 +36,9 @@ def drawTiming(files): fig = figure() ax = fig.add_subplot(111) for fileId in files: - time = ReadTime(fileId) - weightini, Esource = ReadExtraFile(fileId,[2,5]) - nbPhotonsEmitted=ReadProfile(fileId,[3]) - - for powerlaw_index in [1.5,2,2.5]: - # apply source spectrum - weight = weightini*(Esource/min(Esource))**(1-powerlaw_index) - delta_t,dNdt = timing(time,weight,nbPhotonsEmitted) + i=1 + for powerlaw_index in [1,1.5,2,2.5]: + delta_t,dNdt = ReadTiming(fileId,[0,i]) ax.plot(delta_t,dNdt,linestyle="steps-mid",label="p=%.1f"%powerlaw_index) ax.set_xscale('log') @@ -51,7 +46,7 @@ def drawTiming(files): ax.grid(b=True,which='major') ax.legend(loc="best") ax.set_xlabel("Time delay [s]") - ax.set_ylabel("$t\ dN/dt$") + ax.set_ylabel("$t\ dN/dt$ [arbitrary units]") show() diff --git a/Modules/Timing.pyc b/Modules/Timing.pyc index 5d8404f..f0b9fb5 100644 Binary files a/Modules/Timing.pyc and b/Modules/Timing.pyc differ -- libgit2 0.21.2