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