ATLAS Offline Software
PowhegControl_ttHplus_NLO.py
Go to the documentation of this file.
1 #--------------------------------------------------------------
2 # EVGEN configuration
3 #--------------------------------------------------------------
4 evgenConfig.description = 'POWHEG+Madspin+Pythia8 ttbar production with Powheg hdamp equal 1.5*top mass, A14 tune, ME NNPDF30 NLO, A14 NNPDF23 LO, Higgs plus decays for masses from 80-160'
5 evgenConfig.keywords =['Higgs','MSSM','BSMHiggs','chargedHiggs','top','ttbar']
6 evgenConfig.contact = [ 'anna.ivina@cern.ch']
7 evgenConfig.minevents = 10000
8 
9 
10 include('PowhegControl/PowhegControl_tt_Common.py')
11 # Initial settings
12 if hasattr(PowhegConfig, "topdecaymode"):
13  # Use PowhegControl-00-02-XY (and earlier) syntax
14  PowhegConfig.topdecaymode = 22222 # inclusive top decays
15 else:
16  # Use PowhegControl-00-03-XY (and later) syntax
17  PowhegConfig.decay_mode = "t t~ > undecayed"
18  PowhegConfig.MadSpin_enabled = False
19 
20 PowhegConfig.hdamp = 258.75 # 1.5 * mtop
21 DoSingleWeight = False
22 if DoSingleWeight:
23  PowhegConfig.mu_F = 1.0
24  PowhegConfig.mu_R = 1.0
25  PowhegConfig.PDF = 260000
26 else:
27  PowhegConfig.mu_F = [1.0, 2.0, 0.5, 1.0, 1.0, 0.5, 2.0, 0.5, 2.0] # List of factorisation scales which pairs with renormalisation scale below
28  PowhegConfig.mu_R = [1.0, 1.0, 1.0, 2.0, 0.5, 0.5, 2.0, 2.0, 0.5] # List of renormalisation scales
29  PowhegConfig.PDF = [260000, 25200, 13165, 90900] # NNPDF30, MMHT, CT14, PDF4LHC - PDF variations with nominal scale variation
30  PowhegConfig.PDF.extend(range(260001, 260101)) # Include the NNPDF error set
31  PowhegConfig.PDF.extend(range(90901 , 90931 )) # Include the PDF4LHC error set
32 
33 #PowhegConfig.nEvents *= 1.1 # compensate filter efficiency, this will probably is not needed for Hplus
34 PowhegConfig.generate()
35 
36 #--------------------------------------------------------------
37 # Now preparing for MadSpin
38 #--------------------------------------------------------------
39 # Get the DSID and also make the strings
40 thisDSID = runArgs.runNumber #you can use it if you want to decay tops or W to something else
41 model = '2HDMtypeII'
42 madspin_card_rep = 'madspin_card.dat'
43 madspin_in = 'import run.lhe'
44 madspin_rep = 'set ms_dir MadSpin'
45 madspin_seed = runArgs.randomSeed
46 
47 #open the mscard and write the parameters, set seed,and decay
48 #included into the define p and j the b quark
49 mscard = open(madspin_card_rep,'w')
50 mscard.write("""
51 set Nevents_for_max_weigth 250 # number of events for the estimate of the max. weight (default: 75)
52 set max_weight_ps_point 1000 # number of PS to estimate the maximum for each event (default: 400)
53 set seed %i
54 %s
55 %s
56 define l+ = e+ mu+ ta+
57 define l- = e- mu- ta-
58 define vl = ve vm vt
59 define vl~ = ve~ vm~ vt~
60 define p = g u c d b s u~ c~ d~ s~ b~
61 define j = g u c d b s u~ c~ d~ s~ b~
62 \n
63 """%(madspin_seed,madspin_in,madspin_rep))
64 
65 
66 #make the mscard
67 #W decay
68 wstr = ', w- > l- vl~'
69 
70 t2str = 'decay t~ > b~ w-'
71 t1str = 'decay t > b h+'
72 mscard.write("""%s\n%s%s\nlaunch"""%(t1str,t2str,wstr))
73 mscard.close()
74 
75 #--------------------------------------------------------------
76 # Charge Higgs, and all other masses in GeV
77 #--------------------------------------------------------------
78 
79 mhc_str = str(runArgs.jobConfig[0]) # JOB OPTION NAME MUST CONTAIN THE MASS WE WANT TO SIMULATE IN FORMAT LIKE: *_H160_*
80 mhc=0
81 int(s)
82 for s in mhc_str.split("_"):
83  ss=s.replace("H","")
84  if ss.isdigit():
85  mhc = int(ss)
86  print "Charged Higgs mass set to %i"%mhc
87 if mhc==0:
88  raise RuntimeError("Charged Higgs mass not set, mhc=0, check joOption name %s"%mhc_str)
89 
90 import math
91 mh1=1.250e+02
92 mh2=math.sqrt(math.pow(mhc,2)+math.pow(8.0399e+01,2))
93 mh3=mh2
94 
95 
96 masses = {'25':str(mh1)+' # mh1',
97  '35':str(mh2)+' # mh2',
98  '36':str(mh2)+' # mh2',
99  '37':str(mhc)+' # mhc'}
100 decayss={'37': 'DECAY 37 1.300000e+02 # whc'}
101 
102 
103 #write_param_card is part of the 2HDMTypeII model
104 import shutil,os,subprocess
105 
106 paramFileName = 'MadGraph_param_card_ttHplus_NLO.dat'
107 paramFileNameN = 'param_card.dat'
108 
109 from MadGraphControl.MadGraphUtils import *
110 build_param_card(param_card_old=paramFileName,param_card_new=paramFileNameN,masses=masses,decays=decayss)
111 
112 
113 
114 fin = open('PowhegOTF._1.events','r')
115 fout = open('run.lhe','w')
116 line = fin.readline()
117 while line != "-->\n":
118  fout.write(line)
119  line = fin.readline()
120 fout.write(line)
121 
122 
123 # add the process
124 #to use 5 FS you need to set b quark mass to zero at model level, use restrict_nobmass provided in the model folder
125 fout.write('<MG5ProcCard>\n')
126 fout.write('import model /cvmfs/atlas.cern.ch/repo/sw/Generators/madgraph/models/latest/%s-nobmass\n'%(model))
127 fout.write('generate p p > t t~ [QCD] \n')
128 fout.write('</MG5ProcCard>\n')
129 
130 # add run parameters
131 eline = str(PowhegConfig.nEvents)+' = nevents\n'
132 fout.write('<MGRunCard>\n')
133 fout.write('#0.01 = req_acc_FO\n') #Need because of a new check introduced in 2.5.0 : if not there, RunCardLO will be considered, which wants a line like that "x = nhel ! x = 0,1"
134 fout.write(eline)
135 fout.write('50.0 = bwcutoff\n')
136 fout.write('</MGRunCard>\n')
137 
138 
139 
140 #add model parameters
141 fout.write('<slha>\n')
142 shutil.copyfileobj(open(paramFileNameN, 'r'), fout)
143 fout.write('</slha>\n')
144 fout.write('<montecarlomasses>\n')
145 fout.write(' 1 0.330000E+00\n')
146 fout.write(' 2 0.330000E+00\n')
147 fout.write(' 3 0.500000E+00\n')
148 fout.write(' 4 0.1550000E+01\n')
149 fout.write(' 5 0.470000E+01\n')
150 fout.write(' 11 0.510999E-03\n')
151 fout.write(' 13 0.105658E+00\n')
152 fout.write(' 15 0.177682E+01\n')
153 fout.write(' 21 0.000000E+00\n')
154 fout.write('</montecarlomasses>\n')
155 
156 # add the events !
157 line = fin.readline()
158 while 'LesHouchesEvents' not in line:
159  fout.write(line)
160  line = fin.readline()
161 fout.write(line)
162 fout.close()
163 
164 
165 # run MadSpin
166 os.system('$MADPATH/MadSpin/madspin < madspin_card.dat')
167 
168 # This is a Powheg job, so expect lhe file to be called PowhegOTF._1.events
169 unzip = subprocess.Popen(['gunzip','run_decayed.lhe.gz'])
170 unzip.wait()
171 os.system('cp run_decayed.lhe PowhegOTF._1.events')
172 
173 # Does MadSpin add some spurious lines that make the weight names unreadable ?
174 # see https://bugs.launchpad.net/mg5amcnlo/+bug/1720979
175 with open('PowhegOTF._1.events', 'r+') as f:
176  t = f.read()
177  to_delete = [ '<![CDATA[', ']]>' ]
178  f.seek(0)
179  for line in t.split('\n'):
180  if line not in to_delete:
181  f.write(line + '\n')
182  f.truncate()
183 
184 #--------------------------------------------------------------
185 # Pythia8 showering
186 #--------------------------------------------------------------
187 include('Pythia8_i/Pythia8_A14_NNPDF23LO_EvtGen_Common.py')
188 include("Pythia8_i/Pythia8_Powheg_Main31.py")
189 genSeq.Pythia8.Commands += [ 'POWHEG:pThard = 0' ]
190 genSeq.Pythia8.Commands += [ 'POWHEG:nFinal = 2' ]
191 genSeq.Pythia8.Commands += [ 'POWHEG:pTdef = 2' ]
192 genSeq.Pythia8.Commands += [ 'POWHEG:vetoCount = 3' ]
193 genSeq.Pythia8.Commands += [ 'POWHEG:pTemt = 0' ]
194 genSeq.Pythia8.Commands += [ 'POWHEG:emitted = 0' ]
195 genSeq.Pythia8.Commands += [ 'POWHEG:MPIveto = 0' ]
MadGraphUtils
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
python.Include.include
include
Definition: Include.py:318
Trk::open
@ open
Definition: BinningType.h:40
str
Definition: BTagTrackIpAccessor.cxx:11