/[PAMELA software]/gpamela/gpfield/gprbi.F
ViewVC logotype

Annotation of /gpamela/gpfield/gprbi.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3.2 - (hide annotations) (download)
Thu Dec 22 16:26:02 2005 UTC (18 years, 11 months ago) by bottai
Branch: MAIN
CVS Tags: v4r4, v4r5, v4r6, v4r7, v4r8, v4r9, v4r14, v4r12, v4r13, v4r10, v4r11, HEAD
Changes since 3.1: +126 -5 lines
 read the magn. field map for external field

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

  ViewVC Help
Powered by ViewVC 1.1.23