/[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.10 by pam-fi, Mon May 14 11:03:06 2007 UTC revision 1.18 by pam-fi, Mon Aug 20 16:07:16 2007 UTC
# Line 3  Line 3 
3        subroutine idtoc(ipfa,cpfa)        subroutine idtoc(ipfa,cpfa)
4                
5        integer ipfa        integer ipfa
6        character*4 cpfa        character*10 cpfa
7    
8        CPFA='COG4'        CPFA='COG4'
9        if(ipfa.eq.0)CPFA='ETA'        if(ipfa.eq.0)CPFA='ETA'
10        if(ipfa.eq.2)CPFA='ETA2'        if(ipfa.eq.2)CPFA='ETA2'
11        if(ipfa.eq.3)CPFA='ETA3'        if(ipfa.eq.3)CPFA='ETA3'
12        if(ipfa.eq.4)CPFA='ETA4'        if(ipfa.eq.4)CPFA='ETA4'
13          if(ipfa.eq.5)CPFA='ETAL'
14        if(ipfa.eq.10)CPFA='COG'        if(ipfa.eq.10)CPFA='COG'
15        if(ipfa.eq.11)CPFA='COG1'        if(ipfa.eq.11)CPFA='COG1'
16        if(ipfa.eq.12)CPFA='COG2'        if(ipfa.eq.12)CPFA='COG2'
# Line 164  c      print*,pfastrips Line 165  c      print*,pfastrips
165        end        end
166    
167  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
168        real function ris_eta(ic,angle)        real function pfaetal(ic,angle)
169    *--------------------------------------------------------------
170    *     this function returns the position (in strip units)
171    *     it calls:
172    *     - pfaeta2(ic,angle)+pfcorr(ic,angle)
173    *     - pfaeta3(ic,angle)+pfcorr(ic,angle)
174    *     - pfaeta4(ic,angle)+pfcorr(ic,angle)
175    *     according to the angle
176    *--------------------------------------------------------------
177          include 'commontracker.f'
178          include 'level1.f'
179          include 'calib.f'
180          
181          pfaeta = 0
182    
183          if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
184          
185             if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
186                pfaeta = pfaeta2(ic,angle)+pfacorr(ic,angle)
187             elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
188                pfaeta = pfaeta3(ic,angle)+pfacorr(ic,angle)
189             elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
190                pfaeta = pfaeta4(ic,angle)+pfacorr(ic,angle)
191             else
192                pfaeta = cog(4,ic)
193             endif            
194    
195          else                      !X-view
196    
197             if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
198                pfaeta = pfaeta2(ic,angle)+pfacorr(ic,angle)
199             elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
200                pfaeta = pfaeta3(ic,angle)+pfacorr(ic,angle)
201             elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
202                pfaeta = pfaeta4(ic,angle)+pfacorr(ic,angle)
203             else
204                pfaeta = cog(4,ic)
205             endif            
206                
207          endif
208          
209     100  return
210          end
211    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
212    c      real function riseta(ic,angle)
213          real function riseta(iview,angle)
214  *--------------------------------------------------------------  *--------------------------------------------------------------
215  *     this function returns the average spatial resolution  *     this function returns the average spatial resolution
216  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))  *     (in cm) for the ETA algorithm (function pfaeta(ic,angle))
217  *     it calls:  *     it calls:
218  *     - risx_eta2(angle)  *     - risxeta2(angle)
219  *     - risy_eta2(angle)  *     - risyeta2(angle)
220  *     - risx_eta3(angle)  *     - risxeta3(angle)
221  *     - risx_eta4(angle)  *     - risxeta4(angle)
222  *     according to the angle  *     according to the angle
223  *--------------------------------------------------------------  *--------------------------------------------------------------
224        include 'commontracker.f'        include 'commontracker.f'
225        include 'level1.f'        include 'level1.f'
226        include 'calib.f'        include 'calib.f'
227    
228        ris_eta = 0        riseta = 0
229    
230        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  c      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
231          if(mod(iview,2).eq.1)then !Y-view
232                
233    
234           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then           if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
235              ris_eta = risy_eta2(angle)              riseta = risyeta2(angle)
236           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then           elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
237              ris_eta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
238           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then           elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
239              ris_eta = risy_cog(angle) !ATTENZIONE!!              riseta = risy_cog(angle) !ATTENZIONE!!
240           else           else
241              ris_eta = risy_cog(angle)              riseta = risy_cog(angle)
242           endif                       endif            
243    
244        else                      !X-view        else                      !X-view
245    
246           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then           if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
247              ris_eta = risx_eta2(angle)              riseta = risxeta2(angle)
248           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then           elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
249              ris_eta = risx_eta3(angle)              riseta = risxeta3(angle)
250           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then           elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
251              ris_eta = risx_eta4(angle)              riseta = risxeta4(angle)
252           else           else
253              ris_eta = risx_cog(angle)              riseta = risx_cog(angle)
254           endif                       endif            
255                            
256        endif        endif
257    
258    
259   100  return   100  return
260        end        end
261    
# Line 290  c      print*,pfastrips Line 338  c      print*,pfastrips
338              goto 98              goto 98
339           endif           endif
340        enddo        enddo
341        if(DEBUG)        if(DEBUG.EQ.1)
342       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle       $     print*,'pfaeta2 *** warning *** angle out of range: ',angle
343        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
344        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
# Line 366  c$$$         pfaeta2=pfaeta2+1.   !temp Line 414  c$$$         pfaeta2=pfaeta2+1.   !temp
414  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
415  c$$$      endif  c$$$      endif
416    
417        if(DEBUG)print*,'ETA2  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA2  (ic ',ic,' ang',angle,')'
418       $     ,cog2-iadd,' -->',pfaeta2       $     ,cog2-iadd,' -->',pfaeta2
419    
420    
# Line 408  c         print*,'~~~~~~~~~~~~ ',iang,an Line 456  c         print*,'~~~~~~~~~~~~ ',iang,an
456              goto 98              goto 98
457           endif           endif
458        enddo        enddo
459        if(DEBUG)        if(DEBUG.EQ.1)
460       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle       $     print*,'pfaeta3 *** warning *** angle out of range: ',angle
461        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
462        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
# Line 483  c$$$         pfaeta2=pfaeta2+1.   !temp Line 531  c$$$         pfaeta2=pfaeta2+1.   !temp
531  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
532  c$$$      endif  c$$$      endif
533    
534        if(DEBUG)print*,'ETA3  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA3  (ic ',ic,' ang',angle,')'
535       $     ,cog3-iadd,' -->',pfaeta3       $     ,cog3-iadd,' -->',pfaeta3
536    
537   100  return   100  return
# Line 524  c         print*,'~~~~~~~~~~~~ ',iang,an Line 572  c         print*,'~~~~~~~~~~~~ ',iang,an
572              goto 98              goto 98
573           endif           endif
574        enddo        enddo
575        if(DEBUG)        if(DEBUG.EQ.1)
576       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle       $     print*,'pfaeta4 *** warning *** angle out of range: ',angle
577        if(angle.lt.angL(1))iang=1        if(angle.lt.angL(1))iang=1
578        if(angle.gt.angR(nangbin))iang=nangbin        if(angle.gt.angR(nangbin))iang=nangbin
# Line 599  c$$$         pfaeta2=pfaeta2+1.   !temp Line 647  c$$$         pfaeta2=pfaeta2+1.   !temp
647  c$$$         cog2=cog2+1.           !temp  c$$$         cog2=cog2+1.           !temp
648  c$$$      endif  c$$$      endif
649    
650        if(DEBUG)print*,'ETA4  (ic ',ic,' ang',angle,')'        if(DEBUG.EQ.1)print*,'ETA4  (ic ',ic,' ang',angle,')'
651       $     ,cog4-iadd,' -->',pfaeta4       $     ,cog4-iadd,' -->',pfaeta4
652    
653   100  return   100  return
# Line 608  c$$$      endif Line 656  c$$$      endif
656    
657    
658  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
659        real function cog0(ncog,ic)  c$$$      real function cog0(ncog,ic)
660  *-------------------------------------------------  c$$$*-------------------------------------------------
661  *     this function returns  c$$$*     this function returns
662  *  c$$$*
663  *     - the Center-Of-Gravity of the cluster IC  c$$$*     - the Center-Of-Gravity of the cluster IC
664  *     evaluated using NCOG strips,  c$$$*     evaluated using NCOG strips,
665  *     calculated relative to MAXS(IC)  c$$$*     calculated relative to MAXS(IC)
666  *      c$$$*    
667  *     - zero in case that not  enough strips  c$$$*     - zero in case that not  enough strips
668  *     have a positive signal  c$$$*     have a positive signal
669  *      c$$$*    
670  *     NOTE:  c$$$*     NOTE:
671  *     This is the old definition, used by Straulino.  c$$$*     This is the old definition, used by Straulino.
672  *     The new routine, according to Landi,  c$$$*     The new routine, according to Landi,
673  *     is COG(NCOG,IC)  c$$$*     is COG(NCOG,IC)
674  *-------------------------------------------------  c$$$*-------------------------------------------------
675    c$$$
676    c$$$
677        include 'commontracker.f'  c$$$      include 'commontracker.f'
678        include 'level1.f'  c$$$      include 'level1.f'
679          c$$$      
680  *     --> signal of the central strip  c$$$*     --> signal of the central strip
681        sc = CLSIGNAL(INDMAX(ic)) !center  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center
682    c$$$
683  *     signal of adjacent strips  c$$$*     signal of adjacent strips
684  *     --> left  c$$$*     --> left
685        sl1 = 0                  !left 1  c$$$      sl1 = 0                  !left 1
686        if(  c$$$      if(
687       $     (INDMAX(ic)-1).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)
688       $     )  c$$$     $     )
689       $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
690    c$$$
691        sl2 = 0                  !left 2  c$$$      sl2 = 0                  !left 2
692        if(  c$$$      if(
693       $     (INDMAX(ic)-2).ge.INDSTART(ic)  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)
694       $     )  c$$$     $     )
695       $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
696    c$$$
697  *     --> right  c$$$*     --> right
698        sr1 = 0                  !right 1  c$$$      sr1 = 0                  !right 1
699        if(  c$$$      if(
700       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
701       $     .or.  c$$$     $     .or.
702       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
703       $     )  c$$$     $     )
704       $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
705    c$$$
706        sr2 = 0                  !right 2  c$$$      sr2 = 0                  !right 2
707        if(  c$$$      if(
708       $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
709       $     .or.  c$$$     $     .or.
710       $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
711       $     )  c$$$     $     )
712       $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
713          c$$$      
714  ************************************************************  c$$$************************************************************
715  *     COG computation  c$$$*     COG computation
716  ************************************************************  c$$$************************************************************
717    c$$$
718  c      print*,sl2,sl1,sc,sr1,sr2  c$$$c      print*,sl2,sl1,sc,sr1,sr2
719    c$$$
720        COG = 0.  c$$$      COG = 0.
721            c$$$        
722        if(sl1.gt.sr1.and.sl1.gt.0.)then  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then
723            c$$$        
724           if(ncog.eq.2.and.sl1.ne.0)then  c$$$         if(ncog.eq.2.and.sl1.ne.0)then
725              COG = -sl1/(sl1+sc)          c$$$            COG = -sl1/(sl1+sc)        
726           elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
727              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
728           elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
729              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)  c$$$            COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
730           else  c$$$         else
731              COG = 0.  c$$$            COG = 0.
732           endif  c$$$         endif
733            c$$$        
734        elseif(sl1.le.sr1.and.sr1.gt.0.)then  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then
735    c$$$
736           if(ncog.eq.2.and.sr1.ne.0)then  c$$$         if(ncog.eq.2.and.sr1.ne.0)then
737              COG = sr1/(sc+sr1)              c$$$            COG = sr1/(sc+sr1)            
738           elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
739              COG = (sr1-sl1)/(sl1+sc+sr1)  c$$$            COG = (sr1-sl1)/(sl1+sc+sr1)
740           elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
741              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)  c$$$            COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)
742           else  c$$$         else
743              COG = 0.  c$$$            COG = 0.
744           endif  c$$$         endif
745    c$$$
746        endif  c$$$      endif
747    c$$$
748        COG0 = COG  c$$$      COG0 = COG
749    c$$$
750  c      print *,ncog,ic,cog,'/////////////'  c$$$c      print *,ncog,ic,cog,'/////////////'
751    c$$$
752        return  c$$$      return
753        end  c$$$      end
754    
755  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
756        real function cog(ncog,ic)        real function cog(ncog,ic)
# Line 734  c      print *,ncog,ic,cog,'//////////// Line 782  c      print *,ncog,ic,cog,'////////////
782  *     --> signal of the central strip  *     --> signal of the central strip
783           sc = CLSIGNAL(INDMAX(ic)) !center           sc = CLSIGNAL(INDMAX(ic)) !center
784  *     signal of adjacent strips  *     signal of adjacent strips
785           sl1 = 0                !left 1           sl1 = -9999.           !left 1
786           if(           if(
787       $        (INDMAX(ic)-1).ge.INDSTART(ic)       $        (INDMAX(ic)-1).ge.INDSTART(ic)
788       $        )       $        )
789       $        sl1 = CLSIGNAL(INDMAX(ic)-1)       $        sl1 = CLSIGNAL(INDMAX(ic)-1)
790                    
791           sl2 = 0                !left 2           sl2 = -9999.           !left 2
792           if(           if(
793       $        (INDMAX(ic)-2).ge.INDSTART(ic)       $        (INDMAX(ic)-2).ge.INDSTART(ic)
794       $        )       $        )
795       $        sl2 = CLSIGNAL(INDMAX(ic)-2)       $        sl2 = CLSIGNAL(INDMAX(ic)-2)
796                    
797           sr1 = 0                !right 1           sr1 = -9999.           !right 1
798           if(           if(
799       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
800       $        .or.       $        .or.
# Line 754  c      print *,ncog,ic,cog,'//////////// Line 802  c      print *,ncog,ic,cog,'////////////
802       $        )       $        )
803       $        sr1 = CLSIGNAL(INDMAX(ic)+1)       $        sr1 = CLSIGNAL(INDMAX(ic)+1)
804                    
805           sr2 = 0                !right 2           sr2 = -9999.           !right 2
806           if(           if(
807       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
808       $        .or.       $        .or.
# Line 766  c      print *,ncog,ic,cog,'//////////// Line 814  c      print *,ncog,ic,cog,'////////////
814                    
815  c         print*,'## ',sl2,sl1,sc,sr1,sr2  c         print*,'## ',sl2,sl1,sc,sr1,sr2
816    
817    c     ==============================================================
818           if(ncog.eq.1)then           if(ncog.eq.1)then
819              COG = 0.              COG = 0.
820                if(sr1.gt.sc)cog=1. !NEW
821                if(sl1.gt.sc.and.sl1.gt.sr1)cog=-1. !NEW
822    c     ==============================================================
823           elseif(ncog.eq.2)then           elseif(ncog.eq.2)then
824              if(sl1.gt.sr1)then              if(sl1.gt.sr1)then
825                 if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)                         if((sl1+sc).ne.0)COG = -sl1/(sl1+sc)        
826              elseif(sl1.le.sr1)then              elseif(sl1.lt.sr1)then
827                 if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                             if((sc+sr1).ne.0)COG = sr1/(sc+sr1)                        
828              endif              elseif( sl1.eq.sr1.and.sl1.ne.-9999.)then !NEW
829                   if( clsigma(indmax(ic)-1).lt.clsigma(indmax(ic)+1)
830         $              .and.(sl1+sc).ne.0 )cog = -sl1/(sl1+sc) !NEW
831                   if( clsigma(indmax(ic)-1).gt.clsigma(indmax(ic)+1)
832         $              .and.(sc+sr1).ne.0 )cog = sr1/(sc+sr1) !NEW
833                endif
834    c            if(cog==0)print*,'Strange cluster (2) - @maxs ',MAXS(ic)
835    c     $           ,' : ',sl2,sl1,sc,sr1,sr2
836    c     ==============================================================
837           elseif(ncog.eq.3)then           elseif(ncog.eq.3)then
838               if((sl1+sc+sr1).ne.0)COG = (sr1-sl1)/(sl1+sc+sr1)               if( (sl1+sc+sr1).ne.0 )COG = (sr1-sl1)/(sl1+sc+sr1)
839    c             if(cog==0)print*,'Strange cluster (3) - @maxs ',MAXS(ic)
840    c     $            ,' : ',sl2,sl1,sc,sr1,sr2
841    c     ==============================================================
842           elseif(ncog.eq.4)then           elseif(ncog.eq.4)then
843              if(sl2.gt.sr2)then              if(sl2.gt.sr2)then
844                 if((sl2+sl1+sc+sr1).ne.0)                 if((sl2+sl1+sc+sr1).ne.0)
845       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)       $              COG = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1)
846              elseif(sl2.le.sr2)then              elseif(sl2.lt.sr2)then
847                  if((sl2+sl1+sc+sr1).ne.0)                 if((sr2+sl1+sc+sr1).ne.0)
848       $              COG = (2*sr2+sr1-sl1)/(sl2+sl1+sc+sr1)       $              COG = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)
849                elseif(sl2.eq.sr2.and.sl2.ne.-9999.)then !NEW
850                   if( clsigma(indmax(ic)-2).lt.clsigma(indmax(ic)+2)
851         $              .and.(sl2+sl1+sc+sr1).ne.0 )
852         $              cog = (sr1-sl1-2*sl2)/(sl2+sl1+sc+sr1) !NEW
853                   if( clsigma(indmax(ic)-2).gt.clsigma(indmax(ic)+2)
854         $              .and.(sr2+sl1+sc+sr1).ne.0 )
855         $              cog = (2*sr2+sr1-sl1)/(sr2+sl1+sc+sr1)  !NEW              
856              endif              endif
857           else           else
858              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG              print*,'function COG(NCOG,IC) ==> WARNING!! NCOG=',NCOG
# Line 814  c         print*,'-------' Line 884  c         print*,'-------'
884                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
885                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
886                 mu = mu + 1                 mu = mu + 1
887                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
888              else              else
889                 goto 10                 goto 10
890              endif              endif
# Line 827  c         print*,'-------' Line 897  c         print*,'-------'
897                 COG = COG + ipos*CLSIGNAL(i)                 COG = COG + ipos*CLSIGNAL(i)
898                 SGN = SGN + CLSIGNAL(i)                 SGN = SGN + CLSIGNAL(i)
899                 mu = mu + 1                 mu = mu + 1
900                 print*,ipos,CLSIGNAL(i)  c               print*,ipos,CLSIGNAL(i)
901              else              else
902                 goto 20                 goto 20
903              endif              endif
904           enddo           enddo
905   20      continue   20      continue
906           if(SGN.le.0)then           if(SGN.le.0)then
907  c            print*,'cog(0,ic) --> ic, dedx ',ic,SGN              print*,'cog(0,ic) --> ic, dedx ',ic,SGN
908              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)              print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
909              print*,(CLSIGNAL(i),i=istart,istop)              print*,(CLSIGNAL(i),i=istart,istop)
910  c            print*,'cog(0,ic) --> NOT EVALUATED '  c            print*,'cog(0,ic) --> NOT EVALUATED '
# Line 1189  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1259  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1259    
1260  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1261    
1262        FUNCTION risx_eta2(x)        FUNCTION risxeta2(x)
1263    
1264        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1265        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
# Line 1276  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1346  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1346     20 CONTINUE     20 CONTINUE
1347        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1348                
1349        risx_eta2=HQUADF* 1e-4        risxeta2=HQUADF* 1e-4
1350    
1351        END        END
1352    
1353  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1354        FUNCTION risx_eta3(x)        FUNCTION risxeta3(x)
1355        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1356        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1357        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1367  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1437  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1437     20 CONTINUE     20 CONTINUE
1438        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1439                
1440        risx_eta3 = HQUADF* 1e-4        risxeta3 = HQUADF* 1e-4
1441    
1442        END        END
1443  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1444        FUNCTION risx_eta4(x)        FUNCTION risxeta4(x)
1445        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1446        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1447        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1457  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1527  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1527     20 CONTINUE     20 CONTINUE
1528        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1529                
1530        risx_eta4=HQUADF* 1e-4        risxeta4=HQUADF* 1e-4
1531    
1532        END        END
1533  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1534        FUNCTION risy_eta2(x)        FUNCTION risyeta2(x)
1535        DOUBLE PRECISION V( 1)        DOUBLE PRECISION V( 1)
1536        INTEGER NPAR, NDIM, IMQFUN, I, J        INTEGER NPAR, NDIM, IMQFUN, I, J
1537        DOUBLE PRECISION HQDJ, VV, VCONST        DOUBLE PRECISION HQDJ, VV, VCONST
# Line 1529  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1599  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1599     20 CONTINUE     20 CONTINUE
1600        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)        IF (IMQFUN .EQ. 2) HQUADF = VCONST * EXP (HQUADF)
1601    
1602        risy_eta2=HQUADF* 1e-4        risyeta2=HQUADF* 1e-4
1603    
1604        END        END
1605  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
# Line 1680  c         if(BAD(VIEW(ic),nvk(MAXS(ic)), Line 1750  c         if(BAD(VIEW(ic),nvk(MAXS(ic)),
1750        risx_cog = HQUADF * 1e-4        risx_cog = HQUADF * 1e-4
1751    
1752        END        END
1753    
1754    
1755    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1756          real function pfacorr(ic,angle) !(1)
1757    *--------------------------------------------------------------
1758    *     this function returns the landi correction for this cluster
1759    *--------------------------------------------------------------
1760          include 'commontracker.f'
1761          include 'calib.f'
1762          include 'level1.f'
1763    
1764          real angle
1765          integer iview,lad
1766    
1767          iview = VIEW(ic)            
1768          lad = nld(MAXS(ic),VIEW(ic))
1769    
1770    *     find angular bin
1771    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1772          do iang=1,nangbin
1773             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1774                iangle=iang
1775                goto 98
1776             endif
1777          enddo
1778          if(DEBUG.eq.1)
1779         $     print*,'pfacorr *** warning *** angle out of range: ',angle
1780          if(angle.lt.angL(1))iang=1
1781          if(angle.gt.angR(nangbin))iang=nangbin
1782     98   continue                  !jump here if ok
1783    
1784          pfacorr = fcorr(iview,lad,iang)
1785    
1786          if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr
1787    
1788    
1789     100  return
1790          end
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.23