alignmentValidation.py

Go to the documentation of this file.

#! /usr/bin/env python

import re,os,sys,shutil
import optparse

from mutypes import *

execfile("plotscripts.py")

ROOT.gROOT.SetBatch(1);

######################################################

######################################################
# To parse commandline args

usage='%prog [options]\n'+\
  'This script dumps muon alignment validation plots '+\
  '(see https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideMuonAlignValidationPlots) '+\
  'in web-friendly png format into a predefined directory structure '+\
  'and also (to be implemented) does some quality checks for these plots.\n'+\
  'Script uses output of the first and last iterations of muon alignment runs '+\
  'that should be located in inputDir+i1 and inputDir+iN directories (see options descriptions). '+\
  'Each of these directories should contain\ni#prefix+"_plotting.root"\ni#prefix+".root"\ni#prefix+"_report.py"\n'+\
  'files, where # = 1 or N.\n'+\
  'Output will be located in outputDir, which must exist. Plots from the first and last iterations '+\
  'will be in outputDir+"iter1" and outputDir+"iterN" respectively, and plots relevant for all iterations will '+\
  'reside in outputDir+"common/".\n'+\
  'For a first run a --createDirSructure option must be present. If this option is present for subsequent runs '+\
  'with the same outputDir, all previous results in "iter1", "iterN" and and "common" would be deleted!\n'+\
  'If neither --dt or --csc is present, plots for both systems will be created.\n'+\
  'Options must include either -a or any of the following: --map, --segdiff, --fit, --median'

parser=optparse.OptionParser(usage)

parser.add_option("-l", "--runLabel",
  help="[REQUIRED] label to use for a run",
  type="string",
  default='',
  dest="runLabel")

parser.add_option("-i", "--inputDir",
  help="[REQUIRED] input directory: should contain directories with the first and the final iterations' results",
  type="string",
  default='',
  dest="inputDir")

parser.add_option("--i1",
  help="[REQUIRED] directory with the alignment 1st iteration's results relative to inputDir",
  type="string",
  default='',
  dest="i1")

parser.add_option("--iN",
  help="[REQUIRED] directory with the alignment last iteration's results relative to inputDir",
  type="string",
  default='',
  dest="iN")

parser.add_option("--i1prefix",
  help="filename prefix for the alignment 1st iteration's results. If not provided, i1prefix = i1",
  type="string",
  default='',
  dest="i1prefix")

parser.add_option("--iNprefix",
  help="filename prefix for the alignment last iteration's results. If not provided, iNprefix = iN",
  type="string",
  default='',
  dest="iNprefix")

parser.add_option("-o", "--outputDir",
  help="output directory: all plots will be saved with relation to outputDir. If not provided, consider outputDir = inputDir",
  type="string",
  default='',
  dest="outputDir")

parser.add_option("--createDirSructure",
  help="If present, new directory structure for storing plots will be first created for each iteration at outputDir+i1 and outputDir+iN. WARNING: this will delete any existing results!",
  action="store_true",
  default=False,
  dest="createDirSructure")

parser.add_option("--dt",
  help="If it is present, but not --csc, DT only plots will be created",
  action="store_true",
  default=False,
  dest="dt")

parser.add_option("--csc",
  help="If this is present, but not --dt, CSC only plots will be created",
  action="store_true",
  default=False,
  dest="csc")

parser.add_option("-a","--all",
  help="If present, all types of plots will be created",
  action="store_true",
  default=False,
  dest="all")

parser.add_option("--map",
  help="If present, map plots will be created",
  action="store_true",
  default=False,
  dest="map")

parser.add_option("--segdiff",
  help="If present, segdiff plots will be created",
  action="store_true",
  default=False,
  dest="segdiff")

parser.add_option("--curvature",
  help="If present, curvature plots will be created",
  action="store_true",
  default=False,
  dest="curvature")

parser.add_option("--fit",
  help="If present, fit functions plots will be created",
  action="store_true",
  default=False,
  dest="fit")

parser.add_option("--median",
  help="If present, median plots will be created",
  action="store_true",
  default=False,
  dest="median")

parser.add_option("--diagnostic",
  help="If present, will run diagnostic checks",
  action="store_true",
  default=False,
  dest="diagnostic")

parser.add_option("-v", "--verbose",
  help="Degree of debug info verbosity",
  type="int",
  default=0,
  dest="verbose")

options,args=parser.parse_args()

if options.runLabel=='' or options.inputDir=='' or options.i1=='' or options.iN=='':
  print "\nOne or more of REQUIRED options is missing!\n"
  parser.print_help()
  # See \n"+sys.argv[0]+" --help"
  sys.exit()

outdir = options.outputDir
if outdir=='': outdir = options.inputDir

i1prefix = options.i1prefix
if i1prefix=='' : i1prefix = options.i1

iNprefix = options.iNprefix
if iNprefix=='' : iNprefix = options.iN

if not os.access(outdir,os.F_OK):
  print "\noutDir = "+outdir+"\ndoes not exist! Exiting..."
  sys.exit()

# If neither --dt or --csc is present, plots for both systems will be created
DO_DT  = False
DO_CSC = False
if options.dt or not ( options.dt or options.csc):
  DO_DT = True
if options.csc or not ( options.dt or options.csc):
  DO_CSC = True

if not (options.all or options.map or options.curvature or options.segdiff or options.fit or options.median or options.diagnostic):
  print "\nOptions must include either -a or any of the following: --map, --segdiff, --fit, --median, --diagnostic. Exiting..."
  sys.exit()

SINGLE_ITERATION = False
if i1prefix == iNprefix: SINGLE_ITERATION = True

DO_MAP = False
DO_SEGDIFF = False
DO_CURVATURE = False
DO_FIT = False
DO_MEDIAN = False
if options.map or options.all:
  DO_MAP = True
if options.segdiff or options.all:
  DO_SEGDIFF = True
if options.curvature or options.all:
  DO_CURVATURE = True
if options.fit or options.all:
  DO_FIT = True
if options.median or options.all:
  DO_MEDIAN = True

DO_DIAGNOSTIC = options.diagnostic

allOptions = "-l "+options.runLabel+" -i "+options.inputDir+" --i1 "+options.i1+" --iN "+options.iN
if options.i1prefix !='': allOptions += " --i1prefix " + options.i1prefix
if options.iNprefix !='': allOptions += " --iNprefix " + options.iNprefix
allOptions += " -o "+options.outputDir
if options.createDirSructure: allOptions += " --createDirSructure"
if DO_DT: allOptions += " --dt"
if DO_CSC: allOptions += " --csc"
if options.all: allOptions += " -a"
if options.map: allOptions += " --map"
if options.segdiff: allOptions += " --segdiff"
if options.curvature: allOptions += " --curvature"
if options.fit: allOptions += " --fit"
if options.median: allOptions += " --median"
if options.diagnostic: allOptions += " --diagnostic"
print sys.argv[0]+" "+allOptions


QUICKTESTN=10000



######################################################

# template for canvases list
CANVASES_LIST_TEMPLATE = [
['Common',' ',
 ['medians distribution','medians.png']
],
['DT',' ',
 ['Wheel&Station: map of dxdz residual vs phi','map_DTvsphi_dxdz.png'],
 ['Wheel&Station: map of dydz residual vs phi','map_DTvsphi_dydz.png'],
 ['Wheel&Station: map of x residual vs phi','map_DTvsphi_x.png'],
 ['Wheel&Station: map of y residual vs phi','map_DTvsphi_y.png'],
 ['Station&Sector: map of dxdz residual vs z','map_DTvsz_dxdz.png'],
 ['Station&Sector: map of dydz residual vs z','map_DTvsz_dydz.png'],
 ['Station&Sector: map of x residual vs z','map_DTvsz_x.png'],
 ['Station&Sector: map of y residual vs z','map_DTvsz_y.png'],
 ['Station: map of dxdz residual vs z','map_DTvsz_all_dxdz.png'],
 ['Station: map of dydz residual vs z','map_DTvsz_all_dydz.png'],
 ['Station: map of x residual vs z','map_DTvsz_all_x.png'],
 ['Station: map of y residual vs z','map_DTvsz_all_y.png'],
 ['Wheel: segdiff in x residuals vs phi','segdifphi_dt13_resid.png'],
 ['Wheel: segdiff in dxdz residuals vs phi','segdifphi_dt13_slope.png'],
 ['Wheel: segdiff in y residuals vs phi','segdifphi_dt2_resid.png'],
 ['Wheel: segdiff in dydz residuals vs phi','segdifphi_dt2_slope.png'],
 ['Wheel: segdiff DT-CSC in x residuals vs phi','segdifphi_x_dt_csc_resid.png'],
 ['Chamber: segdiff in x residuals','segdif_dt13_resid.png'],
 ['Chamber: segdiff in dxdz residuals','segdif_dt13_slope.png'],
 ['Chamber: segdiff in y residuals','segdif_dt2_resid.png'],
 ['Chamber: segdiff in dydz residuals','segdif_dt2_slope.png'],
 ['Chamber: segdiff DT-CSC in x residuals','segdif_x_dt_csc_resid.png'],
 ['Chamber: residuals distributions','dt_bellcurves.png'],
 ['Chamber: residuals relations to misalignments','dt_polynomials.png'],
 ['Chamber: Delta x residuals vs. curvature','dt_curvature_deltax.png'],
 ['Chamber: Delta dxdz residuals vs. curvature','dt_curvature_deltadxdz.png'],
 ['Extras: Extra plots in a separate window','dt_extras.php']
],
['CSC',' ',
 ['Station&Ring: map of d(rphi)/dz residual vs phi','map_CSCvsphi_dxdz.png'],
 ['Station&Ring: map of rphi residual vs phi','map_CSCvsphi_x.png'],
 ['Station&Chamber: map of d(rphi)/dz residual vs r','map_CSCvsr_dxdz.png'],
 ['Station&Chamber: map of rphi residual vs r','map_CSCvsr_x.png'],
 ['Station: map of d(rphi)/dz residual vs r','map_CSCvsr_all_dxdz.png'],
 ['Station: map of rphi residual vs r','map_CSCvsr_all_x.png'],
 ['Station: segdiff in rphi residuals vs phi','segdifphi_csc_resid.png'],
 ['Station: segdiff in d(rphi)/dz residuals vs phi','segdifphi_csc_slope.png'],
 ['Chamber: segdiff in rphi residuals','segdif_csc_resid.png'],
 ['Chamber: segdiff in d(rphi)/dz residuals','segdif_csc_slope.png'],
 ['Chamber: residuals distributions','csc_bellcurves.png'],
 ['Chamber: residuals relations to misalignments','csc_polynomials.png'],
 #['Chamber: Delta rphi residuals vs. curvature','csc_curvature_x.png'],
 #['Chamber: Delta d(rphi)/dz residuals vs. curvature','csc_curvature_dxdz.png'],
 ['Extras: Extra plots in a separate window','csc_extras.php']
]
]


######################################################
# functions definitions


def isFileUnderDir(dir_name, file_name):
  '''Recursively looks for file named file_name under dir_name directory
  '''
  if not DO_DT and dir_name.find("MB")>-1:
    return False
  if not DO_CSC and dir_name.find("ME")>-1:
    return False
  
  for f in os.listdir(dir_name):
    dirfile = os.path.join(dir_name, f)
    if os.path.isfile(dirfile) and f==file_name:
      return True
    # recursively access file names in subdirectories
    elif os.path.isdir(dirfile):
      #print "Accessing directory:", dirfile
      if isFileUnderDir(dirfile, file_name): return True
  return False


# to time saving of plots
def saveAs(nm): 
  t1 = time.time()
  ddt[15] += 1
  c1.SaveAs(nm)
  tn = time.time()
  ddt[16] = 1./ddt[15]*((ddt[15]-1)*ddt[16] + tn-t1)


def createDirectoryStructure(iteration_name):
  
  if not os.access(iteration_name,os.F_OK):
    os.mkdir(iteration_name)

  csc_basedir = iteration_name+'/'
  for endcap in CSC_TYPES:
    #print csc_basedir+endcap[0]
    shutil.rmtree(csc_basedir+endcap[0],True)
    os.mkdir(csc_basedir+endcap[0])
    for station in endcap[2]:
      #print csc_basedir+endcap[0]+'/'+station[1]
      os.mkdir(csc_basedir+endcap[0]+'/'+station[1])
      for ring in station[2]:
        #print csc_basedir+endcap[0]+'/'+station[1]+'/'+ring[1]
        os.mkdir(csc_basedir+endcap[0]+'/'+station[1]+'/'+ring[1])
        for chamber in range(1,ring[2]+1):
          schamber = "%02d" % chamber
          #print csc_basedir+endcap[0]+'/'+station[1]+'/'+ring[1]+'/'+schamber
          os.mkdir(csc_basedir+endcap[0]+'/'+station[1]+'/'+ring[1]+'/'+schamber)

  dt_basedir = iteration_name+'/MB/'
  #print dt_basedir
  shutil.rmtree(dt_basedir,True)
  os.mkdir(dt_basedir)
  for wheel in DT_TYPES:
    #print dt_basedir+wheel[0]
    os.mkdir(dt_basedir+wheel[0])
    for station in wheel[2]:
      #print dt_basedir+wheel[0]+'/'+station[1]
      os.mkdir(dt_basedir+wheel[0]+'/'+station[1])
      for sector in range(1,station[2]+1):
        ssector = "%02d" % sector
        #print dt_basedir+wheel[0]+'/'+station[1]+'/'+ssector
        os.mkdir(dt_basedir+wheel[0]+'/'+station[1]+'/'+ssector)

  print os.getcwd()

######################################################

def doMapPlotsDT(dt_basedir, tfiles_plotting):
  """write DT map plots

 "DTvsphi_st%dwh%s" % (station, wheelletter):

  plots "integrated" over ALL SECTORS:
  of x, y, dxdz, dydz vs. phi (y and dydz only for stations 1-3)
  made for all (station,wheel) combinations

  Access interface may be arranged into station(1 .. 4) vs. wheel(-2 .. +2) map.
  It could be incorporated into a general DT chambers map (column1: wheel, column2: station,
  columns3-16 correspond to sector #) by making station numbers in column 2 clickable.
  

 "DTvsz_st%dsec%02d" % (station, sector)

  plots "integrated" over ALL WHEELS:
  of x, y, dxdz, dydz vs. z (y and dydz only for stations 1-3)
  made for all (station,sector) combinations

  Interface: may be arranged into station(1 .. 4) vs. sector(1 .. 14) map with sector range
  (1 .. 12) for stations 1-3.
  It could be incorporated into an EXTENDED general DT chambers map (extended by adding an
  identifier "ALL" in column1 for wheel number).


 "DTvsz_st%dsecALL" % (station)

  plots spanning in z over ALL WHEELS and "integrated" over all sectors:
  of x, y, dxdz, dydz vs. z (y and dydz only for stations 1-3)
  made for all stations

  Interface: may be arranged into station(1 .. 4) map 
  It could be incorporated into an EXTENDED general DT chambers map (extended by adding an
  identifier "ALL" in column1 for wheel number)."""
  
  
  for wheel in DT_TYPES:
    if wheel[1]=="ALL": continue
    for station in wheel[2]:
      pdir = dt_basedir+'/'+wheel[0]+'/'+station[1]+'/'
      label = "DTvsphi_st%dwh%s" % (int(station[1]), wheelLetter(int(wheel[1])))
      htitle = "wheel %+d, station %s" % (int(wheel[1]), station[1])
      #mapplot(tfiles_plotting, label, "x", window=25., title=htitle, fitsawteeth=True,fitsine=True)
      mapplot(tfiles_plotting, label, "x", window=10., title=htitle, fitsine=True,fitpeaks=True,peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsphi_x.png')
      #mapplot(tfiles_plotting, label, "dxdz", window=25., title=htitle, fitsawteeth=True,fitsine=True)
      mapplot(tfiles_plotting, label, "dxdz", window=10., title=htitle,peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsphi_dxdz.png')

      if station[1]=='4': continue
      #mapplot(tfiles_plotting, label, "y", window=25., title=htitle, fitsawteeth=True,fitsine=True)
      mapplot(tfiles_plotting, label, "y", window=10., title=htitle,peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsphi_y.png')
      #mapplot(tfiles_plotting, label, "dydz", window=25., title=htitle, fitsawteeth=True,fitsine=True)
      mapplot(tfiles_plotting, label, "dydz", window=10., title=htitle,peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsphi_dydz.png')

  
  qcount=0
  for wheel in DT_TYPES:
    if wheel[1]!="ALL": continue
    for station in wheel[2]:
      for sector in range(1,station[2]+1):
        if qcount>QUICKTESTN: break
        qcount += 1
        ssector = "%02d" % sector
        pdir = dt_basedir+'/'+wheel[0]+'/'+station[1]+'/'+ssector+'/'
        label = "DTvsz_st%ssec%s" % (station[1], ssector)
        htitle = "station %s, sector %d" % (station[1], sector)
        mapplot(tfiles_plotting, label, "x", window=10., title=htitle, peaksbins=2)
        c1.SaveAs(pdir+'map_DTvsz_x.png')
        mapplot(tfiles_plotting, label, "dxdz", window=10., title=htitle, peaksbins=2)
        c1.SaveAs(pdir+'map_DTvsz_dxdz.png')
        
        if station[1]=='4': continue
    mapplot(tfiles_plotting, label, "y", window=10., title=htitle, peaksbins=2)
        c1.SaveAs(pdir+'map_DTvsz_y.png')
        mapplot(tfiles_plotting, label, "dydz", window=10., title=htitle, peaksbins=2)
        c1.SaveAs(pdir+'map_DTvsz_dydz.png')
  
  qcount=0
  for wheel in DT_TYPES:
    if wheel[1]!="ALL": continue
    for station in wheel[2]:
      if qcount>QUICKTESTN: break
      qcount += 1
      pdir = dt_basedir+'/'+wheel[0]+'/'+station[1]+'/'
      label = "DTvsz_st%ssecALL" % (station[1])
      htitle = "station %s" % (station[1])
      
      print label, 
      mapplot(tfiles_plotting, label, "x", window=10., title=htitle, peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsz_all_x.png')
      mapplot(tfiles_plotting, label, "dxdz", window=10., title=htitle, peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsz_all_dxdz.png')
      
      if station[1]=='4': continue
      mapplot(tfiles_plotting, label, "y", window=10., title=htitle, peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsz_all_y.png')
      mapplot(tfiles_plotting, label, "dydz", window=10., title=htitle, peaksbins=2)
      c1.SaveAs(pdir+'map_DTvsz_all_dydz.png')

  saveTestResultsMap(options.runLabel)


def doMapPlotsCSC(csc_basedir, tfiles_plotting):
  """write CSC map plots

 "CSCvsphi_me%s%d%d" % (endcap, station, ring)

  plots "integrated" over ALL SECTORS:
  of rphi, drphi/dz vs. phi
  made for all (endcap,station,ring) combinations

  Interface: may be arranged into two station(1 .. 4) vs. R(1 .. 4) maps for both endcaps
  with R range (1 .. 4) for stations 2-4
 It could be incorporated into a general CSC chambers map (column1: endcap, column2: station,
  column3: ring, columns4-40 correspond to chamber #) by making ring numbers in column 3
  clickable.


 "CSCvsr_me%s%dch%02d" % (endcap, station, chamberNumber)

  plots "integrated" over ALL RINGS:
  of rphi, drphi/dz vs. z
  made for all (endcap,station,chamber) combinations

  Interface: may be arranged into two station(1 .. 4) vs. chamber(1 .. 36) maps for both endcaps
  It could be incorporated into an EXTENDED general CSC chambers map (extended by adding an
  identifier "ALL" in column3 for ring number).
  
 "CSCvsr_me%s%dchALL" % (endcap, station)

  plots spanning over ALL RINGS along r and integrated over all SECTORS:
  of rphi, drphi/dz vs. z
  made for all (endcap,station) combinations

  Interface: may be arranged into two station(1 .. 4) maps for both endcaps
  It could be incorporated into an EXTENDED general CSC chambers map (extended by adding an
  identifier "ALL" in column3 for ring number)."""

  for endcap in CSC_TYPES:
    for station in endcap[2]:
      for ring in station[2]:
        if ring[1]=="ALL": continue
        pdir = csc_basedir+'/'+endcap[0]+'/'+station[1]+'/'+ring[1]+'/'
        label = "CSCvsphi_me%s%s%s" % (endcap[1], station[1], ring[1])
        htitle = "%s%s/%s" % (endcap[0], station[1],ring[1])
        mapplot(tfiles_plotting, label, "x", window=15., title=htitle, fitsine=True,fitpeaks=True, peaksbins=2)
        #mapplot(tfiles_plotting, label, "x", window=15., title=htitle, fitsine=True)
        c1.SaveAs(pdir+'map_CSCvsphi_x.png')
        mapplot(tfiles_plotting, label, "dxdz", window=15., title=htitle, peaksbins=2)
        c1.SaveAs(pdir+'map_CSCvsphi_dxdz.png')

  saveTestResultsMap(options.runLabel)

  qcount = 0
  for endcap in CSC_TYPES:
    for station in endcap[2]:
      for ring in station[2]:
        if ring[1]!="ALL": continue
        for chamber in range(1,ring[2]+1):
          if qcount>QUICKTESTN: break
          qcount += 1
          schamber = "%02d" % chamber
          pdir = csc_basedir+'/'+endcap[0]+'/'+station[1]+'/'+ring[1]+'/'+schamber+'/'
          label = "CSCvsr_me%s%sch%s" % (endcap[1], station[1], schamber)
          htitle = "%s%s/ALL/%d" % (endcap[0], station[1],chamber)
          mapplot(tfiles_plotting, label, "x", window=15., title=htitle, peaksbins=2)
          c1.SaveAs(pdir+'map_CSCvsr_x.png')
          mapplot(tfiles_plotting, label, "dxdz", window=15., title=htitle, peaksbins=2)
          c1.SaveAs(pdir+'map_CSCvsr_dxdz.png')

  qcount = 0
  for endcap in CSC_TYPES:
    for station in endcap[2]:
      for ring in station[2]:
        if ring[1]!="ALL": continue
        if qcount>QUICKTESTN: break
        qcount += 1
        schamber = "%02d" % chamber
        pdir = csc_basedir+'/'+endcap[0]+'/'+station[1]+'/'+ring[1]+'/'
        label = "CSCvsr_me%s%schALL" % (endcap[1], station[1])
        htitle = "%s%s" % (endcap[0], station[1])
        mapplot(tfiles_plotting, label, "x", window=10., title=htitle, peaksbins=2)
        c1.SaveAs(pdir+'map_CSCvsr_all_x.png')
        mapplot(tfiles_plotting, label, "dxdz", window=10., title=htitle, peaksbins=2)
        c1.SaveAs(pdir+'map_CSCvsr_all_dxdz.png')



def doCurvaturePlotsDT(dt_basedir, tfiles_plotting):
  """write DT curvature plots

  "wheel%s_sector%s" % (wheel, sector)

  wheel in "m2", "m1", "z", "p1", "p2"
  station 1 only!
  sector in "01", ..., "12"

  "param" may be one of 
    "deltax" (Delta x position residuals),
    "deltadxdz" (Delta (dx/dz) angular residuals),
    "curverr" (Delta x * d(Delta q/pT)/d(Delta x) = Delta q/pT in the absence of misalignment) - not necessary

  made for all (wheel,station=1,sector) combinations

  Interface: could be accesses through a general DT chambers map for station=1 chambers."""
  
  w_dict = {'-2':'m2', '-1':'m1', '0':'z', '1':'p1', '2':'p2'}
  qcount = 0
  for wheel in DT_TYPES:
    if wheel[1]=="ALL": continue
    #station = 1
    #station = wheel[2][0]
    for station in wheel[2]:
      print "curv in ", wheel[0]+'/'+station[1]
      for sector in range(1,station[2]+1):
        #if sector>12: break
        if qcount>QUICKTESTN: break
        qcount += 1
        ssector = "%02d" % sector
        pdir = dt_basedir+'/'+wheel[0]+'/'+station[1]+'/'+ssector+'/'
        label = "wheel%s_st%s_sector%s" % (w_dict[wheel[1]], station[1], ssector)
        thetitle ="Wheel %s, station %s, sector %s" % (wheel[1], station[1], ssector)
        curvatureplot(tfiles_plotting, label, "deltax", title=thetitle, window=10., fitline=True)
        saveAs(pdir+'dt_curvature_deltax.png')
        curvatureplot(tfiles_plotting, label, "deltadxdz", title=thetitle, window=10., fitline=True)
        saveAs(pdir+'dt_curvature_deltadxdz.png')



def doSegDiffPlotsDT(dt_basedir, tfiles_plotting, iter_reports):
  """write segment-difference plots for DT

 segdiff "dt13_resid" and "dt13_slope"

  set of plots of
  x vs qpt, x for positive, x for negative ("dt13_resid")
  dxdz vs qpt, dxdz for positive, dxdz for negative ("dt13_slope")
  done for MB1-MB2, MB2-MB3, and MB3-MB4 stations combinations with all possible (wheel, sector)

  Interface: could be accessed through a general DT chambers map, but only for chambers in
  stations 2-4 (e.g., station 2 would provide MB1-MB2 plots).

 segdiff "dt2_resid" and "dt2_slope"

  set of plots of
  y vs q/pt, y for positive, y for negative ("dt2_resid")
  dydz vs q/pt, dydz for positive, dydz for negative ("dt2_slope")
  done for MB1-MB2, MB2-MB3 stations combinations with all possible (wheel, sector)

  Interface: then the interface would still be a general DT map,
  but the info only available from station 2 & 3 chambers."""

  qcount = 0
  for iwheel in DT_TYPES:
    if iwheel[1]=="ALL": continue
    for istation in iwheel[2]:
      if istation[1]=="1": continue
      dstations = (int(istation[1])-1)*10 + int(istation[1])
      #print dstations
      for isector in range(1, istation[2] + 1):
        if isector > 12: continue
        if qcount>QUICKTESTN: break
        qcount += 1
        ssector = "%02d" % isector
        pdir = dt_basedir + '/' + iwheel[0] + '/' + istation[1] + '/' + ssector + '/'
        
        segdiff(tfiles_plotting, "dt13_resid", dstations, wheel=int(iwheel[1]), sector=isector, window=15.)
        c1.SaveAs(pdir + 'segdif_dt13_resid.png')
        segdiff(tfiles_plotting, "dt13_slope", dstations, wheel=int(iwheel[1]), sector=isector, window=15.)
        c1.SaveAs(pdir + 'segdif_dt13_slope.png')
        
        if istation[1] != '4':
          segdiff(tfiles_plotting, "dt2_resid", dstations, wheel=int(iwheel[1]), sector=isector, window=15.)
          c1.SaveAs(pdir + 'segdif_dt2_resid.png')
          segdiff(tfiles_plotting, "dt2_slope", dstations, wheel=int(iwheel[1]), sector=isector, window=15.)
          c1.SaveAs(pdir + 'segdif_dt2_slope.png')
        
  qcount = 0
  for iwheel in DT_TYPES:
    if iwheel[1]=="ALL": continue
    if abs(int(iwheel[1])) != 2: continue
    for istation in iwheel[2]:
      if istation[1]=="3": continue
      #print dstations
      for isector in range(1, istation[2] + 1):
        ssector = "%02d" % isector
        pdir = dt_basedir + '/' + iwheel[0] + '/' + istation[1] + '/' + ssector + '/'
        if istation[1]=="1":
          segdiff_xalign(tfiles_plotting, "x_dt1_csc", wheel=int(iwheel[1]), sector=isector, cscstations = "12")
          c1.SaveAs(pdir + 'segdif_x_dt_csc_resid.png')
        if istation[1]=="2":
          segdiff_xalign(tfiles_plotting, "x_dt2_csc", wheel=int(iwheel[1]), sector=isector, cscstations = "1")
          c1.SaveAs(pdir + 'segdif_x_dt_csc_resid.png')

  """segdiffvsphi "dt13_resid" and "dt13_slope"

  plot for a specific wheel #:
  x vs phi of pair ("dt13_resid")
  dxdz vs phi of pair ("dt13_slope")
  contains all three combinations of neighboring stations
  made for all possible wheel values

  Interface: could be accessed by clicking on wheel number under the "wheel" column
  in a general DT map

 segdiffvsphi "dt2_resid" and "dt2_slope"

  plot for a specific wheel #:
  y vs phi of pair ("dt2_resid")
  dydz vs phi of pair ("dt2_slope")
  contains both MB1-MB2 and MB2-MB3 combinations
  made for all possible wheel values

  Interface: could be accessed by clicking on wheel number under the "wheel" column
  in a general DT map"""

  if len(iter_reports)==0: return
  
  for iwheel in DT_TYPES:
    if iwheel[1]=="ALL": continue
    pdir = dt_basedir + '/' + iwheel[0] + '/'
    segdiffvsphi(tfiles_plotting, iter_reports, "dt13_resid", int(iwheel[1]), window=10.)#, excludesectors=(1,7))
    c1.SaveAs(pdir + 'segdifphi_dt13_resid.png')
    segdiffvsphi(tfiles_plotting, iter_reports, "dt13_slope", int(iwheel[1]), window=10.)#, excludesectors=(1,7))
    c1.SaveAs(pdir + 'segdifphi_dt13_slope.png')
    segdiffvsphi(tfiles_plotting, iter_reports, "dt2_resid", int(iwheel[1]), window=10.)#, excludesectors=(1,7))
    c1.SaveAs(pdir + 'segdifphi_dt2_resid.png')
    segdiffvsphi(tfiles_plotting, iter_reports, "dt2_slope", int(iwheel[1]), window=15.)#, excludesectors=(1,7))
    c1.SaveAs(pdir + 'segdifphi_dt2_slope.png')

  for iwheel in DT_TYPES:
    if iwheel[1]=="ALL": continue
    if abs(int(iwheel[1])) != 2: continue
    pdir = dt_basedir + '/' + iwheel[0] + '/'
    segdiffvsphi_xalign(tfiles_plotting, int(iwheel[1]), window=10.)
    c1.SaveAs(pdir + 'segdifphi_x_dt_csc_resid.png')


def doSegDiffPlotsCSC(csc_basedir, tfiles_plotting, iter_reports):
  """write segment-difference plots for CSC
 
 segdiff "csc_resid" and "csc_slope"

  set of plots of
  rphi vs qpt, rphi for positive, rphi for negative ("csc_resid")
  drphidz vs qpt, drphidz for positive, drphidz for negative ("csc_slope")
  done for ME1-ME2, ME2-ME3, and ME3-ME4 stations combinations with
    endcap "m" or "p" 
    ring 1 or 2
    chamber 1-18 (r1) or 1-36 (r2)
  note: there's no ME3-ME4 plots for R2
  
  Interface: could be accessed through a general CSC chambers map, but only for chambers in
  stations 2-4 (e.g., station 2 would provide ME1-ME2 plots)."""

  qcount = 0
  for iendcap in CSC_TYPES:
    for istation in iendcap[2]:
      if istation[1]=="1": continue
      dstations = (int(istation[1])-1)*10 + int(istation[1])
      for iring in istation[2]:
        if iring[1]=="ALL": continue
        if istation[1]=="4" and iring[1]=="2": continue
        for ichamber in range(1,iring[2]+1):
          if qcount>QUICKTESTN: break
          qcount += 1
          schamber = "%02d" % ichamber
          pdir = csc_basedir+'/'+iendcap[0]+'/'+istation[1]+'/'+iring[1]+'/'+schamber+'/'
          segdiff(tfiles_plotting, "csc_resid", dstations, 
                  endcap=iendcap[1], ring=int(iring[1]), chamber=ichamber, window=15.)
          c1.SaveAs(pdir + 'segdif_csc_resid.png')
          segdiff(tfiles_plotting, "csc_slope", dstations, 
                  endcap=iendcap[1], ring=int(iring[1]), chamber=ichamber, window=15.)
          c1.SaveAs(pdir + 'segdif_csc_slope.png')

  """segdiffvsphicsc "csc_resid" and "csc_slope"

  plot for a specific deltaME station differences:
  rphi vs phi of pair ("csc_resid")
  dxdz vs phi of pair ("csc_slope")
  contains plots for two (or one for ME4-ME3) rings
  done for ME1-ME2, ME2-ME3, and ME3-ME4 stations combinations with
    endcap "m" or "p" 
  
  Interface: could be accessed by clicking on ME station boxes, but only for stations 2-4 
  (e.g., station 2 would provide ME1-ME2 plots)."""

  qcount = 0
  for iendcap in CSC_TYPES:
    for istation in iendcap[2]:
      if istation[1]=="1": continue
      dstations = (int(istation[1])-1)*10 + int(istation[1])
      if qcount>QUICKTESTN: break
      qcount += 1
      pdir = csc_basedir+'/'+iendcap[0]+'/'+istation[1]+'/'
      segdiffvsphicsc(tfiles_plotting, "csc_resid", dstations, window=10., endcap=iendcap[1])
      c1.SaveAs(pdir + 'segdifphi_csc_resid.png')
      segdiffvsphicsc(tfiles_plotting, "csc_slope", dstations, window=10., endcap=iendcap[1])
      c1.SaveAs(pdir + 'segdifphi_csc_slope.png')


def doFitFunctionsPlotsDT(dt_basedir, iter_tfile, iter_reports):
  """write fit functions plots for DT

 DT bellcurves and polynomials

  set of plots of bellcurves
    x, dxdz, x vs. dxdz (for all 4 stations)
    y, dydz, x vs. dxdz (only for stations 1-3?)

  set of plots of polynomials -- for stations 1-3 only??
    x vs. xpos,    x vs ypos,    x vs dxdz angle,    x vs dydz angle
    y vs. xpos,    y vs ypos,    y vs dxdz angle,    y vs dydz angle
    dxdz vs. xpos, dxdz vs ypos, dxdz vs dxdz angle, dxdz vs dydz angle
    dydz vs. xpos, dydz vs ypos, dydz vs dxdz angle, dydz vs dydz angle

  set of plots of polynomials -- for station 4 only??
    x vs. xpos,    x vs dxdz angle
    dxdz vs. xpos, dxdz vs dxdz angle

  made for all (wheel,station,sector) combinations

  Interface: could be accesses through a general DT chambers map."""

  qcount = 0
  clearDDT()
  for wheel in DT_TYPES:
    if wheel[1]=="ALL": continue
    for station in wheel[2]:
      print wheel[0]+'/'+station[1]
      for sector in range(1,station[2]+1):
        if qcount>QUICKTESTN: break
        qcount += 1
        ssector = "%02d" % sector
        pdir = dt_basedir+'/'+wheel[0]+'/'+station[1]+'/'+ssector+'/'
        label = "MBwh%sst%ssec%s" % (wheelLetter(int(wheel[1])),station[1],ssector)
        bellcurves(iter_tfile, iter_reports, label, False)
        #c1.SaveAs(pdir+'dt_bellcurves.png')
        saveAs(pdir+'dt_bellcurves.png')
        polynomials(iter_tfile, iter_reports, label, False)
        #c1.SaveAs(pdir+'dt_polynomials.png')
        saveAs(pdir+'dt_polynomials.png')
  #printDeltaTs()


def doFitFunctionsPlotsCSC(csc_basedir, iter_tfile, iter_reports):
  """write fit functions plots for CSC

 CSC bellcurves and polynomials

  set of plots of bellcurves
    rphi, drphidz, rphi vs. drphidz

  set of plots of polynomials
    rphi vs. rphi pos,    rphi vs drphidz angle
    drphidz vs. rphi pos, drphidz vs drphidz angle

  made for all (endcap,station,ring,chamber) combinations

  Interface: could be accesses through a general CSC chambers map."""

  qcount = 0
  clearDDT()
  for endcap in CSC_TYPES:
    for station in endcap[2]:
      for ring in station[2]:
        if ring[1]=="ALL": continue
        print endcap[0]+'/'+station[1]+'/'+ring[1]
        for chamber in range(1,ring[2]+1):
          if qcount>QUICKTESTN: break
          qcount += 1
          schamber = "%02d" % chamber
          pdir = csc_basedir+'/'+endcap[0]+'/'+station[1]+'/'+ring[1]+'/'+schamber+'/'
          label = "ME%s%s%s_%s" % (endcap[1], station[1], ring[1], schamber)
          bellcurves(iter_tfile, iter_reports, label, False)
          #c1.SaveAs(pdir+'csc_bellcurves.png')
          saveAs(pdir+'csc_bellcurves.png')
          polynomials(iter_tfile, iter_reports, label, False)
          #c1.SaveAs(pdir+'csc_polynomials.png')
          saveAs(pdir+'csc_polynomials.png')
  #printDeltaTs()



def doIterationPlots(iteration_directory, tfiles_plotting, iter_tfile, iter_reports):
  dt_basedir = iteration_directory+'/'+'MB'
  csc_basedir = iteration_directory+'/'

  if DO_DT and DO_MAP:
    doMapPlotsDT(dt_basedir, tfiles_plotting)
  if DO_CSC and DO_MAP:
    doMapPlotsCSC(csc_basedir, tfiles_plotting)

  if DO_DT and DO_CURVATURE:
    doCurvaturePlotsDT(dt_basedir, tfiles_plotting)
  #if DO_CSC and DO_CURVATURE:
  #  doCurvaturePlotsCSC(csc_basedir, tfiles_plotting)

  if DO_DT and DO_SEGDIFF:
    doSegDiffPlotsDT(dt_basedir, tfiles_plotting, iter_reports)
  if DO_CSC and DO_SEGDIFF:
    doSegDiffPlotsCSC(csc_basedir, tfiles_plotting, iter_reports)

  if DO_DT and DO_FIT:
    doFitFunctionsPlotsDT(dt_basedir, iter_tfile, iter_reports)
  if DO_CSC and DO_FIT:
    doFitFunctionsPlotsCSC(csc_basedir, iter_tfile, iter_reports)


def createCanvasesList(fname="canvases_list.js"):
  '''Use CANVASES_LIST_TEMPLATE as a template to create a canvases list include for the browser.
     Write out only those canvases which have existing filename.png plots.
  '''
  CANVASES_LIST = []
  for scope in CANVASES_LIST_TEMPLATE:
    scope_entry = []
    if len(scope)>2:
      scope_entry = [scope[0],scope[1]]
      for canvas_entry in scope[2:]:
        if isFileUnderDir("./",canvas_entry[1]):
          scope_entry.append(canvas_entry)
    CANVASES_LIST.append(scope_entry)

  ff = open(fname,mode="w")
  print >>ff, "var CANVASES_LIST = "
  json.dump(CANVASES_LIST,ff)
  ff.close()


def createCanvasToIDList(fname="canvas2id_list.js"):
  '''Writes out a canvas-2-ids list include for the browser.
     Write out only those canvases which have existing filename.png plots.
     Returns: list of unique IDs that have existing filename.png plots.
  '''
  CANVAS2ID_LIST = []
  ID_LIST = []
  for scope in CANVASES_LIST_TEMPLATE:
    if len(scope)>2:
      for canvas_entry in scope[2:]:
        ids = idsForFile("./",canvas_entry[1])
        #  scope_entry.append(canvas_entry)
        # uniquify:
        set_ids = set(ids)
        uids = list(set_ids)
        ID_LIST.extend(uids)
        print canvas_entry, ":", len(uids), "ids"
        if (len(uids)>0):
          CANVAS2ID_LIST.append( (canvas_entry[1],uids) )
  #print CANVAS2ID_LIST
  CANVAS2ID_LIST_DICT = dict(CANVAS2ID_LIST)
  #print CANVAS2ID_LIST_DICT
  ff = open(fname,mode="w")
  print >>ff, "var CANVAS2ID_LIST = "
  json.dump(CANVAS2ID_LIST_DICT,ff)
  ff.close()
  set_ids = set(ID_LIST)
  return list(set_ids)

def idsForFile(dir_name, file_name):
  '''Recursively looks for file named file_name under dir_name directory
  and fill the list with dir names converted to IDs 
  '''
  id_list = []
  for f in os.listdir(dir_name):
    dirfile = os.path.join(dir_name, f)
    if os.path.isfile(dirfile) and f==file_name:
      if file_name[-4:]=='.php': id_list.append(dir_name+'/'+file_name)
      else:    id_list.append(dirToID(dir_name))
    # recursively access file names in subdirectories
    elif os.path.isdir(dirfile):
      #print "Accessing directory:", dirfile
      ids = idsForFile(dirfile, file_name)
      if (len(ids)>0): 
        id_list.extend(ids)
  return id_list


def dirToID(d):
  if d[-1]!='/': d += '/'
  dtn = d.find("/MB/")
  if dtn!=-1:
    return d[dtn+4:-1]
  cscn = d.find("/ME-/")
  if cscn!=-1:
    return 'ME-'+d[cscn+5:-1]
  cscn = d.find("/ME+/")
  if cscn!=-1:
    return 'ME+'+d[cscn+5:-1]
  return ''


############################################################################################################
############################################################################################################
# main script

# open input files:

#basedir='/disks/sdb5/home_reloc/khotilov/db/cms/alignment'
#iteration1 = "iteration_01"
#iteration3 = "iteration_03"
#iteration1 = "NOV4DT_PASS3noweight_TkHIP_01"
#iteration3 = "NOV4DT_PASS3noweight_TkHIP_05"
#os.chdir(options.inputDir)
#iteration1 = options.i1
#iteration3 = options.iN

fname = options.inputDir+'/'+options.i1+'/'+i1prefix
tfiles1_plotting = []
iter1_tfile = None
iter1_reports = []
if not SINGLE_ITERATION:
  if (DO_MAP or DO_SEGDIFF or DO_CURVATURE) and not os.access(fname+"_plotting.root",os.F_OK):
    print "no file "+fname+"_plotting.root"
    sys.exit()
  if DO_FIT and not os.access(fname+".root",os.F_OK):
    print "no file "+fname+".root"
    sys.exit()
  if DO_MEDIAN and not os.access(fname+"_report.py",os.F_OK):
    print "no file "+fname+"_report.py"
    sys.exit()
  if DO_MAP or DO_SEGDIFF or DO_CURVATURE: tfiles1_plotting.append(ROOT.TFile(fname+"_plotting.root"))
  if os.access(fname+".root",os.F_OK):
    iter1_tfile = ROOT.TFile(fname+".root")
  if os.access(fname+"_report.py",os.F_OK):
    execfile(fname+"_report.py")
    iter1_reports = reports

fname = options.inputDir+'/'+options.iN+'/'+iNprefix
tfilesN_plotting = []
iterN_tfile = None
iterN_reports = []
if (DO_MAP or DO_SEGDIFF or DO_CURVATURE) and not os.access(fname+"_plotting.root",os.F_OK):
  print "no file "+fname+"_plotting.root"
  sys.exit()
if DO_FIT and not os.access(fname+".root",os.F_OK):
  print "no file "+fname+".root"
  sys.exit()
if DO_MEDIAN and not os.access(fname+"_report.py",os.F_OK):
  print "no file "+fname+"_report.py"
  sys.exit()
if DO_MAP or DO_SEGDIFF or DO_CURVATURE: tfilesN_plotting.append(ROOT.TFile(fname+"_plotting.root"))
if os.access(fname+".root",os.F_OK):
  iterN_tfile = ROOT.TFile(fname+".root")
if os.access(fname+"_report.py",os.F_OK):
  execfile(fname+"_report.py")
  iterN_reports = reports

if DO_MAP:
  os.chdir(options.inputDir)
  phiedges2c()

######################################################
# setup output:

# cd to outputDIr
os.chdir(options.outputDir)

comdir = "common/"
iteration1 = "iter1"
iterationN = "iterN"

# create directory structure
if options.createDirSructure:
  print "WARNING: all existing results in "+options.outputDir+" will be deleted!"
  if not SINGLE_ITERATION: createDirectoryStructure(iteration1)
  createDirectoryStructure(iterationN)
  if not os.access(comdir,os.F_OK):
    os.mkdir(comdir)

######################################################
# do drawing

c1 = ROOT.TCanvas("c1","c1",800,600)

set_palette("blues")

print "--- ITERATION 1 ---"
if not SINGLE_ITERATION: doIterationPlots(iteration1, tfiles1_plotting, iter1_tfile, iter1_reports)
print "--- ITERATION N ---"
doIterationPlots(iterationN, tfilesN_plotting, iterN_tfile, iterN_reports)

if CPP_LOADED: ROOT.cleanUpHeap()

# write distributions of medians plots

if DO_MEDIAN:
  #plotmedians(iter1_reports, iterN_reports)
  plotmedians(iter1_reports, iterN_reports,binsx=100, windowx=10., binsy=100, windowy=10., binsdxdz=100, windowdxdz=10., binsdydz=100, windowdydz=10.)
  c1.SaveAs(comdir+'medians.png')

# perform diagnostic
if DO_DIAGNOSTIC:
  #if not SINGLE_ITERATION: doTests(iter1_reports,"mu_list_1.js","dqm_report_1.js",options.runLabel)
  createCanvasesList("canvases_list.js")
  pic_ids = createCanvasToIDList("canvas2id_list.js")
  doTests(iterN_reports, pic_ids, "mu_list.js","dqm_report.js",options.runLabel)