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