/[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.11 by pam-fi, Tue May 15 16:22:18 2007 UTC revision 1.22 by pam-fi, Tue Sep 4 09:47:49 2007 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             tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001
74    
75          elseif(mod(iview,2).eq.1)then
76    c     =================================================
77    c     Y view
78    c     =================================================        
79    c     here bbb is the x component of the m.filed
80             angy = ang
81             bx   = bbb
82             tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001        
83    
84          endif      
85          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
86    
87          return
88          end
89  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
90  *     this file contains all subroutines and functions        real function fieldcorr(iview,bbb)
91  *     that are needed for position finding algorithms  
92  *            include 'commontracker.f'
93  *  
94          fieldcorr = 0.
95    
96          if(mod(iview,2).eq.0)then
97    
98    c     =================================================
99    c     X view
100    c     =================================================
101    c     here bbb is the y component of the m.field
102             by   = bbb
103             if(iview.eq.12) by = -1. * bbb
104             fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
105    
106          elseif(mod(iview,2).eq.1)then
107    c     =================================================
108    c     Y view
109    c     =================================================        
110    c     here bbb is the x component of the m.filed
111             bx   = bbb
112             fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
113    
114          endif      
115          
116          return
117          end
118  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
119    
120          subroutine applypfa(PFAtt,ic,ang,corr,res)
121    *---------------------------------------------------------------
122    *     this subroutine calculate the coordinate of cluster ic (in
123    *     strip units), relative to the strip with the maximum signal,
124    *     and its spatial resolution (in cm), applying PFAtt.
125    *     ang is the effective angle, relative to the sensor
126    *---------------------------------------------------------------
127    
128          character*4 PFAtt
129          include 'commontracker.f'
130          include 'level1.f'
131    
132          corr = 0
133          res  = 0
134    
135          if(ic.le.0)return
136    
137          iview   = VIEW(ic)
138    
139          if(mod(iview,2).eq.0)then
140    c     =================================================
141    c     X view
142    c     =================================================
143    
144             res = RESXAV
145    
146             if(PFAtt.eq.'COG1')then
147    
148                corr   = 0
149                res = 1e-4*pitchX/sqrt(12.)!!res
150    
151             elseif(PFAtt.eq.'COG2')then
152    
153                corr    = cog(2,ic)            
154                res = risx_cog(abs(ang))!TEMPORANEO              
155                res = res*fbad_cog(2,ic)
156    
157             elseif(PFAtt.eq.'COG3')then
158    
159                corr    = cog(3,ic)            
160                res = risx_cog(abs(ang))!TEMPORANEO                      
161                res = res*fbad_cog(3,ic)
162    
163             elseif(PFAtt.eq.'COG4')then
164    
165                corr    = cog(4,ic)            
166                res = risx_cog(abs(ang))!TEMPORANEO                      
167                res = res*fbad_cog(4,ic)
168    
169             elseif(PFAtt.eq.'ETA2')then
170    
171                corr  = pfaeta2(ic,ang)          
172                res = risxeta2(abs(ang))
173                res = res*fbad_cog(2,ic)
174    
175             elseif(PFAtt.eq.'ETA3')then                        
176    
177                corr  = pfaeta3(ic,ang)          
178                res = risxeta3(abs(ang))                      
179                res = res*fbad_cog(3,ic)              
180    
181             elseif(PFAtt.eq.'ETA4')then                        
182    
183                corr  = pfaeta4(ic,ang)            
184                res = risxeta4(abs(ang))                      
185                res = res*fbad_cog(4,ic)              
186    
187             elseif(PFAtt.eq.'ETA')then  
188    
189                corr  = pfaeta(ic,ang)            
190    c            res = riseta(ic,ang)                    
191                res = riseta(iview,ang)                    
192                res = res*fbad_eta(ic,ang)            
193    
194             elseif(PFAtt.eq.'ETAL')then  
195    
196                corr  = pfaetal(ic,ang)            
197                res = riseta(iview,ang)                    
198                res = res*fbad_eta(ic,ang)            
199    
200             elseif(PFAtt.eq.'COG')then          
201    
202                corr  = cog(0,ic)            
203                res = risx_cog(abs(ang))                    
204                res = res*fbad_cog(0,ic)
205    
206             else
207                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
208             endif
209    
210    
211    *     ======================================
212    *     temporary patch for saturated clusters
213    *     ======================================
214             if( nsatstrips(ic).gt.0 )then
215                corr  = cog(4,ic)            
216                res = pitchX*1e-4/sqrt(12.)
217    cc            cc=cog(4,ic)
218    c$$$            print*,ic,' *** ',cc
219    c$$$            print*,ic,' *** ',res
220             endif
221    
222    
223          elseif(mod(iview,2).eq.1)then
224    c     =================================================
225    c     Y view
226    c     =================================================
227    
228             res = RESYAV
229    
230             if(PFAtt.eq.'COG1')then
231    
232                corr  = 0  
233                res = 1e-4*pitchY/sqrt(12.)!res  
234    
235             elseif(PFAtt.eq.'COG2')then
236    
237                corr    = cog(2,ic)
238                res = risy_cog(abs(ang))!TEMPORANEO
239                res = res*fbad_cog(2,ic)
240    
241             elseif(PFAtt.eq.'COG3')then
242    
243                corr    = cog(3,ic)
244                res = risy_cog(abs(ang))!TEMPORANEO
245                res = res*fbad_cog(3,ic)
246    
247             elseif(PFAtt.eq.'COG4')then
248    
249                corr    = cog(4,ic)
250                res = risy_cog(abs(ang))!TEMPORANEO
251                res = res*fbad_cog(4,ic)
252    
253             elseif(PFAtt.eq.'ETA2')then
254    
255        integer function npfastrips(ic,PFA,angle)              corr    = pfaeta2(ic,ang)          
256                res = risyeta2(abs(ang))              
257                res = res*fbad_cog(2,ic)
258    
259             elseif(PFAtt.eq.'ETA3')then                      
260    
261                corr    = pfaeta3(ic,ang)
262                res = res*fbad_cog(3,ic)  
263    
264             elseif(PFAtt.eq.'ETA4')then  
265    
266                corr    = pfaeta4(ic,ang)
267                res = res*fbad_cog(4,ic)
268    
269             elseif(PFAtt.eq.'ETA')then
270    
271                corr    = pfaeta(ic,ang)
272    c            res = riseta(ic,ang)  
273                res = riseta(iview,ang)  
274                res = res*fbad_eta(ic,ang)
275    
276             elseif(PFAtt.eq.'ETAL')then
277    
278                corr    = pfaetal(ic,ang)
279                res = riseta(iview,ang)  
280                res = res*fbad_eta(ic,ang)
281    
282             elseif(PFAtt.eq.'COG')then
283    
284                corr    = cog(0,ic)            
285                res = risy_cog(abs(ang))
286                res = res*fbad_cog(0,ic)
287    
288             else
289                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
290             endif
291    
292    
293    *     ======================================
294    *     temporary patch for saturated clusters
295    *     ======================================
296             if( nsatstrips(ic).gt.0 )then
297                corr    = cog(4,ic)            
298                res = pitchY*1e-4/sqrt(12.)
299    cc            cc=cog(4,ic)
300    c$$$            print*,ic,' *** ',cc
301    c$$$            print*,ic,' *** ',res
302             endif
303            
304          endif
305          end
306    
307    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
308          integer function npfastrips(ic,angle)
309  *--------------------------------------------------------------  *--------------------------------------------------------------
310  *     thid function returns the number of strips used  *     thid function returns the number of strips used
311  *     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 314 
314        include 'level1.f'        include 'level1.f'
315        include 'calib.f'        include 'calib.f'
316    
317        character*4 usedPFA,PFA        character*4 usedPFA
318          
319    
320    
321        usedPFA=PFA        call idtoc(pfaid,usedPFA)
322    
323        npfastrips=0        npfastrips=-1
324    
325        if(usedPFA.eq.'COG1')npfastrips=1        if(usedPFA.eq.'COG1')npfastrips=1
326        if(usedPFA.eq.'COG2')npfastrips=2        if(usedPFA.eq.'COG2')npfastrips=2
# Line 50  Line 330 
330        if(usedPFA.eq.'ETA3')npfastrips=3        if(usedPFA.eq.'ETA3')npfastrips=3
331        if(usedPFA.eq.'ETA4')npfastrips=4        if(usedPFA.eq.'ETA4')npfastrips=4
332  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
333        if(usedPFA.eq.'ETA')then        if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
334  c         print*,VIEW(ic),angle  c         print*,VIEW(ic),angle
335           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
336              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 340  c         print*,VIEW(ic),angle
340              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
341                 npfastrips=4                 npfastrips=4
342              else              else
343                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
344              endif                                      endif                        
345           else                   !X-view           else                   !X-view
346              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 350  c               usedPFA='COG'
350              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
351                 npfastrips=4                 npfastrips=4
352              else              else
353                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
354              endif                                      endif                        
355           endif           endif
356        endif        endif
357  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
358        if(usedPFA.eq.'COG')then        if(usedPFA.eq.'COG')then
359    
360           iv=VIEW(ic)           npfastrips=0
361           if(mod(iv,2).eq.1)incut=incuty  
362           if(mod(iv,2).eq.0)incut=incutx  c$$$         iv=VIEW(ic)
363           istart = INDSTART(IC)  c$$$         if(mod(iv,2).eq.1)incut=incuty
364           istop  = TOTCLLENGTH  c$$$         if(mod(iv,2).eq.0)incut=incutx
365           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1  c$$$         istart = INDSTART(IC)
366           mu  = 0  c$$$         istop  = TOTCLLENGTH
367           do i = INDMAX(IC),istart,-1  c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
368              ipos = i-INDMAX(ic)  c$$$         mu  = 0
369              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC),istart,-1
370              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
371                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
372                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
373              else  c$$$               mu = mu + 1
374                 goto 10  c$$$               print*,i,mu
375              endif  c$$$            else
376           enddo  c$$$               goto 10
377   10      continue  c$$$            endif
378           do i = INDMAX(IC)+1,istop  c$$$         enddo
379              ipos = i-INDMAX(ic)  c$$$ 10      continue
380              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC)+1,istop
381              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
382                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
383                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
384              else  c$$$               mu = mu + 1
385                 goto 20  c$$$               print*,i,mu
386              endif  c$$$            else
387           enddo  c$$$               goto 20
388   20      continue  c$$$            endif
389           npfastrips=mu  c$$$         enddo
390    c$$$ 20      continue
391    c$$$         npfastrips=mu
392    
393        endif        endif
394  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
395    
396  c      print*,pfastrips  c      print*,pfaid,usedPFA,angle,npfastrips
397    
398        return        return
399        end        end
# Line 136  c      print*,pfastrips Line 416  c      print*,pfastrips
416    
417        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
418                
419           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
420              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
421           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then  cc            print*,pfaeta2(ic,angle)
422             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
423              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
424           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
425              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
426           else           else
427              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 148  c      print*,pfastrips Line 429  c      print*,pfastrips
429    
430        else                      !X-view        else                      !X-view
431    
432           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
433              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
434           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
435              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
436           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
437              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
438           else           else
439              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 164  c      print*,pfastrips Line 445  c      print*,pfastrips
445        end        end
446    
447  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
448          real function pfaetal(ic,angle)
449    *--------------------------------------------------------------
450    *     this function returns the position (in strip units)
451    *     it calls:
452    *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
453    *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
454    *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
455    *     according to the angle
456    *--------------------------------------------------------------
457          include 'commontracker.f'
458          include 'level1.f'
459          include 'calib.f'
460          
461          pfaetal = 0
462    
463          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
464          
465             if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
466                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
467    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
468             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
469                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
470             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
471                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
472             else
473                pfaetal = cog(4,ic)
474             endif            
475    
476          else                      !X-view
477    
478             if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
479                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
480    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
481             elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
482                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
483             elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
484                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
485             else
486                pfaetal = cog(4,ic)
487             endif            
488                
489          endif
490          
491     100  return
492          end
493    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
494  c      real function riseta(ic,angle)  c      real function riseta(ic,angle)
495        real function riseta(iview,angle)        real function riseta(iview,angle)
496  *--------------------------------------------------------------  *--------------------------------------------------------------
497  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
498  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
499  *     it calls:  *     it calls:
500  *     - risx_eta2(angle)  *     - risxeta2(angle)
501  *     - risy_eta2(angle)  *     - risyeta2(angle)
502  *     - risx_eta3(angle)  *     - risxeta3(angle)
503  *     - risx_eta4(angle)  *     - risxeta4(angle)
504  *     according to the angle  *     according to the angle
505  *--------------------------------------------------------------  *--------------------------------------------------------------
506        include 'commontracker.f'        include 'commontracker.f'
# Line 187  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 514  c      if(mod(int(VIEW(ic)),2).eq.1)then
514                
515    
516           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
517              riseta = risy_eta2(angle)              riseta = risyeta2(angle)
518           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
519              riseta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
520           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
# Line 199  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 526  c      if(mod(int(VIEW(ic)),2).eq.1)then
526        else                      !X-view        else                      !X-view
527    
528           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
529              riseta = risx_eta2(angle)              riseta = risxeta2(angle)
530           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
531              riseta = risx_eta3(angle)              riseta = risxeta3(angle)
532           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
533              riseta = risx_eta4(angle)              riseta = risxeta4(angle)
534           else           else
535              riseta = risx_cog(angle)              riseta = risx_cog(angle)
536           endif                       endif            
537                            
538        endif        endif
539    
       print*,'---- ',riseta,iview,angle  
540    
541   100  return   100  return
542        end        end
# Line 294  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 620  c      if(mod(int(VIEW(ic)),2).eq.1)then
620              goto 98              goto 98
621           endif           endif
622        enddo        enddo
623        if(DEBUG)        if(DEBUG.EQ.1)
624       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
625        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
626        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
627   98   continue                  !jump here if ok   98   continue                  !jump here if ok
628    
629    
# Line 370  c$$$         pfaeta2=pfaeta2+1.   !temp Line 696  c$$$         pfaeta2=pfaeta2+1.   !temp
696  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
697  c$$$      endif  c$$$      endif
698    
699        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
700       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
701    
702    
# Line 412  c         print*,'~~~~~~~~~~~~ ',iang,an Line 738  c         print*,'~~~~~~~~~~~~ ',iang,an
738              goto 98              goto 98
739           endif           endif
740        enddo        enddo
741        if(DEBUG)        if(DEBUG.EQ.1)
742       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
743        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
744        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
745   98   continue                  !jump here if ok   98   continue                  !jump here if ok
746    
747    
# Line 487  c$$$         pfaeta2=pfaeta2+1.   !temp Line 813  c$$$         pfaeta2=pfaeta2+1.   !temp
813  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
814  c$$$      endif  c$$$      endif
815    
816        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
817       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
818    
819   100  return   100  return
# Line 528  c         print*,'~~~~~~~~~~~~ ',iang,an Line 854  c         print*,'~~~~~~~~~~~~ ',iang,an
854              goto 98              goto 98
855           endif           endif
856        enddo        enddo
857        if(DEBUG)        if(DEBUG.EQ.1)
858       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
859        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
860        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
861   98   continue                  !jump here if ok   98   continue                  !jump here if ok
862    
863    
# Line 603  c$$$         pfaeta2=pfaeta2+1.   !temp Line 929  c$$$         pfaeta2=pfaeta2+1.   !temp
929  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
930  c$$$      endif  c$$$      endif
931    
932        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
933       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
934    
935   100  return   100  return
# Line 612  c$$$      endif Line 938  c$$$      endif
938    
939    
940  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
       real function cog0(ncog,ic)  
 *-------------------------------------------------  
 *     this function returns  
 *  
 *     - the Center-Of-Gravity of the cluster IC  
 *     evaluated using NCOG strips,  
 *     calculated relative to MAXS(IC)  
 *      
 *     - zero in case that not  enough strips  
 *     have a positive signal  
 *      
 *     NOTE:  
 *     This is the old definition, used by Straulino.  
 *     The new routine, according to Landi,  
 *     is COG(NCOG,IC)  
 *-------------------------------------------------  
   
   
       include 'commontracker.f'  
       include 'level1.f'  
         
 *     --> signal of the central strip  
       sc = CLSIGNAL(INDMAX(ic)) !center  
   
 *     signal of adjacent strips  
 *     --> left  
       sl1 = 0                  !left 1  
       if(  
      $     (INDMAX(ic)-1).ge.INDSTART(ic)  
      $     )  
      $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
   
       sl2 = 0                  !left 2  
       if(  
      $     (INDMAX(ic)-2).ge.INDSTART(ic)  
      $     )  
      $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
   
 *     --> right  
       sr1 = 0                  !right 1  
       if(  
      $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
      $     .or.  
      $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
      $     )  
      $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
   
       sr2 = 0                  !right 2  
       if(  
      $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
      $     .or.  
      $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
      $     )  
      $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
         
 ************************************************************  
 *     COG computation  
 ************************************************************  
   
 c      print*,sl2,sl1,sc,sr1,sr2  
   
       COG = 0.  
           
       if(sl1.gt.sr1.and.sl1.gt.0.)then  
           
          if(ncog.eq.2.and.sl1.ne.0)then  
             COG = -sl1/(sl1+sc)          
          elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
             COG = (sr1-sl1)/(sl1+sc+sr1)  
          elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
             COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  
          else  
             COG = 0.  
          endif  
           
       elseif(sl1.le.sr1.and.sr1.gt.0.)then  
   
          if(ncog.eq.2.and.sr1.ne.0)then  
             COG = sr1/(sc+sr1)              
          elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
             COG = (sr1-sl1)/(sl1+sc+sr1)  
          elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
             COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  
          else  
             COG = 0.  
          endif  
   
       endif  
   
       COG0 = COG  
   
 c      print *,ncog,ic,cog,'/////////////'  
   
       return  
       end  
   
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
941        real function cog(ncog,ic)        real function cog(ncog,ic)
942  *-------------------------------------------------  *-------------------------------------------------
943  *     this function returns  *     this function returns
# Line 738  c      print *,ncog,ic,cog,'//////////// Line 967  c      print *,ncog,ic,cog,'////////////
967  *     --> signal of the central strip  *     --> signal of the central strip
968           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
969  *     signal of adjacent strips  *     signal of adjacent strips
970           sl1 = 0                !left 1           sl1 = -9999.           !left 1
971           if(           if(
972       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
973       $        )       $        )
974       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
975                    
976           sl2 = 0                !left 2           sl2 = -9999.           !left 2
977           if(           if(
978       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
979       $        )       $        )
980       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
981                    
982           sr1 = 0                !right 1           sr1 = -9999.           !right 1
983           if(           if(
984       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
985       $        .or.       $        .or.
# Line 758  c      print *,ncog,ic,cog,'//////////// Line 987  c      print *,ncog,ic,cog,'////////////
987       $        )       $        )
988       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
989                    
990           sr2 = 0                !right 2           sr2 = -9999.           !right 2
991           if(           if(
992       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
993       $        .or.       $        .or.
# Line 770  c      print *,ncog,ic,cog,'//////////// Line 999  c      print *,ncog,ic,cog,'////////////
999                    
1000  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c         print*,'## ',sl2,sl1,sc,sr1,sr2
1001    
1002    c     ==============================================================
1003           if(ncog.eq.1)then           if(ncog.eq.1)then
1004              COG = 0.              COG = 0.
1005                if(sr1.gt.sc)cog=1.
1006                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1007    c     ==============================================================
1008           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
1009                COG = 0.
1010              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1011                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
1012              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
1013                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
1014              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1015                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1016         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1017                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1018         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1019                endif
1020    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1021    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
1022    c     ==============================================================
1023           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1024               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
1025                sss = sc
1026                if( sl1.ne.-9999. )COG = COG-sl1
1027                if( sl1.ne.-9999. )sss = sss+sl1
1028                if( sr1.ne.-9999. )COG = COG+sr1
1029                if( sr1.ne.-9999. )sss = sss+sr1
1030                if(sss.ne.0)COG=COG/sss
1031    
1032    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1033    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1034    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1035    c     ==============================================================
1036           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1037    
1038                COG = 0
1039                sss = sc
1040                if( sl1.ne.-9999. )COG = COG-sl1
1041                if( sl1.ne.-9999. )sss = sss+sl1
1042                if( sr1.ne.-9999. )COG = COG+sr1
1043                if( sr1.ne.-9999. )sss = sss+sr1
1044              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1045                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1046       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1047              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
1048                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1049       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1050                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1051                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1052         $              .and.(sl2+sss).ne.0 )
1053         $              cog = (cog-2*sl2)/(sl2+sss)
1054                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1055         $              .and.(sr2+sss).ne.0 )
1056         $              cog = (2*sr2+cog)/(sr2+sss)              
1057              endif              endif
1058    c     ==============================================================
1059             elseif(ncog.eq.5)then
1060                COG = 0
1061                sss = sc
1062                if( sl1.ne.-9999. )COG = COG-sl1
1063                if( sl1.ne.-9999. )sss = sss+sl1
1064                if( sr1.ne.-9999. )COG = COG+sr1
1065                if( sr1.ne.-9999. )sss = sss+sr1
1066                if( sl2.ne.-9999. )COG = COG-2*sl2
1067                if( sl2.ne.-9999. )sss = sss+sl2
1068                if( sr2.ne.-9999. )COG = COG+2*sr2
1069                if( sr2.ne.-9999. )sss = sss+sr2
1070                if(sss.ne.0)COG=COG/sss
1071           else           else
1072              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1073       $           ,' not implemented'       $           ,' not implemented'
# Line 818  c         print*,'-------' Line 1098  c         print*,'-------'
1098                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1099                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
1100                 mu = mu + 1                 mu = mu + 1
1101                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
1102              else              else
1103                 goto 10                 goto 10
1104              endif              endif
# Line 831  c         print*,'-------' Line 1111  c         print*,'-------'
1111                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1112                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
1113                 mu = mu + 1                 mu = mu + 1
1114                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
1115              else              else
1116                 goto 20                 goto 20
1117              endif              endif
1118           enddo           enddo
1119   20      continue   20      continue
1120           if(SGN.le.0)then           if(SGN.le.0)then
1121  c            print*,'cog(0,ic) --> ic, dedx ',ic,SGN              print*,'cog(0,ic) --> ic, dedx ',ic,SGN
1122              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1123              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
1124  c            print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
# Line 1003  c            COG = 0. Line 1283  c            COG = 0.
1283           SGN = 0.           SGN = 0.
1284           SNU = 0.           SNU = 0.
1285           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  
1286    
1287           do i=INDMAX(IC),istart,-1           do i=INDMAX(IC),istart,-1
1288              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
# Line 1080  c$$$         endif Line 1330  c$$$         endif
1330        end        end
1331    
1332    
1333  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1334        real function fbad_cog0(ncog,ic)  c$$$      real function fbad_cog0(ncog,ic)
1335  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1336  *     this function returns a factor that takes into  c$$$*     this function returns a factor that takes into
1337  *     account deterioration of the spatial resolution  c$$$*     account deterioration of the spatial resolution
1338  *     in the case BAD strips are included in the cluster.  c$$$*     in the case BAD strips are included in the cluster.
1339  *     This factor should multiply the nominal spatial  c$$$*     This factor should multiply the nominal spatial
1340  *     resolution.  c$$$*     resolution.
1341  *  c$$$*
1342  *     NB!!!  c$$$*     NB!!!
1343  *     (this is the old version. It consider only the two  c$$$*     (this is the old version. It consider only the two
1344  *     strips with the greatest signal. The new one is  c$$$*     strips with the greatest signal. The new one is
1345  *     fbad_cog(ncog,ic) )  c$$$*     fbad_cog(ncog,ic) )
1346  *      c$$$*    
1347  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1348    c$$$
1349        include 'commontracker.f'  c$$$      include 'commontracker.f'
1350        include 'level1.f'  c$$$      include 'level1.f'
1351        include 'calib.f'  c$$$      include 'calib.f'
1352    c$$$
1353  *     --> signal of the central strip  c$$$*     --> signal of the central strip
1354        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
1355    c$$$
1356  *     signal of adjacent strips  c$$$*     signal of adjacent strips
1357  *     --> left  c$$$*     --> left
1358        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
1359        if(  c$$$      if(
1360       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
1361       $     )  c$$$     $     )
1362       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
1363    c$$$
1364        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
1365        if(  c$$$      if(
1366       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
1367       $     )  c$$$     $     )
1368       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
1369    c$$$
1370  *     --> right  c$$$*     --> right
1371        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
1372        if(  c$$$      if(
1373       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1374       $     .or.  c$$$     $     .or.
1375       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1376       $     )  c$$$     $     )
1377       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
1378    c$$$
1379        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
1380        if(  c$$$      if(
1381       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1382       $     .or.  c$$$     $     .or.
1383       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1384       $     )  c$$$     $     )
1385       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
1386    c$$$
1387    c$$$
1388        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1389           f  = 4.  c$$$         f  = 4.
1390           si = 8.4  c$$$         si = 8.4
1391        else                              !X-view  c$$$      else                              !X-view
1392           f  = 6.  c$$$         f  = 6.
1393           si = 3.9  c$$$         si = 3.9
1394        endif  c$$$      endif
1395    c$$$
1396        fbad_cog = 1.  c$$$      fbad_cog = 1.
1397        f0 = 1  c$$$      f0 = 1
1398        f1 = 1  c$$$      f1 = 1
1399        f2 = 1  c$$$      f2 = 1
1400        f3 = 1    c$$$      f3 = 1  
1401        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
1402            c$$$        
1403           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
1404           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
1405  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
1406    c$$$
1407           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
1408              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.)
1409           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
1410              fbad_cog = 1.  c$$$            fbad_cog = 1.
1411           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
1412              fbad_cog = 1.  c$$$            fbad_cog = 1.
1413           else  c$$$         else
1414              fbad_cog = 1.  c$$$            fbad_cog = 1.
1415           endif  c$$$         endif
1416            c$$$        
1417        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
1418    c$$$
1419    c$$$
1420           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
1421           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
1422  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
1423    c$$$
1424           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
1425              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.)
1426           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
1427              fbad_cog = 1.  c$$$            fbad_cog = 1.
1428           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
1429              fbad_cog = 1.  c$$$            fbad_cog = 1.
1430           else  c$$$         else
1431              fbad_cog = 1.  c$$$            fbad_cog = 1.
1432           endif  c$$$         endif
1433    c$$$
1434        endif  c$$$      endif
1435    c$$$
1436        fbad_cog0 = sqrt(fbad_cog)  c$$$      fbad_cog0 = sqrt(fbad_cog)
1437    c$$$
1438        return  c$$$      return
1439        end  c$$$      end
1440    c$$$
1441    c$$$
1442    c$$$
1443    
1444  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1445    
1446        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1447    
1448        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1449        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1280  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1530  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1530     20 CONTINUE     20 CONTINUE
1531        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1532                
1533        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1534    
1535        END        END
1536    
1537  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1538        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1539        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1540        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1541        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1371  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1621  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1621     20 CONTINUE     20 CONTINUE
1622        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1623                
1624        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1625    
1626        END        END
1627  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1628        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1629        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1630        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1631        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1461  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1711  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1711     20 CONTINUE     20 CONTINUE
1712        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1713                
1714        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1715    
1716        END        END
1717  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1718        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1719        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1720        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1721        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1533  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1783  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1783     20 CONTINUE     20 CONTINUE
1784        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1785    
1786        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1787    
1788        END        END
1789  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1684  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1934  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1934        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1935    
1936        END        END
1937    
1938    
1939    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1940          real function pfacorr(ic,angle)
1941    *--------------------------------------------------------------
1942    *     this function returns the landi correction for this cluster
1943    *--------------------------------------------------------------
1944          include 'commontracker.f'
1945          include 'calib.f'
1946          include 'level1.f'
1947    
1948          real angle
1949          integer iview,lad
1950    
1951          iview = VIEW(ic)            
1952          lad = nld(MAXS(ic),VIEW(ic))
1953    
1954    *     find angular bin
1955    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1956          do iang=1,nangbin
1957             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1958                iangle=iang
1959                goto 98
1960             endif
1961          enddo
1962          if(DEBUG.eq.1)
1963         $     print*,'pfacorr *** warning *** angle out of range: ',angle
1964          if(angle.le.angL(1))iang=1
1965          if(angle.ge.angR(nangbin))iang=nangbin
1966     98   continue                  !jump here if ok
1967    
1968          pfacorr = fcorr(iview,lad,iang)
1969    
1970          if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr
1971    
1972    
1973     100  return
1974          end
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.23