In This Package:

IBDn-sim-rootio.py

Go to the documentation of this file.

#!/usr/bin/env python

import sys, os
import GaudiPython as gp
from GaudiKernel import SystemOfUnits as units

try:
    io = sys.argv[1]
except IndexError:
    io = "output"

input = False
if "input" == io: 
    input = True
output = not input
use_aes = True

gen_location = "/Event/Gen/GenHeader"
sim_location = "/Event/Sim/SimHeader"

if output:
    seed = sys.argv[2]
    nEvents = sys.argv[3]
    GaudiIOfile = "AD-neutron-Gaudi-"+nEvents+"evts-ws"+seed+".root"
    RootOutputfile = "AD-neutron-Root-"+nEvents+"evts-ws"+seed+".root"
    print "\tSeed of Neutron is:", seed
    print "\tNumber of events is:", nEvents
    print "\tGaudi output file is:", GaudiIOfile
    print "\tROOT output file is:", RootOutputfile

if input:
    GaudiIOfile = sys.argv[2]
    print "\tGaudi input file is:", GaudiIOfile
    nEvents = sys.argv[3]
    print "\tTake ", nEvents, " from:", GaudiIOfile
    RootOutputfile = sys.argv[4]
    print "\tROOT output file is:", RootOutputfile

if use_aes:
    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",
        ]

if not input:
    import xmldetdesc
    print "\tBuilding detector.\n"
    xmldetdesc.config()

app = gp.AppMgr(outputlevel=4)

msg = app.service("MessageSvc")
#msg.Format = "\--%-F%40W%S%7W%R%T %23W%t\n \-> %0W%M"
msg.Format = "% F%25W%S%7W%R%T %0W%M"
msg.useColors = False
msg.fatalColorCode=['red','white']
msg.errorColorCode=['red']
msg.warningColorCode=['yellow']
msg.debugColorCode=['blue']
msg.verboseColorCode=['cyan']

app.EvtMax = int(nEvents)

if input:
    app.EvtSel =""              # yeah, I know, it's weird.
else:
    app.EvtSel = "NONE"

if input:
    #app.ExtSvc += [ "RootIOEvtSelector/EventSelector" ]
    app.ExtSvc += [ "DybEvtSelector/EventSelector" ]
if output:
    app.ExtSvc += [ "DybStorageSvc" ]
    dss = app.service("DybStorageSvc")
    dss.OutputLevel = 4
    
app.ExtSvc += [ "RootIOCnvSvc" ]

per = app.service("EventPersistencySvc")
per.CnvServices = [ "RootIOCnvSvc" ]

eds = app.service("EventDataService")
eds.OutputLevel = 4

rio = app.property("RootIOCnvSvc")
rio.OutputLevel = 4

iomap = { "default": GaudiIOfile }

if input:
    rio.InputStreams = iomap
if output:
    rio.OutputStreams = iomap

app.TopAlg = [ ]

if not input:
    
    volume = "/dd/Structure/Sites/db-rock/db-ows/db-curtain/db-iws/db-ade1/db-sst1/db-oil1/db-oav1/db-lso1"
    
    # Set up timerator
    tim = app.property("ToolSvc.GtTimeratorTool")
    tim.LifeTime = int(1*units.second)

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

    # Set up gun (Or a generator)
    # gun = app.property("ToolSvc.GtGunGenTool")
    gun = app.property("ToolSvc.GtHepEvtGenTool")
    # gun.Volume = volume
    gun.OutputLevel = 4
    import os
    path=os.getenv("PATH")
    for dir in path.split(":"):
        if os.path.exists(dir + "/InverseBeta.exe"):
            gun.HepEvtDataSource = dir + "/InverseBeta.exe -seed " + seed +" -neutron_only -n " + nEvents + " |"
            break
        continue
    # 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" ]

    # set up gentools
    gen = app.algorithm("GenAlg")
    gen.OutputLevel = 4
    gen.GenTools = [ "GtHepEvtGenTool", "GtPositionerTool", 
                     "GtTimeratorTool", "GtTransformTool" ]
    gen.GenName = "IBD neutron"
    # 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 = 4
    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
    
    gggeo = app.service("GiGaGeo")
    gggeo.OutputLevel = 4
    gggeo.XsizeOfWorldVolume = 2.4*units.kilometer
    gggeo.YsizeOfWorldVolume = 2.4*units.kilometer
    gggeo.ZsizeOfWorldVolume = 2.4*units.kilometer

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

    giga.SteppingAction = "GiGaStepActionSequence"
    stepseq = app.property("GiGa.GiGaStepActionSequence")
    stepseq.OutputLevel = 4
    stepseq.Members = ["HistorianStepAction","UnObserverStepAction"]
    
    TH2DE="TH2DE"       # TH2DE is the fastest
    # TH2DE="TouchableToDetectorElementFast"
    historian = app.property("GiGa.GiGaStepActionSequence.HistorianStepAction")
    historian.TouchableToDetelem = TH2DE
    historian.TrackSelection = "(pdg == 2112)"
    # historian.TrackSelection = "any"
    historian.VertexSelection = "(pdg == 2112)"
    # historian.VertexSelection = "any"

    params = {
        'start' :"(start > 0)", # Same as 'firstvtx'
        'track1':"(id==1)",
        'track2':"(id==2)",
        'GD':    "MaterialName == '/dd/Materials/GdDopedLS'",
        'LS':    "MaterialName == '/dd/Materials/LiquidScintillator'",
        'MO':   "MaterialName == '/dd/Materials/MineralOil'",
        'IAV':   "DetectorElementName == 'db-iav1'",
        'OAV':   "DetectorElementName == 'db-oav1'",
        'IWS': "MaterialName == '/dd/Materials/IwsWater'",
        'OWS': "MaterialName == '/dd/Materials/OwsWater'",
        'lastvtx': "IsStopping == 1",
        'firstvtx': "IsStarting == 1",
        'NeutronTrk': "pdg == 2112",
        'NeutronMom': "creator == 2112",
        'Gamma': "pdg == 22",
        'Muon': "(pdg == 13 or pdg == -13)"
        }
    
    unobs = app.property("GiGa.GiGaStepActionSequence.UnObserverStepAction")
    unobs.TouchableToDetelem = TH2DE
    
    unobs.Stats=[
        
        ["MuonTrkLengthInOws", "dx", "%(Muon)s and %(OWS)s"%params],
        ["MuonTrkLengthInIws", "dx", "%(Muon)s and %(IWS)s"%params],
        ["MuonTrkLengthInLS", "dx", "%(Muon)s and %(LS)s"%params],
        ["MuonTrkLengthInGdLS","dx", "%(Muon)s and %(GD)s"%params],
        ["dEInn","dE", "(pdg!=20022) and %(IWS)s"%params],
        ["dEOut","dE", "(pdg!=20022) and %(OWS)s"%params],
        ["MuonStop", "dx",  "%(Muon)s and %(lastvtx)s"%params],
        
        ["EDepInGdLS", "dE", "%(GD)s"%params],
        ["EDepInLS", "dE", "%(LS)s"%params],
        ["EDepInIAV", "dE", "%(IAV)s"%params],
        ["EDepInOAV", "dE", "%(OAV)s"%params],
        ["EDepInOIL", "dE", "%(MO)s"%params],

        ["QEDepInGdLS", "qdE", "%(GD)s"%params],
        ["QEDepInLS", "qdE", "%(LS)s"%params],
        ["QEDepInIAV", "qdE", "%(IAV)s"%params],
        ["QEDepInOAV", "qdE", "%(OAV)s"%params],
        ["QEDepInOIL", "qdE", "%(MO)s"%params],

        ["tQESumGdLS", "qEt", "%(GD)s"%params],
        ["xQESumGdLS", "qEx", "%(GD)s"%params],
        ["yQESumGdLS", "qEy", "%(GD)s"%params],
        ["zQESumGdLS", "qEz", "%(GD)s"%params],

        ["tQESumLS", "qEt", "%(LS)s"%params],
        ["xQESumLS", "qEx", "%(LS)s"%params],
        ["yQESumLS", "qEy", "%(LS)s"%params],
        ["zQESumLS", "qEz", "%(LS)s"%params],

        ["tQESumMO", "qEt", "%(MO)s"%params],
        ["xQESumMO", "qEx", "%(MO)s"%params],
        ["yQESumMO", "qEy", "%(MO)s"%params],
        ["zQESumMO", "qEz", "%(MO)s"%params],

        ["tGen",   "t","%(NeutronTrk)s and %(firstvtx)s"%params],
        ["xGen",   "x","%(NeutronTrk)s and %(firstvtx)s"%params],
        ["yGen",   "y","%(NeutronTrk)s and %(firstvtx)s"%params],
        ["zGen",   "z","%(NeutronTrk)s and %(firstvtx)s"%params],

        ["tCap",   "t","%(NeutronTrk)s and %(lastvtx)s"%params],
        ["xCap",   "x","%(NeutronTrk)s and %(lastvtx)s"%params],
        ["yCap",   "y","%(NeutronTrk)s and %(lastvtx)s"%params],
        ["zCap",   "z","%(NeutronTrk)s and %(lastvtx)s"%params],

        ["capTarget", "capTargetZ","%(track1)s and %(lastvtx)s"%params],
        
        # 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],
        ["tEnd_Trk1",    "t" , "%(track1)s and %(lastvtx)s"%params],
        ["xEnd_Trk1",    "x", "%(track1)s and %(lastvtx)s"%params],
        ["yEnd_Trk1",    "y", "%(track1)s and %(lastvtx)s"%params],
        ["zEnd_Trk1",    "z", "%(track1)s and %(lastvtx)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],
        ["TrkLength_LS_Trk1",  "dx","%(track1)s and %(LS)s"%params],
        ["TrkLength_oAV_Trk1", "dx","%(track1)s and %(OAV)s"%params],
        ["TrkLength_Oil_Trk1", "dx","%(track1)s and %(MO)s"%params],
        # track 2
        ["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],
        ["tEnd_Trk2",    "t" , "%(track2)s and %(lastvtx)s"%params],
        ["xEnd_Trk2",    "x", "%(track2)s and %(lastvtx)s"%params],
        ["yEnd_Trk2",    "y", "%(track2)s and %(lastvtx)s"%params],
        ["zEnd_Trk2",    "z", "%(track2)s and %(lastvtx)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],
        ["TrkLength_LS_Trk2",  "dx","%(track2)s and %(LS)s"%params],
        ["TrkLength_oAV_Trk2", "dx","%(track2)s and %(OAV)s"%params],
        ["TrkLength_Oil_Trk2", "dx","%(track2)s and %(MO)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/GenSeq" ]
    app.TopAlg += [ "GaudiSequencer/SimSeq" ]
    simseq = app.algorithm("SimSeq")
    simseq.Members = [ "GiGaInputStream/GGInStream" ]
    
    ggin = app.algorithm("GGInStream")
    ggin.OutputLevel = 4
    ggin.ExecuteOnce = True
    ggin.ConversionSvcName = "GiGaGeo"
    ggin.DataProviderSvcName = "DetectorDataSvc"
    ggin.StreamItems = [ "/dd/Structure/Sites/db-rock",
                         "/dd/Geometry/AdDetails/AdSurfacesAll",
                         "/dd/Geometry/AdDetails/AdSurfacesNear"
                         ]
    
    simseq.Members = [ "GiGaInputStream/GGInStream", "DsPushKine/PushKine", 
                       "DsPullEvent/PullEvent" 
                       ]
    
    pull = app.algorithm("PullEvent")
    pull.OutputLevel = 4 
    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

    app.TopAlg += [ 'DetSimVali/dsv' ]
    dsv = app.algorithm("dsv")
    dsv.Volume = volume
    dsv.OutputLevel = 4
    histsvc = app.service("THistSvc")
    histsvc.Output =["file1 DATAFILE='../mc_ntuple/"+RootOutputfile+"' OPT='RECREATE' TYP='ROOT' "]

    pass


if input:
    rioes = app.service("RootIOEvtSelector")
    rioes.OutputLevel = 4
    app.TopAlg += [ 'DetSimVali/dsv' ]
    dsv = app.algorithm("dsv")
    dsv.Volume = volume
    dsv.OutputLevel = 4
    histsvc = app.service("THistSvc")
    histsvc.Output =["file1 DATAFILE='../mc_ntuple/"+RootOutputfile+"' OPT='RECREATE' TYP='ROOT' "]
    
if output:
    app.TopAlg += [ 'DybStoreAlg/dsa' ]
    dsa = app.algorithm("dsa")
    dsa.OutputLevel = 1
    dybstrgsvs = app.service("DybStoreSvc")

# Run...
app.initialize()
app.run(app.EvtMax)

---