MuonAlignment/python/MCScenario_CRAFT1_22X.py

0001 # To use this script:
0002 #     it's an ordinary Python script, not a CMSSW configuration file (though it writes and runs CMSSW configuration files)
0003 #         * you MUST have an inertGlobalPositionRcd.db in your working directory
0004 #         * you MUST NOT have an MCScenario_CRAFT1_22X.db
0005 #
0006 #     to make the inertGlobalPositionRcd:
0007 #         cmsRun Alignment/MuonAlignment/python/makeGlobalPositionRcd_cfg.py
0008 #
0009 #     to get rid of the MCScenario_CRAFT1_22X.db:
0010 #         rm MCScenario_CRAFT1_22X.db    (naturally)
0011 #
0012 #     to run this script:
0013 #         python MCScenario_CRAFT1_22X.py
0014 #
0015 #     it will create
0016 #         * MCScenario_CRAFT1_22X.xml            the XML file with randomly-distributed values, created directly by define_scenario()
0017 #         * convert_cfg.py                       the conversion configuration file
0018 #         * MCScenario_CRAFT1_22X.db             the SQLite database created from the XML
0019 #         * check_cfg.py                         configuration file that converts the SQLite file back into XML
0020 #         * MCScenario_CRAFT1_22X_CHECKME.xml    converted back, so that we can check the values that were saved to the database
0021 #
0022 #     to check the output in Excel, do this
0023 #         ./Alignment/MuonAlignment/python/geometryXMLtoCSV.py < MCScenario_CRAFT1_22X_CHECKME.xml > MCScenario_CRAFT1_22X_CHECKME.csv
0024 #         and then open MCScenario_CRAFT1_22X_CHECKME.csv in Excel
0025
0026 from builtins import range
0027 import random, os
0028 from math import *
0029
0030 # set the initial seed for reproducibility!
0031 random.seed(123456)
0032
0033 #### called once at the end of this script
0034 def make_scenario_sqlite():
0035     scenario = define_scenario()
0036     write_xml(scenario, "MCScenario_CRAFT1_22X.xml")
0037     write_conversion_cfg("convert_cfg.py", "MCScenario_CRAFT1_22X.xml", "MCScenario_CRAFT1_22X.db")
0038     cmsRun("convert_cfg.py")
0039     write_check_cfg("check_cfg.py", "MCScenario_CRAFT1_22X.db", "MCScenario_CRAFT1_22X_CHECKME.xml")
0040     cmsRun("check_cfg.py")
0041 #### that's it!  everything this uses is defined below
0042
0043 def write_conversion_cfg(fileName, xmlFileName, dbFileName):
0044     outfile = file(fileName, "w")
0045     outfile.write("""
0046 from Alignment.MuonAlignment.convertXMLtoSQLite_cfg import *
0047 process.MuonGeometryDBConverter.fileName = "%(xmlFileName)s"
0048 process.PoolDBOutputService.connect = "sqlite_file:%(dbFileName)s"
0049 """ % vars())
0050
0051 def write_check_cfg(fileName, dbFileName, xmlFileName):
0052     outfile = file(fileName, "w")
0053     outfile.write("""
0054 from Alignment.MuonAlignment.convertSQLitetoXML_cfg import *
0055 process.PoolDBESSource.connect = "sqlite_file:%(dbFileName)s"
0056 process.MuonGeometryDBConverter.outputXML.fileName = "%(xmlFileName)s"
0057 process.MuonGeometryDBConverter.outputXML.relativeto = "ideal"
0058 process.MuonGeometryDBConverter.outputXML.suppressDTChambers = False
0059 process.MuonGeometryDBConverter.outputXML.suppressDTSuperLayers = False
0060 process.MuonGeometryDBConverter.outputXML.suppressDTLayers = True
0061 process.MuonGeometryDBConverter.outputXML.suppressCSCChambers = False
0062 process.MuonGeometryDBConverter.outputXML.suppressCSCLayers = False
0063 """ % vars())
0064
0065 def cmsRun(fileName):
0066     os.system("cmsRun %(fileName)s" % vars())
0067
0068 ########### writing a scenario in XML ##############################################################
0069
0070 # only needed to make the output XML readable
0071 DTpreferred_order = {"wheel":1, "station":2, "sector":3, "superlayer":4, "layer":5}
0072 CSCpreferred_order = {"endcap":1, "station":2, "ring":3, "chamber":4, "layer":5}
0073 def DTsorter(a, b): return cmp(DTpreferred_order[a], DTpreferred_order[b])
0074 def CSCsorter(a, b): return cmp(CSCpreferred_order[a], CSCpreferred_order[b])
0075
0076 # an instance of this class corresponds to one <DTChamber ... /> or <CSCStation ... />, etc.
0077 class Alignable:
0078     def __init__(self, alignabletype, **location):
0079         self.alignabletype = alignabletype
0080         self.location = location
0081
0082     def writeXML(self):
0083         parameters = self.location.keys()
0084         if self.alignabletype[0:2] == "DT":
0085             parameters.sort(DTsorter)
0086         else:
0087             parameters.sort(CSCsorter)
0088
0089         output = ["<", self.alignabletype, " "]
0090         for parameter in parameters:
0091             output.extend([parameter, "=\"", str(self.location[parameter]), "\" "])
0092         output.append("/>")
0093
0094         return "".join(output)
0095
0096 preferred_order = {"x":1, "y":2, "z":3, "phix":4, "phiy":5, "phiz":6}
0097 def sorter(a, b): return cmp(preferred_order[a], preferred_order[b])
0098
0099 # an instance of this class corresponds to one <setposition ... />
0100 class Position:
0101     def __init__(self, **location):
0102         self.location = location
0103
0104     def writeXML(self):
0105         parameters = self.location.keys()
0106         parameters.sort(sorter)
0107
0108         output = ["<setposition relativeto=\"ideal\" "]
0109         for parameter in parameters:
0110             output.extend([parameter, "=\"", str(self.location[parameter]), "\" "])
0111         output.append("/>")
0112
0113         return "".join(output)
0114
0115 # an instance of this class corresponds to one <operation> ... </operation> in the XML file
0116 class Operation:
0117     def __init__(self, alignable, position):
0118         self.alignable = alignable
0119         self.position = position
0120
0121     def writeXML(self):
0122         output = ["<operation> ", self.alignable.writeXML(), " ", self.position.writeXML(), " </operation>\n"]
0123         return "".join(output)
0124
0125 def write_xml(scenario, fileName):
0126     # a scenario is an ordered list of Operations
0127     XMLlist = ["<MuonAlignment>\n"]
0128     for operation in scenario:
0129         XMLlist.append(operation.writeXML())
0130     XMLlist.append("</MuonAlignment>\n")
0131     XMLstring = "".join(XMLlist)
0132
0133     outfile = file(fileName, "w")
0134     outfile.write(XMLstring)
0135
0136 class DTChamber:
0137     def __init__(self, **location):
0138         self.__dict__.update(location)
0139
0140 class CSCChamber:
0141     def __init__(self, **location):
0142         self.__dict__.update(location)
0143
0144 ########### defining the actual scenario ##############################################################
0145
0146 # this is the interesting part: where we define a scenario for CRAFT1 MC
0147 def define_scenario():
0148     # this will be a list of operations to write to an XML file
0149     scenario = []
0150
0151     # Uncertainty in DT chamber positions comes in two parts:
0152     #    1. positions within sectors
0153     #    2. positions of the sector-groups
0154
0155     # Aligned chambers (wheels -1, 0, +1 except sectors 1 and 7)
0156     # uncertainty within sectors:
0157     #    x: 0.08 cm (from segment-matching)  phix: 0.0007 rad (from MC)
0158     #    y: 0.10 cm (from MC)                phiy: 0.0007 rad (from segment-matching)
0159     #    z: 0.10 cm (from MC)                phiz: 0.0003 rad (from MC)
0160     # uncertainty of sector-groups (depends on choice of pT cut, not well understood):
0161     #    x: 0.05 cm
0162
0163     # Unaligned chambers uncertainty within sectors:
0164     #    x: 0.08 cm (same as above)          phix: 0.0016 rad
0165     #    y: 0.24 cm                          phiy: 0.0021 rad
0166     #    z: 0.42 cm with a -0.35 cm bias     phiz: 0.0010 rad
0167     # uncertainty of sector-groups:
0168     #    x: 0.65 cm
0169     # These come from actual alignments measured in the aligned
0170     # chambers (we assume that the unaligned chambers have
0171     # misalignments on the same scale)
0172
0173     # Also, superlayer z uncertainty is 0.054 cm
0174
0175     # Before starting, let's build a list of chambers
0176     DTchambers = []
0177     for wheel in -2, -1, 0, 1, 2:
0178         for station in 1, 2, 3, 4:
0179             if station == 4: nsectors = 14
0180             else: nsectors = 12
0181             for sector in range(1, nsectors+1):
0182                 DTchambers.append(DTChamber(wheel = wheel, station = station, sector = sector))
0183
0184     # the superlayers
0185     for dtchamber in DTchambers:
0186         for superlayer in 1, 2, 3:
0187             if superlayer == 2 and dtchamber.station == 4: continue
0188
0189             alignable = Alignable("DTSuperLayer", wheel = dtchamber.wheel, station = dtchamber.station, sector = dtchamber.sector, superlayer = superlayer)
0190             position = Position(x = 0, y = 0, z = random.gauss(0, 0.054), phix = 0, phiy = 0, phiz = 0)
0191             scenario.append(Operation(alignable, position))
0192
0193     sector_errx = {}
0194
0195     # sector-groups for aligned chambers:
0196     for wheel in -1, 0, 1:
0197         for sector in 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14:
0198             sector_errx[wheel, sector] = random.gauss(0., 0.05)
0199
0200     # sector-groups for unaligned chambers:
0201     for wheel in -1, 0, 1:
0202         for sector in 1, 7:
0203             sector_errx[wheel, sector] = random.gauss(0., 0.65)
0204     for wheel in -2, 2:
0205         for sector in 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14:
0206             sector_errx[wheel, sector] = random.gauss(0., 0.65)
0207
0208     for dtchamber in DTchambers:
0209         # within sectors for aligned chambers:
0210         if dtchamber.wheel in (-1, 0, 1) and dtchamber.sector in (2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14):
0211             errx = random.gauss(0, 0.08)
0212             erry = random.gauss(0, 0.10)
0213             errz = random.gauss(0, 0.10)
0214             errphix = random.gauss(0, 0.0007)
0215             errphiy = random.gauss(0, 0.0007)
0216             errphiz = random.gauss(0, 0.0003)
0217
0218         # within sectors for unaligned chambers:
0219         else:
0220             errx = random.gauss(0, 0.08)
0221             erry = random.gauss(0, 0.24)
0222             errz = random.gauss(-0.35, 0.42)
0223             errphix = random.gauss(0, 0.0016)
0224             errphiy = random.gauss(0, 0.0021)
0225             errphiz = random.gauss(0, 0.0010)
0226
0227         errx += sector_errx[dtchamber.wheel, dtchamber.sector]
0228
0229         # now turn this into an operation
0230         alignable = Alignable("DTChamber", wheel = dtchamber.wheel, station = dtchamber.station, sector = dtchamber.sector)
0231         position = Position(x = errx, y = erry, z = errz, phix = errphix, phiy = errphiy, phiz = errphiz)
0232         scenario.append(Operation(alignable, position))
0233
0234     # Uncertainty in CSC chamber positions comes in 5 parts:
0235     #    1. 0.0092 cm layer x misalignments observed with beam-halo tracks
0236     #    2. isotropic photogrammetry uncertainty of 0.03 cm (x, y, z) and 0.00015 rad in phiz
0237     #    3. 0.0023 rad phiy misalignment observed with beam-halo tracks
0238     #    4. 0.1438 cm z and 0.00057 rad phix uncertainty between rings from SLM (from comparison in 0T data with PG)
0239     #    5. 0.05 cm (x, y, z) disk misalignments and 0.0001 rad rotation around beamline
0240
0241     # Before starting, let's build a list of chambers
0242     CSCchambers = []
0243     for endcap in 1, 2:
0244         for station, ring in (1, 1), (1, 2), (1, 3), (1, 4), (2, 1), (2, 2), (3, 1), (3, 2), (4, 1):
0245             if station > 1 and ring == 1:
0246                 nchambers = 18
0247             else:
0248                 nchambers = 36
0249
0250             for chamber in range(1, nchambers+1):
0251                 CSCchambers.append(CSCChamber(endcap = endcap, station = station, ring = ring, chamber = chamber))
0252
0253     # First, the layer uncertainties: x only for simplicity, observed 0.0092 cm in overlaps alignment test
0254     for chamber in CSCchambers:
0255         for layer in 1, 2, 3, 4, 5, 6:
0256             alignable = Alignable("CSCLayer", endcap = chamber.endcap, station = chamber.station, ring = chamber.ring, chamber = chamber.chamber, layer = layer)
0257             position = Position(x = random.gauss(0, 0.0092), y = 0, z = 0, phix = 0, phiy = 0, phiz = 0)
0258             scenario.append(Operation(alignable, position))
0259
0260     # Next, the ring errors from DCOPS (derived from comparison with photogrammetry)
0261     CSCrings = []
0262     for endcap in 1, 2:
0263         for station, ring in (1, 1), (1, 2), (1, 3), (1, 4), (2, 1), (2, 2), (3, 1), (3, 2), (4, 1):
0264             CSCrings.append(CSCChamber(endcap = endcap, station = station, ring = ring, z = random.gauss(0, 0.1438), phix = random.gauss(0, 0.00057)))
0265
0266     # Next, the chamber errors
0267     for chamber in CSCchambers:
0268         errx = random.gauss(0, 0.03)
0269         erry = random.gauss(0, 0.03)
0270         errz = random.gauss(0, 0.03)
0271         errphix = random.gauss(0, 0.00057)
0272         errphiy = random.gauss(0, 0.0023)
0273         errphiz = random.gauss(0, 0.00015)
0274
0275         for ring in CSCrings:
0276             if ring.endcap == chamber.endcap and ring.station == chamber.station and ring.ring == chamber.ring:
0277                 errz += ring.z
0278                 errphix += ring.phix
0279                 break
0280
0281         alignable = Alignable("CSCChamber", endcap = chamber.endcap, station = chamber.station, ring = chamber.ring, chamber = chamber.chamber)
0282         position = Position(x = errx, y = erry, z = errz, phix = errphix, phiy = errphiy, phiz = errphiz)
0283         scenario.append(Operation(alignable, position))
0284
0285     # Finally, the disk errors
0286     for endcap in 1, 2:
0287         for station in 1, 2, 3, 4:
0288             alignable = Alignable("CSCStation", endcap = endcap, station = station)
0289             position = Position(x = random.gauss(0, 0.05), y = random.gauss(0, 0.05), z = random.gauss(0, 0.05), phix = 0., phiy = 0., phiz = random.gauss(0, 0.0001))
0290             scenario.append(Operation(alignable, position))
0291
0292     return scenario
0293
0294 # run it all!
0295 make_scenario_sqlite()