博文
Linux系统下处理Chandra卫星的X射线数据ACIS（IDL程序）
|||
Pro acis_preparation_work, objID=objID, primary_path=primary_path, $
secondary_path=secondary_path
;NAME:
; acis_preparation_work
;PURPOSE:
; produce a batch job file for ACIS Data Preparation in CIAO
;CALLING SEQUENCE:
; acis_preparation
;INPUT:
; none
;OUTPUT:
; a batch job file
;REVISIONHISTORY:
; Original by DL.Wang,Apr-15-2008
;-
print,'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
print,'C Explanation:'
print,'C This code run under the dir secondary before run CIAO'
print,'C It is only used by ACIS !'
print,'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
if n_elements(objID) eq 0 then begin
print,'! Input Observation ID'
objID=''
read,objID
endif
if n_elements(primary_path) eq 0 then begin
print,'! Input primary path'
primary_path=''
read,primary_path
endif
if n_elements(secondary_path) eq 0 then begin
print,'! Input secondary path'
secondary_path=''
read,secondary_path
endif
;------------------------------
; Input level=1 event file
;------------------------------
print,'Input level=1 event file!'
evt1=''
read,evt1
head=headfits(secondary_path+evt1,exten=1)
readmode=sxpar(head,'READMODE')
print,readmode
datamode=sxpar(head,'DATAMODE')
print,datamode
tgain=sxpar(head,'TGAINCOR')
print,tgain
ascdsver=sxpar(head,'ASCDSVER')
print,'ASCDSVER:',ascdsver
print,'==================================================================='
print,' Attention!!:'
print,' if ASCDSVER is s lower than DS 7.6.7, it is necessary to run'
print,' Remove the acis_detect_afterglow Correction!'
print,'==================================================================='
;----------------------
; Check file
;----------------------
;bpix file
bpix=findfile(primary_path+'*_bpix1.fits')
;aso file
aso0=findfile(primary_path+'*_aso*.fits')
if n_elements(aso0) eq 1 then begin
aso=aso0
endif else begin
aso='"@aso.lis"'
openw,lun2,'aso.lis',/get_lun
for i=0L,n_elements(aso0)-1 do begin
printf,lun2,aso0[i]
endfor
endelse
;bias file
bias0=findfile(secondary_path+'*_bias*.fits')
if n_elements(bias0) eq 1 then begin
bias=bias0
endif else begin
bias='"@bias_file.lis"'
openw,lun1,'bias_file.lis',/get_lun
for i=0L,n_elements(bias0)-1 do begin
printf,lun1,bias0[i]
endfor
free_lun,lun1
endelse
;mask file
msk=findfile(secondary_path+'*_msk1.fits')
;pbk file
pbk=findfile(secondary_path+'*_pbk0.fits')
;flt file
flt=findfile(secondary_path+'*_flt1.fits')
;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
;CC
;CC Produce batch file for CIAO
;CC
;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
openw,lun,'acis_preparation'+objID,/get_lun
;###############################################
;### Examine the Afterglow Events
;###############################################
printf,lun,'dmkeypar '+secondary_path+evt1+' ASCDSVER echo+'
printf,lun,'dmcopy "'+secondary_path+evt1+'[exclude status=xxxxxxxxxxxx0000xxxxxxxxxxxxxxxx]" '+'acis_'+objID+'_afterglows_evt1.fits'
;###############################################
;### Reset the Status Bits
;###############################################
printf,lun,'punlearn dmtcalc'
printf,lun,'pset dmtcalc infile='+secondary_path+evt1
printf,lun,'pset dmtcalc outfile='+'acis_'+objID+'_reset_evt1.fits'
printf,lun,'pset dmtcalc expression="status=status,status=X15F,status=X14F,status=X13F,status=X12F"'
printf,lun,'dmtcalc'
;###############################################
;### Create a New Bad Pixel File
;###############################################
printf,lun,'dmkeypar '+secondary_path+evt1+' READMODE echo+'
printf,lun,'dmkeypar '+secondary_path+evt1+' DATAMODE echo+'
printf,lun,'punlearn acis_run_hotpix'
printf,lun,'pset acis_run_hotpix infile='+'acis_'+objID+'_reset_evt1.fits'
printf,lun,'pset acis_run_hotpix outfile='+'acis_'+objID+'_new_bpix1.fits'
printf,lun,'pset acis_run_hotpix badpixfile='+bpix
printf,lun,'pset acis_run_hotpix biasfile='+bias
printf,lun,'pset acis_run_hotpix maskfile='+msk
printf,lun,'pset acis_run_hotpix pbkfile='+pbk
printf,lun,'acis_run_hotpix'
;###############################################
;### Create a New Level=1 File
;###############################################
if readmode eq 'TIMED ' and datamode eq 'FAINT ' then begin
printf,lun,'punlearn acis_process_events'
printf,lun,'pset acis_process_events infile='+'acis_'+objID+'_reset_evt1.fits'
printf,lun,'pset acis_process_events outfile='+'acis_'+objID+'_new_evt1.fits'
printf,lun,'pset acis_process_events badpixfile='+'acis_'+objID+'_new_bpix1.fits'
printf,lun,'pset acis_process_events acaofffile='+aso
printf,lun,'pset acis_process_events eventdef=")stdlev1"'
printf,lun,'pset acis_process_events apply_tgain=yes'
printf,lun,'pset acis_process_events apply_cti=yes'
printf,lun,'pset acis_process_events rand_pix_size=0.5'
printf,lun,'pset acis_process_events rand_pha=yes'
printf,lun,'pset acis_process_events calc_cc_times=yes'
printf,lun,'acis_process_events'
endif
if readmode eq 'TIMED ' and datamode eq 'VFAINT ' then begin
printf,lun,'punlearn acis_process_events'
printf,lun,'pset acis_process_events infile='+'acis_'+objID+'_reset_evt1.fits'
printf,lun,'pset acis_process_events outfile='+'acis_'+objID+'_new_evt1.fits'
printf,lun,'pset acis_process_events badpixfile='+'acis_'+objID+'_new_bpix1.fits'
printf,lun,'pset acis_process_events acaofffile='+aso
printf,lun,'pset acis_process_events eventdef=")stdlev1"'
printf,lun,'pset acis_process_events apply_tgain=no'
printf,lun,'pset acis_process_events apply_cti=no'
printf,lun,'pset acis_process_events doevtgrade=no'
printf,lun,'pset acis_process_events check_vf_pha=yes'
printf,lun,'acis_process_events'
endif
;###############################################
;### Create a New Level=2 File
;###############################################
printf,lun,'punlearn dmcopy'
printf,lun,'dmcopy "'+'acis_'+objID+'_new_evt1.fits'+'[EVENTS][grade=0,2,3,4,6,status=0]" '+'acis_'+objID+'_flt_evt1.fits'
printf,lun,'dmcopy "'+'acis_'+objID+'_flt_evt1.fits'+'[EVENTS][@'+flt+'][cols -phas]" '+'acis_'+objID+'_evt2.fits'
printf,lun,'punlearn destreak'
printf,lun,'pset destreak infile=acis_'+objID+'_evt2.fits'
printf,lun,'pset destreak outfile=acis_'+objID+'_dstrk_evt2.fits'
printf,lun,'destreak'
;###############################################
;### acis_set_ardlib
;###############################################
printf,lun,'grep Id `which acis_set_ardlib`'
printf,lun,'punlearn ardlib'
printf,lun,'acis_set_ardlib '+'acis_'+objID+'_new_bpix1.fits'
free_lun,lun
End
;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Pro acis_extractspec, objID=objID, source_reg=source_reg, $
background_reg=background_reg
;+
;NAME:
; acis_extractspec
;PURPOSE:
; prouduce a batch file for extract spectrum used in ACIS for CIAO
;CALLING SEQUENCE:
; acis_extractspec
;INPUT:
; none
;OUTPUT:
; a batch file
;Optional keyword input:
; objID --- Observation ID
; source_reg --- the file name for source region (String)
; background_reg --- the file name for source background region
; (String)
;REVISION HISTORY:
; Original by DL.Wang,Mar-19-2008
;-
if n_elements(objID) eq 0 then begin
print,'! Input Observation ID'
objID=''
read,objID
endif
if n_elements(source_reg) eq 0 then begin
print,'! Input Source region file'
source_reg=''
read,source_reg
endif
if n_elements(background_reg) eq 0 then begin
print,'! Input Source background region file'
background_reg=''
read,background_reg
endif
;###############################################
;### Extracting Spectra
;###############################################
openw,lun,'acis_extractspec'+objID,/get_lun
printf,lun,'punlearn dmextract'
printf,lun,'pset dmextract infile="'+'acis_'+objID+'_dstrk_evt2.fits'+'[sky=region('+source_reg+')][bin pi]"'
printf,lun,'pset dmextract outfile='+objID+'.pi'
printf,lun,'pset dmextract wmap="det=8"'
printf,lun,'dmextract'
printf,lun,'pset dmextract infile="'+'acis_'+objID+'_dstrk_evt2.fits'+'[sky=region('+background_reg+')][bin pi]"'
printf,lun,'pset dmextract outfile='+objID+'bg.pi'
printf,lun,'dmextract'
printf,lun,'dmstat "'+'acis_'+objID+'_dstrk_evt2.fits'+'[sky=region('+source_reg+')]'+'[cols chipx,chipy,ccd_id,x,y]"'
printf,lun,'dmstat "'+'acis_'+objID+'_dstrk_evt2.fits'+'[sky=region('+background_reg+')]'+'[cols chipx,chipy,ccd_id,x,y]"'
free_lun,lun
End
;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Pro acis_rmfarf_make_work, objID=objID, primary_path=primary_path, $
secondary_path=secondary_path, ccd_id=ccd_id,$
pixelx_s=pixelx_s, pixely_s=pixely_s, $
pixelx_b=pixelx_b,pixely_b=pixely_b
;NAME:
; acis_rmfarf_make_work
;PURPOSE:
; produce a batch job file for make RMF and ARF for CIAO
;CALLING SEQUENCE:
; acis_rmfarf_make_work
;INPUT:
; none
;OUTPUT:
; a batch job file
;OPTIONAL KEYWORD INPUT:
; objID ---- Observation ID
; ccd_id ---- ID for CCD
; pixelx_s ---- source sky x pixel
; pixely_s ---- source sky y pixel
; pixelx_b ---- source background sky x pixel
; pixely_b ---- source background sky y pixel
;REVISION HISTORY:
; Original by DL.Wang,Mar-19-2008
;-
if n_elements(objID) eq 0 then begin
print,'! Input Observation ID'
objID=''
read,objID
endif
if n_elements(primary_path) eq 0 then begin
print,'! Input primary path'
primary_path=''
read,primary_path
endif
if n_elements(secondary_path) eq 0 then begin
print,'! Input secondary path'
secondary_path=''
read,secondary_path
endif
;----------------------
; Check file
;----------------------
;aso file
aso0=findfile(primary_path+'*_aso*.fits')
if n_elements(aso0) eq 0 then begin
print,'Please copy aso*.fits from primary dirtry!!'
return
endif
if n_elements(aso0) eq 1 then begin
aso=aso0
endif else begin
aso='"@aso.lis"'
openw,lun2,'aso.lis',/get_lun
for i=0L,n_elements(aso0)-1 do begin
printf,lun2,aso0[i]
endfor
endelse
;mask file
msk=findfile(secondary_path+'*_msk1.fits')
;pbk file
pbk=findfile(secondary_path+'*_pbk0.fits')
if n_elements(objID) eq 0 then begin
print,'! Input Observation ID'
objID=''
read,objID
endif
if n_elements(ccd_id) eq 0 then begin
print,'! Input CCD_ID'
ccd_id=''
read,ccd_id
endif
if n_elements(pixelx_s) eq 0 then begin
print,'! Input source sky x pixel'
pixelx_s=''
read,pixelx_s
endif
if n_elements(pixely_s) eq 0 then begin
print,'! Input source sky y pixel'
pixely_s=''
read,pixely_s
endif
if n_elements(pixelx_b) eq 0 then begin
print,'! Input source background sky x pixel'
pixelx_b=''
read,pixelx_b
endif
if n_elements(pixely_b) eq 0 then begin
print,'! Input source background sky y pixel'
pixely_b=''
read,pixely_b
endif
;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
;CC
;CC Produce batch file for CIAO
;CC
;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
openw,lun,'acis_rmfarf_make'+objID,/get_lun
;###############################################
;### Creating an RMF at a specific location
;###############################################
printf,lun,'punlearn mkacisrmf'
printf,lun,'pset mkacisrmf infile=CALDB'
printf,lun,'pset mkacisrmf outfile='+objID+'_mkacisrmf.rmf'
printf,lun,'pset mkacisrmf energy=0.1:11.0:0.01'
printf,lun,'pset mkacisrmf channel=1:1024:1'
printf,lun,'pset mkacisrmf chantype=PI'
printf,lun,'pset mkacisrmf wmap='+objID+'.pi'
printf,lun,'pset mkacisrmf gain=CALDB'
printf,lun,'mkacisrmf'
printf,lun,'pset mkacisrmf outfile='+objID+'bg_mkacisrmf.rmf'
printf,lun,'pset mkacisrmf wmap='+objID+'bg.pi'
printf,lun,'mkacisrmf'
;###############################################
;### Creating an ARF
;###############################################
printf,lun,'punlearn asphist'
printf,lun,'pset asphist infile='+aso
printf,lun,'pset asphist outfile='+objID+'.asphist'
printf,lun,'pset asphist evtfile="'+'acis_'+objID+'_dstrk_evt2.fits'+'[ccd_id='+ccd_id+']"'
printf,lun,'asphist'
printf,lun,'dmkeypar '+'acis_'+objID+'_dstrk_evt2.fits'+' GRATING echo+'
printf,lun,'dmkeypar '+'acis_'+objID+'_dstrk_evt2.fits'+' DETNAM echo+'
printf,lun,'punlearn mkarf'
printf,lun,'pset mkarf grating='+'NONE'
printf,lun,'pset mkarf detsubsys='+'ACIS-'+ccd_id
printf,lun,'pset mkarf outfile='+objID+'.arf'
printf,lun,'pset mkarf asphistfile="'+objID+'.asphist[ASPHIST]"'
printf,lun,'pset mkarf obsfile="'+'acis_'+objID+'_dstrk_evt2.fits'+'[EVENTS]"'
printf,lun,'pset mkarf maskfile='+msk
printf,lun,'pset mkarf pbkfile='+pbk
printf,lun,'pset mkarf engrid="grid('+objID+'_mkacisrmf.rmf[cols ENERG_LO,ENERG_HI])"'
printf,lun,'pset mkarf sourcepixelx='+pixelx_s
printf,lun,'pset mkarf sourcepixely='+pixely_s
printf,lun,'mkarf'
printf,lun,'pset mkarf outfile='+objID+'bg.arf'
printf,lun,'pset mkarf engrid="grid('+objID+'bg_mkacisrmf.rmf[cols ENERG_LO,ENERG_HI])"'
printf,lun,'pset mkarf sourcepixelx='+pixelx_b
printf,lun,'pset mkarf sourcepixely='+pixely_b
printf,lun,'mkarf'
free_lun,lun
End
;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Pro acis_lightcure_work, objID=objID, primary_path=primary_path, $
secondary_path=secondary_path, evt2=evt2, $
ccd_id=ccd_id, source_reg=source_reg, $
backgroud_reg=backgroud_reg
;+
;NAME:
; acis_lightcure
;PURPOSE:
; produce a batch job file for extract lightcure in ACIS for CIAO
;CALLING SEQUENCE:
; acis_lightcure
;INPUT:
; none
;OUTPUT:
; a batch job file
;REVISION HISTORY:
; Original by DL.Wang,Mar-19-2008
;-
if n_elements(objID) eq 0 then begin
print,'! Input Observation ID'
objID=''
read,objID
endif
if n_elements(primary_path) eq 0 then begin
print,'! Input primary path'
primary_path=''
read,primary_path
endif
if n_elements(secondary_path) eq 0 then begin
print,'! Input secondary path'
secondary_path=''
read,secondary_path
endif
if n_elements(evt2) eq 0 then begin
print,'Level=2 file is set by user? (y/n)'
select=''
read,select
case select of
'y':begin
evt2=''
read,evt2
end
'n':evt2='acis_'+objID+'_dstrk_evt2.fits'
endcase
endif
if n_elements(source_reg) eq 0 then begin
print,'Source region file is set by user? (y/n)'
select0=''
read,select0
case select0 of
'y':begin
source_reg=''
read,source_reg
end
'n':source_reg='source'+objID+'.reg'
endcase
endif
if n_elements(backgroud_reg) eq 0 then begin
print,'Source backgroud region file is set by user? (y/n)'
select1=''
read,select1
case select1 of
'y':begin
backgroud_reg=''
read,backgroud_reg
end
'n':backgroud_reg='source'+objID+'bg.reg'
endcase
endif
if n_elements(time_binsize) eq 0 then begin
print,'Time bin size is used by default? (y/n)'
select2=''
read,select2
case select2 of
'n':begin
time_binsize=''
read,time_binsize
end
'y':time_binsize='2000'
endcase
endif
eph=findfile(primary_path+'orbit*eph*.fits')
if n_elements(ccd_id) eq 0 then begin
print,'! Input CCD_ID'
ccd_id=''
read,ccd_id
endif
;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
;C
;C Produce CIAO batch job
;C
;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
openw,lun,'acis_lightcure'+objID,/get_lun
;##########################################
;## Apply Barycenter Correction
;##########################################
printf,lun,'dmkeypar '+evt2+' TSTART echo+'
printf,lun,'punlearn axbary'
printf,lun,'pset axbary infile='+evt2
printf,lun,'pset axbary orbitfile='+eph
printf,lun,'pset axbary outfile=acis_'+objID+'_bary_evt2.fits'
printf,lun,'axbary'
printf,lun,'dmlist '+evt2+'"[cols time]" data rows=1:5'
printf,lun,'dmlist acis_'+objID+'_bary_evt2.fits'+'"[cols time]" data rows=1:5'
;################################################
;## Determine which chips are being used
;################################################
printf,lun,'punlearn dmstat'
printf,lun,'dmstat "'+'acis_'+objID+'_bary_evt2.fits'+'[sky=region('+source_reg+')][cols ccd_id]"'
printf,lun,'dmstat "'+'acis_'+objID+'_bary_evt2.fits'+'[sky=region('+backgroud_reg+')][cols ccd_id]"'
;####################################################
;## Create a background-subtracted lightcurve
;####################################################
printf,lun,'punlearn dmextract'
printf,lun,'pset dmextract infile="'+'acis_'+objID+'_bary_evt2.fits'+'[ccd_id='+ccd_id+',sky=region('+source_reg+')][bin time=::'+time_binsize+']"'
printf,lun,'pset dmextract outfile=acis_'+objID+'_'+time_binsize+'_sub_lc.fits'
printf,lun,'pset dmextract bkg="'+'acis_'+objID+'_bary_evt2.fits'+'[ccd_id='+ccd_id+',sky=region('+backgroud_reg+')]"'
printf,lun,'pset dmextract opt="ltc1"'
printf,lun,'dmextract'
free_lun,lun
End
;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Pro acis_lightcure_plot
print,'!! Input acis light cure file'
file=''
read,file
data=mrdfits(file,1,h)
time_min=data.time_min
time_max=data.time_max
time=data.time
counts=data.counts
stat_err=data.stat_err
count_rate=data.count_rate
count_rate_err=data.count_rate_err
net_counts=data.net_counts
net_err=data.net_err
nonzero=where(counts ne 0.0,non0)
print,non0
plotcure:
print,'!! Plot select:'
print,' 0: counts'
print,' 1: count_rate'
print,' 2: net_counts'
select=''
read,select
case select of
'0':ploterror,time,counts,stat_err,psym=6,xtitle='Time(s)',ytitle='Count',charsize=1.2
'1':ploterror,time,count_rate,count_rate_err,psym=6,xtitle='Time(s)',ytitle='Count Rate',charsize=1.2
'2':ploterror,time,net_counts,net_err,psym=6,xtitle='Time(s)',ytitle='Net Count',charsize=1.2
endcase
print,'!! plot continue? (y/n)'
continsel:
contin=''
read,contin
case contin of
'y':goto,plotcure
'n':return
'':goto,continsel
endcase
End
;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Pro acis_batchjob
;+
;NAME:
; acis_batchjob
;PURPOSE:
; batch ACIS for CIAO
;INPUT:
;OUTPUT:
;HISTORY:
; Orignal by DL.Wang,Oct-15-2008
;-
spawn,'source data/software/chandra/ciao-4.0/bin/ciao.csh'
acis_preparation_work
spawn,'ls'
print,'!! Input acis_preparation_work batchjob file:'
batch1=''
read,batch1
spawn,'source '+batch1
cleanimage=''
read,cleanimage
spawn,'ds9 '+cleanimage+'&'
acis_extractspec
spawn,'ls'
print,'!! Input acis_extractspec batchjob file:'
batch2=''
read,batch2
spawn,'source '+batch2
acis_rmfarf_make_work
spawn,'ls'
print,'!! Input acis_rmfarf_make_work batchjob file:'
batch3=''
read,batch3
spawn,'source '+batch3
spawn,'ls'
acis_lightcure_work
print,'!! Input acis_lightcure_work batchjob file:'
batch4=''
read,batch4
spawn,'source '+batch4
spawn,'ls'
acis_lightcure_plot
End
https://blog.sciencenet.cn/blog-456360-658966.html
上一篇：计算光谱的4000埃break（IDL程序）
下一篇：如何检测Linux系统CPU数目
