← Back to branch summary

~vpec/maus/tof_calib_read

« back to all changes in this revision

Viewing changes to workers/.svn/text-base/TimeVSMeanCalib.py.svn-base

  • Committer: tunnell
  • Date: 2010-09-20 10:52:02 UTC
  • Revision ID: tunnell@itchy-20100920105202-ce4w9jm59zvgnsxq
moving stuff from tucs

Show diffs side-by-side

added added

removed removed

Lines of Context:
workers/.svn/text-base/TimeVSMeanCalib.py.svn-base
 
1
# Author: Christopher Tunnell <tunnell@hep.uchicago.edu>
 
2
#
 
3
# March 04, 2009
 
4
#
 
5
 
 
6
from src.GenericWorker import *
 
7
from src.oscalls import *
 
8
import src.MakeCanvas
 
9
import time
 
10
import datetime
 
11
from math import *
 
12
 
 
13
class TimeVSMeanCalib(GenericWorker):
 
14
    "Compute the time stability of a calibration constant"
 
15
 
 
16
    def __init__(self, runType = 'CIS', parameter='calibration', example='EBC_m10_c14', syst=0.007):
 
17
        self.runType = runType
 
18
        self.parameter = parameter
 
19
        self.example = example
 
20
        self.syst= syst
 
21
 
 
22
        self.dir = getPlotDirectory() + '/cis'
 
23
        createDir(self.dir)
 
24
 
 
25
        self.c1 = src.MakeCanvas.MakeCanvas()
 
26
 
 
27
        self.factor = 60*60*24
 
28
 
 
29
    def ProcessStart(self):
 
30
        #self.glo         = ROOT.TGraph()
 
31
        #self.glo_example = ROOT.TGraph()
 
32
        #self.ghi         = ROOT.TGraph()
 
33
        #self.ghi_example = ROOT.TGraph()
 
34
 
 
35
        self.lo = {}
 
36
        self.hi = {}
 
37
 
 
38
        self.lo_example = {}
 
39
        self.hi_example = {}
 
40
 
 
41
    def ProcessStop(self):
 
42
        #da = ROOT.TDatime(1970,01,01,12,00,00);
 
43
        #ROOT.gStyle.SetTimeOffset(da.Convert())
 
44
        ROOT.gStyle.SetNdivisions(304, "y")
 
45
        self.c1.Modified()
 
46
        
 
47
        for all, example, gain in [(self.lo, self.lo_example, 0), (self.hi, self.hi_example, 1)]:
 
48
            gall     = ROOT.TGraphErrors()
 
49
            gall_axis = ROOT.TGraphErrors()
 
50
            gall_sys = ROOT.TGraphErrors()
 
51
            
 
52
            gall_sys.SetFillStyle(3004)
 
53
 
 
54
            i = 0
 
55
            a = all.keys()
 
56
            a.sort()
 
57
            for time in a:
 
58
                calibs = all[time]
 
59
                mean = ROOT.TMath.Mean(len(calibs), array('f', calibs))
 
60
                rms = ROOT.TMath.RMS(len(calibs), array('f', calibs))
 
61
 
 
62
                gall.SetPoint(i, time, mean)
 
63
                gall_axis.SetPoint(i, time, mean)
 
64
                gall_sys.SetPoint(i, time, mean)
 
65
                gall.SetPointError(i, 0, rms/sqrt(len(calibs)))
 
66
                gall_sys.SetPointError(i, 0, self.syst*mean)
 
67
                i += 1
 
68
 
 
69
 
 
70
            gexample = ROOT.TGraphErrors()
 
71
            gexample.SetMarkerStyle(23)
 
72
            gexample.SetMarkerColor(2)
 
73
            gexample_sys = ROOT.TGraphErrors()
 
74
            gexample_sys.SetLineStyle(2)
 
75
            gexample_sys.SetLineWidth(4)
 
76
            #gexample_sys.SetFillColor(ROOT.kGray)
 
77
            #gexample_sys.SetFillStyle(1001)
 
78
 
 
79
            i=0
 
80
            a = example.keys()
 
81
            a.sort()
 
82
 
 
83
            i=0
 
84
            for time in a:
 
85
                calibs = example[time]
 
86
                mean = ROOT.TMath.Mean(len(calibs), array('f', calibs))
 
87
                rms = ROOT.TMath.RMS(len(calibs), array('f', calibs))
 
88
                
 
89
                gexample.SetPoint(i, time, mean)
 
90
                gexample.SetPointError(i, 0, rms/sqrt(len(calibs)))
 
91
                i += 1
 
92
 
 
93
            #gall_axis.Draw("ALP")
 
94
            #self.c1.Update()
 
95
            #self.c1.Modified()
 
96
            #i=0
 
97
            #print self.c1.GetUxmin(), self.c1.GetUxmax()
 
98
            #gexample_sys.SetPoint(0, self.c1.GetUxmin(), gexample.GetMean(2))
 
99
            #gexample_sys.SetPointError(0, 0, self.syst*gexample.GetMean(2))
 
100
            #gexample_sys.SetPoint(1, self.c1.GetUxmax(), gexample.GetMean(2))
 
101
            #gexample_sys.SetPointError(1, 0, self.syst*gexample.GetMean(2))
 
102
            #for time in a:
 
103
            #    gexample_sys.SetPoint(i, time, gexample.GetMean(2))
 
104
            #    gexample_sys.SetPointError(i, 0, self.syst*gexample.GetMean(2))
 
105
            #    i += 1
 
106
                
 
107
            gall_axis.Draw("AP")
 
108
            gall_sys.Draw("A3")
 
109
            gexample_sys.Draw("A3")
 
110
            gall.Draw("P")
 
111
            gexample.Draw("P")
 
112
 
 
113
            if gain ==  0:
 
114
                gall_axis.GetHistogram().GetYaxis().SetTitle("Low gain Calibration (ADC counts/pC)")
 
115
            else:
 
116
                gall_axis.GetHistogram().GetYaxis().SetTitle("High gain Calibration (ADC counts/pC)")
 
117
            gall_axis.GetHistogram().GetXaxis().SetTitle("Time (months)")
 
118
            gall_axis.SetMaximum(gall.GetMean(2)*1.035)
 
119
            gall_axis.SetMinimum(gall.GetMean(2)*0.975)
 
120
            gall_axis.GetHistogram().GetXaxis().SetLabelOffset(1.)
 
121
            gall_axis.Draw("AP")
 
122
            haxis = gall_axis.GetHistogram()
 
123
 
 
124
            ax = gall_sys.GetHistogram().GetXaxis()
 
125
 
 
126
            da = ROOT.TDatime(2008, 8,01,12,00,00).Convert()/self.factor
 
127
            ax.SetBinLabel(ax.FindBin(da), "1 Aug" )
 
128
                
 
129
            da = ROOT.TDatime(2008, 9,01,12,00,00).Convert()/self.factor
 
130
            ax.SetBinLabel(ax.FindBin(da), "1 Sept" )
 
131
 
 
132
            da = ROOT.TDatime(2008,10,01,12,00,00).Convert()/self.factor
 
133
            ax.SetBinLabel(ax.FindBin(da), "1 Oct" )
 
134
 
 
135
            da = ROOT.TDatime(2008,11,01,12,00,00).Convert()/self.factor
 
136
            ax.SetBinLabel(ax.FindBin(da), "1 Nov" )
 
137
 
 
138
            ax.LabelsOption("h")
 
139
            ax.SetLabelSize(0.05)
 
140
 
 
141
            #ax = gall_sys.GetHistogram().GetYaxis()
 
142
             
 
143
            #ax.SetBinLabel(ax.FindBin(80), "80.0" )
 
144
            #ax.SetBinLabel(ax.FindBin(82), "82.0" )
 
145
            #ax.SetBinLabel(ax.FindBin(84), "84.0" )
 
146
 
 
147
            gall_sys.Draw("ALP")
 
148
 
 
149
            self.c1.Clear()
 
150
            gall_sys.GetHistogram().GetYaxis().SetDecimals()
 
151
            haxis.GetYaxis().SetDecimals()
 
152
            haxis.GetXaxis().SetRangeUser(ROOT.TDatime(2008, 8,01,12,00,00).Convert()/self.factor,\
 
153
                                          ROOT.TDatime(2008,11,01,12,00,00).Convert()/self.factor)
 
154
            haxis.GetXaxis().SetNdivisions(503,0)
 
155
            haxis.Draw("AXIS")
 
156
            gall_sys.GetHistogram().Draw("AXISSAME")
 
157
            self.c1.Modified()
 
158
            ##gall_sys.Draw("3")  # too busy with this
 
159
            #tpave2 = ROOT.TBox()
 
160
            #tpave2.SetFillColor(ROOT.kGray)
 
161
            #tpave2.SetFillStyle(1001)
 
162
            #tpave2.SetLineColor(0)
 
163
            #tpave2.SetBorderSize(0)
 
164
            #print self.c1.GetUxmin(),(1-self.syst)*gexample.GetMean(2),self.c1.GetUxmax(),(1+self.syst)*gexample.GetMean(2)
 
165
            #14091.476 81.4252918254 14184.392 
 
166
            #tpave2.DrawBox(14091.476, (1-self.syst)*gexample.GetMean(2),14184.392,(1+self.syst)*gexample.GetMean(2))
 
167
 
 
168
            line = ROOT.TLine()
 
169
            line.SetLineStyle(2)
 
170
            line.SetLineWidth(4)
 
171
            line.DrawLine(14091.476, (1+self.syst)*gexample.GetMean(2), 14184.392, (1+self.syst)*gexample.GetMean(2))
 
172
            line.DrawLine(14091.476, (1-self.syst)*gexample.GetMean(2), 14184.392, (1-self.syst)*gexample.GetMean(2))
 
173
            
 
174
            
 
175
            #gexample_sys.SetPoint(1, self.c1.GetUxmax(), gexample.GetMean(2))
 
176
                        #gexample_sys.SetPointError(1, 0, self.syst*gexample.GetMean(2))   
 
177
            #gexample_sys.Draw("3")
 
178
            self.c1.Modified()
 
179
            haxis.Draw("AXISSAME")
 
180
            self.c1.Modified()
 
181
            gall.Draw("P")
 
182
            gexample.Draw("P")
 
183
 
 
184
            print 'all rms', gall.GetRMS(2)/gall.GetMean(2)
 
185
            print 'example rms', gexample.GetRMS(2)/gexample.GetMean(2)
 
186
 
 
187
            l = ROOT.TLatex()
 
188
            l.SetNDC()
 
189
            l.SetTextFont(72)
 
190
            l.DrawLatex(0.1922,0.867,"ATLAS")
 
191
            
 
192
            l2 = ROOT.TLatex()
 
193
            l2.SetNDC()
 
194
            l2.DrawLatex(0.1922,0.811,"Preliminary")
 
195
 
 
196
            leg = ROOT.TLegend(0.4309045,0.6765734,0.9572864,0.9213287, "", "brNDC")
 
197
            leg.AddEntry(gall,"TileCal average","P")
 
198
            #leg.AddEntry(gall_sys, "detector mean syst. error", "F")
 
199
            leg.AddEntry(gexample,"Typical channel","P")
 
200
            leg.AddEntry(gexample_sys, "\pm 0.7% channel systematic uncertainty", "L") 
 
201
            leg.Draw()
 
202
 
 
203
            tpave = ROOT.TPave()
 
204
            tpave.SetFillColor(0)
 
205
            tpave.SetLineColor(0)            
 
206
            tpave.DrawPave(0.1268844,0.9562937,0.9849246,0.9965035, 0, "brNDC")
 
207
 
 
208
 
 
209
            self.c1.Print("%s/pen_%s_%0.2f.eps" % (self.dir, self.example, gain))
 
210
            self.c1.Print("%s/pen_%s_%0.2f.png" % (self.dir, self.example, gain))
 
211
            self.c1.Print("%s/pen_%s_%0.2f.root" % (self.dir,self.example, gain))
 
212
            
 
213
    
 
214
    def ProcessRegion(self, region):
 
215
        if 'gain' not in region.GetHash():
 
216
            return
 
217
 
 
218
        found = False
 
219
        calibrated = False
 
220
        default = False
 
221
 
 
222
        calib = 0 # calibation constant
 
223
        for event in region.GetEvents():
 
224
            if event.runType == self.runType:
 
225
                #print event.time
 
226
                #print event.time, int(event.getTimeSeconds())
 
227
                if event.data.has_key('goodRegion') and event.data['goodRegion'] and\
 
228
                   event.data.has_key('goodEvent') and event.data['goodEvent']:
 
229
 
 
230
                    # unix/runserver clock starts from 1970, ROOT clock starts from 1995
 
231
                    # so then 788923149 is 25 years in seconds
 
232
                    #print event.getTimeSeconds(), event.getTimeSeconds() -  788923149
 
233
                    time = int(event.getTimeSeconds())/self.factor
 
234
                    
 
235
                    if self.example in region.GetHash():
 
236
                        if 'lowgain' in region.GetHash():
 
237
                            self.lo_example[time] = [event.data['calibration']]
 
238
                        else:
 
239
                            self.hi_example[time] = [event.data['calibration']]
 
240
 
 
241
                    if 'lowgain' in region.GetHash():
 
242
                        if not self.lo.has_key(time):
 
243
                            self.lo[time] = []
 
244
                        self.lo[time].append(event.data['calibration'])
 
245
                    else:
 
246
                        if not self.hi.has_key(time):
 
247
                            self.hi[time] = []
 
248
                        self.hi[time].append(event.data['calibration'])
 
249