pro MAKE_BG_FILES, aaaajjjd, hd, aaaajjjf, hf, $
quarter=quarter, nbgbin=nbgbin, $
keep_interval=keep_interval, $
output_path=output_path, $
verbose=verbose, nodb=nodb
message,'Building Freq_index list',/info
fi = LIST_FREQ(aaaajjjd, hd, aaaajjjf, hf, /freq_index, verbose=verbose)
message,'Building Freq_list from Freq_index',/info
nf = n_elements(fi)
xf = FI_FREQ(fi)
if keyword_set(quarter) then nbgbin=1600 else begin
if keyword_set(nbgbin) then begin
if nbgbin mod 80 ne 0 then message,'NBGbin must be a multiple of 80 !'
endif else nbgbin = 1600
endelse
h = lonarr(nf,nbgbin+1,4)
fon = fltarr(nf,4)
sig = fltarr(nf,4)
fon5 = fltarr(nf,4)
fon10 = fltarr(nf,4)
bt = fltarr(nf,4)
nbt = fltarr(nf,4)
message,'Loading N1 data file List',/info
ptrDataFileList = make_file_list(aaaajjjd, hd, aaaajjjf, hf, level='n1')
if (ptr_valid(ptrDataFileList) eq 0) then begin
nFileList_N1=0
if keyword_set(verbose) then message,"Warning ! No N1 data file selected...",/info
endif else begin
fileList_N1=*ptrDataFileList
nFileList_N1 = n_elements(FileList_N1)
ptr_free, ptrDataFileList
endelse
message,strcompress("found "+string(nFileList_N1)+" N1 data files."),/info
message,'Loading N2 data file List',/info
ptrDataFileList = make_file_list(aaaajjjd, hd, aaaajjjf, hf, level='n2')
if (ptr_valid(ptrDataFileList) eq 0) then begin
nFileList_N2=0
if keyword_set(verbose) then message,"Warning ! No N2 data file selected...",/info
endif else begin
fileList_N2=*ptrDataFileList
nFileList_N2 = n_elements(FileList_N2)
ptr_free, ptrDataFileList
endelse
message,strcompress("found "+string(nFileList_N2)+" N2 data files."),/info
if nFileList_N1 ne nFileList_N2 then begin
if nFileList_N1 gt nFileList_N2 then message,"MISSING N2 DATA FILES ! ABORTING !" $
else message,"MISSING N1 DATA FILES ! ABORTING !"
return
endif
for iFile = 0,nFileList_N1-1l do begin
message,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++",/info
if keyword_set(verbose) then message,fileList_N1(iFile),/info
READ_DATA_BINARY, fileList_N1(iFile),data_n1,level='n1'
message,fileList_N2(iFile),/info
READ_DATA_BINARY, fileList_N2(iFile),data_n2,level='n2'
if keyword_set(verbose) then message,'> Extracting frequency list',/info
t0 = systime(/seconds)
ffi = uniq_list(data_n1.fi,rev=iffi)
if keyword_set(verbose) then message,string(format='(" ",f6.3," seconds elapsed")',systime(/seconds)-t0),/info
nffi = n_elements(ffi)
message,"> Looping on "+string(nffi)+" freq_index",/info
t0 = systime(/seconds)
for i=0,nffi-1 do begin
wfi = where(fi eq ffi(i))
if wfi(0) eq -1 then message,'This frequency index should not exist: '+string(ffi(i))
dataf = data_N2(iffi(iffi(i):iffi(i+1)-1))
for j=1,3 do begin
w=where((dataf.ant mod 10) eq j and dataf.autoX ne 0.)
if w(0) ne -1 then begin
x=(fix((10.*alog10(dataf(w).autoX)+170)*nbgbin/80) < nbgbin) > 0
for k=0L,n_elements(x)-1L do h(wfi,x(k),j)=h(wfi,x(k),j)+1L
bt(wfi,j)=bt(wfi,j)+ total(dataf(w).df*dataf(w).dt)
nbt(wfi,j)=nbt(wfi,j)+ n_elements(w)
endif
endfor
j=0
w=where(dataf.autoZ ne 0.)
if w(0) ne -1 then begin
x=(fix((10.*alog10(dataf(w).autoZ)+170)*nbgbin/80) < nbgbin) > 0
for k=0L,n_elements(x)-1L do h(wfi,x(k),j)=h(wfi,x(k),j)+1L
bt(wfi,j)=bt(wfi,j)+ total(dataf(w).df*dataf(w).dt)
nbt(wfi,j)=nbt(wfi,j)+ n_elements(w)
endif
endfor
if keyword_set(verbose) then message,string(format='(" ",f8.3," seconds elapsed")',systime(/seconds)-t0),/info
print,fileList_N2(iFile), max(h), total(h), nffi
endfor
message,"Analysing histograms",/info
t0 = systime(/seconds)
for i=0,nf-1 do begin
if keyword_set(verbose) then message,'Analysing F='+string(xf(i))+' [FI='+string(fi(i))+']',/info
for j=0,3 do begin
if total(h(i,*,j)) gt 1 then begin
x=[0.]
for k=0,nbgbin do if h(i,k,j) ne 0 then x=[x,(fltarr(h(i,k,j))+k+0.5)*80./nbgbin]
x=x(1:*)
FOND, x, f,s
fon(i,j)=10.^((f-170.)/10.)
sig(i,j)=10.^((f+s-170.)/10.)-fon(i,j)
DYN, x, 0.05,1,f,buf
fon5(i,j)=10.^((f-170.)/10.)
DYN, x, 0.10,1,f,buf
fon10(i,j)=10.^((f-170.)/10.)
endif
endfor
endfor
bt=bt/(nbt > 1)
if keyword_set(verbose) then message,string(format='(" ",f6.3," seconds elapsed")',systime(/seconds)-t0),/info
message,"Saving data_bg",/info
case nbgbin of
320 : data_bg = replicate({data_bg_320},nf)
640 : data_bg = replicate({data_bg_640},nf)
800 : data_bg = replicate({data_bg_800},nf)
1600 : data_bg = replicate({data_bg_1600},nf)
3200 : data_bg = replicate({data_bg_3200},nf)
endcase
data_bg.h = transpose(h,[1,2,0])
data_bg.bt = transpose(bt)
data_bg.nbt = transpose(nbt)
data_bg.xf = xf
data_bg.fi = fi
data_bg.sig = transpose(sig)
data_bg.fon = transpose(fon)
data_bg.fon5 = transpose(fon5)
data_bg.fon10 = transpose(fon10)
if not keyword_set(output_path) then begin
period_path = strmid(fileList_N1(0),0,strlen(fileList_N1(0))-14)
output_path = period_path+'bg/'
endif
if keyword_set(quarter) then begin
pos1=strpos(period_path, '/', /reverse_search)
pos2=strpos(strmid(period_path, 0, pos1-1), '/', /reverse_search)
output_file = output_path+'bg_'+strmid(period_path,pos2+1, pos1-pos2-1)
endif else begin
if keyword_set(keep_interval) then begin
str_ajd = string(format='(I7)',aaaajjjd)
str_hd = string(format='(I2.2)',hd)
str_ajf = string(format='(I7)',aaaajjjf)
str_hf = string(format='(I2.2)',hf)
endif else begin
str_ajd = strmid(FileList_N1(0),strlen(fileList_N1(0))-10,7)
str_hd = string(format='(I2.2)',fix(strmid(FileList_N1(0),strlen(fileList_N1(0))-2)))
str_ajf = strmid(FileList_N1(nFileList_N1-1),strlen(fileList_N1(nFileList_N1-1))-10,7)
str_hf = string(format='(I2.2)',fix(strmid(FileList_N1(nFileList_N1-1),strlen(fileList_N1(nFileList_N1-1))-2))+1)
endelse
output_file = output_path+strjoin(['bg',str_ajd,str_hd,str_ajf,str_hf],'_')+'__'+strtrim(string(nbgbin),2)
endelse
print,'Writing ouput data in :'
print,output_file+' ...'
write_data_binary,output_file,data_bg
print,'done.'
if ~keyword_set(nodb) then begin
upsert_args = getenv('ROOT_RPWS') + '/pro/kronosdb/upsert.sh'
if keyword_set(verbose) then upsert_args = [upsert_args, '-v']
upsert_args = [upsert_args, '-f', output_file]
spawn, /NOSHELL, upsert_args
end
return
end