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