[807] | 1 | # $Id: extractShelve.py,v 1.3 2003/07/04 13:26:31 dressel Exp $ |
---|
| 2 | # ------------------------------------------------------------------- |
---|
| 3 | # GEANT4 tag $Name: $ |
---|
| 4 | # ------------------------------------------------------------------- |
---|
| 5 | # |
---|
| 6 | import string |
---|
| 7 | import math |
---|
| 8 | import os |
---|
| 9 | import shelve |
---|
| 10 | import string |
---|
| 11 | import detector |
---|
| 12 | |
---|
| 13 | |
---|
| 14 | |
---|
| 15 | # some functions for calculating the FOM |
---|
| 16 | |
---|
| 17 | def getR2(measure): |
---|
| 18 | "Error squared of a measure." |
---|
| 19 | v = measure.variance |
---|
| 20 | m2 = math.pow(measure.mean, 2) |
---|
| 21 | n = measure.entries |
---|
| 22 | |
---|
| 23 | r2 = 1e+10 |
---|
| 24 | if n>0 and m2 > 0: |
---|
| 25 | r2 = v/m2/n |
---|
| 26 | |
---|
| 27 | return r2 |
---|
| 28 | |
---|
| 29 | def getFOMbyName(shelveName, tallyName, bin): |
---|
| 30 | "FOM by shelve name, tally name and bin of energy region" |
---|
| 31 | she = shelve.open (shelveName, "r") |
---|
| 32 | measure = she[tallyName].measures[bin] |
---|
| 33 | R2= getR2(measure) |
---|
| 34 | T = she["runTime"] |
---|
| 35 | return 1./(T*R2) |
---|
| 36 | |
---|
| 37 | def getFOM(aShelve, tallyName, bin): |
---|
| 38 | "FOM given a shelve, tally name and energy bin." |
---|
| 39 | measure = aShelve[tallyName].measures[bin] |
---|
| 40 | R2= getR2(measure) |
---|
| 41 | T = aShelve["runTime"] |
---|
| 42 | return 1./(T*R2) |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | |
---|
| 53 | # functions to retrive information from the shelve |
---|
| 54 | |
---|
| 55 | def infoShelve(file): |
---|
| 56 | "Print information stored in a shelve." |
---|
| 57 | she = shelve.open(file, "r") |
---|
| 58 | print file |
---|
| 59 | for k in she.keys(): |
---|
| 60 | value = she[k] |
---|
| 61 | if value.__class__.__module__ == "__builtin__" : |
---|
| 62 | print " ", k, ":", value |
---|
| 63 | |
---|
| 64 | def lsShelve(path): |
---|
| 65 | "List shelve information found in and under the given path. " |
---|
| 66 | if os.path.isfile(path): |
---|
| 67 | if string.find(path,".shelve") > -1: |
---|
| 68 | infoShelve(path) |
---|
| 69 | else: |
---|
| 70 | if os.path.isdir(path): |
---|
| 71 | files = os.listdir(path) |
---|
| 72 | for file in files: |
---|
| 73 | f=path + "/" + file |
---|
| 74 | lsShelve(f) |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | |
---|
| 78 | def getFluxInRegion(she, dist="00", detType="ring", region=1): |
---|
| 79 | "Get the flux in a detector." |
---|
| 80 | tally = she["detector_" + dist + "Tally"] |
---|
| 81 | m = tally.measures[region] |
---|
| 82 | rawFlux = m.sum |
---|
| 83 | mean = m.mean |
---|
| 84 | var = m.getVariance() |
---|
| 85 | n = m.entries |
---|
| 86 | errRel = -1. |
---|
| 87 | if n>0 and var >0 and mean > 0: |
---|
| 88 | errRel = math.sqrt(var/n)/mean |
---|
| 89 | nPeakNeutrons = she["generatorTally"].measures[1].sum |
---|
| 90 | scale = detector.detScale(nPeakNeutrons, she["energy"], 0) |
---|
| 91 | return rawFlux * scale / detector.DetectorVolume[detType][dist], errRel |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | def getFluxes(she, detType = "ring"): |
---|
| 96 | "Get all fluxes and FOM in a shelve." |
---|
| 97 | fluxName = "flux_" + she["energy"]+"_"+\ |
---|
| 98 | she["shieldWidth"] |
---|
| 99 | keys = she.keys() |
---|
| 100 | fluxes = {} |
---|
| 101 | for k in keys: |
---|
| 102 | if string.find(k,"detector_") > -1: |
---|
| 103 | sp = string.split(k,"detector_") |
---|
| 104 | sp = string.split(sp[1],"Tally") |
---|
| 105 | dist = sp[0] |
---|
| 106 | fluxNameD = fluxName + "_" + dist |
---|
| 107 | for i in range(1,3): |
---|
| 108 | fom = getFOM(she, k, i) |
---|
| 109 | f = fluxNameD + "_" + "%(i)d" % vars() |
---|
| 110 | flux, errRel = getFluxInRegion(she, dist, detType, i) |
---|
| 111 | errRel*=100 |
---|
| 112 | sFlux = '%1.2E' % (flux) |
---|
| 113 | sErrRel = '%1.0f' % (errRel) |
---|
| 114 | if errRel < 10.0: |
---|
| 115 | sErrRel = '%1.1f' % (errRel) |
---|
| 116 | print f, sFlux, ", errRel:", sErrRel, " FOM:", fom |
---|
| 117 | # fluxes[f] = sFlux |
---|
| 118 | return fluxes |
---|