C-- SANE_HMS_physics from HMS single arm MC C-- ebeam = (GeV) Energy of incident beam C-- HMS_mom = hsdelta(GeV) Central momentum of HMS C-- HMS_angle = hstheta(radian) Angle between HMS and outgoing beam real function SANE_HMS_pfmc(nrun) c real function SANE_HMS_pfmc implicit none c include 'sane_hms.cmn' include ? logical first data first /.true./ save first real ebeam,beamE,bigcmom,bigcang real ebigc,thbigc real xbSANE,Q2SANE,W2SANE,WSANE real nuSANE,ESANE,thetaSANE,ySANE real xbHMS,Q2HMS,W2HMS,WHMS real nuHMS,EHMS,dEHMS,thetaHMS,phiHMS,yHMS real targxp,targyp,targx,targy,beamz,SRy,SRx real targxpfp,targypfp,targxfp,targyfp real cdt,avei,time,charge,fac,wt real mpr,d2r,r2d,hel integer i,nrun integer n_W,W_index real Wmin,Wmax,Emin,Emax,hsp0,W_int vector trkeff vector trkefffit vector runinfo vector Winfo vector cx vector cwdamc mpr = 0.938272 d2r = 3.14159/180. r2d = 1 / d2r ccc ebeam = 4.734 cc ebeam = 4.731 c ebeam = 5.892 cc ebeam = 5.895 n_W=int(Winfo(1)) Wmin=Winfo(2) Wmax=Winfo(3) Emin=runinfo(2)-0.4 Emax=runinfo(2)+0.4 W_int=0.015 * The MC gives rate assuming 100nA beam current * Need to normalize MC relative to 100nA and computer deadtime. c fac=time*(avei/100.*(1 - cdt)) ebeam=runinfo(1)/1000. hsp0=runinfo(2) if(runinfo(6).eq.1) then cdt=runinfo(7) charge=runinfo(9)*1000. else if(runinfo(6).eq.0) then cdt=runinfo(8) charge=runinfo(10)*1000. endif ccc fac=avei/100.*(1 - cdt)*time fac=charge/100.*(1 - cdt) wt=numer_wt*fac*Winfo(4)!*0.88085!*1.01895!*0.978429!*1.00051 c wt=1. !!! Default value assuming 100nA CC print*,cdt,avei,time,fac if(first) then first=.false. call hbook1(10000,'Q^2! ' ,50,0.5,4.,0) call hbook1(20000,'x?bj! ' ,50,0.,1.,0) call hbook1(300000,'W ' ,100,0.7,3.0,0) cccc 72984 72991 cc call hbook1(30000,'W ' ,50,2.15,2.25,0) call hbook1(30000,'W ' ,21,2.1475,2.2525,0) call hbook1(30001,'W ' ,21,2.1475,2.2525,0) call hbook1(31000,'W ' ,21,2.1475,2.2525,0) call hbook1(31001,'W ' ,21,2.1475,2.2525,0) cccc 72290 call hbook1(30000,'W ' ,27,1.8475,1.9825,0) call hbook1(30001,'W ' ,27,1.8475,1.9825,0) call hbook1(31000,'W ' ,27,1.8475,1.9825,0) call hbook1(31001,'W ' ,27,1.8475,1.9825,0) cccc 73014 73019 cc call hbook1(30000,'W ' ,50,1.2,1.6,0) call hbook1(30000,'W ' ,81,1.1975,1.6025,0) call hbook1(30001,'W ' ,81,1.1975,1.6025,0) call hbook1(31000,'W ' ,81,1.1975,1.6025,0) call hbook1(31001,'W ' ,81,1.1975,1.6025,0) c call hbook1(30000,'W ' ,50,0.7,1.8,0) call hbook1(30000,'W ' ,100,0.9,2.9,0) cc call hbook1(30000,'W ' ,50,0.9,1.8,0)!CHe-C call hbook1(30000,'W ',n_W,Wmin,Wmax,0) call hbook1(30001,'W ',n_W,Wmin,Wmax,0) call hbook1(31000,'W ',n_W,Wmin,Wmax,0) call hbook1(31001,'W ',n_W,Wmin,Wmax,0) cc call hbook1(30000,'W ' ,35,1.080,1.610,0)!4.7_p3200th20_df cc call hbook1(30001,'W ' ,35,1.080,1.610,0) cc call hbook1(31000,'W ' ,35,1.080,1.610,0) cc call hbook1(31001,'W ' ,35,1.080,1.610,0) cc call hbook1(30000,'W ',20,2.15,2.25,0)!4.7_p2200th16 cc call hbook1(30001,'W ',20,2.15,2.25,0) cc call hbook1(31000,'W ',20,2.15,2.25,0) cc call hbook1(31001,'W ',20,2.15,2.25,0) call hbook1(40000,'E?HMS! ' ,50,Emin,Emax,0) cc call hbook1(30000,'W ' ,13,2.1,2.25,0)!CHe-C cc call hbook1(30000,'W ' ,50,1.9,2.5,0)!CHe-C cccc 72984 72991 c call hbook1(40000,'E?HMS! ' ,50,1.7,2.5,0) cccc 72290 c call hbook1(40000,'E?HMS! ' ,50,2.2,3.,0) cccc 72505 cc call hbook1(40000,'E?HMS! ' ,50,2.7,3.5,0) cccc 73014 73019 cc call hbook1(40000,'E?HMS! ' ,50,2.5,4.,0) cccc 72782 c call hbook1(40000,'E?HMS! ' ,50,3.8,5.,0) cccc 72984 72991 72782 call hbook1(50000,'[Q]?HMS! ' ,50,12.,25.,0) cccc 73014 73019 c call hbook1(50000,'[Q]?HMS! ' ,50,15.,25.,0) call hbook1(41000,'[D]E?HMS! ' ,50,-18.,18.,0) call hbook1(41100,'[D]E?HMS! ' ,50,-18.,18.,0) call hbook1(60000,'nu?HMS! ' ,50,0.,4.,0) call hbook1(70000,'y ' ,50,0.,1.,0) call hbook2(120000,'Q^2! vs. x?bj! HMS' ,50,0.,2., > 50,0.,7.0,0) call hbook2(130000,'Q^2! vs. W HMS' ,50,0.5,2., > 50,0.,4.0,0) call hbook2(480000,'E vs. [Q] HMS' ,50,15.,25., > 50,2.,4.0,0) call hbook1(4100,'SR?x!' ,50,-1.5,1.5,0) call hbook1(4200,'SR?y!' ,50,-1.5,1.5,0) call hbook1(1000,'xptarg' ,50,-0.23,0.13,0) call hbook1(1001,'xptarg no [D]P cut' ,50,-0.23,0.13,0) call hbook1(1100,'xptargfp' ,50,-0.08,0.08,0) call hbook1(1200,'xtargfp' ,50,-50.,50.,0) call hbook1(2000,'yptarg' ,50,-0.06,0.06,0) call hbook1(2001,'yptarg no [D]P cut' ,50,-0.06,0.06,0) call hbook1(2100,'yptargfp' ,50,-0.04,0.03,0) call hbook1(2200,'ytargfp' ,50,-30.,25.,0) call hbook1(3000,'beam z' ,50,-10.,10.,0) call hbook1(3001,'beam z no [D]P cut' ,50,-10.,10.,0) call hbook2(12000,'xptar vs. yptar' ,50,-0.06,0.06, > 50,-0.23,0.13,0) call hbook2(141000,'xptar vs. [D]E?HMS!' ,50,-9.,9., > 50,-0.23,0.13,0) call hbook2(111100,'xptar vs. xtar' ,50,-1.5,1.5, > 50,-0.13,0.13,0) call hbook2(211100,'ytar vs. xtar' ,50,-1.5,1.5, > 50,-4.,4.,0) call hbook2(311000,'beam z vs. xtar' ,50,-1.5,1.5, > 50,-10.,10.,0) call hbook2(321000,'beam z vs. ytar' ,50,-4.,4., > 50,-10.,10.,0) call hbook1(1110,'xptarg' ,50,-0.23,0.13,0) call hbook1(2110,'yptarg' ,50,-0.06,0.06,0) call hbook2(412000,'SR?x! vs. ytar' ,50,-2.5,2.5, > 50,-1.5,1.5,0) call hbook2(412100,'SR?x! vs. yptar' ,50,-0.06,0.06, > 50,-1.5,1.5,0) call hbook2(421100,'SR?y! vs. xptar' ,50,-0.23,0.13, > 50,-1.5,1.5,0) call hbook1(1110,'xtar',50,-4.,4.,0) call hbook1(1111,'ytar',50,-4.,4.,0) call hbook1(1112,'[F]?HMS! ' ,100,120,200.,0) call hbook2(1113,'[D]E?HMS vs W',n_W,Wmin,Wmax,50,-18.,18.,0) endif C do i=1,nclust C C hel = i_helicity C ebigc = eclust(i)*1000 C nuSANE = ebeam - ebigc C thetaSANE = r2d*atan(xclust(i)/345.) + 40. C Q2SANE = 4.* ebeam*ebigc * sin(d2r*thetaSANE/2)**2 C xbSANE = Q2SANE / 2. / mpr / nuSANE C W2SANE = mpr**2 + 2.*mpr*nuSANE - Q2SANE C WSANE = sqrt(max(0.,W2SANE)) C ySANE = 1 - ebigc/ebeam C C enddo dEHMS = hsdelta!+0.7 c print *,'hse=',hse ccc EHMS = hse EHMS = hse!hsp0*(1.+dEHMS/100.) c print *,'EHMS=',EHMS c EHMS = (hse/1000) c EHMS = sqrt(dEHMS**2 + mpr**2) nuHMS = ebeam - EHMS * thetaHMS = hstheta*r2d thetaHMS = hstheta phiHMS = hsphi c Q2HMS = 4.* ebeam*dEHMS * sin(d2r*thetaHMS/2)**2 c xbHMS = Q2HMS / 2. / mpr / nuHMS c W2HMS = mpr**2 + 2.*mpr*nuHMS - Q2HMS c WHMS = sqrt(max(0.,W2HMS)) Q2HMS = Q2 xbHMS = Q2HMS / 2. / mpr / nuHMS WHMS = W !-0.006!- 0.01 !to match in elastic region also (run 73014) yHMS = (1 - EHMS/ebeam) c Q2HMS = 4.* ebeam*EHMS * sin((thetaHMS*d2r)/2)**2 c xbHMS = Q2HMS / 2. / mpr / nuHMS c W2HMS = mpr**2 + 2.*mpr*nuHMS - Q2HMS c WHMS = sqrt(max(0.,W2HMS)) targxp = hsxptar targyp = hsyptar targx = hsxtar targy = hsytar beamz = hszbeam SRy = ybeami SRx = xbeami targxpfp = hsxpfp targypfp = hsypfp targxfp = hsxfp targyfp = hsyfp c print*,Q2 ccc If you want to turn on tracking efficiency on MC, you need to turn off data weight if(abs(hsxfp).le.100) then wt=wt*(trkefffit(4)*hsxfp**3+ & trkefffit(3)*hsxfp**2+trkefffit(2)*hsxfp+trkefffit(1)) endif CCC if(dEHMS.gt.-0.4.and.dEHMS.lt.1.1)then CCC wt=wt*0.9 CCC endif CCC if(A.eq.4)then CCC wt=wt*0.86 CCC endif c if(WHMS.ge.(cx(1)-W_int/2.).and.WHMS.le.(cx(n_W)+W_int/2.))then c W_index=0 c do i=1,n_W c if(abs(WHMS-cx(i)).le.(W_int/2.))then c W_index=i c endif c enddo c if(abs(1.-cwdamc(W_index)).lt.0.99)then c wt=wt*cwdamc(W_index) c endif c else c wt=0. c endif if(WHMS.ge.Wmin.and.WHMS.le.Wmax)then ccc if(dEHMS.lt.-1.and.dEHMS.gt.-30)then c if(abs(hsxfp).le.60)then c if(thetaHMS.gt.16.and.thetaHMS.lt.20.and.abs(dEHMS).lt.10) then c if(abs(dEHMS).lt.30.and.abs(hsxfp).lt.37.5)then !.and.abs(dEHMS).gt.1.!.and.abs(hsxfp).lt.30 ccc if(abs(hsxfp).lt.50)then ccc if(dEHMS.gt.-8.and.dEHMS.lt.12)then if(dEHMS.gt.-8.and.dEHMS.lt.8)then ccc if(EHMS.lt.runinfo(5)*1.08.and.EHMS.gt.runinfo(5)*0.92)then ccc & .and.SRy.gt.0.)then call hf1(3001,beamz,wt) call hf1(41000,dEHMS,wt) call hf1(41100,dEHMS,1.) call hf2(141000,dEHMS,targxp,wt) c print *,'EHMS=',EHMS, 'hse=',hse c if(targxp.gt.-0.2.and.targxp.lt.-0.1)then call hf1(1000,targxp,wt) call hf1(2000,targyp,wt) call hf1(1100,targxpfp,wt) call hf1(2100,targypfp,wt) call hf1(1200,hsxfp,wt) call hf1(2200,hsyfp,wt) call hf1(3000,beamz,wt) call hf1(4100,SRx,wt) call hf1(4200,SRy,wt) call hf2(412000,targy,SRx,wt) call hf2(412100,targyp,SRx,wt) call hf2(421100,targxp,SRy,wt) call hf1(10000,Q2HMS,wt) call hf1(20000,xbHMS,wt) call hf1(300000,WHMS,wt) call hf1(30000,WHMS,wt) call hf1(30001,WHMS,1.) if(A.eq.1)then call hf1(31000,WHMS,wt) call hf1(31001,WHMS,1.) endif call hf1(40000,EHMS,wt) call hf1(50000,thetaHMS,wt) call hf1(60000,nuHMS,wt) call hf1(70000,yHMS,wt) call hf2(12000,targyp,targxp,wt) call hf2(111100,targx,targxp,wt) call hf2(211100,targx,targy,wt) call hf2(311000,targx,beamz,wt) call hf2(321000,targy,beamz,wt) call hf2(120000,xbHMS,Q2HMS,wt) call hf2(130000,WHMS,Q2HMS,wt) call hf2(480000,thetaHMS,EHMS,wt) call hf1(1110,targx,wt) call hf1(1111,targy,wt) call hf1(1112,phiHMS,wt) call hf2(1113,WHMS,dEHMS,wt) c endif c endif endif endif end