Alexander Ziegler / ctools

Commit ad205fb2a5a936f6a3d86cd7def439b25f138b1a

Authored by knodlseder
1 parent b421dd14

transforms all .root file in the specified directory

Inline Side-by-side
Showing 1 changed file with 289 additions and 0 deletions
scripts/cta_root2irf.py 0 → 100755
  View file @ad205fb
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates ASCII performance files and file functions for
  4 +# background modelling from MC ROOT files. On input it takes the filename
  5 +# of the ROOT performance file while on output it produced 2 ASCII files.
  6 +# The file "performance.dat" contains the performance data in the usual
  7 +# ASCII format, the file "background.txt" contains a file function that can
  8 +# be used within GammaLib for spectral background modelling.
  9 +#
  10 +# Note that the effective area given in the ASCII file is the full effective
  11 +# area (i.e. the area before applying any theta cut). This is required by
  12 +# GammaLib as the maximum likelihood analysis does not imply any theta cut.
  13 +#
  14 +# Copyright (C) 2011 Jurgen Knodlseder
  15 +#
  16 +# This program is free software: you can redistribute it and/or modify
  17 +# it under the terms of the GNU General Public License as published by
  18 +# the Free Software Foundation, either version 3 of the License, or
  19 +# (at your option) any later version.
  20 +#
  21 +# This program is distributed in the hope that it will be useful,
  22 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  23 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  24 +# GNU General Public License for more details.
  25 +#
  26 +# You should have received a copy of the GNU General Public License
  27 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  28 +#
  29 +# ==========================================================================
  30 +from ROOT import TFile, TH1F, TH2F
  31 +import math
  32 +import sys
  33 +import glob
  34 +import os
  35 +
  36 +
  37 +# ================================== #
  38 +# Compute Eq. (17) of Li & Ma (1983) #
  39 +# ================================== #
  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
  61 +
  62 +
  63 +# ======================================== #
  64 +# Solve Eq. (17) of Li & Ma (1983) for Non #
  65 +# ======================================== #
  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
  93 +
  94 +
  95 +# =============================================== #
  96 +# Compute effective area assuming a Crab spectrum #
  97 +# =============================================== #
  98 +def get_aeff(bgrate, emin, emax, \
  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
  131 +
  132 +
  133 +# ============== #
  134 +# Generate files #
  135 +# ============== #
  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 + print "log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)"
  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
  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
  217 +
  218 + # Write results in file
  219 + line = "%.1f %.1f %.4f %.4f %.4f %.7f %.5e" % \
  220 + (logE, area, r68, r80, dE, bgd, diffs)
  221 + fprf.write(line+"\n")
  222 + fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
  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)
  227 +
  228 + # Write trailer
  229 + fprf.write("---------------------------------------------\n")
  230 + fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n")
  231 + fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")
  232 + fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")
  233 + fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")
  234 + fprf.write(" 5) Fractional Energy Resolution (rms)\n")
  235 + fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")
  236 + fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")
  237 + fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")
  238 + fprf.write(" systematics and statistics and at least 10 photons.\n")
  239 +
  240 + # Close files
  241 + fprf.close()
  242 + fbkg.close()
  243 +
  244 + # Return
  245 + return
  246 +
  247 +
  248 +#==========================#
  249 +# Main routine entry point #
  250 +#==========================#
  251 +if __name__ == '__main__':
  252 + """
  253 + """
  254 + # Get ROOT filename from command line argument
  255 + if (len(sys.argv) > 1):
  256 + filename = sys.argv[1]
  257 + else:
  258 + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010"
  259 + caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_20deg_root"
  260 +
  261 + # Get list of all performance files
  262 + files = glob.glob(caldb+"/*.root")
  263 + files.sort()
  264 +
  265 + # Loop over all files
  266 + for file in files:
  267 +
  268 + # Get filename
  269 + head, tail = os.path.split(file)
  270 +
  271 + # Build name
  272 + name = tail.strip(".root")
  273 + if name.find("kb_") != -1:
  274 + name = name.replace("_20deg", "")
  275 + elif name.find("IFAE_") != -1:
  276 + name = name.replace("Subarray", "")
  277 + name = name.replace("IFAE_", "")
  278 + name = "IFAE_"+name
  279 + i = name.find("hours")
  280 + if i != -1:
  281 + name = name[0:i+1]
  282 +
  283 + # Build output filenames
  284 + prfname = name + ".dat"
  285 + bkgname = "bkg_"+name+".txt"
  286 + print prfname, bkgname
  287 +
  288 + # Create files
  289 + make_files(file, prfname, bkgname)
... ...