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