FAQ:How do I find if a point is in a particular volume?
Jump to navigation
Jump to search
| Offline Documentation: [Offline Category] [FAQ] [Howto] [Reference] [Manual] |
| Offline Documentation |
| This article is part of the Offline Documentation |
| Other pages... |
|
Offline Category |
| See also... |
Short answer
You ask it.
Long answer
First, by "volume" we must mean Detector Element (DE), ie, the oil volume in AD#1 at the far site, but see here, too.
Asking if a volume contains a point
Once you get the desired DE you can ask it if it contains a point like:
// C++
Gaudi::XYZPoint thePoint = ...;
DetectorElement *theDE = ...;
if (theDE->isInside(thePoint)) {
do_something();
}
# Python
import PyCintex
Gaudi = PyCintex.makeNamespace('Gaudi')
point = Gaudi.XYZPoint(pt[0],pt[1],pt[2])
de = ...;
if de.isInside(point):
do_something()
Asking what volume contains a point
Maybe a volume contains a point but it may be further contained by a daughter volume. You can ask what volume contains a point by going through the DE's Geometry Info (GI) like:
// C++ int theLevel = -1; std::string path = theDE->geometry()->belongsToPath(thePoint,theLevel
# Python level = -1; path = de.geometry().belongsToPath(point,level)
In both cases the "level" argument is optional with a default value of -1. From the IGeometryInfo.h comments:
* - if level = 0 - no search, return the name of current level * - if level < 0 - perform search up to the most deepest level * - if level > 0 - perform search up to not more then "level" levels;
The return value is suitable for looking up the corresponding DE in the detector data store.
What is the material of the volume that contains a point ?
To accurately evaluate the material, logical volume or placed volume that contains a given point, you need to delve beyond the DetectorElements into the placed volume path using the belongsTo() method of ILVolume as shown in this lengthy Python example from sumMaterial.py . Detector elements (/dd/Structure) do not, by design, have the full geometry information (/dd/Geometry) of the placed volumes.
evt = self.evtSvc()
simhdr = evt['/Event/Sim/SimHeader']
if debugme>3 : print "SimHeader: ", simhdr
if simhdr == None:
print "No SimHeader in this ReadOut. Skip."
return SUCCESS
execn = simhdr.execNumber()
print "++++++++++++++++++++++++++++= Executing ",self.name()," exec#",execn
# geometry info
self.target_de_name = '/dd/Structure/AD/db-ade1'
# Get underlying DE object
de = self.getDet(self.target_de_name)
if not de :
print 'Failed to get DE ',self.target_de_name
return FAILURE
Gaudi = PyCintex.makeNamespace("Gaudi")
nstep = 10000
# loop over vertices in particle history (historian info)
phis = simhdr.particleHistory()
vtx = phis.vertexVector()
nvtx= vtx.size()
print " found a total of ",nvtx," vertices. Will use ",nstep," steps between vertices."
lastPos = None
mNames = [] # material names
mDX = [] # distance travelled in material
for i in range(0,nvtx-1):
j = i + 1
vsec = vtx[i].secondaries()
vproc = vtx[i].process().name()
pdg = vtx[i].track().track().particle()
energy = vtx[i].totalEnergy()
gp1 = Gaudi.XYZPoint( vtx[i].position() )
gp2 = Gaudi.XYZPoint( vtx[j].position() )
dist = math.sqrt( math.pow(gp2.x()-gp1.x(),2) + math.pow(gp2.y()-gp1.y(),2) + math.pow(gp2.z()-gp1.z(),2) )
dx = (gp2.x()-gp1.x())/dist
dy = (gp2.y()-gp1.y())/dist
dz = (gp2.z()-gp1.z())/dist
x = gp1.x()
y = gp1.y()
z = gp1.z()
delta = dist/float(nstep)
### get global position and then find local position in 'target' detector element
### compute last step size
gloPos = Gaudi.XYZPoint( vtx[i].position() )
for istep in range( 0, nstep+1 ) :
r = float(istep)
gloPos.SetXYZ( x + dx*delta*r, y + dy*delta*r, z + dz*delta*r)
#gloPos = Gaudi.XYZPoint( vtx[i].position() )
locPos = de.geometry().toLocal( gloPos )
# get path of innermost detector element
level= -1
path = de.geometry().belongsToPath( gloPos, level )
# get detector element
what = self.getDet(path)
whatname = what.name()
if debugme>1: print "vtx#",i,"innermost DE",whatname
# get GeometryInfoPlus object for this detector element
gip = what.geometry()
# get pointer to assocated logical volume of detector element
lvName = "NO LOGICAL VOLUME"
lvMaterial = "NONE"
if gip.hasLVolume() :
lvp = gip.lvolume()
lvName = lvp.name()
lvMaterial = lvp.materialName()
if debugme>1: print "vtx#",i,"lVol",lvName,"mat",lvMaterial
# local position in current DetElem
whatPos = what.geometry().toLocal( gloPos )
# get vector of the placed volume path in current DetElem,
# the report the innermost placed volume (pVol)
pvpath = gbl.ILVolume.PVolumePath()
level = -1
sc = lvp.belongsTo( whatPos, level, pvpath)
if sc.isSuccess() :
if debugme>1:
print "vtx#",i,"successfully got pVol path vector. vector size is",pvpath.size()," Here is vector content"
for ipv in range(0, pvpath.size() ):
print "vtx#",i,"ipv",ipv,"path",pvpath[ipv].name(),\
"lVol",pvpath[ipv].lvolume().name(),"mat",pvpath[ipv].lvolume().materialName()
if pvpath.size() > 0 :
ipv = pvpath.size()-1
ll = pvpath[ipv].lvolume()
path = pvpath[ipv].name()
lvName = ll.name()
lvMaterial = ll.materialName()
if debugme>0:
print "vtx#",i,"istep",istep,"pdg",pdg,"energy",energy,"DE/pVol ",path,"lVol",lvName,"material",lvMaterial,\
"XYZglobal(mm)=",gloPos.x()/units.mm," ",gloPos.y()/units.mm," ",gloPos.z()/units.mm,\
"XYZlocal(mm)=",locPos.x()/units.mm," ",locPos.y()/units.mm," ",locPos.z()/units.mm,\
"step(mm)",delta/units.mm
if mNames.count(lvMaterial) == 0 :
mNames.append(lvMaterial)
mDX.append(0.)
im = mNames.index(lvMaterial)
mDX[im] += delta
for im in range( 0, len(mNames) ) :
print "im",im,"material",mNames[im],"distance",mDX[im]
# info from unobservatble stats
statshdr = simhdr.unobservableStatistics()
if statshdr :
stats = statshdr.stats()
if stats :
print "pathlen_acrylic", stats["pathlen_acrylic"].sum()
print "pathlen_MO", stats["pathlen_MO"].sum()
print "pathlen_LS", stats["pathlen_LS"].sum()
print "pathlen_GdLS", stats["pathlen_GdLS"].sum()
print "pathlen_tot", stats["pathlen_tot"].sum()
| Offline Software Documentation: [Offline Categories] [FAQ] [Offline Faq Category] |