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

In This Package:

python » Calibration

SimGe68wBGnoSRC.py

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