00001
00002
00003 '''
00004 usage example:
00005
00006 nuwa.py -A -n 2 -o output.root -m "MDC09a.runPositron -z 4.0 -a 0"
00007
00008 '''
00009
00010 DATETIME_FORMAT = '%Y-%m-%dT%H:%M:%S'
00011
00012 import os, math
00013
00014 def configure(argv = []):
00015 """Configure this module with positron position"""
00016
00017 import sys, getopt
00018 from time import gmtime, mktime, strftime, strptime, timezone
00019 opts,args = getopt.getopt(argv,"p:w:z:a:k:M:s:n:")
00020 wallTime = 0
00021 axis = 0
00022 xpos = 0.0
00023 ypos = 0.0
00024 zpos = 0.0
00025 pmtDataPath = None
00026 mode = ""
00027 energy = 1.0
00028 seed = "42"
00029 nevts = "10000"
00030 volume = "/dd/Structure/AD/db-oil1"
00031
00032 for opt,arg in opts:
00033 if opt == "-p":
00034 pmtDataPath = arg
00035 if opt == "-w":
00036 if -1 != arg.find('T'):
00037 wallTime = int(mktime(strptime(arg,
00038 DATETIME_FORMAT)) - timezone)
00039 else:
00040 wallTime = int(arg)
00041 print "======================================================"
00042 print "Wall clock start time = ", strftime(DATETIME_FORMAT,
00043 gmtime(float(wallTime)))
00044 print "======================================================"
00045 if opt == "-z":
00046 zpos = float(arg)
00047 print "======================================================"
00048 print "Source Z position = ", zpos, " cm"
00049 print "======================================================"
00050 if opt == "-a":
00051 axis = int(arg)
00052 print "======================================================"
00053 print "Source Axis = ", axis, " cm"
00054 print "======================================================"
00055 if axis == 1:
00056 xpos = 135.0 * math.cos(112.5 * math.pi / 180.)
00057 ypos = 135.0 * math.sin(112.5 * math.pi / 180.)
00058 elif axis == 2:
00059 xpos = 177.25 * math.cos( (112.5 + 180) * math.pi / 180.)
00060 ypos = 177.25 * math.sin( (112.5 + 180) * math.pi / 180.)
00061 if opt == "-k":
00062 kinetic = float(arg)
00063 if opt == "-M":
00064 mode = arg
00065 if opt == "-s":
00066 seed = arg
00067 if opt == "-n":
00068 nevts = arg
00069
00070 import GenTools
00071 from GenTools.Helpers import Gun
00072 from GenTools.Helpers import HepEVT
00073 import GaudiKernel.SystemOfUnits as units
00074 gtc = GenTools.Configure()
00075 gtc.generator.TimeStamp = int(wallTime)
00076
00077 if mode == "" or mode == "Cal":
00078 print "Simulating in calibration mode."
00079 mygun = Gun()
00080 mygun.gun.ParticleName = 'e+'
00081 mygun.timerator.LifeTime = 0.020*units.second
00082 mygun.setVolume(volume)
00083 mygun.gun.Momentum = 1.0*units.eV
00084 mygun.positioner.Position = [xpos*units.cm, ypos*units.cm, zpos*units.cm]
00085 gtc.register(mygun)
00086 elif mode == "Mono":
00087 print "Simulating in monoenegetic mode."
00088 mygun = Gun()
00089 mygun.gun.ParticleName = 'e+'
00090 mygun.timerator.LifeTime = 0.020*units.second
00091 mygun.setVolume(volume)
00092 energy = (kinetic + 0.5109989)*units.MeV
00093 momentum = math.sqrt((energy)**2 - (0.5109989*units.MeV)**2)
00094 mygun.gun.Momentum = momentum
00095 mygun.positioner.Strategy = "FullVolume"
00096 mygun.positioner.Mode = "Uniform"
00097 mygun.positioner.Spread = 2.6*units.meter
00098 mygun.positioner.Position = [0,0,0*units.meter]
00099 gtc.register(mygun)
00100 elif mode == "IBD":
00101 print "Simulating in IBD mode."
00102 ibd = "InverseBeta.exe -eplus_only -n " + nevts + " -seed " + seed + " |"
00103 he = HepEVT(hepEvtDataSource = ibd)
00104 he.positioner.Volume = volume
00105 he.positioner.Strategy = "FullVolume"
00106 he.positioner.Mode = "Uniform"
00107 he.positioner.Position = [0,0,0]
00108 he.positioner.Spread = 2.6*units.meter
00109 he.transformer.Volume = volume
00110 gtc.register(he)
00111
00112 import DetSim
00113 detsim = DetSim.Configure(physlist = DetSim.physics_list_basic)
00114 detsim.historian(trackSelection="(pdg == -11)",\
00115 vertexSelection="(pdg == -11)")
00116 params = {
00117 'start' :"(start > 0)",
00118 'track1':"(id==1)",
00119 'track2':"(id==2)",
00120 'GD': "MaterialName == '/dd/Materials/GdDopedLS'",
00121 'LS': "MaterialName == '/dd/Materials/LiquidScintillator'",
00122 'MO': "MaterialName == '/dd/Materials/MineralOil'",
00123 'IAV': "DetectorElementName == 'db-iav1'",
00124 'OAV': "DetectorElementName == 'db-oav1'",
00125 'IWS': "MaterialName == '/dd/Materials/IwsWater'",
00126 'OWS': "MaterialName == '/dd/Materials/OwsWater'",
00127 'lastvtx': "IsStopping == 1",
00128 'firstvtx': "IsStarting == 1",
00129 'Neutron': "pdg == 2112",
00130 'NeutronMom': "creator == 2112",
00131 'Gamma': "pdg == 22",
00132 'Positron': "pdg == -11",
00133 'Muon': "(pdg == 13 or pdg == -13)"
00134 }
00135
00136 detsim.unobserver(stats=[
00137 ["EDepInGdLS", "dE", "%(GD)s"%params],
00138 ["EDepInLS", "dE", "%(LS)s"%params],
00139 ["EDepInIAV", "dE", "%(IAV)s"%params],
00140 ["EDepInOAV", "dE", "%(OAV)s"%params],
00141 ["EDepInOIL", "dE", "%(MO)s"%params],
00142
00143 ["QEDepInGdLS", "qdE", "%(GD)s"%params],
00144 ["QEDepInLS", "qdE", "%(LS)s"%params],
00145 ["QEDepInIAV", "qdE", "%(IAV)s"%params],
00146 ["QEDepInOAV", "qdE", "%(OAV)s"%params],
00147 ["QEDepInOIL", "qdE", "%(MO)s"%params],
00148
00149 ["tQESumGdLS", "qEt", "%(GD)s"%params],
00150 ["xQESumGdLS", "qEx", "%(GD)s"%params],
00151 ["yQESumGdLS", "qEy", "%(GD)s"%params],
00152 ["zQESumGdLS", "qEz", "%(GD)s"%params],
00153
00154 ["tQESumLS", "qEt", "%(LS)s"%params],
00155 ["xQESumLS", "qEx", "%(LS)s"%params],
00156 ["yQESumLS", "qEy", "%(LS)s"%params],
00157 ["zQESumLS", "qEz", "%(LS)s"%params],
00158
00159 ["tQESumMO", "qEt", "%(MO)s"%params],
00160 ["xQESumMO", "qEx", "%(MO)s"%params],
00161 ["yQESumMO", "qEy", "%(MO)s"%params],
00162 ["zQESumMO", "qEz", "%(MO)s"%params],
00163
00164
00165 ["pdgId_Trk1","pdg","%(track1)s and %(start)s"%params],
00166 ["t_Trk1", "t" , "%(track1)s and %(start)s"%params],
00167 ["x_Trk1", "x", "%(track1)s and %(start)s"%params],
00168 ["y_Trk1", "y", "%(track1)s and %(start)s"%params],
00169 ["z_Trk1", "z", "%(track1)s and %(start)s"%params],
00170 ["tEnd_Trk1", "t" , "%(track1)s and %(lastvtx)s"%params],
00171 ["xEnd_Trk1", "x", "%(track1)s and %(lastvtx)s"%params],
00172 ["yEnd_Trk1", "y", "%(track1)s and %(lastvtx)s"%params],
00173 ["zEnd_Trk1", "z", "%(track1)s and %(lastvtx)s"%params],
00174 ["e_Trk1", "E", "%(track1)s and %(start)s"%params],
00175 ["p_Trk1", "p", "%(track1)s and %(start)s"%params],
00176 ["ke_Trk1", "KE", "%(track1)s and %(start)s"%params],
00177 ["vx_Trk1", "lvx","%(track1)s and %(start)s"%params],
00178 ["vy_Trk1", "lvy","%(track1)s and %(start)s"%params],
00179 ["vz_Trk1", "lvz","%(track1)s and %(start)s"%params],
00180 ["TrkLength_GD_Trk1", "dx","%(track1)s and %(GD)s"%params],
00181 ["TrkLength_iAV_Trk1", "dx","%(track1)s and %(IAV)s"%params],
00182 ["TrkLength_LS_Trk1", "dx","%(track1)s and %(LS)s"%params],
00183 ["TrkLength_oAV_Trk1", "dx","%(track1)s and %(OAV)s"%params],
00184 ["TrkLength_Oil_Trk1", "dx","%(track1)s and %(MO)s"%params]
00185 ])
00186
00187 import ElecSim
00188 elecsim = ElecSim.Configure()
00189 if pmtDataPath != None:
00190
00191 elecsim.dataSvc.setPmtSimData( pmtDataPath )
00192
00193 import TrigSim
00194 trigsim = TrigSim.Configure()
00195
00196 import ReadoutSim
00197 rosim = ReadoutSim.Configure()
00198
00199 def run(app):
00200 pass
