/[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.7 by pam-fi, Thu Jan 11 10:20:58 2007 UTC revision 1.15 by pam-fi, Fri Aug 17 13:28:02 2007 UTC
# Line 1  Line 1 
1    
2    
3          subroutine idtoc(ipfa,cpfa)
4          
5          integer ipfa
6          character*4 cpfa
7    
8          CPFA='COG4'
9          if(ipfa.eq.0)CPFA='ETA'
10          if(ipfa.eq.2)CPFA='ETA2'
11          if(ipfa.eq.3)CPFA='ETA3'
12          if(ipfa.eq.4)CPFA='ETA4'
13          if(ipfa.eq.10)CPFA='COG'
14          if(ipfa.eq.11)CPFA='COG1'
15          if(ipfa.eq.12)CPFA='COG2'
16          if(ipfa.eq.13)CPFA='COG3'
17          if(ipfa.eq.14)CPFA='COG4'
18          
19          end
20    
21  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
22  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
23  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms
# Line 6  Line 26 
26  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
27    
28    
29          integer function npfastrips(ic,PFA,angle)
30    *--------------------------------------------------------------
31    *     thid function returns the number of strips used
32    *     to evaluate the position of a cluster, according to the p.f.a.
33    *--------------------------------------------------------------
34          include 'commontracker.f'
35          include 'level1.f'
36          include 'calib.f'
37    
38          character*4 usedPFA,PFA
39    
40    
41          usedPFA=PFA
42    
43          npfastrips=0
44    
45          if(usedPFA.eq.'COG1')npfastrips=1
46          if(usedPFA.eq.'COG2')npfastrips=2
47          if(usedPFA.eq.'COG3')npfastrips=3
48          if(usedPFA.eq.'COG4')npfastrips=4
49          if(usedPFA.eq.'ETA2')npfastrips=2
50          if(usedPFA.eq.'ETA3')npfastrips=3
51          if(usedPFA.eq.'ETA4')npfastrips=4
52    *     ----------------------------------------------------------------
53          if(usedPFA.eq.'ETA')then
54    c         print*,VIEW(ic),angle
55             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
56                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
57                   npfastrips=2
58                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
59                   npfastrips=3
60                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
61                   npfastrips=4
62                else
63                   npfastrips=4
64    c               usedPFA='COG'
65                endif                        
66             else                   !X-view
67                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
68                   npfastrips=2
69                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
70                   npfastrips=3
71                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
72                   npfastrips=4
73                else
74                   npfastrips=4
75    c               usedPFA='COG'
76                endif                        
77             endif
78          endif
79    *     ----------------------------------------------------------------
80          if(usedPFA.eq.'COG')then
81    
82             iv=VIEW(ic)
83             if(mod(iv,2).eq.1)incut=incuty
84             if(mod(iv,2).eq.0)incut=incutx
85             istart = INDSTART(IC)
86             istop  = TOTCLLENGTH
87             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
88             mu  = 0
89             do i = INDMAX(IC),istart,-1
90                ipos = i-INDMAX(ic)
91                cut  = incut*CLSIGMA(i)
92                if(CLSIGNAL(i).ge.cut)then
93                   mu = mu + 1
94                   print*,i,mu
95                else
96                   goto 10
97                endif
98             enddo
99     10      continue
100             do i = INDMAX(IC)+1,istop
101                ipos = i-INDMAX(ic)
102                cut  = incut*CLSIGMA(i)
103                if(CLSIGNAL(i).ge.cut)then
104                   mu = mu + 1
105                   print*,i,mu
106                else
107                   goto 20
108                endif
109             enddo
110     20      continue
111             npfastrips=mu
112    
113          endif
114    *     ----------------------------------------------------------------
115    
116    c      print*,pfastrips
117    
118          return
119          end
120    
121  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
122        real function pfaeta(ic,angle)        real function pfaeta(ic,angle)
123  *--------------------------------------------------------------  *--------------------------------------------------------------
# Line 17  Line 129 
129  *     according to the angle  *     according to the angle
130  *--------------------------------------------------------------  *--------------------------------------------------------------
131        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
132        include 'level1.f'        include 'level1.f'
133          include 'calib.f'
134                
135        pfaeta = 0        pfaeta = 0
136    
137        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
138                
139           pfaeta = pfaeta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
140                  pfaeta = pfaeta2(ic,angle)
141             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
142                pfaeta = pfaeta3(ic,angle)
143             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
144                pfaeta = pfaeta4(ic,angle)
145             else
146                pfaeta = cog(4,ic)
147             endif            
148    
149        else                      !X-view        else                      !X-view
150    
151           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
152              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
153           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
154              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
155           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
156              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
157           endif           else
158                pfaeta = cog(4,ic)
159             endif            
160                            
161        endif        endif
162                
 c      print*,'pfaeta ',pfaeta, angle  
   
163   100  return   100  return
164        end        end
165    
166  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
167        real function ris_eta(ic,angle)  c      real function riseta(ic,angle)
168          real function riseta(iview,angle)
169  *--------------------------------------------------------------  *--------------------------------------------------------------
170  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
171  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
172  *     it calls:  *     it calls:
173  *     - risx_eta2(angle)  *     - risxeta2(angle)
174  *     - risy_eta2(angle)  *     - risyeta2(angle)
175  *     - risx_eta3(angle)  *     - risxeta3(angle)
176  *     - risx_eta4(angle)  *     - risxeta4(angle)
177  *     according to the angle  *     according to the angle
178  *--------------------------------------------------------------  *--------------------------------------------------------------
179        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
180        include 'level1.f'        include 'level1.f'
181          include 'calib.f'
182    
183  c$$$      logical DEBUG        riseta = 0
 c$$$      common/dbg/DEBUG  
         
       ris_eta = 0  
184    
185        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
186          if(mod(iview,2).eq.1)then !Y-view
187                
188           ris_eta = risy_eta2(angle)  
189           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
190                riseta = risyeta2(angle)
191             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
192                riseta = risy_cog(angle) !ATTENZIONE!!
193             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
194                riseta = risy_cog(angle) !ATTENZIONE!!
195             else
196                riseta = risy_cog(angle)
197             endif            
198    
199        else                      !X-view        else                      !X-view
200    
201           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
202              ris_eta = risx_eta2(angle)              riseta = risxeta2(angle)
203           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
204              ris_eta = risx_eta3(angle)              riseta = risxeta3(angle)
205           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
206              ris_eta = risx_eta4(angle)              riseta = risxeta4(angle)
207           elseif(abs(angle).gt.21.)then           else
208              ris_eta = risx_eta4(21.)              riseta = risx_cog(angle)
209           endif           endif            
210                            
211        endif        endif
212    
 c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'  
 c$$$     $     ,' -->',ris_eta  
   
213    
214   100  return   100  return
215        end        end
# Line 104  c$$$     $     ,' -->',ris_eta Line 228  c$$$     $     ,' -->',ris_eta
228    
229        include 'commontracker.f'        include 'commontracker.f'
230        include 'level1.f'        include 'level1.f'
231  *      include 'calib.f'        include 'calib.f'
232        fbad_eta = 0        fbad_eta = 0
233    
234        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
235                
236           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
237                fbad_eta = fbad_cog(2,ic)
238             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
239                fbad_eta = fbad_cog(3,ic)
240             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
241                fbad_eta = fbad_cog(4,ic)
242             else
243                fbad_eta = fbad_cog(4,ic)
244             endif            
245    
246        else                      !X-view        else                      !X-view
247    
248           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
249              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
250           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
251              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
252           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
253              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
254           endif           else
255                fbad_eta = fbad_cog(4,ic)
256             endif            
257                            
258        endif        endif
259    
# Line 127  c$$$     $     ,' -->',ris_eta Line 261  c$$$     $     ,' -->',ris_eta
261        end        end
262    
263  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta2(cog2,view,lad,angle)  
264        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle) !(1)
265  *--------------------------------------------------------------  *--------------------------------------------------------------
266  *     this function returns  *     this function returns
# Line 150  c      real function pfaeta2(cog2,view,l Line 280  c      real function pfaeta2(cog2,view,l
280        real cog2,angle        real cog2,angle
281        integer iview,lad        integer iview,lad
282    
283  c      logical DEBUG        iview = VIEW(ic)            
284  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
285          cog2 = cog(2,ic)          
 c      print*,'## pfaeta2 ',ic,angle  
       iview = VIEW(ic)              !(1)  
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
286        pfaeta2=cog2        pfaeta2=cog2
287    
288  *     find angular bin  *     find angular bin
289  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
290        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
291           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
292              iangle=iang              iangle=iang
293              goto 98              goto 98
# Line 252  c$$$      endif Line 377  c$$$      endif
377        end        end
378    
379  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta3(cog3,view,lad,angle)  
380        real function pfaeta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
381  *--------------------------------------------------------------  *--------------------------------------------------------------
382  *     this function returns  *     this function returns
# Line 275  c      real function pfaeta3(cog3,view,l Line 396  c      real function pfaeta3(cog3,view,l
396        real cog3,angle        real cog3,angle
397        integer iview,lad        integer iview,lad
398    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 c      print*,'## pfaeta3 ',ic,angle  
399    
400        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
401        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
402        cog3 = cog(3,ic)             !(1)        cog3 = cog(3,ic)            
403        pfaeta3=cog3        pfaeta3=cog3
404    
405  *     find angular bin  *     find angular bin
# Line 376  c$$$      endif Line 493  c$$$      endif
493        end        end
494    
495  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
496  c*****************************************************        real function pfaeta4(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta4(cog4,view,lad,angle)  
       real function pfaeta4(ic,angle) !(1)  
497  *--------------------------------------------------------------  *--------------------------------------------------------------
498  *     this function returns  *     this function returns
499  *  *
# Line 399  c      real function pfaeta4(cog4,view,l Line 512  c      real function pfaeta4(cog4,view,l
512        real cog4,angle        real cog4,angle
513        integer iview,lad        integer iview,lad
514    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
515    
516  c      print*,'## pfaeta4 ',ic,angle        iview = VIEW(ic)            
517          lad = nld(MAXS(ic),VIEW(ic))
518        iview = VIEW(ic)             !(1)        cog4=cog(4,ic)
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog4=cog(4,ic)               !(1)  
519        pfaeta4=cog4        pfaeta4=cog4
520    
521  *     find angular bin  *     find angular bin
# Line 502  c$$$      endif Line 611  c$$$      endif
611    
612    
613  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
614        real function cog0(ncog,ic)  c$$$      real function cog0(ncog,ic)
615  *-------------------------------------------------  c$$$*-------------------------------------------------
616  *     this function returns  c$$$*     this function returns
617  *  c$$$*
618  *     - the Center-Of-Gravity of the cluster IC  c$$$*     - the Center-Of-Gravity of the cluster IC
619  *     evaluated using NCOG strips,  c$$$*     evaluated using NCOG strips,
620  *     calculated relative to MAXS(IC)  c$$$*     calculated relative to MAXS(IC)
621  *      c$$$*    
622  *     - zero in case that not  enough strips  c$$$*     - zero in case that not  enough strips
623  *     have a positive signal  c$$$*     have a positive signal
624  *      c$$$*    
625  *     NOTE:  c$$$*     NOTE:
626  *     This is the old definition, used by Straulino.  c$$$*     This is the old definition, used by Straulino.
627  *     The new routine, according to Landi,  c$$$*     The new routine, according to Landi,
628  *     is COG(NCOG,IC)  c$$$*     is COG(NCOG,IC)
629  *-------------------------------------------------  c$$$*-------------------------------------------------
630    c$$$
631    c$$$
632        include 'commontracker.f'  c$$$      include 'commontracker.f'
633        include 'level1.f'  c$$$      include 'level1.f'
634          c$$$      
635  *     --> signal of the central strip  c$$$*     --> signal of the central strip
636        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
637    c$$$
638  *     signal of adjacent strips  c$$$*     signal of adjacent strips
639  *     --> left  c$$$*     --> left
640        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
641        if(  c$$$      if(
642       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
643       $     )  c$$$     $     )
644       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
645    c$$$
646        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
647        if(  c$$$      if(
648       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
649       $     )  c$$$     $     )
650       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
651    c$$$
652  *     --> right  c$$$*     --> right
653        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
654        if(  c$$$      if(
655       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
656       $     .or.  c$$$     $     .or.
657       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
658       $     )  c$$$     $     )
659       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
660    c$$$
661        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
662        if(  c$$$      if(
663       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
664       $     .or.  c$$$     $     .or.
665       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
666       $     )  c$$$     $     )
667       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
668          c$$$      
669  ************************************************************  c$$$************************************************************
670  *     COG computation  c$$$*     COG computation
671  ************************************************************  c$$$************************************************************
672    c$$$
673  c      print*,sl2,sl1,sc,sr1,sr2  c$$$c      print*,sl2,sl1,sc,sr1,sr2
674    c$$$
675        COG = 0.  c$$$      COG = 0.
676            c$$$        
677        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
678            c$$$        
679           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
680              COG = -sl1/(sl1+sc)          c$$$            COG = -sl1/(sl1+sc)        
681           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
682              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
683           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
684              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c$$$            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
685           else  c$$$         else
686              COG = 0.  c$$$            COG = 0.
687           endif  c$$$         endif
688            c$$$        
689        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
690    c$$$
691           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
692              COG = sr1/(sc+sr1)              c$$$            COG = sr1/(sc+sr1)            
693           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
694              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
695           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
696              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c$$$            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
697           else  c$$$         else
698              COG = 0.  c$$$            COG = 0.
699           endif  c$$$         endif
700    c$$$
701        endif  c$$$      endif
702    c$$$
703        COG0 = COG  c$$$      COG0 = COG
704    c$$$
705  c      print *,ncog,ic,cog,'/////////////'  c$$$c      print *,ncog,ic,cog,'/////////////'
706    c$$$
707        return  c$$$      return
708        end  c$$$      end
709    
710  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
711        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 618  c      print *,ncog,ic,cog,'//////////// Line 727  c      print *,ncog,ic,cog,'////////////
727        include 'calib.f'        include 'calib.f'
728        include 'level1.f'        include 'level1.f'
729                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
730    
731    
732        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 630  c      common/dbg/DEBUG Line 737  c      common/dbg/DEBUG
737  *     --> signal of the central strip  *     --> signal of the central strip
738           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
739  *     signal of adjacent strips  *     signal of adjacent strips
740           sl1 = 0                !left 1           sl1 = -9999.           !left 1
741           if(           if(
742       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
743       $        )       $        )
744       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
745                    
746           sl2 = 0                !left 2           sl2 = -9999.           !left 2
747           if(           if(
748       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
749       $        )       $        )
750       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
751                    
752           sr1 = 0                !right 1           sr1 = -9999.           !right 1
753           if(           if(
754       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
755       $        .or.       $        .or.
# Line 650  c      common/dbg/DEBUG Line 757  c      common/dbg/DEBUG
757       $        )       $        )
758       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
759                    
760           sr2 = 0                !right 2           sr2 = -9999.           !right 2
761           if(           if(
762       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
763       $        .or.       $        .or.
# Line 662  c      common/dbg/DEBUG Line 769  c      common/dbg/DEBUG
769                    
770  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c         print*,'## ',sl2,sl1,sc,sr1,sr2
771    
772    c     ==============================================================
773           if(ncog.eq.1)then           if(ncog.eq.1)then
774              COG = 0.              COG = 0.
775                if(sr1.gt.sc)cog=1. !NEW
776                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
777    c     ==============================================================
778           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
779              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
780                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
781              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
782                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
783              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
784                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
785         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
786                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
787         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
788                endif
789    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
790    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
791    c     ==============================================================
792           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
793               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
794    c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
795    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
796    c     ==============================================================
797           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
798              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
799                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sl1+sc+sr1).ne.0)
800       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
801              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
802                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sl1+sc+sr1).ne.0)
803       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
804                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
805                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
806         $              .and.(sl2+sl1+sc+sr1).ne.0 )
807         $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
808                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
809         $              .and.(sr2+sl1+sc+sr1).ne.0 )
810         $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW              
811              endif              endif
812           else           else
813              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
# Line 696  ' Line 825  '
825           iv=VIEW(ic)           iv=VIEW(ic)
826           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
827           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
828           istart = INDSTART(IC)           istart = INDSTART(IC)
829           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
830           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
831           COG = 0             COG = 0  
832             SGN = 0.
833           mu  = 0           mu  = 0
834           do i = istart,istop  c         print*,'-------'
835             do i = INDMAX(IC),istart,-1
836                ipos = i-INDMAX(ic)
837                cut  = incut*CLSIGMA(i)
838                if(CLSIGNAL(i).ge.cut)then              
839                   COG = COG + ipos*CLSIGNAL(i)
840                   SGN = SGN + CLSIGNAL(i)
841                   mu = mu + 1
842    c               print*,ipos,CLSIGNAL(i)
843                else
844                   goto 10
845                endif
846             enddo
847     10      continue
848             do i = INDMAX(IC)+1,istop
849              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
850              cut  = incut*CLSIGMA(i)              cut  = incut*CLSIGMA(i)
851              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
852                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
853                   SGN = SGN + CLSIGNAL(i)
854                 mu = mu + 1                 mu = mu + 1
855    c               print*,ipos,CLSIGNAL(i)
856                else
857                   goto 20
858              endif              endif
859           enddo           enddo
860           if(DEDX(ic).le.0)then   20      continue
861              print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)           if(SGN.le.0)then
862                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
863              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
864              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
865              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
866           else           else
867              COG=COG/DEDX(ic)              COG=COG/SGN
868           endif           endif
869    c         print*,'-------'
870                    
871        else        else
872                    
# Line 734  c      print *,'## cog ',ncog,ic,cog,'// Line 883  c      print *,'## cog ',ncog,ic,cog,'//
883        end        end
884    
885  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
886    
887        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
888  *-------------------------------------------------------  *-------------------------------------------------------
889  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 749  c      print *,'## cog ',ncog,ic,cog,'// Line 899  c      print *,'## cog ',ncog,ic,cog,'//
899        include 'calib.f'        include 'calib.f'
900    
901        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
902           f  = 4.           si = 8.4  !average good-strip noise
903           si = 8.4           f  = 4.   !average bad-strip noise: f*si
904           incut=incuty           incut=incuty
905        else                      !X-view        else                      !X-view
906           f  = 6.           si = 3.9  !average good-strip noise
907           si = 3.9           f  = 6.   !average bad-strip noise: f*si
908           incut=incutx           incut=incutx
909        endif        endif
910                
# Line 765  c      print *,'## cog ',ncog,ic,cog,'// Line 915  c      print *,'## cog ',ncog,ic,cog,'//
915  *     --> signal of the central strip  *     --> signal of the central strip
916           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
917           fsc = 1           fsc = 1
918           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
919             fsc = clsigma(INDMAX(ic))/si
920  *     --> signal of adjacent strips  *     --> signal of adjacent strips
921           sl1 = 0                !left 1           sl1 = 0                !left 1
922           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 773  c      print *,'## cog ',ncog,ic,cog,'// Line 924  c      print *,'## cog ',ncog,ic,cog,'//
924       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
925       $        )then       $        )then
926              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
927              if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f  c            if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f
928  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
929           endif           endif
930    
931           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 784  c            fsl1 = 0 Line 934  c            fsl1 = 0
934       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
935       $        )then       $        )then
936              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
937              if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f  c            if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f
938  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
939           endif           endif
940           sr1 = 0                !right 1           sr1 = 0                !right 1
941           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 796  c            fsl2 = 0 Line 945  c            fsl2 = 0
945       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
946       $        )then       $        )then
947              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
948              if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f  c            if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f
949  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
950           endif               endif    
951           sr2 = 0                !right 2           sr2 = 0                !right 2
952           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 808  c            fsr1 = 0 Line 956  c            fsr1 = 0
956       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
957       $        )then       $        )then
958              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
959              if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f  c            if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f
960  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
961           endif           endif
962    
963    
964    
965  ************************************************************  ************************************************************
966  *     COG computation  *     COG2-3-4 computation
967  ************************************************************  ************************************************************
968    
969  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
970                    
971           COG = 0.           vCOG = cog(ncog,ic)!0.
972                    
973           if(ncog.eq.2)then           if(ncog.eq.2)then
974              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
975                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
976                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
977                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
978              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
979                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
980                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
981                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
982              endif              endif
983           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
984              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
985              fbad_cog =              fbad_cog =
986       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
987              fbad_cog =              fbad_cog =
988       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
989           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
990              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
991                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
992                 fbad_cog =                 fbad_cog =
993       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
994       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
995                 fbad_cog =                 fbad_cog =
996       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
997       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
998              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
999                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1000                 fbad_cog =                 fbad_cog =
1001       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1002       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1003                 fbad_cog =                 fbad_cog =
1004       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1005       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1006              endif              endif
1007           else           else
1008              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1009              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1010              COG = 0.  c            COG = 0.
1011           endif           endif
1012                    
1013        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1014    *     =========================
1015    *     COG computation
1016    *     =========================
1017                    
1018           iv=VIEW(ic)           vCOG = cog(0,ic)
1019           istart=INDSTART(IC)  
1020           istop=TOTCLLENGTH           iv     = VIEW(ic)
1021           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1022           COG=0.           istop  = TOTCLLENGTH
1023           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1024           SDE=0.           SGN = 0.
1025           do i=istart,istop           SNU = 0.
1026              ipos=i-INDMAX(ic)           SDE = 0.
1027              il=nvk(MAXS(ic)+ipos)  c$$$         do i=INDMAX(IC),istart,-1
1028              is=nst(MAXS(ic)+ipos)  c$$$            ipos = i-INDMAX(ic)
1029              cut=incut*SIGMA(iv,il,is)  c$$$            cut  = incut*CLSIGMA(i)
1030    c$$$            if(CLSIGNAL(i).gt.cut)then
1031    c$$$               COG = COG + ipos*CLSIGNAL(i)
1032    c$$$               SGN = SGN + CLSIGNAL(i)
1033    c$$$            else
1034    c$$$               goto 10
1035    c$$$            endif
1036    c$$$         enddo
1037    c$$$ 10      continue
1038    c$$$         do i=INDMAX(IC)+1,istop
1039    c$$$            ipos = i-INDMAX(ic)
1040    c$$$            cut  = incut*CLSIGMA(i)
1041    c$$$            if(CLSIGNAL(i).gt.cut)then
1042    c$$$               COG = COG + ipos*CLSIGNAL(i)
1043    c$$$               SGN = SGN + CLSIGNAL(i)
1044    c$$$            else
1045    c$$$               goto 20
1046    c$$$            endif
1047    c$$$         enddo
1048    c$$$ 20      continue
1049    c$$$         if(SGN.le.0)then
1050    c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN
1051    c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1052    c$$$            print*,(CLSIGNAL(i),i=istart,istop)
1053    c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '
1054    c$$$         else
1055    c$$$            COG=COG/SGN
1056    c$$$         endif
1057    
1058             do i=INDMAX(IC),istart,-1
1059                ipos = i-INDMAX(ic)
1060                cut  = incut*CLSIGMA(i)
1061              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1062                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1063              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1064                   SDE = SDE + (ipos-vCOG)**2
1065                else
1066                   goto 10
1067                endif            
1068           enddo           enddo
1069           COG=COG/DEDX(ic)   10      continue
1070           do i=istart,istop           do i=INDMAX(IC)+1,istop
1071              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1072              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1073              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1074                 fs=1                 fs = clsigma(i)/si
1075                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1076                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1077                 SDE = SDE + (ipos-COG)**2              else
1078                   goto 20
1079              endif                          endif            
1080           enddo           enddo
1081           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1082             if(SDE.ne.0)then
1083                FBAD_COG=SNU/SDE
1084             else
1085                
1086             endif
1087    
1088        else        else
1089                                        
# Line 1026  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1214  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1214    
1215  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1216    
1217        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1218    
1219        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1220        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1098  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1286  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1286       +/       +/
1287    
1288        V(1)= abs(x)        V(1)= abs(x)
1289          if(V(1).gt.20.)V(1)=20.
1290    
1291        HQUADF = 0.        HQUADF = 0.
1292        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1112  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1301  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1301     20 CONTINUE     20 CONTINUE
1302        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1303                
1304        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1305    
1306        END        END
1307    
1308  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1309        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1310        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1311        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1312        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1188  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1377  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1377       +/       +/
1378    
1379        V(1) =  abs(x)        V(1) =  abs(x)
1380          if(V(1).gt.20.)V(1)=20.
1381    
1382        HQUADF = 0.        HQUADF = 0.
1383        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1384           HQDJ = 0.           HQDJ = 0.
# Line 1201  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1392  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1392     20 CONTINUE     20 CONTINUE
1393        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1394                
1395        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1396    
1397        END        END
1398  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1399        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1400        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1401        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1402        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1276  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1467  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1467       +/       +/
1468    
1469        V(1)=abs(x)        V(1)=abs(x)
1470          if(V(1).gt.20.)V(1)=20.
1471    
1472        HQUADF = 0.        HQUADF = 0.
1473        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1290  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1482  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1482     20 CONTINUE     20 CONTINUE
1483        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1484                
1485        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1486    
1487        END        END
1488  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1489        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1490        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1491        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1492        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1347  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1539  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1539       +/       +/
1540    
1541        v(1)= abs(x)        v(1)= abs(x)
1542          if(V(1).gt.20.)V(1)=20.
1543                
1544        HQUADF = 0.        HQUADF = 0.
1545        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1361  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1554  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1554     20 CONTINUE     20 CONTINUE
1555        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1556    
1557        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1558    
1559        END        END
1560  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1413  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1606  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1606       +/       +/
1607                
1608        V(1)=abs(x)        V(1)=abs(x)
1609          if(V(1).gt.20.)V(1)=20.
1610    
1611        HQUADF = 0.        HQUADF = 0.
1612        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1493  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1687  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1687       +/       +/
1688    
1689        V(1)=abs(x)        V(1)=abs(x)
1690          if(V(1).gt.20.)V(1)=20.
1691    
1692        HQUADF = 0.        HQUADF = 0.
1693        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1510  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1705  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1705        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1706    
1707        END        END
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.23