1 |
pam-fi |
1.1 |
************************************************************************* |
2 |
|
|
* |
3 |
|
|
* Subroutine read_B_inner.f |
4 |
|
|
* |
5 |
|
|
* it reads from rz files the two magnetic field maps taken inside the |
6 |
|
|
* spectrometer cavity and fills the variables in common_B_inner.f |
7 |
|
|
* |
8 |
|
|
* needs: |
9 |
|
|
* - ../common/common_B_inner.f common file for the inner magnetic field map |
10 |
|
|
* - .rz map files in ./ containing coordinates of measured points, Bx, By |
11 |
|
|
* and Bz components + errors |
12 |
|
|
* |
13 |
|
|
* output variables: (see common_B_inner.f) |
14 |
|
|
* - px#(nx,3) with #=1,2 for the 2 maps |
15 |
|
|
* - py#(ny,3) |
16 |
|
|
* - pz#(nz,3) |
17 |
|
|
* - b#(nx,ny,nz,3) |
18 |
|
|
* |
19 |
|
|
************************************************************************* |
20 |
|
|
|
21 |
|
|
subroutine read_B_inner |
22 |
|
|
|
23 |
|
|
implicit double precision (a-h,o-z) |
24 |
|
|
include '../common/common_B_inner.f' |
25 |
|
|
|
26 |
|
|
|
27 |
|
|
c------------------------------------------------------------------------ |
28 |
|
|
c |
29 |
|
|
c local variables |
30 |
|
|
c |
31 |
|
|
c------------------------------------------------------------------------ |
32 |
|
|
|
33 |
|
|
character*64 Bmap_file !magnetic field file name |
34 |
|
|
c character*120 cmd1 |
35 |
|
|
c character*120 cmd2 |
36 |
|
|
parameter (lun_Bmap_file=66) !magnetic field map file id number |
37 |
|
|
|
38 |
|
|
parameter (ntpl_Bmap=20) !ntuple identifier |
39 |
|
|
|
40 |
|
|
REAL PFX(3),FX,DFX, !Bx field component coordinates in m, value and error in T |
41 |
|
|
$ PFY(3),FY,DFY |
42 |
|
|
$ ,PFZ(3),FZ,DFZ |
43 |
|
|
INTEGER INDEX(3) !point index |
44 |
|
|
|
45 |
|
|
COMMON /PAWCR4/ INDEX,PFX,FX,DFX,PFY,FY,DFY,PFZ,FZ,DFZ |
46 |
|
|
|
47 |
|
|
|
48 |
|
|
c------------------------------------------------------------------------ |
49 |
|
|
c |
50 |
|
|
c *** FIRST MAP *** |
51 |
|
|
c |
52 |
|
|
c------------------------------------------------------------------------ |
53 |
|
|
|
54 |
|
|
c------------------------------------------------------------------------ |
55 |
|
|
c |
56 |
|
|
c initialization and map file opening |
57 |
|
|
c |
58 |
|
|
c------------------------------------------------------------------------ |
59 |
|
|
|
60 |
|
|
c print*,' ' |
61 |
|
|
c print*,' ' |
62 |
|
|
|
63 |
|
|
Bmap_file='measure_n3_290302.rz' |
64 |
|
|
|
65 |
|
|
c$$$ cmd1='cp $TRK_GRND/source/magnet/' |
66 |
|
|
c$$$ $ //Bmap_file(1:LNBLNK(Bmap_file))//' .' |
67 |
|
|
c$$$ call system(cmd1) |
68 |
|
|
|
69 |
|
|
c opens magnetic field map first file |
70 |
|
|
print *,'Opening file: ',Bmap_file |
71 |
|
|
call HROPEN |
72 |
|
|
$ (lun_Bmap_file,'Bmap','./bin-aux/'//Bmap_file,'P',1024,istat) |
73 |
|
|
if(istat.ne.0) goto 21 |
74 |
|
|
|
75 |
|
|
|
76 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
77 |
|
|
|
78 |
|
|
c call HPRNTU(ntpl_Bmap) |
79 |
|
|
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
80 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
81 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
82 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
83 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
84 |
|
|
|
85 |
|
|
|
86 |
|
|
c------------------------------------------------------------------------ |
87 |
|
|
c |
88 |
|
|
c reads events and fills variables |
89 |
|
|
c |
90 |
|
|
c------------------------------------------------------------------------ |
91 |
|
|
|
92 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
93 |
|
|
|
94 |
|
|
c initializes measurement grid edges |
95 |
|
|
do ic=1,3 |
96 |
|
|
px1max(ic)=0. |
97 |
|
|
px1min(ic)=0. |
98 |
|
|
py1max(ic)=0. |
99 |
|
|
py1min(ic)=0. |
100 |
|
|
pz1max(ic)=0. |
101 |
|
|
pz1min(ic)=0. |
102 |
|
|
enddo |
103 |
|
|
|
104 |
|
|
|
105 |
|
|
do iev=1,iemax !event loop |
106 |
|
|
|
107 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
108 |
|
|
if(ierr.ne.0) goto 22 |
109 |
|
|
|
110 |
|
|
c the output consists of matrices for coordinates, B components values |
111 |
|
|
c and errors: |
112 |
|
|
c e.g. px1(4,2) = X coordinate of the point with index = 4 along X, |
113 |
|
|
c in which By (=2) component has been measured |
114 |
|
|
c e.g. b1(3,23,4,1) = Bx (=1) component value, measured in the point with |
115 |
|
|
c indexes = 3,23,4 along X, Y and Z |
116 |
|
|
|
117 |
|
|
c Bx component |
118 |
|
|
px1(index(1),1) = pfx(1) |
119 |
|
|
if(px1(index(1),1).lt.px1min(1)) px1min(1)=px1(index(1),1) |
120 |
|
|
if(px1(index(1),1).gt.px1max(1)) px1max(1)=px1(index(1),1) |
121 |
|
|
py1(index(2),1) = pfx(2) |
122 |
|
|
if(py1(index(2),1).lt.py1min(1)) py1min(1)=py1(index(2),1) |
123 |
|
|
if(py1(index(2),1).gt.py1max(1)) py1max(1)=py1(index(2),1) |
124 |
|
|
pz1(index(3),1) = pfx(3) |
125 |
|
|
if(pz1(index(3),1).lt.pz1min(1)) pz1min(1)=pz1(index(3),1) |
126 |
|
|
if(pz1(index(3),1).gt.pz1max(1)) pz1max(1)=pz1(index(3),1) |
127 |
|
|
|
128 |
|
|
b1(index(1),index(2),index(3),1) = fx |
129 |
|
|
|
130 |
|
|
|
131 |
|
|
c By component |
132 |
|
|
px1(index(1),2) = pfy(1) |
133 |
|
|
if(px1(index(1),2).lt.px1min(2)) px1min(2)=px1(index(1),2) |
134 |
|
|
if(px1(index(1),2).gt.px1max(2)) px1max(2)=px1(index(1),2) |
135 |
|
|
py1(index(2),2) = pfy(2) |
136 |
|
|
if(py1(index(2),2).lt.py1min(2)) py1min(2)=py1(index(2),2) |
137 |
|
|
if(py1(index(2),2).gt.py1max(2)) py1max(2)=py1(index(2),2) |
138 |
|
|
pz1(index(3),2) = pfy(3) |
139 |
|
|
if(pz1(index(3),2).lt.pz1min(2)) pz1min(2)=pz1(index(3),2) |
140 |
|
|
if(pz1(index(3),2).gt.pz1max(2)) pz1max(2)=pz1(index(3),2) |
141 |
|
|
|
142 |
|
|
b1(index(1),index(2),index(3),2) = fy |
143 |
|
|
|
144 |
|
|
|
145 |
|
|
c Bz component |
146 |
|
|
px1(index(1),3) = pfz(1) |
147 |
|
|
if(px1(index(1),3).lt.px1min(3)) px1min(3)=px1(index(1),3) |
148 |
|
|
if(px1(index(1),3).gt.px1max(3)) px1max(3)=px1(index(1),3) |
149 |
|
|
py1(index(2),3) = pfz(2) |
150 |
|
|
if(py1(index(2),3).lt.py1min(3)) py1min(3)=py1(index(2),3) |
151 |
|
|
if(py1(index(2),3).gt.py1max(3)) py1max(3)=py1(index(2),3) |
152 |
|
|
pz1(index(3),3) = pfz(3) |
153 |
|
|
if(pz1(index(3),3).lt.pz1min(3)) pz1min(3)=pz1(index(3),3) |
154 |
|
|
if(pz1(index(3),3).gt.pz1max(3)) pz1max(3)=pz1(index(3),3) |
155 |
|
|
|
156 |
|
|
b1(index(1),index(2),index(3),3) = fz |
157 |
|
|
|
158 |
|
|
enddo |
159 |
|
|
|
160 |
|
|
|
161 |
|
|
c------------------------------------------------------------------------ |
162 |
|
|
c |
163 |
|
|
c closes files |
164 |
|
|
c |
165 |
|
|
c------------------------------------------------------------------------ |
166 |
|
|
|
167 |
|
|
call HREND('Bmap') |
168 |
|
|
close(lun_Bmap_file) |
169 |
|
|
c$$$ cmd2='rm -f ' |
170 |
|
|
c$$$ $ //Bmap_file(1:LNBLNK(Bmap_file)) |
171 |
|
|
c$$$ call system(cmd2) |
172 |
|
|
c$$$ |
173 |
|
|
|
174 |
|
|
|
175 |
|
|
c------------------------------------------------------------------------ |
176 |
|
|
c |
177 |
|
|
c *** SECOND MAP *** |
178 |
|
|
c |
179 |
|
|
c------------------------------------------------------------------------ |
180 |
|
|
|
181 |
|
|
c------------------------------------------------------------------------ |
182 |
|
|
c |
183 |
|
|
c initialization and map file opening |
184 |
|
|
c |
185 |
|
|
c------------------------------------------------------------------------ |
186 |
|
|
|
187 |
|
|
c print*,' ' |
188 |
|
|
c print*,' ' |
189 |
|
|
|
190 |
|
|
Bmap_file='measure_n4_110402_corrected.rz' |
191 |
|
|
c$$$ cmd1='cp $TRK_GRND/source/magnet/' |
192 |
|
|
c$$$ $ //Bmap_file(1:LNBLNK(Bmap_file))//' .' |
193 |
|
|
c$$$ call system(cmd1) |
194 |
|
|
|
195 |
|
|
c opens magnetic field map first file |
196 |
|
|
print *,'Opening file: ',Bmap_file |
197 |
|
|
call HROPEN |
198 |
|
|
$ (lun_Bmap_file,'Bmap','./bin-aux/'//Bmap_file,'P',1024,istat) |
199 |
|
|
if(istat.ne.0) goto 21 |
200 |
|
|
|
201 |
|
|
|
202 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
203 |
|
|
|
204 |
|
|
c call HPRNTU(ntpl_Bmap) |
205 |
|
|
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
206 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
207 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
208 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
209 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
210 |
|
|
|
211 |
|
|
|
212 |
|
|
c------------------------------------------------------------------------ |
213 |
|
|
c |
214 |
|
|
c reads events and fills variables |
215 |
|
|
c |
216 |
|
|
c------------------------------------------------------------------------ |
217 |
|
|
|
218 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
219 |
|
|
|
220 |
|
|
do ic=1,3 !grid edges |
221 |
|
|
px2max(ic)=0. |
222 |
|
|
px2min(ic)=0. |
223 |
|
|
py2max(ic)=0. |
224 |
|
|
py2min(ic)=0. |
225 |
|
|
pz2max(ic)=0. |
226 |
|
|
pz2min(ic)=0. |
227 |
|
|
enddo |
228 |
|
|
|
229 |
|
|
|
230 |
|
|
do iev=1,iemax !event loop |
231 |
|
|
|
232 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
233 |
|
|
if(ierr.ne.0) goto 22 |
234 |
|
|
|
235 |
|
|
c the output consists of matrices for coordinates, B components values |
236 |
|
|
c and errors: |
237 |
|
|
c e.g. px(4,2) = X coordinate of the point with index = 4 along X, |
238 |
|
|
c in which By (=2) component has been measured |
239 |
|
|
c e.g. b(3,23,4,1) = Bx (=1) component value, measured in the point with |
240 |
|
|
c indexes = 3,23,4 along X, Y and Z |
241 |
|
|
|
242 |
|
|
c Bx component |
243 |
|
|
px2(index(1),1) = pfx(1) |
244 |
|
|
if(px2(index(1),1).lt.px2min(1)) px2min(1)=px2(index(1),1) |
245 |
|
|
if(px2(index(1),1).gt.px2max(1)) px2max(1)=px2(index(1),1) |
246 |
|
|
py2(index(2),1) = pfx(2) |
247 |
|
|
if(py2(index(2),1).lt.py2min(1)) py2min(1)=py2(index(2),1) |
248 |
|
|
if(py2(index(2),1).gt.py2max(1)) py2max(1)=py2(index(2),1) |
249 |
|
|
pz2(index(3),1) = pfx(3) |
250 |
|
|
if(pz2(index(3),1).lt.pz2min(1)) pz2min(1)=pz2(index(3),1) |
251 |
|
|
if(pz2(index(3),1).gt.pz2max(1)) pz2max(1)=pz2(index(3),1) |
252 |
|
|
|
253 |
|
|
b2(index(1),index(2),index(3),1) = fx |
254 |
|
|
|
255 |
|
|
|
256 |
|
|
c By component |
257 |
|
|
px2(index(1),2) = pfy(1) |
258 |
|
|
if(px2(index(1),2).lt.px2min(2)) px2min(2)=px2(index(1),2) |
259 |
|
|
if(px2(index(1),2).gt.px2max(2)) px2max(2)=px2(index(1),2) |
260 |
|
|
py2(index(2),2) = pfy(2) |
261 |
|
|
if(py2(index(2),2).lt.py2min(2)) py2min(2)=py2(index(2),2) |
262 |
|
|
if(py2(index(2),2).gt.py2max(2)) py2max(2)=py2(index(2),2) |
263 |
|
|
pz2(index(3),2) = pfy(3) |
264 |
|
|
if(pz2(index(3),2).lt.pz2min(2)) pz2min(2)=pz2(index(3),2) |
265 |
|
|
if(pz2(index(3),2).gt.pz2max(2)) pz2max(2)=pz2(index(3),2) |
266 |
|
|
|
267 |
|
|
b2(index(1),index(2),index(3),2) = fy |
268 |
|
|
|
269 |
|
|
|
270 |
|
|
c Bz component |
271 |
|
|
px2(index(1),3) = pfz(1) |
272 |
|
|
if(px2(index(1),3).lt.px2min(3)) px2min(3)=px2(index(1),3) |
273 |
|
|
if(px2(index(1),3).gt.px2max(3)) px2max(3)=px2(index(1),3) |
274 |
|
|
py2(index(2),3) = pfz(2) |
275 |
|
|
if(py2(index(2),3).lt.py2min(3)) py2min(3)=py2(index(2),3) |
276 |
|
|
if(py2(index(2),3).gt.py2max(3)) py2max(3)=py2(index(2),3) |
277 |
|
|
pz2(index(3),3) = pfz(3) |
278 |
|
|
if(pz2(index(3),3).lt.pz2min(3)) pz2min(3)=pz2(index(3),3) |
279 |
|
|
if(pz2(index(3),3).gt.pz2max(3)) pz2max(3)=pz2(index(3),3) |
280 |
|
|
|
281 |
|
|
b2(index(1),index(2),index(3),3) = fz |
282 |
|
|
|
283 |
|
|
enddo |
284 |
|
|
|
285 |
|
|
|
286 |
|
|
c------------------------------------------------------------------------ |
287 |
|
|
c |
288 |
|
|
c closes files |
289 |
|
|
c |
290 |
|
|
c------------------------------------------------------------------------ |
291 |
|
|
|
292 |
|
|
call HREND('Bmap') |
293 |
|
|
close(lun_Bmap_file) |
294 |
|
|
c$$$ cmd2='rm -f ' |
295 |
|
|
c$$$ $ //Bmap_file(1:LNBLNK(Bmap_file)) |
296 |
|
|
c$$$ call system(cmd2) |
297 |
|
|
|
298 |
|
|
|
299 |
|
|
c------------------------------------------------------------------------ |
300 |
|
|
c |
301 |
|
|
c no error exit |
302 |
|
|
c |
303 |
|
|
c------------------------------------------------------------------------ |
304 |
|
|
|
305 |
|
|
c$$$ print*,' ' |
306 |
|
|
c$$$ print*,'MAGNETIC FIELD SUCCESSFULLY READ' |
307 |
|
|
c$$$ print*,' ' |
308 |
|
|
c$$$ print*,' ' |
309 |
|
|
|
310 |
|
|
goto 9000 !happy ending |
311 |
|
|
|
312 |
|
|
c------------------------------------------------------------------------ |
313 |
|
|
c |
314 |
|
|
c magnetic field map file opening error |
315 |
|
|
c |
316 |
|
|
c------------------------------------------------------------------------ |
317 |
|
|
|
318 |
|
|
21 continue |
319 |
|
|
|
320 |
|
|
print*,' ' |
321 |
|
|
print*,'read_B_inner: ERROR OPENING MAGNETIC FIELD MAP FILE: ' |
322 |
|
|
$ ,Bmap_file |
323 |
|
|
print*,' ' |
324 |
|
|
print*,' ' |
325 |
|
|
|
326 |
|
|
goto 9000 !the end |
327 |
|
|
|
328 |
|
|
|
329 |
|
|
c------------------------------------------------------------------------ |
330 |
|
|
c |
331 |
|
|
c ntuple event reading error |
332 |
|
|
c |
333 |
|
|
c------------------------------------------------------------------------ |
334 |
|
|
|
335 |
|
|
22 continue |
336 |
|
|
|
337 |
|
|
print*,' ' |
338 |
|
|
print*,'read_B_inner: ERROR WHILE READING NTUPLE, AT EVENT |
339 |
|
|
$ : ',iev |
340 |
|
|
print*,' ' |
341 |
|
|
print*,' ' |
342 |
|
|
|
343 |
|
|
goto 9000 !the end |
344 |
|
|
|
345 |
|
|
|
346 |
|
|
c------------------------------------------------------------------------ |
347 |
|
|
c |
348 |
|
|
c exit |
349 |
|
|
c |
350 |
|
|
c------------------------------------------------------------------------ |
351 |
|
|
|
352 |
|
|
9000 continue |
353 |
|
|
|
354 |
|
|
return |
355 |
|
|
end |