macro sim_yields
* 
*	This Macro calculates simc yields in each -t, phi bin.
* 
* input files:   simc ntuples
* 
* outfiles:  yields/yields.[setting].simc
* 
*   macro/trace on
   set *
   opt *
   close 0
   hi/del 0
   ve/del *
   opt date
   opt utit
   tit    chopt=u
   
   opt liny 
   SELNT 1.0 ; set asiz 0.6; set vsiz 0.4
   opt utit
   set mtyp 27 ; set mscf .7 ; set ylab=0.3
   size 30 30 ; set xsiz 30 ; set ysiz 30
   set *fon -60 ; set gsiz 0.3;set csiz 0.35
   
*   for/file 66 psfiles/sim_yields_fpi.ps
   for/file 66 psfiles/sim_yields_fpiM.ps
   meta 66 -111

* Corners of the diamonds.
* upper left:  x1,y1
* lower left:  x2,y2
* upper right: x3,y3
* lower right: x4,y4

*do ifile=1,3
 ifile=1

 exec simc_yields [ifile] +1 0.30 2.20 .34 +0.000 0.004 0.025 6 0.23 2.232 0.28 2.184 0.32 2.214 0.375 2.164
 exec simc_yields [ifile] +1 0.30 2.20 .66 +0.000 0.004 0.025 6 0.23 2.232 0.28 2.184 0.32 2.214 0.375 2.164
 exec simc_yields [ifile] +1 0.30 2.20 .75 +0.000 0.004 0.025 6 0.23 2.232 0.28 2.184 0.32 2.214 0.375 2.164
 
 exec simc_yields [ifile] +1 1.60 3.00 .41 +0.000 0.02 0.08 6 1.24 3.10 1.47 2.97 1.70 3.03 2.04 2.89
 exec simc_yields [ifile] +1 1.60 3.00 .69 +0.000 0.02 0.08 6 1.24 3.10 1.47 2.97 1.70 3.03 2.04 2.89
 exec simc_yields [ifile] +1 1.60 3.00 .82 +0.000 0.02 0.08 6 1.24 3.10 1.47 2.97 1.70 3.03 2.04 2.89

 exec simc_yields [ifile] +1 2.45 3.20 .38 +0.000 0.03 0.15 6 1.95 3.32 2.30 3.18 2.65 3.22 3.05 3.06
 exec simc_yields [ifile] +1 2.45 3.20 .50 +0.000 0.03 0.15 6 1.95 3.32 2.30 3.18 2.65 3.22 3.05 3.06
 exec simc_yields [ifile] +1 2.45 3.20 .71 +0.000 0.03 0.15 6 1.95 3.32 2.30 3.18 2.65 3.22 3.05 3.06

 exec simc_yields [ifile] +1 3.85 3.07 .30 +0.000 0.09 0.27 6 3.15 3.22 3.70 3.05 3.95 3.09 4.60 2.91
 exec simc_yields [ifile] +1 3.85 3.07 .44 +0.000 0.09 0.27 6 3.15 3.22 3.70 3.05 3.95 3.09 4.60 2.91
 exec simc_yields [ifile] +1 3.85 3.07 .57 +0.000 0.09 0.27 6 3.15 3.22 3.70 3.05 3.95 3.09 4.60 2.91
 exec simc_yields [ifile] +1 3.85 3.07 .67 +0.000 0.09 0.27 6 3.15 3.22 3.70 3.05 3.95 3.09 4.60 2.91

 exec simc_yields [ifile] +1 5.25 3.17 .21 +0.060 0.12 0.39 6 4.40 3.33 5.15 3.15 5.30 3.20 6.20 2.98
 exec simc_yields [ifile] +1 5.25 3.17 .39 +0.000 0.12 0.39 6 4.40 3.33 5.15 3.15 5.30 3.20 6.20 2.98
 exec simc_yields [ifile] +1 5.25 3.17 .52 +0.000 0.12 0.39 6 4.40 3.33 5.15 3.15 5.30 3.20 6.20 2.98

 exec simc_yields [ifile] +1 6.00 3.19 .18 +0.370 0.15 0.42 6 5.10 3.35 5.99 3.15 5.95 3.22 7.10 2.98
 exec simc_yields [ifile] +1 6.00 3.19 .30 +0.000 0.15 0.42 6 5.10 3.35 5.99 3.15 5.95 3.22 7.10 2.98
 exec simc_yields [ifile] +1 6.00 3.19 .45 +0.000 0.15 0.42 6 5.10 3.35 5.99 3.15 5.95 3.22 7.10 2.98

*enddo

   close 66
 return
 
* ****************************************************************

 macro simc_yields ifile pol Q2 W epsil th_pq tmin tmax NBt x1 y1 x2 y2 x3 y3 x4 y4

   mess [pol] [Q2] [W] [epsil] [th_pq] [tmin] [tmax] [NBt]
   
   if ([th_pq]>=0.) then
     theta_pq=+$FORMAT([th_pq],f5.3)
   else
     theta_pq=$FORMAT([th_pq],f6.3)
   endif
   
   q2_set=$FORMAT([Q2]*100.01,i3.3)
   w_set=$FORMAT([W]*100.01,i3.3)
   eps_set=$FORMAT([epsil]*100.01,i2.2)
   if ([theta_pq]>=0.) then
     th_pq_set=+$FORMAT([th_pq]*1000,i4.4)
   else
     th_pq_set=$FORMAT([th_pq]*1000,i4.4)
   endif
     
   set ygti 0.6
   set gsiz 0.5
   
 if ([ifile]=1) then
  tit LH+,Q^2!=[Q2],W=[W],[e]=[epsil],[q]?pq!=[th_pq]
 elseif ([ifile]=2) then
  tit LD+,Q^2!=[Q2],W=[W],[e]=[epsil],[q]?pq!=[th_pq]
 elseif ([ifile]=3) then
  tit LD-,Q^2!=[Q2],W=[W],[e]=[epsil],[q]?pq!=[th_pq]
 endif

   pi=3.14159
   rtd=180.0/[pi]
   hi/del 0
   
   NBphi=16		|# of phi bins

* beamfac=1day*24hr*3600sec*70ua/1mC
   Beamfac=6048         |1 PAC-day of beam at 70uA (SIMC is for 1mC)

* APPLY THE CUTS *********

   cuts $0 -

* HMS

   cuts $11 abs(hsdelta)<8.0
   cuts $12 abs(hsxptar)<0.080   | extra sieve slit fitting should improve things
   cuts $13 abs(hsyptar)<0.035

   cuts $10 $11.and.$12.and.$13

* SHMS

   cuts $21 abs(ssdelta)<15.0      
   cuts $22 abs(ssxptar)<0.040
   cuts $23 abs(ssyptar)<0.024

   cuts $20 $21.and.$22.and.$23

* ************************
*                        *
*   Missing mass Cut     *
*                        *
* ************************

   cuts $45 0.875<missmass<1.025

* ***********************
*                       *
*     DIAMOND CUT       * 
*                       *
* ***********************

* Corners of the diamonds.
* upper left:  x1,y1
* lower left:  x2,y2
* upper right: x3,y3
* lower right: x4,y4

* slopes of four sides
* left:
  m21=([y2]-[y1])/([x2]-[x1])
* top:
  m31=([y3]-[y1])/([x3]-[x1])
* bottom:
  m42=([y4]-[y2])/([x4]-[x2])
* right:
  m43=([y4]-[y3])/([x4]-[x3])

   cuts $51 w>([y1]+[m21]*(Q2-[x1])) 
   cuts $52 w<([y1]+[m31]*(Q2-[x1]))
   cuts $53 w>([y2]+[m42]*(Q2-[x2]))
   cuts $54 w<([y3]+[m43]*(Q2-[x3]))
   
   cuts $50 $51.and.$52.and.$53.and.$54

* limit the -t range to (tmin,tmax)

   cuts $60 (t>[tmin]).and.(t<[tmax])

* *********************************

* sum up all the cuts
       
   cuts $66 $10.and.$20.and.$45.and.$50.and.$60

** Process Simc    *****************************

simc:

 setting=[q2_set]_[w_set]_[eps_set]_[th_pq_set]

 if ([ifile]=1) then
  filename=fpi_[setting].hist
 elseif ([ifile]=2) then
  filename=dpl_[setting].hist
 elseif ([ifile]=3) then
  filename=dmn_[setting].hist
 endif

* Get simc normalization factor.
 
 ve/cr tmp(1) r 0.
 sh rm -f data_simc.tmp
* sh grep -a Ncontribute outfiles/[filename] > data_simc.tmp
 sh grep -a Ncontribute outfilesM/[filename] > data_simc.tmp
 ve/read tmp data_simc.tmp '(28x,f11.0)'
 cntmc=$EVAL(tmp(1))
 sh rm -f data_simc.tmp
 sh grep -a normfac outfiles/[filename]>data_simc.tmp      
 ve/read tmp data_simc.tmp '(28x,e17.6)'  
 normfac=$EVAL(tmp(1))
 norm=[normfac]/[cntmc]
 ve/del tmp
 
 mess outfiles/[filename] :  normfac=[normfac] cntmc=[cntmc] norm=[norm]
* wait
 
* fill histo file

 if ([ifile]=1) then
  filename=fpi_[setting].rzdat
 elseif ([ifile]=2) then
  filename=dpl_[setting].rzdat
 elseif ([ifile]=3) then
  filename=dmn_[setting].rzdat
 endif

 close 1
* hi/file 1 worksim/[filename]
 hi/file 1 worksimM/[filename]
* wait
 ntuid=666
 
* run mask

 nt/mask/file mcmdummy N
 nt/loop [ntuid] $66>>mcmdummy(1)
 nt/mask/close mcmdummy
 
 nt/mask/file mcmdummy
 
 1d 1000  phipq [NBphi] 0 360
 call hbarx(1000)

 zone 1 2
 tit    chopt=u
 1d 2000 MissMass 150 -0.20 0.20  
 1d 3000 MissMass 150 -0.20 0.20 
 nt/pl [ntuid].missmass-0.939565 $66*[norm]*[Beamfac]*Weight  ! ! ! EN 2000
 set hwid 5;set hcol 1
 histo/pl 2000 
 ent1=$hinfo(2000,'max')
 ymax=1.05*[ent1]
 set ltyp 1;set hwid 5;set hcol 2
 histo/pl 3000 s
 set lwid 6;set ltyp 1
 dline 0. 0. 0. [ymax] s
* wait
 
 opt fit
 set hcol 1;set hwid 4
 set fcol 1;set fwid 5
* fit same MM ranges as exp_yields.kumac
  hi/fit 2000(-0.0175:0.0250) g
 atit MissMass
* wait

       zone
       set hcol 1
       set pmci 1
       set hwid 6
** without diamond cut
       nt/plot //lun1/[ntuid].W%Q2 $10.and.$20.and.$45
**       wait
       set pmci 2
       set hcol 2
** with diamond cut
       nt/plot //lun1/[ntuid].W%Q2 $66 ! ! ! S 
**       wait   
*       opt liny 
*
       exec window#push
       exec window#pop

** t-phi polar plot
* without diamond cut
       set hcol 1
       set pmci 1
       nt/plot //lun1/[ntuid].sin(phipq)*t%cos(phipq)*t $10.and.$20.and.$45
* with diamond cut
       set hcol 2
       set pmci 2       
       nt/plot //lun1/[ntuid].sin(phipq)*t%cos(phipq)*t $66 ! ! ! S

       exec drawpolar 		| Draw the polar lines
       exec window#push
       exec window#pop
*wait

 1d 5000 t [NBt] [tmin] [tmax]  
 call hbarx(5000)
 nt/pl [ntuid].t $66*[norm]*[Beamfac]*Weight  ! ! ! EN 5000
 histo/pl 5000
 atit '-t'
 
  cdir //lun1

* get simc yields and errors.

  ve/cr y_s([NBt]) r [NBt]*0
  ve/cr yer_s([NBt]) r [NBt]*0
  ve/cr y1_s([NBt]) r [NBt]*0

  hi/get/contents 5000 y_s
  hi/get/error    5000 yer_s
  sigma y1_s=SUMV(y_s)

  mess simc yield: 
  ve/write y_s ' ' '6f8.0'
  ve/write yer_s ' ' '6f8.0'
*  ve/write y1_s ' ' '6f8.0'

*wait

* Loop through t bins.

 zone 2 3
 do j=1,[NBt]
* simc yields and errors.
  ve/cr y_s_[j]([NBphi]) r [NBphi]*0
  ve/cr yer_s_[j]([NBphi]) r [NBphi]*0
   
* cut t into bins
  cut $71 ([tmin]+([tmax]-[tmin])/[NBt]*([j]-1))<t<([tmin]+([tmax]-[tmin])/[NBt]*[j])
* plot histogram.
  tit [j] chopt=u
  nt/pl [ntuid].(phipq*[rtd]) $71*[norm]*mcmdummy(1)*[Beamfac]*Weight ! ! ! EN 1000
* wait
  nh=1000+[j]
  hi/copy 1000 [nh] [f]?[p]q!
  
  set hwid 5
  set hcol 1
  set pmci 1
  set mtyp 20
*  set msiz 0.3

* gh  16.03.01
* very odd error when running LD2 simulations at low epsilon
* two bins near phi=250deg have -NAN counts
 if ([ifile]=1) then
  histo/pl [nh]
 endif	   
* wait   
* get simc yields and errors.
  hi/get/contents 1000 y_s_[j]
  hi/get/error    1000 yer_s_[j]

*  mess simc yield: [j]
*  ve/write y_s_[j] ' ' '16f6.1'
*  ve/write yer_s_[j] ' ' '16f6.1'

  cut $71 -

enddo	|end t bin loop
atit [f]?[p]q!
*wait

hi/del 1000
hi/del [ntuid]
nt/mask/close mcmdummy
close 1

* save data in file

tw=([tmax]-[tmin])/[NBt]
sigma tc1=array([NBt],([tmin]+[tw]/2)#([tmax]-[tw]/2))

pstp=360/[NBphi]
pmn=[pstp]/2
pmx=360-[pstp]/2
sigma phi=array([NBphi],[pmn]#[pmx])

ve/cr one(1) r 1

* if ([ifile]=1) then
*  filename=yields/fpi_yields.[setting].simc
* elseif ([ifile]=2) then
*  filename=yields/dpl_yields.[setting].simc
* elseif ([ifile]=3) then
*  filename=yields/dmn_yields.[setting].simc
* endif

 if ([ifile]=1) then
  filename=yields/fpi_yieldsM.[setting].simc
 elseif ([ifile]=2) then
  filename=yields/dpl_yieldsM.[setting].simc
 elseif ([ifile]=3) then
  filename=yields/dmn_yieldsM.[setting].simc
 endif

ve/write one [filename] ! O

*  ve/write y_s
*  ve/write yer_s
  ve/write y_s,yer_s,tc1 [filename] '(2e13.5,f8.5)' ' '
  ve/write y1_s(6) [filename] '(1e13.5)' ' '

nh_phi=4000

*do j=1,[NBt]
*
*  tw=([tmax]-[tmin])/[NBt]
*  tc=[tmin]+([j]-0.5)*[tw]
*
*  ve/cr tb([NBphi]) r [NBphi]*[tc]
*  ve/write y_s_[j]
*  ve/write yer_s_[j]
*  ve/write y_s_[j],yer_s_[j],phi,tb [filename] '(2e13.5,f7.2,f8.5)' ' '
*  ve/del tb
*
*enddo |NBt

ve/write one [filename] ! C

ve/del one
ve/del y_s,yer_s,tc1
do j=1,[NBt]
  ve/del y_s_[j],yer_s_[j]
enddo

hi/del 0
cut $0 -
*wait

Return




*************************************************
MACRO drawpolar 
* draw the polar graph for t([f])
*************************************************

*null -[tmax] [tmax]  -[tmax] [tmax]
n=12
delta=0.10 
do i=1,[n]
r=[delta]*[i]
graph/prim/fpoint 0. 0. [r]
enddo

line 0 0 0 1
line 0 0 1 0
line 0 0 0 -1
line 0 0 -1 0
line -1 -1 1 1
line -1  1 1 -1

return