On 03/31/2017 01:58 PM, mitali1203@gmail.com wrote:
> I am a newbie to IDL. I am currently working on event files (which
> are similar to FITS files) in IDL. I have written a code which will sort out
> event files from a list of directories and then compute detector plane
> histograms (DPH) (2D Histograms, basically) for each of those event
> files. In the process, I have also set up a 4-D array (4,25,64,64) and
> run a 'for' loop to save all DPH's to this array.
> However, I keep getting this error "Expression must be an array in
> this context:H", at the line where I have used 'hist_2d' function. I
> heard this is because there may not be any events in that particular
> file; hence IDL does not have any data values to create histogram. I
> guess I have to put a conditional loop (if..else) in order to detect
> such null event files. I am getting confused as to how this problem can
> be solved. It would be of great help if someone can troubleshoot this
> problem.
> Here is the code I have written- The line containing hist_2d
> function is where the error is-
>
> function compute_dph, infile, outpath
>
> dph=fltarr(4,25,64,64)
> for quad=0,3 do begin
> data=mrdfits(infile,quad+1, /unsigned,/dscale)
> for ebin=0,24 do begin
> index=where(data.energy gt 5+ebin*10 and data.energy le 5+(ebin+1)*10)
> qdph=hist_2d(data(index).detx,data(index).dety,MAX1=63,MAX2= 63,MIN1=0,MIN2=0)
> dph(quad,ebin,*,*)=qdph
> endfor
> endfor
>
> names=strsplit(infile,'/', /EXTRACT)
> stem=names(N_ELEMENTS(names)-1)
> outfile=strmid(stem,0,srtlen(stem)-3)+'dph'
> outfile=outpath+outfile
> save,dph,filename=outfile
>
> print,infile,outpath
> return,0
> end
Hi,
the following within the innermost loop should work:
index=where( (data.energy gt 5+ebin*10) and (data.energy le
5+(ebin+1)*10) ,cnt)
if cnt ne 0 then
dph[quad,ebin,*,*]=hist_2d(data(index).detx,data(index).dety ,MAX1=63,MAX2=63,MIN1=0,MIN2=0)
or assuming 0 <= detx,dety <=63 & detx,dety are integers you could put
instead of the inner loop
dat=64L*27L*data.dety+27L*data.detx+(0L>ceil((data.energy-5)/10) <26L)
dat=[reform(dat,n_elements(dat),/overwrite),64L*64L*27L-1L]
dph[0,*,*,*]=( reform(histogram(dat),[27,64,64]) )[1:25,*,*]
The latter I have not debugged.
I hope this helps, Markus
|