I fixed some import statements and if conditions

ctoennis · Apr 7, 2021 · 8a084f6 · 8a084f6
1 parent c37f890
commit 8a084f6
Show file tree

Hide file tree

Showing 10 changed files with 483 additions and 310 deletions.
diff --git a/Acceptance.py b/Acceptance.py
@@ -59,8 +59,6 @@
 
 	        if os.path.isfile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_e.dat"):
 
-#		    print("wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_e.dat " + "exists")
-
 	            spe_e.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_e.dat"))
 		    Nxxx_e.append(spe_e[-1].Integral(spe_e[-1].points[0][0],spe_e[-1].points[-1][0]))
 		    spe_ae.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_ae.dat"))
@@ -119,6 +117,7 @@
 Acc_casc    = ING.H1D.Empty(-1.0*np.pi,np.pi,180)
 Acc_track   = ING.H1D.Empty(-1.0*np.pi,np.pi,180)
 
+#Sun visibility versus devlination
 
 hsun = ING.H1D.FromFile("/home/ctoennis/analyses/standard_analysis_framework/Signal_Ingredients/HSun.txt")
 
@@ -140,14 +139,16 @@
 	if "IC59" in filename:
 
 	    live = 348.138
+	    continue #skip IC59
 
 	elif "IC79" in filename:
 
 	    live = 310.000
+	    continue #skip IC79
 
 	elif "2011" in filename:
 
-	    live = 342.0883333
+	    live = 342.0883333 #livetime in days
 
 	else:
 
@@ -157,7 +158,7 @@
 
 	infile = np.load("/data/ana/analyses/northern_tracks/version-002-p00/"+filename,"r")
 
-	nevents = len(infile["orig_OW"])
+	nevents = len(infile["orig_OW"]) #events in the MC file
 
 #	print nevents
 	for entry in infile:
@@ -195,7 +196,7 @@
 
 					if mctype == 68 or mctype == 14:
 
-					    track[k][j][i][l] += weight*spectrum_mu[k][j][i][l].Eval(mcen/masslist[k])/((N_mu[k][j][i][l]+N_amu[k][j][i][l])*masslist[k])
+					    track[k][j][i][l] += weight*spectrum_mu[k][j][i][l].Eval(mcen/masslist[k])/((N_mu[k][j][i][l]+N_amu[k][j][i][l])*masslist[k]) #separate treatment of neutrinos and antineutrinos
 
 					else:
 				    
@@ -213,7 +214,7 @@
 
         for channel in channellist:
 
-	    outfi.append(open("/home/ctoennis/analyses/standard_analysis_framework/northern_ing/ACC_med" + str(medmass) + "_lifetime" + str(lifetime) + "_channel" + str(channel) + ".txt","w"))
+	    outfi.append(open("/home/ctoennis/analyses/standard_analysis_framework/northern_ing/ACC_med" + str(medmass) + "_lifetime" + str(lifetime) + "_channel" + str(channel) + ".txt","w")) 
 
 	outfil.append(outfi)
 
@@ -229,7 +230,7 @@
 
 	        if not spectrum_mu[k][j][i][l] == "None":
 
-	            outfile[j][i][l].write(str(masslist[k]) + "    " + str(track[k][j][i][l]) + "\n")
+	            outfile[j][i][l].write(str(masslist[k]) + "    " + str(track[k][j][i][l]) + "\n") #frite output to simple text file
 
 
 

diff --git a/Background.py b/Background.py
@@ -2,7 +2,7 @@
 #METAPROJECT icerec/V05-01-00
 
 import numpy as np
-import ingredients as ING
+import src/ING as ING
 import h5py
 import os
 
@@ -12,46 +12,33 @@
 
 hsun.Scale(hsun.nbin*1.0/(hsun.integral))
 
-casc  = ING.H1D.Empty(-1.0,1.0,800)
 track = ING.H1D.Empty(0.998,1.0,400)
 
-#norm = [1.0/(np.sin(np.radians(x))) for x in casc.center]
+#background is generated from time scrambled data by only usng event declination
 
 for filename in os.listdir("/data/ana/analyses/northern_tracks/version-002-p00/"):
 
-    if "exp.npy" in filename:
+    if "exp.npy" in filename and "IC86" in filename: #Only use data files and IC86
 
 	infile = np.load("/data/ana/analyses/northern_tracks/version-002-p00/"+filename,"r")
 
 	for decr,azr,erec in zip(infile["zenith_MPEFit"],infile["azimuth_MPEFit"],infile["logE"]):
 
-#	    print(decr)
 
 	    if np.power(10.0,erec) < 200000.0:
 
 		for azt in map(lambda x: x*np.pi/20.0,range(1,40)):
 
                     for dect in [np.arccos(hsun.GetRandom()) for x in [0]*30]:
 
-
-#		        print("new")
-#		        print(dect)
-#			print(decr)
-
 		        diff = (np.dot(np.array([np.sin(azt)*np.sin(dect),np.cos(azt)*np.sin(dect),np.cos(dect)]),np.array([np.sin(azr)*np.sin(decr),np.cos(azr)*np.sin(decr),np.cos(decr)])))
 
 			if diff > 0.998:
-#			    print(np.arccos(diff))
 
 			    track.Fill(diff,1.0/(40.0*30.0))
 
 
 
-#casc.Write("/home/ctoennis/analyses/standard_analysis_framework/Signal_Ingredients/Background_cscd.txt")
-#casc.Scale(len(casc.content)/casc.integral)
-#casc.Multiply(norm)
-#casc.Write("/home/ctoennis/analyses/standard_analysis_framework/Signal_Ingredients/Background_S_cscd.txt")
 track.Write("/home/ctoennis/analyses/standard_analysis_framework/northern_ing/Background_WIMP.txt")
 track.Scale(len(track.content)/track.integral)
-#track.Multiply(norm)
 track.Write("/home/ctoennis/analyses/standard_analysis_framework/northern_ing/Background_S_WIMP.txt")
diff --git a/Background_ENE.py b/Background_ENE.py
@@ -2,7 +2,7 @@
 #METAPROJECT icerec/V05-01-00
 
 import numpy as np
-import ingredients as ING
+import src/ING as ING
 import h5py
 import os
 
@@ -12,11 +12,10 @@
 
 track = ING.H1D.Empty(0.0,200000.0,2000)
 
-#norm = [1.0/(np.sin(np.radians(x))) for x in casc.center]
 
 for filename in os.listdir("/data/ana/analyses/northern_tracks/version-002-p00/"):
 
-    if "exp.npy" in filename:
+    if "exp.npy" in filename and "IC86" in filename: #only use IC86 data files
 
 	infile = np.load("/data/ana/analyses/northern_tracks/version-002-p00/"+filename,"r")
 

diff --git a/Energy.py b/Energy.py
@@ -4,7 +4,7 @@
 # This code generates histograms for the likelihood that describes the reconstructed signal neutrino energy
 
 import numpy as np
-import ingredients as ING #contains h
+import src/ING as ING #contains histograms
 import h5py
 import os
 import argparse
@@ -24,48 +24,32 @@
 lifetimelist = [1, 10, 100, 1000, 10000]
 channellist  = [5, 8, 11, 13]
 track        = []
-#spectrum_e   = []
-#spectrum_ae  = []
 spectrum_mu  = []
 spectrum_amu = []
-#N_e          = []
-#N_ae         = []
 N_mu         = []
 N_amu        = []
 
 for mass in masslist:
 
     trax      = []
-#    spect_e   = []
-#    spect_ae  = []
     spect_mu  = []
     spect_amu = []
-#    Nx_e      = []
-#    Nx_ae     = []
     Nx_mu     = []
     Nx_amu    = []
 
     for medmass in medmasslist:
 
         traxx    = []
-#        spec_e   = []
-#        spec_ae  = []
         spec_mu  = []
         spec_amu = []
-#        Nxx_e    = []
-#        Nxx_ae   = []
         Nxx_mu   = []
         Nxx_amu  = []
 
         for lifetime in lifetimelist:
 
             traxxx   = []
-#            spe_e    = []
-#            spe_ae   = []
             spe_mu   = []
             spe_amu  = []
-#            Nxxx_e   = []
-#            Nxxx_ae  = []
             Nxxx_mu  = []
             Nxxx_amu = []
 
@@ -74,27 +58,15 @@
                 if os.path.isfile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_e.dat") and\
  mass == masslist[massnumber]:
 
-#                   print("wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_e.dat " + "exists")
-
- #                   spe_e.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_e.dat"))
-#                    Nxxx_e.append(spe_e[-1].Integral(spe_e[-1].points[0][0],spe_e[-1].points[-1][0]))
-#                    spe_ae.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_ae.dat"))
-#                    Nxxx_ae.append(spe_ae[-1].Integral(spe_ae[-1].points[0][0],spe_ae[-1].points[-1][0]))
                     spe_mu.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_mu.dat"))
                     Nxxx_mu.append(spe_mu[-1].Integral(spe_mu[-1].points[0][0],spe_mu[-1].points[-1][0]))
                     spe_amu.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_amu.dat"))
                     Nxxx_amu.append(spe_amu[-1].Integral(spe_amu[-1].points[0][0],spe_amu[-1].points[-1][0]))
 
                     traxxx.append(ING.H1D.Empty(0.0,200000.0,2000))
-#                    print spe_mu[-1].Integral(spe_mu[-1].points[0][0],spe_mu[-1].points[-1][0])
-#                   print [mass,medmass,lifetime,channel]
 
                 else:
 
-#                   spe_e.append("None")
-#                    Nxxx_e.append("None")
-#                    spe_ae.append("None")
-#                    Nxxx_ae.append("None")
                     spe_mu.append("None")
                     Nxxx_mu.append("None")
                     spe_amu.append("None")
@@ -103,32 +75,20 @@
                     traxxx.append("None")
 
 
-#            spec_e.append(spe_e)
-#            spec_ae.append(spe_ae)
             spec_mu.append(spe_mu)
             spec_amu.append(spe_amu)
-#            Nxx_e.append(Nxxx_e)
-#            Nxx_ae.append(Nxxx_ae)
             Nxx_mu.append(Nxxx_mu)
             Nxx_amu.append(Nxxx_amu)
             traxx.append(traxxx)
 
-#        spect_e.append(spec_e)
-#        spect_ae.append(spec_ae)
         spect_mu.append(spec_mu)
         spect_amu.append(spec_amu)
-#        Nx_e.append(Nxx_e)
-#        Nx_ae.append(Nxx_ae)
         Nx_mu.append(Nxx_mu)
         Nx_amu.append(Nxx_amu)
         trax.append(traxx)
 
-#    spectrum_e.append(spect_e)
-#    spectrum_ae.append(spect_ae)
     spectrum_mu.append(spect_mu)
     spectrum_amu.append(spect_amu)
-#    N_e.append(Nx_e)
-#    N_ae.append(Nx_ae)
     N_mu.append(Nx_mu)
     N_amu.append(Nx_amu)
     track.append(trax)
@@ -143,26 +103,13 @@
 
 hsun.Scale(hsun.nbin*1.0/(hsun.integral))
 
-#track  = [[ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000)],[ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000)],[ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000)],[ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000),ING.H1D.Empty(0.0,50000.0,1000)]]
-#casc = [[ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800)],[ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800)],[ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800)],[ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800),ING.H1D.Empty(-1.0,1.0,800)]]
 
 for filename in os.listdir("/data/ana/analyses/northern_tracks/version-002-p00/"): # loop over MC files
 
-#    print len(os.listdir("/data/ana/analyses/northern_tracks/version-002-p00/"))
-
     if "MC.npy" in filename and "IC86" in filename: #Only take MC files
 
-#       with np.load("/data/ana/analyses/northern_tracks/version-002-p00/"+filename,"r") as infile:
-#       count = 0.0
         infile = np.load("/data/ana/analyses/northern_tracks/version-002-p00/"+filename,"r")
 
-#       for entry in infile:
-
-#           print entry["trueAzi"]
-
-        #print len(infile)
-
-#       for azt,dect,azr,decr,nutype,weight,energy in zip(infile["trueAzi"],infile["trueZen"],infile["azi"],infile["zen"],infile["trueType"],infile["orig_OW"],infile["trueE"]) : #loop over events taking the relevant information
 
         for entry in infile:
 
@@ -180,11 +127,6 @@
 
             if diff > 0.998 and energy <  masslist[k]: # only look at events in a small ROI around the Sun
 
-                #count += 1.0
-
-                #for k in range(len(masslist)):
-		#print "event found"
-
                 for j in range(len(medmasslist)):
 
                     for i in range(len(lifetimelist)):

diff --git a/Intervals.py b/Intervals.py
@@ -2,7 +2,7 @@
 #METAPROJECT icerec/V05-01-00
 
 import numpy as np
-import ingredients as ING
+import src/ING as ING
 from scipy.stats import poisson
 import os