ATLAS Offline Software
BandUtils.py
Go to the documentation of this file.
1 #!/bin/python
2 
3 # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
4 # -*- coding: utf-8 -*-
5 
6 import sys
7 sys.argv = [sys.argv[0], '-b']
8 
9 from ROOT import TGraphAsymmErrors, Double, TGraphErrors, TCanvas
10 import atlasStyleMacro
11 from math import fabs, sqrt
12 
13 def makeBand (middleGraph, lowerGraph, upperGraph, noorder) :
14 
15  bandsGraph = TGraphAsymmErrors()
16  bandsGraph.SetName(middleGraph.GetName())
17  bandsGraph.SetMarkerStyle(middleGraph.GetMarkerStyle());
18  bandsGraph.SetMarkerColor(middleGraph.GetMarkerColor());
19  bandsGraph.SetLineColor(middleGraph.GetLineColor());
20 
21  x1 = Double(0)
22  x2 = Double(0)
23  x3 = Double(0)
24  y = []
25  y.append(Double(0))
26  y.append(Double(0))
27  y.append(Double(0))
28 
29  tmp = Double(0)
30 
31  if lowerGraph.GetN() != upperGraph.GetN() or lowerGraph.GetN() != middleGraph.GetN() :
32  print "Plots don't have the same number of points!"
33  print "Lower: ", lowerGraph.GetN()
34  print "Upper: ", upperGraph.GetN()
35  print "Middle: ", middleGraph.GetN()
36  return 0
37 
38  #again a hack to forget about the 1st point
39  for iPoint in xrange(1,lowerGraph.GetN()) :
40  middleGraph.GetPoint(iPoint, x1, y[0])
41  lowerGraph.GetPoint(iPoint, x1, y[1])
42  upperGraph.GetPoint(iPoint, x1, y[2])
43 
44  if (iPoint == lowerGraph.GetN()-1) : x2 = x1
45  else : upperGraph.GetPoint(iPoint+1, x2, buf)
46 
47  if (iPoint == 0) : x3 = x1
48  else : upperGraph.GetPoint(iPoint-1, x3, buf)
49 
50  tmp = Double(0)
51  yup = Double(0)
52  yce = Double(0)
53  ydn = Double(0)
54 
55  if noorder == 0 :
56  yce = y[0]
57  ydn = y[1]
58  yup = y[2]
59 
60  elif noorder == 1 :
61  tmp = y[2]
62  if y[1] < y[2] :
63  y[2] = y[1]
64  y[1] = tmp
65 
66  yce=y[0]
67  ydn=y[1]
68  yup=y[2]
69 
70  else :
71  for p in [0,1,2] :
72  for q in [0,1,2] :
73  if y[q] > y[q+1] :
74  tmp = y[q+1]
75  y[q+1]=y[q]
76  y[q]=tmp
77 
78  ydn=y[0]
79  yce=y[1]
80  yup=y[2]
81 
82  bandsGraph.SetPoint(iPoint, x1, yce)
83 
84  ex0 = Double(lowerGraph.GetErrorX(iPoint))
85 
86  binwl = Double(0)
87  binwh = Double(0)
88 
89  if (ex0==0) :
90  binwl=(x1-x3)/2.;
91  binwh=(x2-x1)/2.;
92  if binwl==0. : binwl= binwh;
93  if binwh==0. : binwh= binwl;
94  else :
95  binwl = ex0
96  binwh = ex0
97 
98  dxl= Double(yce-ydn)
99  dxh= Double(yup-yce)
100 
101  if noorder==0 :
102  if dxl<0 :
103  tmp=-dxl
104  dxl=dxh
105  dxh=tmp
106 
107  bandsGraph.SetPointError(iPoint, binwl, binwh, dxl, dxh)
108 
109  return bandsGraph
110 
111 
112 
113 def getLargestVariationGraph (graph1, graph2) :
114 
115  # Protection
116  if graph1.GetN() != graph2.GetN() :
117 
118  print "Plots don't have the same number of points!"
119  print "Error: ", graph1.GetN()
120  print "Middle: ", graph2.GetN()
121  return 0
122 
123  largestGraph = graph1.Clone()
124 
125  for iPoint in xrange(0, graph1.GetN()) :
126 
127  # Retrieve points
128  graph1PointX = Double(0)
129  graph1PointY = Double(0)
130  graph1ErrorX = graph1.GetErrorX(iPoint)
131  graph1ErrorY = graph1.GetErrorY(iPoint)
132 
133  graph1.GetPoint(iPoint,graph1PointX,graph1PointY)
134 
135  graph2PointX = Double(0)
136  graph2PointY = Double(0)
137  graph2ErrorX = graph2.GetErrorX(iPoint)
138  graph2ErrorY = graph2.GetErrorY(iPoint)
139 
140  graph2.GetPoint(iPoint,graph2PointX,graph2PointY)
141 
142  # More protection - who needs it!! we are young and we listen to malmsteen and we are obviously fearless
143  if graph1PointX - graph2PointX > 0.0001 :
144  print "Points are not ordered!"
145  print "graph1: ", graph1PointX
146  print "graph2: ", graph2PointX
147  return 0
148 
149  largestPoint = Double(0)
150 
151  #worst case scenario: one of the two is _completely_ wrong
152  if fabs(fabs(1-graph1PointY) - fabs(1-graph2PointY)) > 1000 :
153  print "found bad data point: ", graph1PointX, graph2PointY
154  if fabs(1-graph1PointY) > 100 :
155  largestPoint = graph2PointY
156  largestPointErrorX = graph2ErrorX
157  largestPointErrorY = graph2ErrorY
158  if fabs(1-graph2PointY) > 100 :
159  largestPoint = graph1PointX
160  largestPointErrorX = graph1ErrorX
161  largestPointErrorY = graph1ErrorY
162 
163 
164  elif fabs(1-graph1PointY) > fabs(1-graph2PointY) :
165  largestPoint = graph1PointY
166  largestPointErrorX = graph1ErrorX
167  largestPointErrorY = graph1ErrorY
168 
169  else :
170  largestPoint = graph2PointY
171  largestPointErrorX = graph2ErrorX
172  largestPointErrorY = graph2ErrorY
173 
174  print "newpoint:", largestPoint, largestPointErrorY
175 
176  #set the point in the new plot
177  largestGraph.SetPoint(iPoint, graph1PointX, largestPoint)
178  largestGraph.SetPointError(iPoint, largestPointErrorX, largestPointErrorY)
179 
180  return largestGraph
181 
182 #this function returns a plot that contains the variation of a graph wrt 1
183 def getVariationGraph (graph) :
184 
185  varGraph = graph.Clone()
186 
187  for iPoint in xrange(0, graph.GetN()) :
188 
189  # Retrieve points
190  graphPointX = Double(0)
191  graphPointY = Double(0)
192  graphErrorX = graph.GetErrorX(iPoint)
193  graphErrorY = graph.GetErrorY(iPoint)
194 
195  graph.GetPoint(iPoint,graphPointX,graphPointY)
196 
197  newPoint = Double(0)
198 
199  newPoint = fabs(1-graphPointY)
200  newPointErrorX = graphErrorX
201  newPointErrorY = graphErrorY
202 
203  #set the point in the new plot
204  varGraph.SetPoint(iPoint, graphPointX, newPoint)
205  varGraph.SetPointError(iPoint, newPointErrorX, newPointErrorY)
206 
207  return varGraph
208 
209 
210 
211 
212 def makeSingleBandPlot (middleGraph, errorGraphs) :
213 
214  bandsGraph = []
215 
216  for iPoint in xrange(0, middleGraph.GetN()) :
217 
218  # Retrieve middle graph for (x,y) coordinates
219  dataPointX = Double(0)
220  dataPointY = Double(0)
221  middleGraph.GetPoint(iPoint,dataPointX,dataPointY)
222  #print "DATAPOINT: ", dataPointX
223 
224  # Prepare empty arrays
225  errorX = Double(0)
226  errorY = []
227 
228  # Loop over all ErrorGraphs
229  index = 0
230  for errorGraph in errorGraphs:
231  #print "NEW ERRORGRAPH: ", errorGraph.GetName()
232  if iPoint==0 :
233  bandsGraph += [TGraphAsymmErrors()]
234  #print "Point ", iPoint, " Index: ", index
235  thisBandsGraph = bandsGraph[index]
236  thisBandsGraph.SetName(middleGraph.GetName())
237  thisBandsGraph.SetMarkerStyle(middleGraph.GetMarkerStyle());
238  thisBandsGraph.SetMarkerSize(middleGraph.GetMarkerSize());
239  thisBandsGraph.SetMarkerColor(middleGraph.GetMarkerColor());
240  thisBandsGraph.SetLineColor(middleGraph.GetLineColor());
241 
242  index+=1
243 
244  # Protection
245  if errorGraph.GetN() != middleGraph.GetN() :
246  print "Plots don't have the same number of points!"
247  print "Error: ", errorGraph.GetN()
248  print "Middle: ", middleGraph.GetN()
249  return 0
250 
251  # Retrieve errors
252  errorPointX = Double(0)
253  errorPointY = Double(0)
254  errorGraph.GetPoint(iPoint,errorPointX,errorPointY)
255 
256  # More protection - who needs it!!
257  if dataPointX - errorPointX > 0.0001 :
258  print "Points are not ordered!"
259  print "dataPointX: ", dataPointX
260  print "errorX: ", errorPointX
261  return 0
262 
263  errorX = middleGraph.GetErrorX(iPoint)
264  errorY += [fabs(1-errorPointY)]
265 
266  # Add things so far in quadrature
267  sumErrorY = Double(0)
268  for a in errorY:
269  #print "adding ", a
270  sumErrorY += a*a
271 
272  #print '***REMOVED*** ', sumErrorY
273  sumErrorY = sqrt(sumErrorY)
274  #print ***REMOVED***, sumErrorY
275 
276  if( iPoint < middleGraph.GetN()):
277  thisBandsGraph.SetPoint(iPoint, dataPointX, 0)
278  thisBandsGraph.SetPointError(iPoint, errorX, errorX, sumErrorY, sumErrorY)
279  thisBandsGraph.SetPointEYlow(iPoint, 0)
280 
281  return bandsGraph[-1:]
282 
283 def makeSingleBandPlot_error (middleGraph, varGraphs) :
284 
285  bandsGraph = []
286  #contains all the stuff necessary for propagating errors in each plot
287  propagationVarY = []
288 
289  for iPoint in xrange(0, middleGraph.GetN()) :
290  #print "NEW DATA POINT"
291  # Retrieve middle graph for (x,y) coordinates
292  dataPointX = Double(0)
293  dataPointY = Double(0)
294  middleGraph.GetPoint(iPoint,dataPointX,dataPointY)
295 
296  # Prepare empty arrays
297  varX = Double(0)
298  varY = []
299  #contains all the errors
300  errorY = []
301 
302  # Loop over all ErrorGraphs
303  index = 0
304  for iGraph, varGraph in enumerate(varGraphs):
305 
306  #print "NEW GRAPH"
307 
308  if iPoint==0 :
309  bandsGraph += [TGraphAsymmErrors()]
310  #print "Point ", iPoint, " Index: ", index
311  thisBandsGraph = bandsGraph[index]
312  thisBandsGraph.SetName(middleGraph.GetName())
313  thisBandsGraph.SetMarkerStyle(middleGraph.GetMarkerStyle());
314  thisBandsGraph.SetMarkerSize(middleGraph.GetMarkerSize());
315  thisBandsGraph.SetMarkerColor(middleGraph.GetMarkerColor());
316  thisBandsGraph.SetLineColor(middleGraph.GetLineColor());
317 
318  index+=1
319 
320  # Protection
321  if varGraph.GetN() != middleGraph.GetN() :
322  print "Plots don't have the same number of points!"
323  print "Error: ", varGraph.GetN()
324  print "Middle: ", middleGraph.GetN()
325  return 0
326 
327  # Retrieve vars
328  varPointX = Double(0)
329  varPointY = Double(0)
330  varGraph.GetPoint(iPoint,varPointX,varPointY)
331  varErrorY = varGraph.GetErrorY(iPoint)
332  #print "varErrorY:", varErrorY
333 
334  # More protection - who needs it!!
335  if dataPointX - varPointX > 0.0001 :
336  print "Points are not ordered!"
337  print "dataPointX: ", dataPointX
338  print "varX: ", varPointX
339  return 0
340 
341  varX = middleGraph.GetErrorX(iPoint)
342  varY += [fabs(1-varPointY)]
343  errorY += [varErrorY]
344  #print "errorY:", errorY
345 
346  # Add things so far in quadrature
347  sumVarY = Double(0)
348  propagationY = Double(0)
349 
350  for var in varY:
351  sumVarY += var*var
352 
353  sumVarY = sqrt(sumVarY)
354 
355  #the square root of sumVarY is now the denominator of the propagation for all plots
356  #now need to multiply it by the single variation
357  #and store it in a vector for further use
358  for var, error in zip(varY, errorY) :
359  if sumVarY > 0.0000001 : propagationY += var*var / sumVarY * error
360 
361  #print "squared data point: ", sumVarY
362  #print "error: ", propagationY
363 
364  propagationY = sqrt(propagationY)
365 
366  #print "data point: ", sumVarY
367  #print "error: ", propagationY
368 
369  #this will happen once per graph, but only the last time will be the one with all the graphs taken into account - that is the one that will stay!
370  if( iPoint < middleGraph.GetN()):
371  thisBandsGraph.SetPoint(iPoint, dataPointX, 0)
372  thisBandsGraph.SetPointError(iPoint, varX, varX, sumVarY, sumVarY)
373  thisBandsGraph.SetPointEYlow(iPoint, 0)
374 
375  #at the last graph, save the error in array
376  if (iGraph == len(varGraphs)-1) :
377  propagationVarY += [propagationY]
378 
379  #transform the first bandsGraph into a TH1F
380  #print bandsGraph[-1:]#this is a TGraph in an array
381 
382  #pick up the errors calculated previously and assign them to the TH1F
383  #print "varY:", varY
384  #print "errorY:", errorY
385  #print "propagationVarY: ", propagationVarY
386  return (bandsGraph[-1:], propagationVarY)
xrange
void xrange(TH1 *h, bool symmetric)
Definition: computils.cxx:515
BandUtils.getVariationGraph
def getVariationGraph(graph)
Definition: BandUtils.py:183
BandUtils.getLargestVariationGraph
def getLargestVariationGraph(graph1, graph2)
Make a plot with the largest data point out of two ###########.
Definition: BandUtils.py:113
BandUtils.makeSingleBandPlot
def makeSingleBandPlot(middleGraph, errorGraphs)
Make band function ###########.
Definition: BandUtils.py:212
BandUtils.makeBand
def makeBand(middleGraph, lowerGraph, upperGraph, noorder)
Definition: BandUtils.py:13
BandUtils.makeSingleBandPlot_error
def makeSingleBandPlot_error(middleGraph, varGraphs)
Definition: BandUtils.py:283
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
Generated on Thu Nov 7 2024 21:10:38 for ATLAS Offline Software by   doxygen 1.8.18