diff --git a/cscripts/cscaldb.py b/cscripts/cscaldb.py index 018b568..dbdfeb6 100755 --- a/cscripts/cscaldb.py +++ b/cscripts/cscaldb.py @@ -42,7 +42,7 @@ class cscaldb(ctools.cscript): # Set name self.name = "cscaldb" self.version = "1.0.0" - + # Make sure that parfile exists file = self.parfile() @@ -60,7 +60,7 @@ class cscaldb(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -75,35 +75,35 @@ class cscaldb(ctools.cscript): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append_standard() pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile. """ # Return return - + def execute(self): """ Execute the script. """ # Run the script self.run() - + # Return return @@ -117,7 +117,7 @@ class cscaldb(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() @@ -156,7 +156,7 @@ class cscaldb(ctools.cscript): # Loop over instruments for instrument in instruments: - + # Write mission into logger if self.logTerse(): self.log.header3("Response functions in database \""+instrument+"\"") @@ -166,7 +166,7 @@ class cscaldb(ctools.cscript): fits = gammalib.GFits(cifname) cif = fits["CIF"] cals = cif["CAL_CBD"] - + # Extract response names names = [] nrows = cals.length() @@ -181,13 +181,13 @@ class cscaldb(ctools.cscript): if names.count(name) == 0: names.append(name) names.sort() - + # Print response name if self.logTerse(): for name in names: self.log(name+"\n") self.log("\n") - + # Return return @@ -201,10 +201,9 @@ if __name__ == '__main__': """ # Create instance of application app = cscaldb(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - \ No newline at end of file diff --git a/cscripts/cshessobs.py b/cscripts/cshessobs.py index 5ddca56..09bf42c 100755 --- a/cscripts/cshessobs.py +++ b/cscripts/cshessobs.py @@ -52,12 +52,12 @@ class cshessobs(ctools.cscript): except KeyError: raise RuntimeError("Environment variable \"HESSFITS\" not set. " \ "Please set the environment variable \"HESSFITS\" to use this script") - + # Initialise some members self.inmodels = None self.outobs = "obs.xml" self.nodes = 1 - + # Make sure that parfile exists file = self.parfile() @@ -75,7 +75,7 @@ class cshessobs(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -88,64 +88,64 @@ class cshessobs(ctools.cscript): Check if parfile exists. If parfile does not exist then create a default parfile. This kluge avoids shipping the cscript with a parfile. """ - + # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal that parfile was not found print("Parfile \""+parfile+"\" not found. Create default parfile.") - + # Check for folders and possible/available configurations chains = os.listdir(self.datapath) - + # Initialise folder dictionary folders = {} - + # Loop over possible analysis chains for chain in chains: - + # Omit hidden files/folders if chain[0]=="." or not os.path.isdir(self.datapath+"/"+chain): continue - + # Initialise dictionary per chain folders[chain] = {} - + # Read available DST versions dstversions = os.listdir(self.datapath+"/"+chain) - + # Loop over DST versions for dstver in dstversions: - + # Omit hidden files/folders if dstver[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver): continue - + # Initialise available cut configurations folders[chain][dstver] = [] configs = os.listdir(self.datapath+"/"+chain+"/"+dstver) - + # Loop over configurations for config in configs: - + # Omit hidden files/folders if config[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver+"/"+config): continue - + # Append configuration to corresponding dictionary folders[chain][dstver].append(config) - + # Check for appropriate folder structure checked_folders = self.check_structure(folders) if not len(checked_folders): raise KeyError("Appropriate folder structure not found in $HESSFITS. The correct structure can be found on HESS internal documentation") - + # Create default parfile pars = gammalib.GApplicationPars() - + # Create parameters for each available config avchains = "" for chain in checked_folders: @@ -168,38 +168,38 @@ class cshessobs(ctools.cscript): pars.append(gammalib.GApplicationPar("usetrig","b","h","yes","yes|no","","Pre-estimate background normalisation? Recommended for binned analysis")) pars.append_standard() pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Get parameters self.runlistfile = self["runlistfile"].filename() - + # Read models if user input self.inmodel = self["inmodel"].filename() if not (self.inmodel == "NONE" or self.inmodel == ""): self.inmodels = gammalib.GModels(self.inmodel) - + # Observation outfile self.outobs = self["outobs"].filename() - + # Output model file self.outmodel = self["outmodel"].filename() - + # Read hidden parameters self.nodes = self["nodes"].integer() self.usetrig = self["usetrig"].boolean() self.datapath = self["datapath"].string() - + # Read run configuration parameters self.chain = self["chain"].string() self.dstver = self[self.chain+"_dstver"].string() self.config = self[self.chain+"_"+self.dstver+"_config"].string() - + # Build directory name prefix self.directory = self.datapath+"/" self.directory += self.chain +"/" @@ -208,22 +208,22 @@ class cshessobs(ctools.cscript): if not os.path.isdir(self.directory): raise gammalib.Exception.runtime_error("Folder structure not found. Please check your HESSFITS environment variable") - + self.m_log = False # Logging in client tools self.m_debug = False # Debugging in client tools self.m_clobber = self["clobber"].boolean() - + # Return return - + def check_structure(self,folders): - + # Initialise valid folders newfolders = {} - + # Check if we have found possible chains if len(folders): - + # Loop over chains for chain in folders: @@ -232,7 +232,7 @@ class cshessobs(ctools.cscript): # Loop over configs: for config in folders[chain][dst]: - + # check if config is valid path = self.datapath+"/"+chain+"/"+dst+"/"+config if os.path.isdir(path) and len(glob.glob(path+"/run*")): @@ -241,25 +241,25 @@ class cshessobs(ctools.cscript): if not newfolders[chain].has_key(dst): newfolders[chain][dst] = [] newfolders[chain][dst].append(config) - + # Return valid configs return newfolders - + # Create instrument background model def hess_inst_background(self,run,aefffile,trgrate=0.0): - + # Initialise return spectrum spec = 0 - + # Initialise default value value = 1.0 - + if self.usetrig and trgrate > 0.0: # Use value from lookup value = 0.00451521513379 * trgrate + 0.162316854264 - + if self.nodes <= 1: - + # Use a constant scaling function spec = gammalib.GModelSpectralConst() spec["Value"].min(0.1) @@ -269,9 +269,9 @@ class cshessobs(ctools.cscript): spec["Value"].fix() else: spec["Value"].free() - + elif self.nodes == 2: - + # User a power law scaling function e = gammalib.GEnergy(1.0,"TeV") spec = gammalib.GModelSpectralPlaw(1.0,0.0,e) @@ -285,7 +285,7 @@ class cshessobs(ctools.cscript): spec[1].max(2.0) spec[1].scale(1.0) else: - + # Use several energy-dependent scaling factors fits = gammalib.GFits(aefffile) emin = fits["EFFECTIVE AREA"].real("LO_THRES") @@ -301,39 +301,39 @@ class cshessobs(ctools.cscript): elif "Intensity" in par.name(): par.min(0.01) par.max(10.0) - + # Create background model instance bck = gammalib.GCTAModelIrfBackground(spec) - + # Assign specific run id bck.ids(str(int(run))) - + # Assign instrument bck.instruments("HESS") - + # Set name (arbitrary) bck.name("bkg_"+run) - + # Turn off TS calculation for background model bck.tscalc(False) - + # Return model return bck - - + + def execute(self): """ Execute the script. """ # Run the script self.run() - + # Save residual map self.save(self.outobs, self.m_clobber) - + # Return return - + # Get runfolder from runnumber (subject to defined folder structure) def get_runfolder(self,runnr): @@ -346,13 +346,13 @@ class cshessobs(ctools.cscript): end = str(int_end) if len(end) == 5: end = "0"+ end - + runfold = runnr if len(runfold) == 5: runfold = "0"+runfold - + return self.directory+"/run"+start+"-"+end+"/run"+runfold+"/" - + def run(self): """ Run the script. @@ -360,7 +360,7 @@ class cshessobs(ctools.cscript): # Switch screen logging on in debug mode if self.logDebug(): self.log.cout(True) - + # Get parameters self.get_parameters() @@ -368,7 +368,7 @@ class cshessobs(ctools.cscript): if self.logTerse(): self.log_parameters() self.log("\n") - + # Initialise empty models self.models = gammalib.GModels() @@ -376,7 +376,7 @@ class cshessobs(ctools.cscript): self.xml = gammalib.GXml() self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\"")) lib = self.xml.element("observation_list", 0) - + # Read runlist from file runlist = [] runfile = open(self.runlistfile) @@ -390,30 +390,30 @@ class cshessobs(ctools.cscript): else: runlist.append('') runfile.close() - + # Loop over runs # Log output if self.logTerse(): self.log("\n") self.log.header1("Looping over runs") - + for runnr in runlist: - + # Unify run number string length run = runnr if len(run) == 5: run = "0"+run - + # Get folder in which run files should be located runfolder = self.get_runfolder(run) - + # construct filenames eventfile = runfolder + "hess_events_"+run+".fits.gz" aefffile = runfolder + "hess_aeff_2d_"+run+".fits.gz" psffile = runfolder + "hess_psf_king_"+run+".fits.gz" edispfile = ""#runfolder + "hess_edisp_2d_"+run+".fits.gz" bgfile = runfolder + "hess_bkg_offruns_"+run+".fits.gz" - + # Check for existence of files msg = "" skip = False @@ -434,7 +434,7 @@ class cshessobs(ctools.cscript): skip = True else: msg = "Adding run "+str(int(run)) - + # Log ouput and skip run if not available if self.logTerse(): self.log(msg+"\n") @@ -446,11 +446,11 @@ class cshessobs(ctools.cscript): self.log(" Point spread function: "+psffile+"\n") self.log(" Energy dispersion: "+edispfile+"\n") self.log(" Background file: "+bgfile+"\n") - + # Open fits file to determine the observation name fits = gammalib.GFits(eventfile) object_name = fits["EVENTS"].string("OBJECT") - + # Retrieve trigger rate if available trgrate = 0.0 if fits["EVENTS"].has_card("ZTRGRATE"): @@ -458,7 +458,7 @@ class cshessobs(ctools.cscript): # Close FITS file fits.close() - + # Append observation to XML and set attributes obs = lib.append("observation"); obs.attribute("name", object_name); @@ -469,29 +469,29 @@ class cshessobs(ctools.cscript): par = gammalib.GXmlElement("parameter name=\"EventList\"") par.attribute("file", eventfile); obs.append(par) - + # Append effective area aeff = gammalib.GXmlElement("parameter name=\"EffectiveArea\"") aeff.attribute("file",aefffile) - + # Append PSF psf = gammalib.GXmlElement("parameter name=\"PointSpreadFunction\"") psf.attribute("file",psffile) - + # Append energy dispersion edisp = gammalib.GXmlElement("parameter name=\"EnergyDispersion\"") edisp.attribute("file",edispfile) - + # Append background bck = gammalib.GXmlElement("parameter name=\"Background\"") bck.attribute("file",bgfile) - + # assign events and IRFs to observation obs.append(aeff) obs.append(psf) obs.append(edisp) obs.append(bck) - + # Append instrumental background model self.models.append(self.hess_inst_background(run,aefffile,trgrate)) @@ -500,7 +500,7 @@ class cshessobs(ctools.cscript): if self.logTerse(): self.log("\n") self.log.header1("Appending models") - + # Loop over input models for model in self.inmodels: if self.logTerse(): @@ -509,12 +509,12 @@ class cshessobs(ctools.cscript): # Return return - + def save(self,outfile,clobber): - + # Save observation XML file self.xml.save(outfile) - + # Save model XML file self.models.save(self.outmodel) @@ -528,10 +528,9 @@ if __name__ == '__main__': """ # Create instance of application app = cshessobs(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index 6252230..94000a3 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.py @@ -36,14 +36,14 @@ class cslightcrv(ctools.cscript): """ Constructor. """ - + # Set name self.name = "cslightcrv" self.version = "1.0.0" - + # Initialise some members self.obs = None - + # Initialise some members if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): self.obs = argv[0] @@ -52,7 +52,7 @@ class cslightcrv(ctools.cscript): self.obs = gammalib.GObservations() self.obs.clear() self.m_outfile = "" - + # Make sure that parfile exists file = self.parfile() @@ -70,7 +70,7 @@ class cslightcrv(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -78,7 +78,7 @@ class cslightcrv(ctools.cscript): # Write separator into logger if self.logTerse(): self.log("\n") - + # Return return @@ -89,13 +89,13 @@ class cslightcrv(ctools.cscript): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file")) @@ -129,12 +129,12 @@ class cslightcrv(ctools.cscript): pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","cslightcrv.log","","","Log filename")) pars.save(parfile) - + # Return return def create_tbounds(self): - + # Create time bin container self.m_tbins = gammalib.GGti() @@ -172,7 +172,7 @@ class cslightcrv(ctools.cscript): # Return return - + def get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -181,14 +181,14 @@ class cslightcrv(ctools.cscript): if self.obs == None or self.obs.size() == 0: self.require_inobs("cslightcrv::get_parameters()") self.obs = self.get_observations() - + # Set models if we have none if self.obs.models().size() == 0: self.obs.models(self["inmodel"].filename()) - + # Get source name self.m_srcname = self["srcname"].string() - + # Get time boundaries self.create_tbounds() @@ -215,32 +215,32 @@ class cslightcrv(ctools.cscript): # Read other parameters self.m_outfile = self["outfile"].filename() - + # Get other parameeters self.m_calc_ulimit = self["calc_ulim"].boolean() self.m_calc_ts = self["calc_ts"].boolean() self.m_fix_bkg = self["fix_bkg"].boolean() self.m_fix_srcs = self["fix_srcs"].boolean() - + # Set some fixed parameters self.m_log = False # Logging in client tools self.m_chatter = self["chatter"].integer() self.m_clobber = self["clobber"].boolean() self.m_debug = self["debug"].boolean() - + # Return return - + def models(self, models): """ Set model. """ # Copy models self.obs.models(models.clone()) - + # Return return - + def execute(self): """ Execute the script. @@ -250,7 +250,7 @@ class cslightcrv(ctools.cscript): # Save residual map self.fits.saveto(self.m_outfile, self.m_clobber) - + # Return return @@ -264,12 +264,12 @@ class cslightcrv(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() self.log("\n") - + # Write observation into logger if self.logTerse(): self.log("\n") @@ -284,51 +284,51 @@ class cslightcrv(ctools.cscript): # Adjust model parameters dependent on input user parameters for model in self.obs.models(): - + # Set TS flag for all models to false. # Source of interest will be set to true later model.tscalc(False) - + # Log model name if self.logExplicit(): self.log.header3(model.name()) - + # Deal with the source of interest if model.name() == self.m_srcname: if self.m_calc_ts: model.tscalc(True) - + elif self.m_fix_bkg and not model.classname() == "GModelSky": for par in model: if par.is_free() and self.logExplicit(): self.log(" Fixing \""+par.name()+"\"\n") par.fix() - + elif self.m_fix_srcs and model.classname() == "GModelSky": for par in model: if par.is_free() and self.logExplicit(): self.log(" Fixing \""+par.name()+"\"\n") par.fix() - + # Write header if self.logTerse(): self.log("\n") self.log.header1("Generate lightcurve") - + # Initialise FITS Table with extension "LIGHTCURVE" table = gammalib.GFitsBinTable(self.m_tbins.size()) table.extname("LIGHTCURVE") - + # Add Header for compatibility with gammalib.GMWLSpectrum table.card("INSTRUME", "CTA", "Name of Instrument") table.card("TELESCOP", "CTA", "Name of Telescope") - + # Create FITS table columns MJD = gammalib.GFitsTableDoubleCol("MJD", self.m_tbins.size()) MJD.unit("days") e_MJD = gammalib.GFitsTableDoubleCol("e_MJD", self.m_tbins.size()) e_MJD.unit("days") - + # Create a FITS column for every free parameter columns = [] for par in self.obs.models()[self.m_srcname]: @@ -339,7 +339,7 @@ class cslightcrv(ctools.cscript): e_col = gammalib.GFitsTableDoubleCol("e_"+par.name(), self.m_tbins.size()) e_col.unit(par.unit()) columns.append(e_col) - + # Create TS and upper limit columns TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_tbins.size()) ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_tbins.size()) @@ -356,7 +356,7 @@ class cslightcrv(ctools.cscript): # Get time boundaries tmin = self.m_tbins.tstart(i) tmax = self.m_tbins.tstop(i) - + # Compute time bin center and time width tmean = (tmin + tmax) tmean *= 0.5 @@ -366,11 +366,11 @@ class cslightcrv(ctools.cscript): # Store time as MJD MJD[i] = tmean.mjd() e_MJD[i] = twidth.days() - + # Log information if self.logExplicit(): self.log.header3("Selecting events") - + # Select events select = ctools.ctselect(self.obs) select["emin"] = self.m_emin @@ -384,14 +384,14 @@ class cslightcrv(ctools.cscript): # Retrieve observation obs = select.obs() - + # Binned analysis if self.m_binned: # Header if self.logTerse(): self.log.header3("Binning events") - + # Bin events bin = ctools.ctbin(select.obs()) bin["usepnt"] = False @@ -407,11 +407,11 @@ class cslightcrv(ctools.cscript): bin["coordsys"] = self.m_coordsys bin["proj"] = self.m_proj bin.run() - + # Header if self.logTerse(): self.log.header3("Creating exposure cube") - + # Create exposure cube expcube = ctools.ctexpcube(select.obs()) expcube["incube"] = "NONE" @@ -428,11 +428,11 @@ class cslightcrv(ctools.cscript): expcube["coordsys"] = self.m_coordsys expcube["proj"] = self.m_proj expcube.run() - + # Header if self.logTerse(): self.log.header3("Creating PSF cube") - + # Create psf cube psfcube = ctools.ctpsfcube(select.obs()) psfcube["incube"] = "NONE" @@ -449,11 +449,11 @@ class cslightcrv(ctools.cscript): psfcube["coordsys"] = self.m_coordsys psfcube["proj"] = self.m_proj psfcube.run() - + # Header if self.logTerse(): self.log.header3("Creating background cube") - + # Create background cube bkgcube = ctools.ctbkgcube(select.obs()) bkgcube["incube"] = "NONE" @@ -470,37 +470,37 @@ class cslightcrv(ctools.cscript): bkgcube["coordsys"] = self.m_coordsys bkgcube["proj"] = self.m_proj bkgcube.run() - + # Set new binned observation obs = bin.obs() - + # Set precomputed binned response obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube()) # Get new models models = bkgcube.models() - + # Fix background models if required if self.m_fix_bkg: for model in models: if not model.classname() == "GModelSky": for par in model: par.fix() - + # Set new models to binned observation obs.models(models) - + # Header if self.logTerse(): self.log.header3("Performing fit") - + # Likelihood like = ctools.ctlike(obs) like.run() - + # Skip bin if no event was present if like.obs().logL() == 0.0: - + # Log information if self.logTerse(): self.log("No event in this time bin. Bin is skipped\n") @@ -511,7 +511,7 @@ class cslightcrv(ctools.cscript): TSvalues[i] = 0.0 ulim_values[i] = 0.0 continue - + # Get results fitted_models = like.obs().models() source = fitted_models[self.m_srcname] @@ -519,16 +519,16 @@ class cslightcrv(ctools.cscript): # Calculate Upper Limit ulimit_value = -1.0 if self.m_calc_ulimit: - + # Logging information if self.logTerse(): self.log.header3("Computing upper limit") - + # Create upper limit object ulimit = ctools.ctulimit(like.obs()) ulimit["srcname"] = self.m_srcname ulimit["eref"] = 1.0 - + # Try to run upper limit and catch exceptions try: ulimit.run() @@ -537,26 +537,26 @@ class cslightcrv(ctools.cscript): if self.logTerse(): self.log("Upper limit calculation failed\n") ulimit_value = -1.0 - + # Get TS value TS = -1.0 if self.m_calc_ts: TS = source.ts() - + # Set values for storage TSvalues[i] = TS - + # Set FITS column values for col in columns: if "e_" == col.name()[:2]: col[i] = source.spectral()[col.name()[2:]].error() else: col[i] = source.spectral()[col.name()].value() - + # Store upper limit value if available if ulimit_value > 0.0: ulim_values[i] = ulimit_value - + # Log information if self.logExplicit(): self.log.header3("Results of bin "+str(i)+": MJD "+str(tmin.mjd())+"-"+str(tmax.mjd())) @@ -580,14 +580,14 @@ class cslightcrv(ctools.cscript): table.append(col) table.append(TSvalues) table.append(ulim_values) - + # Create the FITS file now self.fits = gammalib.GFits() self.fits.append(table) - + # Return return - + # ======================== # # Main routine entry point # @@ -598,9 +598,9 @@ if __name__ == '__main__': """ # Create instance of application app = cslightcrv(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() diff --git a/cscripts/csobsdef.py b/cscripts/csobsdef.py index 69c919e..01d8ad7 100755 --- a/cscripts/csobsdef.py +++ b/cscripts/csobsdef.py @@ -33,7 +33,7 @@ class csobsdef(ctools.cscript): file with header keywords in the first row followed by a list of pointings (one pointing per row). The following header keywords are supported (case sensitive, column order irrelevant): - + name - Observation name string id - Unique observation identifier string ra - Right Ascension of pointing (deg) @@ -47,14 +47,14 @@ class csobsdef(ctools.cscript): deadc - Deadtime correction factor [0-1] caldb - Calibration database irf - Response function name - + Only the pairs (ra,dec) or (lon,lat) are mandatory header keywords. All other keywords are optional and can be specified when calling csobsdef as user parameters. The only exception is the "duration" keyword that will automatically queried. - + Here some usage examples: - + csobsdef Creates minimal observation definition file. @@ -85,7 +85,7 @@ class csobsdef(ctools.cscript): self.clobber = True self.debug = False self.read_ahead = False - + # Make sure that parfile exists file = self.parfile() @@ -103,7 +103,7 @@ class csobsdef(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -125,7 +125,7 @@ class csobsdef(ctools.cscript): except: # Signal that parfile was not found print("Parfile \""+parfile+"\" not found. Create default parfile.") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inpnt","f","a","NONE","","","Pointing definition file")) @@ -140,10 +140,10 @@ class csobsdef(ctools.cscript): pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","csobsdef.log","","","Log filename")) pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile. @@ -162,20 +162,20 @@ class csobsdef(ctools.cscript): # Return return - + def execute(self): """ Execute the script. """ # Set read ahead flag self.read_ahead = True - + # Run the script self.run() # Save observation definition file self.obs.save(self.outfile) - + # Return return @@ -189,7 +189,7 @@ class csobsdef(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() @@ -325,13 +325,13 @@ class csobsdef(ctools.cscript): # Create an empty event list list = gammalib.GCTAEventList() list.gti(gti) - + # Set optional information if has_ebounds: list.ebounds(ebounds) if has_roi: list.roi(roi) - + # Attach event list to CTA observation obs.events(list) @@ -358,10 +358,10 @@ class csobsdef(ctools.cscript): # Create CTA response response = gammalib.GCTAResponseIrf() response.read(xml) - + # Attach response to observation obs.response(response) - + # Return observation return obs @@ -375,10 +375,9 @@ if __name__ == '__main__': """ # Create instance of application app = csobsdef(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - \ No newline at end of file diff --git a/cscripts/cspull.py b/cscripts/cspull.py index 9c4a686..3430e3d 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -42,7 +42,7 @@ class cspull(ctools.cscript): # Set name self.name = "cspull" self.version = "1.0.0" - + # Initialise some members self.obs = None self.model = None @@ -53,7 +53,7 @@ class cspull(ctools.cscript): self.m_psfcube = None self.m_bckcube = None self.m_stackmodels = None - + # Make sure that parfile exists file = self.parfile() @@ -71,7 +71,7 @@ class cspull(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -79,7 +79,7 @@ class cspull(ctools.cscript): # Write separator into logger if self.logTerse(): self.log("\n") - + # Return return @@ -90,13 +90,13 @@ class cspull(ctools.cscript): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found print("Parfile "+parfile+" not found. Create default parfile.") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file")) @@ -126,29 +126,29 @@ class cspull(ctools.cscript): pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","cspull.log","","","Log filename")) pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Get parameters - + # Set observation if not done before if self.obs == None or self.obs.size() == 0: self.obs = self.get_observations() - + # Check for requested pattern and use above # observation parameters to set wobble pattern self.m_pattern = self["pattern"].string() if self.m_pattern == "four": self.obs = self.set_obs() - + # Get number of energy bins self.m_enumbins = self["enumbins"].integer() - + # Read parameters for binned if requested if not self.m_enumbins == 0: self.m_npix = self["npix"].integer() @@ -157,11 +157,11 @@ class cspull(ctools.cscript): # Set dummy values (required by obsutils) self.m_npix = 0 self.m_binsz = 0.0 - + # Set models if we have none if self.obs.models().size() == 0: self.obs.models(self["inmodel"].filename()) - + # Read other parameters self.m_outfile = self["outfile"].filename() self.m_ntrials = self["ntrials"].integer() @@ -176,7 +176,7 @@ class cspull(ctools.cscript): # Return return - + def models(self, models): """ @@ -184,10 +184,10 @@ class cspull(ctools.cscript): """ # Copy models self.obs.models(models.copy()) - + # Return return - + def execute(self): """ @@ -195,11 +195,11 @@ class cspull(ctools.cscript): """ # Run the script self.run() - + # Return return - + def set_stacked_irf(self): """ For stacked analysis prepare stacked irfs @@ -227,7 +227,7 @@ class cspull(ctools.cscript): # Notify exposure computation if self.logTerse(): self.log("Computed exposure cube\n") - + # Get stacked Psf psfcube = ctools.ctpsfcube(self.obs) psfcube["incube"] = "NONE" @@ -246,7 +246,7 @@ class cspull(ctools.cscript): # Notify Psf computation if self.logTerse(): self.log("Computed Psf cube\n") - + # Get stacked background bkgcube = ctools.ctbkgcube(self.obs) bkgcube["incube"] = "NONE" @@ -265,17 +265,17 @@ class cspull(ctools.cscript): # Notify background cube computation if self.logTerse(): self.log("Computed background cube\n") - + # Store results self.m_exposure = expcube.expcube().copy() self.m_psfcube = psfcube.psfcube().copy() self.m_bckcube = bkgcube.bkgcube().copy() self.m_stackmodels = bkgcube.models().copy() - + # Return return - - + + def run(self): """ Run the script. @@ -286,12 +286,12 @@ class cspull(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() self.log("\n") - + # Write observation into logger if self.logTerse(): self.log("\n") @@ -310,10 +310,10 @@ class cspull(ctools.cscript): # Loop over trials for seed in range(self.m_ntrials): - + # Make a trial result = self.trial(seed) - + # Write out result immediately if seed == 0: file = open(self.m_outfile, 'w') @@ -324,23 +324,23 @@ class cspull(ctools.cscript): writer = csv.DictWriter(file, result['colnames']) writer.writerow(result['values']) file.close() - + # Return return - + def set_obs(self): """ Returns an observation container with a set of CTA observations. - + Keywords: """ - + # Setup observation definition list obsdeflist = obsutils.set_obs_patterns(self.m_pattern, ra=self["ra"].real(), dec=self["dec"].real(), offset=self["offset"].real()) - + # Create list of observations obs = obsutils.set_obs_list(obsdeflist, tstart=self["tmin"].real(), duration=self["tmax"].real()-self["tmin"].real(), @@ -348,7 +348,7 @@ class cspull(ctools.cscript): emin=self["emin"].real(), emax=self["emax"].real(), rad=self["rad"].real(), irf=self["irf"].string(), caldb=self["caldb"].string()) - + # Return observation container return obs @@ -356,11 +356,11 @@ class cspull(ctools.cscript): def trial(self, seed): """ Create the pull for a single trial. - + Parameters: seed - Random number generator seed """ - + # Write header if self.logNormal(): self.log.header2("Trial "+str(seed+1)) @@ -375,12 +375,12 @@ class cspull(ctools.cscript): log=self.m_log, debug=self.m_debug, chatter=self.m_chatter) - + # If stacked, add stacked responses and model if self.obs.size() > 0 and self.m_enumbins > 0: obs[0].response(self.m_exposure, self.m_psfcube, self.m_bckcube) obs.models(self.m_stackmodels) - + # Determine number of events in simulation nevents = 0.0 for run in obs: @@ -416,7 +416,7 @@ class cspull(ctools.cscript): # Write result header if self.logNormal(): self.log.header3("Pulls") - + # Gather results colnames = [] values = {} @@ -432,15 +432,15 @@ class cspull(ctools.cscript): for k in range(model.size()): par = model[k] if par.is_free(): - + # Set parameter name name = model_name+"_"+par.name() - + # Append parameter, Pull_parameter and Unc_parameter colnames.append(name) colnames.append("Pull_"+name) colnames.append("Unc_"+name) - + # Compute pull fitted_value = par.value() real_value = self.obs.models()[i][k].value() @@ -449,7 +449,7 @@ class cspull(ctools.cscript): pull = (fitted_value - real_value) / error else: pull = 99.0 - + # Store results values[name] = fitted_value values["Pull_"+name] = pull @@ -464,10 +464,10 @@ class cspull(ctools.cscript): self.log(" +/- ") self.log(error) self.log(")\n") - + # Bundle together results result = {'colnames': colnames, 'values': values} - + # Return return result @@ -481,10 +481,9 @@ if __name__ == '__main__': """ # Create instance of application app = cspull(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index c64adbe..9561544 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -36,11 +36,11 @@ class csresmap(ctools.cscript): """ Constructor. """ - + # Set name self.name = "csresmap" self.version = "1.0.0" - + # Initialise some members self.obs = None self.algorithm = "SUB" @@ -48,7 +48,7 @@ class csresmap(ctools.cscript): self.m_modcube = "None" self.m_use_maps = False self.m_skip_binning = False - + # Initialise some members if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): self.obs = argv[0] @@ -57,7 +57,7 @@ class csresmap(ctools.cscript): self.obs = gammalib.GObservations() self.obs.clear() self.outfile = "" - + # Make sure that parfile exists file = self.parfile() @@ -75,7 +75,7 @@ class csresmap(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -83,7 +83,7 @@ class csresmap(ctools.cscript): # Write separator into logger if self.logTerse(): self.log("\n") - + # Return return @@ -94,13 +94,13 @@ class csresmap(ctools.cscript): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file")) @@ -128,10 +128,10 @@ class csresmap(ctools.cscript): pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","csresmap.log","","","Log filename")) pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -153,28 +153,28 @@ class csresmap(ctools.cscript): self.m_modcube = self["modcube"].filename() if self.m_modcube != "NONE": self.m_use_maps = True - + # If not two maps are given, proceed to set up observation if not self.m_use_maps: - + # Set observation if not done before if self.obs.size() == 0: self.require_inobs("csresmap.get_parameters()") self.obs = self.get_observations() - + # Check if we have exactly one binned CTA observation if self.obs.size() == 1: - + if self.obs[0].classname() == "GCTAObservation": if self.obs[0].eventtype() == "CountsCube": # Skip ctbin step later on self.m_skip_binning = True - - + + # Set models if we have none if self.obs.models().size() == 0: self.obs.models(self["inmodel"].filename()) - + # Skip query for spatial parameters if a binning is provided in the observation if not self.m_skip_binning: # Read other parameters @@ -189,30 +189,30 @@ class csresmap(ctools.cscript): self.m_nxpix = self["nxpix"].integer() self.m_nypix = self["nypix"].integer() self.m_binsz = self["binsz"].real() - + # Read necessary parameters self.m_outfile = self["outmap"].filename() self.m_algorithm = self["algorithm"].string() - + # Set some fixed parameters self.m_log = False # Logging in client tools self.m_chatter = self["chatter"].integer() self.m_clobber = self["clobber"].boolean() self.m_debug = self["debug"].boolean() - + # Return return - + def models(self, models): """ Set model. """ # Copy models self.obs.models(models.clone()) - + # Return return - + def execute(self): """ Execute the script. @@ -222,7 +222,7 @@ class csresmap(ctools.cscript): # Save residual map self.resmap.save(self.m_outfile, self.m_clobber) - + # Return return @@ -236,12 +236,12 @@ class csresmap(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() self.log("\n") - + # Write observation into logger if self.logTerse(): self.log("\n") @@ -253,19 +253,19 @@ class csresmap(ctools.cscript): if self.m_use_maps: countmap = gammalib.GSkymap(self["inobs"].filename()) modelmap = gammalib.GSkymap(self.m_modcube) - + else: - + if self.m_skip_binning: cta_counts_cube = gammalib.GCTAEventCube(self.obs[0].events().clone()) else: - + # Write header if self.logTerse(): self.log("\n") self.log.header1("Generate binned map (ctbin)") - + # Create countsmap bin = ctools.ctbin(self.obs) bin["nxpix"].integer(self.m_nxpix) @@ -283,18 +283,18 @@ class csresmap(ctools.cscript): bin["clobber"].boolean(self.m_clobber) bin["debug"].boolean(self.m_debug) bin.run() - + # Retrieve counts cube cta_counts_cube = bin.cube() - + # Assign GCTAEventCube to skymap countmap = cta_counts_cube.map() - + # Write header if self.logTerse(): self.log("\n") self.log.header1("Generate model map (ctmodel)") - + # Create model map model = ctools.ctmodel(self.obs) model.cube(cta_counts_cube) @@ -302,11 +302,11 @@ class csresmap(ctools.cscript): model["clobber"].boolean(self.m_clobber) model["debug"].boolean(self.m_debug) model.run() - + # Get model map into GSkymap object modelmap = model.cube().map().copy() - - + + # Store counts map as residual map. Note that we need a # special construct here to avoid memory leaks. This seems # to be a SWIG feature as SWIG creates a new object when @@ -315,15 +315,15 @@ class csresmap(ctools.cscript): self.resmap = countmap.copy() self.resmap.stack_maps() modelmap.stack_maps() - + # Continue calculations depending on given algorithm if self.m_algorithm == "SUB": - + # Subtract maps self.resmap -= modelmap - + elif self.m_algorithm == "SUBDIV": - + # Subtract and divide by model map self.resmap -= modelmap self.resmap /= modelmap @@ -331,7 +331,7 @@ class csresmap(ctools.cscript): # if pixel != 0.0: # pixel = 1.0/pixel #self.resmap *= modelmap - + elif self.m_algorithm == "SUBDIVSQRT": # subtract and divide by sqrt of model map @@ -341,15 +341,15 @@ class csresmap(ctools.cscript): # if pixel != 0.0: # pixel = 1.0/math.sqrt(pixel) #self.resmap *= modelmap - + else: - + # Raise error if algorithm is unkown raise TypeError("Algorithm \""+self.m_algorithm+"\" not known") - + # Return return - + # ======================== # # Main routine entry point # @@ -360,10 +360,9 @@ if __name__ == '__main__': """ # Create instance of application app = csresmap(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - \ No newline at end of file diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 11edebc..7d8d357 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -50,7 +50,7 @@ class cssens(ctools.cscript): self.obs = None self.m_ra = None self.m_dec = None - + # Make sure that parfile exists file = self.parfile() @@ -69,7 +69,7 @@ class cssens(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -77,7 +77,7 @@ class cssens(ctools.cscript): # Write separator into logger if self.logTerse(): self.log("\n") - + # Return return @@ -88,13 +88,13 @@ class cssens(ctools.cscript): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found print("Parfile "+parfile+" not found. Create default parfile.") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file")) @@ -120,10 +120,10 @@ class cssens(ctools.cscript): pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","cssens.log","","","Log filename")) pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get user parameters from parfile. @@ -160,7 +160,7 @@ class cssens(ctools.cscript): self.m_max_iter = self["max_iter"].integer() self.m_num_avg = self["num_avg"].integer() self.m_type = self["type"].string() - + # Set some fixed parameters self.m_log = False # Logging in client tools self.m_debug = self["debug"].boolean() # Debugging in client tools @@ -169,17 +169,17 @@ class cssens(ctools.cscript): self.m_ebounds = gammalib.GEbounds(self.m_bins, gammalib.GEnergy(self.m_emin, "TeV"), gammalib.GEnergy(self.m_emax, "TeV")) - + # Return return - + def execute(self): """ Execute the script. """ # Run the script self.run() - + # Return return @@ -193,12 +193,12 @@ class cssens(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() self.log("\n") - + # Initialise script colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', 'energy_flux', 'sensitivity'] @@ -228,7 +228,7 @@ class cssens(ctools.cscript): # Loop over energy bins for ieng in range(self.m_ebounds.size()): - + # Set energies if self.m_type == "Differential": emin = self.m_ebounds.emin(ieng) @@ -240,11 +240,11 @@ class cssens(ctools.cscript): msg = "Invalid sensitivity type \""+self.m_type+"\" encountered."+ \ " Either use \"Differential\" or \"Integral\"." raise gammalib.GException.invalid_value("cssens", msg) - + # Determine sensitivity result = self.get_sensitivity(self.obs, emin, emax, bkg_model, full_model) - + # Write results if ieng == 0: f = open(self.m_outfile, 'w') @@ -257,17 +257,17 @@ class cssens(ctools.cscript): writer = csv.DictWriter(f, colnames) writer.writerow(result) f.close() - + # Append results results.append(result) - + # Return return def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): """ Returns an observation container. - + Keywords: lpnt - Galactic longitude of pointing [deg] (default: 0.0) bpnt - Galactic latitude of pointing [deg] (default: 0.0) @@ -291,16 +291,16 @@ class cssens(ctools.cscript): # Allocate observation container obs = gammalib.GObservations() - + # Set single pointing pntdir = gammalib.GSkyDir() pntdir.lb_deg(lpnt, bpnt) - + # Create CTA observation run = obsutils.set_obs(pntdir, caldb=caldb, irf=irf, duration=duration, deadc=deadc, emin=emin, emax=emax, rad=rad) - + # Append observation to container obs.append(run) @@ -309,25 +309,25 @@ class cssens(ctools.cscript): pntdir.lb_deg(lpnt, bpnt+offset) self.m_ra = pntdir.ra_deg() self.m_dec = pntdir.dec_deg() - + # Return observation container return obs def set_obs_ebounds(self, emin, emax): """ Set energy boundaries for observation in container. - + Parameters: emin - Minimum energy emax - Maximum energy """ # Loop over all observations in container for obs in self.obs: - + # Set energy boundaries ebounds = gammalib.GEbounds(emin, emax) obs.events().ebounds(ebounds) - + # Return return @@ -340,7 +340,7 @@ class cssens(ctools.cscript): # Get source model model = full_model[self.m_srcname] - + # Check that model has a Prefactor if not model.has_par("Prefactor"): msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \ @@ -399,7 +399,7 @@ class cssens(ctools.cscript): def get_sensitivity(self, obs, emin, emax, bkg_model, full_model): """ Determine sensitivity for given observations. - + Parameters: obs - Observation container emin - Minimum energy for fitting and flux computation @@ -412,7 +412,7 @@ class cssens(ctools.cscript): # Set energy boundaries self.set_obs_ebounds(emin, emax) - + # Determine energy boundaries from first observation in the container loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) e_mean = math.pow(10.0, loge) @@ -517,7 +517,7 @@ class cssens(ctools.cscript): self.log("\n") for par in model: self.log(str(par)+"\n") - + # Fit background and test source sim.models(src_model) like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) @@ -580,7 +580,7 @@ class cssens(ctools.cscript): correct = 1.0 if ts > 0: correct = math.sqrt(self.m_ts_thres/ts) - + # Compute extrapolated fluxes crab_flux = correct * crab_flux photon_flux = correct * photon_flux @@ -590,7 +590,7 @@ class cssens(ctools.cscript): photon_flux_value.append(photon_flux) energy_flux_value.append(energy_flux) sensitivity_value.append(sensitivity) - + # Write background and test source fit results if self.logExplicit(): self.log.parformat("Photon flux") @@ -627,7 +627,7 @@ class cssens(ctools.cscript): self.log(" mCrab = ") self.log(sensitivity) self.log(" erg/cm2/s\n") - + # Compute sliding average of extrapolated fitted prefactor, # photon and energy flux. This damps out fluctuations and # improves convergence @@ -648,12 +648,12 @@ class cssens(ctools.cscript): photon_flux /= num energy_flux /= num sensitivity /= num - + # Compare average flux to last average if iter > self.m_num_avg: if test_crab_flux > 0: ratio = crab_flux/test_crab_flux - + # We have 2 convergence criteria: # 1. The average flux does not change # 2. The flux correction factor is small @@ -666,10 +666,10 @@ class cssens(ctools.cscript): if self.logTerse(): self.log(" Flux is zero.\n") break - + # Use average for next iteration test_crab_flux = crab_flux - + # Exit loop if number of trials exhausted if (iter >= self.m_max_iter): if self.logTerse(): @@ -723,13 +723,13 @@ class cssens(ctools.cscript): self.log("\n") for par in model: self.log(str(par)+"\n") - + # Store result result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ 'crab_flux': crab_flux, 'photon_flux': photon_flux, \ 'energy_flux': energy_flux, \ 'sensitivity': sensitivity} - + # Return result return result @@ -756,7 +756,6 @@ if __name__ == '__main__': """ # Create instance of application app = cssens(sys.argv) - + # Run application app.execute() - \ No newline at end of file diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 65ad310..41dd06a 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -37,14 +37,14 @@ class csspec(ctools.cscript): """ Constructor. """ - + # Set name self.name = "csspec" self.version = "1.0.0" - + # Initialise some members self.obs = None - + # Initialise some members if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): self.obs = argv[0] @@ -53,7 +53,7 @@ class csspec(ctools.cscript): self.obs = gammalib.GObservations() self.obs.clear() self.m_outfile = "" - + # Make sure that parfile exists file = self.parfile() @@ -71,7 +71,7 @@ class csspec(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -79,7 +79,7 @@ class csspec(ctools.cscript): # Write separator into logger if self.logTerse(): self.log("\n") - + # Return return @@ -90,13 +90,13 @@ class csspec(ctools.cscript): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file")) @@ -128,10 +128,10 @@ class csspec(ctools.cscript): pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","csspec.log","","","Log filename")) pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -144,10 +144,10 @@ class csspec(ctools.cscript): # Set models if we have none if self.obs.models().size() == 0: self.obs.models(self["inmodel"].filename()) - + # Get source name self.m_srcname = self["srcname"].string() - + # Get ebounds self.m_ebounds = self.create_ebounds() @@ -165,32 +165,32 @@ class csspec(ctools.cscript): # Read other parameters self.m_outfile = self["outfile"].filename() - + # Get other parameeters self.m_calc_ulimit = self["calc_ulim"].boolean() self.m_calc_ts = self["calc_ts"].boolean() self.m_fix_bkg = self["fix_bkg"].boolean() self.m_fix_srcs = self["fix_srcs"].boolean() - + # Set some fixed parameters self.m_log = False # Logging in client tools self.m_chatter = self["chatter"].integer() self.m_clobber = self["clobber"].boolean() self.m_debug = self["debug"].boolean() - + # Return return - + def models(self, models): """ Set model. """ # Copy models self.obs.models(models.clone()) - + # Return return - + def execute(self): """ Execute the script. @@ -200,7 +200,7 @@ class csspec(ctools.cscript): # Save residual map self.fits.saveto(self.m_outfile, self.m_clobber) - + # Return return @@ -214,12 +214,12 @@ class csspec(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() self.log("\n") - + # Write observation into logger if self.logTerse(): self.log("\n") @@ -234,15 +234,15 @@ class csspec(ctools.cscript): # Adjust model parameters dependent on input user parameters for model in self.obs.models(): - + # Set TS flag for all models to false. # Source of interest will be set to true later model.tscalc(False) - + # Log model name if self.logExplicit(): self.log.header3(model.name()) - + # Deal with the source of interest if model.name() == self.m_srcname: for par in model: @@ -254,32 +254,32 @@ class csspec(ctools.cscript): model.spectral()[0].free() if self.m_calc_ts: model.tscalc(True) - + elif self.m_fix_bkg and not model.classname() == "GModelSky": for par in model: if par.is_free() and self.logExplicit(): self.log(" Fixing \""+par.name()+"\"\n") par.fix() - + elif self.m_fix_srcs and model.classname() == "GModelSky": for par in model: if par.is_free() and self.logExplicit(): self.log(" Fixing \""+par.name()+"\"\n") par.fix() - + # Write header if self.logTerse(): self.log("\n") self.log.header1("Generate spectrum") - + # Initialise FITS Table with extension "SPECTRUM" table = gammalib.GFitsBinTable(self.m_ebounds.size()) table.extname("SPECTRUM") - + # Add Header for compatibility with gammalib.GMWLSpectrum table.card("INSTRUME", "CTA", "Name of Instrument") table.card("TELESCOP", "CTA", "Name of Telescope") - + # Create FITS table columns energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size()) energy.unit("TeV") @@ -295,7 +295,7 @@ class csspec(ctools.cscript): ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size()) ulim_values.unit("erg/cm2/s") Npred_values = gammalib.GFitsTableDoubleCol("Npred", self.m_ebounds.size()) - + # Loop over energy bins for i in range(self.m_ebounds.size()): @@ -310,18 +310,18 @@ class csspec(ctools.cscript): emax = self.m_ebounds.emax(i) elogmean = self.m_ebounds.elogmean(i) elogmean2 = elogmean.MeV() * elogmean.MeV() - + # Store energy as TeV energy[i] = elogmean.TeV() - + # Store energy errors energy_low[i] = (elogmean - emin).TeV() energy_high[i] = (emax - elogmean).TeV() - + # Log information if self.logExplicit(): - self.log.header3("Selecting events") - + self.log.header3("Selecting events") + # Select events select = ctools.ctselect(self.obs) select["emin"] = emin.TeV() @@ -342,7 +342,7 @@ class csspec(ctools.cscript): # Header if self.logTerse(): self.log.header3("Binning events") - + # Bin events bin = ctools.ctbin(select.obs()) bin["usepnt"] = False @@ -358,11 +358,11 @@ class csspec(ctools.cscript): bin["coordsys"] = self.m_coordsys bin["proj"] = self.m_proj bin.run() - + # Header if self.logTerse(): self.log.header3("Creating exposure cube") - + # Create exposure cube expcube = ctools.ctexpcube(select.obs()) expcube["incube"] = "NONE" @@ -379,11 +379,11 @@ class csspec(ctools.cscript): expcube["coordsys"] = self.m_coordsys expcube["proj"] = self.m_proj expcube.run() - + # Header if self.logTerse(): self.log.header3("Creating PSF cube") - + # Create psf cube psfcube = ctools.ctpsfcube(select.obs()) psfcube["incube"] = "NONE" @@ -400,11 +400,11 @@ class csspec(ctools.cscript): psfcube["coordsys"] = self.m_coordsys psfcube["proj"] = self.m_proj psfcube.run() - + # Header if self.logTerse(): self.log.header3("Creating background cube") - + # Create background cube bkgcube = ctools.ctbkgcube(select.obs()) bkgcube["incube"] = "NONE" @@ -421,37 +421,37 @@ class csspec(ctools.cscript): bkgcube["coordsys"] = self.m_coordsys bkgcube["proj"] = self.m_proj bkgcube.run() - + # Set new binned observation obs = bin.obs() - + # Set precomputed binned response obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube()) # Get new models models = bkgcube.models() - + # Fix background models if required if self.m_fix_bkg: for model in models: if not model.classname() == "GModelSky": for par in model: par.fix() - + # Set new models to binned observation obs.models(models) - + # Header if self.logTerse(): self.log.header3("Performing fit") - + # Likelihood like = ctools.ctlike(obs) like.run() - + # Skip bin if no event was present if like.obs().logL() == 0.0: - + # Log information if self.logTerse(): self.log("No event in this bin. Bin is skipped\n") @@ -463,7 +463,7 @@ class csspec(ctools.cscript): ulim_values[i] = 0.0 Npred_values[i] = 0.0 continue - + # Get results fitted_models = like.obs().models() source = fitted_models[self.m_srcname] @@ -471,16 +471,16 @@ class csspec(ctools.cscript): # Calculate Upper Limit ulimit_value = -1.0 if self.m_calc_ulimit: - + # Logging information if self.logTerse(): self.log.header3("Computing upper limit") - + # Create upper limit object ulimit = ctools.ctulimit(like.obs()) ulimit["srcname"] = self.m_srcname ulimit["eref"] = elogmean.TeV() - + # Try to run upper limit and catch exceptions try: ulimit.run() @@ -489,37 +489,37 @@ class csspec(ctools.cscript): if self.logTerse(): self.log("Upper limit calculation failed\n") ulimit_value = -1.0 - + # Get TS value TS = -1.0 if self.m_calc_ts: TS = source.ts() - + # Compute Npred value Npred = 0.0 for observation in like.obs(): Npred += observation.npred(source) - + # Get differential flux fitted_flux = source.spectral().eval(elogmean,gammalib.GTime()) - + # Compute flux error parvalue = source.spectral()[0].value() rel_error = source.spectral()[0].error()/parvalue e_flux = fitted_flux*rel_error - + # Set values for storage TSvalues[i] = TS - + # Set npred values Npred_values[i] = Npred - + # Convert fluxes to nuFnu flux[i] = fitted_flux * elogmean2 * gammalib.MeV2erg flux_err[i] = e_flux * elogmean2 * gammalib.MeV2erg if ulimit_value > 0.0: ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg - + # Log information if self.logExplicit(): self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ") @@ -540,14 +540,14 @@ class csspec(ctools.cscript): table.append(TSvalues) table.append(ulim_values) table.append(Npred_values) - + # Create the FITS file now self.fits = gammalib.GFits() self.fits.append(table) - + # Return return - + # ======================== # # Main routine entry point # @@ -558,10 +558,9 @@ if __name__ == '__main__': """ # Create instance of application app = csspec(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - \ No newline at end of file diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py index af633fa..70b02c1 100755 --- a/cscripts/cstsdist.py +++ b/cscripts/cstsdist.py @@ -43,12 +43,12 @@ class cstsdist(ctools.cscript): # Set name self.name = "cstsdist" self.version = "1.0.0" - + # Initialise some members self.obs = None self.bkg_model = None self.full_model = None - + # Make sure that parfile exists file = self.parfile() @@ -66,7 +66,7 @@ class cstsdist(ctools.cscript): # Return return - + def __del__(self): """ Destructor. @@ -74,7 +74,7 @@ class cstsdist(ctools.cscript): # Write separator into logger if self.logTerse(): self.log("\n") - + # Return return @@ -85,13 +85,13 @@ class cstsdist(ctools.cscript): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file")) @@ -121,15 +121,15 @@ class cstsdist(ctools.cscript): pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","cstsdist.log","","","Log filename")) pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile and setup the observation. """ - + # Set observation if not done before if self.obs == None or self.obs.size() == 0: self.obs = self.get_observations() @@ -142,10 +142,10 @@ class cstsdist(ctools.cscript): # Get source name self.m_srcname = self["srcname"].string() - + # Get number of energy bins self.m_enumbins = self["enumbins"].integer() - + # Read parameters for binned if requested if not self.m_enumbins == 0: self.m_npix = self["npix"].integer() @@ -158,7 +158,7 @@ class cstsdist(ctools.cscript): # Set models if we have none if self.obs.models().size() == 0: self.obs.models(self["inmodel"].filename()) - + # Get other parameters self.m_outfile = self["outfile"].filename() self.m_ntrials = self["ntrials"].integer() @@ -171,24 +171,24 @@ class cstsdist(ctools.cscript): # Return return - + def models(self, models): """ Set model. """ # Copy models self.obs.models(models.copy()) - + # Return return - + def execute(self): """ Execute the script. """ # Run the script self.run() - + # Return return @@ -202,7 +202,7 @@ class cstsdist(ctools.cscript): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() @@ -221,7 +221,7 @@ class cstsdist(ctools.cscript): self.log.header2("Full model") self.log(str(full_model)) self.log("\n") - + # Write observation into logger if self.logTerse(): self.log("\n") @@ -243,10 +243,10 @@ class cstsdist(ctools.cscript): # Loop over trials for seed in range(self.m_ntrials): - + # Make a trial result = self.trial(seed, full_model, bkg_model) - + # Write out result immediately if seed == 0: file = open(self.m_outfile, 'w') @@ -257,10 +257,10 @@ class cstsdist(ctools.cscript): writer = csv.DictWriter(file, result['colnames']) writer.writerow(result['values']) file.close() - + # Return return - + def set_models(self, fitpos=False, fitspec=False): """ Set full and background model. @@ -270,7 +270,7 @@ class cstsdist(ctools.cscript): # Get source model model = full_model[self.m_srcname] - + # Check that model has a Prefactor if not model.has_par("Prefactor"): msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \ @@ -323,16 +323,16 @@ class cstsdist(ctools.cscript): def set_obs(self): """ Returns an observation container with a set of CTA observations. - + Keywords: """ - + # Setup observation definition list obsdeflist = obsutils.set_obs_patterns(self.m_pattern, ra=self["ra"].real(), dec=self["dec"].real(), offset=self["offset"].real()) - + # Create list of observations obs = obsutils.set_obs_list(obsdeflist, tstart=self["tmin"].real(), @@ -342,14 +342,14 @@ class cstsdist(ctools.cscript): emax=self["emax"].real(), rad=self["rad"].real(), irf=self["irf"].string(), caldb=self["caldb"].string()) - + # Return observation container return obs def trial(self, seed, full_model, bkg_model): """ Create the TS for a single trial. - + Parameters: seed - Random number generator seed """ @@ -364,7 +364,7 @@ class cstsdist(ctools.cscript): binsz=self.m_binsz, npix=self.m_npix, log=self.m_log, debug=self.m_debug) - + # Determine number of events in simulation nevents = 0.0 for run in sim: @@ -376,7 +376,7 @@ class cstsdist(ctools.cscript): self.log.parformat("Number of simulated events") self.log(nevents) self.log("\n") - + # Fit background only sim.models(bkg_model) like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) @@ -407,7 +407,7 @@ class cstsdist(ctools.cscript): LogL_all = like_all.opt().value() npred_all = like_all.obs().npred() ts = 2.0*(LogL_bgm-LogL_all) - + # Write background and test source fit results if self.logExplicit(): self.log.header3("Background and test source model fit") @@ -426,7 +426,7 @@ class cstsdist(ctools.cscript): self.log("\n") for par in model: self.log(str(par)+"\n") - + # Write result elif self.logTerse(): self.log.parformat("Trial "+str(seed)) @@ -437,11 +437,11 @@ class cstsdist(ctools.cscript): self.log("+/-") self.log(result_all[self.m_srcname]["Prefactor"].error()) self.log("\n") - + # Initialise results colnames = [] values = {} - + # Set TS value colnames.append("TS") values["TS"] = ts @@ -465,7 +465,7 @@ class cstsdist(ctools.cscript): # Set Npred for full fit colnames.append("Npred_all") values["Npred_all"] = npred_all - + # Gather free full fit parameters for i in range(result_all.size()): model = result_all[i] @@ -473,22 +473,22 @@ class cstsdist(ctools.cscript): for k in range(model.size()): par = model[k] if par.is_free(): - + # Set parameter name name = model_name+"_"+par.name() - + # Append value colnames.append(name) values[name] = par.value() - + # Append error name = "Unc_"+name colnames.append(name) values[name] = par.error() - + # Bundle together results result = {'colnames': colnames, 'values': values} - + # Return return result @@ -502,10 +502,9 @@ if __name__ == '__main__': """ # Create instance of application app = cstsdist(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - \ No newline at end of file diff --git a/cscripts/cta_kb_bkgrates.py b/cscripts/cta_kb_bkgrates.py index 4c473e0..adebc0a 100755 --- a/cscripts/cta_kb_bkgrates.py +++ b/cscripts/cta_kb_bkgrates.py @@ -32,7 +32,7 @@ import sys def make_file_function(irfname, filename): """ Generate file function from model. - + Parameters: irfname - IRF filename filename - File function filename @@ -40,16 +40,16 @@ def make_file_function(irfname, filename): # Open performance file and file function irf = open(irfname, "r") file = open(filename, "w") - + # Read data for row in irf: - + # Skip header and trailer if row.find("log(E)") != -1: continue if row.find("--------") != -1: break - + # Extract data fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",") logE = float(fields[0]) @@ -59,16 +59,16 @@ def make_file_function(irfname, filename): eres = float(fields[4]) bgd = float(fields[5]) sens = float(fields[6]) - + # Compute energy (in MeV) energy = math.pow(10.0, logE)*1.0e6 emin = math.pow(10.0, logE-0.1)*1.0e6 emax = math.pow(10.0, logE+0.1)*1.0e6 ewidth = emax-emin - + # Compute solid angle of r80 region in sr omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) - + # Compute background rate per steradian and MeV if omega > 0.0: bkg_rate = bgd/omega/ewidth @@ -77,11 +77,11 @@ def make_file_function(irfname, filename): # Write file.write(str(energy)+" "+str(bkg_rate)+"\n") - + # Close files irf.close() file.close() - + # Return return @@ -111,6 +111,6 @@ if __name__ == '__main__': irf = tail.strip(".dat") filename = "bkg_"+irf+".txt" print(filename) - + # Make background file make_file_function(irfname, filename) diff --git a/cscripts/cta_kb_root2caldb.py b/cscripts/cta_kb_root2caldb.py index 37c8244..2ae96ab 100755 --- a/cscripts/cta_kb_root2caldb.py +++ b/cscripts/cta_kb_root2caldb.py @@ -50,7 +50,7 @@ class caldb(): """ # Store root path self.path = path - + # Initialise some members self.cif = None self.irf = None @@ -62,7 +62,7 @@ class caldb(): # Store UTC date of job execution self.utc = datetime.utcnow().isoformat()[:19] - + # Set CALDB information self.cal_tel = "CTA" self.cal_inst = inst.upper() @@ -99,13 +99,13 @@ class caldb(): self.bgd_doc = "???" self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] self.bgd_desc = "CTA background" - + # Create directory structure self.make_dirs() - + # Return return - + def __del__(self): """ Destructor. Close all FITS files. @@ -115,18 +115,18 @@ class caldb(): # Return return - + def make_dirs(self): """ Generate CALDB directory structure for one observation identifier. The structure is given by - + data///bcf/ - + where is "cta" and is the instrument specified in the CALDB constructor (the instrument may be used for different array configurations). - + Parameters: None Keywords: @@ -139,7 +139,7 @@ class caldb(): # Set directory names self.bcf_dir = self.base_dir+"/bcf" - + # Optionally prefix path if len(self.path) > 0: self.base_path = self.path+"/"+self.base_dir @@ -151,17 +151,17 @@ class caldb(): # Create directory structure if not os.path.isdir(self.bcf_path): os.makedirs(self.bcf_path) - + # Return return - + def add(self, rspname, split=False, clobber=True): """ Add new calibration. The actual version will put all calibrations in the same file, although each part of the response function will have its own logical name. We can thus easily modify the script to put each calibration information in a separate file. - + Parameters: rspname - Response name Keywords: @@ -179,21 +179,21 @@ class caldb(): self.psf_file = "irf_"+rspname+".dat" self.edisp_file = "irf_"+rspname+".dat" self.bgd_file = "irf_"+rspname+".dat" - + # Open calibration database index if self.cif == None: self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True) - + # If file has no CIF extension than create it now try: self.hdu_cif = self.cif.table("CIF") except: self.create_cif() self.hdu_cif = self.cif.table("CIF") - + # Set response name self.cal_name = "NAME("+rspname+")" - + # Return return @@ -210,7 +210,7 @@ class caldb(): # necessarily have all the information at hand (in particular # about the boundaries) #self.add_cif_info() - + # Close CIF if self.cif != None: self.cif.save(True) @@ -222,7 +222,7 @@ class caldb(): # Return return - + def create_cif(self): """ Create Calibration Index File (CIF) extension in CIF FITS file. @@ -234,9 +234,9 @@ class caldb(): """ # Create binary table table = gammalib.GFitsBinTable() - + # Set boundary - + # Attach columns. Reference: CAL/GEN/92-008 table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10)) table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10)) @@ -256,20 +256,20 @@ class caldb(): table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0)) table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8)) table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70)) - + # Set keywords. Reference: CAL/GEN/92-008 table.extname("CIF") table.card("CIFVERSN", "1992a", "Version of CIF format") - + # Attach table to FITS file. Note that at this moment the # FITS table gets connected to the FITS file. Yet since # nothing was yet written to the FITS file, we cannot read # anything from it. self.cif.append(table) - + # Return return - + def add_cif_info(self): """ Add information to CIF extension. @@ -287,13 +287,13 @@ class caldb(): # Append 4 rows to CIF extension self.hdu_cif.append_rows(4) - + # Add generic information for these 4 rows for i in range(4): - + # Set row number row = i+self.hdu_cif.nrows()-4 - + # Set element self.hdu_cif["TELESCOP"][row] = self.cal_tel self.hdu_cif["INSTRUME"][row] = self.cal_inst @@ -307,7 +307,7 @@ class caldb(): self.hdu_cif["REF_TIME"][row] = self.ref_time self.hdu_cif["CAL_QUAL"][row] = self.cal_qual self.hdu_cif["CAL_DATE"][row] = self.cal_date - + # Add effective area information row = self.hdu_cif.nrows()-4 self.hdu_cif["CAL_DIR"][row] = self.bcf_dir @@ -359,14 +359,14 @@ class caldb(): self.hdu_cif["CAL_CBD"][row, n] = "NONE" else: self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n] - + # Return return - + def root2caldb(self, filename): """ Generate files. - + Parameters: filename - Performance file name """ @@ -383,16 +383,16 @@ class caldb(): if self.ea_file == self.psf_file and \ self.ea_file == self.edisp_file and \ self.ea_file == self.bgd_file : - split = False - fprf = open(ea_filename, "w") - files = [fprf] + split = False + fprf = open(ea_filename, "w") + files = [fprf] else: - split = True - fea = open(ea_filename, "w") - fpsf = open(psf_filename, "w") - fedisp = open(edisp_filename, "w") - fbgd = open(bgd_filename, "w") - files = [fea, fpsf, fedisp, fbgd] + split = True + fea = open(ea_filename, "w") + fpsf = open(psf_filename, "w") + fedisp = open(edisp_filename, "w") + fbgd = open(bgd_filename, "w") + files = [fea, fpsf, fedisp, fbgd] # Write header for file in files: @@ -420,12 +420,12 @@ class caldb(): rootfile.GetObject("AngRes", angres) rootfile.GetObject("AngRes80", angres80) rootfile.GetObject("ERes", eres) - + # Setup vectors energies = bgrate.GetXaxis() nbins_eng = energies.GetNbins() for i in range(nbins_eng): - + # Get logE and background rate per squared degree logE = energies.GetBinCenter(i+1) bgd = bgrate.GetBinContent(i+1,1) @@ -436,7 +436,7 @@ class caldb(): r80 = angres80.GetBinContent(i+1,1) dE = eres.GetBinContent(i+1,1) diffs = sens.GetBinContent(i+1,1) - + # Skip any NaNs if math.isnan(area80): continue @@ -446,10 +446,10 @@ class caldb(): emin = math.pow(10.0, logE-0.1)*1.0e6 emax = math.pow(10.0, logE+0.1)*1.0e6 ewidth = emax-emin - + # Convert into background rate per steradian and MeV bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth - + # Compute control background rate (this only works for KB files # as they are for 80% containment radius). But we don't use this # anyways, it's just for control on display ... @@ -458,7 +458,7 @@ class caldb(): bkg_rate2 = bgd / omega / ewidth else: bkg_rate2 = 0.0 - + # Compute full effective area by dividing area80/0.80 area = area80/0.80 @@ -494,7 +494,7 @@ class caldb(): # necessarily have all the information at hand (in particular # about the boundaries) self.add_cif_info() - + # Return return @@ -521,7 +521,7 @@ if __name__ == '__main__': # Loop over all files for file in files: - + # Get filename head, tail = os.path.split(file) @@ -540,6 +540,6 @@ if __name__ == '__main__': # Open calibration file irf.add(name) - + # Translate ROOT to CALDB information irf.root2caldb(file) diff --git a/cscripts/cta_make_nodes.py b/cscripts/cta_make_nodes.py index 9555a15..24d01c9 100755 --- a/cscripts/cta_make_nodes.py +++ b/cscripts/cta_make_nodes.py @@ -56,7 +56,7 @@ def write_node(file, energy, scale): # Return return - + # ============ # # Create nodes # @@ -72,29 +72,29 @@ def create_nodes(emin, emax, enumbins): elogmin = math.log10(float(emin)) elogmax = math.log10(float(emax)) elogbin = (elogmax - elogmin)/(float(enumbins)-1.0) - + # Fill arrays for i in range(int(enumbins)): - + # Compute energy energy = math.pow(10.0, i*elogbin+elogmin) - + # Compute scale (HESS Crab spectrum) scale = 3.45e-17 * math.pow(energy, -2.63) - + # Write node write_node(file, energy, scale) - + # Debug #sys.stdout.write(energy) #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)±"\n") - + # Close file file.close() - + # Return return - + # ======================= # # Main script entry point # @@ -108,7 +108,7 @@ if __name__ == '__main__': if len(sys.argv) < 3: sys.stdout.write(usage+"\n") sys.exit() - + # Extract parameters emin = sys.argv[1] emax = sys.argv[2] @@ -116,4 +116,3 @@ if __name__ == '__main__': # Create nodes create_nodes(emin, emax, enumbins) - diff --git a/cscripts/cta_root2caldb.py b/cscripts/cta_root2caldb.py index 7c4e062..93c7725 100755 --- a/cscripts/cta_root2caldb.py +++ b/cscripts/cta_root2caldb.py @@ -52,7 +52,7 @@ class caldb(): """ # Store root path self.path = path - + # Initialise some members self.cif = None self.irf = None @@ -64,7 +64,7 @@ class caldb(): # Store UTC date of job execution self.utc = datetime.utcnow().isoformat()[:19] - + # Set CALDB information self.cal_tel = "CTA" self.cal_inst = inst.upper() @@ -102,13 +102,13 @@ class caldb(): self.bgd_doc = "???" self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] self.bgd_desc = "CTA background" - + # Create directory structure self.make_dirs() - + # Return return - + def __del__(self): """ Destructor. Close all FITS files. @@ -118,18 +118,18 @@ class caldb(): # Return return - + def make_dirs(self): """ Generate CALDB directory structure for one observation identifier. The structure is given by - + data///bcf/ - + where is "cta" and is the instrument specified in the CALDB constructor (the instrument may be used for different array configurations). - + Parameters: None Keywords: @@ -147,7 +147,7 @@ class caldb(): self.psf_dir = self.obs_dir self.edisp_dir = self.obs_dir self.bgd_dir = self.obs_dir - + # Optionally prefix path if len(self.path) > 0: self.base_path = self.path+"/"+self.base_dir @@ -175,17 +175,17 @@ class caldb(): os.makedirs(self.edisp_path) if not os.path.isdir(self.bgd_path): os.makedirs(self.bgd_path) - + # Return return - + def open(self, version, split=False, clobber=True): """ Open existing or create new calibration. The actual version will put all calibrations in the same file, although each part of the response function will have its own logical name. We can thus easily modify the script to put each calibration information in a separate file. - + Parameters: version - Filename version Keywords: @@ -203,17 +203,17 @@ class caldb(): self.psf_file = "irf_"+version+".fits" self.edisp_file = "irf_"+version+".fits" self.bgd_file = "irf_"+version+".fits" - + # Open calibration database index self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True) - + # If file has no CIF extension than create it now try: self.hdu_cif = self.cif.table("CIF") except: self.create_cif() self.hdu_cif = self.cif.table("CIF") - + # Set filenames ea_filename = self.ea_path+"/"+self.ea_file psf_filename = self.psf_path+"/"+self.psf_file @@ -232,13 +232,13 @@ class caldb(): self.psf = self.irf self.edisp = self.irf self.bgd = self.irf - + # Open HDUs self.open_ea() self.open_psf() self.open_edisp() self.open_bgd() - + # Return return @@ -255,7 +255,7 @@ class caldb(): # necessarily have all the information at hand (in particular # about the boundaries) self.add_cif_info() - + # Close CIF if self.cif != None: self.cif.save(True) @@ -272,11 +272,11 @@ class caldb(): self.psf = None self.edisp = None self.bgd = None - + # ... otherwise we have split files, so we have to close them # all else: - + # Close effective area file if self.ea != None: self.ea.save(True) @@ -300,10 +300,10 @@ class caldb(): self.bgd.save(True) self.bgd.close() self.bgd = None - + # Return return - + def create_cif(self): """ Create Calibration Index File (CIF) extension in CIF FITS file. @@ -315,9 +315,9 @@ class caldb(): """ # Create binary table table = gammalib.GFitsBinTable() - + # Set boundary - + # Attach columns. Reference: CAL/GEN/92-008 table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10)) table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10)) @@ -337,20 +337,20 @@ class caldb(): table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0)) table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8)) table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70)) - + # Set keywords. Reference: CAL/GEN/92-008 table.extname("CIF") table.card("CIFVERSN", "1992a", "Version of CIF format") - + # Attach table to FITS file. Note that at this moment the # FITS table gets connected to the FITS file. Yet since # nothing was yet written to the FITS file, we cannot read # anything from it. self.cif.append(table) - + # Return return - + def add_cif_info(self): """ Add information to CIF extension. @@ -362,13 +362,13 @@ class caldb(): """ # Append 4 rows to CIF extension self.hdu_cif.append_rows(4) - + # Add generic information for these 4 rows for i in range(4): - + # Set row number row = i+self.hdu_cif.nrows()-4 - + # Set element self.hdu_cif["TELESCOP"][row] = self.cal_tel self.hdu_cif["INSTRUME"][row] = self.cal_inst @@ -382,7 +382,7 @@ class caldb(): self.hdu_cif["REF_TIME"][row] = self.ref_time self.hdu_cif["CAL_QUAL"][row] = self.cal_qual self.hdu_cif["CAL_DATE"][row] = self.cal_date - + # Add effective area information row = self.hdu_cif.nrows()-4 self.hdu_cif["CAL_DIR"][row] = self.ea_dir @@ -434,7 +434,7 @@ class caldb(): self.hdu_cif["CAL_CBD"][row, n] = "NONE" else: self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n] - + # Return return @@ -455,20 +455,20 @@ class caldb(): except: # Create binary table table = gammalib.GFitsBinTable() - + # Set extension name table.extname("EFFECTIVE AREA") - + # Set OGIP keywords self.set_ogip_keywords(table, self.ea_doc, \ ["RESPONSE", self.ea_name]) - + # Append table self.ea.append(table) - + # Get extension self.hdu_ea = self.ea.table("EFFECTIVE AREA") - + # Return return @@ -489,20 +489,20 @@ class caldb(): except: # Create binary table table = gammalib.GFitsBinTable() - + # Set extension name table.extname("POINT SPREAD FUNCTION") - + # Set OGIP keywords self.set_ogip_keywords(table, self.psf_doc, \ ["RESPONSE", self.psf_name]) - + # Append table self.psf.append(table) # Get extension self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") - + # Return return @@ -523,10 +523,10 @@ class caldb(): except: # Create binary table table = gammalib.GFitsBinTable() - + # Set extension name table.extname("ENERGY DISPERSION") - + # Set OGIP keywords self.set_ogip_keywords(table, self.edisp_doc, \ ["RESPONSE", self.edisp_name]) @@ -536,7 +536,7 @@ class caldb(): # Get extension self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") - + # Return return @@ -557,10 +557,10 @@ class caldb(): except: # Create binary table table = gammalib.GFitsBinTable() - + # Set extension name table.extname("BACKGROUND") - + # Set OGIP keywords self.set_ogip_keywords(table, self.bgd_doc, \ ["RESPONSE", self.bgd_name]) @@ -570,7 +570,7 @@ class caldb(): # Get extension self.hdu_bgd = self.bgd.table("BACKGROUND") - + # Return return @@ -597,7 +597,7 @@ class caldb(): key = "HDUCLAS%d" % (i+1) hdu.card(key, item, "") hdu.card("HDUVERS", "1.0.0", "") - + # Return return @@ -608,7 +608,7 @@ class caldb(): offset angle columns, this method will simply add another data column. If name==None, the method will not append any data column. - + Parameters: array - ROOT 2D histogram. hdu - FITS HDU. @@ -690,7 +690,7 @@ class caldb(): offset angle columns, this method will simply add another data column. If name==None, the method will not append any data column. - + Parameters: array - ROOT 2D histogram. hdu - FITS HDU. @@ -701,7 +701,7 @@ class caldb(): """ # Set constants deg2sr = 0.01745329*0.01745329 - + # Extract energy and offset angle vectors energies = array.GetXaxis() neng = energies.GetNbins() @@ -825,7 +825,7 @@ class caldb(): histogram. If the HDU has already the columns, this method will simply add another data column. If name==None, the method will not append any data column. - + Parameters: array - ROOT 3D histogram. hdu - FITS HDU. @@ -925,11 +925,11 @@ class caldb(): # Return boundary information return bounds - + def root2caldb(self, filename, rebin=False, psftype="Gauss", eascale=1.0, bgdscale=1.0): """ Translate ROOT to CALDB information. - + Parameters: filename - ROOT 2D performance filename. Keywords: @@ -943,7 +943,7 @@ class caldb(): # Create effective area self.root2ea(file, rebin=rebin, eascale=eascale) - + # Create point spread function self.root2psf(file, psftype) @@ -953,7 +953,7 @@ class caldb(): # Create background #self.root2bgd(file) self.root2bgd3D(file, bgdscale=bgdscale) - + # Return return @@ -1076,7 +1076,7 @@ class caldb(): for ieng in range(neng): sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma r68.SetBinContent(ieng+1,ioff+1,sigma) - + # Compute scale histogram scale = r68.Clone() for ioff in range(noffset): @@ -1087,7 +1087,7 @@ class caldb(): else: value = 0.0 scale.SetBinContent(ieng+1,ioff+1,value) - + # Set zero histogram zero = r68.Clone() for ioff in range(noffset): @@ -1172,12 +1172,12 @@ class caldb(): # Continue only if both radii are positive if r_68 > 0 and r_80 > 0: - + # Derive constants for equation to solve a = 1.0 - 0.68 b = 1.0 - 0.80 c = r_68*r_68/(r_80*r_80) - + # Solve equation (a^x-1)/(b^x-1)=c for x using secant # method. Stop when we are better than 1e-6. x1 = -0.5 @@ -1202,7 +1202,7 @@ class caldb(): gamma = 1.0 - 1.0/x else: gamma = 1.0 - + # Compute sigma denominator = 2.0 * gamma * (math.pow(a, x) - 1.0) if denominator > 0.0: @@ -1228,7 +1228,7 @@ class caldb(): # Show results on console #print(ioff,ieng,gamma,sigma,x,f,iter) - + # Set bin contents gamma2D.SetBinContent(ieng+1,ioff+1,gamma) sigma2D.SetBinContent(ieng+1,ioff+1,sigma) @@ -1352,7 +1352,7 @@ class caldb(): # BGD_RECO self.make_2D(array, self.hdu_bgd, "BGD_RECO", "") - + # Return return @@ -1384,7 +1384,7 @@ class caldb(): # BGD (reconstructed energy) self.make_3D(array, self.hdu_bgd, "BGD", "1/s/MeV/sr", scale=bgdscale) - + # Return return @@ -1397,7 +1397,7 @@ class caldb(): hdu.card("CCLS0001", self.cal_class, "Calibration class") hdu.card("CDTP0001", self.cal_type, "Calibration type") hdu.card("CCNM0001", name, "Calibration name") - + # Set boundary keywords for n in range(9): keyname = "CBD%d0001" % (n+1) @@ -1406,12 +1406,12 @@ class caldb(): else: value = bounds[n] hdu.card(keyname, value, "Calibration boundary") - + # Set validity keywords hdu.card("CVSD0001", self.val_date, "Calibration validity start date (UTC)") hdu.card("CVST0001", self.val_time, "Calibration validity start time (UTC)") hdu.card("CDEC0001", desc, "Calibration description") - + # Return return @@ -1551,7 +1551,7 @@ if __name__ == '__main__': This script translates a CTA ROOT 2D performance files into a CALDB compliant CTA response. The script will create the response for 3 observations to illustrate the run-wise file structure. - + The CTA ROOT 2D performance files contain the following histograms: - DiffSens_offaxis (2D) - Differential sensitivity - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?) @@ -1592,9 +1592,9 @@ if __name__ == '__main__': # Allocate caldb entry irf = caldb(inst, id) - + # Open calibration file irf.open("file") - + # Translate ROOT to CALDB information irf.root2caldb(filename, rebin=rebin, psftype=psftype, eascale=eascale, bgdscale=bgdscale) diff --git a/cscripts/cta_root2irf.py b/cscripts/cta_root2irf.py index 466750e..e0f6654 100755 --- a/cscripts/cta_root2irf.py +++ b/cscripts/cta_root2irf.py @@ -40,7 +40,7 @@ import os def sigma_lima(Non, Noff, alpha=0.2): """ Compute Eq. (17) of Li & Ma (1983). - + Parameters: Non - Number of on counts Noff - Number of off counts @@ -55,7 +55,7 @@ def sigma_lima(Non, Noff, alpha=0.2): term1 = Non * math.log((alpha1/alpha)*arg1) term2 = Noff * math.log(alpha1*arg2) sigma = math.sqrt(2.0 * (term1 + term2)) - + # Return sensitivity return sigma @@ -66,7 +66,7 @@ def sigma_lima(Non, Noff, alpha=0.2): def Non_lima(sigma, Noff, alpha=0.2): """ Solve Eq. (17) of Li & Ma (1983) for Non - + Parameters: sigma - Required significance Noff - Number of off counts @@ -76,18 +76,18 @@ def Non_lima(sigma, Noff, alpha=0.2): # Set initial guess for Non S = 25.0 Non = S + alpha*Noff - + # Compute initial sigma estimate s = sigma_lima(Non, Noff, alpha=alpha) delta = s - sigma - + # Iterate while abs(delta) > 1.0e-6: S *= sigma/s Non = S + alpha*Noff s = sigma_lima(Non, Noff, alpha=alpha) delta = s - sigma - + # Return return Non @@ -101,7 +101,7 @@ def get_aeff(bgrate, emin, emax, \ Computes effective area assuming a Crab spectrum. This does not include any hadron cut efficiencies, and thus is not very useful in a first place ... - + Parameters: bgrate - Background rate (TeV) emin - Minimum energy (TeV) @@ -118,14 +118,14 @@ def get_aeff(bgrate, emin, emax, \ Noff = B/alpha Non = Non_lima(sigma, Noff) S = (Non - B) - + # Compute Crab flux in specified energy band inx = index+1.0 cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx)) - + # Compute effective area area = S / duration / cu - + # Return effective area return area @@ -136,7 +136,7 @@ def get_aeff(bgrate, emin, emax, \ def make_files(filename, prfname, bkgname): """ Generate files. - + Parameters: filename - Performance file name prfname - Performance file name @@ -146,7 +146,7 @@ def make_files(filename, prfname, bkgname): file = TFile(filename) fprf = open(prfname, "w") fbkg = open(bkgname, "w") - + # Write header fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n") @@ -173,12 +173,12 @@ def make_files(filename, prfname, bkgname): file.GetObject("AngRes", angres) file.GetObject("AngRes80", angres80) file.GetObject("ERes", eres) - + # Setup vectors energies = bgrate.GetXaxis() nbins_eng = energies.GetNbins() for i in range(nbins_eng): - + # Get logE and background rate per squared degree logE = energies.GetBinCenter(i+1) bgd = bgrate.GetBinContent(i+1,1) @@ -189,7 +189,7 @@ def make_files(filename, prfname, bkgname): r80 = angres80.GetBinContent(i+1,1) dE = eres.GetBinContent(i+1,1) diffs = sens.GetBinContent(i+1,1) - + # Skip any NaNs if math.isnan(area80): continue @@ -199,10 +199,10 @@ def make_files(filename, prfname, bkgname): emin = math.pow(10.0, logE-0.1)*1.0e6 emax = math.pow(10.0, logE+0.1)*1.0e6 ewidth = emax-emin - + # Convert into background rate per steradian and MeV bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth - + # Compute control background rate (this only works for KB files # as they are for 80% containment radius). But we don't use this # anyways, it's just for control on display ... @@ -211,7 +211,7 @@ def make_files(filename, prfname, bkgname): bkg_rate2 = bgd / omega / ewidth else: bkg_rate2 = 0.0 - + # Compute full effective area by dividing area80/0.80 area = area80/0.80 @@ -241,7 +241,7 @@ def make_files(filename, prfname, bkgname): # Close files fprf.close() fbkg.close() - + # Return return @@ -263,13 +263,13 @@ if __name__ == '__main__': # Get list of all performance files files = glob.glob(caldb+"/*.root") files.sort() - + # Loop over all files for file in files: - + # Get filename head, tail = os.path.split(file) - + # Build name name = tail.strip(".root") if name.find("kb_") != -1: @@ -281,7 +281,7 @@ if __name__ == '__main__': i = name.find("hours") if i != -1: name = name[0:i+1] - + # Build output filenames prfname = name + ".dat" bkgname = "bkg_"+name+".txt" diff --git a/cscripts/lsmodel.py b/cscripts/lsmodel.py index 6cfd031..0ab020c 100755 --- a/cscripts/lsmodel.py +++ b/cscripts/lsmodel.py @@ -39,7 +39,7 @@ class lsmodel(gammalib.GApplication): # Set name self.name = "lsmodel" self.version = "0.1.0" - + # Make sure that parfile exists file = self.parfile() @@ -57,7 +57,7 @@ class lsmodel(gammalib.GApplication): # Return return - + def __del__(self): """ Destructor. @@ -65,7 +65,7 @@ class lsmodel(gammalib.GApplication): # Write separator into logger if self.logTerse(): self.log("\n") - + # Return return @@ -76,13 +76,13 @@ class lsmodel(gammalib.GApplication): """ # Set parfile name parfile = self.name+".par" - + try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - + # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("cntmap","f","a","cntmap.fits","","","Counts map")) @@ -94,10 +94,10 @@ class lsmodel(gammalib.GApplication): pars.append(gammalib.GApplicationPar("outfile","f","a","model_cube.fits","","","Output file name")) pars.append_standard() pars.save(parfile) - + # Return return - + def get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -114,17 +114,17 @@ class lsmodel(gammalib.GApplication): # Set some fixed parameters self.m_log = False # Logging in client tools self.m_debug = False # Debugging in client tools - + # Return return - + def execute(self): """ Execute the script. """ # Run the script self.run() - + # Return return @@ -138,12 +138,12 @@ class lsmodel(gammalib.GApplication): # Get parameters self.get_parameters() - + # Write input parameters into logger if self.logTerse(): self.log_parameters() self.log("\n") - + # Write header if self.logTerse(): self.log("\n") @@ -194,10 +194,9 @@ if __name__ == '__main__': """ # Create instance of application app = lsmodel(sys.argv) - + # Open logfile app.logFileOpen() - + # Execute application app.execute() - \ No newline at end of file diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index ca125c8..36f995b 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -30,7 +30,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, binsz=0.05, npix=200, proj="TAN", coord="GAL"): """ Simulate events for all observations in the container. - + Parameters: obs - Observation container Keywords: @@ -42,12 +42,12 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, binsz - Pixel size for binned simulation (deg/pixel) npix - Number of pixels in X and Y for binned simulation """ - + # Allocate ctobssim application and set parameters sim = ctools.ctobssim(obs) sim["seed"] = seed sim["edisp"] = edisp - + # Optionally open the log file if log: sim.logFileOpen() @@ -58,16 +58,16 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, # Set chatter level sim["chatter"] = chatter - + # Run ctobssim application. This will loop over all observations in the # container and simulation the events for each observation. Note that # events are not added together, they still apply to each observation # separately. sim.run() - + # Binned option? if nbins > 0: - + # Determine common energy boundaries for observations emin = None emax = None @@ -82,7 +82,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, emax = run_emax elif run_emax > emax: emax = run_emax - + # Allocate ctbin application and set parameters bin = ctools.ctbin(sim.obs()) bin["ebinalg"] = "LOG" @@ -95,7 +95,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, bin["binsz"] = binsz bin["coordsys"] = coord bin["proj"] = proj - + # Optionally open the log file if log: bin.logFileOpen() @@ -110,22 +110,22 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, # Run ctbin application. This will loop over all observations in # the container and bin the events in counts maps bin.run() - + # Make a deep copy of the observation that will be returned # (the ctbin object will go out of scope one the function is # left) obs = bin.obs().copy() else: - + # Make a deep copy of the observation that will be returned # (the ctobssim object will go out of scope one the function is # left) obs = sim.obs().copy() - + # Delete the simulation del sim - + # Return observation container return obs @@ -136,7 +136,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, def fit(obs, log=False, debug=False, chatter=2, edisp=False): """ Perform maximum likelihood fitting of observations in the container. - + Parameters: obs - Observation container Keywords: @@ -147,24 +147,24 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False): """ # Allocate ctlike application like = ctools.ctlike(obs) - + # Optionally open the log file if log: like.logFileOpen() - + # Optionally switch-on debugging model if debug: like["debug"] = True # Set chatter level like["chatter"] = chatter - + # Optionally apply energy dispersion like["edisp"] = edisp # Run ctlike application. like.run() - + # Return observations return like @@ -175,7 +175,7 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False): def cterror(obs, srcname, log=False, debug=False, chatter=2): """ Perform maximum likelihood fitting of observations in the container. - + Parameters: obs - Observation container srcname - Source name @@ -190,11 +190,11 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2): # Set cterror parameters error["srcname"] = srcname - + # Optionally open the log file if log: error.logFileOpen() - + # Optionally switch-on debugging model if debug: error["debug"] = True @@ -204,7 +204,7 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2): # Run cterror application. error.run() - + # Return observations return error @@ -218,7 +218,7 @@ def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ """ Creates a counts map by combining the events of all observations. The counts map will be a summed map over all energies. - + Parameters: obs - Observation container Keywords: @@ -236,25 +236,25 @@ def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ # Set maximum pixel number maxpixel = nxpix * nypix - + # Fill all observations for run in obs: - + # Loop over all events for event in run.events(): - + # Determine sky pixel skydir = event.dir().dir() pixel = map.dir2inx(skydir) - + # Set pixel if pixel < maxpixel: map[pixel] += event.counts() - + # Save sky map. The clobber flag is set to True, so any existing FITS # file will be overwritten. map.save(outname, True) - + # Return counts map return map @@ -269,7 +269,7 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ Make model map for a given reference energy by combining all observations. The model map will be evaluated for a given reference energy 'eref' and will be given in units of [counts/(sr MeV s)]. - + Parameters: obs - Observation container Keywords: @@ -285,7 +285,7 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ """ # Allocate model map map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) - + # Set reference energy, time and direction. The time is not initialised and is # in fact not used (as the IRF is assumed to be time independent for now). # The sky direction is set later using the pixel values. @@ -293,33 +293,33 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ time = gammalib.GTime() instdir = gammalib.GCTAInstDir() energy.TeV(eref) - + # Loop over all map pixels for pixel in range(map.npix()): - + # Get sky direction skydir = map.inx2dir(pixel) instdir.dir(skydir) - + # Create event atom for map pixel atom = gammalib.GCTAEventAtom() atom.dir(instdir) atom.energy(energy) atom.time(time) - + # Initialise model value value = 0.0 - + # Loop over all observations for run in obs: value += obs.models().eval(atom, run) - + # Set map value map[pixel] = value - + # Save sky map map.save(outname, True) - + # Return model map return map @@ -335,15 +335,15 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ """ # Print warning print("Warning: obsutils.set is obsolete, use obsutils.set_obs instead.") - + # Call new function obs_cta = set_obs(pntdir, tstart=tstart, duration=duration, deadc=deadc, \ emin=emin, emax=emax, rad=rad, \ irf=irf, caldb=caldb) - + # Return CTA observation return obs_cta - + # ======================= # # Set one CTA observation # @@ -355,7 +355,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ Returns a single CTA observation containing an empty CTA event list. By looping over this function you can add CTA observations to the observation container. - + Parameters: pntdir - Pointing direction [GSkyDir] Keywords: @@ -379,23 +379,23 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ db.rootdir(caldb) else: db.open("cta", caldb) - + # Set pointing direction pnt = gammalib.GCTAPointing() pnt.dir(pntdir) obs_cta.pointing(pnt) - + # Set ROI roi = gammalib.GCTARoi() instdir = gammalib.GCTAInstDir() instdir.dir(pntdir) roi.centre(instdir) roi.radius(rad) - + # Set GTI gti = gammalib.GGti() gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) - + # Set energy boundaries ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ gammalib.GEnergy(emax, "TeV")) @@ -406,16 +406,16 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ events.gti(gti) events.ebounds(ebounds) obs_cta.events(events) - + # Set instrument response obs_cta.response(irf, db) - + # Set ontime, livetime, and deadtime correction factor obs_cta.ontime(duration) obs_cta.livetime(duration*deadc) obs_cta.deadc(deadc) obs_cta.id(id) - + # Return CTA observation return obs_cta @@ -434,7 +434,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ over other observation proposerties, such as duration, deadtime correction, energy range, etc. If an optional keyword is not specified, the function keyword is used instead. - + Parameters: obsdeflist - Observation definition list [{'ra': x.xx, 'dec': x.xx}] The directory can take the following optional keywords: @@ -470,7 +470,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ # Set pointing direction pntdir = gammalib.GSkyDir() pntdir.radec_deg(obsdef['ra'], obsdef['dec']) - + # Set duration (use default if not found in definition list) if 'duration' in obsdef: obs_duration = obsdef['duration'] @@ -512,7 +512,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ obs_irf = obsdef['irf'] else: obs_irf = irf - + # Generate identifier string id = "%6.6d" % obs_id @@ -520,14 +520,14 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \ deadc=obs_deadc, emin=obs_emin, emax=obs_emax, \ rad=obs_rad, irf=obs_irf, caldb=obs_caldb, id=id) - + # Append to container obs.append(obs_cta) # Update start time and identifier obs_start += obs_duration obs_id += 1 - + # Return observation container return obs @@ -543,7 +543,7 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): pattern - Observation pattern. Possible options are: - "single": single pointing - "four": four pointings 'offset' around pattern centre - + Keywords: ra - Right Ascension of pattern centre [deg] (default: 83.6331) dec - Declination of pattern centre [deg] (default: 22.0145) @@ -551,7 +551,7 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): """ # Initialise observation definition list obsdeflist = [] - + # Add patterns if pattern == "single": obsdef = {'ra': ra, 'dec': dec} @@ -560,7 +560,7 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): # Set pattern centre centre = gammalib.GSkyDir() centre.radec_deg(ra, dec) - + # Append pointings for phi in [0.0, 90.0, 180.0, 270.0]: pntdir = centre.copy() @@ -568,6 +568,6 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): obsdeflist.append({'ra': pntdir.ra_deg(), 'dec': pntdir.dec_deg()}) else: print("Warning: Observation pattern '"+str(pattern)+"' not recognized.") - + # Return observation definition list return obsdeflist diff --git a/cscripts/tsmerge.py b/cscripts/tsmerge.py index 8372c5e..2d1af5d 100644 --- a/cscripts/tsmerge.py +++ b/cscripts/tsmerge.py @@ -25,27 +25,27 @@ def print_help(): # Maps class holding the different sky maps # created by cttsmap class maps: - + # Constructor def __init__(self,fitsfile): - + self.filename = fitsfile if not os.path.isfile(fitsfile): self.error("Maps instance cannot be constructed: File \""+self.filename+"\" does not exist") - + # open FITS file fits = gammalib.GFits(fitsfile) - + # get TS map self.tsmap = gammalib.GSkymap() self.tsmap.read(fits[0]) - + # Get status map if not fits.contains("STATUS MAP"): self.error("No extension \"STATUS MAP\" found in \""+self.filename+"\"") self.statusmap = gammalib.GSkymap() self.statusmap.read(fits["STATUS MAP"]) - + # Get other maps self.maps = [] self.mapnames = [] @@ -72,30 +72,30 @@ class maps: check = False break return check - + # Save maps to outfile def save(self,outfile): - + fits = gammalib.GFits() self.tsmap.write(fits) - + for i in range(len(self.maps)): self.maps[i].write(fits) - + for i in range(len(self.mapnames)): fits[i+1].extname(self.mapnames[i]) - + if not self.check_status(): self.statusmap.write(fits) fits[fits.size()-1].extname("STATUS MAP") fits.saveto(outfile,True) - + # Add map instances to each other # in this way the TS maps get merged def add(self,addmaps): if not len(addmaps) == len(self.maps): self.error("Maps cannot be added: Maps from \""+self.filename+"\" do not match from \""+addmaps.filename+"\"") - + # Loop over bins for i in range(self.tsmap.npix()): if addmaps.statusmap[i] > 0.5: @@ -107,8 +107,8 @@ class maps: def error(self,message): print("*** ERROR ***: "+message) sys.exit(-1) - - + + def pars_from_argv(argv): if argv is None: @@ -139,10 +139,10 @@ def get_kwargs(argv): if s.count('=') == 1: key, value = s.split('=', 1) kwargs[key] = value - + else: print_help() - + return kwargs if __name__ == "__main__": @@ -158,10 +158,10 @@ if __name__ == "__main__": allfiles = open(files).read().splitlines() else: allfiles = glob.glob(files) - + # Initialise files to be merged mergefiles = [] - + # Test files for the entry status map # and add them to the file list for fitsfile in allfiles: @@ -169,14 +169,13 @@ if __name__ == "__main__": if fits.contains("STATUS MAP"): mergefiles.append(fitsfile) fits.close() - + # create a maps instance from the first file finalmaps = maps(mergefiles[0]) - + # add other files to the maps instance for fitsfile in mergefiles[1:]: finalmaps.add(maps(fitsfile)) - + # save the filled maps to file finalmaps.save(outfile=outfile) - diff --git a/examples/benchmark_cta_analysis.py b/examples/benchmark_cta_analysis.py index 96f827a..55fb78e 100755 --- a/examples/benchmark_cta_analysis.py +++ b/examples/benchmark_cta_analysis.py @@ -88,7 +88,7 @@ def unbinned_pipeline(duration): tstop = time.clock() telapsed = tstop - tstart tctlike = tstop - tctlike - + # Return return telapsed, tctlike @@ -161,7 +161,7 @@ def binned_pipeline(duration): tstop = time.clock() telapsed = tstop - tstart tctlike = tstop - tctlike - + # Return return telapsed, tctlike @@ -293,7 +293,7 @@ def stacked_pipeline(duration): tstop = time.clock() telapsed = tstop - tstart tctlike = tstop - tctlike - + # Return return telapsed, tctlike @@ -344,7 +344,7 @@ if __name__ == '__main__': # Increment duration duration *= 2.0 - + # Optionally plot results if has_matplotlib: plt.figure(1) @@ -359,4 +359,3 @@ if __name__ == '__main__': plt.ylabel("CPU time (seconds)") plt.legend(loc="lower right") plt.show() - \ No newline at end of file diff --git a/examples/make_pointings.py b/examples/make_pointings.py index 259e06a..9005284 100755 --- a/examples/make_pointings.py +++ b/examples/make_pointings.py @@ -330,7 +330,7 @@ def write_obsdef(filename, obsdef): # Return return - + # ======================== # # Main routine entry point # @@ -341,7 +341,7 @@ if __name__ == '__main__': """ # Initialise flags need_help = False - + # Test for command line arguments if (len(sys.argv) > 1): if sys.argv[1] == "-h": diff --git a/examples/make_pull_at_sensitivity_limit.py b/examples/make_pull_at_sensitivity_limit.py index 937a20c..9960604 100755 --- a/examples/make_pull_at_sensitivity_limit.py +++ b/examples/make_pull_at_sensitivity_limit.py @@ -33,10 +33,10 @@ from cscripts import cspull # Try importing processing module has_processing = False try: - import processing - has_processing = True + import processing + has_processing = True except: - pass + pass # ======================== # @@ -44,170 +44,170 @@ except: # ======================== # def create_pull(loge, emin, emax, models, ntrials=100, duration=180000.0, enumbins=0, log=False): - """ - Create pull distribution - - Parameters: - emin - Minimum energy (TeV) - emax - Maximum energy (TeV) - Keywords: - log - Create log file(s) - """ - # Generate output filename - outfile = "pull_"+loge+".dat" - - # Setup cspull tool - pull = cspull() - pull.models(models) - pull["outfile"] = outfile - pull["ntrials"] = ntrials - pull["caldb"] = "prod2" - pull["irf"] = "South_50h" - pull["ra"] = 83.6331 - pull["dec"] = 22.0145 - pull["emin"] = emin - pull["emax"] = emax - pull["enumbins"] = enumbins - pull["duration"] = duration - - # Optionally open the log file - if log: - pull.logFileOpen() - - # Run tool - pull.run() - - # Return - return + """ + Create pull distribution + + Parameters: + emin - Minimum energy (TeV) + emax - Maximum energy (TeV) + Keywords: + log - Create log file(s) + """ + # Generate output filename + outfile = "pull_"+loge+".dat" + + # Setup cspull tool + pull = cspull() + pull.models(models) + pull["outfile"] = outfile + pull["ntrials"] = ntrials + pull["caldb"] = "prod2" + pull["irf"] = "South_50h" + pull["ra"] = 83.6331 + pull["dec"] = 22.0145 + pull["emin"] = emin + pull["emax"] = emax + pull["enumbins"] = enumbins + pull["duration"] = duration + + # Optionally open the log file + if log: + pull.logFileOpen() + + # Run tool + pull.run() + + # Return + return # ============ # # Create model # # ============ # def create_model(flux, index=-2.48, fitidx=False): - """ - Add standard CTA background model to observations container. - We use a simple power law here, scaled to Konrad's E configuration - performance table. - - Parameters: - flux - Flux in Crab units - """ - # Define background model - bgd_radial = gammalib.GCTAModelRadialGauss(3.0) - bgd_spectrum = gammalib.GModelSpectralPlaw(61.8, -1.85) - bgd_spectrum["Prefactor"].scale(1.0e-6) - bgd_spectrum["PivotEnergy"].value(1.0) - bgd_spectrum["PivotEnergy"].scale(1.0e6) - if fitidx: - bgd_spectrum["Index"].free() - else: - bgd_spectrum["Index"].fix() - bgd_model = gammalib.GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum) - bgd_model.name("Background") - bgd_model.instruments("CTA") - - # Define source spectrum - location = gammalib.GSkyDir() - location.radec_deg(83.6331, 22.0145) - src_spatial = gammalib.GModelSpatialPtsrc(location) - src_spectrum = gammalib.GModelSpectralPlaw(flux, index) - src_spectrum["Prefactor"].scale(5.7e-16) - src_spectrum["PivotEnergy"].value(0.3) - src_spectrum["PivotEnergy"].scale(1.0e6) - if fitidx: - src_spectrum["Index"].free() - else: - src_spectrum["Index"].fix() - src_model = gammalib.GModelPointSource(src_spatial, src_spectrum) - src_model.name("Test") - - # Add models to container - models = gammalib.GModels() - models.append(bgd_model) - models.append(src_model) - - # Return model container - return models + """ + Add standard CTA background model to observations container. + We use a simple power law here, scaled to Konrad's E configuration + performance table. + + Parameters: + flux - Flux in Crab units + """ + # Define background model + bgd_radial = gammalib.GCTAModelRadialGauss(3.0) + bgd_spectrum = gammalib.GModelSpectralPlaw(61.8, -1.85) + bgd_spectrum["Prefactor"].scale(1.0e-6) + bgd_spectrum["PivotEnergy"].value(1.0) + bgd_spectrum["PivotEnergy"].scale(1.0e6) + if fitidx: + bgd_spectrum["Index"].free() + else: + bgd_spectrum["Index"].fix() + bgd_model = gammalib.GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum) + bgd_model.name("Background") + bgd_model.instruments("CTA") + + # Define source spectrum + location = gammalib.GSkyDir() + location.radec_deg(83.6331, 22.0145) + src_spatial = gammalib.GModelSpatialPtsrc(location) + src_spectrum = gammalib.GModelSpectralPlaw(flux, index) + src_spectrum["Prefactor"].scale(5.7e-16) + src_spectrum["PivotEnergy"].value(0.3) + src_spectrum["PivotEnergy"].scale(1.0e6) + if fitidx: + src_spectrum["Index"].free() + else: + src_spectrum["Index"].fix() + src_model = gammalib.GModelPointSource(src_spatial, src_spectrum) + src_model.name("Test") + + # Add models to container + models = gammalib.GModels() + models.append(bgd_model) + models.append(src_model) + + # Return model container + return models # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Create pull distributions for the energy bands and fluxes given in the - file specified in the argument. The file needs to be created using the - cssens script. - """ - # Get input arguments - usage = """ + """ + Create pull distributions for the energy bands and fluxes given in the + file specified in the argument. The file needs to be created using the +cssens script. + """ + # Get input arguments + usage = """ make_pull_at_sensitivity_limit filename [max_threads] Run cssens to create the input file. """ - if len(sys.argv) < 2 or len(sys.argv) > 3: - print(usage) - sys.exit() - - # Get sensitivity file - filename = sys.argv[1] - - # Set maximum number of threads - if len(sys.argv) > 2: - max_threads = int(sys.argv[2]) - else: - max_threads = 1 - - # Load sensitivity file - file = open(filename, 'r') - reader = csv.DictReader(file) - - # Loop over rows - for row in reader: - - # Get parameters - loge = float(row['loge']) - if loge < 0: - loge = "m"+str(abs(loge)) - else: - loge = "p"+str(abs(loge)) - emin = float(row['emin']) - emax = float(row['emax']) - flux = float(row['crab']) - - # Generate model - models = create_model(flux) - - # Processing support? - if has_processing: - - # Wait until one thread has finished - while len(processing.activeChildren()) >= max_threads: - time.sleep(60) - - # Set arguments - args = (loge, emin, emax, models) - kwargs = {} - - # Generate pull distribution - p = processing.Process(target=create_pull, args=args, kwargs=kwargs) - p.start() - print("Process emin=%.4f emax=%.4f started." % (emin, emax)) - - # Wait some time - time.sleep(1) - - # ... no - else: - create_pull(loge, emin, emax, models) - - # Processing support - if has_processing: - - # Wait until all threads finished - while len(processing.activeChildren()) > 0: - time.sleep(60) - - # Close file - file.close() + if len(sys.argv) < 2 or len(sys.argv) > 3: + print(usage) + sys.exit() + + # Get sensitivity file + filename = sys.argv[1] + + # Set maximum number of threads + if len(sys.argv) > 2: + max_threads = int(sys.argv[2]) + else: + max_threads = 1 + + # Load sensitivity file + file = open(filename, 'r') + reader = csv.DictReader(file) + + # Loop over rows + for row in reader: + + # Get parameters + loge = float(row['loge']) + if loge < 0: + loge = "m"+str(abs(loge)) + else: + loge = "p"+str(abs(loge)) + emin = float(row['emin']) + emax = float(row['emax']) + flux = float(row['crab']) + + # Generate model + models = create_model(flux) + + # Processing support? + if has_processing: + + # Wait until one thread has finished + while len(processing.activeChildren()) >= max_threads: + time.sleep(60) + + # Set arguments + args = (loge, emin, emax, models) + kwargs = {} + + # Generate pull distribution + p = processing.Process(target=create_pull, args=args, kwargs=kwargs) + p.start() + print("Process emin=%.4f emax=%.4f started." % (emin, emax)) + + # Wait some time + time.sleep(1) + + # ... no + else: + create_pull(loge, emin, emax, models) + + # Processing support + if has_processing: + + # Wait until all threads finished + while len(processing.activeChildren()) > 0: + time.sleep(60) + + # Close file + file.close() diff --git a/examples/make_spectrum.py b/examples/make_spectrum.py index bc3820e..885e432 100755 --- a/examples/make_spectrum.py +++ b/examples/make_spectrum.py @@ -22,10 +22,10 @@ import gammalib import ctools from cscripts import obsutils try: - import matplotlib.pyplot as plt - has_matplotlib = True + import matplotlib.pyplot as plt + has_matplotlib = True except: - has_matplotlib = False + has_matplotlib = False # =============== # @@ -67,7 +67,7 @@ def make_spectrum(): # Generate spectral points spectrum = obsutils.spectrum(sim.obs(), "Crab", ebounds) - + # Return spectrum return spectrum @@ -100,7 +100,7 @@ def plot_spectrum(spectrum): # Return return - + #==========================# # Main routine entry point # diff --git a/examples/make_survey.py b/examples/make_survey.py index f5fb536..5af2d9c 100755 --- a/examples/make_survey.py +++ b/examples/make_survey.py @@ -42,33 +42,33 @@ def plot_counts(obs): try: # Import matplotlib import matplotlib.pyplot as plt - + # Set legend fontsize params = {'legend.fontsize': 10} plt.rcParams.update(params) - + # Set plot styles styles = ['b-', 'g-', 'y-', 'n-'] - + # Dump header print("") print("Make plots (using matplotlib):") print("==============================") - + # Create figure 1 plt.figure(1,figsize=(12,6)) plt.subplots_adjust(hspace=.7) - + # Create subplot 1 plt.subplot(121) plt.title("Spectrum (summed over all pixels)") - + # Loop over observations for run in obs: - + # Get event list list = run.events() - + # Create energy axis ebounds = gammalib.GEbounds() emin = gammalib.GEnergy(0.1, "TeV") @@ -77,7 +77,7 @@ def plot_counts(obs): #emax.TeV(100.0) ebounds.setlog(emin, emax, 10) energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())] - + # Create spectrum print("Extract data:") counts = [0.0 for i in range(ebounds.size())] @@ -85,14 +85,14 @@ def plot_counts(obs): for atom in list: index = ebounds.index(atom.energy()) counts[index] = counts[index] + 1.0 - + # Create error bars error = [math.sqrt(c) for c in counts] - + # Plot spectrum plt.loglog(energy, counts, 'ro', label='data') #plt.errorbar(energy, counts, error, fmt=None, ecolor='r') - + # Extract models print("Extract models:") sum_model = [0.0 for i in range(ebounds.size())] @@ -107,26 +107,26 @@ def plot_counts(obs): sum_model[i] = sum_model[i] + model[i] #plt.loglog(energy, model, styles[k], label=m.name()) #plt.loglog(energy, sum_model, 'r-', label='total') - + # Put labels plt.xlabel("Energy (TeV)") plt.ylabel("Counts") plt.legend(loc="lower left") - + # Create subplot 2 plt.subplot(122) plt.title("Offset (summed over all energies)") - + # Set Crab direction crab = gammalib.GSkyDir() crab.radec_deg(83.63, 22.01) - + # Loop over observations for run in obs: - + # Get event list list = run.events() - + # Create offset histogram print("Extract data:") nx = 30 @@ -139,13 +139,13 @@ def plot_counts(obs): index = int(off2/doffset2) if index < nx: counts[index] = counts[index] + 1.0 - + # Create error bars error = [math.sqrt(c) for c in counts] - + # Plot distribution #plt.semilogy(offset2, counts, 'ro', label='data') - + # Extract models print("Extract models:") sum_model = [0.0 for i in range(nx)] @@ -164,15 +164,15 @@ def plot_counts(obs): #plt.plot(offset2, model, styles[k], label=m.name()) #plt.plot(offset2, sum_model, 'r-', label='total') #plt.ylim(ymin=0.1) - + # Put labels plt.xlabel("Offset (deg^2)") plt.ylabel("Counts") #plt.legend(loc="upper right") - + # Show counts spectra plt.show() - + except ImportError: print("Matplotlib is not (correctly) installed on your system.") @@ -186,26 +186,26 @@ def plot_counts(obs): def add_background_model(obs): """ Add standard CTA background model to observations container. - + We use a simple power law here, scaled to Konrad's E configuration performance table. The model needs still to be validated. """ # Recover models from observation models = obs.models() - + # Define background model bgd_radial = gammalib.GCTAModelRadialGauss(3.0) bgd_spectrum = gammalib.GModelSpectralPlaw(61.8e-6, -1.85, gammalib.GEnergy(1.0, "TeV")) bgd_model = gammalib.GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum) bgd_model.name("Background") bgd_model.instruments("CTA") - + # Add background model to container models.append(bgd_model) - + # Put container back in observation container obs.models(models) - + # Return observation container return obs @@ -220,7 +220,7 @@ def crab_spec(): """ # Set parameters spectrum = gammalib.GModelSpectralPlaw(5.7e-16, -2.48, gammalib.GEnergy(0.3, "TeV")) - + # Return spectrum return spectrum @@ -234,13 +234,13 @@ def survey_single(): """ # Allocate observation container obs = gammalib.GObservations() - + # Set single pointing at galactic centre pntdir = gammalib.GSkyDir() pntdir.lb_deg(0.0, 0.0) run = obsutils.set_obs(pntdir) obs.append(run) - + # Define single point source with Crab flux at galactic centre center = gammalib.GSkyDir() center.lb_deg(0.0, 0.0) @@ -248,12 +248,12 @@ def survey_single(): point_spectrum = crab_spec() point = gammalib.GModelSky(point_spatial, point_spectrum) point.name('GC source') - + # Create model container models = gammalib.GModels() models.append(point) obs.models(models) - + # Return observation container return obs @@ -264,24 +264,24 @@ def survey_single(): def survey_gplane(lrange=10, lstep=2): """ Creates a single observation survey for test purposes. - + Keywords: lrange - Longitude range (integer deg) lstep - Longitude step size (integer deg) """ # Allocate observation container obs = gammalib.GObservations() - + # Loop over longitudes for l in range(-lrange,lrange+lstep,lstep): - + # Set pointing pntdir = gammalib.GSkyDir() pntdir.lb_deg(l, 0.0) run = obsutils.set_obs(pntdir) run.id(str(l)) obs.append(run) - + # Define single point source with Crab flux at galactic centre center = gammalib.GSkyDir() center.lb_deg(0.0, 0.0) @@ -289,12 +289,12 @@ def survey_gplane(lrange=10, lstep=2): point_spectrum = crab_spec() point = gammalib.GModelSky(point_spatial, point_spectrum) point.name('GC source') - + # Create model container models = gammalib.GModels() models.append(point) obs.models(models) - + # Return observation container return obs @@ -308,7 +308,7 @@ if __name__ == '__main__': """ # Initialise flags need_help = False - + # Test for command line arguments print(sys.argv[0]) if (len(sys.argv) > 1): @@ -322,42 +322,42 @@ if __name__ == '__main__': print("Usage: example_survey.py [OPTIONS]") print(" -h Display this usage message") sys.exit() - + # Dump header print("***********************") print("* Simulate CTA survey *") print("***********************") - + # Remove any existing result files list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")] for files in list: for file in files: os.remove(file) - + # Setup single observation survey #obs = survey_single() obs = survey_gplane() - + # Add background model obs = add_background_model(obs) - + # Simulate events print("Simulate events") obs = obsutils.sim(obs) - + # Make counts map print("Make counts map") cntmap = obsutils.cntmap(obs) - + # Fit observations print("Fit observations") like = obsutils.fit(obs) #print like.opt() #print like.obs().models() - + # Make model map print("Make model map (this step will take some time)") modmap = obsutils.modmap(obs) - + # Show fit results #plot_counts(like.obs()) diff --git a/examples/make_ts_distributions.py b/examples/make_ts_distributions.py index ca9d0ed..44b7c91 100755 --- a/examples/make_ts_distributions.py +++ b/examples/make_ts_distributions.py @@ -35,155 +35,155 @@ from cscripts import cstsdist # Try importing processing module has_processing = False try: - import processing - has_processing = True + import processing + has_processing = True except: - pass + pass # ====================== # # Create TS distribution # # ====================== # def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, - enumbins=0, log=False): - """ - Create TS distribution. - - Parameters: - loge - Logarithm of mean energy in TeV - emin - Minimum energy (TeV) - emax - Maximum energy (TeV) - Keywords: - log - Create log file(s) - """ - # Generate output filename - outfile = "ts_"+str(loge)+".dat" - - # Setup cstsdist tool - tsdist = cstsdist() - tsdist["inmodel"] = "$CTOOLS/share/models/crab.xml" - tsdist["srcname"] = "Crab" - tsdist["outfile"] = outfile - tsdist["ntrials"] = ntrials - tsdist["caldb"] = "prod2" - tsdist["irf"] = "South_50h" - tsdist["ra"] = 83.63 - tsdist["dec"] = 22.01 - tsdist["emin"] = emin - tsdist["emax"] = emax - tsdist["enumbins"] = enumbins - tsdist["tmin"] = 0 - tsdist["tmax"] = duration - tsdist["rad"] = 5.0 - tsdist["npix"] = 200 - tsdist["binsz"] = 0.05 - #tsdist["debug"].boolean(True) - - # Optionally open the log file - if log: - tsdist.logFileOpen() - - # Run tool - tsdist.run() - - # Return - return + enumbins=0, log=False): + """ + Create TS distribution. + + Parameters: + loge - Logarithm of mean energy in TeV + emin - Minimum energy (TeV) + emax - Maximum energy (TeV) + Keywords: + log - Create log file(s) + """ + # Generate output filename + outfile = "ts_"+str(loge)+".dat" + + # Setup cstsdist tool + tsdist = cstsdist() + tsdist["inmodel"] = "$CTOOLS/share/models/crab.xml" + tsdist["srcname"] = "Crab" + tsdist["outfile"] = outfile + tsdist["ntrials"] = ntrials + tsdist["caldb"] = "prod2" + tsdist["irf"] = "South_50h" + tsdist["ra"] = 83.63 + tsdist["dec"] = 22.01 + tsdist["emin"] = emin + tsdist["emax"] = emax + tsdist["enumbins"] = enumbins + tsdist["tmin"] = 0 + tsdist["tmax"] = duration + tsdist["rad"] = 5.0 + tsdist["npix"] = 200 + tsdist["binsz"] = 0.05 + #tsdist["debug"].boolean(True) + + # Optionally open the log file + if log: + tsdist.logFileOpen() + + # Run tool + tsdist.run() + + # Return + return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Create TS distribution in a number of energy bands. - """ - # Get input arguments - usage = "make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]" - if len(sys.argv) < 1: - print(usage) - sys.exit() - - # Set default parameters - ntrials = 100 - enumbins = 0 - duration = 180000.0 - max_threads = 1 - - # Parameter dictionnary - pars = [{'option': '-n', 'value': ntrials}, - {'option': '-e', 'value': enumbins}, - {'option': '-d', 'value': duration}, - {'option': '-m', 'value': max_threads}] - - # Gather parameters from command line - i = 1 - while i < len(sys.argv): - - # Search for options - for par in pars: - if sys.argv[i] == par['option']: - if len(sys.argv) > i+1: - i += 1 - try: - par['value'] = int(sys.argv[i]) - except: - print(usage) - sys.exit() - else: - print(usage) - sys.exit() - - # Next item - i += 1 - - # Recover parameters - ntrials = pars[0]['value'] - enumbins = pars[1]['value'] - duration = pars[2]['value'] - max_threads = pars[3]['value'] - - # Loop over energy bands. The energy bands are those that are also - # used for sensitivity computation. - for ieng in range(20): - - # Set energies - loge = -1.7 + ieng * 0.2 - emean = pow(10.0, loge) - emin = pow(10.0, loge-0.1) - emax = pow(10.0, loge+0.1) - if loge < 0: - loge = "m"+str(abs(loge)) - else: - loge = "p"+str(abs(loge)) - - # Processing support? - if has_processing: - - # Wait until one thread has finished - while len(processing.activeChildren()) >= max_threads: - time.sleep(10) - - # Set arguments - args = (loge, emin, emax) - kwargs = {'ntrials': ntrials, 'enumbins': enumbins, - 'duration': duration} - - # Generate pull distribution - p = processing.Process(target=create_ts, args=args, kwargs=kwargs) - p.start() - print("Process emin=%.4f emax=%.4f started." % (emin, emax)) - - # Wait a short time to allow process to start - time.sleep(1) - - # ... no - else: - create_ts(loge, emin, emax, ntrials=ntrials, enumbins=enumbins, - duration=duration) - - # Processing support - if has_processing: - - # Wait until all threads finished - while len(processing.activeChildren()) > 0: - time.sleep(10) + """ + Create TS distribution in a number of energy bands. + """ + # Get input arguments + usage = "make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]" + if len(sys.argv) < 1: + print(usage) + sys.exit() + + # Set default parameters + ntrials = 100 + enumbins = 0 + duration = 180000.0 + max_threads = 1 + + # Parameter dictionnary + pars = [{'option': '-n', 'value': ntrials}, + {'option': '-e', 'value': enumbins}, + {'option': '-d', 'value': duration}, + {'option': '-m', 'value': max_threads}] + + # Gather parameters from command line + i = 1 + while i < len(sys.argv): + + # Search for options + for par in pars: + if sys.argv[i] == par['option']: + if len(sys.argv) > i+1: + i += 1 + try: + par['value'] = int(sys.argv[i]) + except: + print(usage) + sys.exit() + else: + print(usage) + sys.exit() + + # Next item + i += 1 + + # Recover parameters + ntrials = pars[0]['value'] + enumbins = pars[1]['value'] + duration = pars[2]['value'] + max_threads = pars[3]['value'] + + # Loop over energy bands. The energy bands are those that are also + # used for sensitivity computation. + for ieng in range(20): + + # Set energies + loge = -1.7 + ieng * 0.2 + emean = pow(10.0, loge) + emin = pow(10.0, loge-0.1) + emax = pow(10.0, loge+0.1) + if loge < 0: + loge = "m"+str(abs(loge)) + else: + loge = "p"+str(abs(loge)) + + # Processing support? + if has_processing: + + # Wait until one thread has finished + while len(processing.activeChildren()) >= max_threads: + time.sleep(10) + + # Set arguments + args = (loge, emin, emax) + kwargs = {'ntrials': ntrials, 'enumbins': enumbins, + 'duration': duration} + + # Generate pull distribution + p = processing.Process(target=create_ts, args=args, kwargs=kwargs) + p.start() + print("Process emin=%.4f emax=%.4f started." % (emin, emax)) + + # Wait a short time to allow process to start + time.sleep(1) + + # ... no + else: + create_ts(loge, emin, emax, ntrials=ntrials, enumbins=enumbins, + duration=duration) + + # Processing support + if has_processing: + + # Wait until all threads finished + while len(processing.activeChildren()) > 0: + time.sleep(10) diff --git a/examples/pipeline_binned_disk.py b/examples/pipeline_binned_disk.py index 7c71c3d..89c3e4d 100755 --- a/examples/pipeline_binned_disk.py +++ b/examples/pipeline_binned_disk.py @@ -39,7 +39,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, caldb="prod2", irf="South_50h"): """ Returns an observation container. - + Keywords: pattern - Pointing pattern, either "single" or "four" ra - RA of pattern centre [deg] (default: 83.6331) @@ -58,7 +58,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, ra=ra, dec=dec, offset=offset) - + # Get observation container obs = obsutils.set_obs_list(obs_def_list, duration=duration, @@ -78,13 +78,13 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): """ Setup model for analysis. - + Keywords: model - Model Xml file """ # Append model from file to observation container obs.models(gammalib.GModels(model)) - + # Return observation container return obs @@ -170,7 +170,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0, like["irf"] = irf like["debug"] = True # Switch this always on for results in console like.execute() - + # Return return diff --git a/examples/pipeline_binned_mem.py b/examples/pipeline_binned_mem.py index 389787b..f85a5b0 100755 --- a/examples/pipeline_binned_mem.py +++ b/examples/pipeline_binned_mem.py @@ -39,7 +39,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, caldb="prod2", irf="South_50h"): """ Returns an observation container. - + Keywords: pattern - Pointing pattern, either "single" or "four" ra - RA of pattern centre [deg] (default: 83.6331) @@ -58,7 +58,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, ra=ra, dec=dec, offset=offset) - + # Get observation container obs = obsutils.set_obs_list(obs_def_list, duration=duration, @@ -78,13 +78,13 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): """ Setup model for analysis. - + Keywords: model - Model Xml file """ # Append model from file to observation container obs.models(gammalib.GModels(model)) - + # Return observation container return obs @@ -144,7 +144,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0, like = ctools.ctlike(obs) like["debug"] = True # Switch this always on for results in console like.run() - + # Return return diff --git a/examples/pipeline_stacked_disk.py b/examples/pipeline_stacked_disk.py index e0876fa..3277caf 100755 --- a/examples/pipeline_stacked_disk.py +++ b/examples/pipeline_stacked_disk.py @@ -39,7 +39,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, caldb="prod2", irf="South_50h"): """ Returns an observation container. - + Keywords: pattern - Pointing pattern, either "single" or "four" ra - RA of pattern centre [deg] (default: 83.6331) @@ -58,7 +58,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, ra=ra, dec=dec, offset=offset) - + # Get observation container obs = obsutils.set_obs_list(obs_def_list, duration=duration, @@ -78,13 +78,13 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): """ Setup model for analysis. - + Keywords: model - Model Xml file """ # Append model from file to observation container obs.models(gammalib.GModels(model)) - + # Return observation container return obs @@ -215,7 +215,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, like["irf"] = irf like["debug"] = True # Switch this always on for results in console like.execute() - + # Return return diff --git a/examples/pipeline_stacked_mem.py b/examples/pipeline_stacked_mem.py index 7e6477d..555055a 100755 --- a/examples/pipeline_stacked_mem.py +++ b/examples/pipeline_stacked_mem.py @@ -39,7 +39,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, caldb="prod2", irf="South_50h"): """ Returns an observation container. - + Keywords: pattern - Pointing pattern, either "single" or "four" ra - RA of pattern centre [deg] (default: 83.6331) @@ -58,7 +58,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, ra=ra, dec=dec, offset=offset) - + # Get observation container obs = obsutils.set_obs_list(obs_def_list, duration=duration, @@ -78,13 +78,13 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): """ Setup model for analysis. - + Keywords: model - Model Xml file """ # Append model from file to observation container obs.models(gammalib.GModels(model)) - + # Return observation container return obs @@ -194,7 +194,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, like = ctools.ctlike(bin.obs()) like["debug"] = True # Switch this always on for results in console like.run() - + # Return return diff --git a/examples/pipeline_unbinned_disk.py b/examples/pipeline_unbinned_disk.py index b81e20d..55fae66 100755 --- a/examples/pipeline_unbinned_disk.py +++ b/examples/pipeline_unbinned_disk.py @@ -39,7 +39,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, caldb="prod2", irf="South_50h"): """ Returns an observation container. - + Keywords: pattern - Pointing pattern, either "single" or "four" ra - RA of pattern centre [deg] (default: 83.6331) @@ -58,7 +58,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, ra=ra, dec=dec, offset=offset) - + # Get observation container obs = obsutils.set_obs_list(obs_def_list, duration=duration, @@ -78,13 +78,13 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): """ Setup model for analysis. - + Keywords: model - Model Xml file """ # Append model from file to observation container obs.models(gammalib.GModels(model)) - + # Return observation container return obs @@ -145,7 +145,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, like["irf"] = irf like["debug"] = True # Switch this always on for results in console like.execute() - + # Return return diff --git a/examples/pipeline_unbinned_mem.py b/examples/pipeline_unbinned_mem.py index d34c09e..392c21c 100755 --- a/examples/pipeline_unbinned_mem.py +++ b/examples/pipeline_unbinned_mem.py @@ -39,7 +39,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, caldb="prod2", irf="South_50h"): """ Returns an observation container. - + Keywords: pattern - Pointing pattern, either "single" or "four" ra - RA of pattern centre [deg] (default: 83.6331) @@ -58,7 +58,7 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, ra=ra, dec=dec, offset=offset) - + # Get observation container obs = obsutils.set_obs_list(obs_def_list, duration=duration, @@ -78,13 +78,13 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): """ Setup model for analysis. - + Keywords: model - Model Xml file """ # Append model from file to observation container obs.models(gammalib.GModels(model)) - + # Return observation container return obs @@ -130,7 +130,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, like = ctools.ctlike(select.obs()) like["debug"] = True # Switch this always on for results in console like.run() - + # Return return diff --git a/examples/show_butterfly.py b/examples/show_butterfly.py index 1f522a7..ad2774b 100755 --- a/examples/show_butterfly.py +++ b/examples/show_butterfly.py @@ -25,46 +25,46 @@ try: import matplotlib.pyplot as plt except: sys.exit("This script needs matplotlib") - + # =============== # # Script entry # # =============== # if __name__ == "__main__": - + # check for given butterfly file if not len(sys.argv) == 2: sys.exit("Usage: show_butterfly.py butterfly.txt") - + # Read given butterfly file filename = sys.argv[1] csv = gammalib.GCsv(filename) - + # initialise arrays to be filled butterfly_x = [] butterfly_y = [] - + line_x = [] line_y = [] - + # loop over rows of the file nrows = csv.nrows() for row in range(nrows): # Compute upper edge of confidence band butterfly_x.append(csv.real(row,0)) butterfly_y.append(csv.real(row,1)+csv.real(row,2)) - + # Set line values line_x.append(csv.real(row,0)) line_y.append(csv.real(row,1)) - + # Loop over the rows backwards to compute the lower edge # of the confidence band for row in range(nrows): index = nrows - 1 - row - + butterfly_x.append(csv.real(index,0)) butterfly_y.append(csv.real(index,1)-csv.real(index,2)) - + # plot the butterfly and spectral line plt.figure() plt.loglog() @@ -73,4 +73,4 @@ if __name__ == "__main__": plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5) plt.xlabel("Energy [MeV]") plt.ylabel(r"dN/dE [MeV$^{-1}$ s$^{-1}$ cm$^{-2}$]") - plt.show() \ No newline at end of file + plt.show() diff --git a/examples/show_model.py b/examples/show_model.py index 5372c09..4fee000 100755 --- a/examples/show_model.py +++ b/examples/show_model.py @@ -35,7 +35,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100): """ # Get spectral component spectrum = model.spectral() - + # Setup energy axis e_min = gammalib.GEnergy(emin, "TeV") e_max = gammalib.GEnergy(emax, "TeV") @@ -61,7 +61,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100): if min < 0.0001*max: min = 0.0001*max min = math.pow(10.0, math.log10(min)-1) - + # Show spectrum plt.figure() @@ -73,11 +73,11 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100): plt.ylabel(r"dN/dE [ph s$^{-1}$ cm$^{-2}$ MeV$^{-1}$]") plt.ylim([min,max]) plt.show() - + # Return return - - + + # ============= # # Script entry # # ============= # @@ -100,6 +100,6 @@ if __name__ == "__main__": # Extract relevant model model = models[index] - + # Plot spectrum plot_spectrum(model) diff --git a/examples/show_pha.py b/examples/show_pha.py index e9a3450..547752f 100755 --- a/examples/show_pha.py +++ b/examples/show_pha.py @@ -53,7 +53,7 @@ def show_pha(pha): for i in range(pha.size()): counts[i] = pha[i] error = [math.sqrt(c) for c in counts] - + # Plot data if channels: plt.semilogy(energy, counts, 'ro') @@ -70,11 +70,11 @@ def show_pha(pha): # Show plot plt.show() - + # Return #return - - + + # ======================== # # Main routine entry point # # ======================== # @@ -107,4 +107,3 @@ if __name__ == '__main__': # Show PHA spectrum if has_matplotlib: show_pha(pha) - \ No newline at end of file diff --git a/examples/show_pull_histogram.py b/examples/show_pull_histogram.py index 042384d..3cf3707 100755 --- a/examples/show_pull_histogram.py +++ b/examples/show_pull_histogram.py @@ -39,85 +39,84 @@ import csv # Read pull distribution # # ====================== # def read_pull(filename, parname): - """ - Read pull distribution. - """ - # Initialise list - values = [] - - # Open reader - reader = csv.reader(open(filename, 'r'), delimiter=',') - - # Read rows - first = True - index = -1 - for row in reader: - - # Get column index if first row - if first: - try: - index = row.index(parname) - except: - sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n') - for p in row: - sys.stdout.write(' "'+p+'"\n') - raise NameError(parname) - - # Handle data rows - else: - values.append(float(row[index])) - - # Flag that first row has been passed - first = False - - # Create numpy array - a = np.array(values) - - # Return array - return a + """ + Read pull distribution. + """ + # Initialise list + values = [] + + # Open reader + reader = csv.reader(open(filename, 'r'), delimiter=',') + + # Read rows + first = True + index = -1 + for row in reader: + + # Get column index if first row + if first: + try: + index = row.index(parname) + except: + sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n') + for p in row: + sys.stdout.write(' "'+p+'"\n') + raise NameError(parname) + + # Handle data rows + else: + values.append(float(row[index])) + + # Flag that first row has been passed + first = False + + # Create numpy array + a = np.array(values) + + # Return array + return a # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Display pull histogram generated by cspull using matplotlib. - """ - # Print usage information - usage = "Usage: show_pull_histogram filename parname [bins]" - if len(sys.argv) < 3: - sys.stdout.write(usage+"\n") - sys.exit() - - # Extract parameters - filename = sys.argv[1] - parname = sys.argv[2] - if len(sys.argv) > 3: - nbins = int(sys.argv[3]) - else: - nbins = 50 - - # Read values from CSV file - values = read_pull(filename, parname) - - # Create histogram - n, bins, patches = plt.hist(values, nbins, range=[-4.0,4.0], - normed=True, facecolor='green') - - # Create expected distribution - y = mlab.normpdf(bins, 0.0, 1.0) - l = plt.plot(bins, y, 'r-', linewidth=2) - - # Set parname for plotting - name = parname - - # Set plot - plt.xlabel('Pull ('+name+')') - plt.ylabel('Arbitrary units') - plt.title(name) - plt.grid(True) - - # Show histogram - plt.show() - \ No newline at end of file + """ + Display pull histogram generated by cspull using matplotlib. + """ + # Print usage information + usage = "Usage: show_pull_histogram filename parname [bins]" + if len(sys.argv) < 3: + sys.stdout.write(usage+"\n") + sys.exit() + + # Extract parameters + filename = sys.argv[1] + parname = sys.argv[2] + if len(sys.argv) > 3: + nbins = int(sys.argv[3]) + else: + nbins = 50 + + # Read values from CSV file + values = read_pull(filename, parname) + + # Create histogram + n, bins, patches = plt.hist(values, nbins, range=[-4.0,4.0], + normed=True, facecolor='green') + + # Create expected distribution + y = mlab.normpdf(bins, 0.0, 1.0) + l = plt.plot(bins, y, 'r-', linewidth=2) + + # Set parname for plotting + name = parname + + # Set plot + plt.xlabel('Pull ('+name+')') + plt.ylabel('Arbitrary units') + plt.title(name) + plt.grid(True) + + # Show histogram + plt.show() diff --git a/examples/show_response.py b/examples/show_response.py index e6a5aaa..3ba9414 100755 --- a/examples/show_response.py +++ b/examples/show_response.py @@ -32,58 +32,58 @@ import matplotlib.pyplot as plt # Compute Eq. (17) of Li & Ma (1983) # # ================================== # def sigma_lima(Non, Noff, alpha=0.2): - """ - Compute Eq. (17) of Li & Ma (1983). - - Parameters: - Non - Number of on counts - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure - """ - # Compute sensitivity - alpha1 = alpha + 1.0 - sum = Non + Noff - arg1 = Non / sum - arg2 = Noff / sum - term1 = Non * math.log((alpha1/alpha)*arg1) - term2 = Noff * math.log(alpha1*arg2) - sigma = math.sqrt(2.0 * (term1 + term2)) - - # Return sensitivity - return sigma + """ + Compute Eq. (17) of Li & Ma (1983). + + Parameters: + Non - Number of on counts + Noff - Number of off counts + Keywords: + alpha - Ratio of on-to-off exposure + """ + # Compute sensitivity + alpha1 = alpha + 1.0 + sum = Non + Noff + arg1 = Non / sum + arg2 = Noff / sum + term1 = Non * math.log((alpha1/alpha)*arg1) + term2 = Noff * math.log(alpha1*arg2) + sigma = math.sqrt(2.0 * (term1 + term2)) + + # Return sensitivity + return sigma # ======================================== # # Solve Eq. (17) of Li & Ma (1983) for Non # # ======================================== # def Non_lima(sigma, Noff, alpha=0.2): - """ - Solve Eq. (17) of Li & Ma (1983) for Non - - Parameters: - sigma - Required significance - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure - """ - # Set initial guess for Non - S = 25.0 - Non = S + alpha*Noff - - # Compute initial sigma estimate - s = sigma_lima(Non, Noff, alpha=alpha) - delta = s - sigma - - # Iterate - while abs(delta) > 1.0e-6: - S *= sigma/s - Non = S + alpha*Noff - s = sigma_lima(Non, Noff, alpha=alpha) - delta = s - sigma - - # Return - return Non + """ + Solve Eq. (17) of Li & Ma (1983) for Non + + Parameters: + sigma - Required significance + Noff - Number of off counts + Keywords: + alpha - Ratio of on-to-off exposure + """ + # Set initial guess for Non + S = 25.0 + Non = S + alpha*Noff + + # Compute initial sigma estimate + s = sigma_lima(Non, Noff, alpha=alpha) + delta = s - sigma + + # Iterate + while abs(delta) > 1.0e-6: + S *= sigma/s + Non = S + alpha*Noff + s = sigma_lima(Non, Noff, alpha=alpha) + delta = s - sigma + + # Return + return Non # ======================= # @@ -95,7 +95,7 @@ def show_one_effective_area(rsp, name, color="r"): """ # Set to figure 1 plt.figure(1) - + # Generate logE(TeV) vector nbins = 20 logE = [-1.7+0.2*i for i in range(nbins)] @@ -130,7 +130,7 @@ def show_one_background_rate(rsp, name, color="r"): """ # Set to figure 2 plt.figure(2) - + # Generate logE(TeV) vector nbins = 20 logE = [-1.7+0.2*i for i in range(nbins)] @@ -142,7 +142,7 @@ def show_one_background_rate(rsp, name, color="r"): bgrate[i] = rsp.background()(logE[i], 0.0, 0.0) if bgrate[i] < 1.0e-10: bgrate[i] = 0.0 - + # Plot data plt.loglog(E, bgrate, color+'-', label=name) plt.loglog(E, bgrate, color+'o') @@ -171,7 +171,7 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig # Set to figure 3 plt.figure(3) - + # Generate logE(TeV) vector nbins = 20 logE = [-1.7+0.2*i for i in range(nbins)] @@ -200,7 +200,7 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig plaw = gammalib.GModelSpectralPlaw(1.0e-6, -2.6, epivot) conv = TeV2erg*E[i]*E[i]/plaw.flux(emin, emax) flux[i] = conv * src_counts / (duration * aeff*0.68) - + # Plot data plt.loglog(E, flux, color+'-', label=name) plt.loglog(E, flux, color+'o') @@ -243,11 +243,11 @@ def show_one_response(rspname, dbname, name, rootdir=None, color="r"): # Show sensitivity show_one_sensitivity(rsp, name, color=color) - + # Return return - - + + # ======================== # # Main routine entry point # # ======================== # @@ -291,7 +291,7 @@ if __name__ == '__main__': rootdir = rsp['rootdir'] else: rootdir = None - + # Show response show_one_response(rspname, dbname, name, rootdir=rootdir, color=color) diff --git a/examples/show_spectrum.py b/examples/show_spectrum.py index d0b3111..8065a37 100755 --- a/examples/show_spectrum.py +++ b/examples/show_spectrum.py @@ -26,18 +26,18 @@ try: except: sys.exit("This script needs matplotlib") - + # ============ # # Script entry # # ============ # if __name__ == "__main__": """ """ - + # Check for given spectrum file if not len(sys.argv) == 2: sys.exit("Usage: show_spectrum.py spectrum.fits") - + # Read spectrum file filename = sys.argv[1] fits = gammalib.GFits(filename) @@ -49,7 +49,7 @@ if __name__ == "__main__": c_eflux = table["e_Flux"] c_ts = table["TS"] c_upper = table["UpperLimit"] - + # Initialise arrays to be filled energies = [] flux = [] @@ -60,7 +60,7 @@ if __name__ == "__main__": ul_ed_engs = [] ul_eu_engs = [] ul_flux = [] - + # Loop over rows of the file nrows = table.nrows() for row in range(nrows): @@ -69,7 +69,7 @@ if __name__ == "__main__": ts = c_ts.real(row) flx = c_flux.real(row) e_flx = c_eflux.real(row) - + # Switch if ts > 9.0 and e_flx < flx: @@ -88,7 +88,7 @@ if __name__ == "__main__": ul_flux.append(c_upper.real(row)) ul_ed_engs.append(c_ed.real(row)) ul_eu_engs.append(c_eu.real(row)) - + # Plot the spectrum plt.figure() plt.loglog() @@ -98,4 +98,3 @@ if __name__ == "__main__": plt.xlabel("Energy [TeV]") plt.ylabel(r"dN/dE [erg cm$^{-2}$ s$^{-1}$]") plt.show() - \ No newline at end of file diff --git a/examples/show_ts_distribution.py b/examples/show_ts_distribution.py index 06de278..027b4c8 100755 --- a/examples/show_ts_distribution.py +++ b/examples/show_ts_distribution.py @@ -34,249 +34,249 @@ import math # Read TS distribution # # ==================== # def read_ts(filename, tsname="TS"): - """ - Read TS distribution values from CSV file. - """ - # Initialise list - values = [] - - # Open reader - reader = csv.reader(open(filename, 'r'), delimiter=',') - - # Read rows - first = True - index = -1 - for row in reader: - - # Get column index if first row - if first: - try: - index = row.index(tsname) - except: - sys.stdout.write('ERROR: Column "TS" not found in file\n') - sys.stdout.write(row+"\n") - raise NameError("TS") - - # Handle data rows - else: - values.append(float(row[index])) - - # Flag that first row has been passed - first = False - - # Create numpy array - a = np.array(values) - - # Return array - return a + """ + Read TS distribution values from CSV file. + """ + # Initialise list + values = [] + + # Open reader + reader = csv.reader(open(filename, 'r'), delimiter=',') + + # Read rows + first = True + index = -1 + for row in reader: + + # Get column index if first row + if first: + try: + index = row.index(tsname) + except: + sys.stdout.write('ERROR: Column "TS" not found in file\n') + sys.stdout.write(row+"\n") + raise NameError("TS") + + # Handle data rows + else: + values.append(float(row[index])) + + # Flag that first row has been passed + first = False + + # Create numpy array + a = np.array(values) + + # Return array + return a # =========== # # Compute erf # # =========== # def erf(x): - """ - Compute error function. - """ - # save the sign of x - sign = 1 - if x < 0: - sign = -1 - x = abs(x) - - # constants - a1 = 0.254829592 - a2 = -0.284496736 - a3 = 1.421413741 - a4 = -1.453152027 - a5 = 1.061405429 - p = 0.3275911 - - # A&S formula 7.1.26 - t = 1.0/(1.0 + p*x) - y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x) - - # Return result - return sign*y # erf(-x) = -erf(x) + """ + Compute error function. + """ + # save the sign of x + sign = 1 + if x < 0: + sign = -1 + x = abs(x) + + # constants + a1 = 0.254829592 + a2 = -0.284496736 + a3 = 1.421413741 + a4 = -1.453152027 + a5 = 1.061405429 + p = 0.3275911 + + # A&S formula 7.1.26 + t = 1.0/(1.0 + p*x) + y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x) + + # Return result + return sign*y # erf(-x) = -erf(x) # =============================== # # Display cumulative distribution # # =============================== # def dist_cdf(values, nbins, title): - """ - Displays cumulative TS distribution and overplots expectation. - """ - # Set range and adapt number of bins. We make sure that we have a - # histogram bin centred on 0 that will capture the TS=0 part. - max_value = max(values) - binsize = max_value/nbins - min_value = 0.5*binsize - while min_value > min(values): - min_value -= binsize - nbins += 1 - - # Create histogram - hist_x = [] - hist_y = [] - hist_expect = [] - norm = 0.5*len(values) - for i in range(nbins): - x_min = min_value+i*binsize - x_max = min_value+(i+1)*binsize - x_val = 0.5*(x_min+x_max) - y_val = 0.0 - for ts in values: - if ts > x_min: - y_val += 1.0 - hist_x.append(x_min) - hist_y.append(y_val) - hist_x.append(x_max) - hist_y.append(y_val) - if x_min > 0.0: - value = norm*erf(math.sqrt(0.5*x_min)) + norm - hist_expect.append(len(values)-value) - hist_expect.append(len(values)-value) - else: - hist_expect.append(norm) - hist_expect.append(norm) - - # Show expected distribution (half) - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label="simulations") - l = plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label="expectation") - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) - - # Set plot - plt.xlabel('TS') - plt.ylabel('Number of trials > TS') - plt.title(title) - plt.grid(True) - plt.legend(loc="upper right") - - # Show histogram - plt.show() - - # Return - return + """ + Displays cumulative TS distribution and overplots expectation. + """ + # Set range and adapt number of bins. We make sure that we have a + # histogram bin centred on 0 that will capture the TS=0 part. + max_value = max(values) + binsize = max_value/nbins + min_value = 0.5*binsize + while min_value > min(values): + min_value -= binsize + nbins += 1 + + # Create histogram + hist_x = [] + hist_y = [] + hist_expect = [] + norm = 0.5*len(values) + for i in range(nbins): + x_min = min_value+i*binsize + x_max = min_value+(i+1)*binsize + x_val = 0.5*(x_min+x_max) + y_val = 0.0 + for ts in values: + if ts > x_min: + y_val += 1.0 + hist_x.append(x_min) + hist_y.append(y_val) + hist_x.append(x_max) + hist_y.append(y_val) + if x_min > 0.0: + value = norm*erf(math.sqrt(0.5*x_min)) + norm + hist_expect.append(len(values)-value) + hist_expect.append(len(values)-value) + else: + hist_expect.append(norm) + hist_expect.append(norm) + + # Show expected distribution (half) + l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label="simulations") + l = plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label="expectation") + l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) + + # Set plot + plt.xlabel('TS') + plt.ylabel('Number of trials > TS') + plt.title(title) + plt.grid(True) + plt.legend(loc="upper right") + + # Show histogram + plt.show() + + # Return + return # ======================================== # # Display probability density distribution # # ======================================== # def dist_pdf(values, nbins, title): - """ - Display probability density TS distribution and overplots - expectation. - """ - # Set range and adapt number of bins. We make sure that we have a - # histogram bin centred on 0 that will capture the TS=0 part. - max_value = max(values) - binsize = max_value/nbins - min_value = 0.5*binsize - while min_value > min(values): - min_value -= binsize - nbins += 1 - - # Create histogram - n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value], - align='mid', facecolor='green', - log=True, label="simulations") - - # Create expected distribution (onyl for positive TS). We compute the - # full and half of the distribution, as for positively constrained - # quantities half of the values have TS=0 - x = [] - y_full = [] - y_half = [] - width = bins[1]-bins[0] - norm = len(values)/(math.sqrt(2.0)*math.sqrt(math.pi))*width - x.append(0.0) - y_full.append(len(values)) - y_half.append(0.5*len(values)) - for bin in bins: - ts = bin + 0.5*width - if ts > 0: - y = norm*math.pow(ts,-0.5)*math.exp(-0.5*ts) - x.append(ts) - y_full.append(y) - y_half.append(0.5*y) - - # Show expected distribution (half) - l = plt.semilogy(x, y_half, 'ro', linewidth=2, label="expectation") - - # Set plot - plt.xlabel('TS') - plt.ylabel('Number of trials') - plt.title(title) - plt.grid(True) - plt.legend(loc="upper right") - - # Show histogram - plt.show() - - + """ + Display probability density TS distribution and overplots + expectation. + """ + # Set range and adapt number of bins. We make sure that we have a + # histogram bin centred on 0 that will capture the TS=0 part. + max_value = max(values) + binsize = max_value/nbins + min_value = 0.5*binsize + while min_value > min(values): + min_value -= binsize + nbins += 1 + + # Create histogram + n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value], + align='mid', facecolor='green', + log=True, label="simulations") + + # Create expected distribution (onyl for positive TS). We compute the + # full and half of the distribution, as for positively constrained + # quantities half of the values have TS=0 + x = [] + y_full = [] + y_half = [] + width = bins[1]-bins[0] + norm = len(values)/(math.sqrt(2.0)*math.sqrt(math.pi))*width + x.append(0.0) + y_full.append(len(values)) + y_half.append(0.5*len(values)) + for bin in bins: + ts = bin + 0.5*width + if ts > 0: + y = norm*math.pow(ts,-0.5)*math.exp(-0.5*ts) + x.append(ts) + y_full.append(y) + y_half.append(0.5*y) + + # Show expected distribution (half) + l = plt.semilogy(x, y_half, 'ro', linewidth=2, label="expectation") + + # Set plot + plt.xlabel('TS') + plt.ylabel('Number of trials') + plt.title(title) + plt.grid(True) + plt.legend(loc="upper right") + + # Show histogram + plt.show() + + # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Display TS distribution generated using cstsdist using matplotlib. - """ - # Print usage information - usage = "Usage: show_ts_distribution filename" \ - " [-n bins] [-c column] [-t title] [-p plot]" - if len(sys.argv) < 2: - sys.stdout.write(usage+"\n") - sys.exit() - - # Extract parameters - filename = sys.argv[1] - - # Set default parameters - nbins = 30 - tsname = "TS" - title = "TS distribution" - plot = "pdf" - - # Parameter dictionnary - pars = [{'option': '-n', 'value': nbins}, - {'option': '-c', 'value': tsname}, - {'option': '-t', 'value': title}, - {'option': '-p', 'value': plot}] - - # Gather parameters from command line - i = 2 - while i < len(sys.argv): - - # Search for options - for par in pars: - if sys.argv[i] == par['option']: - if len(sys.argv) > i+1: - i += 1 - try: - par['value'] = sys.argv[i] - except: - sys.stdout.write(usage+"\n") - sys.exit() - else: - sys.stdout.write(usage+"\n") - sys.exit() - - # Next item - i += 1 - - # Recover parameters - nbins = int(pars[0]['value']) - tsname = pars[1]['value'] - title = pars[2]['value'] - plot = pars[3]['value'] - - # Read values from CSV file - values = read_ts(filename, tsname=tsname) - sys.stdout.write(len(values), "values read.\n") - - # Show histogram - if plot == "pdf": - dist_pdf(values, nbins, title) - else: - dist_cdf(values, nbins, title) + """ + Display TS distribution generated using cstsdist using matplotlib. + """ + # Print usage information + usage = "Usage: show_ts_distribution filename" \ + " [-n bins] [-c column] [-t title] [-p plot]" + if len(sys.argv) < 2: + sys.stdout.write(usage+"\n") + sys.exit() + + # Extract parameters + filename = sys.argv[1] + + # Set default parameters + nbins = 30 + tsname = "TS" + title = "TS distribution" + plot = "pdf" + + # Parameter dictionnary + pars = [{'option': '-n', 'value': nbins}, + {'option': '-c', 'value': tsname}, + {'option': '-t', 'value': title}, + {'option': '-p', 'value': plot}] + + # Gather parameters from command line + i = 2 + while i < len(sys.argv): + + # Search for options + for par in pars: + if sys.argv[i] == par['option']: + if len(sys.argv) > i+1: + i += 1 + try: + par['value'] = sys.argv[i] + except: + sys.stdout.write(usage+"\n") + sys.exit() + else: + sys.stdout.write(usage+"\n") + sys.exit() + + # Next item + i += 1 + + # Recover parameters + nbins = int(pars[0]['value']) + tsname = pars[1]['value'] + title = pars[2]['value'] + plot = pars[3]['value'] + + # Read values from CSV file + values = read_ts(filename, tsname=tsname) + sys.stdout.write(len(values), "values read.\n") + + # Show histogram + if plot == "pdf": + dist_pdf(values, nbins, title) + else: + dist_cdf(values, nbins, title) diff --git a/examples/simulate_events.py b/examples/simulate_events.py index 100a14a..8313e42 100755 --- a/examples/simulate_events.py +++ b/examples/simulate_events.py @@ -86,10 +86,10 @@ def createobs(ra=86.171648, dec=-1.4774586, rad=5.0, events.gti(gti) events.ebounds(ebounds) obs.events(events) - + # Set instrument response obs.response(irf, db) - + # Set ontime, livetime, and deadtime correction factor obs.ontime(duration) obs.livetime(duration*deadc) @@ -155,7 +155,7 @@ def simulate_ctobssim(obs, xmlname, seed=0): # Run simulator sim.run() - + # Retrieve events events = sim.obs()[0].events().copy() @@ -169,7 +169,7 @@ def simulate_ctobssim(obs, xmlname, seed=0): # Return events return events - + # =========== # # Show events # diff --git a/test/test_ctbin.py b/test/test_ctbin.py index b6c1161..3abdd03 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -77,7 +77,7 @@ class Test(gammalib.GPythonTestSuite): bin["proj"] = "CAR" bin["xref"] = 83.63 bin["yref"] = 22.01 - + # Run tool self.test_try("Run ctbin") try: diff --git a/test/test_ctbkgcube.py b/test/test_ctbkgcube.py index bb52469..df17278 100755 --- a/test/test_ctbkgcube.py +++ b/test/test_ctbkgcube.py @@ -85,7 +85,7 @@ class Test(gammalib.GPythonTestSuite): bkgcube["proj"] = "CAR" bkgcube["xref"] = 83.63 bkgcube["yref"] = 22.01 - + # Run tool self.test_try("Run ctbkgcube") try: diff --git a/test/test_ctbutterfly.py b/test/test_ctbutterfly.py index 4fde8f6..ce148ac 100755 --- a/test/test_ctbutterfly.py +++ b/test/test_ctbutterfly.py @@ -75,7 +75,7 @@ class Test(gammalib.GPythonTestSuite): butterfly["emin"] = 0.1 butterfly["emax"] = 100 butterfly["outfile"] = "butterfly.txt" - + # Run tool self.test_try("Run ctbutterfly") try: @@ -91,6 +91,6 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in saving results.") - + # Return return diff --git a/test/test_ctcubemask.py b/test/test_ctcubemask.py index 8d2179f..8b6f8e8 100755 --- a/test/test_ctcubemask.py +++ b/test/test_ctcubemask.py @@ -73,7 +73,7 @@ class Test(gammalib.GPythonTestSuite): mask["rad"] = 2.0 mask["emin"] = 0.1 mask["emax"] = 100 - + # Run tool self.test_try("Run ctcubemask") try: diff --git a/test/test_cterror.py b/test/test_cterror.py index a9aaaa6..2e90bef 100755 --- a/test/test_cterror.py +++ b/test/test_cterror.py @@ -74,7 +74,7 @@ class Test(gammalib.GPythonTestSuite): error["srcname"] = "Crab" error["caldb"] = self.caldb error["irf"] = self.irf - + # Run tool self.test_try("Run cterror") try: @@ -90,6 +90,6 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in saving results.") - + # Return return diff --git a/test/test_ctexpcube.py b/test/test_ctexpcube.py index df5cb23..2e98ff5 100755 --- a/test/test_ctexpcube.py +++ b/test/test_ctexpcube.py @@ -82,7 +82,7 @@ class Test(gammalib.GPythonTestSuite): expcube["proj"] = "CAR" expcube["xref"] = 83.63 expcube["yref"] = 22.01 - + # Run tool self.test_try("Run ctexpcube") try: diff --git a/test/test_ctlike.py b/test/test_ctlike.py index bc6b26b..2912bbc 100755 --- a/test/test_ctlike.py +++ b/test/test_ctlike.py @@ -88,6 +88,6 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in saving results.") - + # Return return diff --git a/test/test_ctmodel.py b/test/test_ctmodel.py index 8e7946c..22581c0 100755 --- a/test/test_ctmodel.py +++ b/test/test_ctmodel.py @@ -90,7 +90,7 @@ class Test(gammalib.GPythonTestSuite): model["proj"] = "CAR" model["xref"] = 83.63 model["yref"] = 22.01 - + # Run tool self.test_try("Run ctmodel") try: diff --git a/test/test_ctobssim.py b/test/test_ctobssim.py index a30917e..bbb91f3 100755 --- a/test/test_ctobssim.py +++ b/test/test_ctobssim.py @@ -78,7 +78,7 @@ class Test(gammalib.GPythonTestSuite): sim["tmax"] = 1800.0 sim["emin"] = 0.1 sim["emax"] = 100.0 - + # Run tool self.test_try("Run ctobssim") try: @@ -98,7 +98,7 @@ class Test(gammalib.GPythonTestSuite): self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "ROI Right Ascension is not 83.63 deg") self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "ROI Declination is not 22.01 deg") self.test_value(evt.size(), 4134, "Number of events is not 4134") - + # Save events self.test_try("Save events") try: @@ -114,7 +114,7 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up observation container obs = self.set_obs(4) - + # Set-up ctobssim sim = ctools.ctobssim(obs) sim["outevents"] = "sim_events.xml" @@ -185,18 +185,18 @@ class Test(gammalib.GPythonTestSuite): pnt = gammalib.GCTAPointing() pnt.dir(pntdir) obs_cta.pointing(pnt) - + # Set ROI roi = gammalib.GCTARoi() instdir = gammalib.GCTAInstDir() instdir.dir(pntdir) roi.centre(instdir) roi.radius(rad) - + # Set GTI gti = gammalib.GGti() gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) - + # Set energy boundaries ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ gammalib.GEnergy(emax, "TeV")) @@ -207,16 +207,15 @@ class Test(gammalib.GPythonTestSuite): events.gti(gti) events.ebounds(ebounds) obs_cta.events(events) - + # Set instrument response obs_cta.response(self.irf, db) - + # Set ontime, livetime, and deadtime correction factor obs_cta.ontime(duration) obs_cta.livetime(duration*deadc) obs_cta.deadc(deadc) obs_cta.id(id) - + # Return CTA observation return obs_cta - \ No newline at end of file diff --git a/test/test_ctpsfcube.py b/test/test_ctpsfcube.py index 07375aa..15bf77f 100755 --- a/test/test_ctpsfcube.py +++ b/test/test_ctpsfcube.py @@ -84,7 +84,7 @@ class Test(gammalib.GPythonTestSuite): psfcube["yref"] = 22.01 psfcube["amax"] = 0.3 psfcube["anumbins"] = 10 - + # Run tool self.test_try("Run ctpsfcube") try: diff --git a/test/test_ctselect.py b/test/test_ctselect.py index f1f57cc..b02237f 100755 --- a/test/test_ctselect.py +++ b/test/test_ctselect.py @@ -89,6 +89,6 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in saving events.") - + # Return return diff --git a/test/test_ctskymap.py b/test/test_ctskymap.py index 811270c..bc557e2 100755 --- a/test/test_ctskymap.py +++ b/test/test_ctskymap.py @@ -75,7 +75,7 @@ class Test(gammalib.GPythonTestSuite): skymap["proj"] = "CAR" skymap["xref"] = 83.63 skymap["yref"] = 22.01 - + # Run tool self.test_try("Run ctskymap") try: diff --git a/test/test_cttsmap.py b/test/test_cttsmap.py index a6382ec..e78439b 100755 --- a/test/test_cttsmap.py +++ b/test/test_cttsmap.py @@ -96,6 +96,6 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in saving results.") - + # Return return diff --git a/test/test_ctulimit.py b/test/test_ctulimit.py index 0684637..8deec9f 100755 --- a/test/test_ctulimit.py +++ b/test/test_ctulimit.py @@ -72,7 +72,7 @@ class Test(gammalib.GPythonTestSuite): ulimit["srcname"] = "Crab" ulimit["caldb"] = self.caldb ulimit["irf"] = self.irf - + # Run tool self.test_try("Run ctulimit") try: @@ -88,6 +88,6 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in saving results.") - + # Return return diff --git a/test/test_pipelines.py b/test/test_pipelines.py index 914c2c1..9c4e366 100755 --- a/test/test_pipelines.py +++ b/test/test_pipelines.py @@ -75,7 +75,7 @@ class Test(gammalib.GPythonTestSuite): emin = 0.1 emax = 100.0 rad_select = 3.0 - + # Simulate events sim = ctools.ctobssim() sim["inmodel"] = model_name @@ -95,7 +95,7 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in ctobssim.") - + # Select events select = ctools.ctselect() select["inobs"] = events_name @@ -113,7 +113,7 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in ctselect.") - + # Perform maximum likelihood fitting like = ctools.ctlike() like["inobs"] = selected_events_name @@ -146,7 +146,7 @@ class Test(gammalib.GPythonTestSuite): emin = 0.1 emax = 100.0 rad_select = 3.0 - + # Simulate events sim = ctools.ctobssim() sim["inmodel"] = model_name @@ -165,7 +165,7 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in ctobssim.") - + # Select events select = ctools.ctselect(sim.obs()) select["ra"] = ra @@ -181,7 +181,7 @@ class Test(gammalib.GPythonTestSuite): self.test_try_success() except: self.test_try_failure("Exception occured in ctselect.") - + # Perform maximum likelihood fitting like = ctools.ctlike(select.obs()) self.test_try("Run ctlike") diff --git a/test/test_python.py b/test/test_python.py index 911178b..4f6daf6 100755 --- a/test/test_python.py +++ b/test/test_python.py @@ -110,7 +110,7 @@ if __name__ == '__main__': # Run test suite success = suites.run() - + # Save test results suites.save("reports/ctools.xml")