In This Package:

pmtCoords.py

Go to the documentation of this file.

#!/usr/bin/env python
# An example script showing how to access the PMT and detector geometry
#
# The script starts up a Gaudi application that does nothing but
# initialize the Gaudi system.
# Once the system is initialized, we have access to Gaudi Services,
# such as the PMT Geometry.
#
# dandwyer@caltech.edu 2008-07-31

# Load ROOT.  Be sure to do this before loading Gaudi.
import ROOT
# Need CLHEP classes for geometry
ROOT.gSystem.Load("libCLHEPRflx")
# GaudiPython overwrites ROOT, so save a reference
_ROOT = ROOT

# Configure the detector geometry
import xmldetdesc
xmldetdesc.config()

# Get units information
from GaudiKernel import SystemOfUnits as units

# Set configuration for the PMT Geometry Service
from Gaudi.Configuration import *
conf = ApplicationMgr()
from DetHelpers.DetHelpersConf import PmtGeomInfoSvc
pgisvc = PmtGeomInfoSvc("PmtGeomInfoSvc")
#pgisvc.OutputLevel = 1
pgisvc.StreamItems = [ "/dd/Structure/DayaBay" ]

# Add a dummy algorithm to initialize gaudi 
from GaudiPython.GaudiAlgs import GaudiAlgo
from GaudiPython import *
class InitAlg(GaudiAlgo):
    """A dummy algorithm that does nothing but provide access to gaudi
    services"""
    def __init__(self,name):
        GaudiAlgo.__init__(self,name)
        print "Making InitAlg",name

    def initialize(self):
        status = GaudiAlgo.initialize(self)
        return status

    def execute(self):
        return SUCCESS

# Prepare Application
from GaudiPython import AppMgr
app = AppMgr()

app.EvtSel = "NONE"
app.ExtSvc += [ "PmtGeomInfoSvc" ]
initalg = InitAlg("InitAlg")
app.setAlgorithms( [ initalg ] )

# Run Application (which does nothing but initialize gaudi)
app.run(1)

# GaudiPython overwrites ROOT, so restore from saved reference
ROOT = _ROOT
del _ROOT

# Get the pmt geometry service
#pmtSvc = initalg.pgis
pmtSvc = initalg.svc("IPmtGeomInfoSvc","PmtGeomInfoSvc")

#################################################################
# From here on out, we treat this like a normal ROOT script.

# PMT ID for DYB Site AD1, Ring 1, Column 1
site = 0x01 # 0x01 (DYB), 0x02 (LA), or 0x04 (FAR)
ad =   1    # 1-2 (near sites), 1-4 (far sites)
ring = 1    # 1-8
col =  1    # 1-24

pmtid =  (site<<24) | (ad<<16) | (ring<<8) | col

# PMT Path for same AD, but Ring 8, Column 1
pmtPath = "/dd/Structure/AdPmtStructure/db-ad/db-ad1-ring8-column1"

# Get the geometry info for the two pmts
pmtGeom_1_1 = pmtSvc.get(pmtid)
pmtGeom_8_1 = pmtSvc.get(pmtPath)

# Get position relative to AD center
pmtPos_1_1 = pmtGeom_1_1.localPosition()
pmtPos_8_1 = pmtGeom_8_1.localPosition()
print "="*60
print "\n\n"
print "Daya Bay Site, AD1"
print ""
print "In AD Coordinate System:"
print "  PMT (Ring 1, Column 1) X Y Z [mm]: ",pmtPos_1_1.x(),pmtPos_1_1.y(),pmtPos_1_1.z()
print "  PMT (Ring 8, Column 1) X Y Z [mm]: ",pmtPos_8_1.x(),pmtPos_8_1.y(),pmtPos_8_1.z()

# Get position relative in global coordinates of the entire experiment
pmtGPos_1_1 = pmtGeom_1_1.globalPosition()
pmtGPos_8_1 = pmtGeom_8_1.globalPosition()
print ""
print "In Global Coordinate System of the entire experiment:"
print "  PMT (Ring 1, Column 1) X Y Z [mm]: ",pmtGPos_1_1.x(),pmtGPos_1_1.y(),pmtGPos_1_1.z()
print "  PMT (Ring 8, Column 1) X Y Z [mm]: ",pmtGPos_8_1.x(),pmtGPos_8_1.y(),pmtGPos_8_1.z()

# Get the AD position in Global Coordinates
adDetector = pmtGeom_1_1.parentDetector()
adName = adDetector.name()
# Detector origin in local coords: 0,0,0 by definition
origin = gbl.Gaudi.XYZPoint(0,0,0)
# Transform to global coordinates
adPos = adDetector.geometry().toGlobal(origin)
print ""
print "Daya Bay Site AD1"
print "  Path:", adName
print "  Global Coordinates X Y Z [mm]:",adPos.x(),adPos.y(),adPos.z()
print "\n\n"
print "="*60


#####################################

# Draw a picture of the AD PMTs in ROOT
nRings = 8
nColumns = 24
nPMTs = nRings*nColumns
pmtHist = ROOT.TH3F("pmtHist",
                    "PMT Positions",
                    200,-2.5,2.5,
                    200,-2.5,2.5,
                    200,-2.5,2.5)

site = 0x01 # 0x01 (DYB), 0x02 (LA), or 0x04 (FAR)
ad =   1    # 1-2 (near sites), 1-4 (far sites)
meter = float(units.meter)
for ring in range(1,nRings+1):
    for col in range(1,nColumns+1):
        pmtid =  (site<<24) | (ad<<16) | (ring<<8) | col
        pos = pmtSvc.get(pmtid).localPosition()
        # Add pmt positions in meters
        point = (ring-1)*nColumns + (col-1)
        [x,y,z] = [pos.x()/meter, pos.y()/meter, pos.z()/meter]
        #print "Point, ring, col: ",point,ring,col,x,y,z
        pmtHist.Fill(x,y,z)
        
#Draw the figure
pmtHist.SetMarkerStyle(8)
pmtHist.SetMarkerSize(1)
pmtHist.SetMarkerColor(4)
pmtHist.GetXaxis().SetTitle("X [m]")
pmtHist.GetYaxis().SetTitle("Y [m]")
pmtHist.GetZaxis().SetTitle("Z [m]")
pmtHist.GetXaxis().SetTitleOffset(1.75)
pmtHist.GetYaxis().SetTitleOffset(1.75)
pmtHist.GetZaxis().SetTitleOffset(1)
pmtHist.Draw()

---