************************************************************************ * ------------------------------------------- SUBROUTINE NENCELSLPV(ENEU,PF,IPAR,MODE,PLEP,PNUC,IERR) * ------------------------------------------- * * (Purpose) * set final lepton momentum vector for elastic scattering * * (Input) * ENE : NEUTRINO ENERGY ( GEV ) * PF(3) : FERMI MOMENTUM ( GEV/C ) * IPAR : NEUTRINO TYPE * 12 : NEU-E * -12 : NEU-E-BAR * 14 : NEU-MU * -14 : NEU-MU-BAR * 16 : NEU-TAU * -16 : NEU-TAU-BAR * MODE : Interaction mode * * (Output) * PLEP(3) : FINAL LEPTON MOMENTUM ( GEV/C ) * IERR : ERROR CODE * * (Creation Date and Author) * 2011.03.15 ; Y.Hayato * ************************************************************************ IMPLICIT NONE #include "necard.h" #include "neutparams.h" REAL ENEU,PF(3),PLEP(3),PNUC(3) INTEGER*4 IPAR,MODE,IERR INTEGER*4 iplep,i,ilast,loop REAL*4 AMIN,AMOUT SAVE ICOUNT,IEVENT INTEGER ICOUNT,IEVENT REAL pfabs,e,beta,gm,enstop,q2maxx,dq2,bmax,blast,q2 real el,pl,cost,amp,dum,hit,sigma REAL fnncq2max,fnnpot,rlu EXTERNAL fnncq2max,fnnpot,rlu REAL*8 DNNCELSQ2 EXTERNAL DNNCELSQ2 REAL PNEU(3),EV(3),DNEU(3),ENF REAL pneulab(3),pnulabi(3),ptotli(3),plepf(3),pnf(3) REAL einnuc,einit,elep,efinal,ptotlis,pnfabs,pnfabs2 REAL phiang,cosphi,ptotliy2 REAL*8 repa,repb,repc,repd,repf,repg,reph REAL*8 DELEP,DTMPY,DENEU REAL*8 DPLEPF(3) REAL*8 DBLE0, DBLQ2 real*8 DSGN REAL XMAP(51,11) REAL FNNUCLPF integer*4 isub,iloop DATA ICOUNT/0/ DATA IEVENT/1/ INTEGER*4 ITARG C IERR=0 DSGN=0.D0 DENEU=ENEU IF (abs(mode).eq.51) then ITARG = 2212 ELSE IF (abs(mode).eq.52) then ITARG = 2112 else write(*,*) 'NENCELSLPV: Invalid mode',MODE STOP ENDIF CALL MCMASS(ITARG,AMIN) AMOUT = AMIN AMIN = AMIN*1.E-3 AMOUT = AMOUT*1.E-3 C C -- CALCUALTE NEUTRINO ENERGY IN NUCELON STOP SYSTEM C C--reduce loop to 1000times. DO 100 LOOP=0,1000 PNEU(1)=ENEU PNEU(2)=0. PNEU(3)=0. C-- set pf (not considering pauli blocking or potential) call nencelsspff(eneu,ipar,itarg,pf,LOOP,ierr) if (ierr.ne.0) goto 9000 PFABS=SQRT(PF(1)**2+PF(2)**2+PF(3)**2) E=SQRT(PFABS**2+AMIN**2) BETA=PFABS/E GM=1./SQRT(1.-BETA**2) IF (PFABS.NE.0.) THEN EV(1)=-PF(1)/PFABS EV(2)=-PF(2)/PFABS EV(3)=-PF(3)/PFABS ELSE EV(1)=1. EV(2)=0. EV(3)=0. ENDIF CALL MCVECBST(PNEU,0.,EV,GM) ENSTOP=SQRT(PNEU(1)**2+PNEU(2)**2+PNEU(3)**2) IF(ENSTOP.LT.0.001 .AND. IABS(IPAR).EQ.12)GO TO 9000 C IF(ENSTOP.LT.0.112 .AND. IABS(IPAR).EQ.14)GO TO 9000 IF(ENSTOP.LT.0.030 .AND. IABS(IPAR).EQ.14)GO TO 9000 C IF(ENSTOP.LT.3.48109 .AND. IABS(IPAR).EQ.16)GO TO 9000 IF(ENSTOP.LT.0.030 .AND. IABS(IPAR).EQ.16)GO TO 9000 DNEU(1)=PNEU(1)/ENSTOP DNEU(2)=PNEU(2)/ENSTOP DNEU(3)=PNEU(3)/ENSTOP C C -- MAKE (Q**2 -- SCATTERING ANGLE) MAP C CALL VZERO(XMAP,561) EV(1)=-EV(1) EV(2)=-EV(2) EV(3)=-EV(3) Q2MAXX=FNNCQ2MAX(ENSTOP,ITARG) DQ2=Q2MAXX/50. BMAX=0. BLAST=0.E+0 FNNUCLPF=FNNPOT(1,PFABS) C C -- SEARCH PEAK ROUGHLY C DO 10 I=0,50 Q2=DQ2*FLOAT(I) EL=ENSTOP-Q2/2./AMIN PL=EL COST=(2.*ENSTOP*EL-Q2)/2./ENSTOP/PL IF(ABS(COST).GT.1.)GO TO 10 DBLE0=ENSTOP DBLQ2=Q2 AMP=DNNCELSQ2(DBLE0,IPAR,ITARG,DBLQ2) IF (BMAX.LT.AMP) BMAX=AMP IF (AMP.LT.BLAST) GOTO 20 BLAST=AMP 10 CONTINUE 20 CONTINUE ILAST=I-1 C C -- SEARCH PEAK CLOSELY C Q2=DQ2*FLOAT(ILAST)+DQ2*0.02 EL=ENSTOP-Q2/2./AMIN PL=EL COST=(2.*ENSTOP*EL-Q2)/2./ENSTOP/PL IF(ABS(COST).GT.1.) THEN AMP=0.0E+0 ELSE DBLE0=ENSTOP DBLQ2=Q2 AMP=DNNCELSQ2(DBLE0,IPAR,ITARG,DBLQ2) ENDIF IF (AMP.GT.BLAST) THEN C C -- FORWARD C BMAX=AMP BLAST=AMP DO 30 I=2,50 Q2=DQ2*FLOAT(ILAST)+DQ2*0.02*FLOAT(I) EL=ENSTOP-Q2/2./AMIN PL=EL COST=(2.*ENSTOP*EL-Q2)/2./ENSTOP/PL IF(ABS(COST).GT.1.)GO TO 30 DBLE0=ENSTOP DBLQ2=Q2 AMP=DNNCELSQ2(DBLE0,IPAR,ITARG,DBLQ2) IF (BMAX.LT.AMP) BMAX=AMP IF (AMP.LT.BLAST) GOTO 50 BLAST=AMP 30 CONTINUE ELSE IF (ILAST.NE.0) THEN C C -- BACKWORD C DO 40 I=1,50 Q2=DQ2*FLOAT(ILAST)-DQ2*0.02*FLOAT(I) EL=ENSTOP-Q2/2./AMIN PL=EL COST=(2.*ENSTOP*EL-Q2)/2./ENSTOP/PL IF(ABS(COST).GT.1.)GO TO 40 DBLE0=ENSTOP DBLQ2=Q2 AMP=DNNCELSQ2(DBLE0,IPAR,ITARG,DBLQ2) IF (BMAX.LT.AMP) BMAX=AMP IF (AMP.LT.BLAST) GOTO 50 BLAST=AMP 40 CONTINUE ENDIF 50 CONTINUE C C -- THROW DICE C EINNUC=sqrt(PFABS**2+AMIN**2) Einit=EINNUC + ENEU #ifdef DEBUG write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') $ 'INLEP =(',ENEU, ',',0, ',',0 , $ ')/Einlep =',eneu write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') $ 'INNUC =(',PF(1), ',',PF(2), ',',PF(3), $ ')/Einnuc =',einnuc write(*,'(A10,F12.8,A1)') $ 'q2=(',q2,')' #endif C DO 100 LOOP=1,1000000 Q2=RLU(DUM)*Q2MAXX DBLQ2=Q2 HIT=RLU(DUM)*BMAX SIGMA=DNNCELSQ2(DBLE0,IPAR,ITARG,DBLQ2) IF (SIGMA.GT.HIT) THEN PHIANG=6.283*RLU(DUM) cosphi=COS(PHIANG) if (COSPHI.gt.0) THEN DSGN=1.D0 ELSE DSGN=-1.D0 ENDIF PNEULAB(1)=ENEU PNEULAB(2)=0. PNEULAB(3)=0. DO 1500 ISUB=1,3 PNULABI(ISUB)=PF(ISUB) 1500 continue C--- for the comparison --- C 1200 EL=ENSTOP-Q2/2./AMIN C PL=EL C COST=(2.*ENSTOP*EL-Q2)/2./ENSTOP/PL C C DLEP(1)=COST*DNEU(1)-sqrt(1.-COST**2)*DNEU(2) C DLEP(2)=sqrt(1.-COST**2)*DNEU(1)+COST*DNEU(2) C DLEP(3)=0. C C do 1210 JJ=1,3 C PLEP(JJ)=PL*DLEP(JJ) C 1210 continue C C CALL MCVECBST(PLEP,0,EV,GM) C el = sqrt(plep(1)**2+plep(2)**2+plep(3)**2) C write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') C $ 'OUTLEPX=(',PLEP(1),',',PLEP(2),',', C $ PLEP(3),')/Eoutlep=',el C write(*,'(A10,F12.8)') 'q^2=', C $ ( (el-eneu)**2 C $ -(plep(1)-eneu)**2-plep(2)**2-plep(3)**2) PNULABI(1)=PF(1) PNULABI(2)=PF(2) PNULABI(3)=0. PTOTLI(1)=PNULABI(1)+ENEU PTOTLI(2)=PNULABI(2) PTOTLI(3)=PNULABI(3) PTOTLIS=(PTOTLI(1)**2+PTOTLI(2)**2+PTOTLI(3)**2) PTOTLIY2=(PTOTLI(2)**2)*COSPHI**2 PNFABS=PFABS DO 1510 ILOOP=1,5 REPA=-1.*q2 C REPB=dble(ENEU)+dble(einnuc) C $ +(FNNUCLPF-FNNUCL(PNFABS)) REPB=dble(ENEU)+dble(einnuc) $ +(FNNUCLPF-FNNPOT(2,PNFABS)) REPC=REPB**2-dble(PTOTLIS)-dble(AMOUT)**2 $ +dble(PTOTLI(1))*REPA/dble(ENEU) REPD=2*(REPB-PTOTLI(1)) REPF= REPC**2+(PTOTLIY2)*(REPA**2)/(ENEU**2) REPG=4*REPA*(PTOTLIY2)/ENEU-2*REPC*REPD REPH=REPG**2-4*REPD**2*REPF C write(*,'(A,F12.7,A,F12.7,A,F12.7,A,F12.7)') C $ 'A=',REPA,', B=',REPB,', C=',REPC,', D=',REPD C write(*,'(A,F12.7,A,F12.7,A,F12.7)') C $ 'F=',REPF,', G=',REPG,', H=',REPH IF (ABS(REPH).lt.1.D-8) REPH=0.D0 IF (REPH.LT.0) THEN if (QUIET.eq.0) then write(*,*) 'neelslpv:ILOOP=',ILOOP, $ 'value in SQRT lt.0', $ '(REPH=',REPH,')' endif IF ((ILOOP.eq.1).or.(ILOOP.eq.5)) then goto 100 endif reph=0. endif 1515 DELEP=(-1.*REPG+DSGN*dsqrt(REPH))/(2*REPD**2) ELEP = REAL(DELEP) DPLEPF(1)=-q2/(2*DENEU)+DELEP DTMPY = DELEP**2-DPLEPF(1)**2 IF (DTMPY.LT.0) THEN if (QUIET.eq.0) then write(*,*) $ 'neelslpv:Insufficient energy of LEPTON', $ '(P_Y^2=',DTMPY,')' endif IF ((ILOOP.eq.1).or.(ILOOP.eq.5)) then goto 100 endif C plepf(1)=0. C plepf(2)=0. C plepf(3)=0. C goto 1575 ENF=REPB-ELEP PNFABS2=ENF**2-AMOUT**2 IF (PNFABS2.lt.0) PNFABS2=0. goto 1585 endif PLEPF(1)=DPLEPF(1) if (abs(PTOTLI(2)).gt.0.000001) then PLEPF(2)=(REPD*ELEP-REPC)/(2*PTOTLI(2)) else PLEPF(2)= (-((REPA/ENEU)**2)/4. $ -(REPA/ENEU)*ELEP)*(COSPHI**2) C PLEPF(2) = (DELEP**2-DPLEPF(1)**2) C $ *(COSPHI**2) IF (PLEPF(2).gt.0.) THEN PLEPF(2)=sqrt(PLEPF(2))*DSGN else PLEPF(2)=0. ENDIF endif IF ( PLEPF(1)**2+PLEPF(2)**2-elep*elep $ .gt.0.0001) THEN IF ((ILOOP.eq.1).or.(ILOOP.eq.5)) then goto 100 endif C plepf(1)=0. C plepf(2)=0. C plepf(3)=0. C goto 1575 ENF=REPB-ELEP PNFABS2=ENF**2-AMOUT**2 IF (PNFABS2.lt.0) PNFABS2=0. goto 1585 endif PLEPF(3)=real(dsqrt(DTMPY))*SIN(PHIANG) 1575 PNF(1)=ENEU+PF(1)-PLEPF(1) PNF(2)=PF(2)-PLEPF(2) PNF(3)=PF(3)-PLEPF(3) PNFABS2=(PNF(1)**2+PNF(2)**2+PNF(3)**2) ENF = sqrt(PNFABS2+AMOUT**2) EFINAL = ELEP + ENF 1585 PNFABS = sqrt(PNFABS2) C write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') C $ 'OUTLEP=(',PLEPF(1),',',PLEPF(2),',',PLEPF(3), C $ ')/Eoutlep=',elep C write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') C $ 'OUTNUC=(',PNF(1), ',',PNF(2), ',',PNF(3), C $ ')/Eoutnuc=',enf C write(*,'(A10,F12.8,A,F12.8)') 'q2(lep)=(', C $ -1*( (ELEP-ENEU)**2-(PLEPF(1)-ENEU)**2 C $ -(PLEPF(2)-0.)**2-(PLEPF(3)-0.)**2) C $ ,') : should be ',q2 C C write(*,'(A10,F12.8,A10,F12.8)') C $ 'Ein : ',Einit, 'Eout : ',EFINAL C C write(*,'(A10,F12.3,A10,F12.3)') C $ 'Ediff : ',(Einit-Efinal)*1000., C $ 'Vnucl: ' ,(FNNPOT(2,PNFABS)-FNNUCLPF)*1000. CC $ 'Vnucl: ' ,(FNNUCL(PNFABS)-FNNUCLPF)*1000. 1510 CONTINUE IF(NEPAUFLG.EQ.0 .AND. PNFABS.LT.PFSURF) GOTO 100 IF (PLEPF(1)/sqrt(ELEP**2).le.-0.998) THEN goto 100 endif #ifdef DEBUG write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') $ 'OUTLEP=(',PLEPF(1),',',PLEPF(2),',',PLEPF(3), $ ')/Eoutlep=',elep write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') $ 'OUTNUC=(',PNF(1), ',',PNF(2), ',',PNF(3), $ ')/Eoutnuc=',enf write(*,'(A10,F12.8,A,F12.8)') 'q2(lep)=(', $ -1*( (ELEP-ENEU)**2-(PLEPF(1)-ENEU)**2 $ -(PLEPF(2)-0.)**2-(PLEPF(3)-0.)**2) $ ,') : should be ',q2 write(*,'(A10,F12.8,A10,F12.8)') $ 'Ein : ',Einit, 'Eout : ',EFINAL write(*,'(A10,F12.3,A10,F12.3)') $ 'Ediff : ',(Einit-Efinal)*1000., $ 'Vnucl: ' ,(FNNPOT(2,PNFABS)-FNNUCLPF)*1000. C $ 'Vnucl: ' ,(FNNUCL(PNFABS)-FNNUCLPF)*1000. write(*,'(A10,F12.3,A10)') $ 'Econserv: ', $ (EINIT+FNNUCLPF-(EFINAL+FNNPOT(2,PNFABS)))*1000, $ 'MeV' C $ (EINIT+FNNUCLPF-(EFINAL+FNNUCL(PNFABS)))*1000, C $ 'MeV' #endif C IF ((abs(EINIT+FNNUCLPF-(EFINAL+FNNUCL(PNFABS)))*1000 IF ((abs(EINIT+FNNUCLPF-(EFINAL+FNNPOT(2,PNFABS)))*1000 $ .GT.2)) THEN write(*,'(A,I8,A)') $ 'EDIFF > 2.0MeV (COUNT=',ICOUNT,')' ICOUNT = ICOUNT+1 ENDIF PLEP(1)=PLEPF(1) PLEP(2)=PLEPF(2) PLEP(3)=PLEPF(3) PNUC(1)=PNF(1) PNUC(2)=PNF(2) PNUC(3)=PNF(3) C write(*,*) '---------------------------------------' C write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') C $ 'OUTLEP=(',PLEP(1),',',PLEP(2),',',PLEP(3), C $ ')/Eoutlep=',elep C write(*,'(A10,F12.8,A,F12.8,A,F12.8,A,F12.8)') C $ 'OUTNUC=(',PNUC(1), ',',PNUC(2), ',',PNUC(3), C $ ')/Eoutnuc=',enf C write(*,*) '---------------------------------------' RETURN ENDIF 100 if (QUIET.eq.0) write(*,*) 'neelslpv: retry..' CONTINUE WRITE(6,*) ' IN NEELSLPV ( TOO MANY TRY ) ' C -Add 97/08/19 _ Y.H. for errortrap IERR=1 RETURN C C ++ ERROR RETURN C 9000 WRITE(6,900) 900 FORMAT(' ***ERROR IN NEELSLPV(NOT ENOUGH NEUTRINO ENERGY)***') IERR=1 RETURN END