/[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.6 by mocchiut, Mon Sep 28 17:06:55 2009 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 392  C       Line 405  C      
405        nrec=0        nrec=0
406        DO i=nsmax-width,nsmax+width        DO i=nsmax-width,nsmax+width
407  C  C
408           IF(i.GT.0.AND.i.LT.97.AND.DEXY(view,nplane,i).GT.0.0) THEN           IF(i.GT.0.AND.i.LT.97) THEN
409  C                          IF(DEXY(view,nplane,i).GT.0.0) THEN
410              nrec=nrec+1  C    
411  C                 nrec=nrec+1
412              CALL STRIP2POS(view,nplane,i,xypos,zpos)  C    
413  C                             CALL STRIP2POS(view,nplane,i,xypos,zpos)
414              s1=s2  C    
415              s2=s2+DEXY(view,nplane,i)                 s1=s2
416              mx=(mx*s1+xypos*DEXY(view,nplane,i))/s2                 s2=s2+DEXY(view,nplane,i)
417              mx2=(mx2*s1+(xypos**2)*DEXY(view,nplane,i))/s2                 mx=(mx*s1+xypos*DEXY(view,nplane,i))/s2
418  C                 mx2=(mx2*s1+(xypos**2)*DEXY(view,nplane,i))/s2
419    C    
420                ENDIF
421           ENDIF           ENDIF
422  C          C        
423        ENDDO        ENDDO
# Line 411  C       Line 426  C      
426           ecnt=mx           ecnt=mx
427           IF (nrec.EQ.1) weight=0.244/sqrt(12*s2)           IF (nrec.EQ.1) weight=0.244/sqrt(12*s2)
428           IF (nrec.NE.1) weight=sqrt(mx2-mx**2)/sqrt(s2)           IF (nrec.NE.1) weight=sqrt(mx2-mx**2)/sqrt(s2)
429    ccc         weight = 1.
430        ELSE        ELSE
431           ecnt=-99.0           ecnt=-99.0
432           weight=100000.0           weight=100000.0
# Line 483  C       Line 499  C      
499        END        END
500    
501  C---------------------------------------------------------------------  C---------------------------------------------------------------------
502        SUBROUTINE STRIP2POS(view,nplane,nstrip,xy,z)        SUBROUTINE STRIP2POS(view,nplane,snstrip,xy,z)
503  c  c
504  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
505  c  c
# Line 500  c z:  z position of the plane in PMAELA Line 516  c z:  z position of the plane in PMAELA
516  C---------------------------------------------------------------------  C---------------------------------------------------------------------
517    
518        IMPLICIT NONE        IMPLICIT NONE
519    C
520        INTEGER view,isy,iseven,nplane,nstrip,nwaf,wstr,npl        INCLUDE 'INTEST.TXT'
521    C
522          INTEGER view,isy,iseven,nplane,snstrip,nwaf,wstr,npl
523        REAL wpos,pos,x,y,z,xy        REAL wpos,pos,x,y,z,xy
524                
525  c     Calculate the x- or y-position in a plane  c     Calculate the x- or y-position in a plane
526    
527        nwaf=int((nstrip-1)/32) !what wafer (0-2)        nwaf=int((snstrip-1)/32) !what wafer (0-2)
528        wstr=mod(nstrip-1,32) !what strip in the wafer (0-31)        wstr=mod(snstrip-1,32) !what strip in the wafer (0-31)
529        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
530  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)  
531            
532  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
533    
534        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  
535        npl=nplane-1        npl=nplane-1
536          iseven=mod(npl,2)! iseven=0 for even and 1 for odd planes
537    c
538        IF (isy.EQ.0) THEN        IF (isy.EQ.0) THEN
539    C
540    C X views
541    C
542           IF (iseven.EQ.0) THEN           IF (iseven.EQ.0) THEN
543  c            print *,'CALCULATING X-ODD'              x = pos + 0.05 - XALIG/10.
544              x = pos+0.05 !x-odd              y = pos + 0.1 - YALIG/10.
545              y = pos-0.2              z = (ZALIG/10.) - 0.581 -0.809*npl -0.2*npl/2.
546              z = -26.762-0.809*npl-0.2*(npl-1)/2  c-26.762-0.809*npl-0.2*(npl-1)/2
547    c            print *,'CALCULATING X-EVEN ',x,y,z ! x0,x2,x4,...
548           ELSE           ELSE
549  c            print *,'CALCULATING X-EVEN'              x = pos - 0.05 - XALIG/10.
550              x = pos-0.05 !x-even              y = pos - 0.1 - YALIG/10.
551              y = pos+0.0              z = (ZALIG/10.) - 0.581 -0.809*npl -0.2*(npl -1)/2.
552              z = -26.762-0.809*npl-0.2*npl/2  c-26.762-0.809*npl-0.2*npl/2
553    c            print *,'CALCULATING X-ODD ',x,y,z !x1,x3,x5
554          ENDIF          ENDIF
555          xy=x          xy=x
556        ELSE        ELSE
557           IF (iseven.EQ.0) THEN           IF (iseven.EQ.0) THEN
558  c            print *,'CALCULATING Y-ODD'              x = pos + 0.1 - XALIG/10.
559              x = pos-0.1 ! y-odd              y = pos + 0.05 - YALIG/10.
560              y = pos-0.15              z = (ZALIG/10.) -0.809*npl-0.2*npl/2.
561              z = -26.181-0.809*npl-0.2*(npl-1)/2              c-26.181-0.809*npl-0.2*(npl-1)/2            
562    c            print *,'CALCULATING Y-EVEN ',x,y,z
563           ELSE           ELSE
564  c            print *,'CALCULATING Y-EVEN'              x = pos - 0.1 - XALIG/10.
565              x =  pos+0.1 ! y-even              y = pos - 0.05 -YALIG/10.
566              y =  pos-0.05              z = (ZALIG/10.) -0.809*npl-0.2*(npl-1)/2.
567              z = -26.181-0.809*npl-0.2*npl/2  c            print *,'CALCULATING Y-ODD ',x,y,z
568    c-26.181-0.809*npl-0.2*npl/2
569          ENDIF          ENDIF
570          xy=y          xy=y
571        ENDIF        ENDIF
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23