diff --git a/cscripts/cspull.par b/cscripts/cspull.par index e70ff62..7ed2397 100644 --- a/cscripts/cspull.par +++ b/cscripts/cspull.par @@ -18,7 +18,7 @@ #============================ inobs, f, h, NONE,,, "Input event list, counts cube, or observation definition XML file" inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file" -onsrc, s, a, ,,, "Name of On source (only for On/Off analysis)" +onsrc, s, a, NONE,,, "Name of On source (only for On/Off analysis)" onrad, r, a, 0.2,,, "Radius of On region (deg)" caldb, s, a, prod2,,, "Calibration database" irf, s, a, South_0.5h,,, "Instrument response function" diff --git a/cscripts/cspull.py b/cscripts/cspull.py index 1b8d943..e025295 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -65,7 +65,7 @@ class cspull(ctools.csobservation): enumbins = self['enumbins'].integer() # Query parameters for On/Off observation - if len(self['onsrc'].string()) > 0: + if gammalib.toupper(self['onsrc'].string()) != 'NONE': self['onrad'].real() # Query parameters for binned if requested @@ -112,7 +112,7 @@ class cspull(ctools.csobservation): # Extract information from observation if obs.classname() == 'GCTAOnOffObservation': emin = obs.on_spec().ebounds().emin().TeV() - emax = obs.on_spec().ebounds().emin().TeV() + emax = obs.on_spec().ebounds().emax().TeV() mode = 'On/Off' else: emin = obs.events().ebounds().emin().TeV() @@ -164,7 +164,7 @@ class cspull(ctools.csobservation): coordsys = 'CEL' # If we have a On source name then set On region radius - if len(onsrc) > 0: + if gammalib.toupper(onsrc) != 'NONE': onrad = self['onrad'].real() emin = self['emin'].real() emax = self['emax'].real() diff --git a/doc/source/users/reference_manual/cspull.rst b/doc/source/users/reference_manual/cspull.rst index 1fbcfe2..559ed0b 100644 --- a/doc/source/users/reference_manual/cspull.rst +++ b/doc/source/users/reference_manual/cspull.rst @@ -41,6 +41,12 @@ General parameters ``inmodel [file]`` Input model XML file. +``onsrc [string]`` + Name of On source (only for On/Off analysis; specify ``NONE`` for other analyses). + +``onrad [real]`` + Radius of On region (deg). + ``caldb [string]`` Calibration database. diff --git a/test/data/sciver/crab_plaw_onoff.xml b/test/data/sciver/crab_plaw_onoff.xml new file mode 100644 index 0000000..76c69c6 --- /dev/null +++ b/test/data/sciver/crab_plaw_onoff.xml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/test/data/sciver/obs_onoff.xml b/test/data/sciver/obs_onoff.xml new file mode 100644 index 0000000..9d1bb78 --- /dev/null +++ b/test/data/sciver/obs_onoff.xml @@ -0,0 +1,18 @@ + + + + + + + + + + + + + + + + + + diff --git a/test/science_verification.py b/test/science_verification.py index 5f61848..15fd381 100755 --- a/test/science_verification.py +++ b/test/science_verification.py @@ -39,8 +39,8 @@ import cscripts # ========================== # # Generate pull distribution # # ========================== # -def generate_pull_distribution(model, obs='NONE', trials=100, \ - caldb='prod2', irf='South_50h', \ +def generate_pull_distribution(model, obs='NONE', onsrc='NONE', onrad=0.2, \ + trials=100, caldb='prod2', irf='South_50h', \ deadc=0.98, edisp=False, \ ra=83.63, dec=22.01, rad=5.0, \ emin=0.1, emax=100.0, enumbins=0, \ @@ -57,6 +57,10 @@ def generate_pull_distribution(model, obs='NONE', trials=100, \ Model XML filename (without .xml extension) obs : str, optional Input observation definition XML filename + onsrc : str, optional + Name of On source for On/Off analysis + onrad : float, optional + Radius of On region trials : int, optional Number of trials caldb : str, optional @@ -108,6 +112,8 @@ def generate_pull_distribution(model, obs='NONE', trials=100, \ pull = cscripts.cspull() pull['inobs'] = obs pull['inmodel'] = inmodel + pull['onsrc'] = onsrc + pull['onrad'] = onrad pull['outfile'] = outfile pull['caldb'] = caldb pull['irf'] = irf @@ -164,7 +170,6 @@ def analyse_pull_distribution(filename): # Read rows first = True - #index = -1 samples = 0.0 for row in reader: @@ -237,6 +242,7 @@ class sciver(gammalib.GPythonTestSuite): self.append(self.spec_plaw, 'Test power law model') self.append(self.spec_plaw_edisp, 'Test power law model with energy dispersion') self.append(self.spec_plaw_stacked, 'Test power law model with stacked analysis') + self.append(self.spec_plaw_onoff, 'Test power law model with On/Off analysis') self.append(self.spec_plaw2, 'Test power law 2 model') self.append(self.spec_smoothbplaw, 'Test smoothly broken power law model') self.append(self.spec_eplaw, 'Test exponentially cut off power law model') @@ -263,7 +269,8 @@ class sciver(gammalib.GPythonTestSuite): return # Generate and analyse pull distributions - def pull(self, model, obs='NONE', trials=100, duration=1800.0, + def pull(self, model, obs='NONE', onsrc='NONE', + trials=100, duration=1800.0, ra=83.63, dec=22.01, rad=5.0, emin=0.1, emax=100.0, enumbins=0, edisp=False, debug=False): @@ -276,6 +283,8 @@ class sciver(gammalib.GPythonTestSuite): Model XML filename (without .xml extension) obs : str, optional Input observation definition XML filename + onsrc : str, optional + Name of On source for On/Off analysis trials : int, optional Number of trials duration : float, optional @@ -298,8 +307,8 @@ class sciver(gammalib.GPythonTestSuite): Enable debugging? """ # Generate pull distribution - outfile = generate_pull_distribution(model, obs=obs, trials=trials, - duration=duration, + outfile = generate_pull_distribution(model, obs=obs, onsrc=onsrc, + trials=trials, duration=duration, ra=ra, dec=dec, rad=rad, emin=emin, emax=emax, enumbins=enumbins, @@ -402,6 +411,20 @@ class sciver(gammalib.GPythonTestSuite): self.test('Pull_BackgroundModel_Index') return + # Test power law model with On/Off analysis + def spec_plaw_onoff(self): + """ + Test power law model with On/Off analysis + """ + self.pull('data/sciver/crab_plaw_onoff', + obs='data/sciver/obs_onoff.xml', onsrc='Crab', + emin=0.1, emax=100.0, enumbins=20) + self.test('Pull_Crab_Prefactor') + self.test('Pull_Crab_Index') + self.test('Pull_Background_Prefactor') + self.test('Pull_Background_Index') + return + # Test power law 2 model def spec_plaw2(self): """