* * $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