MuonAlignmentAlgorithms/scripts/plotscripts.py

0001 from __future__ import print_function
0002 from builtins import range
0003 import ROOT, array, os, re, random
0004 from math import *
0005 import time
0006 import pickle
0007
0008 # python 2.6 has json modue; <2.6 could use simplejson
0009 try:
0010   import json
0011 except ImportError:
0012   import simplejson as json
0013
0014 # sign conventions and some dimensions
0015 from signConventions import *
0016
0017 # common muon types structures
0018 from mutypes import *
0019
0020 CPP_LOADED = False
0021
0022 # containers for test results for map plots
0023 MAP_RESULTS_SAWTOOTH = {}
0024 MAP_RESULTS_FITSIN = {}
0025 MAP_RESULTS_BINS = {}
0026
0027 # general container for all test results
0028 TEST_RESULTS = {}
0029
0030 #############################################################
0031 # Convenience functions
0032
0033 def wheelm2only(dt, wheel, station, sector): return dt == "DT" and wheel == -2
0034 def wheelm1only(dt, wheel, station, sector): return dt == "DT" and wheel == -1
0035 def wheel0only(dt, wheel, station, sector): return dt == "DT" and wheel == 0
0036 def wheelp1only(dt, wheel, station, sector): return dt == "DT" and wheel == 1
0037 def wheelp2only(dt, wheel, station, sector): return dt == "DT" and wheel == 2
0038
0039 def wheelLetter(wheel):
0040   if   wheel == -2: return "A"
0041   elif wheel == -1: return "B"
0042   elif wheel ==  0: return "C"
0043   elif wheel == +1: return "D"
0044   elif wheel == +2: return "E"
0045   else: raise Exception
0046
0047 def wheelNumber(wheell):
0048   if   wheell == "A": return -2
0049   elif wheell == "B": return -1
0050   elif wheell == "C": return 0
0051   elif wheell == "D": return 1
0052   elif wheell == "E": return 2
0053   else: raise Exception
0054
0055 def mean(xlist):
0056   s, n = 0., 0.
0057   for x in xlist:
0058     s += x
0059     n += 1.
0060   return s/n
0061
0062 def rms(xlist):
0063   s2, n = 0., 0.
0064   for x in xlist:
0065     s2 += x**2
0066     n += 1.
0067   return sqrt(s2/n)
0068
0069 def stdev(xlist):
0070   s, s2, n = 0., 0., 0.
0071   for x in xlist:
0072     s += x
0073     s2 += x**2
0074     n += 1.
0075   return sqrt(s2/n - (s/n)**2)
0076
0077 def wmean(xlist):
0078   s, w = 0., 0.
0079   for x, e in xlist:
0080     if e > 0.:
0081       wi = 1./e**2
0082       s += x*wi
0083       w += wi
0084   return s/w, sqrt(1./w)
0085
0086 #############################################################
0087
0088 tdrStyle = None
0089 def setTDRStyle():
0090   global tdrStyle
0091   tdrStyle = ROOT.TStyle("tdrStyle","Style for P-TDR")
0092 # For the canvas:
0093   tdrStyle.SetCanvasBorderMode(0)
0094   tdrStyle.SetCanvasColor(ROOT.kWhite)
0095   tdrStyle.SetCanvasDefH(600) #Height of canvas
0096   tdrStyle.SetCanvasDefW(600) #Width of canvas
0097   tdrStyle.SetCanvasDefX(0)   #POsition on screen
0098   tdrStyle.SetCanvasDefY(0)
0099
0100 # For the Pad:
0101   tdrStyle.SetPadBorderMode(0)
0102   # tdrStyle.SetPadBorderSize(Width_t size = 1)
0103   tdrStyle.SetPadColor(ROOT.kWhite)
0104   tdrStyle.SetPadGridX(False)
0105   tdrStyle.SetPadGridY(False)
0106   tdrStyle.SetGridColor(0)
0107   tdrStyle.SetGridStyle(3)
0108   tdrStyle.SetGridWidth(1)
0109
0110 # For the frame:
0111   tdrStyle.SetFrameBorderMode(0)
0112   tdrStyle.SetFrameBorderSize(1)
0113   tdrStyle.SetFrameFillColor(0)
0114   tdrStyle.SetFrameFillStyle(0)
0115   tdrStyle.SetFrameLineColor(1)
0116   tdrStyle.SetFrameLineStyle(1)
0117   tdrStyle.SetFrameLineWidth(1)
0118
0119 # For the histo:
0120   # tdrStyle.SetHistFillColor(1)
0121   # tdrStyle.SetHistFillStyle(0)
0122   tdrStyle.SetHistLineColor(1)
0123   tdrStyle.SetHistLineStyle(0)
0124   tdrStyle.SetHistLineWidth(1)
0125   # tdrStyle.SetLegoInnerR(Float_t rad = 0.5)
0126   # tdrStyle.SetNumberContours(Int_t number = 20)
0127
0128   tdrStyle.SetEndErrorSize(2)
0129 #  tdrStyle.SetErrorMarker(20)
0130   tdrStyle.SetErrorX(0.)
0131
0132   tdrStyle.SetMarkerStyle(20)
0133
0134 #For the fit/function:
0135   tdrStyle.SetOptFit(1)
0136   tdrStyle.SetFitFormat("5.4g")
0137   tdrStyle.SetFuncColor(2)
0138   tdrStyle.SetFuncStyle(1)
0139   tdrStyle.SetFuncWidth(1)
0140
0141 #For the date:
0142   tdrStyle.SetOptDate(0)
0143   # tdrStyle.SetDateX(Float_t x = 0.01)
0144   # tdrStyle.SetDateY(Float_t y = 0.01)
0145
0146 # For the statistics box:
0147   tdrStyle.SetOptFile(0)
0148   tdrStyle.SetOptStat(0) # To display the mean and RMS:   SetOptStat("mr")
0149   tdrStyle.SetStatColor(ROOT.kWhite)
0150   tdrStyle.SetStatFont(42)
0151   tdrStyle.SetStatFontSize(0.025)
0152   tdrStyle.SetStatTextColor(1)
0153   tdrStyle.SetStatFormat("6.4g")
0154   tdrStyle.SetStatBorderSize(1)
0155   tdrStyle.SetStatH(0.1)
0156   tdrStyle.SetStatW(0.15)
0157   # tdrStyle.SetStatStyle(Style_t style = 1001)
0158   # tdrStyle.SetStatX(Float_t x = 0)
0159   # tdrStyle.SetStatY(Float_t y = 0)
0160
0161 # Margins:
0162   tdrStyle.SetPadTopMargin(0.05)
0163   tdrStyle.SetPadBottomMargin(0.13)
0164   tdrStyle.SetPadLeftMargin(0.13)
0165   tdrStyle.SetPadRightMargin(0.05)
0166
0167 # For the Global title:
0168   tdrStyle.SetOptTitle(0)
0169   tdrStyle.SetTitleFont(42)
0170   tdrStyle.SetTitleColor(1)
0171   tdrStyle.SetTitleTextColor(1)
0172   tdrStyle.SetTitleFillColor(10)
0173   tdrStyle.SetTitleFontSize(0.05)
0174   # tdrStyle.SetTitleH(0) # Set the height of the title box
0175   # tdrStyle.SetTitleW(0) # Set the width of the title box
0176   # tdrStyle.SetTitleX(0) # Set the position of the title box
0177   # tdrStyle.SetTitleY(0.985) # Set the position of the title box
0178   # tdrStyle.SetTitleStyle(Style_t style = 1001)
0179   # tdrStyle.SetTitleBorderSize(2)
0180
0181 # For the axis titles:
0182   tdrStyle.SetTitleColor(1, "XYZ")
0183   tdrStyle.SetTitleFont(42, "XYZ")
0184   tdrStyle.SetTitleSize(0.06, "XYZ")
0185   # tdrStyle.SetTitleXSize(Float_t size = 0.02) # Another way to set the size?
0186   # tdrStyle.SetTitleYSize(Float_t size = 0.02)
0187   tdrStyle.SetTitleXOffset(0.9)
0188   tdrStyle.SetTitleYOffset(1.05)
0189   # tdrStyle.SetTitleOffset(1.1, "Y") # Another way to set the Offset
0190
0191 # For the axis labels:
0192   tdrStyle.SetLabelColor(1, "XYZ")
0193   tdrStyle.SetLabelFont(42, "XYZ")
0194   tdrStyle.SetLabelOffset(0.007, "XYZ")
0195   tdrStyle.SetLabelSize(0.05, "XYZ")
0196
0197 # For the axis:
0198   tdrStyle.SetAxisColor(1, "XYZ")
0199   tdrStyle.SetStripDecimals(True)
0200   tdrStyle.SetTickLength(0.03, "XYZ")
0201   tdrStyle.SetNdivisions(510, "XYZ")
0202   tdrStyle.SetPadTickX(1)  # To get tick marks on the opposite side of the frame
0203   tdrStyle.SetPadTickY(1)
0204
0205 # Change for log plots:
0206   tdrStyle.SetOptLogx(0)
0207   tdrStyle.SetOptLogy(0)
0208   tdrStyle.SetOptLogz(0)
0209
0210 # Postscript options:
0211   tdrStyle.SetPaperSize(20.,20.)
0212   # tdrStyle.SetLineScalePS(Float_t scale = 3)
0213   # tdrStyle.SetLineStyleString(Int_t i, const char* text)
0214   # tdrStyle.SetHeaderPS(const char* header)
0215   # tdrStyle.SetTitlePS(const char* pstitle)
0216
0217   # tdrStyle.SetBarOffset(Float_t baroff = 0.5)
0218   # tdrStyle.SetBarWidth(Float_t barwidth = 0.5)
0219   # tdrStyle.SetPaintTextFormat(const char* format = "g")
0220   # tdrStyle.SetPalette(Int_t ncolors = 0, Int_t* colors = 0)
0221   # tdrStyle.SetTimeOffset(Double_t toffset)
0222   # tdrStyle.SetHistMinimumZero(True)
0223
0224   tdrStyle.cd()
0225
0226 setTDRStyle()
0227
0228 def set_palette(name=None, ncontours=999):
0229     """Set a color palette from a given RGB list
0230     stops, red, green and blue should all be lists of the same length
0231     see set_decent_colors for an example"""
0232
0233     if name == "halfgray":
0234         stops = [0.00, 0.34, 0.61, 0.84, 1.00]
0235         red   = map(lambda x: 1. - (1.-x)/2., [1.00, 0.84, 0.61, 0.34, 0.00])
0236         green = map(lambda x: 1. - (1.-x)/2., [1.00, 0.84, 0.61, 0.34, 0.00])
0237         blue  = map(lambda x: 1. - (1.-x)/2., [1.00, 0.84, 0.61, 0.34, 0.00])
0238     elif name == "gray":
0239         stops = [0.00, 0.34, 0.61, 0.84, 1.00]
0240         red   = [1.00, 0.84, 0.61, 0.34, 0.00]
0241         green = [1.00, 0.84, 0.61, 0.34, 0.00]
0242         blue  = [1.00, 0.84, 0.61, 0.34, 0.00]
0243     elif name == "blues":
0244         stops = [0.00, 0.34, 0.61, 0.84, 1.00]
0245         red   = [1.00, 0.84, 0.61, 0.34, 0.00]
0246         green = [1.00, 0.84, 0.61, 0.34, 0.00]
0247         blue  = [1.00, 1.00, 1.00, 1.00, 1.00]
0248     elif name == "reds":
0249         stops = [0.00, 0.34, 0.61, 0.84, 1.00]
0250         red   = [1.00, 1.00, 1.00, 1.00, 1.00]
0251         green = [1.00, 0.84, 0.61, 0.34, 0.00]
0252         blue  = [1.00, 0.84, 0.61, 0.34, 0.00]
0253     elif name == "antigray":
0254         stops = [0.00, 0.34, 0.61, 0.84, 1.00]
0255         red   = [1.00, 0.84, 0.61, 0.34, 0.00]
0256         green = [1.00, 0.84, 0.61, 0.34, 0.00]
0257         blue  = [1.00, 0.84, 0.61, 0.34, 0.00]
0258         red.reverse()
0259         green.reverse()
0260         blue.reverse()
0261     elif name == "fire":
0262         stops = [0.00, 0.20, 0.80, 1.00]
0263         red   = [1.00, 1.00, 1.00, 0.50]
0264         green = [1.00, 1.00, 0.00, 0.00]
0265         blue  = [0.20, 0.00, 0.00, 0.00]
0266     elif name == "antifire":
0267         stops = [0.00, 0.20, 0.80, 1.00]
0268         red   = [0.50, 1.00, 1.00, 1.00]
0269         green = [0.00, 0.00, 1.00, 1.00]
0270         blue  = [0.00, 0.00, 0.00, 0.20]
0271     else:
0272         # default palette, looks cool
0273         stops = [0.00, 0.34, 0.61, 0.84, 1.00]
0274         red   = [0.00, 0.00, 0.87, 1.00, 0.51]
0275         green = [0.00, 0.81, 1.00, 0.20, 0.00]
0276         blue  = [0.51, 1.00, 0.12, 0.00, 0.00]
0277
0278     s = array.array('d', stops)
0279     r = array.array('d', red)
0280     g = array.array('d', green)
0281     b = array.array('d', blue)
0282
0283     npoints = len(s)
0284     ROOT.TColor.CreateGradientColorTable(npoints, s, r, g, b, ncontours)
0285     ROOT.gStyle.SetNumberContours(ncontours)
0286
0287 set_palette()
0288
0289 ######################################################################################################
0290 ## sector phi edges in: me11 me12 me13 me14 me21 me22 me31 me32 me41 me42 mb1 mb2 mb3 mb4
0291 ## index:               0    1    2    3    4    5    6    7    8    9    10  11  12  13
0292
0293 #phiedgesCSC36 = [pi/180.*(-175. + 10.*i) for i in range(36)]
0294 #phiedgesCSC18 = [pi/180.*(-175. + 20.*i) for i in range(18)]
0295 phiedgesCSC36 = [pi/180.*(-5. + 10.*i) for i in range(36)]
0296 phiedgesCSC18 = [pi/180.*(-5. + 20.*i) for i in range(18)]
0297 phiedges = [
0298    phiedgesCSC36,
0299    phiedgesCSC36,
0300    phiedgesCSC36,
0301    phiedgesCSC36,
0302    phiedgesCSC18,
0303    phiedgesCSC36,
0304    phiedgesCSC18,
0305    phiedgesCSC36,
0306    phiedgesCSC18,
0307    phiedgesCSC36,
0308    [0.35228048120123945, 0.87587781482541827, 1.3994776462193192, 1.923076807996136, 2.4466741416203148, 2.970273973014216,
0309     -2.7893121723885534, -2.2657148387643748, -1.7421150073704739, -1.2185158455936571, -0.69491851196947851, -0.17131868057557731],
0310    [0.22000706229660855, 0.74360690430428489, 1.267204926935573, 1.7908033890915052, 2.3144032310991816, 2.8380012537304697,
0311     -2.9215855912931841, -2.3979857492855081, -1.8743877266542202, -1.3507892644982882, -0.82718942249061178, -0.30359139985932365],
0312    [0.29751957124275596, 0.82111826253905784, 1.3447162969496083, 1.8683158980376524, 2.3919145893339548, 2.915512623744505,
0313     -2.844073082347037, -2.3204743910507353, -1.7968763566401849, -1.2732767555521407, -0.74967806425583894, -0.22608002984528835],
0314    [3.0136655290752188, -2.7530905195097337, -2.2922883025568734, -1.9222915077192773, -1.5707963267948966, -1.2193011458705159,
0315     -0.84930435103291968, -0.38850213408005951, 0.127927124514574, 0.65152597487624719, 1.1322596819239259, 1.5707963267948966,
0316     2.0093329716658674, 2.4900666787135459]]
0317
0318 def phiedges2c():
0319   lines = []
0320   for ed in phiedges[:]:
0321     ed.sort()
0322     #print ed
0323     ed.extend([999 for n in range(0,37-len(ed))])
0324     lines.append('{' + ', '.join(map(str, ed)) + '}')
0325   #print lines
0326   res = ', '.join(lines)
0327   ff = open("phiedges_export.h",mode="w")
0328   print('double phiedges[14][37] = {' + res + '};', file=ff)
0329   ff.close()
0330
0331 class SawTeethFunction:
0332   def __init__(self, name):
0333     self.name = name
0334     self.edges = (phiedges[stationIndex(name)])[:]
0335     self.ed = sorted(self.edges)
0336     # add some padding to the end
0337     self.ed.append(pi+1.)
0338     self.n = len(self.edges)
0339   def __call__(self, xx, par):
0340     # wrap x in the most negative phi sector into positive phi
0341     x = xx[0]
0342     if x < self.ed[0]: x += 2*pi
0343     # locate sector
0344     for i in range(0,self.n):
0345       if x <= self.ed[i]: continue
0346       if x > self.ed[i+1]: continue
0347       return par[i*2] + par[i*2+1]*(x - self.ed[i])
0348     return 0
0349   def pp(self):
0350     print(self.name, self.n)
0351     print(self.edges)
0352     print(self.ed)
0353
0354
0355 def stationIndex(name):
0356   if ("MB" in name or "ME" in name):
0357     # assume the name is ID
0358     pa = idToPostalAddress(name)
0359     if pa is None: return None
0360     if pa[0]=="CSC":
0361       if pa[2]==1 and pa[3]==1: return 0
0362       if pa[2]==1 and pa[3]==2: return 1
0363       if pa[2]==1 and pa[3]==3: return 2
0364       if pa[2]==1 and pa[3]==4: return 3
0365       if pa[2]==2 and pa[3]==1: return 4
0366       if pa[2]==2 and pa[3]==2: return 5
0367       if pa[2]==3 and pa[3]==1: return 6
0368       if pa[2]==3 and pa[3]==2: return 7
0369       if pa[2]==4 and pa[3]==1: return 8
0370       if pa[2]==4 and pa[3]==2: return 9
0371     if pa[0]=="DT":
0372       if pa[2]==1: return 10
0373       if pa[2]==2: return 11
0374       if pa[2]==3: return 12
0375       if pa[2]==4: return 13
0376   else:
0377     if ("mem11" in name or "mep11" in name): return 0
0378     if ("mem12" in name or "mep12" in name): return 1
0379     if ("mem13" in name or "mep13" in name): return 2
0380     if ("mem14" in name or "mep14" in name): return 3
0381     if ("mem21" in name or "mep21" in name): return 4
0382     if ("mem22" in name or "mep22" in name): return 5
0383     if ("mem31" in name or "mep31" in name): return 6
0384     if ("mem32" in name or "mep32" in name): return 7
0385     if ("mem41" in name or "mep41" in name): return 8
0386     if ("mem42" in name or "mep42" in name): return 9
0387     if ("st1" in name): return 10
0388     if ("st2" in name): return 11
0389     if ("st3" in name): return 12
0390     if ("st4" in name): return 13
0391
0392
0393
0394 def philines(name, window, abscissa):
0395     global philine_tlines, philine_labels
0396     philine_tlines = []
0397     edges = phiedges[stationIndex(name)]
0398     #print name, len(edges)
0399     for phi in edges:
0400         if abscissa is None or abscissa[0] < phi < abscissa[1]:
0401             philine_tlines.append(ROOT.TLine(phi, -window, phi, window))
0402             philine_tlines[-1].SetLineStyle(2)
0403             philine_tlines[-1].Draw()
0404     if "st" in name: # DT labels
0405         philine_labels = []
0406         edges = sorted(edges[:])
0407         if "st4" in name:
0408             labels = [" 7", " 8", " 9", "14", "10", "11", "12", " 1", " 2", " 3", "13", " 4", " 5", " 6"]
0409         else:
0410             labels = [" 8", " 9", "10", "11", "12", " 1", " 2", " 3", " 4", " 5", " 6"]
0411             edges = edges[1:]
0412         for phi, label in zip(edges, labels):
0413             littlebit = 0.
0414             if label in (" 7", " 9", "14", "10", "11"): littlebit = 0.05
0415             philine_labels.append(ROOT.TText(phi-0.35+littlebit, -0.9*window, label))
0416             philine_labels[-1].Draw()
0417         philine_labels.append(ROOT.TText(-2.9, -0.75*window, "Sector:"))
0418         philine_labels[-1].Draw()
0419     if "CSC" in name: # DT labels
0420         philine_labels = []
0421         edges = sorted(edges[:])
0422         labels = [" 1", " 2", " 3", " 4", " 5", " 6", " 7", " 8", " 9", "10", "11", "12", "13", "14", "15", "16", "17", "18",
0423                   "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36"]
0424         #else:
0425         #    labels = [" 8", " 9", "10", "11", "12", " 1", " 2", " 3", " 4", " 5", " 6"]
0426         #    edges = edges[1:]
0427         for phi, label in zip(edges, labels):
0428             littlebit = 0.
0429             #if label in (" 7", " 9", "14", "10", "11"): littlebit = 0.05
0430             philine_labels.append(ROOT.TText(phi+littlebit, -0.9*window, label))
0431             philine_labels[-1].SetTextFont(42)
0432             philine_labels[-1].SetTextSize(0.028)
0433             philine_labels[-1].Draw()
0434         philine_labels.append(ROOT.TText(0, -0.78*window, "Chamber:"))
0435         philine_labels[-1].SetTextSize(0.035)
0436         philine_labels[-1].Draw()
0437
0438 def zlines(window, abscissa):
0439     global zline_tlines
0440     zline_tlines = []
0441     for z in -401.625, -133.875, 133.875, 401.625:
0442         if abscissa is None or abscissa[0] < z < abscissa[1]:
0443             zline_tlines.append(ROOT.TLine(z, -window, z, window))
0444             zline_tlines[-1].SetLineStyle(2)
0445             zline_tlines[-1].Draw()
0446     zline_labels = []
0447     zline_labels.append(ROOT.TText(-550, -0.9*window, "-2"))
0448     zline_labels.append(ROOT.TText(-300, -0.9*window, "-1"))
0449     zline_labels.append(ROOT.TText(-10, -0.9*window, "0"))
0450     zline_labels.append(ROOT.TText(250, -0.9*window, "+1"))
0451     zline_labels.append(ROOT.TText(500, -0.9*window, "+2"))
0452     for z in zline_labels: z.Draw()
0453     zline_labels.append(ROOT.TText(-600, -0.75*window, "Wheel:")); zline_labels[-1].Draw()
0454
0455 def rlines(disk, window, abscissa):
0456     global rline_tlines
0457     rline_tlines = []
0458     if disk == 1: rl = [150., 270., 480.]
0459     else: rl = [350.]
0460     for r in rl:
0461         if abscissa is None or abscissa[0] < r < abscissa[1]:
0462             rline_tlines.append(ROOT.TLine(r, -window, r, window))
0463             rline_tlines[-1].SetLineStyle(2)
0464             rline_tlines[-1].Draw()
0465
0466 ######################################################################################################
0467
0468 def getReportByPostalAddress(postal_address, report):
0469   for r in report:
0470     if postal_address == r.postal_address:
0471       return r
0472   return None
0473
0474
0475 ######################################################################################################
0476
0477 def DBMC(database, reports, window=10., windows=None, selection=None, phi=False,
0478          color=ROOT.kBlue-8, style=1, bins=50, normalized=False, getvalues=False, name="", canvas=None, reportdiff=False, inlog=True):
0479     return DBdiff(database, None, reports, None, window, windows, selection, phi, color, style, bins, normalized, getvalues,
0480                   name, canvas, reportdiff, inlog)
0481
0482
0483 def DBdiff(database1, database2, reports1, reports2,
0484            window=10., windows=None, selection=None, phi=False, color=ROOT.kBlue-8,
0485            style=1, bins=50, normalized=False, getvalues=False, name="tmp", canvas=None, reportdiff=False, inlog=False ):
0486
0487     tdrStyle.SetOptStat("emrou")
0488     tdrStyle.SetStatW(0.40)
0489
0490     wnd = [window]*6
0491     if windows is not None:
0492       i=0
0493       for w in windows:
0494         wnd[i] = windows[i]
0495         i+=1
0496
0497     global hx, hy, hz, hphix, hphiy, hphiz
0498
0499     if phi:
0500         hx = ROOT.TH1F("%s_phi" % name, "", bins, -wnd[0], wnd[0])
0501     else:
0502         hx = ROOT.TH1F("%s_x" % name, "", bins, -wnd[0], wnd[0])
0503     hy = ROOT.TH1F("%s_y" % name, "", bins, -wnd[1], wnd[1])
0504     hz = ROOT.TH1F("%s_z" % name, "", bins, -wnd[2], wnd[2])
0505     hphix = ROOT.TH1F("%s_phix" % name, "", bins, -wnd[3], wnd[3])
0506     hphiy = ROOT.TH1F("%s_phiy" % name, "", bins, -wnd[4], wnd[4])
0507     hphiz = ROOT.TH1F("%s_phiz" % name, "", bins, -wnd[5], wnd[5])
0508
0509     for r1 in reports1:
0510         if selection is None or (selection.__code__.co_argcount == len(r1.postal_address) and selection(*r1.postal_address)):
0511             if reports2 is None:
0512                 r2 = Report(r1.chamberId, r1.postal_address, r1.name)
0513                 r2.add_parameters(ValErr(0., 0., 0.), ValErr(0., 0., 0.), ValErr(0., 0., 0.),
0514                     ValErr(0., 0., 0.), ValErr(0., 0., 0.), ValErr(0., 0., 0.), 0., 0., 0., 0.)
0515             else:
0516                 r2 = getReportByPostalAddress(r1.postal_address, reports2)
0517                 if r2 is None: continue
0518
0519             found = False
0520             if r1.postal_address[0] == "DT":
0521                 if r1.postal_address[1:] in database1.dt:
0522                     found = True
0523                     db1 = database1.dt[r1.postal_address[1:]]
0524                     if database2 is None:
0525                         db2 = DTAlignable()
0526                         db2.x = db2.y = db2.z = db2.phix = db2.phiy = db2.phiz = 0.
0527                         db2.xx = db2.xy = db2.xz = db2.yx = db2.yy = db2.yz = db2.zx = db2.zy = db2.zz = 0.
0528                     else:
0529                         db2 = database2.dt[r1.postal_address[1:]]
0530
0531             else:
0532                 # skip ME1/a
0533                 if r1.postal_address[2]==1 and r1.postal_address[3]==4: continue
0534                 if r1.postal_address[1:] in database1.csc:
0535                     found = True
0536                     db1 = database1.csc[r1.postal_address[1:]]
0537                     if database2 is None:
0538                         db2 = CSCAlignable()
0539                         db2.x = db2.y = db2.z = db2.phix = db2.phiy = db2.phiz = 0.
0540                         db2.xx = db2.xy = db2.xz = db2.yx = db2.yy = db2.yz = db2.zx = db2.zy = db2.zz = 0.
0541                     else:
0542                         db2 = database2.csc[r1.postal_address[1:]]
0543
0544             if found and r1.status == "PASS" and r2.status == "PASS":
0545                 if r1.deltax is not None and r2.deltax is not None and r1.deltax.error is not None and \
0546                    r2.deltax.error is not None and (r1.deltax.error**2 + r2.deltax.error**2) > 0.:
0547                     delta = db1.x - db2.x
0548                     if reportdiff: delta -= r1.deltax.value
0549                     if normalized:
0550                         fill = delta/sqrt(r1.deltax.error**2 + r2.deltax.error**2) * signConventions[r1.postal_address][0]
0551                     else:
0552                         if phi:
0553                             fill = delta/signConventions[r1.postal_address][3] * 1000. * signConventions[r1.postal_address][0]
0554                         else:
0555                             fill = delta * 10. * signConventions[r1.postal_address][0]
0556                     hx.Fill(fill)
0557                     if getvalues not in (False, None):
0558                         getvalues["x"].append((fill, 10. * sqrt(r1.deltax.error**2 + r2.deltax.error**2)))
0559
0560                 if r1.deltay is not None and r2.deltay is not None and r1.deltay.error is not None and \
0561                    r2.deltay.error is not None and (r1.deltay.error**2 + r2.deltay.error**2) > 0.:
0562                     delta = db1.y - db2.y
0563                     if reportdiff: delta -= r1.deltay.value
0564                     if normalized:
0565                         fill = delta/sqrt(r1.deltay.error**2 + r2.deltay.error**2) * signConventions[r1.postal_address][1]
0566                     else:
0567                         fill = delta * 10. * signConventions[r1.postal_address][1]
0568                     hy.Fill(fill)
0569                     if getvalues not in (False, None):
0570                         getvalues["y"].append((fill, 10. * sqrt(r1.deltay.error**2 + r2.deltay.error**2)))
0571
0572                 if r1.deltaz is not None and r2.deltaz is not None and r1.deltaz.error is not None and \
0573                    r2.deltaz.error is not None and (r1.deltaz.error**2 + r2.deltaz.error**2) > 0.:
0574                     delta = db1.z - db2.z
0575                     if reportdiff: delta -= r1.deltaz.value
0576                     if normalized:
0577                         fill = delta/sqrt(r1.deltaz.error**2 + r2.deltaz.error**2) * signConventions[r1.postal_address][2]
0578                     else:
0579                         fill = delta * 10. * signConventions[r1.postal_address][2]
0580                     hz.Fill(fill)
0581                     if getvalues not in (False, None):
0582                         getvalues["z"].append((fill, 10. * sqrt(r1.deltaz.error**2 + r2.deltaz.error**2)))
0583
0584                 if r1.deltaphix is not None and r2.deltaphix is not None and r1.deltaphix.error is not None and \
0585                    r2.deltaphix.error is not None and (r1.deltaphix.error**2 + r2.deltaphix.error**2) > 0.:
0586                     delta = db1.phix - db2.phix
0587                     if reportdiff: delta -= r1.deltaphix.value
0588                     if normalized:
0589                         fill = delta/sqrt(r1.deltaphix.error**2 + r2.deltaphix.error**2)
0590                     else:
0591                         fill = delta * 1000.
0592                     hphix.Fill(fill)
0593                     if getvalues not in (False, None):
0594                         getvalues["phix"].append((fill, 10. * sqrt(r1.deltaphix.error**2 + r2.deltaphix.error**2)))
0595
0596                 if r1.deltaphiy is not None and r2.deltaphiy is not None and r1.deltaphiy.error is not None and \
0597                    r2.deltaphiy.error is not None and (r1.deltaphiy.error**2 + r2.deltaphiy.error**2) > 0.:
0598                     delta = db1.phiy - db2.phiy
0599                     if reportdiff:
0600                       delta -= r1.deltaphiy.value
0601                       if abs(delta)>0.02/1000: print(r1.postal_address, 1000*delta, "=", 1000*db1.phiy - 1000*db2.phiy, "-", 1000*r1.deltaphiy.value, "... ",1000*db1.phiy , 1000*db2.phiy)
0602                     if normalized:
0603                         fill = delta/sqrt(r1.deltaphiy.error**2 + r2.deltaphiy.error**2)
0604                     else:
0605                         fill = delta * 1000.
0606                     hphiy.Fill(fill)
0607                     if getvalues not in (False, None):
0608                         getvalues["phiy"].append((fill, 10. * sqrt(r1.deltaphiy.error**2 + r2.deltaphiy.error**2)))
0609
0610                 if r1.deltaphiz is not None and r2.deltaphiz is not None and r1.deltaphiz.error is not None and \
0611                    r2.deltaphiz.error is not None and (r1.deltaphiz.error**2 + r2.deltaphiz.error**2) > 0.:
0612                     delta = db1.phiz - db2.phiz
0613                     if reportdiff: delta -= r1.deltaphiz.value
0614                     if normalized:
0615                         fill = delta/sqrt(r1.deltaphiz.error**2 + r2.deltaphiz.error**2)
0616                     else:
0617                         fill = delta * 1000.
0618                     hphiz.Fill(fill)
0619                     if getvalues not in (False, None):
0620                         getvalues["phiz"].append((fill, 10. * sqrt(r1.deltaphiz.error**2 + r2.deltaphiz.error**2)))
0621
0622     if not normalized:
0623         if phi:
0624             hx.SetXTitle("#delta_{#phi} position (mrad)")
0625         else:
0626             hx.SetXTitle("#delta_{x'} (mm)")
0627         hy.SetXTitle("#delta_{y'} (mm)")
0628         hz.SetXTitle("#delta_{z'} (mm)")
0629         hphix.SetXTitle("#delta_{#phi_{x}} (mrad)")
0630         hphiy.SetXTitle("#delta_{#phi_{y}} (mrad)")
0631         hphiz.SetXTitle("#delta_{#phi_{z}} (mrad)")
0632         if reportdiff:
0633           if phi:
0634               hx.SetXTitle("#delta_{#phi}(XML) - #delta_{#phi}(report) position (mrad)")
0635           else:
0636               hx.SetXTitle("#delta_{x'}(XML) - #delta_{x'}(report) (mm)")
0637           hy.SetXTitle("#delta_{y'}(XML) - #delta_{y'}(report) (mm)")
0638           hz.SetXTitle("#delta_{z'}(XML) - #delta_{z'}(report) (mm)")
0639           hphix.SetXTitle("#delta_{#phi_{x}}(XML) - #delta_{#phi_{x}}(report) (mrad)")
0640           hphiy.SetXTitle("#delta_{#phi_{y}}(XML) - #delta_{#phi_{y}}(report) (mrad)")
0641           hphiz.SetXTitle("#delta_{#phi_{z}}(XML) - #delta_{#phi_{z}}(report) (mrad)")
0642     else:
0643         if phi:
0644             hx.SetXTitle("#delta_{#phi}/#sigma_{#phi} position")
0645         else:
0646             hx.SetXTitle("#delta_{x'}/#sigma_{x'}")
0647         hy.SetXTitle("#delta_{y'}/#sigma_{y'}")
0648         hz.SetXTitle("#delta_{z'}/#sigma_{z'}")
0649         hphix.SetXTitle("#delta_{#phi_{x}}/#sigma_{#phi_{x}}")
0650         hphiy.SetXTitle("#delta_{#phi_{y}}/#sigma_{#phi_{y}}")
0651         hphiz.SetXTitle("#delta_{#phi_{z}}/#sigma_{#phi_{z}}")
0652
0653     for h in hx, hy, hz, hphix, hphiy, hphiz:
0654         h.GetXaxis().CenterTitle()
0655         h.GetYaxis().CenterTitle()
0656         h.SetFillColor(color)
0657         h.SetLineStyle(style)
0658
0659     if canvas is not None: c = canvas
0660     else: c = c1
0661
0662     if normalized:
0663         fx = ROOT.TF1("fx", "%g * exp(-x**2/2.)/sqrt(2.*3.1415926)" % (hx.GetEntries()*2.*window/bins), -window, window)
0664         fy = ROOT.TF1("fy", "%g * exp(-x**2/2.)/sqrt(2.*3.1415926)" % (hy.GetEntries()*2.*window/bins), -window, window)
0665         fz = ROOT.TF1("fz", "%g * exp(-x**2/2.)/sqrt(2.*3.1415926)" % (hz.GetEntries()*2.*window/bins), -window, window)
0666         fphix = ROOT.TF1("fphix", "%g * exp(-x**2/2.)/sqrt(2.*3.1415926)" % (hphix.GetEntries()*2.*window/bins), -window, window)
0667         fphiy = ROOT.TF1("fphiy", "%g * exp(-x**2/2.)/sqrt(2.*3.1415926)" % (hphiy.GetEntries()*2.*window/bins), -window, window)
0668         fphiz = ROOT.TF1("fphiz", "%g * exp(-x**2/2.)/sqrt(2.*3.1415926)" % (hphiz.GetEntries()*2.*window/bins), -window, window)
0669         for f in fx, fy, fz, fphix, fphiy, fphiz:
0670             f.SetLineWidth(2)
0671             f.SetLineColor(ROOT.kBlue)
0672         for h, f in (hx, fx), (hy, fy), (hz, fz), (hphix, fphix), (hphiy, fphiy), (hphiz, fphiz):
0673             h.SetAxisRange(0, 1.1*max(h.GetMaximum(), f.GetMaximum()), "Y")
0674
0675         c.Clear()
0676         c.Divide(3, 2)
0677         c.GetPad(1).cd(); hx.Draw(); fx.Draw("same")
0678         c.GetPad(2).cd(); hy.Draw(); fy.Draw("same")
0679         c.GetPad(3).cd(); hz.Draw(); fz.Draw("same")
0680         c.GetPad(4).cd(); hphix.Draw(); fphix.Draw("same")
0681         c.GetPad(5).cd(); hphiy.Draw(); fphiy.Draw("same")
0682         c.GetPad(6).cd(); hphiz.Draw(); fphiz.Draw("same")
0683         return hx, hy, hz, hphix, hphiy, hphiz, fx, fy, fz, fphix, fphiy, fphiz
0684     else:
0685         nvar = 6
0686
0687         c.Clear()
0688         if nvar == 4: c.Divide(2, 2)
0689         if nvar == 6: c.Divide(3, 2)
0690         c.GetPad(1).cd(); hx.Draw()
0691         c.GetPad(2).cd(); hy.Draw()
0692         if nvar == 4:
0693           c.GetPad(3).cd(); hphiy.Draw()
0694           c.GetPad(4).cd(); hphiz.Draw()
0695         if nvar == 6:
0696           c.GetPad(3).cd(); hz.Draw()
0697           c.GetPad(4).cd(); hphix.Draw()
0698           c.GetPad(5).cd(); hphiy.Draw()
0699           c.GetPad(6).cd(); hphiz.Draw()
0700
0701         if inlog:
0702           if hx.GetEntries()>0:    c.GetPad(1).SetLogy(1)
0703           if hy.GetEntries()>0:    c.GetPad(2).SetLogy(1)
0704           if nvar == 4:
0705             if hphiy.GetEntries()>0: c.GetPad(3).SetLogy(1)
0706             if hphiz.GetEntries()>0: c.GetPad(4).SetLogy(1)
0707           if nvar == 6:
0708             if hz.GetEntries()>0:    c.GetPad(3).SetLogy(1)
0709             if hphix.GetEntries()>0: c.GetPad(4).SetLogy(1)
0710             if hphiy.GetEntries()>0: c.GetPad(5).SetLogy(1)
0711             if hphiz.GetEntries()>0: c.GetPad(6).SetLogy(1)
0712
0713         return hx, hy, hz, hphix, hphiy, hphiz
0714
0715
0716
0717 def DBMCVersus(quantity, versus, database, reports, window=10., selection=None, color=ROOT.kBlack):
0718     return DBdiffVersus(quantity, versus, database, None, reports, None, window, selection, color)
0719
0720 def DBdiffVersus(quantity, versus, database1, database2, reports1, reports2, windwselection=None, color=ROOT.kBlack):
0721     tdrStyle.SetOptStat("")
0722
0723     domain = []
0724     values = []
0725     errors = []
0726
0727     for r1 in reports1:
0728         if selection is None or (selection.__code__.co_argcount == len(r1.postal_address) and selection(*r1.postal_address)):
0729             if reports2 is None:
0730                 r2 = Report(r1.chamberId, r1.postal_address, r1.name)
0731                 r2.add_parameters(ValErr(0., 0., 0.), ValErr(0., 0., 0.), ValErr(0., 0., 0.),
0732                                   ValErr(0., 0., 0.), ValErr(0., 0., 0.), ValErr(0., 0., 0.), 0., 0., 0.)
0733             else:
0734                 found = False
0735                 for r2 in reports2:
0736                     if r1.postal_address == r2.postal_address:
0737                         found = True
0738                         break
0739                 if not found: continue
0740
0741             found = False
0742             if r1.postal_address[0] == "DT":
0743                 if r1.postal_address[1:] in database1.dt:
0744                     found = True
0745                     db1 = database1.dt[r1.postal_address[1:]]
0746                     if database2 is None:
0747                         db2 = DTAlignable()
0748                         db2.x = db2.y = db2.z = db2.phix = db2.phiy = db2.phiz = 0.
0749                         db2.xx = db2.xy = db2.xz = db2.yx = db2.yy = db2.yz = db2.zx = db2.zy = db2.zz = 0.
0750                     else:
0751                         db2 = database2.dt[r1.postal_address[1:]]
0752             else:
0753                 if r1.postal_address[1:] in database1.csc:
0754                     found = True
0755                     db1 = database1.csc[r1.postal_address[1:]]
0756                     if database2 is None:
0757                         db2 = CSCAlignable()
0758                         db2.x = db2.y = db2.z = db2.phix = db2.phiy = db2.phiz = 0.
0759                         db2.xx = db2.xy = db2.xz = db2.yx = db2.yy = db2.yz = db2.zx = db2.zy = db2.zz = 0.
0760                     else:
0761                         db2 = database2.csc[r1.postal_address[1:]]
0762
0763             if found and r1.status == "PASS" and r2.status == "PASS":
0764                 okay = False
0765
0766                 if quantity == "phi":
0767                     if r1.deltax is not None and r2.deltax is not None and r1.deltax.error is not None and \
0768                        r2.deltax.error is not None and (r1.deltax.error**2 + r2.deltax.error**2) > 0.:
0769                         okay = True
0770                         values.append((db1.x - db2.x)/
0771                                       signConventions[r1.postal_address][3] * 1000. * signConventions[r1.postal_address][0])
0772                         errors.append((r1.deltax.error**2 + r2.deltax.error**2)/
0773                                       signConventions[r1.postal_address][3] * 1000. * signConventions[r1.postal_address][0])
0774
0775                 elif quantity == "x":
0776                     if r1.deltax is not None and r2.deltax is not None and r1.deltax.error is not None and \
0777                        r2.deltax.error is not None and (r1.deltax.error**2 + r2.deltax.error**2) > 0.:
0778                         okay = True
0779                         values.append((db1.x - db2.x) * 10. * signConventions[r1.postal_address][0])
0780                         errors.append((r1.deltax.error**2 + r2.deltax.error**2) * 10. * signConventions[r1.postal_address][0])
0781
0782                 elif quantity == "y":
0783                     if r1.deltay is not None and r2.deltay is not None and r1.deltay.error is not None and \
0784                        r2.deltay.error is not None and (r1.deltay.error**2 + r2.deltay.error**2) > 0.:
0785                         okay = True
0786                         values.append((db1.y - db2.y) * 10. * signConventions[r1.postal_address][1])
0787                         errors.append((r1.deltay.error**2 + r2.deltay.error**2) * 10. * signConventions[r1.postal_address][1])
0788
0789                 elif quantity == "z":
0790                     if r1.deltaz is not None and r2.deltaz is not None and r1.deltaz.error is not None and \
0791                        r2.deltaz.error is not None and (r1.deltaz.error**2 + r2.deltaz.error**2) > 0.:
0792                         okay = True
0793                         values.append((db1.z - db2.z) * 10. * signConventions[r1.postal_address][2])
0794                         errors.append((r1.deltaz.error**2 + r2.deltaz.error**2) * 10. * signConventions[r1.postal_address][2])
0795
0796                 elif quantity == "phix":
0797                     if r1.deltaphix is not None and r2.deltaphix is not None and r1.deltaphix.error is not None and \
0798                        r2.deltaphix.error is not None and (r1.deltaphix.error**2 + r2.deltaphix.error**2) > 0.:
0799                         okay = True
0800                         values.append((db1.phix - db2.phix) * 1000.)
0801                         errors.append((r1.deltaphix.error**2 + r2.deltaphix.error**2) * 1000.)
0802
0803                 elif quantity == "phiy":
0804                     if r1.deltaphiy is not None and r2.deltaphiy is not None and r1.deltaphiy.error is not None and \
0805                        r2.deltaphiy.error is not None and (r1.deltaphiy.error**2 + r2.deltaphiy.error**2) > 0.:
0806                         okay = True
0807                         values.append((db1.phiy - db2.phiy) * 1000.)
0808                         errors.append((r1.deltaphiy.error**2 + r2.deltaphiy.error**2) * 1000.)
0809
0810                 elif quantity == "phiz":
0811                     if r1.deltaphiz is not None and r2.deltaphiz is not None and r1.deltaphiz.error is not None and \
0812                        r2.deltaphiz.error is not None and (r1.deltaphiz.error**2 + r2.deltaphiz.error**2) > 0.:
0813                         okay = True
0814                         values.append((db1.phiz - db2.phiz) * 1000.)
0815                         errors.append((r1.deltaphiz.error**2 + r2.deltaphiz.error**2) * 1000.)
0816
0817                 else: raise Exception
0818
0819                 if okay:
0820                     if versus == "r": domain.append(signConventions[r1.postal_address][3])
0821                     elif versus == "phi": domain.append(signConventions[r1.postal_address][4])
0822                     elif versus == "z": domain.append(signConventions[r1.postal_address][5])
0823                     else: raise Exception
0824
0825     if versus == "r":
0826         bkgndhist = ROOT.TH1F("bkgndhist", "", 100, 0., 800.)
0827         bkgndhist.SetXTitle("R (cm)")
0828     elif versus == "phi":
0829         bkgndhist = ROOT.TH1F("bkgndhist", "", 100, -pi, pi)
0830         bkgndhist.SetXTitle("#phi (rad)")
0831     elif versus == "z":
0832         bkgndhist = ROOT.TH1F("bkgndhist", "", 100, -1100., 1100.)
0833         bkgndhist.SetXTitle("z (cm)")
0834     bkgndhist.GetXaxis().CenterTitle()
0835
0836     bkgndhist.SetAxisRange(-window, window, "Y")
0837     if quantity == "phi": bkgndhist.SetYTitle("#delta_{#phi} position (mrad)")
0838     elif quantity == "x": bkgndhist.SetYTitle("#delta_{x'} (mm)")
0839     elif quantity == "y": bkgndhist.SetYTitle("#delta_{y'} (mm)")
0840     elif quantity == "z": bkgndhist.SetYTitle("#delta_{z'} (mm)")
0841     elif quantity == "phix": bkgndhist.SetYTitle("#delta_{#phi_{x}} (mrad)")
0842     elif quantity == "phiy": bkgndhist.SetYTitle("#delta_{#phi_{y}} (mrad)")
0843     elif quantity == "phiz": bkgndhist.SetYTitle("#delta_{#phi_{z}} (mrad)")
0844     else: raise Exception
0845     bkgndhist.GetYaxis().CenterTitle()
0846
0847     if len(domain) == 0:
0848         tgraph = ROOT.TGraphErrors(0)
0849     else:
0850         tgraph = ROOT.TGraphErrors(len(domain), array.array("d", domain), array.array("d", values),
0851                                                 array.array("d", [0.]*len(domain)), array.array("d", errors))
0852     tgraph.SetMarkerColor(color)
0853     tgraph.SetLineColor(color)
0854
0855     bkgndhist.Draw()
0856     if tgraph.GetN() > 0: tgraph.Draw("p")
0857     return bkgndhist, tgraph, domain, values, errors
0858
0859 ######################################################################################################
0860
0861 def idToPostalAddress(id):
0862   # only len==9 ids can correspond to valid postal address
0863   if len(id)!=9: return None
0864   if id[0:2]=="MB":
0865     #print id
0866     pa = ("DT", int(id[2:4]), int(id[5]), int(id[7:9]))
0867     #print pa
0868     if pa[1]<-2 or pa[1]>2: return None
0869     if pa[2]>4: return None
0870     if pa[3]<1 or pa[3]>14 or (pa[3]==4 and pa[3]>12): return None
0871     return pa
0872   elif id[0:2]=="ME":
0873     if id[2]=="+": ec=1
0874     elif id[2]=="-": ec=2
0875     else: return None
0876     pa = ("CSC", ec, int(id[3]), int(id[5]), int(id[7:9]))
0877     if pa[2]<1 or pa[2]>4: return None
0878     if pa[3]<1 or pa[3]>4 or (pa[2]>1 and pa[3]>2): return None
0879     if pa[4]<1 or pa[4]>36 or (pa[2]>1 and pa[3]==1 and pa[4]>18): return None
0880     return pa
0881   else: return None
0882
0883
0884 def postalAddressToId(postal_address):
0885   if postal_address[0] == "DT":
0886     wheel, station, sector = postal_address[1:]
0887     w = "%+d"%wheel
0888     if w=="+0": w = "-0"
0889     return "MB%s/%d/%02d" % (w, station, sector)
0890   elif postal_address[0] == "CSC":
0891     endcap, station, ring, chamber = postal_address[1:]
0892     if endcap != 1: station = -1 * abs(station)
0893     return "ME%+d/%d/%02d" % (station, ring, chamber)
0894
0895
0896 def nameToId(name):
0897   if name[0:2] == "MB":
0898     wh = name[4]
0899     if   wh == "A": w = "-2"
0900     elif wh == "B": w = "-1"
0901     elif wh == "C": w = "-0"
0902     elif wh == "D": w = "+1"
0903     elif wh == "E": w = "+2"
0904     else: return ""
0905     station = name[7]
0906     sector = name[11:13]
0907     return "MB%s/%s/%s" % (w, station, sector)
0908   elif name[0:2] == "ME":
0909     if name[2]=="p": endcap = "+"
0910     elif name[2]=="m": endcap = "-"
0911     else: return ""
0912     station = name[3]
0913     ring = name[4]
0914     chamber = name[6:8]
0915     return "ME%s%s/%s/%s" % (endcap, station, ring, chamber)
0916   return None
0917
0918
0919 def availableCellsDT(reports):
0920   dts = []
0921   # DT wheels
0922   for iwheel in DT_TYPES:
0923     if iwheel[1]=="ALL": continue
0924     dts.append(iwheel[0])
0925   # DT wheel & station
0926   for wheel in DT_TYPES:
0927     if wheel[1]=="ALL": continue
0928     for station in wheel[2]:
0929       dts.append(wheel[0]+'/'+station[1])
0930   # DT station & sector
0931   for wheel in DT_TYPES:
0932     if wheel[1]!="ALL": continue
0933     for station in wheel[2]:
0934       for sector in range(1,station[2]+1):
0935         ssector = "%02d" % sector
0936         dts.append(wheel[0]+'/'+station[1]+'/'+ssector)
0937   # DT station & ALL sectors
0938   for wheel in DT_TYPES:
0939     if wheel[1]!="ALL": continue
0940     for station in wheel[2]:
0941         dts.append(wheel[0]+'/'+station[1])
0942   # DT chambers
0943   for wheel in DT_TYPES:
0944     if wheel[1]=="ALL": continue
0945     for station in wheel[2]:
0946       for sector in range(1,station[2]+1):
0947         ssector = "%02d" % sector
0948         label = "MBwh%sst%ssec%s" % (wheelLetter(int(wheel[1])),station[1],ssector)
0949         if len(reports)==0:
0950           # no reports case: do not include chambers
0951           #dts.append(wheel[0]+'/'+station[1]+'/'+ssector)
0952           continue
0953         found = False
0954         for r in reports:
0955           if r.name == label:
0956             found = True
0957             break
0958         if not found: continue
0959         if r.status == "TOOFEWHITS" and r.posNum+r.negNum==0: continue
0960         if r.status == "NOFIT": continue
0961         dts.append(wheel[0]+'/'+station[1]+'/'+ssector)
0962   return dts
0963
0964
0965 def availableCellsCSC(reports):
0966   cscs = []
0967   # CSC station
0968   for endcap in CSC_TYPES:
0969     for station in endcap[2]:
0970       cscs.append("%s%s" % (endcap[0], station[1]))
0971   # CSC station & ring
0972   for endcap in CSC_TYPES:
0973     for station in endcap[2]:
0974       for ring in station[2]:
0975         if ring[1]=="ALL": continue
0976         #label = "CSCvsphi_me%s%s%s" % (endcap[1], station[1], ring[1])
0977         cscs.append("%s%s/%s" % (endcap[0], station[1],ring[1]))
0978   # CSC station and chamber
0979   for endcap in CSC_TYPES:
0980     for station in endcap[2]:
0981       for ring in station[2]:
0982         if ring[1]!="ALL": continue
0983         for chamber in range(1,ring[2]+1):
0984           #label = "CSCvsr_me%s%sch%02d" % (endcap[1], station[1], chamber)
0985           cscs.append("%s%s/ALL/%02d" % (endcap[0], station[1],chamber))
0986   # CSC station and ALL chambers
0987   for endcap in CSC_TYPES:
0988     for station in endcap[2]:
0989       for ring in station[2]:
0990         if ring[1]!="ALL": continue
0991         #label = "CSCvsr_me%s%schALL" % (endcap[1], station[1])
0992         cscs.append("%s%s/ALL" % (endcap[0], station[1]))
0993   # CSC chambers
0994   for endcap in CSC_TYPES:
0995     for station in endcap[2]:
0996       for ring in station[2]:
0997         if ring[1]=="ALL": continue
0998         for chamber in range(1,ring[2]+1):
0999           # exclude non instrumented ME4/2
1000           if station[1]=="4" and ring[1]=="2":
1001             if endcap[1]=="m": continue
1002             if chamber<9 or chamber>13: continue
1003           schamber = "%02d" % chamber
1004           label = "ME%s%s%s_%s" % (endcap[1], station[1], ring[1], schamber)
1005           if len(reports)==0:
1006             # no reports case: do not include chambers
1007             #cscs.append(endcap[0]+station[1]+'/'+ring[1]+'/'+schamber)
1008             continue
1009           found = False
1010           for r in reports:
1011             if r.name == label:
1012               found = True
1013               break
1014           if not found: continue
1015           if r.status == "TOOFEWHITS" and r.posNum+r.negNum==0: continue
1016           if r.status == "NOFIT": continue
1017           cscs.append(endcap[0]+station[1]+'/'+ring[1]+'/'+schamber)
1018   return cscs
1019
1020
1021 DQM_SEVERITY = [
1022   {"idx":0, "name": "NONE", "color": "lightgreen", "hex":"#90EE90"},
1023   {"idx":1, "name": "LOWSTAT05", "color": "lightgreen", "hex":"#96D953"},
1024   {"idx":2, "name": "LOWSTAT075", "color": "lightgreen", "hex":"#94E26F"},
1025   {"idx":3, "name": "LOWSTAT1", "color": "yellowgreen", "hex":"#9ACD32"},
1026   {"idx":4, "name": "LOWSTAT", "color": "yellow", "hex":"#FFFF00"},
1027   {"idx":5, "name": "TOLERABLE", "color": "lightpink", "hex":"#FFB6C1"},
1028   {"idx":6, "name": "SEVERE", "color": "orange", "hex":"#FFA500"},
1029   {"idx":7, "name": "CRITICAL", "color": "red", "hex":"#FF0000"}];
1030
1031
1032 def addToTestResults(c,res):
1033   if len(res)>0:
1034     if c in TEST_RESULTS: TEST_RESULTS[c].extend(res)
1035     else: TEST_RESULTS[c] = res
1036
1037
1038 def testEntry(testID,scope,descr,severity):
1039   s = 0
1040   for sev in DQM_SEVERITY:
1041     if sev["name"]==severity: s = sev["idx"]
1042   return {"testID":testID,"scope":scope,"descr":descr,"severity":s}
1043
1044
1045 def testZeroWithin5Sigma(x):
1046   if abs(x[1])==0.: return 0.
1047   pull = abs(x[0])/abs(x[1])
1048   if pull <= 5: return 0.
1049   else: return pull
1050
1051
1052 def testDeltaWithin5Sigma(x,sx):
1053   n = len(x)
1054   res = []
1055   dr = []
1056   #print x
1057   #print sx
1058   for i in range(1,n+1):
1059     x1 = x[i-1]
1060     sx1 = sx[i-1]
1061     x2 = x[0]
1062     sx2 = sx[0]
1063     if i<n:
1064       x2 = x[i]
1065       sx2 = sx[i]
1066     sig1 = sqrt( (sx1[0]-sx1[1])**2 + x1[1]**2 )
1067     sig2 = sqrt( (sx2[0]-sx2[1])**2 + x2[1]**2 )
1068     df = abs(x1[0]-x2[0]) - 3*( sig1 + sig2 )
1069     #df = abs(sx1[1]-sx2[0]) - 5*(abs(x1[1]) + abs(x2[1]))
1070     #print i, df, '= abs(',sx1[1],'-',sx2[0],')-5*(abs(',x1[1],')+abs(',x2[1],'))'
1071     dr.append(df)
1072     if df > 0: res.append(i)
1073   #print dr
1074   #print res
1075   return res
1076
1077
1078 def doTestsForReport(cells,reports):
1079   for c in cells:
1080     # can a cell be converted to a chamber postal address?
1081     postal_address = idToPostalAddress(c)
1082     if not postal_address: continue
1083
1084     # is this chamber in _report?
1085     found = False
1086     for r in reports:
1087       if r.postal_address == postal_address:
1088         found = True
1089         break
1090     if not found: continue
1091
1092     # chamber's tests result
1093     res = []
1094     scope = postal_address[0]
1095
1096     # noting could be done if fitting fails
1097     if r.status == "FAIL" or r.status == "MINUITFAIL":
1098       res.append(testEntry("FAILURE",scope,r.status+" failure","CRITICAL"))
1099       addToTestResults(c,res)
1100       continue
1101
1102     # noting could be done if TOOFEWHITS
1103     nseg = r.posNum + r.negNum
1104     if r.status == "TOOFEWHITS" and nseg>0:
1105       res.append(testEntry("LOW_STAT",scope,"low stat, #segments=%d"%nseg,"LOWSTAT"))
1106       addToTestResults(c,res)
1107       continue
1108
1109     # set shades of light green according to sidma(dx)
1110     sdx = 10.*r.deltax.error
1111     if sdx>0.5:
1112       if sdx<0.75: res.append(testEntry("LOW_STAT_DDX05",scope,"low stat, delta(dx)=%f #segments=%d" % (sdx,nseg),"LOWSTAT05"))
1113       elif sdx<1.: res.append(testEntry("LOW_STAT_DDX075",scope,"low stat, delta(dx)=%f #segments=%d" % (sdx,nseg),"LOWSTAT075"))
1114       else: res.append(testEntry("LOW_STAT_DDX1",scope,"low stat, delta(dx)=%f #segments=%d" % (sdx,nseg),"LOWSTAT1"))
1115
1116     # check chi2
1117     if r.redchi2 > 20.: #2.5:
1118       res.append(testEntry("BIG_CHI2",scope,"chi2=%f>20" % r.redchi2,"TOLERABLE"))
1119
1120     # check medians
1121     medx, meddx = 10.*r.median_x, 1000.*r.median_dxdz
1122     #medy, meddy = 10.*r.median_y, 1000.*r.median_dydz
1123     if medx>2:  res.append(testEntry("BIG_MED_X",scope,"median dx=%f>2 mm"%medx,"SEVERE"))
1124     #if medy>3: res.append(testEntry("BIG_MED_Y",scope,"median dy=%f>3 mm"%medy,"SEVERE"))
1125     if meddx>2: res.append(testEntry("BIG_MED_DXDZ",scope,"median d(dx/dz)=%f>2 mrad"%meddx,"SEVERE"))
1126     #if meddy>3: res.append(testEntry("BIG_MED_DYDZ",scope,"median d(dy/dz)=%f>3 mrad"%meddy,"SEVERE"))
1127
1128     # check residuals far from zero
1129     isDTst4 = False
1130     if postal_address[0] == "DT" and postal_address[2]==4: isDTst4 = True
1131     dx, dy, dpy, dpz = 10.*r.deltax.value, 0., 1000.*r.deltaphiy.value, 1000.*r.deltaphiz.value
1132     if not isDTst4: dy = 10.*r.deltay.value
1133     if dx>0.2:   res.append(testEntry("BIG_LAST_ITR_DX",scope,"dx=%f>0.2 mm"%dx,"CRITICAL"))
1134     if dy>0.2:   res.append(testEntry("BIG_LAST_ITR_DY",scope,"dy=%f>0.2 mm"%dy,"CRITICAL"))
1135     if dpy>0.2:   res.append(testEntry("BIG_LAST_ITR_DPHIY",scope,"dphiy=%f>0.2 mrad"%dpy,"CRITICAL"))
1136     if dpz>0.2:   res.append(testEntry("BIG_LAST_ITR_DPHIZ",scope,"dphiz=%f>0.2 mrad"%dpz,"CRITICAL"))
1137     #if ddx>0.03: res.append(testEntry("BIG_DX",scope,"dphix=%f>0.03 mrad"%ddx,"CRITICAL"))
1138
1139     addToTestResults(c,res)
1140
1141
1142 def doTestsForMapPlots(cells):
1143   for c in cells:
1144     res = []
1145
1146     scope = "zzz"
1147     if c[0:2]=="MB": scope = "DT"
1148     if c[0:2]=="ME": scope = "CSC"
1149     if scope == "zzz":
1150       print("strange cell ID: ", c)
1151       return None
1152
1153     if c in MAP_RESULTS_FITSIN:
1154       t = MAP_RESULTS_FITSIN[c]
1155       t_a = testZeroWithin5Sigma(t['a'])
1156       t_s = testZeroWithin5Sigma(t['sin'])
1157       t_c = testZeroWithin5Sigma(t['cos'])
1158       if t_a+t_s+t_c >0:
1159         descr = "map fitsin 5 sigma away from 0; pulls : a=%.2f sin=%.2f, cos=%.2f" % (t_a,t_s,t_c)
1160         res.append(testEntry("MAP_FITSIN",scope,descr,"SEVERE"))
1161
1162     if c in MAP_RESULTS_SAWTOOTH:
1163       t = MAP_RESULTS_SAWTOOTH[c]
1164
1165       t_a = testDeltaWithin5Sigma(t['a'],t['da'])
1166       if len(t_a)>0:
1167         descr = "map discontinuities: %s" % ",".join(map(str,t_a))
1168         res.append(testEntry("MAP_DISCONTIN",scope,descr,"SEVERE"))
1169
1170       t_b = map(testZeroWithin5Sigma, t['b'])
1171       t_bi = []
1172       for i in range(0,len(t_b)):
1173         if t_b[i]>0: t_bi.append(i+1)
1174       if len(t_bi)>0:
1175         descr = "map sawteeth: %s" % ",".join(map(str,t_bi))
1176         res.append(testEntry("MAP_SAWTEETH",scope,descr,"TOLERABLE"))
1177
1178     addToTestResults(c,res)
1179
1180
1181 def saveTestResultsMap(run_name):
1182   if len(MAP_RESULTS_SAWTOOTH)+len(MAP_RESULTS_FITSIN)==0: return None
1183   ff = open("tmp_test_results_map__%s.pkl" % run_name, "wb")
1184   pickle.dump(MAP_RESULTS_SAWTOOTH, ff)
1185   pickle.dump(MAP_RESULTS_FITSIN, ff)
1186   ff.close()
1187
1188
1189 def loadTestResultsMap(run_name):
1190   print("tmp_test_results_map__%s.pkl" % run_name, os.access("tmp_test_results_map__%s.pkl" % run_name,os.F_OK))
1191   if not os.access("tmp_test_results_map__%s.pkl" % run_name,os.F_OK): return None
1192   global MAP_RESULTS_FITSIN, MAP_RESULTS_SAWTOOTH
1193   ff = open("tmp_test_results_map__%s.pkl" % run_name, "rb")
1194   MAP_RESULTS_SAWTOOTH = pickle.load(ff)
1195   MAP_RESULTS_FITSIN = pickle.load(ff)
1196   ff.close()
1197   #execfile("tmp_test_results_map__%s.py" % run_name)
1198   #print 'asasas', MAP_RESULTS_FITSIN
1199   return True
1200
1201 def writeDQMReport(fname_dqm, run_name):
1202   tests = []
1203   for c in TEST_RESULTS:
1204     tests.append({"objID":c, "name":c, "list":TEST_RESULTS[c]})
1205   lt = time.localtime(time.time())
1206   lts = "%04d-%02d-%02d %02d:%02d:%02d %s" % (lt[0], lt[1], lt[2], lt[3], lt[4], lt[5], time.tzname[1])
1207   dqm_report = {"run":run_name, "genDate": lts, "report":tests}
1208   ff = open(fname_dqm,mode="w")
1209   print("var DQM_REPORT = ", file=ff)
1210   json.dump(dqm_report,ff)
1211   #print >>ff, "];"
1212   ff.close()
1213
1214
1215 def doTests(reports, pic_ids, fname_base, fname_dqm, run_name):
1216   # find available baseline
1217   dts = []
1218   cscs = []
1219   if len(reports)>0:
1220     dts  = availableCellsDT(reports)
1221     cscs = availableCellsCSC(reports)
1222   elif len(pic_ids)>0:
1223     dts  = [id for id in pic_ids if 'MB' in id]
1224     cscs = [id for id in pic_ids if 'ME' in id]
1225   mulist = ['Run: '+run_name,['ALL',['MU']],['DT',dts],['CSC',cscs]]
1226   ff = open(fname_base,mode="w")
1227   print("var MU_LIST = [", file=ff)
1228   json.dump(mulist,ff)
1229   print("];", file=ff)
1230   ff.close()
1231
1232   doTestsForReport(dts,reports)
1233   doTestsForReport(cscs,reports)
1234
1235   loadTestResultsMap(run_name)
1236   doTestsForMapPlots(dts)
1237   doTestsForMapPlots(cscs)
1238
1239   writeDQMReport(fname_dqm, run_name)
1240
1241
1242 ######################################################################################################
1243
1244 def plotmedians(reports1, reports2, selection=None, binsx=100, windowx=5., ceilingx=None, binsy=100, windowy=5.,
1245                 ceilingy=None, binsdxdz=100, windowdxdz=5., ceilingdxdz=None, binsdydz=100, windowdydz=5., ceilingdydz=None,
1246                 r1text=" before", r2text=" after", which="median"):
1247     tdrStyle.SetOptStat("emrou")
1248     tdrStyle.SetStatW(0.40)
1249     tdrStyle.SetStatFontSize(0.05)
1250
1251     global hmediandxdz_after, hmediandxdz_before, hmediandxdz_beforecopy, \
1252            hmediandydz_after, hmediandydz_before, hmediandydz_beforecopy, \
1253            hmedianx_after, hmedianx_before, hmedianx_beforecopy, \
1254            hmediany_after, hmediany_before, hmediany_beforecopy, tlegend
1255
1256     hmedianx_before = ROOT.TH1F("hmedianx_before", "", binsx, -windowx, windowx)
1257     hmediany_before = ROOT.TH1F("hmediany_before", "", binsy, -windowy, windowy)
1258     hmediandxdz_before = ROOT.TH1F("hmediandxdz_before", "", binsdxdz, -windowdxdz, windowdxdz)
1259     hmediandydz_before = ROOT.TH1F("hmediandydz_before", "", binsdydz, -windowdydz, windowdydz)
1260     hmedianx_after = ROOT.TH1F("hmedianx_after", "", binsx, -windowx, windowx)
1261     hmediany_after = ROOT.TH1F("hmediany_after", "", binsy, -windowy, windowy)
1262     hmediandxdz_after = ROOT.TH1F("hmediandxdz_after", "", binsdxdz, -windowdxdz, windowdxdz)
1263     hmediandydz_after = ROOT.TH1F("hmediandydz_after", "", binsdydz, -windowdydz, windowdydz)
1264
1265     if which == "median":
1266         whichx = whichy = whichdxdz = whichdydz = "median"
1267     elif which == "bigmean":
1268         whichx = "mean30"
1269         whichy = "mean30"
1270         whichdxdz = "mean20"
1271         whichdydz = "mean50"
1272     elif which == "mean":
1273         whichx = "mean15"
1274         whichy = "mean15"
1275         whichdxdz = "mean10"
1276         whichdydz = "mean25"
1277     elif which == "bigwmean":
1278         whichx = "wmean30"
1279         whichy = "wmean30"
1280         whichdxdz = "wmean20"
1281         whichdydz = "wmean50"
1282     elif which == "wmean":
1283         whichx = "wmean15"
1284         whichy = "wmean15"
1285         whichdxdz = "wmean10"
1286         whichdydz = "wmean25"
1287     elif which == "bigstdev":
1288         whichx = "stdev30"
1289         whichy = "stdev30"
1290         whichdxdz = "stdev20"
1291         whichdydz = "stdev50"
1292     elif which == "stdev":
1293         whichx = "stdev15"
1294         whichy = "stdev15"
1295         whichdxdz = "stdev10"
1296         whichdydz = "stdev25"
1297     else:
1298         raise Exception(which + " not recognized")
1299
1300     for r1 in reports1:
1301         if selection is None or (selection.__code__.co_argcount == len(r1.postal_address) and selection(*r1.postal_address)):
1302             found = False
1303             for r2 in reports2:
1304                 if r1.postal_address == r2.postal_address:
1305                     found = True
1306                     break
1307             if not found: continue
1308
1309             #skip ME1/1a
1310             if r1.postal_address[0]=='CSC':
1311               if r1.postal_address[2]==1 and r1.postal_address[3]==4: continue
1312
1313             if r1.status == "PASS" and r2.status == "PASS":
1314                 hmedianx_before.Fill(10.*eval("r1.%s_x" % whichx))
1315                 hmediandxdz_before.Fill(1000.*eval("r1.%s_dxdz" % whichdxdz))
1316                 hmedianx_after.Fill(10.*eval("r2.%s_x" % whichx))
1317                 hmediandxdz_after.Fill(1000.*eval("r2.%s_dxdz" % whichdxdz))
1318
1319                 if r1.median_y is not None:
1320                     hmediany_before.Fill(10.*eval("r1.%s_y" % whichy))
1321                     hmediandydz_before.Fill(1000.*eval("r1.%s_dydz" % whichdydz))
1322                     hmediany_after.Fill(10.*eval("r2.%s_y" % whichy))
1323                     hmediandydz_after.Fill(1000.*eval("r2.%s_dydz" % whichdydz))
1324
1325     hmedianx_beforecopy = hmedianx_before.Clone()
1326     hmediany_beforecopy = hmediany_before.Clone()
1327     hmediandxdz_beforecopy = hmediandxdz_before.Clone()
1328     hmediandydz_beforecopy = hmediandydz_before.Clone()
1329     hmedianx_beforecopy.SetLineStyle(2)
1330     hmediany_beforecopy.SetLineStyle(2)
1331     hmediandxdz_beforecopy.SetLineStyle(2)
1332     hmediandydz_beforecopy.SetLineStyle(2)
1333
1334     hmedianx_before.SetFillColor(ROOT.kMagenta+2)
1335     hmediany_before.SetFillColor(ROOT.kMagenta+2)
1336     hmediandxdz_before.SetFillColor(ROOT.kMagenta+2)
1337     hmediandydz_before.SetFillColor(ROOT.kMagenta+2)
1338     hmedianx_after.SetFillColor(ROOT.kYellow)
1339     hmediany_after.SetFillColor(ROOT.kYellow)
1340     hmediandxdz_after.SetFillColor(ROOT.kYellow)
1341     hmediandydz_after.SetFillColor(ROOT.kYellow)
1342
1343     hmedianx_after.SetXTitle("median(#Deltax) (mm)")
1344     hmediany_after.SetXTitle("median(#Deltay) (mm)")
1345     hmediandxdz_after.SetXTitle("median(#Deltadx/dz) (mrad)")
1346     hmediandydz_after.SetXTitle("median(#Deltadydz) (mrad)")
1347     hmedianx_after.GetXaxis().CenterTitle()
1348     hmediany_after.GetXaxis().CenterTitle()
1349     hmediandxdz_after.GetXaxis().CenterTitle()
1350     hmediandydz_after.GetXaxis().CenterTitle()
1351
1352     if ceilingx is not None: hmedianx_after.SetAxisRange(0., ceilingx, "Y")
1353     if ceilingy is not None: hmediany_after.SetAxisRange(0., ceilingy, "Y")
1354     if ceilingdxdz is not None: hmediandxdz_after.SetAxisRange(0., ceilingdxdz, "Y")
1355     if ceilingdydz is not None: hmediandydz_after.SetAxisRange(0., ceilingdydz, "Y")
1356
1357     c1.Clear()
1358     c1.Divide(2, 2)
1359
1360     c1.GetPad(1).cd()
1361     hmedianx_after.Draw()
1362     hmedianx_before.Draw("same")
1363     hmedianx_after.Draw("same")
1364     hmedianx_beforecopy.Draw("same")
1365     hmedianx_after.Draw("axissame")
1366
1367     tlegend = ROOT.TLegend(0.17, 0.75-0.05, 0.45+0.05, 0.9)
1368     tlegend.SetFillColor(ROOT.kWhite)
1369     tlegend.SetBorderSize(0)
1370     tlegend.AddEntry(hmedianx_after, r2text, "f")
1371     tlegend.AddEntry(hmedianx_before, r1text, "f")
1372     tlegend.Draw()
1373
1374     c1.GetPad(2).cd()
1375     hmediandxdz_after.Draw()
1376     hmediandxdz_before.Draw("same")
1377     hmediandxdz_after.Draw("same")
1378     hmediandxdz_beforecopy.Draw("same")
1379     hmediandxdz_after.Draw("axissame")
1380
1381     c1.GetPad(3).cd()
1382     hmediany_after.Draw()
1383     hmediany_before.Draw("same")
1384     hmediany_after.Draw("same")
1385     hmediany_beforecopy.Draw("same")
1386     hmediany_after.Draw("axissame")
1387
1388     c1.GetPad(4).cd()
1389     hmediandydz_after.Draw()
1390     hmediandydz_before.Draw("same")
1391     hmediandydz_after.Draw("same")
1392     hmediandydz_beforecopy.Draw("same")
1393     hmediandydz_after.Draw("axissame")
1394
1395     return hmediandxdz_after, hmediandxdz_before, hmediandxdz_beforecopy, \
1396            hmediandydz_after, hmediandydz_before, hmediandydz_beforecopy, \
1397            hmedianx_after,    hmedianx_before, hmedianx_beforecopy, \
1398            hmediany_after,    hmediany_before, hmediany_beforecopy, tlegend
1399
1400
1401 ######################################################################################################
1402
1403 def createPeaksProfile(the2d, rebin=1):
1404   htmp = ROOT.gROOT.FindObject(the2d.GetName()+"_peaks")
1405   if htmp != None: htmp.Delete()
1406
1407   hpeaks = the2d.ProjectionX(the2d.GetName()+"_peaks")
1408   hpeaks.Reset()
1409   hpeaks.Rebin(rebin)
1410   bad_fit_bins = []
1411   for i in range(0, int(the2d.GetNbinsX()), rebin):
1412     tmp = the2d.ProjectionY("tmp", i+1, i + rebin)
1413     nn = tmp.GetEntries()
1414
1415     drange = tmp.GetRMS()
1416     drange = 2.*drange
1417     fgaus = ROOT.TF1("fgaus","gaus", tmp.GetMean() - drange, tmp.GetMean() + drange)
1418     fgaus.SetParameter(0,nn)
1419     fgaus.SetParameter(1,tmp.GetMean())
1420     fgaus.SetParameter(2,tmp.GetRMS())
1421     #print "  ", i, nn, tmp.GetMean() , drange, "[", tmp.GetMean() - drange, tmp.GetMean() + drange, ']'
1422
1423     fitOk = False
1424     if nn > 10:     # good to fit
1425       fr = tmp.Fit("fgaus","RNSQ")
1426       #print "       ", fgaus.GetParameter(1), " +- ", fgaus.GetParError(1), "   fitres = " , fr.Status() , fr.CovMatrixStatus()
1427       hpeaks.SetBinContent(i/rebin+1, fgaus.GetParameter(1))
1428       hpeaks.SetBinError(i/rebin+1, fgaus.GetParError(1))
1429       if fr.Status()==0 and fr.CovMatrixStatus()==3 : fitOk = True
1430     if not fitOk:
1431       bad_fit_bins.append(i/rebin+1)
1432       if nn > 1. and tmp.GetRMS() > 0: # use mean
1433         hpeaks.SetBinContent(i/rebin+1, tmp.GetMean())
1434         hpeaks.SetBinError(i/rebin+1, ROOT.TMath.StudentQuantile(0.841345,nn) * tmp.GetRMS() / sqrt(nn))
1435       else:
1436         hpeaks.SetBinContent(i/rebin+1, 0.)
1437         hpeaks.SetBinError(i/rebin+1, 0.)
1438   if len(bad_fit_bins): print("createPeaksProfile bad fit bins: ", bad_fit_bins)
1439   return hpeaks
1440
1441
1442 ######################################################################################################
1443
1444 def mapplot(tfiles, name, param, mode="from2d", window=10., abscissa=None, title="",
1445             widebins=False, fitsine=False, fitline=False, reset_palette=False, fitsawteeth=False, fitpeaks=False, peaksbins=1, fixfitpars={}, **args):
1446     tdrStyle.SetOptTitle(1)
1447     tdrStyle.SetTitleBorderSize(0)
1448     tdrStyle.SetOptStat(0)
1449     #tdrStyle.SetOptStat("emrou")
1450     tdrStyle.SetOptFit(0)
1451     tdrStyle.SetTitleFontSize(0.05)
1452     tdrStyle.SetPadRightMargin(0.1) # to see the pallete labels on the left
1453
1454     c1.Clear()
1455     c1.ResetAttPad()
1456
1457     if reset_palette: set_palette("blues")
1458     global hist, hist2d, hist2dweight, tline1, tline2, tline3
1459
1460     if fitsine or fitsawteeth:
1461         id = mapNameToId(name)
1462         if not id:
1463             print("bad id for ", name)
1464             raise Exception
1465
1466     hdir = "AlignmentMonitorMuonSystemMap1D/iter1/"
1467     hpref= "%s_%s" % (name, param)
1468     hhh  = hdir+hpref
1469
1470     combine_all = False
1471     if "ALL" in name  and  ("CSCvsr" in name  or  "DTvsz" in name): combine_all = True
1472
1473     add1d =  ("vsphi" in name) and (param == "x")
1474
1475     if "h2d" in args:
1476       hist2d = args["h2d"].Clone(hpref+"_2d_")
1477       if "CSC" in name  and  add1d: hist1d = args["h1d"].Clone(hpref+"_1d_")
1478
1479     elif combine_all:
1480       nch = 12
1481       if  "DT" in name  and  name[6:9]=='st4': nch = 14
1482       if  "CSC" in name: nch = 36
1483       chambers = ["%02d" % ch for ch in range (2,nch+1)]
1484
1485       ch_hhh = hhh.replace('ALL','01')
1486       ch_hpref = hpref.replace('ALL','01')
1487       hist2d = tfiles[0].Get(ch_hhh+"_2d").Clone(ch_hpref+"_2d_")
1488       if "CSC" in name  and  add1d: hist1d = tfiles[0].Get(ch_hhh+"_1d").Clone(ch_hpref+"_1d_")
1489
1490       for ch in chambers:
1491         ch_hhh = hhh.replace('ALL',ch)
1492         ch_hpref = hpref.replace('ALL',ch)
1493         hist2d.Add(tfiles[0].Get(ch_hhh+"_2d"))
1494         if "CSC" in name  and  add1d: hist1d.Add(tfiles[0].Get(ch_hhh+"_1d"))
1495         for tfile in tfiles[1:]:
1496           hist2d.Add(tfile.Get(ch_hhh+"_2d"))
1497           if "CSC" in name   and  add1d: hist1d.Add(tfile.Get(ch_hhh+"_1d"))
1498
1499     else:
1500       hist2d = tfiles[0].Get(hhh+"_2d").Clone(hpref+"_2d_")
1501       if "CSC" in name  and  add1d: hist1d = tfiles[0].Get(hhh+"_1d").Clone(hpref+"_1d_")
1502       for tfile in tfiles[1:]:
1503         hist2d.Add(tfile.Get(hhh+"_2d"))
1504         if "CSC" in name   and  add1d: hist1d.Add(tfile.Get(hhh+"_1d"))
1505
1506
1507     if mode == "from2d":
1508         the2d = hist2d
1509
1510         hist = the2d.ProjectionX()
1511         hist.Reset()
1512
1513         skip = 1
1514         if widebins:
1515             hist.Rebin(10)
1516             skip = 10
1517
1518         #f = ROOT.TF1("g", "gaus", -40., 40)
1519         for i in range(0, int(the2d.GetNbinsX()), skip):
1520             tmp = the2d.ProjectionY("tmp", i+1, i + skip)
1521             if tmp.GetEntries() > 1:
1522                 #tmp.Fit("g","LNq")
1523                 hist.SetBinContent(i/skip+1, tmp.GetMean())
1524                 hist.SetBinError(i/skip+1, ROOT.TMath.StudentQuantile(0.841345,tmp.GetEntries()) * tmp.GetRMS() / sqrt(tmp.GetEntries()))
1525                 #hist.SetBinError(i/skip+1, tmp.GetRMS() / sqrt(tmp.GetEntries()))
1526                 #hist.SetBinError(i/skip+1, f.GetParameter(2))
1527             else:
1528                 #hist.SetBinContent(i/skip+1, 2000.)
1529                 #hist.SetBinError(i/skip+1, 1000.)
1530                 hist.SetBinContent(i/skip+1, 0.)
1531                 hist.SetBinError(i/skip+1, 0.)
1532
1533         hpeaks = createPeaksProfile(the2d, peaksbins)
1534
1535     else:
1536         raise Exception
1537
1538     hist.SetAxisRange(-window, window, "Y")
1539     if abscissa is not None: hist.SetAxisRange(abscissa[0], abscissa[1], "X")
1540     hist.SetMarkerStyle(20)
1541     hist.SetMarkerSize(0.75)
1542     hist.GetXaxis().CenterTitle()
1543     hist.GetYaxis().CenterTitle()
1544     hist.GetYaxis().SetTitleOffset(0.75)
1545     hist.GetXaxis().SetTitleSize(0.05)
1546     hist.GetYaxis().SetTitleSize(0.05)
1547     hist.SetTitle(title)
1548     if "vsphi" in name: hist.SetXTitle("Global #phi position (rad)")
1549     elif "vsz" in name: hist.SetXTitle("Global z position (cm)")
1550     elif "vsr" in name: hist.SetXTitle("Global R position (cm)")
1551     if "DT" in name:
1552         if param == "x": hist.SetYTitle("x' residual (mm)")
1553         if param == "dxdz": hist.SetYTitle("dx'/dz residual (mrad)")
1554         if param == "y": hist.SetYTitle("y' residual (mm)")
1555         if param == "dydz": hist.SetYTitle("dy'/dz residual (mrad)")
1556     if "CSC" in name:
1557         if param == "x": hist.SetYTitle("r#phi residual (mm)")
1558         if param == "dxdz": hist.SetYTitle("d(r#phi)/dz residual (mrad)")
1559     hist.SetMarkerColor(ROOT.kBlack)
1560     hist.SetLineColor(ROOT.kBlack)
1561     hist.Draw()
1562     hist2d.Draw("colzsame")
1563     if widebins: hist.Draw("samee1")
1564     else: hist.Draw("same")
1565
1566     hpeaks.SetMarkerStyle(20)
1567     hpeaks.SetMarkerSize(0.9)
1568     hpeaks.SetMarkerColor(ROOT.kRed)
1569     hpeaks.SetLineColor(ROOT.kRed)
1570     hpeaks.SetLineWidth(2)
1571     #if fitpeaks: hpeaks.Draw("same")
1572     hpeaks.Draw("same")
1573
1574     if fitsine and "vsphi" in name:
1575         global fitsine_const, fitsine_sin, fitsine_cos, fitsine_chi2, fitsine_ndf
1576         if 'CSC' in name:
1577           f = ROOT.TF1("f", "[0] + [1]*sin(x) + [2]*cos(x)", -pi/180.*5., pi*(2.-5./180.))
1578         else:
1579           f = ROOT.TF1("f", "[0] + [1]*sin(x) + [2]*cos(x)", -pi, pi)
1580         f.SetLineColor(ROOT.kRed)
1581         f.SetLineWidth(2)
1582         if len(fixfitpars)>0:
1583           for fpar in fixfitpars.keys():
1584             f.FixParameter(fpar, fixfitpars[fpar])
1585         #hist.Fit(f,"N")
1586         if fitpeaks: hpeaks.Fit(f,"NQ")
1587         else: hist.Fit(f,"NEQ")
1588         if len(fixfitpars)>0:
1589           for fpar in fixfitpars.keys():
1590             f.ReleaseParameter(fpar)
1591         fitsine_const = f.GetParameter(0), f.GetParError(0)
1592         fitsine_sin = f.GetParameter(1), f.GetParError(1)
1593         fitsine_cos = f.GetParameter(2), f.GetParError(2)
1594         fitsine_chi2 = f.GetChisquare()
1595         fitsine_ndf = f.GetNDF()
1596         global MAP_RESULTS_FITSIN
1597         # 'phi' coefficienct will be updated further for CSC
1598         MAP_RESULTS_FITSIN[id] = {'a':fitsine_const, 'phi':fitsine_const, 'sin': fitsine_sin, 'cos': fitsine_cos, 'chi2': fitsine_chi2, 'ndf': fitsine_ndf}
1599         f.Draw("same")
1600         global fitsine_ttext, fitsine_etext
1601         text_xposition = -1.
1602         if 'CSC' in name: text_xposition = 2.
1603         fitsine_ttext = ROOT.TLatex(text_xposition, 0.8*window,
1604                 "%+.3f %+.3f sin#phi %+.3f cos#phi" % (fitsine_const[0], fitsine_sin[0], fitsine_cos[0]))
1605         fitsine_ttext.SetTextColor(ROOT.kRed)
1606         fitsine_ttext.SetTextSize(0.05)
1607         fitsine_ttext.Draw()
1608         fitsine_etext = ROOT.TLatex(text_xposition, 0.70*window,
1609                 " #pm%.3f    #pm%.3f            #pm%.3f" % (fitsine_const[1], fitsine_sin[1], fitsine_cos[1]))
1610         fitsine_etext.SetTextColor(ROOT.kRed)
1611         fitsine_etext.SetTextSize(0.045)
1612         fitsine_etext.Draw()
1613
1614         # additional estimate of phiz ring rotation from 1d distribution
1615         if 'CSC' in name and add1d:
1616           # zero-order rough fit to obtain the fitting range:
1617           f0 = ROOT.TF1("f0", "gaus", hist1d.GetBinLowEdge(1), -hist1d.GetBinLowEdge(1))
1618           fit = hist1d.Fit(f0,"NRQ")
1619           rangea, rangeb = hist1d.GetMean() - hist1d.GetRMS(), hist1d.GetMean() + hist1d.GetRMS()
1620           if fit==0: rangea, rangeb = f0.GetParameter(1) - f0.GetParameter(2), f0.GetParameter(1) + f0.GetParameter(2)
1621           #print rangea, rangeb
1622
1623           # second fit for finding the peak:
1624           f1 = ROOT.TF1("f1", "gaus", rangea, rangeb)
1625           fit = hist1d.Fit(f1,"NRQ")
1626           nn = hist1d.GetEntries()
1627           dphiz, ephiz = 0, 0
1628           if nn>0:  dphiz, ephiz = hist1d.GetMean(), ROOT.TMath.StudentQuantile(0.841345,nn) * hist1d.GetRMS() / sqrt(nn)
1629           if fit==0: dphiz, ephiz = f1.GetParameter(1), f1.GetParError(1)
1630           #print dphiz, ephiz
1631           MAP_RESULTS_FITSIN[id]['phi'] = (dphiz, ephiz)
1632
1633           global ttex_sine_, ttex_sine, ttex_1d_, ttex_1d
1634           postal_address = idToPostalAddress(id+'/01')
1635           ttex_sine_ = ROOT.TLatex(0, 0.8*window,"#Delta#phi_{z}^{sine} (mrad):")
1636           ttex_sine_.SetTextColor(ROOT.kGreen+2); ttex_sine_.SetTextSize(0.04); ttex_sine_.Draw()
1637           ttex_sine = ROOT.TLatex(0, 0.7*window,"   %+.3f#pm%.3f" %
1638                                   (-100*fitsine_const[0]/signConventions[postal_address][3],
1639                                    100*fitsine_const[1]/signConventions[postal_address][3]))
1640           ttex_sine.SetTextColor(ROOT.kGreen+2); ttex_sine.SetTextSize(0.04); ttex_sine.Draw()
1641           ttex_1d_ = ROOT.TLatex(0, 0.6*window,"#Delta#phi_{z}^{phi} (mrad):")
1642           ttex_1d_.SetTextColor(ROOT.kGreen+2); ttex_1d_.SetTextSize(0.04); ttex_1d_.Draw()
1643           ttex_1d = ROOT.TLatex(0, 0.5*window,"   %+.3f#pm%.3f" % (-dphiz, ephiz))
1644           ttex_1d.SetTextColor(ROOT.kGreen+2); ttex_1d.SetTextSize(0.04); ttex_1d.Draw()
1645           ROOT.gPad.Update()
1646
1647     if fitline:
1648         f = ROOT.TF1("f", "[0] + [1]*x", -1000., 1000.)
1649         hist2d.Fit(f, "q")
1650         hist2d.GetFunction("f").SetLineColor(ROOT.kRed)
1651         global fitline_const, fitline_linear, fitline_chi2, fitline_ndf
1652         fitline_const = hist2d.GetFunction("f").GetParameter(0), hist2d.GetFunction("f").GetParError(0)
1653         fitline_linear = hist2d.GetFunction("f").GetParameter(1), hist2d.GetFunction("f").GetParError(1)
1654         fitline_chi2 = hist2d.GetFunction("f").GetChisquare()
1655         fitline_ndf = hist2d.GetFunction("f").GetNDF()
1656         hist2d.GetFunction("f").Draw("same")
1657         global fitline_ttext
1658         if "vsz" in name: which = "Z"
1659         elif "vsr" in name: which = "R"
1660         fitline_ttext = ROOT.TText(hist.GetBinCenter(hist.GetNbinsX()/4),
1661                 0.8*window, "%.3g %+.3g %s" % (fitline_const[0], fitline_linear[0], which))
1662         fitline_ttext.SetTextColor(ROOT.kRed)
1663         fitline_ttext.Draw()
1664
1665     ROOT.gPad.RedrawAxis()
1666
1667     if "vsphi" in name:
1668         if not widebins: philines(name, window, abscissa)
1669         if abscissa is None:
1670           if 'CSC' in name:
1671             tline1 = ROOT.TLine(-pi/180.*5., 0, pi*(2.-5./180.), 0); tline1.Draw()
1672             tline2 = ROOT.TLine(-pi/180.*5., -window, pi*(2.-5./180.), -window); tline2.SetLineWidth(2); tline2.Draw()
1673             tline3 = ROOT.TLine(-pi/180.*5., window, pi*(2.-5./180.), window); tline3.Draw()
1674           else:
1675             tline1 = ROOT.TLine(-pi, 0, pi, 0); tline1.Draw()
1676             tline2 = ROOT.TLine(-pi, -window, pi, -window); tline2.SetLineWidth(2); tline2.Draw()
1677             tline3 = ROOT.TLine(-pi, window, pi, window); tline3.Draw()
1678         else:
1679             tline1 = ROOT.TLine(abscissa[0], 0, abscissa[1], 0); tline1.Draw()
1680             tline2 = ROOT.TLine(abscissa[0], -window, abscissa[1], -window); tline2.SetLineWidth(2); tline2.Draw()
1681             tline3 = ROOT.TLine(abscissa[0], window, abscissa[1], window); tline3.Draw()
1682     elif "vsz" in name:
1683         if not widebins: zlines(window, abscissa)
1684         if abscissa is None:
1685             tline1 = ROOT.TLine(-660, 0, 660, 0); tline1.Draw()
1686             tline2 = ROOT.TLine(-660, -window, 660, -window); tline2.SetLineWidth(2); tline2.Draw()
1687             tline3 = ROOT.TLine(-660, window, 660, window); tline3.Draw()
1688         else:
1689             tline1 = ROOT.TLine(abscissa[0], 0, abscissa[1], 0); tline1.Draw()
1690             tline2 = ROOT.TLine(abscissa[0], -window, abscissa[1], -window); tline2.SetLineWidth(2); tline2.Draw()
1691             tline3 = ROOT.TLine(abscissa[0], window, abscissa[1], window); tline3.Draw()
1692     elif "vsr" in name:
1693         if "mem1" in name or "mep1" in name and not widebins: rlines(1, window, abscissa)
1694         if "mem2" in name or "mep2" in name and not widebins: rlines(2, window, abscissa)
1695         if "mem3" in name or "mep3" in name and not widebins: rlines(3, window, abscissa)
1696         if "mem4" in name or "mep4" in name and not widebins: rlines(4, window, abscissa)
1697         if abscissa is None:
1698             tline1 = ROOT.TLine(100, 0, 700, 0); tline1.Draw()
1699             tline2 = ROOT.TLine(100, -window, 700, -window); tline2.SetLineWidth(2); tline2.Draw()
1700             tline3 = ROOT.TLine(100, window, 700, window); tline3.Draw()
1701         else:
1702             tline1 = ROOT.TLine(abscissa[0], 0, abscissa[1], 0); tline1.Draw()
1703             tline2 = ROOT.TLine(abscissa[0], -window, abscissa[1], -window); tline2.SetLineWidth(2); tline2.Draw()
1704             tline3 = ROOT.TLine(abscissa[0], window, abscissa[1], window); tline3.Draw()
1705
1706     if "vsphi" in name and fitsawteeth:
1707         global CPP_LOADED
1708         if not CPP_LOADED:
1709             phiedges2c()
1710             ROOT.gROOT.ProcessLine(".L phiedges_fitfunctions.C++")
1711             CPP_LOADED = True
1712         fn={0: ROOT.fitf0,
1713             1: ROOT.fitf2,
1714             2: ROOT.fitf2,
1715             3: ROOT.fitf3,
1716             4: ROOT.fitf4,
1717             5: ROOT.fitf5,
1718             6: ROOT.fitf6,
1719             7: ROOT.fitf7,
1720             8: ROOT.fitf8,
1721             9: ROOT.fitf9,
1722             10: ROOT.fitf10,
1723             11: ROOT.fitf11,
1724             12: ROOT.fitf12,
1725             13: ROOT.fitf13
1726         } [stationIndex(name)]
1727         fn.SetNpx(5000)
1728         fn.SetLineColor(ROOT.kYellow)
1729         hist.Fit(fn,"N")
1730         fn.Draw("same")
1731
1732         # get properly arranged phi edges
1733         edges = (phiedges[stationIndex(name)])[:]
1734         ed = sorted(edges)
1735         # add some padding to the end
1736         ed.append(pi+abs(ed[0]))
1737
1738         global sawtooth_a, sawtooth_b
1739         sawtooth_a = []
1740         sawtooth_da = []
1741         sawtooth_b = []
1742         for pr in range(0,fn.GetNpar(),2):
1743             sawtooth_a.append( (fn.GetParameter(pr), fn.GetParError(pr)) )
1744             sawtooth_b.append( (fn.GetParameter(pr+1), fn.GetParError(pr+1)) )
1745             sawtooth_da.append( (fn.Eval(ed[pr/2]+0.01), fn.Eval(ed[pr/2+1]-0.01)) )
1746         global MAP_RESULTS_SAWTOOTH
1747         MAP_RESULTS_SAWTOOTH[id] = {'a': sawtooth_a, 'da': sawtooth_da, 'b': sawtooth_b, 'chi2': fn.GetChisquare(), 'ndf': fn.GetNDF()}
1748
1749     # fill number of contributiong bins
1750
1751
1752     #ROOT.SetOwnership(hist,False)
1753     ROOT.SetOwnership(hist2d,False)
1754     ROOT.SetOwnership(hist,False)
1755     ROOT.SetOwnership(tline1,False)
1756     ROOT.SetOwnership(tline2,False)
1757     ROOT.SetOwnership(tline3,False)
1758     return hist
1759
1760
1761 def mapNameToId(name):
1762   if "DT" in name:
1763     wh = "-ALL"
1764     if name.find('wh')>1: wh = name[name.find('wh')+2]
1765     if   wh == "A": w = "-2"
1766     elif wh == "B": w = "-1"
1767     elif wh == "C": w = "-0"
1768     elif wh == "D": w = "+1"
1769     elif wh == "E": w = "+2"
1770     elif wh == "-ALL": w = "-ALL"
1771     else: return None
1772     station=''
1773     if wh == "-ALL":
1774         if name.find('sec')<0: return None
1775         station = name[name.find('sec')-1]
1776         sector = ''
1777         sector = name[name.find('sec')+3:name.find('sec')+5]
1778         return "MB%s/%s/%s" % (w, station, sector)
1779     if name.find('st')>1: station = name[name.find('st')+2]
1780     else: return None
1781     return "MB%s/%s" % (w, station)
1782   elif "CSC" in name:
1783     p = name.find('me')
1784     if p<0: return None
1785     if name[p+2]=="p": endcap = "+"
1786     elif name[p+2]=="m": endcap = "-"
1787     else: return None
1788     station = name[p+3]
1789     pch = name.find('ch')
1790     if pch<0:
1791         ring = name[p+4]
1792         return "ME%s%s/%s" % (endcap, station, ring)
1793     ring = 'ALL'
1794     chamber = name[pch+2:pch+4]
1795     return "ME%s%s/%s/%s" % (endcap, station, ring, chamber)
1796   return None
1797
1798
1799 ##################################################################################
1800 # "param" may be one of "deltax" (Delta x position residuals),
1801 #      "deltadxdz" (Delta (dx/dz) angular residuals),
1802 #      "curverr" (Delta x * d(Delta q/pT)/d(Delta x) = Delta q/pT in the absence of misalignment)
1803
1804 def curvatureplot(tfiles, name, param, mode="from2d", window=15., widebins=False, title="", fitgauss=False, fitconst=False, fitline=False, fitpeaks=True, reset_palette=False):
1805     tdrStyle.SetOptTitle(1)
1806     tdrStyle.SetTitleBorderSize(0)
1807     tdrStyle.SetOptStat(0)
1808     tdrStyle.SetOptFit(0)
1809     tdrStyle.SetTitleFontSize(0.05)
1810
1811     c1.Clear()
1812     if reset_palette: set_palette("blues")
1813     global hist, histCOPY, hist2d, tline1, tline2, tline3, tline4, tline5
1814
1815     hdir = "AlignmentMonitorMuonVsCurvature/iter1/"
1816
1817     if name not in ("all", "top", "bottom"):
1818         hsuffix = "_%s_%s" % (name, param)
1819         prof = tfiles[0].Get(hdir+"tprofile"+hsuffix).Clone("tprofile_"+hsuffix)
1820         hist2d = tfiles[0].Get(hdir+"th2f"+hsuffix).Clone("th2f_"+hsuffix)
1821         for tfile in tfiles[1:]:
1822             prof.Add(tfile.Get(hdir+"tprofile"+hsuffix))
1823             hist2d.Add(tfile.Get(hdir+"th2f"+hsuffix))
1824     else:
1825         prof = None
1826         hist2d = None
1827         for wheel in "m2", "m1", "z", "p1", "p2":
1828             if name == "all": sectors = "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"
1829             elif name == "top": sectors = "01", "02", "03", "04", "05", "06"
1830             elif name == "bottom": sectors = "07", "08", "09", "10", "11", "12"
1831             else: raise Exception
1832
1833             for sector in sectors:
1834                 hsuffix = "_%s_%s" % ("wheel%s_sector%s" % (wheel, sector), param)
1835                 for tfile in tfiles:
1836                     if prof is None:
1837                         prof = tfiles[0].Get(hdir+"tprofile"+hsuffix).Clone("tprofile_"+hsuffix)
1838                         hist2d = tfiles[0].Get(hdir+"th2f"+hsuffix).Clone("tprofile_"+hsuffix)
1839                     else:
1840                         prof.Add(tfile.Get(hdir+"tprofile"+hsuffix))
1841                         hist2d.Add(tfile.Get(hdir+"th2f"+hsuffix))
1842
1843     hist = ROOT.TH1F("hist", "", prof.GetNbinsX(), prof.GetBinLowEdge(1), -prof.GetBinLowEdge(1))
1844     for i in range(1, prof.GetNbinsX()+1):
1845         hist.SetBinContent(i, prof.GetBinContent(i))
1846         hist.SetBinError(i, prof.GetBinError(i))
1847
1848     if mode == "plain":
1849         hist = prof
1850     elif mode == "from2d":
1851         skip = 1
1852         if widebins:
1853             hist.Rebin(5)
1854             skip = 5
1855         htmp = ROOT.gROOT.FindObject("tmp")
1856         if htmp != None: htmp.Delete()
1857
1858         for i in range(0, int(prof.GetNbinsX()), skip):
1859             tmp = hist2d.ProjectionY("tmp", i+1, i + skip)
1860             if tmp.GetEntries() > 1:
1861                 hist.SetBinContent(i/skip+1, tmp.GetMean())
1862                 hist.SetBinError(i/skip+1, ROOT.TMath.StudentQuantile(0.841345,tmp.GetEntries()) * tmp.GetRMS() / sqrt(tmp.GetEntries()))
1863                 #hist.SetBinError(i/skip+1, tmp.GetRMS() / sqrt(tmp.GetEntries()))
1864             else:
1865                 #hist.SetBinContent(i/skip+1, 2000.)
1866                 #hist.SetBinError(i/skip+1, 1000.)
1867                 hist.SetBinContent(i/skip+1, 0.)
1868                 hist.SetBinError(i/skip+1, 0.)
1869
1870         hpeaks = createPeaksProfile(hist2d, skip)
1871
1872     else:
1873         raise Exception
1874
1875
1876     if fitgauss:
1877         f = ROOT.TF1("f", "[0] + [1]*exp(-x**2/2/0.01**2)", hist.GetBinLowEdge(1), -hist.GetBinLowEdge(1))
1878         f.SetParameters(0, 0., 0.01)
1879         if fitpeaks: hpeaks.Fit(f, "q")
1880         else: hist.Fit(f, "q")
1881         f.SetLineColor(ROOT.kRed)
1882         global fitgauss_diff, fitgauss_chi2, fitgauss_ndf
1883 #         fitter = ROOT.TVirtualFitter.GetFitter()
1884 #         fitgauss_diff = f.GetParameter(0) - f.GetParameter(1), \
1885 #                         sqrt(f.GetParError(0)**2 + f.GetParError(1)**2 + 2.*fitter.GetCovarianceMatrixElement(0, 1))
1886         fitgauss_diff = f.GetParameter(1), f.GetParError(1)
1887         fitgauss_chi2 = f.GetChisquare()
1888         fitgauss_ndf = f.GetNDF()
1889
1890     global fitline_intercept, fitline_slope
1891     if fitconst:
1892         f = ROOT.TF1("f", "[0]", hist.GetBinLowEdge(1), -hist.GetBinLowEdge(1))
1893         if fitpeaks: hpeaks.Fit(f, "q")
1894         else: hist.Fit(f, "q")
1895         f.SetLineColor(ROOT.kRed)
1896         fitline_intercept = f.GetParameter(0), f.GetParError(0)
1897
1898     if fitline:
1899         f = ROOT.TF1("f", "[0] + [1]*x", hist.GetBinLowEdge(1), -hist.GetBinLowEdge(1))
1900         if fitpeaks: hpeaks.Fit(f, "qNE")
1901         else: hist.Fit(f, "qNE")
1902         f.SetLineColor(ROOT.kRed)
1903         global f2, f3
1904         f2 = ROOT.TF1("2", "[0] + [1]*x", hist.GetBinLowEdge(1), -hist.GetBinLowEdge(1))
1905         f3 = ROOT.TF1("2", "[0] + [1]*x", hist.GetBinLowEdge(1), -hist.GetBinLowEdge(1))
1906         f2.SetParameters(f.GetParameter(0), f.GetParameter(1) + f.GetParError(1))
1907         f3.SetParameters(f.GetParameter(0), f.GetParameter(1) - f.GetParError(1))
1908         f2.SetLineColor(ROOT.kRed)
1909         f3.SetLineColor(ROOT.kRed)
1910         f2.SetLineStyle(2)
1911         f3.SetLineStyle(2)
1912         fitline_intercept = f.GetParameter(0), f.GetParError(0)
1913         fitline_slope = f.GetParameter(1), f.GetParError(1)
1914
1915     hist2d.SetAxisRange(-window, window, "Y")
1916     hist2d.SetMarkerStyle(20)
1917     hist2d.SetMarkerSize(0.75)
1918     hist2d.GetXaxis().CenterTitle()
1919     hist2d.GetYaxis().CenterTitle()
1920     if param == "curverr":
1921         hist2d.GetYaxis().SetTitleOffset(1.35)
1922     else:
1923         hist2d.GetYaxis().SetTitleOffset(0.75)
1924     hist2d.GetXaxis().SetTitleOffset(1.2)
1925     hist2d.GetXaxis().SetTitleSize(0.05)
1926     hist2d.GetYaxis().SetTitleSize(0.05)
1927     hist2d.SetTitle(title)
1928     if param == "pterr": hist2d.SetXTitle("qp_{T} (GeV/c)")
1929     else: hist2d.SetXTitle("q/p_{T} (c/GeV)")
1930     if param == "deltax": hist2d.SetYTitle("#Deltax' (mm)")
1931     if param == "deltadxdz": hist2d.SetYTitle("#Deltadx'/dz (mrad)")
1932     if param == "pterr": hist2d.SetYTitle("#Deltap_{T}/p_{T} (%)")
1933     if param == "curverr": hist2d.SetYTitle("#Deltaq/p_{T} (c/GeV)")
1934     hist2d.Draw("colz")
1935     hist.SetMarkerColor(ROOT.kBlack)
1936     hist.SetLineColor(ROOT.kBlack)
1937     hist.Draw("same")
1938     #histCOPY = hist.Clone()
1939     #histCOPY.SetXTitle("")
1940     #histCOPY.SetYTitle("")
1941
1942     #if widebins:
1943     #    histCOPY.Draw("samee1")
1944     #    histCOPY.Draw("sameaxis")
1945     #else:
1946     #    histCOPY.Draw("same")
1947     #    histCOPY.Draw("sameaxis")
1948
1949     if fitline:
1950         f.Draw("same")
1951         #f2.Draw("same")
1952         #f3.Draw("same")
1953
1954     hpeaks.SetMarkerStyle(20)
1955     hpeaks.SetMarkerSize(0.9)
1956     hpeaks.SetMarkerColor(ROOT.kRed)
1957     hpeaks.SetLineColor(ROOT.kRed)
1958     hpeaks.SetLineWidth(2)
1959     #if fitpeaks: hpeaks.Draw("same")
1960     hpeaks.Draw("same")
1961
1962     #tline1 = ROOT.TLine(hist.GetBinLowEdge(1), -window, hist.GetBinLowEdge(1), window)
1963     #tline2 = ROOT.TLine(hist.GetBinLowEdge(1), window, -hist.GetBinLowEdge(1), window)
1964     #tline3 = ROOT.TLine(-hist.GetBinLowEdge(1), window, -hist.GetBinLowEdge(1), -window)
1965     #tline4 = ROOT.TLine(-hist.GetBinLowEdge(1), -window, hist.GetBinLowEdge(1), -window)
1966     tline5 = ROOT.TLine(-hist.GetBinLowEdge(1), 0., hist.GetBinLowEdge(1), 0.)
1967     tline5.Draw()
1968     #for t in tline1, tline2, tline3, tline4, tline5: t.Draw()
1969
1970
1971 def curvatureDTsummary(tfiles, window=15., pdgSfactor=False):
1972     global h, gm2, gm1, gz, gp1, gp2, tlegend
1973
1974     set_palette("blues")
1975     phis = {-2: [], -1: [], 0: [], 1: [], 2: []}
1976     diffs = {-2: [], -1: [], 0: [], 1: [], 2: []}
1977     differrs = {-2: [], -1: [], 0: [], 1: [], 2: []}
1978     for wheelstr, wheel in ("m2", "-2"), ("m1", "-1"), ("z", "0"), ("p1", "+1"), ("p2", "+2"):
1979         for sector in "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12":
1980             curvatureplot(tfiles, "wheel%s_sector%s" % (wheelstr, sector), "deltax",
1981                     title="Wheel %s, sector %s" % (wheel, sector), fitgauss=True, reset_palette=False)
1982             if fitgauss_diff[1] < window:
1983                 uncertainty = fitgauss_diff[1]
1984                 if pdgSfactor and (fitgauss_chi2/fitgauss_ndf) > 1.: uncertainty *= sqrt(fitgauss_chi2/fitgauss_ndf)
1985
1986                 phis[int(wheel)].append(signConventions["DT", int(wheel), 1, int(sector)][4])
1987                 diffs[int(wheel)].append(fitgauss_diff[0])
1988                 differrs[int(wheel)].append(uncertainty)
1989
1990     h = ROOT.TH1F("h", "", 1, -pi, pi)
1991     h.SetAxisRange(-window, window, "Y")
1992     h.SetXTitle("#phi (rad)")
1993     h.SetYTitle("#Deltax(p_{T} #rightarrow #infty) - #Deltax(p_{T} #rightarrow 0) (mm)")
1994     h.GetXaxis().CenterTitle()
1995     h.GetYaxis().CenterTitle()
1996
1997     gm2 = ROOT.TGraphErrors(len(phis[-2]), array.array("d", phis[-2]), array.array("d", diffs[-2]),
1998             array.array("d", [0.]*len(phis[-2])), array.array("d", differrs[-2]))
1999     gm1 = ROOT.TGraphErrors(len(phis[-1]), array.array("d", phis[-1]), array.array("d", diffs[-1]),
2000             array.array("d", [0.]*len(phis[-1])), array.array("d", differrs[-1]))
2001     gz = ROOT.TGraphErrors(len(phis[0]), array.array("d", phis[0]), array.array("d", diffs[0]),
2002             array.array("d", [0.]*len(phis[0])), array.array("d", differrs[0]))
2003     gp1 = ROOT.TGraphErrors(len(phis[1]), array.array("d", phis[1]), array.array("d", diffs[1]),
2004             array.array("d", [0.]*len(phis[1])), array.array("d", differrs[1]))
2005     gp2 = ROOT.TGraphErrors(len(phis[2]), array.array("d", phis[2]), array.array("d", diffs[2]),
2006             array.array("d", [0.]*len(phis[2])), array.array("d", differrs[2]))
2007
2008     gm2.SetMarkerStyle(21); gm2.SetMarkerColor(ROOT.kRed); gm2.SetLineColor(ROOT.kRed)
2009     gm1.SetMarkerStyle(22); gm1.SetMarkerColor(ROOT.kBlue); gm1.SetLineColor(ROOT.kBlue)
2010     gz.SetMarkerStyle(3); gz.SetMarkerColor(ROOT.kBlack); gz.SetLineColor(ROOT.kBlack)
2011     gp1.SetMarkerStyle(26); gp1.SetMarkerColor(ROOT.kBlue); gp1.SetLineColor(ROOT.kBlue)
2012     gp2.SetMarkerStyle(25); gp2.SetMarkerColor(ROOT.kRed); gp2.SetLineColor(ROOT.kRed)
2013
2014     h.Draw()
2015     tlegend = ROOT.TLegend(0.25, 0.2, 0.85, 0.5)
2016     tlegend.SetFillColor(ROOT.kWhite)
2017     tlegend.SetBorderSize(0)
2018     tlegend.AddEntry(gm2, "Wheel -2", "p")
2019     tlegend.AddEntry(gm1, "Wheel -1", "p")
2020     tlegend.AddEntry(gz, "Wheel 0", "p")
2021     tlegend.AddEntry(gp1, "Wheel +1", "p")
2022     tlegend.AddEntry(gp2, "Wheel +2", "p")
2023     tlegend.Draw()
2024
2025     gm2.Draw("p")
2026     gm1.Draw("p")
2027     gz.Draw("p")
2028     gp1.Draw("p")
2029     gp2.Draw("p")
2030
2031
2032 def getname(r):
2033     if r.postal_address[0] == "DT":
2034         wheel, station, sector = r.postal_address[1:]
2035         return "DT wheel %d, station %d, sector %d" % (wheel, station, sector)
2036     elif r.postal_address[0] == "CSC":
2037         endcap, station, ring, chamber = r.postal_address[1:]
2038         if endcap != 1: station = -1 * abs(station)
2039         return "CSC ME%d/%d chamber %d" % (station, ring, chamber)
2040
2041 ddt=[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]
2042 def clearDDT():
2043     for i in range(0,15):
2044         ddt[i]=0.
2045
2046 def printDeltaTs():
2047     n = 0
2048     for t in ddt:
2049         if n==0 or n==7 or n==15: print("%d calls" % t)
2050         else: print("%d : %0.3f ms" % (n,t*1000.0))
2051         n += 1
2052
2053 def bellcurves(tfile, reports, name, twobin=True, suppressblue=False):
2054     t1 = time.time()
2055     ddt[0] += 1
2056     tdrStyle.SetOptTitle(1)
2057     tdrStyle.SetTitleBorderSize(1)
2058     tdrStyle.SetTitleFontSize(0.1)
2059     tdrStyle.SetOptStat(0)
2060     tdrStyle.SetHistMinimumZero()
2061
2062     c1.Clear()
2063     c1.ResetAttPad()
2064
2065     found = False
2066     for r in reports:
2067         if r.name == name:
2068             found = True
2069             break
2070     if not found: raise Exception("Not a valid name")
2071     if r.status == "FAIL":
2072         #raise Exception, "Fit failed"
2073         print("Fit failed")
2074         c1.Clear()
2075         return
2076
2077     Pos = "Pos"; Neg = "Neg"
2078     if not twobin:
2079         Pos = ""; Neg = ""
2080
2081     pdirPos = "MuonAlignmentFromReference/%s%s" % (name, Pos)
2082     pdirNeg = "MuonAlignmentFromReference/%s%s" % (name, Neg)
2083
2084     t2 = time.time()
2085     ddt[1] = 1./ddt[0]*((ddt[0]-1)*ddt[1] + t2-t1)
2086
2087     chamber_x = tfile.Get(pdirPos+"_x")
2088     chamber_x_fit = tfile.Get(pdirPos+"_x_fit")
2089     chamber_y = tfile.Get(pdirPos+"_y")
2090     chamber_y_fit = tfile.Get(pdirPos+"_y_fit")
2091     chamber_dxdz = tfile.Get(pdirPos+"_dxdz")
2092     chamber_dxdz_fit = tfile.Get(pdirPos+"_dxdz_fit")
2093     chamber_dydz = tfile.Get(pdirPos+"_dydz")
2094     chamber_dydz_fit = tfile.Get(pdirPos+"_dydz_fit")
2095     chamber_alphax = tfile.Get(pdirPos+"_alphax")
2096     chamber_alphax_fit = tfile.Get(pdirPos+"_alphax_fit")
2097     chamber_alphay = tfile.Get(pdirPos+"_alphay")
2098     chamber_alphay_fit = tfile.Get(pdirPos+"_alphay_fit")
2099     if twobin:
2100       chamber_x_fit2 = tfile.Get(pdirNeg+"_x_fit")
2101       chamber_y_fit2 = tfile.Get(pdirNeg+"_y_fit")
2102       chamber_dxdz_fit2 = tfile.Get(pdirNeg+"_dxdz_fit")
2103       chamber_dydz_fit2 = tfile.Get(pdirNeg+"_dydz_fit")
2104       chamber_alphax_fit2 = tfile.Get(pdirNeg+"_alphax_fit")
2105       chamber_alphay_fit2 = tfile.Get(pdirNeg+"_alphay_fit")
2106
2107     if not chamber_x:
2108         chamber_x = tfile.Get(pdirPos+"_residual")
2109         chamber_x_fit = tfile.Get(pdirPos+"_residual_fit")
2110         chamber_dxdz = tfile.Get(pdirPos+"_resslope")
2111         chamber_dxdz_fit = tfile.Get(pdirPos+"_resslope_fit")
2112         chamber_alphax = tfile.Get(pdirPos+"_alpha")
2113         chamber_alphax_fit = tfile.Get(pdirPos+"_alpha_fit")
2114         if twobin:
2115           chamber_x_fit2 = tfile.Get(pdirNeg+"_residual_fit")
2116           chamber_dxdz_fit2 = tfile.Get(pdirNeg+"_resslope_fit")
2117           chamber_alphax_fit2 = tfile.Get(pdirNeg+"_alpha_fit")
2118
2119     if not chamber_x:
2120         print("Can't find neither "+pdirPos+"_x  nor "+pdirPos+"_residual")
2121         return
2122
2123     t3 = time.time()
2124     ddt[2] = 1./ddt[0]*((ddt[0]-1)*ddt[2] + t3-t2)
2125
2126     ####
2127     chamber_x.SetAxisRange(-50., 50., "X")
2128     if chamber_x.GetRMS()>15: chamber_x.SetAxisRange(-75., 75., "X")
2129     chamber_dxdz.SetAxisRange(-30., 30., "X")
2130     chamber_alphax.SetAxisRange(-50., 50., "X")
2131     if not not chamber_y:
2132         chamber_y.SetAxisRange(-75., 75., "X")
2133         chamber_dydz.SetAxisRange(-120., 120., "X")
2134         chamber_alphay.SetAxisRange(-120., 120., "X")
2135         chamber_alphay.SetAxisRange(-75., 75., "Y")
2136     ####
2137
2138     chamber_x.SetXTitle("Local x residual (mm)")
2139     chamber_dxdz.SetXTitle("Local dx/dz residual (mrad)")
2140     chamber_alphax.SetXTitle("Local dx/dz residual (mrad)")
2141     chamber_alphax.SetYTitle("Local x residual (mm)")
2142     if not not chamber_y:
2143         chamber_y.SetXTitle("Local y residual (mm)")
2144         chamber_dydz.SetXTitle("Local dy/dz residual (mrad)")
2145         chamber_alphay.SetXTitle("Local dy/dz residual (mrad)")
2146         chamber_alphay.SetYTitle("Local y residual (mm)")
2147     if name[0:2] == "ME":
2148         chamber_x.SetXTitle("Local r#phi residual (mm)")
2149         chamber_dxdz.SetXTitle("Local d(r#phi)/dz residual (mrad)")
2150         chamber_alphax.SetXTitle("Local d(r#phi)/dz residual (mrad)")
2151         chamber_alphax.SetYTitle("Local r#phi residual (mm)")
2152
2153     t4 = time.time()
2154     ddt[3] = 1./ddt[0]*((ddt[0]-1)*ddt[3] + t4-t3)
2155
2156     for h in chamber_x, chamber_dxdz, chamber_alphax, chamber_alphax, \
2157              chamber_y, chamber_dydz, chamber_alphay, chamber_alphay:
2158         if not not h:
2159             h.GetXaxis().CenterTitle()
2160             h.GetYaxis().CenterTitle()
2161             h.GetXaxis().SetLabelSize(0.05)
2162             h.GetYaxis().SetLabelSize(0.05)
2163             h.GetXaxis().SetTitleSize(0.07)
2164             h.GetYaxis().SetTitleSize(0.07)
2165             h.GetXaxis().SetTitleOffset(0.9)
2166             h.GetYaxis().SetTitleOffset(0.9)
2167
2168     if twobin:
2169       for f in chamber_x_fit2, chamber_y_fit2, chamber_dxdz_fit2, chamber_dydz_fit2, \
2170                chamber_alphax_fit2, chamber_alphay_fit2:
2171           if not not f:
2172                f.SetLineColor(4)
2173     if not twobin:
2174         suppressblue = True
2175
2176     t5 = time.time()
2177     ddt[4] = 1./ddt[0]*((ddt[0]-1)*ddt[4] + t5-t4)
2178
2179     global l1, l2, l3, l4
2180     if not not chamber_y:
2181         c1.Clear()
2182         c1.Divide(3, 2)
2183         chamber_x.SetTitle(getname(r))
2184
2185         c1.GetPad(1).cd()
2186         chamber_x.Draw()
2187         if not suppressblue: chamber_x_fit2.Draw("same")
2188         chamber_x_fit.Draw("same")
2189         l1 = ROOT.TLatex(0.67,0.8,"#splitline{#mu: %0.2f#pm%0.2f}{#sigma: %0.1f#pm%0.1f}" % (
2190                          chamber_x_fit.GetParameter(1), chamber_x_fit.GetParError(1),
2191                          chamber_x_fit.GetParameter(2), chamber_x_fit.GetParError(2)))
2192         l1.Draw()
2193
2194         c1.GetPad(2).cd()
2195         chamber_dxdz.Draw()
2196         if not suppressblue: chamber_dxdz_fit2.Draw("same")
2197         chamber_dxdz_fit.Draw("same")
2198         l2 = ROOT.TLatex(0.67,0.8,"#splitline{#mu: %0.2f#pm%0.2f}{#sigma: %0.1f#pm%0.1f}" % (
2199                          chamber_dxdz_fit.GetParameter(1), chamber_dxdz_fit.GetParError(1),
2200                          chamber_dxdz_fit.GetParameter(2), chamber_dxdz_fit.GetParError(2)))
2201         l2.Draw()
2202
2203         c1.GetPad(3).cd()
2204         chamber_alphax.Draw("col")
2205         if not suppressblue: chamber_alphax_fit2.Draw("same")
2206         chamber_alphax_fit.Draw("same")
2207
2208         c1.GetPad(4).cd()
2209         chamber_y.Draw()
2210         if not suppressblue: chamber_y_fit2.Draw("same")
2211         chamber_y_fit.Draw("same")
2212         l3 = ROOT.TLatex(0.67,0.8,"#splitline{#mu: %0.2f#pm%0.2f}{#sigma: %0.1f#pm%0.1f}" % (
2213                          chamber_y_fit.GetParameter(1), chamber_y_fit.GetParError(1),
2214                          chamber_y_fit.GetParameter(2), chamber_y_fit.GetParError(2)))
2215         l3.Draw()
2216
2217         c1.GetPad(5).cd()
2218         chamber_dydz.Draw()
2219         if not suppressblue: chamber_dydz_fit2.Draw("same")
2220         chamber_dydz_fit.Draw("same")
2221         l4 = ROOT.TLatex(0.67,0.8,"#splitline{#mu: %0.2f#pm%0.2f}{#sigma: %0.1f#pm%0.1f}" % (
2222                          chamber_dydz_fit.GetParameter(1), chamber_dydz_fit.GetParError(1),
2223                          chamber_dydz_fit.GetParameter(2), chamber_dydz_fit.GetParError(2)))
2224         l4.Draw()
2225
2226         for lb in l1,l2,l3,l4:
2227           lb.SetNDC(1)
2228           lb.SetTextColor(ROOT.kRed)
2229
2230         c1.GetPad(6).cd()
2231         chamber_alphay.Draw("col")
2232         if not suppressblue: chamber_alphay_fit2.Draw("same")
2233         chamber_alphay_fit.Draw("same")
2234
2235     else:
2236         c1.Clear()
2237         c1.Divide(3, 1)
2238         chamber_x.SetTitle(getname(r))
2239
2240         c1.GetPad(1).cd()
2241         chamber_x.Draw()
2242         if not suppressblue: chamber_x_fit2.Draw("same")
2243         chamber_x_fit.Draw("same")
2244         l1 = ROOT.TLatex(0.67,0.8,"#splitline{#mu: %0.2f#pm%0.2f}{#sigma: %0.1f#pm%0.1f}" % (
2245                          chamber_x_fit.GetParameter(1), chamber_x_fit.GetParError(1),
2246                          chamber_x_fit.GetParameter(2), chamber_x_fit.GetParError(2)))
2247         l1.Draw()
2248
2249         c1.GetPad(2).cd()
2250         chamber_dxdz.Draw()
2251         if not suppressblue: chamber_dxdz_fit2.Draw("same")
2252         chamber_dxdz_fit.Draw("same")
2253         l2 = ROOT.TLatex(0.67,0.8,"#splitline{#mu: %0.2f#pm%0.2f}{#sigma: %0.1f#pm%0.1f}" % (
2254                          chamber_dxdz_fit.GetParameter(1), chamber_dxdz_fit.GetParError(1),
2255                          chamber_dxdz_fit.GetParameter(2), chamber_dxdz_fit.GetParError(2)))
2256         l2.Draw()
2257
2258         c1.GetPad(3).cd()
2259         chamber_alphax.Draw("col")
2260         if not suppressblue: chamber_alphax_fit2.Draw("same")
2261         chamber_alphax_fit.Draw("same")
2262
2263         for lb in l1,l2:
2264           lb.SetNDC(1)
2265           lb.SetTextColor(ROOT.kRed)
2266
2267     t6 = time.time()
2268     ddt[5] = 1./ddt[0]*((ddt[0]-1)*ddt[5] + t6-t5)
2269     ddt[6] = 1./ddt[0]*((ddt[0]-1)*ddt[6] + t6-t1)
2270
2271
2272 def polynomials(tfile, reports, name, twobin=True, suppressblue=False):
2273     t1 = time.time()
2274     ddt[7] += 1
2275     global label1, label2, label3, label4, label5, label6, label7, label8, label9
2276     plotDirectory = "MuonAlignmentFromReference"
2277     tdrStyle.SetOptTitle(1)
2278     tdrStyle.SetTitleBorderSize(1)
2279     tdrStyle.SetTitleFontSize(0.1)
2280     tdrStyle.SetOptStat(0)
2281
2282     c1.Clear()
2283     c1.ResetAttPad()
2284
2285     found = False
2286     for r in reports:
2287         if r.name == name:
2288             found = True
2289             break
2290     if not found: raise Exception("Not a valid name")
2291
2292     if r.status == "FAIL":
2293         #raise Exception, "Fit failed"
2294         print("Fit failed")
2295         c1.Clear()
2296         return
2297
2298     Pos = "Pos"; Neg = "Neg"
2299     if not twobin:
2300         Pos = ""; Neg = ""
2301
2302     pdirPos = "MuonAlignmentFromReference/%s%s" % (name, Pos)
2303     pdirNeg = "MuonAlignmentFromReference/%s%s" % (name, Neg)
2304
2305     global chamber_x_trackx, chamber_x_trackx_fit, chamber_y_trackx, chamber_y_trackx_fit, \
2306         chamber_dxdz_trackx, chamber_dxdz_trackx_fit, chamber_dydz_trackx, chamber_dydz_trackx_fit, \
2307         chamber_x_trackx_fit2, chamber_y_trackx_fit2, chamber_dxdz_trackx_fit2, chamber_dydz_trackx_fit2
2308     global chamber_x_tracky, chamber_x_tracky_fit, chamber_y_tracky, chamber_y_tracky_fit, \
2309         chamber_dxdz_tracky, chamber_dxdz_tracky_fit, chamber_dydz_tracky, chamber_dydz_tracky_fit, \
2310         chamber_x_tracky_fit2, chamber_y_tracky_fit2, chamber_dxdz_tracky_fit2, chamber_dydz_tracky_fit2
2311     global chamber_x_trackdxdz, chamber_x_trackdxdz_fit, chamber_y_trackdxdz, chamber_y_trackdxdz_fit, \
2312         chamber_dxdz_trackdxdz, chamber_dxdz_trackdxdz_fit, chamber_dydz_trackdxdz, chamber_dydz_trackdxdz_fit, \
2313         chamber_x_trackdxdz_fit2, chamber_y_trackdxdz_fit2, chamber_dxdz_trackdxdz_fit2, chamber_dydz_trackdxdz_fit2
2314     global chamber_x_trackdydz, chamber_x_trackdydz_fit, chamber_y_trackdydz, chamber_y_trackdydz_fit, \
2315         chamber_dxdz_trackdydz, chamber_dxdz_trackdydz_fit, chamber_dydz_trackdydz, chamber_dydz_trackdydz_fit, \
2316         chamber_x_trackdydz_fit2, chamber_y_trackdydz_fit2, chamber_dxdz_trackdydz_fit2, chamber_dydz_trackdydz_fit2
2317
2318     chamber_x_trackx = tfile.Get(pdirPos+"_x_trackx")
2319     chamber_x_trackx_fit = tfile.Get(pdirPos+"_x_trackx_fitline")
2320     chamber_y_trackx = tfile.Get(pdirPos+"_y_trackx")
2321     chamber_y_trackx_fit = tfile.Get(pdirPos+"_y_trackx_fitline")
2322     chamber_dxdz_trackx = tfile.Get(pdirPos+"_dxdz_trackx")
2323     chamber_dxdz_trackx_fit = tfile.Get(pdirPos+"_dxdz_trackx_fitline")
2324     chamber_dydz_trackx = tfile.Get(pdirPos+"_dydz_trackx")
2325     chamber_dydz_trackx_fit = tfile.Get(pdirPos+"_dydz_trackx_fitline")
2326     chamber_x_trackx_fit2 = tfile.Get(pdirNeg+"_x_trackx_fitline")
2327     chamber_y_trackx_fit2 = tfile.Get(pdirNeg+"_y_trackx_fitline")
2328     chamber_dxdz_trackx_fit2 = tfile.Get(pdirNeg+"_dxdz_trackx_fitline")
2329     chamber_dydz_trackx_fit2 = tfile.Get(pdirNeg+"_dydz_trackx_fitline")
2330
2331     chamber_x_tracky = tfile.Get(pdirPos+"_x_tracky")
2332     chamber_x_tracky_fit = tfile.Get(pdirPos+"_x_tracky_fitline")
2333     chamber_y_tracky = tfile.Get(pdirPos+"_y_tracky")
2334     chamber_y_tracky_fit = tfile.Get(pdirPos+"_y_tracky_fitline")
2335     chamber_dxdz_tracky = tfile.Get(pdirPos+"_dxdz_tracky")
2336     chamber_dxdz_tracky_fit = tfile.Get(pdirPos+"_dxdz_tracky_fitline")
2337     chamber_dydz_tracky = tfile.Get(pdirPos+"_dydz_tracky")
2338     chamber_dydz_tracky_fit = tfile.Get(pdirPos+"_dydz_tracky_fitline")
2339     chamber_x_tracky_fit2 = tfile.Get(pdirNeg+"_x_tracky_fitline")
2340     chamber_y_tracky_fit2 = tfile.Get(pdirNeg+"_y_tracky_fitline")
2341     chamber_dxdz_tracky_fit2 = tfile.Get(pdirNeg+"_dxdz_tracky_fitline")
2342     chamber_dydz_tracky_fit2 = tfile.Get(pdirNeg+"_dydz_tracky_fitline")
2343
2344     chamber_x_trackdxdz = tfile.Get(pdirPos+"_x_trackdxdz")
2345     chamber_x_trackdxdz_fit = tfile.Get(pdirPos+"_x_trackdxdz_fitline")
2346     chamber_y_trackdxdz = tfile.Get(pdirPos+"_y_trackdxdz")
2347     chamber_y_trackdxdz_fit = tfile.Get(pdirPos+"_y_trackdxdz_fitline")
2348     chamber_dxdz_trackdxdz = tfile.Get(pdirPos+"_dxdz_trackdxdz")
2349     chamber_dxdz_trackdxdz_fit = tfile.Get(pdirPos+"_dxdz_trackdxdz_fitline")
2350     chamber_dydz_trackdxdz = tfile.Get(pdirPos+"_dydz_trackdxdz")
2351     chamber_dydz_trackdxdz_fit = tfile.Get(pdirPos+"_dydz_trackdxdz_fitline")
2352     chamber_x_trackdxdz_fit2 = tfile.Get(pdirNeg+"_x_trackdxdz_fitline")
2353     chamber_y_trackdxdz_fit2 = tfile.Get(pdirNeg+"_y_trackdxdz_fitline")
2354     chamber_dxdz_trackdxdz_fit2 = tfile.Get(pdirNeg+"_dxdz_trackdxdz_fitline")
2355     chamber_dydz_trackdxdz_fit2 = tfile.Get(pdirNeg+"_dydz_trackdxdz_fitline")
2356
2357     chamber_x_trackdydz = tfile.Get(pdirPos+"_x_trackdydz")
2358     chamber_x_trackdydz_fit = tfile.Get(pdirPos+"_x_trackdydz_fitline")
2359     chamber_y_trackdydz = tfile.Get(pdirPos+"_y_trackdydz")
2360     chamber_y_trackdydz_fit = tfile.Get(pdirPos+"_y_trackdydz_fitline")
2361     chamber_dxdz_trackdydz = tfile.Get(pdirPos+"_dxdz_trackdydz")
2362     chamber_dxdz_trackdydz_fit = tfile.Get(pdirPos+"_dxdz_trackdydz_fitline")
2363     chamber_dydz_trackdydz = tfile.Get(pdirPos+"_dydz_trackdydz")
2364     chamber_dydz_trackdydz_fit = tfile.Get(pdirPos+"_dydz_trackdydz_fitline")
2365     chamber_x_trackdydz_fit2 = tfile.Get(pdirNeg+"_x_trackdydz_fitline")
2366     chamber_y_trackdydz_fit2 = tfile.Get(pdirNeg+"_y_trackdydz_fitline")
2367     chamber_dxdz_trackdydz_fit2 = tfile.Get(pdirNeg+"_dxdz_trackdydz_fitline")
2368     chamber_dydz_trackdydz_fit2 = tfile.Get(pdirNeg+"_dydz_trackdydz_fitline")
2369
2370     if not chamber_x_trackx:
2371         chamber_x_trackx = tfile.Get(pdirPos+"_residual_trackx")
2372         chamber_x_trackx_fit = tfile.Get(pdirPos+"_residual_trackx_fitline")
2373         chamber_dxdz_trackx = tfile.Get(pdirPos+"_resslope_trackx")
2374         chamber_dxdz_trackx_fit = tfile.Get(pdirPos+"_resslope_trackx_fitline")
2375         chamber_x_trackx_fit2 = tfile.Get(pdirNeg+"_residual_trackx_fitline")
2376         chamber_dxdz_trackx_fit2 = tfile.Get(pdirNeg+"_resslope_trackx_fitline")
2377
2378         chamber_x_tracky = tfile.Get(pdirPos+"_residual_tracky")
2379         chamber_x_tracky_fit = tfile.Get(pdirPos+"_residual_tracky_fitline")
2380         chamber_dxdz_tracky = tfile.Get(pdirPos+"_resslope_tracky")
2381         chamber_dxdz_tracky_fit = tfile.Get(pdirPos+"_resslope_tracky_fitline")
2382         chamber_x_tracky_fit2 = tfile.Get(pdirNeg+"_residual_tracky_fitline")
2383         chamber_dxdz_tracky_fit2 = tfile.Get(pdirNeg+"_resslope_tracky_fitline")
2384
2385         chamber_x_trackdxdz = tfile.Get(pdirPos+"_residual_trackdxdz")
2386         chamber_x_trackdxdz_fit = tfile.Get(pdirPos+"_residual_trackdxdz_fitline")
2387         chamber_dxdz_trackdxdz = tfile.Get(pdirPos+"_resslope_trackdxdz")
2388         chamber_dxdz_trackdxdz_fit = tfile.Get(pdirPos+"_resslope_trackdxdz_fitline")
2389         chamber_x_trackdxdz_fit2 = tfile.Get(pdirNeg+"_residual_trackdxdz_fitline")
2390         chamber_dxdz_trackdxdz_fit2 = tfile.Get(pdirNeg+"_resslope_trackdxdz_fitline")
2391
2392         chamber_x_trackdydz = tfile.Get(pdirPos+"_residual_trackdydz")
2393         chamber_x_trackdydz_fit = tfile.Get(pdirPos+"_residual_trackdydz_fitline")
2394         chamber_dxdz_trackdydz = tfile.Get(pdirPos+"_resslope_trackdydz")
2395         chamber_dxdz_trackdydz_fit = tfile.Get(pdirPos+"_resslope_trackdydz_fitline")
2396         chamber_x_trackdydz_fit2 = tfile.Get(pdirNeg+"_residual_trackdydz_fitline")
2397         chamber_dxdz_trackdydz_fit2 = tfile.Get(pdirNeg+"_resslope_trackdydz_fitline")
2398
2399     if not chamber_x_trackx:
2400         print("Can't find neither "+pdirPos+"_residual  nor "+pdirPos+"_residual_trackx")
2401         return
2402
2403     chamber_x_trackx = chamber_x_trackx.Clone()
2404     chamber_dxdz_trackx = chamber_dxdz_trackx.Clone()
2405     chamber_x_tracky = chamber_x_tracky.Clone()
2406     chamber_dxdz_tracky = chamber_dxdz_tracky.Clone()
2407     chamber_x_trackdxdz = chamber_x_trackdxdz.Clone()
2408     chamber_dxdz_trackdxdz = chamber_dxdz_trackdxdz.Clone()
2409     chamber_x_trackdydz = chamber_x_trackdydz.Clone()
2410     chamber_dxdz_trackdydz = chamber_dxdz_trackdydz.Clone()
2411
2412     if not not chamber_y_trackx:
2413         chamber_y_trackx = chamber_y_trackx.Clone()
2414         chamber_dydz_trackx = chamber_dydz_trackx.Clone()
2415         chamber_y_tracky = chamber_y_tracky.Clone()
2416         chamber_dydz_tracky = chamber_dydz_tracky.Clone()
2417         chamber_y_trackdxdz = chamber_y_trackdxdz.Clone()
2418         chamber_dydz_trackdxdz = chamber_dydz_trackdxdz.Clone()
2419         chamber_y_trackdydz = chamber_y_trackdydz.Clone()
2420         chamber_dydz_trackdydz = chamber_dydz_trackdydz.Clone()
2421
2422     if not not chamber_y_trackx:
2423         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_x_trackx")); chamber_x_trackx.Merge(tlist)
2424         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dxdz_trackx")); chamber_dxdz_trackx.Merge(tlist)
2425         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_x_tracky")); chamber_x_tracky.Merge(tlist)
2426         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dxdz_tracky")); chamber_dxdz_tracky.Merge(tlist)
2427         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_x_trackdxdz")); chamber_x_trackdxdz.Merge(tlist)
2428         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dxdz_trackdxdz")); chamber_dxdz_trackdxdz.Merge(tlist)
2429         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_x_trackdydz")); chamber_x_trackdydz.Merge(tlist)
2430         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dxdz_trackdydz")); chamber_dxdz_trackdydz.Merge(tlist)
2431         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_y_trackx")); chamber_y_trackx.Merge(tlist)
2432         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dydz_trackx")); chamber_dydz_trackx.Merge(tlist)
2433         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_y_tracky")); chamber_y_tracky.Merge(tlist)
2434         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dydz_tracky")); chamber_dydz_tracky.Merge(tlist)
2435         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_y_trackdxdz")); chamber_y_trackdxdz.Merge(tlist)
2436         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dydz_trackdxdz")); chamber_dydz_trackdxdz.Merge(tlist)
2437         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_y_trackdydz")); chamber_y_trackdydz.Merge(tlist)
2438         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_dydz_trackdydz")); chamber_dydz_trackdydz.Merge(tlist)
2439     else:
2440         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_residual_trackx")); chamber_x_trackx.Merge(tlist)
2441         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_resslope_trackx")); chamber_dxdz_trackx.Merge(tlist)
2442         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_residual_tracky")); chamber_x_tracky.Merge(tlist)
2443         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_resslope_tracky")); chamber_dxdz_tracky.Merge(tlist)
2444         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_residual_trackdxdz")); chamber_x_trackdxdz.Merge(tlist)
2445         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_resslope_trackdxdz")); chamber_dxdz_trackdxdz.Merge(tlist)
2446         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_residual_trackdydz")); chamber_x_trackdydz.Merge(tlist)
2447         tlist = ROOT.TList(); tlist.Add(tfile.Get(pdirNeg+"_resslope_trackdydz")); chamber_dxdz_trackdydz.Merge(tlist)
2448
2449     rr1=10.
2450     rr2=10.
2451     chamber_x_trackx.SetAxisRange(-rr1, rr1, "Y")
2452     chamber_dxdz_trackx.SetAxisRange(-rr2, rr2, "Y")
2453     chamber_x_tracky.SetAxisRange(-rr1, rr1, "Y")
2454     chamber_dxdz_tracky.SetAxisRange(-rr2, rr2, "Y")
2455     chamber_x_trackdxdz.SetAxisRange(-rr1, rr1, "Y")
2456     chamber_dxdz_trackdxdz.SetAxisRange(-rr2, rr2, "Y")
2457     chamber_x_trackdydz.SetAxisRange(-rr1, rr1, "Y")
2458     chamber_dxdz_trackdydz.SetAxisRange(-rr2, rr2, "Y")
2459
2460     rr3=10.
2461     if not not chamber_y_trackx:
2462         chamber_y_trackx.SetAxisRange(-rr3, rr3, "Y")
2463         chamber_dydz_trackx.SetAxisRange(-rr3, rr3, "Y")
2464         chamber_y_tracky.SetAxisRange(-rr3, rr3, "Y")
2465         chamber_dydz_tracky.SetAxisRange(-rr3, rr3, "Y")
2466         chamber_y_trackdxdz.SetAxisRange(-rr3, rr3, "Y")
2467         chamber_dydz_trackdxdz.SetAxisRange(-rr3, rr3, "Y")
2468         chamber_y_trackdydz.SetAxisRange(-rr3, rr3, "Y")
2469         chamber_dydz_trackdydz.SetAxisRange(-rr3, rr3, "Y")
2470
2471     for h in chamber_x_trackx, chamber_y_trackx, chamber_dxdz_trackx, chamber_dydz_trackx, \
2472              chamber_x_tracky, chamber_y_tracky, chamber_dxdz_tracky, chamber_dydz_tracky, \
2473              chamber_x_trackdxdz, chamber_y_trackdxdz, chamber_dxdz_trackdxdz, chamber_dydz_trackdxdz, \
2474              chamber_x_trackdydz, chamber_y_trackdydz, chamber_dxdz_trackdydz, chamber_dydz_trackdydz:
2475         if not not h:
2476             h.SetMarkerStyle(20)
2477             h.SetMarkerSize(0.5)
2478             h.GetXaxis().SetLabelSize(0.12)
2479             h.GetYaxis().SetLabelSize(0.12)
2480             h.GetXaxis().SetNdivisions(505)
2481             h.GetYaxis().SetNdivisions(505)
2482             h.GetXaxis().SetLabelOffset(0.03)
2483             h.GetYaxis().SetLabelOffset(0.03)
2484
2485     trackdxdz_minimum, trackdxdz_maximum = None, None
2486     for h in chamber_x_trackdxdz, chamber_y_trackdxdz, chamber_dxdz_trackdxdz, chamber_dydz_trackdxdz:
2487         if not not h:
2488             for i in range(1, h.GetNbinsX()+1):
2489                 if h.GetBinError(i) > 0.01 and h.GetBinContent(i) - h.GetBinError(i) < 10. and \
2490                    h.GetBinContent(i) + h.GetBinError(i) > -10.:
2491                     if not trackdxdz_minimum or trackdxdz_minimum > h.GetBinCenter(i):
2492                         trackdxdz_minimum = h.GetBinCenter(i)
2493                     if trackdxdz_maximum < h.GetBinCenter(i):
2494                         trackdxdz_maximum = h.GetBinCenter(i)
2495     if not not trackdxdz_minimum and not not trackdxdz_maximum:
2496         for h in chamber_x_trackdxdz, chamber_y_trackdxdz, chamber_dxdz_trackdxdz, chamber_dydz_trackdxdz:
2497             if not not h:
2498                 h.SetAxisRange(trackdxdz_minimum, trackdxdz_maximum, "X")
2499
2500     trackdydz_minimum, trackdydz_maximum = None, None
2501     for h in chamber_x_trackdydz, chamber_y_trackdydz, chamber_dxdz_trackdydz, chamber_dydz_trackdydz:
2502         if not not h:
2503             for i in range(1, h.GetNbinsX()+1):
2504                 if h.GetBinError(i) > 0.01 and h.GetBinContent(i) - h.GetBinError(i) < 10. and \
2505                    h.GetBinContent(i) + h.GetBinError(i) > -10.:
2506                     if not trackdydz_minimum or trackdydz_minimum > h.GetBinCenter(i):
2507                         trackdydz_minimum = h.GetBinCenter(i)
2508                     if trackdydz_maximum < h.GetBinCenter(i):
2509                         trackdydz_maximum = h.GetBinCenter(i)
2510     if not not trackdydz_minimum and not not trackdydz_maximum:
2511         for h in chamber_x_trackdydz, chamber_y_trackdydz, chamber_dxdz_trackdydz, chamber_dydz_trackdydz:
2512             if not not h:
2513                 h.SetAxisRange(trackdydz_minimum, trackdydz_maximum, "X")
2514
2515     for f in chamber_x_trackx_fit2, chamber_y_trackx_fit2, chamber_dxdz_trackx_fit2, chamber_dydz_trackx_fit2, \
2516              chamber_x_tracky_fit2, chamber_y_tracky_fit2, chamber_dxdz_tracky_fit2, chamber_dydz_tracky_fit2, \
2517              chamber_x_trackdxdz_fit2, chamber_y_trackdxdz_fit2, chamber_dxdz_trackdxdz_fit2, chamber_dydz_trackdxdz_fit2, \
2518              chamber_x_trackdydz_fit2, chamber_y_trackdydz_fit2, chamber_dxdz_trackdydz_fit2, chamber_dydz_trackdydz_fit2:
2519         if not not f:
2520             f.SetLineColor(4)
2521
2522     if not not chamber_y_trackx:
2523         c1.Clear()
2524         #c1.Divide(5, 5, 1e-5, 1e-5)
2525         pads = [None]
2526         pads.append(ROOT.TPad("p1" ,"",0.00,0.78,0.07,1.00,0,0,0))
2527         pads.append(ROOT.TPad("p2" ,"",0.07,0.78,0.34,1.00,0,0,0))
2528         pads.append(ROOT.TPad("p3" ,"",0.34,0.78,0.56,1.00,0,0,0))
2529         pads.append(ROOT.TPad("p4" ,"",0.56,0.78,0.78,1.00,0,0,0))
2530         pads.append(ROOT.TPad("p5" ,"",0.78,0.78,1.00,1.00,0,0,0))
2531         pads.append(ROOT.TPad("p6" ,"",0.00,0.56,0.07,0.78,0,0,0))
2532         pads.append(ROOT.TPad("p7" ,"",0.07,0.56,0.34,0.78,0,0,0))
2533         pads.append(ROOT.TPad("p8" ,"",0.34,0.56,0.56,0.78,0,0,0))
2534         pads.append(ROOT.TPad("p9" ,"",0.56,0.56,0.78,0.78,0,0,0))
2535         pads.append(ROOT.TPad("p10","",0.78,0.56,1.00,0.78,0,0,0))
2536         pads.append(ROOT.TPad("p11","",0.00,0.34,0.07,0.56,0,0,0))
2537         pads.append(ROOT.TPad("p12","",0.07,0.34,0.34,0.56,0,0,0))
2538         pads.append(ROOT.TPad("p13","",0.34,0.34,0.56,0.56,0,0,0))
2539         pads.append(ROOT.TPad("p14","",0.56,0.34,0.78,0.56,0,0,0))
2540         pads.append(ROOT.TPad("p15","",0.78,0.34,1.00,0.56,0,0,0))
2541         pads.append(ROOT.TPad("p16","",0.00,0.07,0.07,0.34,0,0,0))
2542         pads.append(ROOT.TPad("p17","",0.07,0.07,0.34,0.34,0,0,0))
2543         pads.append(ROOT.TPad("p18","",0.34,0.07,0.56,0.34,0,0,0))
2544         pads.append(ROOT.TPad("p19","",0.56,0.07,0.78,0.34,0,0,0))
2545         pads.append(ROOT.TPad("p20","",0.78,0.07,1.00,0.34,0,0,0))
2546         pads.append(ROOT.TPad("p21","",0.00,0.00,0.07,0.07,0,0,0))
2547         pads.append(ROOT.TPad("p22","",0.07,0.00,0.34,0.07,0,0,0))
2548         pads.append(ROOT.TPad("p23","",0.34,0.00,0.56,0.07,0,0,0))
2549         pads.append(ROOT.TPad("p24","",0.56,0.00,0.78,0.07,0,0,0))
2550         pads.append(ROOT.TPad("p25","",0.78,0.00,1.00,0.07,0,0,0))
2551         for p in pads:
2552           if not not p:
2553             p.Draw()
2554             ROOT.SetOwnership(p,False)
2555
2556         label1 = ROOT.TPaveLabel(0, 0, 1, 1, "x residuals (mm)","")
2557         label2 = ROOT.TPaveLabel(0, 0, 1, 1, "y residuals (mm)","")
2558         label3 = ROOT.TPaveLabel(0, 0, 1, 1, "dx/dz residuals (mrad)","")
2559         label4 = ROOT.TPaveLabel(0, 0, 1, 1, "dy/dz residuals (mrad)","")
2560         label5 = ROOT.TPaveLabel(0, 0, 1, 1, "x position (cm)","")
2561         label6 = ROOT.TPaveLabel(0, 0, 1, 1, "y position (cm)","")
2562         label7 = ROOT.TPaveLabel(0, 0, 1, 1, "dx/dz angle (rad)","")
2563         label8 = ROOT.TPaveLabel(0, 0, 1, 1, "dy/dz angle (rad)","")
2564         label9 = ROOT.TPaveLabel(0, 0.85, 1, 1, getname(r),"NDC")
2565
2566         for l in label1, label2, label3, label4, label5, label6, label7, label8, label9:
2567             l.SetBorderSize(0)
2568             l.SetFillColor(ROOT.kWhite)
2569
2570         for l in label1, label2, label3, label4:
2571             l.SetTextAngle(90)
2572             l.SetTextSize(0.09)
2573
2574         #label9.SetTextAngle(30)
2575         label9.SetTextSize(0.59)
2576
2577         pads[1].cd(); label1.Draw()
2578         pads[6].cd(); label2.Draw()
2579         pads[11].cd(); label3.Draw()
2580         pads[16].cd(); label4.Draw()
2581         pads[22].cd(); label5.Draw()
2582         pads[23].cd(); label6.Draw()
2583         pads[24].cd(); label7.Draw()
2584         pads[25].cd(); label8.Draw()
2585
2586         pads[2].SetRightMargin(1e-5)
2587         pads[2].SetBottomMargin(1e-5)
2588         pads[2].SetLeftMargin(0.17)
2589         pads[3].SetLeftMargin(1e-5)
2590         pads[3].SetRightMargin(1e-5)
2591         pads[3].SetBottomMargin(1e-5)
2592         pads[4].SetLeftMargin(1e-5)
2593         pads[4].SetRightMargin(1e-5)
2594         pads[4].SetBottomMargin(1e-5)
2595         pads[5].SetLeftMargin(1e-5)
2596         pads[5].SetBottomMargin(1e-5)
2597
2598         pads[7].SetRightMargin(1e-5)
2599         pads[7].SetBottomMargin(1e-5)
2600         pads[7].SetTopMargin(1e-5)
2601         pads[7].SetLeftMargin(0.17)
2602         pads[8].SetLeftMargin(1e-5)
2603         pads[8].SetRightMargin(1e-5)
2604         pads[8].SetBottomMargin(1e-5)
2605         pads[8].SetTopMargin(1e-5)
2606         pads[9].SetLeftMargin(1e-5)
2607         pads[9].SetRightMargin(1e-5)
2608         pads[9].SetBottomMargin(1e-5)
2609         pads[9].SetTopMargin(1e-5)
2610         pads[10].SetLeftMargin(1e-5)
2611         pads[10].SetBottomMargin(1e-5)
2612         pads[10].SetTopMargin(1e-5)
2613
2614         pads[12].SetRightMargin(1e-5)
2615         pads[12].SetBottomMargin(1e-5)
2616         pads[12].SetTopMargin(1e-5)
2617         pads[12].SetLeftMargin(0.17)
2618         pads[13].SetLeftMargin(1e-5)
2619         pads[13].SetRightMargin(1e-5)
2620         pads[13].SetBottomMargin(1e-5)
2621         pads[13].SetTopMargin(1e-5)
2622         pads[14].SetLeftMargin(1e-5)
2623         pads[14].SetRightMargin(1e-5)
2624         pads[14].SetBottomMargin(1e-5)
2625         pads[14].SetTopMargin(1e-5)
2626         pads[15].SetLeftMargin(1e-5)
2627         pads[15].SetBottomMargin(1e-5)
2628         pads[15].SetTopMargin(1e-5)
2629
2630         pads[17].SetRightMargin(1e-5)
2631         pads[17].SetTopMargin(1e-5)
2632         pads[17].SetLeftMargin(0.17)
2633         pads[18].SetLeftMargin(1e-5)
2634         pads[18].SetRightMargin(1e-5)
2635         pads[18].SetTopMargin(1e-5)
2636         pads[19].SetLeftMargin(1e-5)
2637         pads[19].SetRightMargin(1e-5)
2638         pads[19].SetTopMargin(1e-5)
2639         pads[20].SetLeftMargin(1e-5)
2640         pads[20].SetTopMargin(1e-5)
2641
2642         chamber_x_trackx.GetXaxis().SetLabelColor(ROOT.kWhite)
2643         chamber_x_tracky.GetXaxis().SetLabelColor(ROOT.kWhite)
2644         chamber_x_tracky.GetYaxis().SetLabelColor(ROOT.kWhite)
2645         chamber_x_trackdxdz.GetXaxis().SetLabelColor(ROOT.kWhite)
2646         chamber_x_trackdxdz.GetYaxis().SetLabelColor(ROOT.kWhite)
2647         chamber_x_trackdydz.GetXaxis().SetLabelColor(ROOT.kWhite)
2648         chamber_x_trackdydz.GetYaxis().SetLabelColor(ROOT.kWhite)
2649         chamber_y_trackx.GetXaxis().SetLabelColor(ROOT.kWhite)
2650         chamber_y_tracky.GetXaxis().SetLabelColor(ROOT.kWhite)
2651         chamber_y_tracky.GetYaxis().SetLabelColor(ROOT.kWhite)
2652         chamber_y_trackdxdz.GetXaxis().SetLabelColor(ROOT.kWhite)
2653         chamber_y_trackdxdz.GetYaxis().SetLabelColor(ROOT.kWhite)
2654         chamber_y_trackdydz.GetXaxis().SetLabelColor(ROOT.kWhite)
2655         chamber_y_trackdydz.GetYaxis().SetLabelColor(ROOT.kWhite)
2656         chamber_dxdz_trackx.GetXaxis().SetLabelColor(ROOT.kWhite)
2657         chamber_dxdz_tracky.GetXaxis().SetLabelColor(ROOT.kWhite)
2658         chamber_dxdz_tracky.GetYaxis().SetLabelColor(ROOT.kWhite)
2659         chamber_dxdz_trackdxdz.GetXaxis().SetLabelColor(ROOT.kWhite)
2660         chamber_dxdz_trackdxdz.GetYaxis().SetLabelColor(ROOT.kWhite)
2661         chamber_dxdz_trackdydz.GetXaxis().SetLabelColor(ROOT.kWhite)
2662         chamber_dxdz_trackdydz.GetYaxis().SetLabelColor(ROOT.kWhite)
2663
2664         # chamber_dydz_trackx
2665         chamber_dydz_tracky.GetYaxis().SetLabelColor(ROOT.kWhite)
2666         chamber_dydz_trackdxdz.GetYaxis().SetLabelColor(ROOT.kWhite)
2667         chamber_dydz_trackdydz.GetYaxis().SetLabelColor(ROOT.kWhite)
2668
2669         pads[2].cd()
2670         chamber_x_trackx.Draw("e1")
2671         if not suppressblue: chamber_x_trackx_fit2.Draw("samel")
2672         chamber_x_trackx_fit.Draw("samel")
2673         #label99 = ROOT.TPaveLabel(0, 0.8, 1, 1, getname(r),"NDC")
2674         print(getname(r))
2675         #label99 = ROOT.TPaveLabel(0, 0.8, 1, 1, "aaa","NDC")
2676         label9.Draw()
2677         #pads[2].Modified()
2678
2679         pads[3].cd()
2680         chamber_x_tracky.Draw("e1")
2681         if not suppressblue: chamber_x_tracky_fit2.Draw("samel")
2682         chamber_x_tracky_fit.Draw("samel")
2683
2684         pads[4].cd()
2685         chamber_x_trackdxdz.Draw("e1")
2686         if not suppressblue: chamber_x_trackdxdz_fit2.Draw("samel")
2687         chamber_x_trackdxdz_fit.Draw("samel")
2688
2689         pads[5].cd()
2690         chamber_x_trackdydz.Draw("e1")
2691         if not suppressblue: chamber_x_trackdydz_fit2.Draw("samel")
2692         chamber_x_trackdydz_fit.Draw("samel")
2693
2694         pads[7].cd()
2695         chamber_y_trackx.Draw("e1")
2696         if not suppressblue: chamber_y_trackx_fit2.Draw("samel")
2697         chamber_y_trackx_fit.Draw("samel")
2698
2699         pads[8].cd()
2700         chamber_y_tracky.Draw("e1")
2701         if not suppressblue: chamber_y_tracky_fit2.Draw("samel")
2702         chamber_y_tracky_fit.Draw("samel")
2703
2704         pads[9].cd()
2705         chamber_y_trackdxdz.Draw("e1")
2706         if not suppressblue: chamber_y_trackdxdz_fit2.Draw("samel")
2707         chamber_y_trackdxdz_fit.Draw("samel")
2708
2709         pads[10].cd()
2710         chamber_y_trackdydz.Draw("e1")
2711         if not suppressblue: chamber_y_trackdydz_fit2.Draw("samel")
2712         chamber_y_trackdydz_fit.Draw("samel")
2713
2714         pads[12].cd()
2715         chamber_dxdz_trackx.Draw("e1")
2716         if not suppressblue: chamber_dxdz_trackx_fit2.Draw("samel")
2717         chamber_dxdz_trackx_fit.Draw("samel")
2718
2719         pads[13].cd()
2720         chamber_dxdz_tracky.Draw("e1")
2721         if not suppressblue: chamber_dxdz_tracky_fit2.Draw("samel")
2722         chamber_dxdz_tracky_fit.Draw("samel")
2723
2724         pads[14].cd()
2725         chamber_dxdz_trackdxdz.Draw("e1")
2726         if not suppressblue: chamber_dxdz_trackdxdz_fit2.Draw("samel")
2727         chamber_dxdz_trackdxdz_fit.Draw("samel")
2728
2729         pads[15].cd()
2730         chamber_dxdz_trackdydz.Draw("e1")
2731         if not suppressblue: chamber_dxdz_trackdydz_fit2.Draw("samel")
2732         chamber_dxdz_trackdydz_fit.Draw("samel")
2733
2734         pads[17].cd()
2735         chamber_dydz_trackx.Draw("e1")
2736         if not suppressblue: chamber_dydz_trackx_fit2.Draw("samel")
2737         chamber_dydz_trackx_fit.Draw("samel")
2738
2739         pads[18].cd()
2740         chamber_dydz_tracky.Draw("e1")
2741         if not suppressblue: chamber_dydz_tracky_fit2.Draw("samel")
2742         chamber_dydz_tracky_fit.Draw("samel")
2743
2744         pads[19].cd()
2745         chamber_dydz_trackdxdz.Draw("e1")
2746         if not suppressblue: chamber_dydz_trackdxdz_fit2.Draw("samel")
2747         chamber_dydz_trackdxdz_fit.Draw("samel")
2748
2749         pads[20].cd()
2750         chamber_dydz_trackdydz.Draw("e1")
2751         if not suppressblue: chamber_dydz_trackdydz_fit2.Draw("samel")
2752         chamber_dydz_trackdydz_fit.Draw("samel")
2753
2754     else:
2755         c1.Clear()
2756         #c1.Divide(5, 3, 1e-5, 1e-5)
2757         pads = [None]
2758         pads.append(ROOT.TPad("p1" ,"",0.00,0.55,0.07,1.00,0,0,0))
2759         pads.append(ROOT.TPad("p2" ,"",0.07,0.55,0.34,1.00,0,0,0))
2760         pads.append(ROOT.TPad("p3" ,"",0.34,0.55,0.56,1.00,0,0,0))
2761         pads.append(ROOT.TPad("p4" ,"",0.56,0.55,0.78,1.00,0,0,0))
2762         pads.append(ROOT.TPad("p5" ,"",0.78,0.55,1.00,1.00,0,0,0))
2763         pads.append(ROOT.TPad("p6" ,"",0.00,0.1,0.07,0.55,0,0,0))
2764         pads.append(ROOT.TPad("p7" ,"",0.07,0.1,0.34,0.55,0,0,0))
2765         pads.append(ROOT.TPad("p8" ,"",0.34,0.1,0.56,0.55,0,0,0))
2766         pads.append(ROOT.TPad("p9" ,"",0.56,0.1,0.78,0.55,0,0,0))
2767         pads.append(ROOT.TPad("p10","",0.78,0.1,1.00,0.55,0,0,0))
2768         pads.append(ROOT.TPad("p11","",0.00,0.,0.07,0.1,0,0,0))
2769         pads.append(ROOT.TPad("p12","",0.07,0.,0.34,0.1,0,0,0))
2770         pads.append(ROOT.TPad("p13","",0.34,0.,0.56,0.1,0,0,0))
2771         pads.append(ROOT.TPad("p14","",0.56,0.,0.78,0.1,0,0,0))
2772         pads.append(ROOT.TPad("p15","",0.78,0.,1.00,0.1,0,0,0))
2773         for p in pads:
2774           if not not p:
2775             p.Draw()
2776             ROOT.SetOwnership(p,False)
2777
2778         label1 = ROOT.TPaveLabel(0, 0, 1, 1, "x residuals (mm)")
2779         label2 = ROOT.TPaveLabel(0, 0, 1, 1, "dx/dz residuals (mrad)")
2780         label3 = ROOT.TPaveLabel(0, 0.3, 1, 1, "x position (cm)")
2781         label4 = ROOT.TPaveLabel(0, 0.3, 1, 1, "y position (cm)")
2782         label5 = ROOT.TPaveLabel(0, 0.3, 1, 1, "dx/dz angle (rad)")
2783         label6 = ROOT.TPaveLabel(0, 0.3, 1, 1, "dy/dz angle (rad)")
2784         label9 = ROOT.TPaveLabel(0, 0.85, 1, 1, getname(r),"NDC")
2785
2786         if name[0:2] == "ME":
2787             label1 = ROOT.TPaveLabel(0, 0, 1, 1, "r#phi residuals (mm)")
2788             label2 = ROOT.TPaveLabel(0, 0, 1, 1, "d(r#phi)/dz residuals (mrad)")
2789
2790         for l in label1, label2, label3, label4, label5, label6, label9:
2791             l.SetBorderSize(0)
2792             l.SetFillColor(ROOT.kWhite)
2793
2794         for l in label1, label2:
2795             l.SetTextAngle(90)
2796             l.SetTextSize(0.09)
2797
2798         #label9.SetTextAngle(30)
2799         label9.SetTextSize(0.29)
2800
2801         pads[1].cd(); label1.Draw()
2802         pads[6].cd(); label2.Draw()
2803         pads[12].cd(); label3.Draw()
2804         pads[13].cd(); label4.Draw()
2805         pads[14].cd(); label5.Draw()
2806         pads[15].cd(); label6.Draw()
2807         #pads[11].cd(); label9.Draw()
2808
2809         pads[2].SetRightMargin(1e-5)
2810         pads[2].SetBottomMargin(1e-5)
2811         pads[3].SetLeftMargin(1e-5)
2812         pads[3].SetRightMargin(1e-5)
2813         pads[3].SetBottomMargin(1e-5)
2814         pads[4].SetLeftMargin(1e-5)
2815         pads[4].SetRightMargin(1e-5)
2816         pads[4].SetBottomMargin(1e-5)
2817         pads[5].SetLeftMargin(1e-5)
2818         pads[5].SetBottomMargin(1e-5)
2819
2820         pads[7].SetRightMargin(1e-5)
2821         pads[7].SetTopMargin(1e-5)
2822         pads[8].SetLeftMargin(1e-5)
2823         pads[8].SetRightMargin(1e-5)
2824         pads[8].SetTopMargin(1e-5)
2825         pads[9].SetLeftMargin(1e-5)
2826         pads[9].SetRightMargin(1e-5)
2827         pads[9].SetTopMargin(1e-5)
2828         pads[10].SetLeftMargin(1e-5)
2829         pads[10].SetTopMargin(1e-5)
2830
2831         chamber_x_trackx.GetXaxis().SetLabelColor(ROOT.kWhite)
2832         chamber_x_tracky.GetXaxis().SetLabelColor(ROOT.kWhite)
2833         chamber_x_tracky.GetYaxis().SetLabelColor(ROOT.kWhite)
2834         chamber_x_trackdxdz.GetXaxis().SetLabelColor(ROOT.kWhite)
2835         chamber_x_trackdxdz.GetYaxis().SetLabelColor(ROOT.kWhite)
2836         chamber_x_trackdydz.GetXaxis().SetLabelColor(ROOT.kWhite)
2837         chamber_x_trackdydz.GetYaxis().SetLabelColor(ROOT.kWhite)
2838         # chamber_dxdz_trackx
2839         chamber_dxdz_tracky.GetYaxis().SetLabelColor(ROOT.kWhite)
2840         chamber_dxdz_trackdxdz.GetYaxis().SetLabelColor(ROOT.kWhite)
2841         chamber_dxdz_trackdydz.GetYaxis().SetLabelColor(ROOT.kWhite)
2842
2843         pads[2].cd()
2844         chamber_x_trackx.Draw("e1")
2845         if not suppressblue: chamber_x_trackx_fit2.Draw("samel")
2846         chamber_x_trackx_fit.Draw("samel")
2847         label9.Draw()
2848
2849         pads[3].cd()
2850         chamber_x_tracky.Draw("e1")
2851         if not suppressblue: chamber_x_tracky_fit2.Draw("samel")
2852         chamber_x_tracky_fit.Draw("samel")
2853
2854         pads[4].cd()
2855         chamber_x_trackdxdz.Draw("e1")
2856         if not suppressblue: chamber_x_trackdxdz_fit2.Draw("samel")
2857         chamber_x_trackdxdz_fit.Draw("samel")
2858
2859         pads[5].cd()
2860         chamber_x_trackdydz.Draw("e1")
2861         if not suppressblue: chamber_x_trackdydz_fit2.Draw("samel")
2862         chamber_x_trackdydz_fit.Draw("samel")
2863
2864         pads[7].cd()
2865         chamber_dxdz_trackx.Draw("e1")
2866         if not suppressblue: chamber_dxdz_trackx_fit2.Draw("samel")
2867         chamber_dxdz_trackx_fit.Draw("samel")
2868
2869         pads[8].cd()
2870         chamber_dxdz_tracky.Draw("e1")
2871         if not suppressblue: chamber_dxdz_tracky_fit2.Draw("samel")
2872         chamber_dxdz_tracky_fit.Draw("samel")
2873
2874         pads[9].cd()
2875         chamber_dxdz_trackdxdz.Draw("e1")
2876         if not suppressblue: chamber_dxdz_trackdxdz_fit2.Draw("samel")
2877         chamber_dxdz_trackdxdz_fit.Draw("samel")
2878
2879         pads[10].cd()
2880         chamber_dxdz_trackdydz.Draw("e1")
2881         if not suppressblue: chamber_dxdz_trackdydz_fit2.Draw("samel")
2882         chamber_dxdz_trackdydz_fit.Draw("samel")
2883
2884     tn = time.time()
2885     ddt[8] = 1./ddt[7]*((ddt[7]-1)*ddt[8] + tn-t1)
2886
2887 ##################################################################################
2888
2889 def segdiff(tfiles, component, pair, **args):
2890     tdrStyle.SetOptFit(1)
2891     tdrStyle.SetOptTitle(1)
2892     tdrStyle.SetTitleBorderSize(1)
2893     tdrStyle.SetTitleFontSize(0.05)
2894     tdrStyle.SetStatW(0.2)
2895     tdrStyle.SetStatY(0.9)
2896     tdrStyle.SetStatFontSize(0.06)
2897
2898     if component[0:2] == "dt":
2899         wheel = args["wheel"]
2900         wheelletter = wheelLetter(wheel)
2901         sector = args["sector"]
2902         profname = "%s_%s_%02d_%s" % (component, wheelletter, sector, str(pair))
2903         posname = "pos" + profname
2904         negname = "neg" + profname
2905         #print profname
2906
2907         station1 = int(str(pair)[0])
2908         station2 = int(str(pair)[1])
2909         phi1 = signConventions["DT", wheel, station1, sector][4]
2910         phi2 = signConventions["DT", wheel, station2, sector][4]
2911         if abs(phi1 - phi2) > 1.:
2912             if phi1 > phi2: phi1 -= 2.*pi
2913             else: phi1 += 2.*pi
2914         phi = (phi1 + phi2) / 2.
2915         while (phi < -pi): phi += 2.*pi
2916         while (phi > pi): phi -= 2.*pi
2917
2918     elif component[0:3] == "csc":
2919         endcap = args["endcap"]
2920         if endcap=="m":
2921             endcapnum=2
2922             endcapsign="-"
2923         elif endcap=="p":
2924             endcapnum=1
2925             endcapsign="+"
2926         else: raise Exception
2927
2928         ring = args["ring"]
2929         if ring>2 or ring<1: raise Exception
2930         station1 = int(str(pair)[0])
2931         station2 = int(str(pair)[1])
2932         if   ring==1: ringname="inner"
2933         elif ring==2: ringname="outer"
2934         else: raise Exception
2935
2936         chamber = args["chamber"]
2937         if (ring==1 and chamber>18) or (ring==2 and chamber>36): raise Exception
2938
2939         profname = "csc%s_%s_%s_%02d_%s" % (ringname,component[4:], endcap, chamber, str(pair))
2940         posname = "pos" + profname
2941         negname = "neg" + profname
2942         #print profname
2943
2944         station1 = int(str(pair)[0])
2945         station2 = int(str(pair)[1])
2946         phi1 = signConventions["CSC", endcapnum, station1, ring, chamber][4]
2947         phi2 = signConventions["CSC", endcapnum, station2, ring, chamber][4]
2948         if abs(phi1 - phi2) > 1.:
2949             if phi1 > phi2: phi1 -= 2.*pi
2950             else: phi1 += 2.*pi
2951         phi = (phi1 + phi2) / 2.
2952         while (phi < -pi*5./180.): phi += 2.*pi
2953         while (phi > pi*(2.-5./180.)): phi -= 2.*pi
2954
2955     else: raise Exception
2956
2957     if "window" in args: window = args["window"]
2958     else: window = 5.
2959
2960     global tmpprof, tmppos, tmpneg
2961     pdir = "AlignmentMonitorSegmentDifferences/iter1/"
2962     tmpprof = tfiles[0].Get(pdir + profname).Clone()
2963     tmpprof.SetMarkerStyle(8)
2964     tmppos = tfiles[0].Get(pdir + posname).Clone()
2965     tmpneg = tfiles[0].Get(pdir + negname).Clone()
2966     for tfile in tfiles[1:]:
2967         tmpprof.Add(tfile.Get(pdir + profname))
2968         tmppos.Add(tfile.Get(pdir + posname))
2969         tmpneg.Add(tfile.Get(pdir + negname))
2970
2971     for i in range(1, tmpprof.GetNbinsX()+1):
2972         if tmpprof.GetBinError(i) < 1e-5:
2973             tmpprof.SetBinError(i, 100.)
2974     tmpprof.SetAxisRange(-window, window, "Y")
2975
2976     f = ROOT.TF1("p1", "[0] + [1]*x", tmpprof.GetBinLowEdge(1), -tmpprof.GetBinLowEdge(1))
2977     f.SetParameters((tmppos.GetMean() + tmpneg.GetMean())/2., 0.)
2978
2979     tmpprof.SetXTitle("q/p_{T} (c/GeV)")
2980     if component == "dt13_resid":
2981         tmpprof.SetYTitle("#Deltax^{local} (mm)")
2982         tmppos.SetXTitle("#Deltax^{local} (mm)")
2983         tmpneg.SetXTitle("#Deltax^{local} (mm)")
2984         f.SetParNames("#Deltax^{local}_{0}", "Slope")
2985     if component == "dt13_slope":
2986         tmpprof.SetYTitle("#Deltadx/dz^{local} (mrad)")
2987         tmppos.SetXTitle("#Deltadx/dz^{local} (mrad)")
2988         tmpneg.SetXTitle("#Deltadx/dz^{local} (mrad)")
2989         f.SetParNames("#Deltadx/dz^{local}_{0}", "Slope")
2990     if component == "dt2_resid":
2991         tmpprof.SetYTitle("#Deltay^{local} (mm)")
2992         tmppos.SetXTitle("#Deltay^{local} (mm)")
2993         tmpneg.SetXTitle("#Deltay^{local} (mm)")
2994         f.SetParNames("#Deltay^{local}_{0}", "Slope")
2995     if component == "dt2_slope":
2996         tmpprof.SetYTitle("#Deltady/dz^{local} (mrad)")
2997         tmppos.SetXTitle("#Deltady/dz^{local} (mrad)")
2998         tmpneg.SetXTitle("#Deltady/dz^{local} (mrad)")
2999         f.SetParNames("#Deltady/dz^{local}_{0}", "Slope")
3000     if component == "csc_resid":
3001         tmpprof.SetXTitle("q/p_{z} (c/GeV)")
3002         tmpprof.SetYTitle("#Delta(r#phi)^{local} (mm)")
3003         tmppos.SetXTitle("#Delta(r#phi)^{local} (mm)")
3004         tmpneg.SetXTitle("#Delta(r#phi)^{local} (mm)")
3005         f.SetParNames("#Delta(r#phi)^{local}_{0}", "Slope")
3006     if component == "csc_slope":
3007         tmpprof.SetXTitle("q/p_{z} (c/GeV)")
3008         tmpprof.SetYTitle("#Deltad(r#phi)/dz^{local} (mrad)")
3009         tmppos.SetXTitle("#Deltad(r#phi)/dz^{local} (mrad)")
3010         tmpneg.SetXTitle("#Deltad(r#phi)/dz^{local} (mrad)")
3011         f.SetParNames("#Deltad(r#phi)/dz^{local}_{0}", "Slope")
3012
3013     tmpprof.GetXaxis().CenterTitle()
3014     tmpprof.GetYaxis().CenterTitle()
3015     tmppos.GetXaxis().CenterTitle()
3016     tmpneg.GetXaxis().CenterTitle()
3017     if component[0:2] == "dt":
3018         tmpprof.SetTitle("MB%d - MB%d, wheel %d, sector %02d" % (station1, station2, int(wheel), int(sector)))
3019     elif component[0:3] == "csc":
3020         tmpprof.SetTitle("ME%d - ME%d, for ME%s%d/%d/%d" % (station1, station2, endcapsign, station2, ring, chamber))
3021     else: raise Exception
3022
3023     tmppos.SetTitle("Positive muons")
3024     tmpneg.SetTitle("Negative muons")
3025
3026     c1.Clear()
3027     c1.Divide(2, 1)
3028     c1.GetPad(1).cd()
3029     fit1 = tmpprof.Fit("p1", "qS")
3030     tmpprof.Draw("e1")
3031     c1.GetPad(2).cd()
3032     c1.GetPad(2).Divide(1, 2)
3033     c1.GetPad(2).GetPad(1).cd()
3034     tmppos.Draw()
3035     f = ROOT.TF1("gausR", "[0]*exp(-(x - [1])**2 / 2. / [2]**2) / sqrt(2.*3.1415926) / [2]",
3036                  tmppos.GetMean() - tmppos.GetRMS(), tmppos.GetMean() + tmppos.GetRMS())
3037     f.SetParameters(tmppos.GetEntries() * ((10. - -10.)/100.), tmppos.GetMean(), tmppos.GetRMS())
3038     f.SetParNames("Constant", "Mean", "Sigma")
3039     fit2 = tmppos.Fit("gausR", "qRS")
3040     c1.GetPad(2).GetPad(2).cd()
3041     tmpneg.Draw()
3042     f = ROOT.TF1("gausR", "[0]*exp(-(x - [1])**2 / 2. / [2]**2) / sqrt(2.*3.1415926) / [2]",
3043                  tmpneg.GetMean() - tmpneg.GetRMS(), tmpneg.GetMean() + tmpneg.GetRMS())
3044     f.SetParameters(tmpneg.GetEntries() * ((10. - -10.)/100.), tmpneg.GetMean(), tmpneg.GetRMS())
3045     f.SetParNames("Constant", "Mean", "Sigma")
3046     fit3 = tmpneg.Fit("gausR", "qRS")
3047
3048     fit1ok = fit1.Status()==0 and fit1.CovMatrixStatus()==3
3049     fit2ok = fit2.Status()==0 and fit2.CovMatrixStatus()==3
3050     fit3ok = fit3.Status()==0 and fit3.CovMatrixStatus()==3
3051
3052     fitresult1 = None, None
3053     if fit1ok:
3054         fitresult1 = tmpprof.GetFunction("p1").GetParameter(0), tmpprof.GetFunction("p1").GetParError(0)
3055     fitresult2 = None, None
3056     if fit2ok and fit3ok:
3057         fitresult2 = (tmppos.GetFunction("gausR").GetParameter(1) + tmpneg.GetFunction("gausR").GetParameter(1)) / 2., \
3058                      sqrt(tmppos.GetFunction("gausR").GetParError(1)**2 + tmpneg.GetFunction("gausR").GetParError(1)**2) / 2.
3059     return phi, fitresult1[0], fitresult1[1], fitresult2[0], fitresult2[1], fit1ok, fit2ok, fit3ok
3060
3061
3062
3063 ##################################################################################
3064
3065 def segdiff_xalign(tfiles, component, **args):
3066     tdrStyle.SetOptFit(1)
3067     tdrStyle.SetOptTitle(1)
3068     tdrStyle.SetTitleBorderSize(1)
3069     tdrStyle.SetTitleFontSize(0.05)
3070     tdrStyle.SetStatW(0.2)
3071     tdrStyle.SetStatY(0.9)
3072     tdrStyle.SetStatFontSize(0.06)
3073
3074     if component[0:4] == "x_dt":
3075         wheel = int(args["wheel"])
3076         if int(wheel)<0:
3077           wheell = "m%d" % abs(wheel)
3078           endcapsign="-"
3079         else:
3080           wheell = "p%d" % abs(wheel)
3081           endcapsign="+"
3082         station_dt = component[4]
3083         station_csc_1 = args["cscstations"][0]
3084         if station_csc_1=='1':  ring_1 = 3
3085         else:                   ring_1 = 2
3086         sector = args["sector"]
3087         profname = "%s%s_W%sS%02d" % (component, station_csc_1, wheell, sector)
3088         posname_1 = "pos_" + profname
3089         negname_1 = "neg_" + profname
3090         if len(args["cscstations"]) > 1:
3091           station_csc_2 = args["cscstations"][1]
3092           if station_csc_2=='1':  ring_2 = 3
3093           else:                   ring_2 = 2
3094           profname = "%s%s_W%sS%02d" % (component, station_csc_2, wheell, sector)
3095           posname_2 = "pos_" + profname
3096           negname_2 = "neg_" + profname
3097
3098         phi = signConventions["DT", wheel, int(station_dt), sector][4]
3099         while (phi < -pi): phi += 2.*pi
3100         while (phi > pi): phi -= 2.*pi
3101
3102     else: raise Exception
3103
3104     if "window" in args: window = args["window"]
3105     else: window = 5.
3106
3107     global tmppos, tmpneg, tmppos_2, tmpneg_2
3108     pdir = "AlignmentMonitorSegmentDifferences/iter1/"
3109     tmppos = tfiles[0].Get(pdir + posname_1).Clone()
3110     tmpneg = tfiles[0].Get(pdir + negname_1).Clone()
3111     if len(args["cscstations"]) > 1:
3112       tmppos_2 = tfiles[0].Get(pdir + posname_2).Clone()
3113       tmpneg_2 = tfiles[0].Get(pdir + negname_2).Clone()
3114     tmpneg.Rebin(2); tmppos.Rebin(2)
3115     for tfile in tfiles[1:]:
3116       tmppos.Add(tfile.Get(pdir + posname_1))
3117       tmpneg.Add(tfile.Get(pdir + negname_1))
3118       if len(args["cscstations"]) > 1:
3119         tmppos_2.Add(tfile.Get(pdir + posname_2))
3120         tmpneg_2.Add(tfile.Get(pdir + negname_2))
3121       tmpneg_2.Rebin(2); tmppos_2.Rebin(2)
3122
3123     result = {}
3124     result["fit_ok"] = False
3125     result["phi"] = phi
3126     ntot = tmppos.GetEntries() + tmpneg.GetEntries()
3127     if ntot == 0:
3128       return result
3129
3130     tmppos.SetXTitle("#Deltax^{loc}_{MB} - r_{DT}/r_{CSC}#times#Deltax^{loc}_{ME} (mm)")
3131     tmpneg.SetXTitle("#Deltax^{loc}_{MB} - r_{DT}/r_{CSC}#times#Deltax^{loc}_{ME} (mm)")
3132     title1 =  "MB(W%+d St%s Sec%d) - ME%s%s/%d" % (wheel, station_dt, int(sector), endcapsign, station_csc_1, ring_1)
3133     tmppos.SetTitle("Positive #mu:  %s" % title1);
3134     tmpneg.SetTitle("Negative #mu:  %s" % title1);
3135     tmppos.GetXaxis().CenterTitle()
3136     tmpneg.GetXaxis().CenterTitle()
3137     if len(args["cscstations"]) > 1:
3138       tmppos.SetXTitle("#Deltax^{loc}_{DT} - r_{DT}/r_{CSC}#times#Deltax^{loc}_{CSC} (mm)")
3139       tmpneg.SetXTitle("#Deltax^{loc}_{DT} - r_{DT}/r_{CSC}#times#Deltax^{loc}_{CSC} (mm)")
3140       title2 =  "MB(W%+d St%s Sec%d) - ME%s%s/%d" % (wheel, station_dt, int(sector), endcapsign, station_csc_2, ring_2)
3141       tmppos_2.SetTitle("Positive #mu:  %s" % title2);
3142       tmpneg_2.SetTitle("Negative #mu:  %s" % title2);
3143       tmppos_2.GetXaxis().CenterTitle()
3144       tmpneg_2.GetXaxis().CenterTitle()
3145
3146     c1.Clear()
3147     c1.Divide(2, 2)
3148
3149     c1.GetPad(1).cd()
3150     tmppos.Draw()
3151     fpos = ROOT.TF1("gpos", "gaus", tmppos.GetMean() - tmppos.GetRMS(), tmppos.GetMean() + tmppos.GetRMS())
3152     fpos.SetParameters(tmppos.GetEntries() * 2.5, tmppos.GetMean(), tmppos.GetRMS())
3153     fit_pos = tmppos.Fit("gpos", "qRS")
3154
3155     c1.GetPad(3).cd()
3156     tmpneg.Draw()
3157     fneg = ROOT.TF1("gneg", "gaus", tmpneg.GetMean() - tmpneg.GetRMS(), tmpneg.GetMean() + tmpneg.GetRMS())
3158     fneg.SetParameters(tmpneg.GetEntries() * 2.5, tmpneg.GetMean(), tmpneg.GetRMS())
3159     fit_neg = tmpneg.Fit("gneg", "qRS")
3160
3161     result["fit_ok"] = (fit_pos.Status()==0 and fit_pos.CovMatrixStatus()==3 and fit_neg.Status()==0 and fit_neg.CovMatrixStatus()==3)
3162     result["fit_peak"] = (fpos.GetParameter(1)*tmppos.GetEntries() + fneg.GetParameter(1)*tmpneg.GetEntries()) / ntot
3163     result["fit_peak_error"] = sqrt( (fpos.GetParError(1)*tmppos.GetEntries())**2 + (fneg.GetParError(1)*tmpneg.GetEntries())**2) / ntot
3164
3165     if len(args["cscstations"]) > 1:
3166       c1.GetPad(2).cd()
3167       tmppos_2.Draw()
3168       fpos_2 = ROOT.TF1("gpos2", "gaus", tmppos_2.GetMean() - tmppos_2.GetRMS(), tmppos_2.GetMean() + tmppos_2.GetRMS())
3169       fpos_2.SetParameters(tmppos_2.GetEntries() * 2.5, tmppos_2.GetMean(), tmppos_2.GetRMS())
3170       fit_pos_2 = tmppos_2.Fit("gpos2", "qRS")
3171
3172       c1.GetPad(4).cd()
3173       tmpneg_2.Draw()
3174       fneg_2 = ROOT.TF1("gneg2", "gaus", tmpneg_2.GetMean() - tmpneg_2.GetRMS(), tmpneg_2.GetMean() + tmpneg_2.GetRMS())
3175       fneg_2.SetParameters(tmpneg_2.GetEntries() * 2.5, tmpneg_2.GetMean(), tmpneg_2.GetRMS())
3176       fit_neg_2 = tmpneg_2.Fit("gneg2", "qRS")
3177
3178       result["fit_ok_2"] = (fit_pos_2.Status()==0 and fit_pos_2.CovMatrixStatus()==3 and fit_neg_2.Status()==0 and fit_neg_2.CovMatrixStatus()==3)
3179       ntot = tmppos_2.GetEntries() + tmpneg_2.GetEntries()
3180       result["fit_peak_2"] = (fpos_2.GetParameter(1)*tmppos_2.GetEntries() + fneg_2.GetParameter(1)*tmpneg_2.GetEntries()) / ntot
3181       result["fit_peak_error_2"] = sqrt( (fpos_2.GetParError(1)*tmppos_2.GetEntries())**2 + (fneg_2.GetParError(1)*tmpneg_2.GetEntries())**2) / ntot
3182
3183     return result
3184
3185 ##################################################################################
3186
3187 def segdiffvsphi_xalign(tfiles, wheel, window=10.):
3188     tdrStyle.SetOptTitle(1)
3189     tdrStyle.SetTitleBorderSize(1)
3190     tdrStyle.SetTitleFontSize(0.05)
3191
3192     global htemp, gtemp_12, gtemp_21, gtemp_11, tlegend
3193     htemp = ROOT.TH1F("htemp", "", 1, -pi, pi)
3194     gtemp_11_phi, gtemp_11_val, gtemp_11_err = [], [], []
3195     gtemp_12_phi, gtemp_12_val, gtemp_12_err = [], [], []
3196     gtemp_21_phi, gtemp_21_val, gtemp_21_err = [], [], []
3197     for sector in range(1, 12+1):
3198       #print "sect", sector
3199       r1 = segdiff_xalign(tfiles, "x_dt1_csc", wheel=wheel, sector=sector, cscstations = "12")
3200       r2 = segdiff_xalign(tfiles, "x_dt2_csc", wheel=wheel, sector=sector, cscstations = "1")
3201
3202       if r1["fit_ok"]:
3203         gtemp_11_phi.append(r1["phi"])
3204         gtemp_11_val.append(r1["fit_peak"])
3205         gtemp_11_err.append(r1["fit_peak_error"])
3206       if r1["fit_ok_2"]:
3207         gtemp_12_phi.append(r1["phi"])
3208         gtemp_12_val.append(r1["fit_peak_2"])
3209         gtemp_12_err.append(r1["fit_peak_error_2"])
3210       if r2["fit_ok"]:
3211         gtemp_21_phi.append(r2["phi"])
3212         gtemp_21_val.append(r2["fit_peak"])
3213         gtemp_21_err.append(r2["fit_peak_error"])
3214
3215     #print "len(gtemp_11_phi) ",len(gtemp_11_phi)
3216     #print "len(gtemp_12_phi) ",len(gtemp_12_phi)
3217     #print "len(gtemp_21_phi) ",len(gtemp_21_phi)
3218     if len(gtemp_11_phi) > 0:
3219         gtemp_11 = ROOT.TGraphErrors(len(gtemp_11_phi), array.array("d", gtemp_11_phi), array.array("d", gtemp_11_val),
3220                                      array.array("d", [0.] * len(gtemp_11_phi)), array.array("d", gtemp_11_err))
3221     if len(gtemp_12_phi) > 0:
3222         gtemp_12 = ROOT.TGraphErrors(len(gtemp_12_phi), array.array("d", gtemp_12_phi), array.array("d", gtemp_12_val),
3223                                      array.array("d", [0.] * len(gtemp_12_phi)), array.array("d", gtemp_12_err))
3224     if len(gtemp_11_phi) > 0:
3225         gtemp_21 = ROOT.TGraphErrors(len(gtemp_21_phi), array.array("d", gtemp_21_phi), array.array("d", gtemp_21_val),
3226                                      array.array("d", [0.] * len(gtemp_21_phi)), array.array("d", gtemp_21_err))
3227
3228     if len(gtemp_11_phi) > 0:
3229         gtemp_11.SetMarkerStyle(20);  gtemp_11.SetMarkerSize(1.5);
3230         gtemp_11.SetMarkerColor(ROOT.kRed);  gtemp_11.SetLineColor(ROOT.kRed)
3231     if len(gtemp_12_phi) > 0:
3232         gtemp_12.SetMarkerStyle(22);  gtemp_12.SetMarkerSize(1.);
3233         gtemp_12.SetMarkerColor(ROOT.kGreen+2);  gtemp_12.SetLineColor(ROOT.kGreen+2)
3234     if len(gtemp_21_phi) > 0:
3235         gtemp_21.SetMarkerStyle(21);  gtemp_21.SetMarkerSize(1.5);
3236         gtemp_21.SetMarkerColor(ROOT.kBlue);  gtemp_21.SetLineColor(ROOT.kBlue)
3237
3238     htemp.SetTitle("Wheel %+d" % wheel)
3239     htemp.SetAxisRange(-window, window, "Y")
3240     htemp.SetXTitle("#phi of MB")
3241     htemp.SetYTitle("#Deltax^{loc}_{DT} - r_{DT}/r_{CSC}#times#Deltax^{loc}_{CSC} (mm)")
3242     htemp.GetXaxis().CenterTitle()
3243     htemp.GetYaxis().CenterTitle()
3244     htemp.GetYaxis().SetTitleOffset(0.75)
3245
3246     c1.Clear()
3247     htemp.Draw()
3248     if len(gtemp_12_phi) > 0:
3249         gtemp_12.Draw("p")
3250     if len(gtemp_21_phi) > 0:
3251         gtemp_21.Draw("p")
3252     if len(gtemp_11_phi) > 0:
3253         gtemp_11.Draw("p")
3254
3255     tlegend = ROOT.TLegend(0.59, 0.75, 0.99, 0.92)
3256     tlegend.SetBorderSize(0)
3257     tlegend.SetFillColor(ROOT.kWhite)
3258     if len(gtemp_11_phi) > 0:
3259         tlegend.AddEntry(gtemp_11, "MB1 - ME1/3 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_11_val), stdev(gtemp_11_val)), "pl")
3260     if len(gtemp_21_phi) > 0:
3261         tlegend.AddEntry(gtemp_21, "MB2 - ME1/3 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_21_val), stdev(gtemp_21_val)), "pl")
3262     if len(gtemp_12_phi) > 0:
3263         tlegend.AddEntry(gtemp_12, "MB1 - ME2/2 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_12_val), stdev(gtemp_12_val)), "pl")
3264     #if len(gtemp_12_phi) > 0:
3265     #    tlegend.AddEntry(gtemp_12, "total mean: %4.2f, total RMS: %4.2f" % \
3266     #                               (mean(gtemp_11_val + gtemp_12_val + gtemp_21_val),
3267     #                               stdev(gtemp_11_val + gtemp_12_val + gtemp_21_val)), "")
3268     tlegend.Draw()
3269
3270     f_11 = ROOT.TF1("f11", "[0] + [1]*sin(x) + [2]*cos(x)", -pi, pi)
3271     f_11.SetLineColor(ROOT.kRed)
3272     f_11.SetLineWidth(2)
3273     f_21 = ROOT.TF1("f21", "[0] + [1]*sin(x) + [2]*cos(x)", -pi, pi)
3274     f_21.SetLineColor(ROOT.kBlue)
3275     f_21.SetLineWidth(2)
3276     if len(gtemp_11_phi) > 0:
3277       gtemp_11.Fit(f_11,"")
3278     if len(gtemp_21_phi) > 0:
3279       gtemp_21.Fit(f_21,"")
3280
3281     global f_txt,f_11_txt, f_21_txt
3282     f_txt = ROOT.TLatex(-2.9, -0.7*window, "ME1/3 ring corrections equivalents:")
3283     f_txt.SetTextSize(0.028)
3284     f_txt.Draw()
3285     if len(gtemp_11_phi) > 0:
3286       rdt = signConventions[("DT", 2, 1, 1)][3]*10
3287       f_11_txt = ROOT.TLatex(-2.9, -0.8*window, "#Deltax=%.2f#pm%.2f mm   #Deltay=%.2f#pm%.2f mm   #Delta#phi_{z}=%.2f#pm%.2f mrad" % (
3288                              -f_11.GetParameter(1), f_11.GetParError(1), f_11.GetParameter(2), f_11.GetParError(2), -f_11.GetParameter(0)/rdt*1000, f_11.GetParError(0)/rdt*1000))
3289       f_11_txt.SetTextSize(0.028)
3290       f_11_txt.SetTextColor(ROOT.kRed)
3291       f_11_txt.Draw()
3292     if len(gtemp_11_phi) > 0:
3293       rdt = signConventions[("DT", 2, 2, 1)][3]*10
3294       f_21_txt = ROOT.TLatex(-2.9, -0.9*window, "#Deltax=%.2f#pm%.2f mm   #Deltay=%.2f#pm%.2f mm   #Delta#phi_{z}=%.2f#pm%.2f mrad" % (
3295                              -f_21.GetParameter(1), f_21.GetParError(1), f_21.GetParameter(2), f_21.GetParError(2), -f_21.GetParameter(0)/rdt*1000, f_21.GetParError(0)/rdt*1000))
3296       f_21_txt.SetTextSize(0.028)
3297       f_21_txt.SetTextColor(ROOT.kBlue)
3298       f_21_txt.Draw()
3299
3300 ##################################################################################
3301
3302 def segdiffvsphi(tfiles, reports, component, wheel, window=5., excludesectors=()):
3303     tdrStyle.SetOptTitle(1)
3304     tdrStyle.SetTitleBorderSize(1)
3305     tdrStyle.SetTitleFontSize(0.05)
3306
3307     global htemp, gtemp_12, gtemp2_12, gtemp_23, gtemp2_23, gtemp_34, gtemp2_34, tlegend
3308     htemp = ROOT.TH1F("htemp", "", 1, -pi, pi)
3309     gtemp_12_phi, gtemp_12_val, gtemp_12_err, gtemp_12_val2, gtemp_12_err2 = [], [], [], [], []
3310     gtemp_23_phi, gtemp_23_val, gtemp_23_err, gtemp_23_val2, gtemp_23_err2 = [], [], [], [], []
3311     gtemp_34_phi, gtemp_34_val, gtemp_34_err, gtemp_34_val2, gtemp_34_err2 = [], [], [], [], []
3312     for sector in range(1, 12+1):
3313         #print "sect", sector
3314         r1_found, r2_found, r3_found, r4_found = False, False, False, False
3315         for r1 in reports:
3316             if r1.postal_address == ("DT", wheel, 1, sector):
3317                 r1_found = True
3318                 break
3319         for r2 in reports:
3320             if r2.postal_address == ("DT", wheel, 2, sector):
3321                 r2_found = True
3322                 break
3323         for r3 in reports:
3324             if r3.postal_address == ("DT", wheel, 3, sector):
3325                 r3_found = True
3326                 break
3327         for r4 in reports:
3328             if r4.postal_address == ("DT", wheel, 4, sector):
3329                 r4_found = True
3330                 break
3331         #print "rfounds: ", r1_found, r2_found, r3_found, r4_found
3332
3333         if sector not in excludesectors:
3334             if r1_found and r2_found and r1.status == "PASS" and r2.status == "PASS":
3335                 phi, val, err, val2, err2, fit1, fit2, fit3 = segdiff(tfiles, component, 12, wheel=wheel, sector=sector)
3336                 #print "segdif 12", phi, val, err, val2, err2, fit1, fit2, fit3
3337                 if fit1 and fit2 and fit3:
3338                     gtemp_12_phi.append(phi)
3339                     gtemp_12_val.append(val)
3340                     gtemp_12_err.append(err)
3341                     gtemp_12_val2.append(val2)
3342                     gtemp_12_err2.append(err2)
3343             if r2_found and r3_found and r2.status == "PASS" and r3.status == "PASS":
3344                 phi, val, err, val2, err2, fit1, fit2, fit3 = segdiff(tfiles, component, 23, wheel=wheel, sector=sector)
3345                 #print "segdif 23", phi, val, err, val2, err2, fit1, fit2, fit3
3346                 if fit1 and fit2 and fit3:
3347                     gtemp_23_phi.append(phi)
3348                     gtemp_23_val.append(val)
3349                     gtemp_23_err.append(err)
3350                     gtemp_23_val2.append(val2)
3351                     gtemp_23_err2.append(err2)
3352             if component[:4] == "dt13":
3353                 if r3_found and r4_found and r3.status == "PASS" and r4.status == "PASS":
3354                     phi, val, err, val2, err2, fit1, fit2, fit3 = segdiff(tfiles, component, 34, wheel=wheel, sector=sector)
3355                     #print "segdif 34", phi, val, err, val2, err2, fit1, fit2, fit3
3356                     if fit1 and fit2 and fit3:
3357                         gtemp_34_phi.append(phi)
3358                         gtemp_34_val.append(val)
3359                         gtemp_34_err.append(err)
3360                         gtemp_34_val2.append(val2)
3361                         gtemp_34_err2.append(err2)
3362
3363     #print "len(gtemp_12_phi) ", len(gtemp_12_phi)
3364     #print "len(gtemp_23_phi) ",len(gtemp_23_phi)
3365     #print "len(gtemp_34_phi) ",len(gtemp_34_phi)
3366     if len(gtemp_12_phi) > 0:
3367         gtemp_12 = ROOT.TGraphErrors(len(gtemp_12_phi), array.array("d", gtemp_12_phi), array.array("d", gtemp_12_val),
3368                                      array.array("d", [0.] * len(gtemp_12_phi)), array.array("d", gtemp_12_err))
3369         gtemp2_12 = ROOT.TGraphErrors(len(gtemp_12_phi), array.array("d", gtemp_12_phi), array.array("d", gtemp_12_val2),
3370                                       array.array("d", [0.] * len(gtemp_12_phi)), array.array("d", gtemp_12_err2))
3371     if len(gtemp_23_phi) > 0:
3372         gtemp_23 = ROOT.TGraphErrors(len(gtemp_23_phi), array.array("d", gtemp_23_phi), array.array("d", gtemp_23_val),
3373                                      array.array("d", [0.] * len(gtemp_23_phi)), array.array("d", gtemp_23_err))
3374         gtemp2_23 = ROOT.TGraphErrors(len(gtemp_23_phi), array.array("d", gtemp_23_phi), array.array("d", gtemp_23_val2),
3375                                       array.array("d", [0.] * len(gtemp_23_phi)), array.array("d", gtemp_23_err2))
3376     if len(gtemp_34_phi) > 0:
3377         gtemp_34 = ROOT.TGraphErrors(len(gtemp_34_phi), array.array("d", gtemp_34_phi), array.array("d", gtemp_34_val),
3378                                      array.array("d", [0.] * len(gtemp_34_phi)), array.array("d", gtemp_34_err))
3379         gtemp2_34 = ROOT.TGraphErrors(len(gtemp_34_phi), array.array("d", gtemp_34_phi), array.array("d", gtemp_34_val2),
3380                                       array.array("d", [0.] * len(gtemp_34_phi)), array.array("d", gtemp_34_err2))
3381
3382     if len(gtemp_12_phi) > 0:
3383         gtemp_12.SetMarkerStyle(20);  gtemp_12.SetMarkerSize(1.);
3384         gtemp_12.SetMarkerColor(ROOT.kBlue);  gtemp_12.SetLineColor(ROOT.kBlue)
3385         gtemp2_12.SetMarkerStyle(24); gtemp2_12.SetMarkerSize(1.);
3386         gtemp2_12.SetMarkerColor(ROOT.kBlue); gtemp2_12.SetLineColor(ROOT.kBlue)
3387     if len(gtemp_23_phi) > 0:
3388         gtemp_23.SetMarkerStyle(21);  gtemp_23.SetMarkerSize(1.);
3389         gtemp_23.SetMarkerColor(ROOT.kRed);   gtemp_23.SetLineColor(ROOT.kRed)
3390         gtemp2_23.SetMarkerStyle(25); gtemp2_23.SetMarkerSize(1.);
3391         gtemp2_23.SetMarkerColor(ROOT.kRed);  gtemp2_23.SetLineColor(ROOT.kRed)
3392     if len(gtemp_34_phi) > 0 and component[:4] == "dt13":
3393         gtemp_34.SetMarkerStyle(22);  gtemp_34.SetMarkerSize(1.25);
3394         gtemp_34.SetMarkerColor(ROOT.kGreen+2);  gtemp_34.SetLineColor(ROOT.kGreen+2)
3395         gtemp2_34.SetMarkerStyle(26); gtemp2_34.SetMarkerSize(1.25);
3396         gtemp2_34.SetMarkerColor(ROOT.kGreen+2); gtemp2_34.SetLineColor(ROOT.kGreen+2)
3397
3398     if wheel == 0: htemp.SetTitle("Wheel %d" % wheel)
3399     else: htemp.SetTitle("Wheel %+d" % wheel)
3400     htemp.SetAxisRange(-window, window, "Y")
3401     htemp.SetXTitle("Average #phi of pair (rad)")
3402     if component == "dt13_resid": htemp.SetYTitle("#Deltax^{local} (mm)")
3403     if component == "dt13_slope": htemp.SetYTitle("#Deltadx/dz^{local} (mrad)")
3404     if component == "dt2_resid": htemp.SetYTitle("#Deltay^{local} (mm)")
3405     if component == "dt2_slope": htemp.SetYTitle("#Deltady/dz^{local} (mrad)")
3406     htemp.GetXaxis().CenterTitle()
3407     htemp.GetYaxis().CenterTitle()
3408     htemp.GetYaxis().SetTitleOffset(0.75)
3409
3410     c1.Clear()
3411     htemp.Draw()
3412     if len(gtemp_12_phi) > 0:
3413         gtemp_12.Draw("p")
3414         gtemp2_12.Draw("p")
3415     if len(gtemp_23_phi) > 0:
3416         gtemp_23.Draw("p")
3417         gtemp2_23.Draw("p")
3418     if len(gtemp_34_phi) > 0:
3419         gtemp_34.Draw("p")
3420         gtemp2_34.Draw("p")
3421
3422     tlegend = ROOT.TLegend(0.5, 0.72, 0.9, 0.92)
3423     tlegend.SetBorderSize(0)
3424     tlegend.SetFillColor(ROOT.kWhite)
3425     if len(gtemp_12_phi) > 0:
3426         tlegend.AddEntry(gtemp_12, "MB1 - MB2 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_12_val), stdev(gtemp_12_val)), "pl")
3427     if len(gtemp_23_phi) > 0:
3428         tlegend.AddEntry(gtemp_23, "MB2 - MB3 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_23_val), stdev(gtemp_23_val)), "pl")
3429     if len(gtemp_34_phi) > 0:
3430         tlegend.AddEntry(gtemp_34, "MB3 - MB4 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_34_val), stdev(gtemp_34_val)), "pl")
3431     if len(gtemp_12_phi) > 0:
3432         tlegend.AddEntry(gtemp_12, "total mean: %4.2f, total RMS: %4.2f" % \
3433                                    (mean(gtemp_12_val + gtemp_23_val + gtemp_34_val),
3434                                    stdev(gtemp_12_val + gtemp_23_val + gtemp_34_val)), "")
3435     tlegend.Draw()
3436
3437
3438 ##################################################################################
3439
3440 def segdiffvsphicsc(tfiles, component, pair, window=5., **args):
3441     tdrStyle.SetOptTitle(1)
3442     tdrStyle.SetTitleBorderSize(1)
3443     tdrStyle.SetTitleFontSize(0.05)
3444
3445     if not component[0:3] == "csc": Exception
3446
3447     endcap = args["endcap"]
3448     if endcap=="m":
3449       endcapnum=2
3450       endcapsign="-"
3451     elif endcap=="p":
3452       endcapnum=1
3453       endcapsign="+"
3454     else: raise Exception
3455
3456     station1 = int(str(pair)[0])
3457     station2 = int(str(pair)[1])
3458     if not station2-station1==1: raise Exception
3459
3460     rings = [1,2]
3461     if station2==4: rings = [1]
3462
3463
3464     global htemp, gtemp_1, gtemp2_1, gtemp_2, gtemp2_2, tlegend
3465     htemp = ROOT.TH1F("htemp", "", 1, -pi*5./180., pi*(2.-5./180.))
3466     gtemp_1_phi, gtemp_1_val, gtemp_1_err, gtemp_1_val2, gtemp_1_err2 = [], [], [], [], []
3467     gtemp_2_phi, gtemp_2_val, gtemp_2_err, gtemp_2_val2, gtemp_2_err2 = [], [], [], [], []
3468
3469     for ring in rings:
3470       chambers = range(1,37)
3471       if ring == 1: chambers = range(1,19)
3472
3473       for chamber in chambers:
3474         phi, val, err, val2, err2, fit1, fit2, fit3 = segdiff(tfiles, component, pair, endcap=endcap, ring=ring, chamber=chamber)
3475         if fit1 and fit2 and fit3:
3476           if ring==1:
3477             gtemp_1_phi.append(phi)
3478             gtemp_1_val.append(val)
3479             gtemp_1_err.append(err)
3480             gtemp_1_val2.append(val2)
3481             gtemp_1_err2.append(err2)
3482           if ring==2:
3483             gtemp_2_phi.append(phi)
3484             gtemp_2_val.append(val)
3485             gtemp_2_err.append(err)
3486             gtemp_2_val2.append(val2)
3487             gtemp_2_err2.append(err2)
3488
3489     #print "len(gtemp_12_phi) ", len(gtemp_12_phi)
3490     #print "len(gtemp_23_phi) ",len(gtemp_23_phi)
3491     #print "len(gtemp_34_phi) ",len(gtemp_34_phi)
3492     if len(gtemp_1_phi) > 0:
3493         gtemp_1 = ROOT.TGraphErrors(len(gtemp_1_phi), array.array("d", gtemp_1_phi), array.array("d", gtemp_1_val),
3494                                      array.array("d", [0.] * len(gtemp_1_phi)), array.array("d", gtemp_1_err))
3495         gtemp2_1 = ROOT.TGraphErrors(len(gtemp_1_phi), array.array("d", gtemp_1_phi), array.array("d", gtemp_1_val2),
3496                                       array.array("d", [0.] * len(gtemp_1_phi)), array.array("d", gtemp_1_err2))
3497     if len(gtemp_2_phi) > 0:
3498         gtemp_2 = ROOT.TGraphErrors(len(gtemp_2_phi), array.array("d", gtemp_2_phi), array.array("d", gtemp_2_val),
3499                                      array.array("d", [0.] * len(gtemp_2_phi)), array.array("d", gtemp_2_err))
3500         gtemp2_2 = ROOT.TGraphErrors(len(gtemp_2_phi), array.array("d", gtemp_2_phi), array.array("d", gtemp_2_val2),
3501                                       array.array("d", [0.] * len(gtemp_2_phi)), array.array("d", gtemp_2_err2))
3502
3503     if len(gtemp_1_phi) > 0:
3504         gtemp_1.SetMarkerStyle(20);  gtemp_1.SetMarkerSize(1.);
3505         gtemp_1.SetMarkerColor(ROOT.kBlue);  gtemp_1.SetLineColor(ROOT.kBlue)
3506         gtemp2_1.SetMarkerStyle(24); gtemp2_1.SetMarkerSize(1.);
3507         gtemp2_1.SetMarkerColor(ROOT.kBlue); gtemp2_1.SetLineColor(ROOT.kBlue)
3508     if len(gtemp_2_phi) > 0:
3509         gtemp_2.SetMarkerStyle(21);  gtemp_2.SetMarkerSize(1.);
3510         gtemp_2.SetMarkerColor(ROOT.kRed);  gtemp_2.SetLineColor(ROOT.kRed)
3511         gtemp2_2.SetMarkerStyle(25); gtemp2_2.SetMarkerSize(1.);
3512         gtemp2_2.SetMarkerColor(ROOT.kRed); gtemp2_2.SetLineColor(ROOT.kRed)
3513
3514     htemp.SetTitle("ME%s%d - ME%s%d" % (endcapsign,station2,endcapsign,station1))
3515     htemp.SetAxisRange(-window, window, "Y")
3516     htemp.SetXTitle("Average #phi of pair (rad)")
3517     if component == "csc_resid": htemp.SetYTitle("#Delta(r#phi)^{local} (mm)")
3518     if component == "csc_slope": htemp.SetYTitle("#Deltad(r#phi)/dz^{local} (mrad)")
3519     htemp.GetXaxis().CenterTitle()
3520     htemp.GetYaxis().CenterTitle()
3521     htemp.GetYaxis().SetTitleOffset(0.75)
3522
3523     c1.Clear()
3524     htemp.Draw()
3525     if len(gtemp_1_phi) > 0:
3526         gtemp_1.Draw("p")
3527         gtemp2_1.Draw("p")
3528     if len(gtemp_2_phi) > 0:
3529         gtemp_2.Draw("p")
3530         gtemp2_2.Draw("p")
3531
3532     tlegend = ROOT.TLegend(0.5, 0.72, 0.9, 0.92)
3533     tlegend.SetBorderSize(0)
3534     tlegend.SetFillColor(ROOT.kWhite)
3535     if len(gtemp_1_phi) > 0:
3536         tlegend.AddEntry(gtemp_1, "ring 1 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_1_val), stdev(gtemp_1_val)), "pl")
3537     if len(gtemp_2_phi) > 0:
3538         tlegend.AddEntry(gtemp_2, "ring 2 (mean: %4.2f, RMS: %4.2f)" % (mean(gtemp_2_val), stdev(gtemp_2_val)), "pl")
3539     #if len(gtemp_12_phi) > 0:
3540     #    tlegend.AddEntry(gtemp_12, "total mean: %4.2f, total RMS: %4.2f" % \
3541     #                               (mean(gtemp_12_val + gtemp_23_val + gtemp_34_val),
3542     #                               stdev(gtemp_12_val + gtemp_23_val + gtemp_34_val)), "")
3543     tlegend.Draw()
3544
3545
3546
3547 ##################################################################################
3548 # makes a scatterplot of corrections coming either from reports (if xml geometries are None)
3549 # or from geometryX and geometryY (WRT the common initial geometry0)
3550
3551 def corrections2D(reportsX=None, reportsY=None, geometry0=None, geometryX=None, geometryY=None,
3552                   window=25., selection=None, name="tmp", canvas=None, pre_title_x=None, pre_title_y=None,
3553                   which="110011"):
3554
3555   tdrStyle.SetOptStat(0)
3556   tdrStyle.SetStatW(0.40)
3557
3558   # determine what are we plotting: report vs report  or  xml vs xml
3559   mode = None
3560   check_reports = False
3561   if reportsX is not None  and  reportsY is not None:
3562     mode = "reports"
3563     check_reports = True
3564   if geometry0 is not None  and  geometryX is not None  and  geometryY is not None:
3565     mode = "xmls"
3566   if mode is None:
3567     print("Either couple of reports or three geometries have to be given as input. Exiting...")
3568     return
3569
3570   # setup ranges with the maximum [-window,window] that later will be optimized to [-wnd_adaptive,wnd_adaptive]
3571   wnd = [window]*6
3572   wnd_adaptive = [.1]*6
3573
3574   global hx, hy, hz, hphix, hphiy, hphiz
3575   bins=2000
3576   hx    = ROOT.TH2F("%s_x" % name, "", bins, -wnd[0], wnd[0], bins, -wnd[0], wnd[0])
3577   hy    = ROOT.TH2F("%s_y" % name, "", bins, -wnd[1], wnd[1], bins, -wnd[1], wnd[1])
3578   hz    = ROOT.TH2F("%s_z" % name, "", bins, -wnd[2], wnd[2], bins, -wnd[2], wnd[2])
3579   hphix = ROOT.TH2F("%s_phix" % name, "", bins, -wnd[3], wnd[3], bins, -wnd[3], wnd[3])
3580   hphiy = ROOT.TH2F("%s_phiy" % name, "", bins, -wnd[4], wnd[4], bins, -wnd[4], wnd[4])
3581   hphiz = ROOT.TH2F("%s_phiz" % name, "", bins, -wnd[5], wnd[5], bins, -wnd[5], wnd[5])
3582   hhh = [hx, hy, hz, hphix, hphiy, hphiz]
3583
3584   # initialize PCA objects
3585   global pca_x, pca_y, pca_z, pca_phix, pca_phiy, pca_phiz
3586   pca_x = ROOT.TPrincipal(2,"D")
3587   pca_y = ROOT.TPrincipal(2,"D")
3588   pca_z = ROOT.TPrincipal(2,"D")
3589   pca_phix = ROOT.TPrincipal(2,"D")
3590   pca_phiy = ROOT.TPrincipal(2,"D")
3591   pca_phiz = ROOT.TPrincipal(2,"D")
3592   pcas = [pca_x, pca_y, pca_z, pca_phix, pca_phiy, pca_phiz]
3593
3594   # arrays to later fill graphs with
3595   ax=[]; ay=[]; az=[]; aphix=[]; aphiy=[]; aphiz=[]
3596   aaa = [ax, ay, az, aphix, aphiy, aphiz]
3597
3598   # list of postal addresses
3599   postal_addresses = []
3600
3601   # if reports are given, use them to fill addresses and do extra checks
3602   if check_reports:
3603     for r1 in reportsX:
3604       # skip ME1/a
3605       if r1.postal_address[0]=='CSC'  and  r1.postal_address[2]==1 and r1.postal_address[3]==4: continue
3606       if selection is None or (selection.__code__.co_argcount == len(r1.postal_address) and selection(*r1.postal_address)):
3607         r2 = getReportByPostalAddress(r1.postal_address, reportsY)
3608         if r2 is None:
3609           print("bad r2 in ",r1.postal_address)
3610           continue
3611
3612         if r1.status != "PASS" or r2.status != "PASS":
3613           print("bad status", r1.postal_address, r1.status, r2.status)
3614           continue
3615       postal_addresses.append(r1.postal_address)
3616   # otherwise, use chamber addresses from xmls
3617   else:
3618     for key in geometry0.dt.keys():
3619       if len(key)==3  and  key in geometryX.dt  and  key in geometryY.dt:
3620         postal_addresses.append( tuple(['DT'] + list(key)) )
3621     for key in geometry0.csc.keys():
3622       # skip ME1/a
3623       if key[2]==1 and key[3]==4: continue
3624       if len(key)==4  and  key in geometryX.csc  and  key in geometryY.csc:
3625         postal_addresses.append( tuple(['CSC'] + list(key)) )
3626
3627   # fill the values
3628   for addr in postal_addresses:
3629
3630     # checks the selection function
3631     if not (selection is None or (selection.__code__.co_argcount == len(addr) and selection(*addr)) ): continue
3632
3633     factors = [10. * signConventions[addr][0], 10. * signConventions[addr][1], 10. * signConventions[addr][2],
3634                1000., 1000., 1000. ]
3635
3636     if check_reports:
3637       rX = getReportByPostalAddress(addr, reportsX)
3638       rY = getReportByPostalAddress(addr, reportsY)
3639       deltasX = [rX.deltax, rX.deltay, rX.deltaz, rX.deltaphix, rX.deltaphiy, rX.deltaphiz]
3640       deltasY = [rY.deltax, rY.deltay, rY.deltaz, rY.deltaphix, rY.deltaphiy, rY.deltaphiz]
3641
3642     if mode == "reports":
3643
3644       checks = map( lambda d1, d2: d1 is not None  and  d2 is not None  and  d1.error is not None   \
3645                                    and  d2.error is not None and (d1.error**2 + d2.error**2) > 0. , \
3646                     deltasX, deltasY)
3647
3648       for i in range(len(checks)):
3649         if not checks[i]: continue
3650         fillX = deltasX[i].value * factors[i]
3651         fillY = deltasY[i].value * factors[i]
3652         aaa[i].append([fillX,fillY])
3653         pcas[i].AddRow(array.array('d',[fillX,fillY]))
3654         mx = max(abs(fillX), abs(fillY))
3655         if mx > wnd_adaptive[i]: wnd_adaptive[i] = mx
3656
3657     if mode == "xmls":
3658
3659       db0 = dbX = dbY = None
3660       if addr[0] == "DT":
3661         db0, dbX, dbY  = geometry0.dt[addr[1:]], geometryX.dt[addr[1:]], geometryY.dt[addr[1:]]
3662       if addr[0] == 'CSC':
3663         db0, dbX, dbY  = geometry0.csc[addr[1:]], geometryX.csc[addr[1:]], geometryY.csc[addr[1:]]
3664
3665       checks = [True]*6
3666       if check_reports:
3667         checks = map( lambda d1, d2: d1 is not None  and  d2 is not None ,  deltasX, deltasY)
3668
3669       gdeltas0 = [db0.x, db0.y, db0.z, db0.phix, db0.phiy, db0.phiz]
3670       gdeltasX = [dbX.x, dbX.y, dbX.z, dbX.phix, dbX.phiy, dbX.phiz]
3671       gdeltasY = [dbY.x, dbY.y, dbY.z, dbY.phix, dbY.phiy, dbY.phiz]
3672
3673       for i in range(len(checks)):
3674         if not checks[i]: continue
3675         fillX = (gdeltasX[i] - gdeltas0[i]) * factors[i]
3676         fillY = (gdeltasY[i] - gdeltas0[i]) * factors[i]
3677         aaa[i].append([fillX,fillY])
3678         pcas[i].AddRow(array.array('d',[fillX,fillY]))
3679         mx = max(abs(fillX), abs(fillY))
3680         if mx > wnd_adaptive[i]: wnd_adaptive[i] = mx
3681         #if addr[0] == 'CSC' and i==1 and (abs(fillX)>0.01 or abs(fillY)>0.01): print addr, ": hugeCSC i=%d dx=%.03g dy=%.03g"%(i,fillX,fillY)
3682         #if addr[0] == 'CSC' and i==2 and (abs(fillX)>0.02 or abs(fillY)>0.02): print addr, ": hugeCSC i=%d dx=%.03g dy=%.03g"%(i,fillX,fillY)
3683         #if addr[0] == 'CSC' and i==3 and (abs(fillX)>0.05 or abs(fillY)>0.05): print addr, ": hugeCSC i=%d dx=%.03g dy=%.03g"%(i,fillX,fillY)
3684
3685   if mode == "xmls":
3686     if pre_title_x is None: pre_title_x = "geometry 1 "
3687     if pre_title_y is None: pre_title_y = "geometry 2 "
3688   if mode == "reports":
3689     if pre_title_x is None: pre_title_x = "iteration's "
3690     if pre_title_y is None: pre_title_y = "other iteration's "
3691   tmptitles =  ["#Deltax (mm)", "#Deltay (mm)", "#Deltaz (mm)",
3692                 "#Delta#phi_{x} (mrad)", "#Delta#phi_{y} (mrad)", "#Delta#phi_{z} (mrad)"]
3693   htitles = []
3694   for t in tmptitles: htitles.append([pre_title_x + t, pre_title_y + t])
3695
3696   if canvas is not None: c = canvas
3697   else: c = c1
3698   c.Clear()
3699   ndraw = which.count('1')
3700   if ndraw > 4: c.Divide(3, 2)
3701   elif ndraw > 2: c.Divide(2, 2)
3702   elif ndraw > 1: c.Divide(2, 1)
3703
3704   global lines, graphs, texs
3705   lines = [];  graphs = []; texs = []
3706
3707   ipad = 0
3708   for i in range(6):
3709
3710     # decode 'which' binary mask
3711     if ( int(which,2) & (1<<i) ) == 0: continue
3712
3713     ipad += 1
3714     c.GetPad(ipad).cd()
3715     c.GetPad(ipad).SetGridx(1)
3716     c.GetPad(ipad).SetGridy(1)
3717
3718     wn = 1.08 * wnd_adaptive[i]
3719     hhh[i].GetXaxis().SetRangeUser(-wn, wn)
3720     hhh[i].GetYaxis().SetRangeUser(-wn, wn)
3721     hhh[i].SetXTitle(htitles[i][0])
3722     hhh[i].SetYTitle(htitles[i][1])
3723     hhh[i].GetXaxis().CenterTitle()
3724     hhh[i].GetYaxis().CenterTitle()
3725     hhh[i].Draw()
3726
3727     if len(aaa[i]) == 0: continue
3728
3729     a1, a2 = map( lambda x: array.array('d',x),  list(zip(*aaa[i])) )
3730     g = ROOT.TGraph(len(a1), a1, a2)
3731     g.SetMarkerStyle(5)
3732     g.SetMarkerSize(0.3)
3733     g.SetMarkerColor(ROOT.kBlue)
3734     graphs.append(g)
3735
3736     pcas[i].MakePrincipals()
3737     #pcas[i].Print()
3738     #pcas[i].MakeHistograms()
3739     b = pcas[i].GetEigenVectors()(1,0) / pcas[i].GetEigenVectors()(0,0)
3740     a = pcas[i].GetMeanValues()[1] - b * pcas[i].GetMeanValues()[0]
3741     #print a, b, "   ", pcas[i].GetEigenValues()[0], pcas[i].GetEigenValues()[1]
3742
3743     cov = pcas[i].GetCovarianceMatrix()
3744     r = cov(0,1)/sqrt(cov(1,1)*cov(0,0))
3745     print("r, RMSx, RMSy =", r, g.GetRMS(1), g.GetRMS(2))
3746     texrms = ROOT.TLatex(0.17,0.87, "RMS x,y = %.02g, %.02g" % (g.GetRMS(1),g.GetRMS(2)))
3747     texr = ROOT.TLatex(0.17,0.80, "r = %.02g" % r)
3748     for t in texr, texrms:
3749       t.SetNDC(1)
3750       t.SetTextColor(ROOT.kBlue)
3751       t.SetTextSize(0.053)
3752       t.Draw()
3753       texs.append(t)
3754
3755     g.Draw("p")
3756
3757     if not isnan(b):
3758       wn = wnd_adaptive[i]
3759       line = ROOT.TLine(-wn, a - b*wn, wn, a + b*wn)
3760       line.SetLineColor(ROOT.kRed)
3761       line.Draw()
3762       lines.append(line)
3763
3764   #return hx, hy, hphiy, hphiz, pca_x, pca_y, pca_phiy, pca_phiz
3765   return aaa