#! /usr/bin/env python # ========================================================================== # This script performs the ctools science verification. It creates and # analyses the pull distributions for a variety of models in unbinned # analysis mode. At the end the script produces a JUnit compliant # science verification report. # # Usage: # ./science_verification.py # # -------------------------------------------------------------------------- # # Copyright (C) 2015 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # ========================================================================== import gammalib import ctools import cscripts import os import csv import math import sys # ========================== # # Generate pull distribution # # ========================== # def generate_pull_distribution(model, trials=100, \ caldb="prod2", irf="South_50h", \ deadc=0.95, edisp=False, \ ra=83.63, dec=22.01, rad=5.0, \ emin=0.1, emax=100.0, enumbins=0, \ duration=1800.0, \ npix=200, binsz=0.02, \ coordsys="CEL", proj="TAN", \ pattern="single", offset=1.5, \ debug=False, chatter=2): """ Generates pull distribution for a given model. Parameters: model - Model XML filename (without .xml extension) """ # Derive parameters head, tail = os.path.split(model) inmodel = model + ".xml" outfile = "cspull_" + tail + ".dat" # Setup pull distribution generation pull = cscripts.cspull() pull["inobs"] = "NONE" pull["inmodel"] = inmodel pull["outfile"] = outfile pull["caldb"] = caldb pull["irf"] = irf pull["edisp"] = edisp pull["ra"] = ra pull["dec"] = dec pull["rad"] = rad pull["emin"] = emin pull["emax"] = emax pull["tmin"] = 0.0 pull["tmax"] = duration pull["enumbins"] = enumbins pull["npix"] = npix pull["binsz"] = binsz pull["coordsys"] = coordsys pull["proj"] = proj pull["deadc"] = deadc pull["rad"] = rad pull["pattern"] = pattern pull["offset"] = offset pull["ntrials"] = trials pull["debug"] = debug pull["chatter"] = chatter # Generate pull distributions pull.run() # Return return outfile # ========================= # # Analyse pull distribution # # ========================= # def analyse_pull_distribution(filename): """ Compute mean and standard deviation of pull distribution. Parameter: filename - Pull distribution ASCII file to analyse. """ # Initialise column names, means and standard deviations colnames = [] means = [] stds = [] # Open reader reader = csv.reader(open(filename, 'r'), delimiter=',') # Read rows first = True index = -1 samples = 0.0 for row in reader: # Get column names if first row if first: for element in row: colnames.append(element) means.append(0.0) stds.append(0.0) # Handle data rows else: for i, element in enumerate(row): means[i] += float(element) stds[i] += float(element)*float(element) samples += 1.0 # Flag that first row has been passed first = False # Compute mean and standard deviations for i in range(len(stds)): std = math.sqrt(stds[i]/samples - means[i]*means[i]/(samples*samples)) stds[i] = std means[i] /= samples # Setup results results = {} for i in range(len(colnames)): results[colnames[i]] = {'mean': means[i], 'std': stds[i]} # Return results return results # =================================== # # Test class for science verification # # =================================== # class sciver(gammalib.GPythonTestSuite): """ Test class for science verification """ # Constructor def __init__(self): """ Constructor. """ # Call base class constructor gammalib.GPythonTestSuite.__init__(self) # Initialise results self.results = None # Return return # Set test functions def set(self): """ Set all test functions. """ # Set test name self.name("Science Verification") # Append spectral tests self.append(self.spec_plaw, "Test power law model") self.append(self.spec_plaw2, "Test power law 2 model") self.append(self.spec_eplaw, "Test exponentially cut off power law model") self.append(self.spec_supeplaw, "Test super exponentially cut off power law model") self.append(self.spec_logparabola, "Test log parabola model") self.append(self.spec_gauss, "Test Gaussian model") self.append(self.spec_filefct, "Test file function model") self.append(self.spec_nodes, "Test nodes model") # Append spatial tests self.append(self.spat_ptsrc, "Test point source model") self.append(self.spat_rdisk, "Test radial disk model") self.append(self.spat_rgauss, "Test radial Gaussian model") self.append(self.spat_rshell, "Test radial shell model") self.append(self.spat_edisk, "Test elliptical disk model") self.append(self.spat_egauss, "Test elliptical Gaussian model") self.append(self.spat_const, "Test diffuse isotropic model") self.append(self.spat_map, "Test diffuse map model") self.append(self.spat_map_roi, "Test diffuse map model (small ROI)") self.append(self.spat_map_nn, "Test diffuse map model (not normalized and scaled)") self.append(self.spat_cube, "Test diffuse cube model") # Return return # Generate and analyse pull distributions def pull(self, model, trials=100, ra=83.63, dec=22.01, rad=5.0): """ Generate and analyse pull distributions. """ # Generate pull distribution outfile = generate_pull_distribution(model, trials=trials, ra=ra, dec=dec, rad=rad) # Analyse pull distribution self.results = analyse_pull_distribution(outfile) # Return return # Test parameter result def test(self, name): """ Test one parameter. """ # Test mean mean = self.results[name]["mean"] valid = (mean >= -0.40) and (mean <= +0.40) text = "Mean %.5f of %s should be within [-0.4,0.4] range" % (mean, name) self.test_assert(valid, text) # Test standard deviation std = self.results[name]["std"] valid = (std >= 0.80) and (std <= 1.20) text = "Standard deviation %.5f of %s should be within [0.8,1.2] range" % (std, name) self.test_assert(valid, text) # Return return # Test power law model def spec_plaw(self): """ Test power law model. """ self.pull("data/sciver/crab_plaw") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test power law 2 model def spec_plaw2(self): """ Test power law 2 model. """ self.pull("data/sciver/crab_plaw2") self.test("Pull_Crab_Integral") self.test("Pull_Crab_Index") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test exponentially cut off power law model def spec_eplaw(self): """ Test exponentially cut off power law model. """ self.pull("data/sciver/crab_eplaw") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_Cutoff") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test super exponentially cut off power law model def spec_supeplaw(self): """ Test super exponentially cut off power law model. """ self.pull("data/sciver/crab_supeplaw") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index1") self.test("Pull_Crab_Index2") self.test("Pull_Crab_Cutoff") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test log parabola model def spec_logparabola(self): """ Test log parabola model. """ self.pull("data/sciver/crab_logparabola") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_Curvature") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test Gaussian model def spec_gauss(self): """ Test Gaussian model. """ self.pull("data/sciver/crab_gauss") self.test("Pull_Crab_Normalization") self.test("Pull_Crab_Mean") self.test("Pull_Crab_Sigma") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test file function model def spec_filefct(self): """ Test file function model. """ self.pull("data/sciver/crab_filefct") self.test("Pull_Crab_Normalization") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test nodes model def spec_nodes(self): """ Test nodes model. """ self.pull("data/sciver/crab_nodes") self.test("Pull_Crab_Intensity0") self.test("Pull_Crab_Intensity1") self.test("Pull_Crab_Intensity2") self.test("Pull_Crab_Intensity3") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test point source model def spat_ptsrc(self): """ Test point source model. """ self.pull("data/sciver/crab_ptsrc") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_RA") self.test("Pull_Crab_DEC") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test radial disk model def spat_rdisk(self): """ Test radial disk model. """ self.pull("data/sciver/crab_rdisk") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_RA") self.test("Pull_Crab_DEC") self.test("Pull_Crab_Radius") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test radial Gaussian model def spat_rgauss(self): """ Test radial Gaussian model. """ self.pull("data/sciver/crab_rgauss") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_RA") self.test("Pull_Crab_DEC") self.test("Pull_Crab_Sigma") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test radial shell model def spat_rshell(self): """ Test radial shell model. """ self.pull("data/sciver/crab_rshell") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_RA") self.test("Pull_Crab_DEC") self.test("Pull_Crab_Radius") self.test("Pull_Crab_Width") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test elliptical disk model def spat_edisk(self): """ Test elliptical disk model. """ self.pull("data/sciver/crab_edisk") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_RA") self.test("Pull_Crab_DEC") self.test("Pull_Crab_PA") self.test("Pull_Crab_MinorRadius") self.test("Pull_Crab_MajorRadius") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test elliptical Gaussian model def spat_egauss(self): """ Test elliptical Gaussian model. """ self.pull("data/sciver/crab_egauss") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Crab_RA") self.test("Pull_Crab_DEC") self.test("Pull_Crab_PA") self.test("Pull_Crab_MinorRadius") self.test("Pull_Crab_MajorRadius") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test diffuse isotropic model def spat_const(self): """ Test diffuse isotropic model. """ self.pull("data/sciver/crab_const") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test diffuse map model def spat_map(self): """ Test diffuse map model. """ self.pull("data/sciver/crab_map", ra=201.3651, dec=-43.0191) self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test diffuse map model (small ROI) def spat_map_roi(self): """ Test diffuse map model (small ROI). """ self.pull("data/sciver/crab_map_roi", ra=201.3651, dec=-43.0191, rad=1.5) self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test diffuse map model (not normalized and scaled) def spat_map_nn(self): """ Test diffuse map model (not normalized and scaled). """ self.pull("data/sciver/crab_map_nn", ra=201.3651, dec=-43.0191, rad=1.5) self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # Test diffuse cube model def spat_cube(self): """ Test diffuse cube model. """ self.pull("data/sciver/crab_cube") self.test("Pull_Crab_Prefactor") self.test("Pull_Crab_Index") self.test("Pull_Background_Prefactor") self.test("Pull_Background_Index") return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': """ Perform ctools science verification """ # Allocate test suite container suites = gammalib.GTestSuites("ctools science verification") # Allocate test suite and append it to the container suite_sciver = sciver() # Setup test suit suite_sciver.set() # Append test suite to container suites.append(suite_sciver) # Set PFILES environment variable try: os.mkdir("pfiles") except: pass os.system("cp -r ../src/*/*.par pfiles/") os.environ['PFILES'] = "pfiles" # Run test suite success = suites.run() # Save test results suites.save("reports/sciver.xml") # Set return code if success: rc = 0 else: rc = 1 # Exit with return code sys.exit(rc)