In This Package:

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

print "\tBuilding detector\n"
import xmldetdesc
xmldetdesc.config()

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

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

volume = "/dd/Structure/AD/far-lso2"

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

# Set up positioner
poser = app.property("ToolSvc.GtPositionerTool")
#poser.OutputLevel = 1
poser.Strategy = "FullVolume"
poser.Volume = volume
poser.Mode = "Uniform"
#poser.Mode = "Smeared"
poser.Spread = 2.6*units.meter
poser.Position = [0,0,0*units.meter]

# Set up gun
gun = app.property("ToolSvc.GtGunGenTool")
#gun.OutputLevel = 1
gun.ParticlesPerEvent = 1
gun.ParticleName = "n0"
#gun.Momentum = 6*units.eV   # try to enhance rayleigh scattering
gun.Momentum = 5*units.MeV
gun.MomentumMode = "Uniform"
gun.MomentumSpread = 4*units.MeV

gun.DirectionMode = "Uniform"
#from math import sin, cos, pi
#pmt_column_number = 9
#angle = (2*pmt_column_number - 1)*pi/24.0;
#gun.Direction = [ cos(angle),sin(angle),0 ] # aim for PMT 
gun.Direction = [ 1, 0, 0 ] # aim for a PMT
gun.DirectionSpread = 3
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", "GtHepMCDumper/GenDump" ]


gen = app.algorithm("GenAlg")
#gen.OutputLevel = 1
gen.GenTools = [ "GtGunGenTool", "GtPositionerTool", "GtTimeratorTool", "GtTransformTool" ]
gen.GenName = "Bang Bang"
#gen.Location = "/Event/Gen/GenHeader" # this is default anyways

#print " GtDumper"
gendump = app.algorithm("GenDump")
#gendump.Location = "/Event/Gen/GenHeader"  # this is default anyways.

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 = 1
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.OutputLevel = 3
stepseq.Members = ["HistorianStepAction","UnObserverStepAction"]


TH2DE="TH2DE"       # TH2DE is the fastest
#TH2DE="TouchableToDetectorElementFast"
historian = app.property("GiGa.GiGaStepActionSequence.HistorianStepAction")
historian.TouchableToDetelem = TH2DE

params = {
    'start' :"(start > 0)",
    'track1':"(id==1 and ProcessType==1)",
    'track2':"(id==2 and ProcessType==1)",
    'GD':    "MaterialName == '/dd/Materials/GdDopedLS'",
    'LS':    "MaterialName == '/dd/Materials/LiquidScintillator'",
    'oil':   "MaterialName == '/dd/Materials/MineralOil'",
    'iAV':   "MaterialName == '/dd/Materials/Acrylic'",
    'oAV':   "MaterialName == '/dd/Materials/Acrylic'"
    }

Historian = app.property("GiGa.GiGaStepActionSequence.HistorianStepAction");
Historian.TrackSelection = "Process=='HadronCapture'"
#Historian.VertexSelection = "Process=='HadronCapture'"
                        

unobs = app.property("GiGa.GiGaStepActionSequence.UnObserverStepAction")
unobs.TouchableToDetelem = TH2DE
unobs.Stats=[
   # for track 1
    ["pdgId_Trk1","pdg","%(track1)s and %(start)s"%params],
    ["t_Trk1",    "t" , "%(track1)s and %(start)s"%params],
    ["x_Trk1",    "x", "%(track1)s and %(start)s"%params],
    ["y_Trk1",    "y", "%(track1)s and %(start)s"%params],
    ["z_Trk1",    "z", "%(track1)s and %(start)s"%params],
    ["e_Trk1",    "E",  "%(track1)s and %(start)s"%params],
    ["p_Trk1",    "p",  "%(track1)s and %(start)s"%params],
    ["ke_Trk1",   "KE", "%(track1)s and %(start)s"%params],
    ["vx_Trk1",   "lvx","%(track1)s and %(start)s"%params],
    ["vy_Trk1",   "lvy","%(track1)s and %(start)s"%params],
    ["vz_Trk1",   "lvz","%(track1)s and %(start)s"%params],
    ["TrkLength_GD_Trk1",  "dx","%(track1)s and %(GD)s"%params],
    ["TrkLength_iAV_Trk1", "dx","%(track1)s and %(iAV)s"%params],    # not fully implemented
    ["TrkLength_LS_Trk1",  "dx","%(track1)s and %(LS)s"%params],
    ["TrkLength_oAV_Trk1", "dx","%(track1)s and %(oAV)s"%params],    # not fully implemented
    ["TrkLength_Oil_Trk1", "dx","%(track1)s and %(oil)s"%params],
   # for track 2 now
    ["pdgId_Trk2","pdg","%(track2)s and %(start)s"%params],
    ["t_Trk2",    "t" , "%(track2)s and %(start)s"%params],
    ["x_Trk2",    "x", "%(track2)s and %(start)s"%params],
    ["y_Trk2",    "y", "%(track2)s and %(start)s"%params],
    ["z_Trk2",    "z", "%(track2)s and %(start)s"%params],
    ["e_Trk2",    "E",  "%(track2)s and %(start)s"%params],
    ["p_Trk2",    "p",  "%(track2)s and %(start)s"%params],
    ["ke_Trk2",   "KE", "%(track2)s and %(start)s"%params],
    ["vx_Trk2",   "lvx","%(track2)s and %(start)s"%params],
    ["vy_Trk2",   "lvy","%(track2)s and %(start)s"%params],
    ["vz_Trk2",   "lvz","%(track2)s and %(start)s"%params],
    ["TrkLength_GD_Trk2",  "dx","%(track2)s and %(GD)s"%params],
    ["TrkLength_iAV_Trk2", "dx","%(track2)s and %(iAV)s"%params],    # not fully implemented
    ["TrkLength_LS_Trk2",  "dx","%(track2)s and %(LS)s"%params],
    ["TrkLength_oAV_Trk2", "dx","%(track2)s and %(oAV)s"%params],    # not fully implemented
    ["TrkLength_Oil_Trk2", "dx","%(track2)s and %(oil)s"%params]
    ]
                                                

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


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/DayaBay", ]
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" ]
pull = app.algorithm("PullEvent")
pull.OutputLevel = 1 
push = app.algorithm("PushKine")
push.Converter = "HepMCtoG4"
#push.Location = "/Event/Gen/GenHeader" # default anyway

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

# Add the GenHist algorithm
app.TopAlg += ["DetSimVali"]
vali= app.algorithm("DetSimVali")
vali.OutPutLevel=3

histsvc = app.service("THistSvc")
histsvc.Output =["file1 DATAFILE='../mc_ntuple/neutron.root' OPT='RECREATE' TYP='ROOT' "]
    
app.run(10)
    

---