/[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.1 by mocchiut, Fri May 19 13:15:55 2006 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  *     this file contains all subroutines and functions
3  *     that are needed for position finding algorithms  *     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)
40          
41          integer ipfa
42          character*10 cpfa
43    
44          CPFA='COG4'
45          if(ipfa.eq.0)CPFA='ETA'
46          if(ipfa.eq.2)CPFA='ETA2'
47          if(ipfa.eq.3)CPFA='ETA3'
48          if(ipfa.eq.4)CPFA='ETA4'
49          if(ipfa.eq.5)CPFA='ETAL'
50          if(ipfa.eq.10)CPFA='COG'
51          if(ipfa.eq.11)CPFA='COG1'
52          if(ipfa.eq.12)CPFA='COG2'
53          if(ipfa.eq.13)CPFA='COG3'
54          if(ipfa.eq.14)CPFA='COG4'
55          
56          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          real function fieldcorr(iview,bbb)
91    
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        real function pfa_eta(ic,angle)  
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                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
311    *     to evaluate the position of a cluster, according to the p.f.a.
312    *--------------------------------------------------------------
313          include 'commontracker.f'
314          include 'level1.f'
315          include 'calib.f'
316    
317          character*4 usedPFA
318          
319    
320    
321          call idtoc(pfaid,usedPFA)
322    
323          npfastrips=-1
324    
325          if(usedPFA.eq.'COG1')npfastrips=1
326          if(usedPFA.eq.'COG2')npfastrips=2
327          if(usedPFA.eq.'COG3')npfastrips=3
328          if(usedPFA.eq.'COG4')npfastrips=4
329          if(usedPFA.eq.'ETA2')npfastrips=2
330          if(usedPFA.eq.'ETA3')npfastrips=3
331          if(usedPFA.eq.'ETA4')npfastrips=4
332    *     ----------------------------------------------------------------
333          if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
334    c         print*,VIEW(ic),angle
335             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
336                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
337                   npfastrips=2
338                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
339                   npfastrips=3
340                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
341                   npfastrips=4
342                else
343                   npfastrips=4     !COG4
344                endif                        
345             else                   !X-view
346                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
347                   npfastrips=2
348                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
349                   npfastrips=3
350                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
351                   npfastrips=4
352                else
353                   npfastrips=4     !COG4
354                endif                        
355             endif
356          endif
357    *     ----------------------------------------------------------------
358          if(usedPFA.eq.'COG')then
359    
360             npfastrips=0
361    
362    c$$$         iv=VIEW(ic)
363    c$$$         if(mod(iv,2).eq.1)incut=incuty
364    c$$$         if(mod(iv,2).eq.0)incut=incutx
365    c$$$         istart = INDSTART(IC)
366    c$$$         istop  = TOTCLLENGTH
367    c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
368    c$$$         mu  = 0
369    c$$$         do i = INDMAX(IC),istart,-1
370    c$$$            ipos = i-INDMAX(ic)
371    c$$$            cut  = incut*CLSIGMA(i)
372    c$$$            if(CLSIGNAL(i).ge.cut)then
373    c$$$               mu = mu + 1
374    c$$$               print*,i,mu
375    c$$$            else
376    c$$$               goto 10
377    c$$$            endif
378    c$$$         enddo
379    c$$$ 10      continue
380    c$$$         do i = INDMAX(IC)+1,istop
381    c$$$            ipos = i-INDMAX(ic)
382    c$$$            cut  = incut*CLSIGMA(i)
383    c$$$            if(CLSIGNAL(i).ge.cut)then
384    c$$$               mu = mu + 1
385    c$$$               print*,i,mu
386    c$$$            else
387    c$$$               goto 20
388    c$$$            endif
389    c$$$         enddo
390    c$$$ 20      continue
391    c$$$         npfastrips=mu
392    
393          endif
394    *     ----------------------------------------------------------------
395    
396    c      print*,pfaid,usedPFA,angle,npfastrips
397    
398          return
399          end
400    
401    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
402          real function pfaeta(ic,angle)
403  *--------------------------------------------------------------  *--------------------------------------------------------------
404  *     this function returns the position (in strip units)  *     this function returns the position (in strip units)
405  *     it calls:  *     it calls:
406  *     - pfa_eta2(ic,angle)  *     - pfaeta2(ic,angle)
407  *     - pfa_eta3(ic,angle)  *     - pfaeta3(ic,angle)
408  *     - pfa_eta4(ic,angle)  *     - pfaeta4(ic,angle)
409  *     according to the angle  *     according to the angle
410  *--------------------------------------------------------------  *--------------------------------------------------------------
411        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
412        include 'level1.f'        include 'level1.f'
413          include 'calib.f'
414                
415        pfa_eta = 0        pfaeta = 0
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           pfa_eta = pfa_eta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
420                  pfaeta = pfaeta2(ic,angle)
421    cc            print*,pfaeta2(ic,angle)
422             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
423                pfaeta = pfaeta3(ic,angle)
424             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
425                pfaeta = pfaeta4(ic,angle)
426             else
427                pfaeta = cog(4,ic)
428             endif            
429    
430        else                      !X-view        else                      !X-view
431    
432           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
433              pfa_eta = pfa_eta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
434           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
435              pfa_eta = pfa_eta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
436           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
437              pfa_eta = pfa_eta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
438           endif           else
439                pfaeta = cog(4,ic)
440             endif            
441                            
442        endif        endif
443          
   
444   100  return   100  return
445        end        end
446    
447  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
448        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
449  *--------------------------------------------------------------  *--------------------------------------------------------------
450  *     this function returns the average spatial resolution  *     this function returns the position (in strip units)
 *     (in cm) for the ETA algorithm (function pfa_eta(ic,angle))  
451  *     it calls:  *     it calls:
452  *     - risx_eta2(angle)  *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
453  *     - risy_eta2(angle)  *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
454  *     - risx_eta3(angle)  *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
 *     - risx_eta4(angle)  
455  *     according to the angle  *     according to the angle
456  *--------------------------------------------------------------  *--------------------------------------------------------------
457        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
458        include 'level1.f'        include 'level1.f'
459          include 'calib.f'
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
460                
461        ris_eta = 0        pfaetal = 0
462    
463        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
464                
465           ris_eta = risy_eta2(angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
466           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)              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        else                      !X-view
477    
478           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
479              ris_eta = risx_eta2(angle)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
480           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then  cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
481              ris_eta = risx_eta3(angle)           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
482           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then              pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
483              ris_eta = risx_eta4(angle)           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
484           elseif(abs(angle).gt.21.)then              pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
485              ris_eta = risx_eta4(21.)           else
486           endif              pfaetal = cog(4,ic)
487             endif            
488                            
489        endif        endif
490          
491     100  return
492          end
493    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
494    c      real function riseta(ic,angle)
495          real function riseta(iview,angle)
496    *--------------------------------------------------------------
497    *     this function returns the average spatial resolution
498    *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
499    *     it calls:
500    *     - risxeta2(angle)
501    *     - risyeta2(angle)
502    *     - risxeta3(angle)
503    *     - risxeta4(angle)
504    *     according to the angle
505    *--------------------------------------------------------------
506          include 'commontracker.f'
507          include 'level1.f'
508          include 'calib.f'
509    
510  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'        riseta = 0
511  c$$$     $     ,' -->',ris_eta  
512    c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
513          if(mod(iview,2).eq.1)then !Y-view
514          
515    
516             if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
517                riseta = risyeta2(angle)
518             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
519                riseta = risy_cog(angle) !ATTENZIONE!!
520             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
521                riseta = risy_cog(angle) !ATTENZIONE!!
522             else
523                riseta = risy_cog(angle)
524             endif            
525    
526          else                      !X-view
527    
528             if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
529                riseta = risxeta2(angle)
530             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
531                riseta = risxeta3(angle)
532             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
533                riseta = risxeta4(angle)
534             else
535                riseta = risx_cog(angle)
536             endif            
537                
538          endif
539    
540    
541   100  return   100  return
# Line 103  c$$$     $     ,' -->',ris_eta Line 555  c$$$     $     ,' -->',ris_eta
555    
556        include 'commontracker.f'        include 'commontracker.f'
557        include 'level1.f'        include 'level1.f'
558  *      include 'calib.f'        include 'calib.f'
559        fbad_eta = 0        fbad_eta = 0
560    
561        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
562                
563           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
564                fbad_eta = fbad_cog(2,ic)
565             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
566                fbad_eta = fbad_cog(3,ic)
567             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
568                fbad_eta = fbad_cog(4,ic)
569             else
570                fbad_eta = fbad_cog(4,ic)
571             endif            
572    
573        else                      !X-view        else                      !X-view
574    
575           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
576              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
577           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
578              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
579           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
580              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
581           endif           else
582                fbad_eta = fbad_cog(4,ic)
583             endif            
584                            
585        endif        endif
586    
# Line 126  c$$$     $     ,' -->',ris_eta Line 588  c$$$     $     ,' -->',ris_eta
588        end        end
589    
590  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
591  c*****************************************************        real function pfaeta2(ic,angle) !(1)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta2(cog2,view,lad,angle)  
       real function pfa_eta2(ic,angle) !(1)  
592  *--------------------------------------------------------------  *--------------------------------------------------------------
593  *     this function returns  *     this function returns
594  *  *
# Line 149  c      real function pfa_eta2(cog2,view, Line 607  c      real function pfa_eta2(cog2,view,
607        real cog2,angle        real cog2,angle
608        integer iview,lad        integer iview,lad
609    
610        logical DEBUG        iview = VIEW(ic)            
611        common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
612          cog2 = cog(2,ic)          
613        iview = VIEW(ic)              !(1)        pfaeta2=cog2
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
       pfa_eta2=cog2  
614    
615  *     find angular bin  *     find angular bin
616  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
617        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
618           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
619              iangle=iang              iangle=iang
620              goto 98              goto 98
621           endif           endif
622        enddo        enddo
623        if(DEBUG)        if(DEBUG.EQ.1)
624       $     print*,'pfa_eta2 *** 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 177  c$$$*     find extremes of interpolation Line 631  c$$$*     find extremes of interpolation
631  c$$$      iflag=0  c$$$      iflag=0
632  c$$$*     --------------------------------  c$$$*     --------------------------------
633  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then
634  c$$$c         print*,'pfa_eta2 *** warning *** argument out of range: ',cog2  c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2
635  c$$$*         goto 100  c$$$*         goto 100
636  c$$$*     ----------------------------------------------  c$$$*     ----------------------------------------------
637  c$$$*     non salto piu`, ma scalo di 1 o -1  c$$$*     non salto piu`, ma scalo di 1 o -1
# Line 230  c            print*,'-----',x1,x2,y1,y2 Line 684  c            print*,'-----',x1,x2,y1,y2
684        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
685        BB=y1-AA*x1        BB=y1-AA*x1
686    
687        pfa_eta2 = AA*cog2+BB        pfaeta2 = AA*cog2+BB
688        pfa_eta2 = pfa_eta2 - iadd        pfaeta2 = pfaeta2 - iadd
689    
690  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
691  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
692  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
693  c$$$      endif  c$$$      endif
694  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
695  c$$$         pfa_eta2=pfa_eta2+1.   !temp  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,' -->',pfa_eta2       $     ,cog2-iadd,' -->',pfaeta2
701    
702    
703   100  return   100  return
704        end        end
705    
706  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
707  c*****************************************************        real function pfaeta3(ic,angle) !(1)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta3(cog3,view,lad,angle)  
       real function pfa_eta3(ic,angle) !(1)  
708  *--------------------------------------------------------------  *--------------------------------------------------------------
709  *     this function returns  *     this function returns
710  *  *
# Line 273  c      real function pfa_eta3(cog3,view, Line 723  c      real function pfa_eta3(cog3,view,
723        real cog3,angle        real cog3,angle
724        integer iview,lad        integer iview,lad
725    
       logical DEBUG  
       common/dbg/DEBUG  
   
726    
727        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
728        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
729        cog3 = cog(3,ic)             !(1)        cog3 = cog(3,ic)            
730        pfa_eta3=cog3        pfaeta3=cog3
731    
732  *     find angular bin  *     find angular bin
733  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 291  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*,'pfa_eta3 *** 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 354  c            print*,'-----',x1,x2,y1,y2 Line 801  c            print*,'-----',x1,x2,y1,y2
801        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
802        BB=y1-AA*x1        BB=y1-AA*x1
803    
804        pfa_eta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
805        pfa_eta3 = pfa_eta3 - iadd        pfaeta3 = pfaeta3 - iadd
806    
807  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
808  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
809  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
810  c$$$      endif  c$$$      endif
811  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
812  c$$$         pfa_eta2=pfa_eta2+1.   !temp  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,' -->',pfa_eta3       $     ,cog3-iadd,' -->',pfaeta3
818    
819   100  return   100  return
820        end        end
821    
822  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
823  c*****************************************************        real function pfaeta4(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfa_eta4(cog4,view,lad,angle)  
       real function pfa_eta4(ic,angle) !(1)  
824  *--------------------------------------------------------------  *--------------------------------------------------------------
825  *     this function returns  *     this function returns
826  *  *
# Line 396  c      real function pfa_eta4(cog4,view, Line 839  c      real function pfa_eta4(cog4,view,
839        real cog4,angle        real cog4,angle
840        integer iview,lad        integer iview,lad
841    
       logical DEBUG  
       common/dbg/DEBUG  
842    
843        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
844        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
845        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
846        pfa_eta4=cog4        pfaeta4=cog4
847    
848  *     find angular bin  *     find angular bin
849  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
# Line 413  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*,'pfa_eta4 *** 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 476  c            print*,'-----',x1,x2,y1,y2 Line 917  c            print*,'-----',x1,x2,y1,y2
917        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
918        BB=y1-AA*x1        BB=y1-AA*x1
919    
920        pfa_eta4 = AA*cog4+BB        pfaeta4 = AA*cog4+BB
921        pfa_eta4 = pfa_eta4 - iadd        pfaeta4 = pfaeta4 - iadd
922    
923  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
924  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
925  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
926  c$$$      endif  c$$$      endif
927  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
928  c$$$         pfa_eta2=pfa_eta2+1.   !temp  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,' -->',pfa_eta4       $     ,cog4-iadd,' -->',pfaeta4
934    
935   100  return   100  return
936        end        end
# Line 497  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 613  c      print *,ncog,ic,cog,'//////////// Line 957  c      print *,ncog,ic,cog,'////////////
957        include 'calib.f'        include 'calib.f'
958        include 'level1.f'        include 'level1.f'
959                
       logical DEBUG  
       common/dbg/DEBUG  
960    
961    
962        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 625  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 645  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 655  c      print *,ncog,ic,cog,'//////////// Line 997  c      print *,ncog,ic,cog,'////////////
997                    
998           COG = 0.           COG = 0.
999                    
1000           if(ncog.eq.2)then  c         print*,'## ',sl2,sl1,sc,sr1,sr2
1001    
1002    c     ==============================================================
1003             if(ncog.eq.1)then
1004                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
1009                COG = 0.
1010              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1011                 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                 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              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                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sss).ne.0)
1046              elseif(sl2.le.sr2)then       $              COG = (COG-2*sl2)/(sl2+sss)
1047                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)              elseif(sl2.lt.sr2)then
1048                   if((sr2+sss).ne.0)
1049         $              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              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
1074              COG = 0.              COG = 0.
1075           endif           endif
1076    
# Line 685  ' Line 1084  '
1084           iv=VIEW(ic)           iv=VIEW(ic)
1085           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
1086           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
1087             istart = INDSTART(IC)
1088           istart=INDSTART(IC)           istop  = TOTCLLENGTH
          istop=TOTCLLENGTH  
1089           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
1090           COG=0             COG = 0  
1091           mu=0           SGN = 0.
1092           do i=istart,istop           mu  = 0
1093              ipos=i-INDMAX(ic)  c         print*,'-------'
1094              ivk=nvk(MAXS(ic)+ipos)           do i = INDMAX(IC),istart,-1
1095              is=nst(MAXS(ic)+ipos)              ipos = i-INDMAX(ic)
1096  *            print*,'******************',istart,istop,ipos              cut  = incut*CLSIGMA(i)
1097  *     $           ,MAXS(ic)+ipos,iv,ivk,is              if(CLSIGNAL(i).ge.cut)then              
1098              cut=incut*SIGMA(iv,ivk,is)                 COG = COG + ipos*CLSIGNAL(i)
1099                   SGN = SGN + CLSIGNAL(i)
1100                   mu = mu + 1
1101    c               print*,ipos,CLSIGNAL(i)
1102                else
1103                   goto 10
1104                endif
1105             enddo
1106     10      continue
1107             do i = INDMAX(IC)+1,istop
1108                ipos = i-INDMAX(ic)
1109                cut  = incut*CLSIGMA(i)
1110              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
1111                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1112                   SGN = SGN + CLSIGNAL(i)
1113                 mu = mu + 1                 mu = mu + 1
1114  c               print*,ipos,CLSIGNAL(i),incut,cut  c               print*,ipos,CLSIGNAL(i)
1115                else
1116                   goto 20
1117              endif              endif
1118           enddo           enddo
1119           COG=COG/DEDX(ic)   20      continue
1120  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'           if(SGN.le.0)then
1121  c     $        ,cog              print*,'cog(0,ic) --> ic, dedx ',ic,SGN
1122                print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1123                print*,(CLSIGNAL(i),i=istart,istop)
1124    c            print*,'cog(0,ic) --> NOT EVALUATED '
1125             else
1126                COG=COG/SGN
1127             endif
1128    c         print*,'-------'
1129                    
1130        else        else
1131                    
# Line 717  c     $        ,cog Line 1136  c     $        ,cog
1136    
1137        endif        endif
1138    
1139  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1140    
1141        return        return
1142        end        end
1143    
1144  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1145    
1146        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
1147  *-------------------------------------------------------  *-------------------------------------------------------
1148  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 738  c      print *,ncog,ic,cog,'//////////// Line 1158  c      print *,ncog,ic,cog,'////////////
1158        include 'calib.f'        include 'calib.f'
1159    
1160        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1161           f  = 4.           si = 8.4  !average good-strip noise
1162           si = 8.4           f  = 4.   !average bad-strip noise: f*si
1163           incut=incuty           incut=incuty
1164        else                      !X-view        else                      !X-view
1165           f  = 6.           si = 3.9  !average good-strip noise
1166           si = 3.9           f  = 6.   !average bad-strip noise: f*si
1167           incut=incutx           incut=incutx
1168        endif        endif
1169                
# Line 754  c      print *,ncog,ic,cog,'//////////// Line 1174  c      print *,ncog,ic,cog,'////////////
1174  *     --> signal of the central strip  *     --> signal of the central strip
1175           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1176           fsc = 1           fsc = 1
1177           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
1178             fsc = clsigma(INDMAX(ic))/si
1179  *     --> signal of adjacent strips  *     --> signal of adjacent strips
1180           sl1 = 0                !left 1           sl1 = 0                !left 1
1181           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 762  c      print *,ncog,ic,cog,'//////////// Line 1183  c      print *,ncog,ic,cog,'////////////
1183       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1184       $        )then       $        )then
1185              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
1186              if(BAD(VIEW(ic),nvk(MAXS(ic)-1),nst(MAXS(ic)-1)).eq.0)fsl1=f  c            if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f
1187  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
1188           endif           endif
1189    
1190           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 773  c            fsl1 = 0 Line 1193  c            fsl1 = 0
1193       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1194       $        )then       $        )then
1195              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
1196              if(BAD(VIEW(ic),nvk(MAXS(ic)-2),nst(MAXS(ic)-2)).eq.0)fsl2=f  c            if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f
1197  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
1198           endif           endif
1199           sr1 = 0                !right 1           sr1 = 0                !right 1
1200           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 785  c            fsl2 = 0 Line 1204  c            fsl2 = 0
1204       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1205       $        )then       $        )then
1206              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
1207              if(BAD(VIEW(ic),nvk(MAXS(ic)+1),nst(MAXS(ic)+1)).eq.0)fsr1=f  c            if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f
1208  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
1209           endif               endif    
1210           sr2 = 0                !right 2           sr2 = 0                !right 2
1211           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 797  c            fsr1 = 0 Line 1215  c            fsr1 = 0
1215       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1216       $        )then       $        )then
1217              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
1218              if(BAD(VIEW(ic),nvk(MAXS(ic)+2),nst(MAXS(ic)+2)).eq.0)fsr2=f  c            if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f
1219  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
1220           endif           endif
1221    
1222    
1223    
1224  ************************************************************  ************************************************************
1225  *     COG computation  *     COG2-3-4 computation
1226  ************************************************************  ************************************************************
1227    
1228  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
1229                    
1230           COG = 0.           vCOG = cog(ncog,ic)!0.
1231                    
1232           if(ncog.eq.2)then           if(ncog.eq.2)then
1233              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1234                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
1235                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1236                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1237              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
1238                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
1239                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1240                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1241              endif              endif
1242           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1243              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
1244              fbad_cog =              fbad_cog =
1245       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1246              fbad_cog =              fbad_cog =
1247       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1248           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1249              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1250                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1251                 fbad_cog =                 fbad_cog =
1252       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1253       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1254                 fbad_cog =                 fbad_cog =
1255       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1256       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
1257              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
1258                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1259                 fbad_cog =                 fbad_cog =
1260       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1261       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1262                 fbad_cog =                 fbad_cog =
1263       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1264       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1265              endif              endif
1266           else           else
1267              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1268              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1269              COG = 0.  c            COG = 0.
1270           endif           endif
1271                    
1272        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1273    *     =========================
1274    *     COG computation
1275    *     =========================
1276                    
1277           iv=VIEW(ic)           vCOG = cog(0,ic)
1278           istart=INDSTART(IC)  
1279           istop=TOTCLLENGTH           iv     = VIEW(ic)
1280           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1281           COG=0.           istop  = TOTCLLENGTH
1282           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1283           SDE=0.           SGN = 0.
1284           do i=istart,istop           SNU = 0.
1285              ipos=i-INDMAX(ic)           SDE = 0.
1286              il=nvk(MAXS(ic)+ipos)  
1287              is=nst(MAXS(ic)+ipos)           do i=INDMAX(IC),istart,-1
1288              cut=incut*SIGMA(iv,il,is)              ipos = i-INDMAX(ic)
1289                cut  = incut*CLSIGMA(i)
1290              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1291                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1292              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1293                   SDE = SDE + (ipos-vCOG)**2
1294                else
1295                   goto 10
1296                endif            
1297           enddo           enddo
1298           COG=COG/DEDX(ic)   10      continue
1299           do i=istart,istop           do i=INDMAX(IC)+1,istop
1300              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1301              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1302              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1303                 fs=1                 fs = clsigma(i)/si
1304                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1305                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1306                 SDE = SDE + (ipos-COG)**2              else
1307                   goto 20
1308              endif                          endif            
1309           enddo           enddo
1310           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1311             if(SDE.ne.0)then
1312                FBAD_COG=SNU/SDE
1313             else
1314                
1315             endif
1316    
1317        else        else
1318                                        
# Line 902  c      print*,sl2,sl1,sc,sr1,sr2 Line 1330  c      print*,sl2,sl1,sc,sr1,sr2
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 1087  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1515  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1515       +/       +/
1516    
1517        V(1)= abs(x)        V(1)= abs(x)
1518          if(V(1).gt.20.)V(1)=20.
1519    
1520        HQUADF = 0.        HQUADF = 0.
1521        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1101  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 1177  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1606  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1606       +/       +/
1607    
1608        V(1) =  abs(x)        V(1) =  abs(x)
1609          if(V(1).gt.20.)V(1)=20.
1610    
1611        HQUADF = 0.        HQUADF = 0.
1612        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1613           HQDJ = 0.           HQDJ = 0.
# Line 1190  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 1265  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1696  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1696       +/       +/
1697    
1698        V(1)=abs(x)        V(1)=abs(x)
1699          if(V(1).gt.20.)V(1)=20.
1700    
1701        HQUADF = 0.        HQUADF = 0.
1702        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1279  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 1336  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1768  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1768       +/       +/
1769    
1770        v(1)= abs(x)        v(1)= abs(x)
1771          if(V(1).gt.20.)V(1)=20.
1772                
1773        HQUADF = 0.        HQUADF = 0.
1774        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1350  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 1402  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1835  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1835       +/       +/
1836                
1837        V(1)=abs(x)        V(1)=abs(x)
1838          if(V(1).gt.20.)V(1)=20.
1839    
1840        HQUADF = 0.        HQUADF = 0.
1841        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1482  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1916  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1916       +/       +/
1917    
1918        V(1)=abs(x)        V(1)=abs(x)
1919          if(V(1).gt.20.)V(1)=20.
1920    
1921        HQUADF = 0.        HQUADF = 0.
1922        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1499  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.1  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.23