s3tools/src/tabulate-cross-section.py

0001 #!/usr/bin/env python3
0002
0003 # SPDX-License-Identifier: LGPL-3.0-or-later
0004 # Copyright (C) 2023 Christopher Dilks
0005
0006 # make a table of cross section values, for each path to xsec files;
0007 # given a path, the average cross section of the files will be calculated
0008
0009 import os, re, pprint
0010 from numpy import mean
0011
0012 # main directory tree root #############
0013 xsecDir = "datarec/xsec"
0014 ########################################
0015
0016 # tree to hold cross sections and errors
0017 # - structure: { path : { 'val':[list], 'err':[list] } }
0018 # - one list element = one file's numbers, for each file in the path
0019 # - later we will take list averages for each path, and build the table
0020 xsecDict = {}
0021
0022 # traverse directories
0023 for (path,dirs,files) in os.walk(xsecDir):
0024
0025     # truncated path, to be written to table, used as a key in `xsecDict`
0026     pathTrun = re.sub(xsecDir,'',path)
0027     pathTrun = re.sub('^/','',pathTrun)
0028
0029     # loop over xsec files, assume the last line is the most accurate
0030     # calculation of the cross section, and add the value and error
0031     # to `xsecDict` in the appropriate lists
0032     for xsecFile in files:
0033         if re.search(r'\.xsec$',xsecFile):
0034
0035             if not pathTrun in xsecDict:
0036                 print(f'\nPATH = {pathTrun}')
0037                 xsecDict[pathTrun] = {'val':[],'err':[]}
0038
0039             with open(path+'/'+xsecFile) as xf:
0040                 print(f'  FILE = {xsecFile}')
0041                 firstLine = ""
0042                 for line in xf:
0043                   if firstLine=="": firstLine=line
0044                   pass # sets `line` to be the last line
0045                 xsecVal = float(line.split()[3])
0046                 xsecErr = float(line.split()[4])
0047                 xsecDelta = xsecVal - float(firstLine.split()[3]) # difference between first and last line's xsec
0048                 print(f'     xsec = {xsecVal} pb')
0049                 print(f'      err = {xsecErr} pb')
0050                 print(f' delta/xsec = {xsecDelta/xsecVal}')
0051                 xsecDict[pathTrun]['val'].append(xsecVal)
0052                 xsecDict[pathTrun]['err'].append(xsecErr)
0053
0054 print('\n')
0055 pprint.pprint(xsecDict)
0056 print('\n')
0057
0058 # compute average cross sections and tabulate
0059 tableFileName = xsecDir+"/xsec.dat"
0060 print(f'\nCROSS SECTION TABLE:  {tableFileName}\n'+'-'*60)
0061 tableFileTmpName = tableFileName+".tmp"
0062 tableFile = open(tableFileTmpName,'w+')
0063 tableFile.write('#label cross_section_[pb] relative_uncertainty\n')
0064 for path,nums in xsecDict.items():
0065   xsecVal = mean(nums['val'])
0066   xsecErr = mean(nums['err'])
0067   tableFile.write(f'{path} {xsecVal:.5g} {xsecErr/xsecVal:.3g}\n')
0068 tableFile.close()
0069 os.system(f'column -t {tableFileTmpName} | sort -n | tee {tableFileName}')
0070 os.remove(tableFileTmpName)