/[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.8 by pam-fi, Fri Feb 16 14:56:02 2007 UTC revision 1.24 by pam-fi, Sat Mar 22 08:32:50 2008 UTC
# Line 1  Line 1 
1  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2  *     this file contains all subroutines and functions  *     this file contains all subroutines and functions
3  *     that are needed for position finding algorithms  *     that are needed for position finding algorithms:
4  *      *          
5    *     subroutine idtoc(ipfa,cpfa)
6    *
7    *     subroutine applypfa(PFAtt,ic,ang,corr,res)
8    *
9    *     integer function npfastrips(ic,angle)
10    *
11    *     -----------------------------------------------------------------
12    *     p.f.a.
13    *     -----------------------------------------------------------------
14    *     real function pfaeta(ic,angle)
15    *     real function pfaetal(ic,angle)
16    *     real function pfaeta2(ic,angle)
17    *     real function pfaeta3(ic,angle)
18    *     real function pfaeta4(ic,angle)
19    *     real function cog(ncog,ic)
20    *
21    *     -----------------------------------------------------------------
22    *     risoluzione spaziale media, stimata dalla simulazione (samuele)
23    *     -----------------------------------------------------------------
24    *     FUNCTION risxeta2(angle)
25    *     FUNCTION risxeta3(angle)
26    *     FUNCTION risxeta4(angle)
27    *     FUNCTION risyeta2(angle)
28    *     FUNCTION risy_cog(angle)
29    *     FUNCTION risx_cog(angle)
30    *     real function riseta(iview,angle)
31    *     -----------------------------------------------------------------
32    *     fattore moltiplicativo per tenere conto della dipendenza della
33    *     risoluzione dal rumore delle strip
34    *     -----------------------------------------------------------------
35    *     real function fbad_cog(ncog,ic)
36    *     real function fbad_eta(ic,angle)
37    *
38    *     -----------------------------------------------------------------
39    *     NUOVO APPROCCIO PER LA STIMA DELLA RISOLUZIONE
40    *     -----------------------------------------------------------------
41    *     real function riscogtheor(ncog,ic)
42    *     real function risetatheor(ncog,ic,angle)
43    *
44    *     -----------------------------------------------------------------
45    *     correzione landi
46    *     -----------------------------------------------------------------
47    *     real function pfacorr(ic,angle)
48  *  *
49    *     real function effectiveangle(ang,iview,bbb)
50    *     real function fieldcorr(iview,bbb)
51    *
52    *     NB - The angle is the "effective angle", which is relative
53    *          to the sensor and it takes into account the magnetic field
54    *
55    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
56    
57          subroutine idtoc(ipfa,cpfa)
58          
59          integer ipfa
60          character*10 cpfa
61    
62          CPFA='COG4'
63          if(ipfa.eq.0)CPFA='ETA'
64          if(ipfa.eq.2)CPFA='ETA2'
65          if(ipfa.eq.3)CPFA='ETA3'
66          if(ipfa.eq.4)CPFA='ETA4'
67          if(ipfa.eq.5)CPFA='ETAL'
68          if(ipfa.eq.10)CPFA='COG'
69          if(ipfa.eq.11)CPFA='COG1'
70          if(ipfa.eq.12)CPFA='COG2'
71          if(ipfa.eq.13)CPFA='COG3'
72          if(ipfa.eq.14)CPFA='COG4'
73          
74          end
75    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
76          real function effectiveangle(ang,iview,bbb)
77    
78          include 'commontracker.f'
79    
80          effectiveangle = 0.
81    
82          if(mod(iview,2).eq.0)then
83    c     =================================================
84    c     X view
85    c     =================================================
86    c     here bbb is the y component of the m.field
87             angx = ang
88             by   = bbb
89             if(iview.eq.12) angx = -1. * ang
90             if(iview.eq.12) by   = -1. * bbb
91    cc         tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001 !ORRORE!!
92             tgtemp  = tan(angx*acos(-1.)/180.) + pmuH_h*by*0.00001
93    
94          elseif(mod(iview,2).eq.1)then
95    c     =================================================
96    c     Y view
97    c     =================================================        
98    c     here bbb is the x component of the m.filed
99             angy = ang
100             bx   = bbb
101             tgtemp  = tan(angy*acos(-1.)/180.)+pmuH_e*bx*0.00001        
102    
103          endif      
104          effectiveangle = 180.*atan(tgtemp)/acos(-1.)
105    
106          return
107          end
108    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
109          real function fieldcorr(iview,bbb)
110    
111          include 'commontracker.f'
112    
113          fieldcorr = 0.
114    
115          if(mod(iview,2).eq.0)then
116    
117    c     =================================================
118    c     X view
119    c     =================================================
120    c     here bbb is the y component of the m.field
121             by   = bbb
122             if(iview.eq.12) by = -1. * bbb
123             fieldcorr     = -1. * 0.5*pmuH_h*by*0.00001*SiDimZ/pitchX
124    
125          elseif(mod(iview,2).eq.1)then
126    c     =================================================
127    c     Y view
128    c     =================================================        
129    c     here bbb is the x component of the m.filed
130             bx   = bbb
131             fieldcorr     = 0.5*pmuH_e*bx*0.00001*SiDimZ/pitchY
132    
133          endif      
134          
135          return
136          end
137  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
138    
139          subroutine applypfa(PFAtt,ic,ang,corr,res)
140    *---------------------------------------------------------------
141    *     this subroutine calculate the coordinate of cluster ic (in
142    *     strip units), relative to the strip with the maximum signal,
143    *     and its spatial resolution (in cm), applying PFAtt.
144    *     ang is the effective angle, relative to the sensor
145    *---------------------------------------------------------------
146    
147          character*4 PFAtt
148          include 'commontracker.f'
149          include 'level1.f'
150    
151          corr = 0
152          res  = 0
153    
154          if(ic.le.0)return
155    
156          iview   = VIEW(ic)
157    
158          if(mod(iview,2).eq.0)then
159    c     =================================================
160    c     X view
161    c     =================================================
162    
163             res = RESXAV
164    
165             if(PFAtt.eq.'COG1')then
166    
167                corr   = 0
168                res = 1e-4*pitchX/sqrt(12.)!!res
169    
170             elseif(PFAtt.eq.'COG2')then
171    
172                corr    = cog(2,ic)            
173                res = risx_cog(abs(ang))!TEMPORANEO              
174                res = res*fbad_cog(2,ic)
175    
176             elseif(PFAtt.eq.'COG3')then
177    
178                corr    = cog(3,ic)            
179                res = risx_cog(abs(ang))!TEMPORANEO                      
180                res = res*fbad_cog(3,ic)
181    
182             elseif(PFAtt.eq.'COG4')then
183    
184                corr    = cog(4,ic)            
185                res = risx_cog(abs(ang))!TEMPORANEO                      
186                res = res*fbad_cog(4,ic)
187    
188             elseif(PFAtt.eq.'ETA2')then
189    
190                corr  = pfaeta2(ic,ang)          
191                res = risxeta2(abs(ang))
192                res = res*fbad_cog(2,ic)
193    
194             elseif(PFAtt.eq.'ETA3')then                        
195    
196                corr  = pfaeta3(ic,ang)          
197                res = risxeta3(abs(ang))                      
198                res = res*fbad_cog(3,ic)              
199    
200             elseif(PFAtt.eq.'ETA4')then                        
201    
202                corr  = pfaeta4(ic,ang)            
203                res = risxeta4(abs(ang))                      
204                res = res*fbad_cog(4,ic)              
205    
206             elseif(PFAtt.eq.'ETA')then  
207    
208                corr  = pfaeta(ic,ang)            
209    c            res = riseta(ic,ang)                    
210                res = riseta(iview,ang)                    
211                res = res*fbad_eta(ic,ang)            
212    
213             elseif(PFAtt.eq.'ETAL')then  
214    
215                corr  = pfaetal(ic,ang)            
216                res = riseta(iview,ang)                    
217                res = res*fbad_eta(ic,ang)            
218    
219             elseif(PFAtt.eq.'COG')then          
220    
221                corr  = cog(0,ic)            
222                res = risx_cog(abs(ang))                    
223                res = res*fbad_cog(0,ic)
224    
225             else
226                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAtt
227             endif
228    
229    
230    *     ======================================
231    *     temporary patch for saturated clusters
232    *     ======================================
233             if( nsatstrips(ic).gt.0 )then
234                corr  = cog(4,ic)            
235                res = pitchX*1e-4/sqrt(12.)
236    cc            cc=cog(4,ic)
237    c$$$            print*,ic,' *** ',cc
238    c$$$            print*,ic,' *** ',res
239             endif
240    
241    
242          elseif(mod(iview,2).eq.1)then
243    c     =================================================
244    c     Y view
245    c     =================================================
246    
247             res = RESYAV
248    
249             if(PFAtt.eq.'COG1')then
250    
251                corr  = 0  
252                res = 1e-4*pitchY/sqrt(12.)!res  
253    
254             elseif(PFAtt.eq.'COG2')then
255    
256                corr    = cog(2,ic)
257                res = risy_cog(abs(ang))!TEMPORANEO
258                res = res*fbad_cog(2,ic)
259    
260             elseif(PFAtt.eq.'COG3')then
261    
262                corr    = cog(3,ic)
263                res = risy_cog(abs(ang))!TEMPORANEO
264                res = res*fbad_cog(3,ic)
265    
266             elseif(PFAtt.eq.'COG4')then
267    
268                corr    = cog(4,ic)
269                res = risy_cog(abs(ang))!TEMPORANEO
270                res = res*fbad_cog(4,ic)
271    
272             elseif(PFAtt.eq.'ETA2')then
273    
274                corr    = pfaeta2(ic,ang)          
275                res = risyeta2(abs(ang))              
276                res = res*fbad_cog(2,ic)
277    
278             elseif(PFAtt.eq.'ETA3')then                      
279    
280                corr    = pfaeta3(ic,ang)
281                res = res*fbad_cog(3,ic)  
282    
283             elseif(PFAtt.eq.'ETA4')then  
284    
285                corr    = pfaeta4(ic,ang)
286                res = res*fbad_cog(4,ic)
287    
288             elseif(PFAtt.eq.'ETA')then
289    
290                corr    = pfaeta(ic,ang)
291    c            res = riseta(ic,ang)  
292                res = riseta(iview,ang)  
293                res = res*fbad_eta(ic,ang)
294    
295             elseif(PFAtt.eq.'ETAL')then
296    
297                corr    = pfaetal(ic,ang)
298                res = riseta(iview,ang)  
299                res = res*fbad_eta(ic,ang)
300    
301             elseif(PFAtt.eq.'COG')then
302    
303                corr    = cog(0,ic)            
304                res = risy_cog(abs(ang))
305                res = res*fbad_cog(0,ic)
306    
307             else
308                if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (y) --> ',PFAtt
309             endif
310    
311    
312    *     ======================================
313    *     temporary patch for saturated clusters
314    *     ======================================
315             if( nsatstrips(ic).gt.0 )then
316                corr    = cog(4,ic)            
317                res = pitchY*1e-4/sqrt(12.)
318    cc            cc=cog(4,ic)
319    c$$$            print*,ic,' *** ',cc
320    c$$$            print*,ic,' *** ',res
321             endif
322            
323          endif
324          end
325    
326    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
327          integer function npfastrips(ic,angle)
328    *--------------------------------------------------------------
329    *     thid function returns the number of strips used
330    *     to evaluate the position of a cluster, according to the p.f.a.
331    *--------------------------------------------------------------
332          include 'commontracker.f'
333          include 'level1.f'
334          include 'calib.f'
335    
336          character*4 usedPFA
337          
338    
339    
340          call idtoc(pfaid,usedPFA)
341    
342          npfastrips=-1
343    
344          if(usedPFA.eq.'COG1')npfastrips=1
345          if(usedPFA.eq.'COG2')npfastrips=2
346          if(usedPFA.eq.'COG3')npfastrips=3
347          if(usedPFA.eq.'COG4')npfastrips=4
348          if(usedPFA.eq.'ETA2')npfastrips=2
349          if(usedPFA.eq.'ETA3')npfastrips=3
350          if(usedPFA.eq.'ETA4')npfastrips=4
351    *     ----------------------------------------------------------------
352          if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
353    c         print*,VIEW(ic),angle
354             if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
355                if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
356                   npfastrips=2
357                elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
358                   npfastrips=3
359                elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
360                   npfastrips=4
361                else
362                   npfastrips=4     !COG4
363                endif                        
364             else                   !X-view
365                if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
366                   npfastrips=2
367                elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
368                   npfastrips=3
369                elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
370                   npfastrips=4
371                else
372                   npfastrips=4     !COG4
373                endif                        
374             endif
375          endif
376    *     ----------------------------------------------------------------
377          if(usedPFA.eq.'COG')then
378    
379             npfastrips=0
380    
381    c$$$         iv=VIEW(ic)
382    c$$$         if(mod(iv,2).eq.1)incut=incuty
383    c$$$         if(mod(iv,2).eq.0)incut=incutx
384    c$$$         istart = INDSTART(IC)
385    c$$$         istop  = TOTCLLENGTH
386    c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
387    c$$$         mu  = 0
388    c$$$         do i = INDMAX(IC),istart,-1
389    c$$$            ipos = i-INDMAX(ic)
390    c$$$            cut  = incut*CLSIGMA(i)
391    c$$$            if(CLSIGNAL(i).ge.cut)then
392    c$$$               mu = mu + 1
393    c$$$               print*,i,mu
394    c$$$            else
395    c$$$               goto 10
396    c$$$            endif
397    c$$$         enddo
398    c$$$ 10      continue
399    c$$$         do i = INDMAX(IC)+1,istop
400    c$$$            ipos = i-INDMAX(ic)
401    c$$$            cut  = incut*CLSIGMA(i)
402    c$$$            if(CLSIGNAL(i).ge.cut)then
403    c$$$               mu = mu + 1
404    c$$$               print*,i,mu
405    c$$$            else
406    c$$$               goto 20
407    c$$$            endif
408    c$$$         enddo
409    c$$$ 20      continue
410    c$$$         npfastrips=mu
411    
412          endif
413    *     ----------------------------------------------------------------
414    
415    c      print*,pfaid,usedPFA,angle,npfastrips
416    
417          return
418          end
419    
420  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
421        real function pfaeta(ic,angle)        real function pfaeta(ic,angle)
# Line 17  Line 428 
428  *     according to the angle  *     according to the angle
429  *--------------------------------------------------------------  *--------------------------------------------------------------
430        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
431        include 'level1.f'        include 'level1.f'
432          include 'calib.f'
433                
434        pfaeta = 0        pfaeta = 0
435    
436        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
437                
438           pfaeta = pfaeta2(ic,angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
439                  pfaeta = pfaeta2(ic,angle)
440    cc            print*,pfaeta2(ic,angle)
441             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
442                pfaeta = pfaeta3(ic,angle)
443             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
444                pfaeta = pfaeta4(ic,angle)
445             else
446                pfaeta = cog(4,ic)
447             endif            
448    
449        else                      !X-view        else                      !X-view
450    
451           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
452              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
453           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
454              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
455           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
456              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
457           endif           else
458                pfaeta = cog(4,ic)
459             endif            
460                            
461        endif        endif
462                
 c      print*,'pfaeta ',pfaeta, angle  
   
463   100  return   100  return
464        end        end
465    
466  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
467        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
468  *--------------------------------------------------------------  *--------------------------------------------------------------
469  *     this function returns the average spatial resolution  *     this function returns the position (in strip units)
 *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  
470  *     it calls:  *     it calls:
471  *     - risx_eta2(angle)  *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
472  *     - risy_eta2(angle)  *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
473  *     - risx_eta3(angle)  *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
 *     - risx_eta4(angle)  
474  *     according to the angle  *     according to the angle
475  *--------------------------------------------------------------  *--------------------------------------------------------------
476        include 'commontracker.f'        include 'commontracker.f'
 c      include 'calib.f'  
477        include 'level1.f'        include 'level1.f'
478          include 'calib.f'
 c$$$      logical DEBUG  
 c$$$      common/dbg/DEBUG  
479                
480        ris_eta = 0        pfaetal = 0
481    
482        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
483                
484           ris_eta = risy_eta2(angle)           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
485           if(abs(angle).gt.21.)ris_eta = risy_eta2(21.)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
486    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
487             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
488                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
489             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
490                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
491             else
492                pfaetal = cog(4,ic)
493             endif            
494    
495        else                      !X-view        else                      !X-view
496    
497           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
498              ris_eta = risx_eta2(angle)              pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
499           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then  cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
500              ris_eta = risx_eta3(angle)           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
501           elseif(abs(angle).gt.15..and.abs(angle).le.21.)then              pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
502              ris_eta = risx_eta4(angle)           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
503           elseif(abs(angle).gt.21.)then              pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
504              ris_eta = risx_eta4(21.)           else
505           endif              pfaetal = cog(4,ic)
506             endif            
507                            
508        endif        endif
509          
510     100  return
511          end
512    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
513    c      real function riseta(ic,angle)
514          real function riseta(iview,angle)
515    *--------------------------------------------------------------
516    *     this function returns the average spatial resolution
517    *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
518    *     it calls:
519    *     - risxeta2(angle)
520    *     - risyeta2(angle)
521    *     - risxeta3(angle)
522    *     - risxeta4(angle)
523    *     according to the angle
524    *--------------------------------------------------------------
525          include 'commontracker.f'
526          include 'level1.f'
527          include 'calib.f'
528    
529  c$$$      if(DEBUG)print*,'ris   (ic ',ic,' ang',angle,')'        riseta = 0
530  c$$$     $     ,' -->',ris_eta  
531    c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
532          if(mod(iview,2).eq.1)then !Y-view
533          
534    
535             if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
536                riseta = risyeta2(angle)
537             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
538                riseta = risy_cog(angle) !ATTENZIONE!!
539             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
540                riseta = risy_cog(angle) !ATTENZIONE!!
541             else
542                riseta = risy_cog(angle)
543             endif            
544    
545          else                      !X-view
546    
547             if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
548                riseta = risxeta2(angle)
549             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
550                riseta = risxeta3(angle)
551             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
552                riseta = risxeta4(angle)
553             else
554                riseta = risx_cog(angle)
555             endif            
556                
557          endif
558    
559    
560   100  return   100  return
# Line 100  c$$$     $     ,' -->',ris_eta Line 570  c$$$     $     ,' -->',ris_eta
570  *     resolution.  *     resolution.
571  *     It calls the function FBAD_COG(NCOG,IC),  *     It calls the function FBAD_COG(NCOG,IC),
572  *     accordingto the angle  *     accordingto the angle
573    *
574    *     >>> cosi` non e` corretto!!
575    *     >>> l'errore sulla coordinata eta si ottiene moltiplicando
576    *     >>> l'errore sulla coordinata cog per la derivata della
577    *     >>> distribuzione eta... pur sapendolo l'ho sempre ignorato...
578    *     >>> deve essere modificato!!!!
579    *
580  *-------------------------------------------------------  *-------------------------------------------------------
581    
582        include 'commontracker.f'        include 'commontracker.f'
583        include 'level1.f'        include 'level1.f'
584  *      include 'calib.f'        include 'calib.f'
585        fbad_eta = 0        fbad_eta = 0
586    
587        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
588                
589           fbad_eta = fbad_cog(2,ic)           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
590                fbad_eta = fbad_cog(2,ic)
591             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
592                fbad_eta = fbad_cog(3,ic)
593             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
594                fbad_eta = fbad_cog(4,ic)
595             else
596                fbad_eta = fbad_cog(4,ic)
597             endif            
598    
599        else                      !X-view        else                      !X-view
600    
601           if(abs(angle).le.10.)then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
602              fbad_eta = fbad_cog(2,ic)              fbad_eta = fbad_cog(2,ic)
603           elseif(abs(angle).gt.10..and.abs(angle).le.15.)then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
604              fbad_eta = fbad_cog(3,ic)              fbad_eta = fbad_cog(3,ic)
605           elseif(abs(angle).gt.15.)then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
606              fbad_eta = fbad_cog(4,ic)              fbad_eta = fbad_cog(4,ic)
607           endif           else
608                fbad_eta = fbad_cog(4,ic)
609             endif            
610                            
611        endif        endif
612    
# Line 127  c$$$     $     ,' -->',ris_eta Line 614  c$$$     $     ,' -->',ris_eta
614        end        end
615    
616  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
617  c*****************************************************        real function pfaeta2(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta2(cog2,view,lad,angle)  
       real function pfaeta2(ic,angle) !(1)  
618  *--------------------------------------------------------------  *--------------------------------------------------------------
619  *     this function returns  *     this function returns
620  *  *
# Line 150  c      real function pfaeta2(cog2,view,l Line 633  c      real function pfaeta2(cog2,view,l
633        real cog2,angle        real cog2,angle
634        integer iview,lad        integer iview,lad
635    
636  c      logical DEBUG        iview   = VIEW(ic)            
637  c      common/dbg/DEBUG        lad     = nld(MAXS(ic),VIEW(ic))
638          cog2    = cog(2,ic)          
639  c      print*,'## pfaeta2 ',ic,angle        pfaeta2 = cog2
       iview = VIEW(ic)              !(1)  
       lad = nld(MAXS(ic),VIEW(ic)) !(1)  
       cog2 = cog(2,ic)             !(1)  
       pfaeta2=cog2  
640    
641    *     ----------------
642  *     find angular bin  *     find angular bin
643    *     ----------------
644  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
645        do iang=1,nangbin        do iang=1,nangbin
 c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  
646           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
647              iangle=iang              iangle=iang
648              goto 98              goto 98
649           endif           endif
650        enddo        enddo
651        if(DEBUG)        if(DEBUG.EQ.1)
652       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
653        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
654        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
655   98   continue                  !jump here if ok   98   continue                  !jump here if ok
656    
657    *     -------------
658    *     within +/-0.5
659    *     -------------
660    
661  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  
 c$$$c         print*,'pfaeta2 *** warning *** argument out of range: ',cog2  
 c$$$*         goto 100  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
662        iadd=0        iadd=0
663   10   continue   10   continue
664        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
665           cog2 = cog2 + 1           cog2 = cog2 + 1
666           iadd = iadd + 1           iadd = iadd + 1
667             if(iadd>iaddmax)goto 111
668           goto 10           goto 10
669        endif        endif
670   20   continue   20   continue
671        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
672           cog2 = cog2 - 1           cog2 = cog2 - 1
673           iadd = iadd - 1           iadd = iadd - 1
674             if(iadd<-1*iaddmax)goto 111
675           goto 20           goto 20
676        endif        endif
677          goto 1111
678     111  continue
679          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
680          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
681          cog2=0
682     1111 continue
683    
684  *     --------------------------------  *     --------------------------------
685  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 244  c$$$         pfaeta2=pfaeta2+1.   !temp Line 714  c$$$         pfaeta2=pfaeta2+1.   !temp
714  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
715  c$$$      endif  c$$$      endif
716    
717        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
718       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
719    
720    
# Line 252  c$$$      endif Line 722  c$$$      endif
722        end        end
723    
724  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
 c*****************************************************  
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta3(cog3,view,lad,angle)  
725        real function pfaeta3(ic,angle) !(1)        real function pfaeta3(ic,angle) !(1)
726  *--------------------------------------------------------------  *--------------------------------------------------------------
727  *     this function returns  *     this function returns
# Line 275  c      real function pfaeta3(cog3,view,l Line 741  c      real function pfaeta3(cog3,view,l
741        real cog3,angle        real cog3,angle
742        integer iview,lad        integer iview,lad
743    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
744    
745  c      print*,'## pfaeta3 ',ic,angle        iview = VIEW(ic)            
746          lad = nld(MAXS(ic),VIEW(ic))
747        iview = VIEW(ic)              !(1)        cog3 = cog(3,ic)
748        lad = nld(MAXS(ic),VIEW(ic)) !(1)        cc = cog3
749        cog3 = cog(3,ic)             !(1)        cog3 = cc
750        pfaeta3=cog3        pfaeta3=cog3
751    
752    *     ----------------
753  *     find angular bin  *     find angular bin
754    *     ----------------
755  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
756        do iang=1,nangbin        do iang=1,nangbin
757  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 294  c         print*,'~~~~~~~~~~~~ ',iang,an Line 760  c         print*,'~~~~~~~~~~~~ ',iang,an
760              goto 98              goto 98
761           endif           endif
762        enddo        enddo
763        if(DEBUG)        if(DEBUG.EQ.1)
764       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
765        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
766        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
767   98   continue                  !jump here if ok   98   continue                  !jump here if ok
768    
769    *     -------------
770    *     within +/-0.5
771    *     -------------
772    
773  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  
   
         
774        iadd=0        iadd=0
775   10   continue   10   continue
776        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
777           cog3 = cog3 + 1           cog3   = cog3 + 1.
778           iadd = iadd + 1           iadd = iadd + 1
779             if(iadd>iaddmax) goto 111
780           goto 10           goto 10
781        endif        endif
782   20   continue   20   continue
783        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
784           cog3 = cog3 - 1           cog3 = cog3 - 1.
785           iadd = iadd - 1           iadd = iadd - 1
786             if(iadd<-1*iaddmax) goto 111
787           goto 20           goto 20
788        endif        endif
789          goto 1111
790     111  continue
791          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
792          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
793          cog3=0      
794     1111 continue
795    
796  *     --------------------------------  *     --------------------------------
797  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 360  c            print*,'-----',x1,x2,y1,y2 Line 817  c            print*,'-----',x1,x2,y1,y2
817        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
818        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
819    
 c$$$      if(iflag.eq.1)then  
 c$$$         pfaeta2=pfaeta2-1.   !temp  
 c$$$         cog2=cog2-1.           !temp  
 c$$$      endif  
 c$$$      if(iflag.eq.-1)then  
 c$$$         pfaeta2=pfaeta2+1.   !temp  
 c$$$         cog2=cog2+1.           !temp  
 c$$$      endif  
820    
821        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
822       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
823    
824   100  return   100  return
825        end        end
826    
827  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
828  c*****************************************************        real function pfaeta4(ic,angle)
 cccccc 02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
 c      real function pfaeta4(cog4,view,lad,angle)  
       real function pfaeta4(ic,angle) !(1)  
829  *--------------------------------------------------------------  *--------------------------------------------------------------
830  *     this function returns  *     this function returns
831  *  *
# Line 399  c      real function pfaeta4(cog4,view,l Line 844  c      real function pfaeta4(cog4,view,l
844        real cog4,angle        real cog4,angle
845        integer iview,lad        integer iview,lad
846    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
 c      print*,'## pfaeta4 ',ic,angle  
847    
848        iview = VIEW(ic)             !(1)        iview = VIEW(ic)            
849        lad = nld(MAXS(ic),VIEW(ic)) !(1)        lad = nld(MAXS(ic),VIEW(ic))
850        cog4=cog(4,ic)               !(1)        cog4=cog(4,ic)
851        pfaeta4=cog4        pfaeta4=cog4
852    
853    *     ----------------
854  *     find angular bin  *     find angular bin
855    *     ----------------
856  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
857        do iang=1,nangbin        do iang=1,nangbin
858  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 418  c         print*,'~~~~~~~~~~~~ ',iang,an Line 861  c         print*,'~~~~~~~~~~~~ ',iang,an
861              goto 98              goto 98
862           endif           endif
863        enddo        enddo
864        if(DEBUG)        if(DEBUG.EQ.1)
865       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
866        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
867        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
868   98   continue                  !jump here if ok   98   continue                  !jump here if ok
869    
870    *     -------------
871    *     within +/-0.5
872    *     -------------
873    
874  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  
   
         
875        iadd=0        iadd=0
876   10   continue   10   continue
877        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
878           cog4 = cog4 + 1           cog4 = cog4 + 1
879           iadd = iadd + 1           iadd = iadd + 1
880             if(iadd>iaddmax)goto 111
881           goto 10           goto 10
882        endif        endif
883   20   continue   20   continue
884        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
885           cog4 = cog4 - 1           cog4 = cog4 - 1
886           iadd = iadd - 1           iadd = iadd - 1
887             if(iadd<-1*iaddmax)goto 111
888           goto 20           goto 20
889        endif        endif
890          goto 1111
891     111  continue
892          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
893          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
894          cog4=0
895     1111 continue
896    
897  *     --------------------------------  *     --------------------------------
898  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 493  c$$$         pfaeta2=pfaeta2+1.   !temp Line 927  c$$$         pfaeta2=pfaeta2+1.   !temp
927  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
928  c$$$      endif  c$$$      endif
929    
930        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
931       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
932    
933   100  return   100  return
# Line 502  c$$$      endif Line 936  c$$$      endif
936    
937    
938  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
       real function cog0(ncog,ic)  
 *-------------------------------------------------  
 *     this function returns  
 *  
 *     - the Center-Of-Gravity of the cluster IC  
 *     evaluated using NCOG strips,  
 *     calculated relative to MAXS(IC)  
 *      
 *     - zero in case that not  enough strips  
 *     have a positive signal  
 *      
 *     NOTE:  
 *     This is the old definition, used by Straulino.  
 *     The new routine, according to Landi,  
 *     is COG(NCOG,IC)  
 *-------------------------------------------------  
   
   
       include 'commontracker.f'  
       include 'level1.f'  
         
 *     --> signal of the central strip  
       sc = CLSIGNAL(INDMAX(ic)) !center  
   
 *     signal of adjacent strips  
 *     --> left  
       sl1 = 0                  !left 1  
       if(  
      $     (INDMAX(ic)-1).ge.INDSTART(ic)  
      $     )  
      $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
   
       sl2 = 0                  !left 2  
       if(  
      $     (INDMAX(ic)-2).ge.INDSTART(ic)  
      $     )  
      $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
   
 *     --> right  
       sr1 = 0                  !right 1  
       if(  
      $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
      $     .or.  
      $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
      $     )  
      $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
   
       sr2 = 0                  !right 2  
       if(  
      $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
      $     .or.  
      $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
      $     )  
      $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
         
 ************************************************************  
 *     COG computation  
 ************************************************************  
   
 c      print*,sl2,sl1,sc,sr1,sr2  
   
       COG = 0.  
           
       if(sl1.gt.sr1.and.sl1.gt.0.)then  
           
          if(ncog.eq.2.and.sl1.ne.0)then  
             COG = -sl1/(sl1+sc)          
          elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
             COG = (sr1-sl1)/(sl1+sc+sr1)  
          elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
             COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  
          else  
             COG = 0.  
          endif  
           
       elseif(sl1.le.sr1.and.sr1.gt.0.)then  
   
          if(ncog.eq.2.and.sr1.ne.0)then  
             COG = sr1/(sc+sr1)              
          elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
             COG = (sr1-sl1)/(sl1+sc+sr1)  
          elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
             COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  
          else  
             COG = 0.  
          endif  
   
       endif  
   
       COG0 = COG  
   
 c      print *,ncog,ic,cog,'/////////////'  
   
       return  
       end  
   
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
939        real function cog(ncog,ic)        real function cog(ncog,ic)
940  *-------------------------------------------------  *-------------------------------------------------
941  *     this function returns  *     this function returns
# Line 618  c      print *,ncog,ic,cog,'//////////// Line 955  c      print *,ncog,ic,cog,'////////////
955        include 'calib.f'        include 'calib.f'
956        include 'level1.f'        include 'level1.f'
957                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
958    
959    
960        if (ncog.gt.0) then        if (ncog.gt.0) then
# Line 630  c      common/dbg/DEBUG Line 965  c      common/dbg/DEBUG
965  *     --> signal of the central strip  *     --> signal of the central strip
966           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
967  *     signal of adjacent strips  *     signal of adjacent strips
968           sl1 = 0                !left 1           sl1 = -9999.           !left 1
969           if(           if(
970       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
971       $        )       $        )
972       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
973                    
974           sl2 = 0                !left 2           sl2 = -9999.           !left 2
975           if(           if(
976       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
977       $        )       $        )
978       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
979                    
980           sr1 = 0                !right 1           sr1 = -9999.           !right 1
981           if(           if(
982       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
983       $        .or.       $        .or.
# Line 650  c      common/dbg/DEBUG Line 985  c      common/dbg/DEBUG
985       $        )       $        )
986       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
987                    
988           sr2 = 0                !right 2           sr2 = -9999.           !right 2
989           if(           if(
990       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
991       $        .or.       $        .or.
# Line 660  c      common/dbg/DEBUG Line 995  c      common/dbg/DEBUG
995                    
996           COG = 0.           COG = 0.
997                    
998  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
999    
1000    c     ==============================================================
1001           if(ncog.eq.1)then           if(ncog.eq.1)then
1002              COG = 0.              COG = 0.
1003                if(sr1.gt.sc)cog=1.
1004                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1005    c     ==============================================================
1006           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
1007                COG = 0.
1008              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1009                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
1010              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
1011                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
1012              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1013                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1014         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1015                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1016         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1017                endif
1018    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1019    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
1020    c     ==============================================================
1021           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1022               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
1023                sss = sc
1024                if( sl1.ne.-9999. )COG = COG-sl1
1025                if( sl1.ne.-9999. )sss = sss+sl1
1026                if( sr1.ne.-9999. )COG = COG+sr1
1027                if( sr1.ne.-9999. )sss = sss+sr1
1028                if(sss.ne.0)COG=COG/sss
1029    
1030    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1031    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1032    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1033    c     ==============================================================
1034           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1035    
1036                COG = 0
1037                sss = sc
1038                if( sl1.ne.-9999. )COG = COG-sl1
1039                if( sl1.ne.-9999. )sss = sss+sl1
1040                if( sr1.ne.-9999. )COG = COG+sr1
1041                if( sr1.ne.-9999. )sss = sss+sr1
1042              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1043                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1044       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1045              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
1046                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1047       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1048                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1049                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1050         $              .and.(sl2+sss).ne.0 )
1051         $              cog = (cog-2*sl2)/(sl2+sss)
1052                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1053         $              .and.(sr2+sss).ne.0 )
1054         $              cog = (2*sr2+cog)/(sr2+sss)              
1055              endif              endif
1056    c     ==============================================================
1057             elseif(ncog.eq.5)then
1058                COG = 0
1059                sss = sc
1060                if( sl1.ne.-9999. )COG = COG-sl1
1061                if( sl1.ne.-9999. )sss = sss+sl1
1062                if( sr1.ne.-9999. )COG = COG+sr1
1063                if( sr1.ne.-9999. )sss = sss+sr1
1064                if( sl2.ne.-9999. )COG = COG-2*sl2
1065                if( sl2.ne.-9999. )sss = sss+sl2
1066                if( sr2.ne.-9999. )COG = COG+2*sr2
1067                if( sr2.ne.-9999. )sss = sss+sr2
1068                if(sss.ne.0)COG=COG/sss
1069           else           else
1070              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1071       $           ,' not implemented'       $           ,' not implemented'
# Line 696  ' Line 1082  '
1082           iv=VIEW(ic)           iv=VIEW(ic)
1083           if(mod(iv,2).eq.1)incut=incuty           if(mod(iv,2).eq.1)incut=incuty
1084           if(mod(iv,2).eq.0)incut=incutx           if(mod(iv,2).eq.0)incut=incutx
   
1085           istart = INDSTART(IC)           istart = INDSTART(IC)
1086           istop  = TOTCLLENGTH           istop  = TOTCLLENGTH
1087           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
1088           COG = 0             COG = 0  
1089             SGN = 0.
1090           mu  = 0           mu  = 0
1091           do i = istart,istop  c         print*,'-------'
1092             do i = INDMAX(IC),istart,-1
1093                ipos = i-INDMAX(ic)
1094                cut  = incut*CLSIGMA(i)
1095                if(CLSIGNAL(i).ge.cut)then              
1096                   COG = COG + ipos*CLSIGNAL(i)
1097                   SGN = SGN + CLSIGNAL(i)
1098                   mu = mu + 1
1099    c               print*,ipos,CLSIGNAL(i)
1100                else
1101                   goto 10
1102                endif
1103             enddo
1104     10      continue
1105             do i = INDMAX(IC)+1,istop
1106              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
1107              cut  = incut*CLSIGMA(i)              cut  = incut*CLSIGMA(i)
1108              if(CLSIGNAL(i).ge.cut)then              if(CLSIGNAL(i).ge.cut)then
1109                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1110                   SGN = SGN + CLSIGNAL(i)
1111                 mu = mu + 1                 mu = mu + 1
1112    c               print*,ipos,CLSIGNAL(i)
1113                else
1114                   goto 20
1115              endif              endif
1116           enddo           enddo
1117           if(SGNL(ic).le.0)then   20      continue
1118              print*,'cog(0,ic) --> ic, dedx ',ic,SGNL(ic)           if(SGN.le.0)then
1119                print*,'cog(0,ic) --> ic, dedx ',ic,SGN
1120              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
1121              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
1122              print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
1123           else           else
1124              COG=COG/SGNL(ic)              COG=COG/SGN
1125           endif           endif
1126    c         print*,'-------'
1127                    
1128        else        else
1129                    
# Line 730  ' Line 1136  '
1136    
1137  c      print *,'## cog ',ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1138    
1139          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1140             if(DEBUG.eq.1)
1141         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1142             if(DEBUG.eq.1)
1143         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1144          endif
1145    
1146    
1147        return        return
1148        end        end
1149    
1150  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1151    
1152        real function fbad_cog(ncog,ic)        real function fbad_cog(ncog,ic)
1153  *-------------------------------------------------------  *-------------------------------------------------------
1154  *     this function returns a factor that takes into  *     this function returns a factor that takes into
# Line 749  c      print *,'## cog ',ncog,ic,cog,'// Line 1164  c      print *,'## cog ',ncog,ic,cog,'//
1164        include 'calib.f'        include 'calib.f'
1165    
1166        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1167           f  = 4.           si = 8.4  !average good-strip noise
1168           si = 8.4           f  = 4.   !average bad-strip noise: f*si
1169           incut=incuty           incut=incuty
1170        else                      !X-view        else                      !X-view
1171           f  = 6.           si = 3.9  !average good-strip noise
1172           si = 3.9           f  = 6.   !average bad-strip noise: f*si
1173           incut=incutx           incut=incutx
1174        endif        endif
1175                
# Line 765  c      print *,'## cog ',ncog,ic,cog,'// Line 1180  c      print *,'## cog ',ncog,ic,cog,'//
1180  *     --> signal of the central strip  *     --> signal of the central strip
1181           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1182           fsc = 1           fsc = 1
1183  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)fsc=f  c         if( CLBAD(INDMAX(ic)).eq.0 )fsc=f
1184           if( CLBAD(INDMAX(ic)).eq.0 )fsc=f           fsc = clsigma(INDMAX(ic))/si
1185  *     --> signal of adjacent strips  *     --> signal of adjacent strips
1186           sl1 = 0                !left 1           sl1 = 0                !left 1
1187           fsl1 = 1               !left 1           fsl1 = 1               !left 1
# Line 774  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1189  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1189       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1190       $        )then       $        )then
1191              sl1 = CLSIGNAL(INDMAX(ic)-1)              sl1 = CLSIGNAL(INDMAX(ic)-1)
1192  c            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
1193              if( CLBAD(INDMAX(ic)-1).eq.0)fsl1=f              fsl1 = clsigma(INDMAX(ic)-1)/si
 c         else  
 c            fsl1 = 0  
1194           endif           endif
1195    
1196           sl2 = 0                !left 2           sl2 = 0                !left 2
# Line 786  c            fsl1 = 0 Line 1199  c            fsl1 = 0
1199       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1200       $        )then       $        )then
1201              sl2 = CLSIGNAL(INDMAX(ic)-2)              sl2 = CLSIGNAL(INDMAX(ic)-2)
1202  c            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
1203              if(CLBAD(INDMAX(ic)-2).eq.0)fsl2=f              fsl2 = clsigma(INDMAX(ic)-2)/si
 c         else  
 c            fsl2 = 0  
1204           endif           endif
1205           sr1 = 0                !right 1           sr1 = 0                !right 1
1206           fsr1 = 1               !right 1           fsr1 = 1               !right 1
# Line 799  c            fsl2 = 0 Line 1210  c            fsl2 = 0
1210       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1211       $        )then       $        )then
1212              sr1 = CLSIGNAL(INDMAX(ic)+1)              sr1 = CLSIGNAL(INDMAX(ic)+1)
1213  c            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
1214              if(CLBAD(INDMAX(ic)+1).eq.0)fsr1=f              fsr1 = clsigma(INDMAX(ic)+1)/si
 c         else  
 c            fsr1 = 0  
1215           endif               endif    
1216           sr2 = 0                !right 2           sr2 = 0                !right 2
1217           fsr2 = 1               !right 2           fsr2 = 1               !right 2
# Line 812  c            fsr1 = 0 Line 1221  c            fsr1 = 0
1221       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1222       $        )then       $        )then
1223              sr2 = CLSIGNAL(INDMAX(ic)+2)              sr2 = CLSIGNAL(INDMAX(ic)+2)
1224  c            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
1225              if(CLBAD(INDMAX(ic)+2).eq.0)fsr2=f              fsr2 = clsigma(INDMAX(ic)+2)/si
 c         else  
 c            fsr2 = 0  
1226           endif           endif
1227    
1228    
1229    
1230  ************************************************************  ************************************************************
1231  *     COG computation  *     COG2-3-4 computation
1232  ************************************************************  ************************************************************
1233    
1234  c      print*,sl2,sl1,sc,sr1,sr2  c      print*,sl2,sl1,sc,sr1,sr2
1235                    
1236           COG = 0.           vCOG = cog(ncog,ic)!0.
1237                    
1238           if(ncog.eq.2)then           if(ncog.eq.2)then
1239              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1240                 COG = -sl1/(sl1+sc)          c               COG = -sl1/(sl1+sc)        
1241                 fbad_cog = (fsl1*(-1-COG)**2+fsc*(-COG)**2)                 fbad_cog = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1242                 fbad_cog = fbad_cog / ((-1-COG)**2+(-COG)**2)                 fbad_cog = fbad_cog / ((-1-vCOG)**2+(-vCOG)**2)
1243              elseif(sl1.le.sr1)then              elseif(sl1.le.sr1)then
1244                 COG = sr1/(sc+sr1)              c               COG = sr1/(sc+sr1)            
1245                 fbad_cog = (fsc*(-COG)**2+fsr1*(1-COG)**2)                 fbad_cog = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1246                 fbad_cog = fbad_cog / ((-COG)**2+(1-COG)**2)                 fbad_cog = fbad_cog / ((-vCOG)**2+(1-vCOG)**2)
1247              endif              endif
1248           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1249              COG = (sr1-sl1)/(sl1+sc+sr1)  c            COG = (sr1-sl1)/(sl1+sc+sr1)
1250              fbad_cog =              fbad_cog =
1251       $           (fsl1*(-1-COG)**2+fsc*(-COG)**2+fsr1*(1-COG)**2)       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1252              fbad_cog =              fbad_cog =
1253       $           fbad_cog / ((-1-COG)**2+(-COG)**2+(1-COG)**2)       $           fbad_cog / ((-1-vCOG)**2+(-vCOG)**2+(1-vCOG)**2)
1254           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1255              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1256                 COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c               COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
1257                 fbad_cog =                 fbad_cog =
1258       $              (fsl2*(-2-COG)**2+fsl1*(-1-COG)**2       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1259       $              +fsc*(-COG)**2+fsr1*(1-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1260                 fbad_cog =                 fbad_cog =
1261       $              fbad_cog / ((-2-COG)**2+(-1-COG)**2       $              fbad_cog / ((-2-vCOG)**2+(-1-vCOG)**2
1262       $              +(-COG)**2+(1-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2)
1263              elseif(sl2.le.sr2)then              elseif(sl2.le.sr2)then
1264                 COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c               COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
1265                 fbad_cog =                 fbad_cog =
1266       $              (fsl1*(-1-COG)**2       $              (fsl1*(-1-vCOG)**2
1267       $              +fsc*(-COG)**2+fsr1*(1-COG)**2+fsr2*(2-COG)**2)       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1268                 fbad_cog =                 fbad_cog =
1269       $              fbad_cog / ((-1-COG)**2       $              fbad_cog / ((-1-vCOG)**2
1270       $              +(-COG)**2+(1-COG)**2+(2-COG)**2)       $              +(-vCOG)**2+(1-vCOG)**2+(2-vCOG)**2)
1271              endif              endif
1272           else           else
1273              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function FBAD_COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1274              print*,'                               (NCOG must be <= 4)'              print*,'                               (NCOG must be <= 4)'
1275              COG = 0.  c            COG = 0.
1276           endif           endif
1277                    
1278        elseif(ncog.eq.0)then        elseif(ncog.eq.0)then
1279    *     =========================
1280    *     COG computation
1281    *     =========================
1282                    
1283           iv=VIEW(ic)           vCOG = cog(0,ic)
1284           istart=INDSTART(IC)  
1285           istop=TOTCLLENGTH           iv     = VIEW(ic)
1286           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1           istart = INDSTART(IC)
1287           COG=0.           istop  = TOTCLLENGTH
1288           SNU=0.           if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1289           SDE=0.           SGN = 0.
1290           do i=istart,istop           SNU = 0.
1291              ipos=i-INDMAX(ic)           SDE = 0.
1292              il=nvk(MAXS(ic)+ipos)  
1293              is=nst(MAXS(ic)+ipos)           do i=INDMAX(IC),istart,-1
1294              cut=incut*SIGMA(iv,il,is)              ipos = i-INDMAX(ic)
1295                cut  = incut*CLSIGMA(i)
1296              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1297                 COG = COG + ipos*CLSIGNAL(i)                 fs = clsigma(i)/si
1298              endif                 SNU = SNU + fs*(ipos-vCOG)**2
1299                   SDE = SDE + (ipos-vCOG)**2
1300                else
1301                   goto 10
1302                endif            
1303           enddo           enddo
1304           COG=COG/SGNL(ic)   10      continue
1305           do i=istart,istop           do i=INDMAX(IC)+1,istop
1306              ipos=i-INDMAX(ic)              ipos = i-INDMAX(ic)
1307              il=nvk(MAXS(ic)+ipos)              cut  = incut*CLSIGMA(i)
             is=nst(MAXS(ic)+ipos)  
             cut=incut*SIGMA(iv,il,is)  
1308              if(CLSIGNAL(i).gt.cut)then              if(CLSIGNAL(i).gt.cut)then
1309                 fs=1                 fs = clsigma(i)/si
1310                 if(BAD(iv,il,is).eq.0)fs=f                 SNU = SNU + fs*(ipos-vCOG)**2
1311                 SNU = SNU + fs*(ipos-COG)**2                 SDE = SDE + (ipos-vCOG)**2
1312                 SDE = SDE + (ipos-COG)**2              else
1313                   goto 20
1314              endif                          endif            
1315           enddo           enddo
1316           if(SDE.ne.0)FBAD_COG=SNU/SDE   20      continue
1317             if(SDE.ne.0)then
1318                FBAD_COG=SNU/SDE
1319             else
1320                
1321             endif
1322    
1323        else        else
1324                                        
# Line 919  c      print*,sl2,sl1,sc,sr1,sr2 Line 1337  c      print*,sl2,sl1,sc,sr1,sr2
1337    
1338    
1339  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1340        real function fbad_cog0(ncog,ic)  
1341          real function riscogtheor(ncog,ic)
1342  *-------------------------------------------------------  *-------------------------------------------------------
 *     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.  
1343  *  *
1344  *     NB!!!  *     this function returns the expected resolution
1345  *     (this is the old version. It consider only the two  *     obtained by propagating the strip noise
1346  *     strips with the greatest signal. The new one is  *     to the center-of-gravity coordinate
1347  *     fbad_cog(ncog,ic) )  *
1348  *      *     ncog = n.strip used in the coordinate evaluation
1349    *     (ncog=0 => all strips above threshold)
1350    *
1351  *-------------------------------------------------------  *-------------------------------------------------------
1352    
1353        include 'commontracker.f'        include 'commontracker.f'
1354        include 'level1.f'        include 'level1.f'
1355        include 'calib.f'        include 'calib.f'
1356    
1357          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1358             incut = incuty
1359             pitch = pitchY / 1.e4
1360          else                      !X-view
1361             incut = incutx
1362             pitch = pitchX / 1.e4
1363          endif
1364          
1365          func = 100000.
1366          stot = 0.
1367    
1368          if (ncog.gt.0) then
1369    
1370  *     --> signal of the central strip  *     --> signal of the central strip
1371        sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1372             fsc = clsigma(INDMAX(ic))
1373    *     --> signal of adjacent strips
1374             sl1 = 0                !left 1
1375             fsl1 = 1               !left 1
1376             if(
1377         $        (INDMAX(ic)-1).ge.INDSTART(ic)
1378         $        )then
1379                sl1 = CLSIGNAL(INDMAX(ic)-1)
1380                fsl1 = clsigma(INDMAX(ic)-1)
1381             endif
1382    
1383  *     signal of adjacent strips           sl2 = 0                !left 2
1384  *     --> left           fsl2 = 1               !left 2
1385        sl1 = 0                  !left 1           if(
1386        if(       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1387       $     (INDMAX(ic)-1).ge.INDSTART(ic)       $        )then
1388       $     )              sl2 = CLSIGNAL(INDMAX(ic)-2)
1389       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))              fsl2 = clsigma(INDMAX(ic)-2)
1390             endif
1391        sl2 = 0                  !left 2           sr1 = 0                !right 1
1392        if(           fsr1 = 1               !right 1
1393       $     (INDMAX(ic)-2).ge.INDSTART(ic)           if(
1394       $     )       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1395       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))       $        .or.
1396         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1397  *     --> right       $        )then
1398        sr1 = 0                  !right 1              sr1 = CLSIGNAL(INDMAX(ic)+1)
1399        if(              fsr1 = clsigma(INDMAX(ic)+1)
1400       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))           endif    
1401       $     .or.           sr2 = 0                !right 2
1402       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)           fsr2 = 1               !right 2
1403       $     )           if(
1404       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1405         $        .or.
1406        sr2 = 0                  !right 2       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1407        if(       $        )then
1408       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))              sr2 = CLSIGNAL(INDMAX(ic)+2)
1409       $     .or.              fsr2 = clsigma(INDMAX(ic)+2)
1410       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)           endif
      $     )  
      $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
1411    
1412    
       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  
1413    
1414        fbad_cog = 1.  ************************************************************
1415        f0 = 1  *     COG2-3-4 computation
1416        f1 = 1  ************************************************************
1417        f2 = 1  
1418        f3 = 1    c      print*,sl2,sl1,sc,sr1,sr2
1419        if(sl1.gt.sr1.and.sl1.gt.0.)then          
1420             vCOG = cog(ncog,ic)!0.
1421                    
1422           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f           if(ncog.eq.1)then
1423           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f              func = 1./12.
1424  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f              stot = 1.
1425             elseif(ncog.eq.2)then
1426           if(ncog.eq.2.and.sl1.ne.0)then              if(sl1.gt.sr1)then
1427              fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)                 func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1428           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then                 stot = sl1+sc
1429              fbad_cog = 1.              elseif(sl1.le.sr1)then
1430           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then                 func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1431              fbad_cog = 1.                 stot = sc+sr1
1432                endif
1433             elseif(ncog.eq.3)then
1434                func =
1435         $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1436                stot = sl1+sc+sr1
1437             elseif(ncog.eq.4)then
1438                if(sl2.gt.sr2)then
1439                   func =
1440         $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1441         $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1442                   stot = sl2+sl1+sc+sr1
1443                elseif(sl2.le.sr2)then
1444                   func =
1445         $              (fsl1*(-1-vCOG)**2
1446         $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1447                   stot = sl2+sl1+sc+sr1
1448                endif
1449           else           else
1450              fbad_cog = 1.              print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG
1451         $            ,' not implemented'
1452           endif           endif
1453                    
1454        elseif(sl1.le.sr1.and.sr1.gt.0.)then        elseif(ncog.eq.0)then
1455    *     =========================
1456    *     COG computation
1457    *     =========================
1458            
1459             vCOG = cog(0,ic)
1460    
1461             iv     = VIEW(ic)
1462             istart = INDSTART(IC)
1463             istop  = TOTCLLENGTH
1464             if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1465    ccc         SGN = 0.
1466             SNU = 0.
1467    ccc         SDE = 0.
1468    
1469           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f           do i=INDMAX(IC),istart,-1
1470           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f              ipos = i-INDMAX(ic)
1471  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f              cut  = incut*CLSIGMA(i)
1472                if(CLSIGNAL(i).gt.cut)then
1473           if(ncog.eq.2.and.sr1.ne.0)then                 fs = clsigma(i)
1474              fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)                 SNU  = SNU + fs*(ipos-vCOG)**2
1475           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then                 stot = stot + CLSIGNAL(i)
1476              fbad_cog = 1.              else
1477           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then                 goto 10
1478              fbad_cog = 1.              endif            
1479             enddo
1480     10      continue
1481             do i=INDMAX(IC)+1,istop
1482                ipos = i-INDMAX(ic)
1483                cut  = incut*CLSIGMA(i)
1484                if(CLSIGNAL(i).gt.cut)then
1485                   fs = clsigma(i)
1486                   SNU  = SNU + fs*(ipos-vCOG)**2
1487                   stot = stot + CLSIGNAL(i)
1488                else
1489                   goto 20
1490                endif            
1491             enddo
1492     20      continue
1493             if(SDE.ne.0)then
1494                FUNC=SNU
1495           else           else
1496              fbad_cog = 1.              
1497           endif           endif
1498    
1499          else
1500                      
1501             FUNC=0
1502             print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1503             print*,'                               (NCOG must be >= 0)'
1504            
1505    
1506          endif
1507    
1508    
1509          if(stot.gt.0..and.func.gt.0.)then
1510             func = sqrt(func)
1511             func = pitch * func/stot
1512        endif        endif
1513    
1514        fbad_cog0 = sqrt(fbad_cog)        riscogtheor = func
1515    
1516        return        return
1517        end        end
1518    
1519    
1520    
1521    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1522    
1523          real function risetatheor(ncog,ic,angle)
1524    *-------------------------------------------------------
1525    *
1526    *     this function returns the expected resolution
1527    *     obtained by propagating the strip noise
1528    *     to the coordinate evaluated with non-linear eta-function
1529    *
1530    *     ncog = n.strip used in the coordinate evaluation
1531    *     (ncog=0 => ncog=2,3,4 according to angle)
1532    *
1533    *-------------------------------------------------------
1534    
1535          include 'commontracker.f'
1536          include 'level1.f'
1537          include 'calib.f'
1538    
1539    
1540          func = 1.
1541    
1542          iview   = VIEW(ic)            
1543          lad     = nld(MAXS(ic),VIEW(ic))
1544    
1545    *     ------------------------------------------------
1546    *     number of strip to be used (in case of ncog = 0)
1547    *     ------------------------------------------------
1548    
1549          inoeta = 0
1550    
1551          if(mod(int(iview),2).eq.1)then !Y-view
1552    
1553             pitch = pitchY / 1.e4
1554    
1555             if(ncog.eq.0)then
1556                if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
1557                   ncog = 2
1558                elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
1559                   ncog = 3
1560                elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
1561                   ncog = 4
1562                else
1563                   ncog   = 4
1564                   inoeta = 1
1565                endif            
1566             endif
1567    
1568          else                      !X-view
1569    
1570             pitch = pitchX / 1.e4
1571    
1572             if(ncog.eq.0)then
1573                if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
1574                   ncog = 2
1575                elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
1576                   ncog = 3
1577                elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
1578                   ncog = 4
1579                else
1580                   ncog = 4
1581                   inoeta = 1
1582                endif            
1583             endif
1584    
1585          endif
1586          
1587          func = riscogtheor(ncog,ic)
1588    
1589          risetatheor = func
1590          
1591          if(inoeta.eq.1)return ! no eta correction is applied --> exit
1592          if(ncog.lt.1.or.ncog.gt.4)return
1593    
1594    *     ----------------
1595    *     find angular bin
1596    *     ----------------
1597    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1598          do iang=1,nangbin
1599             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1600                iangle=iang
1601                goto 98
1602             endif
1603          enddo
1604          if(DEBUG.EQ.1)print*
1605         $     ,'risetatheor *** warning *** angle out of range: ',angle
1606          if(angle.le.angL(1))iang=1
1607          if(angle.ge.angR(nangbin))iang=nangbin
1608     98   continue                  !jump here if ok
1609    
1610    *     -------------
1611    *     within +/-0.5
1612    *     -------------
1613    
1614          vcog = cog(ncog,ic)          
1615    
1616          etamin = eta2(1,iang)
1617          etamax = eta2(netaval,iang)
1618    
1619          iaddmax=10
1620          iadd=0
1621     10   continue
1622          if(vcog.lt.etamin)then
1623             vcog = vcog + 1
1624             iadd = iadd + 1
1625             if(iadd>iaddmax)goto 111
1626             goto 10
1627          endif
1628     20   continue
1629          if(vcog.gt.etamax)then
1630             vcog = vcog - 1
1631             iadd = iadd - 1
1632             if(iadd<-1*iaddmax)goto 111
1633             goto 20
1634          endif
1635          goto 1111
1636     111  continue
1637          if(DEBUG.eq.1)
1638         $     print*,'risetatheor *** warning *** anomalous cluster'
1639          if(DEBUG.eq.1)
1640         $     print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)'
1641          vcog=0
1642     1111 continue
1643    
1644    *     ------------------------------------------------
1645    *     interpolation
1646    *     ------------------------------------------------
1647    
1648    
1649          ibin = netaval
1650          do i=2,netaval        
1651             if(ncog.eq.2)eta=eta2(i,iang)
1652             if(ncog.eq.3)eta=eta3(i,iang)
1653             if(ncog.eq.4)eta=eta4(i,iang)
1654             if(eta.ge.vcog)then
1655                ibin = i
1656                goto 99
1657             endif
1658          enddo
1659     99   continue
1660    
1661          if(ncog.eq.2)then
1662             x1 = eta2(ibin-1,iang)
1663             x2 = eta2(ibin,iang)
1664             y1 = feta2(ibin-1,iview,lad,iang)
1665             y2 = feta2(ibin,iview,lad,iang)
1666          elseif(ncog.eq.3)then
1667             x1 = eta3(ibin-1,iang)
1668             x2 = eta3(ibin,iang)
1669             y1 = feta3(ibin-1,iview,lad,iang)
1670             y2 = feta3(ibin,iview,lad,iang)
1671          elseif(ncog.eq.4)then
1672             x1 = eta4(ibin-1,iang)
1673             x2 = eta4(ibin,iang)
1674             y1 = feta4(ibin-1,iview,lad,iang)
1675             y2 = feta4(ibin,iview,lad,iang)
1676          endif
1677          
1678          func = func * (y2-y1)/(x2-x1)
1679    
1680          risetatheor = func
1681    
1682          return
1683          end
1684    
1685  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1686    
1687        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1688    
1689        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1690        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1103  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1756  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1756       +/       +/
1757    
1758        V(1)= abs(x)        V(1)= abs(x)
1759          if(V(1).gt.20.)V(1)=20.
1760    
1761        HQUADF = 0.        HQUADF = 0.
1762        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1117  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1771  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1771     20 CONTINUE     20 CONTINUE
1772        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1773                
1774        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1775    
1776        END        END
1777    
1778  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1779        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1780        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1781        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1782        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1193  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1847  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1847       +/       +/
1848    
1849        V(1) =  abs(x)        V(1) =  abs(x)
1850          if(V(1).gt.20.)V(1)=20.
1851    
1852        HQUADF = 0.        HQUADF = 0.
1853        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
1854           HQDJ = 0.           HQDJ = 0.
# Line 1206  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1862  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1862     20 CONTINUE     20 CONTINUE
1863        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1864                
1865        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1866    
1867        END        END
1868  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1869        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1870        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1871        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1872        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1281  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1937  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1937       +/       +/
1938    
1939        V(1)=abs(x)        V(1)=abs(x)
1940          if(V(1).gt.20.)V(1)=20.
1941    
1942        HQUADF = 0.        HQUADF = 0.
1943        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1295  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1952  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1952     20 CONTINUE     20 CONTINUE
1953        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1954                
1955        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1956    
1957        END        END
1958  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1959        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1960        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1961        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1962        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1352  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2009  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2009       +/       +/
2010    
2011        v(1)= abs(x)        v(1)= abs(x)
2012          if(V(1).gt.20.)V(1)=20.
2013                
2014        HQUADF = 0.        HQUADF = 0.
2015        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1366  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2024  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2024     20 CONTINUE     20 CONTINUE
2025        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2026    
2027        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
2028    
2029        END        END
2030  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1418  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2076  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2076       +/       +/
2077                
2078        V(1)=abs(x)        V(1)=abs(x)
2079          if(V(1).gt.20.)V(1)=20.
2080    
2081        HQUADF = 0.        HQUADF = 0.
2082        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1498  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2157  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2157       +/       +/
2158    
2159        V(1)=abs(x)        V(1)=abs(x)
2160          if(V(1).gt.20.)V(1)=20.
2161    
2162        HQUADF = 0.        HQUADF = 0.
2163        DO 20 J = 1, NPAR        DO 20 J = 1, NPAR
# Line 1515  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2175  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2175        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
2176    
2177        END        END
2178    
2179    
2180  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2181          real function pfacorr(ic,angle)
2182    *--------------------------------------------------------------
2183    *     this function returns the landi correction for this cluster
2184    *--------------------------------------------------------------
2185          include 'commontracker.f'
2186          include 'calib.f'
2187          include 'level1.f'
2188    
2189          real angle
2190          integer iview,lad
2191    
2192          iview = VIEW(ic)            
2193          lad = nld(MAXS(ic),VIEW(ic))
2194    
2195    *     find angular bin
2196    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
2197          do iang=1,nangbin
2198             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
2199                iangle=iang
2200                goto 98
2201             endif
2202          enddo
2203          if(DEBUG.eq.1)
2204         $     print*,'pfacorr *** warning *** angle out of range: ',angle
2205          if(angle.le.angL(1))iang=1
2206          if(angle.ge.angR(nangbin))iang=nangbin
2207     98   continue                  !jump here if ok
2208    
2209          pfacorr = fcorr(iview,lad,iang)
2210    
2211          if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr
2212    
2213          
2214     100  return
2215          end
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.23