MuonAlignmentAlgorithms/scripts/plotscripts.py

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