*
* $Id: ghesig.F,v 3.1 2006/11/03 11:19:54 bottai Exp $
*
* $Log: ghesig.F,v $
* Revision 3.1  2006/11/03 11:19:54  bottai
*  new antiproton cross sections
*
* Revision 1.1.1.1  1995/10/24 10:21:15  cernlib
* Geant
*
*
*CMZ :  3.21/02 29/03/94  15.41.37  by  S.Giani
*-- Author :
      FUNCTION GHESIG(PPART,EKIN,AVER,A,Z,W,KK,DENS,QCOR,LPART)
C
C *** CALCULATION OF THE PROBABILITIES FOR (IN)ELASTIC INTERACTIONS ***
C *** OF STABLE PARTICLES ON PROTON AND NEUTRON                     ***
C *** NVE 07-APR-1988 ***
C
C CALLED BY : GPGHEI
C ORIGIN : F.CARMINATI, H.FESEFELDT (ROUTINE INTACT 06-OCT-1987)
C
C OCTOBER 2006 - MODIFIED BY S.Bottai TO CHANGE ANTIPROTON CROSS SECTIONS 
C
C *** IPART DENOTES THE GHEISHA PARTICLE INDEX ***
C
C CONVENTION :
C
C   PARTICLE                 IPART
C   ------------------------------
C   GAMMA                    1
C   NEUTRINO                 2
C   POSITRON                 3
C   ELECTRON                 4
C   MUON +                   5
C   MUON -                   6
C   PION +                   7
C   PION 0                   8
C   PION -                   9
C   KAON +                  10
C   KAON 0 S                11
C   KAON 0 L                12
C   KAON -                  13
C   PROTON                  14
C   PROTON BAR              15
C   NEUTRON                 16
C   NEUTRON BAR             17
C   LAMBDA                  18
C   LAMBDA BAR              19
C   SIGMA +                 20
C   SIGMA 0                 21
C   SIGMA -                 22
C   SIGMA + BAR             23
C   SIGMA 0 BAR             24
C   SIGMA - BAR             25
C   XSI 0                   26
C   XSI -                   27
C   XSI 0 BAR               28
C   XSI - BAR               29
C   DEUTERON                30
C   TRITON                  31
C   ALPHA                   32
C   OMEGA -                 33
C   OMEGA - BAR             34
C   NEW PARTICLES           35
C
#include "gcunit.inc"
#include "gconsp.inc"
C
C --- FOR CROSS-SECTION INFORMATION SEE "PCSDATA" ---
C
#include "gsecti.inc"
C --- GHEISHA COMMONS ---
#include "s_result.inc"
#include "s_prntfl.inc"
C
      DIMENSION A(KK),Z(KK),W(KK)
C
      DIMENSION ALPHA(35),ALPHAC(41),IPART2(7),CSA(4)
      DIMENSION PARTEL(35),PARTIN(35),ICORR(35),INTRC(35)
C
      PARAMETER (ONETHR=1./3.)
      SAVE ALPHA,ALPHAC,PARTEL,PARTIN,CSA,IPART2,ICORR,INTRC
C
#include "s_csdim.inc"
#include "pcodim.inc"
#include "s_csdat.inc"
#include "pcodat.inc"
C
      DATA ALPHA    / 6*0.7,
     +                0.75 ,0.75 ,0.75 ,
     +                0.76,0.76 ,0.76 ,0.76 ,
     +                0.685,0.63 ,0.685,0.63,0.685,0.63,
     +                3*0.685,3*0.63,2*0.685,2*0.63,
     +                3*0.7,0.685,0.63,0.7/
      DATA ALPHAC    /1.2,1.2,1.2,1.15,0.90,0.91,0.98,1.06,1.10,1.11,
     +                1.10,1.08,1.05,1.01,0.985,0.962,0.945,0.932,
     +                0.925,0.920,0.920,0.921,0.922,0.923,0.928,0.931,
     +                0.940,0.945,0.950,0.955,0.958,0.962,0.965,0.976,
     +                0.982,0.988,0.992,1.010,1.020,1.030,1.040/
      DATA PARTEL/6*0.,29*1./
      DATA PARTIN/6*0.,1.00,0.00,1.05,1.20,1.35,1.30,1.20,1.00,1.30,
     +            1.00,1.30,1.00,1.30,1.00,1.00,1.00,1.30,1.30,1.30,
     +            1.00,1.00,1.30,1.30,1.00,1.,1.,1.,1.3,1./
      DATA ICORR /14*1, 0, 1, 0, 1, 0, 3*1, 3*0, 2*1, 2*0, 4*1, 2*0/
      DATA INTRC /6*0, 1, 0, 12*1, 0, 2*2, 0, 10*1, 0/
C
C CROSS SECTIONS ON NUCLEUS ARE KNOWN ONLY FOR PIONS AND PROTONS.
C THE GENERAL LAW SIGMA(A)=1.25*SIGMA(TOT,PROTON)*A**ALPHA IS VALID
C ONLY FOR MOMENTA > 2 GEV.THE PARAMETRIZATION DONE HERE GIVES ONLY
C A BEHAVIOUR AVERAGED OVER MOMENTA AND PARTICLE TYPES.
C FOR A DETECTOR WITH ONLY A FEW MATERIALS IT'S OF COURSE MUCHBETTER
C TO USE TABLES OF THE MEASURED CROSS SECTIONS .
C FOR ELEMENTS WITH THE FOLLOWING ATOMIC NUMBERS MEASURED CROSS
C SECTIONS ARE AVAILABLE (SEE "PCSDATA").
C
C                 H   AL     CU     PB
      DATA  CSA  /1. ,27.00 ,63.54 ,207.19 /
      DATA IPART2/9,8,7,11,10,13,12/
C
C

C --- Initialise GHESIG and switch to GHEISHA particle code ---
      GHESIG=0.0
      IF(LPART.GT.48)GO TO 160
      IPART=KIPART(LPART)

C
C October 2006 - set flag for new antiproton cross section
C
C IF NEWGHEI=1 NEW CROSS SECTIONS FOR ANTIPROTONS ARE USED
C
C
      NEWGHEI=1
C
C --- No interaction for gammas, neutrinos, electrons, positrons, muons,
C --- neutral pions, neutral sigmas and antisigmas and new particles.
      IF(INTRC(IPART).EQ. 0) GO TO 160
      P=PPART
      EK=EKIN
C
C --- Initialise the cross-sections with 0.0 ---
      DO 10  K=1,KK
         AIEL(K)=0.0
         AIIN(K)=0.0
         AIFI(K)=0.0
         AICA(K)=0.0
   10 CONTINUE
C
      IF ((IPART .GE. 30) .AND. (IPART .LE. 32)) THEN
C
C --- Take geometrical cross sections for inelastic scattering ---
C --- of deuterons, tritons and alphas ---
         IF (IPART .EQ. 30) THEN
            APART=2.0**ONETHR
         ELSEIF (IPART .EQ. 31) THEN
            APART=3.0**ONETHR
         ELSEIF (IPART .EQ. 32) THEN
            APART=4.0**ONETHR
         END IF
         DO 20 K=1,KK
            AIIN(K)=49.0*(APART+A(K)**ONETHR)**2
   20    CONTINUE
         IF (NPRT(9)) THEN
            WRITE(CHMAIL,10000)
            CALL GMAIL(0,0)
         END IF
C
      ELSE IF ((IPART .EQ. 16) .AND. (EK .LE. 0.0327)) THEN
C
C --- Use tables for low energy neutrons ---
C --- get energy bin ---
         JE2=17
         DO 30 J=2,17
            IF (EK .LT. ELAB(J)) THEN
               JE2=J
               GO TO 40
            END IF
   30    CONTINUE
C
   40    JE1=JE2-1
         EKX=MAX(EK,1.0E-9)
         DELAB=ELAB(JE2)-ELAB(JE1)
         DO 70 K=1,KK
C
C --- Get A bin ---
            JA2=15
            DO 50 J=2,15
               IF (A(K) .LT. CNLWAT(J)) THEN
                  JA2=J
                  GO TO 60
               END IF
   50       CONTINUE
C
   60       JA1=JA2-1
            DNLWAT=CNLWAT(JA2)-CNLWAT(JA1)
C
C --- Use linear interpolation or extrapolation by Y=RCE*X+RCA*X+B ---
C
C --- Elastic cross section ---
C --- E interpolation or extrapolation at JA1 ---
            DY=CNLWEL(JA1,JE2)-CNLWEL(JA1,JE1)
            RCE=DY/DELAB
C --- A interpolation or extrapolation at JE1 ---
            DY=CNLWEL(JA2,JE1)-CNLWEL(JA1,JE1)
            RCA=DY/DNLWAT
            B=CNLWEL(JA1,JE1)-RCE*ELAB(JE1)-RCA*CNLWAT(JA1)
            AIEL(K)=RCE*EK+RCA*A(K)+B
C
C --- Inelastic cross section ---
C --- E interpolation or extrapolation at JA1 ---
            DY=CNLWIN(JA1,JE2)-CNLWIN(JA1,JE1)
            RCE=DY/DELAB
C --- A interpolation or extrapolation at JE1 ---
            DY=CNLWIN(JA2,JE1)-CNLWIN(JA1,JE1)
            RCA=DY/DNLWAT
            B=CNLWIN(JA1,JE1)-RCE*ELAB(JE1)-RCA*CNLWAT(JA1)
            AIIN(K)=RCE*EK+RCA*A(K)+B
C
            IZNO=Z(K)+0.01
            AICA(K)=11.12*CSCAP(IZNO)/(EKX*1.0E6)**0.577
   70    CONTINUE
         IF (NPRT(9)) THEN
            WRITE(CHMAIL,10100)
            CALL GMAIL(0,0)
         END IF

            
      ELSE
C
C --- Use parametrization of cross section data for all other cases ---
C
         IF (NPRT(9)) THEN
            WRITE(CHMAIL,10200)
            CALL GMAIL(0,0)
         END IF
C
C --- Get momentum bin ---
         J=40
         DO 80 I=2,41
            IF (P .LT. PLAB(I)) THEN
               J=I-1
               GO TO 90
            END IF
   80    CONTINUE
C
C --- Start with  cross sections for scattering on free protons ---
C --- use linear interpolation or extrapolation by Y=RC*X+B     ---
   90    DX=PLAB(J+1)-PLAB(J)
C --- Elastic cross section ---
         DY=CSEL(IPART,J+1)-CSEL(IPART,J)
         RC=DY/DX
         B=CSEL(IPART,J)-RC*PLAB(J)
         AIELIN=RC*P+B
C --- Inelastic cross section ---
         DY=CSIN(IPART,J+1)-CSIN(IPART,J)
         RC=DY/DX
         B=CSIN(IPART,J)-RC*PLAB(J)
         AIININ=RC*P+B
         ALPH=ALPHA(IPART)
         IF(IPART.LT.14) THEN
            DY=ALPHAC(J+1)-ALPHAC(J)
            RC=DY/DX
            B=ALPHAC(J)-RC*PLAB(J)
            CORFAC=RC*P+B
            ALPH=ALPH*CORFAC
C
            IPART3=IPART2(IPART-6)
C
C --- Elastic cross section ---
            DY=CSEL(IPART3,J+1)-CSEL(IPART3,J)
            RC=DY/DX
            B=CSEL(IPART3,J)-RC*PLAB(J)
            XSECEL=RC*P+B
C --- Inelastic cross section ---
            DY=CSIN(IPART3,J+1)-CSIN(IPART3,J)
            RC=DY/DX
            B=CSIN(IPART3,J)-RC*PLAB(J)
            XSECIN=RC*P+B
C
         END IF
C
         DO 100 K=1,KK
            AIEL(K)=AIELIN
            AIIN(K)=AIININ
C
            IF (A(K) .GE. 1.5) THEN

               IF((IPART.EQ.15).AND.(NEWGHEI.EQ.1) ) THEN
C
C October 2006 - new cross section for antiprotons
C      

                 AIEL(K)=ANTIP_AIEL1(A(K),Z(K),P)
                 AIIN(K)=ANTIP_AII1(A(K),Z(K),P)
                 GO TO 100
	       ENDIF


C
C --- A-dependence from parametrization ---
               CREL=1.0
               CRIN=1.0
C --- Get medium bin  1=Hydr.  2=Al  3=Cu  4=Pb ---
               I=3
               IF (A(K) .LT. 50.0) I=2
               IF (A(K) .GT. 100.0) I=4
               IF ((IPART .EQ. 14) .OR. (IPART .EQ. 16)) THEN
C
C --- Protons and neutrons ---
C
C --- Elastic cross section ---
                  DY=CSPNEL(I-1,J+1)-CSPNEL(I-1,J)
                  RC=DY/DX
                  B=CSPNEL(I-1,J)-RC*PLAB(J)
                  XSECEL=RC*P+B
C --- Inelastic cross section ---
                  DY=CSPNIN(I-1,J+1)-CSPNIN(I-1,J)
                  RC=DY/DX
                  B=CSPNIN(I-1,J)-RC*PLAB(J)
                  XSECIN=RC*P+B
                  IF (AIEL(K) .GE. 0.001) CREL=XSECEL/(0.36*AIEL(K)*
     +            CSA(I)**1.17)
                  AITOT=AIEL(K)+AIIN(K)
                  IF (AITOT .GE. 0.001) CRIN=XSECIN/(AITOT*CSA(I)**
     +            ALPH)
C
               ELSEIF (IPART .LT. 15) THEN
C
C --- Calculate correction factors from values on Al,Cu,Pb for all ---
C --- mesons use linear interpolation or extrapolation by Y=RC*X+B ---
C --- Note that data is only available for pions and protons
                  WGCH=0.5
                  IF (A(K) .LT. 20.0) WGCH=0.5+0.5*EXP(-(A(K)-1.0))
                  AIEL(K)=WGCH*AIEL(K)+(1.0-WGCH)*XSECEL
                  AIIN(K)=WGCH*AIIN(K)+(1.0-WGCH)*XSECIN
C
C --- This section not for kaons ---
                  IF (IPART .LT. 10) THEN
C
C --- Elastic cross section ---
                     DY=CSPIEL(I-1,J+1)-CSPIEL(I-1,J)
                     RC=DY/DX
                     B=CSPIEL(I-1,J)-RC*PLAB(J)
                     XSPIEL=RC*P+B
C --- Inelastic cross section ---
                     DY=CSPIIN(I-1,J+1)-CSPIIN(I-1,J)
                     RC=DY/DX
                     B=CSPIIN(I-1,J)-RC*PLAB(J)
                     XSPIIN=RC*P+B
C
                     IF (AIEL(K) .GE. 0.001) CREL=XSPIEL/(0.36* AIEL(K)
     +               *CSA(I)**1.17)
                     AITOT=AIEL(K)+AIIN(K)
                     IF (AITOT .GE. 0.001) CRIN=XSPIIN/(AITOT*CSA(I)
     +               **ALPH)
                  END IF
               END IF
               AIIN(K)=CRIN*(AIIN(K)+AIEL(K))*A(K)**ALPH
cv               AIIN(K)=CRIN*(AIIN(K))*A(K)**ALPH
               AIEL(K)=CREL*0.36*AIEL(K)*A(K)**1.17
               AIEL(K)=AIEL(K)*PARTEL(IPART)
               AIIN(K)=AIIN(K)*PARTIN(IPART)
            END IF
  100    CONTINUE
C
C --- Fission cross sections ---
C --- A-dependence given by  sigma(3 MeV)=-67.0+38.7*Z**(4/3)/A ---
      END IF
      I=21
      DO 110 II=1,21
         IF (EK .LT. EKFISS(II)) THEN
            I=II
            GO TO 120
         END IF
  110 CONTINUE
C
  120 CONTINUE
      DO 130 K=1,KK
C
C --- No fission for materials with A < 230 ---
         IF (A(K) .GE. 230.) THEN
C
C --- Only data for U(233), U(235) and Pu(239) for Ek .le. 0.01 GeV ---
            J=4
            IF (EK .LE. 0.01) THEN
C
C --- Distinguish U(233), U(235), Pu(239) and U(238) and rest by J=1,4
               IF ((Z(K) .EQ. 92.0).AND.(ABS(A(K)-233.0).LT.0.5))THEN
                  J=1
               ELSEIF((Z(K).EQ.92.0).AND.(ABS(A(K)-235.0).LT.0.5))THEN
                  J=2
               ELSEIF((Z(K).EQ.94.0).AND.(ABS(A(K)-239.0).LT.0.5))THEN
                  J=3
               END IF
            END IF
            IF(J.EQ.4) THEN
C
               Z43BA=Z(K)**(4.0*ONETHR)/A(K)
               Z43BA=MAX(-67.0+38.7*Z43BA,0.)
            ELSE
               Z43BA=1.
            ENDIF
C
C --- Energy dependence taken from U(238) ---
C --- approximated as step-function ---
            AIFI(K)=CSFISS(J,I)*Z43BA
         END IF
  130 CONTINUE
C
C --- Corrections for compounds ---
C --- These corrections should only be applied to inorganic scintill. ---
C --- apply the correction only if user selected it within GEANT ---
      IF (QCOR .GT. 0.0.AND.KK.GT.1) THEN
C --- Do not apply corrections for anti-baryons ---
         IF (ICORR(IPART).EQ.1) THEN
C
C --- ACC40 between 0.3 and 0.5 for Pi+ and Pi0, 0.2 for other mesons ---
            ACC40=0.325
C           IF (IPART .GE. 9) ACC40=0.20
C --- ACC40 = 0.15 for baryons ---
            IF (IPART .GE. 14) ACC40=0.15
C --- ACCA=0.08 for all pions, 0.02 for all other particles ---
            ACCA= 0.08
            IF (IPART .GE. 10) ACCA=0.02
C
            ACCB=0.32*(ACC40-ACCA)
            ACC=ACCA-ACCB*LOG(EK)
            IF (ACC .GT. 0.5) ACC=0.5
            IF (ACC .GT. 0.0) THEN
C
               CAVER=AVER**ACC
               CAVER=1.+(CAVER-1.)*QCOR
               DO 140 K=1,KK
                  AIEL(K)=AIEL(K)*CAVER
                  AIIN(K)=AIIN(K)*CAVER
                  AIFI(K)=AIFI(K)*CAVER
                  AICA(K)=AICA(K)*CAVER
  140          CONTINUE
            END IF
         END IF
      END IF
C
C --- Calculate interaction probability ---
C
C --- Correction factor for high (P > 100 GeV/C) energies ---
      CORH=1.0
C --- Assume a LOG(P) dependence of the correction factor with values ---
C P = 100 GeV/C  ==> CORH = 1.
C P =   1 TeV/C  ==> CORH = 1.25
C     DX=LOG(1000.)-LOG(100.)
C     DY=1.25-1.
C     RC=DY/DX
C     B=1.-RC*LOG(100.)
C     CORH=RC*LOG(P)+B

      IF (P .GT. 100.) CORH=0.1085736156*LOG(P)+0.5
      ALAM=0.0
      DO 150 K=1,KK
         AFACT=AVO*1E-3*DENS*W(K)/A(K)
         AIEL(K)=MAX(CORH*AIEL(K)*AFACT,0.)
         AIIN(K)=MAX(CORH*AIIN(K)*AFACT,0.)
         AIFI(K)=MAX(CORH*AIFI(K)*AFACT,0.)
         AICA(K)=MAX(CORH*AICA(K)*AFACT,0.)
C
         ALAM=ALAM+AIEL(K)+AIIN(K)+AIFI(K)+AICA(K)
  150 CONTINUE
C
C --- Pass the interaction probability to GEANT ---
      GHESIG=ALAM
C
      GO TO 999
C
C --- Printout of skipped particles in case of interface debug ---
  160 IF (NPRT(9)) THEN
         WRITE(CHMAIL,10300) IPART
         CALL GMAIL(0,0)
      END IF
10000 FORMAT(' *GHESIG* GEOM X-SECT. FOR INEL. SCAT. OF D,T AND ALPHA')
10100 FORMAT(' *GHESIG* X-SECT. FROM LOW ENERGY NEUTRON TABLES')
10200 FORMAT(' *GHESIG* X-SECT. FROM PARAMETRIZATION OF DATA')
10300 FORMAT('0*GHESIG* GHEISHA PARTICLE ',I3,' SKIPPED')
  999 END


      FUNCTION ANTIP_AIEL1(A,Z,P)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C AUTHOR : Sergio Bottai , October 2006
C
C
C ANTIPROTON ELASTIC CROSS SECTION(mb) 
C
C Bendiscioli and Kharzeev - Nuovo Cimento Vol. 17 N 6
C
C The parametrization for antiproton-Deuterum elastic cross 
C section given in (Bendiscioli and Kharzeev - 
C  Nuovo Cimento Vol. 17 N 6) is extrapolated to other
C nuclei following the A**(2/3) scaling. 
C
C MUST BE IN  millibarn FOR COMPATIBILITY WITH GHESIG
C
C
C CALLED BY : GHESIG
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      
      REAL A,P,Z 
      
      P1=P
      IF(P1.LT.0.27) P1=0.27
      IF(P1.GT.100.) P1=100.
      
      
      IF(P1.LT.2.) THEN
       SIGT=129.+78./P1
       SIGR=105.+22.5/(P1**1.58)-14.9*LOG10(P1)
       ANTIP_AIEL1=SIGT-SIGR
      ELSE IF(P1.GE.2..AND.P1.LT.5.) THEN
        ANTIP_AIEL1=1.8*(37.5-2.63/(P1**2.58)-12.6*LOG10(P1))
      ELSE IF(P1.GE.5.) THEN
         ANTIP_AIEL1=(112.+125./P1**1.08+1.15*LOG10(P1)**2
     +    -12.6*LOG10(P1))/2.33   
      ENDIF
      FACTO=A**(2./3.)/2.**(2./3.)
      ANTIP_AIEL1=ANTIP_AIEL1*FACTO
          

      END
      
      FUNCTION ANTIP_AII1(A,Z,P)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C AUTHOR : Sergio Bottai , October 2006
C
C
C
C ANTIPROTON INELASTIC CROSS SECTION(mb)
C
C SIG=SIG0*A**ALPHA
C
C WE USE A PARAMETRIZATION OF SIG0(P),ALPHA(P)
C WHICH APPROXIMATE BETTER THAN 10% THE VALUE
C TABULATED IN Bendiscioli and Kharzeev - Nuovo Cimento Vol. 17 N 6
C 
C
C MUST BE IN  millibarn FOR COMPATIBILITY WITH GHESIG
C
C
C CALLED BY : GHESIG
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC       

      REAL A,P,Z

      P1=P
      IF(P1.GT.100.) P1=100.

      AL=0.6+0.011*SQRT(P1) 
      IF(P1.LT.0.1) P1=0.1
      SIG0=70.+20.*P1**(-.4)-10.*(LOG10(2.*P1+1.)**(1.3))
      ANTIP_AII1=SIG0*A**AL
      END      
      
