LHeD/scripts/LheSimu.py

0001 """
0002
0003    DD4hep simulation example setup using the python configuration
0004
0005    @author  M.Frank
0006    @version 1.0
0007    modified for LHeC
0008
0009 """
0010 from __future__ import absolute_import, unicode_literals
0011 import logging
0012
0013 logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO)
0014 logger = logging.getLogger(__name__)
0015
0016
0017 def run():
0018   import LHeD
0019   import DDG4
0020   import os
0021   import g4units
0022   from DDG4 import OutputLevel as Output
0023
0024   lhed = LHeD.LHeD()
0025   geant4 = lhed.geant4
0026   kernel = lhed.kernel
0027   lhed.loadGeometry()
0028   geant4.printDetectors()
0029   kernel.UI = "UI"
0030   geant4.setupCshUI()
0031   lhed.setupField(quiet=False)
0032   DDG4.importConstants(kernel.detectorDescription(), debug=False)
0033
0034   dd4hep_dir = os.environ['DD4hep_DIR']
0035   kernel.loadXML("file:" + dd4hep_dir + "/examples/LHeD/scripts/DDG4_field.xml")
0036
0037   geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction')
0038   geant4.printDetectors()
0039
0040   # Configure G4 magnetic field tracking
0041   field = geant4.addConfig('Geant4FieldTrackingSetupAction/MagFieldTrackingSetup')
0042   field.stepper = "HelixGeant4Runge"
0043   field.equation = "Mag_UsualEqRhs"
0044   field.eps_min = 5e-0520 * g4units.mm
0045   field.eps_max = 0.001 * g4units.mm
0046   field.min_chord_step = 0.01 * g4units.mm
0047   field.delta_chord = 0.25 * g4units.mm
0048   field.delta_intersection = 1e-05 * g4units.mm
0049   field.delta_one_step = 0.001 * g4units.mm
0050   logger.info('+++++> %s -> stepper  = %s', field.name, field.stepper)
0051   logger.info('+++++> %s -> equation = %s', field.name, field.equation)
0052   logger.info('+++++> %s -> eps_min  = %s', field.name, field.eps_min)
0053   logger.info('+++++> %s -> eps_max  = %s', field.name, field.eps_max)
0054   logger.info('+++++> %s -> delta_one_step = %s', field.name, field.delta_one_step)
0055
0056   """
0057   # Setup random generator
0058   rndm = DDG4.Action(kernel,'Geant4Random/Random')
0059   rndm.Seed = 987654321
0060   rndm.initialize()
0061   rndm.showStatus()
0062   rndm.Seed = 987654321
0063   """
0064
0065   # Configure Run actions
0066   run1 = DDG4.RunAction(kernel, 'Geant4TestRunAction/RunInit')
0067   """
0068   run1.Property_int    = 12345
0069   run1.Property_double = -5e15*keV
0070   run1.Property_string = 'Startrun: Hello_LHeD'
0071   """
0072   run1.enableUI()
0073   kernel.registerGlobalAction(run1)
0074   kernel.runAction().adopt(run1)
0075
0076   # Configure Event actions
0077   prt = DDG4.EventAction(kernel, 'Geant4ParticlePrint/ParticlePrint')
0078   prt.OutputLevel = Output.INFO
0079   prt.OutputType = 3  # Print both: table and tree
0080   kernel.eventAction().adopt(prt)
0081
0082   # Configure I/O
0083   # evt_lcio = geant4.setupLCIOOutput('LcioOutput','Lhe_dip_sol_circ-higgs-bb')
0084   # evt_lcio.OutputLevel = Output.ERROR
0085
0086   evt_root = geant4.setupROOTOutput('RootOutput', 'Lhe_dip_sol_circ-higgs-bb')
0087   evt_root.OutputLevel = Output.INFO
0088
0089   gen = DDG4.GeneratorAction(kernel, "Geant4GeneratorActionInit/GenerationInit")
0090   kernel.generatorAction().adopt(gen)
0091
0092   """
0093   # First particle generator: e-  non-isotropic generation using Gun:
0094   gun = DDG4.GeneratorAction(kernel,"Geant4ParticleGenerator/Gun");
0095   gun.Particle = 'e-'
0096   gun.Energy = 60 * GeV
0097   gun.Multiplicity = 1
0098   gun.Position = (0.*mm,0.*mm,0.*mm)
0099   gun.Direction = (1.,0.,0.)
0100   gun.Mask = 2
0101   gun.enableUI()
0102   kernel.generatorAction().adopt(gun)
0103   # Install vertex smearing for this primary e-
0104   gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearE-");
0105   gen.Mask = 1
0106   gen.Sigma = (0*mm, 0*mm, 0*mm, 0*ns)
0107   kernel.generatorAction().adopt(gen)
0108   """
0109
0110   # VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
0111   """
0112   Generation of isotrope tracks of a given multiplicity with overlay:
0113   """
0114   # VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
0115   """
0116   # First particle generator: pi+
0117   gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropPi+");
0118   gen.Particle = 'pi+'
0119   gen.Energy = 200*GeV
0120   gen.Multiplicity = 1
0121   gen.Mask = 1
0122   kernel.generatorAction().adopt(gen)
0123   # Install vertex smearing for this interaction
0124   gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearPi+");
0125   gen.Mask = 1
0126   gen.Offset = (20*mm, 10*mm, 10*mm, 0*ns)
0127   gen.Sigma = (4*mm, 1*mm, 1*mm, 0*ns)
0128   kernel.generatorAction().adopt(gen)
0129
0130   # Second particle generator: e-
0131   gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropE-");
0132   gen.Particle = 'e-'
0133   gen.Energy = 60 * GeV
0134   gen.Multiplicity = 2
0135   gen.Mask = 2
0136   kernel.generatorAction().adopt(gen)
0137   # Install vertex smearing for this interaction
0138   gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearE-");
0139   gen.Mask = 2
0140   gen.Offset = (-20*mm, -10*mm, -10*mm, 0*ns)
0141   gen.Sigma = (12*mm, 8*mm, 8*mm, 0*ns)
0142   kernel.generatorAction().adopt(gen)
0143   """
0144   # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0145
0146   # VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
0147   """
0148   Generation of primary particles from LCIO input files
0149   """
0150   # VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
0151
0152   # First particle file reader
0153   gen = DDG4.GeneratorAction(kernel, "LCIOInputAction/LCIO1")
0154   # gen.Input = "LCIOStdHepReader|/afs/.cern.ch/project/lhec/"
0155   #             "software/aidasoft/DD4hep/DD4hep/files/NC_bb_tag_2_pythia_events.hep"
0156   # gen.Input = "LCIOStdHepReader|/opt/DD4hep/files/NC_bb_tag_2_pythia_events.hep"
0157   gen.Input = "LCIOStdHepReader|/opt/DD4hep/files/lhec_for_peter/tag_2_pythia_events.hep"
0158   # gen.Input = "LCIOStdHepReader|/opt/DD4hep/files/single-top-tag_1_pythia_events.hep"
0159
0160   # gen.Input = "Geant4EventReaderHepMC|/opt/DD4hep/files/ePb-q2-0-i.mc2"
0161   # gen.Input = "LCIOStdHepReader|/opt/DD4hep/files/single-top-tag_1_pythia_events.hep"
0162   # gen.Input = "LCIOStdHepReader|/opt/DD4hep/files/root.hep"
0163
0164   gen.OutputLevel = 2  # generator_output_level
0165   gen.MomentumScale = 1.0
0166   gen.Mask = 1
0167   gen.enableUI()
0168   kernel.generatorAction().adopt(gen)
0169
0170   # Install vertex smearing for this interaction
0171   gen = DDG4.GeneratorAction(kernel, "Geant4InteractionVertexSmear/Smear1")
0172   gen.OutputLevel = 4  # generator_output_level
0173   gen.Mask = 1
0174   gen.Offset = (-20 * g4units.mm, -10 * g4units.mm, -10 * g4units.mm, 0 * g4units.ns)
0175   gen.Sigma = (12 * g4units.mm, 8 * g4units.mm, 8 * g4units.mm, 0 * g4units.ns)
0176   gen.enableUI()
0177   kernel.generatorAction().adopt(gen)
0178
0179   # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0180
0181   # Merge all existing interaction records
0182   gen = DDG4.GeneratorAction(kernel, "Geant4InteractionMerger/InteractionMerger")
0183   gen.OutputLevel = 4  # generator_output_level
0184   gen.enableUI()
0185   kernel.generatorAction().adopt(gen)
0186
0187   # Finally generate Geant4 primaries
0188   gen = DDG4.GeneratorAction(kernel, "Geant4PrimaryHandler/PrimaryHandler")
0189   gen.OutputLevel = 4  # generator_output_level
0190   gen.enableUI()
0191   kernel.generatorAction().adopt(gen)
0192
0193   # And handle the simulation particles.
0194   part = DDG4.GeneratorAction(kernel, "Geant4ParticleHandler/ParticleHandler")
0195   kernel.generatorAction().adopt(part)
0196   # part.SaveProcesses = ['conv','Decay']
0197   part.SaveProcesses = ['Decay']
0198   part.MinimalKineticEnergy = 10 * g4units.MeV
0199   part.OutputLevel = 5  # generator_output_level
0200   part.enableUI()
0201
0202   user = DDG4.Action(kernel, "Geant4TCUserParticleHandler/UserParticleHandler")
0203   user.TrackingVolume_Zmax = DDG4.EcalEndcap_zmin_fwd
0204   user.TrackingVolume_Rmax = DDG4.EcalBarrel_rmin
0205   user.enableUI()
0206   part.adopt(user)
0207
0208   """
0209   rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader")
0210   rdr.zSpread = 0.0
0211   rdr.lorentzAngle = 0.0
0212   rdr.OutputLevel = DDG4.OutputLevel.INFO
0213   rdr.Input = "LcioEventReader|test.data"
0214   rdr.enableUI()
0215   kernel.generatorAction().adopt(rdr)
0216   """
0217
0218   # Setup global filters fur use in sensntive detectors
0219   f1 = DDG4.Filter(kernel, 'GeantinoRejectFilter/GeantinoRejector')
0220   kernel.registerGlobalFilter(f1)
0221   f2 = DDG4.Filter(kernel, 'ParticleRejectFilter/OpticalPhotonRejector')
0222   f2.particle = 'opticalphoton'
0223   kernel.registerGlobalFilter(f2)
0224   f3 = DDG4.Filter(kernel, 'ParticleSelectFilter/OpticalPhotonSelector')
0225   f3.particle = 'opticalphoton'
0226   kernel.registerGlobalFilter(f3)
0227
0228   f4 = DDG4.Filter(kernel, 'EnergyDepositMinimumCut')
0229   f4.Cut = 0.5 * g4units.MeV
0230   f4.enableUI()
0231   kernel.registerGlobalFilter(f4)
0232
0233   # First the tracking detectors
0234   seq, act = geant4.setupTracker('SiVertexBarrel')
0235   seq.adopt(f1)
0236   act.adopt(f1)
0237
0238   seq, act = geant4.setupTracker('SiTrackerBarrel')
0239   seq.adopt(f1)
0240   act.adopt(f1)
0241   seq, act = geant4.setupTracker('SiTrackerForward')
0242   seq.adopt(f1)
0243   act.adopt(f1)
0244   seq, act = geant4.setupTracker('SiTrackerBackward')
0245   seq.adopt(f1)
0246   act.adopt(f1)
0247
0248   # Now the calorimeters
0249   seq, act = geant4.setupCalorimeter('EcalBarrel')
0250   seq.adopt(f3)
0251   act.adopt(f3)
0252   seq.adopt(f4)
0253   act.adopt(f4)
0254
0255   seq, act = geant4.setupCalorimeter('EcalEndcap_fwd')
0256   seq.adopt(f3)
0257   act.adopt(f3)
0258   seq.adopt(f4)
0259   act.adopt(f4)
0260   seq, act = geant4.setupCalorimeter('EcalEndcap_bwd')
0261   seq.adopt(f3)
0262   act.adopt(f3)
0263   seq.adopt(f4)
0264   act.adopt(f4)
0265
0266   seq, act = geant4.setupCalorimeter('HcalBarrel')
0267   seq.adopt(f3)
0268   act.adopt(f3)
0269   seq.adopt(f4)
0270   act.adopt(f4)
0271   seq, act = geant4.setupCalorimeter('HcalEndcap_fwd')
0272   seq.adopt(f3)
0273   act.adopt(f3)
0274   seq.adopt(f4)
0275   act.adopt(f4)
0276   seq, act = geant4.setupCalorimeter('HcalEndcap_bwd')
0277   seq.adopt(f3)
0278   act.adopt(f3)
0279   seq.adopt(f4)
0280   act.adopt(f4)
0281
0282   seq, act = geant4.setupCalorimeter('HcalPlug_fwd')
0283   seq.adopt(f3)
0284   act.adopt(f3)
0285   seq.adopt(f4)
0286   act.adopt(f4)
0287   seq, act = geant4.setupCalorimeter('HcalPlug_bwd')
0288   seq.adopt(f3)
0289   act.adopt(f3)
0290   seq.adopt(f4)
0291   act.adopt(f4)
0292
0293   seq, act = geant4.setupCalorimeter('MuonBarrel')
0294   seq.adopt(f2)
0295   act.adopt(f2)
0296   seq, act = geant4.setupCalorimeter('MuonEndcap_fwd1')
0297   seq.adopt(f2)
0298   act.adopt(f2)
0299   seq, act = geant4.setupCalorimeter('MuonEndcap_fwd2')
0300   seq.adopt(f2)
0301   act.adopt(f2)
0302   seq, act = geant4.setupCalorimeter('MuonEndcap_bwd1')
0303   seq.adopt(f2)
0304   act.adopt(f2)
0305   seq, act = geant4.setupCalorimeter('MuonEndcap_bwd2')
0306   seq.adopt(f2)
0307   act.adopt(f2)
0308
0309   """
0310   scan = DDG4.SteppingAction(kernel,'Geant4MaterialScanner/MaterialScan')
0311   kernel.steppingAction().adopt(scan)
0312   """
0313
0314   # Now build the physics list:
0315   phys = geant4.setupPhysics('QGSP_BERT')
0316   ph = DDG4.PhysicsList(kernel, 'Geant4PhysicsList/Myphysics')
0317   ph.addParticleConstructor(str('G4Geantino'))
0318   ph.addParticleConstructor(str('G4BosonConstructor'))
0319   ph.addParticleConstructor(str('G4LeptonConstructor'))
0320   ph.addParticleProcess(str('e[+-]'), str('G4eMultipleScattering'), -1, 1, 1)
0321   ph.addPhysicsConstructor(str('G4OpticalPhysics'))
0322   ph.enableUI()
0323   phys.adopt(ph)
0324   phys.dump()
0325
0326   kernel.configure()
0327   kernel.initialize()
0328
0329   # DDG4.setPrintLevel(Output.DEBUG)
0330   kernel.run()
0331   logging.info('End of run. Terminating .......')
0332   kernel.terminate()
0333
0334
0335 if __name__ == "__main__":
0336   run()