Everything that is constant from the point of view of the photon cohort needs to be collected into the “G4DAEStep” object. Although some things could potentially be looked up from material props on GPU no point doing that when are constant for all photons, just do it once and present as parameter to the kernel launch.
Three files to check for each process (Scintillation and Cerenkov):
No longer need to manually copy any files, the GPU generated photons are now “sidesaved” under types opcerenkov and opscintillation and the Geant4 equivalents are “onlycopied” from DsChromaEventAction::EndOfEventAction via metadata control such as ctrl:noreturn:1 and ctrl:onlycopy:1 applied to processing of all species.
In [1]: cf('wavelength', tag=1, typs='gopcerenkov opcerenkovgen')
## improved by shifting low wavelength to match G4 change 60:801:20 to 80:801:20
In [3]: cf('3xyzw', tag=1, typs='gopcerenkov opcerenkovgen', legend=False)
In [1]: cf('3xyzw', tag=1, typs='gopscintillation opscintillationgen', legend=False)
In [2]: cf('wavelength', tag=1, typs='gopscintillation opscintillationgen', log=True)
Check counts from ipython:
In [2]: genconsistency(1)
INFO:env.g4dae:g4c : GOPCERENKOV (612841, 4, 4)
INFO:env.g4dae:chc : OPCERENKOV (612841, 4, 4)
INFO:env.g4dae:stc : CERENKOV (7836, 6, 4) ==> N 612841
INFO:env.g4dae:g4s : GOPSCINTILLATION (2817543, 4, 4)
INFO:env.g4dae:chs : OPSCINTILLATION (2817543, 4, 4)
INFO:env.g4dae:sts : SCINTILLATION (13898, 6, 4) ==> N 2817543
In [1]: chc, g4c, tst = NPY.mget(1,"opcerenkov","gopcerenkov","test")
Initially had very different wavelength, with chroma flat due to use of uniform wavelength draw, fixed by moving to uniform in reciprocal of wavelenth:
In [1]: cf('wavelength', tag=1, typs="gopcerenkov opcerenkov test")
Create some new “test” arrays with modified standard wavelengths:
g4daechroma.sh --wavelengths 80:801:20
npysend.sh -t1 -icerenkov -otest
G4 distrib has step at 200nm, an artifact of RINDEX edge of quoted range:
In [105]: chroma_refractive_index(cg)
[ 0] LiquidScintillator (18, 2) wl 79.99 : 799.90 1.454 : 1.793
[ 3] GdDopedLS (18, 2) wl 79.99 : 799.90 1.454 : 1.793
[ 5] Acrylic (18, 2) wl 79.99 : 799.90 1.462 : 1.793
[10] MineralOil (18, 2) wl 79.99 : 799.90 1.434 : 1.759
[24] IwsWater (34, 2) wl 199.97 : 799.90 1.333 : 1.390
[27] OwsWater (34, 2) wl 199.97 : 799.90 1.333 : 1.390
[28] DeadWater (34, 2) wl 199.97 : 799.90 1.333 : 1.390
In [6]: plot_refractive_index()
In [7]: plot_refractive_index(standard=True)
In [130]: for i in np.unique(im):print "%2d : %5d : %s " % ( i, bc[i], cg.unique_materials[i].name[17:-9] )
0 : 1133 : LiquidScintillator
3 : 4220 : GdDopedLS
5 : 88 : Acrylic
10 : 1046 : MineralOil
24 : 791 : IwsWater
27 : 530 : OwsWater
28 : 28 : DeadWater
In [9]: cf('time', tag=1, typs="gopcerenkov opcerenkov test")
In [9]: cf('3xyz', g4c, chc)
chroma/chroma/cuda/cerenkov.h:
202 __device__ void
203 generate_cerenkov_photon(Photon& p, CerenkovStep& cs, curandState &rng)
204 {
205 float cosTheta ;
206 float sin2Theta ;
207 float wavelength ;
208 float sampledRI ;
209
210 //
211 // sampling to get wavelength and cone angle
212 //
213 // pick random wavelength inside the range,
214 // lookup refractive index
215 // calculate cosTheta and sinTheta for the refractive index
216 //
217 do {
218
219 wavelength = sample_value(cs.material, curand_uniform(&rng));
220
221 sampledRI = interp_property(cs.material, wavelength, cs.material->refractive_index);
222
223 cosTheta = cs.BetaInverse / sampledRI;
224
225 sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
226
227 } while ( curand_uniform(&rng)*cs.maxSin2 > sin2Theta);
228
229
230 p.wavelength = wavelength ;
231
296 G4double Pmin = Rindex->GetMinPhotonEnergy();
297 G4double Pmax = Rindex->GetMaxPhotonEnergy();
298 G4double dp = Pmax - Pmin;
405 for (G4int i = 0; i < NumPhotons; i++) {
406 // Determine photon energy
407 G4double rand=0;
408 G4double sampledEnergy=0, sampledRI=0;
409 G4double cosTheta=0, sin2Theta=0;
410
411 // sample an energy
412 do {
413 rand = G4UniformRand();
414 sampledEnergy = Pmin + rand * dp;
415 sampledRI = Rindex->GetProperty(sampledEnergy);
416 cosTheta = BetaInverse / sampledRI;
417
418 sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
419 rand = G4UniformRand();
420
421 } while (rand*maxSin2 > sin2Theta);
422
In [48]: cls.refractive_index
Out[48]:
array([[ 79.99 , 1.454],
[ 120.023, 1.454],
[ 129.99 , 1.554],
[ 139.984, 1.664],
[ 149.975, 1.783],
[ 159.98 , 1.793],
[ 169.981, 1.554],
[ 179.974, 1.527],
[ 189.985, 1.618],
[ 199.975, 1.618],
[ 300. , 1.526],
[ 404.7 , 1.499],
[ 435.8 , 1.495],
[ 486.001, 1.492],
[ 546.001, 1.486],
[ 589.002, 1.484],
[ 690.701, 1.48 ],
[ 799.898, 1.478]], dtype=float32)
In [49]: cls.name
Out[49]: '__dd__Materials__LiquidScintillator0xc2308d0'
In [50]: ri = cls.refractive_index
In [51]: plt.scatter(ri[:,0],ri[:,1])
Out[51]: <matplotlib.collections.PathCollection at 0x125b76a90>
In [52]: plt.show()
In [53]: _stc = stc(1)
In [56]: BetaInverse = _stc[:,4,0]
Out[56]: array([ 1., 1., 1., ..., 1., 1., 1.], dtype=float32)
In [57]: BetaInverse.min()
Out[57]: 1.0000062
In [58]: BetaInverse.max()
Out[58]: 1.4531251
In [64]: plt.hist(BetaInverse, bins=100,log=True) # mainly 1.000 with small tail out to 1.45
In [107]: _stc[:,0].view(np.int32)
Out[107]:
array([[ -1, 1, 24, 80],
[ -2, 1, 24, 108],
[ -3, 1, 24, 77],
...,
[-7834, 1, 28, 91],
[-7835, 1, 28, 83],
[-7836, 1, 28, 48]], dtype=int32)
In [108]: _stc[:,0,2].view(np.int32)
Out[108]: array([24, 24, 24, ..., 28, 28, 28], dtype=int32)
In [110]: im
Out[110]: array([24, 24, 24, ..., 28, 28, 28], dtype=int32)
In [111]: np.unique(im)
Out[111]: array([ 0, 3, 5, 10, 24, 27, 28], dtype=int32)
In [129]: bc = np.bincount(im)
In [130]: for i in np.unique(im):print "%2d : %5d : %s " % ( i, bc[i], cg.unique_materials[i].name[17:-9] )
0 : 1133 : LiquidScintillator
3 : 4220 : GdDopedLS
5 : 88 : Acrylic
10 : 1046 : MineralOil
24 : 791 : IwsWater
27 : 530 : OwsWater
28 : 28 : DeadWater
G4DAEChroma/G4DAECerenkovStep.hh:
13 enum {
14
15 _Id, // 0
16 _ParentID,
17 _Material,
18 _NumPhotons,
19
20 _x0_x, // 1
21 _x0_y,
22 _x0_z,
23 _t0,
24
25 _DeltaPosition_x, // 2
26 _DeltaPosition_y,
27 _DeltaPosition_z,
28 _step_length,
29
30 _code, // 3
31 _charge,
32 _weight,
33 _MeanVelocity,
34
35 _BetaInverse, // 4
36 _Pmin,
37 _Pmax,
38 _maxCos,
39
40 _maxSin2, // 5
41 _MeanNumberOfPhotons1,
42 _MeanNumberOfPhotons2,
43 _BialkaliMaterialIndex,
In [73]: maxSin2 = _stc[:,5,0]
In [76]: plt.hist(maxSin2, bins=100, log=True) ## mostly flat with few spikes at high end
In [82]: maxSin2.min()
Out[82]: 0.00065323891
In [83]: maxSin2.max()
Out[83]: 0.53214556
BialkaliMaterialIndex:
n [69]: _stc[:,5,3].view(np.int32).min()
Out[69]: 7
In [70]: _stc[:,5,3].view(np.int32).max()
Out[70]: 7
In [71]: cg.unique_materials[7]
Out[71]: <chroma.geometry.Material at 0x125a9a950>
In [72]: cg.unique_materials[7].name
Out[72]: '__dd__Materials__Bialkali0xc2f2428'
TODO: settle on standard wavelenth range to match G4 better:
In [4]: cf('3xyzw', tag=1, typs='gopcerenkov opcerenkov', legend=False, log=True)
In [8]: g4s, chs = NPY.mget(1, "gopscintillation opscintillation")
In [6]: cf('wavelength', typs="gopscintillation opscintillation",tag=1, log=True, range=(100,900) ) ## hmm clear chroma cut at 600nm ???
Scintillation wavelength, chroma distrib is faithfully representing a “histogram” stepping shape with “bins” of about 25nm. Looks like a problem of mismatched histogram ranges in the chroma sampling and the input histogram
chroma/chroma/geometry.py:
25 # all material/surface properties are interpolated at these
26 # wavelengths when they are sent to the gpu
27 standard_wavelengths = np.arange(60, 810, 20).astype(np.float32)
28
In [45]: standard_wavelengths = np.arange(60, 810, 20).astype(np.float32)
In [46]: standard_wavelengths
Out[46]:
array([ 60., 80., 100., 120., 140., 160., 180., 200., 220.,
240., 260., 280., 300., 320., 340., 360., 380., 400.,
420., 440., 460., 480., 500., 520., 540., 560., 580.,
600., 620., 640., 660., 680., 700., 720., 740., 760.,
780., 800.], dtype=float32)
In [47]: len(standard_wavelengths)
Out[47]: 38
The cdf becomes flat above 600nm:
In [41]: np.set_printoptions(precision=10)
In [42]: s_reemission_cdf
Out[42]:
array([[ 60. , 0. ],
[ 80. , 0. ],
[ 100. , 0. ],
[ 120. , 0. ],
[ 140. , 0. ],
[ 160. , 0. ],
[ 180. , 0. ],
[ 200. , 0.0000000205],
[ 220. , 0.0000161953],
[ 240. , 0.00003237 ],
[ 260. , 0.0000485448],
[ 280. , 0.0000647195],
[ 300. , 0.0000808943],
[ 320. , 0.000097069 ],
[ 340. , 0.0009699235],
[ 360. , 0.003526893 ],
[ 380. , 0.0114005441],
[ 400. , 0.1889369488],
[ 420. , 0.5017676353],
[ 440. , 0.7646831274],
[ 460. , 0.9037991762],
[ 480. , 0.9602615237],
[ 500. , 0.9843546748],
[ 520. , 0.9931191802],
[ 540. , 0.9967452288],
[ 560. , 0.9983744621],
[ 580. , 0.99932307 ],
[ 600. , 0.9999999404],
[ 620. , 1. ],
[ 640. , 1. ],
[ 660. , 1. ],
[ 680. , 1. ],
[ 700. , 1. ],
[ 720. , 1. ],
[ 740. , 1. ],
[ 760. , 1. ],
[ 780. , 1. ],
[ 800. , 1. ]], dtype=float32)
In [43]: s_reemission_cdf.shape
Out[43]: (38, 2)
171 __device__ void
172 generate_scintillation_photon(Photon& p, ScintillationStep& ss, curandState &rng)
173 {
174 float ScintillationTime = ss.ScintillationTime ;
175 if(ss.scnt == 2)
176 {
177 ScintillationTime = ss.slowTimeConstant ;
178 if(curand_uniform(&rng) < ss.slowerRatio)
179 {
180 ScintillationTime = ss.slowerTimeConstant ;
181 }
182 }
183
184 p.wavelength = sample_cdf(&rng, ss.material->n,
185 ss.material->wavelength0,
186 ss.material->step,
187 ss.material->reemission_cdf); // reemission_cdf poorly named, intensity_cdf better
188
random.h:
25 // Draw a random sample given a cumulative distribution function
26 // Assumptions: ncdf >= 2, cdf_y[0] is 0.0, and cdf_y[ncdf-1] is 1.0
27 __device__ float
28 sample_cdf(curandState *rng, int ncdf, float *cdf_x, float *cdf_y)
29 {
30 return interp(curand_uniform(rng),ncdf,cdf_y,cdf_x);
31 }
32
33 // Sample from a uniformly-sampled CDF
34 __device__ float
35 sample_cdf(curandState *rng, int ncdf, float x0, float delta, float *cdf_y)
36 {
37 float u = curand_uniform(rng);
38
39 int lower = 0;
40 int upper = ncdf - 1; // far left, right bin numbers
41
42 while(lower < upper-1) // still not settled on a bin
43 {
44 int half = (lower + upper) / 2;
45
46 if (u < cdf_y[half])
// cdf is normalized to 1 at rhs, so this is appropriate
47 upper = half;
48 else
49 lower = half;
50 }
51
52 float delta_cdf_y = cdf_y[upper] - cdf_y[lower];
53
54 return x0 + delta*lower + delta*(u-cdf_y[lower])/delta_cdf_y;
55 }
Bins (ie wavelengths) beyond where the CDF reaches 1 will never get sampled. The source distribution is SLOWCOMPONENT (same as FASTCOMPONENT):
plt_gdls()
Wavelength comes from sampling:
In [31]: np.allclose( cg.unique_materials[3].reemission_cdf, cg.unique_materials[0].reemission_cdf )
Out[31]: True
In [32]: reemission_cdf = cg.unique_materials[3].reemission_cdf
In [33]: s_reemission_cdf = standardize(reemission_cdf)
catplot(g4s, log=True, range=(100,900))
catplot(g4s, val='wavelength', cat='pdg', log=True, range=(100,900)) ## edges at 200, 800 nm
In [7]: g4s.wavelength.min()
Out[7]: G4ScintillationPhoton(199.97593688964844, dtype=float32)
In [8]: g4s.wavelength.max()
Out[8]: G4ScintillationPhoton(799.8797607421875, dtype=float32)
Where did those edges come from:
In [11]: gdls.extra.properties['SLOWCOMPONENT']
Out[11]:
array([[ 79.99 , 0. ],
[ 120.023, 0. ],
[ 199.975, 0. ], ## note 130nm jump to zero bin
[ 330. , 0.006],
[ 331. , 0.006],
[ 332. , 0.005],
[ 333. , 0.005],
...
[ 598.001, 0.002],
[ 599.001, 0.002],
[ 600.001, 0.002],
[ 799.898, 0. ]]) ## note 200nm jump to zero bin
catplot(g4s, val='time', cat='scnt', log=True, range=(0,50))
[blyth@ntugrid5 geant4.9.2.p01]$ find $PWD -name 'G4PhysicsVector.hh'
/home/blyth/local/env/dyb/external/build/LCG/geant4.9.2.p01/include/G4PhysicsVector.hh
/home/blyth/local/env/dyb/external/build/LCG/geant4.9.2.p01/source/global/management/include/G4PhysicsVector.hh
[blyth@ntugrid5 geant4.9.2.p01]$ find $PWD -name 'G4PhysicsOrderedFreeVector.hh'
/home/blyth/local/env/dyb/external/build/LCG/geant4.9.2.p01/include/G4PhysicsOrderedFreeVector.hh
/home/blyth/local/env/dyb/external/build/LCG/geant4.9.2.p01/source/global/management/include/G4PhysicsOrderedFreeVector.hh
[blyth@ntugrid5 geant4.9.2.p01]$
Grab the scintillation integrals using G4DAEPropList:
G4DAEArray::Allocate nitems 275 nfloat 550
G4DAEArray::Allocate nitems 275 nfloat 550
G4DAEArray::Allocate nitems 28 nfloat 56
G4DAEArray::SavePath [/home/blyth/local/env/prop/ls_fast.npy] itemcount 275 itemshape 2
G4DAEArray::SavePath [/home/blyth/local/env/prop/ls_slow.npy] itemcount 275 itemshape 2
G4DAEArray::SavePath [/home/blyth/local/env/prop/ls_reem.npy] itemcount 28 itemshape 2
G4DAEArray::Allocate nitems 275 nfloat 550
G4DAEArray::Allocate nitems 275 nfloat 550
G4DAEArray::Allocate nitems 28 nfloat 56
G4DAEArray::SavePath [/home/blyth/local/env/prop/gdls_fast.npy] itemcount 275 itemshape 2
G4DAEArray::SavePath [/home/blyth/local/env/prop/gdls_slow.npy] itemcount 275 itemshape 2
G4DAEArray::SavePath [/home/blyth/local/env/prop/gdls_reem.npy] itemcount 28 itemshape 2
physicsList->setCut() start.
delta:~ blyth$ export-prop-rget | sh
gdls_fast.npy 100% 2280 2.2KB/s 00:00
ls_slow.npy 100% 2280 2.2KB/s 00:00
gdls_slow.npy 100% 2280 2.2KB/s 00:00
gdls_reem.npy 100% 304 0.3KB/s 00:00
ls_reem.npy 100% 304 0.3KB/s 00:00
ls_fast.npy 100% 2280 2.2KB/s 00:00
Energy is xscaled to be in reciprocal wavelength (1/nm) and yscale is 1e9:
ls = pro_("ls_fast")
plt.plot(ls[:,0], ls[:,1])
plt.show()
plt.plot(1./ls[:,0], ls[:,1], "r+")
plt.show()
Establish connection between scintillation step and the transported scintillation integeral:
In [3]: g4s = ScintillationStep.get(1)
In [13]: np.unique(g4s[:,5,1]).item()*1e9
Out[13]: 410.0278374608024
Cheat by using purloined ScintillationIntegral in gdct- test_ScintillationIntegral succeeds to reproduce the scintillation wavelengths. But this is essentially using the same G4 code so no surprise:
int test_ScintillationIntegral()
{
G4DAEPropList* cdf = G4DAEPropList::Load("gdls_fast");
cdf->Print();
G4PhysicsOrderedFreeVector* ScintillationIntegral = G4DAEProp::CreatePOFV(cdf);
G4double MaxValue = ScintillationIntegral->GetMaxValue() ;
//size_t size = 1e6 ;
size_t size = 2817543 ; // match the count to current evt "1"
G4DAEArrayHolder* holder = new G4DAEArrayHolder( size, NULL, "2" );
for(size_t n=0 ; n<size ; n++ )
{
G4double CIIvalue = G4UniformRand()*MaxValue;
G4double sampledEnergy = ScintillationIntegral->GetEnergy(CIIvalue);
float* prop = holder->GetNextPointer();
prop[G4DAEProp::_binEdge] = float(CIIvalue) ;
prop[G4DAEProp::_binValue] = float(sampledEnergy) ;
}
G4DAEPropList dist(holder);
dist.Save("1"); // sampledEnergy
//
// cf('wavelength', typs="gopscintillation opscintillation prop",tag=1, log=True, range=(100,900) )
// succeeds to match G4 Scintillation photon distrib
//
return 0 ;
}
What about numpy level:
In [13]: cdf = pro_("gdls_fast")
In [16]: mx = cdf[:,1].max()
In [17]: mx
Out[17]: 410.02786
In [18]: u = np.random.rand( 2817543 )
In [19]: u.shape
Out[19]: (2817543,)
Need to invert x to have wavelength ordinate, but that makes CDF back to front:
In [32]: plt.plot(1/cdf[:,0],cdf[:,1])
Out[32]: [<matplotlib.lines.Line2D at 0x11645d5d0>]
So interpolate in 1/wavelength land and invert afterwards, this avoids the question of how to deal with infinite wavelength:
In [46]: wi = np.interp( u*cdf[:,1].max(), cdf[:,1], cdf[:,0] ) ## NB x-y flip
In [47]: w = 1/wi
In [51]: plt.hist(w, bins=100, log=True, range=(100,900)) ## looking good
So how does the purloined scintillation integral compare with what have been using.
In [3]: cg = cg_get()
In [7]: ls = cg.unique_materials[0]
In [13]: rcdf = ls.reemission_cdf
In [14]: plt.plot(rcdf[:,0], rcdf[:,1])
Out[14]: [<matplotlib.lines.Line2D at 0x116369050>]
In [15]: plt.show()
In [41]: plt.plot( 1/rcdf[:,0], rcdf[:,1], 'b+')
Out[41]: [<matplotlib.lines.Line2D at 0x10ccc4310>]
In [45]: plt.plot( cdf[:,0], 1 - cdf[:,1]/cdf[:,1].max(), 'r+')
Out[45]: [<matplotlib.lines.Line2D at 0x126e99650>]
In [48]: plt.plot( cdf[:,0], 1 - cdf[:,1]/cdf[:,1].max(), 'r+', 1/rcdf[:,0], rcdf[:,1], 'b+' )
In [50]: plt.plot( 1/cdf[:,0], 1 - cdf[:,1]/cdf[:,1].max(), 'r+', rcdf[:,0], rcdf[:,1], 'b+' )
In [64]: plt.plot( 1/cdf[:,0], 1 - cdf[:,1]/cdf[:,1].max(), 'r+-', rcdf[:,0], rcdf[:,1], 'b+-' )
In [52]: ls = get_ls()
In [57]: fast = ls.extra.properties['FASTCOMPONENT'].astype(np.float64)
cy = np.cumsum(fast[:,1], dtype=np.float64) ## cumulative in wavelength land
fcdf = np.vstack([fast[:,0],cy/cy[-1]]).T ## cdf in wavelength
In [112]: np.allclose(fcdf, rcdf)
Out[112]: True
Try duplicating BuildPhysicsTable, fiddly bin averaging:
In [129]: rfast = fast[::-1] # reverse order to be in ascending energy
In [130]: rfast[0]
Out[130]: array([ 799.898, 0. ])
In [146]: x = 1/rfast[:,0] # work in inverse wavelength 1/nm
In [147]: y = rfast[:,1]
(y[1:]+y[:-1])/2 # sum of bins
np.cumsum( (y[1:]+y[:-1])/2 * np.diff(x) )*1e6
In [168]: bcdf = np.vstack( [x[1:], cy/cy[-1]] ).T
In [175]: xcdf = cdf.copy()
In [176]: xcdf[:,1] = xcdf[:,1]/xcdf[:,1].max()
In [180]: np.allclose( xcdf[1:], bcdf )
Out[180]: True
Avoid loosing the bin:
In [185]: bcdf = np.empty( fast.shape )
In [194]: bcdf[0] = 1/fast[-1,0], 0
bcdf[:,0] = x
np.cumsum(ymid*xdif, out=bcdf[1:,1])
bcdf[1:,1] = bcdf[1:,1]/bcdf[1:,1].max()
In [216]: np.allclose(bcdf, xcdf)
Out[216]: True
Lay this down in collada_to_chroma:construct_cdf_energywise
g4daechroma.sh --nogeocache
npysend.sh -t1 -iscintillation -otest
Getting the energywise CDF onto GPU is complicated by chroma wavelength standardization, which does an interpolation to that standard wavelengths. As interpolation requires ascending “x” need to flip order.
Some success with handling energywise cdf, but suspect getting back to front wavelength distrib:
67 def interp_material_property(wavelengths, prop):
68 # note that it is essential that the material properties be
69 # interpolated linearly. this fact is used in the propagation
70 # code to guarantee that probabilities still sum to one.
71
72 ascending = np.all(np.diff(prop[:,0]) >= 0)
73 descending = np.all(np.diff(prop[:,0]) <= 0)
74
75 if ascending:
76 return np.interp(wavelengths, prop[:,0], prop[:,1]).astype(np.float32)
77 elif descending:
78 # the interpolation needs ascending so reverse here, then reverse back after
79 iprop = np.interp(wavelengths, prop[::-1,0], prop[::-1,1]).astype(np.float32)
80 return iprop[::-1].copy()
81 else:
82 assert 0, "needs to be all ascending or all descending "
83 return None
Access test wavelengths:
In [1]: t = ttt_(1)
In [4]: w = t[:,1,3]
Getting some infinites, probably LS material index shift:
In [26]: w[w==np.inf].shape
Out[26]: (598018,)
In [27]: w[w!=np.inf].shape
Out[27]: (2219525,)
The non infinities look like a wavelength distrib:
In [7]: ww=w[w!=np.Inf]
In [9]: plt.hist(ww, bins=100) ## wavelength flipped distribution, maybe need to "1 - cdf"
Succeed to get rid of infinities by establishing order of chroma materials and surfaces to be based on names with pointer address excluded.
After adjusting to use 1/wavelength[::-1] domain reemission_cdf and using kernel sampling that takes that into account,
chroma/cuda/scintillation.h:
194 p.wavelength = sample_reciprocal_cdf(&rng, ss.material->n,
195 ss.material->wavelength0,
196 ss.material->step,
197 ss.material->reemission_cdf);
cf('wavelength', tag=1, typs='gopscintillation opscintillation', log=True )
Scintillation time, almost perfect close match:
In [7]: cf('time', g4s , chs ) ## very long tail
In [30]: cf('time', g4s , chs, range=(0,100))
Position, direction and polarization all almost perfect matches, wavelength needs some attention:
In [32]: cf('3xyzw', g4s, chs, legend=False)
Looking good:
In [3]: cf('3xyzw', tag=1, typs='gopscintillation opscintillation', legend=False, log=True)
delta:~ blyth$ export-
delta:~ blyth$ export-export
delta:~ blyth$ find $DAE_NAME_DYB_CHROMACACHE -name reemission_cdf.npy | grep Gd
/usr/local/env/geant4/geometry/export/DayaBay_VGDX_20140414-1300/g4_00.dae.29c299d81706c62884caf5c3dbdea5c1/chroma_geometry/chroma.detector:Detector:0x11ca48510/unique_materials/003/chroma.geometry:Material:__dd__Materials__GdDopedLS0xc2a8ed0/reemission_cdf.npy
delta:~ blyth$
In [1]: ri = np.load("./chroma.detector:Detector:0x11ca48510/unique_materials/000/chroma.geometry:Material:__dd__Materials__LiquidScintillator0xc2308d0/refractive_index.npy")
In [2]: ri
Out[2]:
array([[ 79.99 , 1.454],
[ 120.023, 1.454],
[ 129.99 , 1.554],
[ 139.984, 1.664],
[ 149.975, 1.783],
[ 159.98 , 1.793],
[ 169.981, 1.554],
[ 179.974, 1.527],
[ 189.985, 1.618],
[ 199.975, 1.618],
[ 300. , 1.526],
[ 404.7 , 1.499],
[ 435.8 , 1.495],
[ 486.001, 1.492],
[ 546.001, 1.486],
[ 589.002, 1.484],
[ 690.701, 1.48 ],
[ 799.898, 1.478]], dtype=float32)
delta:~ blyth$ collada_to_chroma.sh
INFO:env.geant4.geometry.collada.idmap:np.genfromtxt /usr/local/env/geant4/geometry/export/DayaBay_VGDX_20140414-1300/g4_00.idmap
INFO:env.geant4.geometry.collada.idmap:found 685 unique ids
INFO:env.geant4.geometry.collada.g4daenode:idmap exists /usr/local/env/geant4/geometry/export/DayaBay_VGDX_20140414-1300/g4_00.idmap entries 12230
INFO:env.geant4.geometry.collada.g4daenode:index linking DAENode with boundgeom 12230 volumes
INFO:env.geant4.geometry.collada.g4daenode:linking DAENode with idmap 12230 identifiers
INFO:env.geant4.geometry.collada.g4daenode:add_sensitive_surfaces matid __dd__Materials__Bialkali qeprop EFFICIENCY
INFO:env.geant4.geometry.collada.g4daenode:sensitize 684 nodes with matid __dd__Materials__Bialkali and channel_id > 0, uniques 684
INFO:env.geant4.geometry.collada.collada_to_chroma:convert_opticalsurfaces
INFO:env.geant4.geometry.collada.collada_to_chroma:convert_opticalsurfaces creates 44 from 726
WARNING:env.geant4.geometry.collada.collada_to_chroma:setting parent_material to __dd__Materials__Vacuum0xbf9fcc0 as parent is None for node top.0
INFO:env.geant4.geometry.collada.collada_to_chroma:channel_count (nodes with channel_id > 0) : 6888 uniques 684
INFO:env.geant4.geometry.collada.collada_to_chroma:convert_geometry DONE timing_report:
INFO:env.base.timing:timing_report
ColladaToChroma
__init__ : 0.000 1 0.000
convert_flatten : 2.429 1 2.429
convert_geometry_traverse : 4.475 1 4.475
convert_make_maps : 0.000 1 0.000
convert_materials : 0.009 1 0.009
convert_opticalsurfaces : 0.233 1 0.233
INFO:env.geant4.geometry.collada.collada_to_chroma:dropping into IPython.embed() try: cg.<TAB>
Python 2.7.8 (default, Jul 13 2014, 17:11:32)
Type "copyright", "credits" or "license" for more information.
IPython 1.2.1 -- An enhanced Interactive Python.
? -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help -> Python's own help system.
object? -> Details about 'object', use 'object??' for extra details.
In [1]: gdls
Out[1]: <chroma.geometry.Material at 0x10dd0cc50>
In [3]: self = cc
In [5]: collada = self.nodecls.orig
In [6]: collada.materials
Out[6]:
[<Material id=__dd__Materials__PPE0xc12f008 effect=__dd__Materials__PPE_fx_0xc12f008>,
<Material id=__dd__Materials__MixGas0xc21d930 effect=__dd__Materials__MixGas_fx_0xc21d930>,
<Material id=__dd__Materials__Air0xc032550 effect=__dd__Materials__Air_fx_0xc032550>,
<Material id=__dd__Materials__Bakelite0xc2bc240 effect=__dd__Materials__Bakelite_fx_0xc2bc240>,
<Material id=__dd__Materials__Foam0xc558e28 effect=__dd__Materials__Foam_fx_0xc558e28>,
<Material id=__dd__Materials__Aluminium0xc542070 effect=__dd__Materials__Aluminium_fx_0xc542070>,
<Material id=__dd__Materials__Iron0xc542700 effect=__dd__Materials__Iron_fx_0xc542700>,
<Material id=__dd__Materials__GdDopedLS0xc2a8ed0 effect=__dd__Materials__GdDopedLS_fx_0xc2a8ed0>,
<Material id=__dd__Materials__Acrylic0xc02ab98 effect=__dd__Materials__Acrylic_fx_0xc02ab98>,
<Material id=__dd__Materials__Teflon0xc129f90 effect=__dd__Materials__Teflon_fx_0xc129f90>,
<Material id=__dd__Materials__LiquidScintillator0xc2308d0 effect=__dd__Materials__LiquidScintillator_fx_0xc2308d0>,
<Material id=__dd__Materials__Bialkali0xc2f2428 effect=__dd__Materials__Bialkali_fx_0xc2f2428>,
<Material id=__dd__Materials__OpaqueVacuum0xbf5d600 effect=__dd__Materials__OpaqueVacuum_fx_0xbf5d600>,
<Material id=__dd__Materials__Vacuum0xbf9fcc0 effect=__dd__Materials__Vacuum_fx_0xbf9fcc0>,
<Material id=__dd__Materials__Pyrex0xc1005e0 effect=__dd__Materials__Pyrex_fx_0xc1005e0>,
<Material id=__dd__Materials__UnstStainlessSteel0xc5c11e8 effect=__dd__Materials__UnstStainlessSteel_fx_0xc5c11e8>,
<Material id=__dd__Materials__PVC0xc25cfe8 effect=__dd__Materials__PVC_fx_0xc25cfe8>,
<Material id=__dd__Materials__StainlessSteel0xc2adc00 effect=__dd__Materials__StainlessSteel_fx_0xc2adc00>,
<Material id=__dd__Materials__ESR0xbf9f438 effect=__dd__Materials__ESR_fx_0xbf9f438>,
<Material id=__dd__Materials__Nylon0xc3aa360 effect=__dd__Materials__Nylon_fx_0xc3aa360>,
<Material id=__dd__Materials__MineralOil0xbf5c830 effect=__dd__Materials__MineralOil_fx_0xbf5c830>,
<Material id=__dd__Materials__BPE0xc0ad360 effect=__dd__Materials__BPE_fx_0xc0ad360>,
<Material id=__dd__Materials__Ge_680xc2d7e60 effect=__dd__Materials__Ge_68_fx_0xc2d7e60>,
<Material id=__dd__Materials__Co_600xc3cf0c0 effect=__dd__Materials__Co_60_fx_0xc3cf0c0>,
<Material id=__dd__Materials__C_130xc3d0ab0 effect=__dd__Materials__C_13_fx_0xc3d0ab0>,
<Material id=__dd__Materials__Silver0xc3d1370 effect=__dd__Materials__Silver_fx_0xc3d1370>,
<Material id=__dd__Materials__Nitrogen0xc031fd0 effect=__dd__Materials__Nitrogen_fx_0xc031fd0>,
<Material id=__dd__Materials__Water0xc176e30 effect=__dd__Materials__Water_fx_0xc176e30>,
<Material id=__dd__Materials__NitrogenGas0xc17d300 effect=__dd__Materials__NitrogenGas_fx_0xc17d300>,
<Material id=__dd__Materials__IwsWater0xc288f98 effect=__dd__Materials__IwsWater_fx_0xc288f98>,
<Material id=__dd__Materials__ADTableStainlessSteel0xc177178 effect=__dd__Materials__ADTableStainlessSteel_fx_0xc177178>,
<Material id=__dd__Materials__Tyvek0xc246ca0 effect=__dd__Materials__Tyvek_fx_0xc246ca0>,
<Material id=__dd__Materials__OwsWater0xbf90c10 effect=__dd__Materials__OwsWater_fx_0xbf90c10>,
<Material id=__dd__Materials__DeadWater0xbf8a548 effect=__dd__Materials__DeadWater_fx_0xbf8a548>,
<Material id=__dd__Materials__RadRock0xcd2f508 effect=__dd__Materials__RadRock_fx_0xcd2f508>,
<Material id=__dd__Materials__Rock0xc0300c8 effect=__dd__Materials__Rock_fx_0xc0300c8>]
In [7]: collada.materials[7]
Out[7]: <Material id=__dd__Materials__GdDopedLS0xc2a8ed0 effect=__dd__Materials__GdDopedLS_fx_0xc2a8ed0>
In [8]: collada.materials[7].extra
Out[8]: <MaterialProperties keys=['SLOWTIMECONSTANT', 'GammaFASTTIMECONSTANT', 'ReemissionSLOWTIMECONSTANT', 'REEMISSIONPROB', 'AlphaFASTTIMECONSTANT', 'ReemissionFASTTIMECONSTANT', 'SLOWCOMPONENT', 'YIELDRATIO', 'FASTCOMPONENT', 'RINDEX', 'NeutronFASTTIMECONSTANT', 'ReemissionYIELDRATIO', 'RAYLEIGH', 'NeutronYIELDRATIO', 'GammaYIELDRATIO', 'SCINTILLATIONYIELD', 'AlphaYIELDRATIO', 'RESOLUTIONSCALE', 'GammaSLOWTIMECONSTANT', 'AlphaSLOWTIMECONSTANT', 'NeutronSLOWTIMECONSTANT', 'ABSLENGTH', 'FASTTIMECONSTANT'] >
In [9]:
In [11]: collada.materials[7].extra.properties
Out[11]:
{'ABSLENGTH': array([[ 79.9898, 0.001 ],
[ 120.0235, 0.001 ],
[ 199.9746, 0.001 ],
...,
[ 897.916 , 328.4 ],
[ 898.8925, 306.2 ],
[ 899.8711, 299.6 ]]),
'AlphaFASTTIMECONSTANT': array([[ 0.0012, 1. ],
[-0.0012, 1. ]]),
'AlphaSLOWTIMECONSTANT': array([[ 0.0012, 35. ],
[ -0.0012, 35. ]]),
'AlphaYIELDRATIO': array([[ 0.0012, 0.65 ],
[-0.0012, 0.65 ]]),
'FASTCOMPONENT': array([[ 79.9898, 0. ],
[ 120.0235, 0. ],
[ 199.9746, 0. ],
...,
[ 599.0011, 0.0017],
[ 600.0012, 0.0018],
[ 799.8984, 0. ]]),
'FASTTIMECONSTANT': array([[ 0.0012, 3.64 ],
[-0.0012, 3.64 ]]),
'GammaFASTTIMECONSTANT': array([[ 0.0012, 7. ],
[-0.0012, 7. ]]),
'GammaSLOWTIMECONSTANT': array([[ 0.0012, 31. ],
[ -0.0012, 31. ]]),
'GammaYIELDRATIO': array([[ 0.0012, 0.805 ],
[-0.0012, 0.805 ]]),
'NeutronFASTTIMECONSTANT': array([[ 0.0012, 1. ],
[-0.0012, 1. ]]),
'NeutronSLOWTIMECONSTANT': array([[ 0.0012, 34. ],
[ -0.0012, 34. ]]),
'NeutronYIELDRATIO': array([[ 0.0012, 0.65 ],
[-0.0012, 0.65 ]]),
'RAYLEIGH': array([[ 79.9898, 850. ],
[ 120.0235, 850. ],
[ 199.9746, 850. ],
...,
[ 589.8394, 170000. ],
[ 699.9223, 300000. ],
[ 799.8984, 500000. ]]),
'REEMISSIONPROB': array([[ 79.9898, 0.4 ],
[ 120.0235, 0.4 ],
[ 159.9797, 0.4 ],
...,
[ 575.8273, 0.0587],
[ 712.6064, 0. ],
[ 799.8984, 0. ]]),
'RESOLUTIONSCALE': array([[ 0.0012, 1. ],
[-0.0012, 1. ]]),
'RINDEX': array([[ 79.9898, 1.4536],
[ 120.0235, 1.4536],
[ 129.9898, 1.5545],
...,
[ 589.0016, 1.4842],
[ 690.7008, 1.48 ],
[ 799.8984, 1.4781]]),
'ReemissionFASTTIMECONSTANT': array([[ 0.0012, 1.5 ],
[-0.0012, 1.5 ]]),
'ReemissionSLOWTIMECONSTANT': array([[ 0.0012, 1.5 ],
[-0.0012, 1.5 ]]),
'ReemissionYIELDRATIO': array([[ 0.0012, 1. ],
[-0.0012, 1. ]]),
'SCINTILLATIONYIELD': array([[ 0.0012, 11522. ],
[ -0.0012, 11522. ]]),
'SLOWCOMPONENT': array([[ 79.9898, 0. ],
[ 120.0235, 0. ],
[ 199.9746, 0. ],
...,
[ 599.0011, 0.0017],
[ 600.0012, 0.0018],
[ 799.8984, 0. ]]),
'SLOWTIMECONSTANT': array([[ 0.0012, 12.2 ],
[ -0.0012, 12.2 ]]),
'YIELDRATIO': array([[ 0.0012, 0.86 ],
[-0.0012, 0.86 ]])}
In [12]:
In [12]: collada.materials[7].extra.properties['SLOWCOMPONENT']
Out[12]:
array([[ 79.9898, 0. ],
[ 120.0235, 0. ],
[ 199.9746, 0. ],
...,
[ 599.0011, 0.0017],
[ 600.0012, 0.0018],
[ 799.8984, 0. ]])
In [13]: collada.materials[7].extra.properties['FASTCOMPONENT']
Out[13]:
array([[ 79.9898, 0. ],
[ 120.0235, 0. ],
[ 199.9746, 0. ],
...,
[ 599.0011, 0.0017],
[ 600.0012, 0.0018],
[ 799.8984, 0. ]])
In [14]: collada.materials[7].extra.properties['REEMISSIONPROB']
Out[14]:
array([[ 79.9898, 0.4 ],
[ 120.0235, 0.4 ],
[ 159.9797, 0.4 ],
...,
[ 575.8273, 0.0587],
[ 712.6064, 0. ],
[ 799.8984, 0. ]])
In [15]:
In [15]: np.allclose( collada.materials[7].extra.properties['SLOWCOMPONENT'], collada.materials[7].extra.properties['FASTCOMPONENT'] )
Out[15]: True
In [15]: _gdls = gdls()
In [18]: _gdls.__class__
Out[18]: collada.material.Material
In [21]: slow = _gdls.extra.properties['SLOWCOMPONENT']
In [22]: plt.scatter(slow[:,0],slow[:,1])
Out[22]: <matplotlib.collections.PathCollection at 0x115e406d0>
In [23]: plt.show()
Wide range, but very few entries at extremes and near zero anyhow, all action in middle:
In [20]: _gdls.extra.properties['SLOWCOMPONENT']
Out[20]:
array([[ 79.99 , 0. ],
[ 120.023, 0. ],
[ 199.975, 0. ],
[ 330. , 0.006],
[ 331. , 0.006],
[ 332. , 0.005],
[ 333. , 0.005],
...
[ 598.001, 0.002],
[ 599.001, 0.002],
[ 600.001, 0.002],
[ 799.898, 0. ]])
In [24]: slow[:,0].min()
Out[24]: 79.989835277575907
In [25]: slow[:,0].max()
Out[25]: 799.89835277575912
Chopping the extremes:
In [28]: plt.scatter(slow[10:-10,0],slow[10:-10,1])
Out[28]: <matplotlib.collections.PathCollection at 0x124b8f110>
In [29]: plt.show()
The wide range feeds forward into chroma:
In [33]: cg = chroma_geometry()
In [37]: cg.unique_materials[0].name
Out[37]: '__dd__Materials__LiquidScintillator0xc2308d0'
In [38]: cls = cg.unique_materials[0]
In [40]: cls.reemission_cdf.shape
Out[40]: (275, 2)
In [41]: slow.shape
Out[41]: (275, 2)
In [44]: np.allclose( cls.reemission_cdf[:,0], slow[:,0] )
Out[44]: True
chroma/chroma/geometry.py:
25 # all material/surface properties are interpolated at these
26 # wavelengths when they are sent to the gpu
27 standard_wavelengths = np.arange(60, 810, 20).astype(np.float32)
28
Hmm thats pretty coarse, this explains the generated scintillation wavelength distrib.