*
* $Id: gpxtr.F,v 3.1.1.1 2002/07/11 16:02:01 cafagna Exp $
*
* $Log: gpxtr.F,v $
* Revision 3.1.1.1  2002/07/11 16:02:01  cafagna
* First GPAMELA release on CVS
*
*
*CMZ :  3.00/00 05/02/2002  12.51.38  by  Unknown
*-- Author :    Marialuigia Ambriola   11/05/2001
      SUBROUTINE GPXTR
**************************************************************************
*                                                                        *
*       Generates and absorbs transition radiatixon spectrum in TRD       *
*                                                                        *
*                                                                        *
*                                                                        *
* Called by: GUSTEP                                                      *
* Author:  Marialuigia Ambriola   11/05/2001                             *
*                                                                        *
**************************************************************************
#include "gpgene.inc"
#include "gctrak.inc"
#include "gctmed.inc"
#include "gpmed.inc"
#include "gpsed.inc"
#include "gcsets.inc"
#include "gpaltr.inc"
#include "gckine.inc"
#include "gcnum.inc"
#include "gcflag.inc"
#include "gptotr.inc"
#include "gcunit.inc"
	EXTERNAL XTRYIELD
	EXTERNAL NPOISS
        EXTERNAL XTRINTER
        EXTERNAL XTRINTEG
        REAL CAM,RL,G,NF,WR,CR,SPH,CATT,W,E
        REAL GRNDM0(1)
        INTEGER INDEX,INDEXAIR,INDEXC,LAY
c centimeter and GeV:
        DATA RL1/7.0E-4/,RL2/204.E-4/,N/71/,w1/31.59E-9/,w2/.7E-9/
        DATA IEVOLD/-1/,NTMOLD/-1/
        PARAMETER (ALPHA=1/137.,PG=3.1415926536,HTC=197.3269E-16)
        PARAMETER (C=3.E+10)
**************
*particle inside kapton, radiator or xenon?
*if yes, calcolates the length of the track.
*NUMED is the number of the tracking medium (common GCTMED)
*MTRAD, MKAP and MXE are the numbers of radiator, kapton and xenon
*tracking media, like defined by GPAMELA.
        IF(CHARGE.EQ.0) THEN
           RETURN
        ENDIF
        IF(AMASS.LE.0) THEN
           WRITE(CHMAIL,10000) AMASS
           CALL GMAIL(1,0)
           RETURN
        ENDIF
        G=GETOT/AMASS
	IF(G.LE.100.) THEN
           RETURN
	ENDIF
          IF((NUMED.EQ.MKAP.AND.IDET.EQ.IDTRSO).OR.NUMED.EQ.
     +       MTRAD.OR.NUMED.EQ.MXE) THEN
	   IF(INWVOL.EQ.1.) THEN
             CAM=SLENG
	     RETURN
           ELSEIF(INWVOL.EQ.2) THEN
C #
C # Calculate the layer number
C #
              IF(NUMED.EQ.MXE) THEN
               IF(NUMBV(NVNAME-1).GE.29.AND.NUMBV(NVNAME-1).LE.32)LAY=2
               IF(NUMBV(NVNAME-1).GE.25.AND.NUMBV(NVNAME-1).LE.28)LAY=4
               IF(NUMBV(NVNAME-1).GE.21.AND.NUMBV(NVNAME-1).LE.24)LAY=6
               IF(NUMBV(NVNAME-1).GE.17.AND.NUMBV(NVNAME-1).LE.20)LAY=8
               IF(NUMBV(NVNAME-1).GE.13.AND.NUMBV(NVNAME-1).LE.16)LAY=10
               IF(NUMBV(NVNAME-1).GE.10.AND.NUMBV(NVNAME-1).LE.12)LAY=12
               IF(NUMBV(NVNAME-1).GE.7.AND.NUMBV(NVNAME-1).LE.9)LAY=14
               IF(NUMBV(NVNAME-1).GE.4.AND.NUMBV(NVNAME-1).LE.6)LAY=16
               IF(NUMBV(NVNAME-1).GE.1.AND.NUMBV(NVNAME-1).LE.3)LAY=18
               IF(NUMBV(NVNAME).GT.16) THEN
                  LAY = LAY - 1
               ENDIF
              ENDIF
             CAM=SLENG-CAM
	   ELSE
             RETURN
	   ENDIF
	  ELSE
	   RETURN
          ENDIF
****************
	  INDEX=0
	  INDEXAIR=0
	  INDEXC=0
          IF(INWVOL.NE.2) RETURN
	  IF(IEVOLD.NE.IEVENT.OR.NTMOLD.NE.NTMULT) THEN
             CALL VZERO(EY,115)
	  ENDIF
          IEVOLD=IEVENT
          NTMOLD=NTMULT
          SPH=0.
          NPHTR=0
          ENPHTR=0.
          IF(NUMED.EQ.MKAP.AND.IDET.EQ.IDTRSO) INDEX=4
          IF(NUMED.EQ.MXE) INDEX=3
          IF(NUMED.EQ.MTRAD) THEN
             INDEXAIR=1
             INDEXC=2
             RL=RL1+RL2
          ENDIF
****************
c radiator:
          IF(INDEX.EQ.0) THEN
	     NF=CAM/RL
             DO I=1,115
              XTREN=XTRYIELD(I,G)
              WR=0.
              CR=NF*(ATTTRD(I,INDEXAIR)*RL2+ATTTRD(I,INDEXC)*RL1)
              IF(CR.GT.0.) WR=EXP(-CR)
c considering the whole radiator and the absorption.
c when dividing for enatt I consider the number of photons per unit energy,
c which are created and transmitted by the radiator.
C remember that enatt is in keV unity.
              IF(CR.GT.0.AND.NF.GT.1) THEN
                 XTREN=2*NF*XTREN*(1.-WR)/(CR*ENATT(I))
                 EY(I)=EY(I)*WR+XTREN
                 SPH=SPH+EY(I)
	      ENDIF
             ENDDO
          ENDIF
          goto 999
          print*,'index,IEVENT,gamma,spettro'
          print*,index,IEVENT,g
          print*,(ey(i),i=1,4)
          print*,(ey(i),i=5,8)
          print*,(ey(i),i=9,12)
          print*,(ey(i),i=13,16)
          print*,(ey(i),i=17,20)
          print*,(ey(i),i=21,24)
          print*,(ey(i),i=25,28)
          print*,(ey(i),i=29,32)
          print*,(ey(i),i=33,36)
          print*,(ey(i),i=37,40)
          print*,(ey(i),i=41,44)
          print*,(ey(i),i=45,48)
          print*,(ey(i),i=49,52)
          print*,(ey(i),i=53,56)
          print*,(ey(i),i=57,60)
          print*,(ey(i),i=61,64)
          print*,(ey(i),i=65,68)
          print*,(ey(i),i=69,72)
          print*,(ey(i),i=73,76)
          print*,(ey(i),i=77,80)
          print*,(ey(i),i=81,84)
          print*,(ey(i),i=85,88)
          print*,(ey(i),i=89,92)
          print*,(ey(i),i=93,96)
          print*,(ey(i),i=97,100)
          print*,(ey(i),i=101,104)
          print*,(ey(i),i=105,108)
          print*,(ey(i),i=109,112)
          print*,(ey(i),i=113,115)
 999      continue
******************
C now the spectrum is propagated in the other absorbent media in the TRD:
c the absorber could be the kapton or the Xe-CO2 or the gas in between (here
c not considered):
          IF(INDEX.NE.0) THEN
             CALL VZERO(EA,115)
             CALL VZERO(ER,115)
             DO I=1,115
                CATT=ATTTRD(I,INDEX)*CAM
                W=0
                W=EXP(-CATT)
C               PRINT*,'EY(I),I,CATT,ATT',EY(I),I,CATT,ATTTRD(I,INDEX)
c absorbed spectrum:
                EA(I)=EY(I)*(1-W)
c transmitted spectrum:
                EY(I)=EY(I)*W
C               PRINT*,'W,EA(I),EY(I),I,INDEX',W,EA(I),EY(I),I,INDEX
             ENDDO
          ENDIF
C finally the sensitive gas:
          IF(INDEX.EQ.3) THEN
             goto 998
             print*,'index,IEVENT,gamma,spettro assorbito'
             print*,index,IEVENT,g
             print*,(ea(i),i=1,4)
             print*,(ea(i),i=5,8)
             print*,(ea(i),i=9,12)
             print*,(ea(i),i=13,16)
             print*,(ea(i),i=17,20)
             print*,(ea(i),i=21,24)
             print*,(ea(i),i=25,28)
             print*,(ea(i),i=29,32)
             print*,(ea(i),i=33,36)
             print*,(ea(i),i=37,40)
             print*,(ea(i),i=41,44)
             print*,(ea(i),i=45,48)
             print*,(ea(i),i=49,52)
             print*,(ea(i),i=53,56)
             print*,(ea(i),i=57,60)
             print*,(ea(i),i=61,64)
             print*,(ea(i),i=65,68)
             print*,(ea(i),i=69,72)
             print*,(ea(i),i=73,76)
             print*,(ea(i),i=77,80)
             print*,(ea(i),i=81,84)
             print*,(ea(i),i=85,88)
             print*,(ea(i),i=89,92)
             print*,(ea(i),i=93,96)
             print*,(ea(i),i=97,100)
             print*,(ea(i),i=101,104)
             print*,(ea(i),i=105,108)
             print*,(ea(i),i=109,112)
             print*,(ea(i),i=113,115)
 998         continue
c considering the escape peak in Xe-CO2;
             CALL XESCAPE(ENATT,EA,ER,115)
	     SPH=0.
             goto 997
             print*,'index,IEVENT,gamma,spettro con escape'
             print*,index,IEVENT,g
             print*,(er(i),i=1,4)
             print*,(er(i),i=5,8)
             print*,(er(i),i=9,12)
             print*,(er(i),i=13,16)
             print*,(er(i),i=17,20)
             print*,(er(i),i=21,24)
             print*,(er(i),i=25,28)
             print*,(er(i),i=29,32)
             print*,(er(i),i=33,36)
             print*,(er(i),i=37,40)
             print*,(er(i),i=41,44)
             print*,(er(i),i=45,48)
             print*,(er(i),i=49,52)
             print*,(er(i),i=53,56)
             print*,(er(i),i=57,60)
             print*,(er(i),i=61,64)
             print*,(er(i),i=65,68)
             print*,(er(i),i=69,72)
             print*,(er(i),i=73,76)
             print*,(er(i),i=77,80)
             print*,(er(i),i=81,84)
             print*,(er(i),i=85,88)
             print*,(er(i),i=89,92)
             print*,(er(i),i=93,96)
             print*,(er(i),i=97,100)
             print*,(er(i),i=101,104)
             print*,(er(i),i=105,108)
             print*,(er(i),i=109,112)
             print*,(er(i),i=113,115)
 997         continue
c yes escape peak:
ccc             SPH=XTRINTEG(ENATT,ER,EINT,115)
             SPH=XTRINTEG(ENATT,ER,EINT,75)
c if you don't consider the escape-peak, don't comment the next row:
C              SPH=XTRINTEG(ENATT,EA,EINT,115)
C 28/9/2001: FINE TUNING OF TR:
c for PS ALFA=0.65 (momentum.lt.40. GeV/c)
c per SPS ALFA=0.69 (momentum.ge.40. GeV/c):
             IF(ALFATR.NE.-9999.) THEN
                SPH=ALFATR*SPH
             ELSE
                TROK=.FALSE.
                WRITE(CHMAIL,10100) ALFATR
                CALL GMAIL(1,0)
                RETURN
             ENDIF
c end 28/9/2001.
             NPHTR=NPOISS(SPH)
c             print*,'lay=',lay,'sph,nph=',sph,nph
c             PRINT*,'Xe-CO2: XTRINTEG,NPOISS:',SPH,NPHTR
             IF(NPHTR.GT.0) THEN
                DO I=1,NPHTR
ccc 10                E=XTRINTER(RNDM(1)*SPH,EINT,ENATT,115)
c ml 26/03/04
c 10                E=XTRINTER(RNDM(1)*SPH,EINT,ENATT,75)
 10                CALL GRNDM(GRNDM0,1)
                   E=XTRINTER(GRNDM0(1)*SPH,EINT,ENATT,75)
C ml 26/03/04
C 28/9/2001: TENTATIVO DI FINE TUNING DELLA TR:
******** NEI DATI REALI NON VEDO FOTONI CON ENERGIA MAGGIORE DI 60 keV:
c fine modifica 28/9/2001
                   IF(E.GT.45.) THEN
                      PRINT*,'Extracted energy'
                      PRINT*,E
                      GOTO 10
                   ENDIF
                   ENPHTR=ENPHTR+E
                ENDDO
c                PRINT*,'Energy (keV) and photons absorbed, for event:'
c     +               ,IEVENT,ENPHTR,NPHTR ,'layer=',lay
C               NTRHIT=NPHTR
C               ENEHIT=ENPHTR*1.E-6
C* ML: SOTTRAE LA ENERGIA PERSA NEL TRD.
                GEKINT=GEKIN-ENPHTR*1.E-6
                GEKIN=GEKINT
                GETOT=GEKIN +AMASS
                VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
C Routine to find bin number in kinetic energy table stored in ELOW(NEKBIN)
                CALL GEKBIN
C* END ML.
             ELSE
c                PRINT*,'no photons converted in the straw'
             ENDIF
             IF (LAY.GT.0.AND.LAY.LE.18) THEN
                ENTRTOT = ENTRTOT + ENPHTR
                ENLAY(LAY) = ENLAY(LAY) + ENPHTR
                NPHTOTR = NPHTOTR + NPHTR
                NPHLAY(LAY)= NPHLAY(LAY) + NPHTR
                ILI=NUMBV(NVNAME-1)
                IHI=NUMBV(NVNAME)
                ENTRSTRAW(ILI,IHI) = ENTRSTRAW(ILI,IHI)+ ENPHTR
                NTRSTRAW(ILI,IHI) = NTRSTRAW(ILI,IHI)+ NPHTR
c                PRINT*,'NPHTOTR,LAY,ENPHTR,ENTRTOT',NPHTOTR,LAY,ENPHTR,ENTRTOT
C                PRINT*,'NTRSTRAW',NTRSTRAW(ILI,IHI),ILI,IHI
C                PRINT*,'ENTRSTRAW,ENPHTR',ENTRSTRAW(ILI,IHI),ENPHTR
             ELSE
                WRITE(CHMAIL,10200)LAY
                CALL GMAIL(1,0)
C                PRINT*,' ERROR IN LAYER CALCULATION',LAY
             ENDIF
          ENDIF
10000     FORMAT('GPXTR error: negative mass =',F10.3)
10100     FORMAT('GPXTR error: ALFATR not set, TR process ',
     +  'and his tuning will be ignored, ALFATR =',F10.3)
10200     FORMAT('GPXTR error in layer calculation =',F10.3)
          RETURN
          END
