#!/usr/bin/env python

# ----------------------------------------------------------------------------------- #
# 
#   Author: Troels C. Petersen (NBI)
#   Email:  petersen@nbi.dk
#   Date:   10th of September 2012
#
# ----------------------------------------------------------------------------------- #

from ROOT import *
import math

Hist_L30cm = TH1F("Hist_L30cm", "Lengths estimates by 30cm ruler", 1000, 0.0, 5.0);

files = ["data_TableMeasurements2012.txt", "data_TableMeasurements2013.txt"]

List30cm = []

# Loop over files and open them
for file in files :
    with open( file, 'r' ) as f : 

        # Read and ignore header lines:
        header1 = f.readline()
        header2 = f.readline()

        # Loop over lines and extract variables of interest
        for line in f:
            line = line.strip()
            columns = line.split()
            if (len(columns) == 4) :
                print columns[0]
                List30cm.append(float(columns[0]))
                Hist_L30cm.Fill(float(columns[0]))

    f.close()



# ---------------------------------------------------------------------------------- 
# Plot result:
# ---------------------------------------------------------------------------------- 
canvas = TCanvas( "canvas", "canvas", 50, 50, 1200, 600 )

Hist_L30cm.SetTitle("Table lengths of Auditorium A")
Hist_L30cm.GetXaxis().SetRangeUser(0.0, 5.0)
Hist_L30cm.GetXaxis().SetTitle("Table lengths (m)")
Hist_L30cm.GetYaxis().SetTitle("Frequency")
Hist_L30cm.SetLineColor(kBlue)
Hist_L30cm.SetLineWidth(2)
Hist_L30cm.Draw("")                     # The option "e" shows errors (Poisson!)

print "  The mean of the histogram is: ", Hist_L30cm.GetMean()


# Fitting histogram (with predefined function):
fit_L30cm = TF1("fit_L30cm", "gaus", Hist_L30cm.GetMean()-0.1, Hist_L30cm.GetMean()+0.1)
fit_L30cm.SetLineColor(kBlue-8)
fit_L30cm.SetLineWidth(3)
Hist_L30cm.Fit("fit_L30cm", "R")



# ---------------------------------------------------------------------------------- 
# Define range we want values from:
# ---------------------------------------------------------------------------------- 

minOK = fit_L30cm.GetParameter(1) - 3.0 * fit_L30cm.GetParameter(2)
maxOK = fit_L30cm.GetParameter(1) + 3.0 * fit_L30cm.GetParameter(2)
print "  I only want value in the range from %6.4f to %6.4f"%(minOK, maxOK)

sum0 = 0.0
sum1 = 0.0
sum2 = 0.0

for l in List30cm :
    if (minOK < l < maxOK) :
        sum0 += 1
        sum1 += l
        sum2 += l*l

print "  The number of measurements that passed the requirement is ", int(sum0)
mean = sum1/sum0
rms = sqrt(sum2/sum0 - mean*mean)
sigma_mean = rms / sqrt(sum0)

print "  The mean is %7.4f +- %7.4f"%(mean, sigma_mean)


# ---------------------------------------------------------------------------------- 
raw_input('Press Enter to exit')