macro twalk
	
* ..... derivative of the twalk function .....
application COMIS quit
	real function dtw(x)
        common /pawpar/ p(2)
        dtw = -p(1)/(p(2)+x)**2
	return
	end
quit

* ..... twalk function itself (without offset) .....
application COMIS quit
	real function tw(x)
	vector p(2)	
	tw = p(1) / (p(2)+x)
	return
	end
quit

* ..... error on twalk function .....
application COMIS quit
	real function etw(x)
	
	vector p(2)			! fit parameters
	v11 = 28.0			! from MINUIT output
	v12 = 3.24			! from MINUIT output
	v22 = 0.657			! from MINUIT output

	etw = v11/(p(2)+x)**2 + v22*p(1)**2/(p(2)+x)**4
     $      - 2*v12/(p(2)+x)*p(1)/(p(2)+x)**2
	etw = sqrt(etw)
	return
	end
quit

exec sum 0 201
ve/copy par p0

exitm
return

macro sector
| ..... study systematics - twalk function for different sectors
do i=1,4
  hid = 300+[i]
  exec sum [i] [hid]
  ve/copy par p[i]
enddo

opt utit
title 'twalk correction for different sectors' u
null 0. 200. 0. 10.
set fcol 2; ve/copy p0 p; fu/pl tw 0. 200. s; set fcol 1
do i=1,4
  ve/copy p[i] p
  fu/pl tw 0. 200. s
enddo
exitm
return

macro layer
| ..... study systematics - twalk function for different layers .....
do i=5,9
  hid = 300+[i]
  exec sum [i] [hid]
  ve/copy par p[i]
enddo


title 'twalk correction for different layers' u
null 0. 200. 0. 10.

set fcol 2
ve/copy p0 p
fu/pl tw 0. 200 s
set fcol 1
do i=5,9
  ve/copy p[i] p
  fu/pl tw 0. 200 s
enddo
exitm
opt htit

return


macro sum option=0 hid=201

nt/pl 1.twalk.for([option]) 

len  = 20
mrun = 20
ve/cre  d([len]) r
ve/cre ed([len]) r
do iamp=1,[len]
  t  = 0.
  t2 = 0.
  do irun=1,[mrun]
    idx = [iamp] + ([irun]-1)*20
    a  = dtda([idx])
    t  = [t]  + [a]
    t2 = [t2] + [a]*[a]
  enddo
  ve/inp  d([iamp]) $sigma([t]/[mrun])
  ve/inp ed([iamp]) $sigma(sqrt([t2]/[mrun]-([t]/[mrun])**2))
enddo

1d [hid] 'twalk derivative option='//[option] 20 0. 200.
hi/put/co  [hid]  d
hi/put/err [hid] ed
ve/cre  par(2) r 90 10
ve/cre epar(2) r

hi/fit [hid] dtw ! 2 par ! ! ! epar

exitm
return

macro error

title 'Twalk derivative' g
zone 1 1
set fit 111
opt fit
hi/pl 201

| ve/copy par p
| set fcol 1.
| fu/pl tw(x) 0. 200.
| set fcol 4.
| ve/inp p(1) $sigma(par(1)+epar(1))
| fu/pl tw(x) 0. 200. s
| ve/inp p(1) $sigma(par(1)-epar(1))
| fu/pl tw(x) 0. 200. s
| ve/inp p(1) par(1)
| ve/inp p(2) $sigma(par(2)+epar(2))
| fu/pl tw(x) 0. 200. s
| ve/inp p(2) $sigma(par(2)-epar(2))
| fu/pl tw(x) 0. 200. s

title 'Twalk function and its error (cosmic trigger)' g
zone 1 2
opt utit
title 'twalk function' u
null 0. 200. 0. 10.
ve/copy par p
set fcol 2.
set dmod 1; fu/pl tw     0. 200. s
set fcol 1.
set dmod 2; fu/pl tw+etw 0. 200. s
set dmod 2; fu/pl tw-etw 0. 200. s
set dmod 0

set fcol 1.
title 'twalk error' u
null 0. 200. 0. 0.4
set dmod 1; fu/pl etw 0. 200. s
opt htit

exitm
return