#!/usr/bin/env python

# ----------------------------------------------------------------------------------- #
#
#  ROOT macro for reading text (ASCII) file from LabView for ball rolling down incline.
#
#  Author: Troels C. Petersen (NBI)
#  Email:  petersen@nbi.dk
#  Date:   19th of November 2016
#
# ----------------------------------------------------------------------------------- #

from ROOT import *
import math


# ----------------------------------------------------------------------------------- #
# Define a small simple function:
def sqr( a ) : return a*a
# ----------------------------------------------------------------------------------- #

# Setting of general plotting style:
gStyle.SetCanvasColor(0)
gStyle.SetFillColor(kWhite)

# Setting what to be shown in statistics box:
gStyle.SetOptStat("e")
gStyle.SetOptFit(1111)

SavePlots = False


# Set constants and define variables and histograms:
Nevents    = 0
dt = 1.0 / 40000.0   # The time between measurements (for 40 kHz running)

vol_high = 3.619
vol_low  = 0.234
# limit_vol  = (vol_high+vol_low)/2.0     # Define SOME LEVEL (your choice!) at which "the ball passes"
limit_vol  = 0.1*vol_high + 0.9*vol_low

# Make (empty) lists for the time and voltage measurements:
t = []
vol = []



"""
  t1 = 0.46455    t2 = 0.49395    <t> = 0.47925
  t1 = 0.74392    t2 = 0.76348    <t> = 0.75370
  t1 = 0.95888    t2 = 0.97640    <t> = 0.96764
  t1 = 1.13272    t2 = 1.14702    <t> = 1.13987
  t1 = 1.28242    t2 = 1.29608    <t> = 1.28925

  t1 = 0.46542    t2 = 0.49307    <t> = 0.47925
  t1 = 0.74445    t2 = 0.76295    <t> = 0.75370
  t1 = 0.95955    t2 = 0.97575    <t> = 0.96765
  t1 = 1.13315    t2 = 1.14657    <t> = 1.13986
  t1 = 1.28300    t2 = 1.29563    <t> = 1.28931

  t1 = 0.46377    t2 = 0.49455    <t> = 0.47916
  t1 = 0.74352    t2 = 0.76390    <t> = 0.75371
  t1 = 0.95815    t2 = 0.97715    <t> = 0.96765
  t1 = 1.13243    t2 = 1.14735    <t> = 1.13989
  t1 = 1.28200    t2 = 1.29655    <t> = 1.28927
"""



#---------------------------------------------------------------------------------- 
# File to read:
#---------------------------------------------------------------------------------- 

with open("data_RollingBall_example.txt", "rU") as file : 
    for ievent, line in enumerate(file) :         # Loop over lines in file
        line = line.strip()           # This removes all "codes" that may be put into the file.
        line = line.split()           # This splits the line at all spaces, and makes a list.
        t.append(float(line[0]))      # Measured time (s)
        vol.append(float(line[1]))    # Measured voltage (V)

        # Register when the voltage rises and falls:
        # ------------------------------------------
        # Put your own code here...

        # Write out the first 10 measurements:
        if (ievent < 10) : print"  %2d. events read:  t = %9.6f   vol = %12.8f"%(ievent, t[-1], vol[-1])

    Nevents = len(t)
    print "  Number of entries in total: %6d \n"%Nevents


#---------------------------------------------------------------------------------- 
# Define and fill histograms:
#---------------------------------------------------------------------------------- 

t_min = 0.0
t_max = dt * float(Nevents-1)

# The voltage happens to be measured on a "strange" scale, carefully included in the histogram definition:
Hist_vol  = TH1F("Hist_vol" , "Hist_vol;Voltage (arbitrary scale);Number of entries", 6400, 0.0, 4.0)

# Define the histogram to include all data:
Hist_tvol = TH1F("Hist_tvol", "Hist_tvol;time (s);Voltage (V)", Nevents+1, t_min-dt/2, t_max+dt/2)

for ievent in xrange(Nevents) :
    Hist_vol.Fill(vol[ievent])                        # Voltage measured.
    Hist_tvol.SetBinContent(ievent+1, vol[ievent])    # Voltage measured vs. time.
    Hist_tvol.SetBinError(ievent+1, 0.002)            # Voltage uncertainty on that measurement.
                                                      # Note that bin 0 is the "underflow" bin!


#---------------------------------------------------------------------------------- 
# Determine places of ball passage:
#---------------------------------------------------------------------------------- 

lowhigh = True

t_lowtohigh = []
t_hightolow = []


for ievent in xrange(Nevents) :
    if (vol[ievent] > limit_vol and lowhigh) :
        print t[ievent]
        t_lowtohigh.append( t[ievent] )
        lowhigh = False
    if (vol[ievent] < limit_vol and not lowhigh) :
        print t[ievent]
        t_hightolow.append( t[ievent] )
        lowhigh = True


print
for ipass in xrange( len(t_lowtohigh) ) :
    print "  t1 = %7.5f    t2 = %7.5f    <t> = %7.5f"%(t_lowtohigh[ipass], t_hightolow[ipass], (t_hightolow[ipass]+t_lowtohigh[ipass])/2.0)

print


#---------------------------------------------------------------------------------- 
# Fit and Plot:
#---------------------------------------------------------------------------------- 

# Draw the voltage distribution (to see the precision and to define "low" and "high" voltage):
c0 = TCanvas("c0", "", 50, 20, 600, 400)

Hist_vol.GetXaxis().SetRangeUser(0.225, 0.245)   # Zoom in on the "low voltage" part.
fitGauss = TF1("fitGauss", "gaus", 0.23, 0.24)   # Gaussian fit (NOTE: range).
fitGauss.SetLineColor(kBlue)
fitGauss.SetLineWidth(2)
Hist_vol.Fit("fitGauss", "R")
Hist_vol.Draw()

c0.Update()
if (SavePlots) : c1.SaveAs("Fig_RollingBall_VoltageDistribution.pdf")


#---------------------------------------------------------------------------------- 

gStyle.SetStatX(0.32);    # Position of top right corner.
gStyle.SetStatY(0.85);

# Drawing time series:
c1 = TCanvas("c1", "", 100, 50, 1200, 500)

Hist_tvol.Draw()

c1.Update()
if (SavePlots) : c1.SaveAs("Fig_RollingBall_TimeSeries.pdf")


#---------------------------------------------------------------------------------- 
# Print out the results here:
#---------------------------------------------------------------------------------- 





# ----------------------------------------------------------------------------------- #

# Hold script until you want to exit
raw_input( ' ... Press enter to exit ... ' )