#!/usr/bin/env python

import libMausCpp

from ROOT import *

import sys

if len(sys.argv) > 1:

    fname = sys.argv[1]

else:

    fname = "maus_output.root"

print "Using file: " + fname

# Set up the ROOT file and data pointer

f1 = TFile(fname)

t1 = f1.Get("Spill")

data_ptr = MAUS.Data()

t1.SetBranchAddress("data", data_ptr)

# Create the histograms

c1 = TCanvas("c1", "SciFi Data", 1200, 800)

c1.Divide(5, 2)

c2 = TCanvas("c2", "SciFi Spacepoint x - y Profile", 800, 800)

h200 = TH1D("h200", "Cluster NPE", 100, 0, 50);

h200.GetXaxis().SetTitle("Number of photoelectrons")

h200.SetLineColor(kBlue)

h201 = TH1D("h201", "Clusters per Plane", 3, 0, 3);

h201.GetXaxis().SetTitle("Plane Number")

h201.SetLineColor(kBlue)

h300 = TH1D("h300", "Spacepoint NPE", 100, 0, 150);

h300.GetXaxis().SetTitle("Number of photoelectrons")

h300.SetLineColor(kBlue)

h301 = TH1D("h301", "Spacepoints per Tracker", 2, 0, 2);

h301.GetXaxis().SetTitle("Tracker Number")

h301.SetLineColor(kBlue)

h302 = TH1D("h302", "Spacepoints per Station", 5, 1, 6);

h302.GetXaxis().SetTitle("Station Number")

h302.SetLineColor(kBlue)

h310 = TH1D("h310", "Seed Spacepoint NPE", 100, 0, 150);

h310.GetXaxis().SetTitle("Number of photoelectrons")

h310.SetLineColor(kBlue)

h400 = TH1D("h400", "PR Track Number of Points", 3, 3, 6);

h400.GetXaxis().SetTitle("Number of Points")

h400.SetLineColor(kBlue)

h500 = TH1D("h500", "Final Track p-value", 100, 0, 1);

h500.GetXaxis().SetTitle("p-value")

h500.SetLineColor(kBlue)

h501 = TH1D("h501", "Final Track p-value (zoomed)", 100, 0.01, 1);

h501.GetXaxis().SetTitle("p-value")

h501.SetLineColor(kBlue)

h600 = TH1D("h600", "Trackpoint Pull", 100, -1.5, 1.5);

h600.GetXaxis().SetTitle("Trackpoint Pull")

h600.SetLineColor(kBlue)

h700 = TH2D("h700", "Spacepoint Beam Profile", 75, -180, 180, 75, -180, 180);

h700.GetXaxis().SetTitle("x (mm)")

h700.GetYaxis().SetTitle("y (mm)")

# Loop over all spills

for i in range(1, t1.GetEntries()):

    t1.GetEntry(i) # Update the spill pointed to by data_ptr

    spill = data_ptr.GetSpill()

    if spill.GetDaqEventType() == "physics_event":

        # Loop over recon events in spill

        for j in range(spill.GetReconEvents().size()):

            sfevt = spill.GetReconEvents()[j].GetSciFiEvent()

            # Loop over clusters

            for k in range(sfevt.clusters().size()):

                clus = sfevt.clusters()[k]

                h200.Fill(clus.get_npe())

                h201.Fill(clus.get_plane())

            # Loop over spacepoints

            for k in range(sfevt.spacepoints().size()):

                spoint = sfevt.spacepoints()[k]

                h300.Fill(spoint.get_npe())

                h301.Fill(spoint.get_tracker())

                h302.Fill(spoint.get_station())

                h700.Fill(spoint.get_position().x(), spoint.get_position().y())

            # Loop over pat rec straight tracks

            for k in range(sfevt.straightprtracks().size()):

                trk = sfevt.straightprtracks()[k]

                h400.Fill(trk.get_num_points())

                # Loop over seed spacepoints

                for l in range(trk.get_spacepoints_pointers().size()):

                    spoint = trk.get_spacepoints_pointers()[l]

                    h310.Fill(spoint.get_npe())

            # Loop over pat rec helical tracks

            for k in range(sfevt.helicalprtracks().size()):

                trk = sfevt.helicalprtracks()[k]

                h400.Fill(trk.get_num_points())

                # Loop over seed spacepoints

                for l in range(trk.get_spacepoints_pointers().size()):

                    spoint = trk.get_spacepoints_pointers()[l]

                    h310.Fill(spoint.get_npe())

            for k in range(sfevt.scifitracks().size()):

                trk = sfevt.scifitracks()[k]

                h500.Fill(trk.P_value())

                h501.Fill(trk.P_value())

                # Loop over trackpoints

                for l in range(trk.scifitrackpoints().size()):

                    tpoint = trk.scifitrackpoints()[l]

                    h600.Fill(tpoint.pull())

# Draw the histograms

c1.cd(1)

h200.Draw()

c1.cd(2)

h201.Draw()

c1.cd(3)

h300.Draw()

c1.cd(4)

h301.Draw()

c1.cd(5)

h302.Draw()

c1.cd(6)

h310.Draw()

c1.cd(7)

h400.Draw()

c1.cd(8)

h500.Draw()

c1.cd(9)

h501.Draw()

c1.cd(10)

h600.Draw()

c2.cd()

h700.Draw("colz")

c1.Update()

c2.Update()

# Pause before exit

raw_input("Press Enter to finish...")