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