 |
Millepede-II
V04-00-00
|
|
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-03-16 Time: 11:08:55 00004 00037 00039 MODULE mptest2 00040 IMPLICIT NONE 00041 SAVE 00042 00043 INTEGER, PARAMETER :: nlyr=10 !< number of detector layers 00044 INTEGER, PARAMETER :: nmlyr=14 !< number of measurement layers 00045 INTEGER, PARAMETER :: nmx=10 !< number of modules in x direction 00046 INTEGER, PARAMETER :: nmy=5 !< number of modules in y direction 00047 INTEGER, PARAMETER :: ntot=nlyr*nmx*nmy !< total number of modules 00048 ! define detector geometry 00049 REAL, PARAMETER :: dets= 10.0 !< arclength of first plane 00050 REAL, PARAMETER :: diss= 10.0 !< distance between planes 00051 REAL, PARAMETER :: thck= 0.02 !< thickness of plane (X0) 00052 REAL, PARAMETER :: offs= 0.5 !< offset of stereo modules 00053 REAL, PARAMETER :: stereo=0.08727 !< stereo angle 00054 REAL, PARAMETER :: sizel= 20.0 !< size of layers 00055 REAL, PARAMETER :: sigl =0.002 ! <resolution 00056 00057 INTEGER :: nhits !< number of hits 00058 REAL :: the0 !< multiple scattering error 00059 INTEGER, DIMENSION(nmlyr) :: islyr !< (detector) layer 00060 INTEGER, DIMENSION(nmlyr) :: ihits !< module number 00061 REAL, DIMENSION(ntot) :: sdevx !< shift in x (alignment parameter) 00062 REAL, DIMENSION(ntot) :: sdevy !< shift in y (alignment parameter) 00063 REAL, DIMENSION(nmlyr) :: sarc !< arc length 00064 REAL, DIMENSION(nmlyr) :: ssig !< resolution 00065 REAL, DIMENSION(2,nmlyr) :: spro !< projection of measurent direction in (XY) 00066 REAL, DIMENSION(nmlyr) :: xhits !< position perp. to plane (hit) 00067 REAL, DIMENSION(nmlyr) :: yhits !< measured position in plane (hit) 00068 REAL, DIMENSION(nmlyr) :: sigma !< measurement sigma (hit) 00069 00070 END MODULE mptest2 00071 00089 00090 SUBROUTINE mptst2(imodel) ! generate test files 00091 USE mptest2 00092 IMPLICIT NONE 00093 REAL :: cmbbrl 00094 REAL :: dispxm 00095 REAL :: dispym 00096 REAL :: dn 00097 REAL :: dp 00098 REAL :: gran 00099 REAL :: one 00100 REAL :: p 00101 REAL :: s 00102 REAL :: sgn 00103 REAL :: sbrl 00104 REAL :: sold 00105 REAL :: uran 00106 REAL :: wbrl 00107 INTEGER :: i 00108 INTEGER :: ibrl 00109 INTEGER :: icount 00110 INTEGER :: im 00111 INTEGER :: ios 00112 INTEGER :: ip 00113 INTEGER :: j 00114 INTEGER :: k 00115 INTEGER :: l 00116 INTEGER :: labelt 00117 INTEGER :: layer 00118 INTEGER :: lb 00119 INTEGER :: lbrl 00120 INTEGER :: luns 00121 INTEGER :: lunt 00122 INTEGER :: lyr 00123 INTEGER :: nalc 00124 INTEGER :: nbrl 00125 INTEGER :: ncount 00126 INTEGER :: ncx 00127 INTEGER :: nmxy 00128 INTEGER :: nrecds 00129 INTEGER :: nthits 00130 00131 INTEGER, INTENT(IN) :: imodel 00132 00133 REAL :: derlc(nmlyr*2+3) 00134 REAL :: dergl(nmlyr*2+3) 00135 INTEGER :: label(2) 00136 LOGICAL :: ex1 00137 LOGICAL :: ex2 00138 LOGICAL :: ex3 00139 ! for broken lines: 1=fine, 2=coarse 00140 DIMENSION nbrl(2),lbrl(nmlyr,2),sbrl(nmlyr,2),wbrl(nmlyr,2), cmbbrl(2) 00141 DATA cmbbrl / 0.0, 1.0 / ! cut for combining layers 00142 ! ... 00143 !CC CALL RNTIME 00144 INQUIRE(FILE='mp2str.txt',IOSTAT=ios,EXIST=ex1) ! keep, if existing 00145 INQUIRE(FILE='mp2con.txt',IOSTAT=ios,EXIST=ex2) ! keep, if existing 00146 00147 INQUIRE(FILE='mp2tst.bin',IOSTAT=ios,EXIST=ex3) ! remove, if existing 00148 00149 WRITE(*,*) ' ' 00150 WRITE(*,*) 'Generating test data for mp II...' 00151 WRITE(*,*) ' ' 00152 ! file management 00153 IF(ex3) CALL system('rm mp2tst.bin') ! remove old file 00154 00155 IF(.NOT.ex1) OPEN(UNIT=7,ACCESS='SEQUENTIAL',FORM='FORMATTED', & 00156 FILE='mp2str.txt') 00157 IF(.NOT.ex2) OPEN(UNIT=9,ACCESS='SEQUENTIAL',FORM='FORMATTED', & 00158 FILE='mp2con.txt') 00159 OPEN(UNIT=51,ACCESS='SEQUENTIAL',FORM='UNFORMATTED', FILE='mp2tst.bin') 00160 00161 s=dets 00162 i=0 00163 sgn=1.0 00164 DO layer=1,10 00165 i=i+1 00166 islyr(i)=layer ! layer 00167 sarc(i)=s ! arclength 00168 ssig(i)=sigl ! resolution 00169 spro(1,i)=1.0 ! module measures 'X' 00170 spro(2,i)=0.0 00171 IF (MOD(layer,3) == 1) THEN 00172 i=i+1 00173 islyr(i)=layer ! layer 00174 sarc(i)=s+offs ! arclength stereo module 00175 ssig(i)=sigl ! resolution 00176 spro(1,i)=SQRT(1.0-stereo**2) 00177 spro(2,i)=stereo*sgn ! module measures both 'X' and 'Y' 00178 sgn=-sgn ! stereo orientation 00179 END IF 00180 s=s+diss 00181 END DO 00182 00183 ! define broken lines 00184 sold=-1000. 00185 nbrl(1)=0 00186 nbrl(2)=0 00187 DO k=1,2 00188 DO i=1, nmlyr 00189 IF (ABS(sarc(i)-sold) > cmbbrl(k)) nbrl(k)=nbrl(k)+1 00190 lb=nbrl(k) 00191 lbrl(i,k)=lb 00192 sbrl(lb,k)=sbrl(lb,k)+sarc(i) 00193 wbrl(lb,k)=wbrl(lb,k)+1.0 00194 sold=sarc(i) 00195 END DO 00196 DO i=1,nbrl(k) 00197 sbrl(i,k)=sbrl(i,k)/wbrl(i,k) 00198 wbrl(i,k)=SQRT(wbrl(i,k)) 00199 END DO 00200 END DO 00201 ibrl=imodel-2 00202 00203 ! misalign detector modules ----------------------------------------- 00204 00205 dispxm=0.01 ! module displacement in X .05 mm * N(0,1) 00206 dispym=0.01 ! module displacement in Y .05 mm * N(0,1) 00207 00208 DO i=0,nlyr-1 00209 DO k=0,nmy-1 00210 DO l=1,nmx 00211 sdevx(((i*nmy+k)*nmx+l))=dispxm*gran() ! shift in x 00212 sdevy(((i*nmy+k)*nmx+l))=dispym*gran() ! shift in y 00213 END DO 00214 END DO 00215 END DO 00216 ! write text files ------------------------------------------------- 00217 00218 IF(.NOT.ex1) THEN 00219 luns=7 ! steerfile 00220 WRITE(luns,101) '* Default test steering file' 00221 WRITE(luns,101) 'fortranfiles ! following bin files are fortran' 00222 WRITE(luns,101) 'mp2con.txt ! constraints text file ' 00223 WRITE(luns,101) 'mp2tst.bin ! binary data file' 00224 WRITE(luns,101) 'Cfiles ! following bin files are Cfiles' 00225 ! WRITE(LUNS,101) '*outlierrejection 100.0 ! reject if Chi^2/Ndf >' 00226 ! WRITE(LUNS,101) '*outliersuppression 3 ! 3 local_fit iterations' 00227 00228 WRITE(luns,101) '*hugecut 50.0 !cut factor in iteration 0' 00229 WRITE(luns,101) '*chisqcut 1.0 1.0 ! cut factor in iterations 1 and 2' 00230 WRITE(luns,101) '*entries 10 ! lower limit on number of entries/parameter' 00231 WRITE(luns,101) & 00232 '*pairentries 10 ! lower limit on number of parameter pairs', & 00233 ' ! (not yet!)' 00234 WRITE(luns,101) '*printrecord 1 2 ! debug printout for records' 00235 WRITE(luns,101) & 00236 '*printrecord -1 -1 ! debug printout for bad data records' 00237 WRITE(luns,101) & 00238 '*outlierdownweighting 2 ! number of internal iterations (> 1)' 00239 WRITE(luns,101) '*dwfractioncut 0.2 ! 0 < value < 0.5' 00240 WRITE(luns,101) '*presigma 0.01 ! default value for presigma' 00241 WRITE(luns,101) '*regularisation 1.0 ! regularisation factor' 00242 WRITE(luns,101) '*regularisation 1.0 0.01 ! regularisation factor, pre-sigma' 00243 00244 WRITE(luns,101) ' ' 00245 WRITE(luns,101) '*bandwidth 0 ! width of precond. band matrix' 00246 WRITE(luns,101) 'method diagonalization 3 0.001 ! diagonalization ' 00247 WRITE(luns,101) 'method fullMINRES 3 0.01 ! minimal residual ' 00248 WRITE(luns,101) 'method sparseMINRES 3 0.01 ! minimal residual ' 00249 WRITE(luns,101) '*mrestol 1.0D-8 ! epsilon for MINRES' 00250 WRITE(luns,101) 'method inversion 3 0.001 ! Gauss matrix inversion' 00251 WRITE(luns,101) '* last method is applied' 00252 WRITE(luns,101) '*matiter 3 ! recalculate matrix in iterations' 00253 WRITE(luns,101) ' ' 00254 WRITE(luns,101) 'end ! optional for end-of-data' 00255 END IF 00256 00257 ! constraints: fix center modules in first/last layer 00258 00259 ncx=(nmx+1)/2 00260 nmxy=nmx*nmy 00261 lunt=9 00262 one=1.0 00263 DO i=1,nlyr,nlyr-1 00264 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0' 00265 DO k=0,nmy-1 00266 labelt=(i*nmy+k)*nmx+ncx-1 00267 IF(.NOT.ex2) WRITE(lunt,103) labelt,one 00268 sdevx(((i-1)*nmy+k)*nmx+ncx)=0.0 ! fix center modules at 0. 00269 END DO 00270 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0' 00271 DO k=0,nmy-1 00272 labelt=(i*nmy+k)*nmx+ncx+1000-1 00273 IF(.NOT.ex2) WRITE(lunt,103) labelt,one 00274 sdevy(((i-1)*nmy+k)*nmx+ncx)=0.0 ! fix center modules at 0. 00275 END DO 00276 END DO 00277 00278 ! record loop ------------------------------------------------------ 00279 00280 ncount=10000 00281 nthits=0 00282 nrecds=0 00283 00284 DO icount=1,ncount 00285 ! 10..100 GeV 00286 p=10.0**(1.+uran()) 00287 the0=SQRT(thck)*0.014/p 00288 ip=0 00289 ! IF (ICOUNT.LE.3) IP=1 00290 CALL genln2(ip) ! generate hits 00291 00292 00293 DO i=1,nhits 00294 ! simple straight line 00295 lyr=ihits(i)/nmxy+1 00296 im =MOD(ihits(i),nmxy) 00297 nalc=4 00298 derlc(1)=spro(1,lyr) 00299 derlc(2)=spro(2,lyr) 00300 derlc(3)=xhits(i)*spro(1,lyr) 00301 derlc(4)=xhits(i)*spro(2,lyr) 00302 dergl(1)=spro(1,lyr) 00303 dergl(2)=spro(2,lyr) 00304 label(1)=im+nmxy*islyr(lyr) 00305 label(2)=im+nmxy*islyr(lyr)+1000 00306 ! add multiple scattering errors (no correlations) 00307 IF (imodel == 1) THEN 00308 DO j=i,nhits 00309 sigma(j)=SQRT(sigma(j)**2+((xhits(j)-xhits(i))*the0)**2) 00310 END DO 00311 END IF 00312 ! add 'break points' for multiple scattering 00313 IF (imodel == 2.AND.i > 1) THEN 00314 DO j=1,i-1 00315 ! 2 scattering angles from each layer in front of current 00316 nalc=nalc+1 00317 derlc(nalc)=(xhits(i)-xhits(j))*spro(1,lyr) 00318 nalc=nalc+1 00319 derlc(nalc)=(xhits(i)-xhits(j))*spro(2,lyr) 00320 END DO 00321 END IF 00322 ! add 'broken lines' offsets for multiple scattering 00323 IF (imodel >= 3) THEN 00324 nalc=2*nbrl(ibrl) 00325 DO k=1, nalc 00326 derlc(k)=0.0 00327 END DO 00328 ! 2 offsets 00329 lb=lbrl(lyr,ibrl) 00330 derlc(lb*2-1)=spro(1,lyr) 00331 derlc(lb*2 )=spro(2,lyr) 00332 END IF 00333 00334 CALL mille(nalc,derlc,2,dergl,label,yhits(i),sigma(i)) 00335 nthits=nthits+1 ! count hits 00336 END DO 00337 ! additional measurements from MS 00338 IF (imodel == 2) THEN 00339 DO i=1,(nhits-1)*2 00340 nalc=i+4 00341 DO k=1,nalc 00342 derlc(k)=0.0 00343 END DO 00344 derlc(nalc)=1.0 00345 CALL mille(nalc,derlc,0,dergl,label,0.0,the0) 00346 END DO 00347 END IF 00348 00349 IF (imodel >= 3) THEN 00350 DO i=2,nbrl(ibrl)-1 00351 dp=1.0/(sbrl(i,ibrl)-sbrl(i-1,ibrl)) 00352 dn=1.0/(sbrl(i+1,ibrl)-sbrl(i,ibrl)) 00353 nalc=(i+1)*2 00354 DO l=-1,0 00355 DO k=1,nalc 00356 derlc(k)=0.0 00357 END DO 00358 derlc(2*(i-1)+l)= dp 00359 derlc(2* i +l)=-dp-dn 00360 derlc(2*(i+1)+l)= dn 00361 CALL mille(nalc,derlc,0,dergl,label,0.0,the0*wbrl(i,ibrl)) 00362 END DO 00363 END DO 00364 END IF 00365 00366 CALL endle 00367 nrecds=nrecds+1 ! count records 00368 END DO 00369 00370 ! ------------------------------------------------------------------ 00371 IF(.NOT.ex1) THEN 00372 REWIND (7) 00373 CLOSE (7) 00374 END IF 00375 IF(.NOT.ex2) THEN 00376 REWIND (9) 00377 CLOSE (9) 00378 END IF 00379 REWIND (51) 00380 CLOSE (51) 00381 00382 ! WRITE(*,*) ' ' 00383 ! WRITE(*,*) 'Shifts and drift velocity deviations:' 00384 ! DO I=1,NPLAN 00385 ! WRITE(*,102) I,DEL(I),DVD(I) 00386 ! END DO 00387 00388 00389 WRITE(*,*) ' ' 00390 WRITE(*,*) ' ' 00391 WRITE(*,*) ncount,' tracks generated with ',nthits,' hits.' 00392 WRITE(*,*) nrecds,' records written.' 00393 WRITE(*,*) ' ' 00394 101 FORMAT(a) 00395 ! 102 FORMAT(I6,2F10.5) 00396 103 FORMAT(i8,f10.5) 00397 END SUBROUTINE mptst2 00398 00402 00403 SUBROUTINE genln2(ip) 00404 USE mptest2 00405 00406 IMPLICIT NONE 00407 REAL:: ds 00408 REAL:: dx 00409 REAL:: dy 00410 REAL:: gran 00411 INTEGER:: i 00412 INTEGER:: ihit 00413 INTEGER:: imx 00414 INTEGER:: imy 00415 INTEGER:: ioff 00416 REAL:: sold 00417 REAL:: uran 00418 REAL:: x 00419 REAL:: xexit 00420 REAL:: xl 00421 REAL:: xnull 00422 REAL:: xs 00423 REAL:: xslop 00424 REAL:: y 00425 REAL:: yexit 00426 REAL:: yl 00427 REAL:: ynull 00428 REAL:: ys 00429 REAL:: yslop 00430 00431 00432 INTEGER, INTENT(IN) :: ip 00433 00434 ! track parameters 00435 xnull=sizel*(uran()-0.5) ! uniform vertex 00436 ynull=sizel*(uran()-0.5) ! uniform vertex 00437 xexit=sizel*(uran()-0.5) ! uniform exit point 00438 yexit=sizel*(uran()-0.5) ! uniform exit point 00439 xslop=(xexit-xnull)/sarc(nmlyr) 00440 yslop=(yexit-ynull)/sarc(nmlyr) 00441 IF(ip /= 0) THEN 00442 WRITE(*,*) ' ' 00443 WRITE(*,*) ' Track ', xnull, ynull, xslop, yslop 00444 END IF 00445 00446 nhits=0 00447 x=xnull 00448 y=ynull 00449 dx=xslop 00450 dy=yslop 00451 sold=0.0 00452 00453 DO i=1,nmlyr 00454 ds=sarc(i)-sold 00455 sold=sarc(i) 00456 ! position with parameters 1. hit 00457 xs=xnull+sarc(i)*xslop 00458 ys=ynull+sarc(i)*yslop 00459 ! true track position 00460 x=x+dx*ds 00461 y=y+dy*ds 00462 ! multiple scattering 00463 dx=dx+gran()*the0 00464 dy=dy+gran()*the0 00465 00466 imx=IFIX((x+sizel*0.5)/sizel*FLOAT(nmx)) 00467 IF (imx < 0.OR.imx >= nmx) CYCLE 00468 imy=IFIX((y+sizel*0.5)/sizel*FLOAT(nmy)) 00469 IF (imy < 0.OR.imy >= nmy) CYCLE 00470 00471 ihit=((i-1)*nmy+imy)*nmx+imx 00472 ioff=((islyr(i)-1)*nmy+imy)*nmx+imx+1 00473 nhits=nhits+1 00474 ihits(nhits)=ihit 00475 xl=x-sdevx(ioff) 00476 yl=y-sdevy(ioff) 00477 xhits(nhits)=sarc(i) 00478 yhits(nhits)=(xl-xs)*spro(1,i)+(yl-ys)*spro(2,i)+gran()*ssig(i) 00479 sigma(nhits)=ssig(i) 00480 00481 IF(ip /= 0) THEN 00482 WRITE(*,101) nhits,i,ihit,x,y,xhits(nhits), yhits(nhits),sigma(nhits) 00483 END IF 00484 END DO 00485 101 FORMAT(3I3,5F8.4) 00486 END SUBROUTINE genln2 00487
1.8.0
