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