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 |
|
|
call HPRNTU(ntpl_Bmap) |
60 |
|
|
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 |
|
|
|
153 |
|
|
c------------------------------------------------------------------------ |
154 |
|
|
c |
155 |
|
|
c *** SECOND MAP *** |
156 |
|
|
c |
157 |
|
|
c------------------------------------------------------------------------ |
158 |
|
|
|
159 |
|
|
|
160 |
|
|
CALL HCDIR('//FIELD2',' ') |
161 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
162 |
|
|
|
163 |
|
|
call HPRNTU(ntpl_Bmap) |
164 |
|
|
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
165 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
166 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
167 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
168 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
169 |
|
|
|
170 |
|
|
|
171 |
|
|
c------------------------------------------------------------------------ |
172 |
|
|
c |
173 |
|
|
c reads events and fills variables |
174 |
|
|
c |
175 |
|
|
c------------------------------------------------------------------------ |
176 |
|
|
|
177 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
178 |
|
|
|
179 |
|
|
do ic=1,3 !grid edges |
180 |
|
|
px2max(ic)=0. |
181 |
|
|
px2min(ic)=0. |
182 |
|
|
py2max(ic)=0. |
183 |
|
|
py2min(ic)=0. |
184 |
|
|
pz2max(ic)=0. |
185 |
|
|
pz2min(ic)=0. |
186 |
|
|
enddo |
187 |
|
|
|
188 |
|
|
|
189 |
|
|
do iev=1,iemax !event loop |
190 |
|
|
|
191 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
192 |
|
|
if(ierr.ne.0) goto 22 |
193 |
|
|
|
194 |
|
|
c the output consists of matrices for coordinates, B components values |
195 |
|
|
c and errors: |
196 |
|
|
c e.g. px(4,2) = X coordinate of the point with index = 4 along X, |
197 |
|
|
c in which By (=2) component has been measured |
198 |
|
|
c e.g. b(3,23,4,1) = Bx (=1) component value, measured in the point with |
199 |
|
|
c indexes = 3,23,4 along X, Y and Z |
200 |
|
|
|
201 |
|
|
c Bx component |
202 |
|
|
px2(index(1),1) = DBLE(pfx(1)) |
203 |
|
|
if(px2(index(1),1).lt.px2min(1)) px2min(1)=px2(index(1),1) |
204 |
|
|
if(px2(index(1),1).gt.px2max(1)) px2max(1)=px2(index(1),1) |
205 |
|
|
py2(index(2),1) = DBLE(pfx(2)) |
206 |
|
|
if(py2(index(2),1).lt.py2min(1)) py2min(1)=py2(index(2),1) |
207 |
|
|
if(py2(index(2),1).gt.py2max(1)) py2max(1)=py2(index(2),1) |
208 |
|
|
pz2(index(3),1) = DBLE(pfx(3)) |
209 |
|
|
if(pz2(index(3),1).lt.pz2min(1)) pz2min(1)=pz2(index(3),1) |
210 |
|
|
if(pz2(index(3),1).gt.pz2max(1)) pz2max(1)=pz2(index(3),1) |
211 |
|
|
|
212 |
|
|
b2(index(1),index(2),index(3),1) = DBLE(ffx) |
213 |
|
|
|
214 |
|
|
|
215 |
|
|
c By component |
216 |
|
|
px2(index(1),2) = DBLE(pfy(1)) |
217 |
|
|
if(px2(index(1),2).lt.px2min(2)) px2min(2)=px2(index(1),2) |
218 |
|
|
if(px2(index(1),2).gt.px2max(2)) px2max(2)=px2(index(1),2) |
219 |
|
|
py2(index(2),2) = DBLE(pfy(2)) |
220 |
|
|
if(py2(index(2),2).lt.py2min(2)) py2min(2)=py2(index(2),2) |
221 |
|
|
if(py2(index(2),2).gt.py2max(2)) py2max(2)=py2(index(2),2) |
222 |
|
|
pz2(index(3),2) = DBLE(pfy(3)) |
223 |
|
|
if(pz2(index(3),2).lt.pz2min(2)) pz2min(2)=pz2(index(3),2) |
224 |
|
|
if(pz2(index(3),2).gt.pz2max(2)) pz2max(2)=pz2(index(3),2) |
225 |
|
|
|
226 |
|
|
b2(index(1),index(2),index(3),2) = DBLE(ffy) |
227 |
|
|
|
228 |
|
|
|
229 |
|
|
c Bz component |
230 |
|
|
px2(index(1),3) = DBLE(pfz(1)) |
231 |
|
|
if(px2(index(1),3).lt.px2min(3)) px2min(3)=px2(index(1),3) |
232 |
|
|
if(px2(index(1),3).gt.px2max(3)) px2max(3)=px2(index(1),3) |
233 |
|
|
py2(index(2),3) = DBLE(pfz(2)) |
234 |
|
|
if(py2(index(2),3).lt.py2min(3)) py2min(3)=py2(index(2),3) |
235 |
|
|
if(py2(index(2),3).gt.py2max(3)) py2max(3)=py2(index(2),3) |
236 |
|
|
pz2(index(3),3) = DBLE(pfz(3)) |
237 |
|
|
if(pz2(index(3),3).lt.pz2min(3)) pz2min(3)=pz2(index(3),3) |
238 |
|
|
if(pz2(index(3),3).gt.pz2max(3)) pz2max(3)=pz2(index(3),3) |
239 |
|
|
|
240 |
|
|
b2(index(1),index(2),index(3),3) = DBLE(ffz) |
241 |
|
|
|
242 |
|
|
enddo |
243 |
|
|
|
244 |
|
|
c------------------------------------------------------------------------ |
245 |
|
|
c |
246 |
|
|
c closes files |
247 |
|
|
c |
248 |
|
|
c------------------------------------------------------------------------ |
249 |
|
|
|
250 |
|
|
call HREND('FIELD2') |
251 |
|
|
|
252 |
|
|
|
253 |
|
|
c------------------------------------------------------------------------ |
254 |
|
|
c |
255 |
|
|
c no error exit |
256 |
|
|
c |
257 |
|
|
c------------------------------------------------------------------------ |
258 |
|
|
|
259 |
|
|
print*,' ' |
260 |
|
|
print*,'MAGNETIC FIELD SUCCESSFULLY READ' |
261 |
|
|
print*,' ' |
262 |
|
|
print*,' ' |
263 |
|
|
|
264 |
|
|
goto 9000 !happy ending |
265 |
|
|
|
266 |
|
|
c------------------------------------------------------------------------ |
267 |
|
|
c |
268 |
|
|
c magnetic field map file opening error |
269 |
|
|
c |
270 |
|
|
c------------------------------------------------------------------------ |
271 |
|
|
|
272 |
|
|
|
273 |
|
|
|
274 |
|
|
c------------------------------------------------------------------------ |
275 |
|
|
c |
276 |
|
|
c ntuple event reading error |
277 |
|
|
c |
278 |
|
|
c------------------------------------------------------------------------ |
279 |
|
|
|
280 |
|
|
22 continue |
281 |
|
|
|
282 |
|
|
print*,' ' |
283 |
|
|
print*,'read_B_inner: ERROR WHILE READING NTUPLE, AT EVENT |
284 |
|
|
$ : ',iev |
285 |
|
|
print*,' ' |
286 |
|
|
print*,' ' |
287 |
|
|
|
288 |
|
|
goto 9000 !the end |
289 |
|
|
|
290 |
|
|
|
291 |
|
|
c------------------------------------------------------------------------ |
292 |
|
|
c |
293 |
|
|
c exit |
294 |
|
|
c |
295 |
|
|
c------------------------------------------------------------------------ |
296 |
|
|
|
297 |
|
|
9000 continue |
298 |
|
|
|
299 |
|
|
return |
300 |
|
|
end |