1 |
mocchiut |
1.1 |
************************************************************************* |
2 |
|
|
* |
3 |
|
|
* Subroutine read_B.f |
4 |
|
|
* |
5 |
|
|
* it calls the subroutines which read the magnetic field maps for |
6 |
|
|
* the PAMELA spectrometer |
7 |
|
|
* |
8 |
|
|
* needs: |
9 |
|
|
* - ./read_B_inner.f inner map reading subroutine |
10 |
|
|
* - ./read_B_outer.f outer map reading subroutine |
11 |
|
|
* |
12 |
|
|
* to be called before ./inter_B.f (interpolation subroutine) |
13 |
|
|
* |
14 |
|
|
************************************************************************* |
15 |
|
|
|
16 |
|
|
subroutine readB(cpath) |
17 |
|
|
include 'common_c2f.f' |
18 |
pam-fi |
1.2 |
|
19 |
|
|
LOGICAL DEBUG |
20 |
|
|
LOGICAL VERBOSE |
21 |
|
|
LOGICAL WARNING |
22 |
|
|
COMMON/DBG/DEBUG,VERBOSE,WARNING |
23 |
|
|
SAVE/DBG/ |
24 |
|
|
|
25 |
mocchiut |
1.1 |
character*80 cpath |
26 |
|
|
character*80 ppath |
27 |
|
|
|
28 |
|
|
|
29 |
|
|
b_error=0 |
30 |
pam-fi |
1.3 |
ppath='' |
31 |
mocchiut |
1.1 |
|
32 |
|
|
la=80 |
33 |
|
|
do i=1,80 |
34 |
|
|
if(cpath(i:i).eq.'/')la=i |
35 |
|
|
enddo |
36 |
|
|
ppath=cpath(1:la) |
37 |
|
|
|
38 |
|
|
b_path = ppath |
39 |
|
|
b_pathlen = la |
40 |
|
|
|
41 |
|
|
c FUNZIONAAAAAAAAAAAAAAA!!!!!!!!!!!!!!!!!! |
42 |
|
|
|
43 |
pam-fi |
1.2 |
c if(b_debug.eq.1)print*,'Field loaded: ',b_loaded |
44 |
|
|
if(DEBUG)print*,'Field loaded: ',b_loaded |
45 |
mocchiut |
1.1 |
if(b_loaded.eq.0)then |
46 |
|
|
|
47 |
|
|
c call the subroutine which reads the maps of the measurements taken |
48 |
|
|
c inside the magnetic cavity |
49 |
|
|
call readBinner(ppath) |
50 |
|
|
if(b_error.eq.1)return |
51 |
|
|
|
52 |
|
|
c call the subroutine which reads the maps of the measurements taken |
53 |
|
|
c outside the magnetic cavity |
54 |
|
|
call readBouter(ppath) |
55 |
|
|
if(b_error.eq.1)return |
56 |
|
|
|
57 |
|
|
b_loaded = 1 |
58 |
|
|
endif |
59 |
|
|
|
60 |
|
|
return |
61 |
|
|
end |
62 |
|
|
|
63 |
|
|
************************************************************ |
64 |
|
|
|
65 |
|
|
|
66 |
|
|
************************************************************************* |
67 |
|
|
* |
68 |
|
|
* Subroutine readBinner |
69 |
|
|
* |
70 |
|
|
* it reads from rz files the two magnetic field maps taken inside the |
71 |
|
|
* spectrometer cavity and fills the variables in common_B_inner.f |
72 |
|
|
* |
73 |
|
|
* needs: |
74 |
|
|
* - common_B.f common file for the inner magnetic field map |
75 |
|
|
* - .rz map files in ./ containing coordinates of measured points, Bx, By |
76 |
|
|
* and Bz components + errors |
77 |
|
|
* |
78 |
|
|
* output variables: (see common_B_inner.f) |
79 |
|
|
* - px#(nx,3) with #=1,2 for the 2 maps |
80 |
|
|
* - py#(ny,3) |
81 |
|
|
* - pz#(nz,3) |
82 |
|
|
* - b#(nx,ny,nz,3) |
83 |
|
|
* |
84 |
|
|
************************************************************************* |
85 |
|
|
|
86 |
|
|
subroutine readBinner(path) |
87 |
|
|
|
88 |
|
|
c implicit double precision (a-h,o-z) |
89 |
|
|
include 'common_B.f' |
90 |
|
|
include 'common_c2f.f' |
91 |
pam-fi |
1.2 |
|
92 |
|
|
LOGICAL DEBUG |
93 |
|
|
LOGICAL VERBOSE |
94 |
|
|
LOGICAL WARNING |
95 |
|
|
COMMON/DBG/DEBUG,VERBOSE,WARNING |
96 |
|
|
SAVE/DBG/ |
97 |
|
|
|
98 |
mocchiut |
1.1 |
character*80 path |
99 |
|
|
|
100 |
|
|
C |
101 |
|
|
REAL hmemor(10000000) |
102 |
|
|
integer Iquest(100) |
103 |
|
|
COMMON /pawc/hmemor |
104 |
|
|
save /pawc/ |
105 |
|
|
C |
106 |
|
|
Common /QUEST/ Iquest |
107 |
|
|
save /quest/ |
108 |
|
|
|
109 |
|
|
|
110 |
|
|
c------------------------------------------------------------------------ |
111 |
|
|
c |
112 |
|
|
c local variables |
113 |
|
|
c |
114 |
|
|
c------------------------------------------------------------------------ |
115 |
|
|
|
116 |
|
|
character*64 Bmap_file !magnetic field file name |
117 |
|
|
parameter (lun_Bmap_file=66) !magnetic field map file id number |
118 |
|
|
|
119 |
|
|
parameter (ntpl_Bmap=20) !ntuple identifier |
120 |
|
|
|
121 |
|
|
REAL PFX(3),FX,DFX, !Bx field component coordinates in m, value and error in T |
122 |
|
|
$ PFY(3),FY,DFY |
123 |
|
|
$ ,PFZ(3),FZ,DFZ |
124 |
|
|
INTEGER INDEX(3) !point index |
125 |
|
|
|
126 |
|
|
COMMON /PAWCR4/ INDEX,PFX,FX,DFX,PFY,FY,DFY,PFZ,FZ,DFZ |
127 |
|
|
|
128 |
|
|
|
129 |
|
|
CALL HLIMIT(10000000) |
130 |
|
|
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
131 |
|
|
C largest RZ file: IQUEST(10) records x LREC words x 4 byte |
132 |
|
|
C with the following settings: 65000 x 4096 x 4 = 1G |
133 |
|
|
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
134 |
|
|
|
135 |
|
|
c permette di ottenere ntuple funzionanti nonostante |
136 |
|
|
c il messaggio dei 64K di RZOUT |
137 |
|
|
Iquest(10) = 512000 |
138 |
|
|
c------------------------------------------------------------------------ |
139 |
|
|
c |
140 |
|
|
c *** FIRST MAP *** |
141 |
|
|
c |
142 |
|
|
c------------------------------------------------------------------------ |
143 |
|
|
|
144 |
|
|
c------------------------------------------------------------------------ |
145 |
|
|
c |
146 |
|
|
c initialization and map file opening |
147 |
|
|
c |
148 |
|
|
c------------------------------------------------------------------------ |
149 |
|
|
|
150 |
|
|
c print*,' ' |
151 |
|
|
c print*,' ' |
152 |
|
|
|
153 |
|
|
Bmap_file='measure_n3_290302.rz' |
154 |
|
|
|
155 |
|
|
c opens magnetic field map first file |
156 |
pam-fi |
1.2 |
c if(b_debug.eq.1)print * |
157 |
|
|
if(VERBOSE)print * |
158 |
mocchiut |
1.1 |
$ ,path(1:LNBLNK(path))//Bmap_file |
159 |
|
|
call HROPEN |
160 |
|
|
$ (lun_Bmap_file,'Bmap' |
161 |
|
|
$ ,path(1:LNBLNK(path))//Bmap_file |
162 |
|
|
$ ,'P',1024,istat) |
163 |
|
|
if(istat.ne.0) goto 21 |
164 |
|
|
|
165 |
|
|
|
166 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
167 |
|
|
|
168 |
|
|
c call HPRNTU(ntpl_Bmap) |
169 |
|
|
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
170 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
171 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
172 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
173 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
174 |
|
|
|
175 |
|
|
|
176 |
|
|
c------------------------------------------------------------------------ |
177 |
|
|
c |
178 |
|
|
c reads events and fills variables |
179 |
|
|
c |
180 |
|
|
c------------------------------------------------------------------------ |
181 |
|
|
|
182 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
183 |
|
|
|
184 |
|
|
c initializes measurement grid edges |
185 |
|
|
do ic=1,3 |
186 |
|
|
px1max(ic)=0. |
187 |
|
|
px1min(ic)=0. |
188 |
|
|
py1max(ic)=0. |
189 |
|
|
py1min(ic)=0. |
190 |
|
|
pz1max(ic)=0. |
191 |
|
|
pz1min(ic)=0. |
192 |
|
|
enddo |
193 |
|
|
|
194 |
|
|
|
195 |
|
|
do iev=1,iemax !event loop |
196 |
|
|
|
197 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
198 |
|
|
if(ierr.ne.0) goto 22 |
199 |
|
|
|
200 |
|
|
c the output consists of matrices for coordinates, B components values |
201 |
|
|
c and errors: |
202 |
|
|
c e.g. px1(4,2) = X coordinate of the point with index = 4 along X, |
203 |
|
|
c in which By (=2) component has been measured |
204 |
|
|
c e.g. b1(3,23,4,1) = Bx (=1) component value, measured in the point with |
205 |
|
|
c indexes = 3,23,4 along X, Y and Z |
206 |
|
|
|
207 |
|
|
c Bx component |
208 |
|
|
px1(index(1),1) = pfx(1) |
209 |
|
|
if(px1(index(1),1).lt.px1min(1)) px1min(1)=px1(index(1),1) |
210 |
|
|
if(px1(index(1),1).gt.px1max(1)) px1max(1)=px1(index(1),1) |
211 |
|
|
py1(index(2),1) = pfx(2) |
212 |
|
|
if(py1(index(2),1).lt.py1min(1)) py1min(1)=py1(index(2),1) |
213 |
|
|
if(py1(index(2),1).gt.py1max(1)) py1max(1)=py1(index(2),1) |
214 |
|
|
pz1(index(3),1) = pfx(3) |
215 |
|
|
if(pz1(index(3),1).lt.pz1min(1)) pz1min(1)=pz1(index(3),1) |
216 |
|
|
if(pz1(index(3),1).gt.pz1max(1)) pz1max(1)=pz1(index(3),1) |
217 |
|
|
|
218 |
|
|
b1(index(1),index(2),index(3),1) = fx |
219 |
|
|
|
220 |
|
|
|
221 |
|
|
c By component |
222 |
|
|
px1(index(1),2) = pfy(1) |
223 |
|
|
if(px1(index(1),2).lt.px1min(2)) px1min(2)=px1(index(1),2) |
224 |
|
|
if(px1(index(1),2).gt.px1max(2)) px1max(2)=px1(index(1),2) |
225 |
|
|
py1(index(2),2) = pfy(2) |
226 |
|
|
if(py1(index(2),2).lt.py1min(2)) py1min(2)=py1(index(2),2) |
227 |
|
|
if(py1(index(2),2).gt.py1max(2)) py1max(2)=py1(index(2),2) |
228 |
|
|
pz1(index(3),2) = pfy(3) |
229 |
|
|
if(pz1(index(3),2).lt.pz1min(2)) pz1min(2)=pz1(index(3),2) |
230 |
|
|
if(pz1(index(3),2).gt.pz1max(2)) pz1max(2)=pz1(index(3),2) |
231 |
|
|
|
232 |
|
|
b1(index(1),index(2),index(3),2) = fy |
233 |
|
|
|
234 |
|
|
|
235 |
|
|
c Bz component |
236 |
|
|
px1(index(1),3) = pfz(1) |
237 |
|
|
if(px1(index(1),3).lt.px1min(3)) px1min(3)=px1(index(1),3) |
238 |
|
|
if(px1(index(1),3).gt.px1max(3)) px1max(3)=px1(index(1),3) |
239 |
|
|
py1(index(2),3) = pfz(2) |
240 |
|
|
if(py1(index(2),3).lt.py1min(3)) py1min(3)=py1(index(2),3) |
241 |
|
|
if(py1(index(2),3).gt.py1max(3)) py1max(3)=py1(index(2),3) |
242 |
|
|
pz1(index(3),3) = pfz(3) |
243 |
|
|
if(pz1(index(3),3).lt.pz1min(3)) pz1min(3)=pz1(index(3),3) |
244 |
|
|
if(pz1(index(3),3).gt.pz1max(3)) pz1max(3)=pz1(index(3),3) |
245 |
|
|
|
246 |
|
|
b1(index(1),index(2),index(3),3) = fz |
247 |
|
|
|
248 |
|
|
enddo |
249 |
|
|
|
250 |
|
|
c------------------------------------------------------------------------ |
251 |
|
|
c |
252 |
|
|
c closes files |
253 |
|
|
c |
254 |
|
|
c------------------------------------------------------------------------ |
255 |
|
|
|
256 |
|
|
call HREND('Bmap') |
257 |
|
|
close(lun_Bmap_file) |
258 |
|
|
c$$$ cmd2='rm -f ' |
259 |
|
|
c$$$ $ //Bmap_file(1:LNBLNK(Bmap_file)) |
260 |
|
|
c$$$ call system(cmd2) |
261 |
|
|
c$$$ |
262 |
|
|
|
263 |
|
|
|
264 |
|
|
c------------------------------------------------------------------------ |
265 |
|
|
c |
266 |
|
|
c *** SECOND MAP *** |
267 |
|
|
c |
268 |
|
|
c------------------------------------------------------------------------ |
269 |
|
|
|
270 |
|
|
c------------------------------------------------------------------------ |
271 |
|
|
c |
272 |
|
|
c initialization and map file opening |
273 |
|
|
c |
274 |
|
|
c------------------------------------------------------------------------ |
275 |
|
|
|
276 |
|
|
c print*,' ' |
277 |
|
|
c print*,' ' |
278 |
|
|
|
279 |
|
|
Bmap_file='measure_n4_110402_corrected.rz' |
280 |
|
|
c opens magnetic field map first file |
281 |
pam-fi |
1.2 |
c if(b_debug.eq.1)print * !,'Opening file: ' |
282 |
|
|
if(VERBOSE)print * !,'Opening file: ' |
283 |
mocchiut |
1.1 |
$ ,path(1:LNBLNK(path))//Bmap_file |
284 |
|
|
call HROPEN |
285 |
|
|
$ (lun_Bmap_file,'Bmap' |
286 |
|
|
$ ,path(1:LNBLNK(path))//Bmap_file |
287 |
|
|
$ ,'P',1024,istat) |
288 |
|
|
if(istat.ne.0) goto 21 |
289 |
|
|
|
290 |
|
|
|
291 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
292 |
|
|
|
293 |
|
|
c call HPRNTU(ntpl_Bmap) |
294 |
|
|
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
295 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
296 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
297 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
298 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
299 |
|
|
|
300 |
|
|
|
301 |
|
|
c------------------------------------------------------------------------ |
302 |
|
|
c |
303 |
|
|
c reads events and fills variables |
304 |
|
|
c |
305 |
|
|
c------------------------------------------------------------------------ |
306 |
|
|
|
307 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
308 |
|
|
|
309 |
|
|
c print*,'iemax ',iemax |
310 |
|
|
|
311 |
|
|
do ic=1,3 !grid edges |
312 |
|
|
px2max(ic)=0. |
313 |
|
|
px2min(ic)=0. |
314 |
|
|
py2max(ic)=0. |
315 |
|
|
py2min(ic)=0. |
316 |
|
|
pz2max(ic)=0. |
317 |
|
|
pz2min(ic)=0. |
318 |
|
|
enddo |
319 |
|
|
|
320 |
|
|
|
321 |
|
|
do iev=1,iemax !event loop |
322 |
|
|
|
323 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
324 |
|
|
if(ierr.ne.0) goto 22 |
325 |
|
|
|
326 |
|
|
c the output consists of matrices for coordinates, B components values |
327 |
|
|
c and errors: |
328 |
|
|
c e.g. px(4,2) = X coordinate of the point with index = 4 along X, |
329 |
|
|
c in which By (=2) component has been measured |
330 |
|
|
c e.g. b(3,23,4,1) = Bx (=1) component value, measured in the point with |
331 |
|
|
c indexes = 3,23,4 along X, Y and Z |
332 |
|
|
|
333 |
|
|
c Bx component |
334 |
|
|
px2(index(1),1) = pfx(1) |
335 |
|
|
if(px2(index(1),1).lt.px2min(1)) px2min(1)=px2(index(1),1) |
336 |
|
|
if(px2(index(1),1).gt.px2max(1)) px2max(1)=px2(index(1),1) |
337 |
|
|
py2(index(2),1) = pfx(2) |
338 |
|
|
if(py2(index(2),1).lt.py2min(1)) py2min(1)=py2(index(2),1) |
339 |
|
|
if(py2(index(2),1).gt.py2max(1)) py2max(1)=py2(index(2),1) |
340 |
|
|
pz2(index(3),1) = pfx(3) |
341 |
|
|
if(pz2(index(3),1).lt.pz2min(1)) pz2min(1)=pz2(index(3),1) |
342 |
|
|
if(pz2(index(3),1).gt.pz2max(1)) pz2max(1)=pz2(index(3),1) |
343 |
|
|
|
344 |
|
|
b2(index(1),index(2),index(3),1) = fx |
345 |
|
|
|
346 |
|
|
|
347 |
|
|
c By component |
348 |
|
|
px2(index(1),2) = pfy(1) |
349 |
|
|
if(px2(index(1),2).lt.px2min(2)) px2min(2)=px2(index(1),2) |
350 |
|
|
if(px2(index(1),2).gt.px2max(2)) px2max(2)=px2(index(1),2) |
351 |
|
|
py2(index(2),2) = pfy(2) |
352 |
|
|
if(py2(index(2),2).lt.py2min(2)) py2min(2)=py2(index(2),2) |
353 |
|
|
if(py2(index(2),2).gt.py2max(2)) py2max(2)=py2(index(2),2) |
354 |
|
|
pz2(index(3),2) = pfy(3) |
355 |
|
|
if(pz2(index(3),2).lt.pz2min(2)) pz2min(2)=pz2(index(3),2) |
356 |
|
|
if(pz2(index(3),2).gt.pz2max(2)) pz2max(2)=pz2(index(3),2) |
357 |
|
|
|
358 |
|
|
b2(index(1),index(2),index(3),2) = fy |
359 |
|
|
|
360 |
|
|
|
361 |
|
|
c Bz component |
362 |
|
|
px2(index(1),3) = pfz(1) |
363 |
|
|
if(px2(index(1),3).lt.px2min(3)) px2min(3)=px2(index(1),3) |
364 |
|
|
if(px2(index(1),3).gt.px2max(3)) px2max(3)=px2(index(1),3) |
365 |
|
|
py2(index(2),3) = pfz(2) |
366 |
|
|
if(py2(index(2),3).lt.py2min(3)) py2min(3)=py2(index(2),3) |
367 |
|
|
if(py2(index(2),3).gt.py2max(3)) py2max(3)=py2(index(2),3) |
368 |
|
|
pz2(index(3),3) = pfz(3) |
369 |
|
|
if(pz2(index(3),3).lt.pz2min(3)) pz2min(3)=pz2(index(3),3) |
370 |
|
|
if(pz2(index(3),3).gt.pz2max(3)) pz2max(3)=pz2(index(3),3) |
371 |
|
|
|
372 |
|
|
b2(index(1),index(2),index(3),3) = fz |
373 |
|
|
|
374 |
|
|
enddo |
375 |
|
|
|
376 |
|
|
c------------------------------------------------------------------------ |
377 |
|
|
c |
378 |
|
|
c closes files |
379 |
|
|
c |
380 |
|
|
c------------------------------------------------------------------------ |
381 |
|
|
|
382 |
|
|
call HREND('Bmap') |
383 |
|
|
close(lun_Bmap_file) |
384 |
|
|
|
385 |
|
|
c------------------------------------------------------------------------ |
386 |
|
|
c |
387 |
|
|
c no error exit |
388 |
|
|
c |
389 |
|
|
c------------------------------------------------------------------------ |
390 |
|
|
|
391 |
|
|
goto 9000 !happy ending |
392 |
|
|
|
393 |
|
|
c------------------------------------------------------------------------ |
394 |
|
|
c |
395 |
|
|
c magnetic field map file opening error |
396 |
|
|
c |
397 |
|
|
c------------------------------------------------------------------------ |
398 |
|
|
|
399 |
|
|
21 continue |
400 |
|
|
|
401 |
|
|
b_error = 1 |
402 |
pam-fi |
1.2 |
c if(b_debug.eq.1) |
403 |
|
|
if(DEBUG) |
404 |
mocchiut |
1.1 |
$ print* |
405 |
|
|
$ ,'read_B_inner: ERROR OPENING MAGNETIC FIELD MAP FILE: ' |
406 |
|
|
$ ,Bmap_file |
407 |
|
|
|
408 |
|
|
goto 9000 !the end |
409 |
|
|
|
410 |
|
|
|
411 |
|
|
c------------------------------------------------------------------------ |
412 |
|
|
c |
413 |
|
|
c ntuple event reading error |
414 |
|
|
c |
415 |
|
|
c------------------------------------------------------------------------ |
416 |
|
|
|
417 |
|
|
22 continue |
418 |
|
|
|
419 |
|
|
b_error = 1 |
420 |
pam-fi |
1.2 |
c if(b_debug.eq.1) |
421 |
|
|
if(DEBUG) |
422 |
mocchiut |
1.1 |
$ print*,'read_B_inner: ERROR WHILE READING NTUPLE, at entry |
423 |
|
|
$ : ',iev |
424 |
|
|
|
425 |
|
|
goto 9000 !the end |
426 |
|
|
|
427 |
|
|
|
428 |
|
|
c------------------------------------------------------------------------ |
429 |
|
|
c |
430 |
|
|
c exit |
431 |
|
|
c |
432 |
|
|
c------------------------------------------------------------------------ |
433 |
|
|
|
434 |
|
|
9000 continue |
435 |
|
|
|
436 |
|
|
return |
437 |
|
|
end |
438 |
|
|
************************************************************************* |
439 |
|
|
* |
440 |
|
|
* Subroutine readBouter |
441 |
|
|
* |
442 |
|
|
* it reads from rz files the two magnetic field maps taken inside the |
443 |
|
|
* spectrometer cavity and fills the variables in common_B_inner.f |
444 |
|
|
* |
445 |
|
|
* needs: |
446 |
|
|
* - common_B_outer.f common file for the outer magnetic field map |
447 |
|
|
* - .rz map files in ./ containing coordinates of measured points, Bx, By |
448 |
|
|
* and Bz components + errors |
449 |
|
|
* |
450 |
|
|
* output variables: (see common_B_outer.f) |
451 |
|
|
* - pxo(nx,3) |
452 |
|
|
* - pyo(ny,3) |
453 |
|
|
* - pzo(nz,3) |
454 |
|
|
* - bo(nx,ny,nz,3) |
455 |
|
|
* |
456 |
|
|
************************************************************************* |
457 |
|
|
|
458 |
|
|
subroutine readBouter(path) |
459 |
|
|
|
460 |
|
|
c implicit double precision (a-h,o-z) |
461 |
|
|
include 'common_B.f' |
462 |
|
|
include 'common_c2f.f' |
463 |
|
|
C |
464 |
pam-fi |
1.2 |
LOGICAL DEBUG |
465 |
|
|
LOGICAL VERBOSE |
466 |
|
|
LOGICAL WARNING |
467 |
|
|
COMMON/DBG/DEBUG,VERBOSE,WARNING |
468 |
|
|
SAVE/DBG/ |
469 |
|
|
|
470 |
mocchiut |
1.1 |
character*80 path |
471 |
|
|
|
472 |
|
|
REAL hmemor(10000000) |
473 |
|
|
integer Iquest(100) |
474 |
|
|
COMMON /pawc/hmemor |
475 |
|
|
save /pawc/ |
476 |
|
|
C |
477 |
|
|
Common /QUEST/ Iquest |
478 |
|
|
save /quest/ |
479 |
|
|
|
480 |
|
|
|
481 |
|
|
|
482 |
|
|
c------------------------------------------------------------------------ |
483 |
|
|
c |
484 |
|
|
c local variables |
485 |
|
|
c |
486 |
|
|
c------------------------------------------------------------------------ |
487 |
|
|
|
488 |
|
|
character*64 Bmap_file !magnetic field file name |
489 |
|
|
parameter (lun_Bmap_file=66) !magnetic field map file id number |
490 |
|
|
|
491 |
|
|
parameter (ntpl_Bmap=20) !ntuple identifier |
492 |
|
|
|
493 |
|
|
REAL PFX(3),FX,DFX, !Bx field component coordinates in m, value and error in T |
494 |
|
|
$ PFY(3),FY,DFY |
495 |
|
|
$ ,PFZ(3),FZ,DFZ |
496 |
|
|
INTEGER INDEX(3) !point index |
497 |
|
|
|
498 |
|
|
COMMON /PAWCR4/ INDEX,PFX,FX,DFX,PFY,FY,DFY,PFZ,FZ,DFZ |
499 |
|
|
|
500 |
|
|
|
501 |
|
|
c print*,'Calling HLIMIT' |
502 |
|
|
CALL HLIMIT(10000000) |
503 |
|
|
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
504 |
|
|
C largest RZ file: IQUEST(10) records x LREC words x 4 byte |
505 |
|
|
C with the following settings: 65000 x 4096 x 4 = 1G |
506 |
|
|
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
507 |
|
|
|
508 |
|
|
c permette di ottenere ntuple funzionanti nonostante |
509 |
|
|
c il messaggio dei 64K di RZOUT |
510 |
|
|
Iquest(10) = 512000 |
511 |
|
|
c------------------------------------------------------------------------ |
512 |
|
|
c |
513 |
|
|
c *** FIRST MAP *** |
514 |
|
|
c |
515 |
|
|
c------------------------------------------------------------------------ |
516 |
|
|
|
517 |
|
|
c------------------------------------------------------------------------ |
518 |
|
|
c |
519 |
|
|
c initialization and map file opening |
520 |
|
|
c |
521 |
|
|
c------------------------------------------------------------------------ |
522 |
|
|
|
523 |
|
|
c print*,' ' |
524 |
|
|
c print*,' ' |
525 |
|
|
|
526 |
|
|
Bmap_file='External_top_map_n4_150402.rz' |
527 |
|
|
|
528 |
|
|
c opens magnetic field map first file |
529 |
pam-fi |
1.2 |
c if(b_debug.eq.1)print * |
530 |
|
|
if(VERBOSE)print * |
531 |
mocchiut |
1.1 |
$ ,path(1:LNBLNK(path))//Bmap_file |
532 |
|
|
call HROPEN |
533 |
|
|
$ (lun_Bmap_file,'Bmap' |
534 |
|
|
$ ,path(1:LNBLNK(path))//Bmap_file |
535 |
|
|
$ ,'P',1024,istat) |
536 |
|
|
if(istat.ne.0) goto 21 |
537 |
|
|
|
538 |
|
|
|
539 |
|
|
call HRIN(ntpl_Bmap,9999,0) !puts B map ntuple in memory |
540 |
|
|
|
541 |
|
|
c call HPRNTU(ntpl_Bmap) |
542 |
|
|
call HBNAME(ntpl_Bmap,' ',0,'$CLEAR') |
543 |
|
|
call HBNAME(ntpl_Bmap,'INDEX',index,'$SET') |
544 |
|
|
call HBNAME(ntpl_Bmap,'BX',pfx,'$SET') |
545 |
|
|
call HBNAME(ntpl_Bmap,'BY',pfy,'$SET') |
546 |
|
|
call HBNAME(ntpl_Bmap,'BZ',pfz,'$SET') |
547 |
|
|
|
548 |
|
|
|
549 |
|
|
c------------------------------------------------------------------------ |
550 |
|
|
c |
551 |
|
|
c reads events and fills variables |
552 |
|
|
c |
553 |
|
|
c------------------------------------------------------------------------ |
554 |
|
|
|
555 |
|
|
call HNOENT(ntpl_Bmap,iemax) !number of events |
556 |
|
|
|
557 |
|
|
c initializes measurement grid edges |
558 |
|
|
do ic=1,3 |
559 |
|
|
poxmax(ic)=0. |
560 |
|
|
poxmin(ic)=0. |
561 |
|
|
poymax(ic)=0. |
562 |
|
|
poymin(ic)=0. |
563 |
|
|
pozmax(ic)=0. |
564 |
|
|
pozmin(ic)=0. |
565 |
|
|
enddo |
566 |
|
|
|
567 |
|
|
|
568 |
|
|
do iev=1,iemax !event loop |
569 |
|
|
|
570 |
|
|
call HGNT(ntpl_Bmap,iev,ierr) !reads event |
571 |
|
|
if(ierr.ne.0) goto 22 |
572 |
|
|
|
573 |
|
|
c the output consists of matrices for coordinates, B components values |
574 |
|
|
c and errors: |
575 |
|
|
c e.g. px1(4,2) = X coordinate of the point with index = 4 along X, |
576 |
|
|
c in which By (=2) component has been measured |
577 |
|
|
c e.g. b1(3,23,4,1) = Bx (=1) component value, measured in the point with |
578 |
|
|
c indexes = 3,23,4 along X, Y and Z |
579 |
|
|
|
580 |
|
|
c Bx component |
581 |
|
|
pox(index(1),1) = pfx(1) |
582 |
|
|
if(pox(index(1),1).lt.poxmin(1)) poxmin(1)=pox(index(1),1) |
583 |
|
|
if(pox(index(1),1).gt.poxmax(1)) poxmax(1)=pox(index(1),1) |
584 |
|
|
poy(index(2),1) = pfx(2) |
585 |
|
|
if(poy(index(2),1).lt.poymin(1)) poymin(1)=poy(index(2),1) |
586 |
|
|
if(poy(index(2),1).gt.poymax(1)) poymax(1)=poy(index(2),1) |
587 |
|
|
poz(index(3),1) = pfx(3) |
588 |
|
|
if(poz(index(3),1).lt.pozmin(1)) pozmin(1)=poz(index(3),1) |
589 |
|
|
if(poz(index(3),1).gt.pozmax(1)) pozmax(1)=poz(index(3),1) |
590 |
|
|
|
591 |
|
|
bo(index(1),index(2),index(3),1) = fx |
592 |
|
|
|
593 |
|
|
|
594 |
|
|
c By component |
595 |
|
|
pox(index(1),2) = pfy(1) |
596 |
|
|
if(pox(index(1),2).lt.poxmin(2)) poxmin(2)=pox(index(1),2) |
597 |
|
|
if(pox(index(1),2).gt.poxmax(2)) poxmax(2)=pox(index(1),2) |
598 |
|
|
poy(index(2),2) = pfy(2) |
599 |
|
|
if(poy(index(2),2).lt.poymin(2)) poymin(2)=poy(index(2),2) |
600 |
|
|
if(poy(index(2),2).gt.poymax(2)) poymax(2)=poy(index(2),2) |
601 |
|
|
poz(index(3),2) = pfy(3) |
602 |
|
|
if(poz(index(3),2).lt.pozmin(2)) pozmin(2)=poz(index(3),2) |
603 |
|
|
if(poz(index(3),2).gt.pozmax(2)) pozmax(2)=poz(index(3),2) |
604 |
|
|
|
605 |
|
|
bo(index(1),index(2),index(3),2) = fy |
606 |
|
|
|
607 |
|
|
|
608 |
|
|
c Bz component |
609 |
|
|
pox(index(1),3) = pfz(1) |
610 |
|
|
if(pox(index(1),3).lt.poxmin(3)) poxmin(3)=pox(index(1),3) |
611 |
|
|
if(pox(index(1),3).gt.poxmax(3)) poxmax(3)=pox(index(1),3) |
612 |
|
|
poy(index(2),3) = pfz(2) |
613 |
|
|
if(poy(index(2),3).lt.poymin(3)) poymin(3)=poy(index(2),3) |
614 |
|
|
if(poy(index(2),3).gt.poymax(3)) poymax(3)=poy(index(2),3) |
615 |
|
|
poz(index(3),3) = pfz(3) |
616 |
|
|
if(poz(index(3),3).lt.pozmin(3)) pozmin(3)=poz(index(3),3) |
617 |
|
|
if(poz(index(3),3).gt.pozmax(3)) pozmax(3)=poz(index(3),3) |
618 |
|
|
|
619 |
|
|
bo(index(1),index(2),index(3),3) = fz |
620 |
|
|
|
621 |
|
|
enddo |
622 |
|
|
|
623 |
|
|
|
624 |
|
|
c------------------------------------------------------------------------ |
625 |
|
|
c |
626 |
|
|
c closes files |
627 |
|
|
c |
628 |
|
|
c------------------------------------------------------------------------ |
629 |
|
|
|
630 |
|
|
call HREND('Bmap') |
631 |
|
|
close(lun_Bmap_file) |
632 |
|
|
|
633 |
|
|
|
634 |
|
|
c------------------------------------------------------------------------ |
635 |
|
|
c |
636 |
|
|
c no error exit |
637 |
|
|
c |
638 |
|
|
c------------------------------------------------------------------------ |
639 |
|
|
|
640 |
|
|
goto 9000 !happy ending |
641 |
|
|
|
642 |
|
|
c------------------------------------------------------------------------ |
643 |
|
|
c |
644 |
|
|
c magnetic field map file opening error |
645 |
|
|
c |
646 |
|
|
c------------------------------------------------------------------------ |
647 |
|
|
|
648 |
|
|
21 continue |
649 |
|
|
|
650 |
|
|
b_error = 1 |
651 |
pam-fi |
1.2 |
c if(b_debug.eq.1) |
652 |
|
|
if(DEBUG) |
653 |
mocchiut |
1.1 |
$ print* |
654 |
|
|
$ ,'read_B_inner: ERROR OPENING MAGNETIC FIELD MAP FILE: ' |
655 |
|
|
$ ,Bmap_file |
656 |
|
|
|
657 |
|
|
goto 9000 !the end |
658 |
|
|
|
659 |
|
|
|
660 |
|
|
c------------------------------------------------------------------------ |
661 |
|
|
c |
662 |
|
|
c ntuple event reading error |
663 |
|
|
c |
664 |
|
|
c------------------------------------------------------------------------ |
665 |
|
|
|
666 |
|
|
22 continue |
667 |
|
|
|
668 |
|
|
b_error = 1 |
669 |
pam-fi |
1.2 |
c if(b_debug.eq.1) |
670 |
|
|
if(DEBUG) |
671 |
mocchiut |
1.1 |
$ print* |
672 |
|
|
$ ,'read_B_inner: ERROR WHILE READING NTUPLE, at event |
673 |
|
|
$ : ',iev |
674 |
|
|
|
675 |
|
|
goto 9000 !the end |
676 |
|
|
|
677 |
|
|
|
678 |
|
|
c------------------------------------------------------------------------ |
679 |
|
|
c |
680 |
|
|
c exit |
681 |
|
|
c |
682 |
|
|
c------------------------------------------------------------------------ |
683 |
|
|
|
684 |
|
|
9000 continue |
685 |
|
|
|
686 |
|
|
return |
687 |
|
|
end |
688 |
|
|
|