From 631c85a94f95a402e462ad6a1e0eaf4ec7ab3d86 Mon Sep 17 00:00:00 2001 From: Thomas Fitoussi Date: Fri, 17 Jul 2015 18:54:47 +0200 Subject: [PATCH] Ad generation computing to Analysis.py --- Analysis.py | 16 +++++++++++++++- Modules/Generation.py | 24 ++++++------------------ Modules/Generation.pyc | Bin 1678 -> 0 bytes Modules/Read.py | 18 +++--------------- Modules/Read.pyc | Bin 5588 -> 0 bytes Modules/Timing.py | 4 ---- 6 files changed, 24 insertions(+), 38 deletions(-) diff --git a/Analysis.py b/Analysis.py index ebe4354..291ea55 100755 --- a/Analysis.py +++ b/Analysis.py @@ -1,7 +1,7 @@ #!/bin/python from sys import argv -from numpy import append,savetxt, shape, array, newaxis, empty +from numpy import append, savetxt, shape, array, newaxis, zeros, arange from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile from Modules.Spectrum import spectrum from Modules.Angle import radial, angle_vs_energy @@ -23,6 +23,8 @@ for fileId in argv[1:]: weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5]) nbPhotonsEmitted=ReadProfile(fileId,[3]) + Gen_contrib = [] + Gen_cont = zeros((int(max(generation))+1)) Source = [] Spectrum = [] Radial = [] @@ -34,6 +36,17 @@ for fileId in argv[1:]: weight = weightini*(Esource)**(1-powerlaw_index) #=============================================================================# + # GENERATION CONTRIBUTION + #=============================================================================# + NbTotEvents = sum(weight) + if Gen_contrib==[]: + Gen_contrib = arange(0,max(generation)+1,1)[:,newaxis] + + for gen in list(set(generation)): + Gen_cont[int(gen)] = sum(weight[generation==gen])/NbTotEvents *100 + Gen_contrib = append(Gen_contrib,Gen_cont[:,newaxis],axis=1) + + #=============================================================================# # SPECTRUM (SOURCE AND MEASURED) #=============================================================================# nbBins = 100 @@ -100,6 +113,7 @@ for fileId in argv[1:]: else: Timing = append(Timing,dNdt,axis=1) + savetxt(folder+fileId+"/Generation.txt",Gen_contrib) savetxt(folder+fileId+"/Spectrum.txt",Spectrum) savetxt(folder+fileId+"/Source_spectrum.txt",Source) savetxt(folder+fileId+"/Angle_vs_Energy.txt",Angle_Energy) diff --git a/Modules/Generation.py b/Modules/Generation.py index 2c303f0..6ff8b4f 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 ReadExtraFile +from Read import ReadGeneration def drawGeneration(files): @@ -11,29 +11,17 @@ def drawGeneration(files): ''' fig = figure() ax = fig.add_subplot(111) - nbBins = 7 for fileId in files: - 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) + 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) xtick=arange(9) ax.legend(loc='best') ax.set_xticks(xtick+0.5) ax.set_xticklabels(xtick) ax.set_xlabel("Generation") - ax.set_ylabel("Number of events") + ax.set_ylabel("Ratio 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 39fb30c..0bbe0d1 100644 Binary files a/Modules/Generation.pyc and b/Modules/Generation.pyc differ diff --git a/Modules/Read.py b/Modules/Read.py index 85566d7..a43c725 100644 --- a/Modules/Read.py +++ b/Modules/Read.py @@ -38,21 +38,6 @@ def ReadMomentum(fileName): def ReadMomentumAngle(fileName): return loadtxt(resultDirectory+fileName+momentumFile,unpack=True,usecols=[4,5]) -def ReadCharge(fileName): - return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[0]) - -def ReadRedshift(fileName): - return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[1]) - -def ReadWeight(fileName): - return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[2]) - -def ReadGeneration(fileName): - return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[3]) - -def ReadArrivalAngle(fileName): - return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[4]) - def ReadElmag(fileName): Elmag=resultDirectory+fileName+ElmagFile energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2]) @@ -74,3 +59,6 @@ def ReadRadial(fileName,cols=[0,1,2,3,4]): def ReadTiming(fileName,cols=[0,1,2,3,4]): return loadtxt(resultDirectory+fileName+"/Timing.txt",unpack=True,usecols=cols) + +def ReadGeneration(fileName,cols=[0,1,2,3,4]): + return loadtxt(resultDirectory+fileName+"/Generation.txt",unpack=True,usecols=cols) diff --git a/Modules/Read.pyc b/Modules/Read.pyc index 34b3354..88498f4 100644 Binary files a/Modules/Read.pyc and b/Modules/Read.pyc differ diff --git a/Modules/Timing.py b/Modules/Timing.py index 1eb4a5d..385a1ef 100644 --- a/Modules/Timing.py +++ b/Modules/Timing.py @@ -11,10 +11,6 @@ def timing(time,weight,nbBins=100): Input (optional): generation (default = all) Output: time delay (sec), flux ''' - - prop = float(shape(time[time<0])[0])/shape(time)[0] - print shape(time)[0], "events", "negative time:",shape(time[time<0])[0], "~", prop - time=time[time>0] weight=weight[time>0] -- libgit2 0.21.2