DDG4/examples/SiDSim.py

0001 # ==========================================================================
0002 #  AIDA Detector description implementation
0003 # --------------------------------------------------------------------------
0004 # Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 # All rights reserved.
0006 #
0007 # For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 # For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 #
0010 # ==========================================================================
0011 #
0012 #
0013 from __future__ import absolute_import, unicode_literals
0014 import os
0015 import sys
0016 import time
0017 import logging
0018 import DDG4
0019 from DDG4 import OutputLevel as Output
0020 from g4units import keV, GeV, mm, ns, MeV
0021 #
0022 #
0023 logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO)
0024 logger = logging.getLogger(__name__)
0025
0026 """
0027
0028    dd4hep simulation example setup using the python configuration
0029
0030    @author  M.Frank
0031    @version 1.0
0032
0033 """
0034
0035
0036 def show_help():
0037   logging.info("SiDSim.py -option [-option]                           ")
0038   logging.info("       -vis   <file>            Enable visualization  ")
0039   logging.info("                                Macro file is optional")
0040   logging.info("       -macro <file>            Start G4 macro        ")
0041   logging.info("       -batch                   Batch execution       ")
0042   logging.info("       -events <number>         If batch: number of events to be executed")
0043
0044
0045 def run():
0046   args = DDG4.CommandLine()
0047   #
0048   if args.help or args.h:
0049     show_help()
0050     sys.exit(1)
0051
0052   kernel = DDG4.Kernel()
0053   description = kernel.detectorDescription()
0054   install_dir = os.environ['DD4hepINSTALL']
0055   kernel.loadGeometry(str("file:" + install_dir + "/DDDetectors/compact/SiD.xml"))
0056
0057   if args.smartless:
0058     description.worldVolume().setSmartlessValue(int(args.smartless))
0059   DDG4.importConstants(description)
0060
0061   geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction')
0062   geant4.printDetectors()
0063   logger.info("#  Configure UI")
0064   ui = geant4.setupCshUI(macro=args.macro, vis=args.vis)
0065   kernel.UI = 'UI'
0066
0067   import pdb
0068   pdb.set_trace()
0069
0070   cmds = []
0071   if args.verbose:
0072     cmds.append('/run/verbose ' + str(args.verbose))
0073   if args.batch:
0074     if not args.events:
0075       args.events = '5'
0076     cmds.append('/run/beamOn ' + str(args.events))
0077
0078   # Terminate sequence
0079   if args.batch:
0080     cmds.append('/ddg4/UI/terminate')
0081
0082   if len(cmds) > 0:
0083     ui.Commands = cmds
0084
0085   logger.info("#  Configure G4 magnetic field tracking")
0086   geant4.setupTrackingField()
0087
0088   logger.info("#  Setup random generator")
0089   rndm = DDG4.Action(kernel, 'Geant4Random/Random')
0090   rndm.Seed = 987654321
0091   if args.seed_time:
0092     rndm.Seed = int(time.time())
0093   rndm.initialize()
0094   # rndm.showStatus()
0095
0096   logger.info("#  Configure Run actions")
0097   run1 = DDG4.RunAction(kernel, 'Geant4TestRunAction/RunInit')
0098   run1.Property_int = 12345
0099   run1.Property_double = -5e15 * keV
0100   run1.Property_string = 'Startrun: Hello_2'
0101   logger.info("%s %s %s", run1.Property_string, str(run1.Property_double), str(run1.Property_int))
0102   run1.enableUI()
0103   kernel.registerGlobalAction(run1)
0104   kernel.runAction().adopt(run1)
0105
0106   logger.info("#  Configure Event actions")
0107   prt = DDG4.EventAction(kernel, 'Geant4ParticlePrint/ParticlePrint')
0108   prt.OutputLevel = Output.INFO
0109   prt.OutputType = 3  # Print both: table and tree
0110   kernel.eventAction().adopt(prt)
0111
0112   logger.info("""
0113   Configure I/O
0114   """)
0115   # evt_lcio = geant4.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
0116   # evt_lcio.OutputLevel = Output.ERROR
0117
0118   geant4.setupROOTOutput('RootOutput', 'CLICSiD_' + time.strftime('%Y-%m-%d_%H-%M'))
0119
0120   gen = DDG4.GeneratorAction(kernel, "Geant4GeneratorActionInit/GenerationInit")
0121   kernel.generatorAction().adopt(gen)
0122
0123   # VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
0124   logger.info("""
0125   Generation of isotrope tracks of a given multiplicity with overlay:
0126   """)
0127   logger.info("#  First particle generator: pi+")
0128   gen = DDG4.GeneratorAction(kernel, "Geant4IsotropeGenerator/IsotropPi+")
0129   gen.Mask = 1
0130   gen.Particle = 'pi+'
0131   gen.Energy = 100 * GeV
0132   gen.Multiplicity = 2
0133   gen.Distribution = 'cos(theta)'
0134   kernel.generatorAction().adopt(gen)
0135   logger.info("#  Install vertex smearing for this interaction")
0136   gen = DDG4.GeneratorAction(kernel, "Geant4InteractionVertexSmear/SmearPi+")
0137   gen.Mask = 1
0138   gen.Offset = (20 * mm, 10 * mm, 10 * mm, 0 * ns)
0139   gen.Sigma = (4 * mm, 1 * mm, 1 * mm, 0 * ns)
0140   kernel.generatorAction().adopt(gen)
0141
0142   logger.info("#  Second particle generator: e-")
0143   gen = DDG4.GeneratorAction(kernel, "Geant4IsotropeGenerator/IsotropE-")
0144   gen.Mask = 2
0145   gen.Particle = 'e+'
0146   gen.Energy = 25 * GeV
0147   gen.Multiplicity = 2
0148   gen.Distribution = 'uniform'
0149   kernel.generatorAction().adopt(gen)
0150   logger.info("  Install vertex smearing for this interaction")
0151   gen = DDG4.GeneratorAction(kernel, "Geant4InteractionVertexSmear/SmearE-")
0152   gen.Mask = 2
0153   gen.Offset = (-20 * mm, -10 * mm, -10 * mm, 0 * ns)
0154   gen.Sigma = (12 * mm, 8 * mm, 8 * mm, 0 * ns)
0155   kernel.generatorAction().adopt(gen)
0156
0157   logger.info("#  Second particle generator: mu+")
0158   gen = DDG4.GeneratorAction(kernel, "Geant4IsotropeGenerator/IsotropMu+")
0159   gen.Mask = 3
0160   gen.Particle = 'mu+'
0161   gen.Energy = 100 * GeV
0162   gen.Multiplicity = 3
0163   gen.Distribution = 'uniform'
0164   kernel.generatorAction().adopt(gen)
0165   # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0166
0167   logger.info("#  Merge all existing interaction records")
0168   gen = DDG4.GeneratorAction(kernel, "Geant4InteractionMerger/InteractionMerger")
0169   gen.OutputLevel = 4  # generator_output_level
0170   gen.enableUI()
0171   kernel.generatorAction().adopt(gen)
0172   #
0173   logger.info("#  Finally generate Geant4 primaries")
0174   gen = DDG4.GeneratorAction(kernel, "Geant4PrimaryHandler/PrimaryHandler")
0175   gen.OutputLevel = 4  # generator_output_level
0176   gen.enableUI()
0177   kernel.generatorAction().adopt(gen)
0178   #
0179   logger.info("#  ....and handle the simulation particles.")
0180   part = DDG4.GeneratorAction(kernel, "Geant4ParticleHandler/ParticleHandler")
0181   kernel.generatorAction().adopt(part)
0182   # part.SaveProcesses = ['conv','Decay']
0183   part.SaveProcesses = ['Decay']
0184   part.MinimalKineticEnergy = 100 * MeV
0185   part.OutputLevel = 5  # generator_output_level
0186   part.enableUI()
0187   user = DDG4.Action(kernel, "Geant4TCUserParticleHandler/UserParticleHandler")
0188   user.TrackingVolume_Zmax = DDG4.EcalEndcap_zmin
0189   user.TrackingVolume_Rmax = DDG4.EcalBarrel_rmin
0190   user.enableUI()
0191   part.adopt(user)
0192   #
0193   logger.info("#  Setup global filters fur use in sensitive detectors")
0194   f1 = DDG4.Filter(kernel, 'GeantinoRejectFilter/GeantinoRejector')
0195   f2 = DDG4.Filter(kernel, 'ParticleRejectFilter/OpticalPhotonRejector')
0196   f2.particle = 'opticalphoton'
0197   f3 = DDG4.Filter(kernel, 'ParticleSelectFilter/OpticalPhotonSelector')
0198   f3.particle = 'opticalphoton'
0199   f4 = DDG4.Filter(kernel, 'EnergyDepositMinimumCut')
0200   f4.Cut = 10 * MeV
0201   f4.enableUI()
0202   kernel.registerGlobalFilter(f1)
0203   kernel.registerGlobalFilter(f2)
0204   kernel.registerGlobalFilter(f3)
0205   kernel.registerGlobalFilter(f4)
0206   #
0207   logger.info("#  First the tracking detectors")
0208   seq, act = geant4.setupTracker('SiVertexBarrel')
0209   seq.adopt(f1)
0210   act.adopt(f1)
0211   #
0212   seq, act = geant4.setupTracker('SiVertexEndcap')
0213   seq.adopt(f1)
0214   #
0215   seq, act = geant4.setupTracker('SiTrackerBarrel')
0216   seq, act = geant4.setupTracker('SiTrackerEndcap')
0217   seq, act = geant4.setupTracker('SiTrackerForward')
0218   logger.info("#  Now setup the calorimeters")
0219   seq, act = geant4.setupCalorimeter('EcalBarrel')
0220   seq, act = geant4.setupCalorimeter('EcalEndcap')
0221   seq, act = geant4.setupCalorimeter('HcalBarrel')
0222   seq, act = geant4.setupCalorimeter('HcalEndcap')
0223   seq, act = geant4.setupCalorimeter('HcalPlug')
0224   seq, act = geant4.setupCalorimeter('MuonBarrel')
0225   seq, act = geant4.setupCalorimeter('MuonEndcap')
0226   seq, act = geant4.setupCalorimeter('LumiCal')
0227   seq, act = geant4.setupCalorimeter('BeamCal')
0228   #
0229   logger.info("#  Now build the physics list:")
0230   phys = geant4.setupPhysics('QGSP_BERT')
0231   ph = geant4.addPhysics(str('Geant4PhysicsList/Myphysics'))
0232   ph.addPhysicsConstructor(str('G4StepLimiterPhysics'))
0233   #
0234   # Add special particle types from specialized physics constructor
0235   part = geant4.addPhysics(str('Geant4ExtraParticles/ExtraParticles'))
0236   part.pdgfile = os.path.join(install_dir, 'examples/DDG4/examples/particle.tbl')
0237   #
0238   # Add global range cut
0239   rg = geant4.addPhysics(str('Geant4DefaultRangeCut/GlobalRangeCut'))
0240   rg.RangeCut = 0.7 * mm
0241   #
0242   phys.dump()
0243   #
0244   #
0245   if ui and args.vis:
0246     cmds = []
0247     cmds.append('/control/verbose 2')
0248     cmds.append('/run/initialize')
0249     cmds.append('/vis/open OGL')
0250     cmds.append('/vis/verbose errors')
0251     cmds.append('/vis/drawVolume')
0252     cmds.append('/vis/viewer/set/viewpointThetaPhi 55. 45.')
0253     cmds.append('/vis/scene/add/axes 0 0 0 10 m')
0254     ui.Commands = cmds
0255
0256   kernel.configure()
0257   kernel.initialize()
0258
0259   # DDG4.setPrintLevel(Output.DEBUG)
0260   kernel.run()
0261   kernel.terminate()
0262
0263
0264 if __name__ == "__main__":
0265   run()