forked from legolason/PyQSOFit
-
Notifications
You must be signed in to change notification settings - Fork 0
/
PyQSOFit.py
1625 lines (1363 loc) · 83.1 KB
/
PyQSOFit.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/python3
# A code for quasar spectrum fitting
# Last modified on 3/4/2021
# Auther: Hengxiao Guo AT UCI
# Email: hengxiaoguo AT gmail DOT com
# Co-Auther Shu Wang, Yue Shen, Wenke Ren
# version 1.1
# -------------------------------------------------
# fixed the flux error problem, previous error was underestimated by a factor of 1+z
# fixed QSO PCA
# fixed MC parameter unrenewed issue (Thanks Yuming Fu to point this out)
import glob
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import sfdmap
from scipy import interpolate
from scipy import integrate
from kapteyn import kmpfit
from PyAstronomy import pyasl
from astropy.io import fits
from astropy.cosmology import FlatLambdaCDM
from astropy.modeling.blackbody import blackbody_lambda
from astropy.table import Table
import warnings
warnings.filterwarnings("ignore")
class QSOFit():
def __init__(self, lam, flux, err, z, ra=- 999., dec=-999., plateid=None, mjd=None, fiberid=None, path=None,
and_mask=None, or_mask=None):
"""
Get the input data perpared for the QSO spectral fitting
Parameters:
-----------
lam: 1-D array with Npix
Observed wavelength in unit of Angstrom
flux: 1-D array with Npix
Observed flux density in unit of 10^{-17} erg/s/cm^2/Angstrom
err: 1-D array with Npix
1 sigma err with the same unit of flux
z: float number
redshift
ra, dec: float number, optional
the location of the source, right ascension and declination. The default number is 0
plateid, mjd, fiberid: integer number, optional
If the source is SDSS object, they have the plate ID, MJD and Fiber ID in their file herader.
path: str
the path of the input data
and_mask, or_mask: 1-D array with Npix, optional
the bad pixels defined from SDSS data, which can be got from SDSS datacube.
"""
self.lam = np.asarray(lam, dtype=np.float64)
self.flux = np.asarray(flux, dtype=np.float64)
self.err = np.asarray(err, dtype=np.float64)
self.z = z
self.and_mask = and_mask
self.or_mask = or_mask
self.ra = ra
self.dec = dec
self.plateid = plateid
self.mjd = mjd
self.fiberid = fiberid
self.path = path
def Fit(self, name=None, nsmooth=1, and_or_mask=True, reject_badpix=True, deredden=True, wave_range=None,
wave_mask=None, decomposition_host=True, BC03=False, Mi=None, npca_gal=5, npca_qso=20, Fe_uv_op=True,
Fe_flux_range=None, poly=False, BC=False, rej_abs=False, initial_guess=None, MC=True, n_trails=1,
linefit=True, tie_lambda=True, tie_width=True, tie_flux_1=True, tie_flux_2=True, save_result=True,
plot_fig=True, save_fig=True, plot_line_name=True, plot_legend=True, dustmap_path=None, save_fig_path=None,
save_fits_path=None, save_fits_name=None):
"""
Fit the QSO spectrum and get different decomposed components and corresponding parameters
Parameter:
----------
name: str, optinal
source name, Default is None. If None, it will use plateid+mjd+fiberid as the name. If there are no
such parameters, it will be empty.
nsmooth: integer number, optional
do n-pixel smoothing to the raw input flux and err spectra. The default is set to 1 (no smooth).
It will return the same array size. We note that smooth the raw data is not suggested, this function is in case of some fail-fitted low S/N spectra.
and_or_mask: bool, optional
If True, and and_mask or or_mask is not None, it will delete the masked pixels, and only return the remained pixels. Default: False
reject_badpix: bool, optional
reject 10 most possible outliers by the test of pointDistGESD. One important Caveat here is that this process will also delete narrow emission lines
in some high SN ratio object e.g., [OIII]. Only use it when you are definitely clear about what you are doing. It will return the remained pixels.
deredden: bool, optional
correct the Galactic extinction only if the RA and Dec are available. It will return the corrected flux with the same array size. Default: True.
wave_range: 2-element array, optional
trim input wavelength (lam) according to the min and max range of the input 2-element array, e.g.,
np.array([4000.,6000.]) in Rest frame range. Default: None
wave_mask: 2-D array
mask some absorption lines or sky lines in spectrum, e.g., np.array([[2200.,2300.]]), np.array([[5650.,5750.],[5850.,5900.]])
decomposition_host: bool, optional
If True, the host galaxy-QSO decomposition will be applied. If no more than 100 pixels are negative, the result will be applied. The Decomposition is
based on the PCA method of Yip et al. 2004 (AJ, 128, 585) & (128, 2603). Now the template is only available for redshift < 1.16 in specific absolute
magnitude bins. For galaxy, the global model has 10 PCA components and first 5 will enough to reproduce 98.37% galaxy spectra. For QSO, the global model
has 50, and the first 20 will reproduce 96.89% QSOs. If have i-band absolute magnitude, the Luminosity-redshift binned PCA components are available.
Then the first 10 PCA in each bin is enough to reproduce most QSO spectrum. Default: False
BC03: bool, optional
if True, it will use Bruzual1 & Charlot 2003 host model to fit spectrum, high shift host will be low resolution R ~ 300, the rest is R ~ 2000. Default: False
Mi: float, optional
the absolute magnitude of i band. It only works when decomposition_host is True. If not None, the Luminosity redshift binned PCA will be used to decompose
the spectrum. Default: None
npca_gal: int, optional
the number of galaxy PCA components applied. It only works when decomposition_host is True. The default is 5,
which is already account for 98.37% galaxies.
npca_qso: int, optional
the number of QSO PCA components applied. It only works when decomposition_host is True. The default is 20,
No matter the global or luminosity-redshift binned PCA is used, it can reproduce > 92% QSOs. The binned PCA
is better if have Mi information.
Fe_uv_op: bool, optional
if True, fit continuum with UV and optical FeII template. Default: True
Fe_flux_range: 1-D array, 2-D array or None, optional
if 1-D array was given, it should contain two parameters contain a range of wavelength. FeII flux within this range would be calculate and documented in the result fits file.
if 2-D array was given, it should contain a series of ranges. FeII flux within these ranges would be documented respectively.
if None was given, nothing would be return. Default: None
poly: bool, optional
if True, fit continuum with the polynomial component to account for the dust reddening. Default: False
BC: bool, optional
if True, fit continuum with Balmer continua from 1000 to 3646A. Default: False
rej_abs: bool, optional
if True, it will iterate the continuum fitting for deleting some 3 sigmas out continuum window points
(< 3500A), which might fall into the broad absorption lines. Default: False
initial_gauss: 1*14 array, optional
better initial value will help find a solution faster. Default initial is np.array([0., 3000., 0., 0.,
3000., 0., 1., -1.5, 0., 15000., 0.5, 0., 0., 0.]). First six parameters are flux scale, FWHM, small shift for wavelength for UV and optical FeII template,
respectively. The next two parameters are the power-law slope and intercept. The next three are the norm, Te, tau_BE in Balmer continuum model in
Dietrich et al. 2002. the last three parameters are a,b,c in polynomial function a*(x-3000)+b*x^2+c*x^3.
MC: bool, optional
if True, do the Monte Carlo simulation based on the input error array to produce the MC error array.
if False, the code will not save the error produced by kmpfit since it is biased and can not be trusted.
But it can be still output by in kmpfit attribute. Default: False
n_trails: int, optional
the number of trails of the MC process to produce the error array. The conservative number should be larger than 20. It only works when MC is True. Default: 20
linefit: bool, optional
if True, the emission line will be fitted. Default: True
tie_lambda: bool, optional
if True, line center with the same "vindex" will be tied together in the same line complex, this is always used to tie e.g., NII. Default: False
tie_width: bool, optional
if True, line sigma with the same "windex" will be tied together in the same line complex, Default: False
tie_flux_1: bool, optional
if True, line flux with the flag "findex = 1" will be tied to the ratio of fvalue. To fix the narrow line flux, the user should fix the line width
"tie_width" first. Default: False
tie_flux_2: bool, optional
if True, line flux with the flag "findex = 2" will be tied to the ratio of fvalue. To fix the narrow line flux, the user should fix the line width
"tie_width" first. Default: False
save_result: bool, optional
if True, all the fitting results will be saved to a fits file, Default: True
plot_fig: bool, optional
if True, the fitting results will be plotted. Default: True
save_fig: bool, optional
if True, the figure will be saved, and the path can be set by "save_fig_path". Default: True
plot_line_name: bool, optional
if True, serval main emission lines will be plotted in the first panel of the output figure. Default: False
plot_legend: bool, optional
if True, open legend in the first panel of the output figure. Default: False
dustmap_path: str, optional
if Deredden is True, the dustmap_path must be set. If None, the default "dustmap_path" is set to "path"
save_fig_path: str, optional
the output path of the figure. If None, the default "save_fig_path" is set to "path"
save_fit_path: str, optional
the output path of the result fits. If None, the default "save_fits_path" is set to "path"
save_fit_name: str, optional
the output name of the result fits. Default: "result.fits"
Return:
-----------
.wave: array
the rest wavelength, some pixels have been removed.
.flux: array
the rest flux. Dereddened and *(1+z) flux.
.err: array
the error.
.wave_prereduced: array
the wavelength after removing bad pixels, masking, deredden, spectral trim, and smoothing.
.flux_prereduced: array
the flux after removing bad pixels, masking, deredden, spectral trim, and smoothing.
.err_prereduced: array
the error after removing bad pixels, masking, deredden, spectral trim, and smoothing.
.host: array
the model of host galaxy from PCA method
.qso: array
the model of a quasar from PCA method.
.SN_ratio_conti: float
the mean S/N ratio of 1350, 3000 and 5100A.
.conti_fit.: structure
all the kmpfit continuum fitting results, including best-fit parameters and Chisquare, etc. For details,
see https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html
.f_conti_model: array
the continuum model including power-law, polynomial, optical/UV FeII, Balmer continuum.
.f_bc_model: array
the Balmer continuum model.
.f_fe_uv_model: array
the UV FeII model.
.f_fe_op_model: array
the optical FeII model.
.f_pl_model: array
the power-law model.
.f_poly_model: array
the polynomial model.
.PL_poly_BC: array
The combination of Powerlaw, polynomial and Balmer continuum model.
.line_flux: array
the emission line flux after subtracting the .f_conti_model.
.line_fit: structrue
kmpfit line fitting results for last complexes (From Lya to Ha) , including best-fit parameters, errors (kmpfit derived) and Chisquare, etc. For details,
see https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html
.gauss_result: array
3*n Gaussian parameters for all lines in the format of [scale, centerwave, sigma ], n is number of Gaussians for all complexes.
.conti_result: array
continuum parameters, including widely used continuum parameters and monochromatic flux at 1350, 3000
and5100 Angstrom, etc. The corresponding names are listed in .conti_result_name. For all continuum fitting results,
go to .conti_fit.params.
.conti_result_name: array
the names for .conti_result.
.line_result: array
emission line parameters, including FWHM, sigma, EW, measured from whole model of each main broad emission line covered,
and fitting parameters of each Gaussian component. The corresponding names are listed in .line_result_name.
.line_result_name: array
the names for .line_result.
.uniq_linecomp_sort: array
the sorted complex names.
.all_comp_range: array
the start and end wavelength for each complex. e.g., Hb in [4640. 5100.].
.linelist: array
the information listed in the qsopar.fits.
"""
self.name = name
self.wave_range = wave_range
self.wave_mask = wave_mask
self.BC03 = BC03
self.Mi = Mi
self.npca_gal = npca_gal
self.npca_qso = npca_qso
self.initial_guess = initial_guess
self.Fe_uv_op = Fe_uv_op
self.Fe_flux_range = Fe_flux_range
self.poly = poly
self.BC = BC
self.rej_abs = rej_abs
self.MC = MC
self.n_trails = n_trails
self.tie_lambda = tie_lambda
self.tie_width = tie_width
self.tie_flux_1 = tie_flux_1
self.tie_flux_2 = tie_flux_2
self.plot_line_name = plot_line_name
self.plot_legend = plot_legend
self.save_fig = save_fig
# get the source name in plate-mjd-fiber, if no then None
if name is None:
if np.array([self.plateid, self.mjd, self.fiberid]).any() is not None:
self.sdss_name = str(self.plateid).zfill(4)+'-'+str(self.mjd)+'-'+str(self.fiberid).zfill(4)
else:
if self.plateid is None:
self.plateid = 0
if self.mjd is None:
self.mjd = 0
if self.fiberid is None:
self.fiberid = 0
self.sdss_name = ''
else:
self.sdss_name = name
# set default path for figure and fits
if save_result == True and save_fits_path == None:
save_fits_path = self.path
if save_fig == True and save_fig_path == None:
save_fig_path = self.path
if save_fits_name == None:
if self.sdss_name == '':
save_fits_name = 'result'
else:
save_fits_name = self.sdss_name
else:
save_fits_name = save_fits_name
# deal with pixels with error equal 0 or inifity
ind_gooderror = np.where((self.err != 0) & ~np.isinf(self.err), True, False)
err_good = self.err[ind_gooderror]
flux_good = self.flux[ind_gooderror]
lam_good = self.lam[ind_gooderror]
if (self.and_mask is not None) & (self.or_mask is not None):
and_mask_good = self.and_mask[ind_gooderror]
or_mask_good = self.or_mask[ind_gooderror]
del self.and_mask, self.or_mask
self.and_mask = and_mask_good
self.or_mask = or_mask_good
del self.err, self.flux, self.lam
self.err = err_good
self.flux = flux_good
self.lam = lam_good
if nsmooth is not None:
self.flux = self.Smooth(self.flux, nsmooth)
self.err = self.Smooth(self.err, nsmooth)
if (and_or_mask == True) and (self.and_mask is not None or self.or_mask is not None):
self._MaskSdssAndOr(self.lam, self.flux, self.err, self.and_mask, self.or_mask)
if reject_badpix == True:
self._RejectBadPix(self.lam, self.flux, self.err)
if wave_range is not None:
self._WaveTrim(self.lam, self.flux, self.err, self.z)
if wave_mask is not None:
self._WaveMsk(self.lam, self.flux, self.err, self.z)
if deredden == True and self.ra != -999. and self.dec != -999.:
self._DeRedden(self.lam, self.flux, self.err, self.ra, self.dec, dustmap_path)
self._RestFrame(self.lam, self.flux, self.err, self.z)
self._CalculateSN(self.wave, self.flux)
self._OrignialSpec(self.wave, self.flux, self.err)
# do host decomposition --------------
if self.z < 1.16 and decomposition_host == True:
self._DoDecomposition(self.wave, self.flux, self.err, self.path)
else:
self.decomposed = False
if self.z > 1.16 and decomposition_host == True:
print('redshift larger than 1.16 is not allowed for host '
'decomposion!')
# fit continuum --------------------
self._DoContiFit(self.wave, self.flux, self.err, self.ra, self.dec, self.plateid, self.mjd, self.fiberid)
# fit line
if linefit == True:
self._DoLineFit(self.wave, self.line_flux, self.err, self.conti_fit)
else:
self.ncomp = 0
# save data -------
if save_result == True:
if linefit == False:
self.line_result = np.array([])
self.line_result_type = np.array([])
self.line_result_name = np.array([])
self._SaveResult(self.conti_result, self.conti_result_type, self.conti_result_name, self.line_result,
self.line_result_type, self.line_result_name, save_fits_path, save_fits_name)
# plot fig and save ------
if plot_fig == True:
if linefit == False:
self.gauss_result = np.array([])
self.all_comp_range = np.array([])
self.uniq_linecomp_sort = np.array([])
self._PlotFig(self.ra, self.dec, self.z, self.wave, self.flux, self.err, decomposition_host, linefit,
self.tmp_all, self.gauss_result, self.f_conti_model, self.conti_fit, self.all_comp_range,
self.uniq_linecomp_sort, self.line_flux, save_fig_path)
def _MaskSdssAndOr(self, lam, flux, err, and_mask, or_mask):
"""
Remove SDSS and_mask and or_mask points are not zero
Parameter:
----------
lam: wavelength
flux: flux
err: 1 sigma error
and_mask: SDSS flag "and_mask", mask out all non-zero pixels
or_mask: SDSS flag "or_mask", mask out all npn-zero pixels
Retrun:
---------
return the same size array of wavelength, flux, error
"""
ind_and_or = np.where((and_mask == 0) & (or_mask == 0), True, False)
del self.lam, self.flux, self.err
self.lam, self.flux, self.err = lam[ind_and_or], flux[ind_and_or], err[ind_and_or]
def _RejectBadPix(self, lam, flux, err):
"""
Reject 10 most possiable outliers, input wavelength, flux and error. Return a different size wavelength,
flux, and error.
"""
# -----remove bad pixels, but not for high SN spectrum------------
ind_bad = pyasl.pointDistGESD(flux, 10)
wv = np.asarray([i for j, i in enumerate(lam) if j not in ind_bad[1]], dtype=np.float64)
fx = np.asarray([i for j, i in enumerate(flux) if j not in ind_bad[1]], dtype=np.float64)
er = np.asarray([i for j, i in enumerate(err) if j not in ind_bad[1]], dtype=np.float64)
del self.lam, self.flux, self.err
self.lam, self.flux, self.err = wv, fx, er
return self.lam, self.flux, self.err
def _WaveTrim(self, lam, flux, err, z):
"""
Trim spectrum with a range in the rest frame.
"""
# trim spectrum e.g., local fit emiision lines
ind_trim = np.where((lam/(1.+z) > self.wave_range[0]) & (lam/(1.+z) < self.wave_range[1]), True, False)
del self.lam, self.flux, self.err
self.lam, self.flux, self.err = lam[ind_trim], flux[ind_trim], err[ind_trim]
if len(self.lam) < 100:
raise RuntimeError("No enough pixels in the input wave_range!")
return self.lam, self.flux, self.err
def _WaveMsk(self, lam, flux, err, z):
"""Block the bad pixels or absorption lines in spectrum."""
for msk in range(len(self.wave_mask)):
try:
ind_not_mask = ~np.where((lam/(1.+z) > self.wave_mask[msk, 0]) & (lam/(1.+z) < self.wave_mask[msk, 1]),
True, False)
except IndexError:
raise RuntimeError("Wave_mask should be 2D array,e.g., np.array([[2000,3000],[3100,4000]]).")
del self.lam, self.flux, self.err
self.lam, self.flux, self.err = lam[ind_not_mask], flux[ind_not_mask], err[ind_not_mask]
lam, flux, err = self.lam, self.flux, self.err
return self.lam, self.flux, self.err
def _DeRedden(self, lam, flux, err, ra, dec, dustmap_path):
"""Correct the Galactic extinction"""
m = sfdmap.SFDMap(dustmap_path)
zero_flux = np.where(flux == 0, True, False)
flux[zero_flux] = 1e-10
flux_unred = pyasl.unred(lam, flux, m.ebv(ra, dec))
err_unred = err*flux_unred/flux
flux_unred[zero_flux] = 0
del self.flux, self.err
self.flux = flux_unred
self.err = err_unred
return self.flux
def _RestFrame(self, lam, flux, err, z):
"""Move wavelenth and flux to rest frame"""
self.wave = lam/(1.+z)
self.flux = flux*(1.+z)
self.err = err*(1.+z)
return self.wave, self.flux, self.err
def _OrignialSpec(self, wave, flux, err):
"""save the orignial spectrum before host galaxy decompsition"""
self.wave_prereduced = wave
self.flux_prereduced = flux
self.err_prereduced = err
def _CalculateSN(self, wave, flux):
"""calculate the spectral SN ratio for 1350, 3000, 5100A, return the mean value of Three spots"""
if ((wave.min() < 1350. and wave.max() > 1350.) or (wave.min() < 3000. and wave.max() > 3000.) or (
wave.min() < 5100. and wave.max() > 5100.)):
ind5100 = np.where((wave > 5080.) & (wave < 5130.), True, False)
ind3000 = np.where((wave > 3000.) & (wave < 3050.), True, False)
ind1350 = np.where((wave > 1325.) & (wave < 1375.), True, False)
tmp_SN = np.array([flux[ind5100].mean()/flux[ind5100].std(), flux[ind3000].mean()/flux[ind3000].std(),
flux[ind1350].mean()/flux[ind1350].std()])
tmp_SN = tmp_SN[~np.isnan(tmp_SN)]
self.SN_ratio_conti = tmp_SN.mean()
else:
self.SN_ratio_conti = -1.
return self.SN_ratio_conti
def _DoDecomposition(self, wave, flux, err, path):
"""Decompose the host galaxy from QSO"""
datacube = self._HostDecompose(self.wave, self.flux, self.err, self.z, self.Mi, self.npca_gal, self.npca_qso,
path)
# for some negtive host templete, we do not do the decomposition
if np.sum(np.where(datacube[3, :] < 0., True, False)) > 100:
self.host = np.zeros(len(wave))
self.decomposed = False
print('Get negtive host galaxy flux larger than 100 pixels, '
'decomposition is not applied!')
else:
self.decomposed = True
del self.wave, self.flux, self.err
self.wave = datacube[0, :]
# block OIII, ha,NII,SII,OII,Ha,Hb,Hr,hdelta
line_mask = np.where(
(self.wave < 4970.) & (self.wave > 4950.) | (self.wave < 5020.) & (self.wave > 5000.) | (
self.wave < 6590.) & (self.wave > 6540.) | (self.wave < 6740.) & (self.wave > 6710.) | (
self.wave < 3737.) & (self.wave > 3717.) | (self.wave < 4872.) & (self.wave > 4852.) | (
self.wave < 4350.) & (self.wave > 4330.) | (self.wave < 4111.) & (self.wave > 4091.), True,
False)
f = interpolate.interp1d(self.wave[~line_mask], datacube[3, :][~line_mask], bounds_error=False,
fill_value=0)
masked_host = f(self.wave)
self.flux = datacube[1, :]-masked_host # QSO flux without host
self.err = datacube[2, :]
self.host = datacube[3, :]
self.qso = datacube[4, :]
self.host_data = datacube[1, :]-self.qso
return self.wave, self.flux, self.err
def _HostDecompose(self, wave, flux, err, z, Mi, npca_gal, npca_qso, path):
"""
core function to do host decomposition
#Wave is the obs frame wavelength, n_gal and n_qso are the number of eigenspectra used to fit
#If Mi is None then the qso use the globle ones to fit. If not then use the redshift-luminoisty binded ones to fit
#See details:
#Yip, C. W., Connolly, A. J., Szalay, A. S., et al. 2004a, AJ, 128, 585
#Yip, C. W., Connolly, A. J., Vanden Berk, D. E., et al. 2004b, AJ, 128, 2603
"""
# read galaxy and qso eigenspectra -----------------------------------
if self.BC03 == False:
galaxy = fits.open(path+'pca/Yip_pca_templates/gal_eigenspec_Yip2004.fits')
gal = galaxy[1].data
wave_gal = gal['wave'].flatten()
flux_gal = gal['pca'].reshape(gal['pca'].shape[1], gal['pca'].shape[2])
if self.BC03 == True:
cc = 0
flux03 = np.array([])
for i in glob.glob(path+'/bc03/*.gz'):
cc = cc+1
gal_temp = np.genfromtxt(i)
wave_gal = gal_temp[:, 0]
flux03 = np.concatenate((flux03, gal_temp[:, 1]))
flux_gal = np.array(flux03).reshape(cc, -1)
if Mi is None:
quasar = fits.open(path+'pca/Yip_pca_templates/qso_eigenspec_Yip2004_global.fits')
else:
if -24 < Mi <= -22 and 0.08 <= z < 0.53:
quasar = fits.open(path+'pca/Yip_pca_templates/qso_eigenspec_Yip2004_CZBIN1.fits')
elif -26 < Mi <= -24 and 0.08 <= z < 0.53:
quasar = fits.open(path+'pca/Yip_pca_templates/qso_eigenspec_Yip2004_DZBIN1.fits')
elif -24 < Mi <= -22 and 0.53 <= z < 1.16:
quasar = fits.open(path+'pca/Yip_pca_templates/qso_eigenspec_Yip2004_BZBIN2.fits')
elif -26 < Mi <= -24 and 0.53 <= z < 1.16:
quasar = fits.open(path+'pca/Yip_pca_templates/qso_eigenspec_Yip2004_CZBIN2.fits')
elif -28 < Mi <= -26 and 0.53 <= z < 1.16:
quasar = fits.open(path+'pca/Yip_pca_templates/qso_eigenspec_Yip2004_DZBIN2.fits')
else:
raise RuntimeError('Host galaxy template is not available for this redshift and Magnitude!')
qso = quasar[1].data
wave_qso = qso['wave'].flatten()
flux_qso = qso['pca'].reshape(qso['pca'].shape[1], qso['pca'].shape[2])
# get the shortest wavelength range
wave_min = max(wave.min(), wave_gal.min(), wave_qso.min())
wave_max = min(wave.max(), wave_gal.max(), wave_qso.max())
ind_data = np.where((wave > wave_min) & (wave < wave_max), True, False)
ind_gal = np.where((wave_gal > wave_min-1.) & (wave_gal < wave_max+1.), True, False)
ind_qso = np.where((wave_qso > wave_min-1.) & (wave_qso < wave_max+1.), True, False)
flux_gal_new = np.zeros(flux_gal.shape[0]*flux[ind_data].shape[0]).reshape(flux_gal.shape[0],
flux[ind_data].shape[0])
flux_qso_new = np.zeros(flux_qso.shape[0]*flux[ind_data].shape[0]).reshape(flux_qso.shape[0],
flux[ind_data].shape[0])
for i in range(flux_gal.shape[0]):
fgal = interpolate.interp1d(wave_gal[ind_gal], flux_gal[i, ind_gal], bounds_error=False, fill_value=0)
flux_gal_new[i, :] = fgal(wave[ind_data])
for i in range(flux_qso.shape[0]):
fqso = interpolate.interp1d(wave_qso[ind_qso], flux_qso[i, ind_qso], bounds_error=False, fill_value=0)
flux_qso_new[i, :] = fqso(wave[ind_data])
wave_new = wave[ind_data]
flux_new = flux[ind_data]
err_new = err[ind_data]
flux_temp = np.vstack((flux_gal_new[0:npca_gal, :], flux_qso_new[0:npca_qso, :]))
res = np.linalg.lstsq(flux_temp.T, flux_new)[0]
host_flux = np.dot(res[0:npca_gal], flux_temp[0:npca_gal])
qso_flux = np.dot(res[npca_gal:], flux_temp[npca_gal:])
data_cube = np.vstack((wave_new, flux_new, err_new, host_flux, qso_flux))
ind_f4200 = np.where((wave_new > 4160.) & (wave_new < 4210.), True, False)
frac_host_4200 = np.sum(host_flux[ind_f4200])/np.sum(flux_new[ind_f4200])
ind_f5100 = np.where((wave_new > 5080.) & (wave_new < 5130.), True, False)
frac_host_5100 = np.sum(host_flux[ind_f5100])/np.sum(flux_new[ind_f5100])
return data_cube # ,frac_host_4200,frac_host_5100
def _DoContiFit(self, wave, flux, err, ra, dec, plateid, mjd, fiberid):
"""Fit the continuum with PL, Polynomial, UV/optical FeII, Balmer continuum"""
global fe_uv, fe_op
fe_uv = np.genfromtxt(self.path+'fe_uv.txt')
fe_op = np.genfromtxt(self.path+'fe_optical.txt')
# do continuum fit--------------------------
window_all = np.array(
[[1150., 1170.], [1275., 1290.], [1350., 1360.], [1445., 1465.], [1690., 1705.], [1770., 1810.],
[1970., 2400.], [2480., 2675.], [2925., 3400.], [3775., 3832.], [4000., 4050.], [4200., 4230.],
[4435., 4640.], [5100., 5535.], [6005., 6035.], [6110., 6250.], [6800., 7000.], [7160., 7180.],
[7500., 7800.], [8050., 8150.], ])
tmp_all = np.array([np.repeat(False, len(wave))]).flatten()
for jj in range(len(window_all)):
tmp = np.where((wave > window_all[jj, 0]) & (wave < window_all[jj, 1]), True, False)
tmp_all = np.any([tmp_all, tmp], axis=0)
if wave[tmp_all].shape[0] < 10:
print('Continuum fitting pixel < 10. ')
# set initial paramiters for continuum
if self.initial_guess is not None:
pp0 = self.initial_guess
else:
pp0 = np.array([0., 3000., 0., 0., 3000., 0., 1., -1.5, 0., 15000., 0.5, 0., 0., 0.])
conti_fit = kmpfit.Fitter(residuals=self._residuals, data=(wave[tmp_all], flux[tmp_all], err[tmp_all]))
tmp_parinfo = [{'limits': (0., 10.**10)}, {'limits': (1200., 10000.)}, {'limits': (-0.01, 0.01)},
{'limits': (0., 10.**10)}, {'limits': (1200., 10000.)}, {'limits': (-0.01, 0.01)},
{'limits': (0., 10.**10)}, {'limits': (-5., 3.)}, {'limits': (0., 10.**10)},
{'limits': (10000., 50000.)}, {'limits': (0.1, 2.)}, None, None, None, ]
conti_fit.parinfo = tmp_parinfo
conti_fit.fit(params0=pp0)
# Perform one iteration to remove 3sigma pixel below the first continuum fit
# to avoid the continuum windows falls within a BAL trough
if self.rej_abs == True:
if self.poly == True:
tmp_conti = (conti_fit.params[6]*(wave[tmp_all]/3000.0)**conti_fit.params[7]+self.F_poly_conti(
wave[tmp_all], conti_fit.params[11:]))
else:
tmp_conti = (conti_fit.params[6]*(wave[tmp_all]/3000.0)**conti_fit.params[7])
ind_noBAL = ~np.where(((flux[tmp_all] < tmp_conti-3.*err[tmp_all]) & (wave[tmp_all] < 3500.)), True, False)
f = kmpfit.Fitter(residuals=self._residuals, data=(
wave[tmp_all][ind_noBAL], self.Smooth(flux[tmp_all][ind_noBAL], 10), err[tmp_all][ind_noBAL]))
conti_fit.parinfo = tmp_parinfo
conti_fit.fit(params0=pp0)
# calculate continuum luminoisty
L = self._L_conti(wave, conti_fit.params)
# calculate FeII flux
Fe_flux_result = np.array([])
Fe_flux_type = np.array([])
Fe_flux_name = np.array([])
if self.Fe_flux_range is not None:
Fe_flux_result, Fe_flux_type, Fe_flux_name = self.Get_Fe_flux(self.Fe_flux_range, conti_fit.params[:6])
# get conti result -----------------------------
if self.MC == True and self.n_trails > 0:
# calculate MC err
conti_para_std, all_L_std, Fe_flux_std = self._conti_mc(self.wave[tmp_all], self.flux[tmp_all],
self.err[tmp_all], pp0, conti_fit.parinfo,
self.n_trails)
self.conti_result = np.array(
[ra, dec, str(plateid), str(mjd), str(fiberid), self.z, self.SN_ratio_conti, conti_fit.params[0],
conti_para_std[0], conti_fit.params[1], conti_para_std[1], conti_fit.params[2], conti_para_std[2],
conti_fit.params[3], conti_para_std[3], conti_fit.params[4], conti_para_std[4], conti_fit.params[5],
conti_para_std[5], conti_fit.params[6], conti_para_std[6], conti_fit.params[7], conti_para_std[7],
conti_fit.params[8], conti_para_std[8], conti_fit.params[9], conti_para_std[9], conti_fit.params[10],
conti_para_std[10], conti_fit.params[11], conti_para_std[11], conti_fit.params[12], conti_para_std[12],
conti_fit.params[13], conti_para_std[13], L[0], all_L_std[0], L[1], all_L_std[1], L[2], all_L_std[2]])
self.conti_result_type = np.array(
['float', 'float', 'int', 'int', 'int', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
'float', 'float', 'float', 'float', 'float', 'float', 'float'])
self.conti_result_name = np.array(
['ra', 'dec', 'plateid', 'MJD', 'fiberid', 'redshift', 'SN_ratio_conti', 'Fe_uv_norm', 'Fe_uv_norm_err',
'Fe_uv_FWHM', 'Fe_uv_FWHM_err', 'Fe_uv_shift', 'Fe_uv_shift_err', 'Fe_op_norm', 'Fe_op_norm_err',
'Fe_op_FWHM', 'Fe_op_FWHM_err', 'Fe_op_shift', 'Fe_op_shift_err', 'PL_norm', 'PL_norm_err', 'PL_slope',
'PL_slope_err', 'Blamer_norm', 'Blamer_norm_err', 'Balmer_Te', 'Balmer_Te_err', 'Balmer_Tau',
'Balmer_Tau_err', 'POLY_a', 'POLY_a_err', 'POLY_b', 'POLY_b_err', 'POLY_c', 'POLY_c_err', 'L1350',
'L1350_err', 'L3000', 'L3000_err', 'L5100', 'L5100_err'])
for iii in range(Fe_flux_result.shape[0]):
self.conti_result = np.append(self.conti_result, [Fe_flux_result[iii], Fe_flux_std[iii]])
self.conti_result_type = np.append(self.conti_result_type, [Fe_flux_type[iii], Fe_flux_type[iii]])
self.conti_result_name = np.append(self.conti_result_name,
[Fe_flux_name[iii], Fe_flux_name[iii]+'_err'])
else:
self.conti_result = np.array(
[ra, dec, str(plateid), str(mjd), str(fiberid), self.z, self.SN_ratio_conti, conti_fit.params[0],
conti_fit.params[1], conti_fit.params[2], conti_fit.params[3], conti_fit.params[4],
conti_fit.params[5], conti_fit.params[6], conti_fit.params[7], conti_fit.params[8],
conti_fit.params[9], conti_fit.params[10], conti_fit.params[11], conti_fit.params[12],
conti_fit.params[13], L[0], L[1], L[2]])
self.conti_result_type = np.array(
['float', 'float', 'int', 'int', 'int', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
'float'])
self.conti_result_name = np.array(
['ra', 'dec', 'plateid', 'MJD', 'fiberid', 'redshift', 'SN_ratio_conti', 'Fe_uv_norm', 'Fe_uv_FWHM',
'Fe_uv_shift', 'Fe_op_norm', 'Fe_op_FWHM', 'Fe_op_shift', 'PL_norm', 'PL_slope', 'Blamer_norm',
'Balmer_Te', 'Balmer_Tau_e', 'POLY_a', 'POLY_b', 'POLY_c', 'L1350', 'L3000', 'L5100'])
self.conti_result = np.append(self.conti_result, Fe_flux_result)
self.conti_result_type = np.append(self.conti_result_type, Fe_flux_type)
self.conti_result_name = np.append(self.conti_result_name, Fe_flux_name)
self.conti_fit = conti_fit
self.tmp_all = tmp_all
# save different models--------------------
f_fe_mgii_model = self.Fe_flux_mgii(wave, conti_fit.params[0:3])
f_fe_balmer_model = self.Fe_flux_balmer(wave, conti_fit.params[3:6])
f_pl_model = conti_fit.params[6]*(wave/3000.0)**conti_fit.params[7]
f_bc_model = self.Balmer_conti(wave, conti_fit.params[8:11])
f_poly_model = self.F_poly_conti(wave, conti_fit.params[11:])
f_conti_model = (f_pl_model+f_fe_mgii_model+f_fe_balmer_model+f_poly_model+f_bc_model)
line_flux = flux-f_conti_model
self.f_conti_model = f_conti_model
self.f_bc_model = f_bc_model
self.f_fe_uv_model = f_fe_mgii_model
self.f_fe_op_model = f_fe_balmer_model
self.f_pl_model = f_pl_model
self.f_poly_model = f_poly_model
self.line_flux = line_flux
self.PL_poly_BC = f_pl_model+f_poly_model+f_bc_model
return self.conti_result, self.conti_result_name
def _L_conti(self, wave, pp):
"""Calculate continuum Luminoisity at 1350,3000,5100A"""
conti_flux = pp[6]*(wave/3000.0)**pp[7]+self.F_poly_conti(wave, pp[11:])
# plt.plot(wave,conti_flux)
L = np.array([])
for LL in zip([1350., 3000., 5100.]):
if wave.max() > LL[0] and wave.min() < LL[0]:
L_tmp = np.asarray([np.log10(
LL[0]*self.Flux2L(conti_flux[np.where(abs(wave-LL[0]) < 5., True, False)].mean(), self.z))])
else:
L_tmp = np.array([-1.])
L = np.concatenate([L, L_tmp]) # save log10(L1350,L3000,L5100)
return L
def _f_conti_all(self, xval, pp):
"""
Continuum components described by 14 parameters
pp[0]: norm_factor for the MgII Fe_template
pp[1]: FWHM for the MgII Fe_template
pp[2]: small shift of wavelength for the MgII Fe template
pp[3:5]: same as pp[0:2] but for the Hbeta/Halpha Fe template
pp[6]: norm_factor for continuum f_lambda = (lambda/3000.0)^{-alpha}
pp[7]: slope for the power-law continuum
pp[8:10]: norm, Te and Tau_e for the Balmer continuum at <3646 A
pp[11:13]: polynomial for the continuum
"""
# iron flux for MgII line region
f_Fe_MgII = self.Fe_flux_mgii(xval, pp[0:3])
# iron flux for balmer line region
f_Fe_Balmer = self.Fe_flux_balmer(xval, pp[3:6])
# power-law continuum
f_pl = pp[6]*(xval/3000.0)**pp[7]
# Balmer continuum
f_conti_BC = self.Balmer_conti(xval, pp[8:11])
# polynormal conponent for reddened spectra
f_poly = self.F_poly_conti(xval, pp[11:])
if self.Fe_uv_op == True and self.poly == False and self.BC == False:
yval = f_pl+f_Fe_MgII+f_Fe_Balmer
elif self.Fe_uv_op == True and self.poly == True and self.BC == False:
yval = f_pl+f_Fe_MgII+f_Fe_Balmer+f_poly
elif self.Fe_uv_op == True and self.poly == False and self.BC == True:
yval = f_pl+f_Fe_MgII+f_Fe_Balmer+f_conti_BC
elif self.Fe_uv_op == False and self.poly == True and self.BC == False:
yval = f_pl+f_poly
elif self.Fe_uv_op == False and self.poly == False and self.BC == False:
yval = f_pl
elif self.Fe_uv_op == False and self.poly == False and self.BC == True:
yval = f_pl+f_conti_BC
elif self.Fe_uv_op == True and self.poly == True and self.BC == True:
yval = f_pl+f_Fe_MgII+f_Fe_Balmer+f_poly+f_conti_BC
elif self.Fe_uv_op == False and self.poly == True and self.BC == True:
yval = f_pl+f_Fe_Balmer+f_poly+f_conti_BC
else:
raise RuntimeError('No this option for Fe_uv_op, poly and BC!')
return yval
def _residuals(self, pp, data):
"""Continual residual function used in kmpfit"""
xval, yval, weight = data
return (yval-self._f_conti_all(xval, pp))/weight
# ---------MC error for continuum parameters-------------------
def _conti_mc(self, x, y, err, pp0, pp_limits, n_trails):
"""Calculate the continual parameters' Monte carlo errrors"""
all_para = np.zeros(len(pp0)*n_trails).reshape(len(pp0), n_trails)
all_L = np.zeros(3*n_trails).reshape(3, n_trails)
n_Fe_flux = np.array(self.Fe_flux_range).flatten().shape[0]//2
all_Fe_flux = np.zeros(n_Fe_flux*n_trails).reshape(n_Fe_flux, n_trails)
for tra in range(n_trails):
flux = y+np.random.randn(len(y))*err
conti_fit = kmpfit.Fitter(residuals=self._residuals, data=(x, flux, err), maxiter=50)
conti_fit.parinfo = pp_limits
conti_fit.fit(params0=pp0)
all_para[:, tra] = conti_fit.params
all_L[:, tra] = np.asarray(self._L_conti(x, conti_fit.params))
if self.Fe_flux_range is not None:
Fe_flux_result, Fe_flux_type, Fe_flux_name = self.Get_Fe_flux(self.Fe_flux_range, conti_fit.params[:6])
all_Fe_flux[:, tra] = Fe_flux_result
all_para_std = all_para.std(axis=1)
all_L_std = all_L.std(axis=1)
all_Fe_flux_std = all_Fe_flux.std(axis=1)
return all_para_std, all_L_std, all_Fe_flux_std
# line function-----------
def _DoLineFit(self, wave, line_flux, err, f):
"""Fit the emission lines with Gaussian profile """
# remove abosorbtion line in emission line region
# remove the pixels below continuum
ind_neg_line = ~np.where(((((wave > 2700.) & (wave < 2900.)) | ((wave > 1700.) & (wave < 1970.)) | (
(wave > 1500.) & (wave < 1700.)) | ((wave > 1290.) & (wave < 1450.)) | (
(wave > 1150.) & (wave < 1290.))) & (line_flux < -err)), True, False)
# read line parameter
linepara = fits.open(self.path+'qsopar.fits')
linelist = linepara[1].data
self.linelist = linelist
ind_kind_line = np.where((linelist['lambda'] > wave.min()) & (linelist['lambda'] < wave.max()), True, False)
if ind_kind_line.any() == True:
# sort complex name with line wavelength
uniq_linecomp, uniq_ind = np.unique(linelist['compname'][ind_kind_line], return_index=True)
uniq_linecomp_sort = uniq_linecomp[linelist['lambda'][ind_kind_line][uniq_ind].argsort()]
ncomp = len(uniq_linecomp_sort)
compname = linelist['compname']
allcompcenter = np.sort(linelist['lambda'][ind_kind_line][uniq_ind])
# loop over each complex and fit n lines simutaneously
comp_result = np.array([])
comp_result_type = np.array([])
comp_result_name = np.array([])
gauss_result = np.array([])
gauss_result_type = np.array([])
gauss_result_name = np.array([])
all_comp_range = np.array([])
fur_result = np.array([])
fur_result_type = np.array([])
fur_result_name = np.array([])
for ii in range(ncomp):
compcenter = allcompcenter[ii]
ind_line = np.where(linelist['compname'] == uniq_linecomp_sort[ii], True, False) # get line index
nline_fit = np.sum(ind_line) # n line in one complex
linelist_fit = linelist[ind_line]
# n gauss in each line
ngauss_fit = np.asarray(linelist_fit['ngauss'], dtype=int)
# for iitmp in range(nline_fit): # line fit together
comp_range = [linelist_fit[0]['minwav'], linelist_fit[0]['maxwav']] # read complex range from table
all_comp_range = np.concatenate([all_comp_range, comp_range])
# ----tie lines--------
self._do_tie_line(linelist, ind_line)
# get the pixel index in complex region and remove negtive abs in line region
ind_n = np.where((wave > comp_range[0]) & (wave < comp_range[1]) & (ind_neg_line == True), True, False)
if np.sum(ind_n) > 10:
# call kmpfit for lines
line_fit = self._do_line_kmpfit(linelist, line_flux, ind_line, ind_n, nline_fit, ngauss_fit)
# calculate MC err
if self.MC == True and self.n_trails > 0:
all_para_std, fwhm_std, sigma_std, ew_std, peak_std, area_std = self._line_mc(
np.log(wave[ind_n]), line_flux[ind_n], err[ind_n], self.line_fit_ini, self.line_fit_par,
self.n_trails, compcenter)
# ----------------------get line fitting results----------------------
# complex parameters
# tie lines would reduce the number of parameters increasing the dof
dof_fix = 0
if self.tie_lambda == True:
dof_fix += np.max((len(self.ind_tie_vindex1), 1))-1
dof_fix += np.max((len(self.ind_tie_vindex2), 1))-1
if self.tie_width == True:
dof_fix += np.max((len(self.ind_tie_windex1), 1))-1
dof_fix += np.max((len(self.ind_tie_windex2), 1))-1
if self.tie_flux_1 == True:
dof_fix += np.max((len(self.ind_tie_findex1), 1))-1
dof_fix += np.max((len(self.ind_tie_findex2), 1))-1
comp_result_tmp = np.array(
[[linelist['compname'][ind_line][0]], [line_fit.status], [line_fit.chi2_min],
[line_fit.chi2_min/(line_fit.dof+dof_fix)], [line_fit.niter],
[line_fit.dof+dof_fix]]).flatten()
comp_result_type_tmp = np.array(['str', 'int', 'float', 'float', 'int', 'int'])
comp_result_name_tmp = np.array(
[str(ii+1)+'_complex_name', str(ii+1)+'_line_status', str(ii+1)+'_line_min_chi2',
str(ii+1)+'_line_red_chi2', str(ii+1)+'_niter', str(ii+1)+'_ndof'])
comp_result = np.concatenate([comp_result, comp_result_tmp])
comp_result_name = np.concatenate([comp_result_name, comp_result_name_tmp])
comp_result_type = np.concatenate([comp_result_type, comp_result_type_tmp])
# gauss result -------------
gauss_tmp = np.array([])
gauss_type_tmp = np.array([])
gauss_name_tmp = np.array([])
for gg in range(len(line_fit.params)):
gauss_tmp = np.concatenate([gauss_tmp, np.array([line_fit.params[gg]])])
if self.MC == True and self.n_trails > 0:
gauss_tmp = np.concatenate([gauss_tmp, np.array([all_para_std[gg]])])
gauss_result = np.concatenate([gauss_result, gauss_tmp])
# gauss result name -----------------
for n in range(nline_fit):
for nn in range(int(ngauss_fit[n])):
line_name = linelist['linename'][ind_line][n]+'_'+str(nn+1)
if self.MC == True and self.n_trails > 0:
gauss_type_tmp_tmp = ['float', 'float', 'float', 'float', 'float', 'float']
gauss_name_tmp_tmp = [line_name+'_scale', line_name+'_scale_err',
line_name+'_centerwave', line_name+'_centerwave_err',
line_name+'_sigma', line_name+'_sigma_err']
else:
gauss_type_tmp_tmp = ['float', 'float', 'float']
gauss_name_tmp_tmp = [line_name+'_scale', line_name+'_centerwave', line_name+'_sigma']
gauss_name_tmp = np.concatenate([gauss_name_tmp, gauss_name_tmp_tmp])
gauss_type_tmp = np.concatenate([gauss_type_tmp, gauss_type_tmp_tmp])
gauss_result_type = np.concatenate([gauss_result_type, gauss_type_tmp])
gauss_result_name = np.concatenate([gauss_result_name, gauss_name_tmp])
# further line parameters ----------
fur_result_tmp = np.array([])
fur_result_type_tmp = np.array([])
fur_result_name_tmp = np.array([])
fwhm, sigma, ew, peak, area = self.line_prop(compcenter, line_fit.params, 'broad')
br_name = uniq_linecomp_sort[ii]