Analysis of jumps in a two-level system

// function AnalyzeTLS(inWave, decimation, minLifeTime, verbosity)
//  analyzes jumps in a two-level system
//
//  Given an _inWave_ with the following characteristics:
//      - noisy data with two levels A (ground state) and B > A (excited state) (see #1 in code)
//      - both levels have approximately gaussian distribution (see assumption #2 in code)
//      - the ground state is occupied more often than the excited state (see #3 in code)
//  TwoStateAnalyze finds:
//      - where the jumps between A and B occur (output: EventStart)
//      - how long B is occupied in each case (output: LifeTime)
//      - how high the B level is above the mean A level in each case (output: EventLevel)
//
//  Details and input parameters of AnalyzeTLS:
//      -  first decimates the data to reduce noise
//          - set _decimation_ factor to 1 for no decimation
//      - then analyzes the histogram to find mean levels for A and B
//      - uses these mean levels to find (usually too many) candidate events via FindLevels
//      - uses a lifetime criterion to reject spurious events (too short)
//          - set _minLifeTime_ to the desired value (e.g., 4)
//      - prints results to history (_verbosity_ > 0)
//      - displays histogram and inWave, overlaid with colored events (_verbosity_ > 1)
//
//  Original idea and first version by Xia Liu (Canada)
//  Rewritten by W. Harneit (Berlin, Germany) --- 2010/03/25
//
function AnalyzeTLS(inWave, decimation, minLifeTime, verbosity)
wave inWave         // wave to analyze
variable decimation // factor for downscaling original data (> 1)
variable minLifeTime    // minimum lifetime for event to be counted
variable verbosity      // v = 0: be quiet, v = 1: print results but no graphs
                    // v > 1: also show histogram and colored data
    DFREF saveDF = GetDataFolderDFR()
    NewDataFolder/O/S root:Packages
    NewDataFolder/O/S TLSAnalysis
    DFREF pkgDF = GetDataFolderDFR()
    // step 1: downsample data (eliminates noise)
    string sampWaveName = NameOfWave(inWave)+"_samp"
    Duplicate/O inWave $sampWaveName
    wave sampWave = $sampWaveName
    Resample/DOWN=(decimation) sampWave
    // step 2: create histogram
    string histWaveName = NameOfWave(inWave)+"_hist"
    Make/N=100/O $histWaveName
    wave histWave = $histWaveName
    Histogram/B=1 sampWave, histWave
    if( verbosity > 1 )
        DoWindow/K histGraph            // kill histogram graph if it exists
        Display/N=histGraph histWave
        ModifyGraph log(left)=1
    endif
    // step 3: analyze histogram
    WaveStats/M=1/Q histWave
    variable levelA = V_maxLoc, valA = V_max
    variable wid = levelA - leftx(histWave)     // assume left (#1) peak is at edge (#2), and absolute maximum (#3)
    WaveStats/M=1/Q/R=(levelA + wid, inf) histWave  // 2nd peak should be in this range
    variable levelB = V_maxLoc, valB = V_max
    if( verbosity > 0 )
        printf "mean level(A) = %g (%g), mean level(B) = %g (%g)\r", levelA, valA, levelB, valB
    endif
    // step 4: curve fit histogram  // MAY BE UNNECESSARY
    Duplicate/O histWave, $(histWaveName+"_wt")
    wave histWaveWeight = $(histWaveName+"_wt")
    histWaveWeight = sqrt(histWaveWeight)
    Make/D/O/N=6 histFitPara = {valA, levelA, wid/8, valB, levelB, wid/8}
    FuncFit/NTHR=0/Q DoubleGauss histFitPara histWave /W=histWaveWeight /I=1 /D
    if( verbosity > 1 )
        ModifyGraph rgb[1] = (0,0,44444)    // fit curve in blue
    endif
    // step 5: analyze histogram fit    // MAY BE UNNECESSARY
    wave fitHistWave = $("fit_"+histWaveName)
    WaveStats/M=1/Q/R=(histFitPara[1],histFitPara[4]) fitHistWave
    variable minLevel = V_minLoc, minVal = V_min
    if( verbosity > 0 )
        printf "threshold level = %g (%g)\r", minLevel, minVal
    endif
    if( minVal > 0.1 )
        print "it is recommended to decimate more since the value at threshold is too high"
    endif
    // step 6: use FindLevels to evaluate events
    FindLevels/Q/DEST=Events sampWave, (levelA+levelB)/2
    wave Events
    variable numEvents = numpnts(Events)/2
    if( verbosity > 0 )
        printf "found %g events in total\r", numEvents
    endif
    Redimension/N=(2, numEvents) Events
    SetDataFolder saveDF    // want main output in "current" data folder
    Make/D/O/N=(numEvents) EventStart = Events[0][p], LifeTime = Events[1][p] - EventStart
    Make/D/O/N=(numEvents) EventLevel = mean(inWave, EventStart, EventStart+LifeTime) - levelA
    SetDataFolder pkgDF     // back to our den
    // step 7: reject Events of too short duration (lifetime)
    variable n
    for( n = numEvents - 1; n >= 0; n -= 1 )
        if( LifeTime[n] <= minLifeTime )
            DeletePoints n, 1, LifeTime, EventStart, EventLevel
        endif
    endfor
    variable numValidEvents = numpnts(EventStart)
    if( verbosity > 0 )
        printf "rejected %g events based on lifetime criterion (>%g)\r", numEvents-numValidEvents, minLifeTime
    endif
    // step 8: display original trace colorized (verbosity > 1)
    if( verbosity > 1 )
        DoWindow/K DataGraph
        Display/N=DataGraph/W=(436.5,42.5,831,251) inWave
        string eventWaveName = NameOfWave(inWave)+"_events"
        Duplicate/O inWave $eventWaveName
        wave eventWave = $eventWaveName
        variable p1, p2
        eventWave = 0
        for( n = 0; n < numValidEvents; n += 1 )
            p1 = x2pnt(inWave, EventStart[n])
            p2 = x2pnt(inWave, EventStart[n] + LifeTime[n])
            eventWave[p1,p2] += inWave
        endfor
        eventWave = eventWave > 0 ? eventWave : NaN
        AppendToGraph eventWave
        ModifyGraph rgb[1] = (0,0,44444)    // valid events in blue
    endif
    SetDataFolder saveDF
    // KillDataFolder pkgDF
end

//histogram fitting to double gaussian to estimate the step location and stepsize
Function DoubleGauss(w,x) : FitFunc
    Wave w
    Variable x
     return w[0]*exp(-((x-w[1])^2/(2*(w[2]^2))))+w[3]*exp(-((x-w[4])^2/(2*(w[5]^2))))
End