generate_ramp.pro (JWST Pipeline routines in IDL)

; docformat = 'IDL'
;+
; Generates a NIRCam
;-    
pro GENERATE_RAMP,time,ramp,SEED=SEED
;
;Procedure to generate a ramp with parameters specified internally; 
;SEED is the montecarlo seed that may be passed in input
;
;The variables time,ramp are transmitted in output
;
if KEYWORD_SET(SEED) EQ 0 THEN SEED = -6L
;
;make the ramp parameters available to the public
common ramp_params,N,M,Ms,G,tf,tg,RON,rate
;
;set the ramp parameters
N = 10. ;nr of groups
M = FLOAT(4.) ;frames averaged in group
Ms = 0. ;frames skipped between groups
G = 0;1. ;digitization
tf = 10. ;seconds
tg = tf*(M+Ms); seconds
time = indgen(N)*tf*M+tf*(M+1)/2.;timing of the samples
time = DOUBLE(time)
;
;set the count rate, per second
rate = 9.;   if you want it random, add  +(RANDOMU(SEED)+1)*10  ;counts/second
;

;generate the counts within the intervals between frames, Delta_t => tf, according to Poisson stats.
counts_group = RANDOMN(seed,N*(M+Ms),POISSON=rate*tf)
;
;generate the ramp coadding the counts 
ramp_all = lonarr(N*(M+Ms))
ramp_all[0] = counts_group[0]
for i=0,N*(M+Ms)-1 do ramp_all[i]=ramp_all[i-1]+counts_group[i]; MEAN([counts_group[i*M:(i+1)*M-1]])
;
;add zero point/offset
offset=10000L ;offset counts
ramp_all=ramp_all+offset
;
;add readout noise, in output, per single read. 
RON=15.  ;readout noise
ramp_all+=RANDOMN(seed,N*(M+Ms))*RON
;
;group-average of the ramp, as we often do for JWST
ramp = fltarr(N)
for i=0,N-1 do ramp[i]+=MEAN([ramp_all[i*M:(i+1)*M-1]])
;
;I should add the rounding error, since we average 16 bit values and still store the results in 16 bit
;neglect for the moment...
;this should be it! we have time and ramp as output parameters
end

;+
;
;-
pro add_cosmic_rays,CRMAP, MODE=CRmode, CRPATH=crpath, SEED=SEED
;
;if KEYWORD_SET(SEED) EQ 0 THEN SEED=long(systime() )
;
;print,seed
;read the CR datacube (you will probably need to modify this path)
CRs = lonarr(21,21,10000)
for i=0,9 do begin   & $
  filecounter = string(i,FORMAT='(I02)') & $
  fits_read,'~/FUNCTIONAL/JWST/NIRCAM/Studies/2010/2010.CRMODELS/FITSout/CRs_MCT2.5_SUNMIN_'+filecounter+'_IPC.fits',CRcube   & $
  indexes = i*1000+indgen(1000) & $
  CRs[*,*,indexes] = CRcube   & $
endfor  
;pick up a number of CRs;
;we assume a Poisson distribution with average 706 events (from the TR, page 7)
NCrs = (FIX(randomn(systime_seed,1,POISSON=706)))[0]
;select NCrs from the stack
iCRs = FIX(randomu(systime_seed,NCRs)*1E4)
;
the_seed = SEED
;build CR map
CRmap = fltarr(2048+20,2048+20) ;add 20 pixels edge for convenience
x=LONG(randomu(the_seed,NCRs)*2048.+10)
y=LONG(randomu(the_seed+1,NCRs)*2048.+10)
for i= 0, NCRs-1 do CRmap[x[i]-10:x[i]+10,y[i]-10:y[i]+10]+=CRs[*,*,iCRs[i]] 
CRmap = CRmap[10:2047+10,10:2047+10]
CRmap[0:4,*] = 0
CRmap[*,0:4] = 0
CRmap[2044:*,*] = 0
CRmap[*,2044:*] = 0
end


; docformat = 'IDL'
;+
; Generates a NIRCam
;-    
pro GENERATE_DATACUBE,time,ramp,SEED=SEED, Ngroups=N
;
;Procedure to generate a datacube with parameters specified internally; 
;SEED is the montecarlo seed that may be passed in input
;
;The variables time,ramp are transmitted in output
;
if KEYWORD_SET(SEED) EQ 0 THEN SEED=long(systime(/seconds) )
;
;make the ramp parameters available to the public
;common ramp_params,N,M,Ms,G,tf,tg,RON,rate
;
;set the ramp parameters
IF N_ELEMENTS(Ngroups) EQ 0 THEN N = 10. ;nr of groups
M = FLOAT(4.) ;frames averaged in group
Ms = 0. ;frames skipped between groups
G = 0;1. ;digitization
tf = 10. ;seconds
tg = tf*(M+Ms); seconds
time = indgen(N)*tf*M+tf*(M+1)/2.;timing of the samples
time = DOUBLE(time)
;
;set the count rate, per second
rate = 9.;   if you want it random, add  +(RANDOMU(SEED)+1)*10  ;counts/second
;

;generate the electrons within the intervals between frames, Delta_t => tf, according to Poisson stats.
counts_group = RANDOMN(SEED,2048,2048,N*(M+Ms),POISSON=rate*tf)
;
CR_cube=lonarr(2048,2048,N*(M+Ms))
;
;add cosmic rays
for ix=0,N*(M+Ms)-1 do begin &
   add_cosmic_rays,CRMAP,SEED=ix
   counts_group[*,*,ix]+=LONG(CRMAP)   ; CRMAP is in electrons
   CR_cube[*,*,ix] = LONG(CRMAP)
endfor
;
;generate the ramp coadding the counts 
ramp_all = lonarr(2048,2048,N*(M+Ms))
ramp_all[0] = counts_group[0]
for i=0,N*(M+Ms)-1 do ramp_all[*,*,i]=ramp_all[*,*,i-1]+counts_group[*,*,i]; MEAN([counts_group[i*M:(i+1)*M-1]])
;
;add zero point/offset
offset=10000L ;offset counts
ramp_all=ramp_all+offset
;
;add readout noise, in output, per single read. 
RON=13.  ;readout noise
ramp_all+=RANDOMN(seed,2048,2048,N*(M+Ms))*RON
;
;group-average of the ramp, as we often do for JWST
ramp = lonarr(2048,2048,N)
CR_ramp = lonarr(2048,2048,N)
for i=0,N-1 do begin
    ramp[*,*,i]+=MEAN([ramp_all[*,*,i*M:(i+1)*M-1]],DIMENSION=3)
    CR_ramp[*,*,i]+=MEAN([CR_cube[*,*,i*M:(i+1)*M-1]],DIMENSION=3)
endfor
;
;do the CR ramp coadding them
for i=1,N-1 do CR_ramp[*,*,i]+=CR_ramp[*,*,i-1]

;
;I should add the rounding error, since we average 16 bit values and still store the results in 16 bit
;neglect for the moment...
;this should be it! we have time and ramp as output parameters

;ramp = reform(ramp,2048.*2048.,10)
;r=reform(ramp[ (where(ramp[*,9] GT 20000))[1004] ,*])
;r=r<1.E5
;r = r(WHERE(r LT 1.E5)) & time = time(WHERE(r LT 1.E5))
;print,(max(r)-min(r))/(N_ELEMENTS(r)-1),median((shift(r,-1)-r))
writefits,'~/Desktop/ramp.fits',ramp
stop
writefits,'~/Desktop/CR_ramp.fits',CR_ramp
end

pro FIND_CR,CR_found
common ramp_params,N,M,Ms,G,tf,tg,RON,rate
;
;set the ramp parameters
N = 10. ;nr of groups
M = FLOAT(4.) ;frames averaged in group
Ms = 0. ;frames skipped between groups
G = 0;1. ;digitization
tf = 10. ;seconds
tg = tf*(M+Ms); seconds
time = indgen(N)*tf*M+tf*(M+1)/2.;timing of the samples
time = DOUBLE(time)
;

fits_read,'~/Desktop/CR_ramp.fits',CR_ramp
;CR_ramp = reform(CR_ramp,2048.,2048.,10)
;ix_CR= where(CR_ramp[*,9] GT 0,NCrs)
fits_read,'~/Desktop/ramp.fits',ramp
RIGHT=0L
WRONG=0L
MISSED=0L
;for icr =0, NCRs-1 do begin
ramp = reform(ramp,2048.*2048.,10)
                            ramp=FLOAT(ramp)
 delta_time=[t[0],t[1:*]-t[0:N_ELEMENTS(r)-1]]
                           delta_time = transpose(replicate(1,2048*2048.)##delta_time)
 sigma=sqrt(first_diff+6.5^2)
 rate=first_diff/delta_Time
                           ix=WHERE(ramp GT 1.E5) 
                           ramp(ix) = !VALUES.F_NAN 
 resistant_mean,rate,3,mm,smm,goodvec = gv
 RESISTANT_MEAN,RATE,3,M,DIM=2
 CRs=WHERE(RATE-REPLICATE(1,10)##M GT SIGMA/DELTA_TIME * 4)
                             CR_found=lonarr(2048.*2048.*10)
                             CR_found(CRs)=1
                             CR_found=REFORM(CR_found,2048.,2048.,10)
                             tvscl,total(CR_found,3)<;0.9

analysis,CR_found
end

pro analysis,CR_found
fits_read,'~/Desktop/CR_ramp.fits',CR_ramp
CR_frame=float(CR_ramp)
CR_frame[*,*,indgen(9)+1]=cr_ramp[*,*,indgen(9)+1]-cr_ramp[*,*,indgen(9)]
CR_binary = CR_frame*0
CR_binary(WHERE(CR_frame NE 0))=1
for ix=0,9 do begin & $
   FOUND = WHERE(CR_found[*,*,ix] EQ 1 AND CR_binary[*,*,ix] EQ 1, N_FOUND)   & $
   MISSED = WHERE(CR_found[*,*,ix] EQ 0 AND CR_binary[*,*,ix] EQ 1, N_MISSED) & $
   WRONG = WHERE(CR_found[*,*,ix] EQ 1 AND CR_binary[*,*,ix] EQ 0, N_WRONG) & $
   print,ix,N_FOUND,N_MISSEd,N_WRONG & $ 
   c0=CR_FRAME[*,*,ix] & $
    c0=reform(c0,2048.*2048.) & $
     iCR=WHERE(c0 NE 0) & $
   c1=CR_FOUND[*,*,ix] & $
    c1=reform(c1,2048.*2048.) & $
   plot,indgen(1000),histogram(c0[icr],binsize=1),xrange=[0,500],yrange=[0.1,4000],/ylog,psym=10 & $
   oplot,indgen(1000)-1,histogram(c1[icr]*c0[icr],binsize=1),color=255,psym=10  & $
   stop & $
endfor
end