1 """
2 Various aux functions.
3 There should be more bla bla at places like this ...
4
5 """
6
7 import scipy.signal as spsig
8 import numpy as np
9 import math
10 import sp, li, pl, tem, ca
11 import copy as pythoncopy
12 import matplotlib.pyplot as plt
13 import cPickle
14
15
16
18 """detects the n highest peaks in the array and returns their wavelength and their height in a list of two lists"""
19
20 kernel = 1
21 arr = spsig.medfilt(array, kernel)
22 size = len(arr)
23 min = np.min(arr)
24
25 m = []
26 peaksheight = []
27 i = 1
28 while i <= n:
29
30 index = np.argmax(arr)
31 peakheight = max(arr)
32 if index + 1 < size and index - 1 >= 0 and arr[index+1] != min and arr[index-1] != min:
33 i += 1
34 m.append(index+1)
35 peaksheight.append(peakheight)
36 delta1 = - 1
37 delta2 = - 1
38 j = 0
39 while delta2 < 1 and delta1 < 1:
40 j += 1
41
42 if index - j > 0:
43 delta1 = arr[index-j-1]-arr[index-j]
44 else:
45 delta1 = 1000
46 arr[index] = min
47 if index + j < size -1:
48 delta2 = arr[index+j+1]-arr[index+j]
49 else:
50 delta2 = 1000
51 arr[index] = min
52
53 arr[index-j:index+j+1] = min
54 if np.all(arr == min):
55 i = n + 1
56
57 m = np.array(m)
58 peaksheight = np.array(peaksheight)
59 return np.row_stack((m, peaksheight))
60
61
62
64 """calculates the average flux (with and without continum) of a spectra listand returns a spectra"""
65
66 m = []
67 specnumber = len(spectralist)
68 for j in range(specnumber):
69 m.append(spectralist[j].ztopix(spectralist[j].zsdss))
70
71 flux_cs, flux, noise = [], [], []
72 for i in range(len(spectralist[0].fluxes_cs)):
73 a, b, c = [], [], []
74 for j in range(len(spectralist)):
75 index = i+m[j]
76 if index < len(spectralist[j].fluxes_cs):
77 a.append(spectralist[j].fluxes_cs[index])
78 b.append(spectralist[j].fluxes[index])
79 c.append(spectralist[j].noise[i])
80 flux_cs.append(np.mean(a))
81 flux.append(np.mean(b))
82 noise.append(np.mean(c))
83
84 spec = sp.spectra(np.array(flux_cs),np.array( flux), np.array(noise))
85 name = 'mean spectra from ' + str(specnumber) + ' files'
86 spec.name = name
87 spec.setangst()
88 spec.calc_csc()
89 return spec
90
91
92 -def candidat_search(directory_in, directory_out, file, tuplelist, cand_id, minpeaks, minqual, minborder, maxborder, noisetuples = [], masktuples = [], fluxchoice = 'fluxes_csc', n = 10, ratiolimit = 0.5):
93 """from a fit file, it writes candidats files (if they satisfy the selection criteria). It returns the number of candidats created. This function is quite complicated, it calls many other functions !"""
94 cand_number = 0
95 spec = sp.spectra()
96 spec = spec.read(directory_in, file)
97
98 spec.initmask()
99 spec.mask_borders(minborder, maxborder)
100 spec.mask_angst(noisetuples)
101 spec.mask_angst(masktuples, spec.zsdss)
102
103 template = tem.factory(tuplelist)
104
105 ccsp = sp.crosscorrspec(spec, template)
106 ccsp.crosscorr(fluxchoice)
107
108 peaks = peak_detect(ccsp.corr, n)
109 print len(peaks[0]),
110 if peaks != []:
111 mainpeaks = mainpeaksdet(peaks, spec)
112 print len(mainpeaks[0]),
113 z = mainpeaks[0]
114
115 for i in range(len(mainpeaks[0])):
116 temp = template.copy()
117 temp.setglobalz(z[i])
118 temp.optimize(ccsp.spectra, fluxchoice)
119
120 cand_number += create_cand(spec, temp, z[i], directory_out, ratiolimit, cand_id, cand_number, minpeaks, minqual, minborder, maxborder)
121 return cand_number
122
123
124 -def create_cand(spec, temp, z, directory_out, ratiolimit, cand_id, cand_number, minpeaks, minqual, minborder, maxborder):
125 """This is more precisely the function that creats the candidats files, it applies the criteria of selection"""
126 ratio = sigtonoise(spec, temp)
127 number = 0
128
129 lines = []
130 quality = 0
131
132 for j in range(len(ratio)):
133 if ratio[j] > ratiolimit:
134 temp.linelist[j].strength = ratio[j]
135 lines.append(temp.linelist[j])
136 quality += ratio[j]
137
138 peaksnumber = len(lines)
139 if peaksnumber >= minpeaks and quality >= minqual:
140
141 number += 1
142
143 name = str(quality)[0:5] + '_' + str(peaksnumber).zfill(2) + '_' + str(cand_id+cand_number+number).zfill(5) + '.ca'
144 cand = ca.candidat(name, spec, lines, z, quality)
145 cand.write(directory_out)
146
147 return number
148
149
151 """it returns a list containing the ratio between the signal and the noise for the lines in the template"""
152 ratiolist = []
153
154 for i in range(len(tem.linelist)):
155
156 width = tem.linelist[i].width
157 centerangst = tem.linelist[i].angst
158
159 n = spectra.angsttopix(centerangst)
160 i1 = spectra.angsttopix(centerangst)-spectra.angsttopix(centerangst-width)
161 i2 = spectra.angsttopix(centerangst+width) - spectra.angsttopix(centerangst)
162
163 av_noise = np.mean(spectra.noise[n - i1:n + i2 + 1])
164 if av_noise != 0:
165 ratiolist.append(tem.linelist[i].height/av_noise)
166 else:
167 ratiolist.append(0)
168
169 return ratiolist
170
171 -def mainpeaksdet(peaks, sp):
172 """returns a list of the 'real' peaks in a signal, it checks that two peaks very close to each other are not considered as two different ones ; the minimum distance bewteen to peaks is arbitrary"""
173 z = sp.pixtoz(peaks[0])
174 heights = peaks[1]
175 mainpeaks = [[z[0]],[heights[0]]]
176
177 for i in range(len(z)-2):
178 deltax1 = np.absolute(z[i+1] - z[i])
179 if deltax1 > 0.03:
180 mainpeaks[0].append(z[i+1])
181 mainpeaks[1].append(heights[i+1])
182 return mainpeaks
183
185 """returns a list of candidats sorted by quality"""
186 quality = []
187 for file in filelist:
188 qual = file[len(directory):len(directory)+5]
189 quality.append(float(qual))
190 index = np.argsort(quality)
191 index = index[::-1]
192 filelist = np.array(filelist)
193 return filelist[index]
194
195
196
198 """filters a signal using a mean fliter ; N must be odd [not used]"""
199 n = (N-1)/2
200 size = len(arr)
201 a = []
202 for i in range(size):
203 if i < n :
204 a.append(np.mean(arr[0:i+n+1]))
205 elif i > size - n - 1 :
206 a.append(np.mean(arr[i-n:size]))
207 else:
208 a.append(np.mean(arr[i-n:i+n+1]))
209 return np.array(a)
210
211
213 """it is to check that the wavelengths contained in the tuples are in the range of the spectrometer"""
214 useful_tuples = []
215
216 for tuple in reftuples:
217 if tuple[1] * (1 + z) > minangst and tuple[1] * (1 + z) < maxangst:
218 useful_tuples.append(tuple)
219
220 return useful_tuples
221