/[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.12 by pam-fi, Wed May 23 15:31:19 2007 UTC revision 1.26 by mocchiut, Tue Aug 11 12:58:25 2009 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        character*10 cpfa
61    
62        CPFA='COG4'        CPFA='COG4'
63        if(ipfa.eq.0)CPFA='ETA'        if(ipfa.eq.0)CPFA='ETA'
64        if(ipfa.eq.2)CPFA='ETA2'        if(ipfa.eq.2)CPFA='ETA2'
65        if(ipfa.eq.3)CPFA='ETA3'        if(ipfa.eq.3)CPFA='ETA3'
66        if(ipfa.eq.4)CPFA='ETA4'        if(ipfa.eq.4)CPFA='ETA4'
67          if(ipfa.eq.5)CPFA='ETAL'
68        if(ipfa.eq.10)CPFA='COG'        if(ipfa.eq.10)CPFA='COG'
69        if(ipfa.eq.11)CPFA='COG1'        if(ipfa.eq.11)CPFA='COG1'
70        if(ipfa.eq.12)CPFA='COG2'        if(ipfa.eq.12)CPFA='COG2'
# Line 17  Line 72 
72        if(ipfa.eq.14)CPFA='COG4'        if(ipfa.eq.14)CPFA='COG4'
73                
74        end        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  *     this file contains all subroutines and functions        real function fieldcorr(iview,bbb)
110  *     that are needed for position finding algorithms  
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        integer function npfastrips(ic,PFA,angle)           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  *     thid function returns the number of strips used
330  *     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 333 
333        include 'level1.f'        include 'level1.f'
334        include 'calib.f'        include 'calib.f'
335    
336        character*4 usedPFA,PFA        character*4 usedPFA
337          
338    
339    
340        usedPFA=PFA        call idtoc(pfaid,usedPFA)
341    
342        npfastrips=0        npfastrips=-1
343    
344        if(usedPFA.eq.'COG1')npfastrips=1        if(usedPFA.eq.'COG1')npfastrips=1
345        if(usedPFA.eq.'COG2')npfastrips=2        if(usedPFA.eq.'COG2')npfastrips=2
# Line 50  Line 349 
349        if(usedPFA.eq.'ETA3')npfastrips=3        if(usedPFA.eq.'ETA3')npfastrips=3
350        if(usedPFA.eq.'ETA4')npfastrips=4        if(usedPFA.eq.'ETA4')npfastrips=4
351  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
352        if(usedPFA.eq.'ETA')then        if(usedPFA.eq.'ETA'.or.usedPFA.eq.'ETAL')then
353  c         print*,VIEW(ic),angle  c         print*,VIEW(ic),angle
354           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view           if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
355              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 359  c         print*,VIEW(ic),angle
359              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then              elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
360                 npfastrips=4                 npfastrips=4
361              else              else
362                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
363              endif                                      endif                        
364           else                   !X-view           else                   !X-view
365              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 369  c               usedPFA='COG'
369              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then              elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
370                 npfastrips=4                 npfastrips=4
371              else              else
372                 npfastrips=4                 npfastrips=4     !COG4
 c               usedPFA='COG'  
373              endif                                      endif                        
374           endif           endif
375        endif        endif
376  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
377        if(usedPFA.eq.'COG')then        if(usedPFA.eq.'COG')then
378    
379           iv=VIEW(ic)           npfastrips=0
380           if(mod(iv,2).eq.1)incut=incuty  
381           if(mod(iv,2).eq.0)incut=incutx  c$$$         iv=VIEW(ic)
382           istart = INDSTART(IC)  c$$$         if(mod(iv,2).eq.1)incut=incuty
383           istop  = TOTCLLENGTH  c$$$         if(mod(iv,2).eq.0)incut=incutx
384           if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1  c$$$         istart = INDSTART(IC)
385           mu  = 0  c$$$         istop  = TOTCLLENGTH
386           do i = INDMAX(IC),istart,-1  c$$$         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
387              ipos = i-INDMAX(ic)  c$$$         mu  = 0
388              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC),istart,-1
389              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
390                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
391                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
392              else  c$$$               mu = mu + 1
393                 goto 10  c$$$               print*,i,mu
394              endif  c$$$            else
395           enddo  c$$$               goto 10
396   10      continue  c$$$            endif
397           do i = INDMAX(IC)+1,istop  c$$$         enddo
398              ipos = i-INDMAX(ic)  c$$$ 10      continue
399              cut  = incut*CLSIGMA(i)  c$$$         do i = INDMAX(IC)+1,istop
400              if(CLSIGNAL(i).ge.cut)then  c$$$            ipos = i-INDMAX(ic)
401                 mu = mu + 1  c$$$            cut  = incut*CLSIGMA(i)
402                 print*,i,mu  c$$$            if(CLSIGNAL(i).ge.cut)then
403              else  c$$$               mu = mu + 1
404                 goto 20  c$$$               print*,i,mu
405              endif  c$$$            else
406           enddo  c$$$               goto 20
407   20      continue  c$$$            endif
408           npfastrips=mu  c$$$         enddo
409    c$$$ 20      continue
410    c$$$         npfastrips=mu
411    
412        endif        endif
413  *     ----------------------------------------------------------------  *     ----------------------------------------------------------------
414    
415  c      print*,pfastrips  c      print*,pfaid,usedPFA,angle,npfastrips
416    
417        return        return
418        end        end
# Line 136  c      print*,pfastrips Line 435  c      print*,pfastrips
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           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
439              pfaeta = pfaeta2(ic,angle)              pfaeta = pfaeta2(ic,angle)
440           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then  cc            print*,pfaeta2(ic,angle)
441             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
442              pfaeta = pfaeta3(ic,angle)              pfaeta = pfaeta3(ic,angle)
443           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
444              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
445           else           else
446              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 148  c      print*,pfastrips Line 448  c      print*,pfastrips
448    
449        else                      !X-view        else                      !X-view
450    
451           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )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).ge.e3fax.and.abs(angle).le.e3tax )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).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
456              pfaeta = pfaeta4(ic,angle)              pfaeta = pfaeta4(ic,angle)
457           else           else
458              pfaeta = cog(4,ic)              pfaeta = cog(4,ic)
# Line 160  c      print*,pfastrips Line 460  c      print*,pfastrips
460                            
461        endif        endif
462                
463   100  return  c 100  return
464          return
465        end        end
466    
467  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
468          real function pfaetal(ic,angle)
469    *--------------------------------------------------------------
470    *     this function returns the position (in strip units)
471    *     it calls:
472    *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
473    *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
474    *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
475    *     according to the angle
476    *--------------------------------------------------------------
477          include 'commontracker.f'
478          include 'level1.f'
479          include 'calib.f'
480          
481          pfaetal = 0
482    
483          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
484          
485             if( abs(angle).ge.e2fay.and.abs(angle).lt.e2tay )then
486                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
487    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
488             elseif( abs(angle).ge.e3fay.and.abs(angle).lt.e3tay )then
489                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
490             elseif( abs(angle).ge.e4fay.and.abs(angle).lt.e4tay )then
491                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
492             else
493                pfaetal = cog(4,ic)
494             endif            
495    
496          else                      !X-view
497    
498             if( abs(angle).ge.e2fax.and.abs(angle).lt.e2tax )then
499                pfaetal = pfaeta2(ic,angle)+pfacorr(ic,angle)
500    cc            print*,VIEW(ic),angle,pfaeta2(ic,angle),pfacorr(ic,angle)
501             elseif( abs(angle).ge.e3fax.and.abs(angle).lt.e3tax )then
502                pfaetal = pfaeta3(ic,angle)+pfacorr(ic,angle)
503             elseif( abs(angle).ge.e4fax.and.abs(angle).lt.e4tax )then
504                pfaetal = pfaeta4(ic,angle)+pfacorr(ic,angle)
505             else
506                pfaetal = cog(4,ic)
507             endif            
508                
509          endif
510          
511    c 100  return
512          return
513          end
514    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
515  c      real function riseta(ic,angle)  c      real function riseta(ic,angle)
516        real function riseta(iview,angle)        real function riseta(iview,angle)
517  *--------------------------------------------------------------  *--------------------------------------------------------------
518  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
519  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
520  *     it calls:  *     it calls:
521  *     - risx_eta2(angle)  *     - risxeta2(angle)
522  *     - risy_eta2(angle)  *     - risyeta2(angle)
523  *     - risx_eta3(angle)  *     - risxeta3(angle)
524  *     - risx_eta4(angle)  *     - risxeta4(angle)
525  *     according to the angle  *     according to the angle
526  *--------------------------------------------------------------  *--------------------------------------------------------------
527        include 'commontracker.f'        include 'commontracker.f'
# Line 187  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 535  c      if(mod(int(VIEW(ic)),2).eq.1)then
535                
536    
537           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
538              riseta = risy_eta2(angle)              riseta = risyeta2(angle)
539           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
540              riseta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
541           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 547  c      if(mod(int(VIEW(ic)),2).eq.1)then
547        else                      !X-view        else                      !X-view
548    
549           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
550              riseta = risx_eta2(angle)              riseta = risxeta2(angle)
551           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
552              riseta = risx_eta3(angle)              riseta = risxeta3(angle)
553           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
554              riseta = risx_eta4(angle)              riseta = risxeta4(angle)
555           else           else
556              riseta = risx_cog(angle)              riseta = risx_cog(angle)
557           endif                       endif            
558                            
559        endif        endif
560    
       print*,'---- ',riseta,iview,angle  
561    
562   100  return  c 100  return
563          return
564        end        end
565    
566  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 225  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 573  c      if(mod(int(VIEW(ic)),2).eq.1)then
573  *     resolution.  *     resolution.
574  *     It calls the function FBAD_COG(NCOG,IC),  *     It calls the function FBAD_COG(NCOG,IC),
575  *     accordingto the angle  *     accordingto the angle
576    *
577    *     >>> cosi` non e` corretto!!
578    *     >>> l'errore sulla coordinata eta si ottiene moltiplicando
579    *     >>> l'errore sulla coordinata cog per la derivata della
580    *     >>> distribuzione eta... pur sapendolo l'ho sempre ignorato...
581    *     >>> deve essere modificato!!!!
582    *
583  *-------------------------------------------------------  *-------------------------------------------------------
584    
585        include 'commontracker.f'        include 'commontracker.f'
# Line 262  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 617  c      if(mod(int(VIEW(ic)),2).eq.1)then
617        end        end
618    
619  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
620        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle)
621  *--------------------------------------------------------------  *--------------------------------------------------------------
622  *     this function returns  *     this function returns
623  *  *
# Line 281  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 636  c      if(mod(int(VIEW(ic)),2).eq.1)then
636        real cog2,angle        real cog2,angle
637        integer iview,lad        integer iview,lad
638    
639        iview = VIEW(ic)                    iview   = VIEW(ic)            
640        lad = nld(MAXS(ic),VIEW(ic))        lad     = nld(MAXS(ic),VIEW(ic))
641        cog2 = cog(2,ic)                  cog2    = cog(2,ic)          
642        pfaeta2=cog2        pfaeta2 = cog2
643    
644    *     ----------------
645  *     find angular bin  *     find angular bin
646    *     ----------------
647  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
648        do iang=1,nangbin        do iang=1,nangbin
649           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 651  c      if(mod(int(VIEW(ic)),2).eq.1)then
651              goto 98              goto 98
652           endif           endif
653        enddo        enddo
654        if(DEBUG)        if(DEBUG.EQ.1)
655       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
656        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
657        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
658   98   continue                  !jump here if ok   98   continue                  !jump here if ok
659    
660    *     -------------
661    *     within +/-0.5
662    *     -------------
663    
664  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  
   
665        iadd=0        iadd=0
666   10   continue   10   continue
667        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
668           cog2 = cog2 + 1           cog2 = cog2 + 1
669           iadd = iadd + 1           iadd = iadd + 1
670             if(iadd>iaddmax)goto 111
671           goto 10           goto 10
672        endif        endif
673   20   continue   20   continue
674        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
675           cog2 = cog2 - 1           cog2 = cog2 - 1
676           iadd = iadd - 1           iadd = iadd - 1
677             if(iadd<-1*iaddmax)goto 111
678           goto 20           goto 20
679        endif        endif
680          goto 1111
681     111  continue
682          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
683          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
684          cog2=0
685     1111 continue
686    
687  *     --------------------------------  *     --------------------------------
688  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 370  c$$$         pfaeta2=pfaeta2+1.   !temp Line 717  c$$$         pfaeta2=pfaeta2+1.   !temp
717  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
718  c$$$      endif  c$$$      endif
719    
720        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
721       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
722    
723    
724   100  return  c 100  return
725          return
726        end        end
727    
728  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 400  c$$$      endif Line 748  c$$$      endif
748    
749        iview = VIEW(ic)                    iview = VIEW(ic)            
750        lad = nld(MAXS(ic),VIEW(ic))        lad = nld(MAXS(ic),VIEW(ic))
751        cog3 = cog(3,ic)                    cog3 = cog(3,ic)
752          cc = cog3
753          cog3 = cc
754        pfaeta3=cog3        pfaeta3=cog3
755    
756    *     ----------------
757  *     find angular bin  *     find angular bin
758    *     ----------------
759  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
760        do iang=1,nangbin        do iang=1,nangbin
761  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 412  c         print*,'~~~~~~~~~~~~ ',iang,an Line 764  c         print*,'~~~~~~~~~~~~ ',iang,an
764              goto 98              goto 98
765           endif           endif
766        enddo        enddo
767        if(DEBUG)        if(DEBUG.EQ.1)
768       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
769        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
770        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
771   98   continue                  !jump here if ok   98   continue                  !jump here if ok
772    
773    *     -------------
774    *     within +/-0.5
775    *     -------------
776    
777  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  
   
         
778        iadd=0        iadd=0
779   10   continue   10   continue
780        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
781           cog3 = cog3 + 1           cog3   = cog3 + 1.
782           iadd = iadd + 1           iadd = iadd + 1
783             if(iadd>iaddmax) goto 111
784           goto 10           goto 10
785        endif        endif
786   20   continue   20   continue
787        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
788           cog3 = cog3 - 1           cog3 = cog3 - 1.
789           iadd = iadd - 1           iadd = iadd - 1
790             if(iadd<-1*iaddmax) goto 111
791           goto 20           goto 20
792        endif        endif
793          goto 1111
794     111  continue
795          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
796          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
797          cog3=0      
798     1111 continue
799    
800  *     --------------------------------  *     --------------------------------
801  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 821  c            print*,'-----',x1,x2,y1,y2
821        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
822        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
823    
 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  
824    
825        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
826       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
827    
828   100  return  c 100  return
829          return
830        end        end
831    
832  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 519  c$$$      endif Line 855  c$$$      endif
855        cog4=cog(4,ic)        cog4=cog(4,ic)
856        pfaeta4=cog4        pfaeta4=cog4
857    
858    *     ----------------
859  *     find angular bin  *     find angular bin
860    *     ----------------
861  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
862        do iang=1,nangbin        do iang=1,nangbin
863  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 528  c         print*,'~~~~~~~~~~~~ ',iang,an Line 866  c         print*,'~~~~~~~~~~~~ ',iang,an
866              goto 98              goto 98
867           endif           endif
868        enddo        enddo
869        if(DEBUG)        if(DEBUG.EQ.1)
870       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
871        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
872        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
873   98   continue                  !jump here if ok   98   continue                  !jump here if ok
874    
875    *     -------------
876    *     within +/-0.5
877    *     -------------
878    
879  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  
   
         
880        iadd=0        iadd=0
881   10   continue   10   continue
882        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
883           cog4 = cog4 + 1           cog4 = cog4 + 1
884           iadd = iadd + 1           iadd = iadd + 1
885             if(iadd>iaddmax)goto 111
886           goto 10           goto 10
887        endif        endif
888   20   continue   20   continue
889        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
890           cog4 = cog4 - 1           cog4 = cog4 - 1
891           iadd = iadd - 1           iadd = iadd - 1
892             if(iadd<-1*iaddmax)goto 111
893           goto 20           goto 20
894        endif        endif
895          goto 1111
896     111  continue
897          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
898          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
899          cog4=0
900     1111 continue
901    
902  *     --------------------------------  *     --------------------------------
903  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 603  c$$$         pfaeta2=pfaeta2+1.   !temp Line 932  c$$$         pfaeta2=pfaeta2+1.   !temp
932  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
933  c$$$      endif  c$$$      endif
934    
935        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
936       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
937    
938   100  return  c 100  return
939          return
940        end        end
941    
942    
943    
944  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
       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  
   
 *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
945        real function cog(ncog,ic)        real function cog(ncog,ic)
946  *-------------------------------------------------  *-------------------------------------------------
947  *     this function returns  *     this function returns
# Line 738  c      print *,ncog,ic,cog,'//////////// Line 971  c      print *,ncog,ic,cog,'////////////
971  *     --> signal of the central strip  *     --> signal of the central strip
972           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
973  *     signal of adjacent strips  *     signal of adjacent strips
974           sl1 = 0                !left 1           sl1 = -9999.           !left 1
975           if(           if(
976       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
977       $        )       $        )
978       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
979                    
980           sl2 = 0                !left 2           sl2 = -9999.           !left 2
981           if(           if(
982       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
983       $        )       $        )
984       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
985                    
986           sr1 = 0                !right 1           sr1 = -9999.           !right 1
987           if(           if(
988       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
989       $        .or.       $        .or.
# Line 758  c      print *,ncog,ic,cog,'//////////// Line 991  c      print *,ncog,ic,cog,'////////////
991       $        )       $        )
992       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
993                    
994           sr2 = 0                !right 2           sr2 = -9999.           !right 2
995           if(           if(
996       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
997       $        .or.       $        .or.
# Line 768  c      print *,ncog,ic,cog,'//////////// Line 1001  c      print *,ncog,ic,cog,'////////////
1001                    
1002           COG = 0.           COG = 0.
1003                    
1004  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
1005    
1006    c     ==============================================================
1007           if(ncog.eq.1)then           if(ncog.eq.1)then
1008              COG = 0.              COG = 0.
1009             if(sr1.gt.sc)cog=1.
1010             if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1011    c     ==============================================================
1012           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
1013                COG = 0.
1014              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1015                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
1016              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
1017                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)
1018              endif           elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1019                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1020         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1021                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1022         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1023             endif
1024    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1025    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
1026    c     ==============================================================
1027           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
1028               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)              COG = 0
1029                sss = sc
1030                if( sl1.ne.-9999. )COG = COG-sl1
1031                if( sl1.ne.-9999. )sss = sss+sl1
1032                if( sr1.ne.-9999. )COG = COG+sr1
1033                if( sr1.ne.-9999. )sss = sss+sr1
1034                if(sss.ne.0)COG=COG/sss
1035    
1036    c            if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
1037    c     if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
1038    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
1039    c     ==============================================================
1040           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
1041    
1042                COG = 0
1043                sss = sc
1044                if( sl1.ne.-9999. )COG = COG-sl1
1045                if( sl1.ne.-9999. )sss = sss+sl1
1046                if( sr1.ne.-9999. )COG = COG+sr1
1047                if( sr1.ne.-9999. )sss = sss+sr1
1048              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
1049                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sss).ne.0)
1050       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (COG-2*sl2)/(sl2+sss)
1051              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
1052                  if((sr2+sl1+sc+sr1).ne.0)                 if((sr2+sss).ne.0)
1053       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)       $              COG = (2*sr2+COG)/(sr2+sss)
1054                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then
1055                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
1056         $              .and.(sl2+sss).ne.0 )
1057         $              cog = (cog-2*sl2)/(sl2+sss)
1058                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
1059         $              .and.(sr2+sss).ne.0 )
1060         $              cog = (2*sr2+cog)/(sr2+sss)              
1061              endif              endif
1062    c     ==============================================================
1063             elseif(ncog.eq.5)then
1064                COG = 0
1065                sss = sc
1066                if( sl1.ne.-9999. )COG = COG-sl1
1067                if( sl1.ne.-9999. )sss = sss+sl1
1068                if( sr1.ne.-9999. )COG = COG+sr1
1069                if( sr1.ne.-9999. )sss = sss+sr1
1070                if( sl2.ne.-9999. )COG = COG-2*sl2
1071                if( sl2.ne.-9999. )sss = sss+sl2
1072                if( sr2.ne.-9999. )COG = COG+2*sr2
1073                if( sr2.ne.-9999. )sss = sss+sr2
1074                if(sss.ne.0)COG=COG/sss
1075           else           else
1076              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1077       $           ,' not implemented'       $           ,' not implemented'
# Line 858  c         print*,'-------' Line 1142  c         print*,'-------'
1142    
1143  c      print *,'## cog ',ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1144    
1145          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1146             if(DEBUG.eq.1)
1147         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1148             if(DEBUG.eq.1)
1149         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1150          endif
1151    
1152    
1153        return        return
1154        end        end
1155    
# Line 1003  c            COG = 0. Line 1295  c            COG = 0.
1295           SGN = 0.           SGN = 0.
1296           SNU = 0.           SNU = 0.
1297           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  
1298    
1299           do i=INDMAX(IC),istart,-1           do i=INDMAX(IC),istart,-1
1300              ipos = i-INDMAX(ic)              ipos = i-INDMAX(ic)
# Line 1081  c$$$         endif Line 1343  c$$$         endif
1343    
1344    
1345  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1346        real function fbad_cog0(ncog,ic)  
1347          real function riscogtheor(ncog,ic)
1348  *-------------------------------------------------------  *-------------------------------------------------------
 *     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.  
1349  *  *
1350  *     NB!!!  *     this function returns the expected resolution
1351  *     (this is the old version. It consider only the two  *     obtained by propagating the strip noise
1352  *     strips with the greatest signal. The new one is  *     to the center-of-gravity coordinate
1353  *     fbad_cog(ncog,ic) )  *
1354  *      *     ncog = n.strip used in the coordinate evaluation
1355    *     (ncog=0 => all strips above threshold)
1356    *
1357  *-------------------------------------------------------  *-------------------------------------------------------
1358    
1359        include 'commontracker.f'        include 'commontracker.f'
1360        include 'level1.f'        include 'level1.f'
1361        include 'calib.f'        include 'calib.f'
1362    
1363          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1364             incut = incuty
1365             pitch = pitchY / 1.e4
1366          else                      !X-view
1367             incut = incutx
1368             pitch = pitchX / 1.e4
1369          endif
1370          
1371          func = 100000.
1372          stot = 0.
1373    
1374          if (ncog.gt.0) then
1375    
1376  *     --> signal of the central strip  *     --> signal of the central strip
1377        sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
1378             fsc = clsigma(INDMAX(ic))
1379    *     --> signal of adjacent strips
1380             sl1 = 0                !left 1
1381             fsl1 = 1               !left 1
1382             if(
1383         $        (INDMAX(ic)-1).ge.INDSTART(ic)
1384         $        )then
1385                sl1 = CLSIGNAL(INDMAX(ic)-1)
1386                fsl1 = clsigma(INDMAX(ic)-1)
1387             endif
1388    
1389  *     signal of adjacent strips           sl2 = 0                !left 2
1390  *     --> left           fsl2 = 1               !left 2
1391        sl1 = 0                  !left 1           if(
1392        if(       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1393       $     (INDMAX(ic)-1).ge.INDSTART(ic)       $        )then
1394       $     )              sl2 = CLSIGNAL(INDMAX(ic)-2)
1395       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))              fsl2 = clsigma(INDMAX(ic)-2)
1396             endif
1397        sl2 = 0                  !left 2           sr1 = 0                !right 1
1398        if(           fsr1 = 1               !right 1
1399       $     (INDMAX(ic)-2).ge.INDSTART(ic)           if(
1400       $     )       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1401       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))       $        .or.
1402         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1403  *     --> right       $        )then
1404        sr1 = 0                  !right 1              sr1 = CLSIGNAL(INDMAX(ic)+1)
1405        if(              fsr1 = clsigma(INDMAX(ic)+1)
1406       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))           endif    
1407       $     .or.           sr2 = 0                !right 2
1408       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)           fsr2 = 1               !right 2
1409       $     )           if(
1410       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1411         $        .or.
1412        sr2 = 0                  !right 2       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1413        if(       $        )then
1414       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))              sr2 = CLSIGNAL(INDMAX(ic)+2)
1415       $     .or.              fsr2 = clsigma(INDMAX(ic)+2)
1416       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)           endif
      $     )  
      $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
1417    
1418    
       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  
1419    
1420        fbad_cog = 1.  ************************************************************
1421        f0 = 1  *     COG2-3-4 computation
1422        f1 = 1  ************************************************************
1423        f2 = 1  
1424        f3 = 1    c      print*,sl2,sl1,sc,sr1,sr2
       if(sl1.gt.sr1.and.sl1.gt.0.)then  
1425                    
1426           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f           vCOG = cog(ncog,ic)!0.
1427           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f          
1428  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f           if(ncog.eq.1)then
1429                func = 1./12.
1430           if(ncog.eq.2.and.sl1.ne.0)then              stot = 1.
1431              fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)           elseif(ncog.eq.2)then
1432           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then              if(sl1.gt.sr1)then
1433              fbad_cog = 1.                 func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1434           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then                 stot = sl1+sc
1435              fbad_cog = 1.              elseif(sl1.le.sr1)then
1436                   func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1437                   stot = sc+sr1
1438                endif
1439             elseif(ncog.eq.3)then
1440                func =
1441         $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1442                stot = sl1+sc+sr1
1443             elseif(ncog.eq.4)then
1444                if(sl2.gt.sr2)then
1445                   func =
1446         $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1447         $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1448                   stot = sl2+sl1+sc+sr1
1449                elseif(sl2.le.sr2)then
1450                   func =
1451         $              (fsl1*(-1-vCOG)**2
1452         $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1453                   stot = sl2+sl1+sc+sr1
1454                endif
1455           else           else
1456              fbad_cog = 1.              print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG
1457         $            ,' not implemented'
1458           endif           endif
1459                    
1460        elseif(sl1.le.sr1.and.sr1.gt.0.)then        elseif(ncog.eq.0)then
1461    *     =========================
1462    *     COG computation
1463    *     =========================
1464            
1465             vCOG = cog(0,ic)
1466    
1467             iv     = VIEW(ic)
1468             istart = INDSTART(IC)
1469             istop  = TOTCLLENGTH
1470             if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1471    ccc         SGN = 0.
1472             SNU = 0.
1473    ccc         SDE = 0.
1474    
1475           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f           do i=INDMAX(IC),istart,-1
1476           if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f              ipos = i-INDMAX(ic)
1477  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f              cut  = incut*CLSIGMA(i)
1478                if(CLSIGNAL(i).gt.cut)then
1479           if(ncog.eq.2.and.sr1.ne.0)then                 fs = clsigma(i)
1480              fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)                 SNU  = SNU + fs*(ipos-vCOG)**2
1481           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then                 stot = stot + CLSIGNAL(i)
1482              fbad_cog = 1.              else
1483           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then                 goto 10
1484              fbad_cog = 1.              endif            
1485             enddo
1486     10      continue
1487             do i=INDMAX(IC)+1,istop
1488                ipos = i-INDMAX(ic)
1489                cut  = incut*CLSIGMA(i)
1490                if(CLSIGNAL(i).gt.cut)then
1491                   fs = clsigma(i)
1492                   SNU  = SNU + fs*(ipos-vCOG)**2
1493                   stot = stot + CLSIGNAL(i)
1494                else
1495                   goto 20
1496                endif            
1497             enddo
1498     20      continue
1499             if(SDE.ne.0)then
1500                FUNC=SNU
1501           else           else
1502              fbad_cog = 1.              
1503           endif           endif
1504    
1505          else
1506                      
1507             FUNC=0
1508             print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1509             print*,'                               (NCOG must be >= 0)'
1510            
1511    
1512        endif        endif
1513    
1514        fbad_cog0 = sqrt(fbad_cog)  
1515          if(stot.gt.0..and.func.gt.0.)then
1516             func = sqrt(func)
1517             func = pitch * func/stot
1518          endif
1519    
1520          riscogtheor = func
1521    
1522        return        return
1523        end        end
1524    
1525    
1526    
1527    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1528    
1529          real function risetatheor(ncog,ic,angle)
1530    *-------------------------------------------------------
1531    *
1532    *     this function returns the expected resolution
1533    *     obtained by propagating the strip noise
1534    *     to the coordinate evaluated with non-linear eta-function
1535    *
1536    *     ncog = n.strip used in the coordinate evaluation
1537    *     (ncog=0 => ncog=2,3,4 according to angle)
1538    *
1539    *-------------------------------------------------------
1540    
1541          include 'commontracker.f'
1542          include 'level1.f'
1543          include 'calib.f'
1544    
1545    
1546          func = 1.
1547    
1548          iview   = VIEW(ic)            
1549          lad     = nld(MAXS(ic),VIEW(ic))
1550    
1551    *     ------------------------------------------------
1552    *     number of strip to be used (in case of ncog = 0)
1553    *     ------------------------------------------------
1554    
1555          inoeta = 0
1556    
1557          if(mod(int(iview),2).eq.1)then !Y-view
1558    
1559             pitch = pitchY / 1.e4
1560    
1561             if(ncog.eq.0)then
1562                if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
1563                   ncog = 2
1564                elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
1565                   ncog = 3
1566                elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
1567                   ncog = 4
1568                else
1569                   ncog   = 4
1570                   inoeta = 1
1571                endif            
1572             endif
1573    
1574          else                      !X-view
1575    
1576             pitch = pitchX / 1.e4
1577    
1578             if(ncog.eq.0)then
1579                if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
1580                   ncog = 2
1581                elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
1582                   ncog = 3
1583                elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
1584                   ncog = 4
1585                else
1586                   ncog = 4
1587                   inoeta = 1
1588                endif            
1589             endif
1590    
1591          endif
1592          
1593          func = riscogtheor(ncog,ic)
1594    
1595          risetatheor = func
1596          
1597          if(inoeta.eq.1)return ! no eta correction is applied --> exit
1598          if(ncog.lt.1.or.ncog.gt.4)return
1599    
1600    *     ----------------
1601    *     find angular bin
1602    *     ----------------
1603    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1604          do iang=1,nangbin
1605             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1606                iangle=iang
1607                goto 98
1608             endif
1609          enddo
1610          if(DEBUG.EQ.1)print*
1611         $     ,'risetatheor *** warning *** angle out of range: ',angle
1612          if(angle.le.angL(1))iang=1
1613          if(angle.ge.angR(nangbin))iang=nangbin
1614     98   continue                  !jump here if ok
1615    
1616    *     -------------
1617    *     within +/-0.5
1618    *     -------------
1619    
1620          vcog = cog(ncog,ic)          
1621    
1622          etamin = eta2(1,iang)
1623          etamax = eta2(netaval,iang)
1624    
1625          iaddmax=10
1626          iadd=0
1627     10   continue
1628          if(vcog.lt.etamin)then
1629             vcog = vcog + 1
1630             iadd = iadd + 1
1631             if(iadd>iaddmax)goto 111
1632             goto 10
1633          endif
1634     20   continue
1635          if(vcog.gt.etamax)then
1636             vcog = vcog - 1
1637             iadd = iadd - 1
1638             if(iadd<-1*iaddmax)goto 111
1639             goto 20
1640          endif
1641          goto 1111
1642     111  continue
1643          if(DEBUG.eq.1)
1644         $     print*,'risetatheor *** warning *** anomalous cluster'
1645          if(DEBUG.eq.1)
1646         $     print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)'
1647          vcog=0
1648     1111 continue
1649    
1650    *     ------------------------------------------------
1651    *     interpolation
1652    *     ------------------------------------------------
1653    
1654    
1655          ibin = netaval
1656          do i=2,netaval        
1657             if(ncog.eq.2)eta=eta2(i,iang)
1658             if(ncog.eq.3)eta=eta3(i,iang)
1659             if(ncog.eq.4)eta=eta4(i,iang)
1660             if(eta.ge.vcog)then
1661                ibin = i
1662                goto 99
1663             endif
1664          enddo
1665     99   continue
1666    
1667          if(ncog.eq.2)then
1668             x1 = eta2(ibin-1,iang)
1669             x2 = eta2(ibin,iang)
1670             y1 = feta2(ibin-1,iview,lad,iang)
1671             y2 = feta2(ibin,iview,lad,iang)
1672          elseif(ncog.eq.3)then
1673             x1 = eta3(ibin-1,iang)
1674             x2 = eta3(ibin,iang)
1675             y1 = feta3(ibin-1,iview,lad,iang)
1676             y2 = feta3(ibin,iview,lad,iang)
1677          elseif(ncog.eq.4)then
1678             x1 = eta4(ibin-1,iang)
1679             x2 = eta4(ibin,iang)
1680             y1 = feta4(ibin-1,iview,lad,iang)
1681             y2 = feta4(ibin,iview,lad,iang)
1682          endif
1683          
1684          func = func * (y2-y1)/(x2-x1)
1685    
1686          risetatheor = func
1687    
1688          return
1689          end
1690    
1691  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1692    
1693        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1694    
1695        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1696        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1280  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1777  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1777     20 CONTINUE     20 CONTINUE
1778        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1779                
1780        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1781    
1782        END        END
1783    
1784  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1785        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1786        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1787        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1788        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1371  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1868  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1868     20 CONTINUE     20 CONTINUE
1869        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1870                
1871        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1872    
1873        END        END
1874  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1875        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1876        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1877        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1878        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1461  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1958  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1958     20 CONTINUE     20 CONTINUE
1959        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1960                
1961        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1962    
1963        END        END
1964  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1965        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1966        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1967        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1968        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1533  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2030  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2030     20 CONTINUE     20 CONTINUE
2031        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
2032    
2033        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
2034    
2035        END        END
2036  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1684  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 2181  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
2181        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
2182    
2183        END        END
2184    
2185    
2186    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
2187          real function pfacorr(ic,angle)
2188    *--------------------------------------------------------------
2189    *     this function returns the landi correction for this cluster
2190    *--------------------------------------------------------------
2191          include 'commontracker.f'
2192          include 'calib.f'
2193          include 'level1.f'
2194    
2195          real angle
2196          integer iview,lad
2197    
2198          iview = VIEW(ic)            
2199          lad = nld(MAXS(ic),VIEW(ic))
2200    
2201    *     find angular bin
2202    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
2203          do iang=1,nangbin
2204             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
2205                iangle=iang
2206                goto 98
2207             endif
2208          enddo
2209          if(DEBUG.eq.1)
2210         $     print*,'pfacorr *** warning *** angle out of range: ',angle
2211          if(angle.le.angL(1))iang=1
2212          if(angle.ge.angR(nangbin))iang=nangbin
2213     98   continue                  !jump here if ok
2214    
2215          pfacorr = fcorr(iview,lad,iang)
2216    
2217          if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr
2218    
2219          
2220    c 100  return
2221          return
2222          end
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.23