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