| 1 |
C |
| 2 |
C |
| 3 |
C NOTE: THIS ROUTINE DOES NOT SEEMS TO WORK CORRECTLY, HAS TO BE CHECKED CAREFULLY |
| 4 |
C |
| 5 |
C |
| 6 |
C--------------------------------------------------------------------- |
| 7 |
SUBROUTINE SELFTRIG |
| 8 |
C--------------------------------------------------------------------- |
| 9 |
C |
| 10 |
IMPLICIT NONE |
| 11 |
C |
| 12 |
INCLUDE 'INTEST.TXT' |
| 13 |
C |
| 14 |
REAL PI,calwidth,ztopx,ztopy,zbotx,zboty,MIP2GEV |
| 15 |
real wcorr |
| 16 |
parameter (wcorr=1.0) |
| 17 |
PARAMETER (MIP2GEV=0.0001059994) |
| 18 |
PARAMETER (PI=3.14159265358979324) |
| 19 |
CC PARAMETER (calwidth=24.2) |
| 20 |
PARAMETER (calwidth=24.1) |
| 21 |
PARAMETER (ztopx=-26.18) |
| 22 |
PARAMETER (ztopy=-26.76) |
| 23 |
PARAMETER (zbotx=-45.17) |
| 24 |
PARAMETER (zboty=-45.75) |
| 25 |
C |
| 26 |
INTEGER i,j,it |
| 27 |
INTEGER ifin,finpar,finpar1,plent,plentx,plenty |
| 28 |
INTEGER xncnt(NPLAV),yncnt(NPLAV),Nfitx,Nfity |
| 29 |
INTEGER xncnt2(5,NPLAV),yncnt2(5,NPLAV) |
| 30 |
C |
| 31 |
REAL xecnt2(5,NPLAV),xwght2(5,NPLAV),xcorr2(5,NPLAV) |
| 32 |
REAL yecnt2(5,NPLAV),ywght2(5,NPLAV),ycorr2(5,NPLAV) |
| 33 |
REAL ax2(4),bx2(4),eax2(4),ebx2(4) |
| 34 |
REAL ay2(4),by2(4),eay2(4),eby2(4) |
| 35 |
REAL xEmax3(NPLAV),yEmax3(NPLAV),xmax(NPLAV) |
| 36 |
REAL ymax(NPLAV),zx(NPLAV),zy(NPLAV) |
| 37 |
REAL xecnt(NPLAV),xwght(NPLAV),xcorr(NPLAV) |
| 38 |
REAL yecnt(NPLAV),ywght(NPLAV),ycorr(NPLAV) |
| 39 |
REAL ax,bx,eax,ebx,chi2x |
| 40 |
REAL ay,by,eay,eby,chi2y |
| 41 |
REAL thxm,thym,tmisd,pmisd,pmid |
| 42 |
REAL enet,parze,parz(2,NPLAV),ffla,fflaco,parzen2 |
| 43 |
REAL parzen3,funcor2 |
| 44 |
REAL CORRANG,ENCORR |
| 45 |
C |
| 46 |
COMMON / slftrig / tmisd,ax,bx,eax,ebx,chi2x,Nfitx,ay,by,eay,eby, |
| 47 |
& chi2y,Nfity,parzen3 |
| 48 |
SAVE / slftrig / |
| 49 |
C |
| 50 |
COMMON / debug / xncnt2,yncnt2,xecnt2,xwght2,xcorr2,yecnt2,ywght2, |
| 51 |
& ycorr2,ax2,bx2,eax2,ebx2,ay2,by2,eay2,eby2,zx,zy |
| 52 |
SAVE / debug / |
| 53 |
C |
| 54 |
c Energy calculation |
| 55 |
C |
| 56 |
enet = 0. |
| 57 |
parze = 0. |
| 58 |
finpar = 1 |
| 59 |
finpar1 = 1 |
| 60 |
CALL VZERO(parz,2*NPLAV) |
| 61 |
DO i = 1,NPLA |
| 62 |
DO j = 1,96 |
| 63 |
parz(1,i) = parz(1,i) + DEXY(1,i,j) ! sum up the energy in each x-plane |
| 64 |
parz(2,i) = parz(2,i) + DEXY(2,i,j) ! sum up the energy in each y-plane |
| 65 |
enet = enet + DEXY(1,i,j) + DEXY(2,i,j) ! sum up the total energy |
| 66 |
ENDDO |
| 67 |
IF (parz(1,i).GE.parze) THEN ! find plane with max energy |
| 68 |
parze = parz(1,i) ! the energy |
| 69 |
finpar = i ! the plane number |
| 70 |
finpar1 = 1 ! set x or y indicator |
| 71 |
ENDIF |
| 72 |
IF (parz(2,i).GE.parze) THEN |
| 73 |
parze = parz(2,i) |
| 74 |
finpar = i |
| 75 |
finpar1 = 2 |
| 76 |
ENDIF |
| 77 |
ENDDO |
| 78 |
ffla = FLOAT(finpar) |
| 79 |
IF (finpar1.EQ.1) THEN |
| 80 |
ffla = ffla - 1. |
| 81 |
ENDIF |
| 82 |
C |
| 83 |
c Find the center of the 3 strips with maximum total energy in a plane |
| 84 |
C |
| 85 |
CALL VZERO(xemax3,NPLAV) |
| 86 |
CALL VZERO(xncnt,NPLAV) |
| 87 |
CALL VZERO(yemax3,NPLAV) |
| 88 |
CALL VZERO(yncnt,NPLAV) |
| 89 |
DO i = 1,NPLA |
| 90 |
DO j = 1,94 |
| 91 |
|
| 92 |
IF ((DEXY(1,i,j)+ |
| 93 |
& DEXY(1,i,j+1)+ |
| 94 |
& DEXY(1,i,j+2)).GT.xemax3(i)) THEN |
| 95 |
xemax3(i)=DEXY(1,i,j)+ |
| 96 |
& DEXY(1,i,j+1)+ |
| 97 |
& DEXY(1,i,j+2) |
| 98 |
xncnt(i)=j+1 |
| 99 |
ENDIF |
| 100 |
IF ((DEXY(2,i,j)+ |
| 101 |
& DEXY(2,i,j+1)+ |
| 102 |
& DEXY(2,i,j+2)).GT.yemax3(i)) THEN |
| 103 |
yemax3(i)=DEXY(2,i,j)+ |
| 104 |
& DEXY(2,i,j+1)+ |
| 105 |
& DEXY(2,i,j+2) |
| 106 |
yncnt(i)=j+1 |
| 107 |
ENDIF |
| 108 |
|
| 109 |
ENDDO |
| 110 |
ENDDO |
| 111 |
C |
| 112 |
c Calculate the position of the center strip and the center of |
| 113 |
c energy (CoE) of 4 surrounding strips |
| 114 |
C |
| 115 |
DO i=1,NPLA |
| 116 |
CALL STRIP2POS(1,i,xncnt(i),xmax(i),zx(i)) |
| 117 |
c print *,'x plane ',i,' strip ',xncnt(i),' pos ', |
| 118 |
c + xmax(i),' z ',zx(i) |
| 119 |
CALL STRIP2POS(2,i,yncnt(i),ymax(i),zy(i)) |
| 120 |
CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i)) |
| 121 |
CALL ENCENT2(2,yncnt(i),i,4,yecnt(i),ywght(i)) |
| 122 |
xecnt2(1,i)=xecnt(i) |
| 123 |
xwght2(1,i)=xwght(i) |
| 124 |
xncnt2(1,i)=xncnt(i) |
| 125 |
c print *,'x plane ',i,' strip ',xncnt2(1,i),' ecnt ', |
| 126 |
c + xecnt2(1,i),' xncnt2 ',xncnt2(1,i) |
| 127 |
yecnt2(1,i)=yecnt(i) |
| 128 |
ywght2(1,i)=ywght(i) |
| 129 |
yncnt2(1,i)=yncnt(i) |
| 130 |
ENDDO |
| 131 |
C |
| 132 |
c dtmisd=1.0 |
| 133 |
c dthxm=1.0 |
| 134 |
c dthym=1.0 |
| 135 |
tmisd=0. |
| 136 |
thxm=0. |
| 137 |
thym=0. |
| 138 |
C |
| 139 |
c Iterative procedure starts here |
| 140 |
C |
| 141 |
DO it=1,4 |
| 142 |
cc DO it=1,1 |
| 143 |
C |
| 144 |
c Fit the CoEs with a linear function |
| 145 |
C |
| 146 |
CALL LEASTSQR(NPLA,zx,xecnt,xwght,ax,bx,eax,ebx,chi2x,Nfitx) !<----- LINEAR FIT |
| 147 |
CALL LEASTSQR(NPLA,zy,yecnt,ywght,ay,by,eay,eby,chi2y,Nfity) |
| 148 |
ax2(it)=ax |
| 149 |
ay2(it)=ay |
| 150 |
eax2(it)=eax |
| 151 |
eay2(it)=eay |
| 152 |
bx2(it)=bx |
| 153 |
by2(it)=by |
| 154 |
ebx2(it)=ebx |
| 155 |
eby2(it)=eby |
| 156 |
C |
| 157 |
c Calculate theta and phi |
| 158 |
C |
| 159 |
c dtmisd=ABS(tmisd-ATAN(SQRT((bx*bx)+(by*by)))) |
| 160 |
c dthxm=ABS(thxm-ABS(ATAN(bx))) |
| 161 |
c dthym=ABS(thym-ABS(ATAN(by))) |
| 162 |
tmisd=ATAN(SQRT((bx*bx)+(by*by))) |
| 163 |
thxm=ABS(ATAN(bx)) |
| 164 |
thym=ABS(ATAN(by)) |
| 165 |
c |
| 166 |
IF (bx.EQ.0..AND.by.GT.0.) pmisd = 90.*(PI/180.) |
| 167 |
IF (bx.EQ.0..AND.by.LT.0.) pmisd = 270.*(PI/180.) |
| 168 |
IF (by.EQ.0..AND.bx.GE.0.) pmisd = 0.*(PI/180.) |
| 169 |
IF (by.EQ.0..AND.bx.LT.0.) pmisd = 180.*(PI/180.) |
| 170 |
IF (by.NE.0..AND.bx.NE.0.) THEN |
| 171 |
pmid = ATAN(by/bx) |
| 172 |
IF (by.LT.0..AND.bx.GT.0.) pmisd = pmid + 360.*(PI/180.) |
| 173 |
IF (bx.LT.0.) pmisd = pmid + 180.*(PI/180.) |
| 174 |
IF (by.GT.0..AND.bx.GT.0.) pmisd = pmid |
| 175 |
ENDIF |
| 176 |
C |
| 177 |
c Calculate the position of the strip closest to the fitted line and |
| 178 |
c the CoE of 4 surrounding strips. |
| 179 |
C |
| 180 |
DO i=1,NPLA |
| 181 |
c CALL STRIP2POS(1,i,xncnt(i),xmax(i),zx(i)) |
| 182 |
CALL POS2STRIP(1,i,ax+bx*zx(i),xncnt(i)) |
| 183 |
c print *,'Bx plane ',i,' strip ',ax+bx*zx(i),' pos ', |
| 184 |
c + xncnt(i) |
| 185 |
CALL POS2STRIP(2,i,ay+by*zy(i),yncnt(i)) |
| 186 |
IF (xncnt(i).GE.1.AND.xncnt(i).LE.96) THEN |
| 187 |
c CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i)) |
| 188 |
CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i)) |
| 189 |
ELSE |
| 190 |
xecnt(i)=-99. |
| 191 |
xwght(i)=1.0e5 |
| 192 |
ENDIF |
| 193 |
c print *,'Cx plane ',i,' strip ',xncnt(i),' ecnt ', |
| 194 |
c + xecnt(i) |
| 195 |
C |
| 196 |
IF (yncnt(i).GE.1.AND.yncnt(i).LE.96) THEN |
| 197 |
CALL ENCENT2(2,yncnt(i),i,4,yecnt(i),ywght(i)) |
| 198 |
ELSE |
| 199 |
yecnt(i)=-99. |
| 200 |
ywght(i)=1.0e5 |
| 201 |
ENDIF |
| 202 |
xncnt2(it+1,i)=xncnt(i) |
| 203 |
yncnt2(it+1,i)=yncnt(i) |
| 204 |
xecnt2(it+1,i)=xecnt(i) |
| 205 |
xwght2(it+1,i)=xwght(i) |
| 206 |
yecnt2(it+1,i)=yecnt(i) |
| 207 |
ywght2(it+1,i)=ywght(i) |
| 208 |
ENDDO |
| 209 |
C |
| 210 |
c Calculate at which plane the particle has entered the calorimeter |
| 211 |
c according to the fit |
| 212 |
C |
| 213 |
IF (bx.GT.0.) CALL POS2PLANE((calwidth/2)/bx-ax/bx,plentx) |
| 214 |
IF (bx.LT.0.) CALL POS2PLANE(-(calwidth/2)/bx-ax/bx,plentx) |
| 215 |
IF (bx.EQ.0.) plentx=1 |
| 216 |
IF (by.GT.0.) CALL POS2PLANE((calwidth/2)/by-ay/by,plenty) |
| 217 |
IF (by.LT.0.) CALL POS2PLANE(-(calwidth/2)/by-ay/by,plenty) |
| 218 |
IF (by.EQ.0.) plenty=1 |
| 219 |
plent=MAX(plentx,plenty) |
| 220 |
c print *,' plentx ',plentx,' plenty ',plenty,' plent ',plent |
| 221 |
C |
| 222 |
c Calculate the projection in the bottom plane |
| 223 |
C |
| 224 |
c qxp=ax+bx*ztopx |
| 225 |
c qyp=ay+by*ztopy |
| 226 |
c zetamx=ztopx-zbotx |
| 227 |
c zetamy=ztopy-zboty |
| 228 |
cC |
| 229 |
c IF (qxp.GT.calwidth/2) THEN |
| 230 |
c zetamx = zetamx-(qxp-calwidth/2)/bx |
| 231 |
c qxp = calwidth/2 |
| 232 |
c ENDIF |
| 233 |
c IF (qxp.LT.-calwidth/2) THEN |
| 234 |
c zetamx = zetamx-(qxp+calwidth/2)/bx |
| 235 |
c qxp = -calwidth/2 |
| 236 |
c ENDIF |
| 237 |
c IF (qyp.GT.calwidth/2) THEN |
| 238 |
c zetamy = zetamy-(qyp-calwidth/2)/by |
| 239 |
c qyp = calwidth/2 |
| 240 |
c ENDIF |
| 241 |
c IF (qyp.LT.-calwidth/2) THEN |
| 242 |
c zetamy = zetamy-(qyp+calwidth/2)/by |
| 243 |
c qyp = -calwidth/2 |
| 244 |
c ENDIF |
| 245 |
cC |
| 246 |
c posixmd = qxp - zetamx*TAN(tmisd)*COS(pmisd) |
| 247 |
c posiymd = qyp - zetamy*TAN(tmisd)*SIN(pmisd) |
| 248 |
C |
| 249 |
c Energy correction |
| 250 |
C |
| 251 |
IF ((ABS(ax+bx*zbotx).LE.10.1).AND. |
| 252 |
& (ABS(ay+by*zboty).LE.10.1)) THEN |
| 253 |
ifin = finpar + 3 |
| 254 |
IF (ifin.GT.NPLA) ifin = NPLA |
| 255 |
ELSE |
| 256 |
ifin = finpar |
| 257 |
ENDIF |
| 258 |
|
| 259 |
parzen2 = 0.0 |
| 260 |
DO i=1,ifin |
| 261 |
parzen2=parzen2+parz(1,i)+parz(2,i) !Sum up energy until 'ifin' |
| 262 |
ENDDO |
| 263 |
|
| 264 |
fflaco = ifin-plent |
| 265 |
funcor2=parzen2/ENCORR(fflaco/COS(tmisd)) |
| 266 |
|
| 267 |
IF ((ABS(ax+bx*zbotx).LE.10.1).AND. |
| 268 |
& (ABS(ay+by*zboty).LE.10.1)) THEN |
| 269 |
parzen3 = funcor2*(1.-.01775-funcor2*(.2096E-4)+ |
| 270 |
& funcor2*funcor2*funcor2*(.2865E-9)) |
| 271 |
ELSE |
| 272 |
parzen3 = funcor2*(1.-.124+funcor2*(.24E-3)+ |
| 273 |
& funcor2*funcor2*funcor2*(.25E-9))/.7 |
| 274 |
ENDIF |
| 275 |
C |
| 276 |
c Angle correction |
| 277 |
C |
| 278 |
DO i=1,NPLA |
| 279 |
C |
| 280 |
c x view |
| 281 |
C |
| 282 |
IF (xecnt(i).GT.-97.) THEN |
| 283 |
IF (parzen3.GE.250.) THEN |
| 284 |
xcorr(i) = wcorr*CORRANG(FLOAT(i)/ffla,(thxm*180./PI)+ |
| 285 |
& (parzen3-250)*1.764705E-2) |
| 286 |
ELSE |
| 287 |
xcorr(i) = wcorr*CORRANG(FLOAT(i)/ffla,thxm*180./PI) |
| 288 |
ENDIF |
| 289 |
ELSE |
| 290 |
xcorr(i) = 0.0 |
| 291 |
ENDIF |
| 292 |
xecnt(i) = xecnt(i) + xcorr(i) |
| 293 |
C |
| 294 |
c y view |
| 295 |
C |
| 296 |
IF (yecnt(i).GT.-97.) THEN |
| 297 |
IF (parzen3.GE.250.) THEN |
| 298 |
ycorr(i) = wcorr*CORRANG(FLOAT(i)/ffla,thym*180./PI+ |
| 299 |
& (parzen3-250)*1.764705E-2) |
| 300 |
ELSE |
| 301 |
ycorr(i) = wcorr*CORRANG(FLOAT(i)/ffla,thym*180./PI) |
| 302 |
ENDIF |
| 303 |
ELSE |
| 304 |
ycorr(i) = 0.0 |
| 305 |
ENDIF |
| 306 |
yecnt(i) = yecnt(i) + ycorr(i) |
| 307 |
C |
| 308 |
xcorr2(it,i)=xcorr(i) |
| 309 |
ycorr2(it,i)=ycorr(i) |
| 310 |
C |
| 311 |
ENDDO |
| 312 |
C |
| 313 |
ENDDO |
| 314 |
C |
| 315 |
RETURN |
| 316 |
C |
| 317 |
END |
| 318 |
|
| 319 |
C |
| 320 |
C----------------------------------------------------------------------- |
| 321 |
SUBROUTINE ENCENT(view,nsmax,nplane,nit,ecnt,weight) |
| 322 |
c |
| 323 |
c Calculates the center of energy in a cluster |
| 324 |
c |
| 325 |
C----------------------------------------------------------------------- |
| 326 |
C |
| 327 |
IMPLICIT NONE |
| 328 |
C |
| 329 |
INCLUDE 'INTEST.TXT' |
| 330 |
C |
| 331 |
REAL ecnt,weight,esumw,esum,strip1,strip2,xypos,zpos,xymean, |
| 332 |
& xystd |
| 333 |
INTEGER nsmax,st0,st1,st2,st3,nplane,nit,view,npos,p |
| 334 |
C |
| 335 |
PARAMETER (strip1=17.,strip2=9.) |
| 336 |
C |
| 337 |
st0 = NINT(strip1) ! =17 cluster width for iteration 1 |
| 338 |
st1 = NINT(strip1-strip2) ! =8 |
| 339 |
st2 = NINT((strip1-1.)/2. + 1.) ! =9 cluster width for iteration 2 |
| 340 |
st3 = NINT(FLOAT(st1)/2.) ! =4 |
| 341 |
C |
| 342 |
IF (nsmax.GT.0.) THEN |
| 343 |
esumw = 0. |
| 344 |
esum = 0. |
| 345 |
xymean = 0. |
| 346 |
xystd = 0. |
| 347 |
DO P = 1, (st0-(nit-1)*st1) ! loop to 17(9) for iteration 1(2) |
| 348 |
C |
| 349 |
c Calculate strip number |
| 350 |
C |
| 351 |
IF (nsmax.LT.(st2-st3*(nit-1))) THEN ! if nsmax < 9(5) for iteration 1(2) |
| 352 |
npos = P |
| 353 |
IF (npos.GT.(st0-ABS(nsmax-st2-st3*(nit-1)))) |
| 354 |
& GO TO 7100 ! quit loop after the 8th(4th) strip to the right of nsmax for iteration 1(2) |
| 355 |
ELSE |
| 356 |
npos = nsmax+P-st2+st3*(nit-1) ! npos=nsmax-8,-7, ... +7,+8(nsmax-4,-3, ... +3,+4) for iteration 1(2) |
| 357 |
IF (npos.GT.96) |
| 358 |
& GO TO 7100 ! quit loop after the 96th strip |
| 359 |
ENDIF |
| 360 |
C |
| 361 |
c Get the position for npos |
| 362 |
C |
| 363 |
CALL STRIP2POS(view,nplane,npos,xypos,zpos) |
| 364 |
C |
| 365 |
c Sum up energies |
| 366 |
C |
| 367 |
esumw = esumw + xypos*DEXY(view,nplane,npos) |
| 368 |
esum = esum + DEXY(view,nplane,npos) |
| 369 |
ENDDO |
| 370 |
C |
| 371 |
7100 CONTINUE |
| 372 |
C |
| 373 |
c Calculate CoE and wieghts |
| 374 |
C |
| 375 |
IF (esum.GT.0.) THEN |
| 376 |
ecnt = esumw/esum |
| 377 |
weight = ecnt/(esum**.79) |
| 378 |
ELSE |
| 379 |
ecnt = -98. |
| 380 |
weight = 1E5 |
| 381 |
ENDIF |
| 382 |
ELSE |
| 383 |
ecnt = -97. |
| 384 |
weight = 1E5 |
| 385 |
ENDIF |
| 386 |
C |
| 387 |
RETURN |
| 388 |
C |
| 389 |
END |
| 390 |
C |
| 391 |
C----------------------------------------------------------------------- |
| 392 |
SUBROUTINE ENCENT2(view,nsmax,nplane,width,ecnt,weight) |
| 393 |
C----------------------------------------------------------------------- |
| 394 |
C |
| 395 |
IMPLICIT NONE |
| 396 |
C |
| 397 |
INCLUDE 'INTEST.TXT' |
| 398 |
C |
| 399 |
INTEGER i,nplane,view,nsmax,width,nrec |
| 400 |
REAL mx,mx2,s1,s2,xypos,ecnt,weight,zpos |
| 401 |
C |
| 402 |
mx=0.0 |
| 403 |
mx2=0.0 |
| 404 |
s2=0.0 |
| 405 |
nrec=0 |
| 406 |
DO i=nsmax-width,nsmax+width |
| 407 |
C |
| 408 |
IF(i.GT.0.AND.i.LT.97.AND.DEXY(view,nplane,i).GT.0.0) THEN |
| 409 |
C |
| 410 |
nrec=nrec+1 |
| 411 |
C |
| 412 |
CALL STRIP2POS(view,nplane,i,xypos,zpos) |
| 413 |
C |
| 414 |
s1=s2 |
| 415 |
s2=s2+DEXY(view,nplane,i) |
| 416 |
mx=(mx*s1+xypos*DEXY(view,nplane,i))/s2 |
| 417 |
mx2=(mx2*s1+(xypos**2)*DEXY(view,nplane,i))/s2 |
| 418 |
C |
| 419 |
ENDIF |
| 420 |
C |
| 421 |
ENDDO |
| 422 |
C |
| 423 |
IF (nrec.GT.0) THEN |
| 424 |
ecnt=mx |
| 425 |
IF (nrec.EQ.1) weight=0.244/sqrt(12*s2) |
| 426 |
IF (nrec.NE.1) weight=sqrt(mx2-mx**2)/sqrt(s2) |
| 427 |
ccc weight = 1. |
| 428 |
ELSE |
| 429 |
ecnt=-99.0 |
| 430 |
weight=100000.0 |
| 431 |
ENDIF |
| 432 |
C |
| 433 |
RETURN |
| 434 |
C |
| 435 |
END |
| 436 |
|
| 437 |
C--------------------------------------------------------------------------- |
| 438 |
SUBROUTINE POS2PLANE(pos,plane) |
| 439 |
c |
| 440 |
c Find the plane closest to a certain z position |
| 441 |
C--------------------------------------------------------------------------- |
| 442 |
|
| 443 |
IMPLICIT NONE |
| 444 |
C |
| 445 |
INCLUDE 'INTEST.TXT' |
| 446 |
C |
| 447 |
REAL pos, npos, pdiff, xy |
| 448 |
INTEGER view, plane, nplane |
| 449 |
|
| 450 |
pdiff=1000. |
| 451 |
plane=1 |
| 452 |
DO view=1,2 |
| 453 |
DO nplane=1,NPLA |
| 454 |
CALL STRIP2POS(view,nplane,1,xy,npos) |
| 455 |
IF (ABS(pos-npos).LT.pdiff) THEN |
| 456 |
pdiff=ABS(pos-npos) |
| 457 |
plane=nplane |
| 458 |
ENDIF |
| 459 |
ENDDO |
| 460 |
ENDDO |
| 461 |
C |
| 462 |
RETURN |
| 463 |
C |
| 464 |
END |
| 465 |
|
| 466 |
C--------------------------------------------------------------------------- |
| 467 |
SUBROUTINE POS2STRIP(view,nplane,pos,strip) |
| 468 |
c |
| 469 |
c Find the strip closest to a certain x or y position |
| 470 |
C--------------------------------------------------------------------------- |
| 471 |
C |
| 472 |
IMPLICIT NONE |
| 473 |
C |
| 474 |
REAL pos, npos, minpos, maxpos, pdiff, z |
| 475 |
INTEGER view, nplane, strip, i |
| 476 |
C |
| 477 |
pdiff=1000. |
| 478 |
strip=-1 |
| 479 |
CALL STRIP2POS(view,nplane,1,minpos,z) |
| 480 |
CALL STRIP2POS(view,nplane,96,maxpos,z) |
| 481 |
IF (pos.LT.minpos) THEN |
| 482 |
strip=0 |
| 483 |
ELSEIF (pos.GT.maxpos) THEN |
| 484 |
strip=97 |
| 485 |
ELSE |
| 486 |
DO i=1,96 |
| 487 |
CALL STRIP2POS(view,nplane,i,npos,z) |
| 488 |
IF (ABS(pos-npos).LT.pdiff) THEN |
| 489 |
pdiff=ABS(pos-npos) |
| 490 |
strip=i |
| 491 |
ENDIF |
| 492 |
ENDDO |
| 493 |
ENDIF |
| 494 |
C |
| 495 |
RETURN |
| 496 |
C |
| 497 |
END |
| 498 |
|
| 499 |
C--------------------------------------------------------------------- |
| 500 |
SUBROUTINE STRIP2POS(view,nplane,snstrip,xy,z) |
| 501 |
c |
| 502 |
c Calculates the x (or y) and z position of a strip in PAMELA coordinates |
| 503 |
c |
| 504 |
c Arguments |
| 505 |
c --------- |
| 506 |
c view: x=1 and y=2 |
| 507 |
c nplane: plane number 1-NPLA |
| 508 |
c nstrip: strip number 1-96 |
| 509 |
c |
| 510 |
c Return values |
| 511 |
c ------------- |
| 512 |
c xy: x or y position of the strip-center in PAMELA coordinates |
| 513 |
c z: z position of the plane in PMAELA coordinates |
| 514 |
C--------------------------------------------------------------------- |
| 515 |
|
| 516 |
IMPLICIT NONE |
| 517 |
C |
| 518 |
INCLUDE 'INTEST.TXT' |
| 519 |
C |
| 520 |
INTEGER view,isy,iseven,nplane,snstrip,nwaf,wstr,npl |
| 521 |
REAL wpos,pos,x,y,z,xy |
| 522 |
|
| 523 |
c Calculate the x- or y-position in a plane |
| 524 |
|
| 525 |
nwaf=int((snstrip-1)/32) !what wafer (0-2) |
| 526 |
wstr=mod(snstrip-1,32) !what strip in the wafer (0-31) |
| 527 |
wpos=0.096+0.122+wstr*0.244 !the position of the strip center in the wafer |
| 528 |
pos=wpos+nwaf*8.05 !-12.05 !the position in the plane (oigo shifted to center of plane) |
| 529 |
|
| 530 |
c Calculate z position and add x-y-offset depending on the plane number |
| 531 |
|
| 532 |
isy=view-1 !isy=0 for x and 1 for y |
| 533 |
npl=nplane-1 |
| 534 |
iseven=mod(npl,2)! iseven=0 for even and 1 for odd planes |
| 535 |
c |
| 536 |
IF (isy.EQ.0) THEN |
| 537 |
C |
| 538 |
C X views |
| 539 |
C |
| 540 |
IF (iseven.EQ.0) THEN |
| 541 |
x = pos + 0.05 - XALIG/10. |
| 542 |
y = pos + 0.1 - YALIG/10. |
| 543 |
z = (ZALIG/10.) - 0.581 -0.809*npl -0.2*npl/2. |
| 544 |
c-26.762-0.809*npl-0.2*(npl-1)/2 |
| 545 |
c print *,'CALCULATING X-EVEN ',x,y,z ! x0,x2,x4,... |
| 546 |
ELSE |
| 547 |
x = pos - 0.05 - XALIG/10. |
| 548 |
y = pos - 0.1 - YALIG/10. |
| 549 |
z = (ZALIG/10.) - 0.581 -0.809*npl -0.2*(npl -1)/2. |
| 550 |
c-26.762-0.809*npl-0.2*npl/2 |
| 551 |
c print *,'CALCULATING X-ODD ',x,y,z !x1,x3,x5 |
| 552 |
ENDIF |
| 553 |
xy=x |
| 554 |
ELSE |
| 555 |
IF (iseven.EQ.0) THEN |
| 556 |
x = pos + 0.1 - XALIG/10. |
| 557 |
y = pos + 0.05 - YALIG/10. |
| 558 |
z = (ZALIG/10.) -0.809*npl-0.2*npl/2. |
| 559 |
c-26.181-0.809*npl-0.2*(npl-1)/2 |
| 560 |
c print *,'CALCULATING Y-EVEN ',x,y,z |
| 561 |
ELSE |
| 562 |
x = pos - 0.1 - XALIG/10. |
| 563 |
y = pos - 0.05 -YALIG/10. |
| 564 |
z = (ZALIG/10.) -0.809*npl-0.2*(npl-1)/2. |
| 565 |
c print *,'CALCULATING Y-ODD ',x,y,z |
| 566 |
c-26.181-0.809*npl-0.2*npl/2 |
| 567 |
ENDIF |
| 568 |
xy=y |
| 569 |
ENDIF |
| 570 |
C |
| 571 |
RETURN |
| 572 |
END |
| 573 |
|
| 574 |
C--------------------------------------------------------------------------- |
| 575 |
SUBROUTINE LEASTSQR(N,x,y,sy,a,b,ea,eb,chi2,Nfit) |
| 576 |
c |
| 577 |
c Linear least square fit (method is described in "Taylor" chapter 8) |
| 578 |
c |
| 579 |
c |
| 580 |
c |
| 581 |
C--------------------------------------------------------------------------- |
| 582 |
C |
| 583 |
INTEGER N,Nfit |
| 584 |
REAL x(N),y(N),sy(N),w(N),Swx2,Swx,Swy,Swxy,Sw,a,b,ea,eb,chi2 |
| 585 |
C |
| 586 |
Swx2=0 |
| 587 |
Swx=0 |
| 588 |
Swy=0 |
| 589 |
Swxy=0 |
| 590 |
Sw=0 |
| 591 |
Nfit=0 |
| 592 |
DO i=1,N |
| 593 |
IF (y(i).GT.-97.) THEN |
| 594 |
Nfit=Nfit+1 |
| 595 |
w(i)=1/(sy(i)**2) ! the weight |
| 596 |
Swx2=Swx2+w(i)*(x(i)**2)! some sums |
| 597 |
Swx=Swx+w(i)*x(i) |
| 598 |
Swy=Swy+w(i)*y(i) |
| 599 |
Swxy=Swxy+w(i)*x(i)*y(i) |
| 600 |
Sw=Sw+w(i) |
| 601 |
c print *,'FIT LOOP:',sy(i),w(i),x(i),y(i),Swx2,Swx,Swy,Swxy,Sw |
| 602 |
ENDIF |
| 603 |
ENDDO |
| 604 |
|
| 605 |
delta=Sw*Swx2-(Swx**2) |
| 606 |
a=(Swx2*Swy-Swx*Swxy)/delta ! offset |
| 607 |
b=(Sw*Swxy-Swx*Swy)/delta ! slope |
| 608 |
ea=sqrt(Swx2/delta) ! offset standard deviation |
| 609 |
eb=sqrt(Sw/delta) ! slope standard deviation |
| 610 |
|
| 611 |
chi2=0.0 |
| 612 |
DO i=1,N |
| 613 |
IF (y(i).GT.-97.) THEN |
| 614 |
chi2=chi2+((a+b*x(i))-y(i))**2/sy(i)**2 |
| 615 |
ENDIF |
| 616 |
ENDDO |
| 617 |
|
| 618 |
c print *,'FIT RESULT:',a,b,chi2,Nfit |
| 619 |
RETURN |
| 620 |
END |
| 621 |
|
| 622 |
C--------------------------------------------------------------------------- |
| 623 |
REAL FUNCTION ENCORR(X) |
| 624 |
C--------------------------------------------------------------------------- |
| 625 |
|
| 626 |
REAL FFUN |
| 627 |
PARAMETER (CALIB=0.0001059994) |
| 628 |
PARAMETER (P1=.8993962E-2) |
| 629 |
PARAMETER (P2=-.449717) |
| 630 |
PARAMETER (P3=10.90797) |
| 631 |
PARAMETER (P4=-.6768349E-2) |
| 632 |
|
| 633 |
|
| 634 |
FFUN = P1 + P4 * ATAN((X - P3) * P2) |
| 635 |
|
| 636 |
IF (FFUN.EQ.0.) THEN |
| 637 |
FFUN = 1.E-10 |
| 638 |
ENDIF |
| 639 |
|
| 640 |
ENCORR = FFUN/CALIB |
| 641 |
|
| 642 |
RETURN |
| 643 |
END |
| 644 |
|
| 645 |
C--------------------------------------------------------------------------- |
| 646 |
REAL FUNCTION CORRANG(X,ANGX) |
| 647 |
C--------------------------------------------------------------------------- |
| 648 |
|
| 649 |
REAL A,B,X,ANGX |
| 650 |
|
| 651 |
A = -0.017695 + ANGX * 0.0016963 |
| 652 |
B = -0.049583 + ANGX * 0.0050639 |
| 653 |
CORRANG = 0. |
| 654 |
IF (X.LE..9) CORRANG = A * (X - .9) |
| 655 |
IF (X.GE..9.AND.X.LE.1.1) CORRANG = 0. |
| 656 |
IF (X.GE.1.1) CORRANG = B * (X - 1.1) |
| 657 |
IF (ANGX.LT.10.) CORRANG = 0. |
| 658 |
|
| 659 |
RETURN |
| 660 |
END |
| 661 |
|
Parent Directory
| 
Revision Log