SND@LHC Software
Loading...
Searching...
No Matches
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 :: nmy=5
67 INTEGER(mpi), PARAMETER :: ntot=nlyr*nmx*nmy
68 ! define detector geometry
69 REAL(mps), PARAMETER :: dets= 10.0
70 REAL(mps), PARAMETER :: diss= 10.0
71 REAL(mps), PARAMETER :: thck= 0.02
72 REAL(mps), PARAMETER :: offs= 0.5
73 REAL(mps), PARAMETER :: stereo=0.08727
74 REAL(mps), PARAMETER :: sizel= 20.0
75 REAL(mps), PARAMETER :: sigl =0.002 ! <resolution
76
77 INTEGER(mpi) :: nhits
78 REAL(mps) :: the0
79 INTEGER(mpi), DIMENSION(nmlyr) :: islyr
80 INTEGER(mpi), DIMENSION(nmlyr) :: ihits
81 REAL(mps), DIMENSION(ntot) :: sdevx
82 REAL(mps), DIMENSION(ntot) :: sdevy
83 REAL(mps), DIMENSION(nmlyr) :: sarc
84 REAL(mps), DIMENSION(nmlyr) :: ssig
85 REAL(mps), DIMENSION(2,nmlyr) :: spro
86 REAL(mps), DIMENSION(nmlyr) :: xhits
87 REAL(mps), DIMENSION(nmlyr) :: yhits
88 REAL(mps), DIMENSION(nmlyr) :: sigma
89
90END MODULE mptest2
91
109
110SUBROUTINE mptst2(imodel) ! generate test files
111 USE mptest2
112 IMPLICIT NONE
113 REAL(mps) :: cmbbrl
114 REAL(mps) :: dispxm
115 REAL(mps) :: dispym
116 REAL(mps) :: dn
117 REAL(mps) :: dp
118 REAL(mps) :: gran
119 REAL(mps) :: one
120 REAL(mps) :: p
121 REAL(mps) :: s
122 REAL(mps) :: sgn
123 REAL(mps) :: sbrl
124 REAL(mps) :: sold
125 REAL(mps) :: uran
126 REAL(mps) :: wbrl
127 INTEGER(mpi) :: i
128 INTEGER(mpi) :: ibrl
129 INTEGER(mpi) :: icount
130 INTEGER(mpi) :: im
131 INTEGER(mpi) :: ios
132 INTEGER(mpi) :: ip
133 INTEGER(mpi) :: j
134 INTEGER(mpi) :: k
135 INTEGER(mpi) :: l
136 INTEGER(mpi) :: labelt
137 INTEGER(mpi) :: layer
138 INTEGER(mpi) :: lb
139 INTEGER(mpi) :: lbrl
140 INTEGER(mpi) :: luns
141 INTEGER(mpi) :: lunt
142 INTEGER(mpi) :: lyr
143 INTEGER(mpi) :: nalc
144 INTEGER(mpi) :: nbrl
145 INTEGER(mpi) :: ncount
146 INTEGER(mpi) :: ncx
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...'
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 DO k=1,2
208 DO i=1, nmlyr
209 IF (abs(sarc(i)-sold) > cmbbrl(k)) nbrl(k)=nbrl(k)+1
210 lb=nbrl(k)
211 lbrl(i,k)=lb
212 sbrl(lb,k)=sbrl(lb,k)+sarc(i)
213 wbrl(lb,k)=wbrl(lb,k)+1.0
214 sold=sarc(i)
215 END DO
216 DO i=1,nbrl(k)
217 sbrl(i,k)=sbrl(i,k)/wbrl(i,k)
218 wbrl(i,k)=sqrt(wbrl(i,k))
219 END DO
220 END DO
221 ibrl=imodel-2
222
223 ! misalign detector modules -----------------------------------------
224
225 dispxm=0.01 ! module displacement in X .05 mm * N(0,1)
226 dispym=0.01 ! module displacement in Y .05 mm * N(0,1)
227
228 DO i=0,nlyr-1
229 DO k=0,nmy-1
230 DO l=1,nmx
231 sdevx(((i*nmy+k)*nmx+l))=dispxm*gran() ! shift in x
232 sdevy(((i*nmy+k)*nmx+l))=dispym*gran() ! shift in y
233 END DO
234 END DO
235 END DO
236 ! write text files -------------------------------------------------
237
238 IF(.NOT.ex1) THEN
239 luns=7 ! steerfile
240 WRITE(luns,101) '* Default test steering file'
241 WRITE(luns,101) 'fortranfiles ! following bin files are fortran'
242 WRITE(luns,101) 'mp2con.txt ! constraints text file '
243 WRITE(luns,101) 'mp2tst.bin ! binary data file'
244 WRITE(luns,101) 'Cfiles ! following bin files are Cfiles'
245 ! WRITE(LUNS,101) '*outlierrejection 100.0 ! reject if Chi^2/Ndf >'
246 ! WRITE(LUNS,101) '*outliersuppression 3 ! 3 local_fit iterations'
247
248 WRITE(luns,101) '*hugecut 50.0 !cut factor in iteration 0'
249 WRITE(luns,101) '*chisqcut 1.0 1.0 ! cut factor in iterations 1 and 2'
250 WRITE(luns,101) '*entries 10 ! lower limit on number of entries/parameter'
251 WRITE(luns,101) &
252 '*pairentries 10 ! lower limit on number of parameter pairs', &
253 ' ! (not yet!)'
254 WRITE(luns,101) '*printrecord 1 2 ! debug printout for records'
255 WRITE(luns,101) &
256 '*printrecord -1 -1 ! debug printout for bad data records'
257 WRITE(luns,101) &
258 '*outlierdownweighting 2 ! number of internal iterations (> 1)'
259 WRITE(luns,101) '*dwfractioncut 0.2 ! 0 < value < 0.5'
260 WRITE(luns,101) '*presigma 0.01 ! default value for presigma'
261 WRITE(luns,101) '*regularisation 1.0 ! regularisation factor'
262 WRITE(luns,101) '*regularisation 1.0 0.01 ! regularisation factor, pre-sigma'
263
264 WRITE(luns,101) ' '
265 WRITE(luns,101) '*bandwidth 0 ! width of precond. band matrix'
266 WRITE(luns,101) 'method diagonalization 3 0.001 ! diagonalization '
267 WRITE(luns,101) 'method fullMINRES 3 0.01 ! minimal residual '
268 WRITE(luns,101) 'method sparseMINRES 3 0.01 ! minimal residual '
269 WRITE(luns,101) '*mrestol 1.0D-8 ! epsilon for MINRES'
270 WRITE(luns,101) 'method inversion 3 0.001 ! Gauss matrix inversion'
271 WRITE(luns,101) '* last method is applied'
272 WRITE(luns,101) '*matiter 3 ! recalculate matrix in iterations'
273 WRITE(luns,101) ' '
274 WRITE(luns,101) 'end ! optional for end-of-data'
275 END IF
276
277 ! constraints: fix center modules in first/last layer
278
279 ncx=(nmx+1)/2
280 nmxy=nmx*nmy
281 lunt=9
282 one=1.0
283 DO i=1,nlyr,nlyr-1
284 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0'
285 DO k=0,nmy-1
286 labelt=(i*nmy+k)*nmx+ncx-1
287 IF(.NOT.ex2) WRITE(lunt,103) labelt,one
288 sdevx(((i-1)*nmy+k)*nmx+ncx)=0.0 ! fix center modules at 0.
289 END DO
290 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0'
291 DO k=0,nmy-1
292 labelt=(i*nmy+k)*nmx+ncx+1000-1
293 IF(.NOT.ex2) WRITE(lunt,103) labelt,one
294 sdevy(((i-1)*nmy+k)*nmx+ncx)=0.0 ! fix center modules at 0.
295 END DO
296 END DO
297
298 ! record loop ------------------------------------------------------
299
300 ncount=10000
301 nthits=0
302 nrecds=0
303
304 DO icount=1,ncount
305 ! 10..100 GeV
306 p=10.0**(1.+uran())
307 the0=sqrt(thck)*0.014/p
308 ip=0
309 ! IF (ICOUNT.LE.3) IP=1
310 CALL genln2(ip) ! generate hits
311
312
313 DO i=1,nhits
314 ! simple straight line
315 lyr=ihits(i)/nmxy+1
316 im =mod(ihits(i),nmxy)
317 nalc=4
318 derlc(1)=spro(1,lyr)
319 derlc(2)=spro(2,lyr)
320 derlc(3)=xhits(i)*spro(1,lyr)
321 derlc(4)=xhits(i)*spro(2,lyr)
322 dergl(1)=spro(1,lyr)
323 dergl(2)=spro(2,lyr)
324 label(1)=im+nmxy*islyr(lyr)
325 label(2)=im+nmxy*islyr(lyr)+1000
326 ! add multiple scattering errors (no correlations)
327 IF (imodel == 1) THEN
328 DO j=i,nhits
329 sigma(j)=sqrt(sigma(j)**2+((xhits(j)-xhits(i))*the0)**2)
330 END DO
331 END IF
332 ! add 'break points' for multiple scattering
333 IF (imodel == 2.AND.i > 1) THEN
334 DO j=1,i-1
335 ! 2 scattering angles from each layer in front of current
336 nalc=nalc+1
337 derlc(nalc)=(xhits(i)-xhits(j))*spro(1,lyr)
338 nalc=nalc+1
339 derlc(nalc)=(xhits(i)-xhits(j))*spro(2,lyr)
340 END DO
341 END IF
342 ! add 'broken lines' offsets for multiple scattering
343 IF (imodel >= 3) THEN
344 nalc=2*nbrl(ibrl)
345 DO k=1, nalc
346 derlc(k)=0.0
347 END DO
348 ! 2 offsets
349 lb=lbrl(lyr,ibrl)
350 derlc(lb*2-1)=spro(1,lyr)
351 derlc(lb*2 )=spro(2,lyr)
352 END IF
353
354 CALL mille(nalc,derlc,2,dergl,label,yhits(i),sigma(i))
355 nthits=nthits+1 ! count hits
356 END DO
357 ! additional measurements from MS
358 IF (imodel == 2) THEN
359 DO i=1,(nhits-1)*2
360 nalc=i+4
361 DO k=1,nalc
362 derlc(k)=0.0
363 END DO
364 derlc(nalc)=1.0
365 CALL mille(nalc,derlc,0,dergl,label,0.0,the0)
366 END DO
367 END IF
368
369 IF (imodel >= 3) THEN
370 DO i=2,nbrl(ibrl)-1
371 dp=1.0/(sbrl(i,ibrl)-sbrl(i-1,ibrl))
372 dn=1.0/(sbrl(i+1,ibrl)-sbrl(i,ibrl))
373 nalc=(i+1)*2
374 DO l=-1,0
375 DO k=1,nalc
376 derlc(k)=0.0
377 END DO
378 derlc(2*(i-1)+l)= dp
379 derlc(2* i +l)=-dp-dn
380 derlc(2*(i+1)+l)= dn
381 CALL mille(nalc,derlc,0,dergl,label,0.0,the0*wbrl(i,ibrl))
382 END DO
383 END DO
384 END IF
385
386 CALL endle
387 nrecds=nrecds+1 ! count records
388 END DO
389
390 ! ------------------------------------------------------------------
391 IF(.NOT.ex1) THEN
392 rewind(7)
393 CLOSE (7)
394 END IF
395 IF(.NOT.ex2) THEN
396 rewind(9)
397 CLOSE (9)
398 END IF
399 rewind(51)
400 CLOSE (51)
401
402 ! WRITE(*,*) ' '
403 ! WRITE(*,*) 'Shifts and drift velocity deviations:'
404 ! DO I=1,NPLAN
405 ! WRITE(*,102) I,DEL(I),DVD(I)
406 ! END DO
407
408
409 WRITE(*,*) ' '
410 WRITE(*,*) ' '
411 WRITE(*,*) ncount,' tracks generated with ',nthits,' hits.'
412 WRITE(*,*) nrecds,' records written.'
413 WRITE(*,*) ' '
414101 FORMAT(a)
415 ! 102 FORMAT(I6,2F10.5)
416103 FORMAT(i8,f10.5)
417END SUBROUTINE mptst2
418
422
423SUBROUTINE genln2(ip)
424 USE mptest2
425
426 IMPLICIT NONE
427 REAL(mps) :: ds
428 REAL(mps) :: dx
429 REAL(mps) :: dy
430 REAL(mps) :: gran
431 INTEGER(mpi) :: i
432 INTEGER(mpi) :: ihit
433 INTEGER(mpi) :: imx
434 INTEGER(mpi) :: imy
435 INTEGER(mpi) :: ioff
436 REAL(mps) :: sold
437 REAL(mps) :: uran
438 REAL(mps) :: x
439 REAL(mps) :: xexit
440 REAL(mps) :: xl
441 REAL(mps) :: xnull
442 REAL(mps) :: xs
443 REAL(mps) :: xslop
444 REAL(mps) :: y
445 REAL(mps) :: yexit
446 REAL(mps) :: yl
447 REAL(mps) :: ynull
448 REAL(mps) :: ys
449 REAL(mps) :: yslop
450
451
452 INTEGER(mpi), INTENT(IN) :: ip
453
454 ! track parameters
455 xnull=sizel*(uran()-0.5) ! uniform vertex
456 ynull=sizel*(uran()-0.5) ! uniform vertex
457 xexit=sizel*(uran()-0.5) ! uniform exit point
458 yexit=sizel*(uran()-0.5) ! uniform exit point
459 xslop=(xexit-xnull)/sarc(nmlyr)
460 yslop=(yexit-ynull)/sarc(nmlyr)
461 IF(ip /= 0) THEN
462 WRITE(*,*) ' '
463 WRITE(*,*) ' Track ', xnull, ynull, xslop, yslop
464 END IF
465
466 nhits=0
467 x=xnull
468 y=ynull
469 dx=xslop
470 dy=yslop
471 sold=0.0
472
473 DO i=1,nmlyr
474 ds=sarc(i)-sold
475 sold=sarc(i)
476 ! position with parameters 1. hit
477 xs=xnull+sarc(i)*xslop
478 ys=ynull+sarc(i)*yslop
479 ! true track position
480 x=x+dx*ds
481 y=y+dy*ds
482 ! multiple scattering
483 dx=dx+gran()*the0
484 dy=dy+gran()*the0
485
486 imx=int((x+sizel*0.5)/sizel*real(nmx,mps),mpi)
487 IF (imx < 0.OR.imx >= nmx) cycle
488 imy=int((y+sizel*0.5)/sizel*real(nmy,mps),mpi)
489 IF (imy < 0.OR.imy >= nmy) cycle
490
491 ihit=((i-1)*nmy+imy)*nmx+imx
492 ioff=((islyr(i)-1)*nmy+imy)*nmx+imx+1
493 nhits=nhits+1
494 ihits(nhits)=ihit
495 xl=x-sdevx(ioff)
496 yl=y-sdevy(ioff)
497 xhits(nhits)=sarc(i)
498 yhits(nhits)=(xl-xs)*spro(1,i)+(yl-ys)*spro(2,i)+gran()*ssig(i)
499 sigma(nhits)=ssig(i)
500
501 IF(ip /= 0) THEN
502 WRITE(*,101) nhits,i,ihit,x,y,xhits(nhits), yhits(nhits),sigma(nhits)
503 END IF
504 END DO
505101 FORMAT(3i3,5f8.4)
506END SUBROUTINE genln2
507
mille
subroutine mille(nlc, derlc, ngl, dergl, label, rmeas, sigma)
Add data block to record. Called from user code.
Definition mille.f90:76
mptst2
subroutine mptst2(imodel)
Generate test files.
Definition mptest2.f90:111
genln2
subroutine genln2(ip)
Generate line and measurements.
Definition mptest2.f90:424
mpdef
Definition of constants.
Definition mpdef.f90:24
mpdef::mps
integer, parameter mps
Definition mpdef.f90:37
mptest2
Parameters and data.
Definition mptest2.f90:57
mptest2::ssig
real(mps), dimension(nmlyr) ssig
resolution
Definition mptest2.f90:84
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:74
mptest2::sdevx
real(mps), dimension(ntot) sdevx
shift in x (alignment parameter)
Definition mptest2.f90:81
mptest2::sdevy
real(mps), dimension(ntot) sdevy
shift in y (alignment parameter)
Definition mptest2.f90:82
mptest2::nhits
integer(mpi) nhits
number of hits
Definition mptest2.f90:77
mptest2::ihits
integer(mpi), dimension(nmlyr) ihits
module number
Definition mptest2.f90:80
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:73
mptest2::nmy
integer(mpi), parameter nmy
number of modules in y direction
Definition mptest2.f90:66
mptest2::the0
real(mps) the0
multiple scattering error
Definition mptest2.f90:78
mptest2::sarc
real(mps), dimension(nmlyr) sarc
arc length
Definition mptest2.f90:83
mptest2::islyr
integer(mpi), dimension(nmlyr) islyr
(detector) layer
Definition mptest2.f90:79
mptest2::nlyr
integer(mpi), parameter nlyr
number of detector layers
Definition mptest2.f90:63
mptest2::sigl
real(mps), parameter sigl
Definition mptest2.f90:75
mptest2::offs
real(mps), parameter offs
offset of stereo modules
Definition mptest2.f90:72
mptest2::ntot
integer(mpi), parameter ntot
total number of modules
Definition mptest2.f90:67
mptest2::thck
real(mps), parameter thck
thickness of plane (X0)
Definition mptest2.f90:71
mptest2::spro
real(mps), dimension(2, nmlyr) spro
projection of measurent direction in (XY)
Definition mptest2.f90:85
mptest2::yhits
real(mps), dimension(nmlyr) yhits
measured position in plane (hit)
Definition mptest2.f90:87
mptest2::sigma
real(mps), dimension(nmlyr) sigma
measurement sigma (hit)
Definition mptest2.f90:88
mptest2::xhits
real(mps), dimension(nmlyr) xhits
position perp. to plane (hit)
Definition mptest2.f90:86
mptest2::dets
real(mps), parameter dets
arclength of first plane
Definition mptest2.f90:69
mptest2::diss
real(mps), parameter diss
distance between planes
Definition mptest2.f90:70