In This Package:

pmtHits.py

Go to the documentation of this file.

#############################################################################
###
### Python script for visualizing PMT hits on AD detectors
###
### Kim Boddy, kboddy@caltech.edu
### 31 Jul 2008 Created initial code
### 01 Aug 2008 Implemented code to use Gaudi Services, from: 
###    NuWa-trunk/dybgaudi/Detector/DetHelpers/python/DetHelpers/pmtCoords.py
###
#############################################################################

# 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.

import ROOT

### Input settings
# set output: contour plot (0), exact hit location (1), or exact 2D (2)
output_st = raw_input("Choose to plot as contour (default), exact, exact2d: ")
if    output_st=='contour': output=0
elif  output_st=='exact':   output=1
elif  output_st=='exact2d': output=2
else: output=0
# set input file and associated tree
file_st = raw_input("ROOT file to read (default:electron.root): ")
tree_st = raw_input("Tree file to read (default:DetSimValiTree): ")
if not file_st: file_st = 'electron.root'
if not tree_st: tree_st = 'DetSimValiTree'
f = ROOT.TFile(file_st,'READ') # ROOT will return error if file_st does not exist
tree = f.Get(tree_st)          # ROOT will return error if tree_st does not exist

### Indexing Conventions
# site:  DayaBay=0x01, LingAo=0x02, Far=0x04
# DetID: AD1=1, AD2=2, AD3=3, AD4=4
# pmtID consists of 32 bits, 4 pieces of 8-bit information strung together
#   site, detector ID, ring, column 
#   (i.e., site info is leftmost 8 bits and column info is rightmost 8 bits)

### Create histograms and graphs for 8 detectors
hist = [ROOT.TH2F('hist0','PMT Hits: DB site, AD1',24,1,25,8,1,9),
        ROOT.TH2F('hist1','PMT Hits: DB site, AD2',24,1,25,8,1,9),
        ROOT.TH2F('hist2','PMT Hits: LA site, AD1',24,1,25,8,1,9),
        ROOT.TH2F('hist3','PMT Hits: LA site, AD2',24,1,25,8,1,9),
        ROOT.TH2F('hist4','PMT Hits: far site, AD1',24,1,25,8,1,9),
        ROOT.TH2F('hist5','PMT Hits: far site, AD2',24,1,25,8,1,9),
        ROOT.TH2F('hist6','PMT Hits: far site, AD3',24,1,25,8,1,9),
        ROOT.TH2F('hist7','PMT Hits: far site, AD4',24,1,25,8,1,9)]

if output==1:
    gr = [ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph(),
          ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph(),ROOT.TGraph()]

if output==2:
    gr = [ROOT.TGraph2D(),ROOT.TGraph2D(),ROOT.TGraph2D(),ROOT.TGraph2D(),
          ROOT.TGraph2D(),ROOT.TGraph2D(),ROOT.TGraph2D(),ROOT.TGraph2D()]

if (output==1 or output==2):
    gr[0].SetTitle('PMT Hits: DB site, AD1')
    gr[1].SetTitle('PMT Hits: DB site, AD2')
    gr[2].SetTitle('PMT Hits: LA site, AD1')
    gr[3].SetTitle('PMT Hits: LA site, AD2')
    gr[4].SetTitle('PMT Hits: far site, AD1')
    gr[5].SetTitle('PMT Hits: far site, AD2')
    gr[6].SetTitle('PMT Hits: far site, AD3')
    gr[7].SetTitle('PMT Hits: far site, AD4')

c = []

### Set Units
meter  = float(units.meter)
degree = float(units.degree)

### Read tree
entries = tree.GetEntries()

for j in xrange(entries):
    i = tree.LoadTree(j)
    if i < 0:
        break

    nb = tree.GetEntry(j)
    if nb <=0:
        continue

    SimHits = tree.SimHitsEntries
    pmtID   = tree.sensID
    hitx    = tree.hitx
    hity    = tree.hity
    hitz    = tree.hitz
    if SimHits<0:
        continue

    for k in xrange(SimHits):
        site = pmtID[k]>>24
        det  = ((site<<24)^pmtID[k])>>16
        ring = (((site<<24)|(det<<16))^pmtID[k])>>8
        col  = (((site<<24)|(det<<16)|(ring<<8))^pmtID[k])

        if   (site==1) and (det==1): idx=0
        elif (site==1) and (det==2): idx=1
        elif (site==2) and (det==1): idx=2
        elif (site==2) and (det==2): idx=3
        elif (site==4) and (det==1): idx=4
        elif (site==4) and (det==2): idx=5
        elif (site==4) and (det==3): idx=6
        elif (site==4) and (det==4): idx=7
        else:
            print 'Unexpected site/AdID number - skipping this hit'
            continue
        
        # Get detector information
        pmtGeom     = pmtSvc.get(pmtID[k])      # get PMT geometry
        adDetector  = pmtGeom.parentDetector()  # get AD for this PMT
        pmtDetector = pmtGeom.detectorElement() # get PMT detector
        
        # Find PMT hit points 
        pmtHit = gbl.Gaudi.XYZPoint(hitx[k],hity[k],hitz[k]) # hit on PMT (x,y,z)
        gpHit  = pmtDetector.geometry().toGlobal(pmtHit)     # global hit position
        lpHit  = adDetector.geometry().toLocal(gpHit)        # local hit in AD coords
            
        # Convert (x,y,z) coords to (rho,z,phi) coords 
        # Is there a Gaudi Service for this?
        x = lpHit.x() / meter
        y = lpHit.y() / meter
        z = lpHit.z() / meter
        rho = ROOT.TMath.Sqrt(x*x + y*y) 
        phi = ROOT.TMath.ATan2(y,x) / degree

        hist[idx].Fill(col,ring)
        NPoint = int(hist[idx].GetEntries())-1 #since SetPoint starts at point 0
        if output==1: gr[idx].SetPoint(NPoint,phi,z)
        if output==2: gr[idx].SetPoint(NPoint,x,y,z)

ROOT.gStyle.SetPalette(1)
ROOT.gStyle.SetOptStat(11)

### Draw only non-empty histograms/graphs
for i in xrange(8):
    if hist[i].GetEntries()==0: continue
    c += [ROOT.TCanvas()]   

    if output==0:   
        hist[i].GetXaxis().SetTitle('Column')
        hist[i].GetYaxis().SetTitle('Ring')
        hist[i].Draw('COLZ')
    elif output==1: 
        gr[i].GetXaxis().SetTitle('#phi [deg]')
        gr[i].GetYaxis().SetTitle('Z [m]')
        gr[i].SetMarkerStyle(1)
        gr[i].SetMarkerColor(2)
        gr[i].Draw('AP')
    elif output==2:
        gr[i].GetXaxis().SetTitle('X [m]')
        gr[i].GetYaxis().SetTitle('Y [m]')
        gr[i].GetZaxis().SetTitle('Z [m]')
        gr[i].GetXaxis().SetTitleOffset(1.75)
        gr[i].GetYaxis().SetTitleOffset(1.75)
        gr[i].GetZaxis().SetTitleOffset(1)
        gr[i].SetMarkerStyle(1)
        gr[i].SetMarkerColor(2)
        gr[i].Draw('P')
    else: print 'Invalid draw command: problem with reading user input'

---