Debug ====== Debugging optical photon propagation using NumPy + ipython ------------------------------------------------------------- Using *num_hits* to debug an optical propagation is hopeless. You need to enable photon step-by-step recording and save the corresponding arrays in NumPy .npy files. Then you can examine the parameters of all the photons including history flags at every step point of their propagation (up to a configured maximum number of step points) from interactive ipython sessions.:: ipython> a = np.load("/path/to/photon.npy") You could then make plots drawing the paths of the photons. I recommend pyvista if your want to do that A convenient way to install pyvista is to use anaconda. The more commonly used matplotlib python plotting library is not good with 3D plotting or large data sets. To save the arrays you need to:: export OPTICKS_EVENT_MODE=StandardFullDebug # configure step point recording export G4CXOpticks__simulate_saveEvent=1 # enable saveEvent from G4CXOpticks::simulate # optionally enable logging in relevant classes export G4CXOpticks=INFO export SEventConfig=INFO What OPICKS_EVENT_MODE does ------------------------------- To see how OPTICKS_EVENT_MODE works look at:: sysrap/SEventConfig.hh sysrap/SEventConfig.cc Especially:: 324 int SEventConfig::Initialize() // static 325 { 326 const char* mode = EventMode(); 327 int maxbounce = MaxBounce(); 328 329 if(strcmp(mode, Default) == 0 ) 330 { 331 SetCompMaskAuto() ; // comp set based on Max values 332 } 333 else if(strcmp(mode, StandardFullDebug) == 0 ) 334 { 335 SEventConfig::SetMaxRecord(maxbounce+1); 336 SEventConfig::SetMaxRec(maxbounce+1); 337 SEventConfig::SetMaxSeq(maxbounce+1); 338 SEventConfig::SetMaxPrd(maxbounce+1); 339 340 // since moved to compound sflat/stag so MaxFlat/MaxTag should now either be 0 or 1, nothing else 341 SEventConfig::SetMaxTag(1); 342 SEventConfig::SetMaxFlat(1); 343 SetCompMaskAuto() ; // comp set based on Max values 344 } 345 else 346 { 347 LOG(fatal) << "mode [" << mode << "] IS NOT RECOGNIZED " ; 348 LOG(fatal) << " options : " << Default << "," << StandardFullDebug ; 349 assert(0); 350 } Possible arrays saved by sysrap/SEvt ---------------------------------------- The default is not saving any arrays. Array saving is regarded as purely a debugging activity that should not be done in production, as it will greatly slow down the simulation and write many large files. Configure what arrays to save and their limits with *SEventConfig* via envvars or static methods. +------------+------------------+-------------------------------------------------------------------------------+ | array | shape | notes | +============+==================+===============================================================================+ | *genstep* | (num_gs,6,4) | parameters of Cerenkov/Scintillation/Torch/... generation | +------------+------------------+-------------------------------------------------------------------------------+ | *photon* | (num_ph,4,4) | sysrap/sphoton.h : final photon params : position, time, mom, pol, flags | +------------+------------------+-------------------------------------------------------------------------------+ | *record* | (num_ph,16,4,4) | sysrap/sphoton.h : step point records (configurable num points) | +------------+------------------+-------------------------------------------------------------------------------+ | *hit* | (num_ht,4,4) | sysrap/sphoton.h : selection of *photon* | +------------+------------------+-------------------------------------------------------------------------------+ | *rec* | (num_ph,16,4,4) | compressed rec (no longer used, use *record* instead) | +------------+------------------+-------------------------------------------------------------------------------+ | *seq* | (num_ph,2) | sysrap/sseq.h photon level step-by-step history and material recording | +------------+------------------+-------------------------------------------------------------------------------+ | *prd* | | | +------------+------------------+-------------------------------------------------------------------------------+ SEventConfigTest ------------------ The **SEventConfigTest** binary can be used to dump the configuration. :: epsilon:sysrap blyth$ SEventConfigTest test_EstimateAlloc@20: SEventConfig::Desc OPTICKS_EVENT_MODE EventMode : Default OPTICKS_RUNNING_MODE RunningMode : 0 RunningModeLabel : SRM_DEFAULT OPTICKS_G4STATE_SPEC G4StateSpec : 1000:38 G4StateSpecNotes : 38=2*17+4 is appropriate for MixMaxRng OPTICKS_G4STATE_RERUN G4StateRerun : -1 OPTICKS_MAX_GENSTEP MaxGenstep : 1000000 OPTICKS_MAX_PHOTON MaxPhoton : 1000000 OPTICKS_MAX_SIMTRACE MaxSimtrace : 1000000 MaxCurandState : 1000000 OPTICKS_MAX_BOUNCE MaxBounce : 9 OPTICKS_MAX_RECORD MaxRecord : 0 OPTICKS_MAX_REC MaxRec : 0 OPTICKS_MAX_SEQ MaxSeq : 0 OPTICKS_MAX_PRD MaxPrd : 0 OPTICKS_MAX_TAG MaxTag : 0 OPTICKS_MAX_FLAT MaxFlat : 0 OPTICKS_HIT_MASK HitMask : 64 HitMaskLabel : SD OPTICKS_MAX_EXTENT MaxExtent : 1000 OPTICKS_MAX_TIME MaxTime : 10 OPTICKS_RG_MODE RGMode : 2 RGModeLabel : simulate OPTICKS_COMP_MASK CompMask : 262 CompMaskLabel : genstep,photon,hit OPTICKS_OUT_FOLD OutFold : $DefaultOutputDir OPTICKS_OUT_NAME OutName : - OPTICKS_PROPAGATE_EPSILON PropagateEpsilon : 0.0500 OPTICKS_INPUT_PHOTON InputPhoton : - test_EstimateAlloc@25: al.desc - epsilon:sysrap blyth$ Changing the OPTICKS_EVENT_MODE envvar to "StandardFullDebug" has a large effect on the config:: epsilon:sysrap blyth$ OPTICKS_EVENT_MODE=StandardFullDebug SEventConfigTest test_EstimateAlloc@20: SEventConfig::Desc OPTICKS_EVENT_MODE EventMode : StandardFullDebug OPTICKS_RUNNING_MODE RunningMode : 0 RunningModeLabel : SRM_DEFAULT OPTICKS_G4STATE_SPEC G4StateSpec : 1000:38 G4StateSpecNotes : 38=2*17+4 is appropriate for MixMaxRng OPTICKS_G4STATE_RERUN G4StateRerun : -1 OPTICKS_MAX_GENSTEP MaxGenstep : 1000000 OPTICKS_MAX_PHOTON MaxPhoton : 1000000 OPTICKS_MAX_SIMTRACE MaxSimtrace : 1000000 MaxCurandState : 1000000 OPTICKS_MAX_BOUNCE MaxBounce : 9 OPTICKS_MAX_RECORD MaxRecord : 10 OPTICKS_MAX_REC MaxRec : 10 OPTICKS_MAX_SEQ MaxSeq : 10 OPTICKS_MAX_PRD MaxPrd : 10 OPTICKS_MAX_TAG MaxTag : 1 OPTICKS_MAX_FLAT MaxFlat : 1 OPTICKS_HIT_MASK HitMask : 64 HitMaskLabel : SD OPTICKS_MAX_EXTENT MaxExtent : 1000 OPTICKS_MAX_TIME MaxTime : 10 OPTICKS_RG_MODE RGMode : 2 RGModeLabel : simulate OPTICKS_COMP_MASK CompMask : 274814 CompMaskLabel : genstep,photon,record,rec,seq,prd,hit,tag,flat,aux OPTICKS_OUT_FOLD OutFold : $DefaultOutputDir OPTICKS_OUT_NAME OutName : - OPTICKS_PROPAGATE_EPSILON PropagateEpsilon : 0.0500 OPTICKS_INPUT_PHOTON InputPhoton : - test_EstimateAlloc@25: al.desc - epsilon:sysrap blyth$ Saving Photon Propagations into NumPy arrays ---------------------------------------------- To see what G4CXOpticks__simulate_saveEvent is doing look at g4cx/G4CXOpticks.cc simulate method. The directory where the numpy arrays is saved is based on your executable name and the event index you set with:: SEvt::SetIndex(eventid); Enable logging in SEvt to see what it is:: export SEvt=INFO Opticks has lots of python machinery for loading and presenting such NumPy .npy arrays in the "ana" directory and all over the place. However it is better to examine them manually using ipython to begin with, because most people will need to improve their NumPy+ipython skills to make best use of this debugging info and to be able to understand the python machinery. Debugging Lack of Hits --------------------------- When none of your photons have flagmask SURFACE_DETECT you will get no hits. You will still have gensteps, photons and records without having hits. If your photon propagation histories are as expected and you are still not getting hits then your problem is probably that the geometry translation did not notice your sensitive detectors somehow OR that you do not have any sensitive detectors. A common cause of this is loading a geometry from GDML and not reinstating the association. *hits* are the subset of *photons* with flagmask matching the hitmask (default SURFACE_DETECT) so when you get no hits it means that none of your photons .flagmask has all the bits of the hitmask set. You can of course select the hits array from the photons array using one line of NumPy, but that will just match with NumPy what the C++/CUDA would do. You can learn about the mechanics of hit selection in:: ~/opticks/notes/mechanics_of_hit_selection.rst https://bitbucket.org/simoncblyth/opticks/src/master/notes/mechanics_of_hit_selection.rst Checking photon propagation histories ---------------------------------------- So in ipython:: import numpy as np r = np.load("/tmp/ami/opticks/GEOM/example_pet/ALL/007/record.npy”) r[:,:5,0] # (pos, time) of first 5 step point of all photons r[0, :5, 0] # (pos, time) of first 5 step points of the first photon (slot 0) So to start debugging I would look at the sequence of positions, times, momentum directions and polarizations and step point flags to see if they are what I would expect. ~/opticks/ana/tests/check.sh : setup environment for NumPy debug -------------------------------------------------------------------- :: #!/bin/bash -l usage(){ cat << EOU check.sh ========== Setup environment: PYTHONPATH=$HOME allows python scripts to import opticks python machinery eg with "from opticks.ana.fold import Fold" CFBASE=$HOME/.opticks/GEOM/J004 configures where to load geometry from FOLD=$CFBASE/G4CXSimulateTest/ALL configures the directory to load event arrays from, the directory is up to the user To a large degree the directory positions of geometry and event files are controlled by the user. However the example of versioning a geometry name "J004" and keeping event folders within the corresponding geometry folder is a good one to follow as it is important to retain the connection between event data and the geometry used to create the event data. EOU } export PYTHONPATH=$HOME export CFBASE=$HOME/.opticks/GEOM/J004 export FOLD=$CFBASE/G4CXSimulateTest/ALL ${IPYTHON:-ipython} --pdb -i check.py ~/opticks/ana/tests/check.py : python basic example ------------------------------------------------------ :: epsilon:tests blyth$ cat check.py #!/usr/bin/env python import numpy as np from opticks.ana.fold import Fold from opticks.ana.p import * from opticks.ana.histype import HisType if __name__ == '__main__': f = Fold.Load(symbol="f") print(repr(f)) p = f.photon r = f.record s = f.seq h = f.hit Debugging Geometry ---------------------- To debug geometry issues you need to have some familiarity with the translation. It starts from G4CXOpticks::setGeometry. :: 210 void G4CXOpticks::setGeometry(const G4VPhysicalVolume* world ) 211 { 212 LOG(LEVEL) << " G4VPhysicalVolume world " << world ; 213 assert(world); 214 wd = world ; 215 216 sim = SSim::Create(); 217 stree* st = sim->get_tree(); 219 tr = U4Tree::Create(st, world, SensorIdentifier ) ; 220 221 222 // GGeo creation done when starting from a gdml or live G4, still needs Opticks instance 223 Opticks::Configure("--gparts_transform_offset --allownokey" ); 224 225 GGeo* gg_ = X4Geo::Translate(wd) ; 226 setGeometry(gg_); 227 } X4Geo::Translate old way with loads of code, entire extg4 package U4Tree::Create is a simpler approach to translation that I am starting to develop which is aiming to go directly The translation code is still in flux with both old and new approaches in use and an entire geometry model too many.:: . extg4 CSG_GGeo Geant4 ----> GGeo -------> CSG The CSG_GGeo package translates the GGeo geometry model into CSG which gets upload to GPU. While I dont suggest you try to understand the geometry translation in detail, you need some familiarity with how the sensitive detectors get translated in order to debug your issue. The best way to start debugging geometry is to persist it by rerunning with:: export G4CXOpticks=INFO export G4CXOpticks__setGeometry_saveGeometry=$HOME/.opticks/GEOM/$GEOM The above envvar configures directory to save the geometry. Then you can run small executables/scripts which load various parts of the persisted geometry and run tests. One, of many, of such tests is sysrap/tests/stree_test.sh Build and use that:: cd ~/opticks/sysrap/tests ./stree_test.sh build ## builds stree_test binary ./stree_test.sh run ## these load the geometry into C++ and run tests against it ## one of the tests dumps sensor info ./stree_test.sh ana ## this loads the same geometry into ipython ## and run tests against it