#!/usr/bin/env python

# ----------------------------------------------------------------------------------- #
#
#  ROOT macro for reading text (ASCII) files and fitting data written.
#
#  Author: Troels C. Petersen (NBI/CERN)
#  Email:  Troels.Petersen@cern.ch
#  Date:   8th of September 2012
#
#  Python author: Lars Egholm Pedersen 
#  Email:  egholm@nbi.dk
#  Date:   8th of September 2012
#
#
# ----------------------------------------------------------------------------------- #

from ROOT import *

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

# Setting of general plotting style:
gStyle.SetCanvasColor(0)
gStyle.SetFillColor(0)
# Setting what to be shown in statistics box:
gStyle.SetOptStat("emr")
gStyle.SetOptFit(1111)


# Set constants and define variables and histograms:
NeventsMax = 56001
Nevents    = 0
lowhigh    = 0      # To begin with, the voltage (vol) is low (i.e. no ball passing!)
limit_vol  = 2.0    # Define SOME LEVEL (your choice!) at which "the ball passes".

t_min = 0.0
t_max = 0.000025 * float(NeventsMax-1)
Hist_vol  = TH1F("Hist_vol" , "Hist_vol;Voltage (arbitrary scale);Number of entries", 6400, 0.0, 4.0)
Hist_tvol = TH1F("Hist_tvol", "Hist_tvol;time (s);Voltage (arbitrary scale)", NeventsMax, t_min-0.0000125, t_max+0.0000125)


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

# IMPORTANT NOTE ABOUT DATA FILES:
# LabView writes out files, where the codes are from DOS, which does not go well
# with Python reading. In order to avoid this problem, run the following command
# in your shell (assuming your data file names have the ending ".dat"):
#   for f in `ls *.dat`; do perl -pi -e 's/\r/\n/g' $f; done

with open("data_RollingBall_V1.txt", "r") as file : 
    for ievent, line in enumerate(file) : #u
        t   = float(line.strip().split()[0]) #Measured time
        vol = float(line.strip().split()[1]) #Measured voltage

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

        # Register when the voltage rises and falls:
        # ------------------------------------------
        # Going from low -> high:
        if ( vol > limit_vol and lowhigh == 0 ) : 
            lowhigh = 1   # OK, now we are in a "high voltage" interval!
            print "  Pass low->high:  %8.6f  %7.4f"%(t, vol)

        # Going from high -> low:
        if ( vol < limit_vol and lowhigh == 1 ) : 
            lowhigh = 0   # OK, now we are (back) in a "low voltage" interval!
            print "  Pass high->low:  %8.6f  %7.4f"%(t, vol)

        Hist_vol.Fill( vol )

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

        Nevents = ievent

    print "  Number of entries in total: %6d"%Nevents

#---------------------------------------------------------------------------------- 
# 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.
#Hist_vol.GetXaxis().SetRangeUser(3.60,3.65)
fitGauss = TF1("fitGauss", "gaus", 0.23, 0.24)   # Gaussian fit (NOTE: range).
fitGauss.SetLineColor(4)
fitGauss.SetLineWidth(2)
Hist_vol.Fit("fitGauss", "RL")
Hist_vol.Draw()

c0.Update()

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

# Drawing time series:
c1 = TCanvas("c1", "", 100, 50, 1200, 500)
Hist_tvol.Draw("e")
c1.Update()
# c1.SaveAs("Fig_RollingBall_Fys2Lab_TimeSeries.png")

#---------------------------------------------------------------------------------- 
# Your analysis here:
#---------------------------------------------------------------------------------- 

g      = 9.8   # Get your own number!
stat_g = 0.1
syst_g = 0.1

print "  Result: g = %6.4f +- %6.4f (stat) +- %6.4f (syst) m/s^-2"%(g, stat_g, syst_g)

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


# ----------------------------------------------------------------------------------- #
# 
# Consideration:
# --------------
# Little intro:
#    Get a feel for the precision on the vol-measurement.
#    (They are in fact discrete with 1600 possible values per unit [0,1]).
#    Measure the uncertainty by plotting the distribution of "low" (or "high") level
#    measurements (i.e. before/between peaks), and look at the RMS. Is the guess good?
#    (It turns out that the precision is better at low level, than at high!)
# 
# Main analysis:
#    Consider the time series - where would you define the voltage at which the ball passes?
#    Does your final result depend on this definition? That would be an obvious systematic
#    uncertainty!
# 
# Suggestion:
#    Since this program provides the crossing times for the ball, and also how they
#    change with your definition of crossing, it might be a good idea to include the
#    subsequent calculations/analysis in the program as well (as has already been
#    indicated).
# 
# ----------------------------------------------------------------------------------- #