eic-opticks/ana/s2.py

0001 #/usr/bin/env python
0002 """
0003
0004
0005
0006 In [25]: np.c_[s2.ri, np.arange(len(s2.ri))+1 ]
0007 Out[25]:
0008 array([[ 1.55 ,  1.478,  1.   ],
0009        [ 1.795,  1.48 ,  2.   ],
0010        [ 2.105,  1.484,  3.   ],
0011        [ 2.271,  1.486,  4.   ],
0012        [ 2.551,  1.492,  5.   ],
0013        [ 2.845,  1.496,  6.   ],
0014        [ 3.064,  1.499,  7.   ],
0015        [ 4.133,  1.526,  8.   ],
0016        [ 6.2  ,  1.619,  9.   ],
0017        [ 6.526,  1.618, 10.   ],
0018        [ 6.889,  1.527, 11.   ],
0019        [ 7.294,  1.554, 12.   ],
0020        [ 7.75 ,  1.793, 13.   ],
0021        [ 8.267,  1.783, 14.   ],
0022        [ 8.857,  1.664, 15.   ],
0023
0024        [ 9.538,  1.554, 16.   ],
0025        [10.33 ,  1.454, 17.   ],
0026
0027        [15.5  ,  1.454, 18.   ]])
0028
0029
0030
0031 """
0032
0033 import os, numpy as np, logging
0034 log = logging.getLogger(__name__)
0035
0036 import matplotlib.pyplot as plt
0037 from opticks.ana.main import opticks_main
0038 from opticks.ana.key import keydir
0039
0040 def test_digitize():
0041     edom = np.linspace(-1., 12., 131 )
0042     ebin = np.linspace( 1., 10., 10 )
0043     idom = np.digitize(edom, ebin, right=False)  # right==False (the default) indicates the interval does not include the right edge
0044     # below and above yields 0 and len(ebin) = ebin.shape[0]
0045     #print(np.c_[edom,idom])
0046     for e, i in zip(edom,idom):
0047         mkr = "---------------------------" if e in ebin else ""
0048         print(" %7.4f : %d : %s " % (e, i, mkr))
0049     pass
0050
0051
0052
0053 class S2(object):
0054     def __init__(self):
0055         kd = keydir(os.environ["OPTICKS_KEY"])
0056         ri = np.load(os.path.join(kd, "GScintillatorLib/LS_ori/RINDEX.npy"))
0057         ri[:,0] *= 1e6
0058
0059         ri_ = lambda e:np.interp(e, ri[:,0], ri[:,1] )
0060         self.ri = ri
0061         self.ri_ = ri_
0062         self.s2x = np.repeat( np.nan, len(self.ri) )
0063
0064
0065     def bin(self, idx, lhs=1., rhs=1.):
0066         """
0067         Example, bins array of shape (4,) with 3 bins::
0068
0069          |       +-------------+--------------+-------------+         |
0070               0        1             2               3            4
0071
0072
0073         Hmm under/over bins are dangerous for "inflation"
0074         """
0075         dom = self.ri[:,0]
0076         val = self.ri[:,1]
0077
0078         d = None
0079         v = None
0080         if idx == 0:
0081             d = [dom[0] - lhs, dom[0] ]
0082             v = [val[0], val[0]]
0083         elif idx > 0 and idx < len(dom):
0084             d = [dom[idx-1], dom[idx] ]
0085             v = [val[idx-1], val[idx]]
0086         elif idx == len(dom):
0087             d = [dom[idx-1], dom[idx-1] + rhs]
0088             v = [val[idx-1], val[idx-1]]
0089         pass
0090         return np.array(d), np.array(v)
0091
0092
0093     def bin_dump(self, BetaInverse=1.5):
0094         """
0095         """
0096         for idx in range(len(self.ri)+1):
0097             en, ri = self.bin(idx)
0098             s2i, branch, en_cross, meta = self.bin_integral(idx, BetaInverse)
0099             print( " idx %3d  en %20s   ri %20s  s2i %10.4f  " % (idx, str(en), str(ri), s2i ))
0100         pass
0101
0102     def s2_linear_integral(self, en_, ri_, BetaInverse, fix_en_cross=None):
0103         """
0104         :param en_:
0105         :param ri_:
0106         :param BetaInverse:
0107         :param fix_en_cross:
0108
0109         Consider integrating in a region close to
0110         where s2 goes from +ve to -ve with en_0 with s2_0 > 0
0111
0112         As the en_1 is increased beyond the crossing where s2 is zero
0113         the calculation of en_cross will vary a little
0114         depending on en_1 and s2_1 despite the expectation
0115         that extending the en_1 shouuld not be adding to the integral.
0116
0117         This is problematic for forming the CDF because it may cause
0118         the result of the cumulative integral to become slightly non-monotonic.
0119
0120         The root of the problem is repeated calculation of en_cross
0121         as the en_1 in increased despite there being no additional info
0122         only chance to get different en_cross from linear extrapolation numerical
0123         imprecision.
0124
0125         For each crossing bin need to calulate the en_cross once and reuse that.
0126         Then can use fix_en_cross optional argument to ensure that the
0127         same en_cross is used for all integrals relevant to a bin.
0128
0129
0130              0.
0131              |  .
0132              |     .
0133              +-------X----------
0134                         .
0135                            1
0136
0137
0138          Similar triangles to find the en_cross where ri == BetaInverse.
0139          This has advantage over s2 crossing zero in that it is more linear,
0140          so the crossing should be a bit more precise.
0141
0142
0143              en_cross - en_0            en_1 - en_0
0144              -------------------  =  --------------------
0145              BetaInverse - ri_0         ri_1 - ri_0
0146
0147              en_cross_0 =  en_0 + (BetaInverse - ri_0)*(en_1 - en_0)/(ri_1 - ri_0)
0148
0149
0150              en_1 - en_cross            en_1 - en_0
0151              -------------------  =  --------------------
0152              ri_1 - BetaInverse         ri_1 - ri_0
0153
0154              en_cross_1 =  en_1 - (ri_1 - BetaInverse)*(en_1 - en_0)/(ri_1 - ri_0)
0155              en_cross_1 =  en_1 + (BetaInverse - ri_1)*(en_1 - en_0)/(ri_1 - ri_0)
0156
0157
0158              ## when there is no crossing in the range en_cross will not be between en_0 and en_1
0159              ## note potential infinities for flat ri bins,  ri_1 == ri_0
0160
0161         """
0162
0163         en = np.asarray(en_)
0164         ri = np.asarray(ri_)
0165
0166         ct = BetaInverse/ri
0167         s2 = (1.-ct)*(1.+ct)
0168
0169         en_0, en_1 = en
0170         ri_0, ri_1 = ri
0171         s2_0, s2_1 = s2
0172
0173         branch = 0
0174         if s2_0 <= 0. and s2_1 <= 0.:      # s2 -ve OR 0.      : no CK in bin
0175             branch = 1
0176         elif s2_0 < 0. and s2_1 > 0.:      # s2 goes +ve       : CK on RHS of bin
0177             branch = 2
0178         elif s2_0 >= 0. and s2_1 >= 0.:    # s2 +ve or zero    : CK across full bin
0179             branch = 3
0180         elif s2_0 > 0. and s2_1 <= 0.:     # s2 goes -ve or 0. : CK on LHS of bin
0181             branch = 4
0182         else:
0183             assert 0
0184         pass
0185
0186         if branch == 2 or branch == 4:
0187             en_cross_A = (s2_1*en_0 - s2_0*en_1)/(s2_1 - s2_0)    ## s2 zeros should occur at ~same en_cross, but non-linear so less precision
0188             en_cross_B = en_0 + (BetaInverse - ri_0)*(en_1 - en_0)/(ri_1 - ri_0)
0189             en_cross_C = en_1 + (BetaInverse - ri_1)*(en_1 - en_0)/(ri_1 - ri_0)
0190             en_cross = en_cross_B if fix_en_cross is None else fix_en_cross
0191         else:
0192             en_cross_A = np.nan
0193             en_cross_B = np.nan
0194             en_cross_C = np.nan
0195             en_cross = np.nan
0196         pass
0197
0198         if branch == 1.:
0199             area = 0.
0200         elif branch == 2:                            # s2 goes +ve
0201             area = (en_1 - en_cross)*s2_1*0.5
0202         elif branch == 3:                            # s2 +ve or zero across full range
0203             area = (en_1 - en_0)*(s2_0 + s2_1)*0.5
0204         elif branch == 4:                            # s2 goes -ve or to zero
0205             area = (en_cross - en_0)*s2_0*0.5
0206         else:
0207             assert 0
0208         pass
0209
0210         Rfact = 369.81 / 10.
0211         s2i = Rfact * area
0212
0213         meta = np.zeros(16)
0214
0215         meta[0] = en_0
0216         meta[1] = en_1
0217         meta[2] = s2_0
0218         meta[3] = s2_1
0219
0220         meta[4] = ri_0
0221         meta[5] = ri_1
0222         meta[6] = float(branch)
0223         meta[7] = s2i
0224
0225         meta[8] = en_cross_A
0226         meta[9] = en_cross_B
0227         meta[10] = en_cross_C
0228         meta[11] = en_cross
0229
0230         meta[12] = 0.
0231         meta[13] = 0.
0232         meta[14] = 0.
0233         meta[15] = 0.
0234
0235         assert s2i >= 0.
0236         return s2i, en_cross, branch, meta
0237
0238     def bin_integral(self, idx, BetaInverse, fix_en_cross=None):
0239         """
0240         :param idx: 0-based bin index
0241         :param BetaInverse: scalar
0242         :return s2integral: scalar sin^2 ck_angle integral in energy range
0243         """
0244         en, ri = self.bin(idx)
0245         return self.s2_linear_integral(en, ri, BetaInverse, fix_en_cross=fix_en_cross )
0246
0247     def full_integral(self, BetaInverse):
0248         """
0249         :param idx: 0-based bin index
0250         :param BetaInverse: scalar
0251         :return s2integral: scalar sin^2 ck_angle integral in energy range
0252         """
0253         s2integral = 0.
0254         for idx in range(1, len(self.ri)):
0255             s2i, en_cross, branch, meta  = self.bin_integral(idx, BetaInverse)
0256             s2integral += s2i
0257         pass
0258         return s2integral
0259
0260     def full_integral_scan(self):
0261         bis = np.linspace(1, 2, 11)
0262         for BetaInverse in bis:
0263             s2integral = self.full_integral(BetaInverse)
0264             print(" bi %7.4f s2integral  %7.4f  " % (BetaInverse, s2integral))
0265         pass
0266
0267     def cut_integral(self, ecut, BetaInverse=1.5, dump=False, line=0 ):
0268         """
0269         """
0270         ## find index of the bin containing the ecut value
0271         idx = int(np.digitize( ecut, self.ri[:,0], right=False ))
0272
0273         s2integral = np.zeros(len(self.ri)+1)
0274
0275         ## first add up integrals from full bins to the left of idx bin
0276         for i in range(1,idx):
0277             s2integral[i], _en_cross, _branch, _meta = self.bin_integral(i, BetaInverse)
0278         pass
0279
0280         ## find details of the idx bin
0281         en, ri = self.bin(idx)
0282
0283         ## use full bin integral to obtain en_cross
0284         full, full_en_cross, full_branch,  _ = self.bin_integral(idx, BetaInverse, fix_en_cross=None)
0285
0286         ## add s2integral over the partial bin
0287         en_0, en_1 = en[0], ecut
0288         ri_0, ri_1 = ri[0], self.ri_(ecut)
0289
0290         # fix_en_cross only gets used when appropriate
0291         s2integral[-1], en_cross, branch, meta = self.s2_linear_integral( [en_0, en_1], [ri_0, ri_1], BetaInverse, fix_en_cross=full_en_cross )
0292
0293         tot = s2integral.sum()
0294
0295         meta[-4] = tot
0296         meta[-3] = float(idx)
0297         meta[-2] = ecut
0298         meta[-1] = float(line)
0299
0300         if dump:
0301             print("cut_integral  %10.4f idx %3d  en %15s ri %15s  tot %10.4f  " % (ecut, idx, str(en), str(ri), tot ))
0302         pass
0303         return s2integral, meta
0304
0305     def cs2_integral(self, efin, BetaInverse=1.5):
0306         """
0307         Fixed slight non-monotonic by forcing use of common en_cross obtained from the full bin integral
0308         in the partial bin integrals.  Further improvement using ri-BetaInverse cross rather than s2 cross,
0309         as thats more linear.
0310         """
0311         cs2i = np.zeros( [len(efin), len(self.ri)+1])
0312         meta = np.zeros( [len(efin), 16])
0313
0314         for line,ecut in enumerate(efin):
0315             cs2i[line], meta[line] = self.cut_integral(ecut, BetaInverse, line=line)
0316         pass
0317         s_cs2i = cs2i.sum(axis=1)  # sum over full bins and the partial bin for each cut
0318
0319         monotonic = np.all(np.diff(s_cs2i) >= 0.)
0320         assert monotonic
0321         return cs2i, meta
0322
0323
0324 def test_full_vs_cut(s2, dump=False):
0325     bis = np.linspace(1., 2., 101)
0326     res = np.zeros( (len(bis), 4) )
0327     for i, BetaInverse in enumerate(bis):
0328         full = s2.full_integral(BetaInverse)
0329         cut , _cut_meta = s2.cut_integral(s2.ri[-1,0], BetaInverse, dump=dump)
0330         cut_tot = cut.sum()
0331         res[i] = [BetaInverse, full, cut_tot, np.abs(full-cut_tot)*1e10 ]
0332         assert np.abs(full - cut_tot) < 1e-10
0333     pass
0334     return res
0335
0336
0337
0338 if __name__ == '__main__':
0339     logging.basicConfig(level=logging.INFO)
0340     np.set_printoptions(edgeitems=50, precision=4)
0341
0342     s2 = S2()
0343
0344 if 0:
0345     s2.bin_dump()
0346     s2.full_integral_scan()
0347     res = test_full_vs_cut(s2)
0348
0349 if 1:
0350     nx = 100
0351     efin = np.linspace( s2.ri[0,0], s2.ri[-1,0], nx )
0352     cs2i, meta = s2.cs2_integral(efin, BetaInverse=1.5)
0353
0354     # np.c_[cs2i, np.arange(len(efin)), cs2i.sum(axis=1), efin  ]
0355
0356