Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:15:26

0001 #!/usr/bin/env python3
0002 
0003 """
0004 mixture.py replaces the functionality of mixture.f
0005 
0006 Usage:
0007 
0008 $ python3 mixture.py filename.in
0009 
0010 Where `filename.in` contains a set of mixtures and its components:
0011 
0012 The format for `filename.in` is as follows:
0013 
0014 Each mixture description starts with a hash symbol (#) and contains
0015 the following information in a single line:
0016 
0017  * Name (name of mixture of components)
0018  * GMIX name ()
0019  * MC Volume
0020  * MC Area
0021 
0022 A mixture definition is followed by a set of compounds, which is the
0023 set of components that the mixture is made from.
0024 
0025 Each compound description starts with an asterisk (*) and contains
0026 the following information:
0027 
0028  * Item number (An index for the list)
0029  * Comment (A human readable name for the material)
0030  * Material (Component name)
0031  * Volume
0032  * Multiplicity (How many times the material is in the mixture)
0033  * Type (see below)
0034 
0035 Type stands for a classification of the Material using the following
0036 convention:
0037 
0038  * SUP: Support
0039  * COL: Cooling
0040  * ELE: Electronics
0041  * CAB: Cables
0042  * SEN: Sensitive volumes
0043 
0044 Each Material (component of the mixture) belongs to a single category.
0045 `mixture.py` uses these information to compute how the mixture as a whole
0046 contributes to each of the categories in terms of radiation length (x0)
0047 and nuclear interaction length (l0) \lambda_{0}.
0048 
0049 The radiation length or the nuclear interaction lenght become relevant
0050 depending if we are talking about particles that interact with the detector
0051 mainly throgh electromagnetic forces (x0) or either strong interaction or nuclear
0052 forces (l0). Notice that these quantities are different and therefore the
0053 contribution for each category will depend on them.
0054 
0055 Any line in `filename.in` is therefore ignored unless it begins with
0056 a hash or an asterisk.
0057 
0058 Hardcoded in `mixture.py` are the files `mixed_materials.input` and
0059 `pure_materials.input`. These files act as a database that provides
0060 information about the compounds. From each compound defined in
0061 `filename.in` `mixture.py` matches the content of `Material` with a
0062 single element of the database and extracts the following information
0063 (the format of the mentioned .input files):
0064 
0065  * Name (which should match with Material)
0066  * Weight (Atomic Mass) [g/mol]
0067  * Number (Atomic number)
0068  * Density (g/cm^3)
0069  * RadLen (Radiation lenght) (cm)
0070  * IntLen (Nuclear interaction lenght) (cm)
0071 
0072 With these information `material.py` computes for each mixture:
0073 
0074  * The radiation lenght for the mixture
0075  * The nuclear interaction lenght for the mixture
0076  * The mixture density
0077  * The mixture volume
0078  * The total weight
0079  * The relative contribution in each categorio on x0 and l0
0080  * Normalized quantities based on MC Volume and MC Area values
0081 
0082 As well for the compounds:
0083 
0084  * The relative volume within the mixture
0085  * The relative weight
0086  * The relative radiation lenght
0087  * The relative nuclear interaction length
0088 
0089 """
0090 
0091 import re
0092 import pandas as pd
0093 import sys
0094 from os.path import expandvars
0095 
0096 def getMixture(line) :
0097     mixture = re.search(
0098         r'^#\s*"(?P<MixtureName>[^"]+)"\s+"(?P<GMIXName>[^"]+)"\s+(?P<MCVolume>[0-9.Ee\-+]+)\s+(?P<MCArea>[0-9.Ee\-+]+)',
0099         line
0100     )
0101     return {
0102         'mixture_name' : mixture.group("MixtureName"),
0103         'gmix_name' : mixture.group("GMIXName"),
0104         'mc_volume' : float(mixture.group("MCVolume")),
0105         'mc_area' : float(mixture.group("MCArea")),
0106     }
0107 
0108 
0109 def getCompound(line) :
0110     compound = re.search(
0111         r'^\*\s*(?P<Index>[0-9]+)\s+"(?P<Comment>[^"]+)"\s+"(?P<Material>[^"]+)"\s+(?P<Volume>[0-9.Ee\-+]+)\s+(?P<Mult>[0-9.]+)\s+(?P<Type>[^ ]{3})',
0112         line
0113     )
0114     return {
0115         'item_number' : int(compound.group("Index")),
0116         'comment' : compound.group("Comment"),
0117         'material' : compound.group("Material"),
0118         'volume' : float(compound.group("Volume")),
0119         'multiplicity' : float(compound.group("Mult")),
0120         'category' : compound.group("Type"),
0121     }
0122 
0123 def getMixtures(inputFile) :
0124 
0125     """
0126     Returns a `dict` of Mixtures. Each
0127     mixture contains a list of compounds.
0128     """
0129 
0130     inFile = open(inputFile,"r")
0131 
0132     mixtures = []
0133     mixture = {}
0134 
0135     counter = 0
0136     prevCounter = 0
0137 
0138     for line in inFile.readlines():
0139 
0140         if line[0] == "#":
0141             if mixture:
0142                 mixtures.append(dict(mixture)) #Just to make a copy
0143             mixture.clear()
0144             mixture = getMixture(line)
0145             mixture.update({ 'components' : [] })
0146 
0147         if line[0] == "*":
0148             mixture['components'].append(dict(getCompound(line)))
0149 
0150     inFile.close()
0151 
0152     return mixtures
0153 
0154 def loadMaterialsFile(inputFile):
0155 
0156     mixMatFile = open(inputFile,"r")
0157 
0158     mixMats = []
0159 
0160     for line in mixMatFile.readlines():
0161 
0162         if line[0] == '"':
0163             mixMat = re.search(
0164                 r'^"(?P<name>[^"]+)"\s+(?P<weight>[0-9.Ee-]+)\s+(?P<number>[0-9.Ee-]+)\s+(?P<density>[0-9.Ee-]+)\s+(?P<x0>[0-9.Ee-]+)\s+(?P<l0>[0-9Ee.-]+)',
0165                 line
0166                 )
0167 
0168             mixMats.append({
0169                     "name" : mixMat.group("name"),
0170                     "weight" : float(mixMat.group("weight")),
0171                     "number" : float(mixMat.group("number")),
0172                     "density" : float(mixMat.group("density")),
0173                     "x0" : float(mixMat.group("x0")),
0174                     "l0" : float(mixMat.group("l0")),
0175                     })
0176 
0177     return mixMats
0178 
0179 def main(argv):
0180 
0181     listOfMixtures = getMixtures(str(argv[1]))
0182 
0183     fname = re.search("(?P<filename>[a-zA-Z0-9_]+)",str(argv[1]))
0184 
0185     fname = str(fname.group("filename"))
0186     radLenFile = open(fname + ".x0","w");
0187     intLenFile = open(fname + ".l0","w");
0188 
0189     listOfMixedMaterials = loadMaterialsFile(
0190         inputFile = expandvars(
0191             "$CMSSW_BASE/src/Geometry/TrackerCommonData/data/Materials/mixed_materials.input"
0192         )
0193     )
0194     listOfPureMaterials = loadMaterialsFile(
0195         inputFile = expandvars(
0196             "$CMSSW_BASE/src/Geometry/TrackerCommonData/data/Materials/pure_materials.input"
0197         )
0198     )
0199     listOfMaterials = listOfMixedMaterials + listOfPureMaterials
0200 
0201     dfAllMaterials = pd.DataFrame(listOfMaterials)
0202 
0203     for index in range(len(listOfMixtures)):
0204 
0205         gmixName = listOfMixtures[index]["gmix_name"]
0206         print("================================")
0207         print(gmixName)
0208 
0209         components = pd.DataFrame(listOfMixtures[index]['components'])
0210 
0211         components["volume"] = components["multiplicity"]*components["volume"]
0212         totalVolume = sum(
0213             components["volume"]
0214             )
0215 
0216         components["percentVolume"] = components["volume"] / totalVolume
0217 
0218         components = pd.merge(
0219             components,
0220             dfAllMaterials[["name","density","x0","l0"]],
0221             left_on="material",
0222             right_on="name"
0223             )
0224 
0225         components["weight"] = components["density"] * components["volume"]
0226 
0227         totalDensity = sum(
0228             components["percentVolume"] * components["density"]
0229             )
0230 
0231         totalWeight = totalDensity * totalVolume
0232 
0233         components["percentWeight"] = (
0234             components["density"] * components["volume"] / totalWeight
0235             )
0236 
0237         components["ws"] = (
0238             components["percentWeight"] / (components["density"]*components["x0"])
0239             )
0240 
0241         components["ws2"] = (
0242             components["percentWeight"] / (components["density"]*components["l0"])
0243             )
0244 
0245         totalRadLen = 1 / ( sum(components["ws"]) * totalDensity )
0246 
0247         totalIntLen = 1/ ( sum(components["ws2"]) * totalDensity )
0248 
0249         components["percentRadLen"] = (
0250             components["ws"] * totalRadLen * totalDensity
0251             )
0252 
0253         components["percentIntLen"] = (
0254             components["ws2"] * totalIntLen * totalDensity
0255             )
0256 
0257         displayDf = components[[
0258                     "comment",
0259                     "material","volume",
0260                     "percentVolume","percentWeight",
0261                     "density","weight","x0","percentRadLen",
0262                     "l0","percentIntLen"
0263                     ]].copy()
0264 
0265         displayDf["percentVolume"] *= 100.
0266         displayDf["percentWeight"] *= 100.
0267         displayDf["percentRadLen"] *= 100.
0268         displayDf["percentIntLen"] *= 100.
0269 
0270         displayDf = displayDf.rename(
0271             columns = {
0272                 "comment" : "Component",
0273                 "material" : "Material",
0274                 "volume" : "Volume [cm^3]",
0275                 "percentVolume" : "Volume %",
0276                 "density" : "Density",
0277                 "weight" : "Weight",
0278                 "percentWeight" : "Weight %",
0279                 "x0" : "X_0 [cm]",
0280                 "percentRadLen" : "X_0 %",
0281                 "l0": "lambda_0 [cm]",
0282                 "percentIntLen" : "lambda_0 %",
0283                 }
0284             )
0285 
0286         print(displayDf)
0287 
0288 
0289         # Normalized vars:
0290 
0291         mcVolume = listOfMixtures[index]["mc_volume"]
0292         mcArea = listOfMixtures[index]["mc_area"]
0293         normalizationFactor = -1.
0294 
0295         if mcVolume > 0:
0296             normalizationFactor = totalVolume / mcVolume
0297         else:
0298             normalizationFactor = 1.
0299 
0300         normalizedDensity = totalDensity * normalizationFactor
0301         normalizedRadLen = totalRadLen / normalizationFactor
0302         normalizedIntLen = totalIntLen / normalizationFactor
0303 
0304         percentRadL = -1.
0305         percentIntL = -1.
0306 
0307         if mcArea > 0:
0308             percentRadL = mcVolume / (mcArea*normalizedRadLen)
0309             percentIntL = mcVolume / (mcArea*normalizedIntLen)
0310 
0311         pSupRadLen = 0.
0312         pSenRadLen = 0.
0313         pCabRadLen = 0.
0314         pColRadLen = 0.
0315         pEleRadLen = 0.
0316 
0317         pSupIntLen = 0.
0318         pSenIntLen = 0.
0319         pCabIntLen = 0.
0320         pColIntLen = 0.
0321         pEleIntLen = 0.
0322 
0323         for index, component in components.iterrows():
0324             catg = component["category"]
0325             prl = component["percentRadLen"]
0326             irl = component["percentIntLen"]
0327             if catg.upper() == "SUP":
0328                 pSupRadLen += prl
0329                 pSupIntLen += irl
0330             elif catg.upper() == "SEN":
0331                 pSenRadLen += prl
0332                 pSenIntLen += irl
0333             elif catg.upper() == "CAB":
0334                 pCabRadLen += prl
0335                 pCabIntLen += irl
0336             elif catg.upper() == "COL":
0337                 pColRadLen += prl
0338                 pColIntLen += irl
0339             elif catg.upper() == "ELE":
0340                 pEleRadLen += prl
0341                 pEleIntLen += irl
0342 
0343         print("================================")
0344         print("Mixture Density       [g/cm^3]: " , totalDensity)
0345         print("Norm. mixture density [g/cm^3]: ", normalizedDensity)
0346         print("Mixture Volume        [cm^3]: ", totalVolume)
0347         print("MC Volume             [cm^3]: ", mcVolume)
0348         print("MC Area               [cm^2]: ", mcArea)
0349         print("Normalization Factor:       ",  normalizationFactor)
0350         print("Mixture x0            [cm]: ", totalRadLen)
0351         print("Norm. Mixture x0      [cm]: ", normalizedRadLen)
0352         if mcArea > 0 :
0353             print("Norm. Mixture x0      [%]: ", 100. * percentRadL)
0354         print("Mixture l0            [cm]: ", totalIntLen)
0355         print("Norm. Mixture l0      [cm]: ", normalizedIntLen)
0356         if mcArea > 0 :
0357             print("Norm. Mixture l0      [%]: ", 100. * percentIntL)
0358         print("Total Weight          [g]: ", totalWeight)
0359 
0360         print("================================")
0361         print("X0 Contribution: ")
0362         print("Support     :", pSupRadLen)
0363         print("Sensitive   :", pSenRadLen)
0364         print("Cables      :", pCabRadLen)
0365         print("Cooling     :", pColRadLen)
0366         print("Electronics :", pEleRadLen)
0367 
0368         print("================================")
0369         print("l0 Contribution: ")
0370         print("Support     :", pSupIntLen)
0371         print("Sensitive   :", pSenIntLen)
0372         print("Cables      :", pCabIntLen)
0373         print("Cooling     :", pColIntLen)
0374         print("Electronics :", pEleIntLen)
0375 
0376         radLenFile.write("{0:<50}{1:>8}{2:>8}{3:>8}{4:>8}{5:>8}\n".format(
0377                 gmixName,
0378                 round(pSupRadLen,3), round(pSenRadLen,3),
0379                 round(pCabRadLen,3), round(pColRadLen,3),
0380                 round(pEleRadLen,3)
0381                 ))
0382 
0383         intLenFile.write("{0:<50}{1:>8}{2:>8}{3:>8}{4:>8}{5:>8}\n".format(
0384                 gmixName,
0385                 round(pSupIntLen,3), round(pSenIntLen,3),
0386                 round(pCabIntLen,3), round(pColIntLen,3),
0387                 round(pEleIntLen,3)
0388                 ))
0389 
0390     radLenFile.close()
0391     intLenFile.close()
0392 
0393 if __name__== "__main__":
0394     main(sys.argv)