pro acs_clock_charge, FILE_INPUT, FILE_OUTPUT, $ UNCLOCK=unclock, $ JAVA_DIR=java_dir, $ WELL_DEPTH=well_depth, $ WELL_NOTCH_DEPTH=well_notch_depth, $ WELL_FILL_POWER=well_fill_power, $ TRAP_DENSITY=trap_density, $ TRAP_LIFETIME=trap_lifetime, $ TRAP_LOCATION_SEED=trap_location_seed, $ TRAP_DECAY_SEED=trap_decay_seed, $ TRAP_INITIAL_DENSITY=trap_initial_density, $ TRAP_PERMANENT_RATE=trap_permanent_rate, $ TRAP_REMOVABLE_RATE=trap_removable_rate, $ N_ITERATIONS=n_iterations, $ KEEP_TEMP_FILES=keep_temp_files, $ PID=pid ;+ ; NAME: ; ACS_CLOCK_CHARGE ; ; CATEGORY: ; Pixel-by-pixel correction of trailing due to charge transfer inefficiency in CCD detectors. ; ; PURPOSE: ; Add or remove image trailing due to charge transfer inefficiency (CTI) in CCD detectors. ; To add trails, this simply acts as an IDL wrapper for Chris Stoughton's CCD readout code, ; which is written in Java. However, removal requires iteration within IDL. If a whole set ; of images need correcting, this routine can also farm out the operation to a cluster node ; and then return. A wrapper routine should therefore repeatedly call this one with new ; file names. ; ; NOTES: ; This spawns all sorts of UNIX child processes. It is very unlikely to ever work on Windows ; machines or anything like that. ; ; INPUTS: ; FILE_INPUT - String containing the absolute file name a (single) input FITS image. ; FILE_OUTPUT - String containing the absolute file name for the desired output image. ; ; KEYWORD PARAMETERS: ; UNCLOCK - Undo the CCD readout, rather than doing it. ; CLUSTER_FARM- Run, in the background, on a free processor from a cluster. ; NB: SSH public/private keys need to be set up so that it is possible to ; log into the target nodes without interatively entering a password. ; NB: Once the process starts, this routine returns, and IDL moves on ; to the next task. No reporting is done when the image is ready! ; Clashes between temporary file names should not occur, so long as input ; images are not processed twice simultaneously. ; KEEP_TEMP_FI- Do not delete temporary files an difference images (only relevant during ; unclocking). ; ; OPTIONAL INPUTS: ; JAVA_DIR - String containing the directory with the ClockCharge Java code. ; N_ITERATIONS- Number of iterations used if unclocking the CCD readout. ; Residuals of order N+1 will remain. ; WELL_DEPTH - Depth of CCD pixel [electrons]. ; WELL_NOTCH_D- Depth of supplementary buried channel in each pixel [electrons]. ; WELL_FILL_PO- Describes rate at which pixel well fills up above the notch. ; TRAP_DENSITY- Scalar or array containing 3D density of (each species of) charge trap. ; TRAP_LIFETIM- Scalar or array containing characteristic release timescale for ; (each species of) charge trap. This is the 1/e NOT half-life. ; TRAP_LOCATIO- Random seed used to initally scatter charge traps about the CCD. ; TRAP_DECAY_S- Random seed used to decide when traps release charge. ; ; OUTPUTS: ; A new FITS image is written to FILE_OUTPUT. ; ; OPTIONAL OUTPUTS: ; PID - The ID number of the child process managing a farmed-out job. Note that ; the returned value can be systematically lower than the actual process ID, ; due to any any extra processes started when a new shell is initialised. ; ; PROCEDURES USED: ; Various FITSIO routines from the IDL Astronomy Users' Library are required, including ; fits_read.pro, fits_write.pro and sxaddpar.pro ; ; MODIFICATION HISTORY: ; Mar 10 - V2.0 option to keep temporary files added by RM. ; Oct 09 - V1.1 JAVA_DIR correctly parsed for unclocking, and move of all cluster control code to acs_correct_cte completed by RM. ; Jul 09 - V1.0 used by RM for COSMOS v2.0 images and analysis in Massey et al. (2010). ; Mar 09 - V0.3 but with new parameters and allowing negative pixels by RM. ; Jul 08 - V0.3 used by RM for first pass of the COSMOS correction (which was flawed and never went public). ; Jun 08 - CLUSTER_LOCAL option added by RM. ; Mar 08 - CLUSTER options added by RM. ; Feb 08 - Written by Richard Massey. ;- ; Version numbering version="v2.0 - Mar 2010" ; Parse inputs if not keyword_set(file_input) then message,"Input filename not specified!" & if not file_test(file_input) then message,"Input file does not exist!" if not keyword_set(file_output) then message,"Output filename not specified!" & if file_test(file_output) then message,/INFO,"CAUTION: output file "+file_output+" will be overwritten!" ;if not keyword_set(file_residual) then file_residual=file_dirname(file_output,/mark_directory)+file_basename(file_output,".fits")+"_residual.fits" & if file_test(file_residual) then message,/INFO,"CAUTION: residual image "+file_output+" will be overwritten!" if n_elements(n_iterations) eq 0 then n_iterations=1 if not keyword_set(well_depth) then well_depth=84700; if not keyword_set(well_notch_depth) then well_notch_depth=96.5;99.3;93.0;94.08;94.63;98.65;96.85;97.73;90.68;91.34 if not keyword_set(well_fill_power) then well_fill_power=0.576;0.566;0.588;0.604;0.571;0.5636;0.575;0.563;0.628;0.5961;0.462;0.437;0.46 if not keyword_set(trap_location_seed) then trap_location_seed=12 if not keyword_set(trap_decay_seed) then trap_decay_seed=13 if not keyword_set(java_dir) then java_dir="~/bin/CTE/snapSim-current-20090328.230001/" ;if not keyword_set(java_dir) then java_dir="~/bin/CTE/snapSim/" if keyword_set(unclock) then begin ; Initiliase model image using the observed data message,/INFO,"Removing CTE trails from image "+file_input & pid=0 if keyword_set(keep_temp_files) then message,/INFO,"Keeping temporary files" fits_read,file_input,observed,header model=observed ; Iterate to improve the model for iteration=1,n_iterations do begin message,/INFO,"Iteration "+strtrim(iteration,2)+"/"+strtrim(n_iterations,2) sxaddpar,header,"CTE_ITER",iteration," Number of iterations used during CTE correction" if min(model) lt 0 then message,"Image contains negative pixels!",/INFO temp_input=file_dirname(file_input,/mark_directory)+file_basename(file_input,".fits")+"_model"+strtrim(iteration,2)+".fits" temp_output=file_dirname(file_input,/mark_directory)+file_basename(file_input,".fits")+"_modelread"+strtrim(iteration,2)+".fits" temp_residual=file_dirname(file_input,/mark_directory)+file_basename(file_input,".fits")+"_modelresidual"+strtrim(iteration,2)+".fits" fits_write,temp_input,float(model),header ; No longer truncating at zero acs_clock_charge, temp_input, temp_output, $ WELL_DEPTH=well_depth, $ WELL_NOTCH_DEPTH=well_notch_depth, $ WELL_FILL_POWER=well_fill_power, $ TRAP_DENSITY=trap_density, $ TRAP_INITIAL_DENSITY=trap_initial_density, $ TRAP_PERMANENT_RATE=trap_permanent_rate, $ TRAP_REMOVABLE_RATE=trap_removable_rate, $ TRAP_LIFETIME=trap_lifetime, $ TRAP_LOCATION_SEED=trap_location_seed, $ TRAP_DECAY_SEED=trap_decay_seed, $ JAVA_DIR=java_dir if not keyword_set(keep_temp_files) then file_delete,temp_input if file_test(temp_output) then begin fits_read,temp_output,model_readout,header residual=temporary(model_readout)-observed if not keyword_set(keep_temp_files) then file_delete,temp_output else fits_write,temp_residual,float(residual),header print,"Minmax(old image):",minmax(model) model=(temporary(model)-residual)<((2.^16)-3) ; No longer truncating at zero print,"Minmax(new image):",minmax(model) print,"Minmax(change):",minmax(-1*residual) delvarx,residual endif else message,"Java code did not produce an output file!" endfor ; Store best-fit image that, when read out using software, matches the observed data message,/INFO,"Saving image "+file_output fits_write,file_output,model,header endif else if keyword_set(express) then begin message,/info,"Using Jay Anderson's trick to speed up runtime by a factor of ~1000!" message,"To do" endif else begin ; Extract metaparameters needed for CTE model message,/INFO,"Adding CTE trails to image "+file_input & pid=0 fits_read,file_input,image,header if not keyword_set(xrange) then xrange=[1,(size(image,/DIMENSIONS))[0]]-1 if not keyword_set(yrange) then yrange=[1,(size(image,/DIMENSIONS))[1]]-1 delvarx,image date=sxpar(header,"DATE")+2452334.5 background=sxpar(header,"BACKGROUND") if not keyword_set(trap_lifetime) then trap_lifetime=[0.88,10.4] if not keyword_set(trap_density) then trap_density=acs_trap_density(date,TRAP_INITIAL_DENSITY=trap_initial_density,TRAP_PERMANENT_RATE=trap_permanent_rate,TRAP_REMOVABLE_RATE=trap_removable_rate)*[0.25,0.75] if n_elements(trap_density) ne n_elements(trap_lifetime) then message,"Trap properties incorrectly specified!" if n_elements(trap_density) gt 4 then message,/INFO,"Information about some species of charge trap will be missing from the FITS header!" if n_elements(trap_density) eq 1 then trap_density_string=strtrim(trap_density,2) else trap_density_string=strjoin(strtrim(trap_density,2)+[replicate(",",n_elements(trap_density)-1),""]) if n_elements(trap_lifetime) eq 1 then trap_lifetime_string=strtrim(trap_lifetime,2) else trap_lifetime_string=strjoin(strtrim(trap_lifetime,2)+[replicate(",",n_elements(trap_lifetime)-1),""]) command=java_dir+'eagRunner.sh '+java_dir+' -Xmx2048m gov.fnal.eag.sim.pixel.ClockCharge'+$ ' --gain 1 --skyPerPixel 0 --verbosity 1 --showStopwatch 1'+$ ' --quantizeTraps true'+$ ;' --quantizeCharge false --nLevels 3'+$ ; Used for COSMOS correction ' --quantizeCharge false --nLevels 2'+$ ; This can be increased if the trap density is particularly low ;' --quantizeCharge true --nLevels 10000'+$ ; Continuum version ' --serialTrapsPerPixel "0." --serialTrapLifetime "1."'+$ ' --trapLocationRandomSeed '+strtrim(trap_location_seed,2)+$ ' --trapDecayRandomSeed '+strtrim(trap_decay_seed,2)+$ ' --parallelFullWell '+strtrim(well_depth,2)+$ ' --parallelNotchDepth "'+strtrim(well_notch_depth,2)+'" '+$ ' --parallelPower "'+strtrim(well_fill_power,2)+'" '+$ ' --parallelTrapsPerPixel "'+trap_density_string+'"'+$ ' --parallelTrapLifetime "'+trap_lifetime_string+'" '+$ file_input+' '+file_output print & print,command & print ;file_copy,file_input,file_output ; A cheat! spawn,"time "+command fits_read,file_output,image,header ; Update header to include CTE model parameters message,/INFO,"Updating FITS header in output image "+file_output sxaddpar,header,"CTE_VERS",version," CTE model applied "+systime() sxaddpar,header,"CTE_JAVA",java_dir," Directory of Java ClockCharge code" sxaddpar,header,"CTE_WELD",well_depth," Assumed pixel well depth [electrons]" sxaddpar,header,"CTE_WELN",well_notch_depth," Assumed notch depth [electrons]" sxaddpar,header,"CTE_WELP",well_fill_power," Power law controlling filling of pixel well" sxaddpar,header,"CTE_NTRS",n_elements(trap_density)," Number of charge trap species" for i=0,n_elements(trap_density)-1 do sxaddpar,header,"CTE_TR"+strtrim(i+1,2)+"D", ([trap_density,replicate(0,3)])[i]," Density of "+strtrim(i+1,2)+(["st","nd","rd",replicate("th",n_elements(trap_density))])[i]+" charge trap species [pixel^-1]" for i=0,n_elements(trap_lifetime)-1 do sxaddpar,header,"CTE_TR"+strtrim(i+1,2)+"T",([trap_lifetime,replicate(0,3)])[i]," Decay halflife of "+strtrim(i+1,2)+(["st","nd","rd",replicate("th",n_elements(trap_density))])[i]+" charge trap species" sxaddpar,header,"CTE_RND1",trap_location_seed," Random seed for charge trap locations" sxaddpar,header,"CTE_RND2",trap_decay_seed," Random seed for charge trap decays" fits_write,file_output,image,header endelse end