Millepede-II  V04-00-00_preview
 All Classes Files Functions Variables Enumerator Pages
mptest1.f90
Go to the documentation of this file.
1 
2 ! Code converted using TO_F90 by Alan Miller
3 ! Date: 2012-03-16 Time: 11:08:48
4 
13 
14 
16 MODULE mptest1
17  IMPLICIT NONE
18  SAVE
19 
20  INTEGER, PARAMETER :: nplan=100
21 
22  ! define detector geometry
23  REAL, PARAMETER :: detx= 10.0 !< x-value of first plane
24  REAL, PARAMETER :: disx= 10.0 !< distance between planes
25  REAL, PARAMETER :: thck= 2.0 !< thickness of plane
26  REAL, PARAMETER :: heit=100.0 !< height of detector plane
27  REAL, PARAMETER :: effp=0.90 !< plane efficiency
28  REAL, PARAMETER :: sgmp=0.0150 !< measurement sigma
29 
30  ! misalignment
31  REAL, DIMENSION(nplan) :: del !< shift (position deviation) (alignment parameter)
32  REAL, DIMENSION(nplan) :: dvd !< rel. drift velocity deviation (calibration parameter)
33  ! track parameter
34  REAL :: ynull !< track position at vertex
35  REAL :: slope !< track slope
36 
37  INTEGER :: nhits !< number of hits
38  INTEGER, DIMENSION(nplan) :: ihits !< plane numbers (planes with hits)
39  REAL, DIMENSION(nplan) :: eff !< plane efficiency
40  REAL, DIMENSION(nplan) :: sgm !< measurement sigma (plane)
41  REAL, DIMENSION(nplan) :: ydrft !< signed drift length
42  REAL, DIMENSION(nplan) :: xhits !< position perp. to plane (hit)
43  REAL, DIMENSION(nplan) :: yhits !< measured position in plane (hit)
44  REAL, DIMENSION(nplan) :: sigma !< measurement sigma (hit)
45 
46 END MODULE mptest1
47 
56 
57 SUBROUTINE mptest
58  USE mptest1
59 
60  IMPLICIT NONE
61  REAL :: dbar
62  REAL :: det
63  REAL :: displ
64  REAL :: drift
65  REAL :: eps
66  REAL :: eta
67  REAL :: gran
68  REAL :: one
69  REAL :: ww
70  REAL :: x
71  REAL :: xbar
72  INTEGER :: i
73  INTEGER :: icount
74  INTEGER :: ios
75  INTEGER :: ip
76  INTEGER :: ipl
77  INTEGER :: labelt
78  INTEGER :: luns
79  INTEGER :: lunt
80  INTEGER :: ncount
81  INTEGER :: nrecds
82  INTEGER :: nthits
83 
84  DOUBLE PRECISION :: s1
85  DOUBLE PRECISION :: s2
86  DOUBLE PRECISION :: sw
87  DOUBLE PRECISION :: sv
88  DOUBLE PRECISION :: sum1
89  DOUBLE PRECISION :: sum2
90  REAL :: derlc(2)
91  REAL :: dergl(2)
92  INTEGER :: label(2)
93  LOGICAL :: ex1
94  LOGICAL :: ex2
95  LOGICAL :: ex3
96  ! ...
97  !CC CALL RNTIME
98  INQUIRE(file='mp2str.txt',iostat=ios,exist=ex1) ! keep, if existing
99  INQUIRE(file='mp2con.txt',iostat=ios,exist=ex2) ! keep, if existing
100 
101  INQUIRE(file='mp2tst.bin',iostat=ios,exist=ex3) ! remove, if existing
102 
103  WRITE(*,*) ' '
104  WRITE(*,*) 'Generating test data for mp II...'
105  WRITE(*,*) ' '
106  ! file management
107  IF(ex3) CALL system('rm mp2tst.bin') ! remove old file
108 
109  IF(.NOT.ex1) OPEN(unit=7,access='SEQUENTIAL',form='FORMATTED', &
110  file='mp2str.txt')
111  IF(.NOT.ex2) OPEN(unit=9,access='SEQUENTIAL',form='FORMATTED', &
112  file='mp2con.txt')
113  OPEN(unit=51,access='SEQUENTIAL',form='UNFORMATTED', file='mp2tst.bin')
114 
115  DO i=1,nplan
116  eff(i)=effp ! plane efficiency
117  sgm(i)=sgmp ! measurement sigma
118  del(i)=0.0 ! true shift is zero
119  END DO
120 
121  ipl=7 ! modify one plane (7)
122  eff(ipl)=0.1 ! low efficiency
123  sgm(ipl)=0.0400 ! bad resolution
124 
125  ! misalign detector planes -----------------------------------------
126 
127  displ=0.1 ! displacement 1 mm * N(0,1)
128  drift=0.02 ! Vdrift deviation 2 % * N(0,1)
129  DO i=1,nplan
130  del(i)=displ*gran() ! shift
131  dvd(i)=drift*gran() ! rel. drift velocitu deviation
132  END DO
133  del(10)=0.0 ! no shift
134  del(90)=0.0 ! no shift
135 
136  ! write text files -------------------------------------------------
137 
138  IF(.NOT.ex1) THEN
139  luns=7 ! steerfile
140  WRITE(luns,101) '* Default test steering file'
141  WRITE(luns,101) 'fortranfiles ! following bin files are fortran'
142  WRITE(luns,101) 'mp2con.txt ! constraints text file '
143  WRITE(luns,101) 'mp2tst.bin ! binary data file'
144  WRITE(luns,101) 'Cfiles ! following bin files are Cfiles'
145  ! WRITE(LUNS,101) '*outlierrejection 100.0 ! reject if Chi^2/Ndf >'
146  ! WRITE(LUNS,101) '*outliersuppression 3 ! 3 local_fit iterations'
147 
148  WRITE(luns,101) '*hugecut 50.0 !cut factor in iteration 0'
149  WRITE(luns,101) '*chisqcut 1.0 1.0 ! cut factor in iterations 1 and 2'
150  WRITE(luns,101) '*entries 10 ! lower limit on number of entries/parameter'
151  WRITE(luns,101) &
152  '*pairentries 10 ! lower limit on number of parameter pairs', &
153  ' ! (not yet!)'
154  WRITE(luns,101) '*printrecord 1 2 ! debug printout for records'
155  WRITE(luns,101) &
156  '*printrecord -1 -1 ! debug printout for bad data records'
157  WRITE(luns,101) &
158  '*outlierdownweighting 2 ! number of internal iterations (> 1)'
159  WRITE(luns,101) '*dwfractioncut 0.2 ! 0 < value < 0.5'
160  WRITE(luns,101) '*presigma 0.01 ! default value for presigma'
161  WRITE(luns,101) '*regularisation 1.0 ! regularisation factor'
162  WRITE(luns,101) '*regularisation 1.0 0.01 ! regularisation factor, pre-sigma'
163 
164  WRITE(luns,101) ' '
165  WRITE(luns,101) '*bandwidth 0 ! width of precond. band matrix'
166  WRITE(luns,101) 'method diagonalization 3 0.001 ! diagonalization '
167  WRITE(luns,101) 'method fullMINRES 3 0.01 ! minimal residual '
168  WRITE(luns,101) 'method sparseMINRES 3 0.01 ! minimal residual '
169  WRITE(luns,101) '*mrestol 1.0D-8 ! epsilon for MINRES'
170  WRITE(luns,101) 'method inversion 3 0.001 ! Gauss matrix inversion'
171  WRITE(luns,101) '* last method is applied'
172  WRITE(luns,101) '*matiter 3 ! recalculate matrix in iterations'
173  WRITE(luns,101) ' '
174  WRITE(luns,101) 'end ! optional for end-of-data'
175  ENDIF
176 
177  lunt=9 ! constraint file
178  one=1.0 ! shift constraint
179  IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0'
180  DO i=1,nplan
181  labelt=10+i*2
182  x=detx+float(i-1)*disx+0.5*thck
183  IF(.NOT.ex2) WRITE(lunt,103) labelt,one
184  END DO
185 
186  sw=0.0d0 ! tilt constraint
187  sv=0.0d0
188  s1=0.0d0
189  s2=0.0d0
190  IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0' ! write
191  dbar=0.5*float(nplan-1)*disx
192  xbar=detx+0.5*float(nplan-1)*disx! +0.5*THCK
193  DO i=1,nplan
194  labelt=10+i*2
195  x=detx+float(i-1)*disx !+0.5*THCK
196  ww=(x-xbar)/dbar
197  IF(.NOT.ex2) WRITE(lunt,103) labelt,ww ! write
198  s1=s1+del(i)
199  s2=s2+ww*del(i)
200  sw=sw+ww
201  sv=sv+ww*ww
202  END DO
203 
204 
205  det=REAL(float(nplan)*sv-sw*sw)
206  eps=REAL(sv*s1-sw*s2)/det
207  eta=REAL(float(nplan)*s2-sw*s1)/det
208  DO i=1,nplan
209  x=detx+float(i-1)*disx
210  ww=(x-xbar)/dbar
211  del(i)=del(i)-eps-eta*ww ! correct displacement ...
212  END DO ! ... for constraints
213 
214  sum1=0.0
215  sum2=0.0
216  DO i=1,nplan
217  sum1=sum1+del(i)
218  x=detx+float(i-1)*disx !+0.5*THCK
219  ww=(x-xbar)/dbar
220  sum2=sum2+del(i)*ww
221  END DO
222  ! WRITE(*,*) ' Check for constraints ',SUM1,SUM2
223 
224  ! record loop ------------------------------------------------------
225 
226  ncount=10000
227  nthits=0
228  nrecds=0
229 
230  DO icount=1,ncount
231  ip=0
232  IF(icount == 8759) ip=1
233  ! IF(ICOUNT.EQ.6309) IP=1
234  ! IF(ICOUNT.EQ.7468) IP=1
235  CALL genlin(ip) ! generate hits
236 
237  DO i=1,nhits
238  derlc(1)=1.0
239  derlc(2)=xhits(i)
240  dergl(1)=1.0
241  dergl(2)=ydrft(i)
242  label(1)=10+ihits(i)*2
243  label(2)=500 + ihits(i)
244  CALL mille(2,derlc,2,dergl,label,yhits(i),sigma(i))
245  nthits=nthits+1 ! count hits
246  END DO
247  CALL endle
248  nrecds=nrecds+1 ! count records
249  END DO
250 
251  ! ------------------------------------------------------------------
252  IF(.NOT.ex1) THEN
253  rewind(7)
254  CLOSE (7)
255  END IF
256  IF(.NOT.ex2) THEN
257  rewind(9)
258  CLOSE (9)
259  END IF
260  rewind(51)
261  CLOSE (51)
262 
263  ! WRITE(*,*) ' '
264  ! WRITE(*,*) 'Shifts and drift velocity deviations:'
265  ! DO I=1,NPLAN
266  ! WRITE(*,102) I,DEL(I),DVD(I)
267  ! END DO
268 
269 
270  WRITE(*,*) ' '
271  WRITE(*,*) ' '
272  WRITE(*,*) ncount,' tracks generated with ',nthits,' hits.'
273  WRITE(*,*) nrecds,' records written.'
274  WRITE(*,*) ' '
275 101 format(a)
276  ! 102 FORMAT(I6,2F10.5)
277 103 format(i8,f10.5)
278 END SUBROUTINE mptest ! gener
279 
283 
284 SUBROUTINE genlin(ip)
285  USE mptest1
286 
287  IMPLICIT NONE
288  REAL :: gr
289  REAL :: gran
290  REAL :: uran
291  REAL :: x
292  REAL :: ybias
293  REAL :: ydvds
294  REAL :: ylin
295  REAL :: ymeas
296  REAL :: ywire
297  INTEGER :: i
298  INTEGER :: nwire
299 
300  INTEGER, INTENT(IN) :: ip
301 
302  ! ...
303  ynull=0.5*heit+0.1*heit*(uran()-0.5) ! uniform vertex
304  slope=(uran()-0.5)*heit/(float(nplan-1)*disx)
305  IF(ip /= 0) THEN
306  WRITE(*,*) ' '
307  ! WRITE(*,*) 'YNULL=',YNULL,' SLOPE=',SLOPE
308  END IF
309  nhits=0
310  DO i=1,nplan
311  x=detx+float(i-1)*disx ! +0.5*THCK
312  IF(uran() < eff(i)) THEN
313  ylin =ynull+slope*x ! true y value
314  ybias =ylin-del(i) ! biased value
315  nwire=int(1.0+ybias/4.0) ! wire number
316  IF(nwire <= 0.OR.nwire > 25) exit ! check wire number
317  nhits=nhits+1 ! track hits the plane
318  xhits(nhits)=x
319  ihits(nhits)=i
320  gr=gran()
321  ymeas=sgm(i)*gr
322  ydvds=0.0
323  yhits(nhits)=ybias+ymeas+ydvds ! measured
324  ywire=float(nwire)*4.0-2.0
325  ydrft(nhits)=ybias-ywire ! signed drift length
326  ydvds=ydrft(nhits)*dvd(i)
327  yhits(nhits)=ybias+ymeas-ydvds ! measured
328  sigma(nhits)=sgm(i)
329  IF(ip /= 0) THEN
330  ! WRITE(*,101) NHITS,I,X,YLIN,YBIAS,YMEAS,
331  ! + SGM(I),YHITS(NHITS),GR,DEL(I)
332  END IF
333  END IF
334  END DO
335 ! 101 FORMAT(2I3,F5.0,7F8.4)
336 END SUBROUTINE genlin
Generated on Wed May 23 2012 16:53:05 for Millepede-II by   doxygen 1.8.1