#! /usr/bin/awk -i

BEGIN{
    mV=14000;
    Yw=15;
    Gn=12;
    minV = -100;
    maxV = 4000;
    resV = 4*3300/4096.
    resV = 0.838214
    resA = resV
    resB = resV
    offV = 0
    SPARSE = 1
    prescale = 1
    min_dtime = 300
    pmin = 2
    pmax = 1100
    muonthr = 4.96
    maxthr = 150
    Vref=3.3
    name[0] = "A"
    name[1] = "B"
    name[2] = "C"
    name[3] = "D"
    nH = 4
    ph_0 = 1.0468
    ph_1 = 356.5
    if (!U) U="CHAOS"
    @U()
    if (!mthr[0]) mthr[0]=thr[0]
    if (!mthr[1]) mthr[1]=thr[1]
    if (!mthr[2]) mthr[2]=thr[2]
    if (!mthr[3]) mthr[3]=thr[3]
    if (file_epoch) file_epoch_valid = -1
}

function phase(c) { return ph_1*B[c] + ph_0 }

function isE() {
    if (/^EI/) {
	A[0] = $7>=2  ? $6/mV  : -100
	A[1] = $10>=2 ? $9/mV  : -100
	A[2] = $13>=2 ? $12/mV : -100
	A[3] = $16>=2 ? $15/mV : -100
	if ($6<-50*mV || $9<-50*mV || $12<-50*mV || $15<-50*mV) return 0
	return $5 >= min_dtime
    }
    return 0
}

function isEA() {
    if (/^EI/) {
	A[0] = $6/mV
	A[1] = $9/mV
	A[2] = $12/mV
	A[3] = $15/mV
	return $5 >= min_dtime
    }
    return 0
}

function isEE() {
    if (isE()) {
	B[0] = A[0]>0 ?  $8/$6  : -1
	B[1] = A[1]>0 ? $11/$9  : -1
	B[2] = A[2]>0 ? $14/$12 : -1
	B[3] = A[3]>0 ? $17/$15 : -1
	return 1
    }
    return 0
}

function clock() {
    ts = rshift(and(strtonum($3),0x7ffffff0),4)
    if (ts + 0x04000000 < last_ts) 
	ts_roll += 0x08000000
    last_ts = ts
    ts += ts_roll
    return ts
}

function cA(c) { return $(6+3*c) }
function cT(c) { return $(7+3*c) }
function cB(c) { return $(8+3*c) }

function isS() {
    if (/^E/) {
	isEE()
	nE ++
	nS = 0
	TT[0] = cT(0)
	TT[1] = cT(1)
	TT[2] = cT(2)
	TT[3] = cT(3)
	TS = rshift(and(strtonum($3), 0x7ffffff0), 4)
	return 0
    }
    if (/^S/) {
	if (nE==1 && nS==0) {
	    T0 = $2+0;
	    dTS = and(T0, 0x07ffffff) - TS
	    T = T0
	    S0[0] = $3
	    S0[1] = $4
	    S0[2] = $5
	    S0[3] = $6
	    if (dTS >=-2 && dTS<=2) nS = 1
	    nE = 0
	    return nS
	}
	nE=0
	if (nS && $2+0==T+1) {
	    T = $2+0
	    if (nS<16) {
		    S0[0] = (nS*S0[0]+$3)/(nS+1)
		    S0[1] = (nS*S0[1]+$4)/(nS+1)
		    S0[2] = (nS*S0[2]+$5)/(nS+1)
		    S0[3] = (nS*S0[3]+$6)/(nS+1)
	    }
	    nS ++
	    return 1
	}
    }
    return 0
}

function doS(c) {
    if (A[c]>10) {
	if (T==T0) print ""
	print T-T0, $(c+3), S0[c], A[c], B[c], TT[c], phase(c), and(T0, 0x07ffffff), TS, dTS
    }
}

function doSa(c) {
    if (nS==1) print ""
    printf "%d", T-T0
    for (c=0;c<4;c++) printf " %d", $(c+3)
    for (c=0;c<4;c++) printf " %d", $(c+3)-S0[c]
    for (c=0;c<4;c++) printf " %g", A[c]
    for (c=0;c<4;c++) printf " %g", B[c]
    for (c=0;c<4;c++) printf " %g", TT[c]
    printf "\n"
}

function doSS(c, thr) {
    if (A[c]>thr) {
	if (T==T0) print ""
	print ((T-T0)+phase(c)+TT[c])/3.0, ($(c+3)-S0[c])/A[c]*3.3/4.096
    }
}

function degC(hk, R1, R25, B25,  R) {
    R = R1*hk/(4096-hk)
    if (R>0) return B25/(log(R/R25)+B25/298) - 273
    return 222
}

function isH() {
    if (!/^H/) return 0
    Time = $2
    Vcore = $3*Vref/4096
    Vcc   = $4*Vref*3/4096
    Vio   = $5*Vref*2/4096
    Vrpi  = $6*Vref*2/4096
    Tfpga = degC($7, 3.3, 3.3, 4500)
    Tsh   = degC($8, 3.3, 10.0, 3940)
    Tair  = degC($9, 3.3, 3.3, 4500)
    Ibias = $10*Vref/3.3*0.077
    return 1
}

function doH() {
    if (!didH) {
	print "HK time        Vcore Vcc   Vio   Vrpi  Tfpga Tsh  Tair  Ibias"
	didH=1
    }
    printf "HK %u  %.3f %.3f %.3f %.3f  %.1f %.1f %.1f  %.1f\n", \
	Time+file_epoch,  Vcore, Vcc, Vio, Vrpi,  Tfpga, Tsh, Tair,  Ibias
}

function isP(i) {
    if (!/^P/) return 0;
    psn = substr($1,2,1)+0
    if (psn && i && psn != i) return 0;
    Word1 = strtonum($2)
    Word2 = strtonum($3)
    Word3 = strtonum($4)
    Word4 = strtonum($5)
    Data1 = strtonum($6)
    Data2 = strtonum($7)

    if (Data2<10) next
    
    Cal1 = rshift(Word1, 1)
    Cal2 = or(lshift(and(Word3, 0x3f), 6), and(Word4, 0x3f)) 
    Cal3 = rshift(Word4, 6)
    Cal4 = rshift(Word3, 6)
    Cal5 = or(lshift(and(Word1, 1), 10), rshift(Word2, 6))
    Cal6 = and(Word2, 0x3f)

    UT1 = 8*Cal5+20224
    dT = Data2 - UT1
    Temperature = 0.1*(200 + dT*(Cal6+50)/1024)
    
    OFF = Cal2*4 + ((Cal4-512)*dT)/4096
    SENS = Cal1 + (Cal3*dT)/1024 + 24576
    X = (SENS * (Data1-7168))/16384 - OFF
    Pressure = 0.1*(X*10/32 + 2500)
    if (Pressure >= pmin && Pressure <= pmax) return 2
    return 1
}

function isC() {
    if (!/^C/) return 0;
    live = $NF-lastlive 
    lastlive = $NF
    if (live<0) live += 4294967296
    live /= 3e6
    return 1
}

function isCuM() {
    isH()
    isP(1)
    if (!isC() || Time<1633392000) return 0
    if (Pressure<950 || Pressure>1050) return 0
    if (live<59 || live>61) return 0
    for (i=3; i<16; i++) if ($i > 500) return 0
    return 1
}

function doCuM() {
    pc = 1 + (Pressure - 1013)*0.0021
    for (i=3; i<16; i++) $i *= pc
    print $0, Pressure, Time
}

function isMuon13(thr) { return isE() && A[1]>thr && A[3]>thr }

function isMuonA(thr, i) {
    if (thr==0)	thr=muonthr
    n=0
    if (isE()) for (i=0; i<4; i++) if (A[i]>thr) n++;
    return n
}

function do_HIST(C, V,  x) {
    if (V>=minV && V<=maxV) {
	x = int((V-minV)/resV+offV)
	HIST[x][C] += prescale
    }
}
function do_HISTS( C) {
    nHIST ++
    for (C=0; C<4; C++) {
	do_HIST(C, A[C])
    }
}

function phase_HIST(thr, c) {
    if (NR==1) {
	minV = -0.05
	maxV = 0.02
	resV = 0.0005
	iHIST = 1
    }
    if (isEE()) for (c=0; c<4; c++) if (A[c]>thr) do_HIST(c, B[c])
}

function END_HIST() {
    printf "mV"
    for (i=0; i<nH; i++) printf " %s", name[i]
    printf "\n"
    nx = int((maxV-minV)/resV+offV)+1;
    zeros = 1
    for (x=0; x<nx; x++) {
	if (!SPARSE || x in HIST) {
	    if (zeros>1) {
		printf "%g", (xZ+0.5-offV)*resV+minV;
		print " 0 0 0 0"
	    }
	    zeros = 0
	    printf "%g", (x+0.5-offV)*resV+minV;
	    for (i=0; i<nH; i++) {
		printf " %d", HIST[x][i]+0;
	    }
	    print ""
	} else if (zeros) {
	    zeros ++
	    xZ = x
	} else {
	    printf "%g", (x+0.5-offV)*resV+minV;
	    print " 0 0 0 0"
	    zeros = 1
	}
    }
}

function END_HIST_INT() {
    printf "mV"
    for (i=0; i<nH; i++) printf " %s", name[i]
    printf "\n"
    nx = int((maxV-minV)/resV)+1;
    zeros = 1
    last = -1;
    for (x=0; x<nx; x++) {
	if (x in HIST) {
	    printf "%g", x*resV+minV;
	    for (i=0; i<nH; i++) {
		HIST[x][i] += HIST[last][i]
		printf " %d", HIST[x][i];
	    }
	    last = x
	    print ""
	}
    }
}

function mADAM() {
    D1 = 0
    D2 = 3
    HO = 1
    VE = 2
}
function MDDM() {
    D1 = 1
    D2 = 3
    HO = 0
    VE = 2
    muon1 = D1
    muon2 = D2
    validl2 = "0x04 0x08 0x4d 0x10 0x95 0x18 0xdd 0x20 0x26 0x69 0x6f 0xb1 0xb7 0xf9 0xff"
}
function CsI() {
    AA=2;   name[2] = "AA"
    B1=3;   name[3] = "B1"
    B2=0;   name[0] = "B2"
    CC=1;   name[1] = "CC"
    thr[AA] = 8
    thr[B1] = 8
    thr[B2] = 8
    thr[CC] = 8
    muon1 = AA
    muon2 = CC
    resV /= 1
}
function TANOS() {
    AA=2;   name[2] = "AA"
    B1=3;   name[3] = "B1"
    B2=0;   name[0] = "B2"
    CC=1;   name[1] = "CC"
    thr[AA] = 8
    thr[B1] = 8
    thr[B2] = 8
    thr[CC] = 8
    muon1 = AA
    muon2 = CC
    resV /= 1
}

function CHAOS() {
    B1=2;   name[2] = "B1"
    B2=3;   name[3] = "B2"
    AA=0;   name[0] = "AA"
    CC=1;   name[1] = "CC"
    thr[AA] = 7
    thr[B1] = 4
    thr[B2] = 4
    thr[CC] = 7
    muon1 = AA
    muon2 = CC
    resV /= 4
    ph_0 = 1.08133
    ph_1 = 46.1029
    Bcut=0.8
    Bcof=2
    resA = 1/20.
    resB = 1/20.
    LOG10=log(10)
    isMUONf = "CHAOSmuon"
}

function CSRT() {
    B1=2;   name[2] = "B1"
    B2=3;   name[3] = "B2"
    AA=0;   name[0] = "AA"
    CC=1;   name[1] = "CC"
    thr[AA] = 6
    thr[B1] = 6
    thr[B2] = 6
    thr[CC] = 6
    muon1 = AA
    muon2 = CC
    resV /= 4
    ph_0 = 1.08133
    ph_1 = 46.1029
    Bcut=0.8
    Bcof=2
    resA = 1/20.
    resB = 1/20.
    LOG10=log(10)
    isMUONf = "CSRTmuon"
}

function CHAOSmuon() {
    if (!isE()) return 0
    if (A[AA] < 6) return 0
    if (A[CC] < 6) return 0
    if (A[B1]*Bcut-Bcof > A[B2]) return 0
    if (A[B2]*Bcut-Bcof > A[B1]) return 0
    return 1
}

function CSRTmuon() {
    if (!isE()) return 0
    if (A[AA] < 6 && A[CC] < 6) return 0
    if (A[B1] < 6 && A[B2] < 6) return 0
    return 1
}

function CHAOSfish() {
    if (!CHAOSmuon()) return 0
    BB = (A[B1]+A[B2])*18000/39.
    A1 = A[AA]*394/55.
    A2 = A[CC]*394/55.
    if (A1 >100*A2) return 0
    if (A2 >100*A2) return 0
    if (BB < 1000 || BB > 10e6) return 0
    RR = A1/A2
    return 1
}

function doCHAOSfish() {
    do_2DHIST(log(BB)/LOG10, log(RR)/LOG10)
    if (n2DHIST==1) {
	minA =  3/resA
	maxA =  7/resA+1
	minB = -2/resA
	maxB =  2/resA+1
    }
}

function DuToit() {
    D1=0;   name[D1] = "D1"
    D2=1;   name[D2] = "D2"
            name[2]  = "nn"
    Dx=3;   name[Dx] = "Dx"
    thr[0] = 15
    thr[1] = 15
    thr[2] = 15
    thr[3] = 15
    muon1 = D1
    muon2 = D2
    mthr[0] = 20
    mthr[1] = 38
}

function THREED() {
    thr[0] = 10
    thr[1] = 10
    thr[2] = 10
    thr[3] = 10
    muon1 = 0
    muon2 = 2
    name[0] = "A"
    name[1] = "x"
    name[2] = "B"
    name[3] = "C"
    name[4] = "Bc"
    name[5] = "Cc"
    nH = 6
}

function uM() {
    A1=0; name[A1] = "A1"; thr[A1] = 12
    A2=1; name[A2] = "A2"; thr[A2] = 12
    B1=2; name[B1] = "B1"; thr[B1] = 12
    B2=3; name[B2] = "B2"; thr[B2] = 12

    mipsgain[A1]=20.1587
    mipsgain[A2]=22.6824
    mipsgain[B1]=21.7788
    mipsgain[B2]=26.1791

    mipsCal = 5 # MeV
    scale[A1]=mipsCal/mipsgain[A1]
    scale[A2]=mipsCal/mipsgain[A2]
    scale[B1]=mipsCal/mipsgain[B1]
    scale[B2]=mipsCal/mipsgain[B2]
    uMthres = mipsCal*1.0
    uMratio = 3
}

function do_ThreeD() {
    do_HISTS();
    if (A[0]>thr[0]) {
	do_HIST(4, A[2])
	do_HIST(5, A[3])
    }
}

function uM_HIST() {
    nHIST ++
    if (isuMlight(A1,A2)) {
	 do_HIST(B1, A[B1])
	 do_HIST(B2, A[B2])
    }
    if (isuMlight(B1,B2)) {
	 do_HIST(A1, A[A1])
	 do_HIST(A2, A[A2])
    }
}

function do_2DHIST(A, B) {
    iA = int(A/resA+1000)-1000
    iB = int(B/resB+1000)-1000
    if (!n2DHIST) {
	minA = iA
	maxA = iA
	minB = iB
	maxB = iB
    } else {
	if (iA > maxA) maxA = iA
	if (iA < minA) minA = iA
	if (iB > maxB) maxB = iB
	if (iB < minB) minB = iB
    }
    n2DHIST++;
    HIST2D[iA,iB] += 1
}

function is2D(i,j) {
    if (!isE()) return 0
    if (A[i] <= -100 || A[j] <= -100) return 0
    Ai = A[i]
    Aj = A[j]
    Ei = Ai*scale[i]
    Ej = Aj*scale[j]
    return 1
}

function isuMlight(i,j) {
    if (!is2D(i,j)) return 0
    return Ei+Ej > uMthres && Ei*uMratio > Ej && Ej*uMratio > Ei
}

function print_2DHIST(fn) {
    if (max_A && max_A < maxA*resA) maxA = int(max_A/resA)
    if (min_A && min_A > minA*resA) minA = int(min_A/resA)
    if (max_B && max_B < maxB*resB) maxB = int(max_B/resB)
    if (min_B && min_B > minB*resB) minB = int(min_B/resB)
    printf "%u", maxA-minA+1 > fn
    for (iA=minA; iA<=maxA; iA++) printf " %g", iA*resA > fn
    printf "\n" > fn
    for (iB=minB; iB<=maxB; iB++) {
	printf "%g", iB*resB > fn
	for (iA=minA; iA<=maxA; iA++) printf " %u", HIST2D[iA,iB]+0 > fn
	printf "\n" > fn
    }
    delete HIST2D
    n2DHIST = 0
}

END { if (n2DHIST) print_2DHIST("/dev/stdout") }
    
function isMuon() {
    if (isMUONf) return @isMUONf()
    return isE() && A[muon1]>mthr[muon1] && A[muon2]>mthr[muon2]
}

function isXray(  i) {
    N = 0
    if (!isEA()) return 0
    for (i=0; i<4; i++) if (A[i]>thr[i] || A[i]<-50) { N++; C=i }
    return N==1 && $(7+3*C)>=2
}

function incrRATE() {
    if (A[D1] < muonthr && A[D2] < muonthr) next
    if (A[0] > maxthr || A[1] > maxthr || A[2] > maxthr || A[3] > maxthr) next
    if (A[D1] >= muonthr && A[D2] >= muonthr) R[0] += 1
    if (A[HO] >= muonthr && A[VE] <  muonthr) R[1] += 1
    if (A[VE] >= muonthr && A[HO] <  muonthr) R[2] += 1
}

function isGPS() {
    isH()
    if (!/GPRMC/) return 0
    N = split($NF, GPRMC, /[,$]/)
    i = 1
    while (i<=N && GPRMC[i] != "GPRMC") i++
    if (i+9 > N) return 0
    TT=GPRMC[i+1]
    DD=GPRMC[i+9]
    TT=substr(TT,1,2)" "substr(TT,3,2)" "substr(TT,5,2)
    DD="20"substr(DD,5,2)" "substr(DD,3,2)" "substr(DD,1,2)
    GPSTime=mktime(DD" "TT, "UTC")
    valid = GPRMC[i+2]=="A"
    if (Time != GPSHTime && Time) {
	if (!file_epoch_valid || valid) {
	    file_epoch = GPSTime-Time+1
	    file_epoch_valid = valid
	}
	GPSHTime = Time
    }
    return 1 + valid
}

END {
    if (iHIST) { END_HIST_INT() }
    else if (nHIST) {END_HIST() }
}

function isStrangeL2() {
    if (!isE()) return 0
    if (index(validl2, $4)) return 1
    return 2
}