1 |
cafagna |
3.1 |
SUBROUTINE GPRBI |
2 |
|
|
************************************************************************* |
3 |
|
|
* |
4 |
|
|
* Subroutine gprbi.f (DERIVED FROM read_B_inner routine from |
5 |
|
|
* tracker software analysis) |
6 |
|
|
* CALLED BY gpdat.F |
7 |
|
|
* |
8 |
|
|
* |
9 |
|
|
* it reads from rz files the two magnetic field maps taken inside the |
10 |
|
|
* spectrometer cavity and fills the variables in common_B_inner.f |
11 |
|
|
* |
12 |
|
|
* needs: |
13 |
|
|
* - .rz map files in containing coordinates of measured points, Bx, By |
14 |
|
|
* and Bz components + errors |
15 |
|
|
* |
16 |
|
|
* output variables: (see ....) |
17 |
|
|
* - px#(nx,3) with #=1,2 for the 2 maps |
18 |
|
|
* - py#(ny,3) |
19 |
|
|
* - pz#(nz,3) |
20 |
|
|
* - b#(nx,ny,nz,3) |
21 |
|
|
* |
22 |
|
|
************************************************************************* |
23 |
|
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
24 |
|
|
#include "gpfield.inc" |
25 |
|
|
|
26 |
|
|
|
27 |
|
|
|
28 |
|
|
c------------------------------------------------------------------------ |
29 |
|
|
c |
30 |
|
|
c local variables |
31 |
|
|
c |
32 |
|
|
c------------------------------------------------------------------------ |
33 |
|
|
|
34 |
|
|
parameter (ntpl_Bmap=20) !ntuple identifier |
35 |
|
|
|
36 |
|
|
REAL PFX(3),FFX,DFX, !Bx field component coordinates in m, value and error in T |
37 |
|
|
$ PFY(3),FFY,DFY |
38 |
|
|
$ ,PFZ(3),FFZ,DFZ |
39 |
|
|
INTEGER INDEX(3) !point index |
40 |
|
|
|
41 |
|
|
COMMON /PAWCR4/ INDEX,PFX,FFX,DFX,PFY,FFY,DFY,PFZ,FFZ,DFZ |
42 |
|
|
|
43 |
|
|
|
44 |
|
|
c------------------------------------------------------------------------ |
45 |
|
|
c |
46 |
|
|
c *** FIRST MAP *** |
47 |
|
|
c |
48 |
|
|
c------------------------------------------------------------------------ |
49 |
|
|
|
50 |
|
|
c------------------------------------------------------------------------ |
51 |
|
|
c |
52 |
|
|
c initialization and map file opening |
53 |
|
|
c |
54 |
|
|
c------------------------------------------------------------------------ |
55 |
|
|
|
56 |
|
|
CALL HCDIR('//FIELD1',' ') |
57 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
58 |
|
|
|
59 |
bottai |
3.2 |
|
60 |
cafagna |
3.1 |
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
61 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
62 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
63 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
64 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
65 |
|
|
|
66 |
|
|
|
67 |
|
|
c------------------------------------------------------------------------ |
68 |
|
|
c |
69 |
|
|
c reads events and fills variables |
70 |
|
|
c |
71 |
|
|
c------------------------------------------------------------------------ |
72 |
|
|
|
73 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
74 |
|
|
|
75 |
|
|
c initializes measurement grid edges |
76 |
|
|
do ic=1,3 |
77 |
|
|
px1max(ic)=0. |
78 |
|
|
px1min(ic)=0. |
79 |
|
|
py1max(ic)=0. |
80 |
|
|
py1min(ic)=0. |
81 |
|
|
pz1max(ic)=0. |
82 |
|
|
pz1min(ic)=0. |
83 |
|
|
enddo |
84 |
|
|
|
85 |
|
|
|
86 |
|
|
do iev=1,iemax !event loop |
87 |
|
|
|
88 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
89 |
|
|
if(ierr.ne.0) goto 22 |
90 |
|
|
|
91 |
|
|
c the output consists of matrices for coordinates, B components values |
92 |
|
|
c and errors: |
93 |
|
|
c e.g. px1(4,2) = X coordinate of the point with index = 4 along X, |
94 |
|
|
c in which By (=2) component has been measured |
95 |
|
|
c e.g. b1(3,23,4,1) = Bx (=1) component value, measured in the point with |
96 |
|
|
c indexes = 3,23,4 along X, Y and Z |
97 |
|
|
|
98 |
|
|
c Bx component |
99 |
|
|
px1(index(1),1) = DBLE(pfx(1)) |
100 |
|
|
if(px1(index(1),1).lt.px1min(1)) px1min(1)=px1(index(1),1) |
101 |
|
|
if(px1(index(1),1).gt.px1max(1)) px1max(1)=px1(index(1),1) |
102 |
|
|
py1(index(2),1) = DBLE(pfx(2)) |
103 |
|
|
if(py1(index(2),1).lt.py1min(1)) py1min(1)=py1(index(2),1) |
104 |
|
|
if(py1(index(2),1).gt.py1max(1)) py1max(1)=py1(index(2),1) |
105 |
|
|
pz1(index(3),1) = DBLE(pfx(3)) |
106 |
|
|
if(pz1(index(3),1).lt.pz1min(1)) pz1min(1)=pz1(index(3),1) |
107 |
|
|
if(pz1(index(3),1).gt.pz1max(1)) pz1max(1)=pz1(index(3),1) |
108 |
|
|
|
109 |
|
|
b1(index(1),index(2),index(3),1) = DBLE(ffx) |
110 |
|
|
|
111 |
|
|
|
112 |
|
|
c By component |
113 |
|
|
px1(index(1),2) = DBLE(pfy(1)) |
114 |
|
|
if(px1(index(1),2).lt.px1min(2)) px1min(2)=px1(index(1),2) |
115 |
|
|
if(px1(index(1),2).gt.px1max(2)) px1max(2)=px1(index(1),2) |
116 |
|
|
py1(index(2),2) = DBLE(pfy(2)) |
117 |
|
|
if(py1(index(2),2).lt.py1min(2)) py1min(2)=py1(index(2),2) |
118 |
|
|
if(py1(index(2),2).gt.py1max(2)) py1max(2)=py1(index(2),2) |
119 |
|
|
pz1(index(3),2) = DBLE(pfy(3)) |
120 |
|
|
if(pz1(index(3),2).lt.pz1min(2)) pz1min(2)=pz1(index(3),2) |
121 |
|
|
if(pz1(index(3),2).gt.pz1max(2)) pz1max(2)=pz1(index(3),2) |
122 |
|
|
|
123 |
|
|
b1(index(1),index(2),index(3),2) = DBLE(ffy) |
124 |
|
|
|
125 |
|
|
|
126 |
|
|
c Bz component |
127 |
|
|
px1(index(1),3) = DBLE(pfz(1)) |
128 |
|
|
if(px1(index(1),3).lt.px1min(3)) px1min(3)=px1(index(1),3) |
129 |
|
|
if(px1(index(1),3).gt.px1max(3)) px1max(3)=px1(index(1),3) |
130 |
|
|
py1(index(2),3) = DBLE(pfz(2)) |
131 |
|
|
if(py1(index(2),3).lt.py1min(3)) py1min(3)=py1(index(2),3) |
132 |
|
|
if(py1(index(2),3).gt.py1max(3)) py1max(3)=py1(index(2),3) |
133 |
|
|
pz1(index(3),3) = DBLE(pfz(3)) |
134 |
|
|
if(pz1(index(3),3).lt.pz1min(3)) pz1min(3)=pz1(index(3),3) |
135 |
|
|
if(pz1(index(3),3).gt.pz1max(3)) pz1max(3)=pz1(index(3),3) |
136 |
|
|
|
137 |
|
|
b1(index(1),index(2),index(3),3) = DBLE(ffz) |
138 |
|
|
|
139 |
|
|
enddo |
140 |
|
|
|
141 |
|
|
|
142 |
|
|
c------------------------------------------------------------------------ |
143 |
|
|
c |
144 |
|
|
c closes files |
145 |
|
|
c |
146 |
|
|
c------------------------------------------------------------------------ |
147 |
|
|
|
148 |
|
|
call HREND('FIELD1') |
149 |
|
|
|
150 |
|
|
|
151 |
|
|
|
152 |
|
|
c------------------------------------------------------------------------ |
153 |
|
|
c |
154 |
|
|
c *** SECOND MAP *** |
155 |
|
|
c |
156 |
|
|
c------------------------------------------------------------------------ |
157 |
|
|
|
158 |
|
|
|
159 |
|
|
CALL HCDIR('//FIELD2',' ') |
160 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
161 |
|
|
|
162 |
bottai |
3.2 |
|
163 |
cafagna |
3.1 |
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
164 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
165 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
166 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
167 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
168 |
|
|
|
169 |
|
|
|
170 |
|
|
c------------------------------------------------------------------------ |
171 |
|
|
c |
172 |
|
|
c reads events and fills variables |
173 |
|
|
c |
174 |
|
|
c------------------------------------------------------------------------ |
175 |
|
|
|
176 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
177 |
|
|
|
178 |
|
|
do ic=1,3 !grid edges |
179 |
|
|
px2max(ic)=0. |
180 |
|
|
px2min(ic)=0. |
181 |
|
|
py2max(ic)=0. |
182 |
|
|
py2min(ic)=0. |
183 |
|
|
pz2max(ic)=0. |
184 |
|
|
pz2min(ic)=0. |
185 |
|
|
enddo |
186 |
|
|
|
187 |
|
|
|
188 |
|
|
do iev=1,iemax !event loop |
189 |
|
|
|
190 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
191 |
|
|
if(ierr.ne.0) goto 22 |
192 |
|
|
|
193 |
|
|
c the output consists of matrices for coordinates, B components values |
194 |
|
|
c and errors: |
195 |
|
|
c e.g. px(4,2) = X coordinate of the point with index = 4 along X, |
196 |
|
|
c in which By (=2) component has been measured |
197 |
|
|
c e.g. b(3,23,4,1) = Bx (=1) component value, measured in the point with |
198 |
|
|
c indexes = 3,23,4 along X, Y and Z |
199 |
|
|
|
200 |
|
|
c Bx component |
201 |
|
|
px2(index(1),1) = DBLE(pfx(1)) |
202 |
|
|
if(px2(index(1),1).lt.px2min(1)) px2min(1)=px2(index(1),1) |
203 |
|
|
if(px2(index(1),1).gt.px2max(1)) px2max(1)=px2(index(1),1) |
204 |
|
|
py2(index(2),1) = DBLE(pfx(2)) |
205 |
|
|
if(py2(index(2),1).lt.py2min(1)) py2min(1)=py2(index(2),1) |
206 |
|
|
if(py2(index(2),1).gt.py2max(1)) py2max(1)=py2(index(2),1) |
207 |
|
|
pz2(index(3),1) = DBLE(pfx(3)) |
208 |
|
|
if(pz2(index(3),1).lt.pz2min(1)) pz2min(1)=pz2(index(3),1) |
209 |
|
|
if(pz2(index(3),1).gt.pz2max(1)) pz2max(1)=pz2(index(3),1) |
210 |
|
|
|
211 |
|
|
b2(index(1),index(2),index(3),1) = DBLE(ffx) |
212 |
|
|
|
213 |
|
|
|
214 |
|
|
c By component |
215 |
|
|
px2(index(1),2) = DBLE(pfy(1)) |
216 |
|
|
if(px2(index(1),2).lt.px2min(2)) px2min(2)=px2(index(1),2) |
217 |
|
|
if(px2(index(1),2).gt.px2max(2)) px2max(2)=px2(index(1),2) |
218 |
|
|
py2(index(2),2) = DBLE(pfy(2)) |
219 |
|
|
if(py2(index(2),2).lt.py2min(2)) py2min(2)=py2(index(2),2) |
220 |
|
|
if(py2(index(2),2).gt.py2max(2)) py2max(2)=py2(index(2),2) |
221 |
|
|
pz2(index(3),2) = DBLE(pfy(3)) |
222 |
|
|
if(pz2(index(3),2).lt.pz2min(2)) pz2min(2)=pz2(index(3),2) |
223 |
|
|
if(pz2(index(3),2).gt.pz2max(2)) pz2max(2)=pz2(index(3),2) |
224 |
|
|
|
225 |
|
|
b2(index(1),index(2),index(3),2) = DBLE(ffy) |
226 |
|
|
|
227 |
|
|
|
228 |
|
|
c Bz component |
229 |
|
|
px2(index(1),3) = DBLE(pfz(1)) |
230 |
|
|
if(px2(index(1),3).lt.px2min(3)) px2min(3)=px2(index(1),3) |
231 |
|
|
if(px2(index(1),3).gt.px2max(3)) px2max(3)=px2(index(1),3) |
232 |
|
|
py2(index(2),3) = DBLE(pfz(2)) |
233 |
|
|
if(py2(index(2),3).lt.py2min(3)) py2min(3)=py2(index(2),3) |
234 |
|
|
if(py2(index(2),3).gt.py2max(3)) py2max(3)=py2(index(2),3) |
235 |
|
|
pz2(index(3),3) = DBLE(pfz(3)) |
236 |
|
|
if(pz2(index(3),3).lt.pz2min(3)) pz2min(3)=pz2(index(3),3) |
237 |
|
|
if(pz2(index(3),3).gt.pz2max(3)) pz2max(3)=pz2(index(3),3) |
238 |
|
|
|
239 |
|
|
b2(index(1),index(2),index(3),3) = DBLE(ffz) |
240 |
|
|
|
241 |
|
|
enddo |
242 |
|
|
|
243 |
|
|
c------------------------------------------------------------------------ |
244 |
|
|
c |
245 |
|
|
c closes files |
246 |
|
|
c |
247 |
|
|
c------------------------------------------------------------------------ |
248 |
|
|
|
249 |
|
|
call HREND('FIELD2') |
250 |
|
|
|
251 |
|
|
|
252 |
|
|
c------------------------------------------------------------------------ |
253 |
|
|
c |
254 |
|
|
c no error exit |
255 |
|
|
c |
256 |
|
|
c------------------------------------------------------------------------ |
257 |
|
|
|
258 |
|
|
print*,' ' |
259 |
bottai |
3.2 |
print*,'INTERNAL MAGNETIC FIELD SUCCESSFULLY READ' |
260 |
|
|
print*,' ' |
261 |
|
|
print*,' ' |
262 |
|
|
|
263 |
|
|
|
264 |
|
|
c------------------------------------------------------------------------ |
265 |
|
|
c |
266 |
|
|
c *** EXTERNAL MAP *** |
267 |
|
|
c |
268 |
|
|
c------------------------------------------------------------------------ |
269 |
|
|
|
270 |
|
|
c------------------------------------------------------------------------ |
271 |
|
|
c |
272 |
|
|
c initialization and map file opening |
273 |
|
|
c |
274 |
|
|
c------------------------------------------------------------------------ |
275 |
|
|
|
276 |
|
|
|
277 |
|
|
CALL HCDIR('//FIELD3',' ') |
278 |
|
|
|
279 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
280 |
|
|
|
281 |
|
|
|
282 |
|
|
|
283 |
|
|
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
284 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
285 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
286 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
287 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
288 |
|
|
|
289 |
|
|
|
290 |
|
|
c------------------------------------------------------------------------ |
291 |
|
|
c |
292 |
|
|
c reads events and fills variables |
293 |
|
|
c |
294 |
|
|
c------------------------------------------------------------------------ |
295 |
|
|
|
296 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
297 |
|
|
|
298 |
|
|
c initializes measurement grid edges |
299 |
|
|
do ic=1,3 |
300 |
|
|
poxmax(ic)=0. |
301 |
|
|
poxmin(ic)=0. |
302 |
|
|
poymax(ic)=0. |
303 |
|
|
poymin(ic)=0. |
304 |
|
|
pozmax(ic)=0. |
305 |
|
|
pozmin(ic)=0. |
306 |
|
|
enddo |
307 |
|
|
|
308 |
|
|
|
309 |
|
|
do iev=1,iemax !event loop |
310 |
|
|
|
311 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
312 |
|
|
if(ierr.ne.0) goto 22 |
313 |
|
|
|
314 |
|
|
c the output consists of matrices for coordinates, B components values |
315 |
|
|
c and errors: |
316 |
|
|
c e.g. px1(4,2) = X coordinate of the point with index = 4 along X, |
317 |
|
|
c in which By (=2) component has been measured |
318 |
|
|
c e.g. b1(3,23,4,1) = Bx (=1) component value, measured in the point with |
319 |
|
|
c indexes = 3,23,4 along X, Y and Z |
320 |
|
|
|
321 |
|
|
c Bx component |
322 |
|
|
pox(index(1),1) = pfx(1) |
323 |
|
|
if(pox(index(1),1).lt.poxmin(1)) poxmin(1)=pox(index(1),1) |
324 |
|
|
if(pox(index(1),1).gt.poxmax(1)) poxmax(1)=pox(index(1),1) |
325 |
|
|
poy(index(2),1) = pfx(2) |
326 |
|
|
if(poy(index(2),1).lt.poymin(1)) poymin(1)=poy(index(2),1) |
327 |
|
|
if(poy(index(2),1).gt.poymax(1)) poymax(1)=poy(index(2),1) |
328 |
|
|
poz(index(3),1) = pfx(3) |
329 |
|
|
if(poz(index(3),1).lt.pozmin(1)) pozmin(1)=poz(index(3),1) |
330 |
|
|
if(poz(index(3),1).gt.pozmax(1)) pozmax(1)=poz(index(3),1) |
331 |
|
|
|
332 |
|
|
bo(index(1),index(2),index(3),1) = DBLE(ffx) |
333 |
|
|
|
334 |
|
|
|
335 |
|
|
c By component |
336 |
|
|
pox(index(1),2) = pfy(1) |
337 |
|
|
if(pox(index(1),2).lt.poxmin(2)) poxmin(2)=pox(index(1),2) |
338 |
|
|
if(pox(index(1),2).gt.poxmax(2)) poxmax(2)=pox(index(1),2) |
339 |
|
|
poy(index(2),2) = pfy(2) |
340 |
|
|
if(poy(index(2),2).lt.poymin(2)) poymin(2)=poy(index(2),2) |
341 |
|
|
if(poy(index(2),2).gt.poymax(2)) poymax(2)=poy(index(2),2) |
342 |
|
|
poz(index(3),2) = pfy(3) |
343 |
|
|
if(poz(index(3),2).lt.pozmin(2)) pozmin(2)=poz(index(3),2) |
344 |
|
|
if(poz(index(3),2).gt.pozmax(2)) pozmax(2)=poz(index(3),2) |
345 |
|
|
|
346 |
|
|
bo(index(1),index(2),index(3),2) = DBLE(ffy) |
347 |
|
|
|
348 |
|
|
|
349 |
|
|
c Bz component |
350 |
|
|
pox(index(1),3) = pfz(1) |
351 |
|
|
if(pox(index(1),3).lt.poxmin(3)) poxmin(3)=pox(index(1),3) |
352 |
|
|
if(pox(index(1),3).gt.poxmax(3)) poxmax(3)=pox(index(1),3) |
353 |
|
|
poy(index(2),3) = pfz(2) |
354 |
|
|
if(poy(index(2),3).lt.poymin(3)) poymin(3)=poy(index(2),3) |
355 |
|
|
if(poy(index(2),3).gt.poymax(3)) poymax(3)=poy(index(2),3) |
356 |
|
|
poz(index(3),3) = pfz(3) |
357 |
|
|
if(poz(index(3),3).lt.pozmin(3)) pozmin(3)=poz(index(3),3) |
358 |
|
|
if(poz(index(3),3).gt.pozmax(3)) pozmax(3)=poz(index(3),3) |
359 |
|
|
|
360 |
|
|
bo(index(1),index(2),index(3),3) = DBLE(ffz) |
361 |
|
|
|
362 |
|
|
enddo |
363 |
|
|
|
364 |
|
|
|
365 |
|
|
c------------------------------------------------------------------------ |
366 |
|
|
c |
367 |
|
|
c closes files |
368 |
|
|
c |
369 |
|
|
c------------------------------------------------------------------------ |
370 |
|
|
|
371 |
|
|
call HREND('FIELD3') |
372 |
|
|
|
373 |
|
|
c------------------------------------------------------------------------ |
374 |
|
|
c |
375 |
|
|
c no error exit |
376 |
|
|
c |
377 |
|
|
c------------------------------------------------------------------------ |
378 |
|
|
|
379 |
|
|
print*,' ' |
380 |
|
|
print*,'EXTERNAL MAGNETIC FIELD SUCCESSFULLY READ' |
381 |
cafagna |
3.1 |
print*,' ' |
382 |
|
|
print*,' ' |
383 |
bottai |
3.2 |
|
384 |
cafagna |
3.1 |
|
385 |
|
|
goto 9000 !happy ending |
386 |
|
|
|
387 |
|
|
c------------------------------------------------------------------------ |
388 |
|
|
c |
389 |
|
|
c magnetic field map file opening error |
390 |
|
|
c |
391 |
|
|
c------------------------------------------------------------------------ |
392 |
|
|
|
393 |
|
|
|
394 |
|
|
|
395 |
|
|
c------------------------------------------------------------------------ |
396 |
|
|
c |
397 |
|
|
c ntuple event reading error |
398 |
|
|
c |
399 |
|
|
c------------------------------------------------------------------------ |
400 |
|
|
|
401 |
|
|
22 continue |
402 |
|
|
|
403 |
|
|
print*,' ' |
404 |
bottai |
3.2 |
print*,' ERROR WHILE READING NTUPLE, AT EVENT |
405 |
cafagna |
3.1 |
$ : ',iev |
406 |
|
|
print*,' ' |
407 |
|
|
print*,' ' |
408 |
|
|
|
409 |
|
|
goto 9000 !the end |
410 |
|
|
|
411 |
|
|
|
412 |
|
|
c------------------------------------------------------------------------ |
413 |
|
|
c |
414 |
|
|
c exit |
415 |
|
|
c |
416 |
|
|
c------------------------------------------------------------------------ |
417 |
|
|
|
418 |
|
|
9000 continue |
419 |
|
|
|
420 |
|
|
return |
421 |
|
|
end |