/[PAMELA software]/DarthVader/CalorimeterLevel2/src/selftrig.for
ViewVC logotype

Diff of /DarthVader/CalorimeterLevel2/src/selftrig.for

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by mocchiut, Fri Jul 20 08:24:55 2007 UTC revision 1.5 by mocchiut, Mon Sep 29 12:40:43 2008 UTC
# Line 36  C Line 36  C
36        REAL ymax(NPLAV),zx(NPLAV),zy(NPLAV)        REAL ymax(NPLAV),zx(NPLAV),zy(NPLAV)
37        REAL xecnt(NPLAV),xwght(NPLAV),xcorr(NPLAV)        REAL xecnt(NPLAV),xwght(NPLAV),xcorr(NPLAV)
38        REAL yecnt(NPLAV),ywght(NPLAV),ycorr(NPLAV)        REAL yecnt(NPLAV),ywght(NPLAV),ycorr(NPLAV)
39        REAL ax,bx,eax,ebx,chi2x,qxp,posixmd        REAL ax,bx,eax,ebx,chi2x
40        REAL ay,by,eay,eby,chi2y,qyp,posiymd        REAL ay,by,eay,eby,chi2y
41        REAL thxm,thym,tmisd,dthxm,dthym,dtmisd,pmisd,pmid        REAL thxm,thym,tmisd,pmisd,pmid
42        REAL enet,parze,parz(2,NPLAV),ffla,zetamx,zetamy,fflaco,parzen2        REAL enet,parze,parz(2,NPLAV),ffla,fflaco,parzen2
43        REAL parzen3,funcor2        REAL parzen3,funcor2
44        REAL CORRANG,ENCORR        REAL CORRANG,ENCORR
45  C  C
# Line 114  c     energy (CoE) of 4 surrounding stri Line 114  c     energy (CoE) of 4 surrounding stri
114  C  C
115        DO i=1,NPLA        DO i=1,NPLA
116           CALL STRIP2POS(1,i,xncnt(i),xmax(i),zx(i))           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))           CALL STRIP2POS(2,i,yncnt(i),ymax(i),zy(i))
120           CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i))           CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i))
121           CALL ENCENT2(2,yncnt(i),i,4,yecnt(i),ywght(i))           CALL ENCENT2(2,yncnt(i),i,4,yecnt(i),ywght(i))
122           xecnt2(1,i)=xecnt(i)           xecnt2(1,i)=xecnt(i)
123           xwght2(1,i)=xwght(i)           xwght2(1,i)=xwght(i)
124           xncnt2(1,i)=xncnt(i)           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)           yecnt2(1,i)=yecnt(i)
128           ywght2(1,i)=ywght(i)           ywght2(1,i)=ywght(i)
129           yncnt2(1,i)=yncnt(i)           yncnt2(1,i)=yncnt(i)
130        ENDDO        ENDDO
131  C        C      
132        dtmisd=1.0  c      dtmisd=1.0
133        dthxm=1.0  c      dthxm=1.0
134        dthym=1.0  c      dthym=1.0
135        tmisd=0        tmisd=0.
136        thxm=0        thxm=0.
137        thym=0        thym=0.
138  C  C
139  c     Iterative procedure starts here  c     Iterative procedure starts here
140  C  C
141        DO it=1,4        DO it=1,4
142    cc      DO it=1,1
143  C  C
144  c     Fit the CoEs with a linear function  c     Fit the CoEs with a linear function
145  C  C
# Line 151  C Line 156  C
156  C  C
157  c     Calculate theta and phi  c     Calculate theta and phi
158  C  C
159           dtmisd=ABS(tmisd-ATAN(SQRT((bx*bx)+(by*by))))  c         dtmisd=ABS(tmisd-ATAN(SQRT((bx*bx)+(by*by))))
160           dthxm=ABS(thxm-ABS(ATAN(bx)))  c         dthxm=ABS(thxm-ABS(ATAN(bx)))
161           dthym=ABS(thym-ABS(ATAN(by)))  c         dthym=ABS(thym-ABS(ATAN(by)))
162           tmisd=ATAN(SQRT((bx*bx)+(by*by)))           tmisd=ATAN(SQRT((bx*bx)+(by*by)))
163           thxm=ABS(ATAN(bx))           thxm=ABS(ATAN(bx))
164           thym=ABS(ATAN(by))           thym=ABS(ATAN(by))
165    c
166           IF (bx.EQ.0..AND.by.GT.0.) pmisd = 90.*(PI/180.)           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.)           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.)           IF (by.EQ.0..AND.bx.GE.0.) pmisd = 0.*(PI/180.)
# Line 172  c     Calculate the position of the stri Line 178  c     Calculate the position of the stri
178  c     the CoE of 4 surrounding strips.  c     the CoE of 4 surrounding strips.
179  C  C
180           DO i=1,NPLA           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))              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))              CALL POS2STRIP(2,i,ay+by*zy(i),yncnt(i))
186              IF (xncnt(i).GE.1.AND.xncnt(i).LE.96) THEN              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))                 CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i))
189              ELSE              ELSE
190                 xecnt(i)=-99.                 xecnt(i)=-99.
191                 xwght(i)=1.0e5                 xwght(i)=1.0e5
192              ENDIF              ENDIF
193    c      print *,'Cx plane ',i,' strip ',xncnt(i),' ecnt ',
194    c     +        xecnt(i)
195  C              C            
196              IF (yncnt(i).GE.1.AND.yncnt(i).LE.96) THEN              IF (yncnt(i).GE.1.AND.yncnt(i).LE.96) THEN
197                 CALL ENCENT2(2,yncnt(i),i,4,yecnt(i),ywght(i))                 CALL ENCENT2(2,yncnt(i),i,4,yecnt(i),ywght(i))
# Line 205  C Line 217  C
217           IF (by.LT.0.) CALL POS2PLANE(-(calwidth/2)/by-ay/by,plenty)           IF (by.LT.0.) CALL POS2PLANE(-(calwidth/2)/by-ay/by,plenty)
218           IF (by.EQ.0.) plenty=1           IF (by.EQ.0.) plenty=1
219           plent=MAX(plentx,plenty)                   plent=MAX(plentx,plenty)        
220    c         print *,' plentx ',plentx,' plenty ',plenty,' plent ',plent
221  C  C
222  c     Calculate the projection in the bottom plane  c     Calculate the projection in the bottom plane
223  C  C
224           qxp=ax+bx*ztopx  c         qxp=ax+bx*ztopx
225           qyp=ay+by*ztopy  c         qyp=ay+by*ztopy
226           zetamx=ztopx-zbotx  c         zetamx=ztopx-zbotx
227           zetamy=ztopy-zboty  c         zetamy=ztopy-zboty
228  C          cC        
229          IF (qxp.GT.calwidth/2) THEN  c        IF (qxp.GT.calwidth/2) THEN
230             zetamx = zetamx-(qxp-calwidth/2)/bx  c           zetamx = zetamx-(qxp-calwidth/2)/bx
231             qxp = calwidth/2  c           qxp = calwidth/2
232          ENDIF  c        ENDIF
233          IF (qxp.LT.-calwidth/2) THEN  c        IF (qxp.LT.-calwidth/2) THEN
234             zetamx = zetamx-(qxp+calwidth/2)/bx  c           zetamx = zetamx-(qxp+calwidth/2)/bx
235             qxp = -calwidth/2  c           qxp = -calwidth/2
236          ENDIF  c        ENDIF
237          IF (qyp.GT.calwidth/2) THEN  c        IF (qyp.GT.calwidth/2) THEN
238             zetamy = zetamy-(qyp-calwidth/2)/by  c           zetamy = zetamy-(qyp-calwidth/2)/by
239             qyp = calwidth/2  c           qyp = calwidth/2
240          ENDIF  c        ENDIF
241          IF (qyp.LT.-calwidth/2) THEN  c        IF (qyp.LT.-calwidth/2) THEN
242             zetamy = zetamy-(qyp+calwidth/2)/by  c           zetamy = zetamy-(qyp+calwidth/2)/by
243             qyp = -calwidth/2  c           qyp = -calwidth/2
244          ENDIF  c        ENDIF
245  C  cC
246          posixmd = qxp - zetamx*TAN(tmisd)*COS(pmisd)  c        posixmd = qxp - zetamx*TAN(tmisd)*COS(pmisd)
247          posiymd = qyp - zetamy*TAN(tmisd)*SIN(pmisd)  c        posiymd = qyp - zetamy*TAN(tmisd)*SIN(pmisd)
248  C  C
249  c     Energy correction  c     Energy correction
250  C  C
# Line 411  C       Line 424  C      
424           ecnt=mx           ecnt=mx
425           IF (nrec.EQ.1) weight=0.244/sqrt(12*s2)           IF (nrec.EQ.1) weight=0.244/sqrt(12*s2)
426           IF (nrec.NE.1) weight=sqrt(mx2-mx**2)/sqrt(s2)           IF (nrec.NE.1) weight=sqrt(mx2-mx**2)/sqrt(s2)
427    ccc         weight = 1.
428        ELSE        ELSE
429           ecnt=-99.0           ecnt=-99.0
430           weight=100000.0           weight=100000.0
# Line 483  C       Line 497  C      
497        END        END
498    
499  C---------------------------------------------------------------------  C---------------------------------------------------------------------
500        SUBROUTINE STRIP2POS(view,nplane,nstrip,xy,z)        SUBROUTINE STRIP2POS(view,nplane,snstrip,xy,z)
501  c  c
502  c Calculates the x (or y) and z position of a strip in PAMELA coordinates  c Calculates the x (or y) and z position of a strip in PAMELA coordinates
503  c  c
# Line 500  c z:  z position of the plane in PMAELA Line 514  c z:  z position of the plane in PMAELA
514  C---------------------------------------------------------------------  C---------------------------------------------------------------------
515    
516        IMPLICIT NONE        IMPLICIT NONE
517    C
518        INTEGER view,isy,iseven,nplane,nstrip,nwaf,wstr,npl        INCLUDE 'INTEST.TXT'
519    C
520          INTEGER view,isy,iseven,nplane,snstrip,nwaf,wstr,npl
521        REAL wpos,pos,x,y,z,xy        REAL wpos,pos,x,y,z,xy
522                
523  c     Calculate the x- or y-position in a plane  c     Calculate the x- or y-position in a plane
524    
525        nwaf=int((nstrip-1)/32) !what wafer (0-2)        nwaf=int((snstrip-1)/32) !what wafer (0-2)
526        wstr=mod(nstrip-1,32) !what strip in the wafer (0-31)        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        wpos=0.096+0.122+wstr*0.244 !the position of the strip center in the wafer
528  c      pos=wpos+nwaf*8.0005-12.0005 !the position in the plane (oigo shifted to center of plane)        pos=wpos+nwaf*8.05 !-12.05 !the position in the plane (oigo shifted to center of plane)
       pos=wpos+nwaf*8.0505-12.0005 !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  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        isy=view-1 !isy=0 for x and 1 for y
       iseven=mod(nplane,2)! iseven=0 for odd and 1 for even planes  
533        npl=nplane-1        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        IF (isy.EQ.0) THEN
537    C
538    C X views
539    C
540           IF (iseven.EQ.0) THEN           IF (iseven.EQ.0) THEN
541  c            print *,'CALCULATING X-ODD'              x = pos + 0.05 - XALIG/10.
542              x = pos+0.05 !x-odd              y = pos + 0.1 - YALIG/10.
543              y = pos-0.2              z = (ZALIG/10.) - 0.581 -0.809*npl -0.2*npl/2.
544              z = -26.762-0.809*npl-0.2*(npl-1)/2  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           ELSE
547  c            print *,'CALCULATING X-EVEN'              x = pos - 0.05 - XALIG/10.
548              x = pos-0.05 !x-even              y = pos - 0.1 - YALIG/10.
549              y = pos+0.0              z = (ZALIG/10.) - 0.581 -0.809*npl -0.2*(npl -1)/2.
550              z = -26.762-0.809*npl-0.2*npl/2  c-26.762-0.809*npl-0.2*npl/2
551    c            print *,'CALCULATING X-ODD ',x,y,z !x1,x3,x5
552          ENDIF          ENDIF
553          xy=x          xy=x
554        ELSE        ELSE
555           IF (iseven.EQ.0) THEN           IF (iseven.EQ.0) THEN
556  c            print *,'CALCULATING Y-ODD'              x = pos + 0.1 - XALIG/10.
557              x = pos-0.1 ! y-odd              y = pos + 0.05 - YALIG/10.
558              y = pos-0.15              z = (ZALIG/10.) -0.809*npl-0.2*npl/2.
559              z = -26.181-0.809*npl-0.2*(npl-1)/2              c-26.181-0.809*npl-0.2*(npl-1)/2            
560    c            print *,'CALCULATING Y-EVEN ',x,y,z
561           ELSE           ELSE
562  c            print *,'CALCULATING Y-EVEN'              x = pos - 0.1 - XALIG/10.
563              x =  pos+0.1 ! y-even              y = pos - 0.05 -YALIG/10.
564              y =  pos-0.05              z = (ZALIG/10.) -0.809*npl-0.2*(npl-1)/2.
565              z = -26.181-0.809*npl-0.2*npl/2  c            print *,'CALCULATING Y-ODD ',x,y,z
566    c-26.181-0.809*npl-0.2*npl/2
567          ENDIF          ENDIF
568          xy=y          xy=y
569        ENDIF        ENDIF
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23