In This Package:

poolgun.py

Go to the documentation of this file.

#!/usr/bin/env python


particle = "e+"
#particle = "opticalphoton"
use_correct_pmt_angle = False
hit_pmt_gap = False


class MyGenToolsConfig:
    def __init__(self,volume):
        from  GaudiPython import AppMgr
        from GaudiKernel import SystemOfUnits as units
        from math import sin, cos, pi

        app = AppMgr()
        app.TopAlg = []

        # Set up timerator
        import GaudiKernel.SystemOfUnits as units
        tim = app.property("ToolSvc.GtTimeratorTool")
        tim.LifeTime = int(1*units.second)

        # Set up positioner
        poser = app.property("ToolSvc.GtPositionerTool")
        poser.Volume = volume
        poser.Mode = "Fixed"
        if hit_pmt_gap:
            angle=pi/24.0
            radius=0.0*units.meter
            poser.Position = [ radius*cos(angle),radius*sin(angle),0 ] # aim between columns 1 and 2
        else:
            poser.Position = [0,0,5.1*units.meter]

        poser.OutputLevel = 1

        # Set up gun
        gun = app.property("ToolSvc.GtGunGenTool")
        gun.OutputLevel = 1
        if particle == "opticalphoton":
            gun.ParticlesPerEvent = 10000
            gun.ParticleName = "opticalphoton"
            gun.Momentum = 2.5*units.eV
            gun.PolarizeMode = "Random"
        elif particle == "muon":
            gun.ParticlesPerEvent = 1
            gun.ParticleName = "mu-"
            gun.Momentum = 10*units.GeV
            pass
        else:
            gun.ParticlesPerEvent = 1
            gun.ParticleName = "e+"
            gun.Momentum = 3*units.MeV
            pass
        gun.MomentumMode = "Fixed"
        gun.DirectionMode = "Fixed"
        pmt_column_number = 9
        if use_correct_pmt_angle:
            if hit_pmt_gap:
                angle = (2*pmt_column_number)*pi/24.0;
            else:
                angle = (2*pmt_column_number - 1)*pi/24.0;
            gun.Direction = [ cos(angle),sin(angle),0 ] # aim for PMT 
        else:
            if hit_pmt_gap:
                angle=pi/24.0
                gun.Direction = [ cos(angle),sin(angle),0 ] # aim between columns 1 and 2
            else:
                gun.Direction = [ 0,0,-1 ]
        print 'gun.Direction=',gun.Direction

        trans = app.property("ToolSvc.GtTransformTool")
        trans.Volume = volume

        app.TopAlg += [ "GaudiSequencer/GenSeq" ]
        genseq = app.algorithm("GenSeq")
        genseq.Members = [ "GtGenerator/GenAlg" ]

        gen = app.algorithm("GenAlg")
        gen.OutputLevel = 1
        gen.GenTools = [ "GtGunGenTool", "GtPositionerTool", "GtTimeratorTool", "GtTransformTool" ]
        gen.GenName = "Bang Bang"

        if particle != "opticalphoton":
            print " GtDumper"
            genseq.Members += [ "GtHepMCDumper/GenDump" ]
            gendump = app.algorithm("GenDump")
            gendump.OutputLevel = 1

        app.ExtSvc += ["GiGa"]

        modularPL = app.property("GiGa.GiGaPhysListModular")
        modularPL.OutputLevel = 1
        modularPL.CutForElectron = 100*units.micrometer
        modularPL.CutForPositron = 100*units.micrometer
        modularPL.CutForGamma = 1*units.millimeter
        modularPL.PhysicsConstructors = [ 
            "DsPhysConsGeneral", 
            "DsPhysConsOptical",
            "DsPhysConsEM",
#            "DsPhysConsElectroNu",
#            "DsPhysConsHadron",
#            "DsPhysConsIon"
            ]
        optical = app.property("GiGa.GiGaPhysListModular.DsPhysConsOptical")
        #optical.UseCerenkov = False
        #optical.UseScintillation = False

        giga = app.service("GiGa")
        giga.OutputLevel = 1
        giga.PhysicsList = "GiGaPhysListModular"

        gggeo = app.service("GiGaGeo")
        gggeo.OutputLevel = 1
        gggeo.XsizeOfWorldVolume = 2.4*units.kilometer
        gggeo.YsizeOfWorldVolume = 2.4*units.kilometer
        gggeo.ZsizeOfWorldVolume = 2.4*units.kilometer


        giga.SteppingAction = "GiGaStepActionSequence"
        stepseq = app.property("GiGa.GiGaStepActionSequence")
        stepseq.Members = ["HistorianStepAction","UnObserverStepAction"]

        historian = app.property("GiGa.GiGaStepActionSequence.HistorianStepAction")

        params = {
            'muonid':"(pdg==13 or pdg==-13)",
            'ows':"(MaterialName == '/dd/Materials/Water' and "+\
                "Volume == '/dd/Structure/Pool/far-ows')",
            'iws':"(MaterialName == '/dd/Materials/Water' and "+\
                "Volume == '/dd/Structure/Pool/far-iws')",
            'lso':"MaterialName == '/dd/Materials/LiquidScintillator'",
            'gds':"MaterialName == '/dd/Materials/GdDopedLS'"
            }
        unobs = app.property("GiGa.GiGaStepActionSequence.UnObserverStepAction")
        unobs.Stats=[
            ["MuonTrkLengthInOws","dx","%(muonid)s and %(ows)s"%params],
            ["MuonTrkLengthInIws", "dx", "%(muonid)s and %(iws)s"%params],
            ["MuonTrkLengthInLS", "dx", "%(muonid)s and %(lso)s"%params],
            ["MuonTrkLengthInGdLS","dx", "%(muonid)s and %(gds)s"%params]
            ]



        giga.StackingAction="DsOpStackAction"
        
        eds=app.service("EventDataSvc")
        eds.ForceLeaves=True



        # Make Geant4 sing!
        ggrm = app.property("GiGa.GiGaMgr")
        ggrm.Verbosity = 0
        event_ac_cmds = app.property("GiGa.GiGaEventActionCommand")
        event_ac_cmds.BeginOfEventCommands = [
            "/control/verbose 0",
            "/run/verbose 0",
            "/event/verbose 0",
            "/tracking/verbose 0",
            "/geometry/navigator/verbose 0"
            ]
        giga.EventAction = "GiGaEventActionCommand"


        app.TopAlg += [ "GaudiSequencer/SimSeq" ]
        simseq = app.algorithm("SimSeq")
        simseq.Members = [ "GiGaInputStream/GGInStream" ]

        ggin = app.algorithm("GGInStream")
        ggin.OutputLevel = 1
        ggin.ExecuteOnce = True
        ggin.ConversionSvcName = "GiGaGeo"
        ggin.DataProviderSvcName = "DetectorDataSvc"
        ggin.StreamItems = [ "/dd/Structure/Sites/far-rock",
                             "/dd/Geometry/AdDetails/AdSurfacesAll",
                             "/dd/Geometry/AdDetails/AdSurfacesFar",
                             "/dd/Geometry/PoolDetails/FarPoolSurfaces",
                             "/dd/Geometry/PoolDetails/PoolSurfacesAll",
                             ]

        simseq.Members += [ "DsPushKine/PushKine", "DsPullEvent/PullEvent" ]
        push = app.algorithm("PushKine")
        push.Converter = "HepMCtoG4"

        pull = app.algorithm("PullEvent")
        pull.OutputLevel = 1

        # Class name to use is set in DetDesc xml's "sensdet" attribute.
        pmtsd = app.property("GiGaGeo.DsPmtSensDet")
        pmtsd.OutputLevel = 2

        return

if '__main__' == __name__:
    import os

    print "\tConfiguring geometry\n"
    import xmldetdesc
    xmldetdesc.config()

    from GaudiPython import AppMgr
    app = AppMgr()
    app.EvtSel = "NONE"

    #xmldetdesc.dump()

    gtc = MyGenToolsConfig(volume="/dd/Structure/Pool/far-ows")

    app.run(1)
    

---