*
* $Id: gpdspe.F,v 3.10 2007/05/08 15:05:18 bottai Exp $
*
* $Log: gpdspe.F,v $
* Revision 3.10  2007/05/08 15:05:18  bottai
*  avoid simulating Lorentz angle drift in SPE spread positions
*  since new reconstruction correct it at strip level
*
* Revision 3.9  2006/05/02 12:00:24  bottai
*  call gpucollection for strip simulation, condition on FFIELD repositioned
*
* Revision 3.8  2006/04/09 23:29:19  cafagna
* Several new things, among this: ND and CARD
*
* Revision 3.7  2006/02/07 11:38:52  bottai
*  test of mail delivery
*
* Revision 3.6  2006/02/06 14:17:21  bottai
*  adapt again the resolution simulation to the ITRPB volume identification
*
* Revision 3.5  2006/01/24 13:52:01  bottai
* prepare to store the x and y multiplicity of the hit
*
* Revision 3.4  2006/01/18 09:41:33  bottai
*  Some correction to the spe plane ID due to the new geometry tree
*
* Revision 3.3  2005/12/05 12:15:21  pam-ba
* new spectrometer geometry and internal magnetic field
*
* Revision 3.2  2005/07/25 11:53:21  cafagna
* Several updates. See history for details
*
* Revision 3.1.1.1  2002/07/11 16:02:13  cafagna
* First GPAMELA release on CVS
*
*
*CMZ :  2.01/00 05/04/2000  14.35.18  by  Marialuigia Ambriola
*CMZ :  2.00/00 03/03/2000  15.39.06  by  Francesco Cafagna
*CMZ :  1.02/00 19/01/2000  15.50.14  by  Francesco Cafagna
*CMZ :  1.01/00 16/05/96  17.38.13  by  Francesco Cafagna
*CMZ :  1.00/02 27/03/96  18.35.42  by  Francesco Cafagna
*-- Author :    Francesco Cafagna   23/02/96
      SUBROUTINE GPDSPE
*************************************************************************
*                                                                       *
* To perform SPE digitization                                           *
*                                                                       *
* Called by: GUDIGI                                                     *
* Author: Francesco Cafagna, 23/02/96 12.50.38                          *
*                                                                       *
*************************************************************************
#include "gcunit.inc"
#include "gpphit.inc"
      PARAMETER (RADDEG=57.295779)
      PARAMETER (NH=NHSPE,NHM=NHMSPE,NVOL=3)
#include "gppdigi.inc"
c +SEQ,GPDSPE.
#include "gprhit.inc"
#include "gpspe.inc"
#include "gpsed.inc"
#include "gpkey.inc"
*et. 04/11/2005
#include "gpstripspe.inc"
*end et.
*
      INTEGER IS,ID,NV,NVM,GPNHIT,NTHIT,NHFRST
      CHARACTER*4 CSET
      LOGICAL CONT,HF
*
* Data block for digitalization
*
      INTEGER NTRA,IDIG,LTRA,NPLANE,LKDN,NTDIGI,JJ
      INTEGER KDIGI(4),NUMBVD(NVOL)
      REAL XI,YI,XO,YO,PATH,DE,GPGAUS
      REAL VPOS(3),BMAGNET(3),ANGLEX,ANGLEY,XMULT,YMULT,DXPOS,DYPOS
*
      DATA NTRA,LTRA/1,1/

*et. 04/11/2005
*
*
* Vector zeroed
*
      CALL VZERO(ISTRIPX,MAXSTR)            
      CALL VZERO(NPSTRIPX,MAXSTR)            
      CALL VZERO(NTSTRIPX,MAXSTR)
      CALL VZERO(QSTRIPX,MAXSTR)
      CALL VZERO(XSTRIPX,MAXSTR)
      CALL VZERO(ISTRIPY,MAXSTR)            
      CALL VZERO(NPSTRIPY,MAXSTR)            
      CALL VZERO(NTSTRIPY,MAXSTR)
      CALL VZERO(QSTRIPY,MAXSTR)
      CALL VZERO(YSTRIPY,MAXSTR)
      NSTRPX=0
      NSTRPY=0
      
* 

      do ip=1,nplanes
	 do iip=1,nsensors
	     isignyesx=0
	     isignyesy=0
	     do iiip=1,nstripx
		if(proxtanti(ip,iip,iiip).gt.0) then
		    isignyesx=1
		    GO TO 456
    		endif
	     enddo
	     do iiip=1,nstripy
	         if(proytanti(ip,iip,iiip).gt.0) then
		      isignyesy=1
		      GO TO 456
	         endif
	     enddo
456          if(isignyesx.eq.1.or.isignyesy.eq.1) then
		 call gpucollection(ip,iip)
	     endif
         enddo
      enddo




*
*
* Vector zeroed
*
      CALL VZERO(ITRPB ,NHM)
      CALL VZERO(ITRSL ,NHM)
      CALL VZERO(ITSPA ,NHM)
      CALL VZERO(NXMULT ,NHM)
      CALL VZERO(NYMULT ,NHM)      
      CALL VZERO(XTSPA,NHM)
      CALL VZERO(YTSPA,NHM)
      CALL VZERO(ZTSPA,NHM)
      CALL VZERO(XOTSPA,NHM)
      CALL VZERO(YOTSPA,NHM)
      CALL VZERO(ZOTSPA,NHM)
      CALL VZERO(ERTSPA,NHM)
      CALL VZERO(PATSPA,NHM)
      CALL VZERO(NUMD  ,NVOL)
      CALL VZERO(NDPLAN,NDMSPE)
      CALL VZERO(NDSIL ,NDMSPE)
      CALL VZERO(NDISTR,NDMSPE)
      CALL VZERO(NDIADC,NDMSPE)
*
      CONT = .FALSE.
      HF   = .TRUE.
      NTDIGI = 0
      NDISPE = 0
      NV = NVOL
      NVM= NVMAX
      IS   = ISSPE
      ID   = IDTSPA
      CSET = 'SPE '
      IT = -1
*
* Get the total number of hits stored for TSPA detector
*
      NTHIT = GPNHIT(CSET,CNSED(IS,ID) )
      IF (NTHIT.GT.NHM) CONT=.TRUE.
*
* Check on the volume level
*
      IF (NV.GT.NVM) THEN
         WRITE(CHMAIL,10000) CSET,CNSED(IS,ID),NVOL,NVMAX,
     +                        IDRUN,IDEVT
         CALL GMAIL(1,0)
         GO TO 20
      ENDIF
      CALL VZERO(VHITS,NH*NHM)
      CALL VZERO(ITRA ,NHM)
      CALL VZERO(NUMBV,NVOL*NHM)
      NHITS = -1
      NHFRST = 1
      LKDN = 0
   10 CALL GPRHIT(CSET,CNSED(IS,ID),IS,ID,NVOL,NH,NHM,NHFRST,
     +            IT,NUMD,ITRA,NUMBV,VHITS,NHITS)
*
* If there was a problem in GPRHIT then exit
*
      IF(NHITS.LT.0.AND.(NHITS.NE.-(NHM+1)) ) THEN
         WRITE(CHMAIL,10100) CSET,CNSED(IS,ID),NHITS,NHM,
     +                        IDRUN,IDEVT
         CALL GMAIL(1,0)
         GO TO 20
      ENDIF
      IF( NHITS.EQ. -(NHM+1) ) NHITS = NHM
*
* Fill CWN entuple
*
      NTHSPE = NHITS
      DO I=1,NHITS
         NSPEPLANE=INT((NUMBV(2,I)-1)/6)+1
         NUMSENS=NUMBV(2,I)-(NSPEPLANE-1)*6  
         ITRPB(I)   = NSPEPLANE
         ITRSL(I)   = NUMSENS
         ITSPA(I)   = NUMBV(3,I)
#if defined(GPAMELA_INTER)
         XAVSPE(I) = VHITS(1,I)
         YAVSPE(I) = VHITS(2,I)
         ZAVSPE(I) = VHITS(3,I)
         XOUTSPE(I) = VHITS(4,I)
         YOUTSPE(I) = VHITS(5,I)
         ZOUTSPE(I) = VHITS(6,I)
         XINSPE(I) = 2*XAVSPE(I) - XOUTSPE(I)
         YINSPE(I) = 2*YAVSPE(I) - YOUTSPE(I)
         ZINSPE(I) = 2*ZAVSPE(I) - ZOUTSPE(I)
#endif
#if !defined(GPAMELA_INTER)
         XINSPE(I) = VHITS(1,I)
         YINSPE(I) = VHITS(2,I)
         ZINSPE(I) = VHITS(3,I)
         XOUTSPE(I) = VHITS(4,I)
         YOUTSPE(I) = VHITS(5,I)
         ZOUTSPE(I) = VHITS(6,I)
         XAVSPE(I) = (XINSPE(I)+XOUTSPE(I))/2.
         YAVSPE(I) = (YINSPE(I)+YOUTSPE(I))/2.
         ZAVSPE(I) = (ZINSPE(I)+ZOUTSPE(I))/2.
#endif
         ERELSPE(I) = VHITS(7,I)
         PATHSPE(I) = VHITS(8,I)
         IPARSPE(I) = VHITS(9,I)
         P0SPE (I) = VHITS(10,I)

C
C ADD RESOLUTION SIMULATION TO XAVSPE 
C	 
         ANGLEX=0.
         ANGLEY=0.
         IF(ZOUTSPE(I).NE.ZINSPE(I)) THEN 
            ANGLEX=RADDEG*ATAN( ABS(XOUTSPE(I)-XINSPE(I) )/
     +           ABS(ZOUTSPE(I)-ZINSPE(I) ) ) 
            ANGLEY=RADDEG*ATAN( ABS(YOUTSPE(I)-YINSPE(I) )/
     +           ABS(ZOUTSPE(I)-ZINSPE(I) ) )           
         ENDIF
            
         VPOS(1)=XAVSPE(I)
         VPOS(2)=YAVSPE(I)
         VPOS(3)=ZAVSPE(I)	 
         CALL GUFLD(VPOS,BMAGNET)
            
         CALL GPRSPE(ANGLEX,ANGLEY,BMAGNET,XMULT,YMULT,DXPOS,DYPOS)

         NXMULT(I)=INT(XMULT)
         NYMULT(I)=INT(YMULT)
            
         XAVSPE(I) = XAVSPE(I) + DXPOS
         YAVSPE(I) = YAVSPE(I) + DYPOS

C
C ADD THE SHIFT INDUCED BY THE MAG. FIELD
C
cv
cv 8/05/2007 not anymore needed since the effect is corrected in reco
cv before calculating the positions for fit
cv
cv         IF(FFIELD.NE.0) THEN
cv	 
cv            IF(ITRPB(I).LT.6) THEN
cv               XAVSPE(I) = XAVSPE(I) + 2.64e-4/4.5*BMAGNET(2)
cv            ELSE IF(ITRPB(I).EQ.6) THEN
cv               XAVSPE(I) = XAVSPE(I) - 2.64e-4/4.5*BMAGNET(2)
cv            ENDIF
cv         ENDIF   
      ENDDO

*
* Loop on the hits retrived to fill the CWN common and the DIGI b
*
C fc       IF (HF) NTHSPE = NHITS
C fc       DO I=1,NHITS
C fc *
C fc * Start DIGI action
C fc *
C fc          NPLANE=NUMBV(3,I)
C fc          XI=VHITS(1,I)
C fc          YI=VHITS(2,I)
C fc          ZI=VHITS(3,I)
C fc          XO=VHITS(4,I)
C fc          YO=VHITS(5,I)
C fc          ZO=VHITS(6,I)
C fc +SELF,IF=INTER
C fc * interactive case: we need the track edges
C fc          XI=2.*XI-XO
C fc          YI=2.*YI-YO
C fc          ZI=2.*ZI-ZO
C fc +SELF.
C fc          PATH=VHITS(8,I)
C fc          DE=VHITS(7,I)
C fc          IF(DE.GT.0.) THEN
C fc             DO J=1,6
C fc                NUMBVD(J)=NUMBV(J,I)
C fc             ENDDO
C fc c+SELF,IF=PAOLO
C fc c            CALL GPSPEHD(I,NPLANE,XI,YI,ZI,XO,YO,ZO,PATH,DE,LKDN)
C fc c+SELF
C fc c+SELF,IF=-PAOLO
C fc c            CALL GPSPEHD(I,NUMBVD,NVSPE,XI,YI,ZI,XO,YO,ZO,PATH,DE,LKDN)
C fc c+SELF.
C fc             NTDIGI = LKDN
C fc          ENDIF
C fc *
C fc * Fill CWN common
C fc *
C fc          IF(HF) THEN
C fc             ITRPB(I) = NUMBV((NVSPE-3),I)
C fc             ITRSL(I) = NUMBV((NVSPE-1),I)
C fc             ITSPA(I) = NUMBV(NVSPE,I)
C fc             XTSPA(I) = VHITS(1,I)
C fc             YTSPA(I) = VHITS(2,I)
C fc             ZTSPA(I) = VHITS(3,I)
C fc             XOTSPA(I) = VHITS(4,I)
C fc             YOTSPA(I) = VHITS(5,I)
C fc             ZOTSPA(I) = VHITS(6,I)
C fc             ERTSPA(I) = VHITS(7,I)
C fc             PATSPA(I) = VHITS(8,I)
C fc          ENDIF
C fc          IF(CONT) HF = .FALSE.
C fc       ENDDO
C fc       IF(CONT) THEN
C fc          NHFRST = NHFRST + NHM
C fc          IF( (NTHIT-NHFRST).LE.NHM) CONT=.FALSE.
C fc          GO TO 10
C fc       ENDIF
C fc *
C fc * fill digits with noise included
C fc *
C fc       IF(LKDN.GE.1) THEN
C fc          DO II=1,LKDN
C fc * plane=0 if the signal is lesser than 1 ADC
C fc             IF(INT(RLKD4(II)*ADCX).LE.0.AND.LKD2(II).GE.1.
C fc      +      AND.LKD2(II).LE.3) LKD1(II)=0
C fc             IF(INT(RLKD4(II)*ADCY).LE.0.AND.LKD2(II).GE.4.
C fc      +      AND.LKD2(II).LE.5) LKD1(II)=0
C fc *
C fc             JJ=LKDH(II)+(NHFRST-1)
C fc *            KDIGI(1)=LKD1(II)
C fc *            KDIGI(2)=LKD2(II)
C fc             KDIGI(1)=LKD2(II)
C fc             KDIGI(2)=LKD3(II)
C fc *            R1=RAN(IR)
C fc *            R2=RAN(IR)
C fc *            GRNDM0=SIN(2.*3.14159*R1)*SQRT(-2.*LOG(R2))
C fc             GRNDM0 = GPGAUS(0.)
C fc             IF(LKD2(II).GE.1.AND.LKD2(II).LE.3) THEN
C fc                RLKD4(II)=RLKD4(II)+GRNDM0*ENOISX+PIEDX
C fc                KDIGI(3)=INT(RLKD4(II)*ADCX)
C fc             ENDIF
C fc             IF(LKD2(II).GE.4.AND.LKD2(II).LE.5) THEN
C fc                RLKD4(II)=RLKD4(II)+GRNDM0*ENOISY+PIEDY
C fc                KDIGI(3)=INT(RLKD4(II)*ADCY)
C fc             ENDIF
C fc C+SELF,IF=PAOLO
C fc             DO J=1,NVOL
C fc                NUMBVD(J)=NUMBV(J,JJ)
C fc             ENDDO
C fc C+SELF
C fc             LTRA=ITRA(JJ)
C fc             IF(KDIGI(3).LT.0) KDIGI(3)=0
C fc             IF(KDIGI(3).GE.1024) KDIGI(3)=1023
C fc             IF(KDIGI(3).GT.0.AND.LKD1(II).GT.0) THEN
C fc                CALL GSDIGI(IS,ID,LTRA,NTRA,NUMBVD,KDIGI,IDIG)
C fc                IF (IDIG.EQ.0) THEN
C fc                   WRITE(CHMAIL,10200) CSET,CNSED(IS,ID),NVOL,KDIGI
C fc                   CALL GMAIL(1,0)
C fc                   GO TO 20
C fc                ENDIF
C fc             ENDIF
C fc *
C fc * Fills CWN
C fc *
C fc             IF(HF) THEN
C fc                IF(LKD2(II).GE.1.AND.LKD2(II).LE.3) KDIGIT=INT(RLKD4(II)*
C fc      +         ADCX)
C fc                IF(LKD2(II).GE.4.AND.LKD2(II).LE.5) KDIGIT=INT(RLKD4(II)*
C fc      +         ADCY)
C fc                IF(KDIGIT.GT.0.AND.LKD1(II).GT.0) THEN
C fc                   NDISPE=NDISPE+1
C fc                   IF (NDISPE.LE.NDMSPE) THEN
C fc                      NDPLAN(NDISPE)=LKD1(II)
C fc                      NDSIL (NDISPE)=LKD2(II)
C fc                      NDISTR(NDISPE)=LKD3(II)
C fc                      IF(KDIGIT.LT.1024) THEN
C fc                         NDIADC(NDISPE)=KDIGIT
C fc                      ELSE
C fc                         NDIADC(NDISPE)=1023
C fc                      ENDIF
C fc                   ENDIF
C fc                ENDIF
C fc             ENDIF
C fc          ENDDO
C fc       ENDIF
*
10000  FORMAT(' GPDSPE: VOLUME level too big: SET=',A4,', DET=',A4,',
     +        NVOL=',I3,', NVMAX=',I3,', RUN=',I5,', EVT=',I8)
10100  FORMAT(' GPDSPE: Error in GPRHIT: SET=',A4,
     +        ', DET=',A4,',NHITS=',I5,', NHMAX=',I5,
     +        ', RUN=',I5,', EVT=',I8)
10200  FORMAT(' GPDSPE: digitization ERROR: SET=',A4,', DET=',A4,',
     +        NVOL=',I3,', KDIGI=',4I4)
   20 RETURN
      END
