Xspec Manual
User Manual:
Open the PDF directly: View PDF .
Page Count: 414
Download | ![]() |
Open PDF In Browser | View PDF |
An X-Ray Spectral Fitting Package Users’ Guide for version 12.9.1 Keith Arnaud, Craig Gordon & Ben Dorman HEASARC Astrophysics Science Division NASA/GSFC Greenbelt, MD 20771 Jan 2017 Contents 1 Introduction 1 1.1 New in v12.9.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 How to find out more information . . . . . . . . . . . . . . . . . 4 1.3 History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Spectral Fitting and XSPEC 7 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 The Basics of Spectral Fitting . . . . . . . . . . . . . . . . . . . . 7 2.3 The XSPEC implementation . . . . . . . . . . . . . . . . . . . . 9 2.4 Fits and Confidence Intervals . . . . . . . . . . . . . . . . . . . . 11 2.5 A more abstract and generalized approach . . . . . . . . . . . . . 12 2.6 XSPEC Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . 14 2.7 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3 XSPEC Overview and Helpful Hints i 17 ii CONTENTS 3.1 Syntax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 How to return to the XSPEC12>prompt . . . . . . . . . . . . . . 17 3.3 Getting Help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.4 Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5 Issuing Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.6 Control Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.7 Query, chatter and shutting XSPEC up (somewhat) . . . . . . . 20 3.8 Scripts and the Save command . . . . . . . . . . . . . . . . . . . 21 3.9 Miscellaneous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.10 Data Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.11 Model Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.12 Fitting Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.13 Binning and Grouping data . . . . . . . . . . . . . . . . . . . . . 26 3.14 Plotting Commands . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.15 Setting Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.16 Breaking With Ctrl-C . . . . . . . . . . . . . . . . . . . . . . . . 28 3.17 Customizing XSPEC . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.18 Customizing system-wide . . . . . . . . . . . . . . . . . . . . . . 30 4 Walks through XSPEC 31 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Brief Discussion of XSPEC Files . . . . . . . . . . . . . . . . . . 31 4.3 Fitting Models to Data: An Old Example from EXOSAT . . . . 32 4.4 Simultaneous Fitting . . . . . . . . . . . . . . . . . . . . . . . . . 58 CONTENTS iii 4.5 Multiple Models: a Background Modeling Example . . . . . . . . 62 4.6 Using XSPEC to Simulate Data: an Example for Chandra . . . . 65 4.7 Producing Plots: Modifying the Defaults . . . . . . . . . . . . . . 69 4.8 Markov Chain Monte Carlo Example . . . . . . . . . . . . . . . . 72 4.9 INTEGRAL/SPI: A Walk Through Example . . . . . . . . . . . 74 4.10 INTEGRAL Specific Command Line Scripts . . . . . . . . . . . . 81 5 XSPEC Commands 85 5.1 Summary of Commands . . . . . . . . . . . . . . . . . . . . . . . 85 5.2 Description of Syntax . . . . . . . . . . . . . . . . . . . . . . . . 89 5.3 Control Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.3.1 autosave . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.3.2 chatter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.3.3 exit, quit . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.3.4 help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.3.5 log . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.3.6 parallel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.7 query . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.3.8 save . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.3.9 script . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.3.10 show . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.3.11 syscall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.3.12 tclout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.3.13 tcloutr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 iv CONTENTS 5.3.14 time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3.15 undo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3.16 version . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.4 Data Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.4.1 arf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.4.2 backgrnd . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.4.3 corfile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.4.4 cornorm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.4.5 data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.4.6 diagrsp . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.4.7 fakeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.4.8 ignore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.4.9 notice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.4.10 response . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.5 Fit Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.5.1 bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.5.2 chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.5.3 error (and rerror) . . . . . . . . . . . . . . . . . . . . . . . 128 5.5.4 fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.5.5 freeze (and rfreeze) . . . . . . . . . . . . . . . . . . . . . . 131 5.5.6 ftest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.5.7 goodness 5.5.8 margin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 CONTENTS 5.5.9 v renorm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.5.10 steppar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 5.5.11 thaw (and rthaw) . . . . . . . . . . . . . . . . . . . . . . . 135 5.5.12 weight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.6 Model Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.6.1 addcomp . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.6.2 addline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 5.6.3 delcomp . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 5.6.4 dummyrsp . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 5.6.5 editmod . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.6.6 energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 5.6.7 eqwidth . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 5.6.8 flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 5.6.9 gain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 5.6.10 identify . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 5.6.11 initpackage . . . . . . . . . . . . . . . . . . . . . . . . . . 152 5.6.12 lmod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 5.6.13 lumin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 5.6.14 mdefine . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 5.6.15 model (and rmodel) . . . . . . . . . . . . . . . . . . . . . 158 5.6.16 modid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 5.6.17 newpar (and rnewpar) . . . . . . . . . . . . . . . . . . . . 163 5.6.18 systematic . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 vi CONTENTS 5.6.19 untie (and runtie) . . . . . . . . . . . . . . . . . . . . . . 168 5.7 5.8 5.9 Plot Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 5.7.1 cpd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 5.7.2 hardcopy . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 5.7.3 iplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 5.7.4 plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 5.7.5 setplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 Setting Commands . . . . . . . . . . . . . . . . . . . . . . . . . . 182 5.8.1 abund . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 5.8.2 cosmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 5.8.3 method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 5.8.4 statistic . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 5.8.5 xsect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 5.8.6 xset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 Tcl Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 5.9.1 lrt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 5.9.2 multifake . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 5.9.3 rescalecov . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 5.9.4 simftest . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 5.9.5 writefits . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 6 XSPEC Models 193 6.1 Alphabetical Summary of Models . . . . . . . . . . . . . . . . . . 193 6.2 Additive Model Components . . . . . . . . . . . . . . . . . . . . 197 CONTENTS vii 6.2.1 agauss, zagauss: gaussian line profile in wavelength space 197 6.2.2 apec, vapec, vvapec: APEC emission spectrum . . . . . . 198 6.2.3 atable: tabulated additive model . . . . . . . . . . . . . . 200 6.2.4 bapec, bvapec, bvvapec: velocity broadened APEC thermal plasma model . . . . . . . . . . . . . . . . . . . . . . 200 6.2.5 bbody, zbbody: blackbody 6.2.6 bbodyrad: blackbody spectrum, area normalized . . . . . 203 6.2.7 bexrav: reflected e-folded broken power law, neutral medium203 6.2.8 bexriv: reflected e-folded broken power law, ionized medium204 6.2.9 bknpower: broken power law . . . . . . . . . . . . . . . . 205 . . . . . . . . . . . . . . . . . 202 6.2.10 bkn2pow: broken power law, 2 break energies . . . . . . . 206 6.2.11 bmc: Comptonization by relativistic matter . . . . . . . . 206 6.2.12 bremss, vbremss, zbremss: thermal bremsstrahlung . . . . 207 6.2.13 btapec, bvtapec, bvvtapec: velocity broadened APEC emission spectrum with separate continuum and line temperatures . . . . . . . . . . . . . . . . . . . . . . . . . . . 208 6.2.14 c6mekl, c6vmekl, c6pmekl, c6pvmkl: differential emission measure using Chebyshev representations with multitemperature mekal . . . . . . . . . . . . . . . . . . . . . . 210 6.2.15 carbatm: Nonmagnetic carbon atmosphere of a neutron star . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 6.2.16 cemekl, cevmkl: plasma emission, multi-temperature using mekal . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 6.2.17 cflow: cooling flow . . . . . . . . . . . . . . . . . . . . . . 212 6.2.18 compbb: Comptonization, black body . . . . . . . . . . . 213 6.2.19 compLS: Comptonization, Lamb & Sanford . . . . . . . . 213 viii CONTENTS 6.2.20 compmag: Thermal and bulk Comptonization for cylindrical accretion onto the polar cap of a magnetized neutron star . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 6.2.21 compPS: Comptonization, Poutanen & Svenson . . . . . . 214 6.2.22 compST: Comptonization, Sunyaev & Titarchuk . . . . . 218 6.2.23 comptb: Thermal and bulk Comptonization of a seed blackbody-like spectrum . . . . . . . . . . . . . . . . . . . 218 6.2.24 compTT: Comptonization, Titarchuk . . . . . . . . . . . . 219 6.2.25 cplinear: a non-physical piecewise-linear model for low count background spectra. . . . . . . . . . . . . . . . . . . 219 6.2.26 cutoffpl: power law, high energy exponential cutoff . . . . 220 6.2.27 disk: accretion disk, black body . . . . . . . . . . . . . . . 221 6.2.28 diskbb: accretion disk, multi-black body components . . . 221 6.2.29 diskir: Irradiated inner and outer disk . . . . . . . . . . . 221 6.2.30 diskline: accretion disk line emission, relativistic . . . . . 222 6.2.31 diskm: accretion disk with gas pressure viscosity . . . . . 223 6.2.32 disko: accretion disk, inner, radiation pressure viscosity . 223 6.2.33 diskpbb: accretion disk, power-law dependence for T(r) . 223 6.2.34 diskpn: accretion disk, black hole, black body . . . . . . . 224 6.2.35 eplogpar: log-parabolic blazar model with νFν normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224 6.2.36 Eqpair, eqtherm, compth: Paolo Coppi’s hybrid (thermal/nonthermal) hot plasma emission models. . . . . . . . . . . . 225 6.2.37 equil, vequil: collisional plasma, ionization equilibrium . . 226 6.2.38 expdec: exponential decay . . . . . . . . . . . . . . . . . . 228 6.2.39 ezdiskbb: multiple blackbody disk model with zero-torque inner boundary . . . . . . . . . . . . . . . . . . . . . . . . 228 CONTENTS ix 6.2.40 gadem, vgadem: plasma emission, multi-temperature with gaussian distribution of emission measure. . . . . . . . . . 228 6.2.41 gauss, zgauss: gaussian line profile . . . . . . . . . . . . . 229 6.2.42 gnei, vgnei, vvgnei: collisional plasma, non-equilibrium, temperature evolution . . . . . . . . . . . . . . . . . . . . 230 6.2.43 grad: accretion disk, Schwarzschild black hole . . . . . . . 232 6.2.44 grbm: gamma-ray burst continuum . . . . . . . . . . . . . 232 6.2.45 hatm: Nonmagnetic hydrogen atmosphere of a neutron star233 6.2.46 kerrbb: multi-temperature blackbody model for thin accretion disk around a Kerr black hole . . . . . . . . . . . 234 6.2.47 kerrd: optically thick accretion disk around a Kerr black hole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235 6.2.48 kerrdisk: accretion disk line emission with BH spin as free parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . 235 6.2.49 laor: accretion disk, black hole emission line . . . . . . . . 236 6.2.50 laor2: accretion disk with broken-power law emissivity profile, black hole emission line . . . . . . . . . . . . . . . 236 6.2.51 logpar: log-parabolic blazar model . . . . . . . . . . . . . 236 6.2.52 lorentz: lorentz line profile . . . . . . . . . . . . . . . . . . 237 6.2.53 meka, vmeka: emission, hot diffuse gas (Mewe-Gronenschild)237 6.2.54 mekal, vmekal: emission, hot diffuse gas (Mewe-GronenschildKaastra) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 238 6.2.55 mkcflow, vmcflow: cooling flow, mekal . . . . . . . . . . . 239 6.2.56 nei, vnei, vvnei: collisional plasma, non-equilibrium, constant temperature . . . . . . . . . . . . . . . . . . . . . . 240 6.2.57 nlapec: continuum-only APEC emission spectrum . . . . 242 6.2.58 npshock, vnpshock, vvnpshock: shocked plasma, plane parallel, separate ion, electron temperatures . . . . . . . . 243 x CONTENTS 6.2.59 nsa: neutron star atmosphere . . . . . . . . . . . . . . . . 245 6.2.60 nsagrav: NS H atmosphere model for different g . . . . . 246 6.2.61 nsatmos: NS Hydrogen Atmosphere model with electron conduction and self-irradiation . . . . . . . . . . . . . . . 247 6.2.62 nsmax: Neutron Star Magnetic Atmosphere . . . . . . . . 247 6.2.63 nsmaxg: Neutron Star with a Magnetic Atmosphere . . . 249 6.2.64 nsx: neutron star with a non-magnetic atmosphere . . . . 251 6.2.65 nteea: non-thermal pair plasma . . . . . . . . . . . . . . . 251 6.2.66 Nthcomp: Thermally comptonized continuum . . . . . . . 252 6.2.67 Optxagnf, optxagn: Colour temperature corrected disc and energetically coupled Comptonisation model for AGN. 253 6.2.68 pegpwrlw: power law, pegged normalization . . . . . . . . 255 6.2.69 pexmon: neutral Compton reflection with self-consistent Fe and Ni lines. . . . . . . . . . . . . . . . . . . . . . . . . 256 6.2.70 pexrav: reflected powerlaw, neutral medium . . . . . . . . 257 6.2.71 pexriv: reflected powerlaw, ionized medium . . . . . . . . 257 6.2.72 plcabs: powerlaw observed through dense, cold matter . . 258 6.2.73 posm: positronium continuum . . . . . . . . . . . . . . . 259 6.2.74 powerlaw, zpowerlw: power law photon spectrum . . . . . 259 6.2.75 pshock, vpshock, vvpshock: plane-parallel shocked plasma, constant temperature . . . . . . . . . . . . . . . . . . . . 260 6.2.76 raymond, vraymond: emission, hot diffuse gas, RaymondSmith . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262 6.2.77 redge: emission, recombination edge . . . . . . . . . . . . 263 6.2.78 refsch: reflected power law from ionized accretion disk . . 263 CONTENTS xi 6.2.79 rnei, vrnei, vvrnei: non-equilibrium recombining collisional plasma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264 6.2.80 sedov, vsedov, vvsedov: sedov model, separate ion, electron temperatures . . . . . . . . . . . . . . . . . . . . . . 266 6.2.81 sirf: self-irradiated funnel . . . . . . . . . . . . . . . . . . 268 6.2.82 slimbh: Stationary slim accretion disk . . . . . . . . . . . 268 6.2.83 smaug: optically-thin, spherically-symmetric thermal plasma.269 6.2.84 snapec: galaxy cluster spectrum using SN yields . . . . . 271 6.2.85 srcut: synchrotron spectrum, cutoff power law . . . . . . 274 6.2.86 sresc: synchrotron spectrum, cut off by particle escape . . 275 6.2.87 step: step function convolved with gaussian . . . . . . . . 275 6.2.88 tapec, vtapec, vvtapec: APEC emission spectrum with separate continuum and line temperatures . . . . . . . . . 275 6.2.89 voigt: Voigt line profile . . . . . . . . . . . . . . . . . . . 278 6.3 Multiplicative Model Components . . . . . . . . . . . . . . . . . 278 6.3.1 absori: ionized absorber . . . . . . . . . . . . . . . . . . . 278 6.3.2 acisabs: Chandra ACIS q.e. decay . . . . . . . . . . . . . 279 6.3.3 cabs: Optically-thin Compton scattering. . . . . . . . . . 280 6.3.4 constant: energy-independent factor . . . . . . . . . . . . 280 6.3.5 cyclabs: absorption line, cyclotron . . . . . . . . . . . . . 280 6.3.6 dust: dust scattering . . . . . . . . . . . . . . . . . . . . . 281 6.3.7 edge, zedge: absorption edge . . . . . . . . . . . . . . . . 281 6.3.8 etable: exponential tabular model . . . . . . . . . . . . . 281 6.3.9 expabs: exponential roll-off at low E . . . . . . . . . . . . 282 6.3.10 expfac: exponential modification . . . . . . . . . . . . . . 282 xii CONTENTS 6.3.11 gabs: gaussian absorption line . . . . . . . . . . . . . . . . 283 6.3.12 heilin: Voigt absorption profiles for He I series . . . . . . 283 6.3.13 highecut, zhighect: high-energy cutoff . . . . . . . . . . . 283 6.3.14 hrefl: reflection model . . . . . . . . . . . . . . . . . . . . 284 6.3.15 ismabs: A high resolution ISM absorption model with variable columns for individual ions . . . . . . . . . . . . 285 6.3.16 lyman: Voigt absorption profiles for H I or He II Lyman series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 287 6.3.17 mtable: multiplicative tabular model . . . . . . . . . . . . 287 6.3.18 notch: absorption line, notch . . . . . . . . . . . . . . . . 287 6.3.19 pcfabs, zpcfabs: partial covering fraction absorption . . . 288 6.3.20 phabs, vphabs, zphabs, zvphabs: photoelectric absorption 288 6.3.21 plabs: power law absorption . . . . . . . . . . . . . . . . . 289 6.3.22 pwab: power-law distribution of neutral absorbers . . . . 290 6.3.23 recorn: change correction norm for a spectrum . . . . . . 290 6.3.24 redden: interstellar extinction . . . . . . . . . . . . . . . . 290 6.3.25 smedge: smeared edge . . . . . . . . . . . . . . . . . . . . 291 6.3.26 spexpcut: super-exponential cutoff absorption . . . . . . . 291 6.3.27 spline: spline modification . . . . . . . . . . . . . . . . . . 292 6.3.28 SSSice: Einstein SSS ice absorption . . . . . . . . . . . . 292 6.3.29 swind1: absorption by partially ionized material with large velocity shear . . . . . . . . . . . . . . . . . . . . . . . . . 292 6.3.30 tbabs, ztbabs, tbfeo, tbgas, tbgrain, tbpcf, tbvarabs, tbrel: ISM grain absorption . . . . . . . . . . . . . . . . . . . . 293 6.3.31 uvred: interstellar extinction, Seaton Law . . . . . . . . . 295 CONTENTS xiii 6.3.32 varabs, zvarabs: photoelectric absorption . . . . . . . . . 295 6.3.33 wabs, zwabs: photoelectric absorption, Wisconsin crosssections . . . . . . . . . . . . . . . . . . . . . . . . . . . . 296 6.3.34 wndabs, zwndabs: photo-electric absorption, warm absorber296 6.3.35 xion: reflected spectrum of photo-ionized accretion disk/ring297 6.3.36 xscat: dust scattering . . . . . . . . . . . . . . . . . . . . 299 6.3.37 zbabs: EUV ISM attenuation . . . . . . . . . . . . . . . . 299 6.3.38 zdust: extinction by dust grains . . . . . . . . . . . . . . . 300 6.3.39 zigm: UV/Optical attenuation by the intergalactic medium300 6.3.40 zredden: redshifted version of redden . . . . . . . . . . . . 301 6.3.41 zsmdust: extinction by dust grains in starburst galaxies . 301 6.3.42 zvfeabs: photoelectric absorption with free Fe edge energy 301 6.3.43 zxipcf: partial covering absorption by partially ionized material . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302 6.4 Convolution Model Components . . . . . . . . . . . . . . . . . . 302 6.4.1 cflux: calculate flux . . . . . . . . . . . . . . . . . . . . . 302 6.4.2 clumin: calculate luminosity . . . . . . . . . . . . . . . . . 303 6.4.3 cpflux: calculate photon flux . . . . . . . . . . . . . . . . 304 6.4.4 gsmooth: gaussian smoothing . . . . . . . . . . . . . . . . 305 6.4.5 ireflect: reflection from ionized material . . . . . . . . . . 305 6.4.6 kdblur: convolve with the laor model shape . . . . . . . . 307 6.4.7 kdblur2: convolve with the laor2 model shape . . . . . . . 307 6.4.8 kerrconv: accretion disk line shape with BH spin as free parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . 307 6.4.9 lsmooth: lorentzian smoothing . . . . . . . . . . . . . . . 308 xiv CONTENTS 6.4.10 partcov: partial covering . . . . . . . . . . . . . . . . . . . 308 6.4.11 rdblur: convolve with the diskline model shape . . . . . . 309 6.4.12 reflect: reflection from neutral material . . . . . . . . . . 309 6.4.13 rfxconv: angle-dependent reflection from an ionized disk . 310 6.4.14 simpl: comptonization of a seed spectrum . . . . . . . . . 311 6.4.15 xilconv: angle-dependent reflection from an ionized disk . 312 6.4.16 vashift: velocity shift an additive model . . . . . . . . . . 313 6.4.17 vmshift: velocity shift a multiplicative model . . . . . . . 314 6.4.18 zashift: redshift an additive model . . . . . . . . . . . . . 314 6.4.19 zmshift: redshift a multiplicative model . . . . . . . . . . 315 6.5 Pile-Up Model Components . . . . . . . . . . . . . . . . . . . . . 315 6.5.1 6.6 pileup: CCD pile-up model for Chandra . . . . . . . . . . 315 Mixing Model Components . . . . . . . . . . . . . . . . . . . . . 316 6.6.1 ascac: ASCA surface brightness model . . . . . . . . . . . 317 6.6.2 projct: project 3-D ellipsoidal shells onto 2-D elliptical annuli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 317 6.6.3 suzpsf: suzaku surface brightness model . . . . . . . . . . 319 6.6.4 xmmpsf: xmm surface brightness model . . . . . . . . . . 320 A The User Interface 323 A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323 A.2 XSPEC and tcl/tk . . . . . . . . . . . . . . . . . . . . . . . . . . 323 A.3 A note on command processing . . . . . . . . . . . . . . . . . . . 324 A.4 Command Recall/Editing . . . . . . . . . . . . . . . . . . . . . . 324 CONTENTS xv A.5 Logging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325 A.6 Command Completion . . . . . . . . . . . . . . . . . . . . . . . . 326 A.7 Unknown Procedure . . . . . . . . . . . . . . . . . . . . . . . . . 327 A.8 Aliases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 327 A.9 Initialization Script . . . . . . . . . . . . . . . . . . . . . . . . . . 328 A.10 XSPEC Command Result . . . . . . . . . . . . . . . . . . . . . . 328 A.11 Script Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 328 A.12 Unix Shell Commands . . . . . . . . . . . . . . . . . . . . . . . . 330 A.13 Writing Custom XSPEC commands . . . . . . . . . . . . . . . . 331 A.14 Scripting commands that prompt the user . . . . . . . . . . . . . 333 A.15 Script Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334 B Statistics in XSPEC 337 B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 337 B.2 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . 338 B.3 Parameter confidence regions . . . . . . . . . . . . . . . . . . . . 344 B.4 Goodness-of-fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 346 C Adding Models to XSPEC 351 C.1 Table models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 351 C.2 Loading a new model function . . . . . . . . . . . . . . . . . . . . 352 C.3 Writing a new model function . . . . . . . . . . . . . . . . . . . . 353 C.4 Third-Party Libraries In Local Models Build . . . . . . . . . . . 358 C.5 Writing new mixing models . . . . . . . . . . . . . . . . . . . . . 359 xvi CONTENTS D Overview of PLT 361 D.1 Getting started with PLT . . . . . . . . . . . . . . . . . . . . . . 362 D.2 PLT Command summary . . . . . . . . . . . . . . . . . . . . . . 362 E Associated Programs 365 E.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365 E.2 HEAsoft reading tasks . . . . . . . . . . . . . . . . . . . . . . . . 365 E.3 HEAsoft manipulation tasks . . . . . . . . . . . . . . . . . . . . . 366 E.4 HEAsoft subroutines . . . . . . . . . . . . . . . . . . . . . . . . . 366 F Using the XSPEC Models Library in Other Programs 367 F.1 Calling Model Functions From C And Fortran . . . . . . . . . . . 367 F.2 Interface Routines . . . . . . . . . . . . . . . . . . . . . . . . . . 368 F.3 Initializing the Models Library . . . . . . . . . . . . . . . . . . . 370 F.4 Building with the Models Library . . . . . . . . . . . . . . . . . . 370 G Adding a Custom Chain Proposal Algorithm 371 G.1 The randomize.dat Initialization File . . . . . . . . . . . . . . . 372 G.2 Writing a Chain Proposal Class . . . . . . . . . . . . . . . . . . . 373 G.3 Building and Loading the Proposal Class Library . . . . . . . . . 375 H Changes between v11 and v12 377 H.1 Integral Spectrometer/Coded Mask Instrument Support . . . . . 379 H.2 Current Exclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 379 I Older Release Notes 381 CONTENTS xvii I.1 v12.9.0 Jul 2015 . . . . . . . . . . . . . . . . . . . . . . . . . . . 381 I.2 V12.8.2 Jul 2014 . . . . . . . . . . . . . . . . . . . . . . . . . . . 384 I.3 V12.8.1 Aug 2013 . . . . . . . . . . . . . . . . . . . . . . . . . . . 385 I.4 v12.8.0 Dec 2012 . . . . . . . . . . . . . . . . . . . . . . . . . . . 387 I.5 v12.7.1 March 2012 . . . . . . . . . . . . . . . . . . . . . . . . . . 389 I.6 v12.7.0 May 2011 . . . . . . . . . . . . . . . . . . . . . . . . . . . 391 I.7 v12.6.0 March 2010 . . . . . . . . . . . . . . . . . . . . . . . . . . 393 I.8 v12.5.1 Aug. 2009 . . . . . . . . . . . . . . . . . . . . . . . . . . 395 Chapter 1 Introduction XSPEC is a command-driven, interactive, X-ray spectral-fitting program, designed to be completely detector-independent so that it can be used for any spectrometer. XSPEC has been used to analyze data from HEAO-1 A2, Einstein Observatory, EXOSAT, Ginga, ROSAT, BBXRT, ASCA, CGRO, IUE, RXTE, Chandra, XMM-Newton, Integral/SPI, Swift and Suzaku. There now almost 8000 papers listed on ADS which mention XSPEC. This manual describes XSPEC v12, which is available on Linux (gcc 3.2.2 and later), MacOSX/Darwin (gcc 3.3 and later) , Solaris (2.6-9; WS6.0 and later). New users are advised to read Chapter 2, which introduces spectral fitting and the XSPEC approach, Chapter 3, which gives an overview of the program commands, and Chapter 4, which contains walkthroughs of XSPEC sessions. They should then experiment with XSPEC and, if necessary, look up individual commands in Chapter 5 , or descriptions of the spectral models in use, in Chapter 6. The user interface for XSPEC, which employs a tcl scripting shell and the XSPEC parser are described in Appendix A. The statistical methods used within XSPEC are described in Appendix B. Users interested in adding their own models can read how to do so in Appendix C. Appendix D describes the PLT plotting package which XSPEC currently uses for graphical output. Some of the tools (FTOOLS, programs, scripts) that can operate on XSPEC files are listed in Appendix E. The XSPEC model library can be linked into other programs following the instructions in Appendix F. Appendix G describes how to use your own proposal distribution for Markov Chain Monte Carlo. Appendix H includes some notes on the changes between XSPEC v11 and v12. Finally, Appendix I lists release notes for older versions of v12. 1 2 CHAPTER 1. INTRODUCTION 1.1 New in v12.9.1 New features • New Models: carbatm clumin hatm ismabs rfxconv slimbh snapec (b)(v)(v)tapec tbfeo tbgas tbpcf tbrel vashift vmshift voigt xilconv xscat NS non-magnetic C atmosphere calculates luminosity NS non-magnetic H atmosphere High-resolution ISM absorption reflection from ionized disk slim accretion disk galaxy cluster spectrum using SN yields apec model with different temperatures for continuum and emission lines tbabs model with O and Fe only variable tbabs model with no grains tbabs model with partial covering tbabs model that allows negative columns velocity shift an additive component velocity shift a multiplicative component emission line with a Voigt profile reflection from ionized disk dust scattering out of the beam • Updated Models: – apec models updated to AtomDB 3.0.7 – nei models updated to AtomDB 3.0.4 eigen files – The apec code now has the ability to add new line shapes for broadening although at present only a single Gaussian line shape is operational. – A new xset option APECMINFLUX can be used to specify a minimum line flux below which lines are not broadened. Also the maximum flux in any line in the apec model is saved as apecmaxlineflux and this can be recovered using tcloutr modkeyval apecmaxlineflux. At the moment both this and APECMINFLUX do not include the 1e14 normalization factor and any time dilation correction. – The code for the various CEI models was rationalized and there is a new C++ routine calcMultiTempPlasma which can be used by other models. The old Fortran sumdem routine remains as a wrapper for old models. – tbabs model code has been updated. The old version is available by using xset TBABSVERSION 1. 1.1. NEW IN V12.9.1 3 – Updated nsmaxg model with new tabulated model files. – Increased the speed for accessing values from table models. – Modified table model code to support new keywords LOELIMIT and HIELIMIT to specify model values to be taken below and above the tabulated energies. • Plotting: – Added low and high energy limits to setplot id so lines are only shown for the specified energy range. – Added setplot (no)contimage to turn on (off) the background image in contour plots. – The keV to Angstrom conversion value stored in Numerics.h has been updated to that from CODATA 2014. • Fitting: – Added parallelization to chain walker calculations and to the goodness simulations. – The error command can now process parameters from more than 1 model in a single call. – In the various cstat statistics any correction file spectrum is now taken account. This does implicitly assume that the correction file has Poisson statistics. – The Tcl script lrt.tcl now uses the fakeit nowrite option so intermediate faked files are not written out. Also, now works for multiple data groups. • Input/Output: – Increased the precision of the reported output from the tclout command options. – Increased the allowed flux (counts) to beyond the 4-byte limit when running fakeit with counting statistics selected. – New tclout option fileinfo writes the value of the keyword given as the fileinfo argument. e.g. tclout fileinfo detnam returns the value of the detnam keyword in the spectrum extension. • Xspec.init file: – Added PARALLEL to set default parallelization options – Added NEI VERSION to set default version for NEI models – Added CONTOUR IMAGE to set default for adding image to contour plots – Changed default cross-sections from bcmc to vern 4 CHAPTER 1. INTRODUCTION • The manual has now been rewritten in LaTeX. The appendix on statistics has been expanded. Cross-referencing has been improved. Bug fixes All bug fixes to v12.9.0 released as patches are included in v12.9.1. In addition the following problems have been corrected: • A previous revision to the goodness command removed the requirement that the fit must be in a valid state at the start, however it was incorrectly marking the fit as being in a valid state at the end. • When re-reading a previously loaded table model file, XSPEC was not correctly recognizing whether or not parameter error values existed. • A spectrum’s displayed net count rate variance needs to be updated when the fit weighting method gets changed. 1.2 How to find out more information XSPEC is distributed and maintained under the aegis of the GSFC High Energy Astrophysics Science Archival Research Center (HEASARC). It can be downloaded as part of HEAsoft from http://heasarc.gsfc.nasa.gov/docs/software/lheasoft/download.html XSPEC is available either as binaries or source. We recommend downloading the source and compiling locally to avoid potential system library conflicts and allow installation of any patches we release. If you have any problems compiling or running XSPEC please e-mail us at xspec12@athena.gsfc.nasa.gov The XSPEC Web page comprises links to the anonymous ftp site, a Web version of the manual, and several useful documents, including a list of known bugs. The Web address is http://xspec.gsfc.nasa.gov/ with the list of issues and available patches at http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/bugs.html and additional models which can be added at http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/newmodels.html Further useful information can be found on the XSPEC Wiki at https://astrophysics.gsfc.nasa.gov/XSPECwiki/XSPECPage and via the XSPEC Facebook group. 1.3. HISTORY 1.3 5 History The first version of XSPEC was written in 1983 at the Institute of Astronomy, Cambridge, under VAX/VMS by Rick Shafer. It was written to perform spectral analysis of data from the ESA EXOSAT X-ray observatory, which was launched that year. XSPEC is a descendant of a series of spectral-fitting programs written at NASA/Goddard Space Flight Center for the OSO-8, HEAO-1 and EO missions. The initial design was the fruit of many discussions between Rick Shafer and Andy Szymkowiak. Rick Shafer subsequently joined the EXOSAT group, where he enhanced the VAX/VMS version in collaboration with Frank Haberl. Meanwhile, Keith Arnaud ported XSPEC to a Sun/UNIX operating system. The two implementations of XSPEC diverged for several years until a determined and coordinated effort was made to recover a single version. The results of that effort was XSPECv6, described in the first edition of this manual. Subsequent editions have covered later versions of the program. An extensive re-engineering effort was started in the fall of 1998 to improve longterm maintainability, allow significant new features to be added, and support the analysis of spectra from coded-mask instruments. 1.4 Acknowledgements This project would not have been possible without ignoring the advice of far more people than can be mentioned here. We would like to thank all our colleagues for their suggestions, bug reports, and (especially) source code. The GSFC X-ray astronomy group are particularly thanked for their patience exhibited while functioning as the beta test site for new releases. The initial development of XSPEC was funded by a Royal Society grant to Prof. Andy Fabian and subsequent development was partially funded by the European Space Agency’s EXOSAT project and is now funded by the HEASARC at NASA/GSFC. 1.5 References • Arnaud, K.A., 1996, Astronomical Data Analysis Software and Systems V, eds. G. Jacoby and J. Barnes,p17, ASP Conf. Series volume 101. • Dorman, B., and Arnaud, K. A. 2001, Astronomical Data Analysis Software and Systems X eds. F.R. Harnden, Jr., F.A. Primini, and H. E. Payne, vol. 238, p. 415 6 CHAPTER 1. INTRODUCTION • Dorman, B., Arnaud, K. A., and Gordon, C. A. XSPEC12: ObjectOriented X-Ray Data Analysis, AAS HEAD meeting No. 35, #22.10 Chapter 2 Spectral Fitting and XSPEC 2.1 Introduction This chapter comprises a brief description of the basics of spectral fitting, their application in XSPEC, and some helpful hints on how to approach particular problems. We then provide more details on the way XSPEC provides flexibility in its approach to the minimization problem. We also describe the data formats accepted. 2.2 The Basics of Spectral Fitting Although we use a spectrometer to measure the spectrum of a source, what the spectrometer obtains is not the actual spectrum, but rather photon counts (C) within specific instrument channels, (I). This observed spectrum is related to the actual spectrum of the source (f (E)) by: Z C(I) = f (E)R(I, E)dE (2.1) Where R(I, E) is the instrumental response and is proportional to the probability that an incoming photon of energy E will be detected in channel I. Ideally, then, we would like to determine the actual spectrum of a source, f (E), by 7 8 CHAPTER 2. SPECTRAL FITTING AND XSPEC inverting this equation, thus deriving f (E) for a given set of C(I). Regrettably, this is not possible in general, as such inversions tend to be non-unique and unstable to small changes in C(I). (For examples of attempts to circumvent these problems see Blissett & Cruise 1979; Kahn & Blissett 1980; Loredo & Epstein 1989). The usual alternative is to choose a model spectrum, f (E), that can be described in terms of a few parameters (i.e., f (E, p1, p2, ...)), and match, or “fit” it to the data obtained by the spectrometer. For each f (E), a predicted count spectrum (Cp (I)) is calculated and compared to the observed data (C(I)). Then a “fit statistic” is computed from the comparison and used to judge whether the model spectrum “fits” the data obtained by the spectrometer. The model parameters then are varied to find the parameter values that give the most desirable fit statistic. These values are referred to as the best-fit parameters. The model spectrum, fb (E), made up of the best-fit parameters is considered to be the best-fit model. The most common fit statistic in use for determining the “best-fit” model is χ2 , defined as follows: χ2 = X (C(I) − Cp (I))2 (σ(I))2 (2.2) where σ(I) is the (generally unknown) error p for channel I (e.g., if C(I) are counts then σ(I) is usually estimated by C(I); see e.g. Wheaton et.al. 1995 for other possibilities). Once a “best-fit” model is obtained, one must ask two questions: 1. How confident can one be that the observed C(I) can have been produced by the best-fit model fb (E)? The answer to this question is known as the “goodness-of-fit” of the model. The χ2 statistic provides a well-knowngoodness-of-fit criterion for a given number of degrees of freedom (ν, which is calculated as the number of channels minus the number of model parameters) and for a given confidence level. If χ2 exceeds a critical value (tabulated in many statistics texts) one can conclude that fb (E) is not an adequate model for C(I). As a general rule, one wants the “reduced χ2 ” (∼ χ2 /ν) to be approximately equal to one (i.e. χ2 ∼ ν). A reduced χ2 that is much greater than one indicates a poor fit, while a reduced χ2 that is much less than one indicates that the errors on the data have been over-estimated. Even if the best-fit model (fb (E)) does pass the “goodness-of-fit” test, one still cannot say that fb (E) is the only acceptable model. For example, if the data used in the fit are not particularly 2.3. THE XSPEC IMPLEMENTATION 9 good, one may be able to find many different models for which adequate fits can be found. In such a case, the choice of the correct model to fit is a matter of scientific judgment. 2. For a given best-fit parameter (p1), what is the range of values within which one can be confident the true value of the parameter lies? The answer to this question is the “confidence interval” for the parameter. The confidence interval for a given parameter is computed by varying the parameter value until the χ2 increases by a particular amount above the minimum, or “best-fit” value. The amount that the χ2 is allowed to increase (also referred to as the critical ∆χ2 ) depends on the confidence level one requires, and on the number of parameters whose confidence space is being calculated. The critical for common cases are given in the following table (from Avni, 1976): Confidence 0.68 0.90 0.99 2.3 Parameters 1 2 3 1.00 2.30 3.50 2.71 4.61 6.25 6.63 9.21 11.30 The XSPEC implementation To summarize the preceding section, the main components of spectral fitting are as follows: • A set of one or more observed spectra D(I) with background measurements B(I) where available • The corresponding instrumental responses R(I, E) • A set of model spectra M (E) These components are used in the following manner: • Choose a parameterized model which is thought to represent the actual spectrum of the source. • Choose values for the model parameters. • Based on the parameter values given, predict the count spectrum that would be detected by the spectrometer in a given channel for such a model. • Compare the predicted spectrum to the spectrum actually obtained by the instrument. 10 CHAPTER 2. SPECTRAL FITTING AND XSPEC • Manipulate the values of the parameters of the model until the best fit between the theoretical model and the observed data is found. • Then calculate the “goodness” of the fit to determine how well the model explains the observed data, and calculate the confidence intervals for the model’s parameters. This section describes how XSPEC performs these tasks. C(I): The Observed Spectrum To obtain each observed spectrum, C(I), XSPEC uses two files: the data (spectrum) file, containing D(I), and the background file, containing B(I). The data file tells XSPEC how many total photon counts were detected by the instrument in a given channel. XSPEC then uses the background file to derive the set of background-subtracted spectra C(I) in units of counts per second. The background-subtracted count rate is given by, for each spectrum: C(I) = bD(I) B(I) D(I) − aD(I) tD bB(I) ab(I) tB (2.3) where D(I) and B(I) are the counts in the data and background files; tD and tB are the exposure times in the data and background files; bD(I) and bB(I) , aD(I) and aB(I) are the background and area scaling values from the spectrum and background respectively, which together refer the background flux to the same area as the observation as necessary. When this is done, XSPEC has an observed spectrum to which the model spectrum can be fit. R(I,E): The Instrumental Response Before XSPEC can take a set of parameter values and predict the spectrum that would be detected by a given instrument, XSPEC must know the specific characteristics of the instrument. This information is known as the detector response. Recall that for each spectrum the response R(I, E) is proportional to the probability that an incoming photon of energy E will be detected in channel I. As such, the response is a continuous function of E. This continuous function is converted to a discrete function by the creator of a response matrix who defines the energy ranges Ej such that: R EJ RD (I, J) = EJ−1 R(I, E)dE EJ − EJ−1 (2.4) XSPEC reads both the energy ranges, EJ , and the response matrix RD (I, J) from a response file in a compressed format that only stores non-zero elements. 2.4. FITS AND CONFIDENCE INTERVALS 11 XSPEC also includes an option to use an auxiliary response file, which contains an array AD (J) that is multiplied into RD (I, J) as follows: RD (I, J) → RD (I, J) × AD (J) (2.5) This array is designed to represent the efficiency of the detector with the response file representing a normalized Redistribution Matrix Function, or RMF. Conventionally, the response is in units of cm2. M(E): The Model Spectrum The model spectrum, M (E), is calculated within XSPEC using the energy ranges defined by the response file : Z EJ MD (J) = M (E)dE (2.6) EJ−1 and is in units of photons/cm2 /s. XSPEC allows the construction of composite models consisting of additive components representing X-ray sources (e.g., power-laws, blackbodies, and so forth), multiplicative components, which modify additive components by an energy-dependent factor (e.g., photoelectric absorption, edges, ...). Convolution and mixing models can then perform sophisticated operations on the result. Models are defined in algebraic notation. For example, the following expression: phabs (power + phabs (bbody)) defines an absorbed blackbody, phabs(bbody), added to a power-law, power. The result then is modified by another absorption component, phabs. For a list of available models, see Chapter 6. 2.4 Fits and Confidence Intervals Once data have been read in and a model defined, XSPEC uses a fitting algorithm to find the best-fit values of the model parameter. The default is a modified Levenberg-Marquardt algorithm (based on CURFIT from Bevington, 1969). The algorithm used is local rather than global, so be aware that it is possible for the fitting process to get stuck in a local minimum and not find the 12 CHAPTER 2. SPECTRAL FITTING AND XSPEC global best-fit. The process also goes much faster (and is more likely to find the true minimum) if the initial model parameters are set to sensible values. The Levenberg-Marquardt algorithm relies on XSPEC calculating the 2nd derivatives of the fit statistic with respect to the model parameters. By default these are calculated analytically, with the assumption that the 2nd derivatives of the model itself may be ignored. This can be changed by setting the USE NUMERICAL DIFFERENTIATION flag to “true” in the Xspec.init initialization file, in which case XSPEC will perform numerical calculations of the derivatives (which are slower). At the end of a fit, XSPEC will write out the best-fit parameter values, along with estimated confidence intervals. These confidence intervals are one sigma and are calculated from the second derivatives of the fit statistic with respect to the model parameters at the best-fit. These confidence intervals are not reliable and should be used for indicative purposes only. XSPEC has a separate command (error or uncertain) to derive confidence intervals for one interesting parameter, which it does by fixing the parameter of interest at a particular value and fitting for all the other parameters. New values of the parameter of interest are chosen until the appropriate delta-statistic value is obtained. XSPEC uses a bracketing algorithm followed by an iterative cubic interpolation to find the parameter value at each end of the confidence interval. To compute confidence regions for several parameters at a time, XSPEC can run a grid on these parameters. XSPEC also will display a contour plot of the confidence regions of any two parameters. 2.5 A more abstract and generalized approach The sections above provide a simple characterization of the problem. XSPEC actually operates at a more abstract level and considers the following: Given a set of spectra C(I), each supplied as a function of detector channels, a set of theoretical models {M (E)j } each expressed in terms of a vector of energies together with a set of functions {R(I, E)j } that map channels to energies, minimize an objective function S of C, {R(I, E)j }, {M (E)j } using a fitting algorithm, i.e. X S = S(C, Mj × Rj ) (2.7) j In the default case, this reduces to the specific expression for χ2 fitting of a single 2.5. A MORE ABSTRACT AND GENERALIZED APPROACH 13 source S = χ2 = X (Ci − Ri × Mi )2 (2.8) i where i runs over all of the channels in all of the spectra being fitted, and using the Levenberg-Marquadt algorithm to perform the fitting. This differs from the previous formulation in that the operations that control the fitting process make fewer assumptions about how the data are formatted, what function is being minimized, and which algorithm is being employed. At the calculation level, XSPEC requires spectra, backgrounds, responses and models, but places fewer constraints as to how they are represented on disk and how they are combined to compute the objective function (statistic). This has immediate implications for the extension of XSPEC analysis to future missions. New data formats can be implemented independently of the existing code, so that they may be loaded during program execution. The “data format” includes the specification not only of the files on disk but how they combine with models. Multiple sources may be extracted from a spectrum. For example, it generalizes the fitting problem to minimizing, in the case of the χ2 statistic, χ2 = X X (Ci − Rij × Mj )2 i (2.9) j and j runs over one or more models. This allows the analysis of, for instance, coded aperture data where multiple sources may be spatially resolved. Responses, which abstractly represent a mapping from the theoretical energy space of the model to the detector channel space, may be represented in new ways. For example, the INTEGRAL/SPI responses are implemented as a linear superposition of 3 (fixed) components. Instead of explicitly combining responses and models through convolution XSPEC places no prior constraint on how this combination is implemented. For example, analysis of data collected by future large detectors might take advantage of the form of the instrumental response by decomposing the response into components of different frequency. Other differences of approach are in the selection of the statistic of the techniques used for deriving the solution. Statistics and fitting methods may be added to XSPEC at execution time rather than at installation time, so that the analysis package as a whole may more easily keep apace of new techniques. 14 CHAPTER 2. SPECTRAL FITTING AND XSPEC 2.6 XSPEC Data Analysis XSPEC is designed to support multiple input data formats. Support for the earlier SF and Einstein FITS formats are removed. Support for ASCII data is planned, which will allow XSPEC to analyze spectra from other wavelength regions (optical, radio) transparently to the user. OGIP Data The OGIP data format both for single spectrum files (Type I) and multiple spectrum files (Type II) is fully supported. These files can be created and manipulated with programs described in Appendix E and the provided links. The programs are described more fully in George et al. 1992. (the directories below refer to the HEAsoft distribution). • Spectral Data: callib/src/gen/rdpha2.f, wtpha3.f • Responses: callib/src/gen rdebd4.f, rdrmf5.f, wtebd4.f, wtrmf5.f. The “rmf” programs read and write the RMF extension, while the “ebd” programs write an extension called EBOUNDS that contains nominal energies for the detector channels. • Auxiliary Responses: callib/src/gen rdarf1.f, and wtarf1.f INTEGRAL/SPI Data XSPEC also includes an add-in module to read and simulate INTEGRAL/SPI data, which can be loaded by the user on demand. The INTEGRAL/SPI datasets are similar to OGIP Type II, but contain an additional FITS extension that stores information on the multiple files used to construct the responses. The INTEGRAL Spectrometer (SPI) is a coded-mask telescope, with a 19element Germanium detector array. The Spectral resolution is 500, and the angular resolution is 3deg. Unlike focusing instruments however, the detected photons are not directionally tagged, and a statistical analysis procedure, using for example cross-correlation techniques, must be employed to reconstruct an image. The description of the XSPEC analysis approach which follows assumes that an image reconstruction has already been performed; see the SPIROS utility within the INTEGRAL offline software analysis package (OSA), OR, the positions on the sky of all sources to be analyzed are already known (which is often the case). Those unfamiliar with INTEGRAL data analysis should refer to the OSA documentation. Thus, the INTEGRAL/SPI analysis chain must be run up to the event binning level [if the field of view (FoV) source content is known, e.g. from published catalogs, or from IBIS image analysis], or the image 2.6. XSPEC DATA ANALYSIS 15 reconstruction level. SPIHIST should be run selecting the “PHA” output option, and selecting detectors 0-18. This will produce an OGIP standard type-II PHA spectral file, which contains multiple, detector count spectra. In addition, the SPIARF procedure should be run once for each source to be analyzed, plus one additional time to produce a special response for analysis of the instrumental background. If this is done correctly, and in the proper sequence, SPIARF will create a table in the PHA-II spectral file, which will associate each spectrum with the appropriate set of response matrices. The response matrices are then automatically loaded into XSPEC upon execution of the data command in a manner very transparent to the user. You will also need to run SPIRMF (unless you have opted to use the default energy bins of the template SPI RMFs). Finally, you will need to run the FTOOL SPIBKG INIT. Each of these utilities - SPIHIST, SPIARF, SPIRMF and SPIBKG INIT - are documented elsewhere, either in the INTEGRAL or (for SPIBKG INIT the HEAsoft) software documentation. There are several complications regarding the spectral de-convolution of codedaperture data. One already mentioned is the source confusion issue; there may be multiple sources in the FoV, which are lead to different degrees of shadowing on different detectors. Thus, a separate instrumental response must be applied to a spectral model for each possible source, for each detector. This is further compounded by the fact that INTEGRAL’s typical mode of observation is “dithering.” A single observation may consist of 10’s of individual exposures at raster points separated by ∼ 2 deg. This further enumerates the number of individual response matrices required for the analysis. If there are multiple sources in the FoV, then additional spectral models can be applied to an additional set of response matrices, enumerated as before over detector and dither pointing. This capability - to model more than one source at a time in a given χ2 (or alternative) minimization procedure - did not exist in previous versions of XSPEC. For an observation with the INTEGRAL/SPI instrument, where the apparent detector efficiency is sensitive to the position of the source on the sky relative to the axis of the instrument, the statistic is: P P X X X Dd,p (I) − j E Rd,p;j × Mj (E; xs ) − Bd,p (I; xb ) 2 χ = ( )2 (2.10) σd,p (I) P d I where p, d run over instrument pointings and detectors; I runs over individual detector channels; j enumerates the sources detected in the field at different position (θ, φ); E indexes the energies in the source model; xs are the parameters of the source model, which is combined with the response; and xb are parameters of the background model. Examination of this equation reveals one more complication; the term B represents the background, which, unlike for chopping, scanning or imaging experiments, must be solved for simultaneously with the desired source content. 16 CHAPTER 2. SPECTRAL FITTING AND XSPEC The proportion of background-to-source counts for a bright source such as the Crab is 1%. Furthermore, the background varies as a function of detector, and time (dither-points), making simple subtraction implausible. Thus, a model of the background is applied to a special response matrix, and included in the de-convolution algorithm. 2.7 References • Arnaud, K.A., George, I.M., Tennant, A.F., 1992. Legacy, 2, 65. • Avni, Y., 1976. ApJ, 210, 642. • Bevington, P.R., 2002, 3rd Edition. Data Reduction and Error Analysis for the Physical Sciences, McGraw-Hill. • Blissett, R.J., Cruise, A.M., 1979. MNRAS, 186, 45. • George, I.M., Arnaud, K.A., Pence, W., Ruamsuwan, L., 1992. Legacy, 2, 51. • Kahn, S.M., Blissett, R.J., 1980. ApJ, 238, 417. • Loredo, T.J., Epstein, R.I., 1989. ApJ, 336, 896. • Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes (2nd edition), p687ff, CUP. • Wheaton, W.A. et.al., 1995. ApJ, 438, 322. Chapter 3 XSPEC Overview and Helpful Hints 3.1 Syntax XSPEC is a command-driven, interactive program. You will see a prompt XSPEC12> whenever input is required. Command recall and inline editing are available using the arrow keys. XSPEC uses Tcl as its user interface, providing looping, conditionals, file I/O and so on. For further details of the Tcl syntax, consult the “Description of Syntax” section, appendix A, and links therein. 3.2 How to return to the XSPEC12>prompt The string /* acts as an emergency escape back to the XSPEC prompt. This string in answer to any question should bounce XSPEC out of whatever it is doing and back to the command prompt. 17 18 CHAPTER 3. XSPEC OVERVIEW AND HELPFUL HINTS 3.3 Getting Help Quick help: If you are uncertain about command syntax, typing a command followed by a “?” will print a one-line summary. The help command: XSPEC12> help without arguments will bring up the full XSPEC manual in a PDF document reader (e.g. AdobeTM Acrobat Reader), or will open a browser to the XSPEC manual home page either locally or on the HEASARC site. See “Customizing XSPEC” later in this chapter to see how to select between these options, and how to assign a PDF reader and web browser to XSPEC. Typing XSPEC12> helpwill display the manual section corresponding to ¡command¿. Help for individual model components can be displayed by XSPEC12> help model if all else fails you can e-mail your questions to the XSPEC team at xspec12@athena.gsfc.nasa.gov 3.4 Commands XSPEC commands can be divided into 6 categories: Control, Data, Model, Fitting, Plotting and Setting, as follows: • Control commands include items such as controlling logging, obtaining help, executing scripts, and other miscellaneous items to do with the program control rather than manipulating data or theoretical models. • Data commands load spectral data and calibration data such as backgrounds and responses, and specify channel ranges to be fit. • Model commands define and manipulate theoretical models and their parameters, and compute additional information such as fluxes or line identifications. 3.5. ISSUING COMMANDS 19 • Fit commands initiate the fitting routines, control the parameter set, perform statistical tests and compute confidence levels. • Plot commands generate about 50 different kinds of 2-dimensional plots • Setting commands change a variety of XSPEC internals which control details of models, statistics, and fitting methods. These command types are summarized below. For full details see Chapter 5. 3.5 Issuing Commands In an interactive session, the command interpreter accepts the shortest unambiguous abbreviation for any command. Since the interpreter is built on the Tcl language, some possible XSPEC command abbreviations might coincide with both XSPEC and Tcl commands. In this case, the interpreter will print the multiple possibilities and stop. Command abbreviations may be specified in a start-up script ($HOME/.xspec/xspec.rc), see Appendix A for details. Inside a script, command abbreviations are not recognized. This behavior can be overridden but we do not recommended it: however, specific abbreviations can be defined in the custom startup script. Operating-system commands can be given from within XSPEC simply by typing the command, as at the shell prompt: however, if wild cards are needed (e.g. ls *.pha) the operating system command must be preceded by syscall. Note that an XSPEC abbreviation which corresponds to a system command will run the latter. 3.6 Control Commands Control commands include those for: • controlling parallel operations: parallel • writing program information: log, (save session to an ASCII file) script (record a set of commands), save (save commands needed to restore the current program state), autosave (automatically run the save command at specified intervals); 20 CHAPTER 3. XSPEC OVERVIEW AND HELPFUL HINTS • controlling program output: chatter (control the amount of program output), query (give an automatic response to prompts), tclout and tcloutr (create Tcl variables for manipulation in scripts) • displaying status information: show, time, and version • ending the session: exit (or quit) • displaying online help help and ?. Help can be given either in summary, in specific manual pages, a manual section, or the entire manual. 3.7 Query, chatter and shutting XSPEC up (somewhat) The fit command will run a certain number of iterations and then query the user whether he or she wants to continue. While useful under most circumstances, the constant questioning can be a pain if one leaves XSPEC running and expects to find it finished when one gets back, or if one habitually runs scripts. One way around this problem is to reset the number of iterations before the question is asked by giving a single argument. For example, XSPEC12> fit 100 will run 100 iterations before asking a question. A more drastic solution is to use the query command. XSPEC12> query yes This will suppress all questions and assume that their answer is yes, while XSPEC12> query no will suppress all questions but assume that their answer is no. The amount of output that XSPEC writes is set by the chatter command, which takes two arguments applying to the terminal and to the log file. 3.8. SCRIPTS AND THE SAVE COMMAND 3.8 21 Scripts and the Save command XSPEC commands can be written into a file and then executed by entering XSPEC12> @filename Alternatively, from the shell prompt % xspec - filename & for batch execution (remember to end the script in file filename.xcm with exit if you want the program to terminate after completion!). Note that the default suffix for xspec scripts is .xcm The save command writes the current XSPEC status to a command file, which later can be run to reset XSPEC to the same configuration. XSPEC has a mechanism to save the current status automatically. This is controlled through the autosave command. This command is very useful when reading a large number of data sets and/or fitting complicated models. If autosaving is operating (the default) then the equivalent of XSPEC12> save all xautosav.xcm is run after each command, so if a disaster occurs it is possible to recover. 3.9 Miscellaneous Information on the current XSPEC status can be printed out using the show command. The time command writes out system-timing information, and the version command writes out the version number and the build time and date. Finally, XSPEC can be terminated with the exit or quit commands. 3.10 Data Commands XSPEC is designed to allow complicated, multi-instrument analysis, so most commands can take arguments specifying more than one data set. Arguments 22 CHAPTER 3. XSPEC OVERVIEW AND HELPFUL HINTS in XSPEC are separated by either blanks or commas. A single argument can define a range. The ranges are delimited by a dash (-). A colon (:) is used to separate ranges (e.g., the phrase 1-2:11-24 refers to channels 11-24 in files 1 and 2). Reading data and modifying calibration and auxiliary files XSPEC reads in spectra from spectral files using the data command. Several datasets may be specified in one command. Several datasets may be stored in a single file and accessed separately. A particular dataset in use may be replaced by another or dropped entirely. The input data file contains pointers to background, redistribution and auxiliary response files, but these pointers may be overridden by the backgrnd, response, and arf commands. All these commands have the same syntax as data. An auxiliary background file, called the correction file (an absolute subtraction with zero variance), also can be included using the corfile command. Its use is described in the section on fitting. The current response can be replaced by a diagonal version using diagrsp. A dummy response for testing purposes can be defined using dummyrsp. Controlling channels being fitted PHA channels may be left out of fitting using the ignore command and included again using the notice command. These commands have a syntax allowing the same channels to be specified for more than one input file. The ignored and noticed ranges can be specified either as channels or as energies. If the command setplot wave has been entered, real ranges are interpreted as wavelengths. Simulations The fakeit command is used to generate simulated data. The current response matrix and model (a model must be defined prior to using the fakeit command) are used to create fake data. The user is prompted for various options. To make fake data when only a response matrix is available, give the command XSPEC12> fakeit none. XSPEC will prompt the user for the response and ancillary filenames from which to build the simulated data. It is important to note that a model must be defined prior to issuing this command. Data groups The most common use of XSPEC is to fit one or more data sets with responses to a particular model. However, it is often useful to be able to fit simultaneously several data sets with a model whose parameters can be different for each data 3.11. MODEL COMMANDS 23 set. A simple example would be a number of data sets that we expect to have the same model spectrum shape but different normalizations. XSPEC caters to this need through the use of data groups. When files are read in they can be labeled as belonging to a particular data group. When a model is defined a set of model parameters is allocated for each data group. These parameters can all vary freely or they can be linked together across data groups as required. To set up data groups, the data command should be given as in the following example : XSPEC12> data 1:1 file1 1:2 file2 2:3 file3 which sets up two data groups. The first data group comprises data sets from file1 and file2, and the second data group takes the data set from file3. Now when a model is defined, XSPEC will give two sets of model parameters, one for the first datagroup and one for the second. 3.11 Model Commands XSPEC allows users to fit data with models constructed from individual components. These components may be either additive, multiplicative, mixing, or convolution. Multiplicative components simply multiply the model by an energy-dependent factor. Convolutions apply a transformation to the model component they operated on whereby the output can be affected by a range of input energies, such as in smoothing operations. Mixing components are two dimensional and designed to transform fluxes between different spatial regions (such as in projection). Multiplicative, and convolution components can act on individual components, on groups of components, or on the entire model, whereas mixing transformations apply to the whole model. The model command defines the model to be used and prompts for the starting values of its parameters. The user also can set the allowed ranges of the parameter. Parameters can be linked to an algebraic function of the other parameters, and unlinked using the untie command. The value of an individual parameter can be changed with the command newpar (and the current setting queried with newpar 0). Parameters can be fixed at their current value with the freeze command and allowed to vary freely with the thaw command. Individual components can be added or subtracted from the model using addcomp, delcomp, and editmod. The plasma emission and photoelectric absorption models require an assumption about relative elemental abundances. The flux command calculates the flux from the current model in the given 24 CHAPTER 3. XSPEC OVERVIEW AND HELPFUL HINTS energy range. This energy range must be within that defined by the current response matrix. If a larger energy range is required, then the energies command can be given to compute the model over the desired range. The lumin command calculates the luminosity for the source redshift given. The eqwidth command determines the equivalent width of a model component, usually a line. The user of either of these last two commands should read the help descriptions carefully. The Tcl script addline can be used to automatically add lines to a model. These can be identified using identify and modid. New model components which can be described by a simple algebraic formula can be set up using mdefine and used in the same way as the standard models except they will run slower being interpreted rather than compiled. Models with multiple responses and background models Multiple models and responses can be assigned to a single spectrum. This generalizes and replaces the “/b” technique of specifying background models in v11. In the FITS file format, a single response file can be associated with a spectrum either through a header keyword or a table column entry. XSPEC always assigns this response to a spectrum’s source number 1. The model command by default also creates new models for source number 1. The response command in tandem with model can be used to create additional sources. For example to add a background model to loaded spectrum 1, first load a 2nd response: XSPEC12> response 2:1 resp2.rsp then define a background model to apply to source 2: XSPEC12> model 2:my_background_model_name wa(po) This model will now apply to spectrum 1 and any other spectrum that has a response loaded for source 2. To apply a different background model to spectrum 2, load a response for source 3 rather than 2: XSPEC12> response 3:2 another_response.rsp XSPEC12> model 3:another_background_model ga An arf can also be assigned to a particular source number and spectrum: XSPEC12> arf 2:1 arf_file.pha Source numbers do not need to be entered in consecutive order for a given spectrum, and gaps in numbering are allowed. Please see the individual model 3.12. FITTING COMMANDS 25 and response entries in the “XSPEC Commands” section for more information and examples. 3.12 Fitting Commands The basic fit command is called fit. This command performs a minimization using the currently selected algorithm (default: Levenberg-Marquardt). fit takes arguments that are passed to the fitting method: by default, these are the number of iterations to execute before asking the user whether to continue, and the numerical convergence criterion. A systematic model uncertainty can be included using the systematic command. The error or uncertain command calculates error bounds for one interesting parameter for the specified parameters and confidence levels. To produce multi-dimensional errors the steppar command is used to generate a fit-statistic grid. Two-dimensional grids may be expressed as contour plots (using plot contour). The model normalization can be set using the renorm command. The normalization of the correction file background can be set with cornorm. ftest and the Tcl script simftest can be used to calculate F-test probabilities. Markov Chain Monte Carlo runs can be performed using the chain command. The proposal distribution covariance can be rescaled using chain rescale if the Metropolis-Hastings algorithm is selected. If an MCMC chain has been loaded then error uses this chain. The analog of steppar is the margin command with the results being plotted using plot margin or plot integprob. What to do when you have Poisson data The χ2 statistic assumes that all the spectral channels are Gaussian distributed and that the estimate for the variance is uncorrelated with the observed counts. If the data are Poisson then these are bad assumptions especially if there are small numbers of counts in a channel. An alternative fit statistic, the C-statistic, should be used in this case. The C-statistic can also provide confidence intervals in exactly the same way as χ2 . To use, give the command XSPEC12> statistic cstat and then use the fit and error commands as usual. An alternative (and deprecated) approach is to continue using the χ2 statistic but change the weighting to provide a better estimate of the variance in the small number limit. This can be done using the weight gehrels or weight churazov commands. The 26 CHAPTER 3. XSPEC OVERVIEW AND HELPFUL HINTS latter is to be preferred. The goodness-of-fit statistic can be set using the command statistic test. There are a number of options available. They can be interpreted using the goodness command, which utilizes Monte Carlo methods. 3.13 Binning and Grouping data Often one does not want to use the full resolution of a spectrum, either because the channels over-sample the spectral resolution or because the S/N is low. XSPEC and the associated programs provide a number of ways of handling this. Firstly, the XSPEC command setplot rebin can be used to add channels together in the plot. It is important to realize that this effects only the plot and not the data being fitted. Two FTOOLS are available to bin and group data for fitting purposes. RBNPHA bins up the data in a non-reversible manner and should only be used to ensure that the number of bins in the spectrum is the same as that in the response. GRPPHA is the tool of choice for grouping the data to get adequate S/N or number of counts in each channel. GRPPHA does not actually add together channels, but instead sets a flag which is read by XSPEC and causes XSPEC to sum the appropriate channels. If a data file is read with some grouping then XSPEC will apply the same operation to any background or response files used. 3.14 Plotting Commands XSPEC plotting is currently performed using the PLT interface. There are two basic commands: plot and iplot. The plot command makes a plot and returns the user to the XSPEC prompt, while the iplot command leaves the user in the interactive plotting interface, thus allowing the user to edit the plot. A variety of different quantities may be plotted, including the data and the current model; the integrated counts; the fit residuals; the ratio of data to model; the contributions to the fit statistic; the theoretical model; the unfolded (incident) spectrum; the detector efficiency; the results of the goodness command; and the fit-statistic contours. All data plots can have an x-axis of channels, energy, or wavelength, which are specified with setplot channel, setplot energy, setplot wave respectively. A number of different units are available for energy or wavelength. The plotting device to be used is set using setplot device or cpd. Separate spectra may be added together and channels binned up (for plotting 3.15. SETTING COMMANDS 27 purposes only) using setplot group (and ungrouped with setplot ungroup) and setplot rebin. There is an option to plot individual additive model components on data plots, this option is enabled by setplot add and disabled by setplot noadd. The effective area can be divided out of data plots using setplot area (which option can be turned off using setplot noarea). Line IDs can be plotted using setplot id and turned off by setplot noid. A stack of PLT commands can be created and manipulated with setplot command, setplot delete, and setplot list. This command stack then is applied to every plot. PLT is built on top of the PGPLOT package, which comes with a standard set of device drivers. Any machine running X-windows should support /xs and /xw, while xterm windows should support /xt. PGPLOT supports monochrome and color postscript and both landscape and portrait orientation with the drivers /ps, /cps, /vps, and /vcps. The easiest way to make a hardcopy of an XSPEC plot is to use XSPEC12> iplot command and then at the PLT prompt to enter PLT> hard /ps This will make a file called pgplot.ps which can be printed. Alternatively, the sequence XSPEC12> cpd /ps XSPEC12> ... plot commands ... XSPEC12> cpd none will place the plots in a PostScript file . 3.15 Setting Commands The fit and goodness-of-fit test statistics are set using the statistic command. Other fit-minimization algorithms are available, and can be selected using the method command. The various fit methods require first and in some cases second derivatives of the statistic with respect to the parameters. By default XSPEC calculates these analytically, using an approximation for the second derivatives. This may be changed by setting the USE NUMERICAL DIFFERENTIATION 28 CHAPTER 3. XSPEC OVERVIEW AND HELPFUL HINTS flag in the user’s startup Xspec.init file. The weighting algorithm used to calculate χ2 can be altered by the weight command. Other setting commands modify: • cosmological parameters used to calculate luminosity (cosmo) • solar abundances for 18 elements (abund) • photoionization cross-sections (xsect) The xset command can be used as an interface for abund, cosmo, method, statistic, and xsect. Additionally, xset may set string expressions that are used by models, for example the path to, and version number of AtomDB atomic line calculations, or the coordinate system for surface brightness calculations used in the xmmpsf mixing model. 3.16 Breaking With Ctrl-C Ctrl-C can be used to break out of the data, chain, error, fit, and steppar commands. If a Ctrl-C is entered elsewhere, it will have no effect. When a break is entered during the fitting commands (error, fit, and steppar), the fit will proceed until the end of the current fit iteration (ie. current lambda value when using Levenberg-Marquardt) before breaking. This is to ensure the program remains in a stable well-defined state. Therefore on slower machines, a user may notice a slight delay before the program actually breaks. Ctrl-C breaking is currently only implemented for the Levenberg-Marquardt fitting method. Breaking is implemented for the data command primarily for users who load a large number of Type-II spectra with one data command. So if you enter XSPEC12> data my_data{1-1000} and decide it is taking too long to load, you can break out at any time. However, if you do choose to break, all spectra loaded from that particular data set will be lost. For example, if the command below is entered and a Ctrl-C is sent while the spectra from my data2 are loading, the 50 spectra from my data1 will be retained while none will be from my data2: XSPEC12> data my_data1{1-50} mydata2{1-50} 3.17. CUSTOMIZING XSPEC 3.17 29 Customizing XSPEC The XSPEC environment can be customized using two separate files, both of which are searched for in the directory $HOME/.xspec The first file, Xspec.init contains a number of settings that control items in XSPEC. The settings available are as follows. USE ONLINE HELP LOCAL HELP FORMAT PDF COMMAND HTML COMMAND GUI DUMMY CHAT XSECT ABUND LEVEL STATISTIC WEIGHT USE NUMERICAL DIFFER COSMO GRAPH PLOTDEVICE WAVE PLOT UNITS USER SCRIPT DIRECTORY FIT DELTAS ATOMDB VERSION PARALLEL CONTOUR IMAGE whether to use on-line or local help files local help format (html or pdf) command to read pdf command to read html whether to use GUI mode (false only) argument for dummyrsp command terminal and log chatter levels photoelectric cross-section table solar abundance table fitting method fit statistic weighting for S 2 statistic calculate parameter derivatives numerically cosmology parameters (H0 , q0 , λ0 ) graphics package (plt only) plotting device units for setplot wave mode directory for user scripts choice for parameters’ fit delta values the AtomDB version number Parallel processing settings Put background image on contour plots A copy of this file is placed in the $HOME/.xspec directory on XSPEC12’s first start-up or when it is not found. After this users can modify settings such as default cosmology or the energy range for dummy response for their own requirements. This is also the place where users can select if they want to view PDF help files from the XSPEC distribution or HTML either locally or from the HEASARC site. Setting USE ONLINE HELP to true uses the remote HTML files while false will use either PDF or HTML local files depending on the value of LOCAL HELP FORMAT. The PDF COMMAND and HTML COMMAND entries determine which applications are run for the two viewing cases. The HTML COMMAND value should typically just be the name of a web browser, or “open” for Mac OS X 30 CHAPTER 3. XSPEC OVERVIEW AND HELPFUL HINTS users. The default settings for the PDF COMMAND entry are what work best for launching Adobe Acrobat Reader 7.0.x on Linux/Unix systems. For those launching earlier versions, the “-openInNewWindow” flag should be replaced with “-useFrontEndProgram”. For Mac users, again we recommend the entire entry simply be replaced with “open”. The second file that is searched for is the xspec.rc file. This contains users’ own customizations, for example Tcl or XSPEC command abbreviations, packages to be loaded on startup, or Tcl scripts containing procedures that are to be executed as commands. Please consult Appendix A and references/links therein for details of Tcl commands and scripting. 3.18 Customizing system-wide When an XSPEC build is intended for many users across a system, it is also possible for the installer (or whoever has write access to the distribution and installation areas) to globally customize XSPEC. This is done through the file global customize.tcl, located in the Xspec/src/scripts directory. (This was done in the xspec.tcl file prior to v12.2.1) Any of the customizations mentioned above for the individual’s own xspec.rc file can also be placed in the global customize.tcl file. After making the additions, run “hmake install” out of the Xspec/src/scripts directory in order to copy the modified global customize.tcl file to the installation area. This additional code will be executed for all users upon startup, BEFORE any of their own customizations in their xspec.rc files. Chapter 4 Walks through XSPEC 4.1 Introduction This chapter demonstrates the use of XSPEC. The brief discussion of data and response files is followed by fully worked examples using real data that include all the screen input and output with a variety of plots. The topics covered are as follows: defining models, fitting data, determining errors, fitting more than one set of data simultaneously, simulating data, and producing plots. 4.2 Brief Discussion of XSPEC Files At least two files are necessary for use with XSPEC: a data file and a response file. In some cases, a file containing background may also be used, and, in rare cases, a correction file is needed to adjust the background during fitting. If the response is split between an rmf and an arf then an ancillary response file is also required. However, most of the time the user need only specify the data file, as the names and locations of the correct response and background files should be written in the header of the data file by whatever program created the files. 31 32 4.3 CHAPTER 4. WALKS THROUGH XSPEC Fitting Models to Data: An Old Example from EXOSAT Our first example uses very old data which is much simpler than more modern observations and so can be used to better illustrate the basics of XSPEC analysis. The 6s X-ray pulsar 1E1048.1-5937 was observed by EXOSAT in June 1985 for 20 ks. In this example, we’ll conduct a general investigation of the spectrum from the Medium Energy (ME) instrument, i.e. follow the same sort of steps as the original investigators (Seward, Charles & Smale, 1986). The ME spectrum and corresponding response matrix were obtained from the HEASARC On-line service. Once installed, XSPEC is invoked by typing %xspec as in this example: %xspec XSPEC version: 12.9.0 Build Date/Time: Sun Jun 28 17:58:18 2015 XSPEC12>data s54405.pha 1 spectrum in use Spectral Data File: s54405.pha Spectrum 1 Net count rate (cts/s) for Spectrum:1 3.783e+00 +/- 1.367e-01 Assigned to Data Group 1 and Plot Group 1 Noticed Channels: 1-125 Telescope: EXOSAT Instrument: ME Channel Type: PHA Exposure Time: 2.358e+04 sec Using fit statistic: chi Using test statistic: chi Using Response (RMF) File s54405.rsp for Source 1 The data command tells the program to read the data as well as the response file that is named in the header of the data file. In general, XSPEC commands can be truncated, provided they remain unambiguous. Since the default extension of a data file is .pha, and the abbreviation of data to the first two letters is unambiguous, the above command can be abbreviated to da s54405, if desired. To obtain help on the data command, or on any other command, type help command followed by the name of the command. One of the first things most users will want to do at this stage – even before fitting models – is to look at their data. The plotting device should be set 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT33 first, with the command cpd (change plotting device). Here, we use the pgplot X-Window server, /xs. XSPEC12> cpd /xs There are more than 50 different things that can be plotted, all related in some way to the data, the model, the fit and the instrument. To see them, type: XSPEC12> plot ? plot data/models/fits etc Syntax: plot commands: background chain chisq contour data delchi dem emodel efficiency eufspec eeufspec foldmodel icounts insensitivity lcounts ldata model ratio residuals sensitivity ufspec Multi-panel plots are created by entering multiple commands e.g. "plot data chisq" counts eemodel goodness margin sum The most fundamental is the data plotted against instrument channel (data); next most fundamental, and more informative, is the data plotted against channel energy. To do this plot, use the XSPEC command setplot energy. Figure 4.1 shows the result of the commands: XSPEC12> setplot energy XSPEC12> plot data Note the label on the y-axis. The word “normalized” indicates that this plot has been divided by the value of the EFFAREA keyword in the response file. Usually this is unity so can be ignored. The label also has no cm-2 so the plot is not corrected for the effective area of the detector. We are now ready to fit the data with a model. Models in XSPEC are specified using the model command, followed by an algebraic expression of a combination of model components. There are two basic kinds of model components: additive, which represent X-Ray sources of different kinds. After being convolved with the instrument response, the components prescribe the number of counts per energy bin (e.g., a bremsstrahlung continuum); and multiplicative models components, which represent phenomena that modify the observed X-Radiation (e.g. reddening or an absorption edge). They apply an energy-dependent multiplicative factor to the source radiation before the result is convolved with the instrumental response. 34 CHAPTER 4. WALKS THROUGH XSPEC Figure 4.1: The result of the command plot data when the data file in question is the EXOSAT ME spectrum of the 6s X-ray pulsar 1E1048.1–5937 available from the HEASARC on-line service. 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT35 More generally, XSPEC allows three types of modifying components: convolutions and mixing models in addition to the multiplicative type. Since there must be a source, there must be least one additive component in a model, but there is no restriction on the number of modifying components. To see what components are available, just type model : XSPEC12>model Additive Models: agauss apec bexriv bkn2pow bvvapec c6mekl cevmkl cflow compbb compmag disk diskbb diskpbb diskpn expdec ezdiskbb grbm kerrbb logpar lorentz nlapec npshock nsmaxg nsx pegpwrlw pexmon powerlaw pshock sedov sirf vapec vbremss vmeka vmekal vrnei vsedov vvpshock vvrnei zgauss zpowerlw bapec bknpower c6pmekl compLS comptb diskir eplogpar gadem kerrd meka nsa nteea pexrav raymond smaug vequil vnei vvapec vvsedov bbody bmc c6pvmkl compPS compth diskline eqpair gaussian kerrdisk mekal nsagrav nthComp pexriv redge srcut vgadem vnpshock vvgnei zagauss bbodyrad bremss c6vmekl compST cplinear diskm eqtherm gnei laor mkcflow nsatmos optxagn plcabs refsch sresc vgnei vpshock vvnei zbbody bexrav bvapec cemekl compTT cutoffpl disko equil grad laor2 nei nsmax optxagnf posm rnei step vmcflow vraymond vvnpshock zbremss Multiplicative Models: SSS_ice TBabs cabs constant expfac gabs notch pcfabs redden smedge varabs vphabs zbabs zdust zphabs zredden zwabs zwndabs TBgrain cyclabs heilin phabs spexpcut wabs zedge zsmdust zxipcf TBvarabs dust highecut plabs spline wndabs zhighect zvarabs absori edge hrefl pwab swind1 xion zigm zvfeabs acisabs expabs lyman recorn uvred zTBabs zpcfabs zvphabs Convolution Models: cflux cpflux kerrconv lsmooth simpl zashift gsmooth partcov zmshift ireflect rdblur kdblur reflect kdblur2 rgsxsrc Mixing Models: ascac projct suzpsf xmmpsf Pile-up Models: 36 CHAPTER 4. WALKS THROUGH XSPEC pileup Additional models are available at : legacy.gsfc.nasa.gov/docs/xanadu/xspec/newmodels.html For information about a specific component, type help model followed by the name of the component): XSPEC12>help model apec Given the quality of our data, as shown by the plot, we’ll choose an absorbed power law, specified as follows : XSPEC12> model phabs(powerlaw) Or, abbreviating unambiguously: XSPEC12> mo pha(po) The user is then prompted for the initial values of the parameters. Entering ¡return¿ or / in response to a prompt uses the default values. We could also have set all parameters to their default values by entering /* at the first prompt. As well as the parameter values themselves, users also may specify step sizes and ranges (value, delta, min, bot, top, and max values), but here, we’ll enter the defaults: XSPEC12>mo pha(po) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001( 0.01) 0 0 100000 1:phabs:nH>/* ======================================================================== Model: phabs<1>*powerlaw<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 1.00000 +/- 0.0 2 2 powerlaw PhoIndex 1.00000 +/- 0.0 3 2 powerlaw norm 1.00000 +/- 0.0 ________________________________________________________________________ Fit statistic : Chi-Squared = 4.864244e+08 using 125 PHA bins. 1E+06 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT37 Test statistic : Chi-Squared = 4.864244e+08 using 125 PHA bins. Reduced chi-squared = 3.987085e+06 for 122 degrees of freedom Null hypothesis probability = 0.000000e+00 Current data and model not fit yet. The current statistic is χ2 and is huge for the initial, default values – mostly because the power law normalization is two orders of magnitude too large. This is easily fixed using the renorm command. XSPEC12> renorm Fit statistic : Chi-Squared = 852.19 using 125 PHA bins. Test statistic : Chi-Squared = 852.19 using 125 PHA bins. Reduced chi-squared = 6.9852 for 122 degrees of freedom Null hypothesis probability = 7.320765e-110 Current data and model not fit yet. We are not quite ready to fit the data (and obtain a better χ2 ), because not all of the 125 PHA bins should be included in the fitting: some are below the lower discriminator of the instrument and therefore do not contain valid data; some have imperfect background subtraction at the margins of the pass band; and some may not contain enough counts for χ2 to be strictly meaningful. To find out which channels to discard (ignore in XSPEC terminology), consult mission-specific documentation that will include information about discriminator settings, background subtraction problems and other issues. For the mature missions in the HEASARC archives, this information already has been encoded in the headers of the spectral files as a list of “bad” channels. Simply issue the command: XSPEC12> ignore bad ignore: 40 channels ignored from Fit statistic : Chi-Squared = source number 1 799.74 using 85 PHA bins. Test statistic : Chi-Squared = 799.74 using 85 PHA bins. Reduced chi-squared = 9.7529 for 82 degrees of freedom Null hypothesis probability = 3.545709e-118 Current data and model not fit yet. XSPEC12> plot ldata chi Giving two options for the plot command generates a plot with vertically stacked windows. Up to six options can be given to the plot command at a time. Forty channels were rejected because they were flagged as bad – but do we need to ignore any more? Figure 4.2 shows the result of plotting the data and the model 38 CHAPTER 4. WALKS THROUGH XSPEC Figure 4.2: The result of the command plot ldata chi after the command ignore bad on the EXOSAT ME spectrum 1E1048.1-5937. 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT39 (in the upper window) and the contributions to χ2 (in the lower window). We see that above about 15 keV the S/N becomes small. We also see, comparing Figure 4.2 with Figure 4.1, which bad channels were ignored. Although visual inspection is not the most rigorous method for deciding which channels to ignore (more on this subject later), it’s good enough for now, and will at least prevent us from getting grossly misleading results from the fitting. To ignore energies above 15 keV: XSPEC12> ignore 15.0-** 78 channels (48-125) ignored in spectrum # Fit statistic : Chi-Squared = 1 721.56 using 45 PHA bins. Test statistic : Chi-Squared = 721.56 using 45 PHA bins. Reduced chi-squared = 17.180 for 42 degrees of freedom Null hypothesis probability = 1.253044e-124 Current data and model not fit yet. If the ignore command is handed a real number it assumes energy in keV while if it is handed an integer it will assume channel number. The “**” just means the highest energy. Starting a range with “**” means the lowest energy. The inverse of ignore is notice, which has the same syntax. We are now ready to fit the data. Fitting is initiated by the command fit. As the fit proceeds, the screen displays the status of the fit for each iteration until either the fit converges to the minimum χ2 , or we are asked whether the fit is to go through another set of iterations to find the minimum. The default number of iterations before prompting is ten. XSPEC12>fit Parameters Chi-Squared |beta|/N Lvl 1:nH 2:PhoIndex 451.814 150.854 -3 0.0961968 1.60719 413.575 63220.2 -3 0.264111 2.30580 53.9398 28104.3 -4 0.517263 2.14118 43.816 4617.17 -5 0.551755 2.23947 43.802 139.682 -6 0.538816 2.23680 43.802 0.58082 -7 0.537846 2.23646 ======================================== Variances and Principal Axes 1 2 3 4.7889E-08| -0.0025 -0.0151 0.9999 8.6827E-02| -0.9153 -0.4026 -0.0084 2.2916E-03| -0.4027 0.9153 0.0128 ---------------------------------------==================================== 3:norm 0.00386478 0.00908043 0.0121356 0.0130926 0.0130394 0.0130320 40 CHAPTER 4. WALKS THROUGH XSPEC Covariance Matrix 1 2 3 7.312e-02 3.115e-02 6.565e-04 3.115e-02 1.599e-02 3.208e-04 6.565e-04 3.208e-04 6.562e-06 -----------------------------------======================================================================== Model phabs<1>*powerlaw<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 0.537846 +/- 0.270409 2 2 powerlaw PhoIndex 2.23646 +/- 0.126458 3 2 powerlaw norm 1.30320E-02 +/- 2.56170E-03 ________________________________________________________________________ Fit statistic : Chi-Squared = Test statistic : Chi-Squared = Reduced chi-squared = Null hypothesis probability = 43.80 using 45 PHA bins. 43.80 using 45 PHA bins. 1.043 for 42 degrees of freedom 3.949507e-01 There is a fair amount of information here so we will unpack it a bit at a time. One line is written out after each fit iteration. The columns labeled Chi-Squared and Parameters are obvious. The other two provide additional information on fit convergence. At each step in the fit a numerical derivative of the statistic with respect to the parameters is calculated. We call the vector of these derivatives beta. At the best-fit the norm of beta should be zero so we write out —beta— divided by the number of parameters as a check. The actual default convergence criterion is when the fit statistic does not change significantly between iterations so it is possible for the fit to end while —beta— is still significantly different from zero. The —beta—/N column helps us spot this case. The Lvl column also indicates how the fit is converging and should generally decrease. Note that on the first iteration only the powerlaw norm is varied. While not necessary this simple model, for more complicated models only varying the norms on the first iteration helps the fit proper get started in a reasonable region of parameter space. At the end of the fit XSPEC writes out the Variances and Principal Axes and Covariance Matrix sections. These are both based on the second derivatives of the statistic with respect to the parameters. Generally, the larger these second derivatives, the better determined the parameter (think of the case of a parabola in 1-D). The Covariance Matrix is the inverse of the matrix of second derivatives. The Variances and Principal Axes section is based on an eigenvector decomposition of the matrix of second derivatives and indicates which parameters are correlated. We can see in this case that the first eigenvector depends 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT41 almost entirely on the powerlaw norm while the other two are combinations of the nH and powerlaw PhoIndex. This tells us that the norm is independent but the other two parameters are correlated. The next section shows the best-fit parameters and error estimates. The latter are just the square roots of the diagonal elements of the covariance matrix so implicitly assume that the parameter space is multidimensional Gaussian with all parameters independent. We already know in this case that the parameters are not independent so these error estimates should only be considered guidelines to help us determine the true errors later. The final section shows the statistic values at the end of the fit. XSPEC defines a fit statistic, used to determine the best-fit parameters and errors, and test statistic, used to decide whether this model and parameters provide a good fit to the data. By default, both statistics are χ2 . When the test statistic is χ2 we can also calculate the reduced χ2 and the null hypothesis probability. This latter is the probability of getting a value of χ2 as large or larger than observed if the model is correct. If this probability is small then the model is not a good fit. The null hypothesis probability can be calculated analytically for χ2 but not for some other test statistics so XSPEC provides another way of determining the meaning of the statistic value. The goodness command performs simulations of the data based on the current model and parameters and compares the statistic values calculated with that for the real data. If the observed statistic is larger than the values for the simulated data this implies that the real data do not come from the model. To see how this works we will use the command for this case (where it is not necessary) XSPEC12>goodness 1000 50.40% of realizations are < best fit statistic 43.80 XSPEC12>plot goodness (nosim) Approximately half of the simulations give a statistic value less than that observed, consistent with this being a good fit. Figure 4.3 shows a histogram of the χ2 values from the simulations with the observed value shown by the vertical dotted line. So the statistic implies the fit is good but it is still always a good idea to look at the data and residuals to check for any systematic differences that may not be caught by the test. To see the fit and the residuals, we produce figure 4.4 using the command XSPEC12>plot data resid Now that we think we have the correct model we need to determine how well the parameters are determined. The screen output at the end of the fit shows the 42 CHAPTER 4. WALKS THROUGH XSPEC Figure 4.3: The result of the command plot goodness. The histogram shows the fraction of simulations with a given value of the statistic and the dotted line marks that for the observed data. There is no reason to reject the model. 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT43 Figure 4.4: The result of the command plot data resid with: the ME data file from 1E1048.1-5937; “bad” and negative channels ignored; the best-fitting absorbed power-law model; the residuals of the fit. 44 CHAPTER 4. WALKS THROUGH XSPEC best-fitting parameter values, as well as approximations to their errors. These errors should be regarded as indications of the uncertainties in the parameters and should not be quoted in publications. The true errors, i.e. the confidence ranges, are obtained using the error command. We want to run error on all three parameters which is an intrinsically parallel operation so we can use XSPEC’s support for multiple cores and run the error estimations in parallel: XSPEC12>parallel error 3 XSPEC12>error 1 2 3 Parameter Confidence Range (2.706) 1 0.107599 1.00722 (-0.430231,0.469393) 2 2.03775 2.44916 (-0.198718,0.212699) 3 0.00954178 0.0181617 (-0.00349016,0.00512979) Here, the numbers 1, 2, 3 refer to the parameter numbers in the Model par column of the output at the end of the fit. For the first parameter, the column of absorbing hydrogen atoms, the 90% confidence range is . This corresponds to an excursion in of 2.706. The reason these “better” errors are not given automatically as part of the fit output is that they entail further fitting. When the model is simple, this does not require much CPU, but for complicated models the extra time can be considerable. The error for each parameter is determined allowing the other two parameters to vary freely. If the parameters are uncorrelated this is all the information we need to know. However, we have an indication from the covariance matrix at the end of the fit that the column and photon index are correlated. To investigate this further we can use the command steppar to run a grid over these two parameters: XSPEC12>steppar 1 0.0 1.5 25 2 1.5 3.0 25 Chi-Squared Delta Chi-Squared 162.65 118.84 171.59 127.79 180.87 137.06 190.44 146.64 200.29 156.49 . . . . . . . 316.02 272.22 334.68 290.88 354.2 310.4 374.62 330.82 395.94 352.14 nH 1 PhoIndex 2 0 1 2 3 4 0 0.06 0.12 0.18 0.24 0 0 0 0 0 1.5 1.5 1.5 1.5 1.5 4 3 2 1 0 0.24 0.18 0.12 0.06 0 25 25 25 25 25 3 3 3 3 3 and make the contour plot shown in figure 4.5 using: 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT45 Figure 4.5: The result of the command plot contour. The contours shown are for one, two and three sigma. The cross marks the best-fit position. XSPEC12>plot contour What else can we do with the fit? One thing is to derive the flux of the model. The data by themselves only give the instrument-dependent count rate. The model, on the other hand, is an estimate of the true spectrum emitted. In XSPEC, the model is defined in physical units independent of the instrument. The command flux integrates the current model over the range specified by the user: XSPEC12> flux 2 10 Model Flux 0.003539 photons (2.2321e-11 ergs/cm^2/s) range (2.0000 - 10.000 keV) Here we have chosen the standard X-ray range of 2–10 keV and find that the energy flux is 2.2 × 10−11 ergs/cm2 /s. Note that flux will integrate only within the energy range of the current response matrix. If the model flux outside this range is desired – in effect, an extrapolation beyond the data – then the command energies should be used. This command defines a set of energies on which the model will be calculated. The resulting model is then remapped onto the response energies for convolution with the response matrix. For example, if 46 CHAPTER 4. WALKS THROUGH XSPEC we want to know the flux of our model in the ROSAT PSPC band of 0.2–2 keV, we enter: XSPEC12>energies extend low 0.2 100 Models will use response energies extended to: Low: 0.2 in 100 log bins Fit statistic : Chi-Squared = 43.80 using 45 PHA bins. Test statistic : Chi-Squared = 43.80 using 45 PHA bins. Reduced chi-squared = 1.043 for 42 degrees of freedom Null hypothesis probability = 3.949507e-01 Current data and model not fit yet. XSPEC12>flux 0.2 2. Model Flux 0.0043484 photons (8.8419e-12 ergs/cm^2/s) range (0.20000 - 2.0000 keV) The energy flux, at 8.8 × 10−12 ergs/cm2 /s is lower in this band but the photon flux is higher. The model energies can be reset to the response energies using energies reset. Calculating the flux is not usually enough, we want its uncertainty as well. The best way to do this is to use the cflux model. Suppose further that what we really want is the flux without the absorption then we include the new cflux model by: XSPEC12>editmod pha*cflux(pow) Input parameter value, delta, min, bot, top, and max values for ... 0.5 -0.1( 0.005) 0 0 1e+06 2:cflux:Emin>0.2 10 -0.1( 0.1) 0 0 1e+06 3:cflux:Emax>2.0 -12 0.01( 0.12) -100 -100 100 4:cflux:lg10Flux>-10.3 Fit statistic : Chi-Squared = 52.01 using 45 PHA bins. Test statistic : Chi-Squared = 52.01 using 45 PHA bins. Reduced chi-squared = 1.268 for 41 degrees of freedom Null hypothesis probability = 1.163983e-01 Current data and model not fit yet. ======================================================================== Model phabs<1>*cflux<2>*powerlaw<3> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 0.537843 +/- 0.270399 2 2 cflux Emin keV 0.200000 frozen 1e+06 1e+06 100 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT47 3 4 5 6 2 2 3 3 cflux cflux powerlaw powerlaw Emax lg10Flux PhoIndex norm keV cgs 2.00000 -10.3000 2.23646 1.30320E-02 frozen +/- 0.0 +/- 0.126455 +/- 2.56146E-03 ________________________________________________________________________ The Emin and Emax parameters are set to the energy range over which we want the flux to be calculated. We also have to fix the norm of the powerlaw because the normalization of the model will now be determined by the lg10Flux parameter. This is done using the freeze command: XSPEC12>freeze 6 We now run fit to get the best-fit value of lg10Flux as -10.2903 then: XSPEC12>error 4 Parameter Confidence Range (2.706) 4 -10.458 -10.0789 (-0.167672,0.211462) for a 90% confidence range on the 0.2–2 keV unabsorbed flux of 3.49 × 10−11 – 8.33 × 10−11 ergs/cm2 /s. The fit, as we’ve remarked, is good, and the parameters are constrained. But unless the purpose of our investigation is merely to measure a photon index, it’s a good idea to check whether alternative models can fit the data just as well. We also should derive upper limits to components such as iron emission lines and additional continua, which, although not evident in the data nor required for a good fit, are nevertheless important to constrain. First, let’s try an absorbed black body: XSPEC12>mo pha(bb) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001( 0.01) 0 0 100000 1:phabs:nH>/* ======================================================================== Model phabs<1>*bbody<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 1.00000 +/- 0.0 2 2 bbody kT keV 3.00000 +/- 0.0 3 2 bbody norm 1.00000 +/- 0.0 ________________________________________________________________________ 1e+06 48 Fit statistic : Chi-Squared = CHAPTER 4. WALKS THROUGH XSPEC 3.377094e+09 using 45 PHA bins. Test statistic : Chi-Squared = 3.377094e+09 using 45 PHA bins. Reduced chi-squared = 8.040700e+07 for 42 degrees of freedom Null hypothesis probability = 0.000000e+00 Current data and model not fit yet. XSPEC12>fit Parameters Chi-Squared |beta|/N Lvl 1:nH 2:kT 3:norm 1535.61 63.3168 0 0.334306 3.01647 0.000673086 1523.48 112166 0 0.157480 2.96616 0.000613284 1491.73 170831 0 0.0668724 2.87680 0.000570109 1444.74 204638 0 0.0228530 2.76753 0.000535209 1387.84 226856 0 0.00206016 2.64901 0.000504576 1325.39 243767 0 0.000851460 2.52617 0.000476469 1256.03 258162 0 0.000299777 2.40098 0.000449911 1180.44 271941 0 4.82004e-05 2.27529 0.000425076 1093.77 284819 0 1.99048e-05 2.14630 0.000402842 976.95 291176 0 9.75848e-06 1.99181 0.000380678 Number of trials exceeded: continue fitting? Y ... ... 123.773 25.397 -8 1.87147e-08 0.890295 0.000278599 Number of trials exceeded: continue fitting? ***Warning: Zero alpha-matrix diagonal element for parameter 1 Parameter 1 is pegged at 1.87147e-08 due to zero or negative pivot element, likely caused by the fit being insensitive to the parameter. 123.773 1.92501 -3 1.87147e-08 0.890205 0.000278596 ============================== Variances and Principal Axes 2 3 2.8677E-04| -1.0000 -0.0000 2.2370E-11| 0.0000 -1.0000 -----------------------------======================== Covariance Matrix 1 2 2.867e-04 9.315e-09 9.315e-09 2.267e-11 -----------------------======================================================================== Model phabs<1>*bbody<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 2.24621E-09 +/- -1.00000 2 2 bbody kT keV 0.890178 +/- 1.69312E-02 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT49 3 2 bbody norm 2.78595E-04 +/- 4.76156E-06 ________________________________________________________________________ Fit statistic : Chi-Squared = 123.77 using 45 PHA bins. Test statistic : Chi-Squared = 123.77 using 45 PHA bins. Reduced chi-squared = 2.9470 for 42 degrees of freedom Null hypothesis probability = 5.417154e-10 Note that after each set of 10 iterations you are asked whether you want to continue. Replying no at these prompts is a good idea if the fit is not converging quickly. Conversely, to avoid having to keep answering the question, i.e., to increase the number of iterations before the prompting question appears, begin the fit with, say fit 100. This command will put the fit through 100 iterations before pausing. To automatically answer yes to all such questions use the command query yes. Note that the fit has written out a warning about the first parameter and its estimated error is written as -1. This indicates that the fit is unable to constrain the parameter and it should be considered indeterminate. This usually indicates that the model is not appropriate. One thing to check in this case is that the model component has any contribution within the energy range being calculated. Plotting the data and residuals again we obtain Figure 4.6. The black body fit is obviously not a good one. Not only is χ2 large, but the best-fitting NH is indeterminate. Inspection of the residuals confirms this: the pronounced wave-like shape is indicative of a bad choice of overall continuum. Let’s try thermal bremsstrahlung next: XSPEC12>mo pha(br) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001( 0.01) 0 0 100000 1:phabs:nH>/* ======================================================================== Model phabs<1>*bremss<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 1.00000 +/- 0.0 2 2 bremss kT keV 7.00000 +/- 0.0 3 2 bremss norm 1.00000 +/- 0.0 ________________________________________________________________________ Fit statistic : Chi-Squared = 4.534834e+07 using 45 PHA bins. 1e+06 50 CHAPTER 4. WALKS THROUGH XSPEC Figure 4.6: As for Figure 4.4, but the model is the best-fitting absorbed black body. Note the wave-like shape of the residuals which indicates how poor the fit is, i.e. that the continuum is obviously not a black body. 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT51 Test statistic : Chi-Squared = 4.534834e+07 using 45 PHA bins. Reduced chi-squared = 1.079722e+06 for 42 degrees of freedom Null hypothesis probability = 0.000000e+00 Current data and model not fit yet. XSPEC12>fit Parameters Chi-Squared |beta|/N Lvl 1:nH 2:kT 3:norm 105.28 24.2507 -3 0.273734 6.18714 0.00724161 46.8022 16593.4 -4 0.0371200 5.59937 0.00785588 ... ... 40.0373 270.662 0 7.62629e-05 5.28989 0.00830799 ======================================== Variances and Principal Axes 1 2 3 1.9514E-08| -0.0016 0.0007 1.0000 1.1574E-02| 0.9736 0.2281 0.0014 5.3111E-01| 0.2281 -0.9736 0.0011 ---------------------------------------==================================== Covariance Matrix 1 2 3 3.862e-02 -1.148e-01 1.431e-04 -1.148e-01 5.015e-01 -5.379e-04 1.431e-04 -5.379e-04 6.290e-07 -----------------------------------======================================================================== Model phabs<1>*bremss<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 7.62629E-05 +/- 0.196509 2 2 bremss kT keV 5.28989 +/- 0.708184 3 2 bremss norm 8.30799E-03 +/- 7.93082E-04 ________________________________________________________________________ Fit statistic : Chi-Squared = 40.04 using 45 PHA bins. Test statistic : Chi-Squared = 40.04 using 45 PHA bins. Reduced chi-squared = 0.9533 for 42 degrees of freedom Null hypothesis probability = 5.574343e-01 Bremsstrahlung is a better fit than the black body – and is as good as the power law – although it shares the low NH . With two good fits, the power law and the bremsstrahlung, it’s time to scrutinize their parameters in more detail. 52 CHAPTER 4. WALKS THROUGH XSPEC First, we reset our fit to the powerlaw (output omitted): XSPEC12>mo pha(po) From the EXOSAT database on HEASARC, we know that the target in question, 1E1048.1–5937, has a Galactic latitude of , i.e., almost on the plane of the Galaxy. In fact, the database also provides the value of the Galactic NH based on 21-cm radio observations. At 4 × 1022 cm−2 , it is higher than the 90 percentconfidence upper limit from the power-law fit. Perhaps, then, the power-law fit is not so good after all. What we can do is fix (freeze in XSPEC terminology) the value of NH at the Galactic value and refit the power law. Although we won’t get a good fit, the shape of the residuals might give us a clue to what is missing. To freeze a parameter in XSPEC, use the command freeze followed by the parameter number, like this: XSPEC12> freeze 1 The inverse of freeze is thaw: XSPEC12> thaw 1 Alternatively, parameters can be frozen using the newpar command, which allows all the quantities associated with a parameter to be changed. We can flip between frozen and thawed states by entering 0 after the new parameter value. In our case, we want NH frozen at 4 × 1022 cm−2 , so we go back to the power law best fit and do the following : XSPEC12>newpar 1 Current value, delta, min, bot, top, and max values 0.537843 0.001(0.00537843) 0 1:phabs[1]:nH:1>4,0 Fit statistic : Chi-Squared = 0 100000 823.34 using 45 PHA bins. Test statistic : Chi-Squared = 823.34 using 45 PHA bins. Reduced chi-squared = 19.148 for 43 degrees of freedom Null hypothesis probability = 6.152922e-145 Current data and model not fit yet. The same result can be obtained by putting everything onto the command line, i.e., newpar 1 4, 0, or by issuing the two commands, newpar 1 4 followed by freeze 1. Now, if we fit and plot again, we get the following model (Fig. 4.7). 1e+06 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT53 Figure 4.7: As for Figure 4.4 & 4.6, but the model is the best-fitting power law with the absorption fixed at the Galactic value. Under the assumptions that the absorption really is the same as the 21-cm value and that the continuum really is a power law, this plot provides some indication of what other components might be added to the model to improve the fit. 54 CHAPTER 4. WALKS THROUGH XSPEC XSPEC12>fit ... ======================================================================== Model phabs<1>*powerlaw<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 4.00000 frozen 2 2 powerlaw PhoIndex 3.59784 +/- 6.76670E-02 3 2 powerlaw norm 0.116579 +/- 9.43208E-03 ________________________________________________________________________ Fit statistic : Chi-Squared = 136.04 using 45 PHA bins. The fit is not good. In Figure 4.7 we can see why: there appears to be a surplus of softer photons, perhaps indicating a second continuum component. To investigate this possibility we can add a component to our model. The editmod command lets us do this without having to respecify the model from scratch. Here, we’ll add a black body component. XSPEC12>editmod pha(po+bb) Input parameter value, delta, min, bot, top, and max values for ... 3 0.01( 0.03) 0.0001 0.01 100 4:bbody:kT>2,0 1 0.01( 0.01) 0 0 1e+24 5:bbody:norm>1e-5 Fit statistic : Chi-Squared = 132.76 using 45 PHA bins. Test statistic : Chi-Squared = 132.76 using 45 PHA bins. Reduced chi-squared = 3.1610 for 42 degrees of freedom Null hypothesis probability = 2.387580e-11 Current data and model not fit yet. ======================================================================== Model phabs<1>(powerlaw<2> + bbody<3>) Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 4.00000 frozen 2 2 powerlaw PhoIndex 3.59784 +/- 6.76670E-02 3 2 powerlaw norm 0.116579 +/- 9.43208E-03 4 3 bbody kT keV 2.00000 frozen 5 3 bbody norm 1.00000E-05 +/- 0.0 ________________________________________________________________________ 200 1e+24 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT55 Figure 4.8: The result of the command plot model in the case of the ME data file from 1E1048.1-5937. Here, the model is the best-fitting combination of power law, black body and fixed Galactic absorption. The three lines show the two continuum components (absorbed to the same degree) and their sum. 56 CHAPTER 4. WALKS THROUGH XSPEC Notice that in specifying the initial values of the black body, we have frozen kT at 2 keV (the canonical temperature for nuclear burning on the surface of a neutron star in a low-mass X-ray binary) and started the normalization small. Without these measures, the fit might “lose its way”. Now, if we fit, we get (not showing all the iterations this time): ======================================================================== Model phabs<1>(powerlaw<2> + bbody<3>) Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 4.00000 frozen 2 2 powerlaw PhoIndex 4.89633 +/- 0.158893 3 2 powerlaw norm 0.365391 +/- 5.26695E-02 4 3 bbody kT keV 2.00000 frozen 5 3 bbody norm 2.29745E-04 +/- 2.03915E-05 ________________________________________________________________________ Fit statistic : Chi-Squared = 69.53 using 45 PHA bins. The fit is better than the one with just a power law and the fixed Galactic column, but it is still not good. Thawing the black body temperature and fitting does of course improve the fit but the powerl law index becomes even steeper. Looking at this odd model with the command XSPEC12> plot model We see, in Figure 4.8, that the black body and the power law have changed places, in that the power law provides the soft photons required by the high absorption, while the black body provides the harder photons. We could continue to search for a plausible, well-fitting model, but the data, with their limited signal-to-noise and energy resolution, probably don’t warrant it (the original investigators published only the power law fit). There is, however, one final, useful thing to do with the data: derive an upper limit to the presence of a fluorescent iron emission line. First we delete the black body component using delcomp then thaw NH and refit to recover our original best fit. Now, we add a gaussian emission line of fixed energy and width then fit to get: ======================================================================== Model phabs<1>(powerlaw<2> + gaussian<3>) Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 0.753994 +/- 0.320347 4.3. FITTING MODELS TO DATA: AN OLD EXAMPLE FROM EXOSAT57 2 2 powerlaw PhoIndex 2.38165 +/- 0.166974 3 2 powerlaw norm 1.59131E-02 +/- 3.94933E-03 4 3 gaussian LineE keV 6.40000 frozen 5 3 gaussian Sigma keV 0.100000 frozen 6 3 gaussian norm 7.47374E-05 +/- 4.74253E-05 ________________________________________________________________________ The energy and width have to be frozen because, in the absence of an obvious line in the data, the fit would be completely unable to converge on meaningful values. Besides, our aim is to see how bright a line at 6.4 keV can be and still not ruin the fit. To do this, we fit first and then use the error command to derive the maximum allowable iron line normalization. We then set the normalization at this maximum value with newpar and, finally, derive the equivalent width using the eqwidth command. That is: XSPEC12>err 6 Parameter Confidence Range (2.706) ***Warning: Parameter pegged at hard limit: 0 6 0 0.000151164 (-7.476e-05,7.64036e-05) XSPEC12>new 6 0.000151164 Fit statistic : Chi-Squared = 46.03 using 45 PHA bins. Test statistic : Chi-Squared = 46.03 using 45 PHA bins. Reduced chi-squared = 1.123 for 41 degrees of freedom Null hypothesis probability = 2.717072e-01 Current data and model not fit yet. XSPEC12>eqwidth 3 Data group number: 1 Additive group equiv width for Component 3: 0.784169 keV Things to note: • The true minimum value of the gaussian normalization is less than zero, but the error command stopped searching for a ∆χ2 of 2.706 when the minimum value hit zero, the “hard” lower limit of the parameter. Hard limits can be adjusted with the newpar command, and they correspond to the quantities min and max associated with the parameter values. • The command eqwidth takes the component number as its argument. • The upper limit on the equivalent width of a 6.4 keV emission line is high (784 eV)! 58 CHAPTER 4. WALKS THROUGH XSPEC 4.4 Simultaneous Fitting XSPEC has the very useful facility of allowing models to be fitted simultaneously to more than one data file. It is even possible to group files together and to fit different models simultaneously. Reasons for fitting in this manner include: • The same target is observed at several epochs but, although the source stays constant, the response matrix has changed. When this happens, the data files cannot be added together; they have to be fitted separately. Fitting the data files simultaneously yields tighter constraints. • The same target is observed with different instruments. All the instruments on Suzaku, for example, observe in the same direction simultaneously. As far as XSPEC is concerned, this is just like the previous case: two data files with two responses fitted simultaneously with the same model. • Different targets are observed, but the user wants to fit the same model to each data file with some parameters shared and some allowed to vary separately. For example, if we have a series of spectra from a variable AGN, we might want to fit them simultaneously with a model that has the same, common photon index but separately vary the normalization and absorption. Other scenarios are possible—the important thing is to recognize the flexibility of XSPEC in this regard. As an example we will look at a case of fitting the same model to two different data files but where not all the parameters are identical. Again, this is an older dataset that provides a simpler illustration than more modern data. The massive X-ray binary Centaurus X-3 was observed with the LAC on Ginga in 1989. Its flux level before eclipse was much lower than the level after eclipse. Here, we’ll use XSPEC to see whether spectra from these two phases can be fitted with the same model, which differs only in the amount of absorption. This kind of fitting relies on introducing an extra dimension, the group, to the indexing of the data files. The files in each group share the same model but not necessarily the same parameter values, which may be shared as common to all the groups or varied separately from group to group. Although each group may contain more than one file, there is only one file in each of the two groups in this example. Groups are specified with the data command, with the group number preceding the file number, like this: XSPEC12>data 1:1 losum 2:2 hisum 2 spectra in use 4.4. SIMULTANEOUS FITTING 59 Spectral Data File: losum.pha Spectrum 1 Net count rate (cts/s) for Spectrum:1 1.401e+02 +/- 3.549e-01 Assigned to Data Group 1 and Plot Group 1 Noticed Channels: 1-48 Telescope: GINGA Instrument: LAC Channel Type: PHA Exposure Time: 1 sec Using fit statistic: chi Using test statistic: chi Using Response (RMF) File ginga_lac.rsp for Source 1 Spectral Data File: hisum.pha Spectrum 2 Net count rate (cts/s) for Spectrum:2 1.371e+03 +/- 3.123e+00 Assigned to Data Group 2 and Plot Group 2 Noticed Channels: 1-48 Telescope: GINGA Instrument: LAC Channel Type: PHA Exposure Time: 1 sec Using fit statistic: chi Using test statistic: chi Using Response (RMF) File ginga_lac.rsp for Source 1 Here, the first group makes up the file losum.pha, which contains the spectrum of all the low, pre-eclipse emission. The second group makes up the second file, hisum.pha, which contains all the high, post-eclipse emission. Note that file number is “absolute” in the sense that it is independent of group number. Thus, if there were three files in each of the two groups (lo1.pha, lo2.pha, lo3.pha, hi1.pha, hi2.pha, and hi3.pha, say), rather than one, the six files would be specified as da 1:1 lo1 1:2 lo2 1:3 lo3 2:4 hi1 2:5 hi2 2:6 hi3. The ignore command works on file number, and does not take group number into account. So, to ignore channels 1–3 and 37–48 of both files: XSPEC12> ignore 1-2:1-3 37-48 The model we’ll use at first to fit the two files is an absorbed power law with a high-energy cut-off: XSPEC12> mo phabs * highecut (po) After defining the model, we will be prompted for two sets of parameter values, one for the first group of data files (losum.pha), the other for the second group (hisum.pha). Here, we’ll enter the absorption column of the first group as 1024 cm−2 and enter the default values for all the other parameters in the first group. Now, when it comes to the second group of parameters, we enter a column of 1022 cm− 2 and then enter defaults for the other parameters. The rule being applied here is as follows: to tie parameters in the second group to their equivalents in 60 CHAPTER 4. WALKS THROUGH XSPEC the first group, take the default when entering the second-group parameters; to allow parameters in the second group to vary independently of their equivalents in the first group, enter different values explicitly: XSPEC12>mo phabs*highecut(po) Input parameter value, delta, min, bot, top, and max values for ... Current: 1 0.001 0 0 1E+05 1E+06 DataGroup 1:phabs:nH>100 Current: 10 0.01 0.0001 0.01 1E+06 1E+06 DataGroup 1:highecut:cutoffE> Current: 15 0.01 0.0001 0.01 1E+06 1E+06 DataGroup 1:highecut:foldE> Current: 1 0.01 -3 -2 9 10 DataGroup 1:powerlaw:PhoIndex> Current: 1 0.01 0 0 1E+24 1E+24 DataGroup 1:powerlaw:norm> Current: 100 0.001 0 0 1E+05 1E+06 DataGroup 2:phabs:nH>1 Current: 10 0.01 0.0001 0.01 1E+06 1E+06 DataGroup 2:highecut:cutoffE>/* ======================================================================== Model phabs<1>*highecut<2>*powerlaw<3> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp Data group: 1 1 1 phabs nH 10^22 100.000 +/- 0.0 2 2 highecut cutoffE keV 10.0000 +/- 0.0 3 2 highecut foldE keV 15.0000 +/- 0.0 4 3 powerlaw PhoIndex 1.00000 +/- 0.0 5 3 powerlaw norm 1.00000 +/- 0.0 Data group: 2 6 1 phabs nH 10^22 1.00000 +/- 0.0 7 2 highecut cutoffE keV 10.0000 = 2 8 2 highecut foldE keV 15.0000 = 3 9 3 powerlaw PhoIndex 1.00000 = 4 10 3 powerlaw norm 1.00000 = 5 ________________________________________________________________________ Notice how the summary of the model, displayed immediately above, is different now that we have two groups, as opposed to one (as in all the previous examples). We can see that of the 10 model parameters, 6 are free (i.e., 4 of the second group parameters are tied to their equivalents in the first group). Fitting this model results in a huge χ2 (not shown here), because our assumption that only a change in absorption can account for the spectral variation before and after eclipse is clearly wrong. Perhaps scattering also plays a role in reducing the flux before eclipse. This could be modeled (simply at first) by allowing the normalization of the power law to be smaller before eclipse than after eclipse. 4.4. SIMULTANEOUS FITTING 61 To decouple tied parameters, we change the parameter value in the second group to a value – any value – different from that in the first group (changing the value in the first group has the effect of changing both without decoupling). As usual, the newpar command is used: XSPEC12>newpar 10 1 Fit statistic : Chi-Squared = 2.025975e+07 using 66 PHA bins. Test statistic : Chi-Squared = 2.025975e+07 using 66 PHA bins. Reduced chi-squared = 343385.7 for 59 degrees of freedom Null hypothesis probability = 0.000000e+00 Current data and model not fit yet. XSPEC12>fit ... ======================================================================== Model phabs<1>*highecut<2>*powerlaw<3> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp Data group: 1 1 1 phabs nH 10^22 20.1548 +/- 0.181919 2 2 highecut cutoffE keV 14.6847 +/- 5.59260E-02 3 2 highecut foldE keV 7.41661 +/- 8.99590E-02 4 3 powerlaw PhoIndex 1.18693 +/- 6.33041E-03 5 3 powerlaw norm 5.88350E-02 +/- 9.30517E-04 Data group: 2 6 1 phabs nH 10^22 1.27017 +/- 3.77710E-02 7 2 highecut cutoffE keV 14.6847 = p2 8 2 highecut foldE keV 7.41661 = p3 9 3 powerlaw PhoIndex 1.18693 = p4 10 3 powerlaw norm 0.312138 +/- 4.49104E-03 ________________________________________________________________________ Fit statistic : Chi-Squared = 15423.79 using 66 PHA bins. After fitting, this decoupling reduces χ2 by a factor of six to 15,478, but this is still too high. Indeed, this simple attempt to account for the spectral variability in terms of “blanket” cold absorption and scattering does not work. More sophisticated models, involving additional components and partial absorption, should be tried. 62 CHAPTER 4. WALKS THROUGH XSPEC 4.5 Multiple Models: a Background Modeling Example In the previous section we showed how to fit the same model to multiple datasets. We now demonstrate how to fit multiple models, each with their own response, to the same dataset. There are several reasons why this may be useful, for instance: • We are using data from a coded aperture mask. If there are multiple sources in the field they will all contribute to the spectrum from each detector. However, each source may have a different response due to its position. • it We are observing an extended source using a telescope whose PSF is large enough that the signal from different regions are mixed together. In this case we will want to analyze spectra from all regions of the source simultaneously with each spectrum having a contribution from the model in other regions. • We wish to model the background spectrum that includes a particle component. The particle background will have a different response from the X-ray background because the particles come from all directions, not just down the telescope. We will demonstrate the third example here. Suppose we have a model for the background spectrum that requires a different response to that for the source spectrum. Read in the source and background spectra as separate files: XSPEC12>data 1:1 source.pha 2:2 back.pha The source and background files have their own response matrices: XSPEC12>response 1 source.rsp 2 back.rsp Set up the model for the source. Here we will take the simple case of an absorbed power-law: XSPEC12>model phabs(pow) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001( 0.01) 0 0 100000 1e+06 1:data group 1::phabs:nH> 4.5. MULTIPLE MODELS: A BACKGROUND MODELING EXAMPLE 63 1 0.01( 0.01) 2:data group 1::powerlaw:PhoIndex> 1 0.01( 0.01) 3:data group 1::powerlaw:norm> -3 -2 9 10 0 0 1e+24 1e+24 Input parameter value, delta, min, bot, top, and max values for ... 1 0.001( 0.01) 0 0 100000 4:data group 2::phabs:nH> 1 0.01( 0.01) -3 -2 9 5:data group 2::powerlaw:PhoIndex> 1 0.01( 0.01) 0 0 1e+24 6:data group 2::powerlaw:norm>0 0 1e+06 10 1e+24 Note that we have fixed the normalization of the source model for the background dataset at zero so it doesn’t contribute. Now we need to set up the background model for both datasets with the appropriate response matrices. XSPEC12>response 2:1 source.rsp 2:2 back.rsp This tells XSPEC that both these datasets have a second model which must be multiplied by the source.rsp response matrix for its contribution to the source region and back.rsp for its contribution to the background region. This is likely to be the case for a standard imaging data where the response will only depend on the extraction region. Note that for a coded-aperture mask the situation may be more complicated with a different response for the source and the background even if they are extracted from the same region of the detector. We now define the background model to be used. In this case take the simple example of a single power-law XSPEC12>model 2:myback pow Input parameter value, delta, min, bot, top, and max values for ... 1 0.01( 0.01) -3 -2 9 1:myback:data group 1::powerlaw:PhoIndex> 1 0.01( 0.01) 0 0 1e+24 2:myback:data group 1::powerlaw:norm> Input parameter value, delta, min, bot, top, and max values for ... 1 0.01( 0.01) -3 -2 9 3:myback:data group 2::powerlaw:PhoIndex> 1 0.01( 0.01) 0 0 1e+24 4:myback:data group 2::powerlaw:norm> We have now set up XSPEC so that the source data is compared to a source model multiplied by the source response plus a background model multiplied 10 1e+24 10 1e+24 64 CHAPTER 4. WALKS THROUGH XSPEC by the background response and the background data is compared to the background model multiplied by the background response. The background models fitted to the source and background data are constrained to be the same. XSPEC12>show param Parameters defined: ======================================================================== Model phabs<1>*powerlaw<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp Data group: 1 1 1 phabs nH 10^22 1.00000 +/- 0.0 2 2 powerlaw PhoIndex 1.00000 +/- 0.0 3 2 powerlaw norm 1.00000 +/- 0.0 Data group: 2 4 1 phabs nH 10^22 1.00000 = 1 5 2 powerlaw PhoIndex 1.00000 = 2 6 2 powerlaw norm 0.00000 frozen ________________________________________________________________________ ======================================================================== Model myback:powerlaw<1> Source No.: 2 Active/On Model Model Component Parameter Unit Value par comp Data group: 1 1 1 powerlaw PhoIndex 1.00000 +/- 0.0 2 1 powerlaw norm 1.00000 +/- 0.0 Data group: 2 3 1 powerlaw PhoIndex 1.00000 = myback:1 4 1 powerlaw norm 1.00000 = myback:2 ________________________________________________________________________ It is often the case that the response information is split into an RMF and ARF, where the RMF describes the instrument response and the ARF the telescope effective area. The particle background can then be included by using the RMF but not the ARF: XSPEC12>data 1:1 source.pha 2:2 back.pha XSPEC12>response 1 source.rmf 2 source.rmf XSPEC12>arf 1 source.arf XSPEC12>response 2:1 source.rmf 2:2 source.rmf 4.6. USING XSPEC TO SIMULATE DATA: AN EXAMPLE FOR CHANDRA65 4.6 Using XSPEC to Simulate Data: an Example for Chandra In several cases, analyzing simulated data is a powerful tool to demonstrate feasibility. For example: • To support an observing proposal. That is, to demonstrate what constraints a proposed observation would yield. • To support a hardware proposal. If a response matrix is generated, it can be used to demonstrate what kind of science could be done with a new instrument. • To support a theoretical paper. A theorist could write a paper describing a model, and then show how these model spectra would appear when observed. This, of course, is very like the first case. Here, we will illustrate the first example. The first step is to define a model on which to base the simulation. The way XSPEC creates simulated data is to take the current model, convolve it with the current response matrix, while adding noise appropriate to the integration time specified. Once created, the simulated data can be analyzed in the same way as real data to derive confidence limits. Let us suppose that we want to measure the intrinsic absorption of a faint highredshift source using Chandra. Our model is thus a power-law absorbed both by the local Galactic column and an intrinsic column. First, we set up the model. From the literature we have the Galactic absorption column and redshift so: XSPEC12>mo pha*zpha(zpo) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001( 0.01) 0 0 100000 1:phabs:nH>0.08 1 0.001( 0.01) 0 0 100000 2:zphabs:nH>1.0 0 -0.01( 0.01) -0.999 -0.999 10 3:zphabs:Redshift>5.1 1 0.01( 0.01) -3 -2 9 4:zpowerlw:PhoIndex>1.7 0 -0.01( 0.01) -0.999 -0.999 10 5:zpowerlw:Redshift>5.1 1 0.01( 0.01) 0 0 1e+24 6:zpowerlw:norm> ======================================================================== 1e+06 1e+06 10 10 10 1e+24 66 CHAPTER 4. WALKS THROUGH XSPEC Model phabs<1>*zphabs<2>*zpowerlw<3> Source No.: 1 Active/Off Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 8.00000E-02 +/- 0.0 2 2 zphabs nH 10^22 1.00000 +/- 0.0 3 2 zphabs Redshift 5.10000 frozen 4 3 zpowerlw PhoIndex 1.70000 +/- 0.0 5 3 zpowerlw Redshift 5.10000 frozen 6 3 zpowerlw norm 1.00000 +/- 0.0 Now suppose that we know that the observed 0.5–2.5 keV flux is 1.1 × 10−13 ergs/cm2 /s. We now can derive the correct normalization by using the commands energies, flux and newpar. That is, we’ll determine the flux of the model with the normalization of unity. We then work out the new normalization and reset it: XSPEC12>energies 0.5 2.5 1000 XSPEC12>flux 0.5 2.5 Model Flux 0.052736 photons (1.0017e-10 ergs/cm^2/s) range (0.50000 - 2.5000 keV) XSPEC12> newpar 6 1.1e-3 3 variable fit parameters XSPEC12>flux Model Flux 2.6368e-05 photons (5.0086e-14 ergs/cm^2/s) range (0.50000 - 2.5000 keV) Here, we have changed the value of the normalization (the fifth parameter) from 1.0 to 1.1 × 10−13 / 1.00 × 10−10 = 1.1 × 10−3 to give the required flux. The simulation is initiated with the command fakeit. If the argument none is given, we will be prompted for the name of the response matrix. If no argument is given, the current response will be used. We also need to reset the energies command before the fakeit to ensure that the model is calculated on the entire energy range of the response: XSPEC12>energies reset XSPEC12>fakeit none For fake data, file #1 needs response file: aciss_aimpt_cy15.rmf ... and ancillary response file: aciss_aimpt_cy15.arf There then follows a series of prompts asking the user to specify whether he or she wants counting statistics (yes!), the name of the fake data (file test.fak in our example), and the integration time (40,000 seconds – the correction norm and background exposure time can be left at their default values). Use counting statistics in creating fake data? (y): Input optional fake file prefix: 4.6. USING XSPEC TO SIMULATE DATA: AN EXAMPLE FOR CHANDRA67 Fake data file name (aciss_aimpt_cy15.fak): test.fak Exposure time, correction norm, bkg exposure time (1.00000, 1.00000, 1.00000): 40000.0 No background will be applied to fake spectrum #1 1 spectrum in use Fit statistic : Chi-Squared = 350.95 using 1024 PHA bins. ***Warning: Chi-square may not be valid due to bins with zero variance in spectrum number(s): 1 Test statistic : Chi-Squared = 350.95 using 1024 PHA bins. Reduced chi-squared = 0.34407 for 1020 degrees of freedom Null hypothesis probability = 1.000000e+00 ***Warning: Chi-square may not be valid due to bins with zero variance in spectrum number(s): 1 Current data and model not fit yet. The first thing we should note is that the default statistics are not correct for these data. For Poisson data and no background we should cstat for the fit statistic and pchi for the test statistic: XSPEC12>statistic cstat Default fit statistic is set to: C-Statistic This will apply to all current and newly loaded spectra. Fit statistic : C-Statistic = 513.63 using 1024 PHA bins and 1020 degrees of freedom. Test statistic : Chi-Squared = 350.95 using 1024 PHA bins. Reduced chi-squared = 0.34407 for 1020 degrees of freedom Null hypothesis probability = 1.000000e+00 ***Warning: Chi-square may not be valid due to bins with zero variance in spectrum number(s): 1 Current data and model not fit yet. XSPEC12>statistic test pchi Default test statistic is set to: Pearson Chi-Squared This will apply to all current and newly loaded spectra. Fit statistic : C-Statistic = 513.63 using 1024 PHA bins and 1020 degrees of freedom. Test statistic : Pearson Chi-Squared = Reduced chi-squared = 0.62682 for 639.35 using 1024 PHA bins. 1020 degrees of freedom 68 Null hypothesis probability = CHAPTER 4. WALKS THROUGH XSPEC 1.000000e+00 ***Warning: Pearson Chi-square may not be valid due to bins with zero model value in spectrum number(s): 1 Current data and model not fit yet. As we can see from the warning message we need to ignore some channels where there is no effective response. Looking at a plot of the data and model indicates we should ignore below 0.15 keV and above 10 keV so: XSPEC12>ignore **-0.15 10.0-** 11 channels (1-11) ignored in spectrum # 1 340 channels (685-1024) ignored in spectrum # Fit statistic : C-Statistic = 1 510.55 using 673 PHA bins and 669 degrees of freedom. Test statistic : Pearson Chi-Squared = 635.19 using 673 PHA bins. Reduced chi-squared = 0.94947 for 669 degrees of freedom Null hypothesis probability = 8.217205e-01 Current data and model not fit yet. We assume that the Galactic column is known so freeze the first parameter. We then perform a fit followed by the error command: XSPEC12> freeze 1 ... XSPEC12>fit ... XSPEC12>parallel error 3 XSPEC12>err 2 4 6 Parameter Confidence Range ( 2.706) 2 1.16302 5.64805 (-2.00255,2.48247) 4 1.73345 1.95111 (-0.106137,0.111521) 6 0.00126229 0.00221906 (-0.000397759,0.000559019) Note that our input parameters do not lie within the 90% confidence errors however since this will happen one times in ten (by definition) this should not worry us unduly. For a real observing proposal we would likely repeat this experiment with different input values of the intrinsic absorption to learn how well we could constrain it given a range of possible actual values. 4.7. PRODUCING PLOTS: MODIFYING THE DEFAULTS 4.7 69 Producing Plots: Modifying the Defaults The final results of using XSPEC are usually one or more tables containing confidence ranges and fit statistics, and one or more plots showing the fits themselves. So far, the plots shown have generally used the default settings, but it is possible to edit plots to improve their appearance. The plotting package used by XSPEC is PGPLOT, which is comprised of a library of low-level tasks. At a higher level is QDP/PLT, the interactive program that forms the interface between the XSPEC user and PGPLOT. QDP/PLT has its own manual; it also comes with on-line help. Here, we show how to make some of the most common modifications to plots. In this example, we’ll take the simulated Chandra spectrum and make a better plot. Figure 4.9 shows the basic data and folded model plot. The only additional changes we have made to this plot are to increase the line widths to make them print better. We made this plot as follows: XSPEC12>setplot energy XSPEC12>iplot data PLT>lwidth 3 PLT>lwidth 3 on 1..2 PLT>time off PLT>hard figi.png/ps The first lwidth command increases the line widths on the frame while the second increases it on the data and model. The time off command just removes a username and time stamp from the bottom right of the plot. The hard command makes a “hardcopy”, in this case a PostScript file. Before looking at other PLT commands we can use to modify the plot there is one more thing we can try at the XSPEC level. The current bin sizes are much smaller than the resolution so we may as well combine bins in the plot (but not in the fitting) to make it clearer. XSPEC12>setplot rebin 100 4 Combines four spectral bins into one. The setplot rebin command combines bins till either a S/N of the first argument is reached or the number of bins in the second argument have been combined. We do an iplot again then do the following modifications: PLT> viewport 0.2 0.2 0.8 0.7 70 CHAPTER 4. WALKS THROUGH XSPEC Figure 4.9: The data and folded model for the simulated Chandra ACIS-S spectrum. The first thing we’ll do is change the aspect ratio of the box that contains the plot (viewport in QDP terminology). The viewport is defined as the coordinates of the lower left and upper right corners of the page. The units are such that the full width and height of the page are unity. The labels fall outside the viewport, so if the full viewport were specified, only the plot would appear. The default box has a viewport with corners at (0.1, 0.1) and (0.9, 0.9). For our purposes, we want a viewport with corners at (0.2, 0.2) and (0.8, 0.7): with this size and shape, the hardcopy will fit nicely on the page and not have to be reduced for photocopying. Next we want to change some of the labels: PLT> label top Simulated Spectrum PLT> label file Chandra ACIS-S PLT> label y counts s\u-1\d keV\u-1\d Note the change in the y-axis label is to remove the string “normalized”. The normalization referred to is almost always unity so this label can generally be changed. To get help on a PLT command, just type help followed by the 4.7. PRODUCING PLOTS: MODIFYING THE DEFAULTS 71 Figure 4.10: A simulated Chandra ACIS-S spectrum produced to show how a plot can be modified by the user. name of the command. Note that PLT commands can be abbreviated, just like XSPEC commands. To see the results of changing the viewport and the labels, just enter the command plot. The two changes we want to make next are to rescale the axes and to change the y-axis to a logarithmic scale. The commands for these changes also are straightforward: the rescale command takes the minimum and maximum values as its arguments, while the log command takes x or y as arguments: PLT> PLT> PLT> PLT> rescale x 0.3 6.0 rescale y 1.0e-4 0.03 log y plot To revert to a linear scale, use the command log off y. The only remaining extra change is to reduce the size of the characters. This is done using the csize command, which takes the normalization as its argument. One (1) will not change the size, a number less than one will reduce it and a number bigger than one will increase it. 72 CHAPTER 4. WALKS THROUGH XSPEC PLT> csize 0.8 PLT> plot We make the PostScript file and also save the plot information using the wenv command that, in this case, writes files figj.qdp and figj.pco containing the plot data and commands, respectively. PLT> hardcopy figj.png/ps PLT> wenv figj PLT> quit The result of all this manipulation is shown proudly in Figure 4.10. 4.8 Markov Chain Monte Carlo Example To illustrate MCMC methods we will use the same data as the first walkthrough. XSPEC12> ... XSPEC12> ... XSPEC12> ... XSPEC12> XSPEC12> XSPEC12> data s54405 model phabs(pow) fit chain type gw chain walkers 8 chain length 10000 We use the Goodman-Weare algorithm with 8 walkers and a total length of 10,000. For the G-W algorithm the actual number of steps are 10,000/8 but we combine the results from all 8 walkers into a single output file. We start the chain at the default model parameters except that we use the renorm command to make sure that the model and the data have the same normalization. If we had multiple additive parameters with their own norms then a good starting place would be to use the fit 1 command to initially set the normalizations to something sensible. XSPEC12> chain run test1.fits The first thing to check is what happened to the fit statistic during the run. XSPEC12> plot chain 0 4.8. MARKOV CHAIN MONTE CARLO EXAMPLE 73 Figure 4.11: The statistic from an MCMC run showing the initial burn-in phase. 74 CHAPTER 4. WALKS THROUGH XSPEC The result is shown in Figure 4.11, which plots the statistic value against the chain step. It is clear that after about 2000 steps the chain reached a steady state. We would usually have told XSPEC to discard the first few thousand steps but included them for illustrative purposes. Let us do this again but specifying a burn-in phase that will not be stored. XSPEC12> chain burn 5000 XSPEC12> chain run test1.fits The output chain now comprises 10,000 steady-state samples of the parameter probability distribution. Repeating plot chain 0 will confirm that the chain is in a steady state. The other parameter values can be plotted either singly using eg plot chain 2 for the power-law index or in pairs eg plot chain 1 2 giving a scatter plot as shown in Figure 4.12. Using the error command at this point will generate errors based on the chain values. XSPEC12>error 1 2 3 Errors calculated from chains Parameter Confidence Range (2.706) 1 0.264971 0.919546 2 2.1134 2.41307 3 0.0107304 0.0171814 The 90% confidence ranges are determined by ordering the parameter values in the chain then finding the center 90%. To make confidence regions in two dimensions the margin command takes the same arguments as steppar. XSPEC12>margin 1 0.0 1.5 25 2 1.5 3.0 25 Then use plot integprob to produce a plot of confidence regions. 4.9 INTEGRAL/SPI: A Walk Through Example Consider an observation of the Crab, for which a (standard) 5 × 5 dithering observation strategy was employed. Since the Crab (pulsar and nebular components are of course un-resolvable at INTEGRAL’s spatial resolution) is by 4.9. INTEGRAL/SPI: A WALK THROUGH EXAMPLE Figure 4.12: The scatter plot from a 10,000 step MCMC run. 75 76 CHAPTER 4. WALKS THROUGH XSPEC far the brightest source in it immediate region of the sky, and its position is precisely known, we can opt not to perform SPI or IBIS imaging analysis prior to XSPEC analysis. We thus run the standard INTEGRAL/SPI analysis chain on detectors 0-18 up to the SPIHIST level for (or BIN I level in the terminology of the INTEGRAL documentation), selecting the “PHA” output option. We then run SPIARF, providing the name of the PHA-II file just created, and selecting the “update” option in the spiarf.par parameter file (you should refer to the SPIARF documentation prior to this step if it is unfamiliar). The celestial coordinates for the Crab are provided in decimal degrees (RA,Dec = 83.63,22.01) interactively or by editing the parameter file. This may take a few minutes, depending on the speed of your computer and the length of your observation. Once completed, SPIARF must be run one more time, setting the “bkg resp” option to “y”; this creates the response matrices to be applied to the background model, and updates the PHA-II response database table accordingly. Then SPIRMF, which interpolates the template RMFs to the users desired spectral binning, also writes information to the PHA response database table to be used by XSPEC. Finally, you should run SPIBKG INIT, which will construct a set of bbackground spectral templates to initialize the SPI background model currently installed in XSPEC (read the FTOOLS help for that utility carefully your first time). You are now ready to run XSPEC; a sample session might look like this (some repetitive output has been suppressed): % % xspec XSPEC version: 12.2.1 Build Date/Time: Wed Nov 2 17:14:21 2005 XSPEC12>package require Integral 1.0 1.0 XSPEC12>data ./myDataDir/rev0044_crab.pha{1-19} 19 spectra in use RMF # 1 Using Response (RMF) File RMF # 2 Using Response (RMF) File RMF # 3 Using Response (RMF) File resp/comp1_100x100.rmf resp/comp2_100x100.rmf resp/comp3_100x100.rmf Using Multiple Sources For Source # 1 Using Auxiliary Response (ARF) Files resp/rev0044_100ch_crab_cmp1.arf.fits 4.9. INTEGRAL/SPI: A WALK THROUGH EXAMPLE 77 resp/rev0044_100ch_crab_cmp2.arf.fits resp/rev0044_100ch_crab_cmp3.arf.fits For Source # 2 Using Auxiliary Response (ARF) Files resp/rev0044_100ch_bkg_cmp1.arf.fits resp/rev0044_100ch_bkg_cmp2.arf.fits resp/rev0044_100ch_bkg_cmp3.arf.fits Source File: ./myDataDir/rev0044_crab.pha{1} Net count rate (cts/s) for Spectrum No. 1 3.7011e+01 Assigned to Data Group No. : 1 Assigned to Plot Group No. : 1 Source File: ./myDataDir/rev0044_crab.pha{2} Net count rate (cts/s) for Spectrum No. 2 3.7309e+01 Assigned to Data Group No. : 1 Assigned to Plot Group No. : 2 ... Source File: ./myDataDir/rev0044_crab.pha{19} Net count rate (cts/s) for Spectrum No. 19 3.6913e+01 Assigned to Data Group No. : 1 Assigned to Plot Group No. : 19 +/- 1.2119e-01 +/- 1.2167e-01 +/- 1.2103e-01 XSPEC12>mo 1:crab po Input parameter value, delta, min, bot, top, and max values for ... 1 PhoIndex 1.0000E+00 1.0000E-02 -3.0000E+00 crab::powerlaw:PhoIndex>2.11 0.01 1.5 1.6 2.5 2.6 2 norm 1.0000E+00 1.0000E-02 0.0000E+00 crab::powerlaw:norm>8. 0.1 1. 2. 18. 20. ... XSPEC12>mo 2:bkg spibkg5 -2.0000E+00 9.0000E+00 0.0000E+00 1.0000E+24 Input parameter value, delta, min, bot, top, and max values for ... 1 Par_1 0.0000E+00 1.0000E-02 -2.0000E-01 -1.5000E-01 1.5000E-01 bkg::spibkg5:Par_1>/* ... _____________________________________________________________________________________________________ XSPEC12>ign 1-19:68-80 ... XSPEC12>ign 1-19:90-100 ... XSPEC12>fit Number of trials and critical delta: 10 1.0000000E-02 ... ======================================================================== Model bkg:spibkg5 Source No.: 2 Active/On Model Component Name: spibkg5 Number: 1 78 CHAPTER 4. WALKS THROUGH XSPEC N Name Unit Value Sigma 1 Par_1 9.0650E-03 +/2.8651E-03 2 Par_2 1.6174E-02 +/3.4778E-03 ... 25 Par_25 -1.9537E-02 +/6.1429E-03 26 norm 9.7286E-01 +/1.3527E-03 ________________________________________________________________________ ======================================================================== Model crab:powerlaw Source No.: 1 Active/On Model Component Name: powerlaw Number: 1 N Name Unit Value Sigma 1 PhoIndex 2.1163E+00 +/1.8946E-02 2 norm 1.1390E+01 +/8.1414E-01 ________________________________________________________________________ Chi-Squared = 1.8993005E+03 using 1463 PHA bins. Reduced chi-squared = 1.3235544E+00 for 1435 degrees of freedom Null hypothesis probability = 1.5268098E-15 XSPEC12> Note that the syntax used for the data statement (appended curly bracket, specifying use of spectra 1-19), and the separate model commands, which are indexed and named (in this case simply “crab” for the source of interest and “bkg” for the background model, “spibkg lo”. These commands are described in detail elsewhere in this document, as are the the spibkg lo, spibkg med and spibkg hi models. In this case, 100 logarithmically-spaced energy bins spanning the nominal 20-8000 keV band of the SPI instrument were used. In this example, only one dither-point was used to solve for the Crab spectrum, and the background. The simple assumption of a single background spectrum (i.e. no detector-to-detector variations) was assumed. In general, and particularly for fainter sources, a much larger number of spectra will be needed for a solution (and even for the Crab, the quality of the fit, and the accuracy of the inferred parameters can be improved). Also, detector-to-detector and/or time (i.e. pointing-to-pointing) variations will need to be considered. This can be accomplished using the data-grouping feature of XSPEC, which will be described subsequently. Also notice that channels between about 70 and 80 were ignored; this is because there are detector electronic effects contaminating the single-event data for energies from 1250-1400 keV (refer to the SPI data analysis manual for additional discussion), and that there are a lot of (scientifically uninteresting) background model parameters. Also, the highest energies were ignored, since the source flux becomes insignificant relative to the background. Some results are illustrated below. These plots were generated with the sequence of commands: 4.9. INTEGRAL/SPI: A WALK THROUGH EXAMPLE 79 XSPEC12> setplot group 1-19 XSPEC12> plot ldata res ... XSPEC12> plot ufspec Note that without the “setplot group” command, XSPEC would plot 19 sets of spectral data, models and residuals. The can become confusing, especially as the number of spectra included in an analysis becomes much larger than 19! On the other hand, it can be useful to divide the data into subsets for plotting purposes, e.g. setplot group 1-6 7-12 13-19, to get an idea of relative shadowing effects of the coded-mask. The left hand plot illustrates the source model, the background model, the total model (i.e. source + background), and the data (here in count rates per channel). The right hand plot illustrates the “unfolded model” (blue, power-law curve), the summed model, and the data as a photon flux. A possible source of confusion is the similarity of the background model curves plotted in theses two separate representations. The explanation is that the background, which is dominated by instrumental contributions, is modeled in detector count space (i.e. the background response matrix has unit effective area. Thus, to be strictly correct, the right-hand plot is a hybrid of the photon source model and the detector-rate background model. We further note that at the present time, XSPEC does not have the capability to plot (or store and manipulate) the background subtracted data. This is a feature under consideration for a future release. If we had chosen a observation containing more than a single source, the procedure would have been similar, except that the sequence of model commands would be extended, e.g. XSPEC12>data ./MyDataDir/GCDE_aug_03.pha{1-475} ... XSPEC12> model 1:1e1740 po ... XSPEC12> model 2:gx1_4 po ... XSPEC12> model 3:bkg spibkg_lo ... Here data from the Galactic Center deep exposure campaign were loaded, and two sources are sought. In this case, a much larger number of spectra were loaded (475 spectra corresponds to one full 5 × 5 dither using all 19 detectors. In this case, the simple approach of applying constant background (i.e. no detector-to-detector or pointing-to-pointing variation) to the full data set is 80 CHAPTER 4. WALKS THROUGH XSPEC likely to be a poor approximation. A more realistic approach would be to use the XSPEC grouping capability to handle such variations in the background solution. This can be accomplished in the usual manner (refer to the description of the grouping command in this document), however, it can become tedious in terms of the required command line inputs. For example, to establish a separate data group for each detector for a long (e.g. 5 × 5 dither) observations, a sequence of commands such as this would be required: XSPEC12> XSPEC12> XSPEC12> ... XSPEC12> XSPEC12> XSPEC12> XSPEC12> ... XSPEC12> XSPEC12> XSPEC12> XSPEC12> ... XSPEC12> XSPEC12> data data data 1:1 2:2 3:3 ./MyDataDir/rev0044_Crab.pha.fits{1} ./MyDataDir/rev0044_Crab.pha.fits{2} ./MyDataDir/rev0044_Crab.pha.fits{3} data data data data 19:19 1:20 2:21 3:22 ./MyDataDir/rev0044_Crab.pha.fits{19} ./MyDataDir/rev0044_Crab.pha.fits{20} ./MyDataDir/rev0044_Crab.pha.fits{21} ./MyDataDir/rev0044_Crab.pha.fits{22} data data data data 19:38 1:39 2:40 3:41 ./MyDataDir/rev0044_Crab.pha.fits{38} /MyDataDir/rev0044_Crab.pha.fits{39} ./MyDataDir/rev0044_Crab.pha.fits{40} ./MyDataDir/rev0044_Crab.pha.fits{41} data data 18:474 19:475 ./MyDataDir/rev0044_Crab.pha.fits{474} ./MyDataDir/rev0044_Crab.pha.fits{475} One might then for example, make a first cut attempt by fitting a constant background. Then, as a next step, one might allow the normalization terms of the background model to vary over the groups (i.e. over the detector plane). This is accomplished with the “untie” command, using the following sequence: XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> XSPEC12> untie untie untie untie untie untie untie untie untie untie untie untie untie untie untie untie untie untie bkg:52 bkg:78 bkg:104 bkg:130 bkg:156 bkg:182 bkg:208 bkg:234 bkg:260 bkg:286 bkg:312 bkg:338 bkg:364 bkg:390 bkg:416 bkg:442 bkg:468 bkg:487 4.10. INTEGRAL SPECIFIC COMMAND LINE SCRIPTS 81 Note that use of the “bkg” identifier, which associates the parameters index with the background model. The specific sequence of numbers use here requires some explanation; the particular background model employed has 25 parameters (which simply correspond in rank order to the 25 most variable individual bins), and a normalization term, i.e. parameter 26. Thus, the normalization for the second detector group is parameter 52, for the third parameter 78, and so on. Similar command sequences can be used to untie additional background model parameters. Supposing that we did this and refitted the data. We then might, for example wish to go back and freeze the individual normalization terms with the freeze command: XSPEC12> XSPEC12> ... XSPEC12> freeze freeze bkg:26 bkg:52 freeze bkg:487 By now though, you probably get the idea that this all requires an unreasonable amount of command-line input. To circumvent this problem, a number of INTEGRAL/SPI specific “tcl” scripts are available which greatly streamline this process. 4.10 INTEGRAL Specific Command Line Scripts SPIdata The SPIdata procedure, which when installed can be treated as an XSPEC command, greatly facilitates the data initialization step. For example, the command XSPEC12> SPIdata ./MyData/Dir/rev0044_crab.pha 475 det Y Opens the Crab observation spectral data file, reads the 475 spectra into memory, grouping them by detector. The “Y” then indicates that, yes, I wish to ignore the spectral data channels corresponding to the known detector-electronic noise contamination (this is the default). Instead of “det” as the grouping option I could have selected “time” to group by time (quantized into dither-pointing intervals). A “-” lead to the data being initialzed into a single group. The command: XSPEC12> SPIdata ./MyData/Dir/rev0044_crab.pha 475 Reads the 475 spectra into a single data group, and ignores the undesirable channels. If you forget all this, the command 82 XSPEC12> CHAPTER 4. WALKS THROUGH XSPEC SPIdata -h will remind you. The scripts SPIuntie, and SPIfreeze have similar command-line syntax. SPIuntie and SPIfreeze XSPEC12> SPIuntie bkg 475 19 -1 The SPIuntie command script will accomplish the same result as the sequence of “untie” commands in the INTEGRAL/SPI example presented in this document. In that case, we had loaded 475 spectra associated with a single 5×5-dither pattern centered on the Crab nebula. The spectra were grouped by detector, which is a common approach to SPI analysis given the known detector-to-detector variations in the background rates. Suppose after an initial fitting pass, for which we assumed a single background spectrum, we know wish to untie the individual data group (i.e. detector) background models. This can be accomplished by issuing 25 “untie” commands as previously noted, or in a single command line using the SPIuntie command: XSPEC12> SPIuntie bkg 475 19 -1 untie bkg:52 Chi-Squared = 1.2030200E+04 using 1615 PHA bins. Reduced chi-squared = 7.5852458E+00 for 1586 degrees of freedom Null hypothesis probability = 0.0000000E+00 untie bkg:78 Chi-Squared = 1.2030200E+04 using 1615 PHA bins. Reduced chi-squared = 7.5900314E+00 for 1585 degrees of freedom Null hypothesis probability = 0.0000000E+00 untie bkg:104 renorm: no renormalization necessary Chi-Squared = 1.2030200E+04 using 1615 PHA bins. Reduced chi-squared = 7.5948231E+00 for 1584 degrees of freedom Null hypothesis probability = 0.0000000E+00 ... One might then make a second pass at fitting the data, hopefully leading to improved statistics. Subsequently, additional background model parameters could be untied using the SPIuntie procedure as well. For example, to untie three additional parameters over the full data set, the command syntax is: 4.10. INTEGRAL SPECIFIC COMMAND LINE SCRIPTS 83 XSPEC12> SPIuntie bkg 475 19 1 3 ... This will untie the first 3 parameters of the background model identified by “bkg”, i.e. equivalent to issuing (475 − 1) × 3 individual untie commands. Note that you can always be reminded of the command-line argument definitions by typing “SPIuntie -h” at the XSPEC prompt. Suppose now that you are satisfied with the relative background normalization terms, and wish to freeze them at their current values for subsequent fitting passes. This could be accomplished using the SPIfreeze command script: XSPEC12> SPIfreeze bkg 475 -1 XSPEC12>SPIfreeze bkg 19 1 -1 freeze bkg:52 1 Chi-Squared = 6.6232600E+05 using 1805 PHA bins. Reduced chi-squared = 3.7589444E+02 for 1762 degrees of freedom Null hypothesis probability = 0.0000000E+00 freeze bkg:78 Chi-Squared = 6.5791894E+05 using 1805 PHA bins. Reduced chi-squared = 3.7318148E+02 for 1763 degrees of freedom Null hypothesis probability = 0.0000000E+00 ... As with the SPIuntie command script, typing “SPIfreeze -h” at the XSPEC prompt will scroll the command-line definitions to your screen. Note that the current SPI background models, which are documented elsewhere, are designed so that the parameter list is hierarchically ordered in terms of decreasing criticality. Thus, freeing the first parameter is likely to have the most significant impact on the statistics, the second parameter, the next most significant, and so on. 84 CHAPTER 4. WALKS THROUGH XSPEC Chapter 5 XSPEC Commands 5.1 Summary of Commands The following is a list of the commands available in XSPEC, together with a brief description of the purpose of each. The commands have been categorized under six headings: Control, Data, Fit, Model, Plot, Script, and Setting. The Control commands contain the interface with the operating system: they cause commands to be executed, or user input written to disk, or control how much is output. The Data commands manipulate the data being analyzed, by reading data into the program or replacing spectra or their ancillary detector, background, correction, or efficiency (auxiliary response) arrays. Additionally data commands control the channels under analysis. The fit commands invoke the fitting routines, modify their behavior by interchanging fitting algorithms or statistics in use, fixing parameters, or perform statistical testing. The Model commands create or manipulate the model, adding or editing components, modifying parameters, or alternatively performing analytical calculations from a model. The Plot commands deal with all aspects of plotting. The scripts are auto-loaded Tcl scripts that can be used in the same ways as commands. Finally the Setting commands sets variables that affect theoretical models. Command Category Description abund SETTING Set the abundance table. addcomp MODEL Add a component to the model. addline MODEL Add lines to a model arf DATA Read an auxiliary response file. 85 86 CHAPTER 5. XSPEC COMMANDS autosave CONTROL Periodically save the XSPEC status. backgrnd DATA Reset the files to be used for background subtraction. bayes FIT Set up for Bayesian inference. chain FIT Run a Monte Carlo Markov Chain. chatter CONTROL Control the verbosity of XSPEC. corfile DATA Reset the files to be used for background correction. cornorm DATA Reset the normalization to be used in correcting the background. cosmo SETTING Set H0 , q0 , and Λ0 cpd PLOT Alias for setplot device. data DATA Input one or more PHA data files. delcomp MODEL Delete a component from the model. diagrsp DATA Diagonalize the current response for an ideal response. dummyrsp MODEL Create a dummy response, covering a given energy range. editmod MODEL Add, delete, or replace one component in the model. energies MODEL Specify new energy binning for model fluxes. eqwidth MODEL Calculate a model component2̆019s equivalent width. error (rerror) FIT Determine a single parameter confidence region. rerror is for response parameters. exec CONTROL Execute a shell command from within XSPEC. exit CONTROL Wind up any hardcopy plots and exit from XSPEC. extend MODEL This is now obsolete. See energies command. fakeit DATA Produce simulated data files for sensitivity studies. fit FIT Find the best fit model parameters. flux MODEL Calculate the current model’s flux over an energy range. freeze (rfreeze) FIT Do not allow a model parameter to vary during the fit. rfreeze is for response parameters. ftest FIT Calculate the F-statistic between two model fits gain MODEL Perform a simple modification of the response gain. goodness FIT Monte Carlo calculation of goodness-of-fit. 5.1. SUMMARY OF COMMANDS 87 hardcopy PLOT Spool the current plot to the printer. help CONTROL Obtain help on XSPEC commands. identify MODEL List possible lines in the specified energy range. ignore DATA Ignore a range of PHA channels in future fit operations. initpackage MODEL Compile, build, and initialize a package of local models. iplot PLOT As plot command but interactive using PLT. lmod MODEL Load a package of local models. log CONTROL Open the log file to save output. lrt SCRIPT Likelihood ratio test between two models. lumin MODEL Calculate the current model’s luminosity over a given rest frame energy range and redshift. margin FIT MCMC probability distribution. mdefine MODEL Define a simple expression. method SETTING Set the minimization method. model (rmodel) MODEL Define the model to be used when fitting the data. modid MODEL Guess line IDs in the model. multifake SCRIPT Perform many iterations of fakeit and save the results in a FITS file. newpar (rnewpar) MODEL Modify the model parameters. rnewpar is for response parameters notice DATA Restore a range of PHA channels for future operations. parallel CONTROL Enable parallel processing for particular tasks in XSPEC. plot PLOT Plot various information on the current plot device. query CONTROL Switch on/off prompt to continue fitting. quit CONTROL An alias for exit renorm FIT Adjust the model norms, and/or control automatic renorming. rescalecov SCRIPT Rescale the covariance matrix used in the proposal chain command. model using an arithmetic 88 CHAPTER 5. XSPEC COMMANDS response DATA Reset the files used to determine the detector responses. save CONTROL Save aspects of the current state to a command file. script CONTROL Open the script file to save all commands input. setplot PLOT Modify the plot device and other values used by the plot routines. show CONTROL Display current file and model information. simftest SCRIPT Generate simulated datasets to estimate the F-test probability for adding a model component. source CONTROL Execute a script file. statistic SETTING Change the fit statistic in use. steppar FIT Step through a range of parameter values; perform a fit at each step. syscall CONTROL Run a shell command. systematic MODEL Set the model systematic error. tclout CONTROL write xspec data to a tcl variable tcloutr CONTROL tclout with return value thaw (rthaw) FIT Allow a model parameter to vary during the fit. rthaw is for response parameters time CONTROL Display elapsed information. uncertain FIT Alias for error untie (runtie) MODEL Untie linked parameters. parameters. version CONTROL Print XSPEC version and build date/time weight FIT Change the weighting function used for chi-squared fits. writefits SCRIPT Write information about the current fit and errors to a FITS file. xsect SETTING Change the photoelectric absorption cross-sections in use. xset SETTING Modify a number of XSPEC internal switches time and other statistical runtie is for response 5.2. DESCRIPTION OF SYNTAX 5.2 89 Description of Syntax The individual commands are treated in alphabetical order in the following section. The novice would be well-served by reading the treatments of the data, model, newpar, and fit commands, in t hat order, then the other commands as needed. The write-up for each command includes a brief description of the purpose, an outline of the correct syntax, a more detailed discussion of the command assumptions and purpose, and a series of examples. Some commands have one or more subcommands that are similarly described following the command. In the command description, the syntax uses the following conventions. an argument to the command =:: defines as followed by ... a repeated string of arguments of the same type [ ] is an optional argument | indicates a choice between an argument of or Exceptional responses to the command prompt are : empty line or line containing only spaces and tab characters / nothing performed, prompt repeated any remaining arguments will have the values given on the last invocation of the command (Ctrl-D on Unix) same as quit /* skip input and return to prompt. Defaults for prompted values will be used. ?, ?command, command ? Print list of commands, or summary help for a single command 90 5.3 5.3.1 CHAPTER 5. XSPEC COMMANDS Control Commands autosave set frequency of saving commands Set or disable autosave, which saves the XSPEC environment to a file periodically. Syntax: autosave