注册 |登录

deliangwang的个人博客分享 http://blog.sciencenet.cn/u/deliangwang

博文

Linux系统下处理Chandra卫星的X射线数据ACIS（IDL程序）

已有 3606 次阅读 2013-2-3 00:40 |个人分类:编程笔记|系统分类:科研笔记| Linux, 程序, job, file, produce

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数目
收藏 IP: 123.86.145.*| 热度|

1 杨华磊

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-9-26 00:51

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部