/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/functionspfa.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/functionspfa.f

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

revision 1.21 by pam-fi, Fri Aug 31 14:56:52 2007 UTC revision 1.29 by mocchiut, Fri Jan 17 10:14:39 2014 UTC
# Line 8  Line 8 
8  *  *
9  *     integer function npfastrips(ic,angle)  *     integer function npfastrips(ic,angle)
10  *  *
11    *     -----------------------------------------------------------------
12    *     p.f.a.
13    *     -----------------------------------------------------------------
14  *     real function pfaeta(ic,angle)  *     real function pfaeta(ic,angle)
15  *     real function pfaetal(ic,angle)  *     real function pfaetal(ic,angle)
16  *     real function pfaeta2(ic,angle)  *     real function pfaeta2(ic,angle)
# Line 15  Line 18 
18  *     real function pfaeta4(ic,angle)  *     real function pfaeta4(ic,angle)
19  *     real function cog(ncog,ic)  *     real function cog(ncog,ic)
20  *  *
21    *     -----------------------------------------------------------------
22    *     risoluzione spaziale media, stimata dalla simulazione (samuele)
23    *     -----------------------------------------------------------------
24    *     FUNCTION risxeta2(angle)
25    *     FUNCTION risxeta3(angle)
26    *     FUNCTION risxeta4(angle)
27    *     FUNCTION risyeta2(angle)
28    *     FUNCTION risy_cog(angle)
29    *     FUNCTION risx_cog(angle)
30    *     real function riseta(iview,angle)
31    *     -----------------------------------------------------------------
32    *     fattore moltiplicativo per tenere conto della dipendenza della
33    *     risoluzione dal rumore delle strip
34    *     -----------------------------------------------------------------
35  *     real function fbad_cog(ncog,ic)  *     real function fbad_cog(ncog,ic)
36  *     real function fbad_eta(ic,angle)  *     real function fbad_eta(ic,angle)
37  *  *
38  *     real function riseta(iview,angle)  *     -----------------------------------------------------------------
39  *     FUNCTION risxeta2(x)  *     NUOVO APPROCCIO PER LA STIMA DELLA RISOLUZIONE
40  *     FUNCTION risxeta3(x)  *     -----------------------------------------------------------------
41  *     FUNCTION risxeta4(x)  *     real function riscogtheor(ncog,ic)
42  *     FUNCTION risyeta2(x)  *     real function risetatheor(ncog,ic,angle)
 *     FUNCTION risy_cog(x)  
 *     FUNCTION risx_cog(x)  
43  *  *
44    *     -----------------------------------------------------------------
45    *     correzione landi
46    *     -----------------------------------------------------------------
47  *     real function pfacorr(ic,angle)  *     real function pfacorr(ic,angle)
48  *  *
49  *     real function effectiveangle(ang,iview,bbb)  *     real function effectiveangle(ang,iview,bbb)
# Line 39  Line 57 
57        subroutine idtoc(ipfa,cpfa)        subroutine idtoc(ipfa,cpfa)
58                
59        integer ipfa        integer ipfa
60        character*10 cpfa  c      character*10 cpfa
61          character*4 cpfa ! EM GCC4.7
62    
63        CPFA='COG4'        CPFA='COG4'
64        if(ipfa.eq.0)CPFA='ETA'        if(ipfa.eq.0)CPFA='ETA'
# Line 56  Line 75 
75        end        end
76  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
77        real function effectiveangle(ang,iview,bbb)        real function effectiveangle(ang,iview,bbb)
78          
79        include 'commontracker.f'        include 'commontracker.f'
80          real tgtemp
81    
82        effectiveangle = 0.        effectiveangle = 0.
83    
# Line 70  c     here bbb is the y component of the Line 90  c     here bbb is the y component of the
90           by   = bbb           by   = bbb
91           if(iview.eq.12) angx = -1. * ang           if(iview.eq.12) angx = -1. * ang
92           if(iview.eq.12) by   = -1. * bbb           if(iview.eq.12) by   = -1. * bbb
93           tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001  cc         tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!!
94             tgtemp  = tan(angx*acos(-1.)/180.) + REAL(pmuH_h*by*0.00001)  ! EM GCC4.7 pmuH_h is double precision but all the others are real...
95    
96        elseif(mod(iview,2).eq.1)then        elseif(mod(iview,2).eq.1)then
97  c     =================================================  c     =================================================
# Line 79  c     ================================== Line 100  c     ==================================
100  c     here bbb is the x component of the m.filed  c     here bbb is the x component of the m.filed
101           angy = ang           angy = ang
102           bx   = bbb           bx   = bbb
103           tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001                   tgtemp  = tan(angy*acos(-1.)/180.)+real(pmuH_e*bx*0.00001) ! EM GCC4.7 pmuH_h is double precision but all the others are real...
104    
105        endif              endif      
106        effectiveangle = 180.*atan(tgtemp)/acos(-1.)        effectiveangle = 180.*atan(tgtemp)/acos(-1.)
# Line 101  c     ================================== Line 122  c     ==================================
122  c     here bbb is the y component of the m.field  c     here bbb is the y component of the m.field
123           by   = bbb           by   = bbb
124           if(iview.eq.12) by = -1. * bbb           if(iview.eq.12) by = -1. * bbb
125           fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX           fieldcorr = -1. * 0.5*REAL(pmuH_h*by*0.00001*SiDimZ/pitchX)! EM GCC4.7 pmuH_h is double precision but all the others are real...
126    
127        elseif(mod(iview,2).eq.1)then        elseif(mod(iview,2).eq.1)then
128  c     =================================================  c     =================================================
# Line 109  c     Y view Line 130  c     Y view
130  c     =================================================          c     =================================================        
131  c     here bbb is the x component of the m.filed  c     here bbb is the x component of the m.filed
132           bx   = bbb           bx   = bbb
133           fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY           fieldcorr     = 0.5*real(pmuH_e*bx*0.00001*SiDimZ/pitchY) ! EM GCC4.7 pmuH_h is double precision but all the others are real...
134    
135        endif              endif      
136                
# Line 128  c     here bbb is the x component of the Line 149  c     here bbb is the x component of the
149        character*4 PFAtt        character*4 PFAtt
150        include 'commontracker.f'        include 'commontracker.f'
151        include 'level1.f'        include 'level1.f'
152          real corr, res ! EM GCC4.7
153        corr = 0        corr = 0.
154        res  = 0        res  = 0.
155    
156        if(ic.le.0)return        if(ic.le.0)return
157    
# Line 145  c     ================================== Line 166  c     ==================================
166    
167           if(PFAtt.eq.'COG1')then           if(PFAtt.eq.'COG1')then
168    
169              corr   = 0              corr   = 0.
170              res = 1e-4*pitchX/sqrt(12.)!!res              res = REAL(1e-4*pitchX/sqrt(12.))!!res  EM GCC4.7
171    
172           elseif(PFAtt.eq.'COG2')then           elseif(PFAtt.eq.'COG2')then
173    
# Line 212  c            res = riseta(ic,ang)       Line 233  c            res = riseta(ic,ang)      
233  *     temporary patch for saturated clusters  *     temporary patch for saturated clusters
234  *     ======================================  *     ======================================
235           if( nsatstrips(ic).gt.0 )then           if( nsatstrips(ic).gt.0 )then
236              corr  = cog(4,ic)              c            corr  = cog(4,ic)            
237              res = pitchX*1e-4/sqrt(12.)              corr = digsat(ic)
238                res = REAL(pitchX*1e-4/sqrt(12.)) !EM GCC4.7
239  cc            cc=cog(4,ic)  cc            cc=cog(4,ic)
240  c$$$            print*,ic,' *** ',cc  c$$$            print*,ic,' *** ',cc
241  c$$$            print*,ic,' *** ',res  c$$$            print*,ic,' *** ',res
# Line 230  c     ================================== Line 252  c     ==================================
252           if(PFAtt.eq.'COG1')then           if(PFAtt.eq.'COG1')then
253    
254              corr  = 0                corr  = 0  
255              res = 1e-4*pitchY/sqrt(12.)!res                res = REAL(1e-4*pitchY/sqrt(12.))!res  EM GCC4.7
256    
257           elseif(PFAtt.eq.'COG2')then           elseif(PFAtt.eq.'COG2')then
258    
# Line 294  c            res = riseta(ic,ang)   Line 316  c            res = riseta(ic,ang)  
316  *     temporary patch for saturated clusters  *     temporary patch for saturated clusters
317  *     ======================================  *     ======================================
318           if( nsatstrips(ic).gt.0 )then           if( nsatstrips(ic).gt.0 )then
319              corr    = cog(4,ic)              c            corr    = cog(4,ic)            
320              res = pitchY*1e-4/sqrt(12.)              corr = digsat(ic)
321                res = REAL(pitchY*1e-4/sqrt(12.)) ! EM GCC4.7
322  cc            cc=cog(4,ic)  cc            cc=cog(4,ic)
323  c$$$            print*,ic,' *** ',cc  c$$$            print*,ic,' *** ',cc
324  c$$$            print*,ic,' *** ',res  c$$$            print*,ic,' *** ',res
# Line 441  cc            print*,pfaeta2(ic,angle) Line 464  cc            print*,pfaeta2(ic,angle)
464                            
465        endif        endif
466                
467   100  return  c 100  return
468          return
469        end        end
470    
471  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 488  cc            print*,VIEW(ic),angle,pfae Line 512  cc            print*,VIEW(ic),angle,pfae
512                            
513        endif        endif
514                
515   100  return  c 100  return
516          return
517        end        end
518  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
519  c      real function riseta(ic,angle)  c      real function riseta(ic,angle)
# Line 507  c      real function riseta(ic,angle) Line 532  c      real function riseta(ic,angle)
532        include 'level1.f'        include 'level1.f'
533        include 'calib.f'        include 'calib.f'
534    
535        riseta = 0        riseta = 0.
536    
537  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
538        if(mod(iview,2).eq.1)then !Y-view        if(mod(iview,2).eq.1)then !Y-view
# Line 538  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 563  c      if(mod(int(VIEW(ic)),2).eq.1)then
563        endif        endif
564    
565    
566   100  return  c 100  return
567          return
568        end        end
569    
570  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 551  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 577  c      if(mod(int(VIEW(ic)),2).eq.1)then
577  *     resolution.  *     resolution.
578  *     It calls the function FBAD_COG(NCOG,IC),  *     It calls the function FBAD_COG(NCOG,IC),
579  *     accordingto the angle  *     accordingto the angle
580    *
581    *     >>> cosi` non e` corretto!!
582    *     >>> l'errore sulla coordinata eta si ottiene moltiplicando
583    *     >>> l'errore sulla coordinata cog per la derivata della
584    *     >>> distribuzione eta... pur sapendolo l'ho sempre ignorato...
585    *     >>> deve essere modificato!!!!
586    *
587  *-------------------------------------------------------  *-------------------------------------------------------
588    
589        include 'commontracker.f'        include 'commontracker.f'
# Line 588  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 621  c      if(mod(int(VIEW(ic)),2).eq.1)then
621        end        end
622    
623  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
624        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle)
625  *--------------------------------------------------------------  *--------------------------------------------------------------
626  *     this function returns  *     this function returns
627  *  *
# Line 607  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 640  c      if(mod(int(VIEW(ic)),2).eq.1)then
640        real cog2,angle        real cog2,angle
641        integer iview,lad        integer iview,lad
642    
643        iview = VIEW(ic)                    iview   = VIEW(ic)            
644        lad = nld(MAXS(ic),VIEW(ic))        lad     = nld(MAXS(ic),VIEW(ic))
645        cog2 = cog(2,ic)                  cog2    = cog(2,ic)          
646        pfaeta2=cog2        pfaeta2 = cog2
647    
648    *     ----------------
649  *     find angular bin  *     find angular bin
650    *     ----------------
651  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
652        do iang=1,nangbin        do iang=1,nangbin
653           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
# Line 622  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 657  c      if(mod(int(VIEW(ic)),2).eq.1)then
657        enddo        enddo
658        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
659       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
660        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
661        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
662   98   continue                  !jump here if ok   98   continue                  !jump here if ok
663    
664    *     -------------
665    *     within +/-0.5
666    *     -------------
667    
668  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  
 c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2  
 c$$$*         goto 100  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
669        iadd=0        iadd=0
670   10   continue   10   continue
671        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
672           cog2 = cog2 + 1           cog2 = cog2 + 1
673           iadd = iadd + 1           iadd = iadd + 1
674             if(iadd>iaddmax)goto 111
675           goto 10           goto 10
676        endif        endif
677   20   continue   20   continue
678        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
679           cog2 = cog2 - 1           cog2 = cog2 - 1
680           iadd = iadd - 1           iadd = iadd - 1
681             if(iadd<-1*iaddmax)goto 111
682           goto 20           goto 20
683        endif        endif
684          goto 1111
685     111  continue
686          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
687          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
688          cog2=0
689     1111 continue
690    
691  *     --------------------------------  *     --------------------------------
692  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 700  c$$$      endif Line 725  c$$$      endif
725       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
726    
727    
728   100  return  c 100  return
729          return
730        end        end
731    
732  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 726  c$$$      endif Line 752  c$$$      endif
752    
753        iview = VIEW(ic)                    iview = VIEW(ic)            
754        lad = nld(MAXS(ic),VIEW(ic))        lad = nld(MAXS(ic),VIEW(ic))
755        cog3 = cog(3,ic)                    cog3 = cog(3,ic)
756          cc = cog3
757          cog3 = cc
758        pfaeta3=cog3        pfaeta3=cog3
759    
760    *     ----------------
761  *     find angular bin  *     find angular bin
762    *     ----------------
763  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
764        do iang=1,nangbin        do iang=1,nangbin
765  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 740  c         print*,'~~~~~~~~~~~~ ',iang,an Line 770  c         print*,'~~~~~~~~~~~~ ',iang,an
770        enddo        enddo
771        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
772       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
773        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
774        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
775   98   continue                  !jump here if ok   98   continue                  !jump here if ok
776    
777    *     -------------
778    *     within +/-0.5
779    *     -------------
780    
781  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
782        iadd=0        iadd=0
783   10   continue   10   continue
784        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
785           cog3 = cog3 + 1           cog3   = cog3 + 1.
786           iadd = iadd + 1           iadd = iadd + 1
787             if(iadd>iaddmax) goto 111
788           goto 10           goto 10
789        endif        endif
790   20   continue   20   continue
791        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
792           cog3 = cog3 - 1           cog3 = cog3 - 1.
793           iadd = iadd - 1           iadd = iadd - 1
794             if(iadd<-1*iaddmax) goto 111
795           goto 20           goto 20
796        endif        endif
797          goto 1111
798     111  continue
799          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
800          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
801          cog3=0      
802     1111 continue
803    
804  *     --------------------------------  *     --------------------------------
805  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 804  c            print*,'-----',x1,x2,y1,y2 Line 825  c            print*,'-----',x1,x2,y1,y2
825        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
826        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
827    
 c$$$      if(iflag.eq.1)then  
 c$$$         pfaeta2=pfaeta2-1.   !temp  
 c$$$         cog2=cog2-1.           !temp  
 c$$$      endif  
 c$$$      if(iflag.eq.-1)then  
 c$$$         pfaeta2=pfaeta2+1.   !temp  
 c$$$         cog2=cog2+1.           !temp  
 c$$$      endif  
828    
829        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
830       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
831    
832   100  return  c 100  return
833          return
834        end        end
835    
836  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 845  c$$$      endif Line 859  c$$$      endif
859        cog4=cog(4,ic)        cog4=cog(4,ic)
860        pfaeta4=cog4        pfaeta4=cog4
861    
862    *     ----------------
863  *     find angular bin  *     find angular bin
864    *     ----------------
865  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
866        do iang=1,nangbin        do iang=1,nangbin
867  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 856  c         print*,'~~~~~~~~~~~~ ',iang,an Line 872  c         print*,'~~~~~~~~~~~~ ',iang,an
872        enddo        enddo
873        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
874       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
875        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
876        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
877   98   continue                  !jump here if ok   98   continue                  !jump here if ok
878    
879    *     -------------
880    *     within +/-0.5
881    *     -------------
882    
883  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
884        iadd=0        iadd=0
885   10   continue   10   continue
886        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
887           cog4 = cog4 + 1           cog4 = cog4 + 1
888           iadd = iadd + 1           iadd = iadd + 1
889             if(iadd>iaddmax)goto 111
890           goto 10           goto 10
891        endif        endif
892   20   continue   20   continue
893        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
894           cog4 = cog4 - 1           cog4 = cog4 - 1
895           iadd = iadd - 1           iadd = iadd - 1
896             if(iadd<-1*iaddmax)goto 111
897           goto 20           goto 20
898        endif        endif
899          goto 1111
900     111  continue
901          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
902          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
903          cog4=0
904     1111 continue
905    
906  *     --------------------------------  *     --------------------------------
907  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 932  c$$$      endif Line 939  c$$$      endif
939        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
940       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
941    
942   100  return  c 100  return
943          return
944        end        end
945    
946    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
947          real function digsat(ic)
948    *-------------------------------------------------
949    *
950    *          
951    *-------------------------------------------------
952          include 'commontracker.f'
953          include 'calib.f'
954          include 'level1.f'
955          
956          integer nsat
957          real pitchsat
958          
959          nsat = 0
960          pitchsat = 0.
961          iv=VIEW(ic)              
962          istart = INDSTART(IC)
963          istop  = TOTCLLENGTH
964          if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
965          do i = INDMAX(IC),istart,-1
966             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
967         $        .or.
968         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
969                nsat = nsat + 1
970                pitchsat = pitchsat + i - INDMAX(IC)
971             else
972                goto 10
973             endif
974          enddo
975     10   continue
976          do i = INDMAX(IC)+1,istop
977             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
978         $        .or.
979         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
980                nsat = nsat + 1
981                pitchsat = pitchsat + i - INDMAX(IC)
982             else
983                goto 20
984             endif
985          enddo
986     20   continue
987          
988          digsat = 0
989          if (nsat.gt.0) digsat = pitchsat / nsat
990          
991          return
992          end
993    
994    
995  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
996        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 997  c$$$      endif Line 1052  c$$$      endif
1052                    
1053           COG = 0.           COG = 0.
1054                    
1055  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
1056    
1057  c     ==============================================================  c     ==============================================================
1058           if(ncog.eq.1)then           if(ncog.eq.1)then
1059              COG = 0.              COG = 0.
1060              if(sr1.gt.sc)cog=1.           if(sr1.gt.sc)cog=1.
1061              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.           if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1062  c     ==============================================================  c     ==============================================================
1063           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
1064              COG = 0.              COG = 0.
1065              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1066                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
1067              elseif(sl1.lt.sr1)then              elseif(sl1.lt.sr1)then
1068                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                                         if((sc+sr1).ne.0)COG = sr1/(sc+sr1)
1069              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then           elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1070                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1071       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1072                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1073       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1074              endif           endif
1075  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1076  c     $           ,' : ',sl2,sl1,sc,sr1,sr2  c     $           ,' : ',sl2,sl1,sc,sr1,sr2
1077  c     ==============================================================  c     ==============================================================
# Line 1082  ' Line 1137  '
1137  *     =========================  *     =========================
1138    
1139           iv=VIEW(ic)           iv=VIEW(ic)
1140           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=NINT(incuty) ! incut is implicitly INTEGER, incuty is REAL
1141           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=NINT(incutx) ! incut is implicitly INTEGER, incutx is REAL
1142           istart = INDSTART(IC)           istart = INDSTART(IC)
1143           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
1144           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
# Line 1138  c         print*,'-------' Line 1193  c         print*,'-------'
1193    
1194  c      print *,'## cog ',ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1195    
1196          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1197             if(DEBUG.eq.1)
1198         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1199             if(DEBUG.eq.1)
1200         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1201          endif
1202    
1203    
1204        return        return
1205        end        end
1206    
# Line 1160  c      print *,'## cog ',ncog,ic,cog,'// Line 1223  c      print *,'## cog ',ncog,ic,cog,'//
1223        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1224           si = 8.4  !average good-strip noise           si = 8.4  !average good-strip noise
1225           f  = 4.   !average bad-strip noise: f*si           f  = 4.   !average bad-strip noise: f*si
1226           incut=incuty           incut=NINT(incuty)
1227        else                      !X-view        else                      !X-view
1228           si = 3.9  !average good-strip noise           si = 3.9  !average good-strip noise
1229           f  = 6.   !average bad-strip noise: f*si           f  = 6.   !average bad-strip noise: f*si
1230           incut=incutx           incut=NINT(incutx)
1231        endif        endif
1232                
1233        fbad_cog = 1.        fbad_cog = 1.
# Line 1330  c            COG = 0. Line 1393  c            COG = 0.
1393        end        end
1394    
1395    
1396  c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1397  c$$$      real function fbad_cog0(ncog,ic)  
1398  c$$$*-------------------------------------------------------        real function riscogtheor(ncog,ic)
1399  c$$$*     this function returns a factor that takes into  *-------------------------------------------------------
1400  c$$$*     account deterioration of the spatial resolution  *
1401  c$$$*     in the case BAD strips are included in the cluster.  *     this function returns the expected resolution
1402  c$$$*     This factor should multiply the nominal spatial  *     obtained by propagating the strip noise
1403  c$$$*     resolution.  *     to the center-of-gravity coordinate
1404  c$$$*  *
1405  c$$$*     NB!!!  *     ncog = n.strip used in the coordinate evaluation
1406  c$$$*     (this is the old version. It consider only the two  *     (ncog=0 => all strips above threshold)
1407  c$$$*     strips with the greatest signal. The new one is  *
1408  c$$$*     fbad_cog(ncog,ic) )  *-------------------------------------------------------
1409  c$$$*      
1410  c$$$*-------------------------------------------------------        include 'commontracker.f'
1411  c$$$        include 'level1.f'
1412  c$$$      include 'commontracker.f'        include 'calib.f'
1413  c$$$      include 'level1.f'  
1414  c$$$      include 'calib.f'        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1415  c$$$           incut = NINT(incuty) ! EM GCC4.7
1416  c$$$*     --> signal of the central strip           pitch = REAL(pitchY / 1.e4)
1417  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center        else                      !X-view
1418  c$$$           incut = NINT(incutx) ! EM GCC4.7
1419  c$$$*     signal of adjacent strips           pitch = REAL(pitchX / 1.e4)
1420  c$$$*     --> left        endif
1421  c$$$      sl1 = 0                  !left 1        
1422  c$$$      if(        func = 100000.
1423  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)        stot = 0.
1424  c$$$     $     )  
1425  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))        if (ncog.gt.0) then
1426  c$$$  
1427  c$$$      sl2 = 0                  !left 2  *     --> signal of the central strip
1428  c$$$      if(           sc = CLSIGNAL(INDMAX(ic)) !center
1429  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)           fsc = clsigma(INDMAX(ic))
1430  c$$$     $     )  *     --> signal of adjacent strips
1431  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))           sl1 = 0                !left 1
1432  c$$$           fsl1 = 1               !left 1
1433  c$$$*     --> right           if(
1434  c$$$      sr1 = 0                  !right 1       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1435  c$$$      if(       $        )then
1436  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))              sl1 = CLSIGNAL(INDMAX(ic)-1)
1437  c$$$     $     .or.              fsl1 = clsigma(INDMAX(ic)-1)
1438  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)           endif
1439  c$$$     $     )  
1440  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))           sl2 = 0                !left 2
1441  c$$$           fsl2 = 1               !left 2
1442  c$$$      sr2 = 0                  !right 2           if(
1443  c$$$      if(       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1444  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        )then
1445  c$$$     $     .or.              sl2 = CLSIGNAL(INDMAX(ic)-2)
1446  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)              fsl2 = clsigma(INDMAX(ic)-2)
1447  c$$$     $     )           endif
1448  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))           sr1 = 0                !right 1
1449  c$$$           fsr1 = 1               !right 1
1450  c$$$           if(
1451  c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1452  c$$$         f  = 4.       $        .or.
1453  c$$$         si = 8.4       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1454  c$$$      else                              !X-view       $        )then
1455  c$$$         f  = 6.              sr1 = CLSIGNAL(INDMAX(ic)+1)
1456  c$$$         si = 3.9              fsr1 = clsigma(INDMAX(ic)+1)
1457  c$$$      endif           endif    
1458  c$$$           sr2 = 0                !right 2
1459  c$$$      fbad_cog = 1.           fsr2 = 1               !right 2
1460  c$$$      f0 = 1           if(
1461  c$$$      f1 = 1       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1462  c$$$      f2 = 1       $        .or.
1463  c$$$      f3 = 1         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1464  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then       $        )then
1465  c$$$                      sr2 = CLSIGNAL(INDMAX(ic)+2)
1466  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f              fsr2 = clsigma(INDMAX(ic)+2)
1467  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f           endif
1468  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  
1469  c$$$  
1470  c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
1471  c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  ************************************************************
1472  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  *     COG2-3-4 computation
1473  c$$$            fbad_cog = 1.  ************************************************************
1474  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
1475  c$$$            fbad_cog = 1.  c      print*,sl2,sl1,sc,sr1,sr2
1476  c$$$         else          
1477  c$$$            fbad_cog = 1.           vCOG = cog(ncog,ic)!0.
1478  c$$$         endif          
1479  c$$$                   if(ncog.eq.1)then
1480  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then              func = 1./12.
1481  c$$$              stot = 1.
1482  c$$$           elseif(ncog.eq.2)then
1483  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f              if(sl1.gt.sr1)then
1484  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f                 func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1485  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f                 stot = sl1+sc
1486  c$$$              elseif(sl1.le.sr1)then
1487  c$$$         if(ncog.eq.2.and.sr1.ne.0)then                 func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1488  c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)                 stot = sc+sr1
1489  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then              endif
1490  c$$$            fbad_cog = 1.           elseif(ncog.eq.3)then
1491  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then              func =
1492  c$$$            fbad_cog = 1.       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1493  c$$$         else              stot = sl1+sc+sr1
1494  c$$$            fbad_cog = 1.           elseif(ncog.eq.4)then
1495  c$$$         endif              if(sl2.gt.sr2)then
1496  c$$$                 func =
1497  c$$$      endif       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1498  c$$$       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1499  c$$$      fbad_cog0 = sqrt(fbad_cog)                 stot = sl2+sl1+sc+sr1
1500  c$$$              elseif(sl2.le.sr2)then
1501  c$$$      return                 func =
1502  c$$$      end       $              (fsl1*(-1-vCOG)**2
1503  c$$$       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1504  c$$$                 stot = sl2+sl1+sc+sr1
1505  c$$$              endif
1506             else
1507                print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG
1508         $            ,' not implemented'
1509             endif
1510            
1511          elseif(ncog.eq.0)then
1512    *     =========================
1513    *     COG computation
1514    *     =========================
1515            
1516             vCOG = cog(0,ic)
1517    
1518             iv     = VIEW(ic)
1519             istart = INDSTART(IC)
1520             istop  = TOTCLLENGTH
1521             if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1522    ccc         SGN = 0.
1523             SNU = 0.
1524    ccc         SDE = 0.
1525    
1526             do i=INDMAX(IC),istart,-1
1527                ipos = i-INDMAX(ic)
1528                cut  = incut*CLSIGMA(i)
1529                if(CLSIGNAL(i).gt.cut)then
1530                   fs = clsigma(i)
1531                   SNU  = SNU + fs*(ipos-vCOG)**2
1532                   stot = stot + CLSIGNAL(i)
1533                else
1534                   goto 10
1535                endif            
1536             enddo
1537     10      continue
1538             do i=INDMAX(IC)+1,istop
1539                ipos = i-INDMAX(ic)
1540                cut  = incut*CLSIGMA(i)
1541                if(CLSIGNAL(i).gt.cut)then
1542                   fs = clsigma(i)
1543                   SNU  = SNU + fs*(ipos-vCOG)**2
1544                   stot = stot + CLSIGNAL(i)
1545                else
1546                   goto 20
1547                endif            
1548             enddo
1549     20      continue
1550             if(SDE.ne.0)then
1551                FUNC=SNU
1552             else
1553                
1554             endif
1555    
1556          else
1557                      
1558             FUNC=0
1559             print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1560             print*,'                               (NCOG must be >= 0)'
1561            
1562    
1563          endif
1564    
1565    
1566          if(stot.gt.0..and.func.gt.0.)then
1567             func = sqrt(func)
1568             func = pitch * func/stot
1569          endif
1570    
1571          riscogtheor = func
1572    
1573          return
1574          end
1575    
1576    
1577    
1578    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1579    
1580          real function risetatheor(ncog,ic,angle)
1581    *-------------------------------------------------------
1582    *
1583    *     this function returns the expected resolution
1584    *     obtained by propagating the strip noise
1585    *     to the coordinate evaluated with non-linear eta-function
1586    *
1587    *     ncog = n.strip used in the coordinate evaluation
1588    *     (ncog=0 => ncog=2,3,4 according to angle)
1589    *
1590    *-------------------------------------------------------
1591    
1592          include 'commontracker.f'
1593          include 'level1.f'
1594          include 'calib.f'
1595    
1596    
1597          func = 1.
1598    
1599          iview   = VIEW(ic)            
1600          lad     = nld(MAXS(ic),VIEW(ic))
1601    
1602    *     ------------------------------------------------
1603    *     number of strip to be used (in case of ncog = 0)
1604    *     ------------------------------------------------
1605    
1606          inoeta = 0
1607    
1608          if(mod(int(iview),2).eq.1)then !Y-view
1609    
1610             pitch = REAL(pitchY / 1.e4) !EM GCC 4.7
1611    
1612             if(ncog.eq.0)then
1613                if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
1614                   ncog = 2
1615                elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
1616                   ncog = 3
1617                elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
1618                   ncog = 4
1619                else
1620                   ncog   = 4
1621                   inoeta = 1
1622                endif            
1623             endif
1624    
1625          else                      !X-view
1626    
1627             pitch = REAL(pitchX / 1.e4) ! EM GCC4.7
1628    
1629             if(ncog.eq.0)then
1630                if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
1631                   ncog = 2
1632                elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
1633                   ncog = 3
1634                elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
1635                   ncog = 4
1636                else
1637                   ncog = 4
1638                   inoeta = 1
1639                endif            
1640             endif
1641    
1642          endif
1643          
1644          func = riscogtheor(ncog,ic)
1645    
1646          risetatheor = func
1647          
1648          if(inoeta.eq.1)return ! no eta correction is applied --> exit
1649          if(ncog.lt.1.or.ncog.gt.4)return
1650    
1651    *     ----------------
1652    *     find angular bin
1653    *     ----------------
1654    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1655          do iang=1,nangbin
1656             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1657                iangle=iang
1658                goto 98
1659             endif
1660          enddo
1661          if(DEBUG.EQ.1)print*
1662         $     ,'risetatheor *** warning *** angle out of range: ',angle
1663          if(angle.le.angL(1))iang=1
1664          if(angle.ge.angR(nangbin))iang=nangbin
1665     98   continue                  !jump here if ok
1666    
1667    *     -------------
1668    *     within +/-0.5
1669    *     -------------
1670    
1671          vcog = cog(ncog,ic)          
1672    
1673          etamin = eta2(1,iang)
1674          etamax = eta2(netaval,iang)
1675    
1676          iaddmax=10
1677          iadd=0
1678     10   continue
1679          if(vcog.lt.etamin)then
1680             vcog = vcog + 1
1681             iadd = iadd + 1
1682             if(iadd>iaddmax)goto 111
1683             goto 10
1684          endif
1685     20   continue
1686          if(vcog.gt.etamax)then
1687             vcog = vcog - 1
1688             iadd = iadd - 1
1689             if(iadd<-1*iaddmax)goto 111
1690             goto 20
1691          endif
1692          goto 1111
1693     111  continue
1694          if(DEBUG.eq.1)
1695         $     print*,'risetatheor *** warning *** anomalous cluster'
1696          if(DEBUG.eq.1)
1697         $     print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)'
1698          vcog=0
1699     1111 continue
1700    
1701    *     ------------------------------------------------
1702    *     interpolation
1703    *     ------------------------------------------------
1704    
1705    
1706          ibin = netaval
1707          do i=2,netaval        
1708             if(ncog.eq.2)eta=eta2(i,iang)
1709             if(ncog.eq.3)eta=eta3(i,iang)
1710             if(ncog.eq.4)eta=eta4(i,iang)
1711             if(eta.ge.vcog)then
1712                ibin = i
1713                goto 99
1714             endif
1715          enddo
1716     99   continue
1717    
1718          if(ncog.eq.2)then
1719             x1 = eta2(ibin-1,iang)
1720             x2 = eta2(ibin,iang)
1721             y1 = feta2(ibin-1,iview,lad,iang)
1722             y2 = feta2(ibin,iview,lad,iang)
1723          elseif(ncog.eq.3)then
1724             x1 = eta3(ibin-1,iang)
1725             x2 = eta3(ibin,iang)
1726             y1 = feta3(ibin-1,iview,lad,iang)
1727             y2 = feta3(ibin,iview,lad,iang)
1728          elseif(ncog.eq.4)then
1729             x1 = eta4(ibin-1,iang)
1730             x2 = eta4(ibin,iang)
1731             y1 = feta4(ibin-1,iview,lad,iang)
1732             y2 = feta4(ibin,iview,lad,iang)
1733          endif
1734          
1735          func = func * (y2-y1)/(x2-x1)
1736    
1737          risetatheor = func
1738    
1739          return
1740          end
1741    
1742  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1743    
1744        FUNCTION risxeta2(x)        FUNCTION risxeta2(x)
1745    
1746          DOUBLE PRECISION HQUADF ! EM GCC4.7
1747        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1748        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1749        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1530  c$$$ Line 1829  c$$$
1829     20 CONTINUE     20 CONTINUE
1830        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1831                
1832        risxeta2=HQUADF* 1e-4        risxeta2=REAL(HQUADF* 1e-4) ! EM GCC4.7 all computation here are done in double precision but the function returns REAL since it is undefined and it is used in the code in single precision variables
1833    
1834        END        END
1835    
1836  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1837        FUNCTION risxeta3(x)        FUNCTION risxeta3(x)
1838          DOUBLE PRECISION HQUADF ! EM GCC4.7
1839        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1840        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1841        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1621  c$$$ Line 1921  c$$$
1921     20 CONTINUE     20 CONTINUE
1922        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1923                
1924        risxeta3 = HQUADF* 1e-4        risxeta3 = REAL(HQUADF* 1e-4) ! EM GCC4.7 all computation here are done in double precision but the function returns REAL since it is undefined and it is used in the code in single precision variables
1925    
1926        END        END
1927  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1928        FUNCTION risxeta4(x)        FUNCTION risxeta4(x)
1929          DOUBLE PRECISION HQUADF ! EM GCC4.7
1930        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1931        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1932        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1711  c$$$ Line 2012  c$$$
2012     20 CONTINUE     20 CONTINUE
2013        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2014                
2015        risxeta4=HQUADF* 1e-4        risxeta4=REAL(HQUADF* 1e-4) ! EM GCC4.7 all computation here are done in double precision but the function returns REAL since it is undefined and it is used in the code in single precision variables
2016    
2017        END        END
2018  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2019        FUNCTION risyeta2(x)        FUNCTION risyeta2(x)
2020          DOUBLE PRECISION HQUADF ! EM GCC4.7
2021        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
2022        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
2023        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1783  c$$$ Line 2085  c$$$
2085     20 CONTINUE     20 CONTINUE
2086        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2087    
2088        risyeta2=HQUADF* 1e-4        risyeta2=REAL(HQUADF* 1e-4) ! EM GCC4.7 all computation here are done in double precision but the function returns REAL since it is undefined and it is used in the code in single precision variables
2089    
2090        END        END
2091  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2092    
2093        FUNCTION risy_cog(x)        FUNCTION risy_cog(x)
2094          DOUBLE PRECISION HQUADF ! EM GCC4.7
2095        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
2096        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
2097        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1850  c$$$ Line 2153  c$$$
2153     20 CONTINUE     20 CONTINUE
2154        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2155    
2156        risy_cog=HQUADF* 1e-4        risy_cog=REAL(HQUADF* 1e-4) ! EM GCC4.7 all computation here are done in double precision but the function returns REAL since it is undefined and it is used in the code in single precision variables
2157                
2158        END        END
2159  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2160        FUNCTION risx_cog(x)        FUNCTION risx_cog(x)
2161          DOUBLE PRECISION HQUADF ! EM GCC4.7
2162        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
2163        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
2164        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1931  c$$$ Line 2235  c$$$
2235     20 CONTINUE     20 CONTINUE
2236        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2237    
2238        risx_cog = HQUADF * 1e-4        risx_cog = REAL(HQUADF * 1e-4) ! EM GCC4.7 all computation here are done in double precision but the function returns REAL since it is undefined and it is used in the code in single precision variables
2239    
2240        END        END
2241    
# Line 1961  c$$$ Line 2265  c$$$
2265        enddo        enddo
2266        if(DEBUG.eq.1)        if(DEBUG.eq.1)
2267       $     print*,'pfacorr *** warning *** angle out of range: ',angle       $     print*,'pfacorr *** warning *** angle out of range: ',angle
2268        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
2269        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
2270   98   continue                  !jump here if ok   98   continue                  !jump here if ok
2271    
2272        pfacorr = fcorr(iview,lad,iang)        pfacorr = fcorr(iview,lad,iang)
2273    
2274        if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr        if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr
   
2275    
2276   100  return        
2277    c 100  return
2278          return
2279        end        end
Legend:
Removed from v.1.21  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.23