 |
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:48 00004 00013 00014 00016 MODULE mptest1 00017 IMPLICIT NONE 00018 SAVE 00019 00020 INTEGER, PARAMETER :: nplan=100 00021 00022 ! define detector geometry 00023 REAL, PARAMETER :: detx= 10.0 !< x-value of first plane 00024 REAL, PARAMETER :: disx= 10.0 !< distance between planes 00025 REAL, PARAMETER :: thck= 2.0 !< thickness of plane 00026 REAL, PARAMETER :: heit=100.0 !< height of detector plane 00027 REAL, PARAMETER :: effp=0.90 !< plane efficiency 00028 REAL, PARAMETER :: sgmp=0.0150 !< measurement sigma 00029 00030 ! misalignment 00031 REAL, DIMENSION(nplan) :: del !< shift (position deviation) (alignment parameter) 00032 REAL, DIMENSION(nplan) :: dvd !< rel. drift velocity deviation (calibration parameter) 00033 ! track parameter 00034 REAL :: ynull !< track position at vertex 00035 REAL :: slope !< track slope 00036 00037 INTEGER :: nhits !< number of hits 00038 INTEGER, DIMENSION(nplan) :: ihits !< plane numbers (planes with hits) 00039 REAL, DIMENSION(nplan) :: eff !< plane efficiency 00040 REAL, DIMENSION(nplan) :: sgm !< measurement sigma (plane) 00041 REAL, DIMENSION(nplan) :: ydrft !< signed drift length 00042 REAL, DIMENSION(nplan) :: xhits !< position perp. to plane (hit) 00043 REAL, DIMENSION(nplan) :: yhits !< measured position in plane (hit) 00044 REAL, DIMENSION(nplan) :: sigma !< measurement sigma (hit) 00045 00046 END MODULE mptest1 00047 00056 00057 SUBROUTINE mptest 00058 USE mptest1 00059 00060 IMPLICIT NONE 00061 REAL :: dbar 00062 REAL :: det 00063 REAL :: displ 00064 REAL :: drift 00065 REAL :: eps 00066 REAL :: eta 00067 REAL :: gran 00068 REAL :: one 00069 REAL :: ww 00070 REAL :: x 00071 REAL :: xbar 00072 INTEGER :: i 00073 INTEGER :: icount 00074 INTEGER :: ios 00075 INTEGER :: ip 00076 INTEGER :: ipl 00077 INTEGER :: labelt 00078 INTEGER :: luns 00079 INTEGER :: lunt 00080 INTEGER :: ncount 00081 INTEGER :: nrecds 00082 INTEGER :: nthits 00083 00084 DOUBLE PRECISION :: s1 00085 DOUBLE PRECISION :: s2 00086 DOUBLE PRECISION :: sw 00087 DOUBLE PRECISION :: sv 00088 DOUBLE PRECISION :: sum1 00089 DOUBLE PRECISION :: sum2 00090 REAL :: derlc(2) 00091 REAL :: dergl(2) 00092 INTEGER :: label(2) 00093 LOGICAL :: ex1 00094 LOGICAL :: ex2 00095 LOGICAL :: ex3 00096 ! ... 00097 !CC CALL RNTIME 00098 INQUIRE(FILE='mp2str.txt',IOSTAT=ios,EXIST=ex1) ! keep, if existing 00099 INQUIRE(FILE='mp2con.txt',IOSTAT=ios,EXIST=ex2) ! keep, if existing 00100 00101 INQUIRE(FILE='mp2tst.bin',IOSTAT=ios,EXIST=ex3) ! remove, if existing 00102 00103 WRITE(*,*) ' ' 00104 WRITE(*,*) 'Generating test data for mp II...' 00105 WRITE(*,*) ' ' 00106 ! file management 00107 IF(ex3) CALL system('rm mp2tst.bin') ! remove old file 00108 00109 IF(.NOT.ex1) OPEN(UNIT=7,ACCESS='SEQUENTIAL',FORM='FORMATTED', & 00110 FILE='mp2str.txt') 00111 IF(.NOT.ex2) OPEN(UNIT=9,ACCESS='SEQUENTIAL',FORM='FORMATTED', & 00112 FILE='mp2con.txt') 00113 OPEN(UNIT=51,ACCESS='SEQUENTIAL',FORM='UNFORMATTED', FILE='mp2tst.bin') 00114 00115 DO i=1,nplan 00116 eff(i)=effp ! plane efficiency 00117 sgm(i)=sgmp ! measurement sigma 00118 del(i)=0.0 ! true shift is zero 00119 END DO 00120 00121 ipl=7 ! modify one plane (7) 00122 eff(ipl)=0.1 ! low efficiency 00123 sgm(ipl)=0.0400 ! bad resolution 00124 00125 ! misalign detector planes ----------------------------------------- 00126 00127 displ=0.1 ! displacement 1 mm * N(0,1) 00128 drift=0.02 ! Vdrift deviation 2 % * N(0,1) 00129 DO i=1,nplan 00130 del(i)=displ*gran() ! shift 00131 dvd(i)=drift*gran() ! rel. drift velocitu deviation 00132 END DO 00133 del(10)=0.0 ! no shift 00134 del(90)=0.0 ! no shift 00135 00136 ! write text files ------------------------------------------------- 00137 00138 IF(.NOT.ex1) THEN 00139 luns=7 ! steerfile 00140 WRITE(luns,101) '* Default test steering file' 00141 WRITE(luns,101) 'fortranfiles ! following bin files are fortran' 00142 WRITE(luns,101) 'mp2con.txt ! constraints text file ' 00143 WRITE(luns,101) 'mp2tst.bin ! binary data file' 00144 WRITE(luns,101) 'Cfiles ! following bin files are Cfiles' 00145 ! WRITE(LUNS,101) '*outlierrejection 100.0 ! reject if Chi^2/Ndf >' 00146 ! WRITE(LUNS,101) '*outliersuppression 3 ! 3 local_fit iterations' 00147 00148 WRITE(luns,101) '*hugecut 50.0 !cut factor in iteration 0' 00149 WRITE(luns,101) '*chisqcut 1.0 1.0 ! cut factor in iterations 1 and 2' 00150 WRITE(luns,101) '*entries 10 ! lower limit on number of entries/parameter' 00151 WRITE(luns,101) & 00152 '*pairentries 10 ! lower limit on number of parameter pairs', & 00153 ' ! (not yet!)' 00154 WRITE(luns,101) '*printrecord 1 2 ! debug printout for records' 00155 WRITE(luns,101) & 00156 '*printrecord -1 -1 ! debug printout for bad data records' 00157 WRITE(luns,101) & 00158 '*outlierdownweighting 2 ! number of internal iterations (> 1)' 00159 WRITE(luns,101) '*dwfractioncut 0.2 ! 0 < value < 0.5' 00160 WRITE(luns,101) '*presigma 0.01 ! default value for presigma' 00161 WRITE(luns,101) '*regularisation 1.0 ! regularisation factor' 00162 WRITE(luns,101) '*regularisation 1.0 0.01 ! regularisation factor, pre-sigma' 00163 00164 WRITE(luns,101) ' ' 00165 WRITE(luns,101) '*bandwidth 0 ! width of precond. band matrix' 00166 WRITE(luns,101) 'method diagonalization 3 0.001 ! diagonalization ' 00167 WRITE(luns,101) 'method fullMINRES 3 0.01 ! minimal residual ' 00168 WRITE(luns,101) 'method sparseMINRES 3 0.01 ! minimal residual ' 00169 WRITE(luns,101) '*mrestol 1.0D-8 ! epsilon for MINRES' 00170 WRITE(luns,101) 'method inversion 3 0.001 ! Gauss matrix inversion' 00171 WRITE(luns,101) '* last method is applied' 00172 WRITE(luns,101) '*matiter 3 ! recalculate matrix in iterations' 00173 WRITE(luns,101) ' ' 00174 WRITE(luns,101) 'end ! optional for end-of-data' 00175 ENDIF 00176 00177 lunt=9 ! constraint file 00178 one=1.0 ! shift constraint 00179 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0' 00180 DO i=1,nplan 00181 labelt=10+i*2 00182 x=detx+FLOAT(i-1)*disx+0.5*thck 00183 IF(.NOT.ex2) WRITE(lunt,103) labelt,one 00184 END DO 00185 00186 sw=0.0D0 ! tilt constraint 00187 sv=0.0D0 00188 s1=0.0D0 00189 s2=0.0D0 00190 IF(.NOT.ex2) WRITE(lunt,*) 'Constraint 0.0' ! write 00191 dbar=0.5*FLOAT(nplan-1)*disx 00192 xbar=detx+0.5*FLOAT(nplan-1)*disx! +0.5*THCK 00193 DO i=1,nplan 00194 labelt=10+i*2 00195 x=detx+FLOAT(i-1)*disx !+0.5*THCK 00196 ww=(x-xbar)/dbar 00197 IF(.NOT.ex2) WRITE(lunt,103) labelt,ww ! write 00198 s1=s1+del(i) 00199 s2=s2+ww*del(i) 00200 sw=sw+ww 00201 sv=sv+ww*ww 00202 END DO 00203 00204 00205 det=REAL(FLOAT(nplan)*sv-sw*sw) 00206 eps=REAL(sv*s1-sw*s2)/det 00207 eta=REAL(FLOAT(nplan)*s2-sw*s1)/det 00208 DO i=1,nplan 00209 x=detx+FLOAT(i-1)*disx 00210 ww=(x-xbar)/dbar 00211 del(i)=del(i)-eps-eta*ww ! correct displacement ... 00212 END DO ! ... for constraints 00213 00214 sum1=0.0 00215 sum2=0.0 00216 DO i=1,nplan 00217 sum1=sum1+del(i) 00218 x=detx+FLOAT(i-1)*disx !+0.5*THCK 00219 ww=(x-xbar)/dbar 00220 sum2=sum2+del(i)*ww 00221 END DO 00222 ! WRITE(*,*) ' Check for constraints ',SUM1,SUM2 00223 00224 ! record loop ------------------------------------------------------ 00225 00226 ncount=10000 00227 nthits=0 00228 nrecds=0 00229 00230 DO icount=1,ncount 00231 ip=0 00232 IF(icount == 8759) ip=1 00233 ! IF(ICOUNT.EQ.6309) IP=1 00234 ! IF(ICOUNT.EQ.7468) IP=1 00235 CALL genlin(ip) ! generate hits 00236 00237 DO i=1,nhits 00238 derlc(1)=1.0 00239 derlc(2)=xhits(i) 00240 dergl(1)=1.0 00241 dergl(2)=ydrft(i) 00242 label(1)=10+ihits(i)*2 00243 label(2)=500 + ihits(i) 00244 CALL mille(2,derlc,2,dergl,label,yhits(i),sigma(i)) 00245 nthits=nthits+1 ! count hits 00246 END DO 00247 CALL endle 00248 nrecds=nrecds+1 ! count records 00249 END DO 00250 00251 ! ------------------------------------------------------------------ 00252 IF(.NOT.ex1) THEN 00253 REWIND (7) 00254 CLOSE (7) 00255 END IF 00256 IF(.NOT.ex2) THEN 00257 REWIND (9) 00258 CLOSE (9) 00259 END IF 00260 REWIND (51) 00261 CLOSE (51) 00262 00263 ! WRITE(*,*) ' ' 00264 ! WRITE(*,*) 'Shifts and drift velocity deviations:' 00265 ! DO I=1,NPLAN 00266 ! WRITE(*,102) I,DEL(I),DVD(I) 00267 ! END DO 00268 00269 00270 WRITE(*,*) ' ' 00271 WRITE(*,*) ' ' 00272 WRITE(*,*) ncount,' tracks generated with ',nthits,' hits.' 00273 WRITE(*,*) nrecds,' records written.' 00274 WRITE(*,*) ' ' 00275 101 FORMAT(a) 00276 ! 102 FORMAT(I6,2F10.5) 00277 103 FORMAT(i8,f10.5) 00278 END SUBROUTINE mptest ! gener 00279 00283 00284 SUBROUTINE genlin(ip) 00285 USE mptest1 00286 00287 IMPLICIT NONE 00288 REAL :: gr 00289 REAL :: gran 00290 REAL :: uran 00291 REAL :: x 00292 REAL :: ybias 00293 REAL :: ydvds 00294 REAL :: ylin 00295 REAL :: ymeas 00296 REAL :: ywire 00297 INTEGER :: i 00298 INTEGER :: nwire 00299 00300 INTEGER, INTENT(IN) :: ip 00301 00302 ! ... 00303 ynull=0.5*heit+0.1*heit*(uran()-0.5) ! uniform vertex 00304 slope=(uran()-0.5)*heit/(FLOAT(nplan-1)*disx) 00305 IF(ip /= 0) THEN 00306 WRITE(*,*) ' ' 00307 ! WRITE(*,*) 'YNULL=',YNULL,' SLOPE=',SLOPE 00308 END IF 00309 nhits=0 00310 DO i=1,nplan 00311 x=detx+FLOAT(i-1)*disx ! +0.5*THCK 00312 IF(uran() < eff(i)) THEN 00313 ylin =ynull+slope*x ! true y value 00314 ybias =ylin-del(i) ! biased value 00315 nwire=INT(1.0+ybias/4.0) ! wire number 00316 IF(nwire <= 0.OR.nwire > 25) EXIT ! check wire number 00317 nhits=nhits+1 ! track hits the plane 00318 xhits(nhits)=x 00319 ihits(nhits)=i 00320 gr=gran() 00321 ymeas=sgm(i)*gr 00322 ydvds=0.0 00323 yhits(nhits)=ybias+ymeas+ydvds ! measured 00324 ywire=FLOAT(nwire)*4.0-2.0 00325 ydrft(nhits)=ybias-ywire ! signed drift length 00326 ydvds=ydrft(nhits)*dvd(i) 00327 yhits(nhits)=ybias+ymeas-ydvds ! measured 00328 sigma(nhits)=sgm(i) 00329 IF(ip /= 0) THEN 00330 ! WRITE(*,101) NHITS,I,X,YLIN,YBIAS,YMEAS, 00331 ! + SGM(I),YHITS(NHITS),GR,DEL(I) 00332 END IF 00333 END IF 00334 END DO 00335 ! 101 FORMAT(2I3,F5.0,7F8.4) 00336 END SUBROUTINE genlin
1.8.0
