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 |