In This Package:

TestHist.py

Go to the documentation of this file.

#!/usr/bin/env python

import os, sys, PyCintex
Gaudi = PyCintex.makeNamespace('Gaudi')

import GaudiKernel.SystemOfUnits as units
from DybPython.Util import irange
from DybPython.DybPythonAlg import DybPythonAlg
from GaudiPython import SUCCESS, FAILURE
from GaudiPython import gbl
from GaudiPython import loaddict
loaddict('libHepMCRflx')

# Make shortcuts to ROOT classes
TH1F = gbl.TH1F
TH2F = gbl.TH2F

class GenMuonHists(DybPythonAlg):
    def __init__(self,name = "GenMuonHists", volume='/dd/Structure/AD/db-ade1'):
        DybPythonAlg.__init__(self,name)
        self.lastTime = 0
        self.volume = volume
        return

    def initialize(self):
        status = DybPythonAlg.initialize(self)
        if status.isFailure(): return status

        self.de = self.getDet(self.volume)
        
        self.stats.defaultPath = "/file1/muons"
        self.stats["dT"] = self.dT = TH1F("dT","Time between events (microsec)",
                                          10000,0,10);
        self.stats["xy"] = self.xy = TH2F("xy","Vertex, Y vs. X (m)",
                                          200,-100,100,200,-100,100)
        self.stats["yz"] = self.yz = TH2F("yz","Vertex, Z vs. Y (m)",
                                          200,-100,100,200,-100,100)
        self.stats["zx"] = self.zx = TH2F("zx","Vertex, X vs. Z (m)",
                                          200,-100,100,200,-100,100)
        return SUCCESS

    def execute(self):
        evt = self.evtSvc()
        genHdr = evt["/Event/Gen/GenHeader"]
        genEvt= genHdr.event()

        for vertex in irange(genEvt.vertices_begin(), genEvt.vertices_end()):
            xGen = vertex.position().x()/units.m
            yGen = vertex.position().y()/units.m
            zGen = vertex.position().z()/units.m
            tGen = vertex.position().t()

            point = Gaudi.XYZPoint(xGen,yGen,zGen)
            gpoint = self.de.geometry().toGlobal(point)

            self.info("Vtx: %f %f %f %f"%(xGen,yGen,zGen,tGen))

            self.dT.Fill(tGen-self.lastTime)
            self.lastTime = tGen
            self.xy.Fill(gpoint.x(),gpoint.y())
            self.yz.Fill(gpoint.y(),gpoint.z())
            self.zx.Fill(gpoing.z(),gpoint.x())
 
            # for particle in irange(vertex.particles_out_const_begin(),
            #                        vertex.particles_out_const_end()):
            #     GenKineticEnergy = (particle.momentum().e()
            #                         - particle.momentum().m())
            #     if particle.pdg_id()==11 or particle.pdg_id()==-11:
            #         self.BetaKE[nbeta] = GenKineticEnergy 
            #         nbeta += 1
            #         pass
            #     continue
            continue
        return SUCCESS

def configure(argv=[]):
    '''
    Configure GtMuoneratorTool based kinematics and an algorithm to
    histogram some quantities.
    '''
    from optparse import OptionParser
    op = OptionParser(usage=configure.__doc__)
    op.add_option("-r","--rotate-site",default=False,action="store_true",
                  help="Rotate muons to match site")
    op.add_option("-s","--site",default='SAB',
                  help="Detector site (DYB, LA, Mid, Far,SAB)")
    op.add_option("-M","--muon-file",default="",
                  help="Path to text file of muon energy and directions to sample")
    op.add_option("-R","--ratio-file",default="",
                  help="Path to ROOT file of mu+/mu- ratio histogram")
    op.add_option("-v","--volume",default="rock",
                  help="Volume to use to place vertices (rock, RPC, ADE)")
    (opts,args) = op.parse_args(args=argv)

    from StatisticsSvc.StatisticsSvcConf import StatisticsSvc
    statsSvc = StatisticsSvc()
    statsSvc.Output ={"file1":"MuonTestHist.root"}

    volume = '/dd/Structure/AD/db-ade1'
    name = 'Mu'

    from GaudiKernel import SystemOfUnits as units

    from GenMuon.GenMuonConf import GtMuoneratorTool
    muon = GtMuoneratorTool(name+'onerator')
    muon.Rotation = opts.rotate_site
    muon.WhichSite = opts.site
    muon.MuonFile = opts.muon_file
    muon.RatioFile = opts.ratio_file
    muon.Volume = opts.volume

    from GenTools.GenToolsConf import GtPositionerTool, GtTransformTool, GtTimeratorTool, GtGenerator



    # Set up positioner
    pos = GtPositionerTool(name+'Pos',Volume=volume)
    pos.Mode = "Relative"
    pos.Position = [0,0,0]
        
    # Set up timerator. note high rate of muons
    tim = GtTimeratorTool(name+'Tim')
    muonRate = 1000.   # in Hz
    lifet = 1./float(muonRate)
    tim.LifeTime = (lifet)*units.s
    print "Muon: Rate per second is ",muonRate
        
    # transform
    tra = GtTransformTool(name+'Tra',Volume=volume)
    tra.Offset = [0., 0., (0.042)*units.meter]

    gtgen = GtGenerator("muonerator")
    gtgen.GenName = 'muonerator'
    gtgen.GenTools = [muon,pos,tim,tra]

    from Gaudi.Configuration import ApplicationMgr
    theApp = ApplicationMgr()
    theApp.TopAlg.append(gtgen)

    return

def run(app):
    app.ExtSvc += ["StaticCableSvc", "StatisticsSvc"]
    alg = GenMuonHists()
    app.addAlgorithm(alg)
    pass

---