PROGRAM MUON c --- generate muon decay for a muon beam c --- and store 4-vectors in an ntuple c --- nue -> numu variant implicit none c --- common /ev/ stores output event record integer kk double precision rescale parameter (rescale = 10**24) ! factor to make cross sections reasonable real event real ccgen real throwwt real zdecay real cckin(25) real pe_lab(4), pe_cms(4) real pnue_lab(4), pnue_cms(4), r_nue, sigma_nue real pnum_lab(4), pnum_cms(4), r_num, sigma_num real emu_det, ee_det real emax_nue, emax_num real nueconst ! 2 GF^2 me/pi c real emass parameter (nueconst=17.2E-42) real sin2tw real eynum,eynue real ncenuml,ncenuel real ncenumr,ncenuer real ccenuml,ccenuel real ccenumr,ccenuer integer icycle real factm,facte real snum,snue integer ityp common /ev/ event, zdecay, pe_lab, pe_cms, * pnue_lab, pnue_cms, r_nue, sigma_nue, * pnum_lab, pnum_cms, r_num, sigma_num, * emu_det, ee_det, emax_nue, emax_num c ********************* added for lepto ***************** real cut(14),parl(30),x,y,w2,q2,u,p(4000,5),v(4000,5) integer lst(40),n,k(4000,5) common /leptou/ cut(14),lst(40),parl(30),x,y,w2,q2,u common /lujets/ n,k(4000,5),p(4000,5),v(4000,5) integer lepin, inter, interaction, nucleon real ptot, nuthetax, nuthetay real targa, targz, ptot2 c smeared variables the prefix v denotes the varience of the variable real smrehad, smremu, thetamu, smrenu real smrx, smry, smrq2, u1, u2, u3, u4, u5, u6, u7, u8 real z1, z2, z3, z4 real emcal, hadcal, hadcal2, mucal, vthnu, vthmum, vthmu real ptotmu, smrptotmu, smrhadcal, smremcal, smrmucal, smrthmu real shadcal, semcal, smucal, sthnu, sthmum c ******************************************************* c --- MUSPIN input/output common c q = +1. for mu+ decay and -1. for mu- decay real q, pelec(4), pnue(4), pnumu(4) common/mugen/ q, pelec, pnue, pnumu c --- oscillation info stored in /osc/ real l ! baseline length (cm) real l0 ! characteristic length (cm) real nk ! neutrino energy (eV) real psum(200,100) ! probability sums real mvals(100) ! delta m^2 points real svals(200) ! sin^2 2 theta points integer nm ! number of mval points integer ns ! number of sval points common/osc/ l, l0, nk, nm, ns, mvals, svals, psum c --- local variables real straight ! straight section length (cm) integer ndecays logical output, leptoflag real polarize ! -1 -> spin opposite momentum vector ! 0 -> unpolarized real qmu double precision npbeam,dpbeam,pbeam real beamp real sigma_x real sigma_y real sigma_dx real sigma_dy real m_d real n_mu real r_max ! for r < r_max (km) will use nu to calculate real length real dmsq real ssq2t real n_a double precision movere, thetae real thetaex,thetaey,phie real emass double precision mumass double precision twopi double precision beta, gamma, bgam, ebeam real p4(4) integer idecay real cost,pmax real pt_nue, pt_num real nx(28) integer nvar2 parameter (nvar2=55) real space(1000) real ny(nvar2) real space2(1000) real nu_wt, polar real det_wt real dummy1, dummy2, dummy3, dummy4, dummy5, dummy6 real dummy7, dummy8 real rndm, test, pprob real cos_theta, sin_theta, phi real a_max, n_int_nue, n_int_num, n_nue, n_num cc real n_int_numo, numo_wt cc real n_int_nue_noosc, n_int_nue_nomat cc real n_int_numo_nomat cc real n_int_num_noosc, n_int_num_nomat, nueo_wt,n_int_nueo cc real n_int_nueo_nomat, enueo_detav real nue_flux, num_flux, enue_av, enum_av real emu_detav, ee_detav, enue_detav, enum_detav cc real enumo_detav real y_nu, y_antinu real ptelec, ptelec_new, pzelec_new real ptnue, ptnue_new, pznue_new real ptnumu, ptnumu_new, pznumu_new real eav_nue, eav_num, eav_e real rvec(2), xprime, yprime real sum(4),sumcm(4),masscm,masslab real rg32,dummyg integer ntotin, nto2in real pe_lab_prime(4), pnue_lab_prime(4), pnum_lab_prime(4) integer i, j real r1, r2, wt cc real p_osc, P_osc_nomat, dm2_0, s2t_0 cc real sigma_nue_noosc, sigma_nue_nomat, sigma_numo_nomat cc real sigma_num_noosc, sigma_num_nomat, sigma_nueo_nomat integer nutype real*8 ll real z, enu, delta, ssq_0, p_noosc real dummy13 real eleftmu,erightmu real elefte,erighte #include "param.inc" parameter (eleftmu = ( -.5 + sin2tw)**2) parameter (erightmu = (sin2tw)**2) parameter (elefte = (.5+sin2tw)**2) parameter (erighte = sin2tw**2) c --- set up polarization polar = polarize c --- zero interaction counters n_nue = 0.0 n_num = 0.0 eav_nue = 0.0 eav_num = 0.0 eav_e = 0.0 enue_av = 0.0 enum_av = 0.0 emu_detav = 0.0 ee_detav = 0.0 enue_detav= 0.0 enum_detav= 0.0 cc enumo_detav= 0.0 cc enueo_detav= 0.0 n_int_nue = 0.0 n_int_num = 0.0 cc n_int_numo = 0.0 cc n_int_nueo = 0.0 cc n_int_nue_noosc = 0.0 cc n_int_num_noosc = 0.0 cc n_int_nue_nomat = 0.0 cc n_int_num_nomat = 0.0 cc n_int_numo_nomat = 0.0 cc n_int_nueo_nomat = 0.0 q = qmu ! set up mu charge (+ or -) c --- prepare ntuple c********************** change nutple if lepto is being used ************ if (.not.leptoflag) then call book_ntuple elseif (leptoflag) then call book_ntuple2 endif c ******************************************************************** c --- compute normalization area associated with r_max cc a_max = 3.142 * (r_max*10**5)**2 ! cm**2 a_max = 1.00 c --- compute total number of neutrinos that each generated neutrino c represents for a 1 year (10**7 secs) expt nu_wt = n_mu / float(ndecays) det_wt = (n_a*m_d) c --- load lab beam four-vector c --- get characteristic length in matter L = length * 10.**5 ! baseline length (cm) c call get_L0(L,L0) ! L0 (cm) c c --- Write header c 9010 format(' nu_e -> nu_mu') 9011 format(' nu_e bar -> nu_mu bar') 9012 format(' nu_mu -> nu_e') 9013 format(' nu_mu bar -> nu_e bar') 9014 format(' nu_mu bar -> nu_e bar and nu_e -> nu_mu') 9015 format(' nu_mu -> nu_e and nu_e bar -> nu_mu bar') write(6,9999) qmu, n_mu, polarize, npbeam,dpbeam, * sigma_dx,sigma_dy, * sigma_x, sigma_y, straight, * m_d, nu_wt, * det_wt, sin2tw, ndecays 9999 format( * ' muon charge = ',f6.3,/, * ' no. muons per year = ',e20.10,/, c * ' f_dk = ',g15.7,/, * ' polarization = ',f10.5,/, * ' muon energy = ',g20.10,' GeV',/, * ' muon energy spread = ',g20.10,' %',/, * ' beam divergence sigdx= ',g15.7,/, * ' beam divergence sigdy= ',g15.7,/, * ' beam size(cm) sigx = ',g15.7,/, * ' beam size(cm) sigy = ',g15.7,/, * ' length of straight section in cm = ',g15.7,/, * ' detector thickness = ',g20.10,' g/cm2',/, * ' neutrino wt = ',g20.10,/, * ' detector wt = ',g20.10,/, * ' sin^2 2 theta W = ',g20.10,/, * ' No. decays generated = ',i15,/, * ' ----------------------------------------------------',/) if(output) write(6,9020) 9020 format(' px',4x,'py',4x,'pz',4x,'E',4x,'posc') c ************lepto init ************* if (leptoflag) then c setup c cut(1) = 0.0001 cut(5) = 1.00 cut(7) = 3.00 c cut(11) = 1.0 lst(1)=1 c output switch lst(3)=5 lst(7)=1 lst(8)=12 lst(9)=5 c lst(10)=0 lst(17)=1 lst(19)=-1 c lst(22) selects target 1=proton or 2=neutron lst(22)=nucleon c interaction=2 is charged current, interaction=3 is neutral current lst(23)=interaction c right or left handed neutrinos if (lepin.gt.0) then lst(30)=-1 parl(6)=-1.0 elseif (lepin.lt.0) then lst(30)=1 parl(6)=1.0 endif lst(31)=1 c parl(1)=A of target, parl(2)=Z of target parl(1)=targa parl(2)=targz c parl(11)=0.01 c parl(15)=0.01 beamp = npbeam call linit(0,lepin,beamp,0,interaction) print*, parl(23)," cross section in pb calc. by LEPTO init" endif c************************************** c --- loop over ndecays muon decays do idecay = 1, ndecays * print *, ' event wts ', nu_wt, det_wt pbeam = npbeam*(1. +rg32(dummy3)*dpbeam) p4(1) = 0. p4(2) = 0. p4(3) = pbeam p4(4) = dsqrt(pbeam**2 + mumass**2) ebeam = dsqrt(pbeam**2 + mumass**2) gamma = ebeam / mumass beta = pbeam / ebeam bgam = beta*gamma r_max = 3./gamma if(polarize.eq.999)then if(rndm(dummyg).gt.0.5)then polar = 1 else polar = -1 endif endif event = float(idecay) call muspin(polar) ! generates 4-vectors in mu rest-frame ! for polarization polar pe_cms(1) = pelec(1) pe_cms(2) = pelec(2) pe_cms(3) = pelec(3) pe_cms(4) = pelec(4) pnue_cms(1) = pnue(1) pnue_cms(2) = pnue(2) pnue_cms(3) = pnue(3) pnue_cms(4) = pnue(4) pmax = mumass/2. cost = pnue_cms(3)/pnue_cms(4) emax_nue = gamma*pmax*(1+beta*cost) pnum_cms(1) = pnumu(1) pnum_cms(2) = pnumu(2) pnum_cms(3) = pnumu(3) pnum_cms(4) = pnumu(4) cost = pnum_cms(3)/pnum_cms(4) emax_num = gamma*pmax*(1+beta*cost) c --- set up boost vector so that we can get from cms to lab c --- E' = gamma*E+beta*gamma*pz c --- pz' = beta*gamma*E + gamma*pz, c where primed quantities are in lab, unprimed in mu rest frame c and beta and gamma are for beam in lab. pe_lab(1) = pe_cms(1) pe_lab(2) = pe_cms(2) pe_lab(3) = bgam*pe_cms(4) + gamma*pe_cms(3) pe_lab(4) = gamma*pe_cms(4) + bgam*pe_cms(3) pnue_lab(1) = pnue_cms(1) pnue_lab(2) = pnue_cms(2) pnue_lab(3) = bgam*pnue_cms(4) + gamma*pnue_cms(3) pnue_lab(4) = gamma*pnue_cms(4) + bgam*pnue_cms(3) pnum_lab(1) = pnum_cms(1) pnum_lab(2) = pnum_cms(2) pnum_lab(3) = bgam*pnum_cms(4) + gamma*pnum_cms(3) pnum_lab(4) = gamma*pnum_cms(4) + bgam*pnum_cms(3) c c --- If beam divergence switch on, smear beam direction c xprime = 0. yprime = 0. if(sigma_dx.gt.0..or.sigma_dy.gt.0.) then rvec(1) = rg32(dummyg) rvec(2) = rg32(dummyg) if(sigma_dx.gt.0.) xprime = sigma_dx * rvec(1) if(sigma_dy.gt.0.) yprime = sigma_dy * rvec(2) pnum_lab(1) = pnum_lab(1)+pnum_lab(4)*xprime pnum_lab(2) = pnum_lab(2)+pnum_lab(4)*yprime pnum_lab(3) = sqrt(pnum_lab(4)**2 + - pnum_lab(1)**2 - pnum_lab(2)**2) pnue_lab(1) = pnue_lab(1)+pnue_lab(4)*xprime pnue_lab(2) = pnue_lab(2)+pnue_lab(4)*yprime pnue_lab(3) = sqrt(pnue_lab(4)**2 + - pnue_lab(1)**2 - pnue_lab(2)**2) c call rotate(pnum_lab,xprime,yprime) c call rotate(pnue_lab,xprime,yprime) end if c c --- generate energy fractions given to lepton and antilepton c in charged current interaction 40 y_nu = rndm(dummy3) pprob = 0.833 + 0.167*(1.-y_nu)**2 ! Assumes 20% antiquark contribution test = rndm(dummy4) if(test.gt.pprob) goto 40 y_nu = 1. - y_nu ! Note: y_nu = E_mu/E_nu = 1-y 50 y_antinu = rndm(dummy5) pprob = 0.167 + 0.833*(1.-y_antinu)**2 ! Assumes 20% antiquark contribution test = rndm(dummy6) if(test.gt.pprob) goto 50 y_antinu = 1. - y_antinu ! Note: y_antinu = E_mu/E_nu = 1-y c generate electron scattering eynum = rndm(dummy3) eynue = eynum C C --- Consider electron scattering first C movere = emass/pnum_lab(4) thetae = acos( * (1.+movere)*sqrt(eynum*eynum-movere*movere)/ * (eynum+movere)) phie = twopi*rndm(dummy6) thetaex = thetae*cos(phie) thetaey = thetae*sin(phie) factm = nueconst*rescale*pnum_lab(4) facte = nueconst*rescale*pnue_lab(4) snum = 2.*emass*pnum_lab(4) snue = 2.*emass*pnue_lab(4) if(qmu.gt.0)then ncenuml = factm*(eleftmu*(1-eynum)**2) ncenuel = facte*(erighte*(1-eynum)**2) ncenumr = factm*(erightmu) ncenuer = facte*(elefte) ccenuel = 0.0 ccenuml = 0.0 else ncenuml = factm*(eleftmu) ncenuel = facte*(erighte) ncenumr = factm*(erightmu*(1-eynum)**2) ncenuer = facte*(elefte*(1-eynue)**2) cc ncenuml = factm*(eleftmu) cc ncenuel = facte*(erighte*(1-eynue)**2) cc ncenumr = factm*(erightmu*(1-eynum)**2) cc ncenuer = facte*(elefte) ccenumr = 0.0 ccenuel = 0.0 if(snue.gt.mumass**2)then ccenuer = facte*(eynue*(1.-eynue)- + mumass**2*(snue-mumass**2)/snue) else ccenuer = 0.0 endif if(snum.gt.mumass**2)then ccenuml = factm*(snum-mumass**2)/snum else ccenuml = 0.0 endif endif n_nue = n_nue + 1. enue_av = enue_av + pnue_lab(4) c --- compute nue -> numu oscillation probability c --- compute sigma for event and energy of leptons from CC interactions if(qmu.gt.0.) then ! mu+ decay -> nue sigma_nue = 0.67e-38 * pnue_lab(4) ee_det = y_nu * pnue_lab(4) else ! mu- decay -> anti nue sigma_nue = 0.34e-38 * pnue_lab(4) ee_det = y_antinu * pnue_lab(4) end if n_int_nue = n_int_nue + nu_wt*sigma_nue*det_wt ee_detav = ee_detav + ee_det*nu_wt*sigma_nue*det_wt enue_detav = enue_detav + * pnue_lab(4)*nu_wt*sigma_nue*det_wt 100 continue C C --- Now consider nu_mu C pt_num = sqrt(pnum_lab(1)**2 + pnum_lab(2)**2) r_num = length * pt_num / pnum_lab(4) if(r_num.gt.r_max) goto 200 n_num = n_num + 1. enum_av = enum_av + pnum_lab(4) c --- compute sigma for event and energy of leptons from CC interactions if(qmu.gt.0.) then ! mu+ decay -> anti numu sigma_num = 0.34e-38 * pnum_lab(4) emu_det = y_antinu * pnum_lab(4) else ! mu- decay -> anti nue + numu sigma_num = 0.67e-38 * pnum_lab(4) emu_det = y_nu * pnum_lab(4) end if n_int_num = n_int_num + nu_wt*sigma_num*det_wt emu_detav = emu_detav + emu_det*nu_wt*sigma_num*det_wt enum_detav = enum_detav + * pnum_lab(4)*nu_wt*sigma_num*det_wt 200 continue c ****************** lepto ******************* if (leptoflag) then c open(unit=15, file='other.log') if (abs(lepin).eq.14) then ptot = sqrt(pnum_lab(3)**2+pnum_lab(1)**2+pnum_lab(2)**2) elseif (abs(lepin).eq.12) then ptot = sqrt(pnue_lab(3)**2+pnue_lab(1)**2+pnue_lab(2)**2) endif c nuthetax = pnum_lab(1)/ptot c nuthetay = pnum_lab(2)/ptot c cut on total neutrino momentum, if the cut is too low lepto c will freeze up if (ptot.gt.8.5) then c write(15,101) ptot if (abs(lepin).eq.14) then p(1,1)=pnum_lab(1) p(1,2)=pnum_lab(2) p(1,3)=pnum_lab(3) p(1,4)=ptot p(1,5)=0. p(2,1)=0. p(2,2)=0. p(2,3)=0. p(2,4)=0.938 p(2,5)=0.938 elseif (abs(lepin).eq.12) then p(1,1)=pnue_lab(1) p(1,2)=pnue_lab(2) p(1,3)=pnue_lab(3) p(1,4)=ptot p(1,5)=0. p(2,1)=0. p(2,2)=0. p(2,3)=0. p(2,4)=0.938 p(2,5)=0.938 endif c generate the event... call lepto if (lst(21).ne.0) then goto 45 endif c call lulist(1) call luedit(2) c do i= 1, n c if (k(i,1).eq.1) then c write(15,103) k(i,2),p(i,1),p(i,2),p(i,3),p(i,5) c endif c enddo c elseif (ptot.lt.5.0) then c write(15,102) ptot endif c 101 format('used event total momentum ',f7.3) c 102 format('unused event total momentum ',f7.3) c 103 format('pid ',i4,' xmom ',f7.4,' ymom ',f7.4,' zmom ' c & ,f7.4,' mass ',f7.4) c smear x,y,q2 for events... emcal = 0 hadcal2 = 0 do i=1, n ptot2=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2) if (p(i,5).gt.0.110.and.ptot2.gt.0.5) then hadcal2 = hadcal2 + p(i,4) elseif (p(i,5).lt.0.0001.and.ptot2.gt.0.5) then emcal = emcal + p(i,4) elseif (0.0001.lt.p(i,5).and.p(i,5).lt.0.110) then mucal = p(i,4) ptotmu = sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2) thetamu = atan(sqrt(p(i,1)**2+p(i,2)**2)/p(i,3)) endif enddo u1 = rndm(dummy1) u2 = rndm(dummy2) u3 = rndm(dummy3) u4 = rndm(dummy4) u5 = rndm(dummy5) u6 = rndm(dummy6) u7 = rndm(dummy7) u8 = rndm(dummy8) z1 = sin(2*3.1415926*u1)*sqrt(-2*log(u2)) z2 = sin(2*3.1415926*u3)*sqrt(-2*log(u4)) z3 = sin(2*3.1415926*u5)*sqrt(-2*log(u6)) z4 = sin(2*3.1415926*u7)*sqrt(-2*log(u8)) smrhadcal = hadcal2 + shadcal*sqrt(hadcal2)*z1 smremcal = emcal + semcal*sqrt(emcal)*z2 smrmucal = mucal + smucal*mucal*z3 if (abs(lepin).eq.14) then smrptotmu = sqrt(smrmucal**2 - mumass**2) elseif (abs(lepin).eq.12) then smrptotmu = sqrt(smrmucal**2 - 0.000511**2) endif vthnu = (sthnu*0.001)**2 vthmum = (sthmum*0.001)**2 vthmu = vthnu + vthmum smrthmu = thetamu + sqrt(vthmu)*z4 smrehad = smrhadcal + smremcal smrenu = smrehad + smrmucal smrq2 = mumass**2 + 2*smrenu*(smrmucal - & smrptotmu*cos(smrthmu)) smrx = smrq2/(2*0.938*smrehad) smry = smrehad/smrenu endif c leptoflag c ******************************************** if (.not.leptoflag) then ityp = 2 throwwt = ccgen(ityp,-qmu,pnum_lab(4),cckin) c muon ny(1) = float(idecay) ! event number ny(2) = pbeam ! muon beam energy ny(3) = qmu ! muon charge ny(4) = polar ! muon polarization ny(5) = ityp ! 1 for e, 2 for mu neutrino ny(6) = -qmu ! -1 for anti, 1 for neutrino c cm muon neutrino ny(7) = pnum_cms(4) ! neutrino energy cm ny(8) = pnum_cms(3)/pnum_cms(4) ! cos theta in cm c decay effects ny(9) = -straight*rndm(dummy13) ! decay position ny(10) = xprime ! beam angle in x ny(11) = yprime ! beam angle in y c lab muon neutrino ny(12) = pnum_lab(1) ! lab momentum ny(13) = pnum_lab(2) ny(14) = pnum_lab(3) ny(15) = pnum_lab(4) ny(16) = emax_num c ny(17) = nu_wt*det_wt/rescale !nue lum / gr/cm2 hydrogen ny(18) = eynum ! ybj for nu-e scatter ny(19) = ncenuml ny(20) = ncenumr ny(21) = ccenuml ny(22) = ccenumr ny(23) = rg32(dummy3)*sigma_x ny(24) = rg32(dummy3)*sigma_y ny(25) = thetaex ny(26) = thetaey ! will be theta for nu-e c ny(27) = det_wt*nu_wt/rescale ! lum for dis ny(28) = sigma_num*rescale ! lum*sigma for dis isoscalar do kk = 1,25 ny(kk+28) = cckin(kk) enddo if(r_num.lt.r_max) then cc call hfn(1,nx) ! Fill ntuple call hfn(2,ny) end if elseif (leptoflag.and.ptot.gt.8.5) then ityp = 2 throwwt = ccgen(ityp,-qmu,pnum_lab(4),cckin) do i=1,n if (k(i,1).eq.1) then c muon ny(1) = float(idecay) ! event number ny(2) = pbeam ! muon beam energy ny(3) = qmu ! muon charge ny(4) = polar ! muon polarization c decay effects ny(9) = -straight*rndm(dummy13) ! decay position ny(10) = xprime ! beam angle in x ny(11) = yprime ! beam angle in y if (abs(lepin).eq.14) then c lab muon neutrino ny(5) = 2 ! 1 for e, 2 for mu neutrino ny(6) = -qmu ! -1 for anti, 1 for neutrino c cm muon neutrino ny(7) = pnum_cms(4) ! neutrino energy cm ny(8) = pnum_cms(3)/pnum_cms(4) ! cos theta in cm ny(12) = pnum_lab(1) ! lab momentum ny(13) = pnum_lab(2) ny(14) = pnum_lab(3) ny(15) = pnum_lab(4) elseif (abs(lepin).eq.12) then ny(5) = 1 ! 1 for e, 2 for mu neutrino ny(6) = -qmu ! -1 for anti, 1 for neutrino c cm muon neutrino ny(7) = pnue_cms(4) ! neutrino energy cm ny(8) = pnue_cms(3)/pnue_cms(4) ! cos theta in cm ny(12) = pnue_lab(1) ! lab momentum ny(13) = pnue_lab(2) ny(14) = pnue_lab(3) ny(15) = pnue_lab(4) endif ny(16) = emax_num c ny(17) = nu_wt*det_wt/rescale !nue lum / gr/cm2 hydrogen ny(18) = eynum ! ybj for nu-e scatter ny(19) = ncenuml ny(20) = ncenumr ny(21) = ccenuml ny(22) = ccenumr ny(23) = rg32(dummy3)*sigma_x ny(24) = rg32(dummy3)*sigma_y ny(25) = thetaex ny(26) = thetaey ! will be theta for nu-e c ny(27) = det_wt*nu_wt/rescale ! lum for dis ny(28) = sigma_num*rescale ! lum*sigma for dis isoscalar ny(29) = i ny(30) = p(i,1) ny(31) = p(i,2) ny(32) = p(i,3) ny(33) = p(i,4) ny(34) = p(i,5) ny(35) = x ny(36) = y ny(37) = w2 ny(38) = q2 ny(39) = u ny(40) = smrq2 ny(41) = smrx ny(42) = smry if(r_num.lt.r_max) then cc call hfn(1,nx) ! Fill ntuple call hfn(2,ny) end if endif enddo endif go to 45 ityp = 1 throwwt = ccgen(ityp,qmu,pnue_lab(4),cckin) c electron ny(1) = float(idecay) ! event number ny(2) = pbeam ! muon beam energy ny(3) = qmu ! muon charge ny(4) = polar ! muon polarization ny(5) = ityp ! 1 for e, 2 for mu neutrino ny(6) = qmu ! -1 for anti, 1 for neutrino c cm muon neutrino ny(7) = pnue_cms(4) ! neutrino energy cm ny(8) = pnue_cms(3)/pnue_cms(4) ! cos theta in cm c decay effects ny(9) = -straight*rndm(dummy13) ! decay position ny(10) = xprime ! beam angle in x ny(11) = yprime ! beam angle in y c lab muon neutrino ny(12) = pnue_lab(1) ! lab momentum ny(13) = pnue_lab(2) ny(14) = pnue_lab(3) ny(15) = pnue_lab(4) ny(16) = emax_nue ny(17) = nu_wt*det_wt/rescale ! lum for nu-e on proton/gr/cm2 ny(18) = eynue ! ybj for nu-e scatter ny(19) = ncenuel ny(20) = ncenuer ny(21) = ccenuel ny(22) = ccenuer ny(23) = rg32(dummy3)*sigma_x ny(24) = rg32(dummy3)*sigma_y ny(25) = thetaex ! will be theta for nu-e ny(26) = thetaey ny(27) = det_wt*nu_wt/rescale ny(28) = sigma_nue*rescale ! sigma for dis isoscalar do kk = 1,25 ny(kk+28) = cckin(kk) enddo if(r_nue.lt.r_max) then cc call hfn(1,nx) ! Fill ntuple call hfn(2,ny) end if 45 continue enddo c --- output print*, parl(24)," cross section in pb estimate from MC events" nue_flux = n_nue * nu_wt / a_max ! m^-2 per year num_flux = n_num * nu_wt / a_max if(n_num.gt.0.) then enum_av = enum_av / n_num emu_detav = emu_detav / n_int_num enum_detav = enum_detav / n_int_num cc enueo_detav = enueo_detav / n_int_nueo end if if(n_nue.gt.0.) then enue_av = enue_av / n_nue ee_detav = ee_detav / n_int_nue enue_detav = enue_detav / n_int_nue cc enumo_detav = enumo_detav / n_int_numo end if write(6,1000) n_num, num_flux, n_int_num, cc * n_int_num_noosc, n_int_num_nomat, cc * n_int_nueo, n_int_nueo_nomat, * enum_av, emu_detav, enum_detav 1000 format( * ' ----------------------------------------------------',/, * ' No.generated nu_mu inside r_max = ',g15.7,/, * ' nu_mu flux inside r max (m-2) = ',g15.7,/, * ' No. nu_mu interactions per year = ',g15.7,/, cc * ' No. nu_mu interactns/yr no osc = ',g15.7,/, cc * ' No. nu_mu interactns/yr no mat = ',g15.7,/, cc * ' No. osc. nu_e interactns per year = ',g15.7,/, cc * ' No. osc. nu_e interactns/yr no mat= ',g15.7,/, * ' Av. numu energy inside r_max = ',g15.7,' GeV',/, * ' Av. detected mu energy = ',g15.7,' GeV',/, * ' Av. energy of interacting nu_mu = ',g15.7,' GeV',/) write(6,2000) n_nue, nue_flux, n_int_nue, cc * n_int_nue_noosc, n_int_nue_nomat, cc * n_int_numo, n_int_numo_nomat, * enue_av, ee_detav, enue_detav 2000 format( * ' No.generated nu_e inside r_max = ',g15.7,/, * ' nu_e flux inside r max (m-2) = ',g15.7,/, * ' No. nu_e interactions per year = ',g15.7,/, cc * ' No. nu_e interactns/yr no osc = ',g15.7,/, cc * ' No. nu_e interactns/yr no mat = ',g15.7,/, cc * ' No. osc. nu_mu interactns per year = ',g15.7,/, cc * ' No. osc. nu_mu interactns/yr no mat= ',g15.7,/, * ' Av. nue energy inside r_max = ',g15.7,' GeV',/, * ' Av. detected e energy = ',g15.7,' GeV',/, * ' Av. energy of interacting nu_e = ',g15.7,' GeV',/) CALL HPAGSZ(50) call hcdir('//NUE',' ') CALL HROUT(2,icycle,' ') print *,' icycle was ', icycle CALL HREND('NUE') stop end subroutine book_ntuple c --- common /ev/ stores output event record real event real pe_lab(4), pe_cms(4) real pnue_lab(4), pnue_cms(4), r_nue, sigma_nue real pnum_lab(4), pnum_cms(4), r_num, sigma_num real emu_det, ee_det,emax_nue, emax_num real sigma_numo common /ev/ event, zdecay, pe_lab, pe_cms, * pnue_lab, pnue_cms, r_nue, sigma_nue, * pnum_lab, pnum_cms, r_num, sigma_num, * emu_det, ee_det, emax_nue, emax_num integer status, lun character*2 flstat integer nvar,nvar2 parameter( nvar = 28 ) parameter(nvar2 = 55) character*80 filnam character*8 names( nvar ) character*8 names2(nvar2) character*4 topdir integer npags, ihist, hmem, icycle parameter( npags = 6000) parameter( ihist = 128 * npags ) common /PAWC/ hmem( ihist ) data topdir/ 'NUE' / data filnam/ 'mu.hbook' / DATA NAMES/'event','pxe', 'pye', 'pze', 'ee', 'eecm', + 'pxnue', 'pynue', 'pznue', 'enue', 'enucm', 'rnue', 'wtnue', + 'pxnum', 'pynum', 'pznum', 'enum', 'enumcm', 'rnum', 'wtnum', + 'emudet', 'eedet','zdecay', 'pbeam','coscmnue','emaxnue', + 'coscmnum', + 'emaxnum'/ data names2/'event','pbeam','qmu','pol','ityp','nuanti', + 'ecm','cosc', + 'zdecay', 'xdiv','ydiv', + 'px','py','pz','enu','emaxnu', + 'nuelum', + 'ybje', + 'ncel','ncer','ccel','ccer','beamx','beamy', + 'thetaex','thetaey','dislum','ccisonu', + 'x','y','q2', + 'u','d','s','c','ubar','dbar','sbar','cbar','gluon', + 'sigmap','sigmaiso','weight','xsi','charm','charmiso','ewt' + ,'thetalx','thetaly','plep', + 'dummy3','dummy4','dummy5','dummy6','dummy7'/ INTEGER NSEL, NBYTES CHARACTER*72 FILEOUT COMMON/EVSEL/FILEOUT,NSEL,NBYTES INTEGER IQUEST COMMON/QUEST/IQUEST(100) c --- hbook call hlimit(ihist) C C --- open PAW file C lun = 12 flstat = 'NQ' iquest(10) = 65000 c call hropen( lun, topdir, filnam, flstat, 8192, status ) c flstat = 'N' call hropen( lun, topdir, filnam, flstat, 1024, status ) C C --- book histograms C cc call hbookn(1, 'two body', nvar, topdir, 6000, names ) call hbookn(2, 'data', nvar2, topdir, 6000, names2 ) call hbook1(500,'el energy distribution: r2GeV', *100,0.,50.,0.) call hbook1(503,'nu_tau CC E distribution: r2GeV', *100,0.,50.,0.) return end subroutine book_ntuple2 c creates ntuple that includes lepto events c --- common /ev/ stores output event record real event real pe_lab(4), pe_cms(4) real pnue_lab(4), pnue_cms(4), r_nue, sigma_nue real pnum_lab(4), pnum_cms(4), r_num, sigma_num real emu_det, ee_det,emax_nue, emax_num real sigma_numo common /ev/ event, zdecay, pe_lab, pe_cms, * pnue_lab, pnue_cms, r_nue, sigma_nue, * pnum_lab, pnum_cms, r_num, sigma_num, * emu_det, ee_det, emax_nue, emax_num integer status, lun character*2 flstat integer nvar,nvar2 parameter( nvar = 28) parameter(nvar2 = 42) character*80 filnam character*8 names( nvar ) character*8 names2(nvar2) character*4 topdir integer npags, ihist, hmem, icycle parameter( npags = 6000) parameter( ihist = 128 * npags ) common /PAWC/ hmem( ihist ) data topdir/ 'NUE' / data filnam/ 'mu.hbook' / DATA NAMES/'event','pxe', 'pye', 'pze', 'ee', 'eecm', + 'pxnue', 'pynue', 'pznue', 'enue', 'enucm', 'rnue', 'wtnue', + 'pxnum', 'pynum', 'pznum', 'enum', 'enumcm', 'rnum', 'wtnum', + 'emudet', 'eedet','zdecay', 'pbeam','coscmnue','emaxnue', + 'coscmnum', + 'emaxnum'/ data names2/'event','pbeam','qmu','pol','ityp','nuanti', + 'ecm','cosc', + 'zdecay', 'xdiv','ydiv', + 'px','py','pz','enu','emaxnu', + 'nuelum', + 'ybje', + 'ncel','ncer','ccel','ccer','beamx','beamy', + 'thetaex','thetaey','dislum','ccisonu','i','pxi','pyi', + 'pzi','ei','mi','x','y','w2', + 'q2','nu','smrq2','smrx','smry'/ INTEGER NSEL, NBYTES CHARACTER*72 FILEOUT COMMON/EVSEL/FILEOUT,NSEL,NBYTES INTEGER IQUEST COMMON/QUEST/IQUEST(100) c --- hbook call hlimit(ihist) C C --- open PAW file C lun = 12 flstat = 'NQ' iquest(10) = 65000 c call hropen( lun, topdir, filnam, flstat, 8192, status ) c flstat = 'N' call hropen( lun, topdir, filnam, flstat, 1024, status ) C C --- book histograms C cc call hbookn(1, 'two body', nvar, topdir, 6000, names ) call hbookn(2, 'data', nvar2, topdir, 6000, names2 ) call hbook1(500,'el energy distribution: r2GeV', *100,0.,50.,0.) call hbook1(503,'nu_tau CC E distribution: r2GeV', *100,0.,50.,0.) return end real FUNCTION fnue(polar,nx,cth) C polarized muon decay C see PRD 39,3532 (1989) G.Barr,T.K.Gaisser,T.Stanev C or "Cosmic Rays and Particle Physics" T.K. Gaisser chaper 7 c --- input polarization added ... S.G. implicit none real nx, cth, polar C Functions are dN/(dx dphi d(costh)) C For mu+ relative to spin direction (same as parent direction in mu cm) C FNUMU(X,CTH) = ( 2.*X*X*(3.-2.*X) - 2.*X*X*(1.-2.*X)*CTH )/(4.*PI) C FNUE (X,CTH) = (12.*X*X*(1.-X) - 12.*X*X*(1.-X) *CTH )/(4.*PI) C For mu-, CTH term flips sign relative to spin, and and spin flips C sign relative to parent direction, so can use same function C by using pion in muon rest frame. C now strip them to run faster FNUE = nX*nX*(1.-nX)*(1.-polar*CTH) RETURN END real FUNCTION fnumu(polar,nx,cth) C polarized muon decay C see PRD 39,3532 (1989) G.Barr,T.K.Gaisser,T.Stanev C or "Cosmic Rays and Particle Physics" T.K. Gaisser chaper 7 c --- input polarization added .... S.G. implicit none real nx, cth, polar C Functions are dN/(dx dphi d(costh)) C For mu+ relative to spin direction (same as parent direction in mu cm) C FNUMU(X,CTH) = ( 2.*X*X*(3.-2.*X) - 2.*X*X*(1.-2.*X)*CTH )/(4.*PI) C FNUE (X,CTH) = (12.*X*X*(1.-X) - 12.*X*X*(1.-X) *CTH )/(4.*PI) C For mu-, CTH term flips sign relative to spin, and and spin flips C sign relative to parent direction, so can use same function C by using pion in muon rest frame. C now strip them to run faster FNUMU = nX*nX* ( 3. - polar*CTH - 2.*nX*(1.-polar*CTH) ) RETURN END SUBROUTINE MUSPIN(polar) C polarized muon decay C see PRD 39,3532 (1989) G.Barr,T.K.Gaisser,T.Stanev C or "Cosmic Rays and Particle Physics" T.K. Gaisser chaper 7, p.89 implicit none c --- MUSPIN input/output common c q = +1. for mu+ decay and -1. for mu- decay real q, pelec(4), pnue(4), pnumu(4) common/mugen/ q, pelec, pnue, pnumu real rndm, polar real FNUMU,FNUE ! statement func real emass, mumass, twopi ! parameters real FNUMU_max, FNUE_max, Elep, nX, STH, CTH real CTH_numu, X_numu, phi_numu real CTH_nue, X_nue, phi_nue real CTH_elec, X_elec, phi_elec real STH_elec, STH_nue, STH_numu real e_elec, e_numu, e_nue real dummy1, dummy2, dummy3, dummy4, dummy5, dummy6 real dummy7, dummy8, dummy9, dummy10, dummy11, dummy12 parameter (EMASS=0.0005109990615) parameter (mumass = .105658389) parameter (TWOPI=6.28318530717958648) parameter (FNUMU_max=2.) ! MAX VALUE OF FNUMU IS AT X=1 ,CTH=-1 parameter (FNUE_max= 8./27.) ! MAX VALUE OF FNUE IS AT X=2/3,CTH= 1 c--------------------------------------------------------------------------- C Functions are dN/(dx dphi d(costh)) C For mu+ relative to spin direction (same as parent direction in mu cm) C FNUMU(X,CTH) = ( 2.*X*X*(3.-2.*X) - 2.*X*X*(1.-2.*X)*CTH )/(4.*PI) C FNUE (X,CTH) = (12.*X*X*(1.-X) - 12.*X*X*(1.-X) *CTH )/(4.*PI) C For mu-, CTH term flips sign relative to spin, and and spin flips C sign relative to parent direction, so can use same function C by using pion in muon rest frame.... c we wont assume anything about muon polarization ... so we generate c wrt spin direction in muon restframe ... and need to fold in mu sign C now strip them to run faster c FNUMU(X,CTH) = X*X* ( 3. - CTH - 2.*X*(1.-CTH) ) c FNUE (X,CTH) = X*X*(1.-X)*(1.-CTH) c--------------------------------------------------------------------------- C anti- numu from mu+ decay of numu from mu- decay 20 CTH_numu = 2.*RNDM(dummy1) - 1. X_numu=RNDM(dummy2) IF (FNUMU_max*RNDM(dummy3) * .GT.FNUMU(polar,X_numu,CTH_numu)) GOTO 20 cth_numu = q * cth_numu ! added to be correct for mu- decay e_numu = 0.5 * X_numu * mumass C nu e from mu+ decay or anti-nue from mu- decay 30 CTH_nue = 2.*RNDM(dummy4) - 1. X_nue=RNDM(dummy5) IF (FNUE_max*RNDM(dummy6).GT.FNUE(polar,X_nue,CTH_nue)) GOTO 30 cth_nue = q * cth_nue ! added to be correct for mu- decay e_nue = 0.5 * X_nue * mumass C In approx that electron mass = 0, e+ has same dist as anti-numu 40 CTH_elec = 2.*RNDM(dummy7) - 1. X_elec=RNDM(dummy8) IF (FNUMU_max*RNDM(dummy9) * .GT.FNUMU(polar,X_elec,CTH_elec)) GOTO 40 cth_elec = q * cth_elec ! added to be correct for mu- decay e_elec = 0.5 * X_elec * mumass phi_numu = twopi * rndm(dummy10) phi_nue = twopi * rndm(dummy11) phi_elec = twopi * rndm(dummy12) STH_numu = sqrt(1.-CTH_numu**2) pnumu(1) = e_numu * STH_numu * cos(phi_numu) pnumu(2) = e_numu * STH_numu * sin(phi_numu) pnumu(3) = e_numu * CTH_numu pnumu(4) = e_numu STH_nue = sqrt(1.-CTH_nue**2) pnue(1) = e_nue * STH_nue * cos(phi_nue) pnue(2) = e_nue * STH_nue * sin(phi_nue) pnue(3) = e_nue * CTH_nue pnue(4) = e_nue STH_elec = sqrt(1.-CTH_elec**2) pelec(1) = e_elec * STH_elec * cos(phi_elec) pelec(2) = e_elec * STH_elec * sin(phi_elec) pelec(3) = e_elec * CTH_elec pelec(4) = e_elec RETURN END real function ccgen(ityp,anti,Enu,yy) * HMS 11-2 implicit none real anti integer ityp real rndm real dummy real Enu real yy(25) double precision wt,xwt,wtx double precision pmass double precision vus,vud double precision ny,nq2 double precision space(1000) double precision nx double precision qcall double precision upv,dnv,usea,dsea,str,chm,bot,top,glu double precision space2(1000) double precision xsi double precision upvc,dnvc,useac,dseac,strc,chmc,botc,topc,gluc double precision space3(1000) double precision uq,d,s,c,ubar,dbar,sbar,cbar,gluon double precision un,dn, unbar,dnbar, sn, snbar cc double precision q2max, q2min,lq2max,lq2min double precision strup,strun,struiso double precision yfac double precision ecm2 double precision mc double precision charm double precision charmiso double precision udfac,usfac double precision sin2th,theta,thetax,thetay,plep,phi c double precision xsi double precision q2cteq real cwt double precision gfermi double precision twopi double precision sigmafac double precision gevcm integer first data q2cteq/4./ parameter (pmass = .938) data mc/1.4/ parameter (twopi = 6.28318530717958648) parameter (gfermi = 1.16639 E-5) ! in GeV-2 parameter (vus = .22) parameter (gevcm = .197E-13) parameter (sigmafac = pmass*gfermi*gfermi/twopi*gevcm*gevcm*4.) double precision rescale parameter (rescale = 10E24) ! factor to make cross sections reasonable data first/0/ ccgen = 0.0 if(first.eq.0)then vud = sqrt(1.-vus*vus) udfac = vud*vud usfac = vus*vus c print *,'factors',udfac,usfac first = 1 endif first = first + 1 ecm2 = 2.*Enu*pmass * throw events in nx with a (1-nx)**2 distribution nx = wtx(dummy) ! xbj ny = rndm(dummy) nq2 = ecm2*nx*ny qcall =sqrt( max(nq2,q2cteq)) call structm(nx,qcall,Upv,Dnv,Usea,Dsea,Str,Chm,Bot,Top,Glu) * try to generate charm * check kinematics before bothering with parton distributions xsi=nx+dble(mc*mc/pmass/2/enu/ny ) c c-- now apply the thresholds on xsi. c if(xsi.lt.mc*mc/2/pmass/enu)go to 100 if(xsi.gt.1) go to 100 c c-- and y c if(ny.lt.mc*mc/2/pmass/enu/xsi)go to 100 c if(sqrt(ecm2)*ny.lt.3)go to 100 call structm(xsi,qcall,Upvc,Dnvc,Useac *,Dseac,Strc,Chmc,Botc,Topc,Gluc) * print * ,' strange ',xsi,strc,s,strc/s go to 200 100 continue xsi = -1 upvc = 0 dnvc = 0 useac = 0 dseac = 0 strc = 0 chmc = 0 botc = 0 topc = 0 gluc = 0 c print 'xsi' , xsi,qcall,nx,enu 200 continue * handle electrons and rotation from mass to weak usfac = vus*vus udfac = 1.-usfac * proton uq = upv+usea ubar = usea d = (dnv+dsea)*udfac + str*usfac dbar = dsea*udfac + str*usfac s = strc*udfac + (dnvc+dseac)*usfac sbar = strc*udfac + dseac*usfac c = chm cbar = chm gluon = glu if (first.lt.1000)then c print *, 'quarks ',nx,nq2, upv,usea, uq, dnv,dsea,str, c *d,udfac,usfac endif * neutron dn = upv+usea dnbar = usea un = (dnv+dsea)*udfac + str*usfac unbar = dsea*udfac + str*usfac sn = strc*udfac + (upvc+useac)*usfac snbar = strc*udfac + useac*usfac c print *,'strange2',d,dbar,s,sbar,udfac,usfac yfac = (1-ny)*(1-ny) * make cross section factors if(anti.gt.0)then charm = s charmiso = .5*(s+sn) strup = (d + s) + (ubar+cbar)*yfac strun = (dn + sn) +(unbar+cbar)*yfac struiso = .5*(strup+strun) else charm = sbar charmiso = .5*(sbar+snbar) strup = (dbar + sbar) + ( uq +c)*yfac strun = (dnbar + snbar) + (un +c) *yfac struiso = .5*(strup+strun) endif * make a cross section weight * compensate for (1-nx)**2 nx distribution if(anti.gt.0) then wt = 1/xwt(nx) else wt = 1/xwt(nx) endif if(nx.gt.0.99) wt = 0.0 cc cwt = enu*8.2 cwt = sigmafac*enu*rescale ! factor that goes in front of SF part yy(1) = nx ! xbj yy(2) = ny ! ybj yy(3) = nq2 ! yy(4) = uq ! quarks and gluons in proton yy(5) = d yy(6) = s yy(7) = c yy(8) = ubar yy(9) = dbar yy(10) = sbar yy(11) = cbar yy(12) = gluon yy(13) = strup *cwt ! dsigma dx dy on proton yy(14) = struiso * cwt ! dsgima dx dy on iso yy(15) = wt ! event wt yy(16) = xsi ! slow rescaling x yy(17) = charm *cwt ! charm function * y factors in proton yy(18) = charmiso *cwt ! charm function * y factors on iso plep = enu*(1-ny) sin2th = sqrt(nq2/4./enu/plep) theta = 2.*asin(sin2th) phi = twopi*rndm(dummy) thetax = theta*cos(phi) thetay = theta*sin(phi) yy(20) = thetax yy(21) = thetay yy(22) = plep ccgen = wt ! weight of event, needs to g * print *, ' event wt ', nx, wt, cwt return end double precision function xwt(nx) implicit none real rndm double precision nx,wtx real dum xwt = 3.*(1.-nx)*(1.-nx) return entry wtx(dum) wtx = 1.- rndm(dum)**(.333) return end #include "structm_cteq.f"