Ad generation computing to Analysis.py

Thomas Fitoussi
1 parent 074c1729
Showing 6 changed files with 24 additions and 38 deletions Show diff stats
Analysis.py
Modules/Generation.py
Modules/Generation.pyc
Modules/Read.py
Modules/Read.pyc
Modules/Timing.py
 #!/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)
 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"
@@ -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)
@@ -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]