/[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.3 by pam-fi, Thu Sep 28 14:04:40 2006 UTC revision 1.27 by pam-fi, Tue Mar 12 11:02:02 2013 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    *     -----------------------------------------------------------------
12    *     p.f.a.
13    *     -----------------------------------------------------------------
14    *     real function pfaeta(ic,angle)
15    *     real function pfaetal(ic,angle)
16    *     real function pfaeta2(ic,angle)
17    *     real function pfaeta3(ic,angle)
18    *     real function pfaeta4(ic,angle)
19    *     real function cog(ncog,ic)
20    *
21    *     -----------------------------------------------------------------
22    *     risoluzione spaziale media, stimata dalla simulazione (samuele)
23    *     -----------------------------------------------------------------
24    *     FUNCTION risxeta2(angle)
25    *     FUNCTION risxeta3(angle)
26    *     FUNCTION risxeta4(angle)
27    *     FUNCTION risyeta2(angle)
28    *     FUNCTION risy_cog(angle)
29    *     FUNCTION risx_cog(angle)
30    *     real function riseta(iview,angle)
31    *     -----------------------------------------------------------------
32    *     fattore moltiplicativo per tenere conto della dipendenza della
33    *     risoluzione dal rumore delle strip
34    *     -----------------------------------------------------------------
35    *     real function fbad_cog(ncog,ic)
36    *     real function fbad_eta(ic,angle)
37    *
38    *     -----------------------------------------------------------------
39    *     NUOVO APPROCCIO PER LA STIMA DELLA RISOLUZIONE
40    *     -----------------------------------------------------------------
41    *     real function riscogtheor(ncog,ic)
42    *     real function risetatheor(ncog,ic,angle)
43    *
44    *     -----------------------------------------------------------------
45    *     correzione landi
46    *     -----------------------------------------------------------------
47    *     real function pfacorr(ic,angle)
48    *
49    *     real function effectiveangle(ang,iview,bbb)
50    *     real function fieldcorr(iview,bbb)
51    *
52    *     NB - The angle is the "effective angle", which is relative
53    *          to the sensor and it takes into account the magnetic field
54  *  *
55  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
56    
57          subroutine idtoc(ipfa,cpfa)
58          
59          integer ipfa
60          character*10 cpfa
61    
62          CPFA='COG4'
63          if(ipfa.eq.0)CPFA='ETA'
64          if(ipfa.eq.2)CPFA='ETA2'
65          if(ipfa.eq.3)CPFA='ETA3'
66          if(ipfa.eq.4)CPFA='ETA4'
67          if(ipfa.eq.5)CPFA='ETAL'
68          if(ipfa.eq.10)CPFA='COG'
69          if(ipfa.eq.11)CPFA='COG1'
70          if(ipfa.eq.12)CPFA='COG2'
71          if(ipfa.eq.13)CPFA='COG3'
72          if(ipfa.eq.14)CPFA='COG4'
73          
74          end
75    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
76          real function effectiveangle(ang,iview,bbb)
77    
78          include 'commontracker.f'
79    
80          effectiveangle = 0.
81    
82          if(mod(iview,2).eq.0)then
83    c     =================================================
84    c     X view
85    c     =================================================
86    c     here bbb is the y component of the m.field
87             angx = ang
88             by   = bbb
89             if(iview.eq.12) angx = -1. * ang
90             if(iview.eq.12) by   = -1. * bbb
91    cc         tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!!
92             tgtemp  = tan(angx*acos(-1.)/180.) + pmuH_h*by*0.00001
93    
94          elseif(mod(iview,2).eq.1)then
95    c     =================================================
96    c     Y view
97    c     =================================================        
98    c     here bbb is the x component of the m.filed
99             angy = ang
100             bx   = bbb
101             tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001        
102    
103          endif      
104          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
105    
106          return
107          end
108    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
109          real function fieldcorr(iview,bbb)
110    
111          include 'commontracker.f'
112    
113          fieldcorr = 0.
114    
115          if(mod(iview,2).eq.0)then
116    
117    c     =================================================
118    c     X view
119    c     =================================================
120    c     here bbb is the y component of the m.field
121             by   = bbb
122             if(iview.eq.12) by = -1. * bbb
123             fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
124    
125          elseif(mod(iview,2).eq.1)then
126    c     =================================================
127    c     Y view
128    c     =================================================        
129    c     here bbb is the x component of the m.filed
130             bx   = bbb
131             fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
132    
133          endif      
134          
135          return
136          end
137    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
138    
139          subroutine applypfa(PFAtt,ic,ang,corr,res)
140    *---------------------------------------------------------------
141    *     this subroutine calculate the coordinate of cluster ic (in
142    *     strip units), relative to the strip with the maximum signal,
143    *     and its spatial resolution (in cm), applying PFAtt.
144    *     ang is the effective angle, relative to the sensor
145    *---------------------------------------------------------------
146    
147          character*4 PFAtt
148          include 'commontracker.f'
149          include 'level1.f'
150    
151          corr = 0
152          res  = 0
153    
154          if(ic.le.0)return
155    
156          iview   = VIEW(ic)
157    
158          if(mod(iview,2).eq.0)then
159    c     =================================================
160    c     X view
161    c     =================================================
162    
163             res = RESXAV
164    
165             if(PFAtt.eq.'COG1')then
166    
167                corr   = 0
168                res = 1e-4*pitchX/sqrt(12.)!!res
169    
170             elseif(PFAtt.eq.'COG2')then
171    
172                corr    = cog(2,ic)            
173                res = risx_cog(abs(ang))!TEMPORANEO              
174                res = res*fbad_cog(2,ic)
175    
176             elseif(PFAtt.eq.'COG3')then
177    
178                corr    = cog(3,ic)            
179                res = risx_cog(abs(ang))!TEMPORANEO                      
180                res = res*fbad_cog(3,ic)
181    
182             elseif(PFAtt.eq.'COG4')then
183    
184                corr    = cog(4,ic)            
185                res = risx_cog(abs(ang))!TEMPORANEO                      
186                res = res*fbad_cog(4,ic)
187    
188             elseif(PFAtt.eq.'ETA2')then
189    
190                corr  = pfaeta2(ic,ang)          
191                res = risxeta2(abs(ang))
192                res = res*fbad_cog(2,ic)
193    
194             elseif(PFAtt.eq.'ETA3')then                        
195    
196                corr  = pfaeta3(ic,ang)          
197                res = risxeta3(abs(ang))                      
198                res = res*fbad_cog(3,ic)              
199    
200             elseif(PFAtt.eq.'ETA4')then                        
201    
202                corr  = pfaeta4(ic,ang)            
203                res = risxeta4(abs(ang))                      
204                res = res*fbad_cog(4,ic)              
205    
206             elseif(PFAtt.eq.'ETA')then  
207    
208                corr  = pfaeta(ic,ang)            
209    c            res = riseta(ic,ang)                    
210                res = riseta(iview,ang)                    
211                res = res*fbad_eta(ic,ang)            
212    
213             elseif(PFAtt.eq.'ETAL')then  
214    
215                corr  = pfaetal(ic,ang)            
216                res = riseta(iview,ang)                    
217                res = res*fbad_eta(ic,ang)            
218    
219             elseif(PFAtt.eq.'COG')then          
220    
221                corr  = cog(0,ic)            
222                res = risx_cog(abs(ang))                    
223                res = res*fbad_cog(0,ic)
224    
225             else
226                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
227             endif
228    
229    
230    *     ======================================
231    *     temporary patch for saturated clusters
232    *     ======================================
233             if( nsatstrips(ic).gt.0 )then
234    c            corr  = cog(4,ic)            
235                corr = digsat(ic)
236                res = pitchX*1e-4/sqrt(12.)
237    cc            cc=cog(4,ic)
238    c$$$            print*,ic,' *** ',cc
239    c$$$            print*,ic,' *** ',res
240             endif
241    
242    
243          elseif(mod(iview,2).eq.1)then
244    c     =================================================
245    c     Y view
246    c     =================================================
247    
248             res = RESYAV
249    
250             if(PFAtt.eq.'COG1')then
251    
252                corr  = 0  
253                res = 1e-4*pitchY/sqrt(12.)!res  
254    
255             elseif(PFAtt.eq.'COG2')then
256    
257                corr    = cog(2,ic)
258                res = risy_cog(abs(ang))!TEMPORANEO
259                res = res*fbad_cog(2,ic)
260    
261             elseif(PFAtt.eq.'COG3')then
262    
263                corr    = cog(3,ic)
264                res = risy_cog(abs(ang))!TEMPORANEO
265                res = res*fbad_cog(3,ic)
266    
267             elseif(PFAtt.eq.'COG4')then
268    
269                corr    = cog(4,ic)
270                res = risy_cog(abs(ang))!TEMPORANEO
271                res = res*fbad_cog(4,ic)
272    
273             elseif(PFAtt.eq.'ETA2')then
274    
275                corr    = pfaeta2(ic,ang)          
276                res = risyeta2(abs(ang))              
277                res = res*fbad_cog(2,ic)
278    
279             elseif(PFAtt.eq.'ETA3')then                      
280    
281                corr    = pfaeta3(ic,ang)
282                res = res*fbad_cog(3,ic)  
283    
284             elseif(PFAtt.eq.'ETA4')then  
285    
286                corr    = pfaeta4(ic,ang)
287                res = res*fbad_cog(4,ic)
288    
289             elseif(PFAtt.eq.'ETA')then
290    
291                corr    = pfaeta(ic,ang)
292    c            res = riseta(ic,ang)  
293                res = riseta(iview,ang)  
294                res = res*fbad_eta(ic,ang)
295    
296             elseif(PFAtt.eq.'ETAL')then
297    
298                corr    = pfaetal(ic,ang)
299                res = riseta(iview,ang)  
300                res = res*fbad_eta(ic,ang)
301    
302             elseif(PFAtt.eq.'COG')then
303    
304                corr    = cog(0,ic)            
305                res = risy_cog(abs(ang))
306                res = res*fbad_cog(0,ic)
307    
308             else
309                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
310             endif
311    
312    
313    *     ======================================
314    *     temporary patch for saturated clusters
315    *     ======================================
316             if( nsatstrips(ic).gt.0 )then
317    c            corr    = cog(4,ic)            
318                corr = digsat(ic)
319                res = pitchY*1e-4/sqrt(12.)
320    cc            cc=cog(4,ic)
321    c$$$            print*,ic,' *** ',cc
322    c$$$            print*,ic,' *** ',res
323             endif
324            
325          endif
326          end
327    
328  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
329        real function pfa_eta(ic,angle)        integer function npfastrips(ic,angle)
330    *--------------------------------------------------------------
331    *     thid function returns the number of strips used
332    *     to evaluate the position of a cluster, according to the p.f.a.
333    *--------------------------------------------------------------
334          include 'commontracker.f'
335          include 'level1.f'
336          include 'calib.f'
337    
338          character*4 usedPFA
339          
340    
341    
342          call idtoc(pfaid,usedPFA)
343    
344          npfastrips=-1
345    
346          if(usedPFA.eq.'COG1')npfastrips=1
347          if(usedPFA.eq.'COG2')npfastrips=2
348          if(usedPFA.eq.'COG3')npfastrips=3
349          if(usedPFA.eq.'COG4')npfastrips=4
350          if(usedPFA.eq.'ETA2')npfastrips=2
351          if(usedPFA.eq.'ETA3')npfastrips=3
352          if(usedPFA.eq.'ETA4')npfastrips=4
353    *     ----------------------------------------------------------------
354          if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
355    c         print*,VIEW(ic),angle
356             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
357                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
358                   npfastrips=2
359                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
360                   npfastrips=3
361                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
362                   npfastrips=4
363                else
364                   npfastrips=4     !COG4
365                endif                        
366             else                   !X-view
367                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
368                   npfastrips=2
369                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
370                   npfastrips=3
371                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
372                   npfastrips=4
373                else
374                   npfastrips=4     !COG4
375                endif                        
376             endif
377          endif
378    *     ----------------------------------------------------------------
379          if(usedPFA.eq.'COG')then
380    
381             npfastrips=0
382    
383    c$$$         iv=VIEW(ic)
384    c$$$         if(mod(iv,2).eq.1)incut=incuty
385    c$$$         if(mod(iv,2).eq.0)incut=incutx
386    c$$$         istart = INDSTART(IC)
387    c$$$         istop  = TOTCLLENGTH
388    c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
389    c$$$         mu  = 0
390    c$$$         do i = INDMAX(IC),istart,-1
391    c$$$            ipos = i-INDMAX(ic)
392    c$$$            cut  = incut*CLSIGMA(i)
393    c$$$            if(CLSIGNAL(i).ge.cut)then
394    c$$$               mu = mu + 1
395    c$$$               print*,i,mu
396    c$$$            else
397    c$$$               goto 10
398    c$$$            endif
399    c$$$         enddo
400    c$$$ 10      continue
401    c$$$         do i = INDMAX(IC)+1,istop
402    c$$$            ipos = i-INDMAX(ic)
403    c$$$            cut  = incut*CLSIGMA(i)
404    c$$$            if(CLSIGNAL(i).ge.cut)then
405    c$$$               mu = mu + 1
406    c$$$               print*,i,mu
407    c$$$            else
408    c$$$               goto 20
409    c$$$            endif
410    c$$$         enddo
411    c$$$ 20      continue
412    c$$$         npfastrips=mu
413    
414          endif
415    *     ----------------------------------------------------------------
416    
417    c      print*,pfaid,usedPFA,angle,npfastrips
418    
419          return
420          end
421    
422    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
423          real function pfaeta(ic,angle)
424  *--------------------------------------------------------------  *--------------------------------------------------------------
425  *     this function returns the position (in strip units)  *     this function returns the position (in strip units)
426  *     it calls:  *     it calls:
427  *     - pfa_eta2(ic,angle)  *     - pfaeta2(ic,angle)
428  *     - pfa_eta3(ic,angle)  *     - pfaeta3(ic,angle)
429  *     - pfa_eta4(ic,angle)  *     - pfaeta4(ic,angle)
430  *     according to the angle  *     according to the angle
431  *--------------------------------------------------------------  *--------------------------------------------------------------
432        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
433        include 'level1.f'        include 'level1.f'
434          include 'calib.f'
435                
436        pfa_eta = 0        pfaeta = 0
437    
438        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
439                
440           pfa_eta = pfa_eta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
441                  pfaeta = pfaeta2(ic,angle)
442    cc            print*,pfaeta2(ic,angle)
443             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
444                pfaeta = pfaeta3(ic,angle)
445             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
446                pfaeta = pfaeta4(ic,angle)
447             else
448                pfaeta = cog(4,ic)
449             endif            
450    
451        else                      !X-view        else                      !X-view
452    
453           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
454              pfa_eta = pfa_eta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
455           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
456              pfa_eta = pfa_eta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
457           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
458              pfa_eta = pfa_eta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
459           endif           else
460                pfaeta = cog(4,ic)
461             endif            
462                            
463        endif        endif
464          
465    c 100  return
466   100  return        return
467        end        end
468    
469  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
470        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
471  *--------------------------------------------------------------  *--------------------------------------------------------------
472  *     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))  
473  *     it calls:  *     it calls:
474  *     - risx_eta2(angle)  *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
475  *     - risy_eta2(angle)  *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
476  *     - risx_eta3(angle)  *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
 *     - risx_eta4(angle)  
477  *     according to the angle  *     according to the angle
478  *--------------------------------------------------------------  *--------------------------------------------------------------
479        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
480        include 'level1.f'        include 'level1.f'
481          include 'calib.f'
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
482                
483        ris_eta = 0        pfaetal = 0
484    
485        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
486                
487           ris_eta = risy_eta2(angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
488           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
489    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
490             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
491                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
492             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
493                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
494             else
495                pfaetal = cog(4,ic)
496             endif            
497    
498        else                      !X-view        else                      !X-view
499    
500           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
501              ris_eta = risx_eta2(angle)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
502           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then  cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
503              ris_eta = risx_eta3(angle)           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
504           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then              pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
505              ris_eta = risx_eta4(angle)           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
506           elseif(abs(angle).gt.21.)then              pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
507              ris_eta = risx_eta4(21.)           else
508           endif              pfaetal = cog(4,ic)
509             endif            
510                            
511        endif        endif
512          
513    c 100  return
514          return
515          end
516    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
517    c      real function riseta(ic,angle)
518          real function riseta(iview,angle)
519    *--------------------------------------------------------------
520    *     this function returns the average spatial resolution
521    *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
522    *     it calls:
523    *     - risxeta2(angle)
524    *     - risyeta2(angle)
525    *     - risxeta3(angle)
526    *     - risxeta4(angle)
527    *     according to the angle
528    *--------------------------------------------------------------
529          include 'commontracker.f'
530          include 'level1.f'
531          include 'calib.f'
532    
533          riseta = 0
534    
535    c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
536          if(mod(iview,2).eq.1)then !Y-view
537          
538    
539             if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
540                riseta = risyeta2(angle)
541             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
542                riseta = risy_cog(angle) !ATTENZIONE!!
543             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
544                riseta = risy_cog(angle) !ATTENZIONE!!
545             else
546                riseta = risy_cog(angle)
547             endif            
548    
549          else                      !X-view
550    
551  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
552  c$$$     $     ,' -->',ris_eta              riseta = risxeta2(angle)
553             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
554                riseta = risxeta3(angle)
555             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
556                riseta = risxeta4(angle)
557             else
558                riseta = risx_cog(angle)
559             endif            
560                
561          endif
562    
563    
564   100  return  c 100  return
565          return
566        end        end
567    
568  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 99  c$$$     $     ,' -->',ris_eta Line 575  c$$$     $     ,' -->',ris_eta
575  *     resolution.  *     resolution.
576  *     It calls the function FBAD_COG(NCOG,IC),  *     It calls the function FBAD_COG(NCOG,IC),
577  *     accordingto the angle  *     accordingto the angle
578    *
579    *     >>> cosi` non e` corretto!!
580    *     >>> l'errore sulla coordinata eta si ottiene moltiplicando
581    *     >>> l'errore sulla coordinata cog per la derivata della
582    *     >>> distribuzione eta... pur sapendolo l'ho sempre ignorato...
583    *     >>> deve essere modificato!!!!
584    *
585  *-------------------------------------------------------  *-------------------------------------------------------
586    
587        include 'commontracker.f'        include 'commontracker.f'
588        include 'level1.f'        include 'level1.f'
589  *      include 'calib.f'        include 'calib.f'
590        fbad_eta = 0        fbad_eta = 0
591    
592        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
593                
594           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
595                fbad_eta = fbad_cog(2,ic)
596             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
597                fbad_eta = fbad_cog(3,ic)
598             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
599                fbad_eta = fbad_cog(4,ic)
600             else
601                fbad_eta = fbad_cog(4,ic)
602             endif            
603    
604        else                      !X-view        else                      !X-view
605    
606           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
607              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
608           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
609              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
610           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
611              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
612           endif           else
613                fbad_eta = fbad_cog(4,ic)
614             endif            
615                            
616        endif        endif
617    
# Line 126  c$$$     $     ,' -->',ris_eta Line 619  c$$$     $     ,' -->',ris_eta
619        end        end
620    
621  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
622  c*****************************************************        real function pfaeta2(ic,angle)
 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)  
623  *--------------------------------------------------------------  *--------------------------------------------------------------
624  *     this function returns  *     this function returns
625  *  *
# Line 149  c      real function pfa_eta2(cog2,view, Line 638  c      real function pfa_eta2(cog2,view,
638        real cog2,angle        real cog2,angle
639        integer iview,lad        integer iview,lad
640    
641  c      logical DEBUG        iview   = VIEW(ic)            
642  c      common/dbg/DEBUG        lad     = nld(MAXS(ic),VIEW(ic))
643          cog2    = cog(2,ic)          
644        iview = VIEW(ic)              !(1)        pfaeta2 = cog2
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
       pfa_eta2=cog2  
645    
646    *     ----------------
647  *     find angular bin  *     find angular bin
648    *     ----------------
649  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
650        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
651           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
652              iangle=iang              iangle=iang
653              goto 98              goto 98
654           endif           endif
655        enddo        enddo
656        if(DEBUG)        if(DEBUG.EQ.1)
657       $     print*,'pfa_eta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
658        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
659        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
660   98   continue                  !jump here if ok   98   continue                  !jump here if ok
661    
662    *     -------------
663    *     within +/-0.5
664    *     -------------
665    
666  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*,'pfa_eta2 *** 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  
   
667        iadd=0        iadd=0
668   10   continue   10   continue
669        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
670           cog2 = cog2 + 1           cog2 = cog2 + 1
671           iadd = iadd + 1           iadd = iadd + 1
672             if(iadd>iaddmax)goto 111
673           goto 10           goto 10
674        endif        endif
675   20   continue   20   continue
676        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
677           cog2 = cog2 - 1           cog2 = cog2 - 1
678           iadd = iadd - 1           iadd = iadd - 1
679             if(iadd<-1*iaddmax)goto 111
680           goto 20           goto 20
681        endif        endif
682          goto 1111
683     111  continue
684          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
685          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
686          cog2=0
687     1111 continue
688    
689  *     --------------------------------  *     --------------------------------
690  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 230  c            print*,'-----',x1,x2,y1,y2 Line 707  c            print*,'-----',x1,x2,y1,y2
707        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
708        BB=y1-AA*x1        BB=y1-AA*x1
709    
710        pfa_eta2 = AA*cog2+BB        pfaeta2 = AA*cog2+BB
711        pfa_eta2 = pfa_eta2 - iadd        pfaeta2 = pfaeta2 - iadd
712    
713  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
714  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
715  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
716  c$$$      endif  c$$$      endif
717  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
718  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
719  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
720  c$$$      endif  c$$$      endif
721    
722        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
723       $     ,cog2-iadd,' -->',pfa_eta2       $     ,cog2-iadd,' -->',pfaeta2
724    
725    
726   100  return  c 100  return
727          return
728        end        end
729    
730  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
731  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)  
732  *--------------------------------------------------------------  *--------------------------------------------------------------
733  *     this function returns  *     this function returns
734  *  *
# Line 273  c      real function pfa_eta3(cog3,view, Line 747  c      real function pfa_eta3(cog3,view,
747        real cog3,angle        real cog3,angle
748        integer iview,lad        integer iview,lad
749    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
750    
751          iview = VIEW(ic)            
752          lad = nld(MAXS(ic),VIEW(ic))
753          cog3 = cog(3,ic)
754          cc = cog3
755          cog3 = cc
756          pfaeta3=cog3
757    
758        iview = VIEW(ic)              !(1)  *     ----------------
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog3 = cog(3,ic)             !(1)  
       pfa_eta3=cog3  
   
759  *     find angular bin  *     find angular bin
760    *     ----------------
761  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
762        do iang=1,nangbin        do iang=1,nangbin
763  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 291  c         print*,'~~~~~~~~~~~~ ',iang,an Line 766  c         print*,'~~~~~~~~~~~~ ',iang,an
766              goto 98              goto 98
767           endif           endif
768        enddo        enddo
769        if(DEBUG)        if(DEBUG.EQ.1)
770       $     print*,'pfa_eta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
771        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
772        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
773   98   continue                  !jump here if ok   98   continue                  !jump here if ok
774    
775    *     -------------
776    *     within +/-0.5
777    *     -------------
778    
779  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  
   
         
780        iadd=0        iadd=0
781   10   continue   10   continue
782        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
783           cog3 = cog3 + 1           cog3   = cog3 + 1.
784           iadd = iadd + 1           iadd = iadd + 1
785             if(iadd>iaddmax) goto 111
786           goto 10           goto 10
787        endif        endif
788   20   continue   20   continue
789        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
790           cog3 = cog3 - 1           cog3 = cog3 - 1.
791           iadd = iadd - 1           iadd = iadd - 1
792             if(iadd<-1*iaddmax) goto 111
793           goto 20           goto 20
794        endif        endif
795          goto 1111
796     111  continue
797          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
798          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
799          cog3=0      
800     1111 continue
801    
802  *     --------------------------------  *     --------------------------------
803  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 354  c            print*,'-----',x1,x2,y1,y2 Line 820  c            print*,'-----',x1,x2,y1,y2
820        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
821        BB=y1-AA*x1        BB=y1-AA*x1
822    
823        pfa_eta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
824        pfa_eta3 = pfa_eta3 - iadd        pfaeta3 = pfaeta3 - iadd
825    
 c$$$      if(iflag.eq.1)then  
 c$$$         pfa_eta2=pfa_eta2-1.   !temp  
 c$$$         cog2=cog2-1.           !temp  
 c$$$      endif  
 c$$$      if(iflag.eq.-1)then  
 c$$$         pfa_eta2=pfa_eta2+1.   !temp  
 c$$$         cog2=cog2+1.           !temp  
 c$$$      endif  
826    
827        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
828       $     ,cog3-iadd,' -->',pfa_eta3       $     ,cog3-iadd,' -->',pfaeta3
829    
830   100  return  c 100  return
831          return
832        end        end
833    
834  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
835  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)  
836  *--------------------------------------------------------------  *--------------------------------------------------------------
837  *     this function returns  *     this function returns
838  *  *
# Line 396  c      real function pfa_eta4(cog4,view, Line 851  c      real function pfa_eta4(cog4,view,
851        real cog4,angle        real cog4,angle
852        integer iview,lad        integer iview,lad
853    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
854    
855        iview = VIEW(ic)              !(1)        iview = VIEW(ic)            
856        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
857        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
858        pfa_eta4=cog4        pfaeta4=cog4
859    
860    *     ----------------
861  *     find angular bin  *     find angular bin
862    *     ----------------
863  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
864        do iang=1,nangbin        do iang=1,nangbin
865  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 413  c         print*,'~~~~~~~~~~~~ ',iang,an Line 868  c         print*,'~~~~~~~~~~~~ ',iang,an
868              goto 98              goto 98
869           endif           endif
870        enddo        enddo
871        if(DEBUG)        if(DEBUG.EQ.1)
872       $     print*,'pfa_eta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
873        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
874        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
875   98   continue                  !jump here if ok   98   continue                  !jump here if ok
876    
877    *     -------------
878    *     within +/-0.5
879    *     -------------
880    
881  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  
   
         
882        iadd=0        iadd=0
883   10   continue   10   continue
884        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
885           cog4 = cog4 + 1           cog4 = cog4 + 1
886           iadd = iadd + 1           iadd = iadd + 1
887             if(iadd>iaddmax)goto 111
888           goto 10           goto 10
889        endif        endif
890   20   continue   20   continue
891        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
892           cog4 = cog4 - 1           cog4 = cog4 - 1
893           iadd = iadd - 1           iadd = iadd - 1
894             if(iadd<-1*iaddmax)goto 111
895           goto 20           goto 20
896        endif        endif
897          goto 1111
898     111  continue
899          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
900          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
901          cog4=0
902     1111 continue
903    
904  *     --------------------------------  *     --------------------------------
905  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 476  c            print*,'-----',x1,x2,y1,y2 Line 922  c            print*,'-----',x1,x2,y1,y2
922        AA=(y2-y1)/(x2-x1)        AA=(y2-y1)/(x2-x1)
923        BB=y1-AA*x1        BB=y1-AA*x1
924    
925        pfa_eta4 = AA*cog4+BB        pfaeta4 = AA*cog4+BB
926        pfa_eta4 = pfa_eta4 - iadd        pfaeta4 = pfaeta4 - iadd
927    
928  c$$$      if(iflag.eq.1)then  c$$$      if(iflag.eq.1)then
929  c$$$         pfa_eta2=pfa_eta2-1.   !temp  c$$$         pfaeta2=pfaeta2-1.   !temp
930  c$$$         cog2=cog2-1.           !temp  c$$$         cog2=cog2-1.           !temp
931  c$$$      endif  c$$$      endif
932  c$$$      if(iflag.eq.-1)then  c$$$      if(iflag.eq.-1)then
933  c$$$         pfa_eta2=pfa_eta2+1.   !temp  c$$$         pfaeta2=pfaeta2+1.   !temp
934  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
935  c$$$      endif  c$$$      endif
936    
937        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
938       $     ,cog4-iadd,' -->',pfa_eta4       $     ,cog4-iadd,' -->',pfaeta4
939    
940   100  return  c 100  return
941          return
942        end        end
943    
   
   
944  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
945        real function cog0(ncog,ic)        real function digsat(ic)
946  *-------------------------------------------------  *-------------------------------------------------
947  *     this function returns  *
948  *  *          
 *     - 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)  
949  *-------------------------------------------------  *-------------------------------------------------
   
   
950        include 'commontracker.f'        include 'commontracker.f'
951          include 'calib.f'
952        include 'level1.f'        include 'level1.f'
953                
954  *     --> signal of the central strip        integer nsat
955        sc = CLSIGNAL(INDMAX(ic)) !center        real pitchsat
   
 *     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))  
956                
957  ************************************************************        nsat = 0
958  *     COG computation        pitchsat = 0.
959  ************************************************************        iv=VIEW(ic)              
960          istart = INDSTART(IC)
961  c      print*,sl2,sl1,sc,sr1,sr2        istop  = TOTCLLENGTH
962          if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
963        COG = 0.        do i = INDMAX(IC),istart,-1
964                     if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
965        if(sl1.gt.sr1.and.sl1.gt.0.)then       $        .or.
966                 $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
967           if(ncog.eq.2.and.sl1.ne.0)then              nsat = nsat + 1
968              COG = -sl1/(sl1+sc)                      pitchsat = pitchsat + i - INDMAX(IC)
          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)  
969           else           else
970              COG = 0.              goto 10
971           endif           endif
972                  enddo
973        elseif(sl1.le.sr1.and.sr1.gt.0.)then   10   continue
974          do i = INDMAX(IC)+1,istop
975           if(ncog.eq.2.and.sr1.ne.0)then           if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
976              COG = sr1/(sc+sr1)                   $        .or.
977           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then       $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
978              COG = (sr1-sl1)/(sl1+sc+sr1)              nsat = nsat + 1
979           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then              pitchsat = pitchsat + i - INDMAX(IC)
             COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  
980           else           else
981              COG = 0.              goto 20
982           endif           endif
983          enddo
984        endif   20   continue
985          
986        COG0 = COG        digsat = 0
987          if (nsat.gt.0) digsat = pitchsat / nsat
988  c      print *,ncog,ic,cog,'/////////////'        
   
989        return        return
990        end        end
991    
992    
993  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
994        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 613  c      print *,ncog,ic,cog,'//////////// Line 1010  c      print *,ncog,ic,cog,'////////////
1010        include 'calib.f'        include 'calib.f'
1011        include 'level1.f'        include 'level1.f'
1012                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1013    
1014    
1015        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 625  c      common/dbg/DEBUG Line 1020  c      common/dbg/DEBUG
1020  *     --> signal of the central strip  *     --> signal of the central strip
1021           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1022  *     signal of adjacent strips  *     signal of adjacent strips
1023           sl1 = 0                !left 1           sl1 = -9999.           !left 1
1024           if(           if(
1025       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1026       $        )       $        )
1027       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
1028                    
1029           sl2 = 0                !left 2           sl2 = -9999.           !left 2
1030           if(           if(
1031       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1032       $        )       $        )
1033       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
1034                    
1035           sr1 = 0                !right 1           sr1 = -9999.           !right 1
1036           if(           if(
1037       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1038       $        .or.       $        .or.
# Line 645  c      common/dbg/DEBUG Line 1040  c      common/dbg/DEBUG
1040       $        )       $        )
1041       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
1042                    
1043           sr2 = 0                !right 2           sr2 = -9999.           !right 2
1044           if(           if(
1045       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1046       $        .or.       $        .or.
# Line 655  c      common/dbg/DEBUG Line 1050  c      common/dbg/DEBUG
1050                    
1051           COG = 0.           COG = 0.
1052                    
1053           if(ncog.eq.2)then  c      print *,'## ',sl2,sl1,sc,sr1,sr2
1054    
1055    c     ==============================================================
1056             if(ncog.eq.1)then
1057                COG = 0.
1058             if(sr1.gt.sc)cog=1.
1059             if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1060    c     ==============================================================
1061             elseif(ncog.eq.2)then
1062                COG = 0.
1063              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1064                 COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
1065              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
1066                 COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)
1067              endif           elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1068                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1069         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1070                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1071         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1072             endif
1073    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1074    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
1075    c     ==============================================================
1076           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1077              COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
1078                sss = sc
1079                if( sl1.ne.-9999. )COG = COG-sl1
1080                if( sl1.ne.-9999. )sss = sss+sl1
1081                if( sr1.ne.-9999. )COG = COG+sr1
1082                if( sr1.ne.-9999. )sss = sss+sr1
1083                if(sss.ne.0)COG=COG/sss
1084    
1085    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1086    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1087    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1088    c     ==============================================================
1089           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1090    
1091                COG = 0
1092                sss = sc
1093                if( sl1.ne.-9999. )COG = COG-sl1
1094                if( sl1.ne.-9999. )sss = sss+sl1
1095                if( sr1.ne.-9999. )COG = COG+sr1
1096                if( sr1.ne.-9999. )sss = sss+sr1
1097              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1098                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)                 if((sl2+sss).ne.0)
1099              elseif(sl2.le.sr2)then       $              COG = (COG-2*sl2)/(sl2+sss)
1100                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)              elseif(sl2.lt.sr2)then
1101                   if((sr2+sss).ne.0)
1102         $              COG = (2*sr2+COG)/(sr2+sss)
1103                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1104                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1105         $              .and.(sl2+sss).ne.0 )
1106         $              cog = (cog-2*sl2)/(sl2+sss)
1107                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1108         $              .and.(sr2+sss).ne.0 )
1109         $              cog = (2*sr2+cog)/(sr2+sss)              
1110              endif              endif
1111    c     ==============================================================
1112             elseif(ncog.eq.5)then
1113                COG = 0
1114                sss = sc
1115                if( sl1.ne.-9999. )COG = COG-sl1
1116                if( sl1.ne.-9999. )sss = sss+sl1
1117                if( sr1.ne.-9999. )COG = COG+sr1
1118                if( sr1.ne.-9999. )sss = sss+sr1
1119                if( sl2.ne.-9999. )COG = COG-2*sl2
1120                if( sl2.ne.-9999. )sss = sss+sl2
1121                if( sr2.ne.-9999. )COG = COG+2*sr2
1122                if( sr2.ne.-9999. )sss = sss+sr2
1123                if(sss.ne.0)COG=COG/sss
1124           else           else
1125              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1126              print*,'                 (NCOG must be <= 4)'       $           ,' not implemented'
1127              COG = 0.              COG = 0.
1128           endif           endif
1129    
# Line 685  ' Line 1137  '
1137           iv=VIEW(ic)           iv=VIEW(ic)
1138           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
1139           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
1140           istart = INDSTART(IC)           istart = INDSTART(IC)
1141           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
1142           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
1143           COG = 0             COG = 0  
1144             SGN = 0.
1145           mu  = 0           mu  = 0
1146           do i = istart,istop  c         print*,'-------'
1147             do i = INDMAX(IC),istart,-1
1148              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
1149              ivk  = nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
1150              is   = nst(MAXS(ic)+ipos)              if(CLSIGNAL(i).ge.cut)then              
1151  *            print*,'******************',istart,istop,ipos                 COG = COG + ipos*CLSIGNAL(i)
1152  *     $           ,MAXS(ic)+ipos,iv,ivk,is                 SGN = SGN + CLSIGNAL(i)
1153              cut  = incut*SIGMA(iv,ivk,is)                 mu = mu + 1
1154              if(SIGMA(iv,ivk,is).ne.CLSIGMA(i))  c               print*,ipos,CLSIGNAL(i)
1155       $           print*,'cog(0,ic) --> hai fatto qualche cazzata'              else
1156                   goto 10
1157                endif
1158             enddo
1159     10      continue
1160             do i = INDMAX(IC)+1,istop
1161                ipos = i-INDMAX(ic)
1162                cut  = incut*CLSIGMA(i)
1163              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
1164                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1165                   SGN = SGN + CLSIGNAL(i)
1166                 mu = mu + 1                 mu = mu + 1
1167  c               print*,ipos,CLSIGNAL(i),incut,cut  c               print*,ipos,CLSIGNAL(i)
1168                else
1169                   goto 20
1170              endif              endif
1171           enddo           enddo
1172           if(DEDX(ic).le.0)then   20      continue
1173              print*,'cog(0,ic) --> ic, dedx ',ic,DEDX(ic)           if(SGN.le.0)then
1174                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
1175              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1176              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
1177              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
1178           else           else
1179              COG=COG/DEDX(ic)              COG=COG/SGN
1180           endif           endif
1181  c         if(DEBUG)print*,'COG  (ic ',ic,' m',mu,')'  c         print*,'-------'
 c     $        ,cog  
1182                    
1183        else        else
1184                    
# Line 726  c     $        ,cog Line 1189  c     $        ,cog
1189    
1190        endif        endif
1191    
1192  c      print *,ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1193    
1194          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1195             if(DEBUG.eq.1)
1196         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1197             if(DEBUG.eq.1)
1198         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1199          endif
1200    
1201    
1202        return        return
1203        end        end
1204    
1205  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1206    
1207        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
1208  *-------------------------------------------------------  *-------------------------------------------------------
1209  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 747  c      print *,ncog,ic,cog,'//////////// Line 1219  c      print *,ncog,ic,cog,'////////////
1219        include 'calib.f'        include 'calib.f'
1220    
1221        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1222           f  = 4.           si = 8.4  !average good-strip noise
1223           si = 8.4           f  = 4.   !average bad-strip noise: f*si
1224           incut=incuty           incut=incuty
1225        else                      !X-view        else                      !X-view
1226           f  = 6.           si = 3.9  !average good-strip noise
1227           si = 3.9           f  = 6.   !average bad-strip noise: f*si
1228           incut=incutx           incut=incutx
1229        endif        endif
1230                
# Line 763  c      print *,ncog,ic,cog,'//////////// Line 1235  c      print *,ncog,ic,cog,'////////////
1235  *     --> signal of the central strip  *     --> signal of the central strip
1236           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1237           fsc = 1           fsc = 1
1238           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
1239             fsc = clsigma(INDMAX(ic))/si
1240  *     --> signal of adjacent strips  *     --> signal of adjacent strips
1241           sl1 = 0                !left 1           sl1 = 0                !left 1
1242           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 771  c      print *,ncog,ic,cog,'//////////// Line 1244  c      print *,ncog,ic,cog,'////////////
1244       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1245       $        )then       $        )then
1246              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
1247              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
1248  c         else              fsl1 = clsigma(INDMAX(ic)-1)/si
 c            fsl1 = 0  
1249           endif           endif
1250    
1251           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 782  c            fsl1 = 0 Line 1254  c            fsl1 = 0
1254       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1255       $        )then       $        )then
1256              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
1257              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
1258  c         else              fsl2 = clsigma(INDMAX(ic)-2)/si
 c            fsl2 = 0  
1259           endif           endif
1260           sr1 = 0                !right 1           sr1 = 0                !right 1
1261           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 794  c            fsl2 = 0 Line 1265  c            fsl2 = 0
1265       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1266       $        )then       $        )then
1267              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
1268              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
1269  c         else              fsr1 = clsigma(INDMAX(ic)+1)/si
 c            fsr1 = 0  
1270           endif               endif    
1271           sr2 = 0                !right 2           sr2 = 0                !right 2
1272           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 806  c            fsr1 = 0 Line 1276  c            fsr1 = 0
1276       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1277       $        )then       $        )then
1278              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
1279              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
1280  c         else              fsr2 = clsigma(INDMAX(ic)+2)/si
 c            fsr2 = 0  
1281           endif           endif
1282    
1283    
1284    
1285  ************************************************************  ************************************************************
1286  *     COG computation  *     COG2-3-4 computation
1287  ************************************************************  ************************************************************
1288    
1289  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
1290                    
1291           COG = 0.           vCOG = cog(ncog,ic)!0.
1292                    
1293           if(ncog.eq.2)then           if(ncog.eq.2)then
1294              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1295                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
1296                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1297                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1298              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
1299                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
1300                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1301                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1302              endif              endif
1303           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1304              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
1305              fbad_cog =              fbad_cog =
1306       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1307              fbad_cog =              fbad_cog =
1308       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1309           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1310              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1311                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1312                 fbad_cog =                 fbad_cog =
1313       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1314       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1315                 fbad_cog =                 fbad_cog =
1316       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1317       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
1318              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
1319                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1320                 fbad_cog =                 fbad_cog =
1321       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1322       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1323                 fbad_cog =                 fbad_cog =
1324       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1325       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1326              endif              endif
1327           else           else
1328              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1329              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1330              COG = 0.  c            COG = 0.
1331           endif           endif
1332                    
1333        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1334    *     =========================
1335    *     COG computation
1336    *     =========================
1337                    
1338           iv=VIEW(ic)           vCOG = cog(0,ic)
1339           istart=INDSTART(IC)  
1340           istop=TOTCLLENGTH           iv     = VIEW(ic)
1341           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1342           COG=0.           istop  = TOTCLLENGTH
1343           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1344           SDE=0.           SGN = 0.
1345           do i=istart,istop           SNU = 0.
1346              ipos=i-INDMAX(ic)           SDE = 0.
1347              il=nvk(MAXS(ic)+ipos)  
1348              is=nst(MAXS(ic)+ipos)           do i=INDMAX(IC),istart,-1
1349              cut=incut*SIGMA(iv,il,is)              ipos = i-INDMAX(ic)
1350                cut  = incut*CLSIGMA(i)
1351              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1352                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1353              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1354                   SDE = SDE + (ipos-vCOG)**2
1355                else
1356                   goto 10
1357                endif            
1358           enddo           enddo
1359           COG=COG/DEDX(ic)   10      continue
1360           do i=istart,istop           do i=INDMAX(IC)+1,istop
1361              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1362              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1363              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1364                 fs=1                 fs = clsigma(i)/si
1365                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1366                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1367                 SDE = SDE + (ipos-COG)**2              else
1368                   goto 20
1369              endif                          endif            
1370           enddo           enddo
1371           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1372             if(SDE.ne.0)then
1373                FBAD_COG=SNU/SDE
1374             else
1375                
1376             endif
1377    
1378        else        else
1379                                        
# Line 912  c      print*,sl2,sl1,sc,sr1,sr2 Line 1392  c      print*,sl2,sl1,sc,sr1,sr2
1392    
1393    
1394  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1395        real function fbad_cog0(ncog,ic)  
1396          real function riscogtheor(ncog,ic)
1397  *-------------------------------------------------------  *-------------------------------------------------------
 *     this function returns a factor that takes into  
 *     account deterioration of the spatial resolution  
 *     in the case BAD strips are included in the cluster.  
 *     This factor should multiply the nominal spatial  
 *     resolution.  
1398  *  *
1399  *     NB!!!  *     this function returns the expected resolution
1400  *     (this is the old version. It consider only the two  *     obtained by propagating the strip noise
1401  *     strips with the greatest signal. The new one is  *     to the center-of-gravity coordinate
1402  *     fbad_cog(ncog,ic) )  *
1403  *      *     ncog = n.strip used in the coordinate evaluation
1404    *     (ncog=0 => all strips above threshold)
1405    *
1406  *-------------------------------------------------------  *-------------------------------------------------------
1407    
1408        include 'commontracker.f'        include 'commontracker.f'
1409        include 'level1.f'        include 'level1.f'
1410        include 'calib.f'        include 'calib.f'
1411    
1412          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1413             incut = incuty
1414             pitch = pitchY / 1.e4
1415          else                      !X-view
1416             incut = incutx
1417             pitch = pitchX / 1.e4
1418          endif
1419          
1420          func = 100000.
1421          stot = 0.
1422    
1423          if (ncog.gt.0) then
1424    
1425  *     --> signal of the central strip  *     --> signal of the central strip
1426        sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1427             fsc = clsigma(INDMAX(ic))
1428    *     --> signal of adjacent strips
1429             sl1 = 0                !left 1
1430             fsl1 = 1               !left 1
1431             if(
1432         $        (INDMAX(ic)-1).ge.INDSTART(ic)
1433         $        )then
1434                sl1 = CLSIGNAL(INDMAX(ic)-1)
1435                fsl1 = clsigma(INDMAX(ic)-1)
1436             endif
1437    
1438  *     signal of adjacent strips           sl2 = 0                !left 2
1439  *     --> left           fsl2 = 1               !left 2
1440        sl1 = 0                  !left 1           if(
1441        if(       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1442       $     (INDMAX(ic)-1).ge.INDSTART(ic)       $        )then
1443       $     )              sl2 = CLSIGNAL(INDMAX(ic)-2)
1444       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))              fsl2 = clsigma(INDMAX(ic)-2)
1445             endif
1446        sl2 = 0                  !left 2           sr1 = 0                !right 1
1447        if(           fsr1 = 1               !right 1
1448       $     (INDMAX(ic)-2).ge.INDSTART(ic)           if(
1449       $     )       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1450       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))       $        .or.
1451         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1452  *     --> right       $        )then
1453        sr1 = 0                  !right 1              sr1 = CLSIGNAL(INDMAX(ic)+1)
1454        if(              fsr1 = clsigma(INDMAX(ic)+1)
1455       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))           endif    
1456       $     .or.           sr2 = 0                !right 2
1457       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)           fsr2 = 1               !right 2
1458       $     )           if(
1459       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1460         $        .or.
1461        sr2 = 0                  !right 2       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1462        if(       $        )then
1463       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))              sr2 = CLSIGNAL(INDMAX(ic)+2)
1464       $     .or.              fsr2 = clsigma(INDMAX(ic)+2)
1465       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)           endif
      $     )  
      $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
1466    
1467    
       if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  
          f  = 4.  
          si = 8.4  
       else                              !X-view  
          f  = 6.  
          si = 3.9  
       endif  
1468    
1469        fbad_cog = 1.  ************************************************************
1470        f0 = 1  *     COG2-3-4 computation
1471        f1 = 1  ************************************************************
1472        f2 = 1  
1473        f3 = 1    c      print*,sl2,sl1,sc,sr1,sr2
1474        if(sl1.gt.sr1.and.sl1.gt.0.)then          
1475             vCOG = cog(ncog,ic)!0.
1476                    
1477           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f           if(ncog.eq.1)then
1478           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f              func = 1./12.
1479  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f              stot = 1.
1480             elseif(ncog.eq.2)then
1481           if(ncog.eq.2.and.sl1.ne.0)then              if(sl1.gt.sr1)then
1482              fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)                 func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1483           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then                 stot = sl1+sc
1484              fbad_cog = 1.              elseif(sl1.le.sr1)then
1485           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then                 func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1486              fbad_cog = 1.                 stot = sc+sr1
1487                endif
1488             elseif(ncog.eq.3)then
1489                func =
1490         $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1491                stot = sl1+sc+sr1
1492             elseif(ncog.eq.4)then
1493                if(sl2.gt.sr2)then
1494                   func =
1495         $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1496         $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1497                   stot = sl2+sl1+sc+sr1
1498                elseif(sl2.le.sr2)then
1499                   func =
1500         $              (fsl1*(-1-vCOG)**2
1501         $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1502                   stot = sl2+sl1+sc+sr1
1503                endif
1504           else           else
1505              fbad_cog = 1.              print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG
1506         $            ,' not implemented'
1507           endif           endif
1508                    
1509        elseif(sl1.le.sr1.and.sr1.gt.0.)then        elseif(ncog.eq.0)then
1510    *     =========================
1511    *     COG computation
1512    *     =========================
1513            
1514             vCOG = cog(0,ic)
1515    
1516             iv     = VIEW(ic)
1517             istart = INDSTART(IC)
1518             istop  = TOTCLLENGTH
1519             if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1520    ccc         SGN = 0.
1521             SNU = 0.
1522    ccc         SDE = 0.
1523    
1524           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f           do i=INDMAX(IC),istart,-1
1525           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f              ipos = i-INDMAX(ic)
1526  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f              cut  = incut*CLSIGMA(i)
1527                if(CLSIGNAL(i).gt.cut)then
1528           if(ncog.eq.2.and.sr1.ne.0)then                 fs = clsigma(i)
1529              fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)                 SNU  = SNU + fs*(ipos-vCOG)**2
1530           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then                 stot = stot + CLSIGNAL(i)
1531              fbad_cog = 1.              else
1532           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then                 goto 10
1533              fbad_cog = 1.              endif            
1534             enddo
1535     10      continue
1536             do i=INDMAX(IC)+1,istop
1537                ipos = i-INDMAX(ic)
1538                cut  = incut*CLSIGMA(i)
1539                if(CLSIGNAL(i).gt.cut)then
1540                   fs = clsigma(i)
1541                   SNU  = SNU + fs*(ipos-vCOG)**2
1542                   stot = stot + CLSIGNAL(i)
1543                else
1544                   goto 20
1545                endif            
1546             enddo
1547     20      continue
1548             if(SDE.ne.0)then
1549                FUNC=SNU
1550           else           else
1551              fbad_cog = 1.              
1552           endif           endif
1553    
1554          else
1555                      
1556             FUNC=0
1557             print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1558             print*,'                               (NCOG must be >= 0)'
1559            
1560    
1561          endif
1562    
1563    
1564          if(stot.gt.0..and.func.gt.0.)then
1565             func = sqrt(func)
1566             func = pitch * func/stot
1567        endif        endif
1568    
1569        fbad_cog0 = sqrt(fbad_cog)        riscogtheor = func
1570    
1571        return        return
1572        end        end
1573    
1574    
1575    
1576    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1577    
1578          real function risetatheor(ncog,ic,angle)
1579    *-------------------------------------------------------
1580    *
1581    *     this function returns the expected resolution
1582    *     obtained by propagating the strip noise
1583    *     to the coordinate evaluated with non-linear eta-function
1584    *
1585    *     ncog = n.strip used in the coordinate evaluation
1586    *     (ncog=0 => ncog=2,3,4 according to angle)
1587    *
1588    *-------------------------------------------------------
1589    
1590          include 'commontracker.f'
1591          include 'level1.f'
1592          include 'calib.f'
1593    
1594    
1595          func = 1.
1596    
1597          iview   = VIEW(ic)            
1598          lad     = nld(MAXS(ic),VIEW(ic))
1599    
1600    *     ------------------------------------------------
1601    *     number of strip to be used (in case of ncog = 0)
1602    *     ------------------------------------------------
1603    
1604          inoeta = 0
1605    
1606          if(mod(int(iview),2).eq.1)then !Y-view
1607    
1608             pitch = pitchY / 1.e4
1609    
1610             if(ncog.eq.0)then
1611                if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
1612                   ncog = 2
1613                elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
1614                   ncog = 3
1615                elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
1616                   ncog = 4
1617                else
1618                   ncog   = 4
1619                   inoeta = 1
1620                endif            
1621             endif
1622    
1623          else                      !X-view
1624    
1625             pitch = pitchX / 1.e4
1626    
1627             if(ncog.eq.0)then
1628                if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
1629                   ncog = 2
1630                elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
1631                   ncog = 3
1632                elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
1633                   ncog = 4
1634                else
1635                   ncog = 4
1636                   inoeta = 1
1637                endif            
1638             endif
1639    
1640          endif
1641          
1642          func = riscogtheor(ncog,ic)
1643    
1644          risetatheor = func
1645          
1646          if(inoeta.eq.1)return ! no eta correction is applied --> exit
1647          if(ncog.lt.1.or.ncog.gt.4)return
1648    
1649    *     ----------------
1650    *     find angular bin
1651    *     ----------------
1652    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1653          do iang=1,nangbin
1654             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1655                iangle=iang
1656                goto 98
1657             endif
1658          enddo
1659          if(DEBUG.EQ.1)print*
1660         $     ,'risetatheor *** warning *** angle out of range: ',angle
1661          if(angle.le.angL(1))iang=1
1662          if(angle.ge.angR(nangbin))iang=nangbin
1663     98   continue                  !jump here if ok
1664    
1665    *     -------------
1666    *     within +/-0.5
1667    *     -------------
1668    
1669          vcog = cog(ncog,ic)          
1670    
1671          etamin = eta2(1,iang)
1672          etamax = eta2(netaval,iang)
1673    
1674          iaddmax=10
1675          iadd=0
1676     10   continue
1677          if(vcog.lt.etamin)then
1678             vcog = vcog + 1
1679             iadd = iadd + 1
1680             if(iadd>iaddmax)goto 111
1681             goto 10
1682          endif
1683     20   continue
1684          if(vcog.gt.etamax)then
1685             vcog = vcog - 1
1686             iadd = iadd - 1
1687             if(iadd<-1*iaddmax)goto 111
1688             goto 20
1689          endif
1690          goto 1111
1691     111  continue
1692          if(DEBUG.eq.1)
1693         $     print*,'risetatheor *** warning *** anomalous cluster'
1694          if(DEBUG.eq.1)
1695         $     print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)'
1696          vcog=0
1697     1111 continue
1698    
1699    *     ------------------------------------------------
1700    *     interpolation
1701    *     ------------------------------------------------
1702    
1703    
1704          ibin = netaval
1705          do i=2,netaval        
1706             if(ncog.eq.2)eta=eta2(i,iang)
1707             if(ncog.eq.3)eta=eta3(i,iang)
1708             if(ncog.eq.4)eta=eta4(i,iang)
1709             if(eta.ge.vcog)then
1710                ibin = i
1711                goto 99
1712             endif
1713          enddo
1714     99   continue
1715    
1716          if(ncog.eq.2)then
1717             x1 = eta2(ibin-1,iang)
1718             x2 = eta2(ibin,iang)
1719             y1 = feta2(ibin-1,iview,lad,iang)
1720             y2 = feta2(ibin,iview,lad,iang)
1721          elseif(ncog.eq.3)then
1722             x1 = eta3(ibin-1,iang)
1723             x2 = eta3(ibin,iang)
1724             y1 = feta3(ibin-1,iview,lad,iang)
1725             y2 = feta3(ibin,iview,lad,iang)
1726          elseif(ncog.eq.4)then
1727             x1 = eta4(ibin-1,iang)
1728             x2 = eta4(ibin,iang)
1729             y1 = feta4(ibin-1,iview,lad,iang)
1730             y2 = feta4(ibin,iview,lad,iang)
1731          endif
1732          
1733          func = func * (y2-y1)/(x2-x1)
1734    
1735          risetatheor = func
1736    
1737          return
1738          end
1739    
1740  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1741    
1742        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1743    
1744        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1745        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1096  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1811  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1811       +/       +/
1812    
1813        V(1)= abs(x)        V(1)= abs(x)
1814          if(V(1).gt.20.)V(1)=20.
1815    
1816        HQUADF = 0.        HQUADF = 0.
1817        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1110  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1826  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1826     20 CONTINUE     20 CONTINUE
1827        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1828                
1829        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1830    
1831        END        END
1832    
1833  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1834        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1835        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1836        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1837        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1186  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1902  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1902       +/       +/
1903    
1904        V(1) =  abs(x)        V(1) =  abs(x)
1905          if(V(1).gt.20.)V(1)=20.
1906    
1907        HQUADF = 0.        HQUADF = 0.
1908        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1909           HQDJ = 0.           HQDJ = 0.
# Line 1199  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1917  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1917     20 CONTINUE     20 CONTINUE
1918        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1919                
1920        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1921    
1922        END        END
1923  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1924        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1925        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1926        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1927        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1274  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1992  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1992       +/       +/
1993    
1994        V(1)=abs(x)        V(1)=abs(x)
1995          if(V(1).gt.20.)V(1)=20.
1996    
1997        HQUADF = 0.        HQUADF = 0.
1998        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1288  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2007  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2007     20 CONTINUE     20 CONTINUE
2008        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2009                
2010        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
2011    
2012        END        END
2013  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2014        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
2015        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
2016        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
2017        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1345  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2064  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2064       +/       +/
2065    
2066        v(1)= abs(x)        v(1)= abs(x)
2067          if(V(1).gt.20.)V(1)=20.
2068                
2069        HQUADF = 0.        HQUADF = 0.
2070        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1359  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2079  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2079     20 CONTINUE     20 CONTINUE
2080        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2081    
2082        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
2083    
2084        END        END
2085  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1411  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2131  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2131       +/       +/
2132                
2133        V(1)=abs(x)        V(1)=abs(x)
2134          if(V(1).gt.20.)V(1)=20.
2135    
2136        HQUADF = 0.        HQUADF = 0.
2137        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1491  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2212  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2212       +/       +/
2213    
2214        V(1)=abs(x)        V(1)=abs(x)
2215          if(V(1).gt.20.)V(1)=20.
2216    
2217        HQUADF = 0.        HQUADF = 0.
2218        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1508  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2230  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2230        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
2231    
2232        END        END
2233    
2234    
2235  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2236          real function pfacorr(ic,angle)
2237    *--------------------------------------------------------------
2238    *     this function returns the landi correction for this cluster
2239    *--------------------------------------------------------------
2240          include 'commontracker.f'
2241          include 'calib.f'
2242          include 'level1.f'
2243    
2244          real angle
2245          integer iview,lad
2246    
2247          iview = VIEW(ic)            
2248          lad = nld(MAXS(ic),VIEW(ic))
2249    
2250    *     find angular bin
2251    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
2252          do iang=1,nangbin
2253             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
2254                iangle=iang
2255                goto 98
2256             endif
2257          enddo
2258          if(DEBUG.eq.1)
2259         $     print*,'pfacorr *** warning *** angle out of range: ',angle
2260          if(angle.le.angL(1))iang=1
2261          if(angle.ge.angR(nangbin))iang=nangbin
2262     98   continue                  !jump here if ok
2263    
2264          pfacorr = fcorr(iview,lad,iang)
2265    
2266          if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr
2267    
2268          
2269    c 100  return
2270          return
2271          end
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.23