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

Annotation of /gpamela/gpfield/gprbi.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3.1 - (hide annotations) (download)
Tue Dec 6 01:07:35 2005 UTC (18 years, 11 months ago) by cafagna
Branch: MAIN
Adding new magnetic field routines

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

  ViewVC Help
Powered by ViewVC 1.1.23