Hi,
It seems that the CODA data files opened by espace for reading are not
closed. ( at least I could not find the place where they are closed, if
anybody knows please let me know). This is not a problem most of the time
becase only a couple of data files are analyzed at a time.
But for optics optimization sometimes we have to analyze hundreds of data
files together. In such cases the fortran limit for the maximum
number of opened files is reached and espace crashes.
So I changes the routines espace_lib/process_data.f and
espace_halla/rawdata_til_eof.f to close data files after reading. Please
see the attached routines with the comments followed by NL:
If noone has any objections I would ask Ole to incorporate these changes
in to the next release.
Nilanga
PS: I also increased the parameter MAXRUNS to 150
c---------------------------------------------------------------------------
c Copyright (c) 1996 Southeastern Universities Research Association,
c Thomas Jefferson National Accelerator Facility
c
c This software was developed under a United States Government
c license described in the NOTICE file included as part of this
c distribution.
c
c Eddy A.J.M. Offermann, 12000 Jefferson Ave., Newport News, VA 23606
c Email: offerman@cebaf.gov Tel: (757) 249-7544 Fax: (757) 249-5800
c---------------------------------------------------------------------------
c
subroutine process_data(scanfile,scan_it,review_it,loop,
> first_interrupt,last_interrupt)
implicit none
include 'parameter.h'
include 'coolhands/cool_type.h'
include 'coolhands/coolhands.h'
include 'espace_type.h'
include 'transport.h'
include 'detector.h'
include 'espace_user.h'
include 'experiment.h'
include 'constant.h'
include 'spectra/communications.h'
character scanfile*(*)
logical scan_it,review_it,loop
integer last_interrupt,first_interrupt
logical filter_til_eof,rawdata_til_eof,geterr
real*8 reaction_z,dot3
integer nbytes,isp,jbit,jspec,indx,ifo,i
logical open_it,scaler_event,data_event,filter_indx
real*8 distance,distance_min,center,num,denom,lambda,
> react_z,dreact_z,weight_sum,p3v(3,2)
real*8 xdiff, zdiff, ratdiff, x0(2), z0(2), rat(2)
c
c April 22, 98, M. Kuss
c initialization of the different reaction points
c
do i = 1, 3
beam.react.v(i) = rbig
do isp = 1, 2
beam.spec(isp).v(i) = rbig
enddo
enddo
twoarm_x = rbig
twoarm_z = rbig
if (loop) then
open_it = .false.
else
open_it = .true.
endif
c
c NL: Thu Mar 25 1999
c changed nr_interrupts.lt.last_interrupt to
c nr_interrupts.le.last_interrupt below to make sure that the
c CODA file is closed when the desired number of events are analyzed.
c
nr_interrupts = 0
do while(last_interrupt.lt.0.or.
> nr_interrupts.le.last_interrupt)
if (rawdata_input) then
c
c NL: Thu Mar 25 1999
c pass the last_interrupt to rawdata_til_eof to close the
c CODA file when this number is reached.
c
if (rawdata_til_eof(open_it,scanfile,nbytes,
> scaler_event,first_interrupt,data_event,
> last_interrupt))
> return
else
if (filter_til_eof(open_it,scanfile,nbytes,
> scaler_event,first_interrupt,data_event))
> return
endif
if (geterr(1)) then
call derror('PROCESS_DATA: Error in the EVENT decoding')
return
endif
open_it = .false.
c If it is a run-number event, make sure that it is written to the
c filtered file even when the filter-condition is not fulfilled and
c the event_number is not in the required range.
if (event_type_espace.eq.17) then
do jspec = 1,nspec+nfilter+noptimize
filter_indx = ((jspec.gt.nspec).and.
> .not.(jspec.gt.nspec+nfilter))
if (filter_indx) then
indx = jspec-nspec
call filter_event(filter(indx).unit)
endif
enddo
endif
if (nr_interrupts.ge.first_interrupt) then
if (data_event) then
do jbit = 1,NHITBIT
hits(jbit) = .false.
enddo
c If the user has indicated that he wants to analyze the
c particular spectrometer, the appropriate subroutine is called.
c In this routine is checked first whether any data has been
c registered in this spectrometer for this interrupt.
weight_sum = 0.d0
beam.react.z = 0.d0
do isp = 1,2
beam.spec(isp).x = rbig
beam.spec(isp).y = rbig
beam.spec(isp).z = rbig
if (spstat.sp(isp).analyze) then
call process_spectrometer(spdetector.sp(isp),
> sptransport.sp(isp),
> spstat.sp(isp),hits)
c Calculate the z-coordinate of the reaction point using spectrometer
c ISP and the beam raster information if the variable RECONSTRUCT > 0
c If not, assume scattering took place in the center of the target.
if (beam.reconstruct.gt.0) then
if (tar.model.gt.1.or.
> (tar.model.eq.1.and.
> tar.thin.nfoil.gt.1)) then
react_z =reaction_z(tar,sptransport.sp(isp),
> beam,dreact_z)
c
c Jan 15, 98, M. Kuss
c beam.spec(isp).z stores react_z for one spectrometer
c
beam.spec(isp).z = react_z
beam.react.z = beam.react.z+
> react_z/dreact_z**2
weight_sum = weight_sum+1./dreact_z**2
endif
endif
endif
enddo
c Average z over all the spectrometers that are being analyzed
if (weight_sum.gt.0.d0) then
beam.react.z = beam.react.z/weight_sum
beam.react.x = beam.pos.x+beam.react.z*beam.angle.th
beam.react.y = beam.pos.y+beam.react.z*beam.angle.ph
c
c Jan 15, 98 M. Kuss
c beam.spec(isp).x and y store the reaction point for one spectrometer.
c
do isp = 1,2
if (spstat.sp(isp).analyze) then
beam.spec(isp).x = beam.pos.x
* + beam.spec(isp).z * beam.angle.th
beam.spec(isp).y = beam.pos.y
* + beam.spec(isp).z * beam.angle.ph
endif
enddo
else
beam.react.x = rbig
beam.react.y = rbig
beam.react.z = rbig
endif
c write ( 6, * ) 'process_data ', beam.react.v
c
c March 11, 98, M. Kuss, with code of St. Jaminion
c calculate the intersection of the two vertices
c
if( spstat.sp(1).analyze .AND. spstat.sp(2).analyze .AND.
* sptransport.sp(1).ok .AND. sptransport.sp(2).ok ) then
c write ( 6, * ) '*********************************'
do isp = 1, 2
call transport_to_beam(sptransport.sp(isp).index,
* sptransport.sp(isp).particle.pmom,
* sptransport.sp(isp).particle.tg.th,
* sptransport.sp(isp).particle.tg.ph,
* p3v(1,isp))
c write ( 6, * ) sptransport.sp(isp).particle.tg.y
c write ( 6, * ) (p3v(i,isp),i=1,3)
x0(isp) = sptransport.sp(isp).x_off
* + sptransport.sp(isp).particle.tg.y
* * sptransport.sp(isp).cs_th_geo
z0(isp) = sptransport.sp(isp).z_off
* - sptransport.sp(isp).particle.tg.y
* * sptransport.sp(isp).sn_th_geo
rat(isp) = p3v(1,isp) / p3v(3,isp)
enddo
c
xdiff = x0(2) - x0(1)
zdiff = z0(2) - z0(1)
ratdiff = rat(2) - rat(1)
twoarm_z = (rat(2)*z0(2)-rat(1)*z0(1)-xdiff) / ratdiff
twoarm_x = x0(1) + rat(1) * ( twoarm_z - z0(1) )
c write ( 6, * ) twoarm_x, twoarm_z
endif
c
c Use the reconstructed reaction point to determine in which foil the
c electron scattered (of course, trivial for one foil).
if (tar.model.eq.1) then
tar.thin.ifoil = (tar.thin.nfoil+1)/2
if (beam.reconstruct.gt.0) then
if (tar.thin.nfoil.gt.1) then
distance_min = rbig
do ifo = 1,tar.thin.nfoil
center = 0.01*dot3(tar.thin.foil(ifo).
> center_hall.v,tar.thin.norm.v)
distance = abs(dot3(beam.react.v,
> tar.thin.norm.v)-center)
if (distance.lt.distance_min) then
tar.thin.ifoil = ifo
distance_min = distance
endif
enddo
endif
c Force the reaction point to be in the selected foil by determining
c the intersection of beam and foil.
ifo = tar.thin.ifoil
num = beam.pos.x-tar.thin.foil(ifo).center_hall.x/
> 100.+beam.angle.th*
> (tar.thin.foil(ifo).center_hall.z/100.-
> beam.pos.z)
denom = tar.thin.surface.x-
> beam.angle.th*tar.thin.surface.z
lambda = num/denom
beam.react.z =tar.thin.foil(ifo).center_hall.z/100.
> +lambda*tar.thin.surface.z
beam.react.x = beam.pos.x+
> beam.react.z*beam.angle.th
beam.react.y = beam.pos.y+
> beam.react.z*beam.angle.ph
endif
endif
c Convert the reconstructed vertex of each analyzed spectrometer
c to the beam coordinate system. Also the energy loss and
c extended target effects on the momentum are calculated.
c This part has to be done AFTER processing the information of
c each spectrometer because the z-coordinate of the reaction point
c at the target needs information input from all spectrometers.
do isp = 1,2
if (spstat.sp(isp).analyze) then
call beam_vertex(mom_e0,momloss_e0,beam,tar,
> ilepton,ihadron,
> sptransport.sp(isp),
> spdetector.sp(isp))
endif
enddo
c Calculate the kinematic variables.
c There are two scenarios:
c 1) SPSTAT.COIN.ANALYZE=.TRUE., so coincidence analysis has been
c requested and therefore the electron kinematics and the
c coincidence related ones are calculated.
c 2) SPSTAT.COIN.ANALYZE=.FALSE.,SPSTAT.E.ANALYZE=.TRUE.
c only the electron kinematic variables are calculated.
call calc_kinematics(sptransport,spstat,ilepton,ihadron)
if (spstat.coin.analyze) then
call correct_cointime(sptransport,spdetector,
> spstat.coin)
endif
nwie = nbytes
if (review_it) call revsub
if (scan_it) call scansub
if (spstat.sp(1).analyze) call shcal_events
> (sptransport.sp(1),spdetector.sp(1))
else if (scaler_event) then
c If it is a scaler event, make sure that it is written to the
c filtered file even when the filter-condition is not fulfilled.
do jspec = 1,nspec+nfilter+noptimize
filter_indx = ((jspec.gt.nspec).and.
> .not.(jspec.gt.nspec+nfilter))
if (filter_indx) then
indx = jspec-nspec
call filter_event(filter(indx).unit)
endif
enddo
endif
endif
enddo
return
end
c---------------------------------------------------------------------------
c Copyright (c) 1996 Southeastern Universities Research Association,
c Thomas Jefferson National Accelerator Facility
c
c This software was developed under a United States Government
c license described in the NOTICE file included as part of this
c distribution.
c
c Eddy A.J.M. Offermann, 12000 Jefferson Ave., Newport News, VA 23606
c Email: offerman@cebaf.gov Tel: (757) 249-7544 Fax: (757) 249-5800
c---------------------------------------------------------------------------
c
logical function rawdata_til_eof(open,scanfile,nbytes,
> scaler_event,first_interrupt,
> data_event,last_interrupt)
c
c----------------------------------------------------------------------
c
c Regarding the trigger types. We have 12 inputs to the TS:
c First 4 can be prescaled by 2^24, next 4 by 2^16 (65536) and
c 4 cannot be prescaled.For the hall A data we have now the following
c event_type definition scheme:
c
c T1 = Golden electron-arm singles
c T2 = Junk electron-arm
c T3 = Golden hadron-arm singles
c T4 = Junk hadron-arm
c T5 = Coincidence electron-hadron
c T6,T7 = open for creative people
c T8 = pulser trigger
c T9,T10,T11,T12 = open (cannot prescale these)
implicit none
include 'spectra/communications.h'
include 'espace_lib/parameter.h'
include 'espace_lib/espace_type.h'
include 'espace_lib/detector.h'
include 'espace_lib/transport.h'
include 'espace_lib/experiment.h'
include 'espace_lib/espace_user.h'
include 'raw_type.h'
include 'rawdata.h'
include 'decoder.h'
include 'run_info.h'
include 'coda.h'
c arguments
character scanfile*(*)
integer nbytes,first_interrupt,last_interrupt
logical open,scaler_event,data_event
c functions
logical interpret_event_coin,interpret_event_e,interpret_event_h,
> corrupt_synch,superperiod_check,interpret_bpm_rast
integer evopen,evread
c local variables
integer i,iw,istatus,ihandle,maxlen,event_length,
& i1,i2,jdev,jflag
logical statusok,found
integer ntrk,nclu1,nclu2
integer MAXIARRAY,MAXCSTRING
parameter(MAXIARRAY=500,MAXCSTRING=4*MAXIARRAY)
structure /charint/
union
map
integer int(MAXIARRAY)
end map
map
character*(MAXCSTRING) char
end map
end union
end structure
record /charint/ cfigstring
integer ev_offset
parameter(ev_offset=30)
integer bout,bevcnt
data bout,bevcnt / 1, 0 /
save bout,bevcnt
save ihandle
c--------------------------------
rawdata_til_eof = .false.
scaler_event = .false.
data_event = .false.
statusok = .true.
nbytes = 0
maxlen = MAXLENGTH
c Opening of the raw data file is done by the CODA routine EVOPEN.
if (open) then
call nospac(scanfile,i1,i2)
istatus = evopen(scanfile(i1:i2),'r',ihandle)
if (istatus.ne.0) then
call derror('RAWDATA_TIL_EOF: error opening raw data file')
rawdata_til_eof = .true.
return
endif
endif
c Loop around till uncorrupted data has been found.
c Keep track of the number of synch problems found
1 continue
istatus = evread(ihandle,evbuffer,maxlen)
if (istatus.ne.0) then
rawdata_til_eof = .true.
c
c NL: Thu Mar 25 1999
c
c No more data in the file, so close the CODA data file
c
call evclose(ihandle)
return
endif
call nospac(scanfile,i1,i2)
if (corrupt_synch(scanfile(i1:i2),debug_sync,evbuffer,
> exp_scaler.data.syn_extr.run.value,
> exp_scaler.data.syn_miss.run.value)) then
c 3 May 98 D.Prout
c Do not skip 'missing gate' events, but analyze them using an offset
c event type.c
c goto 1
event_type_espace=event_type+ev_offset
endif
c
c Get event length and type
event_length = evbuffer(1)
event_type = ishft(evbuffer(2),-16)
event_type_espace = event_type
nbytes = 4*event_length
c Now go down the list of all the different EVENT_TYPES. Only
c for physics events we have to continue the analysis !!
if (event_type.eq.17) then
c Get run number and run type
run_number = evbuffer(4) ! get run number from PRESTART event
run_number_espace = run_number
run_type = evbuffer(5) ! get run type
else if (event_type.eq.133) then
c A pre-scaler event.
do i = 1,min(MAXCSTRING,evbuffer(1)+1)
cfigstring.int(i) = evbuffer(i)
end do
call config_parse(cfigstring)
else if (event_type.eq.spstat.scaler.flag(1).or.
& event_type.eq.spstat.scaler.flag(2)) then
c A scaler event; nr_interrupts of the previous routine entry is
c used here !!
scaler_event = .true.
if (nr_interrupts.ge.first_interrupt) then
call scaler_update(evbuffer,event_length,istatus)
if (istatus.ne.0)
> call derror('RAWDATA_TIL_EOF: error in scaler_update')
endif
else if (event_type.eq.131) then
call get_bcm
endif
c Continue if event type is < 16 ==> physics event, see scheme above.
if (event_type.lt.16) then
c Get event number and increment the event counter
event_num = evbuffer(5) ! store event number
event_ctr = event_ctr+1 ! increment event counter
nr_interrupts = event_num
if ( mod(event_num,1000) .EQ. 0 .OR. debug ) then
write ( 6, * ) ' event number: ', event_num
endif
c Return if event length is < 0
if (event_length.le.8) return ! no data
c Return if interrupt number < first_interrupt
if (nr_interrupts.lt.first_interrupt) return
if (debug_coda) then
write(6,'(2x,a,i4)')'===> event number = ',evbuffer(5)
do i = 1, event_length+1
write(6,'(3x,i4,a,z8)') i,' : ',evbuffer(i)
enddo
endif
c Now we get serious, it was a data event and ESPACE has to analyze
c it.
data_event = .true.
c Check whether an electron, hadron and/or coincidence event has been
c registered.
do jdev = 1,3
jflag = 1
found = .false.
spstat.dev(jdev).registered = .false.
do while (.not.found.and.jflag.le.spstat.dev(jdev).nr_flag)
if (event_type.eq.spstat.dev(jdev).flag(jflag)) then
spstat.dev(jdev).registered = .true.
found = .true.
else
jflag = jflag+1
endif
enddo
enddo
call init_eventcounters
call event_decoder(evbuffer,event_length,istatus)
if (istatus.gt.1) then
call derror('RAWDATA_TIL_EOF:error in EVENT_DECODER')
return
endif
c
c Interpret the BPM and Raster Data
statusok = statusok.and.interpret_bpm_rast()
c
c
c Unpack digitized data from the CODA banks - only if we want to
c analyze it later.
if (spstat.coin.analyze) then
statusok = statusok.and.
> interpret_event_coin(raw,wdctr,sptransport,
> spdetector,
> tdc_trig_cor,
> vdc_tdc_cut)
else
if (spstat.e.analyze)
> statusok = statusok.and.
> interpret_event_e(raw,wdctr,sptransport.e,
> spdetector.e,
> tdc_trig_cor,
> vdc_tdc_cut)
if (spstat.h.analyze)
> statusok = statusok.and.
> interpret_event_h(raw,wdctr,sptransport.h,
> spdetector.h,
> tdc_trig_cor,
> vdc_tdc_cut)
endif
endif
c
c Februar 18, 98, M. Kuss
c read the latched trigger pattern
c
call read_spare
c
c NL: If the required number of events has been reached
c 03/25/99 close the coda file
c
if (nr_interrupts.ge.last_interrupt) then
call evclose(ihandle)
endif
return
end