;;
;; .run readalldata2  ; reads all data
;; .run cuts          ; makes various cuts
;
; (zfindcurves there somewhere, and findafe and such)
;
;; .run zassignr     ; "current best" attempt to assign charges
;;

; in the end this should assign two arrays, Za and Zb for "above" and
; "below" C0 cutoff. Maybe also "Zc" if/when there's a cross-hit
; assignment independent of Zb.
;
; Zb is assigned for all events, whether they're above or below C0
; Za is only assigned for events above C0 (from C1/C0)
;

tt = systime(1)

restore,'par.idl'
print,'Assigning Zb to ALL events:',(tcount=n_elements(data[0,*]))

s = (data[1,*]*par.s2f + data[2,*] )[*]
;; the number s2f is the (apx) ratio of the average s2/s1 value.

c1 = (data[5,*])[*]
afe = par.afe ;/ 10000.d

; from findafe
;afe=[6.132879e+03,-9.573490e+01,8.700731e+00,-1.383504e+00,6.935768e-02]

; self-modified
;afe=[6.176558e+03,-9.479262e+01,8.808393e+00,-1.422378e+00,6.760263e-02]

afe = afe/1d4

; The following is described in Jason's thesis chapter 5.5.2 
; The fit has two components, one which is Z-idependent and deals wit
; the energy-dependence of S, the other half is a Voltz model (with
; weakly energy dependent parameters)
;

zb = c1*0.

p= par.pza

startz = 7.5 ; this is the lowest Z that will ever be assigned.

endz = 35 ; this is only where the histogram ends for plotting; charges 
 ; are assigned until I run out of events, i.e. beyond this number

zstep = 0.01 ; in 0.01e steps

nstep = fix((endz-startz)/zstep) ; in 0.01e steps

h=(hx=dblarr(nstep))

window,retain=2
plot,hx,h,psym=10,xrange=[5,30],/xst,yrange=[0,1.2e4]
        
k1 = where(s gt 0,ko)

i=0l
repeat begin

    z=i*zstep+startz
    z2=z^2
    
    mx = c1/z2
    A=p[0]+mx*p[1]
    B=p[2]+mx*p[3]
    C=p[4]+mx*p[5]
    y3 = poly(mx,afe) * A*z^( B*z + c)

    print,"Z/ko=",z,ko,'  '+string(bytarr(75)+8b),format='($,(a),f5.2,i8,(a))'

    if i ge 1 then begin
        kk = where(s gt y1 and s le y3 and c1 lt 450*z-1500,nk)

        if nk gt 0 then $
          zb[kk] = z + zstep* ((s[kk]-y1[kk])/(y3[kk]-y1[kk])-1)

        if i le nstep then begin
            hx[i-1] = z+zstep/2
            h[i-1] = nk
            oplot,hx[0:i-1],h[0:i-1],psym=10
        endif
                
    endif  

    k1 = where(s gt y3,ko)

    y1 = y3
    i=i+1
endrep  until ko eq 0

plot,hx,h,psym=10,xrange=[8,30],/xst
print

print,'Assigning Za to above-cO events: ',format='($,(A))'
xs = where((rflag and 55) eq 16,tcount)
print,tcount

za = zb*0.
za[xs] = sqrt( c1[xs] - par.k * data[10,xs] ) / par.f

print,'Finished '+systime(0)
print,'Total time elapsed: ',fix(systime(1)-tt)/60,' minutes'
print,'See variables Za and Zb for results.'


end