![]() ![]() |
![]() |
File: [HallC] / pol_hms_single / hms / mc_hms_recon.f.simc
(download)
Revision: 1.1.1.1 (vendor branch), Tue Dec 16 13:39:55 2003 UTC (20 years, 9 months ago) by jones Branch: poltar, MAIN CVS Tags: start, HEAD Changes since 1.1: +0 -0 lines Starting Poltar HMS single MC |
subroutine mc_hms_recon (delta_p,delta_t,delta_phi,y_tgt,fry) C+______________________________________________________________________________ ! ! MC_HMS_RECON : Reconstruct target quantities from tracks. ! This subroutine is part of the MC_HMS program. ! ! Right-handed coordinates are assumed: X=down, Z=downstream, Y = (Z cross X) ! ! Inputs are from common block in track_*.inc (except for fry): ! xs, ys, fry are in cm. ! dxdzs, dydzs are unitless slopes (we say "radians" to contrast to "mr"). ! Matrix Elements want meters and radians, so have to convert cm-->m. ! Output: ! delta_p is deviation from central momentum(unitless). Convert to % for outpu+ ! delta_t, delta_phi are in "radians". ! y_tgt is in meters, convert to cm for output. ! ! Author: D. Potterveld, ANL, 18-Mar-1993 ! ! Modification History: ! ! 2-August-1995 (RMM) Hardcoded in an extra index on the coeff, expon, and ! n_terms variables to specify HMS. (HMS = 1) ! ! 19-AUG-1993 (DHP) Modified to use COSY INFINITY reconstruction coefficients. C-______________________________________________________________________________ implicit none include 'spectrometers.inc' integer*4 specnum parameter (specnum = 1) !this is the HMS routine C Argument definitions. real*8 delta_p,delta_t,delta_phi,y_tgt real*8 fry !vertical position at target (+y=down) C Cosy reconstruction matrix elements. integer*4 max_elements parameter (max_elements = 1000) real*8 coeff(nspectr,4,max_elements) integer*2 expon(nspectr,5,max_elements) integer*4 n_terms(nspectr),max_order real*8 sum(4),hut(5),term C Misc. variables. integer*4 i,j integer*4 chan logical firsttime /.true./ character*132 line C Functions. logical locforunt C No amnesia, please... save C ============================= Executable Code ================================ C First time through, read in coefficients from data file. if (firsttime) then if (.not.locforunt(chan)) stop 'MC_HMS_RECON: No I/O channels!' open (unit=chan,status='old',readonly,file='hms/recon_cosy.dat') ! Skip past header. line = '!' do while (line(1:1).eq.'!') read (chan,1001) line enddo ! Read in coefficients and exponents. n_terms(specnum) = 0 max_order = 0 do while (line(1:4).ne.' ---') n_terms(specnum) = n_terms(specnum) + 1 if (n_terms(specnum).gt.max_elements) > stop 'WCRECON: too many COSY terms!' read (line,1200) (coeff(specnum,i,n_terms(specnum)),i=1,4), > (expon(specnum,j,n_terms(specnum)),j=1,5) read (chan,1001) line max_order = max(max_order, expon(specnum,1,n_terms(specnum)) + > expon(specnum,2,n_terms(specnum)) + > expon(specnum,3,n_terms(specnum)) + > expon(specnum,4,n_terms(specnum)) + > expon(specnum,5,n_terms(specnum))) enddo c write(6,*) 'HMS: N_TERMS, MAX_ORDER = ',n_terms(specnum),max_order close (unit=chan) firsttime = .false. endif C Reset COSY sums. do i = 1,4 sum(i) = 0. enddo C Convert hut quantities to right-handed coordinates, in meters and "radians". C Make sure hut(5) is non-zero, to avoid taking 0.0**0 (which crashes) hut(1) = xs/100. !cm --> m hut(2) = dxdzs !slope ("radians") hut(3) = ys/100. !cm --> m hut(4) = dydzs !slope ("radians") hut(5) = fry/100. !vert. position at target(cm-->m) if (abs(hut(5)).le.1.d-30) hut(5)=1.d-30 C Compute COSY sums. do i = 1,n_terms(specnum) term = hut(1)**expon(specnum,1,i) * hut(2)**expon(specnum,2,i) > * hut(3)**expon(specnum,3,i) * hut(4)**expon(specnum,4,i) > * hut(5)**expon(specnum,5,i) sum(1) = sum(1) + term*coeff(specnum,1,i) sum(2) = sum(2) + term*coeff(specnum,2,i) sum(3) = sum(3) + term*coeff(specnum,3,i) sum(4) = sum(4) + term*coeff(specnum,4,i) enddo C Load output values. delta_phi = sum(1) !slope ("radians") y_tgt = sum(2)*100. !m --> cm delta_t = sum(3) !slope ("radians") delta_p = sum(4)*100. !percent deviation return C ============================ Format Statements =============================== 1001 format(a) 1200 format(1x,4g16.9,1x,5i1) END
Analyzer/Replay: Mark Jones, Documents: Stephen Wood |
Powered by ViewCVS 0.9.2-cvsgraph-1.4.0 |