Project

General

Profile

Analysis Code » scifi_run_summary_plots.py

Dobbs, Adam, 29 October 2015 17:18

 
1
#!/usr/bin/env python
2

    
3
import libMausCpp
4
from ROOT import *
5
import sys
6

    
7
if len(sys.argv) > 1:
8
    fname = sys.argv[1]
9
else:
10
    fname = "maus_output.root"
11
print "Using file: " + fname
12

    
13
# Set up the ROOT file and data pointer
14
f1 = TFile(fname)
15
t1 = f1.Get("Spill")
16
data_ptr = MAUS.Data()
17
t1.SetBranchAddress("data", data_ptr)
18

    
19
# Create the histograms
20

    
21
c1 = TCanvas("c1", "SciFi Data", 1200, 800)
22
c1.Divide(5, 2)
23

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

    
26
h200 = TH1D("h200", "Cluster NPE", 100, 0, 50);
27
h200.GetXaxis().SetTitle("Number of photoelectrons")
28
h200.SetLineColor(kBlue)
29

    
30
h201 = TH1D("h201", "Clusters per Plane", 3, 0, 3);
31
h201.GetXaxis().SetTitle("Plane Number")
32
h201.SetLineColor(kBlue)
33

    
34
h300 = TH1D("h300", "Spacepoint NPE", 100, 0, 150);
35
h300.GetXaxis().SetTitle("Number of photoelectrons")
36
h300.SetLineColor(kBlue)
37

    
38
h301 = TH1D("h301", "Spacepoints per Tracker", 2, 0, 2);
39
h301.GetXaxis().SetTitle("Tracker Number")
40
h301.SetLineColor(kBlue)
41

    
42
h302 = TH1D("h302", "Spacepoints per Station", 5, 1, 6);
43
h302.GetXaxis().SetTitle("Station Number")
44
h302.SetLineColor(kBlue)
45

    
46
h310 = TH1D("h310", "Seed Spacepoint NPE", 100, 0, 150);
47
h310.GetXaxis().SetTitle("Number of photoelectrons")
48
h310.SetLineColor(kBlue)
49

    
50
h400 = TH1D("h400", "PR Track Number of Points", 3, 3, 6);
51
h400.GetXaxis().SetTitle("Number of Points")
52
h400.SetLineColor(kBlue)
53

    
54
h500 = TH1D("h500", "Final Track p-value", 100, 0, 1);
55
h500.GetXaxis().SetTitle("p-value")
56
h500.SetLineColor(kBlue)
57

    
58
h501 = TH1D("h501", "Final Track p-value (zoomed)", 100, 0.01, 1);
59
h501.GetXaxis().SetTitle("p-value")
60
h501.SetLineColor(kBlue)
61

    
62
h600 = TH1D("h600", "Trackpoint Pull", 100, -1.5, 1.5);
63
h600.GetXaxis().SetTitle("Trackpoint Pull")
64
h600.SetLineColor(kBlue)
65

    
66
h700 = TH2D("h700", "Spacepoint Beam Profile", 75, -180, 180, 75, -180, 180);
67
h700.GetXaxis().SetTitle("x (mm)")
68
h700.GetYaxis().SetTitle("y (mm)")
69

    
70

    
71
# Loop over all spills
72
for i in range(1, t1.GetEntries()):
73
    t1.GetEntry(i) # Update the spill pointed to by data_ptr
74
    spill = data_ptr.GetSpill()
75
    if spill.GetDaqEventType() == "physics_event":
76
        # Loop over recon events in spill
77
        for j in range(spill.GetReconEvents().size()):
78
            sfevt = spill.GetReconEvents()[j].GetSciFiEvent()
79
            # Loop over clusters
80
            for k in range(sfevt.clusters().size()):
81
                clus = sfevt.clusters()[k]
82
                h200.Fill(clus.get_npe())
83
                h201.Fill(clus.get_plane())
84
            # Loop over spacepoints
85
            for k in range(sfevt.spacepoints().size()):
86
                spoint = sfevt.spacepoints()[k]
87
                h300.Fill(spoint.get_npe())
88
                h301.Fill(spoint.get_tracker())
89
                h302.Fill(spoint.get_station())
90
                h700.Fill(spoint.get_position().x(), spoint.get_position().y())
91
            # Loop over pat rec straight tracks
92
            for k in range(sfevt.straightprtracks().size()):
93
                trk = sfevt.straightprtracks()[k]
94
                h400.Fill(trk.get_num_points())
95
                # Loop over seed spacepoints
96
                for l in range(trk.get_spacepoints_pointers().size()):
97
                    spoint = trk.get_spacepoints_pointers()[l]
98
                    h310.Fill(spoint.get_npe())
99
            # Loop over pat rec helical tracks
100
            for k in range(sfevt.helicalprtracks().size()):
101
                trk = sfevt.helicalprtracks()[k]
102
                h400.Fill(trk.get_num_points())
103
                # Loop over seed spacepoints
104
                for l in range(trk.get_spacepoints_pointers().size()):
105
                    spoint = trk.get_spacepoints_pointers()[l]
106
                    h310.Fill(spoint.get_npe())
107
            for k in range(sfevt.scifitracks().size()):
108
                trk = sfevt.scifitracks()[k]
109
                h500.Fill(trk.P_value())
110
                h501.Fill(trk.P_value())
111
                # Loop over trackpoints
112
                for l in range(trk.scifitrackpoints().size()):
113
                    tpoint = trk.scifitrackpoints()[l]
114
                    h600.Fill(tpoint.pull())
115

    
116
# Draw the histograms
117
c1.cd(1)
118
h200.Draw()
119
c1.cd(2)
120
h201.Draw()
121
c1.cd(3)
122
h300.Draw()
123
c1.cd(4)
124
h301.Draw()
125
c1.cd(5)
126
h302.Draw()
127
c1.cd(6)
128
h310.Draw()
129
c1.cd(7)
130
h400.Draw()
131
c1.cd(8)
132
h500.Draw()
133
c1.cd(9)
134
h501.Draw()
135
c1.cd(10)
136
h600.Draw()
137

    
138
c2.cd()
139
h700.Draw("colz")
140

    
141
c1.Update()
142
c2.Update()
143

    
144
# Pause before exit
145
raw_input("Press Enter to finish...")
(2-2/3)