File indexing completed on 2024-04-06 12:15:26
0001
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))
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
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)
