Alexander Ziegler / ctools

Commit 8d74cdabb5363a466b1b96e09d3ba2a503509e5d

Authored by Jürgen Knödlseder
1 parent 62d75f67

Convert tabs to spaces

Inline Side-by-side
Showing 1 changed file with 217 additions and 217 deletions
scripts/cta_root2irf.py
  View file @8d74cda
@@ -38,58 +38,58 @@ import os @@ -38,58 +38,58 @@ import os
38 # Compute Eq. (17) of Li & Ma (1983) # 38 # Compute Eq. (17) of Li & Ma (1983) #
39 # ================================== # 39 # ================================== #
40 def sigma_lima(Non, Noff, alpha=0.2): 40 def sigma_lima(Non, Noff, alpha=0.2):
41 - """  
42 - Compute Eq. (17) of Li & Ma (1983).  
43 -  
44 - Parameters:  
45 - Non - Number of on counts  
46 - Noff - Number of off counts  
47 - Keywords:  
48 - alpha - Ratio of on-to-off exposure  
49 - """  
50 - # Compute sensitivity  
51 - alpha1 = alpha + 1.0  
52 - sum = Non + Noff  
53 - arg1 = Non / sum  
54 - arg2 = Noff / sum  
55 - term1 = Non * math.log((alpha1/alpha)*arg1)  
56 - term2 = Noff * math.log(alpha1*arg2)  
57 - sigma = math.sqrt(2.0 * (term1 + term2))  
58 -  
59 - # Return sensitivity  
60 - return sigma 41 + """
  42 + Compute Eq. (17) of Li & Ma (1983).
  43 +
  44 + Parameters:
  45 + Non - Number of on counts
  46 + Noff - Number of off counts
  47 + Keywords:
  48 + alpha - Ratio of on-to-off exposure
  49 + """
  50 + # Compute sensitivity
  51 + alpha1 = alpha + 1.0
  52 + sum = Non + Noff
  53 + arg1 = Non / sum
  54 + arg2 = Noff / sum
  55 + term1 = Non * math.log((alpha1/alpha)*arg1)
  56 + term2 = Noff * math.log(alpha1*arg2)
  57 + sigma = math.sqrt(2.0 * (term1 + term2))
  58 +
  59 + # Return sensitivity
  60 + return sigma
61 61
62 62
63 # ======================================== # 63 # ======================================== #
64 # Solve Eq. (17) of Li & Ma (1983) for Non # 64 # Solve Eq. (17) of Li & Ma (1983) for Non #
65 # ======================================== # 65 # ======================================== #
66 def Non_lima(sigma, Noff, alpha=0.2): 66 def Non_lima(sigma, Noff, alpha=0.2):
67 - """  
68 - Solve Eq. (17) of Li & Ma (1983) for Non  
69 -  
70 - Parameters:  
71 - sigma - Required significance  
72 - Noff - Number of off counts  
73 - Keywords:  
74 - alpha - Ratio of on-to-off exposure  
75 - """  
76 - # Set initial guess for Non  
77 - S = 25.0  
78 - Non = S + alpha*Noff  
79 -  
80 - # Compute initial sigma estimate  
81 - s = sigma_lima(Non, Noff, alpha=alpha)  
82 - delta = s - sigma  
83 -  
84 - # Iterate  
85 - while abs(delta) > 1.0e-6:  
86 - S *= sigma/s  
87 - Non = S + alpha*Noff  
88 - s = sigma_lima(Non, Noff, alpha=alpha)  
89 - delta = s - sigma  
90 -  
91 - # Return  
92 - return Non 67 + """
  68 + Solve Eq. (17) of Li & Ma (1983) for Non
  69 +
  70 + Parameters:
  71 + sigma - Required significance
  72 + Noff - Number of off counts
  73 + Keywords:
  74 + alpha - Ratio of on-to-off exposure
  75 + """
  76 + # Set initial guess for Non
  77 + S = 25.0
  78 + Non = S + alpha*Noff
  79 +
  80 + # Compute initial sigma estimate
  81 + s = sigma_lima(Non, Noff, alpha=alpha)
  82 + delta = s - sigma
  83 +
  84 + # Iterate
  85 + while abs(delta) > 1.0e-6:
  86 + S *= sigma/s
  87 + Non = S + alpha*Noff
  88 + s = sigma_lima(Non, Noff, alpha=alpha)
  89 + delta = s - sigma
  90 +
  91 + # Return
  92 + return Non
93 93
94 94
95 # =============================================== # 95 # =============================================== #
@@ -97,195 +97,195 @@ def Non_lima(sigma, Noff, alpha=0.2): @@ -97,195 +97,195 @@ def Non_lima(sigma, Noff, alpha=0.2):
97 # =============================================== # 97 # =============================================== #
98 def get_aeff(bgrate, emin, emax, \ 98 def get_aeff(bgrate, emin, emax, \
99 duration=180000.0, alpha=0.2, sigma=5.0, norm=2.79e-7, index=-2.57): 99 duration=180000.0, alpha=0.2, sigma=5.0, norm=2.79e-7, index=-2.57):
100 - """  
101 - Computes effective area assuming a Crab spectrum. This does not  
102 - include any hadron cut efficiencies, and thus is not very useful in a first  
103 - place ...  
104 -  
105 - Parameters:  
106 - bgrate - Background rate (TeV)  
107 - emin - Minimum energy (TeV)  
108 - emax - Maximum energy (TeV)  
109 - Keywords:  
110 - duration - Observing time in seconds (default: 50h)  
111 - alpha - Non/Noff ratio (default: 0.2)  
112 - sigma - Significance (default: 5.0)  
113 - norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7)  
114 - index - Crab spectral index (default: -2.57)  
115 - """  
116 - # Compute number of source events  
117 - B = bgrate*duration  
118 - Noff = B/alpha  
119 - Non = Non_lima(sigma, Noff)  
120 - S = (Non - B)  
121 -  
122 - # Compute Crab flux in specified energy band  
123 - inx = index+1.0  
124 - cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx))  
125 -  
126 - # Compute effective area  
127 - area = S / duration / cu  
128 -  
129 - # Return effective area  
130 - return area 100 + """
  101 + Computes effective area assuming a Crab spectrum. This does not
  102 + include any hadron cut efficiencies, and thus is not very useful in a first
  103 + place ...
  104 +
  105 + Parameters:
  106 + bgrate - Background rate (TeV)
  107 + emin - Minimum energy (TeV)
  108 + emax - Maximum energy (TeV)
  109 + Keywords:
  110 + duration - Observing time in seconds (default: 50h)
  111 + alpha - Non/Noff ratio (default: 0.2)
  112 + sigma - Significance (default: 5.0)
  113 + norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7)
  114 + index - Crab spectral index (default: -2.57)
  115 + """
  116 + # Compute number of source events
  117 + B = bgrate*duration
  118 + Noff = B/alpha
  119 + Non = Non_lima(sigma, Noff)
  120 + S = (Non - B)
  121 +
  122 + # Compute Crab flux in specified energy band
  123 + inx = index+1.0
  124 + cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx))
  125 +
  126 + # Compute effective area
  127 + area = S / duration / cu
  128 +
  129 + # Return effective area
  130 + return area
131 131
132 132
133 # ============== # 133 # ============== #
134 # Generate files # 134 # Generate files #
135 # ============== # 135 # ============== #
136 def make_files(filename, prfname, bkgname): 136 def make_files(filename, prfname, bkgname):
137 - """  
138 - Generate files.  
139 -  
140 - Parameters:  
141 - filename - Performance file name  
142 - prfname - Performance file name  
143 - bkgname - File function filename  
144 - """  
145 - # Open performance file and file function  
146 - file = TFile(filename)  
147 - fprf = open(prfname, "w")  
148 - fbkg = open(bkgname, "w")  
149 -  
150 - # Write header  
151 - fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")  
152 - #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n") 137 + """
  138 + Generate files.
  139 +
  140 + Parameters:
  141 + filename - Performance file name
  142 + prfname - Performance file name
  143 + bkgname - File function filename
  144 + """
  145 + # Open performance file and file function
  146 + file = TFile(filename)
  147 + fprf = open(prfname, "w")
  148 + fbkg = open(bkgname, "w")
  149 +
  150 + # Write header
  151 + fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
  152 + #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n")
153 153
154 - # Get histograms. Konrad's root files do not have an EffectiveArea80  
155 - # histogram, but his effective area is for 80% containment radius  
156 - # anyways, so we simply read this histogram into aeff80.  
157 - sens = TH1F()  
158 - bgrate = TH1F()  
159 - bgratesqdeg = TH1F()  
160 - aeff = TH1F()  
161 - aeff80 = TH1F()  
162 - angres = TH1F()  
163 - angres80 = TH1F()  
164 - eres = TH1F()  
165 - file.GetObject("DiffSens", sens)  
166 - file.GetObject("BGRate", bgrate)  
167 - file.GetObject("BGRatePerSqDeg", bgratesqdeg)  
168 - file.GetObject("EffectiveArea", aeff)  
169 - try:  
170 - file.GetObject("EffectiveArea80", aeff80)  
171 - except:  
172 - file.GetObject("EffectiveArea", aeff80)  
173 - file.GetObject("AngRes", angres)  
174 - file.GetObject("AngRes80", angres80)  
175 - file.GetObject("ERes", eres)  
176 -  
177 - # Setup vectors  
178 - energies = bgrate.GetXaxis()  
179 - nbins_eng = energies.GetNbins()  
180 - for i in range(nbins_eng):  
181 -  
182 - # Get logE and background rate per squared degree  
183 - logE = energies.GetBinCenter(i+1)  
184 - bgd = bgrate.GetBinContent(i+1,1)  
185 - bgdsq = bgratesqdeg.GetBinContent(i+1,1)  
186 - area = aeff.GetBinContent(i+1,1)  
187 - area80 = aeff80.GetBinContent(i+1,1)  
188 - r68 = angres.GetBinContent(i+1,1)  
189 - r80 = angres80.GetBinContent(i+1,1)  
190 - dE = eres.GetBinContent(i+1,1)  
191 - diffs = sens.GetBinContent(i+1,1)  
192 -  
193 - # Skip any NaNs  
194 - if math.isnan(area80):  
195 - continue 154 + # Get histograms. Konrad's root files do not have an EffectiveArea80
  155 + # histogram, but his effective area is for 80% containment radius
  156 + # anyways, so we simply read this histogram into aeff80.
  157 + sens = TH1F()
  158 + bgrate = TH1F()
  159 + bgratesqdeg = TH1F()
  160 + aeff = TH1F()
  161 + aeff80 = TH1F()
  162 + angres = TH1F()
  163 + angres80 = TH1F()
  164 + eres = TH1F()
  165 + file.GetObject("DiffSens", sens)
  166 + file.GetObject("BGRate", bgrate)
  167 + file.GetObject("BGRatePerSqDeg", bgratesqdeg)
  168 + file.GetObject("EffectiveArea", aeff)
  169 + try:
  170 + file.GetObject("EffectiveArea80", aeff80)
  171 + except:
  172 + file.GetObject("EffectiveArea", aeff80)
  173 + file.GetObject("AngRes", angres)
  174 + file.GetObject("AngRes80", angres80)
  175 + file.GetObject("ERes", eres)
  176 +
  177 + # Setup vectors
  178 + energies = bgrate.GetXaxis()
  179 + nbins_eng = energies.GetNbins()
  180 + for i in range(nbins_eng):
  181 +
  182 + # Get logE and background rate per squared degree
  183 + logE = energies.GetBinCenter(i+1)
  184 + bgd = bgrate.GetBinContent(i+1,1)
  185 + bgdsq = bgratesqdeg.GetBinContent(i+1,1)
  186 + area = aeff.GetBinContent(i+1,1)
  187 + area80 = aeff80.GetBinContent(i+1,1)
  188 + r68 = angres.GetBinContent(i+1,1)
  189 + r80 = angres80.GetBinContent(i+1,1)
  190 + dE = eres.GetBinContent(i+1,1)
  191 + diffs = sens.GetBinContent(i+1,1)
  192 +
  193 + # Skip any NaNs
  194 + if math.isnan(area80):
  195 + continue
196 196
197 - # Compute energy (in MeV)  
198 - energy = math.pow(10.0, logE)*1.0e6  
199 - emin = math.pow(10.0, logE-0.1)*1.0e6  
200 - emax = math.pow(10.0, logE+0.1)*1.0e6  
201 - ewidth = emax-emin  
202 -  
203 - # Convert into background rate per steradian and MeV  
204 - bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth  
205 -  
206 - # Compute control background rate (this only works for KB files  
207 - # as they are for 80% containment radius). But we don't use this  
208 - # anyways, it's just for control on display ...  
209 - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))  
210 - if omega > 0.0:  
211 - bkg_rate2 = bgd / omega / ewidth  
212 - else:  
213 - bkg_rate2 = 0.0  
214 -  
215 - # Compute full effective area by dividing area80/0.80  
216 - area = area80/0.80 197 + # Compute energy (in MeV)
  198 + energy = math.pow(10.0, logE)*1.0e6
  199 + emin = math.pow(10.0, logE-0.1)*1.0e6
  200 + emax = math.pow(10.0, logE+0.1)*1.0e6
  201 + ewidth = emax-emin
  202 +
  203 + # Convert into background rate per steradian and MeV
  204 + bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
  205 +
  206 + # Compute control background rate (this only works for KB files
  207 + # as they are for 80% containment radius). But we don't use this
  208 + # anyways, it's just for control on display ...
  209 + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
  210 + if omega > 0.0:
  211 + bkg_rate2 = bgd / omega / ewidth
  212 + else:
  213 + bkg_rate2 = 0.0
  214 +
  215 + # Compute full effective area by dividing area80/0.80
  216 + area = area80/0.80
217 217
218 - # Write results in file  
219 - line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \  
220 - (logE, area, r68, r80, dE, bgd, diffs)  
221 - fprf.write(line+"\n")  
222 - fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") 218 + # Write results in file
  219 + line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \
  220 + (logE, area, r68, r80, dE, bgd, diffs)
  221 + fprf.write(line+"\n")
  222 + fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
223 223
224 - # Show performance file  
225 - #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \  
226 - # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2) 224 + # Show performance file
  225 + #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \
  226 + # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2)
227 227
228 - # Write trailer  
229 - fprf.write("---------------------------------------------\n")  
230 - fprf.write("Notes\n")  
231 - fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n")  
232 - fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")  
233 - fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")  
234 - fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")  
235 - fprf.write(" 5) Fractional Energy Resolution (rms)\n")  
236 - fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")  
237 - fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")  
238 - fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")  
239 - fprf.write(" systematics and statistics and at least 10 photons.\n") 228 + # Write trailer
  229 + fprf.write("---------------------------------------------\n")
  230 + fprf.write("Notes\n")
  231 + fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n")
  232 + fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")
  233 + fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")
  234 + fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")
  235 + fprf.write(" 5) Fractional Energy Resolution (rms)\n")
  236 + fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")
  237 + fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")
  238 + fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")
  239 + fprf.write(" systematics and statistics and at least 10 photons.\n")
240 240
241 - # Close files  
242 - fprf.close()  
243 - fbkg.close() 241 + # Close files
  242 + fprf.close()
  243 + fbkg.close()
244 244
245 - # Return  
246 - return 245 + # Return
  246 + return
247 247
248 248
249 #==========================# 249 #==========================#
250 # Main routine entry point # 250 # Main routine entry point #
251 #==========================# 251 #==========================#
252 if __name__ == '__main__': 252 if __name__ == '__main__':
253 - """  
254 - """  
255 - # Get ROOT filename from command line argument  
256 - if (len(sys.argv) > 1):  
257 - filename = sys.argv[1]  
258 - else:  
259 - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010"  
260 - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root"  
261 - caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" 253 + """
  254 + """
  255 + # Get ROOT filename from command line argument
  256 + if (len(sys.argv) > 1):
  257 + filename = sys.argv[1]
  258 + else:
  259 + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010"
  260 + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root"
  261 + caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root"
262 262
263 - # Get list of all performance files  
264 - files = glob.glob(caldb+"/*.root")  
265 - files.sort()  
266 -  
267 - # Loop over all files  
268 - for file in files:  
269 -  
270 - # Get filename  
271 - head, tail = os.path.split(file)  
272 -  
273 - # Build name  
274 - name = tail.strip(".root")  
275 - if name.find("kb_") != -1:  
276 - name = name.replace("_20deg", "")  
277 - elif name.find("IFAE_") != -1:  
278 - name = name.replace("Subarray", "")  
279 - name = name.replace("IFAE_", "")  
280 - name = "IFAE_"+name  
281 - i = name.find("hours")  
282 - if i != -1:  
283 - name = name[0:i+1]  
284 -  
285 - # Build output filenames  
286 - prfname = name + ".dat"  
287 - bkgname = "bkg_"+name+".txt"  
288 - #print prfname, bkgname 263 + # Get list of all performance files
  264 + files = glob.glob(caldb+"/*.root")
  265 + files.sort()
  266 +
  267 + # Loop over all files
  268 + for file in files:
  269 +
  270 + # Get filename
  271 + head, tail = os.path.split(file)
  272 +
  273 + # Build name
  274 + name = tail.strip(".root")
  275 + if name.find("kb_") != -1:
  276 + name = name.replace("_20deg", "")
  277 + elif name.find("IFAE_") != -1:
  278 + name = name.replace("Subarray", "")
  279 + name = name.replace("IFAE_", "")
  280 + name = "IFAE_"+name
  281 + i = name.find("hours")
  282 + if i != -1:
  283 + name = name[0:i+1]
  284 +
  285 + # Build output filenames
  286 + prfname = name + ".dat"
  287 + bkgname = "bkg_"+name+".txt"
  288 + #print prfname, bkgname
289 289
290 - # Create files  
291 - make_files(file, prfname, bkgname) 290 + # Create files
  291 + make_files(file, prfname, bkgname)