PRO pt, chip ;get interactive parameters if not entered if n_elements(chip) ne 1 then begin chip=" " read, 'ENTER THE CHIP NAME :',chip endif filename=" " chipname=chip ;-------------------------------------- - end parameter ;definisce e apre il file log filelog='pt_'+chip+'.logo' openw,1,filelog ;controllo dei file nel direttorio e tempi di esposizione ;---------------------------- listafi=findfile('*.hhh') ;read all the .FITS files in the dir n_file=n_elements(listafi) ;count them if n_file le 0 then print,"THERE ARE NO FILES .HHH IN THIS DIRECTORY: " if (fix(n_file /2) - n_file/2.0) ne 0 then begin print,'!!!! WARNING: you should have 2 images for each exposure time:' print,'check for missing files:' print,listafile endif ;-------------------------------------------------------------------------- exptime=fltarr(n_file) ;define a 1-D vector with the exposure time print,'...please,wait: I am checking the exposure time' for h=0,n_file-1 do begin ;cicle on images strd,ima,he,listafi(h) stime=sxpar(he,'SAMPTIME') print,listafi(h)+" sec "+ strtrim(string(stime),2) printf,1,listafi(h)+" sec "+ strtrim(string(stime),2) exptime(h)=stime endfor onlytime=exptime(uniq(exptime,sort(exptime))) ;contains only exptime values print,"your images have the following exposure time (sec)" print,onlytime ;----------------------- fine lettura tempi esposizione ;definizione parametri per rejection------------------------- print,"define a threshold in stddev unit to reject deviant pixel (3-5)" read,nsig if nsig eq 0 then nsig=4 print,"how many iterations (2-4)" read,niter if niter eq 0 then niter=2 ; ----------------------------- fine definizione parametri rejection ;display the image in the middle of the list to the the user to chose ;the boxes for the statistics: sizbox=40 ;define the length of the box-side to use for the statistics filename=listafi(n_file/2) strd,ima,hea,filename ;read the image xdim=sxpar(hea,'NAXIS1') ;get the x dim from the header ydim=sxpar(hea,'NAXIS2') ;get the y dim from the header IF xdim LT 1550 THEN bfa=2 ELSE bfa=4 ;set the rebin factor to 4 print,xdim/bfa,ydim/bfa ; create the window and plot the image window,0,xsize=xdim/bfa,ysize=ydim/bfa,title=filename wset,0 & tvscl,bytscl(congrid(ima,xdim/bfa,ydim/bfa),min=(median(ima)-median(ima)/8.),max=(median(ima)+median(ima)/8.)) ;define and plot a frame x1=0+sizbox/bfa y1=0+sizbox/bfa x2=(xdim-sizbox)/bfa y2=(ydim-sizbox)/bfa plots,[x1,x2],[y1,y1],/dev plots,[x2,y2],/continue,/dev plots,[x1,y2],/continue,/dev plots,[x1,y1],/continue,/dev ; get the boxes print," " print,"Mark the center OF A BOX within the marked frame - TO EXIT CLICK THE RIGHT BUTTON TWICE" print," " !MOUSE.button=2 ax=intarr(10000) ay=intarr(10000) ai=0 ;index of marked points WHILE (!MOUSE.button NE 4) DO BEGIN ;Repeat until the right button is pressed (twice). cursor,x,y,/device,/down ;get the coordinates xx1=x-(sizbox/(2*bfa)) xx2=x+(sizbox/(2*bfa)) yy1=y-(sizbox/(2*bfa)) yy2=y+(sizbox/(2*bfa)) plots,[xx1,xx2],[yy1,yy1],/dev plots,[xx2,yy2],/continue,/dev plots,[xx1,yy2],/continue,/dev plots,[xx1,yy1],/continue,/dev ;riconverte i centri alla misura originale del chip ax(ai)=x*bfa ay(ai)=y*bfa ai=ai+1 ENDWHILE ax=ax(where(ax ne 0)) ;keep only the marked boxes ay=ay(where(ax ne 0)) ;--------------------------------------- fine acquisizione quadratini ma=fltarr(4,n_elements(ax),n_elements(onlytime)) posi=fltarr(2,n_elements(ax)) posi(0,*)=ax posi(1,*)=ay for gh=0,n_elements(onlytime)-1 do begin ;cycle on exptime ma(0:1,*,gh)=posi filename=listafi(where(exptime eq onlytime(gh))) ;controlla i due file con stessa exptime strd,ima,h,filename(0) ;apre il primo file strd,imb,h,filename(1) ;apre il secondo file print,filename(0)+" "+ filename(1)+"exptime (sec)"+string(onlytime(gh)) sz=size(ima) ;------------------------------------------------------------------------------------ overa=fltarr(3) staf1=moment(ima(2:15,30:sz(2)-30)) ; overscan f1 staf2=moment(ima(sz(1)-18:sz(1)-2,30:sz(2)-30)) ;overscan f2 stav1=moment(ima(40:sz(1)-40,sz(2)-18:sz(2)-2)) ;virtual ov if abs(staf1(0)-staf2(0))/staf1(0) gt 8 then print,"WARNING BIAS LEVEL" overa=[staf1(0),staf2(0),stav1(0)] ;values from overscan v1,v2,f1 rnsa=sqrt([staf1(1),staf2(1),stav1(1)]) ;readnoise in DN from the 3 overscan print,"image _1: physical overscan values: "+strtrim(staf1(0),2)+" / "+strtrim(staf2(0),2)+" virtual overscan: "+strtrim(stav1(0),2) printf,1,"image _1: physical overscan values: "+strtrim(staf1(0),2)+" / "+strtrim(staf2(0),2)+" virtual overscan: "+strtrim(stav1(0),2) printf,1,"------------------------------------------------------------------------------------" overb=fltarr(3) stbf1=moment(imb(2:15,30:sz(2)-30)) ; overscan v1 stbf2=moment(imb(sz(1)-18:sz(1)-2,30:sz(2)-30)) ;overscan v2 stbv1=moment(imb(40:sz(1)-40,sz(2)-18:sz(2)-2)) if abs(stbf1(0)-stbf2(0))/stbf1(0) gt 8 then print,"WARNING BIAS LEVEL" overb=[stbf1(0),stbf2(0),stbv1(0)] ;values from overscan v1,v2,f1 rnsb=sqrt([stbf1(1),stbf2(1),stbv1(1)]) ;readnoise in DN from the 3 overscan print,"image _2: physical overscan values: "+strtrim(stbf1(0),2)+" / "+strtrim(stbf2(0),2)+" virtual overscan: "+strtrim(stbv1(0),2) printf,1,"image _2: physical overscan values: "+strtrim(stbf1(0),2)+" / "+strtrim(stbf2(0),2)+" virtual overscan: "+strtrim(stbv1(0),2) printf,1,"------------------------------------------------------------------------------------" !p.multi=[0,2,2] plot,rebin(ima(2:15,30:sz(2)-30),1,sz(2)-60+1),yr=[staf1(0)-1.5*rnsa(0),staf1(0)+1.5*rnsa(0)],title="-A- physical overscan" oplot,rebin(ima(sz(1)-18:sz(1)-2,30:sz(2)-30),1,sz(2)-60+1),color=80 plot,rebin(ima(40:sz(1)-40,sz(2)-18:sz(2)-2),sz(1)-80+1,1),yr=[stav1(0)-1.5*rnsa(2),stav1(0)+1.5*rnsa(2)],title="-A- virtual ov." plot,rebin(imb(2:15,30:sz(2)-30),1,sz(2)-60+1),yr=[stbf1(0)-1.5*rnsb(0),stbf1(0)+1.5*rnsb(0)],title="-B- physical overscan" oplot,rebin(imb(sz(1)-18:sz(1)-2,30:sz(2)-30),1,sz(2)-60+1),color=80 plot,rebin(imb(40:sz(1)-40,sz(2)-18:sz(2)-2),sz(1)-80+1,1),yr=[stbv1(0)-1.5*rnsb(2),stbv1(0)+1.5*rnsb(2)],title="-B- virtual ov." !p.multi=0 ;------------------------------------- bias subtraction ima=ima-median(ima(40:sz(1)-40,sz(2)-18:sz(2)-2)) imb=imb-median(imb(40:sz(1)-40,sz(2)-18:sz(2)-2)) print,".. performing bias subtraction with virtual overscan median value" ;----------------------------------------------------------------- ;------------------------------------------------------- print,"... and statistics on the average and on the difference of the 2 images" imd=ima-imb ;create the image from the difference for the variance imm=(ima+imb)/2. ;create the average of the image for the signal for k=0,n_elements(posi(1,*))-1 do begin ;ciclo sui quadratini sezd=imd(posi(0,k)-20:posi(0,k)+20,posi(1,k)-20:posi(1,k)+20) sezm=imm(posi(0,k)-20:posi(0,k)+20,posi(1,k)-20:posi(1,k)+20) for vf=0,niter-1 do BEGIN ;ciclo sul nomero di iterazioni per la immagine sottratta stad1=moment(sezd) vdmin=median(sezd)-nsig*sqrt(stad1(1)) vdmax=median(sezd)+nsig*sqrt(stad1(1)) sezdtmp=sezd if min(sezd) lt vdmin or max(sezd) gt vdmax then begin sezdtmp(where(sezdtmp lt vdmin or sezdtmp gt vdmax))=median(sezd) endif sezd=sezdtmp ENDFOR ;fine ciclo interazioni stad2=moment(sezd) ma(3,k,gh)=stad2(1)/2. ;write the vriance for vf=0,niter-1 do BEGIN ;ciclo sul nomero di iterazioni per la immagine mediata stam1=moment(sezm) vmmin=median(sezm)-nsig*sqrt(stam1(1)) vmmax=median(sezm)+nsig*sqrt(stam1(1)) sezmtmp=sezm if min(sezm) lt vmmin or max(sezm) gt vmmax then begin sezmtmp(where(sezmtmp lt vmmin or sezmtmp gt vmmax))=median(sezm) endif sezm=sezmtmp endfor stam2=moment(sezmtmp) ma(2,k,gh)=stam2(0) ;write the mean signal endfor ;fine ciclo quadratini print,ma(0:3,*,gh) endfor ;end cycle on exptime ;------------------------------------------ gain calculation ga=fltarr(2,n_elements(posi(0,*))) ; definisce un vettore per i gain for l=0,n_elements(ga(0,*))-1 do begin gafit=poly_fit(ma(2,l,*),ma(3,l,*),1,yfi,yba,ysig) ga(0,l)=1.0/gafit(1) ga(1,l)=ysig endfor ;-------------------------------------plot section ;--------- 1 fg=fltarr(2,n_elements(onlytime)) ; calculate the mean signal and variance for each image for pk=0,n_elements(onlytime)-1 do begin fg(0:1,pk)=[median(ma(2,*,pk)),median(ma(3,*,pk))] endfor nlims=n_elements(onlytime)-1 ;plot all the values for all the boxes plot,ma(2,*,*),ma(3,*,*),psym=3,yra=[0,median(ma(2,*,nlims))+1000],xtit="SIGNAL (DN)",ytitl="VARIANCE" ;overplot the mean values and fit the slope gw=findgen(16)*(!pi*2/16.) ;define a new symbol usersym,cos(gw),sin(gw),/fill oplot,fg(0,*),fg(1,*),psym=8 curva2=poly_fit(fg(0,*),fg(1,*),1,yfi) ;fit the slope on mean points oplot,fg(0,*),yfi xyouts,fg(0,0)+1000,fg(1,nlims),"gain= "+strtrim(1/curva2(1),2)+" e/DN" print,"y=mx+q"+" with m="+curva2(1)+" and q="+curva2(0) print," gain from mean points: "+1/curva2(1)+" e/DN" print,'FIGURE 1: this figure shows all the measurements ' print,' for each exposure time (points), the ' print,' mean values (black circles) and the fit' print,' to the mean values' ;save the plot in a file answ1="" print,"do you want to create a postscript file (yY/nN)?" & read,answ1 if answ1 eq "y" or answ1 eq "Y" then begin set_plot,'ps device,filename="gain_slope.ps" plot,ma(2,*,*),ma(3,*,*),psym=3,yra=[0,median(ma(2,*,nlims))+1000],xtit="SIGNAL (DN)",ytitl="VARIANCE" oplot,fg(0,*),fg(1,*),psym=8 oplot,fg(0,*),yfi xyouts,fg(0,0)+1000,fg(1,nlims),"gain= "+strtrim(1/curva2(1),2)+" e/DN" device,/close set_plot,'x' print," the file gain_slope.ps has been created" endif ;-------2 ;histogram of gain values mngain=(fix((1.0/curva2(1)*100)))/100.0 xva=mngain/2+indgen(40)*0.05 plot,xva,histogram(ga(0,*),min=min(xva),max=max(xva),binsize=(max(xva)-min(xva))/n_elements(xva)),psym=10,xtitle="gain e/DN",ytitle="#" print,'FIGURE 2: this figure shows the histogram ' print,' for the gain estimate in each ' print,' section' ;save the plot in a file answ2="" print,"do you want to create a postscript file (yY/nN)?" & read,answ2 if answ2 eq "y" or answ2 eq "Y" then begin set_plot,'ps device,filename="gain_hist.ps" ;plot,xva,histogram(ga(0,*),min=min(xva),max=max(xva),binsize=(max(xva)-min(xva))/n_elements(xva)),psym=10,xtitle="gain e/DN",ytitle="#" plot,xva,histogram(ga(0,*),min=min(xva),max=max(xva),binsize=(max(xva)-min(xva))/n_elements(xva)),psym=10,xtitle="gain e/DN",ytitle="#" device,/close set_plot,'x' print," the file gain_hist.ps has been created" endif print,'NOTE: if the result is weird and the plot show points ' print,' with unusual variance check the images and run again' print,' this script using a lower value for the sigma-threshold' print,' and/or an higher number of iterations' print,' ' print,' you just used niter= '+strtrim(niter,2)+' and threshold= '+strtrim(nsig,2) printf,1,'you just used niter= '+strtrim(niter,2)+' and threshold= '+strtrim(nsig,2) printf,1,"GAIN [e/DN] =" + strtrim(1/curva2(1),2) close,1 END