C C ++ GAMMA,PI- SCATTERING SIMULATION FOR O16 C IMPLICIT NONE REAL PI PARAMETER (PI=3.141593) REAL XIN(3),PIN(3),PFI(3) REAL H(50000) REAL RMAX,TIN,PPP,PINABS,XZIG,PFIABS,TPI REAL THETA,WEIGHT,COST,SINT,WEIGT2 REAL AREA,DTHETA INTEGER KIN,ITOT,IABS,IFREE,ICHEX,I INTEGER NABS,NZIG,KFI,IINEL REAL EGAM,RADIUS,S,PGAM(3),PNI(3),PF,ABSPNI REAL ABSCM,GAMMA,BETA,BSTVEC(3),ROTPHI,ROTCOS REAL PGAMCM(3),PNICM(3),PPIABS,PITTA,PIPHI,PPICM(3) REAL DUM,APN,AP REAL COSPI,PGAMABS,DIFCRS INTEGER IDUM C REAL RNDM C EXTERNAL RNDM REAL RLU EXTERNAL RLU REAL EFABRHO EXTERNAL EFABRHO REAL EFFRMMOM EXTERNAL EFFRMMOM REAL EFGAMMN EXTERNAL EFGAMMN INTEGER LTINEL,LTABS,LTCOH INTEGER ISEED COMMON /PAWC/H COMMON /RNSEED/ISEED #include COMMON /YHDEBUG/LTINEL,LTABS,LTCOH C CALL HLIMIT(50000) CALL HBOOK1(1,'SCATTERING ANGLE DISTRIBUTION COS(THETA)$',100, & -1.,1.,0.) CALL HBOOK1(2,'SCATTERING ANGLE DISTRIBUTION THETA$',90, & 0.,180.,0.) CALL HBOOK1(3,'FINAL MOMENTUM DISTRIBUTION$',100,0.,300.,0.) CALL HBOOK1(4,'FINAL KINETIC ENERGY DISTRIBUTION$',100,0.,200.,0.) CALL HBOOK1(5,'PPI (THETA = 44.2)$',100,0.,500.,0.) CALL HBOOK1(6,'PPI (THETA = 44.2)$',100,0.,500.,0.) CALL HBOOK1(10,'NZIG$',10,0.,10.,0.) CALL HBOOK1(30,'NZIG BEFORE ABSORPTION$',50,0.,50.,0.) CALL HBOOK1(20,'SCATTERING ANGLE DISTRIBUTION THETA(WEIGHTED)$', &90,0.,180.,0.) CALL HBOOK1(21,'DIFFERENTIAL CROSS-SECTION$', &90,0.,180.,0.) CALL HBOOK1(31,'KINETIC ENERGY BEFORE ABSORPTION$',100, &0.,200.,0.) CALL HBOOK2(32,'PLACE WHERE ABSORBED(X .VS. INPACT PARAMETER)$', &100,-7.,7.,50,0.,7.,0.) CALL HBOOK1(33,'COS(THETA) BEFORE ABSORPTION$',100, &-1.,1.,0.) CALL HBOOK1(40,'SCATTERING ANGLE DISTRIBUTION THETA(WEIGHTED)' & //' CHARGE EXCHANGE PI+,PI0$',90,0.,180.,0.) CALL RNSET(121761) CALL HBOOK1(99,'RANDOM NUMBER DISTRIBUTION',100,0.,1.,0.) C C C READ PARAMETERS C C WRITE(*,*) "RANDOM SEED" READ(*,*) ISEED WRITE(*,*) "EVENT #" READ(*,*) ITOT WRITE(*,*) "E(GAMMA)" READ(*,*) EGAM WRITE(*,*) "EVENT : ",ITOT WRITE(*,*) "E(GAMMA) : ",EGAM C C C C FOR DEBUG C LTINEL=0 LTABS=0 LTCOH=0 C C SET PARAMETERS C RMAX=2.5*C C (RMAX: MAXIMUM RADIUS,C:WOODS-SAXON CONST.) KIN=-211 C (KIN:KIND OF INCOMING PION,PI+ -> 211,PI- -> -211,PI0 -> 111) C AREA=PI*RMAX*RMAX*10 C (AREA:AREA OF CIRCLE RADIUS RMAX FM, C 1FM^2 = 10MB,SO X10 WAS ADDED) C DTHETA=PI/180.*2. C (DTHETA IS FOR CALC. CROSS-SECTION, C CURRENTRY 2 DEGREE / BIN.) C IABS=0 IFREE=0 ICHEX=0 IINEL=0 DO 100 I=1,ITOT C C SET INTERACTION POINT C C WRITE(*,*) "EVENT NO. ",I 105 CONTINUE C XIN(1)=-RMAX+RNDM(DUM)*RMAX*2. C XIN(2)=-RMAX+RNDM(DUM)*RMAX*2. C XIN(3)=-RMAX+RNDM(DUM)*RMAX*2. XIN(1)=-RMAX+RLU(IDUM)*RMAX*2. XIN(2)=-RMAX+RLU(IDUM)*RMAX*2. XIN(3)=-RMAX+RLU(IDUM)*RMAX*2. RADIUS = SQRT(XIN(1)**2 + XIN(2)**2 + XIN(3)**2) IF (RADIUS.GT.RMAX) GOTO 105 C IF(EFABRHO(RADIUS).LT.RNDM(DUM)) GOTO 105 IF(EFABRHO(RADIUS).LT.RLU(IDUM)) GOTO 105 C C SET INITIAL NUCLEON AND GAMMA MOMENTUM C C 110 PF = 275.*RNDM(DUM) 110 PF = 275.*RLU(IDUM) C IF (EFFRMMOM(PF).GT.RNDM(DUM)) GOTO 110 IF (EFFRMMOM(PF).GT.RLU(IDUM)) GOTO 110 C 101 PNI(1)=(-1.+RNDM(DUM)*2)*PF C PNI(2)=(-1.+RNDM(DUM)*2)*PF 101 PNI(1)=(-1.+RLU(IDUM)*2)*PF PNI(2)=(-1.+RLU(IDUM)*2)*PF PNI(3)=(PF**2-PNI(1)**2-PNI(2)**2) IF (PNI(3).LT.0.) GOTO 101 PNI(3)=SQRT(PNI(3)) C IF (RNDM(DUM).GT.0.5) THEN IF (RLU(IDUM).GT.0.5) THEN PNI(3)=PNI(3)*(-1.) ENDIF PGAM(1)=0. PGAM(2)=0. PGAM(3)=EGAM C C CALC BETA AND BOOST VECTOR(UNIT VECTOR) C ABSPNI = PF ABSCM = SQRT(PNI(1)**2 + PNI(2)**2 + (PNI(3)+PGAM(3))**2) BETA = ABSCM/(SQRT(939.**2+ABSPNI)+EGAM) GAMMA = 1/SQRT((1-BETA**2)) BSTVEC(1) = -PNI(1)/ $ SQRT(PNI(1)**2 + PNI(2)**2 + (PNI(3)+PGAM(3))**2) BSTVEC(2) = -PNI(2)/ $ SQRT(PNI(1)**2 + PNI(2)**2 + (PNI(3)+PGAM(3))**2) BSTVEC(3) = -(PNI(3)+PGAM(3))/ $ SQRT(PNI(1)**2 + PNI(2)**2 + (PNI(3)+PGAM(3))**2) CALL MCVECBST(PNI,939.,BSTVEC,GAMMA) C CALL VECBST(PGAM,0.,BSTVEC,GAMMA) ABSPNI = SQRT(PNI(1)**2 + PNI(2)**2 + PNI(3)**2) ROTPHI = ATAN(-PNI(2)/PNI(1)) ROTCOS = ACOS(-PNI(3)/ABSPNI) PNICM(1)=COS(ROTPHI)*PNI(1)-SIN(ROTPHI)*PNI(2) PNICM(2)=SIN(ROTPHI)*PNI(1)+COS(ROTPHI)*PNI(2) PNICM(3)=COS(ROTCOS)*PNI(3)-SIN(ROTCOS)*PNICM(1) PNICM(1)=SIN(ROTCOS)*PNI(3)+COS(ROTCOS)*PNICM(1) PGAMCM(1)=-PNICM(1) PGAMCM(2)=-PNICM(2) PGAMCM(3)=-PNICM(3) PGAMABS=SQRT(PGAMCM(1)**2+PGAMCM(2)**2+PGAMCM(3)**2) C COSPI=RNDM(DUM) COSPI=RLU(IDUM) 102 DIFCRS=EFGAMMN(PGAMABS,COSPI) IF (DIFCRS.LT.0.0001) GOTO 100 C IF (DIFCRS/30.0.LT.RNDM(DUM)) THEN IF (DIFCRS/30.0.LT.RLU(IDUM)) THEN C COSPI=RNDM(DUM) COSPI=RLU(IDUM) GOTO 102 ENDIF PITTA=ACOS(COSPI) S=(ABSPNI+SQRT(ABSPNI**2+939.**2)) PINABS=SQRT((S**4+139.**4+939.**4-2*S**2*(939.**2+139.**2) $ -2*139.*2*939.*2)/4/S**2) PPICM(3)=PINABS*COS(PITTA) PPICM(1)=PINABS*SIN(PITTA) C PIPHI=RNDM(DUM)*2*3.14 PIPHI=RLU(IDUM)*2*3.14 PIN(3)=PPICM(3) PIN(1)=PPICM(1)*COS(PIPHI) PIN(2)=PPICM(1)*SIN(PIPHI) PNI(1)=PIN(1) PNI(2)=PIN(2) PNI(3)=PIN(3) BSTVEC(1)=-BSTVEC(1) BSTVEC(2)=-BSTVEC(2) BSTVEC(3)=-BSTVEC(3) CALL MCVECBST(PIN,139.,BSTVEC,GAMMA) CALL MCVECBST(PNI,939.,BSTVEC,GAMMA) IF ((PNI(1)**2+PNI(2)**2+PNI(3)**2).LT.275.**2) GOTO 100 C C C CALL EFTRACE(APN,AP,KIN,XIN,PIN,NABS,NZIG,KFI,PFI) XZIG=FLOAT(NZIG)+0.5 CALL HF1(10,XZIG,1.) IF(NABS.EQ.0)GO TO 200 IF(KIN.NE.KFI) GO TO 100 PFIABS=SQRT(PFI(1)**2+PFI(2)**2+PFI(3)**2) TPI=SQRT(PFIABS**2+139.**2)-139. COST=PGAM(3)*PFI(3)/EGAM/PFIABS THETA=ACOS(COST) THETA=THETA/3.141593*180. SINT=SQRT(1.-COST**2) WEIGHT=1./SINT WEIGT2=WEIGHT*AREA/FLOAT(ITOT)/(2 * PI * DTHETA) CALL HF1(1,COST,1.) CALL HF1(2,THETA,1.) CALL HF1(20,THETA,WEIGHT) CALL HF1(21,THETA,WEIGT2) CALL HF1(3,PFIABS,1.) CALL HF1(4,TPI,1.) IF(ABS(COST-0.717) .LT. 0.1)CALL HF1(5,PFIABS,WEIGT2) IINEL=IINEL+1 GO TO 100 200 IABS=IABS+1 PFIABS=SQRT(PFI(1)**2+PFI(2)**2+PFI(3)**2) TPI=SQRT(PFIABS**2+139.**2)-139. COST=PIN(1)*PFI(1)/PINABS/PFIABS C PACT=SQRT(XOUT(2)**2+XOUT(3)**2) CALL HF1(30,XZIG,1.) CALL HF1(31,TPI,1.) C CALL HF2(32,XOUT(1),PACT,1.) CALL HF1(33,COST,1.) GO TO 100 100 CONTINUE C C IINEL=ITOT-IABS-IFREE-ICHEX WRITE(100,600)ITOT,IFREE,IABS,IINEL,ICHEX,RMAX 600 FORMAT(///5X,'TOTAL # OF EVENT :',I6/ &5X,'FREE ESCAPE :',I6/ &5X,'ABSORPTION :',I6/ &5X,'INELASTIC SCATTER:',I6/ &5X,'CHARGE EXCHANGE :',I6/ &5X,'RMAX :',G15.7,' FERMI') WRITE(200,*) 'INEL=',LTINEL,' :ABS=',LTABS,' :COH=',LTCOH C CALL HISTDO C CALL HSTORE(0,10) CALL HRPUT(1,'COS-DIST.HB','N') CALL HRPUT(2,'ANGL-DIST.HB','N') CALL HRPUT(5,'MOM-DIST.HB','N') CALL HRPUT(20,'WTDANGL-DIST.HB','N') CALL HRPUT(21,'CRSSECT-DIST.HB','N') CALL HRPUT(3,'FINMOM-DIST.HB','N') CALL HRPUT(10,'ZIG-DIST.HB','N') C CALL HRPUT(50,'FMMMOM-DIST.HB','N') CALL HRPUT(99,'RANDOM-DIST.HB','N') STOP END