/[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.20 by pam-fi, Mon Aug 27 12:57:15 2007 UTC revision 1.28 by mocchiut, Thu Jan 16 15:29:54 2014 UTC
# Line 4  Line 4 
4  *            *          
5  *     subroutine idtoc(ipfa,cpfa)  *     subroutine idtoc(ipfa,cpfa)
6  *  *
7    *     subroutine applypfa(PFAtt,ic,ang,corr,res)
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 13  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)
50    *     real function fieldcorr(iview,bbb)
51    *
52  *     NB - The angle is the "effective angle", which is relative  *     NB - The angle is the "effective angle", which is relative
53  *          to the sensor and it takes into account the magnetic field  *          to the sensor and it takes into account the magnetic field
54  *  *
# Line 34  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 49  Line 73 
73        if(ipfa.eq.14)CPFA='COG4'        if(ipfa.eq.14)CPFA='COG4'
74                
75        end        end
76    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
77          real function effectiveangle(ang,iview,bbb)
78          
79          include 'commontracker.f'
80          real tgtemp
81    
82          effectiveangle = 0.
83    
84          if(mod(iview,2).eq.0)then
85    c     =================================================
86    c     X view
87    c     =================================================
88    c     here bbb is the y component of the m.field
89             angx = ang
90             by   = bbb
91             if(iview.eq.12) angx = -1. * ang
92             if(iview.eq.12) by   = -1. * bbb
93    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
97    c     =================================================
98    c     Y view
99    c     =================================================        
100    c     here bbb is the x component of the m.filed
101             angy = ang
102             bx   = bbb
103             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      
106          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
107    
108          return
109          end
110    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
111          real function fieldcorr(iview,bbb)
112    
113          include 'commontracker.f'
114    
115          fieldcorr = 0.
116    
117          if(mod(iview,2).eq.0)then
118    
119    c     =================================================
120    c     X view
121    c     =================================================
122    c     here bbb is the y component of the m.field
123             by   = bbb
124             if(iview.eq.12) by = -1. * bbb
125             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
128    c     =================================================
129    c     Y view
130    c     =================================================        
131    c     here bbb is the x component of the m.filed
132             bx   = bbb
133             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      
136          
137          return
138          end
139    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
140    
141          subroutine applypfa(PFAtt,ic,ang,corr,res)
142    *---------------------------------------------------------------
143    *     this subroutine calculate the coordinate of cluster ic (in
144    *     strip units), relative to the strip with the maximum signal,
145    *     and its spatial resolution (in cm), applying PFAtt.
146    *     ang is the effective angle, relative to the sensor
147    *---------------------------------------------------------------
148    
149          character*4 PFAtt
150          include 'commontracker.f'
151          include 'level1.f'
152          real corr, res ! EM GCC4.7
153          corr = 0.
154          res  = 0.
155    
156          if(ic.le.0)return
157    
158          iview   = VIEW(ic)
159    
160          if(mod(iview,2).eq.0)then
161    c     =================================================
162    c     X view
163    c     =================================================
164    
165             res = RESXAV
166    
167             if(PFAtt.eq.'COG1')then
168    
169                corr   = 0.
170                res = REAL(1e-4*pitchX/sqrt(12.))!!res  EM GCC4.7
171    
172             elseif(PFAtt.eq.'COG2')then
173    
174                corr    = cog(2,ic)            
175                res = risx_cog(abs(ang))!TEMPORANEO              
176                res = res*fbad_cog(2,ic)
177    
178             elseif(PFAtt.eq.'COG3')then
179    
180                corr    = cog(3,ic)            
181                res = risx_cog(abs(ang))!TEMPORANEO                      
182                res = res*fbad_cog(3,ic)
183    
184             elseif(PFAtt.eq.'COG4')then
185    
186                corr    = cog(4,ic)            
187                res = risx_cog(abs(ang))!TEMPORANEO                      
188                res = res*fbad_cog(4,ic)
189    
190             elseif(PFAtt.eq.'ETA2')then
191    
192                corr  = pfaeta2(ic,ang)          
193                res = risxeta2(abs(ang))
194                res = res*fbad_cog(2,ic)
195    
196             elseif(PFAtt.eq.'ETA3')then                        
197    
198                corr  = pfaeta3(ic,ang)          
199                res = risxeta3(abs(ang))                      
200                res = res*fbad_cog(3,ic)              
201    
202             elseif(PFAtt.eq.'ETA4')then                        
203    
204                corr  = pfaeta4(ic,ang)            
205                res = risxeta4(abs(ang))                      
206                res = res*fbad_cog(4,ic)              
207    
208             elseif(PFAtt.eq.'ETA')then  
209    
210                corr  = pfaeta(ic,ang)            
211    c            res = riseta(ic,ang)                    
212                res = riseta(iview,ang)                    
213                res = res*fbad_eta(ic,ang)            
214    
215             elseif(PFAtt.eq.'ETAL')then  
216    
217                corr  = pfaetal(ic,ang)            
218                res = riseta(iview,ang)                    
219                res = res*fbad_eta(ic,ang)            
220    
221             elseif(PFAtt.eq.'COG')then          
222    
223                corr  = cog(0,ic)            
224                res = risx_cog(abs(ang))                    
225                res = res*fbad_cog(0,ic)
226    
227             else
228                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
229             endif
230    
231    
232    *     ======================================
233    *     temporary patch for saturated clusters
234    *     ======================================
235             if( nsatstrips(ic).gt.0 )then
236    c            corr  = cog(4,ic)            
237                corr = digsat(ic)
238                res = REAL(pitchX*1e-4/sqrt(12.)) !EM GCC4.7
239    cc            cc=cog(4,ic)
240    c$$$            print*,ic,' *** ',cc
241    c$$$            print*,ic,' *** ',res
242             endif
243    
244    
245          elseif(mod(iview,2).eq.1)then
246    c     =================================================
247    c     Y view
248    c     =================================================
249    
250             res = RESYAV
251    
252             if(PFAtt.eq.'COG1')then
253    
254                corr  = 0  
255                res = REAL(1e-4*pitchY/sqrt(12.))!res  EM GCC4.7
256    
257             elseif(PFAtt.eq.'COG2')then
258    
259                corr    = cog(2,ic)
260                res = risy_cog(abs(ang))!TEMPORANEO
261                res = res*fbad_cog(2,ic)
262    
263             elseif(PFAtt.eq.'COG3')then
264    
265                corr    = cog(3,ic)
266                res = risy_cog(abs(ang))!TEMPORANEO
267                res = res*fbad_cog(3,ic)
268    
269             elseif(PFAtt.eq.'COG4')then
270    
271                corr    = cog(4,ic)
272                res = risy_cog(abs(ang))!TEMPORANEO
273                res = res*fbad_cog(4,ic)
274    
275             elseif(PFAtt.eq.'ETA2')then
276    
277                corr    = pfaeta2(ic,ang)          
278                res = risyeta2(abs(ang))              
279                res = res*fbad_cog(2,ic)
280    
281             elseif(PFAtt.eq.'ETA3')then                      
282    
283                corr    = pfaeta3(ic,ang)
284                res = res*fbad_cog(3,ic)  
285    
286             elseif(PFAtt.eq.'ETA4')then  
287    
288                corr    = pfaeta4(ic,ang)
289                res = res*fbad_cog(4,ic)
290    
291             elseif(PFAtt.eq.'ETA')then
292    
293                corr    = pfaeta(ic,ang)
294    c            res = riseta(ic,ang)  
295                res = riseta(iview,ang)  
296                res = res*fbad_eta(ic,ang)
297    
298             elseif(PFAtt.eq.'ETAL')then
299    
300                corr    = pfaetal(ic,ang)
301                res = riseta(iview,ang)  
302                res = res*fbad_eta(ic,ang)
303    
304             elseif(PFAtt.eq.'COG')then
305    
306                corr    = cog(0,ic)            
307                res = risy_cog(abs(ang))
308                res = res*fbad_cog(0,ic)
309    
310             else
311                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
312             endif
313    
314    
315    *     ======================================
316    *     temporary patch for saturated clusters
317    *     ======================================
318             if( nsatstrips(ic).gt.0 )then
319    c            corr    = cog(4,ic)            
320                corr = digsat(ic)
321                res = REAL(pitchY*1e-4/sqrt(12.)) ! EM GCC4.7
322    cc            cc=cog(4,ic)
323    c$$$            print*,ic,' *** ',cc
324    c$$$            print*,ic,' *** ',res
325             endif
326            
327          endif
328          end
329    
330    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
331        integer function npfastrips(ic,angle)        integer function npfastrips(ic,angle)
332  *--------------------------------------------------------------  *--------------------------------------------------------------
333  *     thid function returns the number of strips used  *     thid function returns the number of strips used
# Line 188  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 235  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 254  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 285  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 298  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 335  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 354  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 369  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 447  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 473  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 487  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 551  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 592  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 603  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 679  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 744  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. !NEW           if(sr1.gt.sc)cog=1.
1061              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW           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.
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 !NEW           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) !NEW       $              .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) !NEW       $              .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     ==============================================================
1078           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1079               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
1080  c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)              sss = sc
1081                if( sl1.ne.-9999. )COG = COG-sl1
1082                if( sl1.ne.-9999. )sss = sss+sl1
1083                if( sr1.ne.-9999. )COG = COG+sr1
1084                if( sr1.ne.-9999. )sss = sss+sr1
1085                if(sss.ne.0)COG=COG/sss
1086    
1087    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1088    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1089  c     $            ,' : ',sl2,sl1,sc,sr1,sr2  c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1090  c     ==============================================================  c     ==============================================================
1091           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1092    
1093                COG = 0
1094                sss = sc
1095                if( sl1.ne.-9999. )COG = COG-sl1
1096                if( sl1.ne.-9999. )sss = sss+sl1
1097                if( sr1.ne.-9999. )COG = COG+sr1
1098                if( sr1.ne.-9999. )sss = sss+sr1
1099              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1100                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1101       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1102              elseif(sl2.lt.sr2)then              elseif(sl2.lt.sr2)then
1103                 if((sr2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1104       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1105              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1106                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1107       $              .and.(sl2+sl1+sc+sr1).ne.0 )       $              .and.(sl2+sss).ne.0 )
1108       $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW       $              cog = (cog-2*sl2)/(sl2+sss)
1109                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1110       $              .and.(sr2+sl1+sc+sr1).ne.0 )       $              .and.(sr2+sss).ne.0 )
1111       $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW                     $              cog = (2*sr2+cog)/(sr2+sss)              
1112              endif              endif
1113    c     ==============================================================
1114             elseif(ncog.eq.5)then
1115                COG = 0
1116                sss = sc
1117                if( sl1.ne.-9999. )COG = COG-sl1
1118                if( sl1.ne.-9999. )sss = sss+sl1
1119                if( sr1.ne.-9999. )COG = COG+sr1
1120                if( sr1.ne.-9999. )sss = sss+sr1
1121                if( sl2.ne.-9999. )COG = COG-2*sl2
1122                if( sl2.ne.-9999. )sss = sss+sl2
1123                if( sr2.ne.-9999. )COG = COG+2*sr2
1124                if( sr2.ne.-9999. )sss = sss+sr2
1125                if(sss.ne.0)COG=COG/sss
1126           else           else
1127              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1128       $           ,' not implemented'       $           ,' not implemented'
# Line 800  ' 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 856  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 878  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 1048  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 1248  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 1339  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 1429  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 1501  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 1568  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 1649  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 1679  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.20  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.23