/[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.2 by mocchiut, Tue Jan 16 11:06:05 2007 UTC revision 1.4 by mocchiut, Fri Jul 20 08:24:55 2007 UTC
# Line 1  Line 1 
1  C  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---------------------------------------------------------------------  C---------------------------------------------------------------------
7        SUBROUTINE SELFTRIG        SUBROUTINE SELFTRIG
8  C---------------------------------------------------------------------  C---------------------------------------------------------------------
# Line 21  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 yecnt(NPLAV),ywght(NPLAV),ycorr(NPLAV)
39        REAL ax,bx,eax,ebx,chi2x,qxp,posixmd        REAL ax,bx,eax,ebx,chi2x,qxp,posixmd
40        REAL ay,by,eay,eby,chi2y,qyp,posiymd        REAL ay,by,eay,eby,chi2y,qyp,posiymd
41        REAL thxm,thym,tmisd,dthxm,dthym,dtmisd,pmisd,pmid        REAL thxm,thym,tmisd,dthxm,dthym,dtmisd,pmisd,pmid
42        REAL enet,parze,parz(2,22),ffla,zetamx,zetamy,fflaco,parzen2        REAL enet,parze,parz(2,NPLAV),ffla,zetamx,zetamy,fflaco,parzen2
43        REAL parzen3,funcor2        REAL parzen3,funcor2
44        REAL CORRANG,ENCORR        REAL CORRANG,ENCORR
45  C  C
# Line 52  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) + ESTRIP(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) + ESTRIP(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
65              enet = enet + ESTRIP(1,i,j) + ESTRIP(2,i,j) ! sum up the total energy              enet = enet + DEXY(1,i,j) + DEXY(2,i,j) ! sum up the total energy
66           ENDDO           ENDDO
67           IF (parz(1,i).GE.parze) THEN ! find plane with max energy           IF (parz(1,i).GE.parze) THEN ! find plane with max energy
68              parze = parz(1,i)   ! the energy              parze = parz(1,i)   ! the energy
# Line 77  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 ((ESTRIP(1,i,j)+              IF ((DEXY(1,i,j)+
93       &           ESTRIP(1,i,j+1)+       &           DEXY(1,i,j+1)+
94       &           ESTRIP(1,i,j+2)).GT.xemax3(i)) THEN       &           DEXY(1,i,j+2)).GT.xemax3(i)) THEN
95                 xemax3(i)=ESTRIP(1,i,j)+                 xemax3(i)=DEXY(1,i,j)+
96       &                   ESTRIP(1,i,j+1)+       &                   DEXY(1,i,j+1)+
97       &                   ESTRIP(1,i,j+2)       &                   DEXY(1,i,j+2)
98                 xncnt(i)=j+1                 xncnt(i)=j+1
99              ENDIF              ENDIF
100              IF ((ESTRIP(2,i,j)+              IF ((DEXY(2,i,j)+
101       &           ESTRIP(2,i,j+1)+       &           DEXY(2,i,j+1)+
102       &           ESTRIP(2,i,j+2)).GT.yemax3(i)) THEN       &           DEXY(2,i,j+2)).GT.yemax3(i)) THEN
103                 yemax3(i)=ESTRIP(2,i,j)+                 yemax3(i)=DEXY(2,i,j)+
104       &                   ESTRIP(2,i,j+1)+       &                   DEXY(2,i,j+1)+
105       &                   ESTRIP(2,i,j+2)       &                   DEXY(2,i,j+2)
106                 yncnt(i)=j+1                 yncnt(i)=j+1
107              ENDIF              ENDIF
108    
# Line 107  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           CALL STRIP2POS(2,i,yncnt(i),ymax(i),zy(i))           CALL STRIP2POS(2,i,yncnt(i),ymax(i),zy(i))
118           CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i))           CALL ENCENT2(1,xncnt(i),i,4,xecnt(i),xwght(i))
# Line 133  C Line 138  C
138  C  C
139  c     Fit the CoEs with a linear function  c     Fit the CoEs with a linear function
140  C  C
141           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
142           CALL LEASTSQR(22,zy,yecnt,ywght,ay,by,eay,eby,chi2y,Nfity)           CALL LEASTSQR(NPLA,zy,yecnt,ywght,ay,by,eay,eby,chi2y,Nfity)
143           ax2(it)=ax           ax2(it)=ax
144           ay2(it)=ay           ay2(it)=ay
145           eax2(it)=eax           eax2(it)=eax
# Line 166  C Line 171  C
171  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
172  c     the CoE of 4 surrounding strips.  c     the CoE of 4 surrounding strips.
173  C  C
174           DO i=1,22           DO i=1,NPLA
175              CALL POS2STRIP(1,i,ax+bx*zx(i),xncnt(i))              CALL POS2STRIP(1,i,ax+bx*zx(i),xncnt(i))
176              CALL POS2STRIP(2,i,ay+by*zy(i),yncnt(i))              CALL POS2STRIP(2,i,ay+by*zy(i),yncnt(i))
177              IF (xncnt(i).GE.1.AND.xncnt(i).LE.96) THEN              IF (xncnt(i).GE.1.AND.xncnt(i).LE.96) THEN
# Line 233  C Line 238  C
238           IF ((ABS(ax+bx*zbotx).LE.10.1).AND.           IF ((ABS(ax+bx*zbotx).LE.10.1).AND.
239       &        (ABS(ay+by*zboty).LE.10.1)) THEN       &        (ABS(ay+by*zboty).LE.10.1)) THEN
240              ifin = finpar + 3              ifin = finpar + 3
241              IF (ifin.GT.22) ifin = 22              IF (ifin.GT.NPLA) ifin = NPLA
242           ELSE           ELSE
243              ifin = finpar              ifin = finpar
244           ENDIF           ENDIF
# Line 257  C Line 262  C
262  C  C
263  c     Angle correction  c     Angle correction
264  C  C
265           DO i=1,22                   DO i=1,NPLA        
266  C            C          
267  c           x view              c           x view            
268  C  C
# Line 346  C Line 351  C
351  C  C
352  c     Sum up energies  c     Sum up energies
353  C  C
354              esumw = esumw + xypos*ESTRIP(view,nplane,npos)              esumw = esumw + xypos*DEXY(view,nplane,npos)
355              esum = esum + ESTRIP(view,nplane,npos)              esum = esum + DEXY(view,nplane,npos)
356           ENDDO           ENDDO
357  C  C
358   7100    CONTINUE   7100    CONTINUE
# Line 387  C       Line 392  C      
392        nrec=0        nrec=0
393        DO i=nsmax-width,nsmax+width        DO i=nsmax-width,nsmax+width
394  C  C
395           IF(i.GT.0.AND.i.LT.97.AND.ESTRIP(view,nplane,i).GT.0.0) THEN           IF(i.GT.0.AND.i.LT.97.AND.DEXY(view,nplane,i).GT.0.0) THEN
396  C              C            
397              nrec=nrec+1              nrec=nrec+1
398  C  C
399              CALL STRIP2POS(view,nplane,i,xypos,zpos)              CALL STRIP2POS(view,nplane,i,xypos,zpos)
400  C              C            
401              s1=s2              s1=s2
402              s2=s2+ESTRIP(view,nplane,i)              s2=s2+DEXY(view,nplane,i)
403              mx=(mx*s1+xypos*ESTRIP(view,nplane,i))/s2              mx=(mx*s1+xypos*DEXY(view,nplane,i))/s2
404              mx2=(mx2*s1+(xypos**2)*ESTRIP(view,nplane,i))/s2              mx2=(mx2*s1+(xypos**2)*DEXY(view,nplane,i))/s2
405  C  C
406           ENDIF           ENDIF
407  C          C        
# Line 422  c Find the plane closest to a certain z Line 427  c Find the plane closest to a certain z
427  C---------------------------------------------------------------------------  C---------------------------------------------------------------------------
428    
429        IMPLICIT NONE        IMPLICIT NONE
430    C
431          INCLUDE 'INTEST.TXT'
432    C
433        REAL pos, npos, pdiff, xy        REAL pos, npos, pdiff, xy
434        INTEGER view, plane, nplane        INTEGER view, plane, nplane
435    
436        pdiff=1000.        pdiff=1000.
437        plane=1        plane=1
438        DO view=1,2        DO view=1,2
439           DO nplane=1,22           DO nplane=1,NPLA
440              CALL STRIP2POS(view,nplane,1,xy,npos)              CALL STRIP2POS(view,nplane,1,xy,npos)
441              IF (ABS(pos-npos).LT.pdiff) THEN              IF (ABS(pos-npos).LT.pdiff) THEN
442                 pdiff=ABS(pos-npos)                 pdiff=ABS(pos-npos)
# Line 483  c Line 490  c
490  c Arguments  c Arguments
491  c ---------  c ---------
492  c view:   x=1 and y=2  c view:   x=1 and y=2
493  c nplane: plane number 1-22  c nplane: plane number 1-NPLA
494  c nstrip: strip number 1-96  c nstrip: strip number 1-96
495  c  c
496  c Return values  c Return values
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23