| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

python » Calibration

SimGe68wBG.py

Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 
00003 '''
00004 
00005 Script to simulate Ge68 events in the geometry with background sources.
00006 
00007 usage example:
00008 
00009   nuwa.py -n 2 -o output.root -m "Calibration.SimGe68wBG -z 4.0 -a A"
00010 
00011 '''
00012 
00013 from DybPython.DybPythonAlg import DybPythonAlg
00014 from DetHelpers.DetHelpersConf import AutoPositionerTool
00015 
00016 import os, math
00017 
00018 class PositionerAlg(DybPythonAlg):
00019     # An algorithm that places a volume in the detector geometry
00020     def _init_(self,name):
00021         DybPython._init_(self,name)
00022         print "Making PositionerAlg",name
00023         self.posToolConf = None
00024 
00025     def initialize(self):
00026         status = DybPythonAlg.initialize(self)
00027         print "getFullName = ",self.posToolConf.getFullName()
00028         self.posTool = self.tool('IPositionerTool',self.posToolConf.getFullName())
00029 
00030         # Place volume
00031         status = self.posTool.placeVolume()
00032         return status
00033 
00034 
00035 def configure(argv = []):
00036     """Configure this module with source position"""
00037     
00038     import sys, getopt
00039     import GaudiKernel.SystemOfUnits as units
00040     from time import gmtime, mktime, strftime, strptime, timezone
00041     opts,args = getopt.getopt(argv,
00042                               "z:a:")
00043     from DybPython.Control import nuwa
00044     wallTime = nuwa.opts.time
00045     runNumber = nuwa.opts.run
00046     #wallTime = 0
00047     axis = 'A'
00048     acuName = 'ACU_A_Center'
00049     xpos = 0.0
00050     ypos = 0.0
00051     zpos = 0.0
00052     pmtDataPath = None
00053     for opt,arg in opts:
00054         if opt == "-z":
00055             zpos = float(arg) * units.cm
00056             print "======================================================"
00057             print "Source Z position = ", zpos / units.cm, " cm"
00058             print "======================================================"
00059         if opt == "-a":
00060             axis = arg
00061             print "======================================================"
00062             print "ACU Axis = ", axis
00063             print "======================================================"
00064             if axis == 'B':
00065                 xpos = 135.0 * math.cos(112.5 * math.pi / 180.) * units.cm
00066                 ypos = 135.0 * math.sin(112.5 * math.pi / 180.) * units.cm
00067                 acuName = 'ACU_B_GdlsEdge'
00068             elif axis == 'C':
00069                 xpos = 177.25*math.cos( (112.5 + 180) * math.pi / 180.)*units.cm
00070                 ypos = 177.25*math.sin( (112.5 + 180) * math.pi / 180.)*units.cm
00071                 acuName = 'ACU_C_GammaCatcher'
00072     
00073     # Set run info
00074     from RunDataSvc.RunDataSvcConf import RunDataSvc
00075     runDataSvc = RunDataSvc()
00076     runDataSvc.SimRunType = "Calibration"
00077     sourceName = "DayaBayAD1_"+acuName+"_Germanium_68"
00078     runDataSvc.SimCalibSources = [ sourceName ]
00079     runDataSvc.SimCalibZPosition = { sourceName : zpos }    
00080 
00081     # Add pull-mixing simulation
00082     from Stage import Configure as StageConfigure
00083     stageCfg = StageConfigure()
00084     stageCfg.addStages( ['Kinematic','Detector','Electronic','TrigRead',
00085                          'SingleLoader'] )
00086 
00087     # Add kinematic generators
00088     from Gnrtr.GnrtrConf import Gnrtr
00089     from GenTools.Helpers import Gun
00090     from GenTools.Helpers import HepEVT
00091 
00092     SimTime=nuwa.opts.executions * 2 * 0.010
00093 
00094     # Ge68
00095     
00096     seed = runNumber
00097     Ge68lifetime=0.010
00098     Ge68events= SimTime // Ge68lifetime
00099    
00100     mygenGe68 = HepEVT("Ge68.exe -n %d -seed %d |" % (Ge68events,seed), name = "Ge68")
00101     mygenGe68.positioner.Volume = "/dd/Structure/CalibrationSources/db-ad1-Ge68SourceActive"
00102     mygenGe68.transformer.Volume = "/dd/Structure/CalibrationSources/db-ad1-Ge68SourceActive"
00103     mygenGe68.positioner.Position = [0, 0, 0]
00104     mygenGe68.positioner.Strategy = "FullVolume"
00105     mygenGe68.positioner.Mode = "Uniform"
00106     mygenGe68.positioner.Spread = 3.0 * units.cm
00107     mygenGe68.timerator.LifeTime = Ge68lifetime*units.second
00108 
00109     gnrtrGe68 = Gnrtr("gnrtrGe68");
00110     gnrtrGe68.GenTools = mygenGe68.tools()
00111     gnrtrGe68.ThisStageName = "Kinematic"
00112     gnrtrGe68.TimeStamp = int(wallTime)
00113     stageCfg.KinematicSequence.Members.append(gnrtrGe68)
00114     
00115 
00116     # K40 in PMTs
00117     seed = runNumber
00118     K40PMTlifetime=1/689.472
00119     K40PMTevents= SimTime // K40PMTlifetime
00120     mygenK40PMT = HepEVT("K40.exe -n %d -seed %d |" % (K40PMTevents,seed), name = "K40PMT")
00121     mygenK40PMT.positioner.Volume = "/dd/Structure/AD/db-oil1"
00122     mygenK40PMT.transformer.Volume = "/dd/Structure/AD/db-oil1"
00123     mygenK40PMT.positioner.Position = [0, 0, 0]
00124     mygenK40PMT.positioner.Strategy = "VolumeType"
00125     mygenK40PMT.positioner.FillVolumes = ["lvPmtHemiVacuum"]
00126     mygenK40PMT.positioner.Mode = "Uniform"
00127     mygenK40PMT.positioner.Spread = 3. * units.m
00128     mygenK40PMT.timerator.LifeTime = K40PMTlifetime*units.second
00129 
00130     gnrtrK40PMT = Gnrtr("gnrtrK40PMT");
00131     gnrtrK40PMT.GenTools = mygenK40PMT.tools()
00132     gnrtrK40PMT.ThisStageName = "Kinematic"
00133     gnrtrK40PMT.TimeStamp = int(wallTime)
00134     stageCfg.KinematicSequence.Members.append(gnrtrK40PMT)
00135 
00136     # K40 in Steel
00137     seed = runNumber
00138     K40STLlifetime=(1.0/260.0)*(20.0/24.875)
00139     K40STLevents= SimTime // K40STLlifetime
00140     mygenK40STL = HepEVT("K40.exe -n %d -seed %d |" % (K40STLevents,seed), name = "K40STL")
00141     mygenK40STL.positioner.Volume = "/dd/Structure/AD/db-ade1"
00142     mygenK40STL.transformer.Volume = "/dd/Structure/AD/db-ade1"
00143     mygenK40STL.positioner.Position = [0, 0, 0]
00144     #mygenK40STL.positioner.Strategy = "VolumeType"
00145     #mygenK40STL.positioner.FillVolumes = ["lvSST"]
00146     mygenK40STL.positioner.Strategy = "Material"
00147     mygenK40STL.positioner.FillMaterials = ["StainlessSteel"]
00148     mygenK40STL.positioner.Mode = "Uniform"
00149     mygenK40STL.positioner.Spread = 3.0 * units.m
00150     mygenK40STL.timerator.LifeTime = K40STLlifetime*units.second
00151 
00152     gnrtrK40STL = Gnrtr("gnrtrK40STL");
00153     gnrtrK40STL.GenTools = mygenK40STL.tools()
00154     gnrtrK40STL.ThisStageName = "Kinematic"
00155     gnrtrK40STL.TimeStamp = int(wallTime)
00156     stageCfg.KinematicSequence.Members.append(gnrtrK40STL)
00157 
00158     # Thorium in PMTs
00159     thPMTlifetime=1.0/810.929
00160     thPMTMaxEvents = SimTime // thPMTlifetime
00161     thSTLlifetime=1.0/480.0*20.0/24.875
00162     thSTLMaxEvents = SimTime // thSTLlifetime
00163     thPMTFirstEvent = (thPMTMaxEvents + thSTLMaxEvents)*(runNumber - 1001)
00164     while (thPMTFirstEvent + thPMTMaxEvents) > 1000000:
00165         thPMTFirstEvent = thPMTFirstEvent // 5.0
00166     thPMTHepEvtFile = "${UNDERSTANDINGENERGYROOT}/share/th232.asc"
00167     thoriumPMTGen = "${UNDERSTANDINGENERGYROOT}/share/chooseHepEVT.py -f %s -n %d -s %d |" % (thPMTHepEvtFile,thPMTMaxEvents,thPMTFirstEvent)
00168     mygenThPMT = HepEVT(thoriumPMTGen,name = "ThoriumPMT")
00169     mygenThPMT.positioner.Volume = "/dd/Structure/AD/db-oil1"
00170     mygenThPMT.transformer.Volume = "/dd/Structure/AD/db-oil1"
00171     mygenThPMT.positioner.Position = [0, 0, 0]
00172     mygenThPMT.positioner.Strategy = "VolumeType"
00173     mygenThPMT.positioner.FillVolumes = ["lvPmtHemiVacuum"]
00174     mygenThPMT.positioner.Mode = "Uniform"
00175     mygenThPMT.positioner.Spread = 3. * units.m
00176     mygenThPMT.timerator.LifeTime = thPMTlifetime*units.second
00177 
00178     gnrtrTh232PMT = Gnrtr("gnrtrTh232PMT");
00179     gnrtrTh232PMT.GenTools = mygenThPMT.tools()
00180     gnrtrTh232PMT.ThisStageName = "Kinematic"
00181     gnrtrTh232PMT.TimeStamp = int(wallTime)
00182     stageCfg.KinematicSequence.Members.append(gnrtrTh232PMT)
00183 
00184 
00185     # Thorium in Steel
00186     thSTLFirstEvent = (thPMTMaxEvents + thSTLMaxEvents)*(runNumber - 1001) + thPMTMaxEvents
00187     while (thSTLFirstEvent + thSTLMaxEvents) > 1000000:
00188         thSTLFirstEvent = thSTLFirstEvent // 5.0
00189     thHepEvtFile = "${UNDERSTANDINGENERGYROOT}/share/th232.asc"
00190     thoriumSTLGen = "${UNDERSTANDINGENERGYROOT}/share/chooseHepEVT.py -f %s -n %d -s %d |" % (thHepEvtFile,thSTLMaxEvents,thSTLFirstEvent)
00191     mygenThSTL = HepEVT(thoriumSTLGen,name = "ThoriumSTL")
00192     mygenThSTL.positioner.Volume = "/dd/Structure/AD/db-ade1"
00193     mygenThSTL.transformer.Volume = "/dd/Structure/AD/db-ade1"
00194     mygenThSTL.positioner.Position = [0, 0, 0]
00195     mygenThSTL.positioner.Strategy = "Material"
00196     mygenThSTL.positioner.FillMaterials = ["StainlessSteel"]
00197     mygenThSTL.positioner.Mode = "Uniform"
00198     mygenThSTL.positioner.Spread = 3. * units.m
00199     mygenThSTL.timerator.LifeTime = thSTLlifetime*units.second
00200 
00201     gnrtrTh232STL = Gnrtr("gnrtrTh232STL");
00202     gnrtrTh232STL.GenTools = mygenThSTL.tools()
00203     gnrtrTh232STL.ThisStageName = "Kinematic"
00204     gnrtrTh232STL.TimeStamp = int(wallTime)
00205     stageCfg.KinematicSequence.Members.append(gnrtrTh232STL)
00206 
00207     # Uranium in PMTs
00208     uPMTlifetime=1.0/2823.805
00209     uPMTMaxEvents = SimTime // uPMTlifetime
00210     uSTLlifetime=1.0/216.0*20.0/24.875
00211     uSTLMaxEvents = SimTime // uSTLlifetime
00212     uPMTFirstEvent = (uPMTMaxEvents + uSTLMaxEvents)*(runNumber - 1001)
00213     while (uPMTFirstEvent + uPMTMaxEvents) > 1000000:
00214         uPMTFirstEvent = uPMTFirstEvent // 5.0
00215     uPMTHepEvtFile = "${UNDERSTANDINGENERGYROOT}/share/u238.asc"
00216     uraniumPMTGen = "${UNDERSTANDINGENERGYROOT}/share/chooseHepEVT.py -f %s -n %d -s %d |" % (uPMTHepEvtFile,uPMTMaxEvents,uPMTFirstEvent)
00217     mygenUPMT = HepEVT(uraniumPMTGen,name = "UraniumPMT")
00218     mygenUPMT.positioner.Volume = "/dd/Structure/AD/db-oil1"
00219     mygenUPMT.transformer.Volume = "/dd/Structure/AD/db-oil1"
00220     mygenUPMT.positioner.Position = [0, 0, 0]
00221     mygenUPMT.positioner.Strategy = "VolumeType"
00222     mygenUPMT.positioner.FillVolumes = ["lvPmtHemiVacuum"]
00223     mygenUPMT.positioner.Mode = "Uniform"
00224     mygenUPMT.positioner.Spread = 3. * units.m
00225     mygenUPMT.timerator.LifeTime = uPMTlifetime*units.second
00226 
00227     gnrtrU238PMT = Gnrtr("gnrtrU232PMT");
00228     gnrtrU238PMT.GenTools = mygenUPMT.tools()
00229     gnrtrU238PMT.ThisStageName = "Kinematic"
00230     gnrtrU238PMT.TimeStamp = int(wallTime)
00231     stageCfg.KinematicSequence.Members.append(gnrtrU238PMT)
00232 
00233 
00234     # Uranium in Steel
00235     uSTLFirstEvent = (uPMTMaxEvents + uSTLMaxEvents)*(runNumber - 1001) + uPMTMaxEvents
00236     while (uSTLFirstEvent + uSTLMaxEvents) > 1000000:
00237         uSTLFirstEvent = uSTLFirstEvent // 5.0
00238     uSTLHepEvtFile = "${UNDERSTANDINGENERGYROOT}/share/u238.asc"
00239     uraniumSTLGen = "${UNDERSTANDINGENERGYROOT}/share/chooseHepEVT.py -f %s -n %d -s %d |" % (uSTLHepEvtFile,uSTLMaxEvents,uSTLFirstEvent)
00240     mygenUSTL = HepEVT(uraniumSTLGen,name = "UraniumSTL")
00241     mygenUSTL.positioner.Volume = "/dd/Structure/AD/db-ade1"
00242     mygenUSTL.transformer.Volume = "/dd/Structure/AD/db-ade1"
00243     mygenUSTL.positioner.Position = [0, 0, 0]
00244     #mygenUSTL.positioner.Strategy = "VolumeType"
00245     #mygenUSTL.positioner.FillVolumes = ["lvSST"]
00246     mygenUSTL.positioner.Strategy = "Material"
00247     mygenUSTL.positioner.FillMaterials = ["StainlessSteel"]
00248     mygenUSTL.positioner.Mode = "Uniform"
00249     mygenUSTL.positioner.Spread = 3. * units.m
00250     mygenUSTL.timerator.LifeTime = uSTLlifetime*units.second
00251 
00252     gnrtrU238STL = Gnrtr("gnrtrU232STL");
00253     gnrtrU238STL.GenTools = mygenUSTL.tools()
00254     gnrtrU238STL.ThisStageName = "Kinematic"
00255     gnrtrU238STL.TimeStamp = int(wallTime)
00256     stageCfg.KinematicSequence.Members.append(gnrtrU238STL)
00257 
00258     # Rest of Simulation stages
00259     import DetSim
00260     detsim = DetSim.Configure(use_push_algs = False)
00261     detsim.historian(trackSelection="(pdg == -11)",\
00262                      vertexSelection="(pdg == -11)")
00263     params = {
00264         'start' :"(start > 0)",
00265         'track1':"(id==1)",
00266         'track2':"(id==2)",
00267         'GD':    "MaterialName == '/dd/Materials/GdDopedLS'",
00268         'LS':    "MaterialName == '/dd/Materials/LiquidScintillator'",
00269         'MO':   "MaterialName == '/dd/Materials/MineralOil'",
00270         'IAV':   "DetectorElementName == 'db-iav1'",
00271         'OAV':   "DetectorElementName == 'db-oav1'",
00272         'IWS': "MaterialName == '/dd/Materials/IwsWater'",
00273         'OWS': "MaterialName == '/dd/Materials/OwsWater'",
00274         'lastvtx': "IsStopping == 1",
00275         'firstvtx': "IsStarting == 1",
00276         'Neutron': "pdg == 2112",
00277         'NeutronMom': "creator == 2112",
00278         'Gamma': "pdg == 22",
00279         'Positron': "pdg == -11",
00280         'Muon': "(pdg == 13 or pdg == -13)"
00281         }
00282 
00283     detsim.unobserver(stats=[
00284             ["EDepInGdLS", "dE", "%(GD)s"%params],
00285             ["EDepInLS", "dE", "%(LS)s"%params],
00286             ["EDepInIAV", "dE", "%(IAV)s"%params],
00287             ["EDepInOAV", "dE", "%(OAV)s"%params],
00288             ["EDepInOIL", "dE", "%(MO)s"%params],
00289 
00290             ["QEDepInGdLS", "qdE", "%(GD)s"%params],
00291             ["QEDepInLS", "qdE", "%(LS)s"%params],
00292             ["QEDepInIAV", "qdE", "%(IAV)s"%params],
00293             ["QEDepInOAV", "qdE", "%(OAV)s"%params],
00294             ["QEDepInOIL", "qdE", "%(MO)s"%params],
00295 
00296             ["tQESumGdLS", "qEt", "%(GD)s"%params],
00297             ["xQESumGdLS", "qEx", "%(GD)s"%params],
00298             ["yQESumGdLS", "qEy", "%(GD)s"%params],
00299             ["zQESumGdLS", "qEz", "%(GD)s"%params],
00300 
00301             ["tQESumLS", "qEt", "%(LS)s"%params],
00302             ["xQESumLS", "qEx", "%(LS)s"%params],
00303             ["yQESumLS", "qEy", "%(LS)s"%params],
00304             ["zQESumLS", "qEz", "%(LS)s"%params],
00305 
00306             ["tQESumMO", "qEt", "%(MO)s"%params],
00307             ["xQESumMO", "qEx", "%(MO)s"%params],
00308             ["yQESumMO", "qEy", "%(MO)s"%params],
00309             ["zQESumMO", "qEz", "%(MO)s"%params],
00310             
00311             # track 1
00312             ["pdgId_Trk1","pdg","%(track1)s and %(start)s"%params],
00313             ["t_Trk1",    "t" , "%(track1)s and %(start)s"%params],
00314             ["x_Trk1",    "x", "%(track1)s and %(start)s"%params],
00315             ["y_Trk1",    "y", "%(track1)s and %(start)s"%params],
00316             ["z_Trk1",    "z", "%(track1)s and %(start)s"%params],
00317             ["tEnd_Trk1",    "t" , "%(track1)s and %(lastvtx)s"%params],
00318             ["xEnd_Trk1",    "x", "%(track1)s and %(lastvtx)s"%params],
00319             ["yEnd_Trk1",    "y", "%(track1)s and %(lastvtx)s"%params],
00320             ["zEnd_Trk1",    "z", "%(track1)s and %(lastvtx)s"%params],
00321             ["e_Trk1",    "E",  "%(track1)s and %(start)s"%params],
00322             ["p_Trk1",    "p",  "%(track1)s and %(start)s"%params],
00323             ["ke_Trk1",   "KE", "%(track1)s and %(start)s"%params],
00324             ["vx_Trk1",   "vx","%(track1)s and %(start)s"%params],
00325             ["vy_Trk1",   "vy","%(track1)s and %(start)s"%params],
00326             ["vz_Trk1",   "vz","%(track1)s and %(start)s"%params],
00327             ["TrkLength_GD_Trk1",  "dx","%(track1)s and %(GD)s"%params],
00328             ["TrkLength_iAV_Trk1", "dx","%(track1)s and %(IAV)s"%params],
00329             ["TrkLength_LS_Trk1",  "dx","%(track1)s and %(LS)s"%params],
00330             ["TrkLength_oAV_Trk1", "dx","%(track1)s and %(OAV)s"%params],
00331             ["TrkLength_Oil_Trk1", "dx","%(track1)s and %(MO)s"%params]
00332             ])
00333 
00334     from DetSimProc.DetSimProcConf import DetSimProc
00335     dsp = DetSimProc()
00336     dsp.ThisStageName = "Detector"
00337     dsp.LowerStageName = "Kinematic"
00338     stageCfg.DetectorSequence.Members.append(dsp)
00339 
00340     import ElecSim
00341     elecsim = ElecSim.Configure(use_push_algs = False)
00342     from ElecSimProc.ElecSimProcConf import ElecSimProc
00343     esp = ElecSimProc()
00344     esp.ThisStageName = "Electronic"
00345     esp.LowerStageName = "Detector"
00346     stageCfg.ElectronicSequence.Members.append(esp)
00347     
00348     #import TrigSim
00349     #trigsim = TrigSim.Configure()
00350     from TrigReadProc.TrigReadProcConf import TrigReadProc
00351     tsp = TrigReadProc()
00352     tsp.ThisStageName = "TrigRead"
00353     tsp.LowerStageName = "Electronic"
00354     stageCfg.TrigReadSequence.Members.append(tsp)
00355 
00356     #import ReadoutSim
00357     #rosim = ReadoutSim.Configure()
00358     from SingleLoader.SingleLoaderConf import SingleLoader
00359     sll = SingleLoader()
00360     sll.ThisStageName = "SingleLoader"
00361     sll.LowerStageName = "TrigRead"
00362     stageCfg.SingleLoaderSequence.Members.append(sll)
00363 
00364     from Stage.StageConf import Sim15
00365     sim15=Sim15()
00366     sim15.TopStage="SingleLoader"
00367         
00368     from Gaudi.Configuration import ApplicationMgr
00369     theApp = ApplicationMgr()
00370     theApp.TopAlg.append(sim15)
00371     
00372     global autoPositioner
00373     autoPositioner = AutoPositionerTool("posAlg.AutoPositioner")
00374     autoPositioner.PhysicalVolume = "pvGe68SourceAssy"
00375     autoPositioner.LogicalVolume = "/dd/Geometry/CalibrationSources/lvGe68SourceAssy"
00376     autoPositioner.CoordinateDetElem = "/dd/Structure/AD/db-oil1"
00377     #autoPositioner.Position = [0,0,0]
00378     autoPositioner.Position = [xpos, ypos, zpos]
00379     autoPositioner.Element = "db-ad1-Ge68SourceAssy"
00380     autoPositioner.ElementPath = "/dd/Structure/CalibrationSources"
00381     
00382     autoPositioner.SubDetectorElements = [["db-ad1-Ge68SourceActive",
00383                                            "/dd/Structure/CalibrationSources/db-ad1-Ge68SourceAssy",
00384                                            "/dd/Geometry/CalibrationSources/lvGe68SourceActive",
00385                                            "pvGe68SourceShell/pvGe68Source/pvGe68SourceInterior/pvGe68SourceActive"]]                                           
00386     
00387     
00388 def run(app):
00389      
00390     global autoPositioner
00391     posAlg = PositionerAlg("posAlg")
00392     posAlg.posToolConf = autoPositioner
00393     print "dir(app)= ",dir(app)
00394     print "dir(app._appMgr)= ",dir(app._appmgr)
00395     app.addAlgorithm(posAlg)
00396     print "Algorithms:"
00397     for alg in app.algorithms():
00398         print "    ",alg
00399     topAlg = app.topAlg
00400     app.topAlg = topAlg[-1:] + topAlg[:-1]
00401     print "Algorithms (fixed):"
00402     for alg in app.algorithms():
00403         print "    ",alg
00404     
00405     pass
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |
Generated on Mon Apr 11 20:11:24 2011 for Calibration by doxygen 1.4.7