REAL FUNCTION eeppb(hid0,jeloss,ikin) ********************************************************* * * * This file was generated by HUWFUN. * * * ********************************************************* * * Ntuple Id: 20 * Ntuple Title: eep coin * Creation: 28/02/2005 16.02.58 * ********************************************************* * LOGICAL CHAIN CHARACTER*128 CFILE INTEGER IDNEVT,NCHEVT,ICHEVT REAL OBS(13) * COMMON /PAWIDN/ IDNEVT,OBS COMMON /PAWCHN/ CHAIN, NCHEVT, ICHEVT COMMON /PAWCHC/ CFILE * *-- Ntuple Variable Declarations * real xvert,yvert,zvert,eex,eey,eez,ee,ppx,ppy, + ppz,pproton,xe7,ye7 + ,ze7,xp8,yp8,zp8,ppxv,ppyv,ppzv,ppv,eexv,eeyv,eezv,eev,xrast0 + ,yrast0 * COMMON /PAWCR4/ xvert,yvert,zvert,eex,eey,eez,ee,ppx,ppy,ppz + ,pproton,xe7,ye7,ze7,xp8,yp8,zp8,ppxv,ppyv,ppzv,ppv,eexv,eeyv + ,eezv,eev,xrast0,yrast0 * * *-- Enter user code here c change the function name to refer to the variable being plotted integer hid0, hid1, hidp, jloss, iep, i12,jeloss,ikin real ztgtex,ztgtey,ztgtpx,ztgtpy,ztgte,ztgtp,ztgtep real e3s,p3s, xgs,xdummy,p1s,p2s double precision dex,dey,dpx,dpy,de,dp,dpe,p2bite double precision ea,eb,e1,e2,e3,ma,mb,m1,m2,m3,pa,p1,p2,p3 double precision pa0,pax,pay,paz,p1x,p1y,p1z,p2x,p2y,p2z,p3x, + p3y,p3z,xg,p1m,p2m,q,qz,csthq,csth2,pfact,p2set,p2lo,p2hi double precision aa,ab,a1,a2,a3,da,db,d1,d2,d3,amu, dpa2,p32,eloss double precision zlo1,zhi1,zlo2,zhi2,zlo3,zhi3,delEa,delE1,delE2 double precision the1,thp1,the2,thp2,the3,thp3,thtgt,rad, + snthe,snthp,theta2(7) double precision tLH1,tLH2,tLH3,tLHx,cap,wall,pb,kap,ni,pbx double precision ededxH,ededxAl,ededxPb,ededxC,ededxNi double precision pdedxH,pdedxAl,pdedxPb,pdedxC,pdedxNi,delta c c this file is setup for 208Pb(e,e'p)X, diamond two foils, TGT2, kin01, pmiss=200MeV/c data amu/931.49432/, aa/0./,ab/208./,a1/0./,a2/1./,a3/207./ data da/0.511/,db/-21.764/,d1/0.511/,d2/7.289/,d3/-21.044/ data pa0/2445.0/,rad/57.2958/,p2set/990.0/,p2bite/0.045/ data the1/19.80/,thp1/60.67/,the2/23.45/, + the3/23.19/,thp3/58.27/,thtgt/0./ data zlo1/9.54/,zhi1/8.54/,zlo2/-3./,zhi2/3.0/, + zlo3/-9.0/,zhi3/-6.0/ data tlH1/0.0/,tLH2/0./,tLH3/0.0/,thHx/0./, + cap/0./ data wall/0.0/,kap/0.02528/,ni/0.18614/,pb/0.227/ data ededxH/6.3/,ededxAl/2.4/,ededxPb/1.9/,ededxC/2.5/, + ededxNi/2.2/ data pdedxH/6.3/,pdedxAl/2.4/,pdedxPb/1.6/,pdedxC/2.8/, + pdedxNi/2.2/ data theta2/53.19,58.92,47.46,64.72,41.67,70.53,35.85/ c *********************************** c first call specf.f to convert the geant variables to spectrometer variables xdummy = specf(hid0) c**** this section calculates ztgt in various ways **** c1 hid1=hid0 dex = (xe7-xrast0)*ee/eex ztgtex=ze7-eez/ee*dex call hfill(hid1,ztgtex,1.,1.) c2 hid1=hid1+1 dey= (ye7-yrast0)*ee/eey ztgtey=ze7-eez/ee*dey call hfill(hid1,ztgtey,1.,1.) c3 hid1=hid1+1 dpx= (xp8-xrast0)*pproton/ppx ztgtpx=zp8-ppz/pproton*dpx call hfill(hid1,ztgtpx,1.,1.) c4 hid1=hid1+1 dpy= (yp8-yrast0)*pproton/ppy ztgtpy=zp8-ppz/pproton*dpy call hfill(hid1,ztgtpy,1.,1.) c5 hid1=hid1+1 de= (dex+dey)/2. ztgte=ze7-eez/ee*de call hfill(hid1,ztgte,1.,1.) c6 hid1=hid1+1 dp = (dpx+dpy)/2. ztgtp=zp8-ppz/pproton*dp call hfill(hid1,ztgtp,1.,1.) c7 hid1=hid1+1 ztgtep = (ztgte+ztgtp)/2. call hfill(hid1,ztgtep,1.,1.) iep = hid0 + 1005 ! this is the ID for the 2dim emiss vs pmiss spectrum i12 = hid0 + 1006 ! ID for 2dim emom vs pmom p2lo = (1.- p2bite)*p2set ! low momentum acceptance of proton spectrometer p2hi = (1.+ p2bite)*p2set !high momentum acceptance of proton spectrometer thp2 = theta2(ikin) c hid1 = hid1 + 1 c hidp = hid1 c************************************************** c c ******* start calculating the missing mass spectrum c convert GeV to MeV -- make sure pa0 in the data statement is in MeV/c p1x=1000.*eex p1y=1000.*eey p1z=1000.*eez p2x=1000.*ppx p2y=1000.*ppy p2z=1000.*ppz eloss = jeloss c determine masses ma = (aa*amu + da) mb = (ab*amu + db) m1 = (a1*amu + d1) m2 = (a2*amu + d2) m3 = (a3*amu + d3) eb = mb ! stationary target pbx = pb/cos(thtgt/rad) ! target is tilted wrt beam so path length is greater than pb c c determine which target is the vertex site c pa0 is the incident beam energy, pa is the incident energy at the vertex if((ztgtex.ge.zlo1).and.(ztgtex.le.zhi1))then delta = eloss*(cap*ededxAl+tlH1*ededxH+pbx*ededxPb/2.) ! assume scatter from mid point c write(6,*)'pa0, delta ',pa0,delta pa = pa0 - delta ! incident energy at vertex ea = sqrt(pa*pa + ma*ma) c calculate p1 electron momenta at vertex due to ionization energy loss p1 = sqrt(p1x*p1x+p1y*p1y+p1z*p1z) ! this is the measured momentum at sp07 c xg = xgauss() ! get dp/p to simulate spectrometer resolution c xgs = xg ! single precision value c call hfill(hidp,xgs,1.,1.) c p1 = p1*(1. + xg) ! smear p1 to simulate spectrometer resolution snthe = sin(the1/rad) delta = eloss*(tLHx/snthe*ededxH+ni/snthe*ededxNi+kap*ededxC+ + pb/2./cos((the1+thtgt)/rad)*ededxPb) ! energy loss of electron out of target p1 = p1 + delta ! p1 vertex momentum e1 = sqrt(p1*p1 + m1*m1) c calculate proton momentum, p2, at vertex correcting for ionization loss p2 = sqrt(p2x*p2x+p2y*p2y+p2z*p2z) ! this is the measured momentum at sp08 c p2 = p2*(1. + xgauss()) ! smear p2 to simulate spectrometer resolution snthp = sin(thp1/rad) delE2 = (tlHx/snthp*pdedxH + ni/snthp*pdedxNi + kap*pdedxC + + pb/2./cos((thp1-thtgt)/rad)*pdedxPb) ! this is the energy loss of the proton e2 = sqrt(p2*p2+m2*m2) ! proton energy at sp08 dp2 = e2*delE2/p2*eloss ! momentum loss of proton p2 = p2 + dp2 ! this is the proton momentum at the vertex e2 = sqrt(p2*p2 + m2*m2) ! proton energy at vertex c calculate missing momentum and energy at the vertex p32 = pa*pa-2.*pa*p1*cos(the1/rad)-2.*pa*p2*cos(thp1/rad)+p1*p1 + +p2*p2 +2.*p1*p2*cos((the1+thp1)/rad) e3 = sqrt((ea+eb-e1-e2)*(ea+eb-e1-e2) - p32) - m3 !excitation energy in GeV e3s = e3 ! single precision value hid1 = hid1+1 c write(6,*)'ea eb e1 e2 e3 p32',ea,eb,e1,e2,e3,p32 call hfill(hid1,e3s,1.,1.) c this section for making a cut on ee c hid1 = hid1+1 c if(ee.gt.1.9)then c call hfill(hid1,e3,1.,1.) c endif elseif((ztgtex.ge.zlo2).and.(ztgtex.le.zhi2))then pa = pa0 -eloss*(cap*ededxAl + tlH2*ededxH + 0.5*pbx*ededxPb) ! assume scatter from mid point ea = sqrt(pa*pa + ma*ma) c calculate p1 electron momenta at vertex due to ionization energy loss p1m = sqrt(p1x*p1x+p1y*p1y+p1z*p1z) ! this is the measured mometum at sp07 c xg = xgauss() ! get dp/p to simulate spectrometer resolution c xgs = xg c call hfill(hidp,xgs,1.,1.) c p1 = p1m*(1. + xg) ! smear p1 to simulate spectrometer resolution snthe = sin(the2/rad) p1 = p1m+eloss*(tLHx/snthe*ededxH+ni/snthe*ededxNi+kap*ededxC+ + pb/2./cos((the2+thtgt)/rad)*ededxPb) ! p1 vertex momentum c write(6,*)'p1,p1m',p1,p1m e1 = sqrt(p1*p1 + m1*m1) c calculate proton momentum, p2, at vertex correcting for ionization loss p2m = sqrt(p2x*p2x+p2y*p2y+p2z*p2z) ! this is the measured momentum at sp08 c p2 = p2m*(1. + xgauss()) ! smear p2 to simulate spectrometer resolution snthp = sin(thp2/rad) delE2 = (tlHx/snthp*pdedxH + ni/snthp*pdedxNi + kap*pdedxC + + pb/2./cos((thp2-thtgt)/rad)*pdedxPb) ! this is the energy loss of the proton e2 = sqrt(p2m*p2m+m2*m2) ! proton energy at sp08 dp2 = e2*delE2/p2m*eloss ! momentum loss of proton p2 = p2m + dp2 ! this is the proton momentum at the vertex c write(6,*)'p2,p2m',p2,p2m e2 = sqrt(p2*p2 + m2*m2) ! proton energy at vertex c calculate missing momentum and energy at the vertex c recalculate the p1x,p1y,p1z,p2x,p2y,p2z components assuming the energy loss did not change directions p1x=p1x/p1m*p1 p1y=p1y/p1m*p1 p1z=p1z/p1m*p1 p2x=p2x/p2m*p2 p2y=p2y/p2m*p2 p2z=p2z/p2m*p2 p32=pa*pa+p1*p1+p2*p2+2.*(p1x*p2x+p1y*p2y+p1z*p2z)-2.*pa*(p1z+p2z) p3 = sqrt(p32) c p32 = pa*pa -2.*pa*p1*cos(the2/rad)-2.*pa*p2*cos(thp2/rad)+p1*p1 c + +p2*p2 +2.*p1*p2*cos((the2+thp2)/rad) e3 = sqrt((ea+eb-e1-e2)*(ea+eb-e1-e2) - p32) - m3 ! missing energy e3s = e3 p3s = p3 p1s = p1 p2s = p2 hid1 = hid1+2 c write(6,*)'ea eb e1 e2 e3 p32',ea,eb,e1,e2,e3,p32 if((p2.ge.p2lo).and.(p2.le.p2hi))then call hfill(hid1,e3s,1.,1.) ! only fill emiss spectrum if p2 is in proton acceptance q = sqrt(pa*pa+p1*p1-2.*pa*p1z) ! three momentum transfer q qz = pa - p1z ! z componemt of q csthq = qz/q csth2 = p2z/p2 pfact = 1. if(csth2.gt.csthq)then ! if csth2>csthq the proton is emitted forward of q pfact = -1. ! and we will define the missing momentum p3s as negative endif p3s = pfact*p3s call hfill(iep,e3s,p3s,1.) !emiss vs pmiss endif call hfill(i12,p1s,p2s,1.) !emom vs pmom elseif((ztgtex.ge.zlo3).and.(ztgtex.le.zhi3))then pa = pa0 -eloss*(cap*ededxAl + tlH3*ededxH + 2.5*pbx*ededxPb) ! assume scatter from mid point ea = sqrt(pa*pa + ma*ma) c calculate p1 electron momenta at vertex due to ionization energy loss p1 = sqrt(p1x*p1x+p1y*p1y+p1z*p1z) ! this is the measured mometum at sp07 c xg = xgauss() ! get dp/p to simulate spectrometer resolution c xgs = xg c call hfill(hidp,xgs,1.,1.) c p1 = p1*(1. + xg) ! smear p1 to simulate spectrometer resolution snthe = sin(the3/rad) p1 = p1+eloss*(tLHx/snthe*ededxH+ni/snthe*ededxNi+kap*ededxC+ + pb/2./cos((the3+thtgt)/rad)*ededxPb) ! p1 vertex momentum e1 = sqrt(p1*p1 + m1*m1) c calculate proton momentum, p2, at vertex correcting for ionization loss p2 = sqrt(p2x*p2x+p2y*p2y+p2z*p2z) ! this is the measured momentum at sp08 c p2 = p2*(1. + xgauss()) ! smear p2 to simulate spectrometer resolution snthp = sin(thp3/rad) delE2 = (tlHx/snthp*pdedxH + ni/snthp*pdedxNi + kap*pdedxC + + pb/2./cos((thp3-thtgt)/rad)*pdedxPb) ! this is the energy loss of the proton e2 = sqrt(p2*p2+m2*m2) ! proton energy at sp08 dp2 = e2*delE2/p2*eloss ! momentum loss of proton p2 = p2 + dp2 ! this is the proton momentum at the vertex e2 = sqrt(p2*p2 + m2*m2) ! proton energy at vertex c calculate missing momentum and energy at the vertex p32 = pa*pa -2.*pa*p1*cos(the3/rad)-2.*pa*p2*cos(thp3/rad)+p1*p1 + +p2*p2 +2.*p1*p2*cos((the3+thp3)/rad) e3 = sqrt((ea+eb-e1-e2)*(ea+eb-e1-e2) - p32) - m3 ! missing energy e3s = e3 hid1 = hid1+3 c write(6,*)'ea eb e1 e2 e3 p32',ea,eb,e1,e2,e3,p32 call hfill(hid1,e3s,1.,1.) endif eepmm = 1. * * END c ************************************************************************* REAL FUNCTION specf(hid0) ********************************************************* * * * This file was generated by HUWFUN. * * * ********************************************************* * * Ntuple Id: 20 * Ntuple Title: eep coin * Creation: 11/03/2005 14.42.29 * ********************************************************* * LOGICAL CHAIN CHARACTER*128 CFILE INTEGER IDNEVT,NCHEVT,ICHEVT REAL OBS(13) * COMMON /PAWIDN/ IDNEVT,OBS COMMON /PAWCHN/ CHAIN, NCHEVT, ICHEVT COMMON /PAWCHC/ CFILE * *-- Ntuple Variable Declarations * REAL xvert,yvert,zvert,eex,eey,eez,ee,ppx,ppy,ppz,pproton,xe7,ye7 + ,ze7,xp8,yp8,zp8,ppxv,ppyv,ppzv,ppv,eexv,eeyv,eezv,eev,xrast0 + ,yrast0 * COMMON /PAWCR4/ xvert,yvert,zvert,eex,eey,eez,ee,ppx,ppy,ppz + ,pproton,xe7,ye7,ze7,xp8,yp8,zp8,ppxv,ppyv,ppzv,ppv,eexv,eeyv + ,eezv,eev,xrast0,yrast0 * * *-- Enter user code here * simulate spectrometer smearing given geant input real eeg(5),ees(5),ppg(5),pps(5),x integer hid0,ip,ih,iv,ie,ipr common /nuts/ eeg,ees,ppg,pps,ip,ih,iv,ie,ipr ip = hid0 + 1000 ih = hid0 + 1001 iv = hid0 + 1002 ie = hid0 + 1003 ipr = hid0 + 1004 eeg(1)=eex ! electron geant values at sp07 eeg(2)=eey eeg(3)=eez eeg(4)=xe7 eeg(5)=ye7 ppg(1)=ppx ! proton geant values at sp08 ppg(2)=ppy ppg(3)=ppz ppg(4)=xp8 ppg(5)=yp8 x=speclr(ispec) c change the sp07 and sp08 variables to spectrometer values eex=ees(1) eey=ees(2) eez=ees(3) xe7=ees(4) ye7=ees(5) ppx=pps(1) ppy=pps(2) ppz=pps(3) xp8=pps(4) yp8=pps(5) specf = 1. * END ********************************************** c broaden the geant variables to simulate the spectrometer effects real function speclr(ispec) real eeg(5),ees(5),ppg(5),pps(5),xh,xv,xp,xee,xpp real dh,dv,dp,thh,thv,th0h,th0v,tnh,tnv,x integer ispec,i,ip,ih,iv,ie,ipr data dh/0.3e-3/,dv/1.0e-3/,dp/1.0e-4/ common /nuts/ eeg,ees,ppg,pps,ip,ih,iv,ie,ipr xee = 0. xpp = 0. do 10 i=1,3 xee=xee+eeg(i)*eeg(i) ! left spectrometer = electron xpp=xpp+ppg(i)*ppg(i) ! right spectrometer = proton 10 continue xee = sqrt(xee) ! geant electron momentum xpp = sqrt(xpp) ! geant proton momentum c first calculate spectrometer variables for electron xp = dp*xgauss() ! momentum uncertainty for electron xee = xee*(1. + xp) !spectrometer electron momentum call hfill(ip,xp,1.,1.) call hfill(ie,xee,1.,1.) tnh = eeg(1)/eeg(3) ! tan(phi_tgt) from geant variables tnv = eeg(2)/eeg(3) ! tan(th_tgt) from geant variables th0h = atan(tnh) th0v = atan(tnv) xh=dh*xgauss() !smear phi_tgt xv=dv*xgauss() !smear th_tgt call hfill(ih,xh,1.,1.) !only fill these histo once per event call hfill(iv,xv,1.,1.) thh = th0h + xh ! phi_tgt from spec variables thv = th0v + xv ! th_tgt from spec variables tnh = tan(thh) tnv = tan(thv) ees(3) = xee/sqrt(1.+tnh*tnh+tnv*tnv) !spectrometer electron z momentum ees(2) = ees(3)*tnv !spectrometer electron y momentum ees(1) = ees(3)*tnh !spectrometer electron x momentum ees(4) = eeg(4) !spectrometer electron xe7 position ees(5) = eeg(5) !spectrometer electron ye7 position c c second calculate spectrometer proton variables xpp = xpp*(1. + dp*xgauss()) !spectrometer proton momentum call hfill(ipr,xpp,1.,1.) tnh = ppg(1)/ppg(3) ! tan(phi_tgt) from geant variables tnv = ppg(2)/ppg(3) ! tan(th_tgt) from geant variables th0h = atan(tnh) th0v = atan(tnv) xh=dh*xgauss() !smear phi_tgt xv=dv*xgauss() !smear th_tgt thh = th0h + xh ! phi_tgt from spec variables thv = th0v + xv ! th_tgt from spec variables tnh = tan(thh) tnv = tan(thv) pps(3) = xpp/sqrt(1.+tnh*tnh+tnv*tnv) !spectrometer proton z momentum pps(2) = pps(3)*tnv !spectrometer proton y momentum pps(1) = pps(3)*tnh !spectrometer proton x momentum pps(4) = ppg(4) !spectrometer proton xp8 position pps(5) = ppg(5) !spectrometer proton yp8 position c speclr = 1. end c ************************************************************************* c this function simulates a gaussian distribution by a sum of 12 random numbers from 0. to 1. real function xgauss() double precision sum integer i sum = 0. do 20 i=1,12 sum = sum + rndm() 20 continue xgauss = (sum-6.0) end