/[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.15 by pam-fi, Fri Aug 17 13:28:02 2007 UTC revision 1.23 by pam-fi, Wed Mar 5 17:00:20 2008 UTC
# Line 1  Line 1 
1    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2    *     this file contains all subroutines and functions
3    *     that are needed for position finding algorithms:
4    *          
5    *     subroutine idtoc(ipfa,cpfa)
6    *
7    *     subroutine applypfa(PFAtt,ic,ang,corr,res)
8    *
9    *     integer function npfastrips(ic,angle)
10    *
11    *     real function pfaeta(ic,angle)
12    *     real function pfaetal(ic,angle)
13    *     real function pfaeta2(ic,angle)
14    *     real function pfaeta3(ic,angle)
15    *     real function pfaeta4(ic,angle)
16    *     real function cog(ncog,ic)
17    *
18    *     real function fbad_cog(ncog,ic)
19    *     real function fbad_eta(ic,angle)
20    *
21    *     real function riseta(iview,angle)
22    *     FUNCTION risxeta2(x)
23    *     FUNCTION risxeta3(x)
24    *     FUNCTION risxeta4(x)
25    *     FUNCTION risyeta2(x)
26    *     FUNCTION risy_cog(x)
27    *     FUNCTION risx_cog(x)
28    *
29    *     real function pfacorr(ic,angle)
30    *
31    *     real function effectiveangle(ang,iview,bbb)
32    *     real function fieldcorr(iview,bbb)
33    *
34    *     NB - The angle is the "effective angle", which is relative
35    *          to the sensor and it takes into account the magnetic field
36    *
37    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
38    
39        subroutine idtoc(ipfa,cpfa)        subroutine idtoc(ipfa,cpfa)
40                
41        integer ipfa        integer ipfa
42        character*4 cpfa        character*10 cpfa
43    
44        CPFA='COG4'        CPFA='COG4'
45        if(ipfa.eq.0)CPFA='ETA'        if(ipfa.eq.0)CPFA='ETA'
46        if(ipfa.eq.2)CPFA='ETA2'        if(ipfa.eq.2)CPFA='ETA2'
47        if(ipfa.eq.3)CPFA='ETA3'        if(ipfa.eq.3)CPFA='ETA3'
48        if(ipfa.eq.4)CPFA='ETA4'        if(ipfa.eq.4)CPFA='ETA4'
49          if(ipfa.eq.5)CPFA='ETAL'
50        if(ipfa.eq.10)CPFA='COG'        if(ipfa.eq.10)CPFA='COG'
51        if(ipfa.eq.11)CPFA='COG1'        if(ipfa.eq.11)CPFA='COG1'
52        if(ipfa.eq.12)CPFA='COG2'        if(ipfa.eq.12)CPFA='COG2'
# Line 17  Line 54 
54        if(ipfa.eq.14)CPFA='COG4'        if(ipfa.eq.14)CPFA='COG4'
55                
56        end        end
57    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
58          real function effectiveangle(ang,iview,bbb)
59    
60          include 'commontracker.f'
61    
62          effectiveangle = 0.
63    
64          if(mod(iview,2).eq.0)then
65    c     =================================================
66    c     X view
67    c     =================================================
68    c     here bbb is the y component of the m.field
69             angx = ang
70             by   = bbb
71             if(iview.eq.12) angx = -1. * ang
72             if(iview.eq.12) by   = -1. * bbb
73    cc         tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!!
74             tgtemp  = tan(angx*acos(-1.)/180.) + pmuH_h*by*0.00001
75    
76          elseif(mod(iview,2).eq.1)then
77    c     =================================================
78    c     Y view
79    c     =================================================        
80    c     here bbb is the x component of the m.filed
81             angy = ang
82             bx   = bbb
83             tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001        
84    
85          endif      
86          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
87    
88          return
89          end
90  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
91  *     this file contains all subroutines and functions        real function fieldcorr(iview,bbb)
92  *     that are needed for position finding algorithms  
93  *            include 'commontracker.f'
94  *  
95          fieldcorr = 0.
96    
97          if(mod(iview,2).eq.0)then
98    
99    c     =================================================
100    c     X view
101    c     =================================================
102    c     here bbb is the y component of the m.field
103             by   = bbb
104             if(iview.eq.12) by = -1. * bbb
105             fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
106    
107          elseif(mod(iview,2).eq.1)then
108    c     =================================================
109    c     Y view
110    c     =================================================        
111    c     here bbb is the x component of the m.filed
112             bx   = bbb
113             fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
114    
115          endif      
116          
117          return
118          end
119  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
120    
121          subroutine applypfa(PFAtt,ic,ang,corr,res)
122    *---------------------------------------------------------------
123    *     this subroutine calculate the coordinate of cluster ic (in
124    *     strip units), relative to the strip with the maximum signal,
125    *     and its spatial resolution (in cm), applying PFAtt.
126    *     ang is the effective angle, relative to the sensor
127    *---------------------------------------------------------------
128    
129          character*4 PFAtt
130          include 'commontracker.f'
131          include 'level1.f'
132    
133          corr = 0
134          res  = 0
135    
136          if(ic.le.0)return
137    
138          iview   = VIEW(ic)
139    
140          if(mod(iview,2).eq.0)then
141    c     =================================================
142    c     X view
143    c     =================================================
144    
145             res = RESXAV
146    
147             if(PFAtt.eq.'COG1')then
148    
149                corr   = 0
150                res = 1e-4*pitchX/sqrt(12.)!!res
151    
152             elseif(PFAtt.eq.'COG2')then
153    
154                corr    = cog(2,ic)            
155                res = risx_cog(abs(ang))!TEMPORANEO              
156                res = res*fbad_cog(2,ic)
157    
158             elseif(PFAtt.eq.'COG3')then
159    
160                corr    = cog(3,ic)            
161                res = risx_cog(abs(ang))!TEMPORANEO                      
162                res = res*fbad_cog(3,ic)
163    
164             elseif(PFAtt.eq.'COG4')then
165    
166                corr    = cog(4,ic)            
167                res = risx_cog(abs(ang))!TEMPORANEO                      
168                res = res*fbad_cog(4,ic)
169    
170             elseif(PFAtt.eq.'ETA2')then
171    
172                corr  = pfaeta2(ic,ang)          
173                res = risxeta2(abs(ang))
174                res = res*fbad_cog(2,ic)
175    
176             elseif(PFAtt.eq.'ETA3')then                        
177    
178        integer function npfastrips(ic,PFA,angle)              corr  = pfaeta3(ic,ang)          
179                res = risxeta3(abs(ang))                      
180                res = res*fbad_cog(3,ic)              
181    
182             elseif(PFAtt.eq.'ETA4')then                        
183    
184                corr  = pfaeta4(ic,ang)            
185                res = risxeta4(abs(ang))                      
186                res = res*fbad_cog(4,ic)              
187    
188             elseif(PFAtt.eq.'ETA')then  
189    
190                corr  = pfaeta(ic,ang)            
191    c            res = riseta(ic,ang)                    
192                res = riseta(iview,ang)                    
193                res = res*fbad_eta(ic,ang)            
194    
195             elseif(PFAtt.eq.'ETAL')then  
196    
197                corr  = pfaetal(ic,ang)            
198                res = riseta(iview,ang)                    
199                res = res*fbad_eta(ic,ang)            
200    
201             elseif(PFAtt.eq.'COG')then          
202    
203                corr  = cog(0,ic)            
204                res = risx_cog(abs(ang))                    
205                res = res*fbad_cog(0,ic)
206    
207             else
208                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
209             endif
210    
211    
212    *     ======================================
213    *     temporary patch for saturated clusters
214    *     ======================================
215             if( nsatstrips(ic).gt.0 )then
216                corr  = cog(4,ic)            
217                res = pitchX*1e-4/sqrt(12.)
218    cc            cc=cog(4,ic)
219    c$$$            print*,ic,' *** ',cc
220    c$$$            print*,ic,' *** ',res
221             endif
222    
223    
224          elseif(mod(iview,2).eq.1)then
225    c     =================================================
226    c     Y view
227    c     =================================================
228    
229             res = RESYAV
230    
231             if(PFAtt.eq.'COG1')then
232    
233                corr  = 0  
234                res = 1e-4*pitchY/sqrt(12.)!res  
235    
236             elseif(PFAtt.eq.'COG2')then
237    
238                corr    = cog(2,ic)
239                res = risy_cog(abs(ang))!TEMPORANEO
240                res = res*fbad_cog(2,ic)
241    
242             elseif(PFAtt.eq.'COG3')then
243    
244                corr    = cog(3,ic)
245                res = risy_cog(abs(ang))!TEMPORANEO
246                res = res*fbad_cog(3,ic)
247    
248             elseif(PFAtt.eq.'COG4')then
249    
250                corr    = cog(4,ic)
251                res = risy_cog(abs(ang))!TEMPORANEO
252                res = res*fbad_cog(4,ic)
253    
254             elseif(PFAtt.eq.'ETA2')then
255    
256                corr    = pfaeta2(ic,ang)          
257                res = risyeta2(abs(ang))              
258                res = res*fbad_cog(2,ic)
259    
260             elseif(PFAtt.eq.'ETA3')then                      
261    
262                corr    = pfaeta3(ic,ang)
263                res = res*fbad_cog(3,ic)  
264    
265             elseif(PFAtt.eq.'ETA4')then  
266    
267                corr    = pfaeta4(ic,ang)
268                res = res*fbad_cog(4,ic)
269    
270             elseif(PFAtt.eq.'ETA')then
271    
272                corr    = pfaeta(ic,ang)
273    c            res = riseta(ic,ang)  
274                res = riseta(iview,ang)  
275                res = res*fbad_eta(ic,ang)
276    
277             elseif(PFAtt.eq.'ETAL')then
278    
279                corr    = pfaetal(ic,ang)
280                res = riseta(iview,ang)  
281                res = res*fbad_eta(ic,ang)
282    
283             elseif(PFAtt.eq.'COG')then
284    
285                corr    = cog(0,ic)            
286                res = risy_cog(abs(ang))
287                res = res*fbad_cog(0,ic)
288    
289             else
290                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
291             endif
292    
293    
294    *     ======================================
295    *     temporary patch for saturated clusters
296    *     ======================================
297             if( nsatstrips(ic).gt.0 )then
298                corr    = cog(4,ic)            
299                res = pitchY*1e-4/sqrt(12.)
300    cc            cc=cog(4,ic)
301    c$$$            print*,ic,' *** ',cc
302    c$$$            print*,ic,' *** ',res
303             endif
304            
305          endif
306          end
307    
308    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
309          integer function npfastrips(ic,angle)
310  *--------------------------------------------------------------  *--------------------------------------------------------------
311  *     thid function returns the number of strips used  *     thid function returns the number of strips used
312  *     to evaluate the position of a cluster, according to the p.f.a.  *     to evaluate the position of a cluster, according to the p.f.a.
# Line 35  Line 315 
315        include 'level1.f'        include 'level1.f'
316        include 'calib.f'        include 'calib.f'
317    
318        character*4 usedPFA,PFA        character*4 usedPFA
319          
320    
321    
322        usedPFA=PFA        call idtoc(pfaid,usedPFA)
323    
324        npfastrips=0        npfastrips=-1
325    
326        if(usedPFA.eq.'COG1')npfastrips=1        if(usedPFA.eq.'COG1')npfastrips=1
327        if(usedPFA.eq.'COG2')npfastrips=2        if(usedPFA.eq.'COG2')npfastrips=2
# Line 50  Line 331 
331        if(usedPFA.eq.'ETA3')npfastrips=3        if(usedPFA.eq.'ETA3')npfastrips=3
332        if(usedPFA.eq.'ETA4')npfastrips=4        if(usedPFA.eq.'ETA4')npfastrips=4
333  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
334        if(usedPFA.eq.'ETA')then        if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
335  c         print*,VIEW(ic),angle  c         print*,VIEW(ic),angle
336           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
337              if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then              if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
# Line 60  c         print*,VIEW(ic),angle Line 341  c         print*,VIEW(ic),angle
341              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
342                 npfastrips=4                 npfastrips=4
343              else              else
344                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
345              endif                                      endif                        
346           else                   !X-view           else                   !X-view
347              if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then              if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
# Line 71  c               usedPFA='COG' Line 351  c               usedPFA='COG'
351              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
352                 npfastrips=4                 npfastrips=4
353              else              else
354                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
355              endif                                      endif                        
356           endif           endif
357        endif        endif
358  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
359        if(usedPFA.eq.'COG')then        if(usedPFA.eq.'COG')then
360    
361           iv=VIEW(ic)           npfastrips=0
362           if(mod(iv,2).eq.1)incut=incuty  
363           if(mod(iv,2).eq.0)incut=incutx  c$$$         iv=VIEW(ic)
364           istart = INDSTART(IC)  c$$$         if(mod(iv,2).eq.1)incut=incuty
365           istop  = TOTCLLENGTH  c$$$         if(mod(iv,2).eq.0)incut=incutx
366           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1  c$$$         istart = INDSTART(IC)
367           mu  = 0  c$$$         istop  = TOTCLLENGTH
368           do i = INDMAX(IC),istart,-1  c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
369              ipos = i-INDMAX(ic)  c$$$         mu  = 0
370              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC),istart,-1
371              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
372                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
373                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
374              else  c$$$               mu = mu + 1
375                 goto 10  c$$$               print*,i,mu
376              endif  c$$$            else
377           enddo  c$$$               goto 10
378   10      continue  c$$$            endif
379           do i = INDMAX(IC)+1,istop  c$$$         enddo
380              ipos = i-INDMAX(ic)  c$$$ 10      continue
381              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC)+1,istop
382              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
383                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
384                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
385              else  c$$$               mu = mu + 1
386                 goto 20  c$$$               print*,i,mu
387              endif  c$$$            else
388           enddo  c$$$               goto 20
389   20      continue  c$$$            endif
390           npfastrips=mu  c$$$         enddo
391    c$$$ 20      continue
392    c$$$         npfastrips=mu
393    
394        endif        endif
395  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
396    
397  c      print*,pfastrips  c      print*,pfaid,usedPFA,angle,npfastrips
398    
399        return        return
400        end        end
# Line 136  c      print*,pfastrips Line 417  c      print*,pfastrips
417    
418        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
419                
420           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
421              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
422           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then  cc            print*,pfaeta2(ic,angle)
423             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
424              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
425           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
426              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
427           else           else
428              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 148  c      print*,pfastrips Line 430  c      print*,pfastrips
430    
431        else                      !X-view        else                      !X-view
432    
433           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
434              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
435           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
436              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
437           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
438              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
439           else           else
440              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 164  c      print*,pfastrips Line 446  c      print*,pfastrips
446        end        end
447    
448  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
449          real function pfaetal(ic,angle)
450    *--------------------------------------------------------------
451    *     this function returns the position (in strip units)
452    *     it calls:
453    *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
454    *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
455    *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
456    *     according to the angle
457    *--------------------------------------------------------------
458          include 'commontracker.f'
459          include 'level1.f'
460          include 'calib.f'
461          
462          pfaetal = 0
463    
464          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
465          
466             if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
467                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
468    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
469             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
470                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
471             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
472                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
473             else
474                pfaetal = cog(4,ic)
475             endif            
476    
477          else                      !X-view
478    
479             if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
480                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
481    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
482             elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
483                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
484             elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
485                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
486             else
487                pfaetal = cog(4,ic)
488             endif            
489                
490          endif
491          
492     100  return
493          end
494    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
495  c      real function riseta(ic,angle)  c      real function riseta(ic,angle)
496        real function riseta(iview,angle)        real function riseta(iview,angle)
497  *--------------------------------------------------------------  *--------------------------------------------------------------
# Line 285  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 613  c      if(mod(int(VIEW(ic)),2).eq.1)then
613        cog2 = cog(2,ic)                  cog2 = cog(2,ic)          
614        pfaeta2=cog2        pfaeta2=cog2
615    
616    *     ----------------
617  *     find angular bin  *     find angular bin
618    *     ----------------
619  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
620        do iang=1,nangbin        do iang=1,nangbin
621           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
# Line 293  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 623  c      if(mod(int(VIEW(ic)),2).eq.1)then
623              goto 98              goto 98
624           endif           endif
625        enddo        enddo
626        if(DEBUG)        if(DEBUG.EQ.1)
627       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
628        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
629        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
630   98   continue                  !jump here if ok   98   continue                  !jump here if ok
631    
632    *     -------------
633    *     within +/-0.5
634    *     -------------
635    
636  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  
   
637        iadd=0        iadd=0
638   10   continue   10   continue
639        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
640           cog2 = cog2 + 1           cog2 = cog2 + 1
641           iadd = iadd + 1           iadd = iadd + 1
642             if(iadd>iaddmax)goto 111
643           goto 10           goto 10
644        endif        endif
645   20   continue   20   continue
646        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
647           cog2 = cog2 - 1           cog2 = cog2 - 1
648           iadd = iadd - 1           iadd = iadd - 1
649             if(iadd<-1*iaddmax)goto 111
650           goto 20           goto 20
651        endif        endif
652          goto 1111
653     111  continue
654          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
655          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
656          cog2=0
657     1111 continue
658    
659  *     --------------------------------  *     --------------------------------
660  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 369  c$$$         pfaeta2=pfaeta2+1.   !temp Line 689  c$$$         pfaeta2=pfaeta2+1.   !temp
689  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
690  c$$$      endif  c$$$      endif
691    
692        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
693       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
694    
695    
# Line 399  c$$$      endif Line 719  c$$$      endif
719    
720        iview = VIEW(ic)                    iview = VIEW(ic)            
721        lad = nld(MAXS(ic),VIEW(ic))        lad = nld(MAXS(ic),VIEW(ic))
722        cog3 = cog(3,ic)                    cog3 = cog(3,ic)
723          cc = cog3
724          cog3 = cc
725        pfaeta3=cog3        pfaeta3=cog3
726    
727    *     ----------------
728  *     find angular bin  *     find angular bin
729    *     ----------------
730  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
731        do iang=1,nangbin        do iang=1,nangbin
732  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 411  c         print*,'~~~~~~~~~~~~ ',iang,an Line 735  c         print*,'~~~~~~~~~~~~ ',iang,an
735              goto 98              goto 98
736           endif           endif
737        enddo        enddo
738        if(DEBUG)        if(DEBUG.EQ.1)
739       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
740        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
741        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
742   98   continue                  !jump here if ok   98   continue                  !jump here if ok
743    
744    *     -------------
745    *     within +/-0.5
746    *     -------------
747    
748  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  
   
         
749        iadd=0        iadd=0
750   10   continue   10   continue
751        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
752           cog3 = cog3 + 1           cog3   = cog3 + 1.
753           iadd = iadd + 1           iadd = iadd + 1
754             if(iadd>iaddmax) goto 111
755           goto 10           goto 10
756        endif        endif
757   20   continue   20   continue
758        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
759           cog3 = cog3 - 1           cog3 = cog3 - 1.
760           iadd = iadd - 1           iadd = iadd - 1
761             if(iadd<-1*iaddmax) goto 111
762           goto 20           goto 20
763        endif        endif
764          goto 1111
765     111  continue
766          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
767          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
768          cog3=0      
769     1111 continue
770    
771  *     --------------------------------  *     --------------------------------
772  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 477  c            print*,'-----',x1,x2,y1,y2 Line 792  c            print*,'-----',x1,x2,y1,y2
792        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
793        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
794    
 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  
795    
796        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
797       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
798    
799   100  return   100  return
# Line 518  c$$$      endif Line 825  c$$$      endif
825        cog4=cog(4,ic)        cog4=cog(4,ic)
826        pfaeta4=cog4        pfaeta4=cog4
827    
828    *     ----------------
829  *     find angular bin  *     find angular bin
830    *     ----------------
831  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
832        do iang=1,nangbin        do iang=1,nangbin
833  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 527  c         print*,'~~~~~~~~~~~~ ',iang,an Line 836  c         print*,'~~~~~~~~~~~~ ',iang,an
836              goto 98              goto 98
837           endif           endif
838        enddo        enddo
839        if(DEBUG)        if(DEBUG.EQ.1)
840       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
841        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
842        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
843   98   continue                  !jump here if ok   98   continue                  !jump here if ok
844    
845    *     -------------
846    *     within +/-0.5
847    *     -------------
848    
849  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  
   
         
850        iadd=0        iadd=0
851   10   continue   10   continue
852        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
853           cog4 = cog4 + 1           cog4 = cog4 + 1
854           iadd = iadd + 1           iadd = iadd + 1
855             if(iadd>iaddmax)goto 111
856           goto 10           goto 10
857        endif        endif
858   20   continue   20   continue
859        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
860           cog4 = cog4 - 1           cog4 = cog4 - 1
861           iadd = iadd - 1           iadd = iadd - 1
862             if(iadd<-1*iaddmax)goto 111
863           goto 20           goto 20
864        endif        endif
865          goto 1111
866     111  continue
867          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
868          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
869          cog4=0
870     1111 continue
871    
872  *     --------------------------------  *     --------------------------------
873  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 602  c$$$         pfaeta2=pfaeta2+1.   !temp Line 902  c$$$         pfaeta2=pfaeta2+1.   !temp
902  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
903  c$$$      endif  c$$$      endif
904    
905        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
906       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
907    
908   100  return   100  return
# Line 611  c$$$      endif Line 911  c$$$      endif
911    
912    
913  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c$$$      real function cog0(ncog,ic)  
 c$$$*-------------------------------------------------  
 c$$$*     this function returns  
 c$$$*  
 c$$$*     - the Center-Of-Gravity of the cluster IC  
 c$$$*     evaluated using NCOG strips,  
 c$$$*     calculated relative to MAXS(IC)  
 c$$$*      
 c$$$*     - zero in case that not  enough strips  
 c$$$*     have a positive signal  
 c$$$*      
 c$$$*     NOTE:  
 c$$$*     This is the old definition, used by Straulino.  
 c$$$*     The new routine, according to Landi,  
 c$$$*     is COG(NCOG,IC)  
 c$$$*-------------------------------------------------  
 c$$$  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$        
 c$$$*     --> signal of the central strip  
 c$$$      sc = CLSIGNAL(INDMAX(ic)) !center  
 c$$$  
 c$$$*     signal of adjacent strips  
 c$$$*     --> left  
 c$$$      sl1 = 0                  !left 1  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
 c$$$  
 c$$$      sl2 = 0                  !left 2  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
 c$$$  
 c$$$*     --> right  
 c$$$      sr1 = 0                  !right 1  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
 c$$$  
 c$$$      sr2 = 0                  !right 2  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
 c$$$        
 c$$$************************************************************  
 c$$$*     COG computation  
 c$$$************************************************************  
 c$$$  
 c$$$c      print*,sl2,sl1,sc,sr1,sr2  
 c$$$  
 c$$$      COG = 0.  
 c$$$          
 c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then  
 c$$$          
 c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
 c$$$            COG = -sl1/(sl1+sc)          
 c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
 c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)  
 c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
 c$$$            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  
 c$$$         else  
 c$$$            COG = 0.  
 c$$$         endif  
 c$$$          
 c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then  
 c$$$  
 c$$$         if(ncog.eq.2.and.sr1.ne.0)then  
 c$$$            COG = sr1/(sc+sr1)              
 c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
 c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)  
 c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
 c$$$            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  
 c$$$         else  
 c$$$            COG = 0.  
 c$$$         endif  
 c$$$  
 c$$$      endif  
 c$$$  
 c$$$      COG0 = COG  
 c$$$  
 c$$$c      print *,ncog,ic,cog,'/////////////'  
 c$$$  
 c$$$      return  
 c$$$      end  
   
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
914        real function cog(ncog,ic)        real function cog(ncog,ic)
915  *-------------------------------------------------  *-------------------------------------------------
916  *     this function returns  *     this function returns
# Line 767  c$$$      end Line 970  c$$$      end
970                    
971           COG = 0.           COG = 0.
972                    
973  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
974    
975  c     ==============================================================  c     ==============================================================
976           if(ncog.eq.1)then           if(ncog.eq.1)then
977              COG = 0.              COG = 0.
978              if(sr1.gt.sc)cog=1. !NEW              if(sr1.gt.sc)cog=1.
979              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
980  c     ==============================================================  c     ==============================================================
981           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
982                COG = 0.
983              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
984                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
985              elseif(sl1.lt.sr1)then              elseif(sl1.lt.sr1)then
986                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                                         if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
987              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
988                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
989       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
990                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
991       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
992              endif              endif
993  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
994  c     $           ,' : ',sl2,sl1,sc,sr1,sr2  c     $           ,' : ',sl2,sl1,sc,sr1,sr2
995  c     ==============================================================  c     ==============================================================
996           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
997               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
998  c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)              sss = sc
999                if( sl1.ne.-9999. )COG = COG-sl1
1000                if( sl1.ne.-9999. )sss = sss+sl1
1001                if( sr1.ne.-9999. )COG = COG+sr1
1002                if( sr1.ne.-9999. )sss = sss+sr1
1003                if(sss.ne.0)COG=COG/sss
1004    
1005    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1006    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1007  c     $            ,' : ',sl2,sl1,sc,sr1,sr2  c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1008  c     ==============================================================  c     ==============================================================
1009           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1010    
1011                COG = 0
1012                sss = sc
1013                if( sl1.ne.-9999. )COG = COG-sl1
1014                if( sl1.ne.-9999. )sss = sss+sl1
1015                if( sr1.ne.-9999. )COG = COG+sr1
1016                if( sr1.ne.-9999. )sss = sss+sr1
1017              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1018                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1019       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1020              elseif(sl2.lt.sr2)then              elseif(sl2.lt.sr2)then
1021                 if((sr2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1022       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1023              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW              elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1024                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1025       $              .and.(sl2+sl1+sc+sr1).ne.0 )       $              .and.(sl2+sss).ne.0 )
1026       $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW       $              cog = (cog-2*sl2)/(sl2+sss)
1027                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)                 if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1028       $              .and.(sr2+sl1+sc+sr1).ne.0 )       $              .and.(sr2+sss).ne.0 )
1029       $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW                     $              cog = (2*sr2+cog)/(sr2+sss)              
1030              endif              endif
1031    c     ==============================================================
1032             elseif(ncog.eq.5)then
1033                COG = 0
1034                sss = sc
1035                if( sl1.ne.-9999. )COG = COG-sl1
1036                if( sl1.ne.-9999. )sss = sss+sl1
1037                if( sr1.ne.-9999. )COG = COG+sr1
1038                if( sr1.ne.-9999. )sss = sss+sr1
1039                if( sl2.ne.-9999. )COG = COG-2*sl2
1040                if( sl2.ne.-9999. )sss = sss+sl2
1041                if( sr2.ne.-9999. )COG = COG+2*sr2
1042                if( sr2.ne.-9999. )sss = sss+sr2
1043                if(sss.ne.0)COG=COG/sss
1044           else           else
1045              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1046       $           ,' not implemented'       $           ,' not implemented'
# Line 879  c         print*,'-------' Line 1111  c         print*,'-------'
1111    
1112  c      print *,'## cog ',ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1113    
1114          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1115             if(DEBUG.eq.1)
1116         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1117             if(DEBUG.eq.1)
1118         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1119          endif
1120    
1121    
1122        return        return
1123        end        end
1124    
# Line 1024  c            COG = 0. Line 1264  c            COG = 0.
1264           SGN = 0.           SGN = 0.
1265           SNU = 0.           SNU = 0.
1266           SDE = 0.           SDE = 0.
 c$$$         do i=INDMAX(IC),istart,-1  
 c$$$            ipos = i-INDMAX(ic)  
 c$$$            cut  = incut*CLSIGMA(i)  
 c$$$            if(CLSIGNAL(i).gt.cut)then  
 c$$$               COG = COG + ipos*CLSIGNAL(i)  
 c$$$               SGN = SGN + CLSIGNAL(i)  
 c$$$            else  
 c$$$               goto 10  
 c$$$            endif  
 c$$$         enddo  
 c$$$ 10      continue  
 c$$$         do i=INDMAX(IC)+1,istop  
 c$$$            ipos = i-INDMAX(ic)  
 c$$$            cut  = incut*CLSIGMA(i)  
 c$$$            if(CLSIGNAL(i).gt.cut)then  
 c$$$               COG = COG + ipos*CLSIGNAL(i)  
 c$$$               SGN = SGN + CLSIGNAL(i)  
 c$$$            else  
 c$$$               goto 20  
 c$$$            endif  
 c$$$         enddo  
 c$$$ 20      continue  
 c$$$         if(SGN.le.0)then  
 c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN  
 c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)  
 c$$$            print*,(CLSIGNAL(i),i=istart,istop)  
 c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '  
 c$$$         else  
 c$$$            COG=COG/SGN  
 c$$$         endif  
1267    
1268           do i=INDMAX(IC),istart,-1           do i=INDMAX(IC),istart,-1
1269              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
# Line 1101  c$$$         endif Line 1311  c$$$         endif
1311        end        end
1312    
1313    
1314  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1315        real function fbad_cog0(ncog,ic)  c$$$      real function fbad_cog0(ncog,ic)
1316  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1317  *     this function returns a factor that takes into  c$$$*     this function returns a factor that takes into
1318  *     account deterioration of the spatial resolution  c$$$*     account deterioration of the spatial resolution
1319  *     in the case BAD strips are included in the cluster.  c$$$*     in the case BAD strips are included in the cluster.
1320  *     This factor should multiply the nominal spatial  c$$$*     This factor should multiply the nominal spatial
1321  *     resolution.  c$$$*     resolution.
1322  *  c$$$*
1323  *     NB!!!  c$$$*     NB!!!
1324  *     (this is the old version. It consider only the two  c$$$*     (this is the old version. It consider only the two
1325  *     strips with the greatest signal. The new one is  c$$$*     strips with the greatest signal. The new one is
1326  *     fbad_cog(ncog,ic) )  c$$$*     fbad_cog(ncog,ic) )
1327  *      c$$$*    
1328  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1329    c$$$
1330        include 'commontracker.f'  c$$$      include 'commontracker.f'
1331        include 'level1.f'  c$$$      include 'level1.f'
1332        include 'calib.f'  c$$$      include 'calib.f'
1333    c$$$
1334  *     --> signal of the central strip  c$$$*     --> signal of the central strip
1335        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
1336    c$$$
1337  *     signal of adjacent strips  c$$$*     signal of adjacent strips
1338  *     --> left  c$$$*     --> left
1339        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
1340        if(  c$$$      if(
1341       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
1342       $     )  c$$$     $     )
1343       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
1344    c$$$
1345        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
1346        if(  c$$$      if(
1347       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
1348       $     )  c$$$     $     )
1349       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
1350    c$$$
1351  *     --> right  c$$$*     --> right
1352        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
1353        if(  c$$$      if(
1354       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1355       $     .or.  c$$$     $     .or.
1356       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1357       $     )  c$$$     $     )
1358       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
1359    c$$$
1360        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
1361        if(  c$$$      if(
1362       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1363       $     .or.  c$$$     $     .or.
1364       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1365       $     )  c$$$     $     )
1366       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
1367    c$$$
1368    c$$$
1369        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1370           f  = 4.  c$$$         f  = 4.
1371           si = 8.4  c$$$         si = 8.4
1372        else                              !X-view  c$$$      else                              !X-view
1373           f  = 6.  c$$$         f  = 6.
1374           si = 3.9  c$$$         si = 3.9
1375        endif  c$$$      endif
1376    c$$$
1377        fbad_cog = 1.  c$$$      fbad_cog = 1.
1378        f0 = 1  c$$$      f0 = 1
1379        f1 = 1  c$$$      f1 = 1
1380        f2 = 1  c$$$      f2 = 1
1381        f3 = 1    c$$$      f3 = 1  
1382        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
1383            c$$$        
1384           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
1385           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
1386  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
1387    c$$$
1388           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
1389              fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
1390           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
1391              fbad_cog = 1.  c$$$            fbad_cog = 1.
1392           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
1393              fbad_cog = 1.  c$$$            fbad_cog = 1.
1394           else  c$$$         else
1395              fbad_cog = 1.  c$$$            fbad_cog = 1.
1396           endif  c$$$         endif
1397            c$$$        
1398        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
1399    c$$$
1400    c$$$
1401           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
1402           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
1403  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
1404    c$$$
1405           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
1406              fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)  c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
1407           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
1408              fbad_cog = 1.  c$$$            fbad_cog = 1.
1409           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
1410              fbad_cog = 1.  c$$$            fbad_cog = 1.
1411           else  c$$$         else
1412              fbad_cog = 1.  c$$$            fbad_cog = 1.
1413           endif  c$$$         endif
1414    c$$$
1415        endif  c$$$      endif
1416    c$$$
1417        fbad_cog0 = sqrt(fbad_cog)  c$$$      fbad_cog0 = sqrt(fbad_cog)
1418    c$$$
1419        return  c$$$      return
1420        end  c$$$      end
1421    c$$$
1422    c$$$
1423    c$$$
1424    
1425  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1426    
# Line 1705  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1915  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1915        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1916    
1917        END        END
1918    
1919    
1920    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1921          real function pfacorr(ic,angle)
1922    *--------------------------------------------------------------
1923    *     this function returns the landi correction for this cluster
1924    *--------------------------------------------------------------
1925          include 'commontracker.f'
1926          include 'calib.f'
1927          include 'level1.f'
1928    
1929          real angle
1930          integer iview,lad
1931    
1932          iview = VIEW(ic)            
1933          lad = nld(MAXS(ic),VIEW(ic))
1934    
1935    *     find angular bin
1936    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1937          do iang=1,nangbin
1938             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1939                iangle=iang
1940                goto 98
1941             endif
1942          enddo
1943          if(DEBUG.eq.1)
1944         $     print*,'pfacorr *** warning *** angle out of range: ',angle
1945          if(angle.le.angL(1))iang=1
1946          if(angle.ge.angR(nangbin))iang=nangbin
1947     98   continue                  !jump here if ok
1948    
1949          pfacorr = fcorr(iview,lad,iang)
1950    
1951          if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr
1952    
1953    
1954     100  return
1955          end
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23