1 gaskelld 1.1 {------------------------------------------------------------------------------
2 HRS-MONTE.FOX
3
4 Generate Monte-Carlo forward and reconstruction matrix elements.
5
6 Hall A HRS spectrometer.
7
8 All output matrix elements are in TRANSPORT coordinates.
9 Note that COSY-7 uses different units than COSY-6 when printing the maps
10 in transport notation (procedure PT.) The reconstruction map is still
11 calculated using the old COSY-6 units. The output format corresponds to what
12 is expected by my program monte_hms.
13
14 X-dependent (raster) reconstruction map generated.
15
16 Version: 0.0-dev
17 Author: D. Potterveld, Argonne National Lab, 26-Jun-2001
18
19 Modification History:
20
21 ------------------------------------------------------------------------------
22 gaskelld 1.1 }
23 INCLUDE 'COSY' ;
24 PROCEDURE RUN ;
25
26 {=========================== Variable Declarations ===========================}
27
28 VARIABLE Z1C 1 ; VARIABLE L1 1 ; VARIABLE L1E 1 ;
29 VARIABLE Z2C 1 ; VARIABLE L2 1 ; VARIABLE L2E 1 ;
30 VARIABLE Z3C 1 ; VARIABLE L3 1 ; VARIABLE L3E 1 ;
31 VARIABLE ZDI 1 ;
32 VARIABLE DR1 1 ;
33 VARIABLE DR2 1 ;
34 VARIABLE DR3 1 ;
35 VARIABLE DRD 1 ;
36 VARIABLE Q1 1 ; { quadrupole strength }
37 VARIABLE Q2 1 ; { quadrupole strength }
38 VARIABLE Q3 1 ; VARIABLE DD3 1 ; {Q3 parameters }
39
40 VARIABLE Q1MID 1 ; VARIABLE Q2MID 1 ; VARIABLE Q3MID 1 ;
41
42 {dipole variables}
43 gaskelld 1.1
44 variable D_RADIUS 1 ; { dipole bend radius }
45 variable D_BEND 1 ; { dipole bend angle }
46 variable D_GAP 1 ; { dipole half gap }
47 variable D_ALPHA 1 ; { entrance shim angle }
48 variable D_BETA 1 ; { exit shim angle }
49 variable D_INDEX 1 ; { Field index }
50 variable D_N 1 5 ; { field radial dependence coefficients }
51 variable DELTA_D_N 1 ; { error on the field index }
52 variable D_S1 1 5 ; { entrance polynomial coefficients }
53 variable D_S2 1 5 ; { exit polynomial coefficients }
54 variable D_H1 1 ; { dipole entrance curvature }
55 variable D_H2 1 ; { dipole exit curvature }
56 variable D_S01 1 5 ;
57 variable D_S2 1 5 ;
58 variable D_S02 1 5 ;
59
60 VARIABLE OBJ 1 ; { function used in evaluating first order focus }
61 VARIABLE FILE 80 ; { filename }
62 VARIABLE LINE 80 ; { dummy for user to answer into }
63 VARIABLE TILT 1 ; { Focal plane tilt angle }
64 gaskelld 1.1 VARIABLE I 1 ;
65 VARIABLE J 1 ;
66 VARIABLE TMP 1 ;
67 VARIABLE SNAME 3 ; {Spectrometer name}
68 VARIABLE FRMODE 1 ; { Fringe field mode }
69 VARIABLE MORDER 1 ; { Maximum order of maps }
70 VARIABLE FORDER 1 ; { Order of forward maps }
71 VARIABLE RORDER 1 ; { Order of reconstruction map }
72 VARIABLE RES 1 4 ; { Resolution array }
73
74 VARIABLE TMPSPOS 1000 ; { Saved value of SPOS }
75 VARIABLE POSITION 1 ; { Position of each output plane }
76 VARIABLE GMAP 2000 8 ; { Global transfer map }
77 VARIABLE RMAP 2000 5 ; { Reconstruction map }
78
79 {========================== Procedure declarations ===========================}
80
81 { Print rmap to IU }
82
83 PROCEDURE PMR IU ;
84 VARIABLE I 1 ; VARIABLE M 2000 5 ;
85 gaskelld 1.1 LOOP I 1 5 ; M(I) := RMAP(I) + 0*DD(1) ; ENDLOOP ;
86 DAPRV M 4 5 TWOND IU ; ENDPROCEDURE ;
87
88 { Reconstruction with X as independent variable }
89
90 PROCEDURE RR1 M X0 X A Y B D PR AR N OR MR RES ;{TRAJECTORY RECONSTRUCTION}
91 PROCEDURE REC ; VARIABLE NN 1 ; VARIABLE S MAX(NM1,N+1) ;
92 VARIABLE I 1 ; VARIABLE M4 NM1 5 ; VARIABLE CM 1 6 ;
93 VARIABLE RI N 6 ; VARIABLE RO N 4 ; VARIABLE RB N 4 ; VARIABLE RE N 4 ;
94 VARIABLE IRAN 1 ; VARIABLE J NM1 ; VARIABLE NOM 1 ; VARIABLE IER 1 ;
95 VARIABLE MR4 NM1 4 ;
96 FUNCTION RAN IRAN ; RAN := 98765*IRAN+.12345 ; RAN := 2*(RAN-NINT(RAN)) ;
97 IRAN := RAN ; ENDFUNCTION ; IRAN := 0 ;
98 IF OR>0 ; LOOP I 1 4 ; CM(I+1) := DD(I) ;
99 ENDLOOP ; NOM := NOC ; CO OR ; CM(1) := DD(5) ; CM(6) := 0*DD(1) ;
100 IF ND=3 ; CM(6) := DD(4) ; CM(5) := 0*DD(1) ; ENDIF ;
101 POLVAL 1 M 4 CM 6 M4 4 ;
102 M4(5) := DD(5) ; MI M4 MR 5 IER ; CO NOM ; ENDIF ;
103 S := 0 ; LOOP I 1 N-1 ; S := S&0 ; ENDLOOP ;
104 LOOP I 1 N ; VELSET S I X0+X*RAN(IRAN) ; ENDLOOP ; RI(1) := S ;
105 LOOP I 1 N ; VELSET S I A*RAN(IRAN) ; ENDLOOP ; RI(2) := S ;
106 gaskelld 1.1 LOOP I 1 N ; VELSET S I Y*RAN(IRAN) ; ENDLOOP ; RI(3) := S ;
107 LOOP I 1 N ; VELSET S I B*RAN(IRAN) ; ENDLOOP ; RI(4) := S ;
108 LOOP I 1 N ; VELSET S I D*RAN(IRAN) ; ENDLOOP ; RI(5) := S ;
109 RI(6) := S ; POLVAL 1 M 4 RI 6 RO 4 ;
110 LOOP I 1 N ; VELSET S I PR*RAN(IRAN) ; ENDLOOP ; RE(1) := S ;
111 LOOP I 1 N ; VELSET S I AR*RAN(IRAN) ; ENDLOOP ; RE(2) := S ;
112 LOOP I 1 N ; VELSET S I PR*RAN(IRAN) ; ENDLOOP ; RE(3) := S ;
113 LOOP I 1 N ; VELSET S I AR*RAN(IRAN) ; ENDLOOP ; RE(4) := S ;
114 RO(1) := RO(1) + RE(1) ; RO(2) := RO(2) + RE(2) ;
115 RO(3) := RO(3) + RE(3) ; RO(4) := RO(4) + RE(4) ;
116 LOOP I 1 4 ; DAPLU MR(I) 5 X0 S ; MR4(I) := S ; ENDLOOP ;
117 POLVAL 1 MR4 4 RO 4 RB 4 ; RE(1) := RI(2)-RB(1) ; RE(2) := RI(3)-RB(2) ;
118 RE(3) := RI(4)-RB(3) ; RE(4) := RI(6)-RB(4) ;
119 LOOP J 1 4 ; RES(J) := 0 ; NN := 0 ; LOOP I 1 N ; VELGET RE(J) I S ;
120 IF S#0 ; RES(J) := RES(J) + ABS(S) ; NN := NN + 1 ; ENDIF ; ENDLOOP ;
121 RES(J) := NN/(RES(J)+1E-20)/4 ; IF IER#0 ; RES(J) := 0 ; ENDIF ;
122 ENDLOOP ; ENDPROCEDURE ;
123 NM1 := MAX(NM1,N) ; SCRLEN 2*NM1 ; REC ;
124 NM1 := NMON(NO+1,NV) ; SCRLEN LSCR ; ENDPROCEDURE ;
125
126 { Ask user for input value }
127 gaskelld 1.1
128 PROCEDURE ASK VAR PROMPT LL UL ;
129 VARIABLE TMP 1 ;
130 VAR := 2*ABS(UL) + 1 ;
131 WHILE (VAR<LL)+(VAR>UL) ; WRITE 6 '$'&PROMPT ; READ 5 VAR ;
132 ENDWHILE ; ENDPROCEDURE ;
133
134 { Output forward transformation data at an aperture and update internal maps }
135
136 PROCEDURE APERTURE LABEL REGION ;
137 WRITE 20 'NAME: '&LABEL ;
138 WRITE 20 'REGION: '®ION ;
139 WRITE 20 'OFFSET: '&SF(TMPSPOS,'(F25.16)')&' (in meters)' ;
140 WRITE 20 'LENGTH: '&SF(SPOS,'(F25.16)')&' (canonical length in meters)' ;
141 PT 20 ; {Output data}
142 TMPSPOS := TMPSPOS + SPOS ;
143 ANM MAP GMAP GMAP ; UM ; {Update global map + reset}
144 WRITE 6 ' Aperture at S = '&S(TMPSPOS) ;
145 ENDPROCEDURE ;
146
147 { Compute and output maps to inner aperture and exit of HRS quad magnet. }
148 gaskelld 1.1
149 PROCEDURE HRSQUAD NAME LEFF Q HEX OCT DEC DDEC A FRAC ;
150 VARIABLE DRIFT 1 ; VARIABLE LMID 1 ; VARIABLE QFLAG 1 ;
151
152 { Compute flag for element type }
153
154 QFLAG := ABS(HEX)+ABS(OCT)+ABS(DEC)+ABS(DDEC);
155 IF QFLAG=0 ; WRITE 6 NAME&': Using normal quad' ;
156 ELSEIF TRUE ; WRITE 6 NAME&': Using M5 quad' ;
157 ENDIF ;
158
159 { Compute the location of the inner aperture }
160
161 LMID := (LEFF*FRAC) ;
162
163 IF FRMODE#0 ; FR -1 ; {Entrance fringe field}
164 IF QFLAG=0 ; MQ LEFF Q A ;
165 ELSEIF TRUE ; M5 LEFF Q HEX OCT DEC DDEC A ; ENDIF ;
166 ENDIF;
167
168 FR 0 ; {First part of quad}
169 gaskelld 1.1 IF QFLAG=0 ; MQ LMID Q A ;
170 ELSEIF TRUE ; M5 LMID Q HEX OCT DEC DDEC A ; ENDIF ;
171
172 APERTURE NAME&'_MID' 'FRONT' ;
173
174 IF QFLAG=0 ; MQ (LEFF-LMID) Q A ; {Second part of quad}
175 ELSEIF TRUE ; M5 LEFF-LMID Q HEX OCT DEC DDEC A ; ENDIF ;
176
177 IF FRMODE#0 ; FR -2 ; {Exit fringe field}
178 IF QFLAG=0 ; MQ LEFF Q A ;
179 ELSEIF TRUE ; M5 LEFF Q HEX OCT DEC DDEC A ; ENDIF ;
180 ENDIF;
181
182 APERTURE NAME&'_EXIT' 'BACK' ;
183
184 FR FRMODE ; ENDPROCEDURE ; {Restore fringe field mode}
185
186
187 {================================= Main Code =================================}
188
189 { Computation order specifications }
190 gaskelld 1.1
191 MORDER := 5 ; {Max order for calculations}
192 ASK FORDER 'Order for maps (1-'&SF(morder,'(i1)')&'): ' 1 morder ;
193
194 OV FORDER 3 0 ; {Order Value <order><phase dim><#par>}
195 CO FORDER ; RORDER := FORDER ;
196
197 { Beam specifications }
198
199 RPE 837.56 ; {REF PARTICLE <Momentum>}
200
201 SB .01 .030 0 .05 .020 0 0 .05 0 0 0 ;
202 { --> PX PA R12 PY PB R34 PT PD R56 PG PZ }
203
204
205 { Spectrometer name }
206
207 SNAME := 'HRS' ;
208
209 { Fringe-field specifications.}
210
211 gaskelld 1.1 ASK frmode 'Fringe-field mode (0-3): ' 0 3 ;
212 FR FRMODE ;
213
214 { From Dave Meekins HRS file }
215
216 FC 1 1 1 0.04725 2.2395 -0.9768 0.7288 -0.1299 0.0222 ;
217 FC 1 2 1 0.04725 2.2395 -0.9768 0.7288 -0.1299 0.0222 ;
218
219 LOOP I 2 6 ; LOOP J 1 2 ;
220 FC I J 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
221 ENDLOOP ; ENDLOOP ;
222
223 { Magnet geometry}
224
225 L1 := 0.9413 ; {Magnetic length of Q1}
226 L2 := 1.8266 ; {Magnetic length of Q2}
227 L3 := 1.8268 ; {Magnetic length of Q3}
228
229 DR1 := 1.5903 ; {Drift from target to Q1}
230 DR2 := 1.172 ; {Drift from Q1 to Q2}
231 DRD := 4.4308 ; {Drift from Q2 to dipole}
232 gaskelld 1.1 DR3 := 1.5925 ; {DRIFT FROM DIPOLE TO Q3}
233
234 { QUAD STRENGTHS }
235
236 Q1 := -0.243838 ;
237 Q2 := 0.1934 ;
238 Q3 := 0.17892 ;
239 DD3 := 0.00245 ;
240
241 { Dipole parameters }
242
243 D_RADIUS := 8.40 ; D_BEND := 45.0 ; D_GAP := 0.125 ;
244 D_ALPHA := -30.0 ; D_BETA := -30.0 ; D_INDEX := -1.25 ;
245
246 D_N(1) := D_INDEX ;
247 loop i 2 5 ; d_n(i) := -d_n(i-1)*D_INDEX ; endloop ;
248
249 D_S1(1) := TAN(-30.*DEGRAD) ; {Entrance pole shape}
250 D_S1(2) := 0.4014526 ;
251 D_S1(3) := 0.21710 ;
252 D_S1(4) := -0.0825 ;
253 gaskelld 1.1 D_S1(5) := -0.3110 ;
254 D_S2(1) := TAN(-30.*DEGRAD) ; {Exit pole shape}
255 D_S2(2) := -0.4433318 ;
256 D_S2(3) := 0.18020 ;
257 D_S2(4) := -0.2443 ;
258 D_S2(5) := 1.1649 ;
259
260
261 { Begin composing forward maps }
262
263 file := SNAME&'_fr'&sf(frmode,'(i1)')&'_for_maps.dat' ;
264 openf 20 file 'NEW' ;
265
266 FR FRMODE ;
267 UM ; SM GMAP ; TMPSPOS := SPOS ; {Initialize global map}
268
269 { FRONT QUADS }
270
271 DL DR1 ; {Drift to Q1 physical aperture}
272 APERTURE 'Q1_ENTRANCE' 'FULL' ;
273 HRSQUAD 'Q1' L1 Q1 0 0 0 0 0.15 2/3 ; {Q1}
274 gaskelld 1.1
275 DL DR2 ; {Drift to Q2}
276 APERTURE 'Q2_ENTRANCE' 'FULL' ;
277 HRSQUAD 'Q2' L2 Q2 0 0 0 0 0.30 2/3 ; {Q2}
278
279 { DIPOLE }
280
281 DL DRD ; {Drift to dipole}
282 APERTURE 'DIPOLE_ENTRANCE' 'FULL' ;
283
284 { Consistency check: we better be at the right Z! }
285
286 {
287 IF ABS(TMPSPOS-ZDI)>1E-9 ;
288 WRITE 6 'Warning: Dipole location is bad:' ;
289 WRITE 6 'POSITION = '&S(TMPSPOS)&', but should be '&S(ZDI) ;
290 ENDIF ;
291 }
292 MC D_RADIUS D_BEND D_GAP D_N D_S1 D_S2 5 ; {B1}
293
294 APERTURE 'DIPOLE_EXIT' 'FULL' ;
295 gaskelld 1.1
296 { REAR QUAD }
297
298 DL DR3 ; {Drift to Q1 physical aperture}
299 APERTURE 'Q3_ENTRANCE' 'FULL' ;
300 HRSQUAD 'Q3' L3 Q3 0 0 0 DD3 0.30 2/3 ; {Q3}
301
302 { Reconstruction plane }
303
304 DL 3.4523 ; {Drift to detector plane}
305 APERTURE 'RECON PLANE' 'FULL' ;
306
307 closef 20 ; {Done outputing forward MAPS}
308 { Convert GMAP from canonical to TRANSPORT form}
309
310 CATR GMAP MAP TMPSPOS ; spos := tmpspos ;
311
312 { Compute focal objective function }
313
314 OBJ := SQRT(ME(1,2)*ME(1,2)+ME(3,4)*ME(3,4)) ;
315 WRITE 6 ' Q1,Q2, 1st order OBJ = ' Q1&Q2&OBJ ;
316 gaskelld 1.1 WRITE 6 ' Total length = '&S(CONS(SPOS))&' meters' ;
317
318 { Compute focal plane tilt angle }
319
320 TILT := (180*atan(ME(1,1)*ME(1,26)/ME(1,6))/pi)+90 ;
321 write 6 ' Focal plane tilt ='&SF(tilt,'(f8.3)')&' degrees' ;
322
323 { Compute RECON map }
324
325 RR1 MAP 0 0 .04 0 .04 .15 0 0 1000 RORDER RMAP RES ;
326 write 6 ' RECON resolutions = ' RES(1)&RES(2)&RES(3)&RES(4) ;
327 file := SNAME&'_fr'&sf(frmode,'(i1)')&'_rec.dat' ;
328 openf 20 file 'NEW' ; pmr 20 ; closef 20 ;
329
330 ENDPROCEDURE ;
331 RUN ;
332 END ;
333
334
335
336
337 gaskelld 1.1
338
339
340
|
Return to
hrs-monte.fox
CVS log
Up to [HallC] / simc_gfortran / hrsr