FAQ:How do I find if a point is in a particular volume?

From Daya Bay
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
How Tos
Getting Started
Software Installation

See also...

General Help on this Wiki

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)) {
# Python
import PyCintex
Gaudi = PyCintex.makeNamespace('Gaudi')
point = Gaudi.XYZPoint(pt[0],pt[1],pt[2])
de = ...;
if de.isInside(point):

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(),\
                    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,\
                if mNames.count(lvMaterial) == 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]