--- DarthVader/TrackerLevel2/src/F77/functions.f	2006/05/19 13:15:55	1.1.1.1
+++ DarthVader/TrackerLevel2/src/F77/functions.f	2007/11/27 15:28:58	1.6
@@ -1,4 +1,4 @@
-*************************************************************************
+************************************************************************
 *     
 *     functions.f
 *     
@@ -132,169 +132,8 @@
 
       end
 
-c$$$      
-c$$$c$$$c------------------------------------------------------------------------
-c$$$c$$$
-c$$$c$$$c     NB: le coordinate in mech_pos.dat sono calcolate a partire da alcuni dati
-c$$$c$$$c     contenuti in commontracker.f. forse si puo' evitare mech_pos.dat e mettere
-c$$$c$$$c     tutto in commontracker.f
-c$$$c$$$
-c$$$c$$$
-c$$$c$$$      subroutine mech_sensor    !it reads sensors coordinates (in PAMELA reference
-c$$$c$$$                                ! frame) from a text file and it uses them to fill
-c$$$c$$$                                ! x/y/z_mech_sensor variables, taking into account
-c$$$c$$$                                ! last plane inversion
-c$$$c$$$
-c$$$c$$$      include './commontracker.f'
-c$$$c$$$      include './common_tracks.f'
-c$$$c$$$
-c$$$c$$$      real xvec(nladders_view),yvec(2),zvec(nplanes)
-c$$$c$$$
-c$$$c$$$      integer id                !file identifier
-c$$$c$$$      logical od                !.true. if the specified unit is connected to a file
-c$$$c$$$
-c$$$c$$$      do id=20,100,1            !opens the file using a free file id
-c$$$c$$$        inquire (id, opened=od)
-c$$$c$$$        if(.not.od) goto 666
-c$$$c$$$      enddo      
-c$$$c$$$ 666  continue
-c$$$c$$$
-c$$$c$$$      open(id,FILE='../common/mech_pos.dat') !sensors centres coordinates in mm in 
-c$$$c$$$                                ! PAMELA reference frame:
-c$$$c$$$                                ! the first plane is the one with lowest Z (the one
-c$$$c$$$                                ! nearest the calorimeter)
-c$$$c$$$                                ! the first ladder is the one with lowest X (the
-c$$$c$$$                                ! one on which the first X strip is)
-c$$$c$$$                                ! the first sensor is the one with lowest Y (the
-c$$$c$$$                                ! one on which the first Y strip is) for planes
-c$$$c$$$                                ! 2..6. for plane 1 the first sensor has higher Y
-c$$$c$$$
-c$$$c$$$      read(20,*) xvec
-c$$$c$$$      read(20,*) yvec
-c$$$c$$$      read(20,*) zvec
-c$$$c$$$
-c$$$c$$$      do i=1,nplanes
-c$$$c$$$        do j=1,nladders_view
-c$$$c$$$          do k=1,2
-c$$$c$$$            x_mech_sensor(i,j,k)=xvec(j)
-c$$$c$$$            y_mech_sensor(i,j,k)=yvec(k)
-c$$$c$$$            z_mech_sensor(i,j,k)=zvec(i)
-c$$$c$$$            if(i.eq.1) then     !y coordinates of first plane (11th view) are
-c$$$c$$$              y_mech_sensor(i,j,k)=-yvec(k) ! exchanged due to last plane inversion
-c$$$c$$$            endif              
-c$$$c$$$          enddo
-c$$$c$$$        enddo
-c$$$c$$$      enddo
-c$$$c$$$
-c$$$c$$$      close(id)
-c$$$c$$$
-c$$$c$$$
-c$$$c$$$c$$$  ! *** INIZIO DEBUG ***
-c$$$c$$$c$$$  do i=1,6
-c$$$c$$$c$$$  do j=1,3
-c$$$c$$$c$$$  do k=1,2
-c$$$c$$$c$$$  c            print*,x_mech_sensor(1,j,k)
-c$$$c$$$c$$$  print*,y_mech_sensor(i,j,k)
-c$$$c$$$c$$$  c            print*,z_mech_sensor(i,j,k)
-c$$$c$$$c$$$  enddo
-c$$$c$$$c$$$  enddo
-c$$$c$$$c$$$  print*,' '
-c$$$c$$$c$$$  enddo
-c$$$c$$$c$$$  ! *** FINE DEBUG ***
-c$$$c$$$
-c$$$c$$$      
-c$$$c$$$      return
-c$$$c$$$      end
-
-
-c$$$c------------------------------------------------------------------------
-c$$$
-c$$$c     NB: le coordinate in mech_pos.dat sono calcolate a partire da alcuni dati
-c$$$c     contenuti in commontracker.f. forse si puo' evitare mech_pos.dat e mettere
-c$$$c     tutto in commontracker.f
-c$$$
-c$$$
-c$$$      subroutine mech_sensor    
-c$$$c     !it reads sensors coordinates (in PAMELA reference
-c$$$c     ! frame) from a text file and it uses them to fill
-c$$$c     ! x/y/z_mech_sensor variables, taking into account
-c$$$c     ! last plane inversion
-c$$$
-c$$$      include './commontracker.f'
-c$$$      include './common_tracks.f'
-c$$$
-c$$$      real xvec(nladders_view),yvec(2),zvec(nplanes)
-c$$$
-c$$$      integer id                !file identifier
-c$$$      logical od                !.true. if the specified unit is connected to a file
-c$$$
-c$$$      do id=20,100,1            !opens the file using a free file id
-c$$$        inquire (id, opened=od)
-c$$$        if(.not.od) goto 666
-c$$$      enddo      
-c$$$ 666  continue
-c$$$
-c$$$c      open(id,FILE='../common/mech_pos.dat') !sensors centres coordinates in mm in 
-c$$$c      open(id,FILE='source/common/mech_pos.dat') 
-c$$$c      call system('cp $TRK_GRND/source/common/mech_pos.dat .')
-c$$$      print *,'Opening file: mech_pos.dat'
-c$$$      open(id,FILE='./bin-aux/mech_pos.dat',IOSTAT=iostat) 
-c$$$c     !sensors centres coordinates in mm in 
-c$$$c     ! PAMELA reference frame:
-c$$$c     ! the first plane is the one with lowest Z (the one
-c$$$c     ! nearest the calorimeter)
-c$$$c     ! the first ladder is the one with lowest X (the
-c$$$c     ! one on which the first X strip is)
-c$$$c     ! the first sensor is the one with lowest Y (the
-c$$$c     ! one on which the first Y strip is) for planes
-c$$$c     ! 2..6. for plane 1 the first sensor has higher Y
-c$$$      
-c$$$      if(iostat.ne.0)then
-c$$$         print*,'MECH_SENSOR: *** Error in opening file ***'
-c$$$         return
-c$$$      endif
-c$$$
-c$$$      read(id,*) xvec
-c$$$      read(id,*) yvec
-c$$$      read(id,*) zvec
-c$$$
-c$$$      do i=1,nplanes
-c$$$        do j=1,nladders_view
-c$$$          do k=1,2
-c$$$            x_mech_sensor(i,j,k)=xvec(j)
-c$$$            y_mech_sensor(i,j,k)=yvec(k)
-c$$$            z_mech_sensor(i,j,k)=zvec(i)
-c$$$            if(i.eq.1) then     !y coordinates of first plane (11th view) are
-c$$$              y_mech_sensor(i,j,k)=-yvec(k) ! exchanged due to last plane inversion
-c$$$            endif              
-c$$$          enddo
-c$$$        enddo
-c$$$      enddo
-c$$$
-c$$$      close(id)
-c$$$c      call system('rm -f mech_pos.dat')
-c$$$
-c$$$c$$$                                ! *** INIZIO DEBUG ***
-c$$$c$$$      do i=1,6
-c$$$c$$$c        do j=1,3
-c$$$c$$$          do k=1,2
-c$$$c$$$            j=1
-c$$$c$$$c            print*,x_mech_sensor(1,j,k)
-c$$$c$$$            print*,y_mech_sensor(i,j,k)
-c$$$c$$$c            print*,z_mech_sensor(i,j,k)
-c$$$c$$$          enddo
-c$$$c$$$c        enddo
-c$$$c$$$        print*,' '
-c$$$c$$$      enddo
-c$$$c$$$                                ! *** FINE DEBUG ***
-c$$$
-c$$$      
-c$$$      return
-c$$$      end
-c$$$
-c$$$
-c$$$c------------------------------------------------------------------------
 
+c------------------------------------------------------------------------
 
       function coordsi(istrip,view) 
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
@@ -325,9 +164,9 @@
       
       if(mod(view,2).eq.0) then !X view
 
-        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
-          print*,'functions: WARNING: false X strip: strip ',is
-        endif
+c        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
+c          print*,'functions: WARNING: false X strip: strip ',is
+c        endif
 
         is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
 
@@ -396,9 +235,9 @@
       
       if(mod(view,2).eq.0) then !X view
 
-        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
-          print*,'functions: WARNING: false X strip: strip ',is
-        endif
+c        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
+c          print*,'functions: WARNING: false X strip: strip ',is
+c        endif
 
         is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
 
@@ -436,6 +275,85 @@
 c------------------------------------------------------------------------
 
 
+      double precision function dcoordsi(rstrip,view) 
+*
+*     same as COORDSI, but accept a real value of strip!!!
+*     and gives a double precision output
+*
+* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+*     it gives the strip coordinate in micrometers,
+*     knowing the strip number (1..3072) and the view
+*     number. the origin of the coordinate is on the
+*     centre of the sensor the strip belongs to.
+*     the axes directions are the same as in the PAMELA
+*     reference frame (i.e.: the 11th view coordinate
+*     direction has to be inverted here)
+* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+c      integer is,view,istrip
+
+      integer view,is,istrip
+      real rstrip
+      double precision strip,stripladder,p
+      double precision edge,dim
+      double precision coord1
+
+
+      include 'commontracker.f'
+
+c     NB mettere il 1024 nel commontracker...!???
+
+      strip=DBLE(rstrip)
+
+      istrip = int(strip+0.5)   !istrip stores the closest integer to strip      
+
+      is=istrip                 !it stores istrip number
+      is=mod(is-1,1024)+1       !it puts all clusters on a single ladder
+      
+      dcoordsi=0.
+      
+      if(mod(view,2).eq.0) then !X view
+
+c        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
+c          print*,'functions: WARNING: false X strip: strip ',is
+c        endif
+
+        is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
+
+        edge=edgeX
+        dim=SiDimX
+
+      elseif(mod(view,2).eq.1) then !Y view
+
+        edge=edgeY
+        dim=SiDimY
+
+c$$$        if(view.eq.11) then     !INVERSIONE!???
+c$$$          is=1025-is
+c$$$        endif
+
+      endif
+
+
+      stripladder = DBLE(is)+(strip-DBLE(istrip))!cluster position relative to ladder
+      p=pitch(view)
+
+ccccc coord1=(is-1)*p           !referred to 1st sensor strip
+      coord1=(stripladder-1)*p  !referred to 1st sensor strip
+      coord1=coord1+edge        !referred to sensor edge
+      dcoordsi=coord1-dim/2     !referred to the centre of the sensor
+
+      if(view.eq.11) then       !INVERSION: it puts y axis in the same direction for all views
+        dcoordsi=-dcoordsi
+      endif
+      
+      end
+
+
+
+c------------------------------------------------------------------------
+
+
       function coord(coordsi,view,ladder,sen) 
 *     it gives the coordinate in 
 *     micrometers, knowing the coordinate in the sensor
@@ -520,10 +438,25 @@
 
       dcoord=0.
       
+      if(
+     $     ladder.lt.1.or.
+     $     ladder.gt.nladders_view.or.
+     $     sen.lt.1.or.
+     $     sen.gt.2.or.
+     $     view.lt.1.or.
+     $     view.gt.nviews.or.
+     $     .false.)then
+         print*,'dcoord ---> wrong input: ladder ',ladder
+     $        ,' sensor ',sen
+     $        ,' view ',view
+         return
+      endif
+
       sx=ladder
       sy=sen
       sz=npl(view)
 
+
       if(mod(view,2).eq.0) then !X view
 
         trasl=x_mech_sensor(sz,sx,sy) !in mm
@@ -544,6 +477,185 @@
 
 
 c------------------------------------------------------------------------
+      integer function nsatstrips(ic)
+*--------------------------------------------------------------
+*     this function returns the number of saturated strips 
+*     inside a cluster
+*--------------------------------------------------------------
+      include 'commontracker.f'
+      include 'level1.f'
+      include 'calib.f'
 
+      
 
+      integer nsat 
+      nsat = 0 
+      iv=VIEW(ic)               
+      if(mod(iv,2).eq.1)incut=incuty
+      if(mod(iv,2).eq.0)incut=incutx
+      istart = INDSTART(IC)
+      istop  = TOTCLLENGTH
+      if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
+      do i = INDMAX(IC),istart,-1
+c         cut  = incut*CLSIGMA(i)
+c         if(CLSIGNAL(i).ge.cut)then
+         if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
+     $        .or.
+     $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
+            nsat = nsat +1 
+         else
+            goto 10
+         endif
+      enddo
+ 10   continue
+      do i = INDMAX(IC)+1,istop
+c         cut  = incut*CLSIGMA(i)
+c         if(CLSIGNAL(i).ge.cut)then
+         if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
+     $        .or. 
+     $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
+            nsat = nsat +1 
+         else
+            goto 20
+         endif
+      enddo
+ 20   continue
+            
+      nsatstrips = nsat
+      return 
+      end
 
+c------------------------------------------------------------------------
+      integer function nbadstrips(ncog,ic)
+*--------------------------------------------------------------
+*     this function returns the number of BAD strips 
+*     inside a cluster:
+*     - if NCOG=0, the number BAD strips inside the whole cluster
+*     are given, according to the cluster multiplicity
+*     
+*     - if NCOG>0, the number BAD strips is evaluated using NCOG 
+*     strips, even if they have a negative signal (according to Landi)
+*--------------------------------------------------------------
+      include 'commontracker.f'
+      include 'level1.f'
+      include 'calib.f'
+
+      integer nbad 
+      nbad = 0 
+
+      if (ncog.gt.0) then
+
+*     --> signal of the central strip
+         sc = CLSIGNAL(INDMAX(ic)) !center
+*     signal of adjacent strips
+         sl1 = -100000                !left 1
+         if(
+     $        (INDMAX(ic)-1).ge.INDSTART(ic)
+     $        )
+     $        sl1 = CLSIGNAL(INDMAX(ic)-1)
+         
+         sl2 = -100000                !left 2
+         if(
+     $        (INDMAX(ic)-2).ge.INDSTART(ic)
+     $        )
+     $        sl2 = CLSIGNAL(INDMAX(ic)-2)
+         
+         sr1 = -100000                !right 1
+         if(
+     $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
+     $        .or.
+     $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
+     $        )
+     $        sr1 = CLSIGNAL(INDMAX(ic)+1)
+         
+         sr2 = -100000                !right 2
+         if(
+     $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
+     $        .or.
+     $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
+     $        )
+     $        sr2 = CLSIGNAL(INDMAX(ic)+2)
+         
+         if(ncog.ge.1)nbad = nbad+1-CLBAD(INDMAX(ic))
+         
+         if(ncog.ge.2)then
+            if(sl1.gt.sr1.and.(INDMAX(ic)-1).ge.1)then
+               nbad=nbad+1-CLBAD(INDMAX(ic)-1)      
+            elseif(sl1.le.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
+               nbad=nbad+1-CLBAD(INDMAX(ic)+1)           
+            endif
+         endif
+         
+         if(ncog.ge.3)then
+            if(sl1.gt.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
+               nbad=nbad+1-CLBAD(INDMAX(ic)+1)
+            elseif(sl1.le.sr1.and.(INDMAX(ic)-1).ge.1)then
+c               if(INDMAX(ic)-1.eq.0)
+c     $              print*,' ======= ',sl2,sl1,sc,sr1,sr2
+               nbad=nbad+1-CLBAD(INDMAX(ic)-1)
+            endif
+         endif 
+         
+         if(ncog.ge.4)then
+            if(sl2.gt.sr2.and.(INDMAX(ic)-2).ge.1)then
+               nbad=nbad+1-CLBAD(INDMAX(ic)-2)      
+            elseif(sl2.le.sr2.and.(INDMAX(ic)+2).le.NCLSTR1)then
+               nbad=nbad+1-CLBAD(INDMAX(ic)+2)           
+            endif
+         endif
+         
+c         if(ncog.ge.5)then
+c            print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
+c     $           ,' not implemented'
+c         endif
+            
+      elseif(ncog.eq.0)then
+*     =========================
+*     COG computation
+*     =========================
+
+
+         
+         iv=VIEW(ic)         
+         if(mod(iv,2).eq.1)incut=incuty
+         if(mod(iv,2).eq.0)incut=incutx
+         
+         istart = INDSTART(IC)
+         istop  = TOTCLLENGTH
+         if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
+         nbad = 0
+c$$$         do i=istart,istop
+c$$$            cut  = incut*CLSIGMA(i)            
+c$$$            if(CLSIGNAL(i).ge.cut)nbad = nbad +1 -CLBAD(i)
+c$$$         enddo
+         do i = INDMAX(IC),istart,-1
+            cut  = incut*CLSIGMA(i)
+            if(CLSIGNAL(i).ge.cut)then
+               nbad = nbad +1 -CLBAD(i)
+            else
+               goto 10
+            endif
+         enddo
+ 10      continue
+         do i = INDMAX(IC)+1,istop
+            cut  = incut*CLSIGMA(i)
+            if(CLSIGNAL(i).ge.cut)then
+               nbad = nbad +1 -CLBAD(i)
+            else
+               goto 20
+            endif
+         enddo
+ 20      continue
+         
+      else
+         
+c         print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
+c     $        ,' not implemented'
+         
+
+      endif
+
+      nbadstrips = nbad 
+
+      return
+      end