| 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
|
|---|
