macro deadtime run=1

*Read chargetemp.dat including tot_scaler=gscaler_change(510) and fill ntuple 9
nt/cr 9 'Scaler_information' 14 ! ! event ch2s tch tch_p ch2s_p tch_m ch2s_m p_scaler p_trig dtime_p n_scaler n_trig dtime_n t_scaler

nt/read 9 replay/[run].dat

zone 2 2
opt stat
nt/pl 9.n_trig%n_scaler n_trig.gt.0
nt/pl 9.n_trig n_trig.gt.0
nt/pl 9.p_trig%p_scaler n_trig.gt.0
nt/pl 9.p_trig n_trig.gt.0
pic/prin [run]_trig_scalers.gif
zone 1 1
wait
zone 1 2
nt/pl 9.ch2s_m
nt/pl 9.ch2s_p
pic/prin [run]_ch2s.gif
zone 1 1
wait
opt nstat

*Scan coefficients to find best fit close to negative helicity case
nstep=200
stepsize=0.0001
coeff=0.980
first_coeff=[coeff]
final_coeff=[coeff]+[stepsize]*([nstep]-1)

min_trig=0
max_trig=400
min_scal=300.5
max_scal=350.5
*72795
*min_trig=0
*max_trig=1200
*min_scal=1080.5
*max_scal=1200.5
*73014
*min_trig=250
*max_trig=350
*min_trig=0
*max_trig=600
*min_scal=300.5
*max_scal=350.5

*min_scal=1000.5
*max_scal=1200.5
nbin=[max_scal]-[min_scal]
mess Scanning coefficients from [coeff] to [final_coeff]
*mess Fitting [min_trig] < Observed_event_for_2s < [max_trig]
mess Fitting [min_scal] < Recorded_event_for_2s < [max_scal]

opt fit
prof 91 'n_trig vs n_scaler' [nbin] [min_scal] [max_scal] [min_trig] [max_trig]
cut $1 n_scaler.gt.0.and.(n_trig/n_scaler).gt.0.8
nt/proj 91 9.n_trig%n_scaler $1
*abs(ch2s_m-0.070).lt.0.005
*n_trig.gt.[min_trig].and.n_trig.lt.[max_trig]
*abs(n_scaler*0.87218+14.754-n_trig).lt.20
v/cr n_pattern(2) r
hi/fit 91 p1 ! 2 n_pattern
pic/prin [run]_n_fit.gif

v/cr n_x([nbin]) r
get/abs 91 n_x
m=n_pattern(2)
n=n_pattern(1)
sigma fitline = [m]*n_x+[n]

   v/cr p_y([nbin]) r
   v/cr p_yer([nbin]) r
   get/con 91 p_y
   get/err 91 p_yer
   chi2=0.
   nsample=0
   do j=1,[nbin]

      y_obs=p_y([j])
      y_exp=fitline([j])
      sig=p_yer([j])
      if ([y_obs].ne.0) then
         chi2=[chi2]+([y_obs]-[y_exp])*([y_obs]-[y_exp])/([sig]*[sig])
         nsample=[nsample]+1
      endif

   enddo
   chi2_N=[chi2]/([nsample]-2)
   mess chi2 = [chi2], N+2=[nsample], chi2/N=[chi2_N]

*sample=340
*refvalue=[sample]*n_pattern(2)+n_pattern(1)
*mess n_trig = [refvalue] at n_scaler = [sample]
wait

v/cr coeff_list([nstep]) r
v/cr par1([nstep]) r
v/cr par2([nstep]) r
v/cr chi2_ndf([nstep]) r
do i=1,[nstep]

   *Get fitting parameters of reconstructed positive helicity scaler
   prof 92 'p_trig vs p_scaler(recon.)' [nbin] [min_scal] [max_scal] [min_trig] [max_trig]
   cut $2 ([coeff]*t_scaler-n_scaler).gt.0.and.(n_trig/([coeff]*t_scaler-n_scaler)).gt.0.8
   nt/proj 92 9.p_trig%([coeff]*t_scaler-n_scaler) $2
*abs(ch2s_p-0.070).lt.0.005
* p_trig.gt.[min_trig].and.p_trig.lt.[max_trig]
   v/cr p_pattern(2) r
   hi/fit 92 p1 'Q' 2 p_pattern

   v/input coeff_list([i]) [coeff]
   v/input par1([i]) p_pattern(1)
   v/input par2([i]) p_pattern(2)

   v/cr p_y([nbin]) r
   v/cr p_y_er([nbin]) r
   get/con 92 p_y
   get/err 92 p_y_er
   chi2=0.
   nsample=0
   do j=1,[nbin]

      y_obs=p_y([j])
      y_exp=fitline([j])
      sig=p_y_er([j])
      if ([y_obs].ne.0) then
         chi2=[chi2]+([y_obs]-[y_exp])*([y_obs]-[y_exp])/([sig]*[sig])
         nsample=[nsample]+1
      endif

   enddo
   chi2_N=[chi2]/([nsample]-2)
   *v/prin p_y
   *mess nsample=[nsample] chi2_N=[chi2_N]
   v/input chi2_ndf([i]) [chi2_N]

   coeff=[coeff]+[stepsize]

enddo

*Integrate |hel_p_fitting_line-hel_n_fitting_line| to get area between them
m=n_pattern(2)
n=n_pattern(1)
v/cr area_ns([nstep]) r
v/cr chi2Dndf([nstep]) r

do i=1,[nstep]

   a=par2([i])
   b=par1([i])
   interx=([n]-[b])/([a]-[m])
   if ([interx].gt.[min_scal]).and.([interx].lt.[max_scal]) then
      area1=0.5*([a]-[m])*([interx]*[interx]-[min_scal]*[min_scal])+([b]-[n])*([interx]-[min_scal])
      area2=0.5*([a]-[m])*([max_scal]*[max_scal]-[interx]*[interx])+([b]-[n])*([max_scal]-[interx])
      area=[area1]-[area2]
   else
      area=0.5*([a]-[m])*([max_scal]*[max_scal]-[min_scal]*[min_scal])+([b]-[n])*([max_scal]-[min_scal])
      *area=1000.
   endif
   v/input area_ns([i]) [area]
   *mess [min_scal],[interx],[max_scal]
   *mess [a]x+[b] and [m]x+[n]
   *mess area = [area]

   chi2=0.
   do j=[min_scal],[max_scal]
      p_fit=[a]*[j]+[b]
      n_fit=[m]*[j]+[n]
      chi2=[chi2]+([p_fit]-[n_fit])*([p_fit]-[n_fit])/[n_fit]
   enddo
   chi2_N=[chi2]/([nbin]-2)
   v/input chi2Dndf([i]) [chi2_N]
enddo

sigma area_s = abs(area_ns)
*v/prin area_s
*wait

*sigma diff_np = par1 + par2*[sample] - [refvalue]
graph/hpl/null [first_coeff] [final_coeff] 0. 1000
graph/hpl/sym coeff_list area_s [nstep] 24 0.2
graph/hpl/atitle 'Coefficient' 'Area between two fitting lines'
pic/prin [run]_area.gif
wait

graph/hpl/null [first_coeff] [final_coeff] 0. 0.01
graph/hpl/sym coeff_list chi2Dndf [nstep] 24 0.1
graph/hpl/atitle 'Coefficient' 'Chi2/ndf of two fitting lines'
pic/prin [run]_chi2_np_fit.gif
wait

graph/hpl/null [first_coeff] [final_coeff] 0. 10.
graph/hpl/sym coeff_list chi2_ndf [nstep] 24 0.1
graph/hpl/atitle 'Coefficient' 'Chi2/ndf of p data with n fit'
pic/prin [run]_chi2_n_fit_with_p.gif