Millepede-II V04-17-03
mptest2.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:55
4
55
57MODULE mptest2
58 USE mpdef
59
60 IMPLICIT NONE
61 SAVE
62
63 INTEGER(mpi), PARAMETER :: nlyr=10
64 INTEGER(mpi), PARAMETER :: nmlyr=14
65 INTEGER(mpi), PARAMETER :: nmx=10
66 INTEGER(mpi), PARAMETER :: ncx=5
67 INTEGER(mpi), PARAMETER :: nmy=5
68 INTEGER(mpi), PARAMETER :: ntot=nlyr*nmx*nmy
69 ! define detector geometry
70 REAL(mps), PARAMETER :: dets= 10.0
71 REAL(mps), PARAMETER :: diss= 10.0
72 REAL(mps), PARAMETER :: thck= 0.02
73 REAL(mps), PARAMETER :: offs= 0.5
74 REAL(mps), PARAMETER :: stereo=0.08727
75 REAL(mps), PARAMETER :: sizel= 20.0
76 REAL(mps), PARAMETER :: sigl =0.002 ! <resolution
77
78 INTEGER(mpi) :: nhits
79 REAL(mps) :: the0
80 INTEGER(mpi), DIMENSION(nmlyr) :: islyr
81 INTEGER(mpi), DIMENSION(nmlyr) :: ihits
82 REAL(mps), DIMENSION(ntot) :: sdevx
83 REAL(mps), DIMENSION(ntot) :: sdevy
84 REAL(mps), DIMENSION(nmlyr) :: sarc
85 REAL(mps), DIMENSION(nmlyr) :: ssig
86 REAL(mps), DIMENSION(2,nmlyr) :: spro
87 REAL(mps), DIMENSION(nmlyr) :: xhits
88 REAL(mps), DIMENSION(nmlyr) :: yhits
89 REAL(mps), DIMENSION(nmlyr) :: sigma
90
91END MODULE mptest2
92
110
111SUBROUTINE mptst2(imodel) ! generate test files
112 USE mptest2
113 IMPLICIT NONE
114 REAL(mps) :: cmbbrl
115 REAL(mps) :: dispxm
116 REAL(mps) :: dispym
117 REAL(mps) :: dn
118 REAL(mps) :: dp
119 REAL(mps) :: gran
120 REAL(mps) :: one
121 REAL(mps) :: p
122 REAL(mps) :: s
123 REAL(mps) :: sgn
124 REAL(mps) :: sbrl
125 REAL(mps) :: sold
126 REAL(mps) :: uran
127 REAL(mps) :: wbrl
128 INTEGER(mpi) :: i
129 INTEGER(mpi) :: ibrl
130 INTEGER(mpi) :: icount
131 INTEGER(mpi) :: im
132 INTEGER(mpi) :: ios
133 INTEGER(mpi) :: ip
134 INTEGER(mpi) :: j
135 INTEGER(mpi) :: k
136 INTEGER(mpi) :: l
137 INTEGER(mpi) :: labelt
138 INTEGER(mpi) :: layer
139 INTEGER(mpi) :: lb
140 INTEGER(mpi) :: lbrl
141 INTEGER(mpi) :: luns
142 INTEGER(mpi) :: lunt
143 INTEGER(mpi) :: lyr
144 INTEGER(mpi) :: nalc
145 INTEGER(mpi) :: nbrl
146 INTEGER(mpi) :: ncount
147 INTEGER(mpi) :: nmxy
148 INTEGER(mpi) :: nrecds
149 INTEGER(mpi) :: nthits
150
151 INTEGER(mpi), INTENT(IN) :: imodel
152
153 REAL(mps) :: derlc(nmlyr*2+3)
154 REAL(mps) :: dergl(nmlyr*2+3)
155 INTEGER(mpi) :: label(2)
156 LOGICAL :: ex1
157 LOGICAL :: ex2
158 LOGICAL :: ex3
159 ! for broken lines: 1=fine, 2=coarse
160 dimension nbrl(2),lbrl(nmlyr,2),sbrl(nmlyr,2),wbrl(nmlyr,2), cmbbrl(2)
161 DATA cmbbrl / 0.0, 1.0 / ! cut for combining layers
162 ! ...
163 !CC CALL RNTIME
164 INQUIRE(file='mp2str.txt',iostat=ios,exist=ex1) ! keep, if existing
165 INQUIRE(file='mp2con.txt',iostat=ios,exist=ex2) ! keep, if existing
166
167 INQUIRE(file='mp2tst.bin',iostat=ios,exist=ex3) ! remove, if existing
168
169 WRITE(*,*) ' '
170 WRITE(*,*) 'Generating test data for MP-II ...', imodel
171 WRITE(*,*) ' '
172 ! file management
173 IF(ex3) CALL system('rm mp2tst.bin') ! remove old file
174
175 IF(.NOT.ex1) OPEN(unit=7,access='SEQUENTIAL',form='FORMATTED', &
176 file='mp2str.txt')
177 IF(.NOT.ex2) OPEN(unit=9,access='SEQUENTIAL',form='FORMATTED', &
178 file='mp2con.txt')
179 OPEN(unit=51,access='SEQUENTIAL',form='UNFORMATTED', file='mp2tst.bin')
180
181 s=dets
182 i=0
183 sgn=1.0
184 DO layer=1,10
185 i=i+1
186 islyr(i)=layer ! layer
187 sarc(i)=s ! arclength
188 ssig(i)=sigl ! resolution
189 spro(1,i)=1.0 ! module measures 'X'
190 spro(2,i)=0.0
191 IF (mod(layer,3) == 1) THEN
192 i=i+1
193 islyr(i)=layer ! layer
194 sarc(i)=s+offs ! arclength stereo module
195 ssig(i)=sigl ! resolution
196 spro(1,i)=sqrt(1.0-stereo**2)
197 spro(2,i)=stereo*sgn ! module measures both 'X' and 'Y'
198 sgn=-sgn ! stereo orientation
199 END IF
200 s=s+diss
201 END DO
202
203 ! define broken lines
204 sold=-1000.
205 nbrl(1)=0
206 nbrl(2)=0
207 sbrl=0.
208 wbrl=0.
209 DO k=1,2
210 DO i=1, nmlyr
211 IF (abs(sarc(i)-sold) > cmbbrl(k)) nbrl(k)=nbrl(k)+1
212 lb=nbrl(k)
213 lbrl(i,k)=lb
214 sbrl(lb,k)=sbrl(lb,k)+sarc(i)
215 wbrl(lb,k)=wbrl(lb,k)+1.0
216 sold=sarc(i)
217 END DO
218 DO i=1,nbrl(k)
219 sbrl(i,k)=sbrl(i,k)/wbrl(i,k)
220 wbrl(i,k)=sqrt(wbrl(i,k))
221 END DO
222 END DO
223 ibrl=imodel-2
224
225 ! misalign detector modules -----------------------------------------
226
227 dispxm=0.01 ! module displacement in X 0.1 mm * N(0,1)
228 dispym=0.01 ! module displacement in Y 0.1 mm * N(0,1)
229
230 DO i=0,nlyr-1
231 DO k=0,nmy-1
232 DO l=1,nmx
233 sdevx(((i*nmy+k)*nmx+l))=dispxm*gran() ! shift in x
234 sdevy(((i*nmy+k)*nmx+l))=dispym*gran() ! shift in y
235 END DO
236 END DO
237 END DO
238 ! write text files -------------------------------------------------
239
240 IF(.NOT.ex1) THEN
241 luns=7 ! steerfile
242 WRITE(luns,101) '* *** Default test steering file ***'
243 WRITE(luns,101) ' '
244 WRITE(luns,101) '* Additinal files'
245 WRITE(luns,101) 'fortranfiles ! following bin files are fortran'
246 WRITE(luns,101) 'mp2con.txt ! constraints text file '
247 WRITE(luns,101) 'mp2tst.bin ! binary data file'
248 WRITE(luns,101) 'Cfiles ! following bin files are Cfiles'
249 WRITE(luns,101) ' '
250 WRITE(luns,101) '* Selection of variable global parameters'
251 WRITE(luns,101) '*entries 10 ! lower limit on number of entries/parameter'
252 WRITE(luns,101) '*entries 25 ! default limit on number of entries/parameter'
253 WRITE(luns,101) '*entries 50 ! higher limit on number of entries/parameter'
254 WRITE(luns,101) ' '
255 WRITE(luns,101) '* Initial parameter values, presigmas'
256 WRITE(luns,101) '*Parameter ! define parameter attributes (start of list)'
257 WRITE(luns,101) '*205 -0.01 0. ! start value'
258 WRITE(luns,101) '*206 0.01 0. ! start value'
259 WRITE(luns,101) '*215 0.01 -1. ! fix parameter at value'
260 WRITE(luns,101) '*216 0. -1. ! fix parameter at value'
261 WRITE(luns,101) ' '
262 WRITE(luns,101) '* Handling of outliers, tails etc'
263 WRITE(luns,101) '*hugecut 50.0 !cut factor in iteration 0'
264 WRITE(luns,101) '*chisqcut 30.0 6.0 ! cut factor in iterations 1 and 2'
265 WRITE(luns,101) '*outlierdownweighting 2 ! number of internal iterations (> 1)'
266 WRITE(luns,101) '*dwfractioncut 0.2 ! 0 < value < 0.5'
267 WRITE(luns,101) '*presigma 0.01 ! default value for presigma'
268 WRITE(luns,101) '*regularisation 1.0 ! regularisation factor'
269 WRITE(luns,101) '*regularisation 1.0 0.01 ! regularisation factor, pre-sigma'
270 WRITE(luns,101) ' '
271 WRITE(luns,101) '* Solution methods'
272 WRITE(luns,101) 'method fullMINRES 3 0.01 ! (approximate) MINRES, full storage'
273 WRITE(luns,101) 'method sparseMINRES 3 0.01 ! (approximate) MINRES, sparse storage'
274 WRITE(luns,101) '*mrestol 1.0D-8 ! epsilon for MINRES convergence'
275 WRITE(luns,101) '*bandwidth 0 ! width of MINRES precond. band matrix'
276 WRITE(luns,101) 'method diagonalization 3 0.001 ! diagonalization (-> millepede.eve)'
277 WRITE(luns,101) 'method decomposition 3 0.001 ! Cholesky decomposition'
278 WRITE(luns,101) 'method inversion 3 0.001 ! Gauss matrix inversion'
279 WRITE(luns,101) '* last method is applied'
280 WRITE(luns,101) ' '
281 WRITE(luns,101) '* Additional output. monitoring'
282 WRITE(luns,101) 'printcounts ! print number of entries'
283 WRITE(luns,101) '*monitorresiduals ! poor man''s DMR (-> millepede.mon)'
284 WRITE(luns,101) '*printrecord 1 2 ! debug printout for records'
285 WRITE(luns,101) '*printrecord -1 -1 ! debug printout for bad data records'
286 WRITE(luns,101) ' '
287 WRITE(luns,101) 'end ! optional for end-of-data'
288 END IF
289
290 ! constraints: fix center modules in first/last layer
291
292 nmxy=nmx*nmy
293 lunt=9
294 one=1.0
295 DO i=1,nlyr,nlyr-1
296 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0'
297 DO k=0,nmy-1
298 labelt=(i*nmy+k)*nmx+ncx-1
299 IF(.NOT.ex2) WRITE(lunt,103) labelt,one
300 sdevx(((i-1)*nmy+k)*nmx+ncx)=0.0 ! fix center modules at 0.
301 END DO
302 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0'
303 DO k=0,nmy-1
304 labelt=(i*nmy+k)*nmx+ncx+1000-1
305 IF(.NOT.ex2) WRITE(lunt,103) labelt,one
306 sdevy(((i-1)*nmy+k)*nmx+ncx)=0.0 ! fix center modules at 0.
307 END DO
308 END DO
309
310 ! record loop ------------------------------------------------------
311
312 ncount=10000
313 nthits=0
314 nrecds=0
315
316 DO icount=1,ncount
317 ! 10..100 GeV
318 p=10.0**(1.+uran())
319 the0=sqrt(thck)*0.014/p
320 ip=0
321 ! IF (ICOUNT.LE.3) IP=1
322 CALL genln2(ip) ! generate hits
323
324
325 DO i=1,nhits
326 ! simple straight line
327 lyr=ihits(i)/nmxy+1
328 im =mod(ihits(i),nmxy)
329 nalc=4
330 derlc(1)=spro(1,lyr)
331 derlc(2)=spro(2,lyr)
332 derlc(3)=xhits(i)*spro(1,lyr)
333 derlc(4)=xhits(i)*spro(2,lyr)
334 dergl(1)=spro(1,lyr)
335 dergl(2)=spro(2,lyr)
336 label(1)=im+nmxy*islyr(lyr)
337 label(2)=im+nmxy*islyr(lyr)+1000
338 ! add multiple scattering errors (no correlations)
339 IF (imodel == 1) THEN
340 DO j=i,nhits
341 sigma(j)=sqrt(sigma(j)**2+((xhits(j)-xhits(i))*the0)**2)
342 END DO
343 END IF
344 ! add 'break points' for multiple scattering
345 IF (imodel == 2.AND.i > 1) THEN
346 DO j=1,i-1
347 ! 2 scattering angles from each layer in front of current
348 nalc=nalc+1
349 derlc(nalc)=(xhits(i)-xhits(j))*spro(1,lyr)
350 nalc=nalc+1
351 derlc(nalc)=(xhits(i)-xhits(j))*spro(2,lyr)
352 END DO
353 END IF
354 ! add 'broken lines' offsets for multiple scattering
355 IF (imodel >= 3) THEN
356 nalc=2*nbrl(ibrl)
357 DO k=1, nalc
358 derlc(k)=0.0
359 END DO
360 ! 2 offsets
361 lb=lbrl(lyr,ibrl)
362 derlc(lb*2-1)=spro(1,lyr)
363 derlc(lb*2 )=spro(2,lyr)
364 END IF
365
366 CALL mille(nalc,derlc,2,dergl,label,yhits(i),sigma(i))
367 nthits=nthits+1 ! count hits
368 END DO
369 ! additional measurements from MS
370 IF (imodel == 2) THEN
371 DO i=1,(nhits-1)*2
372 nalc=i+4
373 DO k=1,nalc
374 derlc(k)=0.0
375 END DO
376 derlc(nalc)=1.0
377 CALL mille(nalc,derlc,0,dergl,label,0.0,the0)
378 END DO
379 END IF
380
381 IF (imodel >= 3) THEN
382 DO i=2,nbrl(ibrl)-1
383 dp=1.0/(sbrl(i,ibrl)-sbrl(i-1,ibrl))
384 dn=1.0/(sbrl(i+1,ibrl)-sbrl(i,ibrl))
385 nalc=(i+1)*2
386 DO l=-1,0
387 DO k=1,nalc
388 derlc(k)=0.0
389 END DO
390 derlc(2*(i-1)+l)= dp
391 derlc(2* i +l)=-dp-dn
392 derlc(2*(i+1)+l)= dn
393 CALL mille(nalc,derlc,0,dergl,label,0.0,the0*wbrl(i,ibrl))
394 END DO
395 END DO
396 END IF
397
398 CALL endle
399 nrecds=nrecds+1 ! count records
400 END DO
401
402 ! ------------------------------------------------------------------
403 IF(.NOT.ex1) THEN
404 rewind(7)
405 CLOSE (7)
406 END IF
407 IF(.NOT.ex2) THEN
408 rewind(9)
409 CLOSE (9)
410 END IF
411 rewind(51)
412 CLOSE (51)
413
414 ! WRITE(*,*) ' '
415 ! WRITE(*,*) 'Shifts and drift velocity deviations:'
416 ! DO I=1,NPLAN
417 ! WRITE(*,102) I,DEL(I),DVD(I)
418 ! END DO
419
420
421 WRITE(*,*) ' '
422 WRITE(*,*) ' '
423 WRITE(*,*) ncount,' tracks generated with ',nthits,' hits.'
424 WRITE(*,*) nrecds,' records written.'
425 WRITE(*,*) ' '
426101 FORMAT(a)
427 ! 102 FORMAT(I6,2F10.5)
428103 FORMAT(i8,f10.5)
429END SUBROUTINE mptst2
430
434
435SUBROUTINE genln2(ip)
436 USE mptest2
437
438 IMPLICIT NONE
439 REAL(mps) :: ds
440 REAL(mps) :: dx
441 REAL(mps) :: dy
442 REAL(mps) :: gran
443 INTEGER(mpi) :: i
444 INTEGER(mpi) :: ihit
445 INTEGER(mpi) :: imx
446 INTEGER(mpi) :: imy
447 INTEGER(mpi) :: ioff
448 REAL(mps) :: sold
449 REAL(mps) :: uran
450 REAL(mps) :: x
451 REAL(mps) :: xexit
452 REAL(mps) :: xl
453 REAL(mps) :: xnull
454 REAL(mps) :: xs
455 REAL(mps) :: xslop
456 REAL(mps) :: y
457 REAL(mps) :: yexit
458 REAL(mps) :: yl
459 REAL(mps) :: ynull
460 REAL(mps) :: ys
461 REAL(mps) :: yslop
462
463
464 INTEGER(mpi), INTENT(IN) :: ip
465
466 ! track parameters
467 xnull=sizel*(uran()-0.5) ! uniform vertex
468 ynull=sizel*(uran()-0.5) ! uniform vertex
469 xexit=sizel*(uran()-0.5) ! uniform exit point
470 yexit=sizel*(uran()-0.5) ! uniform exit point
471 xslop=(xexit-xnull)/sarc(nmlyr)
472 yslop=(yexit-ynull)/sarc(nmlyr)
473 IF(ip /= 0) THEN
474 WRITE(*,*) ' '
475 WRITE(*,*) ' Track ', xnull, ynull, xslop, yslop
476 END IF
477
478 nhits=0
479 x=xnull
480 y=ynull
481 dx=xslop
482 dy=yslop
483 sold=0.0
484
485 DO i=1,nmlyr
486 ds=sarc(i)-sold
487 sold=sarc(i)
488 ! position with parameters 1. hit
489 xs=xnull+sarc(i)*xslop
490 ys=ynull+sarc(i)*yslop
491 ! true track position
492 x=x+dx*ds
493 y=y+dy*ds
494 ! multiple scattering
495 dx=dx+gran()*the0
496 dy=dy+gran()*the0
497
498 imx=int((x+sizel*0.5)/sizel*real(nmx,mps),mpi)
499 IF (imx < 0.OR.imx >= nmx) cycle
500 imy=int((y+sizel*0.5)/sizel*real(nmy,mps),mpi)
501 IF (imy < 0.OR.imy >= nmy) cycle
502
503 ihit=((i-1)*nmy+imy)*nmx+imx
504 ioff=((islyr(i)-1)*nmy+imy)*nmx+imx+1
505 nhits=nhits+1
506 ihits(nhits)=ihit
507 xl=x-sdevx(ioff)
508 yl=y-sdevy(ioff)
509 xhits(nhits)=sarc(i)
510 yhits(nhits)=(xl-xs)*spro(1,i)+(yl-ys)*spro(2,i)+gran()*ssig(i)
511 sigma(nhits)=ssig(i)
512
513 IF(ip /= 0) THEN
514 WRITE(*,101) nhits,i,ihit,x,y,xhits(nhits), yhits(nhits),sigma(nhits)
515 END IF
516 END DO
517101 FORMAT(3i3,5f8.4)
518END SUBROUTINE genln2
519
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
mptst2
subroutine mptst2(imodel)
Generate test files.
Definition: mptest2.f90:112
genln2
subroutine genln2(ip)
Generate line and measurements.
Definition: mptest2.f90:436
mpdef
Definition of constants.
Definition: mpdef.f90:24
mpdef::mps
integer, parameter mps
single precision
Definition: mpdef.f90:37
mptest2
Parameters and data.
Definition: mptest2.f90:57
mptest2::ssig
real(mps), dimension(nmlyr) ssig
resolution
Definition: mptest2.f90:85
mptest2::nmx
integer(mpi), parameter nmx
number of modules in x direction
Definition: mptest2.f90:65
mptest2::sizel
real(mps), parameter sizel
size of layers
Definition: mptest2.f90:75
mptest2::sdevx
real(mps), dimension(ntot) sdevx
shift in x (alignment parameter)
Definition: mptest2.f90:82
mptest2::sdevy
real(mps), dimension(ntot) sdevy
shift in y (alignment parameter)
Definition: mptest2.f90:83
mptest2::ncx
integer(mpi), parameter ncx
center (vertical/x) row (used in cons.)
Definition: mptest2.f90:66
mptest2::nhits
integer(mpi) nhits
number of hits
Definition: mptest2.f90:78
mptest2::ihits
integer(mpi), dimension(nmlyr) ihits
module number
Definition: mptest2.f90:81
mptest2::nmlyr
integer(mpi), parameter nmlyr
number of measurement layers
Definition: mptest2.f90:64
mptest2::stereo
real(mps), parameter stereo
stereo angle
Definition: mptest2.f90:74
mptest2::nmy
integer(mpi), parameter nmy
number of modules in y direction
Definition: mptest2.f90:67
mptest2::the0
real(mps) the0
multiple scattering error
Definition: mptest2.f90:79
mptest2::sarc
real(mps), dimension(nmlyr) sarc
arc length
Definition: mptest2.f90:84
mptest2::islyr
integer(mpi), dimension(nmlyr) islyr
(detector) layer
Definition: mptest2.f90:80
mptest2::nlyr
integer(mpi), parameter nlyr
number of detector layers
Definition: mptest2.f90:63
mptest2::sigl
real(mps), parameter sigl
Definition: mptest2.f90:76
mptest2::offs
real(mps), parameter offs
offset of stereo modules
Definition: mptest2.f90:73
mptest2::ntot
integer(mpi), parameter ntot
total number of modules
Definition: mptest2.f90:68
mptest2::thck
real(mps), parameter thck
thickness of plane (X0)
Definition: mptest2.f90:72
mptest2::spro
real(mps), dimension(2, nmlyr) spro
projection of measurent direction in (XY)
Definition: mptest2.f90:86
mptest2::yhits
real(mps), dimension(nmlyr) yhits
measured position in plane (hit)
Definition: mptest2.f90:88
mptest2::sigma
real(mps), dimension(nmlyr) sigma
measurement sigma (hit)
Definition: mptest2.f90:89
mptest2::xhits
real(mps), dimension(nmlyr) xhits
position perp.
Definition: mptest2.f90:87
mptest2::dets
real(mps), parameter dets
arclength of first plane
Definition: mptest2.f90:70
mptest2::diss
real(mps), parameter diss
distance between planes
Definition: mptest2.f90:71
Generated on Wed Sep 18 2024 14:36:52 for Millepede-II by doxygen 1.9.3