u4/tests/U4Scintillation_Debug.py

0001 #!/usr/bin/env python
0002 """
0003 U4Scintillation_Debug.py
0004 ==========================
0005
0006 DsG4Scintillation.cc::
0007
0008      358         //////////////////////////////////// Birks' law ////////////////////////
0009      359         // J.B.Birks. The theory and practice of Scintillation Counting.
0010      360         // Pergamon Press, 1964.
0011      361         // For particles with energy much smaller than minimum ionization
0012      362         // energy, the scintillation response is non-linear because of quenching
0013      363         // effect. The light output is reduced by a parametric factor:
0014      364         // 1/(1 + birk1*delta + birk2* delta^2).
0015      365         // Delta is the energy loss per unit mass thickness. birk1 and birk2
0016      366         // were measured for several organic scintillators.
0017      367         // Here we use birk1 = 0.0125*g/cm2/MeV and ignore birk2.
0018      368         // R.L.Craun and D.L.Smith. Nucl. Inst. and Meth., 80:239-244, 1970.
0019      369         // Liang Zhan  01/27/2006
0020      370         // /////////////////////////////////////////////////////////////////////
0021      371
0022
0023 Light output reduced by parametric factor::
0024
0025     1/(1 + birk1*delta + birk2* delta^2)
0026
0027
0028 * birk1 : 0.0125*g/cm2/MeV   : density*length/energy
0029 * delta :                    : energy/length/density
0030 * birk1*delta                : unitless
0031
0032
0033 G4double dE = TotalEnergyDeposit;
0034 G4double dx = aStep.GetStepLength();
0035 G4double dE_dx = dE/dx;             // energy over length
0036 G4double delta = dE_dx/aMaterial->GetDensity();//get scintillator density
0037
0038 In [10]: d["dE_dx"]
0039 Out[10]: 9.293368125518313
0040
0041 In [11]: d["dE_dx"]/d["Density"]
0042 Out[11]: 1.7333662276764504e-18
0043
0044 In [12]: d["delta"]
0045 Out[12]: 1.7333662276764504e-18
0046
0047 In [13]: d["TotalEnergyDeposit"]/d["dx"]
0048 Out[13]: 7.887351897652502e-08
0049
0050
0051
0052
0053
0054
0055
0056
0057 """
0058
0059 import os, re, numpy as np
0060
0061 class _U4Scintillation_Debug(object):
0062     PKGDIR = os.path.dirname(os.path.realpath(os.path.dirname(__file__)))
0063     PTN = re.compile("^    double (\S*) ;$")
0064     @classmethod
0065     def Parse(cls):
0066         path = os.path.join(cls.PKGDIR, "U4Scintillation_Debug.hh")
0067         lines = open(path).read().splitlines()
0068         fields = []
0069         for line in lines:
0070             match = cls.PTN.match(line)
0071             if match:
0072                 g = match.groups()[0]
0073                 fields.append(g)
0074             pass
0075         pass
0076         return fields
0077
0078
0079 class U4Scintillation_Debug(object):
0080     BASE = "/tmp/u4debug"
0081     NAME = "U4Scintillation_Debug.npy"
0082     FIELDS = _U4Scintillation_Debug.Parse()
0083
0084     def __init__(self, rel, symbol):
0085         path = os.path.join(self.BASE, rel, self.NAME)
0086         self.a = np.load(path)
0087         self.rel = rel
0088         self.symbol = symbol
0089         self.path = path
0090
0091     def __repr__(self):
0092         return "%s %s %s %s" % (self.symbol, self.rel, str(self.a.shape), self.path)
0093
0094
0095     def QuenchedTotalEnergyDeposit(self, idx):
0096         d = self[idx]
0097         return d["TotalEnergyDeposit"]/(1.+d["birk1"]*d["delta"])
0098
0099     def QuenchedTotalEnergyDeposit_untypo(self, idx):
0100         d = self[idx]
0101         return d["TotalEnergyDeposit"]/(1.+d["birk1"]*1e-6*d["delta"])
0102
0103     def MeanNumberOfPhotons(self, idx):
0104         d = self[idx]
0105         return d["ScintillationYield"]*d["QuenchedTotalEnergyDeposit"]
0106
0107     def MeanNumberOfPhotons_untypo(self, idx):
0108         d = self[idx]
0109         return d["ScintillationYield"]*self.QuenchedTotalEnergyDeposit_untypo(idx)
0110
0111     def __getitem__(self, idx):
0112         vals = self.a[idx].ravel()
0113         keys = self.FIELDS
0114         return dict(zip(keys, vals))
0115
0116     def desc_MeanNumberOfPhotons(self, idx):
0117         d = self[idx]
0118         al = "d[\"MeanNumberOfPhotons\"]"
0119         av = d["MeanNumberOfPhotons"]
0120
0121         bl = "%s.MeanNumberOfPhotons(%3d)" % (self.symbol, idx)
0122         bv = self.MeanNumberOfPhotons(idx)
0123
0124         cl = "%s.MeanNumberOfPhotons_untypo(%3d)" % (self.symbol, idx)
0125         cv = self.MeanNumberOfPhotons_untypo(idx)
0126
0127         line = "%4d : %s:%10.4f %s:%10.4f %s:%10.4f " % (idx, al,av,bl,bv,cl,cv)
0128         return line
0129
0130     def __str__(self):
0131         lines = []
0132         for idx in list(range(len(self.a))):
0133             line = self.desc_MeanNumberOfPhotons(idx)
0134             lines.append(line)
0135         pass
0136         return "\n".join(lines)
0137
0138
0139
0140 if __name__ == '__main__':
0141     s30 = U4Scintillation_Debug("ntds3/000", "s30" )
0142     s31 = U4Scintillation_Debug("ntds3/001", "s31" )
0143
0144     print(s30)
0145     print(s31)
0146     print(U4Scintillation_Debug.FIELDS)
0147
0148     print(str(s30))
0149
0150