/[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.21 by pam-fi, Fri Aug 31 14:56:52 2007 UTC revision 1.27 by pam-fi, Tue Mar 12 11:02:02 2013 UTC
# Line 8  Line 8 
8  *  *
9  *     integer function npfastrips(ic,angle)  *     integer function npfastrips(ic,angle)
10  *  *
11    *     -----------------------------------------------------------------
12    *     p.f.a.
13    *     -----------------------------------------------------------------
14  *     real function pfaeta(ic,angle)  *     real function pfaeta(ic,angle)
15  *     real function pfaetal(ic,angle)  *     real function pfaetal(ic,angle)
16  *     real function pfaeta2(ic,angle)  *     real function pfaeta2(ic,angle)
# Line 15  Line 18 
18  *     real function pfaeta4(ic,angle)  *     real function pfaeta4(ic,angle)
19  *     real function cog(ncog,ic)  *     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)  *     real function fbad_cog(ncog,ic)
36  *     real function fbad_eta(ic,angle)  *     real function fbad_eta(ic,angle)
37  *  *
38  *     real function riseta(iview,angle)  *     -----------------------------------------------------------------
39  *     FUNCTION risxeta2(x)  *     NUOVO APPROCCIO PER LA STIMA DELLA RISOLUZIONE
40  *     FUNCTION risxeta3(x)  *     -----------------------------------------------------------------
41  *     FUNCTION risxeta4(x)  *     real function riscogtheor(ncog,ic)
42  *     FUNCTION risyeta2(x)  *     real function risetatheor(ncog,ic,angle)
 *     FUNCTION risy_cog(x)  
 *     FUNCTION risx_cog(x)  
43  *  *
44    *     -----------------------------------------------------------------
45    *     correzione landi
46    *     -----------------------------------------------------------------
47  *     real function pfacorr(ic,angle)  *     real function pfacorr(ic,angle)
48  *  *
49  *     real function effectiveangle(ang,iview,bbb)  *     real function effectiveangle(ang,iview,bbb)
# Line 70  c     here bbb is the y component of the Line 88  c     here bbb is the y component of the
88           by   = bbb           by   = bbb
89           if(iview.eq.12) angx = -1. * ang           if(iview.eq.12) angx = -1. * ang
90           if(iview.eq.12) by   = -1. * bbb           if(iview.eq.12) by   = -1. * bbb
91           tgtemp  = tan(ang*acos(-1.)/180.) + pmuH_h*by*0.00001  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        elseif(mod(iview,2).eq.1)then
95  c     =================================================  c     =================================================
# Line 212  c            res = riseta(ic,ang)       Line 231  c            res = riseta(ic,ang)      
231  *     temporary patch for saturated clusters  *     temporary patch for saturated clusters
232  *     ======================================  *     ======================================
233           if( nsatstrips(ic).gt.0 )then           if( nsatstrips(ic).gt.0 )then
234              corr  = cog(4,ic)              c            corr  = cog(4,ic)            
235                corr = digsat(ic)
236              res = pitchX*1e-4/sqrt(12.)              res = pitchX*1e-4/sqrt(12.)
237  cc            cc=cog(4,ic)  cc            cc=cog(4,ic)
238  c$$$            print*,ic,' *** ',cc  c$$$            print*,ic,' *** ',cc
# Line 294  c            res = riseta(ic,ang)   Line 314  c            res = riseta(ic,ang)  
314  *     temporary patch for saturated clusters  *     temporary patch for saturated clusters
315  *     ======================================  *     ======================================
316           if( nsatstrips(ic).gt.0 )then           if( nsatstrips(ic).gt.0 )then
317              corr    = cog(4,ic)              c            corr    = cog(4,ic)            
318                corr = digsat(ic)
319              res = pitchY*1e-4/sqrt(12.)              res = pitchY*1e-4/sqrt(12.)
320  cc            cc=cog(4,ic)  cc            cc=cog(4,ic)
321  c$$$            print*,ic,' *** ',cc  c$$$            print*,ic,' *** ',cc
# Line 441  cc            print*,pfaeta2(ic,angle) Line 462  cc            print*,pfaeta2(ic,angle)
462                            
463        endif        endif
464                
465   100  return  c 100  return
466          return
467        end        end
468    
469  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 488  cc            print*,VIEW(ic),angle,pfae Line 510  cc            print*,VIEW(ic),angle,pfae
510                            
511        endif        endif
512                
513   100  return  c 100  return
514          return
515        end        end
516  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
517  c      real function riseta(ic,angle)  c      real function riseta(ic,angle)
# Line 538  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 561  c      if(mod(int(VIEW(ic)),2).eq.1)then
561        endif        endif
562    
563    
564   100  return  c 100  return
565          return
566        end        end
567    
568  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 551  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 575  c      if(mod(int(VIEW(ic)),2).eq.1)then
575  *     resolution.  *     resolution.
576  *     It calls the function FBAD_COG(NCOG,IC),  *     It calls the function FBAD_COG(NCOG,IC),
577  *     accordingto the angle  *     accordingto the angle
578    *
579    *     >>> cosi` non e` corretto!!
580    *     >>> l'errore sulla coordinata eta si ottiene moltiplicando
581    *     >>> l'errore sulla coordinata cog per la derivata della
582    *     >>> distribuzione eta... pur sapendolo l'ho sempre ignorato...
583    *     >>> deve essere modificato!!!!
584    *
585  *-------------------------------------------------------  *-------------------------------------------------------
586    
587        include 'commontracker.f'        include 'commontracker.f'
# Line 588  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 619  c      if(mod(int(VIEW(ic)),2).eq.1)then
619        end        end
620    
621  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
622        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle)
623  *--------------------------------------------------------------  *--------------------------------------------------------------
624  *     this function returns  *     this function returns
625  *  *
# Line 607  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 638  c      if(mod(int(VIEW(ic)),2).eq.1)then
638        real cog2,angle        real cog2,angle
639        integer iview,lad        integer iview,lad
640    
641        iview = VIEW(ic)                    iview   = VIEW(ic)            
642        lad = nld(MAXS(ic),VIEW(ic))        lad     = nld(MAXS(ic),VIEW(ic))
643        cog2 = cog(2,ic)                  cog2    = cog(2,ic)          
644        pfaeta2=cog2        pfaeta2 = cog2
645    
646    *     ----------------
647  *     find angular bin  *     find angular bin
648    *     ----------------
649  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
650        do iang=1,nangbin        do iang=1,nangbin
651           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then           if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
# Line 622  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 655  c      if(mod(int(VIEW(ic)),2).eq.1)then
655        enddo        enddo
656        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
657       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
658        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
659        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
660   98   continue                  !jump here if ok   98   continue                  !jump here if ok
661    
662    *     -------------
663    *     within +/-0.5
664    *     -------------
665    
666  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog2.lt.eta2(1,iang).or.cog2.gt.eta2(netaval,iang))then  
 c$$$c         print*,'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  
   
667        iadd=0        iadd=0
668   10   continue   10   continue
669        if(cog2.lt.eta2(1,iang))then        if(cog2.lt.eta2(1,iang))then
670           cog2 = cog2 + 1           cog2 = cog2 + 1
671           iadd = iadd + 1           iadd = iadd + 1
672             if(iadd>iaddmax)goto 111
673           goto 10           goto 10
674        endif        endif
675   20   continue   20   continue
676        if(cog2.gt.eta2(netaval,iang))then        if(cog2.gt.eta2(netaval,iang))then
677           cog2 = cog2 - 1           cog2 = cog2 - 1
678           iadd = iadd - 1           iadd = iadd - 1
679             if(iadd<-1*iaddmax)goto 111
680           goto 20           goto 20
681        endif        endif
682          goto 1111
683     111  continue
684          if(DEBUG.eq.1)print*,'pfaeta2 *** warning *** anomalous cluster'
685          if(DEBUG.eq.1)print*,'--> COG(2) = ',cog2-iadd,' (set to zero)'
686          cog2=0
687     1111 continue
688    
689  *     --------------------------------  *     --------------------------------
690  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 700  c$$$      endif Line 723  c$$$      endif
723       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
724    
725    
726   100  return  c 100  return
727          return
728        end        end
729    
730  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 726  c$$$      endif Line 750  c$$$      endif
750    
751        iview = VIEW(ic)                    iview = VIEW(ic)            
752        lad = nld(MAXS(ic),VIEW(ic))        lad = nld(MAXS(ic),VIEW(ic))
753        cog3 = cog(3,ic)                    cog3 = cog(3,ic)
754          cc = cog3
755          cog3 = cc
756        pfaeta3=cog3        pfaeta3=cog3
757    
758    *     ----------------
759  *     find angular bin  *     find angular bin
760    *     ----------------
761  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
762        do iang=1,nangbin        do iang=1,nangbin
763  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 740  c         print*,'~~~~~~~~~~~~ ',iang,an Line 768  c         print*,'~~~~~~~~~~~~ ',iang,an
768        enddo        enddo
769        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
770       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
771        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
772        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
773   98   continue                  !jump here if ok   98   continue                  !jump here if ok
774    
775    *     -------------
776    *     within +/-0.5
777    *     -------------
778    
779  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
780        iadd=0        iadd=0
781   10   continue   10   continue
782        if(cog3.lt.eta3(1,iang))then        if(cog3.lt.eta3(1,iang))then
783           cog3 = cog3 + 1           cog3   = cog3 + 1.
784           iadd = iadd + 1           iadd = iadd + 1
785             if(iadd>iaddmax) goto 111
786           goto 10           goto 10
787        endif        endif
788   20   continue   20   continue
789        if(cog3.gt.eta3(netaval,iang))then        if(cog3.gt.eta3(netaval,iang))then
790           cog3 = cog3 - 1           cog3 = cog3 - 1.
791           iadd = iadd - 1           iadd = iadd - 1
792             if(iadd<-1*iaddmax) goto 111
793           goto 20           goto 20
794        endif        endif
795          goto 1111
796     111  continue
797          if(DEBUG.eq.1)print*,'pfaeta3 *** warning *** anomalous cluster'
798          if(DEBUG.eq.1)print*,'--> COG(3) = ',cog3-iadd,' (set to zero)'
799          cog3=0      
800     1111 continue
801    
802  *     --------------------------------  *     --------------------------------
803  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 804  c            print*,'-----',x1,x2,y1,y2 Line 823  c            print*,'-----',x1,x2,y1,y2
823        pfaeta3 = AA*cog3+BB        pfaeta3 = AA*cog3+BB
824        pfaeta3 = pfaeta3 - iadd        pfaeta3 = pfaeta3 - iadd
825    
 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  
826    
827        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
828       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
829    
830   100  return  c 100  return
831          return
832        end        end
833    
834  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 845  c$$$      endif Line 857  c$$$      endif
857        cog4=cog(4,ic)        cog4=cog(4,ic)
858        pfaeta4=cog4        pfaeta4=cog4
859    
860    *     ----------------
861  *     find angular bin  *     find angular bin
862    *     ----------------
863  *     (in futuro possiamo pensare di interpolare anche sull'angolo)  *     (in futuro possiamo pensare di interpolare anche sull'angolo)
864        do iang=1,nangbin        do iang=1,nangbin
865  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle  c         print*,'~~~~~~~~~~~~ ',iang,angL(iang),angR(iang),angle
# Line 856  c         print*,'~~~~~~~~~~~~ ',iang,an Line 870  c         print*,'~~~~~~~~~~~~ ',iang,an
870        enddo        enddo
871        if(DEBUG.EQ.1)        if(DEBUG.EQ.1)
872       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
873        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
874        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
875   98   continue                  !jump here if ok   98   continue                  !jump here if ok
876    
877    *     -------------
878    *     within +/-0.5
879    *     -------------
880    
881  c$$$*     find extremes of interpolation        iaddmax=10
 c$$$      iflag=0  
 c$$$*     --------------------------------  
 c$$$      if(cog3.lt.eta3(1,iang).or.cog3.gt.eta3(netaval,iang))then  
 c$$$*     ----------------------------------------------  
 c$$$*     non salto piu`, ma scalo di 1 o -1  
 c$$$*     nel caso si tratti di un cluster  
 c$$$*     in cui la strip con il segnale massimo non coincide  
 c$$$*     con la strip con il rapposto s/n massimo!!!  
 c$$$*     ----------------------------------------------  
 c$$$         if(cog2.lt.eta2(1,iang))then !temp  
 c$$$            cog2=cog2+1.        !temp  
 c$$$            iflag=1             !temp  
 c$$$         else                   !temp  
 c$$$            cog2=cog2-1.        !temp  
 c$$$            iflag=-1            !temp  
 c$$$         endif                  !temp  
 c$$$c         print*,'shifted >>> ',cog2  
 c$$$      endif  
   
         
882        iadd=0        iadd=0
883   10   continue   10   continue
884        if(cog4.lt.eta4(1,iang))then        if(cog4.lt.eta4(1,iang))then
885           cog4 = cog4 + 1           cog4 = cog4 + 1
886           iadd = iadd + 1           iadd = iadd + 1
887             if(iadd>iaddmax)goto 111
888           goto 10           goto 10
889        endif        endif
890   20   continue   20   continue
891        if(cog4.gt.eta4(netaval,iang))then        if(cog4.gt.eta4(netaval,iang))then
892           cog4 = cog4 - 1           cog4 = cog4 - 1
893           iadd = iadd - 1           iadd = iadd - 1
894             if(iadd<-1*iaddmax)goto 111
895           goto 20           goto 20
896        endif        endif
897          goto 1111
898     111  continue
899          if(DEBUG.eq.1)print*,'pfaeta4 *** warning *** anomalous cluster'
900          if(DEBUG.eq.1)print*,'--> COG(4) = ',cog4-iadd,' (set to zero)'
901          cog4=0
902     1111 continue
903    
904  *     --------------------------------  *     --------------------------------
905  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2  c      print*,'*****',i,view,lad,iang,'------> cog2 ',cog2
# Line 932  c$$$      endif Line 937  c$$$      endif
937        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
938       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
939    
940   100  return  c 100  return
941          return
942        end        end
943    
944    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
945          real function digsat(ic)
946    *-------------------------------------------------
947    *
948    *          
949    *-------------------------------------------------
950          include 'commontracker.f'
951          include 'calib.f'
952          include 'level1.f'
953          
954          integer nsat
955          real pitchsat
956          
957          nsat = 0
958          pitchsat = 0.
959          iv=VIEW(ic)              
960          istart = INDSTART(IC)
961          istop  = TOTCLLENGTH
962          if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
963          do i = INDMAX(IC),istart,-1
964             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
965         $        .or.
966         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
967                nsat = nsat + 1
968                pitchsat = pitchsat + i - INDMAX(IC)
969             else
970                goto 10
971             endif
972          enddo
973     10   continue
974          do i = INDMAX(IC)+1,istop
975             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
976         $        .or.
977         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
978                nsat = nsat + 1
979                pitchsat = pitchsat + i - INDMAX(IC)
980             else
981                goto 20
982             endif
983          enddo
984     20   continue
985          
986          digsat = 0
987          if (nsat.gt.0) digsat = pitchsat / nsat
988          
989          return
990          end
991    
992    
993  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
994        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 997  c$$$      endif Line 1050  c$$$      endif
1050                    
1051           COG = 0.           COG = 0.
1052                    
1053  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c      print *,'## ',sl2,sl1,sc,sr1,sr2
1054    
1055  c     ==============================================================  c     ==============================================================
1056           if(ncog.eq.1)then           if(ncog.eq.1)then
1057              COG = 0.              COG = 0.
1058              if(sr1.gt.sc)cog=1.           if(sr1.gt.sc)cog=1.
1059              if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.           if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1.
1060  c     ==============================================================  c     ==============================================================
1061           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
1062              COG = 0.              COG = 0.
1063              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
1064                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
1065              elseif(sl1.lt.sr1)then              elseif(sl1.lt.sr1)then
1066                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                                         if((sc+sr1).ne.0)COG = sr1/(sc+sr1)
1067              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then           elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then
1068                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
1069       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)       $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc)
1070                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)                 if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
1071       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)       $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1)
1072              endif           endif
1073  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)  c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
1074  c     $           ,' : ',sl2,sl1,sc,sr1,sr2  c     $           ,' : ',sl2,sl1,sc,sr1,sr2
1075  c     ==============================================================  c     ==============================================================
# Line 1138  c         print*,'-------' Line 1191  c         print*,'-------'
1191    
1192  c      print *,'## cog ',ncog,ic,cog,'/////////////'  c      print *,'## cog ',ncog,ic,cog,'/////////////'
1193    
1194          if(COG.lt.-0.75.or.COG.gt.+0.75)then
1195             if(DEBUG.eq.1)
1196         $        print*,'cog *** warning *** anomalous cluster ??? --> '
1197             if(DEBUG.eq.1)
1198         $        print*,sl2,sl1,sc,sr1,sr2,' --> COG(',ncog,') = ',COG
1199          endif
1200    
1201    
1202        return        return
1203        end        end
1204    
# Line 1330  c            COG = 0. Line 1391  c            COG = 0.
1391        end        end
1392    
1393    
1394  c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1395  c$$$      real function fbad_cog0(ncog,ic)  
1396  c$$$*-------------------------------------------------------        real function riscogtheor(ncog,ic)
1397  c$$$*     this function returns a factor that takes into  *-------------------------------------------------------
1398  c$$$*     account deterioration of the spatial resolution  *
1399  c$$$*     in the case BAD strips are included in the cluster.  *     this function returns the expected resolution
1400  c$$$*     This factor should multiply the nominal spatial  *     obtained by propagating the strip noise
1401  c$$$*     resolution.  *     to the center-of-gravity coordinate
1402  c$$$*  *
1403  c$$$*     NB!!!  *     ncog = n.strip used in the coordinate evaluation
1404  c$$$*     (this is the old version. It consider only the two  *     (ncog=0 => all strips above threshold)
1405  c$$$*     strips with the greatest signal. The new one is  *
1406  c$$$*     fbad_cog(ncog,ic) )  *-------------------------------------------------------
1407  c$$$*      
1408  c$$$*-------------------------------------------------------        include 'commontracker.f'
1409  c$$$        include 'level1.f'
1410  c$$$      include 'commontracker.f'        include 'calib.f'
1411  c$$$      include 'level1.f'  
1412  c$$$      include 'calib.f'        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1413  c$$$           incut = incuty
1414  c$$$*     --> signal of the central strip           pitch = pitchY / 1.e4
1415  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center        else                      !X-view
1416  c$$$           incut = incutx
1417  c$$$*     signal of adjacent strips           pitch = pitchX / 1.e4
1418  c$$$*     --> left        endif
1419  c$$$      sl1 = 0                  !left 1        
1420  c$$$      if(        func = 100000.
1421  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)        stot = 0.
1422  c$$$     $     )  
1423  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))        if (ncog.gt.0) then
1424  c$$$  
1425  c$$$      sl2 = 0                  !left 2  *     --> signal of the central strip
1426  c$$$      if(           sc = CLSIGNAL(INDMAX(ic)) !center
1427  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)           fsc = clsigma(INDMAX(ic))
1428  c$$$     $     )  *     --> signal of adjacent strips
1429  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))           sl1 = 0                !left 1
1430  c$$$           fsl1 = 1               !left 1
1431  c$$$*     --> right           if(
1432  c$$$      sr1 = 0                  !right 1       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1433  c$$$      if(       $        )then
1434  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))              sl1 = CLSIGNAL(INDMAX(ic)-1)
1435  c$$$     $     .or.              fsl1 = clsigma(INDMAX(ic)-1)
1436  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)           endif
1437  c$$$     $     )  
1438  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))           sl2 = 0                !left 2
1439  c$$$           fsl2 = 1               !left 2
1440  c$$$      sr2 = 0                  !right 2           if(
1441  c$$$      if(       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1442  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        )then
1443  c$$$     $     .or.              sl2 = CLSIGNAL(INDMAX(ic)-2)
1444  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)              fsl2 = clsigma(INDMAX(ic)-2)
1445  c$$$     $     )           endif
1446  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))           sr1 = 0                !right 1
1447  c$$$           fsr1 = 1               !right 1
1448  c$$$           if(
1449  c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1450  c$$$         f  = 4.       $        .or.
1451  c$$$         si = 8.4       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1452  c$$$      else                              !X-view       $        )then
1453  c$$$         f  = 6.              sr1 = CLSIGNAL(INDMAX(ic)+1)
1454  c$$$         si = 3.9              fsr1 = clsigma(INDMAX(ic)+1)
1455  c$$$      endif           endif    
1456  c$$$           sr2 = 0                !right 2
1457  c$$$      fbad_cog = 1.           fsr2 = 1               !right 2
1458  c$$$      f0 = 1           if(
1459  c$$$      f1 = 1       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1460  c$$$      f2 = 1       $        .or.
1461  c$$$      f3 = 1         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1462  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then       $        )then
1463  c$$$                      sr2 = CLSIGNAL(INDMAX(ic)+2)
1464  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f              fsr2 = clsigma(INDMAX(ic)+2)
1465  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f           endif
1466  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  
1467  c$$$  
1468  c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
1469  c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  ************************************************************
1470  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  *     COG2-3-4 computation
1471  c$$$            fbad_cog = 1.  ************************************************************
1472  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
1473  c$$$            fbad_cog = 1.  c      print*,sl2,sl1,sc,sr1,sr2
1474  c$$$         else          
1475  c$$$            fbad_cog = 1.           vCOG = cog(ncog,ic)!0.
1476  c$$$         endif          
1477  c$$$                   if(ncog.eq.1)then
1478  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then              func = 1./12.
1479  c$$$              stot = 1.
1480  c$$$           elseif(ncog.eq.2)then
1481  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f              if(sl1.gt.sr1)then
1482  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f                 func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1483  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f                 stot = sl1+sc
1484  c$$$              elseif(sl1.le.sr1)then
1485  c$$$         if(ncog.eq.2.and.sr1.ne.0)then                 func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1486  c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)                 stot = sc+sr1
1487  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then              endif
1488  c$$$            fbad_cog = 1.           elseif(ncog.eq.3)then
1489  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then              func =
1490  c$$$            fbad_cog = 1.       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1491  c$$$         else              stot = sl1+sc+sr1
1492  c$$$            fbad_cog = 1.           elseif(ncog.eq.4)then
1493  c$$$         endif              if(sl2.gt.sr2)then
1494  c$$$                 func =
1495  c$$$      endif       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1496  c$$$       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1497  c$$$      fbad_cog0 = sqrt(fbad_cog)                 stot = sl2+sl1+sc+sr1
1498  c$$$              elseif(sl2.le.sr2)then
1499  c$$$      return                 func =
1500  c$$$      end       $              (fsl1*(-1-vCOG)**2
1501  c$$$       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1502  c$$$                 stot = sl2+sl1+sc+sr1
1503  c$$$              endif
1504             else
1505                print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG
1506         $            ,' not implemented'
1507             endif
1508            
1509          elseif(ncog.eq.0)then
1510    *     =========================
1511    *     COG computation
1512    *     =========================
1513            
1514             vCOG = cog(0,ic)
1515    
1516             iv     = VIEW(ic)
1517             istart = INDSTART(IC)
1518             istop  = TOTCLLENGTH
1519             if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1520    ccc         SGN = 0.
1521             SNU = 0.
1522    ccc         SDE = 0.
1523    
1524             do i=INDMAX(IC),istart,-1
1525                ipos = i-INDMAX(ic)
1526                cut  = incut*CLSIGMA(i)
1527                if(CLSIGNAL(i).gt.cut)then
1528                   fs = clsigma(i)
1529                   SNU  = SNU + fs*(ipos-vCOG)**2
1530                   stot = stot + CLSIGNAL(i)
1531                else
1532                   goto 10
1533                endif            
1534             enddo
1535     10      continue
1536             do i=INDMAX(IC)+1,istop
1537                ipos = i-INDMAX(ic)
1538                cut  = incut*CLSIGMA(i)
1539                if(CLSIGNAL(i).gt.cut)then
1540                   fs = clsigma(i)
1541                   SNU  = SNU + fs*(ipos-vCOG)**2
1542                   stot = stot + CLSIGNAL(i)
1543                else
1544                   goto 20
1545                endif            
1546             enddo
1547     20      continue
1548             if(SDE.ne.0)then
1549                FUNC=SNU
1550             else
1551                
1552             endif
1553    
1554          else
1555                      
1556             FUNC=0
1557             print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1558             print*,'                               (NCOG must be >= 0)'
1559            
1560    
1561          endif
1562    
1563    
1564          if(stot.gt.0..and.func.gt.0.)then
1565             func = sqrt(func)
1566             func = pitch * func/stot
1567          endif
1568    
1569          riscogtheor = func
1570    
1571          return
1572          end
1573    
1574    
1575    
1576    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1577    
1578          real function risetatheor(ncog,ic,angle)
1579    *-------------------------------------------------------
1580    *
1581    *     this function returns the expected resolution
1582    *     obtained by propagating the strip noise
1583    *     to the coordinate evaluated with non-linear eta-function
1584    *
1585    *     ncog = n.strip used in the coordinate evaluation
1586    *     (ncog=0 => ncog=2,3,4 according to angle)
1587    *
1588    *-------------------------------------------------------
1589    
1590          include 'commontracker.f'
1591          include 'level1.f'
1592          include 'calib.f'
1593    
1594    
1595          func = 1.
1596    
1597          iview   = VIEW(ic)            
1598          lad     = nld(MAXS(ic),VIEW(ic))
1599    
1600    *     ------------------------------------------------
1601    *     number of strip to be used (in case of ncog = 0)
1602    *     ------------------------------------------------
1603    
1604          inoeta = 0
1605    
1606          if(mod(int(iview),2).eq.1)then !Y-view
1607    
1608             pitch = pitchY / 1.e4
1609    
1610             if(ncog.eq.0)then
1611                if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
1612                   ncog = 2
1613                elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
1614                   ncog = 3
1615                elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
1616                   ncog = 4
1617                else
1618                   ncog   = 4
1619                   inoeta = 1
1620                endif            
1621             endif
1622    
1623          else                      !X-view
1624    
1625             pitch = pitchX / 1.e4
1626    
1627             if(ncog.eq.0)then
1628                if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
1629                   ncog = 2
1630                elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
1631                   ncog = 3
1632                elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
1633                   ncog = 4
1634                else
1635                   ncog = 4
1636                   inoeta = 1
1637                endif            
1638             endif
1639    
1640          endif
1641          
1642          func = riscogtheor(ncog,ic)
1643    
1644          risetatheor = func
1645          
1646          if(inoeta.eq.1)return ! no eta correction is applied --> exit
1647          if(ncog.lt.1.or.ncog.gt.4)return
1648    
1649    *     ----------------
1650    *     find angular bin
1651    *     ----------------
1652    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1653          do iang=1,nangbin
1654             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1655                iangle=iang
1656                goto 98
1657             endif
1658          enddo
1659          if(DEBUG.EQ.1)print*
1660         $     ,'risetatheor *** warning *** angle out of range: ',angle
1661          if(angle.le.angL(1))iang=1
1662          if(angle.ge.angR(nangbin))iang=nangbin
1663     98   continue                  !jump here if ok
1664    
1665    *     -------------
1666    *     within +/-0.5
1667    *     -------------
1668    
1669          vcog = cog(ncog,ic)          
1670    
1671          etamin = eta2(1,iang)
1672          etamax = eta2(netaval,iang)
1673    
1674          iaddmax=10
1675          iadd=0
1676     10   continue
1677          if(vcog.lt.etamin)then
1678             vcog = vcog + 1
1679             iadd = iadd + 1
1680             if(iadd>iaddmax)goto 111
1681             goto 10
1682          endif
1683     20   continue
1684          if(vcog.gt.etamax)then
1685             vcog = vcog - 1
1686             iadd = iadd - 1
1687             if(iadd<-1*iaddmax)goto 111
1688             goto 20
1689          endif
1690          goto 1111
1691     111  continue
1692          if(DEBUG.eq.1)
1693         $     print*,'risetatheor *** warning *** anomalous cluster'
1694          if(DEBUG.eq.1)
1695         $     print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)'
1696          vcog=0
1697     1111 continue
1698    
1699    *     ------------------------------------------------
1700    *     interpolation
1701    *     ------------------------------------------------
1702    
1703    
1704          ibin = netaval
1705          do i=2,netaval        
1706             if(ncog.eq.2)eta=eta2(i,iang)
1707             if(ncog.eq.3)eta=eta3(i,iang)
1708             if(ncog.eq.4)eta=eta4(i,iang)
1709             if(eta.ge.vcog)then
1710                ibin = i
1711                goto 99
1712             endif
1713          enddo
1714     99   continue
1715    
1716          if(ncog.eq.2)then
1717             x1 = eta2(ibin-1,iang)
1718             x2 = eta2(ibin,iang)
1719             y1 = feta2(ibin-1,iview,lad,iang)
1720             y2 = feta2(ibin,iview,lad,iang)
1721          elseif(ncog.eq.3)then
1722             x1 = eta3(ibin-1,iang)
1723             x2 = eta3(ibin,iang)
1724             y1 = feta3(ibin-1,iview,lad,iang)
1725             y2 = feta3(ibin,iview,lad,iang)
1726          elseif(ncog.eq.4)then
1727             x1 = eta4(ibin-1,iang)
1728             x2 = eta4(ibin,iang)
1729             y1 = feta4(ibin-1,iview,lad,iang)
1730             y2 = feta4(ibin,iview,lad,iang)
1731          endif
1732          
1733          func = func * (y2-y1)/(x2-x1)
1734    
1735          risetatheor = func
1736    
1737          return
1738          end
1739    
1740  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1741    
# Line 1961  c$$$ Line 2257  c$$$
2257        enddo        enddo
2258        if(DEBUG.eq.1)        if(DEBUG.eq.1)
2259       $     print*,'pfacorr *** warning *** angle out of range: ',angle       $     print*,'pfacorr *** warning *** angle out of range: ',angle
2260        if(angle.lt.angL(1))iang=1        if(angle.le.angL(1))iang=1
2261        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.ge.angR(nangbin))iang=nangbin
2262   98   continue                  !jump here if ok   98   continue                  !jump here if ok
2263    
2264        pfacorr = fcorr(iview,lad,iang)        pfacorr = fcorr(iview,lad,iang)
2265    
2266        if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr        if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr
   
2267    
2268   100  return        
2269    c 100  return
2270          return
2271        end        end
Legend:
Removed from v.1.21  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.23