#!/usr/bin/env python

# ----------------------------------------------------------------------------------- #
#
#  ROOT macro for reading text (ASCII) files and fitting data written.
#
#  Author: Troels C. Petersen (NBI)
#  Email:  petersen@nbi.dk
#  Date:   16th of November 2015
#
# ----------------------------------------------------------------------------------- #

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)

SavePlots = False

# 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  = 0.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 ".txt"):
#   for f in `ls *.txt`; do perl -pi -e 's/\r/\n/g' $f; done

with open("data_RollingBall_example.txt", "r") 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   = float(line[0])          # Measured time
        vol = float(line[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!
        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


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


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

c0.Update()

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

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

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


#---------------------------------------------------------------------------------- 
# 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?
# 
# 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).
# 
# ----------------------------------------------------------------------------------- #