[u/mrichter/AliRoot.git] / GEANT321 / neutron / thrsel.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:59  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.48  by  S.Giani
*-- Author :
      SUBROUTINE THRSEL(NE,NP7,NB7,E,EOUT,FM,TDK,ETH,ALPHA,BETA,W,
     +                  PMUPS,PMDS,F,AWR,IIN,IPTMD,IPMDS,IOUT)
C      THIS ROUTINE SELECTS THE EXIT ENERGY AND SCATTERING ANGLE
C      FROM S(ALPHA,BETA) DATA TABLES
      DIMENSION ETH(*),ALPHA(*),BETA(*),ABETA2(2),PMDS(*),
     +   IPTMD(NE),W(NE),IPMDS(*),PMUPS(NB7,NE),F(NP7,NB7),AWR(*)
      SAVE
C
      AAWR=AWR(IIN)
      DO 10 IE=1,NE
         IF(E.LT.ETH(IE))GO TO 20
   10 CONTINUE
      INDX=NE
      GO TO 30
   20 INDX=IE
      IF(INDX.EQ.1)GO TO 30
      R=FLTRNF(0)
      DELE=(ETH(INDX)-E)/(ETH(INDX)-ETH(INDX-1))
      DELEN=DELE
      GO TO 40
   30 DELEN=1.0
      INDX=2
   40 PROB=DELEN*W(INDX-1)+(1.0-DELEN)*W(INDX)
      IF(R.LE.PROB)GO TO 120
C       NEUTRON DOWNSCATTERS
      R=FLTRNF(0)
      DO 90 I=1,2
         IP=IPTMD(INDX)
         NB=IPMDS(IP)
         DO 50 IB=1,NB
            IF(R.LE.PMDS(IP+NB+IB))GO TO 60
   50    CONTINUE
         WRITE(IOUT,10000)PMDS(IP+2*NB)
10000 FORMAT(' MICAP: CUMULATIVE DOWNSCATTER DIST. DOES NOT END ',
     +       'IN 1.0 IN THRSEL',E12.4)
         PRINT *,' CALOR: ERROR in MICAP ====> STOP '
         STOP
   60    IF(IB.EQ.1)GO TO 70
         DELE=(PMDS(IP+NB+IB)-R)/(PMDS(IP+NB+IB)-PMDS(IP+NB+IB-1))
         ABETA=DELE*(PMDS(IP+IB-1)-PMDS(IP+IB))+PMDS(IP+IB)
         GO TO 80
   70    ABETA=BETA(IB)
   80    ABETA2(I)=ABETA
         INDX=INDX-1
   90 CONTINUE
      ABETA=DELEN*ABETA2(2)+(1.0-DELEN)*ABETA2(1)
      EOUT=E-TDK*ABETA
      IF(EOUT.LT.1.0E-05)EOUT=1.0E-05
      DO 100 IB=1,NB7
         IF(ABETA.LE.BETA(IB))GO TO 110
  100 CONTINUE
      IB=NB7
  110 DELE=(ABETA-BETA(IB))/(BETA(IB-1)-BETA(IB))
      GO TO 180
C       NEUTRON UPSCATTERS
  120 R=FLTRNF(0)
      DO 170 I=1,2
         DO 130 IB=1,NB7
            IF(R.LE.PMUPS(IB,INDX))GO TO 140
  130    CONTINUE
         WRITE(IOUT,10100)PMUPS(NB7,INDX)
10100 FORMAT(' MICAP: CUMULATIVE UPSCATTER DIST. DOES NOT END ',
     +       'IN 1.0 IN THRSEL',E12.4)
         PRINT *,' CALOR: ERROR in MICAP ====> STOP '
         STOP
  140    IF(IB.EQ.1)GO TO 150
         DELE=(PMUPS(IB,INDX)-R)/(PMUPS(IB,INDX)-PMUPS(IB-1,INDX))
         ABETA=DELE*(BETA(IB-1)-BETA(IB))+BETA(IB)
         GO TO 160
  150    WRITE(IOUT,10200)PMUPS(1,INDX)
10200 FORMAT(' MICAP: CUMULATIVE UPSCATTER DIST. DOES NOT BEGIN ',
     +       'AT 0.0 IN THRSEL',E12.4)
         PRINT *,' CALOR: ERROR in MICAP ====> STOP '
         STOP
  160    ABETA2(I)=ABETA
         INDX=INDX-1
  170 CONTINUE
      ABETA=DELEN*ABETA2(2)+(1.0-DELEN)*ABETA2(1)
      EOUT=E+TDK*ABETA
C       SELECT ANGLE
  180 AMAX=(EOUT+E+2.0*SQRT(E*EOUT))/(AAWR*TDK)
      AMIN=(EOUT+E-2.0*SQRT(E*EOUT))/(AAWR*TDK)
      DO 190 IA=1,NP7
         IF(AMAX.LT.ALPHA(IA))GO TO 200
  190 CONTINUE
      IA=NP7
      DELA=0.0
      GO TO 210
  200 DELA=(ALPHA(IA)-AMAX)/(ALPHA(IA)-ALPHA(IA-1))
  210 F4=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
      F3=DELE*(F(IA-1,IB-1)-F(IA-1,IB))+F(IA-1,IB)
      F2=DELA*(F3-F4)+F4
      DO 220 IA=1,NP7
         IF(AMIN.LT.ALPHA(IA))GO TO 230
  220 CONTINUE
      IA=NP7
      DELA=0.0
      GO TO 240
  230 DELA=(ALPHA(IA)-AMIN)/(ALPHA(IA)-ALPHA(IA-1))
  240 F4=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
      F3=DELE*(F(IA-1,IB-1)-F(IA-1,IB))+F(IA-1,IB)
      F1=DELA*(F3-F4)+F4
      R=FLTRNF(0)
      F0=R*F2+(1.0-R)*F1
      F1=0.0
      DO 250 IA=1,NP7
         F2=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
         IF(F0.LE.F2)GO TO 260
         F1=F2
  250 CONTINUE
  260 IF(F1.EQ.F2)GO TO 270
      DELA=(F2-F0)/(F2-F1)
      GO TO 280
  270 ALP=ALPHA(IA)
      GO TO 290
  280 ALP=DELA*ALPHA(IA-1)+(1.0-DELA)*ALPHA(IA)
  290 FM=(E+EOUT-ALP*AAWR*TDK)/(2.0*SQRT(E*EOUT))
      IF(ABS(FM).LE.1.0)RETURN
      WRITE(IOUT,10300)FM,E,EOUT,R,IA,IB
10300 FORMAT(' MICAP: ERROR IN THRSEL, COSINE OF ANGLE >1'/,
     +' ',1P4E12.4,2I11)
      FM=2.0*FLTRNF(0)-1.0
      RETURN
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:21:59 cernlib
	6	* Geant
	7	*
	8	*
	9	#include "geant321/pilot.h"
	10	*CMZ : 3.21/02 29/03/94 15.41.48 by S.Giani
	11	*-- Author :
	12	SUBROUTINE THRSEL(NE,NP7,NB7,E,EOUT,FM,TDK,ETH,ALPHA,BETA,W,
	13	+ PMUPS,PMDS,F,AWR,IIN,IPTMD,IPMDS,IOUT)
	14	C THIS ROUTINE SELECTS THE EXIT ENERGY AND SCATTERING ANGLE
	15	C FROM S(ALPHA,BETA) DATA TABLES
	16	DIMENSION ETH(),ALPHA(),BETA(),ABETA2(2),PMDS(),
	17	+ IPTMD(NE),W(NE),IPMDS(),PMUPS(NB7,NE),F(NP7,NB7),AWR()
	18	SAVE
	19	C
	20	AAWR=AWR(IIN)
	21	DO 10 IE=1,NE
	22	IF(E.LT.ETH(IE))GO TO 20
	23	10 CONTINUE
	24	INDX=NE
	25	GO TO 30
	26	20 INDX=IE
	27	IF(INDX.EQ.1)GO TO 30
	28	R=FLTRNF(0)
	29	DELE=(ETH(INDX)-E)/(ETH(INDX)-ETH(INDX-1))
	30	DELEN=DELE
	31	GO TO 40
	32	30 DELEN=1.0
	33	INDX=2
	34	40 PROB=DELENW(INDX-1)+(1.0-DELEN)W(INDX)
	35	IF(R.LE.PROB)GO TO 120
	36	C NEUTRON DOWNSCATTERS
	37	R=FLTRNF(0)
	38	DO 90 I=1,2
	39	IP=IPTMD(INDX)
	40	NB=IPMDS(IP)
	41	DO 50 IB=1,NB
	42	IF(R.LE.PMDS(IP+NB+IB))GO TO 60
	43	50 CONTINUE
	44	WRITE(IOUT,10000)PMDS(IP+2*NB)
	45	10000 FORMAT(' MICAP: CUMULATIVE DOWNSCATTER DIST. DOES NOT END ',
	46	+ 'IN 1.0 IN THRSEL',E12.4)
	47	PRINT *,' CALOR: ERROR in MICAP ====> STOP '
	48	STOP
	49	60 IF(IB.EQ.1)GO TO 70
	50	DELE=(PMDS(IP+NB+IB)-R)/(PMDS(IP+NB+IB)-PMDS(IP+NB+IB-1))
	51	ABETA=DELE*(PMDS(IP+IB-1)-PMDS(IP+IB))+PMDS(IP+IB)
	52	GO TO 80
	53	70 ABETA=BETA(IB)
	54	80 ABETA2(I)=ABETA
	55	INDX=INDX-1
	56	90 CONTINUE
	57	ABETA=DELENABETA2(2)+(1.0-DELEN)ABETA2(1)
	58	EOUT=E-TDK*ABETA
	59	IF(EOUT.LT.1.0E-05)EOUT=1.0E-05
	60	DO 100 IB=1,NB7
	61	IF(ABETA.LE.BETA(IB))GO TO 110
	62	100 CONTINUE
	63	IB=NB7
	64	110 DELE=(ABETA-BETA(IB))/(BETA(IB-1)-BETA(IB))
65	GO TO 180
66	C NEUTRON UPSCATTERS
67	120 R=FLTRNF(0)
68	DO 170 I=1,2
69	DO 130 IB=1,NB7
70	IF(R.LE.PMUPS(IB,INDX))GO TO 140
71	130 CONTINUE
72	WRITE(IOUT,10100)PMUPS(NB7,INDX)
73	10100 FORMAT(' MICAP: CUMULATIVE UPSCATTER DIST. DOES NOT END ',
74	+ 'IN 1.0 IN THRSEL',E12.4)
75	PRINT *,' CALOR: ERROR in MICAP ====> STOP '
76	STOP
77	140 IF(IB.EQ.1)GO TO 150
78	DELE=(PMUPS(IB,INDX)-R)/(PMUPS(IB,INDX)-PMUPS(IB-1,INDX))
79	ABETA=DELE*(BETA(IB-1)-BETA(IB))+BETA(IB)
80	GO TO 160
81	150 WRITE(IOUT,10200)PMUPS(1,INDX)
82	10200 FORMAT(' MICAP: CUMULATIVE UPSCATTER DIST. DOES NOT BEGIN ',
83	+ 'AT 0.0 IN THRSEL',E12.4)
84	PRINT *,' CALOR: ERROR in MICAP ====> STOP '
85	STOP
86	160 ABETA2(I)=ABETA
87	INDX=INDX-1
88	170 CONTINUE
89	ABETA=DELENABETA2(2)+(1.0-DELEN)ABETA2(1)
90	EOUT=E+TDK*ABETA
91	C SELECT ANGLE
92	180 AMAX=(EOUT+E+2.0SQRT(EEOUT))/(AAWR*TDK)
93	AMIN=(EOUT+E-2.0SQRT(EEOUT))/(AAWR*TDK)
94	DO 190 IA=1,NP7
95	IF(AMAX.LT.ALPHA(IA))GO TO 200
96	190 CONTINUE
97	IA=NP7
98	DELA=0.0
99	GO TO 210
100	200 DELA=(ALPHA(IA)-AMAX)/(ALPHA(IA)-ALPHA(IA-1))
101	210 F4=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
102	F3=DELE*(F(IA-1,IB-1)-F(IA-1,IB))+F(IA-1,IB)
103	F2=DELA*(F3-F4)+F4
104	DO 220 IA=1,NP7
105	IF(AMIN.LT.ALPHA(IA))GO TO 230
106	220 CONTINUE
107	IA=NP7
108	DELA=0.0
109	GO TO 240
110	230 DELA=(ALPHA(IA)-AMIN)/(ALPHA(IA)-ALPHA(IA-1))
111	240 F4=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
112	F3=DELE*(F(IA-1,IB-1)-F(IA-1,IB))+F(IA-1,IB)
113	F1=DELA*(F3-F4)+F4
114	R=FLTRNF(0)
115	F0=RF2+(1.0-R)F1
116	F1=0.0
117	DO 250 IA=1,NP7
118	F2=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
119	IF(F0.LE.F2)GO TO 260
120	F1=F2
121	250 CONTINUE
122	260 IF(F1.EQ.F2)GO TO 270
123	DELA=(F2-F0)/(F2-F1)
124	GO TO 280
125	270 ALP=ALPHA(IA)
126	GO TO 290
127	280 ALP=DELAALPHA(IA-1)+(1.0-DELA)ALPHA(IA)
128	290 FM=(E+EOUT-ALPAAWRTDK)/(2.0SQRT(EEOUT))
129	IF(ABS(FM).LE.1.0)RETURN
130	WRITE(IOUT,10300)FM,E,EOUT,R,IA,IB
131	10300 FORMAT(' MICAP: ERROR IN THRSEL, COSINE OF ANGLE >1'/,
132	+' ',1P4E12.4,2I11)
133	FM=2.0*FLTRNF(0)-1.0
134	RETURN
135	END