Plot PMT hit time residuals from the RAT DS

time_residuals.py

'''Plot PMT hit time residuals from the RAT DS.

Usage:

    $ python time_residuals.py "filenames_*.root"
    
    With one filename or wildcards (use quotes when using wildcards).
    
The histogram opens in an interactive ROOT window for editing.
'''
import glob
import sys
from rat import ROOT, dsreader

if len(sys.argv) != 2:
    print 'Usage:', sys.argv[0], 'filenames_*.root'

h = ROOT.TH1F('h', ';Time residual (ns);Counts', 500, -50, 250)

for filename in glob.glob(sys.argv[1]):
    dsr = dsreader(filename)
    for ds, run in dsr:
        # Get light path stuff from the Run
        pp = run.GetPMTProp()
        evt = run.GetEffectiveVelocityTime()
        lp = run.GetLightPath()
        dsc, dav, dwa = ROOT.Double(), ROOT.Double(), ROOT.Double()
        for ev in ds.ev:
            # Check for a valid fit
            try:
                fit_valid = ev.GetFitResult('scintFitter').GetValid()
            except Exception:
                fit_valid = False
            if not fit_valid:
                continue

            # Fit info
            fitpos = ev.GetFitResult('scintFitter').GetVertex(0).GetPosition()
            fittime = ev.GetFitResult('scintFitter').GetVertex(0).GetTime()

            # Loop over hits
            for pmt in ev.pmtCal:
                pmtpos = pp.GetPos(pmt.GetID())
                lp.CalcByPosition(fitpos, pmtpos, dsc, dav, dwa)
                tof = evt.CalcByDistance(dsc, dav, dwa)
                tres = pmt.GetsPMTt() - tof - fittime
                h.Fill(tres)

c = ROOT.TCanvas()
h.Draw()
raw_input()