1 gaskelld 1.1 subroutine mc_hrsl_recon (delta_p,delta_t,delta_phi,y_tgt,fry)
2
3 C+______________________________________________________________________________
4 !
5 ! MC_HRSL_RECON : Reconstruct target quantities from tracks.
6 ! This subroutine is part of the MC_HRSL program.
7 !
8 ! Right-handed coordinates are assumed: X=down, Z=downstream, Y = (Z cross X)
9 !
10 ! Inputs are from common block in track_*.inc (except for fry):
11 ! xs, ys, fry are in cm.
12 ! dxdzs, dydzs are unitless slopes (we say "radians" to contrast to "mr").
13 ! Matrix Elements want meters and radians, so have to convert cm-->m.
14 ! Output:
15 ! delta_p is deviation from central momentum(unitless). Convert to % for output.
16 ! delta_t, delta_phi are in "radians".
17 ! y_tgt is in meters, convert to cm for output.
18 !
19 ! Author: D. Potterveld, ANL, 18-Mar-1993
20 !
21 ! Modification History:
22 gaskelld 1.1 !
23 ! 2-August-1995 (RMM) Hardcoded in an extra index on the coeff, expon, and
24 ! n_terms variables to specify HRSL. (HRSL = 4)
25 !
26 ! 19-AUG-1993 (DHP) Modified to use COSY INFINITY reconstruction coefficients.
27 C-______________________________________________________________________________
28
29 implicit none
30
31 include '../track.inc'
32
33 integer*4 specnum
34 parameter (specnum = 4) !this is the HRSL routine
35
36 C Argument definitions.
37
38 real*8 delta_p,delta_t,delta_phi,y_tgt
39 real*8 fry !vertical position at target (+y=down)
40
41 C Cosy reconstruction matrix elements.
42
43 gaskelld 1.1 integer*4 max_elements
44 parameter (max_elements = 1000)
45
46 integer*4 nspectr
47 parameter (nspectr=4)
48
49 real*8 coeff(nspectr,4,max_elements)
50 integer*2 expon(nspectr,4,max_elements)
51 C integer*2 expon(nspectr,5,max_elements)
52 integer*4 n_terms(nspectr),max_order
53 C real*8 sum(4),hut(5),term
54 real*8 sum(4),hut(4),term
55
56 C Misc. variables.
57
58 integer*4 i,j
59 integer*4 chan
60
61 logical firsttime /.true./
62 character*132 line
63
64 gaskelld 1.1 C Functions.
65
66 logical locforunt
67
68 C No amnesia, please...
69
70 save
71
72 C ============================= Executable Code ================================
73
74 C First time through, read in coefficients from data file.
75
76 if (firsttime) then
77 if (.not.locforunt(chan)) stop 'MC_HRSL_RECON: No I/O channels!'
78 open (unit=chan,status='old',readonly,file='hrsl/hrs_recon_cosy.dat')
79
80 ! Skip past header.
81
82 line = '!'
83 do while (line(1:1).eq.'!')
84 read (chan,1001) line
85 gaskelld 1.1 enddo
86
87 ! Read in coefficients and exponents.
88 n_terms(specnum) = 0
89 max_order = 0
90 do while (line(1:4).ne.' ---')
91 n_terms(specnum) = n_terms(specnum) + 1
92 if (n_terms(specnum).gt.max_elements)
93 > stop 'WCRECON: too many COSY terms!'
94 read (line,1200) (coeff(specnum,i,n_terms(specnum)),i=1,4),
95 > (expon(specnum,j,n_terms(specnum)),j=1,4)
96 read (chan,1001) line
97 max_order = max(max_order, expon(specnum,1,n_terms(specnum)) +
98 > expon(specnum,2,n_terms(specnum)) +
99 > expon(specnum,3,n_terms(specnum)) +
100 > expon(specnum,4,n_terms(specnum)))
101 C > expon(specnum,5,n_terms(specnum)))
102 enddo
103 write(6,*) 'HRSL:N_TERMS, MAX_ORDER = ',n_terms(specnum),max_order
104 close (unit=chan)
105 firsttime = .false.
106 gaskelld 1.1 endif
107
108 C Reset COSY sums.
109
110 do i = 1,4
111 sum(i) = 0.
112 enddo
113
114 C Convert hut quantities to right-handed coordinates, in meters and "radians".
115 hut(1) = xs/100. !cm --> m
116 hut(2) = dxdzs !slope ("radians")
117 hut(3) = ys/100. !cm --> m
118 hut(4) = dydzs !slope ("radians")
119
120 C Compute COSY sums.
121
122 do i = 1,n_terms(specnum)
123 term = hut(1)**expon(specnum,1,i)*hut(2)**expon(specnum,2,i)
124 > * hut(3)**expon(specnum,3,i)*hut(4)**expon(specnum,4,i)
125 sum(1) = sum(1) + term*coeff(specnum,1,i)
126 sum(2) = sum(2) + term*coeff(specnum,2,i)
127 gaskelld 1.1 sum(3) = sum(3) + term*coeff(specnum,3,i)
128 sum(4) = sum(4) + term*coeff(specnum,4,i)
129 enddo
130
131 C Load output values.
132
133 delta_phi = sum(1) !slope ("radians")
134 y_tgt = sum(2)*100. !m --> cm
135 delta_t = sum(3) !slope ("radians")
136 delta_p = sum(4)*100. !percent deviation
137
138 return
139
140 C ============================ Format Statements ===============================
141
142 1001 format(a)
143 1200 format(1x,4g16.9,1x,4i1)
144
145 END
|