DetSim Parallel Cerenkov and Scintillation ============================================ * Can the photon loops be parallelized by moving generation onto the GPU ? * This would largely avoid transport overheads. 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. Cerenkov Refs --------------- * http://large.stanford.edu/courses/2014/ph241/alaeian2/ * http://math.ucr.edu/home/baez/physics/Relativity/SpeedOfLight/cherenkov.html * https://thespectrumofriemannium.wordpress.com/tag/tamm-frank-formula/ Rejection Sampling ------------------- * http://www.youtube.com/watch?v=wRdjeCtc8d0&feature=youtube_gdata_player Check Generation Consistency ----------------------------- Three files to check for each process (Scintillation and Cerenkov): * CERENKOV, the step file * GOPCERENKOV, g4 optical photons generated from those steps * OPCERENKOV, chroma optical photons generated from the steps on GPU * SCINTILLATION, the step file * GOPSCINTILLATION, g4 optical photons generated from those steps * OPSCINTILLATION, chroma optical photons generated from the steps on GPU 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. Revisit With New Type Names ----------------------------- :: 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) Count Checking ---------------- * Cerenkov counts are matching only as ApplyWaterQE killing (via continue in photon loop) is still disabled 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 Cerenkov ---------- tee up the arrays ~~~~~~~~~~~~~~~~~~~ :: In [1]: chc, g4c, tst = NPY.mget(1,"opcerenkov","gopcerenkov","test") wavelength ~~~~~~~~~~~~ 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 * need to split test into testcerenkov and testscintillation, otherwise write stomping potential ? 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 time ~~~~~~ :: In [9]: cf('time', tag=1, typs="gopcerenkov opcerenkov test") * with ApplyWaterQE killing enabled * very closely matched up to 18ns, beyond that much less g4 * without ApplyWaterQE * almost perfect match xyz pos,dir,pol ~~~~~~~~~~~~~~~~~~ :: In [9]: cf('3xyz', g4c, chc) * with ApplyWaterQE killing enabled * pos : clear spatial discrepancy, less at extremes of x and y * without ApplyWaterQE * pos : almost perfect * dir : vgood agreement, except that chroma spikes are more spiky * pol : same as dir with chroma spikes more spiky investigate cerenkov wavelength ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ `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]: 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]: In [72]: cg.unique_materials[7].name Out[72]: '__dd__Materials__Bialkali0xc2f2428' cerenkov review ~~~~~~~~~~~~~~~~~ TODO: settle on standard wavelenth range to match G4 better:: In [4]: cf('3xyzw', tag=1, typs='gopcerenkov opcerenkov', legend=False, log=True) Scintillation -------------- tee up the arrays ~~~~~~~~~~~~~~~~~~~ :: In [8]: g4s, chs = NPY.mget(1, "gopscintillation opscintillation") wavelength ~~~~~~~~~~~ :: 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 * not quite, just a case of coarse interpolation `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 * what to do about that ? * tighten the range to a more relevant one, and reduce bin size to keep roughly the same number of bins * reduce bin size * variable bin size ? bad performance impact presumably * could use a coarse and a fine closer look at scintillation wavelength ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * chroma has sharp cutoff at 600nm and less clear drop down at around 350nm 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) investigate G4 scintillation wavelength ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: 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 ; nGetEnergy(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]: [] 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 compare cdfs ~~~~~~~~~~~~~~~~ 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]: [] In [15]: plt.show() :: In [41]: plt.plot( 1/rcdf[:,0], rcdf[:,1], 'b+') Out[41]: [] In [45]: plt.plot( cdf[:,0], 1 - cdf[:,1]/cdf[:,1].max(), 'r+') Out[45]: [] :: 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 test with energywise cdf ~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: 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 ) time ~~~~~~~~ Scintillation time, almost perfect close match:: In [7]: cf('time', g4s , chs ) ## very long tail In [30]: cf('time', g4s , chs, range=(0,100)) xyz pos,dir,pol ~~~~~~~~~~~~~~~~~ Position, direction and polarization all almost perfect matches, wavelength needs some attention:: In [32]: cf('3xyzw', g4s, chs, legend=False) scintillation review ~~~~~~~~~~~~~~~~~~~~~~ Looking good:: In [3]: cf('3xyzw', tag=1, typs='gopscintillation opscintillation', legend=False, log=True) Properties ---------- :: 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$ Lookups for Cerenkov --------------------- :: 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) Material Properties for Scintillation/Cerenkov GPU generation --------------------------------------------------------------- :: 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. 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]: In [3]: self = cc In [5]: collada = self.nodecls.orig In [6]: collada.materials Out[6]: [, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ] In [7]: collada.materials[7] Out[7]: In [8]: collada.materials[7].extra Out[8]: 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 Wavelength Ranges from G4 to Chroma ------------------------------------- :: 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]: 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]: 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.