macro plot_fluka2 e=2000 theta=20
* 
* Plot parameterizations of Bcal dynamic range
* 05/08/07  ES
*
* set options
*
option ndate
option nbox
*set stat 1111111
set stat 111
option stat
set fit 111
option fit
option grid
*
* plotting options
*
set *
set xmgl 4.
set ymgl 4.
set asiz 0.5
set xlab 2.
set ylab 1.
set xsiz 20.
set xmgl 3.
set ymgl 3.
set ysiz 20.
set gsiz 0.4
* set xwin 0.1
*
* set font definitions to bold roman 
*
set CFON -21
set GFON -21
set LFON -21
set TFON -21
set VFON -21
set txfp -21
set SMGU 0.02
set SMGR 0.02
set CSIZ 0.45
set VSIZ 0.3
set TSIZ 0.3
set YHTI 0.9
set HWID 3.0
set BWID 3.0
* set YHTI 2.
*
* open inputfile
*
* e = 2000
* theta = 20
infile = 'fluka_de_phot_' // [theta] // 'deg_' // [e] // 'mev-shift.hbook'
hi/file 2 [infile] 0 -x
hi/create/title_global 'Energy in Bcal Segments E?[g]!='//[e]//' MeV Shift'
*
zone 2 2
*
* open metafile
*
outfile = 'plot_fluka2_' // [theta] // 'deg_' // [e] // 'mev-shift.ps'
mess ' e=' [e] 'theta=' [theta] ' outfile=' [outfile]
for/file 66 [outfile]
meta 66 -111
* outfile = 'plot_fluka2_' // [theta] // 'deg_' // [e] 'mev.eps'
*for/file 66 [outfile]
*meta 66 -113
*
npts1 = 101
*
xmin = 0 
xmax = 10
ymin = 10
ymax = 70
* option logx 
* option logy
*
npts = [npts1] - 1
sigma j=array([npts1],0#[npts1]-1)
sigma PDE = [xmin] + ([xmax]-[xmin])*j/[npts]
*
*
message 'wait'
* wait
*
label 1 1 ' '
csize = 0.20
igset chhe [csize]
set ndvx 505
set ndvy 505
* hplot/null [xmin] [xmax] [ymin] [ymax]
* hplot/atitle  'Energy(MeV)' 'Counts'
set pmci 2
set plci 2
*
*
* project and plot histograms
*
* Plot deposited energy
*
zone 2 2
hi/del 101
hi/create/1d 101 'Total Energy' 100 0 0.2*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(1,100)
npar = 3
mean = 0.12*[e]
sigma = 0.05*0.12*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a0 = [frac]*$sigma(sqrt([e]/1000))
c0 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' c0=' [c0]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Direct sum'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a0 = ' [a0]
  itx 0.1 0.05 'c0 = ' [c0]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Total Energy' 100 0 0.2*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(1,101)
npar = 3
mean = 0.12*[e]
sigma = 0.05*0.12*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a1 = [a0]
c1 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a1=' [a1] ' c1=' [c1]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=2"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a1 = ' [a1]
  itx 0.1 0.05 'c1 = ' [c1]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Total Energy' 100 0 0.2*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(1,102)
npar = 3
mean = 0.12*[e]
sigma = 0.05*0.12*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a2 = [a0]
c2 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a2=' [a2]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a2 = ' [a2]
  itx 0.1 0.05 'c2 = ' [c2]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Total Energy' 100 0 0.2*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(1,103)
npar = 3
mean = 0.12*[e]
sigma = 0.05*0.12*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a3 = [a0]
c3 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a3=' [a3]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum Norm Groups sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a3 = ' [a3]
  itx 0.1 0.05 'c3 = ' [c3]
exe window#pop
*
* now go for sqrt(Pel*Per)
*
zone 2 2
hi/del 101
hi/create/1d 101 'Sqrt(Pel*Per)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(6,100)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a0 = [frac]*$sigma(sqrt([e]/1000))
c0 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' c0=' [c0]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Direct sum'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a0 = ' [a0]
  itx 0.1 0.05 'c0 = ' [c0]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(Pel*Per)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(6,101)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a1 = [a0]
c1 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a1=' [a1] ' c1=' [c1]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=2"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a1 = ' [a1]
  itx 0.1 0.05 'c1 = ' [c1]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(Pel*Per)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(6,102)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a2 = [a0]
c2 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a2=' [a2]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a2 = ' [a2]
  itx 0.1 0.05 'c2 = ' [c2]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(Pel*Per)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(6,103)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a3 = [a0]
c3 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a3=' [a3]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum Norm Groups sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a3 = ' [a3]
  itx 0.1 0.05 'c3 = ' [c3]
exe window#pop
*
* now go for sqrt(firel*firer)
*
zone 2 2
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(7,100)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a0 = [frac]*$sigma(sqrt([e]/1000))
c0 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' c0=' [c0]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Direct sum'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a0 = ' [a0]
  itx 0.1 0.05 'c0 = ' [c0]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(7,101)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a1 = [a0]
c1 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a1=' [a1] ' c1=' [c1]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=2"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a1 = ' [a1]
  itx 0.1 0.05 'c1 = ' [c1]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(7,102)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a2 = [a0]
c2 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a2=' [a2]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a2 = ' [a2]
  itx 0.1 0.05 'c2 = ' [c2]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR)' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(7,103)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a3 = [a0]
c3 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a3=' [a3]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum Norm Groups sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a3 = ' [a3]
  itx 0.1 0.05 'c3 = ' [c3]
exe window#pop
*
* now go for sqrt(firel*firer) with threshold
*
zone 2 2
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR) Npe"G#10' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(8,100)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a0 = [frac]*$sigma(sqrt([e]/1000))
c0 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' c0=' [c0]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Direct sum'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a0 = ' [a0]
  itx 0.1 0.05 'c0 = ' [c0]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR) Npe"G#10' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(8,101)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a1 = [a0]
c1 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a1=' [a1] ' c1=' [c1]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=2"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a1 = ' [a1]
  itx 0.1 0.05 'c1 = ' [c1]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR) Npe"G#10' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(8,102)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a2 = [a0]
c2 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a2=' [a2]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a2 = ' [a2]
  itx 0.1 0.05 'c2 = ' [c2]
exe window#pop
*
hi/del 101
hi/create/1d 101 'Sqrt(FireL*FireR) Npe"G#10' 100 0 5*[e]
cd //lun2
nt/project 101 77.plot_fluka.f(8,103)
npar = 3
mean = 3*[e]
sigma = 0.05*3*[e]
vec/create par([npar]) R 1000 [mean] [sigma]
hi/fit 101 g ! [npar] par
mean = par(2)
sigma = $sigma(abs(par(3)))
frac = [sigma]/[mean]
a3 = [a0]
c3 = $sigma(sqrt([frac]*[frac] - [a0]*[a0]*1000/[e]))
mess ' mean=' [mean] ' sigma=' [sigma] ' frac=' [frac] ' a3=' [a3]
*
exe window#push
igset chhe 0.3
  itx 0.1 0.35 'Sum Norm Groups sigma=10"Y#'
  itx 0.1 0.25 'Frac = ' [frac]
  itx 0.1 0.15 'a3 = ' [a3]
  itx 0.1 0.05 'c3 = ' [c3]
exe window#pop
*
close 66
*
exitm 
return