/[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.3 by mocchiut, Fri Mar 30 11:17:17 2007 UTC revision 1.5 by mocchiut, Mon Sep 29 12:40:43 2008 UTC
# Line 25  CC      PARAMETER (calwidth=24.2) Line 25  CC      PARAMETER (calwidth=24.2)
25  C        C      
26        INTEGER i,j,it        INTEGER i,j,it
27        INTEGER ifin,finpar,finpar1,plent,plentx,plenty        INTEGER ifin,finpar,finpar1,plent,plentx,plenty
28        INTEGER xncnt(22),yncnt(22),Nfitx,Nfity        INTEGER xncnt(NPLAV),yncnt(NPLAV),Nfitx,Nfity
29        INTEGER xncnt2(5,22),yncnt2(5,22)        INTEGER xncnt2(5,NPLAV),yncnt2(5,NPLAV)
30  C  C
31        REAL xecnt2(5,22),xwght2(5,22),xcorr2(5,22)        REAL xecnt2(5,NPLAV),xwght2(5,NPLAV),xcorr2(5,NPLAV)
32        REAL yecnt2(5,22),ywght2(5,22),ycorr2(5,22)        REAL yecnt2(5,NPLAV),ywght2(5,NPLAV),ycorr2(5,NPLAV)
33        REAL ax2(4),bx2(4),eax2(4),ebx2(4)        REAL ax2(4),bx2(4),eax2(4),ebx2(4)
34        REAL ay2(4),by2(4),eay2(4),eby2(4)        REAL ay2(4),by2(4),eay2(4),eby2(4)
35        REAL xEmax3(22),yEmax3(22),xmax(22),ymax(22),zx(22),zy(22)        REAL xEmax3(NPLAV),yEmax3(NPLAV),xmax(NPLAV)
36        REAL xecnt(22),xwght(22),xcorr(22)        REAL ymax(NPLAV),zx(NPLAV),zy(NPLAV)
37        REAL yecnt(22),ywght(22),ycorr(22)        REAL xecnt(NPLAV),xwght(NPLAV),xcorr(NPLAV)
38        REAL ax,bx,eax,ebx,chi2x,qxp,posixmd        REAL yecnt(NPLAV),ywght(NPLAV),ycorr(NPLAV)
39        REAL ay,by,eay,eby,chi2y,qyp,posiymd        REAL ax,bx,eax,ebx,chi2x
40        REAL thxm,thym,tmisd,dthxm,dthym,dtmisd,pmisd,pmid        REAL ay,by,eay,eby,chi2y
41        REAL enet,parze,parz(2,22),ffla,zetamx,zetamy,fflaco,parzen2        REAL thxm,thym,tmisd,pmisd,pmid
42          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 56  C Line 57  C
57        parze = 0.        parze = 0.
58        finpar = 1        finpar = 1
59        finpar1 = 1        finpar1 = 1
60        CALL VZERO(parz,2*22)        CALL VZERO(parz,2*NPLAV)
61        DO i = 1,22        DO i = 1,NPLA
62           DO j = 1,96           DO j = 1,96
63              parz(1,i) = parz(1,i) + DEXY(1,i,j) ! sum up the energy in each x-plane              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              parz(2,i) = parz(2,i) + DEXY(2,i,j) ! sum up the energy in each y-plane
# Line 81  C Line 82  C
82  C        C      
83  c     Find the center of the 3 strips with maximum total energy in a plane  c     Find the center of the 3 strips with maximum total energy in a plane
84  C  C
85        CALL VZERO(xemax3,22)        CALL VZERO(xemax3,NPLAV)
86        CALL VZERO(xncnt,22)        CALL VZERO(xncnt,NPLAV)
87        CALL VZERO(yemax3,22)        CALL VZERO(yemax3,NPLAV)
88        CALL VZERO(yncnt,22)        CALL VZERO(yncnt,NPLAV)
89        DO i = 1,22        DO i = 1,NPLA
90           DO j = 1,94           DO j = 1,94
91    
92              IF ((DEXY(1,i,j)+              IF ((DEXY(1,i,j)+
# Line 111  C Line 112  C
112  c     Calculate the position of the center strip and the center of  c     Calculate the position of the center strip and the center of
113  c     energy (CoE) of 4 surrounding strips  c     energy (CoE) of 4 surrounding strips
114  C  C
115        DO i=1,22        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
146           CALL LEASTSQR(22,zx,xecnt,xwght,ax,bx,eax,ebx,chi2x,Nfitx) !<----- LINEAR FIT           CALL LEASTSQR(NPLA,zx,xecnt,xwght,ax,bx,eax,ebx,chi2x,Nfitx) !<----- LINEAR FIT
147           CALL LEASTSQR(22,zy,yecnt,ywght,ay,by,eay,eby,chi2y,Nfity)           CALL LEASTSQR(NPLA,zy,yecnt,ywght,ay,by,eay,eby,chi2y,Nfity)
148           ax2(it)=ax           ax2(it)=ax
149           ay2(it)=ay           ay2(it)=ay
150           eax2(it)=eax           eax2(it)=eax
# Line 150  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 170  C Line 177  C
177  c     Calculate the position of the strip closest to the fitted line and  c     Calculate the position of the strip closest to the fitted line and
178  c     the CoE of 4 surrounding strips.  c     the CoE of 4 surrounding strips.
179  C  C
180           DO i=1,22           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 204  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
251           IF ((ABS(ax+bx*zbotx).LE.10.1).AND.           IF ((ABS(ax+bx*zbotx).LE.10.1).AND.
252       &        (ABS(ay+by*zboty).LE.10.1)) THEN       &        (ABS(ay+by*zboty).LE.10.1)) THEN
253              ifin = finpar + 3              ifin = finpar + 3
254              IF (ifin.GT.22) ifin = 22              IF (ifin.GT.NPLA) ifin = NPLA
255           ELSE           ELSE
256              ifin = finpar              ifin = finpar
257           ENDIF           ENDIF
# Line 261  C Line 275  C
275  C  C
276  c     Angle correction  c     Angle correction
277  C  C
278           DO i=1,22                   DO i=1,NPLA        
279  C            C          
280  c           x view              c           x view            
281  C  C
# Line 410  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 426  c Find the plane closest to a certain z Line 441  c Find the plane closest to a certain z
441  C---------------------------------------------------------------------------  C---------------------------------------------------------------------------
442    
443        IMPLICIT NONE        IMPLICIT NONE
444    C
445          INCLUDE 'INTEST.TXT'
446    C
447        REAL pos, npos, pdiff, xy        REAL pos, npos, pdiff, xy
448        INTEGER view, plane, nplane        INTEGER view, plane, nplane
449    
450        pdiff=1000.        pdiff=1000.
451        plane=1        plane=1
452        DO view=1,2        DO view=1,2
453           DO nplane=1,22           DO nplane=1,NPLA
454              CALL STRIP2POS(view,nplane,1,xy,npos)              CALL STRIP2POS(view,nplane,1,xy,npos)
455              IF (ABS(pos-npos).LT.pdiff) THEN              IF (ABS(pos-npos).LT.pdiff) THEN
456                 pdiff=ABS(pos-npos)                 pdiff=ABS(pos-npos)
# Line 480  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
504  c Arguments  c Arguments
505  c ---------  c ---------
506  c view:   x=1 and y=2  c view:   x=1 and y=2
507  c nplane: plane number 1-22  c nplane: plane number 1-NPLA
508  c nstrip: strip number 1-96  c nstrip: strip number 1-96
509  c  c
510  c Return values  c Return values
# Line 497  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.3  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23