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