/[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.11 by pam-fi, Tue May 15 16:22:18 2007 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
3    *     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)        subroutine idtoc(ipfa,cpfa)
58                
59        integer ipfa        integer ipfa
60        character*4 cpfa  c      character*10 cpfa
61          character*4 cpfa ! EM GCC4.7
62    
63        CPFA='COG4'        CPFA='COG4'
64        if(ipfa.eq.0)CPFA='ETA'        if(ipfa.eq.0)CPFA='ETA'
65        if(ipfa.eq.2)CPFA='ETA2'        if(ipfa.eq.2)CPFA='ETA2'
66        if(ipfa.eq.3)CPFA='ETA3'        if(ipfa.eq.3)CPFA='ETA3'
67        if(ipfa.eq.4)CPFA='ETA4'        if(ipfa.eq.4)CPFA='ETA4'
68          if(ipfa.eq.5)CPFA='ETAL'
69        if(ipfa.eq.10)CPFA='COG'        if(ipfa.eq.10)CPFA='COG'
70        if(ipfa.eq.11)CPFA='COG1'        if(ipfa.eq.11)CPFA='COG1'
71        if(ipfa.eq.12)CPFA='COG2'        if(ipfa.eq.12)CPFA='COG2'
# Line 17  Line 73 
73        if(ipfa.eq.14)CPFA='COG4'        if(ipfa.eq.14)CPFA='COG4'
74                
75        end        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  *     this file contains all subroutines and functions        real function fieldcorr(iview,bbb)
112  *     that are needed for position finding algorithms  
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        integer function npfastrips(ic,PFA,angle)              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  *     thid function returns the number of strips used
334  *     to evaluate the position of a cluster, according to the p.f.a.  *     to evaluate the position of a cluster, according to the p.f.a.
# Line 35  Line 337 
337        include 'level1.f'        include 'level1.f'
338        include 'calib.f'        include 'calib.f'
339    
340        character*4 usedPFA,PFA        character*4 usedPFA
341          
342    
343    
344        usedPFA=PFA        call idtoc(pfaid,usedPFA)
345    
346        npfastrips=0        npfastrips=-1
347    
348        if(usedPFA.eq.'COG1')npfastrips=1        if(usedPFA.eq.'COG1')npfastrips=1
349        if(usedPFA.eq.'COG2')npfastrips=2        if(usedPFA.eq.'COG2')npfastrips=2
# Line 50  Line 353 
353        if(usedPFA.eq.'ETA3')npfastrips=3        if(usedPFA.eq.'ETA3')npfastrips=3
354        if(usedPFA.eq.'ETA4')npfastrips=4        if(usedPFA.eq.'ETA4')npfastrips=4
355  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
356        if(usedPFA.eq.'ETA')then        if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
357  c         print*,VIEW(ic),angle  c         print*,VIEW(ic),angle
358           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
359              if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then              if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
# Line 60  c         print*,VIEW(ic),angle Line 363  c         print*,VIEW(ic),angle
363              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
364                 npfastrips=4                 npfastrips=4
365              else              else
366                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
367              endif                                      endif                        
368           else                   !X-view           else                   !X-view
369              if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then              if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
# Line 71  c               usedPFA='COG' Line 373  c               usedPFA='COG'
373              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
374                 npfastrips=4                 npfastrips=4
375              else              else
376                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
377              endif                                      endif                        
378           endif           endif
379        endif        endif
380  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
381        if(usedPFA.eq.'COG')then        if(usedPFA.eq.'COG')then
382    
383           iv=VIEW(ic)           npfastrips=0
384           if(mod(iv,2).eq.1)incut=incuty  
385           if(mod(iv,2).eq.0)incut=incutx  c$$$         iv=VIEW(ic)
386           istart = INDSTART(IC)  c$$$         if(mod(iv,2).eq.1)incut=incuty
387           istop  = TOTCLLENGTH  c$$$         if(mod(iv,2).eq.0)incut=incutx
388           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1  c$$$         istart = INDSTART(IC)
389           mu  = 0  c$$$         istop  = TOTCLLENGTH
390           do i = INDMAX(IC),istart,-1  c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
391              ipos = i-INDMAX(ic)  c$$$         mu  = 0
392              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC),istart,-1
393              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
394                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
395                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
396              else  c$$$               mu = mu + 1
397                 goto 10  c$$$               print*,i,mu
398              endif  c$$$            else
399           enddo  c$$$               goto 10
400   10      continue  c$$$            endif
401           do i = INDMAX(IC)+1,istop  c$$$         enddo
402              ipos = i-INDMAX(ic)  c$$$ 10      continue
403              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC)+1,istop
404              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
405                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
406                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
407              else  c$$$               mu = mu + 1
408                 goto 20  c$$$               print*,i,mu
409              endif  c$$$            else
410           enddo  c$$$               goto 20
411   20      continue  c$$$            endif
412           npfastrips=mu  c$$$         enddo
413    c$$$ 20      continue
414    c$$$         npfastrips=mu
415    
416        endif        endif
417  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
418    
419  c      print*,pfastrips  c      print*,pfaid,usedPFA,angle,npfastrips
420    
421        return        return
422        end        end
# Line 136  c      print*,pfastrips Line 439  c      print*,pfastrips
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           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
443              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
444           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then  cc            print*,pfaeta2(ic,angle)
445             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
446              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
447           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
448              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
449           else           else
450              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 148  c      print*,pfastrips Line 452  c      print*,pfastrips
452    
453        else                      !X-view        else                      !X-view
454    
455           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
456              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
457           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
458              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
459           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
460              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
461           else           else
462              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 160  c      print*,pfastrips Line 464  c      print*,pfastrips
464                            
465        endif        endif
466                
467   100  return  c 100  return
468          return
469        end        end
470    
471  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
472          real function pfaetal(ic,angle)
473    *--------------------------------------------------------------
474    *     this function returns the position (in strip units)
475    *     it calls:
476    *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
477    *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
478    *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
479    *     according to the angle
480    *--------------------------------------------------------------
481          include 'commontracker.f'
482          include 'level1.f'
483          include 'calib.f'
484          
485          pfaetal = 0
486    
487          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
488          
489             if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
490                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
501    
502             if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
503                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
504    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
505             elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
506                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
507             elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
508                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
509             else
510                pfaetal = cog(4,ic)
511             endif            
512                
513          endif
514          
515    c 100  return
516          return
517          end
518    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
519  c      real function riseta(ic,angle)  c      real function riseta(ic,angle)
520        real function riseta(iview,angle)        real function riseta(iview,angle)
521  *--------------------------------------------------------------  *--------------------------------------------------------------
522  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
523  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
524  *     it calls:  *     it calls:
525  *     - risx_eta2(angle)  *     - risxeta2(angle)
526  *     - risy_eta2(angle)  *     - risyeta2(angle)
527  *     - risx_eta3(angle)  *     - risxeta3(angle)
528  *     - risx_eta4(angle)  *     - risxeta4(angle)
529  *     according to the angle  *     according to the angle
530  *--------------------------------------------------------------  *--------------------------------------------------------------
531        include 'commontracker.f'        include 'commontracker.f'
532        include 'level1.f'        include 'level1.f'
533        include 'calib.f'        include 'calib.f'
534    
535        riseta = 0        riseta = 0.
536    
537  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
538        if(mod(iview,2).eq.1)then !Y-view        if(mod(iview,2).eq.1)then !Y-view
539                
540    
541           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
542              riseta = risy_eta2(angle)              riseta = risyeta2(angle)
543           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
544              riseta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
545           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
# Line 199  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 551  c      if(mod(int(VIEW(ic)),2).eq.1)then
551        else                      !X-view        else                      !X-view
552    
553           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
554              riseta = risx_eta2(angle)              riseta = risxeta2(angle)
555           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
556              riseta = risx_eta3(angle)              riseta = risxeta3(angle)
557           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
558              riseta = risx_eta4(angle)              riseta = risxeta4(angle)
559           else           else
560              riseta = risx_cog(angle)              riseta = risx_cog(angle)
561           endif                       endif            
562                            
563        endif        endif
564    
       print*,'---- ',riseta,iview,angle  
565    
566   100  return  c 100  return
567          return
568        end        end
569    
570  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 225  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 577  c      if(mod(int(VIEW(ic)),2).eq.1)then
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'
# Line 262  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 621  c      if(mod(int(VIEW(ic)),2).eq.1)then
621        end        end
622    
623  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
624        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle)
625  *--------------------------------------------------------------  *--------------------------------------------------------------
626  *     this function returns  *     this function returns
627  *  *
# Line 281  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 640  c      if(mod(int(VIEW(ic)),2).eq.1)then
640        real cog2,angle        real cog2,angle
641        integer iview,lad        integer iview,lad
642    
643        iview = VIEW(ic)                    iview   = VIEW(ic)            
644        lad = nld(MAXS(ic),VIEW(ic))        lad     = nld(MAXS(ic),VIEW(ic))
645        cog2 = cog(2,ic)                  cog2    = cog(2,ic)          
646        pfaeta2=cog2        pfaeta2 = 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
653           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
# Line 294  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 655  c      if(mod(int(VIEW(ic)),2).eq.1)then
655              goto 98              goto 98
656           endif           endif
657        enddo        enddo
658        if(DEBUG)        if(DEBUG.EQ.1)
659       $     print*,'pfaeta2 *** 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*,'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  
   
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 370  c$$$         pfaeta2=pfaeta2+1.   !temp Line 721  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,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
726    
727    
728   100  return  c 100  return
729          return
730        end        end
731    
732  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 400  c$$$      endif Line 752  c$$$      endif
752    
753        iview = VIEW(ic)                    iview = VIEW(ic)            
754        lad = nld(MAXS(ic),VIEW(ic))        lad = nld(MAXS(ic),VIEW(ic))
755        cog3 = cog(3,ic)                    cog3 = cog(3,ic)
756          cc = cog3
757          cog3 = cc
758        pfaeta3=cog3        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 412  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*,'pfaeta3 *** 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 478  c            print*,'-----',x1,x2,y1,y2 Line 825  c            print*,'-----',x1,x2,y1,y2
825        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
826        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
827    
 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  
828    
829        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
830       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
831    
832   100  return  c 100  return
833          return
834        end        end
835    
836  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 519  c$$$      endif Line 859  c$$$      endif
859        cog4=cog(4,ic)        cog4=cog(4,ic)
860        pfaeta4=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 528  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*,'pfaeta4 *** 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 603  c$$$         pfaeta2=pfaeta2+1.   !temp Line 936  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,' -->',pfaeta4       $     ,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 738  c      print *,ncog,ic,cog,'//////////// Line 1022  c      print *,ncog,ic,cog,'////////////
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 758  c      print *,ncog,ic,cog,'//////////// Line 1042  c      print *,ncog,ic,cog,'////////////
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 768  c      print *,ncog,ic,cog,'//////////// Line 1052  c      print *,ncog,ic,cog,'////////////
1052                    
1053           COG = 0.           COG = 0.
1054                    
1055  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
1056    
1057    c     ==============================================================
1058           if(ncog.eq.1)then           if(ncog.eq.1)then
1059              COG = 0.              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           elseif(ncog.eq.2)then
1064                COG = 0.
1065              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1066                 if((sl1+sc).ne.0)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                 if((sc+sr1).ne.0)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               if((sl1+sc+sr1).ne.0)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                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1101       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1102              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
1103                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1104       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              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       $           ,' not implemented'       $           ,' not implemented'
# Line 802  ' 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
# Line 818  c         print*,'-------' Line 1153  c         print*,'-------'
1153                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1154                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
1155                 mu = mu + 1                 mu = mu + 1
1156                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
1157              else              else
1158                 goto 10                 goto 10
1159              endif              endif
# Line 831  c         print*,'-------' Line 1166  c         print*,'-------'
1166                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
1167                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
1168                 mu = mu + 1                 mu = mu + 1
1169                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
1170              else              else
1171                 goto 20                 goto 20
1172              endif              endif
1173           enddo           enddo
1174   20      continue   20      continue
1175           if(SGN.le.0)then           if(SGN.le.0)then
1176  c            print*,'cog(0,ic) --> ic, dedx ',ic,SGN              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  c            print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
# Line 858  c         print*,'-------' Line 1193  c         print*,'-------'
1193    
1194  c      print *,'## cog ',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    
# Line 880  c      print *,'## cog ',ncog,ic,cog,'// Line 1223  c      print *,'## cog ',ncog,ic,cog,'//
1223        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1224           si = 8.4  !average good-strip noise           si = 8.4  !average good-strip noise
1225           f  = 4.   !average bad-strip noise: f*si           f  = 4.   !average bad-strip noise: f*si
1226           incut=incuty           incut=NINT(incuty)
1227        else                      !X-view        else                      !X-view
1228           si = 3.9  !average good-strip noise           si = 3.9  !average good-strip noise
1229           f  = 6.   !average bad-strip noise: f*si           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 1003  c            COG = 0. Line 1346  c            COG = 0.
1346           SGN = 0.           SGN = 0.
1347           SNU = 0.           SNU = 0.
1348           SDE = 0.           SDE = 0.
 c$$$         do i=INDMAX(IC),istart,-1  
 c$$$            ipos = i-INDMAX(ic)  
 c$$$            cut  = incut*CLSIGMA(i)  
 c$$$            if(CLSIGNAL(i).gt.cut)then  
 c$$$               COG = COG + ipos*CLSIGNAL(i)  
 c$$$               SGN = SGN + CLSIGNAL(i)  
 c$$$            else  
 c$$$               goto 10  
 c$$$            endif  
 c$$$         enddo  
 c$$$ 10      continue  
 c$$$         do i=INDMAX(IC)+1,istop  
 c$$$            ipos = i-INDMAX(ic)  
 c$$$            cut  = incut*CLSIGMA(i)  
 c$$$            if(CLSIGNAL(i).gt.cut)then  
 c$$$               COG = COG + ipos*CLSIGNAL(i)  
 c$$$               SGN = SGN + CLSIGNAL(i)  
 c$$$            else  
 c$$$               goto 20  
 c$$$            endif  
 c$$$         enddo  
 c$$$ 20      continue  
 c$$$         if(SGN.le.0)then  
 c$$$            print*,'fbad_cog(0,ic) --> ic, dedx ',ic,SGN  
 c$$$            print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)  
 c$$$            print*,(CLSIGNAL(i),i=istart,istop)  
 c$$$            print*,'fbad_cog(0,ic) --> NOT EVALUATED '  
 c$$$         else  
 c$$$            COG=COG/SGN  
 c$$$         endif  
1349    
1350           do i=INDMAX(IC),istart,-1           do i=INDMAX(IC),istart,-1
1351              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
# Line 1081  c$$$         endif Line 1394  c$$$         endif
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
1476        if(sl1.gt.sr1.and.sl1.gt.0.)then          
1477             vCOG = cog(ncog,ic)!0.
1478                    
1479           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f           if(ncog.eq.1)then
1480           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f              func = 1./12.
1481  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f              stot = 1.
1482             elseif(ncog.eq.2)then
1483           if(ncog.eq.2.and.sl1.ne.0)then              if(sl1.gt.sr1)then
1484              fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)                 func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1485           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then                 stot = sl1+sc
1486              fbad_cog = 1.              elseif(sl1.le.sr1)then
1487           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)
1488              fbad_cog = 1.                 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 1280  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 1371  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 1461  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 1533  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 1600  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 1681  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.11  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.23