Generating Higgs To 4 Muons at NIKHEF
Revision as of 15:42, 20 April 2005 by Ivov@nikhef.nl (talk | contribs)
This page contains a description on how to generate Higgs production at the LHC, where the Higgs boson decays into 2 Z bosons that each decay into either 2 electrons or 2 muons. It is ment as a starting point. In this example we will use AtlFast so we see what the events look like in Atlas. The next steop will be the full simulation and reconstruction. We will produce a combined ntuple (CBNT) that contains the MC truth and reconstructed AtlFast objects.
Note: We assume you have the CMT and Athena set-up at NIKHEF in ordnung Starting with CMT and Athena at NIKHEF
1. producing 10 events
- a) Go to your favorite area (your project for example) and create a running directory
- cd /project/atlas/users/<your_login_name>
- mkdir MyGeneration
- b) Create your joboptions file
- Each athena job requires a joboptions file as input. Here we will create a joboption file that will:
- Define what 'algorithms to run (in our case Pythia and Atlfast)
- Define the Pythia settings
- Define output parameters/ntuples
- Open a file called joboptions_HiggsGeneration.py and fill it with the following code
###############################################################
#
# Job options file
#
#==============================================================
#--------------------------------------------------------------
# General Application Configuration options
#--------------------------------------------------------------
include("AthenaCommon/Atlas_Gen.UnixStandardJob.py")
theApp.setup( MONTECARLO )
include( "PartPropSvc/PartPropSvc.py" )
include( "AtlfastStandardOptions.py" )
#--------------------------------------------------------------
# Private Application Configuration options
#--------------------------------------------------------------
theApp.DLLs += [ "TruthExamples", "Pythia_i" ]
theApp.Dlls += [ "GaudiAlg", "GaudiAud" ]
theApp.Dlls += [ "AtlfastAlgs" ]
# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL )
MessageSvc = Service( "MessageSvc" )
MessageSvc.OutputLevel = 3
MessageSvc.defaultLimit = 9999999
#--------------------------------------------------------------
# Create the top level algorithm (and the Python equivalent)
#--------------------------------------------------------------
theApp.TopAlg = [ "Sequencer/TopSequencer" ]
TopSequencer = Algorithm( "TopSequencer" )
TopSequencer.StopOverride = TRUE
# Create two paths and populate them
TopSequencer.Members = [ "Sequencer/Generator", "Sequencer/Atlfast" ]
Generator = Algorithm( "Generator" )
Generator.Members = [ "Pythia", ]
#--------------------------------------------------------------
# Event related parameters
#--------------------------------------------------------------
# Number of events to be processed (default is 10)
theApp.EvtMax = 10
#theApp.EvtMax = XNEVENTSX
#--------------------------------------------------------------
# Random numbers
#--------------------------------------------------------------
AtRndmGenSvc = Service( "AtRndmGenSvc" )
AtRndmGenSvc.Seeds = ["PYTHIA 5769791 690419913", "PYTHIA_INIT 690501 4106941"];
#AtRndmGenSvc.Seeds = ["PYTHIA XRNDPYTHIA0X XRNDPYTHIA1X", "PYTHIA_INIT XRNDPYTHIA2X XRNDPYTHIA3X"]
#--------------------------------------------------------------
# Algorithms Private Options (Pythia)
#--------------------------------------------------------------
Pythia = Algorithm( "Pythia" )
#-- Stolen from DC1 production (simple test)
Pythia.PythiaCommand = [
# Higgs production
"pysubs msel 16", # Higgs production
"pydat2 pmas 25 1 150.", # Higgs mass
# Higgs decay
"pydat3 mdme 210 1 0", # H -> dd
"pydat3 mdme 211 1 0", # H -> uu
"pydat3 mdme 212 1 0", # H -> ss
"pydat3 mdme 213 1 0", # H -> cc
"pydat3 mdme 214 1 0", # H -> bb
"pydat3 mdme 215 1 0", # H -> tt
"pydat3 mdme 218 1 0", # H -> e+e-
"pydat3 mdme 219 1 0", # H -> mu+mu-
"pydat3 mdme 220 1 0", # H -> tau+tau-
"pydat3 mdme 222 1 0", # H -> gluon gluon
"pydat3 mdme 223 1 0", # H -> gamma gamma
"pydat3 mdme 224 1 0", # H -> gluon + gamma
"pydat3 mdme 225 1 1", # H -> ZZ (ON)
"pydat3 mdme 226 1 0", # H -> WW
# Z decay
"pydat3 mdme 174 1 0", # Z -> dd
"pydat3 mdme 175 1 0", # Z -> uu
"pydat3 mdme 176 1 0", # Z -> ss
"pydat3 mdme 177 1 0", # Z -> cc
"pydat3 mdme 178 1 0", # Z -> bb
"pydat3 mdme 179 1 0", # Z -> tt
"pydat3 mdme 182 1 1", # Z -> e+e- (ON)
"pydat3 mdme 183 1 0", # Z -> nu_e nu_e
"pydat3 mdme 184 1 1", # Z -> mu+mu- (ON)
"pydat3 mdme 185 1 0", # Z -> nu_mu nu_mu
"pydat3 mdme 186 1 0", # Z -> tau+tau-
"pydat3 mdme 187 1 0" # Z -> nu_tau nu_tau
]
#--------------------------------------------------------------
# What do you want in the CBNT
#--------------------------------------------------------------
include( "CBNT_Athena/CBNT_Athena_jobOptions.py" )
include( "CBNT_Athena/CBNT_EventInfo_jobOptions.py" )
include( "RecExCommon/CBNT_Truth_jobOptions.py" )
CBNT_Truth.MaxNbParticles = 6000 # maximum number of particles in the ntuple
CBNT_Truth.MaxNbVertices = 6000 # maximum number of vertices in the ntuple
CBNT_Athena.NtupleLocID = "/FILE1/CBNT/t3333" # name of the Tree
All.Enable = True # save ALL particles
#--------------------------------------------------------------
# (Root) output service
#--------------------------------------------------------------
# ROOT Output Parameters
theApp.Dlls += [ "RootHistCnv" ]
# Select ROOT persistency
theApp.HistogramPersistency = "ROOT"
#--------------------------------------------------------------
# Output file for the Ntuple
#--------------------------------------------------------------
NTupleSvc = Service( "NTupleSvc" )
NTupleSvc.Output = [ "FILE1 DATAFILE='./HiggsNtuple.root' OPT='NEW'" ]
#NTupleSvc.Output = [ "FILE1 DATAFILE='./XOutputDirCBNTX/HiggsNtuple.JobXJOBNUMBERX.root' OPT='NEW'" ]
#==============================================================
#
# End of job options file
#
###############################################################
- b) Get additional steering files
- There are some requires input files that you have to get. The way to obbtain them is using the command get_files.
- get_files PDGTABLE.MeV
- get_files PartPropSvc.py
- get_files AtlfastStandardOptions.py
- c) Run Athena
- Running Athena is now a single line: athena.py joboptions_HiggsGeneration.py
- You can ask athena to print out al commands it is processing:
- Full output of what athena is doing: athena.py -bs joboptions_HiggsGeneration.py
- d) Check the output
- In the directory now an output file called HiggsNtuple.root has been produced. When you look in the file you'll find the Tree for the CBNT (MC Truth) and that for AtlFast (Event in ATLAS). Now you only have to write a ROOT macro to read the file and do your analysis.
- Finished!
2. producing 10,000 events in 10 sets of 1000
When using the joboptions file in the example above you have produced 10 events, but when you want to have a bit large production you want to automatise everything a bit:
Create a new joboption file for each job each having
- A unique random number seed for Pythia
- A user defined number of events
- An output ntuple that is different for each event
- Store output and Logfiles in a separate directory
To do just this a small script has been created.
- a) Create a BASICS-directory
- Create a directory with the basic information you require.
mkdir HiggsGen_BASICS cp PDGTABLE.MeV ./HiggsGen_BASICS/ cp PartPropSvc.py ./HiggsGen_BASICS/ cp AtlfastStandardOptions.py ./HiggsGen_BASICS/
- We also copy the joboptions file there and rename it
cp joboptions_HiggsGeneration.py ./HiggsGen_BASICS/joboptions_HiggsGeneration.py.BASIC
- b) Edit the standard joboptions file
To allow the script to change the job-dependent settings in the joboptions file we'll now have to change 3 lines in the joboptions file.
- a) change Number of events
theApp.EvtMax = 10
- now becomes
theApp.EvtMax = XNEVENTSX
- b) Change Random number seeds for Pythia
AtRndmGenSvc.Seeds = ["PYTHIA 5769791 690419913", "PYTHIA_INIT 690501 4106941"];
- now becomes
AtRndmGenSvc.Seeds = ["PYTHIA XRNDPYTHIA0X XRNDPYTHIA1X", "PYTHIA_INIT XRNDPYTHIA2X XRNDPYTHIA3X"]
- c) Change Name output file
<NTupleSvc.Output = [ "FILE1 DATAFILE='./HiggsNtuple.root' OPT='NEW'" ]
- now becomes
NTupleSvc.Output = [ "FILE1 DATAFILE='./XOutputDirCBNTX/HiggsNtuple.JobXJOBNUMBERX.root' OPT='NEW'/font>