Support #1891
ctlike runs for a long time
| Status: | Closed | Start date: | 01/03/2017 | |
|---|---|---|---|---|
| Priority: | Normal | Due date: | ||
| Assigned To: |  Knödlseder Jürgen | % Done: | 0% | |
| Category: | - | |||
| Target version: | - | |||
| Duration: |
Description
Hi,
I’m testing out my datafiles on csspec, but it seems to take forever with only ~1000 events, so I’m wondering if I’ve set something up wrong:
My GObservations object:
=== GObservations === Number of observations ....: 22 Number of models ..........: 23 Number of observed events .: 1178 Number of predicted events : 0
First GObservation object:
=== GCTAObservation === Name ......................: Identifier ................: VR65555_CH6 Instrument ................: VERITAS Event file ................: /nv/hp11/nkelleyh3/scratch/temp/temp_dm/pntsrc/2017.01.03./nv/hp11/nkelleyh3/data2/datadm_leviathan/pntsrc/ReconMethodGeom.Cut-NTel2-ExtendedSource-Hard.crab.pntsrc.json/fits_dir/VR65555.ReconMethodGeom.Cut-NTel2-ExtendedSource-Hard.chunk6.fits Event type ................: EventList Statistics ................: Poisson Ontime ....................: 120 s Livetime ..................: 109.285579281195 s Deadtime correction .......: 0.910713160676625 User energy range .........: undefined === GCTAPointing === Pointing direction ........: (RA,Dec)=(83.6340572925004,21.5144713376947) === GCTAResponseIrf === Caldb mission .............: Caldb instrument ..........: Response name .............: Energy dispersion .........: Not used Save energy range .........: undefined === GCTAEventList === Number of events ..........: 56 (disposed in "/nv/hp11/nkelleyh3/scratch/temp/temp_dm/pntsrc/2017.01.03./nv/hp11/nkelleyh3/data2/datadm_leviathan/pntsrc/ReconMethodGeom.Cut-NTel2-ExtendedSource-Hard.crab.pntsrc.json/fits_dir/VR65555.ReconMethodGeom.Cut-NTel2-ExtendedSource-Hard.chunk6.fits") Time interval .............: 56278.4661111111 - 56278.4674074074 days Energy interval ...........: 0.02 - 500 TeV Region of interest ........: RA=83.634057, DEC=21.514471 [0,0] Radius=2.25 deg
I’m noticing its GCTAResponseIrf is empty, but I’m not sure why.
One of the background models:
=== GCTAModelIrfBackground === Name ......................: plaw_VR65555_CH6 Instruments ...............: VERITAS Instrument scale factors ..: unity Observation identifiers ...: VR65555_CH6 Model type ................: "PowerLaw" * "Constant" Number of parameters ......: 4 Number of spectral par's ..: 3 Prefactor ................: 2.83e-17 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient) Index ....................: -2.7 +/- 0 [-10,-0.1] (free,scale=1,gradient) PivotEnergy ..............: 1000000 MeV (fixed,scale=1,gradient) Number of temporal par's ..: 1 Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
My point source:
=== GModelSky === Name ......................: custom_point_source Instruments ...............: all Test Statistic ............: Computation requested Instrument scale factors ..: unity Observation identifiers ...: all Model type ................: PointSource Model components ..........: "PointSource" * "PowerLaw" * "Constant" Number of parameters ......: 6 Number of spatial par's ...: 2 RA .......................: 83.6333315523972 deg (fixed,scale=1) DEC ......................: 22.0144708833273 deg (fixed,scale=1) Number of spectral par's ..: 3 Prefactor ................: 2.83e-17 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient) Index ....................: -2 [-10,10] (fixed,scale=1,gradient) PivotEnergy ..............: 1500000 MeV (fixed,scale=1,gradient) Number of temporal par's ..: 1 Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
A sample of events:
Dir=RA=83.9617156982422, DEC=20.9645671844482 [0.0109378429056319,0.00094791545451527] Energy=2.2598831653595 TeV Time=209301068.645663 s (TT) Dir=RA=83.198860168457, DEC=20.9807472229004 [0.00558848230144966,-0.0102796696899487] Energy=2.87266707420349 TeV Time=209301079.256931 s (TT) Dir=RA=82.6423492431641, DEC=21.2432231903076 [-0.00232960458481034,-0.0166376952749173] Energy=4.44022941589355 TeV Time=209301084.938934 s (TT) Dir=RA=84.2462387084961, DEC=20.1632556915283 [0.0256108191743704,-0.000495546041862671] Energy=2.62396693229675 TeV Time=209301108.688608 s (TT) Dir=RA=82.4786834716797, DEC=20.4782047271729 [0.00870472043391108,-0.0246182902450169] Energy=3.08900237083435 TeV Time=209301117.877033 s (TT)
Each event list is saved in a separate file along with its effective area, psf, and background, so I have a special function that loads it all into one GObservation, with the background model creation and linking done afterwards.
def load_chunk( fname ) :
"""Load fits file 'fname' and return a gammalib.GCTAObservation object
**Args:**
fname : filename of fits file
**Returns:**
gammalib.GCTAObservation() object
"""
# initial
check_if_file_exists( fname )
run = gammalib.GCTAObservation()
# convert unicode to regular string
if isinstance(fname, bytes) :
fname = fname.decode('utf-8')
run.load( fname )
# parse the run and chunk id from the filename,
# and set them in the run
m = re.search( 'VR(\d+).*chunk(\d+)', fname )
obsid = VeritasObsID()
obsid.runnumber = int( m.group(1) )
obsid.chunkid = int( m.group(2) )
run.id( obsid2str( obsid ) )
# setup irfs
cal = gammalib.GCaldb( os.path.dirname( os.path.abspath( fname ) ) )
rsp = gammalib.GCTAResponseIrf()
rsp.caldb( cal )
rsp.load_aeff( fname )
rsp.load_psf( fname )
#rsp.load_edisp( fname )
rsp.load_background( fname )
run.response( rsp )
return run
When I run csspec, I get this:
events in each energy bin: i: 0 e: 2.00- 3.16 n: 161 i: 1 e: 3.16- 5.00 n: 148 2017-01-03T16:05:15: +============+ 2017-01-03T16:05:15: | Parameters | 2017-01-03T16:05:15: +============+ 2017-01-03T16:05:15: inobs .....................: [not queried] 2017-01-03T16:05:15: inmodel ...................: [not queried] 2017-01-03T16:05:15: srcname ...................: custom_point_source 2017-01-03T16:05:15: expcube ...................: [not queried] 2017-01-03T16:05:15: psfcube ...................: [not queried] 2017-01-03T16:05:15: edispcube .................: [not queried] 2017-01-03T16:05:15: bkgcube ...................: [not queried] 2017-01-03T16:05:15: caldb .....................: [not queried] 2017-01-03T16:05:15: irf .......................: [not queried] 2017-01-03T16:05:15: edisp .....................: no 2017-01-03T16:05:15: outfile ...................: [not queried] 2017-01-03T16:05:15: emin ......................: 2 2017-01-03T16:05:15: emax ......................: 5 2017-01-03T16:05:15: enumbins ..................: 2 2017-01-03T16:05:15: ebinalg ...................: LOG 2017-01-03T16:05:15: ebinfile ..................: [not queried] 2017-01-03T16:05:15: calc_ts ...................: yes 2017-01-03T16:05:15: calc_ulim .................: yes 2017-01-03T16:05:15: fix_srcs ..................: yes 2017-01-03T16:05:15: fix_bkg ...................: no 2017-01-03T16:05:15: publish ...................: no 2017-01-03T16:05:15: chatter ...................: 4 2017-01-03T16:05:15: clobber ...................: yes 2017-01-03T16:05:15: debug .....................: yes 2017-01-03T16:05:15: mode ......................: ql 2017-01-03T16:05:15: logfile ...................: csspec.log 2017-01-03T16:05:15: 2017-01-03T16:05:15: +==================+ 2017-01-03T16:05:15: | Spectral binning | 2017-01-03T16:05:15: +==================+ 2017-01-03T16:05:15: Bin 1 .....................: 2 TeV - 3.16227766016838 TeV 2017-01-03T16:05:15: Bin 2 .....................: 3.16227766016838 TeV - 5 TeV 2017-01-03T16:05:15: 2017-01-03T16:05:15: +=============+ 2017-01-03T16:05:15: | Observation | 2017-01-03T16:05:15: +=============+ 2017-01-03T16:05:15: === GObservations === 2017-01-03T16:05:15: Number of observations ....: 22 2017-01-03T16:05:15: Number of models ..........: 23 2017-01-03T16:05:15: Number of observed events .: 1178 2017-01-03T16:05:15: Number of predicted events : 0 2017-01-03T16:05:15: 2017-01-03T16:05:15: +=========================+ 2017-01-03T16:05:15: | Adjust model parameters | 2017-01-03T16:05:15: +=========================+ 2017-01-03T16:05:15: === plaw_VR65555_CH6 === 2017-01-03T16:05:15: === plaw_VR65457_CH5 === 2017-01-03T16:05:15: === plaw_VR65457_CH8 === 2017-01-03T16:05:15: === plaw_VR65555_CH10 === 2017-01-03T16:05:15: === plaw_VR65555_CH1 === 2017-01-03T16:05:15: === plaw_VR65555_CH9 === 2017-01-03T16:05:15: === plaw_VR65555_CH8 === 2017-01-03T16:05:15: === plaw_VR65457_CH11 === 2017-01-03T16:05:15: === plaw_VR65457_CH3 === 2017-01-03T16:05:15: === plaw_VR65457_CH2 === 2017-01-03T16:05:15: === plaw_VR65555_CH3 === 2017-01-03T16:05:15: === plaw_VR65457_CH4 === 2017-01-03T16:05:15: === plaw_VR65457_CH7 === 2017-01-03T16:05:15: === plaw_VR65457_CH9 === 2017-01-03T16:05:15: === plaw_VR65555_CH2 === 2017-01-03T16:05:15: === plaw_VR65555_CH7 === 2017-01-03T16:05:15: === plaw_VR65457_CH10 === 2017-01-03T16:05:15: === plaw_VR65555_CH4 === 2017-01-03T16:05:15: === plaw_VR65555_CH5 === 2017-01-03T16:05:15: === plaw_VR65457_CH6 === 2017-01-03T16:05:15: === plaw_VR65555_CH11 === 2017-01-03T16:05:15: === plaw_VR65457_CH1 === 2017-01-03T16:05:15: === custom_point_source === 2017-01-03T16:05:15: Fixing "Prefactor" 2017-01-03T16:05:15: Freeing "Prefactor" 2017-01-03T16:05:15: 2017-01-03T16:05:15: +===================+ 2017-01-03T16:05:15: | Generate spectrum | 2017-01-03T16:05:15: +===================+ 2017-01-03T16:05:15: === GEbounds === 2017-01-03T16:05:15: Number of intervals .......: 2 2017-01-03T16:05:15: Energy range ..............: 2 TeV - 5 TeV 2017-01-03T16:05:15: +--------------+ 2017-01-03T16:05:15: | Energy bin 1 | 2017-01-03T16:05:15: +--------------+ 2017-01-03T16:05:15: === Selecting events === 2017-01-03T16:05:16: === Performing fit === 2017-01-03T16:08:28: === Computing upper limit === 2017-01-03T16:08:28: Upper limit calculation failed. 2017-01-03T16:08:28: Bin 1 .....................: 9.78507298965837e-11 +/- 1.4816555160430926e-11 erg/cm2/s 2017-01-03T16:08:28: +--------------+ 2017-01-03T16:08:28: | Energy bin 2 | 2017-01-03T16:08:28: +--------------+ 2017-01-03T16:08:28: === Selecting events === 2017-01-03T16:08:29: === Performing fit === 2017-01-03T16:11:46: === Computing upper limit === 2017-01-03T16:11:46: Upper limit calculation failed.
So it looks like for 110 events in 22 GObservations, it took
2017-01-03T16:05:16: === Performing fit === 2017-01-03T16:08:28: === Computing upper limit ===
~3min to fit one energy bin?
This is a shortened analysis, only using 1/100 of the available data, and with no edisp table, so I’m wondering if this is the normal amount of time?
Any suggestions would be appreciated.
Recurrence
No recurrence.
History
#1
Updated by Knödlseder Jürgen over 7 years ago
- Status changed from New to In Progress
- Assigned To set to Knödlseder Jürgen
The time it takes for a model fit depends sensitively on the model. I see that you have many model components in your example, so this may explain the issue. In particular if some of your components are not point sources they will take more time. Not sure why the response information is empty.
Maybe you can make a ctlike run on a single energy bin to see what actually ctlike is doing (I recognize that csspec has not the possibility to switch-on the logging of ctlike; this should probably be changed; you may also modify the csspec script for doing that).
#2
Updated by Mayer Michael over 7 years ago
Not sure why the response information is empty.
I believe the problem is that you specify a GCaldb object before the loading the IRFs:
cal = gammalib.GCaldb( os.path.dirname( os.path.abspath( fname ) ) ) rsp = gammalib.GCTAResponseIrf() rsp.caldb( cal ) rsp.load_aeff( fname )
You should simply omit this by removing the first and third line. I am not sure, but this could affect the printing of the information. Try using
csspec chatter=4 to get more details.
I guess your problem of the long fit comes from the fact that you have very very few events per observation but you intend to fit 2 background parameter and at least 1 source parameter on the handful of events. This could cause problems in the fit convergence leading to the fitter jumping back and forth. My suggestion would be to try to fix the background parameters and rerun the fit (hidden parameter in csspec: fix_bkg=yes).
If the fit convergence is the problem, the upper limit calculation will be affected, too. The UL will start its calculations based on the errors of the precedent fit. If you have very large errors (due to strongly degenerate parameters), the UL computation will need its maximum iterations (default 50). You can avoid the UL computation using the hidden parameter in csspec: calc_ulim=no).
#3
Updated by Kelley-Hoskins Nathan over 7 years ago
Hi again,
After removing the caldb lines and digging around some more, I think I was loading backgrounds full of zeroes. Fixing this seems to get me plausible values for things. Thanks for your help!
-Nathan
#4
Updated by Kelley-Hoskins Nathan over 7 years ago
- Status changed from In Progress to Resolved
#5
Updated by Knödlseder Jürgen about 7 years ago
- Status changed from Resolved to Closed