/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/functionspfa.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/functionspfa.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.7 by pam-fi, Thu Jan 11 10:20:58 2007 UTC revision 1.23 by pam-fi, Wed Mar 5 17:00:20 2008 UTC
# Line 1  Line 1 
1  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2  *     this file contains all subroutines and functions  *     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    cc         tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!!
74             tgtemp  = tan(angx*acos(-1.)/180.) + pmuH_h*by*0.00001
75    
76          elseif(mod(iview,2).eq.1)then
77    c     =================================================
78    c     Y view
79    c     =================================================        
80    c     here bbb is the x component of the m.filed
81             angy = ang
82             bx   = bbb
83             tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001        
84    
85          endif      
86          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
87    
88          return
89          end
90  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
91          real function fieldcorr(iview,bbb)
92    
93          include 'commontracker.f'
94    
95          fieldcorr = 0.
96    
97          if(mod(iview,2).eq.0)then
98    
99    c     =================================================
100    c     X view
101    c     =================================================
102    c     here bbb is the y component of the m.field
103             by   = bbb
104             if(iview.eq.12) by = -1. * bbb
105             fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
106    
107          elseif(mod(iview,2).eq.1)then
108    c     =================================================
109    c     Y view
110    c     =================================================        
111    c     here bbb is the x component of the m.filed
112             bx   = bbb
113             fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
114    
115          endif      
116          
117          return
118          end
119    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
120    
121          subroutine applypfa(PFAtt,ic,ang,corr,res)
122    *---------------------------------------------------------------
123    *     this subroutine calculate the coordinate of cluster ic (in
124    *     strip units), relative to the strip with the maximum signal,
125    *     and its spatial resolution (in cm), applying PFAtt.
126    *     ang is the effective angle, relative to the sensor
127    *---------------------------------------------------------------
128    
129          character*4 PFAtt
130          include 'commontracker.f'
131          include 'level1.f'
132    
133          corr = 0
134          res  = 0
135    
136          if(ic.le.0)return
137    
138          iview   = VIEW(ic)
139    
140          if(mod(iview,2).eq.0)then
141    c     =================================================
142    c     X view
143    c     =================================================
144    
145             res = RESXAV
146    
147             if(PFAtt.eq.'COG1')then
148    
149                corr   = 0
150                res = 1e-4*pitchX/sqrt(12.)!!res
151    
152             elseif(PFAtt.eq.'COG2')then
153    
154                corr    = cog(2,ic)            
155                res = risx_cog(abs(ang))!TEMPORANEO              
156                res = res*fbad_cog(2,ic)
157    
158             elseif(PFAtt.eq.'COG3')then
159    
160                corr    = cog(3,ic)            
161                res = risx_cog(abs(ang))!TEMPORANEO                      
162                res = res*fbad_cog(3,ic)
163    
164             elseif(PFAtt.eq.'COG4')then
165    
166                corr    = cog(4,ic)            
167                res = risx_cog(abs(ang))!TEMPORANEO                      
168                res = res*fbad_cog(4,ic)
169    
170             elseif(PFAtt.eq.'ETA2')then
171    
172                corr  = pfaeta2(ic,ang)          
173                res = risxeta2(abs(ang))
174                res = res*fbad_cog(2,ic)
175    
176             elseif(PFAtt.eq.'ETA3')then                        
177    
178                corr  = pfaeta3(ic,ang)          
179                res = risxeta3(abs(ang))                      
180                res = res*fbad_cog(3,ic)              
181    
182             elseif(PFAtt.eq.'ETA4')then                        
183    
184                corr  = pfaeta4(ic,ang)            
185                res = risxeta4(abs(ang))                      
186                res = res*fbad_cog(4,ic)              
187    
188             elseif(PFAtt.eq.'ETA')then  
189    
190                corr  = pfaeta(ic,ang)            
191    c            res = riseta(ic,ang)                    
192                res = riseta(iview,ang)                    
193                res = res*fbad_eta(ic,ang)            
194    
195             elseif(PFAtt.eq.'ETAL')then  
196    
197                corr  = pfaetal(ic,ang)            
198                res = riseta(iview,ang)                    
199                res = res*fbad_eta(ic,ang)            
200    
201             elseif(PFAtt.eq.'COG')then          
202    
203                corr  = cog(0,ic)            
204                res = risx_cog(abs(ang))                    
205                res = res*fbad_cog(0,ic)
206    
207             else
208                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
209             endif
210    
211    
212    *     ======================================
213    *     temporary patch for saturated clusters
214    *     ======================================
215             if( nsatstrips(ic).gt.0 )then
216                corr  = cog(4,ic)            
217                res = pitchX*1e-4/sqrt(12.)
218    cc            cc=cog(4,ic)
219    c$$$            print*,ic,' *** ',cc
220    c$$$            print*,ic,' *** ',res
221             endif
222    
223    
224          elseif(mod(iview,2).eq.1)then
225    c     =================================================
226    c     Y view
227    c     =================================================
228    
229             res = RESYAV
230    
231             if(PFAtt.eq.'COG1')then
232    
233                corr  = 0  
234                res = 1e-4*pitchY/sqrt(12.)!res  
235    
236             elseif(PFAtt.eq.'COG2')then
237    
238                corr    = cog(2,ic)
239                res = risy_cog(abs(ang))!TEMPORANEO
240                res = res*fbad_cog(2,ic)
241    
242             elseif(PFAtt.eq.'COG3')then
243    
244                corr    = cog(3,ic)
245                res = risy_cog(abs(ang))!TEMPORANEO
246                res = res*fbad_cog(3,ic)
247    
248             elseif(PFAtt.eq.'COG4')then
249    
250                corr    = cog(4,ic)
251                res = risy_cog(abs(ang))!TEMPORANEO
252                res = res*fbad_cog(4,ic)
253    
254             elseif(PFAtt.eq.'ETA2')then
255    
256                corr    = pfaeta2(ic,ang)          
257                res = risyeta2(abs(ang))              
258                res = res*fbad_cog(2,ic)
259    
260             elseif(PFAtt.eq.'ETA3')then                      
261    
262                corr    = pfaeta3(ic,ang)
263                res = res*fbad_cog(3,ic)  
264    
265             elseif(PFAtt.eq.'ETA4')then  
266    
267                corr    = pfaeta4(ic,ang)
268                res = res*fbad_cog(4,ic)
269    
270             elseif(PFAtt.eq.'ETA')then
271    
272                corr    = pfaeta(ic,ang)
273    c            res = riseta(ic,ang)  
274                res = riseta(iview,ang)  
275                res = res*fbad_eta(ic,ang)
276    
277             elseif(PFAtt.eq.'ETAL')then
278    
279                corr    = pfaetal(ic,ang)
280                res = riseta(iview,ang)  
281                res = res*fbad_eta(ic,ang)
282    
283             elseif(PFAtt.eq.'COG')then
284    
285                corr    = cog(0,ic)            
286                res = risy_cog(abs(ang))
287                res = res*fbad_cog(0,ic)
288    
289             else
290                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
291             endif
292    
293    
294    *     ======================================
295    *     temporary patch for saturated clusters
296    *     ======================================
297             if( nsatstrips(ic).gt.0 )then
298                corr    = cog(4,ic)            
299                res = pitchY*1e-4/sqrt(12.)
300    cc            cc=cog(4,ic)
301    c$$$            print*,ic,' *** ',cc
302    c$$$            print*,ic,' *** ',res
303             endif
304            
305          endif
306          end
307    
308    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
309          integer function npfastrips(ic,angle)
310    *--------------------------------------------------------------
311    *     thid function returns the number of strips used
312    *     to evaluate the position of a cluster, according to the p.f.a.
313    *--------------------------------------------------------------
314          include 'commontracker.f'
315          include 'level1.f'
316          include 'calib.f'
317    
318          character*4 usedPFA
319          
320    
321    
322          call idtoc(pfaid,usedPFA)
323    
324          npfastrips=-1
325    
326          if(usedPFA.eq.'COG1')npfastrips=1
327          if(usedPFA.eq.'COG2')npfastrips=2
328          if(usedPFA.eq.'COG3')npfastrips=3
329          if(usedPFA.eq.'COG4')npfastrips=4
330          if(usedPFA.eq.'ETA2')npfastrips=2
331          if(usedPFA.eq.'ETA3')npfastrips=3
332          if(usedPFA.eq.'ETA4')npfastrips=4
333    *     ----------------------------------------------------------------
334          if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
335    c         print*,VIEW(ic),angle
336             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
337                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
338                   npfastrips=2
339                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
340                   npfastrips=3
341                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
342                   npfastrips=4
343                else
344                   npfastrips=4     !COG4
345                endif                        
346             else                   !X-view
347                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
348                   npfastrips=2
349                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
350                   npfastrips=3
351                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
352                   npfastrips=4
353                else
354                   npfastrips=4     !COG4
355                endif                        
356             endif
357          endif
358    *     ----------------------------------------------------------------
359          if(usedPFA.eq.'COG')then
360    
361             npfastrips=0
362    
363    c$$$         iv=VIEW(ic)
364    c$$$         if(mod(iv,2).eq.1)incut=incuty
365    c$$$         if(mod(iv,2).eq.0)incut=incutx
366    c$$$         istart = INDSTART(IC)
367    c$$$         istop  = TOTCLLENGTH
368    c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
369    c$$$         mu  = 0
370    c$$$         do i = INDMAX(IC),istart,-1
371    c$$$            ipos = i-INDMAX(ic)
372    c$$$            cut  = incut*CLSIGMA(i)
373    c$$$            if(CLSIGNAL(i).ge.cut)then
374    c$$$               mu = mu + 1
375    c$$$               print*,i,mu
376    c$$$            else
377    c$$$               goto 10
378    c$$$            endif
379    c$$$         enddo
380    c$$$ 10      continue
381    c$$$         do i = INDMAX(IC)+1,istop
382    c$$$            ipos = i-INDMAX(ic)
383    c$$$            cut  = incut*CLSIGMA(i)
384    c$$$            if(CLSIGNAL(i).ge.cut)then
385    c$$$               mu = mu + 1
386    c$$$               print*,i,mu
387    c$$$            else
388    c$$$               goto 20
389    c$$$            endif
390    c$$$         enddo
391    c$$$ 20      continue
392    c$$$         npfastrips=mu
393    
394          endif
395    *     ----------------------------------------------------------------
396    
397    c      print*,pfaid,usedPFA,angle,npfastrips
398    
399          return
400          end
401    
402  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
403        real function pfaeta(ic,angle)        real function pfaeta(ic,angle)
# Line 17  Line 410 
410  *     according to the angle  *     according to the angle
411  *--------------------------------------------------------------  *--------------------------------------------------------------
412        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
413        include 'level1.f'        include 'level1.f'
414          include 'calib.f'
415                
416        pfaeta = 0        pfaeta = 0
417    
418        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
419                
420           pfaeta = pfaeta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
421                  pfaeta = pfaeta2(ic,angle)
422    cc            print*,pfaeta2(ic,angle)
423             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
424                pfaeta = pfaeta3(ic,angle)
425             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
426                pfaeta = pfaeta4(ic,angle)
427             else
428                pfaeta = cog(4,ic)
429             endif            
430    
431        else                      !X-view        else                      !X-view
432    
433           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
434              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
435           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
436              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
437           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
438              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
439           endif           else
440                pfaeta = cog(4,ic)
441             endif            
442                            
443        endif        endif
444                
 c      print*,'pfaeta ',pfaeta, angle  
   
445   100  return   100  return
446        end        end
447    
448  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
449        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
450  *--------------------------------------------------------------  *--------------------------------------------------------------
451  *     this function returns the average spatial resolution  *     this function returns the position (in strip units)
 *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  
452  *     it calls:  *     it calls:
453  *     - risx_eta2(angle)  *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
454  *     - risy_eta2(angle)  *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
455  *     - risx_eta3(angle)  *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
 *     - risx_eta4(angle)  
456  *     according to the angle  *     according to the angle
457  *--------------------------------------------------------------  *--------------------------------------------------------------
458        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
459        include 'level1.f'        include 'level1.f'
460          include 'calib.f'
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
461                
462        ris_eta = 0        pfaetal = 0
463    
464        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
465                
466           ris_eta = risy_eta2(angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
467           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
468    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
469             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
470                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
471             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
472                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
473             else
474                pfaetal = cog(4,ic)
475             endif            
476    
477        else                      !X-view        else                      !X-view
478    
479           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
480              ris_eta = risx_eta2(angle)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
481           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then  cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
482              ris_eta = risx_eta3(angle)           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
483           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then              pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
484              ris_eta = risx_eta4(angle)           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
485           elseif(abs(angle).gt.21.)then              pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
486              ris_eta = risx_eta4(21.)           else
487           endif              pfaetal = cog(4,ic)
488             endif            
489                            
490        endif        endif
491          
492     100  return
493          end
494    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
495    c      real function riseta(ic,angle)
496          real function riseta(iview,angle)
497    *--------------------------------------------------------------
498    *     this function returns the average spatial resolution
499    *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
500    *     it calls:
501    *     - risxeta2(angle)
502    *     - risyeta2(angle)
503    *     - risxeta3(angle)
504    *     - risxeta4(angle)
505    *     according to the angle
506    *--------------------------------------------------------------
507          include 'commontracker.f'
508          include 'level1.f'
509          include 'calib.f'
510    
511          riseta = 0
512    
513    c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
514          if(mod(iview,2).eq.1)then !Y-view
515          
516    
517             if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
518                riseta = risyeta2(angle)
519             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
520                riseta = risy_cog(angle) !ATTENZIONE!!
521             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
522                riseta = risy_cog(angle) !ATTENZIONE!!
523             else
524                riseta = risy_cog(angle)
525             endif            
526    
527  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'        else                      !X-view
528  c$$$     $     ,' -->',ris_eta  
529             if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
530                riseta = risxeta2(angle)
531             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
532                riseta = risxeta3(angle)
533             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
534                riseta = risxeta4(angle)
535             else
536                riseta = risx_cog(angle)
537             endif            
538                
539          endif
540    
541    
542   100  return   100  return
# Line 104  c$$$     $     ,' -->',ris_eta Line 556  c$$$     $     ,' -->',ris_eta
556    
557        include 'commontracker.f'        include 'commontracker.f'
558        include 'level1.f'        include 'level1.f'
559  *      include 'calib.f'        include 'calib.f'
560        fbad_eta = 0        fbad_eta = 0
561    
562        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
563                
564           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
565                fbad_eta = fbad_cog(2,ic)
566             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
567                fbad_eta = fbad_cog(3,ic)
568             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
569                fbad_eta = fbad_cog(4,ic)
570             else
571                fbad_eta = fbad_cog(4,ic)
572             endif            
573    
574        else                      !X-view        else                      !X-view
575    
576           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
577              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
578           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
579              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
580           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
581              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
582           endif           else
583                fbad_eta = fbad_cog(4,ic)
584             endif            
585                            
586        endif        endif
587    
# Line 127  c$$$     $     ,' -->',ris_eta Line 589  c$$$     $     ,' -->',ris_eta
589        end        end
590    
591  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta2(cog2,view,lad,angle)  
592        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle) !(1)
593  *--------------------------------------------------------------  *--------------------------------------------------------------
594  *     this function returns  *     this function returns
# Line 150  c      real function pfaeta2(cog2,view,l Line 608  c      real function pfaeta2(cog2,view,l
608        real cog2,angle        real cog2,angle
609        integer iview,lad        integer iview,lad
610    
611  c      logical DEBUG        iview = VIEW(ic)            
612  c      common/dbg/DEBUG        lad = nld(MAXS(ic),VIEW(ic))
613          cog2 = cog(2,ic)          
 c      print*,'## pfaeta2 ',ic,angle  
       iview = VIEW(ic)              !(1)  
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
614        pfaeta2=cog2        pfaeta2=cog2
615    
616    *     ----------------
617  *     find angular bin  *     find angular bin
618    *     ----------------
619  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
620        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
621           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
622              iangle=iang              iangle=iang
623              goto 98              goto 98
624           endif           endif
625        enddo        enddo
626        if(DEBUG)        if(DEBUG.EQ.1)
627       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
628        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
629        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
630   98   continue                  !jump here if ok   98   continue                  !jump here if ok
631    
632    *     -------------
633    *     within +/-0.5
634    *     -------------
635    
636  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  
 c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2  
 c$$$*         goto 100  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
637        iadd=0        iadd=0
638   10   continue   10   continue
639        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
640           cog2 = cog2 + 1           cog2 = cog2 + 1
641           iadd = iadd + 1           iadd = iadd + 1
642             if(iadd>iaddmax)goto 111
643           goto 10           goto 10
644        endif        endif
645   20   continue   20   continue
646        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
647           cog2 = cog2 - 1           cog2 = cog2 - 1
648           iadd = iadd - 1           iadd = iadd - 1
649             if(iadd<-1*iaddmax)goto 111
650           goto 20           goto 20
651        endif        endif
652          goto 1111
653     111  continue
654          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
655          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
656          cog2=0
657     1111 continue
658    
659  *     --------------------------------  *     --------------------------------
660  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 244  c$$$         pfaeta2=pfaeta2+1.   !temp Line 689  c$$$         pfaeta2=pfaeta2+1.   !temp
689  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
690  c$$$      endif  c$$$      endif
691    
692        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
693       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
694    
695    
# Line 252  c$$$      endif Line 697  c$$$      endif
697        end        end
698    
699  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta3(cog3,view,lad,angle)  
700        real function pfaeta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
701  *--------------------------------------------------------------  *--------------------------------------------------------------
702  *     this function returns  *     this function returns
# Line 275  c      real function pfaeta3(cog3,view,l Line 716  c      real function pfaeta3(cog3,view,l
716        real cog3,angle        real cog3,angle
717        integer iview,lad        integer iview,lad
718    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
719    
720  c      print*,'## pfaeta3 ',ic,angle        iview = VIEW(ic)            
721          lad = nld(MAXS(ic),VIEW(ic))
722        iview = VIEW(ic)              !(1)        cog3 = cog(3,ic)
723        lad = nld(MAXS(ic),VIEW(ic)) !(1)        cc = cog3
724        cog3 = cog(3,ic)             !(1)        cog3 = cc
725        pfaeta3=cog3        pfaeta3=cog3
726    
727    *     ----------------
728  *     find angular bin  *     find angular bin
729    *     ----------------
730  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
731        do iang=1,nangbin        do iang=1,nangbin
732  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 294  c         print*,'~~~~~~~~~~~~ ',iang,an Line 735  c         print*,'~~~~~~~~~~~~ ',iang,an
735              goto 98              goto 98
736           endif           endif
737        enddo        enddo
738        if(DEBUG)        if(DEBUG.EQ.1)
739       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
740        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
741        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
742   98   continue                  !jump here if ok   98   continue                  !jump here if ok
743    
744    *     -------------
745    *     within +/-0.5
746    *     -------------
747    
748  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
749        iadd=0        iadd=0
750   10   continue   10   continue
751        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
752           cog3 = cog3 + 1           cog3   = cog3 + 1.
753           iadd = iadd + 1           iadd = iadd + 1
754             if(iadd>iaddmax) goto 111
755           goto 10           goto 10
756        endif        endif
757   20   continue   20   continue
758        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
759           cog3 = cog3 - 1           cog3 = cog3 - 1.
760           iadd = iadd - 1           iadd = iadd - 1
761             if(iadd<-1*iaddmax) goto 111
762           goto 20           goto 20
763        endif        endif
764          goto 1111
765     111  continue
766          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
767          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
768          cog3=0      
769     1111 continue
770    
771  *     --------------------------------  *     --------------------------------
772  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 360  c            print*,'-----',x1,x2,y1,y2 Line 792  c            print*,'-----',x1,x2,y1,y2
792        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
793        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
794    
 c$$$      if(iflag.eq.1)then  
 c$$$         pfaeta2=pfaeta2-1.   !temp  
 c$$$         cog2=cog2-1.           !temp  
 c$$$      endif  
 c$$$      if(iflag.eq.-1)then  
 c$$$         pfaeta2=pfaeta2+1.   !temp  
 c$$$         cog2=cog2+1.           !temp  
 c$$$      endif  
795    
796        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
797       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
798    
799   100  return   100  return
800        end        end
801    
802  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
803  c*****************************************************        real function pfaeta4(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta4(cog4,view,lad,angle)  
       real function pfaeta4(ic,angle) !(1)  
804  *--------------------------------------------------------------  *--------------------------------------------------------------
805  *     this function returns  *     this function returns
806  *  *
# Line 399  c      real function pfaeta4(cog4,view,l Line 819  c      real function pfaeta4(cog4,view,l
819        real cog4,angle        real cog4,angle
820        integer iview,lad        integer iview,lad
821    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 c      print*,'## pfaeta4 ',ic,angle  
822    
823        iview = VIEW(ic)             !(1)        iview = VIEW(ic)            
824        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
825        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
826        pfaeta4=cog4        pfaeta4=cog4
827    
828    *     ----------------
829  *     find angular bin  *     find angular bin
830    *     ----------------
831  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
832        do iang=1,nangbin        do iang=1,nangbin
833  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 418  c         print*,'~~~~~~~~~~~~ ',iang,an Line 836  c         print*,'~~~~~~~~~~~~ ',iang,an
836              goto 98              goto 98
837           endif           endif
838        enddo        enddo
839        if(DEBUG)        if(DEBUG.EQ.1)
840       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
841        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
842        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
843   98   continue                  !jump here if ok   98   continue                  !jump here if ok
844    
845    *     -------------
846    *     within +/-0.5
847    *     -------------
848    
849  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
850        iadd=0        iadd=0
851   10   continue   10   continue
852        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
853           cog4 = cog4 + 1           cog4 = cog4 + 1
854           iadd = iadd + 1           iadd = iadd + 1
855             if(iadd>iaddmax)goto 111
856           goto 10           goto 10
857        endif        endif
858   20   continue   20   continue
859        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
860           cog4 = cog4 - 1           cog4 = cog4 - 1
861           iadd = iadd - 1           iadd = iadd - 1
862             if(iadd<-1*iaddmax)goto 111
863           goto 20           goto 20
864        endif        endif
865          goto 1111
866     111  continue
867          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
868          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
869          cog4=0
870     1111 continue
871    
872  *     --------------------------------  *     --------------------------------
873  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 493  c$$$         pfaeta2=pfaeta2+1.   !temp Line 902  c$$$         pfaeta2=pfaeta2+1.   !temp
902  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
903  c$$$      endif  c$$$      endif
904    
905        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
906       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
907    
908   100  return   100  return
# Line 502  c$$$      endif Line 911  c$$$      endif
911    
912    
913  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
       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  
   
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
914        real function cog(ncog,ic)        real function cog(ncog,ic)
915  *-------------------------------------------------  *-------------------------------------------------
916  *     this function returns  *     this function returns
# Line 618  c      print *,ncog,ic,cog,'//////////// Line 930  c      print *,ncog,ic,cog,'////////////
930        include 'calib.f'        include 'calib.f'
931        include 'level1.f'        include 'level1.f'
932                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
933    
934    
935        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 630  c      common/dbg/DEBUG Line 940  c      common/dbg/DEBUG
940  *     --> signal of the central strip  *     --> signal of the central strip
941           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
942  *     signal of adjacent strips  *     signal of adjacent strips
943           sl1 = 0                !left 1           sl1 = -9999.           !left 1
944           if(           if(
945       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
946       $        )       $        )
947       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
948                    
949           sl2 = 0                !left 2           sl2 = -9999.           !left 2
950           if(           if(
951       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
952       $        )       $        )
953       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
954                    
955           sr1 = 0                !right 1           sr1 = -9999.           !right 1
956           if(           if(
957       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
958       $        .or.       $        .or.
# Line 650  c      common/dbg/DEBUG Line 960  c      common/dbg/DEBUG
960       $        )       $        )
961       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
962                    
963           sr2 = 0                !right 2           sr2 = -9999.           !right 2
964           if(           if(
965       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
966       $        .or.       $        .or.
# Line 660  c      common/dbg/DEBUG Line 970  c      common/dbg/DEBUG
970                    
971           COG = 0.           COG = 0.
972                    
973  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
974    
975    c     ==============================================================
976           if(ncog.eq.1)then           if(ncog.eq.1)then
977              COG = 0.              COG = 0.
978                if(sr1.gt.sc)cog=1.
979                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
980    c     ==============================================================
981           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
982                COG = 0.
983              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
984                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
985              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
986                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
987              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
988                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
989         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
990                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
991         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
992                endif
993    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
994    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
995    c     ==============================================================
996           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
997               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
998                sss = sc
999                if( sl1.ne.-9999. )COG = COG-sl1
1000                if( sl1.ne.-9999. )sss = sss+sl1
1001                if( sr1.ne.-9999. )COG = COG+sr1
1002                if( sr1.ne.-9999. )sss = sss+sr1
1003                if(sss.ne.0)COG=COG/sss
1004    
1005    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1006    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1007    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1008    c     ==============================================================
1009           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1010    
1011                COG = 0
1012                sss = sc
1013                if( sl1.ne.-9999. )COG = COG-sl1
1014                if( sl1.ne.-9999. )sss = sss+sl1
1015                if( sr1.ne.-9999. )COG = COG+sr1
1016                if( sr1.ne.-9999. )sss = sss+sr1
1017              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1018                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1019       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1020              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
1021                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1022       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1023                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1024                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1025         $              .and.(sl2+sss).ne.0 )
1026         $              cog = (cog-2*sl2)/(sl2+sss)
1027                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1028         $              .and.(sr2+sss).ne.0 )
1029         $              cog = (2*sr2+cog)/(sr2+sss)              
1030              endif              endif
1031    c     ==============================================================
1032             elseif(ncog.eq.5)then
1033                COG = 0
1034                sss = sc
1035                if( sl1.ne.-9999. )COG = COG-sl1
1036                if( sl1.ne.-9999. )sss = sss+sl1
1037                if( sr1.ne.-9999. )COG = COG+sr1
1038                if( sr1.ne.-9999. )sss = sss+sr1
1039                if( sl2.ne.-9999. )COG = COG-2*sl2
1040                if( sl2.ne.-9999. )sss = sss+sl2
1041                if( sr2.ne.-9999. )COG = COG+2*sr2
1042                if( sr2.ne.-9999. )sss = sss+sr2
1043                if(sss.ne.0)COG=COG/sss
1044           else           else
1045              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1046       $           ,' not implemented'       $           ,' not implemented'
# Line 696  ' Line 1057  '
1057           iv=VIEW(ic)           iv=VIEW(ic)
1058           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
1059           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
1060           istart = INDSTART(IC)           istart = INDSTART(IC)
1061           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
1062           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
1063           COG = 0             COG = 0  
1064             SGN = 0.
1065           mu  = 0           mu  = 0
1066           do i = istart,istop  c         print*,'-------'
1067             do i = INDMAX(IC),istart,-1
1068                ipos = i-INDMAX(ic)
1069                cut  = incut*CLSIGMA(i)
1070                if(CLSIGNAL(i).ge.cut)then              
1071                   COG = COG + ipos*CLSIGNAL(i)
1072                   SGN = SGN + CLSIGNAL(i)
1073                   mu = mu + 1
1074    c               print*,ipos,CLSIGNAL(i)
1075                else
1076                   goto 10
1077                endif
1078             enddo
1079     10      continue
1080             do i = INDMAX(IC)+1,istop
1081              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
1082              cut  = incut*CLSIGMA(i)              cut  = incut*CLSIGMA(i)
1083              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
1084                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1085                   SGN = SGN + CLSIGNAL(i)
1086                 mu = mu + 1                 mu = mu + 1
1087    c               print*,ipos,CLSIGNAL(i)
1088                else
1089                   goto 20
1090              endif              endif
1091           enddo           enddo
1092           if(DEDX(ic).le.0)then   20      continue
1093              print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)           if(SGN.le.0)then
1094                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
1095              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1096              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
1097              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
1098           else           else
1099              COG=COG/DEDX(ic)              COG=COG/SGN
1100           endif           endif
1101    c         print*,'-------'
1102                    
1103        else        else
1104                    
# Line 730  ' Line 1111  '
1111    
1112  c      print *,'## cog ',ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1113    
1114          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1115             if(DEBUG.eq.1)
1116         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1117             if(DEBUG.eq.1)
1118         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1119          endif
1120    
1121    
1122        return        return
1123        end        end
1124    
1125  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1126    
1127        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
1128  *-------------------------------------------------------  *-------------------------------------------------------
1129  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 749  c      print *,'## cog ',ncog,ic,cog,'// Line 1139  c      print *,'## cog ',ncog,ic,cog,'//
1139        include 'calib.f'        include 'calib.f'
1140    
1141        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1142           f  = 4.           si = 8.4  !average good-strip noise
1143           si = 8.4           f  = 4.   !average bad-strip noise: f*si
1144           incut=incuty           incut=incuty
1145        else                      !X-view        else                      !X-view
1146           f  = 6.           si = 3.9  !average good-strip noise
1147           si = 3.9           f  = 6.   !average bad-strip noise: f*si
1148           incut=incutx           incut=incutx
1149        endif        endif
1150                
# Line 765  c      print *,'## cog ',ncog,ic,cog,'// Line 1155  c      print *,'## cog ',ncog,ic,cog,'//
1155  *     --> signal of the central strip  *     --> signal of the central strip
1156           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1157           fsc = 1           fsc = 1
1158           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
1159             fsc = clsigma(INDMAX(ic))/si
1160  *     --> signal of adjacent strips  *     --> signal of adjacent strips
1161           sl1 = 0                !left 1           sl1 = 0                !left 1
1162           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 773  c      print *,'## cog ',ncog,ic,cog,'// Line 1164  c      print *,'## cog ',ncog,ic,cog,'//
1164       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1165       $        )then       $        )then
1166              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
1167              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
1168  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
1169           endif           endif
1170    
1171           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 784  c            fsl1 = 0 Line 1174  c            fsl1 = 0
1174       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1175       $        )then       $        )then
1176              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
1177              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
1178  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
1179           endif           endif
1180           sr1 = 0                !right 1           sr1 = 0                !right 1
1181           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 796  c            fsl2 = 0 Line 1185  c            fsl2 = 0
1185       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1186       $        )then       $        )then
1187              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
1188              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
1189  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
1190           endif               endif    
1191           sr2 = 0                !right 2           sr2 = 0                !right 2
1192           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 808  c            fsr1 = 0 Line 1196  c            fsr1 = 0
1196       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1197       $        )then       $        )then
1198              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
1199              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
1200  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
1201           endif           endif
1202    
1203    
1204    
1205  ************************************************************  ************************************************************
1206  *     COG computation  *     COG2-3-4 computation
1207  ************************************************************  ************************************************************
1208    
1209  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
1210                    
1211           COG = 0.           vCOG = cog(ncog,ic)!0.
1212                    
1213           if(ncog.eq.2)then           if(ncog.eq.2)then
1214              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1215                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
1216                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1217                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1218              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
1219                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
1220                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1221                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1222              endif              endif
1223           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1224              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
1225              fbad_cog =              fbad_cog =
1226       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1227              fbad_cog =              fbad_cog =
1228       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1229           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1230              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1231                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1232                 fbad_cog =                 fbad_cog =
1233       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1234       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1235                 fbad_cog =                 fbad_cog =
1236       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1237       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
1238              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
1239                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1240                 fbad_cog =                 fbad_cog =
1241       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1242       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1243                 fbad_cog =                 fbad_cog =
1244       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1245       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1246              endif              endif
1247           else           else
1248              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1249              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1250              COG = 0.  c            COG = 0.
1251           endif           endif
1252                    
1253        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1254    *     =========================
1255    *     COG computation
1256    *     =========================
1257                    
1258           iv=VIEW(ic)           vCOG = cog(0,ic)
1259           istart=INDSTART(IC)  
1260           istop=TOTCLLENGTH           iv     = VIEW(ic)
1261           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1262           COG=0.           istop  = TOTCLLENGTH
1263           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1264           SDE=0.           SGN = 0.
1265           do i=istart,istop           SNU = 0.
1266              ipos=i-INDMAX(ic)           SDE = 0.
1267              il=nvk(MAXS(ic)+ipos)  
1268              is=nst(MAXS(ic)+ipos)           do i=INDMAX(IC),istart,-1
1269              cut=incut*SIGMA(iv,il,is)              ipos = i-INDMAX(ic)
1270                cut  = incut*CLSIGMA(i)
1271              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1272                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1273              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1274                   SDE = SDE + (ipos-vCOG)**2
1275                else
1276                   goto 10
1277                endif            
1278           enddo           enddo
1279           COG=COG/DEDX(ic)   10      continue
1280           do i=istart,istop           do i=INDMAX(IC)+1,istop
1281              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1282              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1283              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1284                 fs=1                 fs = clsigma(i)/si
1285                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1286                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1287                 SDE = SDE + (ipos-COG)**2              else
1288                   goto 20
1289              endif                          endif            
1290           enddo           enddo
1291           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1292             if(SDE.ne.0)then
1293                FBAD_COG=SNU/SDE
1294             else
1295                
1296             endif
1297    
1298        else        else
1299                                        
# Line 913  c      print*,sl2,sl1,sc,sr1,sr2 Line 1311  c      print*,sl2,sl1,sc,sr1,sr2
1311        end        end
1312    
1313    
1314  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1315        real function fbad_cog0(ncog,ic)  c$$$      real function fbad_cog0(ncog,ic)
1316  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1317  *     this function returns a factor that takes into  c$$$*     this function returns a factor that takes into
1318  *     account deterioration of the spatial resolution  c$$$*     account deterioration of the spatial resolution
1319  *     in the case BAD strips are included in the cluster.  c$$$*     in the case BAD strips are included in the cluster.
1320  *     This factor should multiply the nominal spatial  c$$$*     This factor should multiply the nominal spatial
1321  *     resolution.  c$$$*     resolution.
1322  *  c$$$*
1323  *     NB!!!  c$$$*     NB!!!
1324  *     (this is the old version. It consider only the two  c$$$*     (this is the old version. It consider only the two
1325  *     strips with the greatest signal. The new one is  c$$$*     strips with the greatest signal. The new one is
1326  *     fbad_cog(ncog,ic) )  c$$$*     fbad_cog(ncog,ic) )
1327  *      c$$$*    
1328  *-------------------------------------------------------  c$$$*-------------------------------------------------------
1329    c$$$
1330        include 'commontracker.f'  c$$$      include 'commontracker.f'
1331        include 'level1.f'  c$$$      include 'level1.f'
1332        include 'calib.f'  c$$$      include 'calib.f'
1333    c$$$
1334  *     --> signal of the central strip  c$$$*     --> signal of the central strip
1335        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
1336    c$$$
1337  *     signal of adjacent strips  c$$$*     signal of adjacent strips
1338  *     --> left  c$$$*     --> left
1339        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
1340        if(  c$$$      if(
1341       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
1342       $     )  c$$$     $     )
1343       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
1344    c$$$
1345        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
1346        if(  c$$$      if(
1347       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
1348       $     )  c$$$     $     )
1349       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
1350    c$$$
1351  *     --> right  c$$$*     --> right
1352        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
1353        if(  c$$$      if(
1354       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1355       $     .or.  c$$$     $     .or.
1356       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1357       $     )  c$$$     $     )
1358       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
1359    c$$$
1360        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
1361        if(  c$$$      if(
1362       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1363       $     .or.  c$$$     $     .or.
1364       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1365       $     )  c$$$     $     )
1366       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
1367    c$$$
1368    c$$$
1369        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1370           f  = 4.  c$$$         f  = 4.
1371           si = 8.4  c$$$         si = 8.4
1372        else                              !X-view  c$$$      else                              !X-view
1373           f  = 6.  c$$$         f  = 6.
1374           si = 3.9  c$$$         si = 3.9
1375        endif  c$$$      endif
1376    c$$$
1377        fbad_cog = 1.  c$$$      fbad_cog = 1.
1378        f0 = 1  c$$$      f0 = 1
1379        f1 = 1  c$$$      f1 = 1
1380        f2 = 1  c$$$      f2 = 1
1381        f3 = 1    c$$$      f3 = 1  
1382        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
1383            c$$$        
1384           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
1385           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
1386  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
1387    c$$$
1388           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
1389              fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
1390           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
1391              fbad_cog = 1.  c$$$            fbad_cog = 1.
1392           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
1393              fbad_cog = 1.  c$$$            fbad_cog = 1.
1394           else  c$$$         else
1395              fbad_cog = 1.  c$$$            fbad_cog = 1.
1396           endif  c$$$         endif
1397            c$$$        
1398        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
1399    c$$$
1400    c$$$
1401           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f
1402           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
1403  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
1404    c$$$
1405           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
1406              fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)  c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
1407           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
1408              fbad_cog = 1.  c$$$            fbad_cog = 1.
1409           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
1410              fbad_cog = 1.  c$$$            fbad_cog = 1.
1411           else  c$$$         else
1412              fbad_cog = 1.  c$$$            fbad_cog = 1.
1413           endif  c$$$         endif
1414    c$$$
1415        endif  c$$$      endif
1416    c$$$
1417        fbad_cog0 = sqrt(fbad_cog)  c$$$      fbad_cog0 = sqrt(fbad_cog)
1418    c$$$
1419        return  c$$$      return
1420        end  c$$$      end
1421    c$$$
1422    c$$$
1423    c$$$
1424    
1425  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1426    
1427        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1428    
1429        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1430        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1098  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1496  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1496       +/       +/
1497    
1498        V(1)= abs(x)        V(1)= abs(x)
1499          if(V(1).gt.20.)V(1)=20.
1500    
1501        HQUADF = 0.        HQUADF = 0.
1502        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1112  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1511  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1511     20 CONTINUE     20 CONTINUE
1512        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1513                
1514        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1515    
1516        END        END
1517    
1518  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1519        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1520        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1521        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1522        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1188  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1587  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1587       +/       +/
1588    
1589        V(1) =  abs(x)        V(1) =  abs(x)
1590          if(V(1).gt.20.)V(1)=20.
1591    
1592        HQUADF = 0.        HQUADF = 0.
1593        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1594           HQDJ = 0.           HQDJ = 0.
# Line 1201  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1602  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1602     20 CONTINUE     20 CONTINUE
1603        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1604                
1605        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1606    
1607        END        END
1608  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1609        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1610        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1611        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1612        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1276  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1677  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1677       +/       +/
1678    
1679        V(1)=abs(x)        V(1)=abs(x)
1680          if(V(1).gt.20.)V(1)=20.
1681    
1682        HQUADF = 0.        HQUADF = 0.
1683        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1290  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1692  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1692     20 CONTINUE     20 CONTINUE
1693        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1694                
1695        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1696    
1697        END        END
1698  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1699        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1700        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1701        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1702        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1347  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1749  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1749       +/       +/
1750    
1751        v(1)= abs(x)        v(1)= abs(x)
1752          if(V(1).gt.20.)V(1)=20.
1753                
1754        HQUADF = 0.        HQUADF = 0.
1755        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1361  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1764  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1764     20 CONTINUE     20 CONTINUE
1765        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1766    
1767        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1768    
1769        END        END
1770  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1413  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1816  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1816       +/       +/
1817                
1818        V(1)=abs(x)        V(1)=abs(x)
1819          if(V(1).gt.20.)V(1)=20.
1820    
1821        HQUADF = 0.        HQUADF = 0.
1822        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1493  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1897  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1897       +/       +/
1898    
1899        V(1)=abs(x)        V(1)=abs(x)
1900          if(V(1).gt.20.)V(1)=20.
1901    
1902        HQUADF = 0.        HQUADF = 0.
1903        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1510  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1915  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1915        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1916    
1917        END        END
1918    
1919    
1920  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1921          real function pfacorr(ic,angle)
1922    *--------------------------------------------------------------
1923    *     this function returns the landi correction for this cluster
1924    *--------------------------------------------------------------
1925          include 'commontracker.f'
1926          include 'calib.f'
1927          include 'level1.f'
1928    
1929          real angle
1930          integer iview,lad
1931    
1932          iview = VIEW(ic)            
1933          lad = nld(MAXS(ic),VIEW(ic))
1934    
1935    *     find angular bin
1936    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1937          do iang=1,nangbin
1938             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1939                iangle=iang
1940                goto 98
1941             endif
1942          enddo
1943          if(DEBUG.eq.1)
1944         $     print*,'pfacorr *** warning *** angle out of range: ',angle
1945          if(angle.le.angL(1))iang=1
1946          if(angle.ge.angR(nangbin))iang=nangbin
1947     98   continue                  !jump here if ok
1948    
1949          pfacorr = fcorr(iview,lad,iang)
1950    
1951          if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr
1952    
1953    
1954     100  return
1955          end
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23