C********************************************* C %W% modified on %G% C program Kinematics implicit none C************************************************************************ C @(#)fiducial.cdk 1.1 modified on 12/29/92 real FV common/fidv/FV integer i1 C**** Declare space for any hbook stuff that may be done later integer hmax parameter (hmax=50000) integer hmem(hmax) common /pawc/ hmem C**** Declare a seed value integer iseed data iseed /0/ C**** The File name to read the mc input from. character*132 inpfile C**** The File Name to write the kinematics output to. character*132 outfile C**** The number of arguments on the command line. integer iargc C**** The argument character*132 argv C**** The flux to use integer Use_Flux common /Flux_To_Use/ Use_Flux C**** Flag that only contained events are to be counted. Events outside the C fiducial volume will be generated, but not counted. integer contained common /Contained_Events/ contained C**** The number of files found. integer FileFound fv = 0.0 contained = -1 FileFound = 0 do i1 = 1, iargc() call getarg(i1,argv) if (argv(1:2).eq.'-s') then read(argv(3:132),*) iseed write(6,*) 'seed read from command line:',iseed else if (argv(1:5).eq.'-flux') then read(argv(6:132),*) Use_Flux else if (argv(1:2).eq.'-f') then read(argv(3:132),*) fv write(6,*) 'fiducial distance read from command line:',fv else if (argv(1:2).eq.'-c') then contained = 1 write(6,*) 'only counting contained events.' else FileFound = FileFound + 1 if (FileFound.eq.1) then inpfile = argv write(6,*) 'input file from command line: ',inpfile(1:32) C*** Open up the input file. open(unit=5, $ file = inpfile, $ status='old', $ form='formatted', $ readonly, $ shared) else if (FileFound.eq.2) then outfile = argv write(6,*) 'output file from command line: ',outfile(1:32) C*** Open up the input file. open(unit=6, $ file = outfile, $ status='new', $ form='formatted') endif endif enddo C*** Seed the random number generator if (iseed.eq.0) then call initseed(iseed) type *, 'Using Random Seed',iseed endif call seedranf(iseed) C*** Initialize for hbook. call hlimit(hmax) C*** Generate the events. CALL SPOT close(21) stop end subroutine getarg_parser() implicit none integer i1 C**** The number of arguments in the command line integer nargs C**** The length of the command line integer larg0 C**** The command line character*1024 carg0 C**** The individual arguements integer nchar character*128 cargs(40) common /getarg_common/ nargs, larg0, carg0, cargs data nargs /-1/ data larg0 /0/ C**** The state and states for the parser. integer state, START, HAVEARG, HAVEQUOTE parameter (START=0,HAVEARG=1,HAVEQUOTE=2) if (nargs.ge.0) return C**** Get the command line call lib$get_foreign(carg0,0,larg0) C**** Now Parse it. nargs = 0 do i1 = 1,40 cargs(i1) = ' ' enddo state = START do i1 = 1,larg0 if (state.eq.START) then if (carg0(i1:i1).eq.'''') then nargs = nargs + 1 state = HAVEQUOTE nchar = 0 else if (carg0(i1:i1).ne.' ') then nargs = nargs + 1 state = HAVEARG nchar = 1 cargs(nargs)(nchar:nchar) = carg0(i1:i1) else state = START endif else if (state.eq.HAVEARG) then if (carg0(i1:i1).eq.' ') then state = START else nchar = nchar + 1 cargs(nargs)(nchar:nchar) = carg0(i1:i1) state = HAVEARG endif else if (state.eq.HAVEQUOTE) then if (carg0(i1:i1).eq.'''') then state = START else nchar = nchar + 1 cargs(nargs)(nchar:nchar) = carg0(i1:i1) state = HAVEQUOTE endif endif enddo return end integer function iargc() implicit none C**** The number of arguments in the command line integer nargs C**** The length of the command line integer larg0 C**** The command line character*1024 carg0 C**** The individual arguements integer nchar character*128 cargs(40) common /getarg_common/ nargs, larg0, carg0, cargs call getarg_parser() if (nargs.gt.0) then iargc = nargs else iargc = 0 endif return end subroutine getarg(iarg,carg) implicit none integer i1 C**** The number of arguments in the command line integer nargs C**** The length of the command line integer larg0 C**** The command line character*1024 carg0 C**** The individual arguements integer nchar character*128 cargs(40) common /getarg_common/ nargs, larg0, carg0, cargs C**** The index of the argument to return integer iarg C**** The string of the command argument to return character*(*) carg call getarg_parser() if (iarg.ge.1.and.iarg.le.nargs) then carg = cargs(iarg) do i1 =1,len(carg) if (carg(i1:i1).ge.'A'.and.carg(i1:i1).le.'Z') then carg(i1:i1) = char(ichar(carg(i1:i1)) $ - ichar('A') + ichar('a')) endif enddo else if (iarg.eq.0) then carg = 'COMMAND' else carg = ' ' endif return end C********************************************************************** C Provide a ranf call if there is not a good random number generater C defined. real function ranf() implicit none integer seedranf integer iseed, seed data iseed /987654321/ ranf = ran(iseed) return entry seedranf(seed) iseed = seed return end subroutine INITSEED(iseed) integer ijk(2) call sys$gettim(ijk) iseed = iabs(ijk(1)).or.1 return end C************************************************* C %W% modified on %G% C SUBROUTINE SPOT implicit none integer i1, i2 real r1 C SPOT: C 1) READS THE DATA FROM UNIT=5 C 2) CALLS KINEM FOR EVERY EVENT THAT IS GOING TO BE GENERATED C 3) WRITES THE RESULTS ON UNIT=1 IN BINARY FORM C AND ON UNIT=20 IN FORMATTED FORM C******************************************************************** C C**** The input unit. integer inunit parameter (inunit=5) C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC C**** The common block with the detector wall positions in it. real HT(3) COMMON/GEO/HT character *60 text character *32 file1,file2 C******************************************************************** C @(#)prdec.cdk 1.1 modified on 12/29/92 C pass information about the initail particls and the desired products. C**** The amass and info about the interaction from spot. real amint, amlep, amd, amm(5) integer np(5), npar, kt COMMON/PRDEC/AMINT,AMLEP,AMD,AMM,NP,NPAR,KT C******************************************************************* C @(#)inter.cdk 1.1 modified on 12/29/92 C the interaction point. C real xint(3) common/inter/xint C**************************************************************** C @(#)momentum.cdk 1.2 modified on 1/12/93 C Include the information about the primary particles to be tracked. C**** The number of momenta. integer nmom C**** The particle 4 momentum. real ppr(4,20) C**** The particle masses real ppm(20) C**** The particle charges real chh(20) C**** The particle id. integer pid(20) COMMON/MOM/ nmom,ppr,ppm,chh,pid C***** Gamma ray information real ppg(4,20) integer ngam common /gamm1/ ppg,ngam integer inot, itrans, idec, inel, iel, iabs, icx common /hadron/inot,itrans,idec,inel,iel,iabs,icx logical fileopen data fileopen /.false./ C**** The number of events to generat integer nev, iev C**** Flag for the event type. integer key C**** Flag for if fermi momentum is generated. integer ifer C**** A mode variable that is not used. integer mode integer moda C**** A random seed to start the program with. integer iseed C**** Flag this run as a test integer ktest C**** Flag that only contained events are to be counted. Events outside the C fiducial volume will be generated, but not counted. integer contained common /Contained_Events/ contained 1000 continue call mcinit C*************************************************************** C NEV IS NUMBER OF EVENTS C KEY IS -3 for atmospheric muon c -2 for thru mu C -1 FOR SINGLE TRACKS C 0 FOR PDK C 1 FOR NU+NUCLEON=LEPTON+DELTA C 2 FOR NU+NUCLEON=LEPTON+PRODUCTS (ONEPI MODEL) C IFER=0 FOR FREE NUCLEONS C =1 FOR NUCLEONS WITH FERMI MOTION C MODE FOR KEY=-1 IRRELEVANT C KEY=0 irrelevant c KEY=1 MODE=0 IF A MIXTURE OF NUEL/NUMU C 10 ',5) else if (i1.lt.npar) then call setcat(text,t1,' + ',3) endif enddo call setcat(text,t1,' **********',11) call setcat(file1,f1,'.dk91',5) call setcat(file2,f2,'.kk91',5) return end subroutine setcat(str1,lst1,str2,lst2) implicit none C**** The string to add the characters to. character*(*) str1 C**** The last character in str1 integer lst1 C**** The string to add to str2. character*(*) str2 C**** The number of characters in str2. integer lst2 C**** some counters. integer i1, i2 C**** The character to istert the string after in str1. integer ins do i1 = 1,lst2 lst1 = lst1 + 1 str1(lst1:lst1) = str2(i1:i1) enddo return end C***************************************************** C @(#)ranve.f 1.1 modified on 12/7/92 SUBROUTINE RANVE(VR,V,Q1,Q2) C C ASSIGNS TO V(3) A RANDOM VECTOR OF LENGTH VR DIMENSION V(3) CT=Q1 PHI=Q2 IF(abs(CT).GT.1.0)CT=2.*ranf()-1. IF(PHI.LT.0..or.phi.gt.6.28318)PHI=6.28318*ranf() ST=SQRT(1.-CT**2) V(1)=VR*ST*COS(PHI) V(2)=VR*ST*SIN(PHI) V(3)=VR*CT RETURN END C******************************************************** C @(#)interpolation.f 1.1 modified on 12/7/92 C Interpolate in 2 dimensions. C C Do a two dimensional interpolation that will work for N dimensions C with N greater than NMax. This uses the save interpolation routine C to get the best value of the function. I would recomend copying this C routine when you want to use it and substituting the best interpolation C routine for the one dimensional interpolations. subroutine interpol2(x1a,x2a,ya,d1,d2,x1,x2,y,dy,int1,int2) implicit none integer i1 real r1 C**** The interpolation functions for the x1 and x2 directions. These should C be either polint or ratint depending on the function. external int1, int2 C**** The dimensions of the arrays. integer d1, d2 C**** The independent coordinates of the function. real x1a(d1), x2a(d2) C**** The dependent values of the function. real ya(d1,d2) C**** The x1 and x2 values to find the function at. real x1, x2 C**** The estimated function at the function value. real y C**** The estimated error at the function value. real dy C**** The maximum number of secondary interpolation values integer MaxD2 parameter (MaxD2=3) C**** The intermediate interpolations values real xt(MaxD2), yt(MaxD2) C**** The number of temporary interpolation values to use. integer NumD2 C**** The offset from the first x2a value to the first one to be used. integer OffD2 C**** Find the index of the nearest x2a to the requested x2. NumD2 = d2 OffD2 = 0 if (NumD2.gt.MaxD2) then NumD2 = MaxD2 call HUNT(X2A,D2,X2,i1) OffD2 = i1 - NumD2/2 if (OffD2+NumD2-1.gt.D2) OffD2 = D2 - NumD2 + 1 if (OffD2.lt.1) OffD2 = 1 OffD2 = OffD2 - 1 endif C**** Fill xt with x2a values that bracket x2 and interpolate to find C function values for the x2a values at x1. dy = 0.0 do i1 = 1, NumD2 xt(i1) = x2a(i1+OffD2) call Int1(x1a,ya(1,i1+OffD2),D1,x1,yt(i1),r1) dy = dy + r1*r1 enddo C**** Use the interpolated yt values for x2a values to interpolate the final C y values. call int2(xt,yt,NumD2,x2,y,r1) dy = dy + r1*r1 dy = sqrt(dy) return end C From Numerical Recipes chap 3.1 c Polynomial interpolation and exterpolation: c Given arrays XA and YA each of length N and given a value c of X, this routine returns a value of Y, and an error estimate c DY. If P(x) is the polynomial of degree n-1 such that c P(XA) = YA, then the returned value Y = P(X) C 2/8/89 - Routine modified to handle arrays N larger than NMAX C SUBROUTINE POLINT(XA,YA,N,X,Y,DY) PARAMETER (NMAX=5, HUGE = 1.E+25) REAL XA(*), YA(*), C(NMAX), D(NMAX) integer NElem NTop = N NOff = 0 if (NTop.gt.NMax) then NTop = NMax C Binary search for the element closest to X. call HUNT(XA,N,X,NElem) NBase = NElem - NTop/2 if (NBase+NTop-1.gt.N) NBase = N - NTop + 1 if (NBase.lt.1) NBase = 1 NOff = NBase - 1 endif NS = 1 DIF = ABS(X-XA(1+NOff)) DO I = 1,NTop DIFT = ABS(X-XA(I+NOff)) IF (DIFT.LT.DIF) THEN NS = I DIF = DIFT ENDIF C(I) = YA(I+NOff) D(I) = YA(I+NOff) ENDDO Y = YA(NS+NOff) NS = NS - 1 DO M = 1, NTop-1 DO I = 1, NTop-M HO = XA(I+NOff) - X HP = XA(I+M+NOff) - X W = C(I+1) - D(I) DEN = HO - HP IF (DEN.EQ.0) THEN TYPE *, $ 'POLINT: The input points are to close together:', $ HO, HP Y = 0.0 DY = HUGE GOTO 99 ENDIF DEN = W/DEN C(I) = HO*DEN D(I) = HP*DEN ENDDO IF (2*NS.LT.NTop-M) THEN DY = C(NS+1) ELSE DY = D(NS) NS = NS - 1 ENDIF Y = Y + DY ENDDO 99 continue if (DY.gt.0.05*Y) then C type *, 'POLINT: Big Error',X,Y,' +-',DY C do i1 = 1,ntop C type *, x,xa(i1+noff),ya(i1+noff),y C enddo endif RETURN END c Routine take from Numerical Recipes Chap 3.2 c c Rational exterpolation routine: c Given arrays XA and YA each of length N and given a value X, this c routine returns a value of Y and accuracy estimate DY. The value c returned is that of a tiagonal rational function, evaluated at X c which passes through the N points (XA, YA). C 2/8/89 - modified to handle arrays N larger than nmax. C SUBROUTINE RATINT(XA,YA,N,X,Y,DY) PARAMETER (NMAX=5, TINY=1.E-25, HUGE = 1.E+25) REAL XA(*), YA(*), X, Y, DY REAL C(NMAX), D(NMAX) INTEGER N integer NElem NTop = N NOff = 0 if (NTop.gt.NMax) then NTop = NMax C Binary search for the element closest to X. call HUNT(XA,N,X,NElem) NBase = NElem - NTop/2 if (NBase+NTop-1.gt.N) NBase = N - NTop + 1 if (NBase.lt.1) NBase = 1 NOff = NBase - 1 endif NS = 1 HH = ABS(X-XA(1+NOff)) DO I = 1,NTop H = ABS(X-XA(I+NOff)) IF (H.EQ.0) THEN Y = YA(I+NOff) DY = 0.0 GOTO 99 ELSE IF (H.LT.HH) THEN NS = I HH = H ENDIF C(I) = YA(I+NOff) D(I) = YA(I+NOff) + TINY ENDDO Y = YA(NS+NOff) NS = NS - 1 DO M = 1, NTop-1 DO I = 1, NTop-M W = C(I+1) - D(I) H = XA(I+M+NOff) - X T = (XA(I+NOff) - X)*D(I)/H DD = T-C(I+1) IF (DD.EQ.0) THEN TYPE *, 'RATINT: There is a pole at X =', X Y = 0 DY = HUGE GOTO 99 ENDIF DD = W/DD D(I) = C(I*1)*DD C(I) = T*DD ENDDO IF (2*NS.LT.NTop-M) THEN DY = C(NS+1) ELSE DY = D(NS) NS = NS - 1 ENDIF Y = Y + DY ENDDO 99 CONTINUE if (DY.gt.0.05*Y) then C type *, 'RATINT: Big Error',X,Y,' +-',DY C do i1 = 1,ntop C type *, x,xa(i1+noff),ya(i1+noff),y C enddo endif RETURN END C From Numerical Recipes Chap 3.4 C Given an array XX of length N and a given value X, returns a value C such that X is between XX(J) and XX(J+1). XX must be monotonic. C J = 0 or N is returned if X is out of range. Jon input is C taken as the initial guess for J on output subroutine HUNT(XX,N,X,JL) integer N real XX(N), X C**** THE UPPER AND LOWER BOUND ON THE POSITION BEING SEARCHED FOR. integer jl, ju C**** The last position found. save j integer j data j /0/ C**** flag assending or decending order. LOGICAL ASCND C**** The increment for the lower and upper bound integer inc C**** The middle value for j. integer JM ASCND = XX(N).gt.XX(1) C**** Is the first guess useful? if (j.le.0.or.j.gt.n) then JL = 0 JU = N+1 goto 10 else jl = j endif INC = 1 if (X.ge.XX(JL).eqv.ASCND) then 1 JU = JL + INC if (JU.gt.N) then JU = N + 1 else if (X.ge.XX(JU).eqv.ASCND) then JL = JU INC = INC + INC goto 1 endif else ju = jl 2 jl = ju - inc if (JL.lt.1) then JL = 0 else if (X.lt.XX(JL).eqv.ASCND) then JU = JL INC = INC + INC goto 2 endif endif 10 if (JU-JL.gt.1) then JM = (JU+JL)/2 if (ASCND.eqv.(X.gt.XX(JM))) then JL = JM else JU = JM endif goto 10 endif j = JL return end C************************************************ C C Perform a linear interpolation on the values in XA, and YA. C XA, YA must be ordered. C SUBROUTINE LININT(XA,YA,N,X,Y,DY) PARAMETER (NMAX=10, TINY=1.E-25, HUGE = 1.E+25) REAL XA(N), YA(N), X, Y, DY integer NElem call HUNT(XA,N,X,NElem) if (NElem.ge.N) NElem = N - 1 if (NElem.lt.1) NElem = 1 DYDX = (YA(NElem+1) - YA(NElem)) / (XA(NElem+1) - XA(NElem)) DX = X - XA(NElem) + Tiny Y = YA(NElem) + DYDX*DX DY = DYDX*DX return end C****************************** C Convert the probablity distribution Y into a cumulative probablity C distribution YI. Y does not need to be normalized. YI will be C normalized correctly. SUBROUTINE CFD(YI,N,Y) INTEGER N REAL Y(N),YI(N+1) CF=0.0 DO I=1,N CF=CF+Y(I) YI(I+1)=CF ENDDO YI(1)=0.0 cf=1./cf DO I=2,n+1 YI(I)=YI(I)*cf ENDDO RETURN END C***************************************************** C @(#)prbin.f 1.1 modified on 12/7/92 FUNCTION PRBIN(YI,N,XMIN,XSTEP,YR) DIMENSION YI(1) MIN=1 MAX=N+1 GO TO 3 1 IF((MAX-MIN).EQ.1) GO TO 5 IF(YI(M)-YR) 2,6,4 2 MIN=M 3 M=(MAX+MIN)/2 GO TO 1 4 MAX=M GO TO 3 5 M=MIN FR=(YR-YI(MIN))/(YI(MAX)-YI(MIN)) PRBIN=XMIN+XSTEP*(FLOAT(MIN)+FR-1.0) RETURN 6 FR=0.0 PRBIN=XMIN+XSTEP*(FLOAT(M)-1.0) RETURN END C******************************************* C @(#)cone.f 1.1 modified on 12/7/92 SUBROUTINE CONE(SAa,SB,CTH,PHI) C ASSIGNS TO SB(3) A VECTOR FORMING AN ANGLE ACOS(CTH) WITH SA(3) C IF PHI<0 THEN A RANDOM PHI IS TAKEN (FROM CONICAL SURFACE) C DIMENSION SA(4),saa(3),SB(3),D(3,3),RK(3) r=0. do i=1,3 r=r+saa(i)**2 end do r=sqrt(r) do i=1,3 sa(i)=saa(i)/r end do CALL RANVE(1.,RK,CTH,PHI) R=SQRT(SA(1)**2+SA(2)**2) IF(R.Eq.0.)then SB(1)=RK(1) sb(2)=rk(2) sb(3)=rk(3)*sa(3) else D(1,1)=SA(1)*SA(3)/R D(1,2)=SA(2)*SA(3)/R D(1,3)=-R D(2,1)=-SA(2)/R D(2,2)=SA(1)/R D(2,3)=0. D(3,1)=SA(1) D(3,2)=SA(2) D(3,3)=SA(3) DO 1 I=1,3 SB(I)=0. DO 1 J=1,3 1 SB(I)=SB(I)+D(J,I)*RK(J) end if 10 RETURN END C***************************************************** C @(#)lorentz.f 1.1 modified on 12/7/92 subroutine lorentz(pi,pf,beta) dimension pi(4),pf(4),beta(3) betamag=0. do i=1,3 betamag=betamag+beta(i)**2 end do if(betamag.le.0.) return betamag=sqrt(betamag) if(betamag.ge.1.) betamag=.9999 gamma=1./sqrt(1.-betamag**2) pperp=0. if(betamag.eq.0.) then do i=1,4 pf(i)=pi(i) end do return end if do i=1,3 pperp=pperp+pi(i)*beta(i)/betamag end do do i=1,3 pf(i)=(gamma*pperp-gamma*betamag*pi(4))*beta(i)/betamag+ 2 (pi(i)-pperp*beta(i)/betamag) end do pf(4)=gamma*pi(4)-gamma*betamag*pperp return end C**************************************************** C @(#)lloren.f 1.1 modified on 12/7/92 SUBROUTINE lLOREN(PPRIME,EPRIME,P,E,BETA) DIMENSION PPRIME(3),P(3),BETA(3) BSQ=0. BDP=0. DO 10 I=1,3 BSQ=BSQ+BETA(I)**2 10 BDP=BDP+BETA(I)*P(I) IF(BSQ) 30,30,15 15 GAM=SQRT(1./(1.-BSQ)) SHIFT=(GAM-1.)*BDP/BSQ -GAM*E DO 20 I=1,3 20 PPRIME(I)=P(I)+SHIFT*BETA(I) EPRIME=GAM*(E-BDP) RETURN 30 DO 40 I=1,3 40 PPRIME(I)=P(I) EPRIME=E RETURN END C*********************************************** C @(#)loren.f 1.1 modified on 12/7/92 SUBROUTINE LOREN(C,P,Q,N,K) DIMENSION C(4),P(4,N),Q(4,N) S=0. DO 1 I=1,3 1 S=S+C(I)**2 IF(S.EQ.0.)GO TO 100 S=C(4)**2-S W=SQRT(S) Z=W+C(4) DO 20 J=1,N D=C(4)*P(4,J) DO 21 I=1,3 21 D=D-K*C(I)*P(I,J) E=D/W A=(P(4,J)+E)/Z Q(4,J)=E DO 22 I=1,3 22 Q(I,J)=P(I,J)-A*K*C(I) 20 CONTINUE GO TO 1000 100 DO 101 J=1,N DO 101 I=1,4 101 Q(I,J)=P(I,J) GO TO 1000 1000 RETURN END C**************************************************** C @(#)twob.f 1.2 modified on 1/27/93 SUBROUTINE TWOB(AMD,AM1,AM2,PD,P1) DIMENSION PD(4),P1(4,2),PTEM(4,2) AMD2=AMD**2 A2M1=AM1**2 A2M2=AM2**2 E1=(AMD2+A2M1-A2M2)/AMD/2. E2=(AMD2+A2M2-A2M1)/AMD/2. P=E1**2-A2M1 IF(P.LT.0)THEN P=0 print 100,AMD,AM1,AM2 100 FORMAT(/' error in TWOB - masses:',3G12.4) END IF P=SQRT(P) CALL RANVE(P,PTEM,2.,-1.) PTEM(4,1)=E1 PTEM(4,2)=E2 DO I=1,3 PTEM(I,2)=-PTEM(I,1) enddo CALL LOREN(PD,PTEM,P1,2,-1) RETURN END C******************************************************* C @(#)threeb.f 1.2 modified on 12/29/92 SUBROUTINE THREEB(AMD,AM1,AM2,AM3,PD,PP) common/number/ntry DIMENSION PD(4),PP(4,3) DIMENSION PD1(4),PTEM(4,3) A=(AM1+AM2)**2 A1=(AMD-AM3)**2-A B=(AM1+AM3)**2 B1=(AMD-AM2)**2-B SM1=AM1**2 SM2=AM2**2 SM3=AM3**2 SMD=AMD**2 2 CONTINUE ntry=ntry+1 AM12=A+ranf()*A1 AM13=B+ranf()*B1 DM12=1/(2*SQRT(AM12)) E1=(AM12+SM1-SM2)*DM12 E3=(SMD-AM12-SM3)*DM12 P1=SQRT(E1**2-SM1) P3=SQRT(E3**2-SM3) C=(E1+E3)**2 IF(AM13.LT.C-(P1+P3)**2)GOTO 2 IF(AM13.GT.C-(P1-P3)**2)GOTO 2 E1=(AM12+AM13-SM3-SM2)/2/AMD P1=SQRT(E1**2-SM1) EMR=AMD-E1 AMR=SQRT(EMR**2-P1**2) CALL RANVE(P1,PD1,2.,-1.) PD1(4)=EMR CALL TWOB(AMR,AM2,AM3,PD1,PTEM(1,2)) DO 1 I=1,3 1 PTEM(I,1)=-PD1(I) PTEM(4,1)=E1 CALL LOREN(PD,PTEM,PP,3,-1) RETURN END C******************************************************** C @(#)putgam.f 1.3 modified on 1/4/93 C SUBROUTINE PUTGAM(P) implicit none real P(4) integer k C**************************************************************** C @(#)momentum.cdk 1.2 modified on 1/12/93 C Include the information about the primary particles to be tracked. C**** The number of momenta. integer nmom C**** The particle 4 momentum. real ppr(4,20) C**** The particle masses real ppm(20) C**** The particle charges real chh(20) C**** The particle id. integer pid(20) COMMON/MOM/ nmom,ppr,ppm,chh,pid C***** Gamma ray information real ppg(4,20) integer ngam common /gamm1/ ppg,ngam NGAM = NGAM + 1 DO K=1,4 PPG(K,NGAM)=P(K) enddo RETURN END C************************************************** C @(#)entryexit.f 1.2 modified on 10/20/93 C c ****************** Entryexit c This takes a track reference point and direction and returns the c distance to the point the track enters and exits the detector. c not enter the region entry and exit are equal and arbitrary. subroutine Entryexit(Pos,Dir,PEntry,PExit) implicit none integer i1,i2 real r1,r2 c The Position of the track reference point. real Pos(3) c The direction of the track real Dir(3) C The distance along the track direction to the entray and exit point. real PEntry,PExit C**** The detector box. real Walls(6) data Walls / -1154.0, -843.0, -876.0, 1154.0, 843.0, 876.0 / C high and low corner of the box. real hc(3), lc(3) c the distances to be traveled real t1(6) integer it C**** The fiducial distance. real Fidd, PFidd Fidd = 0.0 goto 10 entry enterfidu(Pos,Dir,PFidd,PEntry,PExit) Fidd = Pfidd 10 continue c make sure HCorner and LCorner are in the right order do i1 = 1,3 LC(i1) = walls(i1) + Fidd HC(i1) = walls(i1+3) - Fidd enddo c find the distances to each plane crossing. it = 0 do i1=1,3 if (Dir(i1).ne.0) then it = it + 1 t1(it) = (LC(i1)-Pos(i1))/Dir(i1) it = it + 1 t1(it) = (HC(i1)-Pos(i1))/Dir(i1) endif enddo C**** Make sure the box was hit. (Direction non zero.) if (it.lt.2) then PExit = PEntry return endif c sort the distances do i1 = 1,it-1 do i2 = i1+1,it if (t1(i1).gt.t1(i2)) then r1 = t1(i1) t1(i1) = t1(i2) t1(i2) = r1 endif enddo enddo c the entry and exit points will be the middle two entries i1 = it/2 PEntry = t1(i1) PExit = t1(i1+1) c make sure the track entered the region. r1 = (PEntry + PExit)/ 2.0 do i1 = 1,3 r2 = Pos(i1) + r1*Dir(i1) if ((r2.lt.LC(i1)).or.(r2.gt.HC(i1))) then PEntry = PExit return endif enddo return end SUBROUTINE bspline(XA,YA,Y2A,N,X,Y) DIMENSION XA(N),YA(N),Y2A(N) KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) IF (H.EQ.0.) then type *, ' BAD XA, must be monotonic' stop endif A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. RETURN END SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2) PARAMETER (NMAX=100) DIMENSION X(N),Y(N),Y2(N),U(NMAX) IF (YP1.GT..99E30) THEN Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN.GT..99E30) THEN QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END C************************************************************** C %W% modified on %G% C SUBROUTINE KINEM(KEY,IFER,MODE,moda) C***************************************** C Generate the kinimatics of the initial particles. The results are placed C into the ppr and pgam arrays. C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC C**************************************************************** C @(#)momentum.cdk 1.2 modified on 1/12/93 C Include the information about the primary particles to be tracked. C**** The number of momenta. integer nmom C**** The particle 4 momentum. real ppr(4,20) C**** The particle masses real ppm(20) C**** The particle charges real chh(20) C**** The particle id. integer pid(20) COMMON/MOM/ nmom,ppr,ppm,chh,pid C***** Gamma ray information real ppg(4,20) integer ngam common /gamm1/ ppg,ngam C******************************************************************** C @(#)prdec.cdk 1.1 modified on 12/29/92 C pass information about the initail particls and the desired products. C**** The amass and info about the interaction from spot. real amint, amlep, amd, amm(5) integer np(5), npar, kt COMMON/PRDEC/AMINT,AMLEP,AMD,AMM,NP,NPAR,KT C*************************************************************** C ASSIGN A GENERATED VERTEX TO XINT(3) C*************************************************************** CALL GEOM(KEY) if (KEY.eq.0) then C*************************************************************** C MAKE A PROTON DECAY C*************************************************************** CALL PROTON(IFER,MODE,MODA) else if(key.eq.1) then C*************************************************************** C MAKE A NEUTRINO INTERACTION C*************************************************************** c old delta production model CALL DELT(IFER,MODE,MODA) else if (key.eq.2) then C*************************************************************** C MAKE A NEUTRINO INTERACTION C*************************************************************** c new single-pi model with the addition of multipips CALL MCPI(ifer,mode,moda) else if (key.eq.-1) then C*************************************************************** C PREPARE PPR ARRAY FOR SINGLE PARTICLE TRACKING C*************************************************************** nmom = 1 call idgeant(np(1),ppm(nmom),chh(nmom)) call ranve(1.0,ppr,2.0,-2.0) ppr(4,nmom) = np(2) if (np(3).gt.np(2)) then ppr(4,nmom) = ppr(4,nmom) + ranf()*(np(3)-ppr(4,nmom)) endif ppr(4,nmom) = ppr(4,nmom)/1000000.0 else if ((key.eq.-2).or.(key.eq.-3)) then C*************************************************************** C CALL STMIU FOR COSMIC MUONS C KEY = -2 Central Muon C KEY = -3 Atmospheric Muon C*************************************************************** CALL STMIU(key) endif RETURN END SUBROUTINE GEOM(KEY) C C GENERATES VERTEX XINT(3) IN THE FIDUCIAL VOLUME C C******************************************************************* C @(#)inter.cdk 1.1 modified on 12/29/92 C the interaction point. C real xint(3) common/inter/xint C**** The common block with the detector wall positions in it. real HT(3) COMMON/GEO/HT C************************************************************************ C @(#)fiducial.cdk 1.1 modified on 12/29/92 real FV common/fidv/FV data fv/0./ DO 1 I=1,3 1 XINT(I)=(2*ranf()-1)*(HT(I)-fv) RETURN END C******************************************* C @(#)onepionmodel.f 1.9 modified on 2/23/93 SUBROUTINE mcpi(IFER,MODE,MODA) implicit none real r1 integer i1, i2, i3 c this subroutine calculates the single pion production model c developed at Irvine c c *** np = 0 for gamma c 1 for pi negative c 2 for pi zero c 3 for pi plus c 4 for K negative c 5 for K zero stable c 6 for K plus c C**** Use fermi 1) momentum 1 or 0) not. integer ifer C**** The input control parameters from kinem that the onepi model C ignores. integer mode, moda c ***************************************************** c dimensions to handle inernuclear interactions c ***************************************************** c C**** Number of input particles for partnuc integer nin C**** The initial energy (ei), mass (ui), momentum (pi) and charge (chi) real ei(20), ui(20), pi(3,20), chi(20) C**** The final energy (ei), mass (ui), momentum (pi) and charge (chi) real eo(20), uo(20), po(3,20), cho(20) real stpt(3,20),start(3),decp(3) C**** Flags to pass information to partnuc C to use or not use fermi momentum. integer imode parameter (imode = 0) C**** integer icont C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC C******************************************************************* C @(#)inter.cdk 1.1 modified on 12/29/92 C the interaction point. C real xint(3) common/inter/xint C******************************************************************** C @(#)prdec.cdk 1.1 modified on 12/29/92 C pass information about the initail particls and the desired products. C**** The amass and info about the interaction from spot. real amint, amlep, amd, amm(5) integer np(5), npar, kt COMMON/PRDEC/AMINT,AMLEP,AMD,AMM,NP,NPAR,KT C**************************************************************** C @(#)momentum.cdk 1.2 modified on 1/12/93 C Include the information about the primary particles to be tracked. C**** The number of momenta. integer nmom C**** The particle 4 momentum. real ppr(4,20) C**** The particle masses real ppm(20) C**** The particle charges real chh(20) C**** The particle id. integer pid(20) COMMON/MOM/ nmom,ppr,ppm,chh,pid C***** Gamma ray information real ppg(4,20) integer ngam common /gamm1/ ppg,ngam C**** God Only Knows, But PARTNUC wants it. integer isetn parameter (isetn = 0) C**** The Atomic number of oxygen. REAL ANuc parameter (ANuc = 16.0) C**** The momentum of the final nucleus, lepton, and initial neutrino. real pnuc(4), plep(4) C**** The momentum of the initial neutrino and nucleus. real pnu(4), pfermi(4) real efnuc C**** The transverse momentum and center of mass energy of nu neucleon. real pt, tecm C**** The initial type of the neutrino. 1) electron 2) muon integer nutype C**** The initial type of the neucleon. 0) neutron 1) proton integer nuctype C**** The initial energy of the neutrino. real enu C**** The initial direction of the neutrino. real dirnu(3) C**** The final type of the lepton and nucleon. integer leptypef,nuctypef C**** The pions from the final state. integer npion real ppi(4,10) integer pitypef(10) C**** The total weight of this event from final state. real wttot, wtmax common/weights/wttot,wtmax C**** The weight of the subprocess from final state. real wt C**** The weight of the neutrino energy and direction. real wtnu C**** the indentity of the interaction from final state. integer ident C**** The pauli blocking energy in GeV. real ppauli parameter (ppauli=.26) C**** Flag if this has been called before. integer icalled data icalled/0/ real ranf C**** The number of trials to get a successful event. integer ntry C**** The number of good events found. integer nnn C**** Flag if event includes fermi momentum 0) has fermi -1) has no fermi. integer nc C**** Geant neutrino type. integer gnutype integer wtblock, pblock real pmag C**** The version of the nukin file if one already exists. character*6 nvers if(icalled.eq.0) then C**** If this is the first call then open up a nukin.b86 kinematics file. icalled=1 ntry = 0 nnn = 0 call hbook2(1000,'Neutrino Energy Spectrum Produced by UCIKine$', $ 50,0.0,5.0, $ 40,1.0,21.0,0.0) call hbook2(1001,'Lepton Energy Spectrum Produced by UCIKine$', $ 50,0.0,5.0, $ 4,1.0,3.0,0.0) call hbook1(1002,'Fermi Momentum of Nucleon$',50,0.0,.500,0.0) end if c***************************************************************** c INITIALIZATION c***************************************************************** do i1 = 1,20 ppm(i1)=0. do i2=1,4 ppr(i2,i1)=0. end do end do do i1=1,3 start(i1)=0. decp(i1)=0. end do 30 continue do i1=1,3 pfermi(i1)=0.0 end do pfermi(4)=amp C**** Generate an initial energy and direction for the neutrino. call nu_energy(enu,dirnu,nutype,wtnu) C type *, enu, nutype, wtnu c**** Decide if the struck nucleon was a neutron or a proton. if(ranf().gt..555555) then nuctype=0 !neutron else nuctype=1 !proton end if c**** Find initial fermi momentum of nucleon. nc=nuctype if(ifer.gt.0) then call fermid(nc,pfermi,efnuc) !nc=0 -> got fermi'd else nc=-1 end if C******************************************************* c decide on final state particles/momenta c c first calculate total energy in the cms c for now let the neutrino direction be along x c and remember the initital nucleon momentum is in pfermi pnu(1)=enu pnu(2)=0. pnu(3)=0. pnu(4)=enu C**** Find the transverse momentum and the center of mass energy of C the collision. pt=(pfermi(1)+enu)**2+pfermi(2)**2+pfermi(3)**2 tecm=sqrt((enu+pfermi(4))**2-(pnu(1)+pfermi(1))**2 $ -(pnu(2)+pfermi(2))**2-(pnu(3)+pfermi(3))**2) C******************************************************************* C Find the final state products from the interaction. call final_state(nuctype,pfermi,nutype,pnu,tecm,nuctypef,pnuc, $ leptypef,plep,npion,pitypef,ppi,wt,ident) ntry=ntry+1 if(wt.le.0.0)then wtblock=wtblock+1 go to 30 ! dont keep event end if if(npion.gt.3) then print *,' error npion too big',npion go to 30 end if c********************************************************** c do pauli blocking of the final nucleon c if(nc.eq.0)then pmag = sqrt(pnuc(1)**2 + pnuc(2)**2 + pnuc(3)**2) if(pmag.lt.ppauli)then pblock = pblock + 1 go to 30 !dont keep event end if end if C************************** C We have a good event. nnn=nnn+1 C**** Write the incoming neutrino to the output if (nutype.eq.1) then gnutype = -3 else if (nutype.eq.-1) then gnutype = -2 else if (nutype.eq.2) then gnutype = -6 else if (nutype.eq.-2) then gnutype = -5 else gnutype = 4 endif i2 = gnutype write(6,'('' $ INCOMING '',I3,F12.6,3F10.6,F10.1)'),I2,enu, $ dirnu(1),dirnu(2),dirnu(3), wtnu C**** Write the incoming nucleon to the output if (nuctype.eq.0) then i2 = 13 else if (nuctype.eq.1) then i2 = 14 else write(6,*) 'ONEPI: bad nuctype', nuctype endif call rotatenu(pnu,dirnu,pfermi,1) write(6,'('' $ INCOMING '',I3,F12.6,3F10.6,F10.1)'),I2, pfermi(4), $ pfermi(1), pfermi(2), pfermi(3), wtnu c************************************** c fill arrays with the final state do i1=1,4 ppr(i1,1)=plep(i1) ppr(i1,2)=pnuc(i1) do i2=1,npion if(abs(pitypef(i2)).le.1) then ppr(i1,2+i2)=ppi(i1,i2) else ppr(i1,2+i2)=0. end if end do end do nmom=2 if(abs(pitypef(1)).le.1) nmom=3 if(npion.ge.1) nmom=2+npion C**** Set the Lepton Type C LepTypeF: 1 - electron C 2 - muon C 0 - neutrino if(abs(leptypef).eq.1)then ppm(1)=amel chh(1)=sign(1,leptypef) else if(abs(leptypef).eq.2)then ppm(1)=ammu chh(1)=sign(1,leptypef) else ppm(1)=0. chh(1)=0. end if amlep=ppm(1) C**** Set the nucleon type C NucTypeF: 1 - Proton C 0 - Neutron if(nuctypef.eq.1) then ppm(2)=amp chh(2)=1. else if (nuctypef.eq.0) then ppm(2)=amn chh(2)=0. else write(6,*) 'ONEPI: invalid nuctypef',nuctypef stop end if C***** Set the Pion Types if(nmom.gt.2) then do i1=1,npion chh(2+i1)=pitypef(i1) if(pitypef(i1).eq.0)then ppm(2+i1)=ampo else ppm(2+i1)=ampc end if end do else chh(3)=0. ppm(3)=0. end if C**** Write the outgoing lepton information to output. if (leptypef.eq.0) then i2 = gnutype else if (abs(leptypef).eq.1) then if (leptypef.gt.0) then i2 = 2 else i2 = 3 endif else if (abs(leptypef).eq.2) then if (leptypef.gt.0) then i2 = 5 else i2 = 6 endif else write(6,*) 'ONEPI:: Bad leptypef', leptypef call exit(1) endif write(6,'('' $ OUTGOING '',I3,F12.6,3F10.6,I4)'),I2, plep(4), $ plep(1), plep(2), plep(3), ident C**** Write the outgoing hadron information to output if (NucTypeF.eq.0) then i2 = 13 else if (NucTypeF.eq.1) then i2 = 14 else write(6,*) 'ONEPI:: bad nuctypef',nuctypef call exit(1) endif write(6,'('' $ OUTGOING '',I3,F12.6,3F10.6,I4)'),I2, pnuc(4), $ pnuc(1), pnuc(2), pnuc(3), ident C**** Write the outgoing pions to the output. do i1 = 1,npion if (pitypef(i1).eq.0) then i2 = 7 else if (pitypef(i1).eq.1) then i2 = 8 else if (pitypef(i1).eq.-1) then i2 = 9 else write(6,*) 'ONEPI: bad pitypef', i1, pitypef(i1) call exit(1) endif write(6,'('' $ OUTGOING '',I3,F12.6,3F10.6,I4)'),I2, ppi(4,i1), $ ppi(1,i1), ppi(2,i1), ppi(3,i1), ident enddo ********************************************************************** c c do nuclear interactions of pion and nucleon c if(nc.eq.-1) go to 40 !no nuclear interactions of free proton icont=0 nin=0 do i1=2,nmom nin=nin+1 ei(nin)=ppr(4,i1) ui(nin)=ppm(i1) do i2=1,3 pi(i2,nin)=ppr(i2,i1) end do chi(nin)=chh(i1) end do if(nin.eq.0) go to 40 ! no pions or nucleons found call partnuc( isetn,anuc,stpt, $ ei,pi,ui,chi, ! input hadrons $ nin, ! number of initial hadrons. $ eo,po,uo,cho, ! output hadrons $ start, decp,imode,icont) i3=1 do i1=1,nin if(abs(uo(i1)).ge.1.e-4) then !pion not absorbed i3=i3+1 do i2=1,3 ppr(i2,i3)=po(i2,i1) end do ppr(4,i3)=eo(i1) ppm(i3)=uo(i1) chh(i3)=cho(i1) if (abs((abs(chh(i3))-1.0)).lt.0.5 $ .and.abs(ppm(i3)-AMPC).lt.0.01) ppm(i3) = AMPC if (abs(chh(i3)).lt.0.5 $ .and.abs(ppm(i3)-AMP).lt.0.05) ppm(i3) = AMN endif end do nmom=i3 do i1=nmom+1,20 ppm(i1)=0 chh(i1)=0 do i2=1,4 ppr(i2,i1)=0 end do end do 40 continue C**** Rotate the neutrino and produces into the detector coordinate system. call rotatenu(pnu,dirnu,ppr,nmom) C************************** C Histogram the results. C******************* Histogram the energy of the incoming neutrino. C r1 will contain the ident of the interaction. if (ident.lt.20) then r1 = ident else if (ident.lt.40) then r1 = 17 else if (ident.lt.60) then r1 = 18 else if (ident.lt.120) then r1 = 19 else r1 = 20 endif r1 = r1 + 0.25 if (abs(nutype).eq.2) r1 = r1 + 0.5 call hfill(1000,enu,r1,1.0) C******************* Histogram the energy of the outgoing lepton. C Ignore Neutrinos. C r1 0.25 -- e- C 0.75 -- e+ C 1.25 -- m- C 1.75 -- m+ if (ident.eq.15.or.ident.eq.16) then r1 = abs(LepTypeF) + 0.25 if (LepTypeF.gt.0) r1 = r1 + 0.5 call hfill(1001,plep(4),r1,1.0) endif C******************** Histogram the Fermi momentum. r1 = pfermi(1)**2 + pfermi(2)**2 + pfermi(3)**2 if (r1.gt.0.0) then r1 = sqrt(r1) else r1 = 0.0 endif call hfill(1002,r1,0.0,1.0) return END subroutine rotatenu(pnu,dirnu,ppr,nmom) c c routine to rotate ppr fron pnu to dirnu c C**** The momentum of the neutrino. real pnu(4) C**** The direction of the neutrino. That things should be rotated to. real dirnu(3) C**** The number of primary particles. integer nmom C**** The momenum of the primary particles in the detector. real ppr(4,nmom) C**** The rotation matrices ru and rv are rotations around the euler axis. real r(3,3), rv(3,3), ru(3,3) C**** the normalized direction of pnu() real dir(3) C**** The normalized direction of dirnu() real dir2(3) C**** A temporary variable used to rotate ppr real xr(3) C**** The basis vectors in the new coordinate system. real khat(3),ahat(3),bhat(3) rr=0. rrr=0. do i=1,3 rr=rr+pnu(i)**2 rrr=rrr+dirnu(i)**2 end do rr=sqrt(rr) rrr=sqrt(rrr) do i=1,3 dir(i)=pnu(i)/rr dir2(i)=dirnu(i)/rrr end do c c pick random khat c 10 call ranve(1.0,khat,2.0,-2.0) udotk=khat(1)*dir(1)+khat(2)*dir(2)+khat(3)*dir(3) if(udotk.eq.0) go to 10 vdotk=khat(1)*dir2(1)+khat(2)*dir2(2)+khat(3)*dir2(3) if(vdotk.eq.0) go to 10 call cross_product(khat,dir2,ahat) call normalize(ahat) call cross_product(dir2,ahat,bhat) call normalize(bhat) do i=1,3 rv(i,1)=dir2(i) rv(i,2)=ahat(i) rv(i,3)=bhat(i) end do call cross_product(khat,dir,ahat) call normalize(ahat) call cross_product(dir,ahat,bhat) call normalize(bhat) do i=1,3 ru(1,i)=dir(i) ru(2,i)=ahat(i) ru(3,i)=bhat(i) end do do i=1,3 do j=1,3 r(i,j)=0. do k=1,3 r(i,j)=r(i,j)+rv(i,k)*ru(k,j) end do end do end do do ii=1,nmom do i=1,3 xr(i)=0. do j=1,3 xr(i)=xr(i)+r(i,j)*ppr(j,ii) end do end do do i=1,3 ppr(i,ii)=xr(i) end do end do return end subroutine cross_product(a,b,c) dimension a(3),b(3),c(3) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) return end subroutine normalize(x) dimension x(3) r=0. do i=1,3 r=r+x(i)**2 end do r=sqrt(r) do i=1,3 x(i)=x(i)/r end do return end C****************************************************************** C Generate the final state of the neutrino and the nucleon. subroutine final_state(nuctype,pnuci,nutype,pnui,tecm,nuctypef, $ pnuc,leptypef,plep,npion, $ pitypef,ppi,wtf,identf) implicit none integer i1, i2 C**** nuctype: The nucleon type. 0) neutron 1) proton integer nuctype C**** pnuci: The initial momentum of the nucleon. real pnuci(4) C**** nutype: The neutrino type. 1) electron 2) muon +) nu -) nubar integer nutype C**** pnui: The initial momentum of the neutrion real pnui(4) C**** tecm: The center of mass energy. real tecm C**** nuctypef:The final nucleon type 0) neutron 1) proton integer nuctypef C**** pnuc: The final nucleon momentum. real pnuc(4) C**** leptypef:The final lepton type 0) neutrino +-1) electron +-2) muon integer leptypef C**** plep: The lepton momentum. real plep(4) C**** npion: The number of produced pions. integer npion, npion2p C**** pitypef: The type of the produced pions 0) pi0 +1) pi+ -1) pi- integer pitypef(*) C**** ppi: the momentum of the pions. real ppi(4,*) C**** wtf: The weight of the final state interaction returned to the C calling routine. If this is zero then no valid interaction was found. real wtf C**** identf: The identity of the interaction. C 1) nu p --> l- p pi+ 4) nubar p --> l+ p pi- C 2) nu n --> l- p pi0 5) nubar p --> l+ n pi0 C 3) nu n --> l- n pi+ 6) nubar n --> l+ n pi- C 7) nu p --> nu p pi0 11) nubar p --> nubar p pi0 C 8) nu n --> nu n pi0 12) nubar n --> nubar n pi0 C 9) nu p --> nu n pi+ 13) nubar p --> nubar n pi+ C 10) nu n --> nu p pi- 14) nubar n --> nubar p pi- C 15) nu n --> l- p 16) nubar p --> l+ n integer identf C**** The phase space weight for the two body charged current interactions. C 15, 16 real wt2 C**** The matrix element for the two body interactions. This is forced to C zero for the neutral current interactions since they produce no C visible interactions in imb3 real wtm2 C**** The phase space weight for 3 body charged current interactions C 1, 2, 3, 4, 5, 6 real wt3 C**** The phase space weight for 3 body neutral current interactions C 7, 8, 9, 10, 11, 12, 13, 14 real wt3nc C**** The interaction weights of the individual onepion modes. The mode for C each of these elements is identified nchanel(1..nfin,nutype,nuctype+1) C where nutype is 1) nubar, 2) nu and nuctype is 0) neutron, 1) proton. real wtt(5) C**** The interaction weight of the 2 pion charged current mode. real wt2pcc C**** The interaction weight of the 2 pion neutral current mode. real wt2pnc C**** The interaction weight of the 3 pion charged current mode. real wt3pcc C**** The interaction weight of the 3 pion neutral current mode. real wt3pnc C**** Dummy variables to hold the initial neutrino or nucleon. real pnuii(4),pnui2(4) real pnui3(4),pnuci3(4),pnucii(4),pnuci2(4) real pnuc2(4),pnuc3(4),plep2(4),plep3(4),ppi3(4) real pnui3nc(4),pnuci3nc(4),pnuc3nc(4),plep3nc(4),ppi3nc(4) integer nuctypef2pcc, leptypefpcc integer ident2pcc, nuctypef2pnc integer ident2pnc, npion3p, nuctypef3pcc, leptypef2pcc integer ident3pcc, nuctypef3pnc, leptypefpnc integer leptypef2pnc, leptypef3pcc, leptypef3pnc, ident3pnc C**** The number of times the total weight has been greater than wtmax. integer ngreat C**** The sum of the partial weights real wtsum C**** The random weight for this event real wtrand C**** The total and maximum weights passed to onepi for debuging. real wttot, wtmax common/weights/wttot,wtmax C************************************************** C Control for matvect, matform, matcoeff, and matdia c**** iel=1 for elastic (fills arrays .mnu) c iel=0 for pi production integer iel c**** ncc=-1 for isoscalar nc (fills arrays a..s) c ncc= 0 for isovector nc (fills arrays ...2) c ncc= 1 for isovector cc (fills arrays ...1) integer ncc common/elastic/iel,ncc C**** The momentum input vectors for matvect. real p1(4),p2(4),pe(4),pnub(4),dumm(9),s,ped,z common/vectors/p1,p2,pe,pnub common/matscal/dumm equivalence (s,dumm(1)),(ped,dumm(4)) C**** The proton mass plus the pion mass. real scut parameter (scut=(0.938+0.139)*1.078) C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC C**** Flag an anti neutrion or a neutrino. 1) anti-neutrino 2) neutrino integer intype C**** The number of possible final states given the initial state. C indexed by (anti-nu/nu,neutron/proton) integer nfin, ifinal integer nfinal(2,2) data nfinal/3,5, $ 5,3/ C**** The identity of the channels available for a particular particle C combination. integer ident integer nchanel(5,2,2) data nchanel/ $ 6,12,14, 0, 0, $ 2, 3, 8,10,15, $ 4, 5,11,13,16, $ 1, 7, 9, 0, 0 $ / C**** The type of the produced pion for various quasi-elastic reactions. integer iptype(16) data iptype/1,0,1,-1,0,-1,0,0,1,-1,0,0,1,-1,3,3/ C**** The type of interacton for twopi and threepi interactions. C 0) charge current interactions 1) neutral current interactions integer itype real pnuc2pcc(4),pnuc2pnc(4),plep2pcc(4),plep2pnc(4), $ ppi2pcc(4,10),ppi2pnc(4,10) integer pitypef2pcc(10),pitypef2pnc(10) real pnuc3pcc(4),plep3pcc(4),ppi3pcc(4,10) integer pitypef3pcc(10) real pnuc3pnc(4),plep3pnc(4),ppi3pnc(4,10) integer pitypef3pnc(10) real ranf, qmat C************************************************ C INITIALIZE FOR THIS EVENT. wtmax = 60.0 wtf=0.0 npion=-1 do i1=1,10 pitypef(i1)=10 end do do i1=1,5 wtt(i1)=0.0 end do wt2pcc = 0.0 wt2pnc = 0.0 wt3pcc = 0.0 wt3pnc = 0.0 C**** Find if this is an anti neutrino or a neutrino if(nutype.lt.0) then intype=1 else intype=2 endif C************************************************ C Find the weights for the zero and one pion final states. nfin=nfinal(intype,nuctype+1) c**************************************** c**** Do phase space for one pion C.C. pnuc2(4)=amp pnuc3(4)=amp ppi3(4)=ampc plep2(4)=amel plep3(4)=amel if(abs(nutype).eq.2) then plep2(4)=ammu plep3(4)=ammu end if c after calls to threebod and twobod pnui has been changed c so save original values do i1=1,4 pnuii(i1)=pnui(i1) pnucii(i1)=pnuci(i1) end do call threebod(pnucii,pnuii,tecm,pnuc3,plep3,ppi3,wt3) do i1=1,4 pnui3(i1)=pnuii(i1) pnuci3(i1)=pnucii(i1) pnuii(i1)=pnui(i1) pnucii(i1)=pnuci(i1) end do call twobod(pnucii,pnuii,tecm,pnuc2,plep2,wt2) do i1=1,4 pnui2(i1)=pnuii(i1) pnuci2(i1)=pnucii(i1) pnuii(i1)=pnui(i1) pnucii(i1)=pnuci(i1) end do c*************************************** c**** Do Phase Space for one pion N.C. pnuc3nc(4)=amp ppi3nc(4)=ampc plep3nc(4)=0. do i1=1,4 pnuii(i1) = pnui(i1) pnucii(i1) = pnuci(i1) enddo call threebod(pnucii,pnuii,tecm,pnuc3nc,plep3nc,ppi3nc,wt3nc) do i1=1,4 pnui3nc(i1)=pnuii(i1) pnuci3nc(i1)=pnucii(i1) end do do i2=1,nfin wtt(i2)=0. ident=nchanel(i2,intype,nuctype+1) C**** For Three Body Final States INELASTIC. if(ident.lt.15) then C**** Set iel for pion production. iel=0 if(ident.lt.7) then C**** For three body Charged Current interactions if(wt3.le.0.0) go to 100 C**** Set the input kinematics do i1=1,4 p1(i1)=pnuci3(i1) p2(i1)=pnuc3(i1) pe(i1)=plep3(i1) pnub(i1)=pnui3(i1) end do C**** Set ncc for isovector charged current. ncc=1 C**** Get the kinematic variables. call matvect if(s.lt.scut.or.ped.gt.0.)go to 100 call matform ! get form factors call matcoeff ! get coeffs call matdia(ident) !get a's and b's C**** Set z for neutrino or anti neutrino. 1) neutrino -1) anti-neutrino if(nutype.gt.0) then z = 1. else z = -1. endif C**** Find the weight for this interaction. wtt(i2)=wt3*qmat(z)*1.02e-3 !x10**-38 cm**2 else C**** For three body Neutral Current interactions if(wt3nc.le.0.0) go to 100 C**** Set the input kinematics do i1=1,4 p1(i1)=pnuci3nc(i1) p2(i1)=pnuc3nc(i1) pe(i1)=plep3nc(i1) pnub(i1)=pnui3nc(i1) end do call matvect if(s.lt.scut.or.ped.gt.0.)go to 100 C**** Set ncc for isovector neutral current. ncc=0 call matform call matcoeff C**** Set ncc for isoscalar neutral current. ncc=-1 call matform call matcoeff call matdia(ident) C**** Set z for neutrino or anti neutrino. 1) neutrino -1) anti-neutrino if(nutype.gt.0) then z = 1. else z = -1. endif wtt(i2)=wt3nc*qmat(z)*1.02e-3 !x10**-38 cm**2 end if else C**** For two body final states. ELASTIC if(wt2.le.0.0) go to 100 C**** Set the input kinematics do i1=1,4 p1(i1)=pnuci2(i1) p2(i1)=pnuc2(i1) pe(i1)=plep2(i1) pnub(i1)=pnui2(i1) end do C**** Set iel for no pion production. iel=1 C**** Set ncc for charged current isovector interactions ncc=1 call matvect call matform call matcoeff if(nutype.gt.0 .and. nuctype.eq.0) then C**** nu n --> l- p z=1. call matdia(15) wtm2=qmat(z) else if (nutype.lt.0 .and. nuctype.eq.1) then C**** nubar p --> l+ n z=-1. call matdia(16) wtm2=qmat(z) else C**** interactions that make no visible particle. C**** nu p --> nu p C**** nubar n --> nubar n wtm2 = 0.0 end if wtt(i2)=wt2*wtm2*.0631 ! x10**-38 cm**2 end if 100 continue end do c***************************** c Find the weights for the two pion final states npion2p=2 C**** Find the charged current weight. itype=0 call twopi(npion2p,itype,nutype,pnui,nuctype,pnuci, $ nuctypef2pcc,pnuc2pcc,leptypef2pcc,plep2pcc, $ pitypef2pcc,ppi2pcc,wt2pcc) if(wt2pcc.gt.0.) then ident2pcc=1 do i2=1,2 ident2pcc=ident2pcc+(pitypef2pcc(i2)+1)*3**(2-i2) end do ident2pcc=ident2pcc+nuctypef2pcc*9+20 else C***** In Original Version These are NOT SET TO ZERO. wt2pcc = 0.0 ident2pcc = 0 end if C**** Find the neutral current weight. itype=1 call twopi(npion2p,itype,nutype,pnui,nuctype,pnuci, $ nuctypef2pnc,pnuc2pnc,leptypef2pnc,plep2pnc, $ pitypef2pnc,ppi2pnc,wt2pnc) if(wt2pnc.gt.0.) then ident2pnc=1 do i2=1,2 ident2pnc=ident2pnc+(pitypef2pnc(i2)+1)*3**(2-i2) end do ident2pnc=ident2pnc+nuctypef2pnc*9+40 else wt2pnc = 0.0 ident2pnc = 0 end if C*************************************************** C Find the weights for the three pion final states. npion3p=3 C**** Find the weight for three pion charged current. itype=0 call threepi(npion3p,itype,nutype,pnui,nuctype,pnuci, $ nuctypef3pcc,pnuc3pcc,leptypef3pcc,plep3pcc, $ pitypef3pcc,ppi3pcc,wt3pcc) if(wt3pcc.gt.0.) then ident3pcc=1 do i2=1,3 ident3pcc=ident3pcc+(pitypef3pcc(i2)+1)*3**(3-i2) end do ident3pcc=ident3pcc+nuctypef3pcc*27 + 60 else wt3pcc = 0.0 ident3pcc = 0 end if C**** Find the weight for three pion neutral current. itype=1 call threepi(npion3p,itype,nutype,pnui,nuctype,pnuci, $ nuctypef3pnc,pnuc3pnc,leptypef3pnc,plep3pnc, $ pitypef3pnc,ppi3pnc,wt3pnc) if(wt3pnc.gt.0.) then ident3pnc=1 do i2=1,3 ident3pnc=ident3pnc+(pitypef3pnc(i2)+1)*3**(3-i2) end do ident3pnc=ident3pnc+nuctypef3pnc*27 + 120 else wt3pnc = 0.0 ident3pnc = 0 end if c******************************************************** c pick particular final state from nfinal choices C**** Find the total weight of having an interaction. wttot=0. do i1=1,nfin if(wtt(i1).lt.0.) wtt(i1) = 0.0 wttot=wttot+wtt(i1) end do wttot=wttot+wt2pcc+wt2pnc wttot=wttot+wt3pcc+wt3pnc C**** Generate a random weight for this event. wtrand = ranf()*wtmax C**** If wttot is larger than wtmax warn incase this is happening to much. if(wttot.gt.wtmax) then ngreat=ngreat+1 type *, 'FINAL_STATE:: wtmax no large enough, rescaling wttot.' type '('' number: '',I4,'' weight:'',F10.6, $ '' energy:'',2F7.3)', $ ngreat,wttot,pnui(4),pnuci(4) wtrand = wtrand*wttot/wtmax end if C**** If the random weight is greater that the total weight C of this event then return with no interaction. if(wtrand.gt.wttot) then wtf=0. return end if C******************************************************** C Find the final state interaction by summing weights until the sum C is greater than the random weight of this event. wtsum=0. C**** Check if this event is a quasi elastic or a one pion interaction. do i1=1,nfin wtsum=wtsum+wtt(i1) if(wtsum.gt.wtrand) then ifinal=i1 go to 500 end if end do C**** Check if this event is a two pion interaction. wtsum=wtsum+wt2pcc if(wtsum.gt.wtrand) go to 600 !two pion cc wtsum=wtsum+wt2pnc if(wtsum.gt.wtrand) go to 700 !two pion nc C**** Check if this event is a three pion interaction. wtsum=wtsum+wt3pcc if(wtsum.gt.wtrand) go to 6003 !three pion cc wtsum=wtsum+wt3pnc if(wtsum.gt.wtrand) go to 7003 !three pion nc C**** If this is reached then wtsum is less than wttot and we have a problem. type *, 'FINAL_STATE: wtsum less than random weight' type *, ' wtsum',wtsum,' wtrand',wtrand,' wttot',wttot wtf = 0.0 return 500 continue !single pion+QE wtf=wtt(ifinal) ident=nchanel(ifinal,intype,nuctype+1) identf=ident if(ident.lt.15) then ! final state single pi prod npion=1 pitypef(1)=iptype(ident) if(ident.lt.7) then ! C.C. leptypef=-nutype nuctypef=nuctype-pitypef(1)-sign(1,leptypef) do i1=1,4 ppi(i1,1)=ppi3(i1) pnuc(i1)=pnuc3(i1) pnuci(i1)=pnuci3(i1) plep(i1)=plep3(i1) pnui(i1)=pnui3(i1) end do if (wt3.gt.149.) then type *, 'FINAL_STATE:: wt3 is big',wt3,ident,nutype endif else ! N.C. leptypef=0 nuctypef=nuctype-pitypef(1) do i1=1,4 ppi(i1,1)=ppi3nc(i1) pnuc(i1)=pnuc3nc(i1) pnuci(i1)=pnuci3nc(i1) plep(i1)=plep3nc(i1) pnui(i1)=pnui3nc(i1) end do end if else npion=0 pitypef(1)=3 ! Q.E. so no pion in f.s. leptypef=-nutype nuctypef=abs(1.-nuctype) do i1=1,4 ppi(i1,1)=0. pnuc(i1)=pnuc2(i1) pnuci(i1)=pnuci2(i1) plep(i1)=plep2(i1) pnui(i1)=pnui2(i1) end do end if return 600 continue !double pi wtf=wt2pcc ident=ident2pcc identf=ident npion=2 do i1=1,npion pitypef(i1)=pitypef2pcc(i1) end do leptypef=leptypef2pcc nuctypef=nuctypef2pcc do i1=1,4 do i2=1,npion ppi(i1,i2)=ppi2pcc(i1,i2) end do pnuc(i1)=pnuc2pcc(i1) pnuci(i1)=pnuci(i1) plep(i1)=plep2pcc(i1) pnui(i1)=pnui(i1) end do return 700 continue wtf=wt2pnc ident=ident2pnc identf=ident npion=2 do i1=1,npion pitypef(i1)=pitypef2pnc(i1) end do leptypef=leptypef2pnc nuctypef=nuctypef2pnc do i1=1,4 do i2=1,npion ppi(i1,i2)=ppi2pnc(i1,i2) end do pnuc(i1)=pnuc2pnc(i1) pnuci(i1)=pnuci(i1) plep(i1)=plep2pnc(i1) pnui(i1)=pnui(i1) end do return 6003 continue !triple pi wtf=wt3pcc ident=ident3pcc identf=ident npion=3 do i1=1,npion pitypef(i1)=pitypef3pcc(i1) end do leptypef=leptypef3pcc nuctypef=nuctypef3pcc do i1=1,4 do i2=1,npion ppi(i1,i2)=ppi3pcc(i1,i2) end do pnuc(i1)=pnuc3pcc(i1) pnuci(i1)=pnuci(i1) plep(i1)=plep3pcc(i1) pnui(i1)=pnui(i1) end do return 7003 continue wtf=wt3pnc ident=ident3pnc identf=ident npion=3 do i1=1,npion pitypef(i1)=pitypef3pnc(i1) end do leptypef=leptypef3pnc nuctypef=nuctypef3pnc do i1=1,4 do i2=1,npion ppi(i1,i2)=ppi3pnc(i1,i2) end do pnuc(i1)=pnuc3pnc(i1) pnuci(i1)=pnuci(i1) plep(i1)=plep3pnc(i1) pnui(i1)=pnui(i1) end do return end C********************************************************** subroutine twobod(pnuci,pnui,tecm,pnuc,plep,wt) dimension pnuci(4),pnui(4),pnuc(4),plep(4),ppi(4) dimension ppiprime(4),pnucprime(4),beta(3),dir(3),pdum(4) C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC real mpn,m0,ksq,ksqmin,ksqmax data slope/2.26/ c c routine to generate three body phase space c in cms and lab frames c pass lepton mass input in plep(4), ditto for pion + nucleon c amlep=plep(4) amnuc=pnuc(4) c c calculate lepton 4-momenta c do i=1,3 beta(i)=(pnui(i)+pnuci(i))/(pnui(4)+pnuci(4)) end do call lorentz(pnui,pdum,beta) do i=1,4 pnui(i)=pdum(i) end do call lorentz(pnuci,pdum,beta) do i=1,4 pnuci(i)=pdum(i) end do tecm2=sqrt((pnui(4)+pnuci(4))**2-(pnui(1)+pnuci(1))**2 2 -(pnui(2)+pnuci(2))**2-(pnui(3)+pnuci(3))**2) tecm=tecm2 if(tecm.lt.(amlep+amnuc)) then wt=0. return end if elep=(tecm**2+amlep**2-amnuc**2)/(2.*tecm) plepm=sqrt(elep**2-amlep**2) c c now do k-squared importance sampling c k-squared usually <0. c pncm=pnui(4) pn=pnui(4) ksq=2.*pn*plepm ksqmin=amlep**2-2.*pn*elep-ksq ksqmax=ksqmin+2.*ksq if(ksqmax.gt.amlep**2) then wt=0. return end if r=ranf() if((slope*ksqmin).gt.-30.) then b=r*exp(slope*ksqmax)+(1.-r)*exp(slope*ksqmin) else b=(r)*exp(slope*ksqmax) end if ksq=log(b)/slope c c calculate k**2 weight c if(slope*ksqmin.gt.-30.) then zmien=exp(slope*ksqmax)-exp(slope*ksqmin) if(abs(zmien).lt.1e-10)then wt=0. return end if aa=slope/zmien else aa=slope/(exp(slope*ksqmax)) end if wtksq=aa*exp(slope*ksq) c c calculate total phase space weight c including the "1/V" factor c wt=3.1416/(wtksq*tecm*pncm)/(tecm*pncm) if(wt.gt.150.) then type *, 'TWOBOD:: wt2 is two big',wt wt=150. endif c c calculate nu-lep angle given k**2 c ct=(ksq-amlep**2+2.*pn*elep)/(2.*pn*plepm) if(abs(ct).gt.1.) then print *,'TWOBOD:: trouble in twobod abs(ct)>1' wt=0. return end if call cone(pnui,dir,ct,-1.) !get nu dir given k**2 pnucm=0. do i=1,3 plep(i)=dir(i)*plepm pnuc(i)=-plep(i) pnucm=pnucm+pnuc(i)**2 end do plep(4)=elep pnuc(4)=sqrt(pnucm+amnuc**2) c c transform everything into lab frame c do i=1,3 beta(i)=-beta(i) end do call lorentz(plep,pdum,beta) do i=1,4 plep(i)=pdum(i) end do call lorentz(pnuc,pdum,beta) do i=1,4 pnuc(I)=pdum(i) end do call lorentz(pnuci,pdum,beta) do i=1,4 pnuci(I)=pdum(i) end do call lorentz(pnui,pdum,beta) do i=1,4 pnui(I)=pdum(i) end do return end C************************************************************* subroutine threebod(pnuci,pnui,tecm,pnuc,plep,ppi,wt) dimension pnuci(4),pnui(4),pnuc(4),plep(4),ppi(4) dimension ppiprime(4),pnucprime(4),beta(3),dir(3),pdum(4) dimension beta2(3) C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC data gm/.25/,m0/1.35/ real mpn,m0,ksq,ksqmin,ksqmax,mpncut data slope/2.26/ data mpncut/1.5/ !cut to fit nu data c c routine to generate three body phase space c in lab frame for nu-nuc->lep nuc pi c pass lepton mass input in plep(4), ditto for pion + nucleon c amlep=plep(4) ampi=ppi(4) amnuc=pnuc(4) do i=1,3 beta(i)=(pnui(i)+pnuci(i))/(pnui(4)+pnuci(4)) end do call lorentz(pnui,pdum,beta) do i=1,4 pnui(i)=pdum(i) end do call lorentz(pnuci,pdum,beta) do i=1,4 pnuci(i)=pdum(i) end do tecm2=sqrt((pnui(4)+pnuci(4))**2-(pnui(1)+pnuci(1))**2 2 -(pnui(2)+pnuci(2))**2-(pnui(3)+pnuci(3))**2) tecm=tecm2 if(tecm.lt.(amlep+ampi+amnuc)) then wt=0. return end if c c importance sample on mpn c ampnmin=ampi+amnuc ampnmax=tecm-amlep if(ampnmin.ge.ampnmax) then wt=0. return end if thmin=atan(2.*(ampnmin-m0)/gm) thmax=atan(2.*(ampnmax-m0)/gm) anorm=gm/(2.*(thmax-thmin)) r=ranf() th=thmin+gm*r/(2.*anorm) mpn=m0+gm*tan(th)/2. if(mpn.gt.mpncut) then wt=0. return end if if(mpn.gt.tecm) then type *,'THREEBOD:: mpn>tecm',mpn,tecm wt=0. return end if c c calculate weight of mpn c wtmpn=anorm/((mpn-m0)**2+gm**2/4.) c c calculate lepton 4-momenta c elep=(tecm**2+amlep**2-mpn**2)/(2.*tecm) plepm=(elep**2-amlep**2) if(plepm.le..000025) then plepm=.000025 elep=sqrt(plepm+amlep**2) end if plepm=sqrt(plepm) c c now do k-squared importance sampling in lab c k-squared is usually < 0. c pncm=pnui(4) pn=pnui(4) ksq=2.*pn*plepm ksqmin=amlep**2-2.*pn*elep-ksq ksqmax=ksqmin+2.*ksq if(ksqmax.gt.amlep**2) then type *,'THREEBOD:: ksqmax too big' wt=0. return end if if(ksqmax.le.ksqmin) then type *,'THREEBOD:: ksqmax.le.ksqmin' wt=0. return end if r=ranf() if((slope*ksqmin).gt.-30.) then b=r*exp(slope*ksqmax)+(1.-r)*exp(slope*ksqmin) else b=(1.-r)*exp(slope*ksqmax) end if ksq=log(b)/slope c c calculate k**2 weight c if(slope*ksqmin.gt.-30.) then cc=(exp(slope*ksqmax)-exp(slope*ksqmin)) if(cc.le.0.) cc=1.e-10 aa=slope/cc else aa=slope/(exp(slope*ksqmax)) end if wtksq=aa*exp(slope*ksq) c c calculate nu-lep angle given k**2 c ct=(ksq-amlep**2+2.*pn*elep)/(2.*pn*plepm) if(abs(ct).gt.1.) then if (abs(ct).gt.1.01) then type *, 'THREEBOD::' type *, 'err in 3b abs(ct)>1.' type *, 'ct,ksq,amlep,pn,elep,plepm,ksqmin,ksqmax' type *, ct,ksq,amlep,pn,elep,plepm,ksqmin,ksqmax endif ct=sign(1.,ct) end if call cone(pnui,dir,ct,-1.) !get nu dir given k**2 c c get new lep dir c do i=1,3 plep(i)=dir(i)*plepm end do plep(4)=elep c c calculate pion and nucleon 4-momenta in the mpn frame c epi=(mpn**2+ampi**2-amnuc**2)/(2.*mpn) enuc=mpn-epi ppim=sqrt(epi**2-ampi**2) ppimprime=ppim call ranve(1.0,dir,2.,-2.) do i=1,3 ppiprime(i)=dir(i)*ppim pnucprime(i)=-ppiprime(i) end do ppiprime(4)=epi pnucprime(4)=enuc c c lorentz pion and nuc momenta into cms frame c do i=1,3 beta2(i)=plep(i)/sqrt(mpn**2+plepm**2) end do call lorentz(pnucprime,pnuc,beta2) call lorentz(ppiprime,ppi,beta2) c c calculate total phase space weight c wt=3.1416**2*ppimprime/(2.*wtksq*wtmpn*tecm*pncm)/(tecm*pncm) if(wt.gt.150.) then type '(''THREEBOD:: wt is two big.'',F8.1,F,F,F,F6.3,F6.3)', $ wt,ppimprime,wtksq,wtmpn,tecm,pncm wt=150. endif c c transform everything into lab frame c do i=1,3 beta(i)=-beta(i) end do call lorentz(plep,pdum,beta) do i=1,4 plep(i)=pdum(i) end do call lorentz(pnuc,pdum,beta) do i=1,4 pnuc(i)=pdum(i) end do call lorentz(ppi,pdum,beta) do i=1,4 ppi(i)=pdum(i) end do call lorentz(pnui,pdum,beta) do i=1,4 pnui(i)=pdum(i) end do call lorentz(pnuci,pdum,beta) do i=1,4 pnuci(i)=pdum(i) end do return end subroutine twopi(npion,itype,nutype,pnui,nuctype,pnuci, $ nuctypef,pnuc,leptypef,plep,pitypef,ppi,wt) c c call once per cc and nc event c itype=0 C.C. c itype=1 N.C. c dimension pnui(4),pnuci(4),pnuc(4),plep(4),ppi(4,10) integer pitypef(10) dimension wtcorrect(18,2) !correction (ident,itype+1) data wtcorrect/36*0./ if(icalled.eq.0) then icalled=1 end if 1234 continue if(npion.ne.2) then wt=0. return end if call multi2pi(npion,itype,nutype,pnui,nuctype,pnuci, $ nuctypef,pnuc,leptypef,plep,pitypef,ppi,wt) ident=1 do ijk=1,npion ident=ident + (pitypef(ijk)+1)*3**(npion-ijk) end do ident=ident+nuctypef*9 return end subroutine threepi(npion,itype,nutype,pnui,nuctype,pnuci, $ nuctypef,pnuc,leptypef,plep,pitypef,ppi,wt) c c call once per cc and nc event c itype=0 C.C. c itype=1 N.C. c dimension pnui(4),pnuci(4),pnuc(4),plep(4),ppi(4,10) integer pitypef(10) dimension wtcorrect(18,2) !correction (ident,itype+1) data wtcorrect/36*0./ if(icalled.eq.0) then icalled=1 end if 1234 continue if(npion.ne.3) then wt=0. return end if call multi3pi(npion,itype,nutype,pnui,nuctype,pnuci, $ nuctypef,pnuc,leptypef,plep,pitypef,ppi,wt) ident=1 do ijk=1,npion ident=ident+(pitypef(ijk)+1)*3**(npion-ijk) end do ident=ident+nuctypef*9 return end subroutine multi2pi(npion,itype,nutype,pnui,nuctype,pnuci $ ,nuctypef,pnuc,leptypef,plep,pitypef,ppi,wt) c c itype=0 C.C. c itype=1 N.C. c dimension pnui(4),pnuci(4),pnuc(4),plep(4),ppi(4,10) integer pitypef(10) dimension constant(2,2) ! const(itype+1,i) i=1 for nubar i=2 for nu data constant/1.,1.,1.,1./ ! to convert to abs. cross sec. c c initialization c if(icalled.eq.0) then icalled=1 constant(1,1)=19500. ! cc nu in 10**-38 cm**2 constant(2,1)=.34*19500. ! nc nu constant(1,2)=19500. ! cc nubar constant(2,2)=.39*19500. ! nc nubar end if if(npion.gt.10) then wt=0. return end if do i=1,10 pitypef(i)=10 end do nupquarki=1 if(nuctype.eq.1) nupquarki=2 initcharge=nuctype wtquark=0. c c pick initial meson type c if(itype.eq.0) then ! do C.C. wtquark=3-nupquarki ! weight from # quarks if(nutype.lt.0) wtquark=nupquarki leptypef=-nutype lepcharge=sign(1,leptypef) iquark=lepcharge !identify struck quark ss=(nuctype-.5)*nutype if(ss.lt.0) then ! nu-n or nubar-p r=ranf() if(r.lt..333333) then nuctypef=0 pitypef(1)=1 pitypef(2)=0 else if(r.gt..66666) then nuctypef=1 pitypef(1)=1 pitypef(2)=-1 else nuctypef=1 pitypef(1)=0 pitypef(2)=0 end if end if if(nutype.lt.0) then nuctypef=1-nuctypef pitypef(1)=-pitypef(1) pitypef(2)=-pitypef(2) end if else ! nu-p or nubar-n r=ranf() if(r.lt..5) then nuctypef=0 pitypef(1)=1 pitypef(2)=1 else nuctypef=1 pitypef(1)=1 pitypef(2)=0 end if if(nutype.lt.0) then nuctypef=1-nuctypef pitypef(1)=-pitypef(1) pitypef(2)=-pitypef(2) end if end if else ! do N.C.s leptypef=0 lepcharge=0 wtquark=1.5 r=ranf() if(r.lt..33333) then pitypef(1)=1 pitypef(2)=-1 else if(r.gt..666666) then pitypef(1)=0 pitypef(2)=0 else pitypef(1)=-1 pitypef(2)=0 end if end if if(nuctype.eq.1) then pitypef(1)=-pitypef(1) pitypef(2)=-pitypef(2) end if nuctypef=nuctype-pitypef(1)-pitypef(2) end if if(nuctype-nuctypef-lepcharge-pitypef(1)-pitypef(2).ne.0) then wt=0. print *,'error bad charges' print *,nuctype,nuctypef,leptypef,pitypef return end if call multipips(npion,itype,nutype,pnui,nuctype,pnuci $ ,nuctypef,pnuc,leptypef,plep,pitypef,ppi,wtps) qq=(pnui(4)-plep(4))**2-(pnui(1)-plep(1))**2 $ -(pnui(2)-plep(2))**2- $ (pnui(3)-plep(3))**2 wtff=1. ijk=1 if(nutype.lt.0) ijk=2 wt=constant(itype+1,ijk)*wtps*wtquark*wtff return end subroutine multi3pi(npion,itype,nutype,pnui,nuctype,pnuci, $ nuctypef,pnuc,leptypef,plep,pitypef,ppi,wt) c c itype=0 C.C. c itype=1 N.C. c dimension pnui(4),pnuci(4),pnuc(4),plep(4),ppi(4,10) integer pitypef(10) dimension constant(2,2) ! const(itype+1,i) i=1 for nubar i=2 for nu data constant/1.,1.,1.,1./ ! to convert to abs. cross sec. c c initialization c if(icalled.eq.0) then icalled=1 constant(1,1)=7.6e5 ! cc nu in 10**-38 cm**2 constant(2,1)=7.6e5*.34 ! nc nu constant(1,2)=7.6e5 ! cc nubar constant(2,2)=7.6e5*.39 ! nc nubar end if if(npion.gt.10) then wt=0. return end if do i=1,10 pitypef(i)=10 end do nupquarki=1 if(nuctype.eq.1) nupquarki=2 initcharge=nuctype wtquark=0. c c pick charges c if(itype.eq.0) then ! do C.C. wtquark=3-nupquarki ! weight from # quarks if(nutype.lt.0) wtquark=nupquarki leptypef=-nutype lepcharge=sign(1,leptypef) iquark=lepcharge !identify struck quark if(nutype.gt.0) then if(nuctype.gt.0) then !nu-p r=ranf() if(r.lt..33333333) then pitypef(1)=1 pitypef(2)=0 pitypef(3)=0 nuctypef=1 else if(r.gt..66666666) then pitypef(1)=1 pitypef(2)=1 pitypef(3)=0 nuctypef=0 else pitypef(1)=1 pitypef(2)=1 pitypef(3)=-1 nuctypef=1 end if end if else !nu-n r=ranf() if(r.lt..25) then pitypef(1)=0 pitypef(2)=0 pitypef(3)=0 nuctypef=1 end if if(r.gt..25.and.r.lt..5) then pitypef(1)=1 pitypef(2)=0 pitypef(3)=0 nuctypef=0 end if if(r.gt..5.and.r.lt..75) then pitypef(1)=1 pitypef(2)=0 pitypef(3)=-1 nuctypef=1 end if if(r.gt..75) then pitypef(1)=1 pitypef(2)=1 pitypef(3)=-1 nuctypef=0 end if end if else !nubar-p if(nuctype.gt.0) then r=ranf() if(r.lt..25) then pitypef(1)=0 pitypef(2)=0 pitypef(3)=0 nuctypef=0 end if if(r.gt..25.and.r.lt..5) then pitypef(1)=-1 pitypef(2)=0 pitypef(3)=0 nuctypef=1 end if if(r.gt..5.and.r.lt..75) then pitypef(1)=1 pitypef(2)=0 pitypef(3)=-1 nuctypef=0 end if if(r.gt..75) then pitypef(1)=1 pitypef(2)=-1 pitypef(3)=-1 nuctypef=1 end if else !nubar-n r=ranf() if(r.lt..3333333333) then pitypef(1)=-1 pitypef(2)=0 pitypef(3)=0 nuctypef=0 else if(r.gt..666666666) then pitypef(1)=-1 pitypef(2)=-1 pitypef(3)=0 nuctypef=1 else pitypef(1)=1 pitypef(2)=-1 pitypef(3)=-1 nuctypef=0 end if end if end if end if else ! do N.C.s leptypef=0 lepcharge=0 nuctypef=nuctype-ich wtquark=3./2. r=ranf() if(r.lt..5) then nuctypef=nuctype r=ranf() if(r.lt..5) then pitypef(1)=0 pitypef(2)=0 pitypef(3)=0 else pitypef(1)=1 pitypef(2)=-1 pitypef(3)=0 end if else nuctypef=1-abs(nuctype) r=ranf() if(r.gt..5) then pitypef(1)=1 pitypef(2)=-1 else pitypef(1)=0 pitypef(2)=0 end if pitypef(3)=nuctype-nuctypef-pitypef(1)-pitypef(2) end if end if 7654 continue if(abs(pitypef(npion)).gt.1) then wt=0. print *,'pi charge wrong',initcharge,lepcharge,nuctypef print *,(pitypef(i),i=1,npion) return end if call multipips(npion,itype,nutype,pnui,nuctype,pnuci $ ,nuctypef,pnuc,leptypef,plep,pitypef,ppi,wtps) ijk=1 if(nutype.lt.0) ijk=2 wt=constant(itype+1,ijk)*wtps*wtquark return end subroutine multipips(npion,itype,nutype,pnui,nuctype,pnuci 2 ,nuctypef,pnuc,leptypef,plep,pitypef,ppi,wt) c c itype=0 C.C. c itype=1 N.C. c dimension pnui(4),pnuci(4),pnuc(4),plep(4),ppi(4,10) integer pitypef(10) dimension pnucm(4),pnuccm(4),beta(3),dir(3),phad(4),pdum(4) dimension betacmtolab(3),constant(2),pdumdum(4),plepcm(4) data constant/1.,1./ common/genin/npgenbod,tecmgenbod,amass(18),kgenev common/genout/pcm(5,18),wtgenbod C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC data twopi/6.2832/ if(npion.gt.10) then wt=0. return end if c c next, pick scaling variables x and y c call pickx(x,nutype,itype) if(x.gt.1.) then print *,' x>1',x wt=0. return end if wtx=1. !should be like # of proper quarks in the struk nucleon??? enu=pnui(4) amnucleon=amp call picky(x,enu,amnucleon,y,nutype,itype) if(y.gt.1.) then print *,' y>1',y wt=0. return end if c c now generate plep c amlep=0. !remember partons are massless elep=enu*(1.-y) plepm=elep amlep=0. if(itype.eq.0) then amlep=amel if(abs(nutype).gt.1) amlep=ammu end if plep(4)=sqrt(elep**2+amlep**2) c anuscal=enu*y c qsqscal=x*2.*amp*anuscal ct=1.- amp*x*y/elep if(abs(ct).gt.1.) then print *,' ct>1',ct wt=0. return end if call cone(pnui,dir,ct,-1.) do i=1,3 plep(i)=dir(i)*plepm end do c c go to cms system c do i=1,3 beta(i)=(pnui(i)+pnuci(i))/(pnui(4)+pnuci(4)) betacmtolab(i)=-beta(i) end do call lorentz(pnui,pnucm,beta) call lorentz(pnuci,pnuccm,beta) call lorentz(plep,plepcm,beta) enucm=pnucm(4) elepcm=plepcm(4) tecm=pnucm(4)+pnuccm(4) s=tecm**2 c c calculate weights for x and y distributions c wtx=1. !should be like # of proper quarks in the struk nucleon??? wty=s*x c c calculate hadronic momenta c do i=1,4 phad(i)=pnui(i)+pnuci(i)-plep(i) end do c c go to hadronic cms c do i=1,3 beta(i)=phad(i)/phad(4) end do call lorentz(phad,pdum,beta) teh=pdum(4) kgenev=2 amass(1)=amp amtot=amass(1) do i=2,npion+1 amass(i)=ampc amtot=amtot+amass(i) end do if(teh.lt.amtot) then wt=0. return end if tecmgenbod=teh npgenbod=npion+1 call genbod wtps=wtgenbod*(twopi**4.) 2 /((twopi**3.)**float(npgenbod))/(tecm*enucm) c c get hadronic momenta back into lab c do i=1,3 beta(i)=-beta(i) end do do i=1,4 pnuc(i)=pcm(i,1) end do call lorentz(pnuc,pdum,beta) do i=1,4 pnuc(i)=pdum(i) end do do i=1,npion do j=1,4 pdumdum(j)=pcm(j,i+1) end do call lorentz(pdumdum,pdum,beta) do j=1,4 ppi(j,i)=pdum(j) end do end do wt=wtx*wty*wtps*constant(itype+1) return end subroutine pickx(x,nutype,itype) external fnubarcc,fnucc,fnubarnc,fnunc dimension probx(100,2,2) c c itype=0 C.C. c itype=1 N.C. c if(icalled.eq.0) then icalled=1 xlow=0. xhigh=1. call funpre(fnubarcc,probx(1,1,1),xlow,xhigh) call funpre(fnucc,probx(1,2,1),xlow,xhigh) call funpre(fnubarnc,probx(1,1,2),xlow,xhigh) call funpre(fnunc,probx(1,2,2),xlow,xhigh) end if inutype=2 !neutrino if(nutype.lt.0) inutype=1 ! anti-nu call funran(probx(1,inutype,itype+1),x) return end subroutine picky(x,enu,amnuc,y,nutype,itype) c c itype=0 C.C. c itype=1 N.C. c inutype=1. !neutrino if(nutype.lt.0) inutype=2 !anti-nu r=ranf() c ymin=0. chi=2.*enu/(amnuc*x) ymax=chi/(1.+chi) if(inutype.eq.1) then y=r*ymax else xxx= 1. - r*(1.-(1.-ymax)**3) y=1.-xxx**.333333333 end if return end real function fnubarcc(xx) real x, xx if (xx.le.0.0) then x = 1.0E-4 else if (xx.ge.1.0) then x = 1.0 - 1.0E-4 else x = xx endif fnubarcc = (3.9*2./3.-1.2)*x**.55 $ * (1.-x)**3.2 + 1.1*2./3.*(1.-x)**8. return end real function fnubarnc(xx) real x, xx if (xx.le.0.0) then x = 1.0E-4 else if (xx.ge.1.0) then x = 1.0 - 1.0E-4 else x = xx endif fnubarnc=2./3.*(.1464+.1874) $ *(3.9*x**.55*(1.-x)**3.2+1.1*(1.-x)**8.) $ -1./3.*(.1033+.1767)*3.6*x**.55*(1.-x)**3.2 return end real function fnucc(xx) real x, xx if (xx.le.0.0) then x = 1.0E-4 else if (xx.ge.1.0) then x = 1.0 - 1.0E-4 else x = xx endif fnucc = (3.9*2./3.+1.2)*x**.55 * (1.-x)**3.2 $ + 1.1*2./3.*(1.-x)**8. return end real function fnunc(xx) real x, xx if (xx.le.0.0) then x = 1.0E-4 else if (xx.ge.1.0) then x = 1.0 - 1.0E-4 else x = xx endif fnunc = 2./3.*(.1464+.1874) $ *(3.9*x**.55*(1.-x)**3.2+1.1*(1.-x)**8.) $ +1./3.*(.1033+.1767)*3.6*x**.55*(1.-x)**3.2 return end C************************************************** C @(#)fermid.f 1.3 modified on 1/12/93 SUBROUTINE FERMId(NC,P,EFNUC) C NUCLEON MOMENTUM INSIDE O16 NUCLEUS IS GENERATED FOLLOWING C THE FIT OF PHYS. REV. C20(1979)744 DIMENSION P(4),PROBP(20),PROBPI(21),A(6) DATA A/.06140,-.38597,1.7283,-.30007,-.05410,-.04974/ DATA ALFA/.59206/,BETA/1.5428/,PI/3.1416/,HC/197./ logical isw/.true./ C*************************************************** C @(#)masses.cdk 1.2 modified on 12/30/92 C The particle masses C PROTON - AMP C NEUTRON - AMN C LAMBDA - AML C SIGMA - AMS C ELECTRON - AMEL C MUON - AMMU C PI 0 - AMPO C PI +/- - AMPC C ETA = AMETA C OMEGA = AMOM C RHO = AMRHO C K 0 = AMKO C K +- = AMKC C K* 0 = AMKSO C K* +- = AMKSC C*** The leptons real AMEL, AMMU C*** The baryons real AMP, AMN, AML, AMS C*** The mesons real AMPO, AMPC, AMETA, AMOMG, AMRHO C*** The strange mesons. real AMKO, AMKC, AMKSO, AMKSC COMMON/MASSES/AMEL,AMMU, $ AMP, AMN, AML, AMS, $ AMPO, AMPC, AMETA, AMRHO, AMOMG, $ AMKO, AMKC, AMKSO, AMKSC C**** The mass of the nucleon. REAL NucMass IF(ISW)then isw=.false. DO I=1,20 Q=((I-1)*20+10)/HC PR=0 DO J=1,6 AN=ALFA*BETA**J PR=PR+A(J)*(AN/PI)**1.5*EXP(-AN*Q**2) END DO PROBP(I)=PR*Q**2 END DO CALL CFD(PROBPI,20,PROBP) end if if (NC.eq.1) then NucMass = amp else NucMass = amn endif IF(NC.EQ.1.AND.ranf().LT..2)THEN C**** If nucleon is proton then it may be in a hydrogen. DO I=1,3 P(I)=0.0 END DO p(4)= NucMass EFNUC = NucMass NC=-1 ELSE C**** This nucleon is in an oxygen so find fermi momentum. PR=ranf() Q=PRBIN(PROBPI,20,0.,20.,PR) Q=Q*1.e-3 CALL RANVE(Q,P,2.,-1.) P(4)=sqrt(NucMass**2+Q**2) EFNUC=SQRT(.931**2-Q**2) NC=0 END IF RETURN END C********************************************* C @(#)nu_energy.f 1.2 modified on 12/29/92 C Generate a neutrino energy from a spectrum. subroutine nu_energy(enu, nudir, nutype, wtnu) implicit none integer i1 real r1, r2 C**** The neutrino energy real enu C**** The direction of the neutrino real nudir(3) C**** The internal definition of the neutrino type is: C 1 - electron C 2 - muon C 3 - anti-electron C 4 - anti-muon C**** The external definition of the neutrino type is: C 1 - electron C 2 - muon C -1 - anti-electron C -2 - anti-muon integer nutype C**** The event weight real wtnu C**** The number of energy bins. integer ebins, ebin parameter (ebins=200) C**** The number of cosine bins. integer cbins, cbin parameter (cbins=10) real cmin, cmax C**** The average cosine of the cosine bins. parameter (cmin=-1.0, cmax=1.0) real cosine(cbins) real dcos parameter (dcos = (cmax-cmin)/cbins) C**** The lowest energy of the energy bin. real emax, emin parameter (emax=5.0, emin=0.050) real energy(ebins+1) C**** The flux of the energy bin. real flux(ebins,cbins,4) C**** Make the conversion table from probability to energy real prob(ebins,cbins,4) C**** Make the sum table from probability to energy integer sbins, sbin parameter (sbins=ebins*cbins*4) real sum(sbins+1) integer called data called/0/ real ranf, nu_flux if (called.eq.0) then called = 1 C**** Set up the energy binning. do ebin = 1, ebins+1 energy(ebin) = emax*exp(-emax)*exp(emax*ebin/ebins) enddo C**** Set up the cosine binning. do cbin = 1, cbins cosine(cbin) = cmin + dcos*cbin - 0.5*dcos enddo C**** Set the flux. do nutype = 1,4 do cbin = 1,cbins do ebin = 1,ebins flux(ebin,cbin,nutype) = $ nu_flux(energy(ebin),cosine(cbin),nutype) enddo enddo enddo C**** Make the probability table. sbin = 1 sum(sbin) = 0 r1 = 0.0 do nutype = 1,4 do cbin = 1,cbins do ebin = 1,ebins sbin = sbin + 1 prob(ebin,cbin,nutype) = flux(ebin,cbin,nutype) sum(sbin) = sum(sbin-1) $ + prob(ebin,cbin,nutype) $ *dcos*(energy(ebin+1)-energy(ebin)) enddo enddo enddo C**** Normalize everything do nutype = 1, 4 do cbin = 1, cbins do ebin = 1, ebins prob(ebin,cbin,nutype) = $ prob(ebin,cbin,nutype)/sum(sbins+1) enddo enddo enddo do sbin = 1,sbins sum(sbin) = sum(sbin)/sum(sbins+1) enddo sum(sbins+1) = 1.0 endif C**** Find the bin of the probability. r1 = ranf() sbin = sbins/2 i1 = sbin 100 continue if (i1.gt.1) i1 = i1 / 2 if (r1.lt.sum(sbin)) then if (sbin.gt.1) then sbin = sbin - i1 goto 100 endif else if (r1.gt.sum(sbin+1)) then if (sbin.lt.sbins) then sbin = sbin + i1 goto 100 endif endif C**** Find the ebin,cbin,nutype index of the prob. nutype = sbin/ebins/cbins sbin = sbin - ebins*cbins*nutype cbin = sbin/ebins ebin = sbin - ebins*cbin nutype = nutype + 1 cbin = cbin + 1 ebin = ebin + 1 C**** Do a linear interpolation to find the actual values of enu and dirnu. r1 = ranf() enu = energy(ebin) + r1*(energy(ebin+1)-energy(ebin)) r2 = ranf() wtnu = cosine(cbin) + (r2-0.5)*dcos call ranve(1.0,nudir,wtnu,-1.0) C**** Return the weight of this event in Number/GeV/Cos/NuType wtnu = flux(ebin,cbin,nutype) if (nutype.gt.2) nutype = 2-nutype return end real function nu_flux(energy,cosine,nutype) implicit none real energy, cosine integer nutype real fluxnaumov real fluxleekoh89 integer last_used data last_used /0/ integer Use_Flux data Use_Flux /2/ common /Flux_To_Use/ Use_Flux if (Use_Flux.eq.1) then if (Last_Used.ne.1) then Last_Used = 1 write(6,*) 'NEUTRINO FLUX: Bugaev and Naumov Solar Minimum Flux.' endif nu_flux = fluxnaumov(energy,cosine,nutype) else if (Use_Flux.eq.2) then if (Last_Used.ne.2) then Last_Used = 2 write(6,*) 'NEUTRINO FLUX: Lee and Koh 89 Solar Average Flux.' endif nu_flux = FluxLeeKoh89(energy,cosine,nutype) else if (Use_Flux.eq.3) then if (Last_Used.ne.3) then Last_Used = 3 write(6,*) 'NEUTRINO FLUX: Uniform Flux.' endif nu_flux = 1.0 endif return end C***************************************************** C @(#)matrixelements.f 1.2 modified on 2/3/94 C C This contains the routines MAT* for calculating matrix elements in the C uci onepi model. C C FUNCTION FOR CALCULATING A FOUR DIMENSIONAL C DOT PRODUCT (A METRIC CONTRACTION OF TWO FOUR VECTORS) C THE METRIC IS gmunu=(1,-1,-1,-1) C 10/2/84 KSG FUNCTION D4(A,B) DIMENSION A(4),B(4) C=0. DO I=1,3 C=C+A(I)*B(I) END DO D4=A(4)*B(4)-C RETURN END FUNCTION EPS(A,B,C,D) dimension a(4),b(4),c(4),d(4) dimension v(24) V(1)=-A(1)*B(2)*C(3)*D(4) V(2)= A(1)*B(2)*C(4)*D(3) V(3)= A(1)*B(3)*C(2)*D(4) V(4)=-A(1)*B(4)*C(2)*D(3) V(5)=-A(1)*B(3)*C(4)*D(2) V(6)= A(1)*B(4)*C(3)*D(2) V(7)= A(2)*B(1)*C(3)*D(4) V(8)=-A(2)*B(1)*C(4)*D(3) V(9)=-A(3)*B(1)*C(2)*D(4) V(10)= A(4)*B(1)*C(2)*D(3) V(11)= A(3)*B(1)*C(4)*D(2) V(12)=-A(4)*B(1)*C(3)*D(2) V(13)=-A(2)*B(3)*C(1)*D(4) V(14)= A(2)*B(4)*C(1)*D(3) V(15)= A(3)*B(2)*C(1)*D(4) V(16)=-A(4)*B(2)*C(1)*D(3) V(17)=-A(3)*B(4)*C(1)*D(2) V(18)= A(4)*B(3)*C(1)*D(2) V(19)= A(2)*B(3)*C(4)*D(1) V(20)=-A(2)*B(4)*C(3)*D(1) V(21)=-A(3)*B(2)*C(4)*D(1) V(22)= A(4)*B(2)*C(3)*D(1) V(23)= A(3)*B(4)*C(2)*D(1) V(24)=-A(4)*B(3)*C(2)*D(1) EPS=0 DO I=1,24 EPS=EPS+V(I) END DO RETURN END subroutine matcoeff c****************************************************** c iel=1 for elastic (fills arrays .mnu) * c iel=0 for pi production * c * c ncc=-1 for isoscalar nc (fills arrays a..s) * c ncc= 0 for isovector nc (fills arrays ...2) * c ncc= 1 for isovector cc (fills arrays ...1) * c****************************************************** implicit real m,k common/elastic/iel,ncc common/coeff/AMpi1(8),BMpi1(8),AMN1(8),BMN1(8),AMNN1(8),BMNN1(8), 1 AMP1(8), BMP1(8), AMs1(8),BMs1(8),AMDL1(8),BMDL1(8), 2 AMD1(8),BMD1(8), 3 AMpi2(8),BMpi2(8),AMN2(8),BMN2(8),AMNN2(8),BMNN2(8), 4 AMP2(8), BMP2(8), AMs2(8),BMs2(8),AMDL2(8),BMDL2(8), 5 AMD2(8),BMD2(8), amnu(8),bmnu(8), 6 amns(8),amnns(8),amss(8),amds(8) COMPLEX AMpi(6),BMpi(6),ampi1,bmpi1,ampi2,bmpi2 COMPLEX AMN(6) ,BMN(6) ,amn1 ,bmn1, amn2 ,bmn2 ,amns COMPLEX AMNN(6),BMNN(6),amnn1,bmnn1,amnn2,bmnn2,amnns COMPLEX AMP(6) ,BMP(6) ,amp1 ,bmp1, amp2 ,bmp2 COMPLEX AMs(6) ,BMs(6) ,ams1 ,bms1 ,ams2 ,bms2 ,amss COMPLEX AMDL(6),BMDL(6),amdl1,bmdl1,amdl2,bmdl2 COMPLEX AMD(6) ,BMD(6) ,amd1 ,bmd1 ,amd2 ,bmd2 ,amds COMPLEX AMNU ,BMNU complex c1,q1,r1,e1,s1 common/matkin/p1p2,p1pe,p2pe,p1pnub,p2pnub,pepnub common/matscal/s,t,u,kk,kq,p1k,p2k,p1q,p2q common/form/fpit,f1vt,f2vt,c3vt,c4vt,g1vt,g2vt,h3vt,h4vt,h5vt, 1 fat,c5at,d1at,g1at,h5at c######################################################################### parameter (PI=3.14159263) c parameter (gNNpi=13.637, g=2.11, fp=4.45, fs=.48, fD=1.56) parameter (gNNpi=13.405, g=2.11, fp=4.45, fs=.48, fD=1.56) parameter (RT2=1.414213562, rt3=1.732050808) parameter (C2=1.5/RT2, C3=RT2, C4=.5/RT2, C5=1/RT2) parameter (C6=2/3., C7=1/3.*RT2, C8=1./RT3, C9=RT2/rt3) parameter (C10=RT3/2., C11=1/3., C12=2/3., C13=RT2/3.) parameter (C14=1./6., C15=3./2.) parameter (M=.938, MDL=1.232, MS=1.505, MP=1.434) parameter (MD=1.514, mpi=.139, mmu=.106) parameter (m2=M*M, mmdl=M+MDL, mmd=M+MD, mpi2=mpi*mpi) parameter (mndl=MDL-M, mnd=MD-M) parameter (mdl6=1./(6.*MDL), md6=1./(6.*MD)) parameter (mdl23=1./(3.*MDL*MDL), md23=1./(3.*MD*MD)) parameter (m2mdl=2.*m*mdl, m2md=2.*m*md) dimension lr(4),gr(4),qr(4) c width of the resonances c n=1 - p33 c n=2 - p11 c n=3 - s11 c n=4 - d13 data lr/1,1,0,2/ data gr/0.114,0.200,0.100,0.130/ data qr/0.2254,0.390,0.446,0.450/ parameter (xsq=0.1225) complex gam GAM(n)=(0.,1.)/2.*Gr(n)*(QP/QR(n))**(2*Lr(n)+1) 1 *((QR(n)**2+xsq)/(QP**2+xsq))**Lr(n) c******************* if(iel.eq.1)then c1=(1.,0.) amnu(5)=c1*fat bmnu(5)=c1*f1vt bmnu(6)=-f2vt/m go to 1000 end if c************************************************************* ep=(s+mpi2-m2)/2./sqrt(s) xxxx=ep**2-mpi2 if(xxxx.gt.0.) then qp=sqrt(xxxx) else open(unit=91,status='new',form='formatted',name='spoty.err') write(91,*) ' error in matcoeff qp=0.' write(91,*) ' xx,ep,s,mpi2',xxxx,ep,s,mpi2 write(91,*) ' p1p2,p1pe,p2pe,p1pnub,p2pnub,pepnub, 2s,t,u,kk,kq,p1k,p2k,p1q,p2q' write(91,*) p1p2,p1pe,p2pe,p1pnub,p2pnub,pepnub, 2 s,t,u,kk,kq,p1k,p2k,p1q,p2q close(unit=91) qp=0.01 end if p1kq=p1q+kq p1kk=p1k+kk p12kk=s-m2 p14kk=4*p1k+kk p1p2k=p1k+p2k p1p22k=p1p2+p2k sc1=s-2*p1p22k sc2=sc1+m2 sc3=sc1-m2+m2md sc4=sc1-m2-m2mdl c******************** FFVt=F1Vt+2*F2Vt C1=(1.,0.)*RT2*gNNpi/(u-M2) q1=c1*F1Vt AMN(1)= q1 AMN(2)=-q1 e1=C1*F2Vt/M AMN(3)=-e1 AMN(4)= e1 AMN(5)= e1*(2*p2k-kk) AMN(6)=-C1*FFVt c AMN(7)= C1*FFVt q1=c1*FAt BMN(1)= q1 BMN(2)=-q1 BMN(5)=-q1*2*M BMN(6)=-q1 c BMN(7)= q1 c******************** C1=(1.,0.)*RT2*gNNpi/(s-M2) e1=C1*F1Vt AMNN(1)= e1 AMNN(2)= e1 r1=C1*F2Vt/M AMNN(3)=-r1 AMNN(4)=-r1 AMNN(5)= r1*p12kk AMNN(6)=-C1*FFVt c AMNN(7)= C1*FFVt q1=C1*FAt BMNN(1)=-q1 BMNN(2)=-q1 BMNN(5)= q1*2*M BMNN(6)= q1 c BMNN(7)=-q1 c*********************** C1=C15*fs/(s-(Ms-gam(3))**2) q1=C1*G1Vt AMs(1)= q1 AMs(2)= q1 e1=C1*G2Vt/M AMs(3)= e1 AMs(4)= e1 AMs(5)=-q1*(M+Ms)-e1*p12kk Ams(6)=-q1-e1*(Ms-M) c Ams(7)= q1+e1*(Ms-M) q1=C1*G1At BMs(1)=-q1 BMs(2)=-q1 BMs(5)=-q1*(Ms-M) BMs(6)= q1 c BMs(7)=-q1 c*********************** C1=3/(2*mpi)*fD/(s-(MD-gam(4))**2) q1=C1*H3Vt/M r1=C1*H4Vt/m2 s1=C1*H5Vt/m2 y2=kq-(p1kq*md23+md6*mnd)*kk y3=mmd*kq/2+(mnd*C14-md6*p1kq)*kk y4=mmd*kq/2-((mmd*md23-md6)*p1kq+md6*sc3+mnd*C14)*kk AMD(1)=q1*y2+r1*y3+s1*y4 AMD(2)=q1*(y2-p12kk)+r1*(y3-mmd*p1kk)+s1*(y4-mmd*p1k) y2=kq/2+md6*p1kq*mnd-C14*sc2 y3=y2-(mnd*md6+md23*p1kq)*kk y4=-m*p1kq*md23-md6*(s-m2)+c12*mnd AMD(3)=q1*y4+r1*y2+s1*y3 AMD(4)=q1*(y4-mnd)+r1*(y2-p1kk)+s1*(y3-p1k) y2=kq-2*p1kq*p1kk*md23 y3=-2*md6*p1kq*mnd+C11*sc2 AMD(5)=q1*(mnd*y2 1 +2*md6*(2*p1kq*p12kk-p1kk*sc2)-C12*mnd*p12kk) 2 +r1*y3*p1kk+s1*y3*p1k y3=2*md6*p1kq-C11*mnd AMD(6)=q1*(-y2 1 +2*md6*(mnd*p1kk-2*m*p1kq)+c12*(sc3-p1kq)) 2 +r1*p1kk*y3+s1*p1k*y3 c y2=-kq/2+md23*p1kq*p12kk c y3=md6*(sc3*p12kk-p1kq*p14kk) c AMD(7)=q1*(-kk*(md23*p1kq+md6*mnd) c 1 +2*md6*(2*mmd*p1kq+md*(3*kq-2*sc3))) c 2 +r1*(kq/2*mmd-md6*p1kq*kk+C14*mnd*kk) c 3 +s1*(mmd*y2+y3+c14*mnd*p14kk) c AMD(8)=q1*(md23*mnd*p1kq+md6*sc2) c 1 +r1*(kq/2+md6*mnd*p1kq-C14*sc2) c 2 +s1*(y2+md6*(mmd*(m2-p1p22k)+md*p12kk)) c*************************************************************** if(ncc.lt.0)go to 100 c*************************************************************** C1=c1*H5At y2=-md6*(s-m2)+md23*mnd*p1kq BMD(1)=C1*(C11*mnd+y2) BMD(2)=C1*(-C12*mnd+y2) BMD(3)=-C1*(md6*mnd+p1kq*md23) BMD(4)=BMD(3)+c1 BMD(5)=C1*(-C11*sc3+p1kq*mmd*2*md6) BMD(6)=C1*(p1kq*2*md6-C11*mnd) c BMD(7)=C1*((mnd*sc2-p12kk*m)*md23/2+mnd*c11) c BMD(8)=C1*(md23*(p1p22k-s)-mnd*md6) c******************* C1=(1.,0.)*RT2*gNNpi*fpit/(t-mpi2) AMpi(2)=-C1*2 c AMpi(7)=C1 c******************* C1=C15*fp*D1At/(s-(Mp-gam(2))**2) BMp(1)=-C1 BMp(2)=-C1 BMp(5)= C1*(M+Mp) BMp(6)= C1 c BMp(7)=-C1 c****************** C1=g/mpi/(s-(MDL-gam(1))**2) q1=C1*C3Vt/M r1=C1*C4Vt/m2 y2=kq-(p1kq*mdl23+mdl6*mmdl)*kk y3=mndl*kq/2+(mmdl*C14-mdl6*p1kq)*kk AMDL(1)=q1*y2+r1*y3 AMDL(2)=q1*(y2-p12kk)+r1*(y3-mndl*p1kk) y2=-m*mdl23*p1kq+mdl6*p12kk-C12*mmdl y3=-kq/2-mdl6*p1kq*mmdl+C14*sc2 AMDL(3)=q1*y2+r1*y3 AMDL(4)=q1*(mmdl+y2)+r1*(y3+p1kk) y3=kq-2*p1kq*p1kk*mdl23 AMDL(5)=q1*(-mmdl*y3-2*mdl6*(2*p1kq*p12kk-sc2*p1kk) 1 +C12*mmdl*p12kk)+ 2 r1*p1kk*(2*mdl6*p1kq*mmdl-C11*sc2) AMDL(6)=q1*(-y3+2*mdl6*mmdl*(p1kk-2*p1p22k)+ 1 4*mdl6*m*s-c12*m*(mmdl+mdl))+ 1 r1*p1kk*(2*mdl6*p1kq-C11*mmdl) c AMDL(7)=q1*(-kk*(mdl23*p1kq+mdl6*mmdl)+ c 1 2*mdl6*(2*mndl*p1kq+mdl*(3*kq-2*sc4)))+ c 2 r1*(kq/2*mndl-mdl6*p1kq*kk+C14*mmdl*kk) c AMDL(8)=q1*(-mdl23*mmdl*p1kq-mdl6*sc2)+ c 1 r1*(-kq/2-mdl6*mmdl*p1kq+C14*sc2) C1=c1*C5At y2=-mdl6*(s-m2)+p1kq*mdl23*mmdl BMDL(1)=C1*( C11*mmdl+y2) BMDL(2)=C1*(-C12*mmdl+y2) BMDL(3)=C1*(mdl6*mmdl+p1kq*mdl23) BMDL(4)=BMDL(3)-c1 BMDL(5)=C1*(C11*sc4-p1kq*mndl*2*mdl6) BMDL(6)=C1*(2*p1kq*mdl6-C11*mmdl) c BMDL(7)=C1*((mmdl*sc2+p12kk*m)*mdl23/2+mmdl*c11) c BMDL(8)=C1*(mdl23*(s-p1p22k)+mmdl*mdl6) c************************************************************ c************************************************************ 100 if(ncc)1,2,3 1 continue do i=1,6 amns(i)=amn(i) amnns(i)=amnn(i) amss(i)=ams(i) amds(i)=amd(i) end do go to 1000 2 continue do i=1,6 ampi2(i)=ampi(i) bmpi2(i)=bmpi(i) amn2(i) =amn(i) bmn2(i) =bmn(i) amnn2(i)=amnn(i) bmnn2(i)=bmnn(i) amp2(i) =amp(i) bmp2(i) =bmp(i) ams2(i) =ams(i) bms2(i) =bms(i) amdl2(i)=amdl(i) bmdl2(i)=bmdl(i) amd2(i) =amd(i) bmd2(i) =bmd(i) end do go to 1000 3 continue do i=1,6 ampi1(i)=ampi(i) bmpi1(i)=bmpi(i) amn1(i) =amn(i) bmn1(i) =bmn(i) amnn1(i)=amnn(i) bmnn1(i)=bmnn(i) amp1(i) =amp(i) bmp1(i) =bmp(i) ams1(i) =ams(i) bms1(i) =bms(i) amdl1(i)=amdl(i) bmdl1(i)=bmdl(i) amd1(i) =amd(i) bmd1(i) =bmd(i) end do 1000 return end subroutine matdia(id) c****************************************************** c iel=1 for elastic (fills arrays .mnu) * c iel=0 for pi production * c * c ncc=-1 for isoscalar nc (fills arrays a..s) * c ncc= 0 for isovector nc (fills arrays ...2) * c ncc= 1 for isovector cc (fills arrays ...1) * c****************************************************** common/elastic/iel,ncc common/coeff/AMpi1(8),BMpi1(8),AMN1(8),BMN1(8),AMNN1(8),BMNN1(8), 1 AMP1(8), BMP1(8), AMs1(8),BMs1(8),AMDL1(8),BMDL1(8), 2 AMD1(8),BMD1(8), 3