In This Package:

U-Th-K.py

Go to the documentation of this file.

#!/usr/bin/env python
from  GaudiPython import AppMgr
from GaudiKernel import SystemOfUnits as units

import os,sys

import xmldetdesc
xmldetdesc.config()

import GaudiPython as gp

iapp = gp.iService("ApplicationMgr")
iapp.SvcMapping = [
    'EvtDataSvc/EventDataArchiveSvc',
    'DybDataSvc/EventDataSvc',
    #"EvtDataSvc/EventDataSvc",
    "DetDataSvc/DetectorDataSvc",
    "HistogramSvc/HistogramDataSvc",
    "HbookCnv::PersSvc/HbookHistSvc",
    "RootHistCnv::PersSvc/RootHistSvc",
    "EvtPersistencySvc/EventPersistencySvc",
    "DetPersistencySvc/DetectorPersistencySvc",
    "HistogramPersistencySvc/HistogramPersistencySvc",
    ]
        
app = gp.AppMgr()

app.EvtSel = "NONE"

trimSvc = app.service('ArchiveTrimSvc')
trimSvc.DefaultWindowSeconds = 0.0001 # not in system-of-units time!

app.TopAlg = ["GiGaInputStream/GGInStream",
              'Sim15',
              'DybTrimIO'
              ]

StageKin=app.property("ToolSvc.Kinematic")
StageKin.Sequencer='GnrtrSeq'
GenSeq=app.property('GnrtrSeq')
GenSeq.Members = [
#    "Gnrtr/GenK40",
    "Gnrtr/GenU",
#    "Gnrtr/GenTh",
    ]

StageDet=app.property("ToolSvc.Detector")
StageDet.Sequencer='DetSimSeq'
DSSeq=app.property('DetSimSeq')
DSSeq.Members = ['DetSimProc']

StageElec=app.property("ToolSvc.Electronic")
StageElec.Sequencer='ElecSimSeq'
ESSeq=app.property('ElecSimSeq')
ESSeq.Members = ['ElecSimProc']

StageTrig=app.property("ToolSvc.Trigger")
StageTrig.Sequencer='TrigSimSeq'
TSSeq=app.property('TrigSimSeq')
TSSeq.Members = ['TrigReadProc']

StageSingle=app.property("ToolSvc.SingleLoader")
StageSingle.Sequencer='SingleSeq'
SingleSeq=app.property('SingleSeq')
SingleSeq.Members = ['SingleLoader']


# Kinematic
genk40=app.algorithm('GenK40')
genk40.ThisStageName='Kinematic'
genk40.LowerStageName=''

genu=app.algorithm('GenU')
genu.ThisStageName='Kinematic'
genu.LowerStageName=''

genth=app.algorithm('GenTh')
genth.ThisStageName='Kinematic'
genth.LowerStageName=''

# Detector
sim=app.algorithm('DetSimProc')
sim.ThisStageName='Detector'
sim.LowerStageName='Kinematic'

# Electronic
elec=app.algorithm('ElecSimProc')
elec.ThisStageName='Electronic'
elec.LowerStageName='Detector'

# Trigger
trig=app.algorithm('TrigReadProc')
trig.ThisStageName='Trigger'
trig.LowerStageName='Electronic'

# SingleLoader
singll=app.algorithm('SingleLoader')
singll.ThisStageName='SingleLoader'
singll.LowerStageName='Trigger'

# Configure generator
# K40
genk40.GenTools = [ "GtHepEvtGenTool/K40Gen",
                    "GtPositionerTool/K40Pos",
                    "GtTimeratorTool/K40Tim",
                    "GtTransformTool/K40Tra"]
genk40.OutputLevel = 2
genk40.GenName = "K40_ticktack"

volume =  "/dd/Structure/AD/far-lso2"
#genk40.Location = "/Event/Gen/HepMCEvents" # this is default anyways
#
# Set up K40Gen
seed=1234567
nevents=1000
hepevt_source="K40.exe -seed %(seed)s -n %(nevents)s|"
# Deal with executable
if hepevt_source[-1] == "|":
    exe = hepevt_source.split(' ')[0]
    if exe[0] != '/':   # Try to find full path
        import os, os.path
        path = os.getenv('PATH')
        for p in path:
            if (os.path.isfile(path+"/"+exe)):
                exe = path+"/"+exe
                break
            continue
        pass
    source = exe + ' ' + ' '.join(hepevt_source.split(' ')[1:])
    if "%" in source:   # Fill in any placemarks
        source = source%{'nevents':str(nevents),
                         'seed':str(seed)}
        pass
    pass


#generator
k40gen=app.property('ToolSvc.K40Gen')
k40gen.OutputLevel = 2
k40gen.HepEvtDataSource = source

# Set up positioner
k40pos=app.property('ToolSvc.K40Pos')
k40pos.OutputLevel = 2
k40pos.Volume = volume
k40pos.Strategy = "FullVolume" # also "AvoidDaughters" and "Surface"
#k40pos.Mode = "Uniform"
k40pos.Mode = "Fixed"
k40pos.Spread = 4*units.m
k40pos.Position = [0,0,2.5*units.m]

# Set up timerator
app.property('ToolSvc.K40Tim').LifeTime = int(30000*units.ns)

# transform
k40tra = app.property("ToolSvc.K40Tra")
k40tra.Volume = volume


#setup U decay
genu.GenTools = [ "GtUDecaysGenTool/UGen",
                  "GtPositionerTool/UPos",
                  "GtUDecaysTimeTool/UTim",
                  "GtTransformTool/UTra" ]
genu.OutputLevel = 2
genu.GenName = "U_ticktack"

# U decay generator
ugen=app.property("ToolSvc.UGen")
ugen.OutputLevel =1
ugen.RandomSeed =450

spray_photons = True

# Set up positioner
upos = app.property("ToolSvc.UPos")
upos.Volume = volume
upos.Mode = "Fixed"
if spray_photons:
    upos.Position = [0,0,2500*units.mm]
else:
    upos.Position = [-3*units.meter,3*units.meter,0]
upos.OutputLevel = 2
                
# Set up timerator
utim = app.property("ToolSvc.UTim")
utim.LifeTime = int(30000*units.ns)
utim.OutputLevel = 2

# set up transform
utrans = app.property("ToolSvc.UTra")
utrans.Volume = volume
utrans.OutputLevel = 2

# set up thorium decay
genth.GenTools = [ "GtThDecaysGenTool/ThGen",
                   "GtPositionerTool/ThPos",
                   "GtTimeratorTool/ThTim",   # ??? no GtThDecaysTimeTool
                   "GtTransformTool/ThTra" ]
genth.OutputLevel = 2
genth.GenName = "Th_ticktack"

# Th decay generator
thgen=app.property("ToolSvc.ThGen")
thgen.OutputLevel =2
thgen.RandomSeed =100

# Set up positioner
thpos = app.property("ToolSvc.ThPos")
thpos.Volume = volume
thpos.Mode = "Fixed"
if spray_photons:
    upos.Position = [0,0,2500*units.mm]
else:
    upos.Position = [-3*units.meter,3*units.meter,0]
thpos.OutputLevel = 2

# Set up timerator
thtim = app.property("ToolSvc.ThTim")
thtim.LifeTime = int(30000*units.ns)
thtim.OutputLevel = 2

# set up transform
thtrans = app.property("ToolSvc.ThTra")
thtrans.Volume = volume
thtrans.OutputLevel = 2

#-------------------------------------------------------------------------
# Configure geant4 simulation
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"
            ]
# in DsPhysConsOptical.cc these properties are declared.
optical = app.property("GiGa.GiGaPhysListModular.DsPhysConsOptical")
#
#optical.UseScintillation = False  # false option will kill all optical photons, should be true all the time
#optical.UseRayleigh = True       # enable rayleigh scattering
optical.CerenPhotonScaleWeight=4.0  # scale photon yield down and scale QE up
optical.ScintPhotonScaleWeight=4.0

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

gggeo = app.service("GiGaGeo")
gggeo.OutputLevel = 3
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"]

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

event_ac_cmds.BeginOfEventCommands = silent_cmds
giga.EventAction = "GiGaEventActionCommand"


ggin = app.algorithm("GGInStream")

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

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

sim.Converter = "HepMCtoG4"
#push.Location = "/Event/Gen/GenHeader" # default anyway
sim.OutputLevel=2

# other properties for ElecSim
#elec.PmtTool = 'EsPmtEffectPulseTool'
elec.PmtTool = 'EsIdealPulseTool'
elec.RpcTool = 'EsIdealPulseTool'
elec.FeeTool = 'EsIdealFeeTool'
elec.FecTool = 'EsIdealFecTool'
elec.OutputLevel=2

# trigger tools
trig.TrigTools=['TsMultTriggerTool']  # default
trig.RoTools=['ROsFeeReadoutTool','ROsFecReadoutTool']
trig.OutputLevel=2
nhits=app.property("ToolSvc.TsMultTriggerTool")
nhits.OutputLevel=2
fee=app.property("ToolSvc.ROsFeeReadoutTool")
fec=app.property("ToolSvc.ROsFecReadoutTool")
fee.OutputLevel=2
fec.OutputLevel=2

# Single Loader
singll.OutputLevel=2

# others
eds=app.service("EventDataSvc")
eds.ForceLeaves=True

#
app.EvtMax = 150
app.EvtSel ="NONE"

#
sim15=app.algorithm('Sim15')
#sim15.TopStage='Kinematic'
#sim15.TopStage='Detector'
sim15.TopStage='Electronic'
#sim15.TopStage='Trigger'
#sim15.TopStage='SingleLoader'
sim15.OutputLevel=2

histsvc = app.service("THistSvc")

msv = app.service('MessageSvc')
msv.OutputLevel = 2
app.initialize()
app.run( app.EvtMax )
app.exit()

---