bug fixed JOs for HZZ offshell ggF samples with dimension-six SMEFT operator cpG*ctp

507xxx/507599/MadGraphControl_Pythia8EvtGen_SMEFTNLO_ZZ.py

+from MadGraphControl.MadGraphUtils import *
+
+from MadGraphControl.MadGraphUtils import *
+from MadGraphControl import MadGraphUtils
+MadGraphUtils.MADGRAPH_CATCH_ERRORS=False
+import math
+
+
+import os, shutil, copy
+### get MC job-options filename
+FIRST_DIR = (os.environ['JOBOPTSEARCHPATH']).split(":")[0]
+jofiles = [f for f in os.listdir(FIRST_DIR) if (f.startswith('mc') and f.endswith('.py'))]
+
+#---------------------------------------------------------------------------
+# General Settings
+#  need to provide the following in the top JO:
+#    EFTorder: order of the calculation, 0 for SM only, 1 for EFT linear term, 2 for EFT quadratic term, 3 for all
+#    EFTop: EFT parameter name
+#    EFTval: value of the EFT parameter, 1 by default 
+#    dogg:  loop-induced or not
+#    doNLO: LO or NLO model
+#---------------------------------------------------------------------------
+mode=0 # 0: single core; 1: cluster; 2: multicore
+cluster_type="condor"
+cluster_queue=None
+if mode!=1:
+  cluster_type=None
+  cluster_queue=None
+
+gridpack_mode=False
+
+## safe factor applied to nevents, to account for the filter efficiency
+safefactor = 1.1
+nevents = runArgs.maxEvents*safefactor if runArgs.maxEvents>0 else safefactor*evgenConfig.nEventsPerJob
+nevents = int(nevents)
+
+## EFT parameters
+dim6_ops_LO = {
+  'dim6'    : ['cpG'],
+  'dim62f'  : ['ctp','cpt']
+}
+
+dim6_ops_NLO = {
+  'dim6'    : ['cpG'],
+  'dim62f'  : ['ctp','cpt']
+}
+
+# the block DIM64F2L not available for NLO, some parameters not available for NLO
+if doNLO:
+  dim6_ops=copy.deepcopy(dim6_ops_NLO)
+else:
+  dim6_ops=copy.deepcopy(dim6_ops_LO)
+
+# check user defined EFT operator names
+l_EFTop=[]
+l_EFTval=[]
+if EFTop:    
+  if "," in EFTop:
+    for op in EFTop.split(","): l_EFTop.append(op)
+  else: l_EFTop.append(EFTop)
+
+  for op in l_EFTop:
+    op_found = 0
+    for bl in dim6_ops.keys():
+      for aop in dim6_ops[bl]:
+        if aop.upper() == op.upper():
+          op_found = 1
+          break
+    if not op_found:
+      raise RuntimeError("EFT parameter %s not avaiable." % op)
+
+if EFTval:    
+  if "," in EFTval:
+    for val in EFTval.split(","): l_EFTval.append(val)
+  else: l_EFTval.append(EFTval)
+else:
+  # use default value 1
+  for ie, op in enumerate(l_EFTop): l_EFTval.append(1)
+
+#---------------------------------------------------------------------------
+# Process
+#---------------------------------------------------------------------------
+
+# model with restrictions: SM, LO, NLO
+model_name_SM = "SMEFTatNLO"
+model_name_LO = model_name_SM+"-cpG_ctp_cpt"
+model_name_NLO = model_name_SM+"-cpG_ctp_cpt"
+model_name = model_name_SM
+# production mode
+prod_str = "p p"
+if dogg:
+  prod_str = "g g"
+# lepton final states
+decay_str = ""
+if doZZ<=0:
+  if do2l2v:
+    decay_str = "l+ l- vl vl~ / l+ l- vl vl~ w+ w- g"
+  else:
+    decay_str = "l+ l- l+ l- / g"
+
+
+if EFTorder==0: ## SM
+  if doZZ==1: ## on-shell ZZ
+    proc_str = "generate %s > z z QCD=2 QED=2 NP=0" % (prod_str)
+  elif doZZ==2: ## HZZ
+    proc_str = "generate p p > h > z z QCD=2 QED=2 NP=0"
+  elif doZZ==3: ## no Higgs
+    proc_str = "generate %s > z z / h QCD=2 QED=2 NP=0" % (prod_str)
+  elif doZZ==-1: ## H4l
+    proc_str = "generate %s > h > %s QCD=2 QED=4 NP=0" % (prod_str, decay_str)
+  else:  ## off-shell ll
+    proc_str = "generate %s > %s QCD=2 QED=4 NP=0" % (prod_str, decay_str)
+else:
+  if EFTorder==1: ## linear
+    if doZZ==1: ## on-shell ZZ
+      proc_str = "generate %s > z z QCD=2 QED=2 NP=2 NP^2==2" % (prod_str)
+    elif doZZ==2: ## HZZ
+      proc_str = "generate %s > h > z z QCD=2 QED=2 NP=2 NP^2==2" % (prod_str)
+    elif doZZ==3: ## no Higgs
+      proc_str = "generate %s > z z / h QCD=2 QED=2 NP=2 NP^2==2" % (prod_str)
+    elif doZZ==-1: ## H4l
+      proc_str = "generate %s > h > %s QCD=2 QED=4 NP=2 NP^2==2" % (prod_str, decay_str)
+    else:
+      proc_str = "generate %s > %s QCD=2 QED=4 NP=2 NP^2==2" % (prod_str, decay_str)
+  elif EFTorder==2: ## quadratic
+    if doZZ==1: ## on-shell ZZ
+      proc_str = "generate %s > z z QCD=2 QED=2 NP=2 NP^2==4" % (prod_str)
+    elif doZZ==2: ## HZZ
+      proc_str = "generate %s > h > z z QCD=2 QED=2 NP=2 NP^2==4" % (prod_str)
+    elif doZZ==3: ## no Higgs
+      proc_str = "generate %s > z z / h QCD=2 QED=2 NP=2 NP^2==4" % (prod_str)
+    elif doZZ==-1: ## H4l
+      proc_str = "generate %s > h > %s QCD=2 QED=4 NP=2 NP^2==4" % (prod_str, decay_str)
+    else:
+      proc_str = "generate %s > %s QCD=2 QED=4 NP=2 NP^2==4" % (prod_str, decay_str)
+  elif EFTorder==3: ## All
+    if doZZ==1: ## on-shell ZZ
+      proc_str = "generate %s > z z QCD=2 QED=2 NP=2" % (prod_str)
+    elif doZZ==2: ## HZZ
+      proc_str = "generate %s > h > z z QCD=2 QED=2 NP=2" % (prod_str)
+    elif doZZ==3: ## no Higgs
+      proc_str = "generate %s > z z / h QCD=2 QED=2 NP=2" % (prod_str)
+    elif doZZ==-1: ## H4l
+      proc_str = "generate %s > h > %s QCD=2 QED=4 NP=2" % (prod_str, decay_str)
+    else:
+      proc_str = "generate %s > %s QCD=2 QED=4 NP=2" % (prod_str, decay_str)
+
+  if doNLO:
+    model_name = model_name_NLO
+  else:
+    model_name = model_name_LO
+
+# QCD order: None, "[QCD]", "[noborn=QCD]"
+if QCDmode:
+  proc_str += " "+QCDmode
+
+# use model files locally, to fix cpG
+#mgmodels='/cvmfs/atlas.cern.ch/repo/sw/Generators/madgraph/models/latest/'
+mgmodels='/cvmfs/atlas.cern.ch/repo/sw/Generators/madgraph/models/latest'
+model_name_LO_org = mgmodels+"/" + model_name_SM+"-cpG_ctp_cpt"
+model_name_NLO_org = mgmodels+"/" + model_name_SM+"-cpG_ctp_cpt"
+model_name_custom = mgmodels+"/"+model_name_SM+"-cpG_ctp_cpt"
+
+restricted_model=model_name
+model_path='/cvmfs/atlas.cern.ch/repo/sw/Generators/madgraph/models/latest'
+if fixModel:
+  ## copy model files to local
+  model_name_def=model_name.split('-')[0]
+  local_dir='mgmodels_local'
+  if os.path.exists(local_dir):
+    shutil.rmtree(local_dir)
+  os.mkdir(local_dir)
+  local_model=local_dir+'/'+model_name_def
+  shutil.copytree(mgmodels+model_name_def,local_model)
+
+  ## fix cpG
+  tmp_cards=['restrict_LO.dat', 'restrict_NLO.dat']
+  for tmp_card in tmp_cards:
+    rst_card=local_model+'/'+tmp_card
+    if os.access(rst_card,os.R_OK):
+      rst_card_old = local_model+"/"+tmp_card.replace(".dat", "Full.dat")
+      mglog.info('Modifying %s to turn off cpG' % tmp_card)
+      os.rename(rst_card, rst_card_old)
+
+      oldCard = open(rst_card_old, 'r')
+      newCard = open(rst_card, 'w')
+      for line in iter(oldCard):
+        if "cpG" in line: line=line.replace('0.8e-04', '0.0e-04')
+        newCard.write(line.strip()+'\n')
+      oldCard.close()
+      newCard.close()
+
+  #model_path="./%s" % local_dir
+  model_path="/cvmfs/atlas.cern.ch/repo/sw/Generators/madgraph/models/latest"
+  restricted_model="%s/%s" % (model_path, model_name)
+  ## update model names which might be used later for reweighting
+  model_path=os.getcwd()+"/" + local_dir
+  model_name_LO = model_path+"/" + model_name_SM+"-LO"
+  model_name_NLO = model_path+"/" + model_name_SM+"-NLO"
+
+  model_name_LO_org = model_path+"/" + model_name_SM+"-LOFull"
+  model_name_NLO_org = model_path+"/" + model_name_SM+"-NLOFull"
+  
+  
+  ## remove cpG
+  if 'dim6' in dim6_ops and 'cpG' in dim6_ops['dim6']:
+    dim6_ops['dim6'].remove('cpG')
+
+restricted_model="%s/%s" % (model_path, model_name)
+
+if not is_gen_from_gridpack():
+    process = """
+    set max_npoint_for_channel 4
+    import model %s
+    define p = g u c d s b u~ c~ d~ s~ b~
+    define j = g u c d s b u~ c~ d~ s~ b~
+    define l+ = e+ mu+
+    define l- = e- mu-
+    %s
+    output -f
+    """ % (restricted_model, proc_str)
+    
+    process_dir = new_process(process)
+else:
+    process_dir = MADGRAPH_GRIDPACK_LOCATION
+
+
+beamEnergy=-999
+if hasattr(runArgs,'ecmEnergy'):
+  beamEnergy = runArgs.ecmEnergy / 2.
+else: 
+  raise RuntimeError("No center of mass energy found.")
+
+#----------------------------------------------------------------------------
+# Random Seed
+#----------------------------------------------------------------------------
+randomSeed = 0
+if hasattr(runArgs,'randomSeed'): randomSeed = runArgs.randomSeed
+
+#---------------------------------------------------------------------------
+# Number of Events
+#---------------------------------------------------------------------------
+skip_events=0
+if hasattr(runArgs,'skipEvents'): skip_events=runArgs.skipEvents
+
+#---------------------------------------------------------------------------
+# MG5 run Card
+#---------------------------------------------------------------------------
+
+extras = {
+    'asrwgtflavor':"5",
+    'lhe_version':"3.0",
+    'ptj':"0",
+    'ptb':"0",
+    'pta':"0",
+    'ptl':"3",
+    'etaj':"-1",
+    'etab':"-1",
+    'etaa':"-1",
+    'etal':"3.0",
+    'drjj':"0",
+    'mmll':"0",
+    'draa':"0",
+    'draj':"0",
+    'use_syst':"False",
+    'scalefact':'1.0',
+    'drjl':"0",
+    'dral':"0",
+    'drll':"0.04",
+    'maxjetflavor':"5" ,
+    'cut_decays'  :'T',
+    'auto_ptj_mjj': 'F',
+    'nevents'     : nevents,
+}
+## m4l cut for 4l
+if (doZZ==0 or doZZ==-1) and (not do2l2v):
+  extras['mmnl'] = "130"
+## settings for NLO
+# call MadGraphUtils.is_NLO_run to make sure it's really an NLO run
+if is_NLO_run(process_dir=process_dir):
+  extras = {
+    'lhe_version':"3.0",
+    'parton_shower' :'PYTHIA8',
+    'ptl':"3",
+    'etal':"3.0",
+    'ptj':"0",
+    'etaj':"-1",
+    'mll':"0",
+    'mll_sf':"0",
+    'drll':"0.04",
+    'drll_sf':"0.04",
+    'jetalgo':'-1',
+    'jetradius' : '0.4',
+    'maxjetflavor':"5" ,
+    'nevents'     : nevents,
+  }
+
+modify_run_card(process_dir=process_dir,runArgs=runArgs,settings=extras)
+
+
+# Set up batch jobs
+if mode==1 and (not is_gen_from_gridpack()):
+  modify_config_card(process_dir=process_dir,settings={'cluster_type':cluster_type,'cluster_queue':cluster_queue})
+else:
+  modify_config_card(process_dir=process_dir)
+
+
+#---------------------------------------------------------------------------
+# MG5 param Card
+#---------------------------------------------------------------------------
+if EFTorder>0: 
+    # Set SMEFT@NLO parameters
+    ## params is a dictionary of dictionaries (each dictionary is a separate block)
+    params = dict() 
+    
+    # initial setting, disable all EFT parameters
+    for bl in dim6_ops.keys():
+      params[bl] = dict()
+      for op in dim6_ops[bl]:
+        params[bl][op] = '1.000000e-20'
+
+    params['dim6']['Lambda'] = '1000.'
+    # set the EFT operator under test
+    for ie, op in enumerate(l_EFTop):
+      val = l_EFTval[ie]
+      for bl in dim6_ops.keys():
+        if op in dim6_ops[bl]:
+          params[bl][op] = val
+    
+    modify_param_card(process_dir=process_dir,params=params)
+
+
+# Build reweight_card.dat
+if doReweight:
+  # Create reweighting card
+  reweight_card_loc=process_dir+'/Cards/reweight_card.dat'
+  rwcard = open(reweight_card_loc,'w')
+
+  ## parameters to reweight
+  # format is [ 'paramName1_paraValue1', ...] for 1D reweighting
+  # or [ 'paramName1_paraValue1--paramName2_paraValue2', ...] for 2D reweighting, etc
+  dict_para_rwt={}
+  for rw_name in reweights:
+    dict_para_rwt[rw_name]=[]
+    for param in rw_name.split('--'):
+      param_name, value = param.split('_')
+      dict_para_rwt[rw_name].append([param_name, value])
+
+  ## model
+  curr_model_name=model_name ## keep track of the model used for the previous reweighting
+
+  # !!! NLO model has to be used for interference (either tree*loop or loop*loop)
+  dim6_ops=copy.deepcopy(dim6_ops_NLO)
+
+  ## linear
+  if doReweight==1 or doReweight==3:
+    if doZZ==1: ## on-shell ZZ
+      proc_str = "%s > z z QCD=2 QED=2 NP=2 NP^2==2" % (prod_str)
+    elif doZZ==2: ## HZZ
+      proc_str = "%s > h > z z QCD=2 QED=2 NP=2 NP^2==2" % (prod_str)
+    elif doZZ==3: ## no Higgs
+      proc_str = "%s > z z / h QCD=2 QED=2 NP=2 NP^2==2" % (prod_str)
+    elif doZZ==-1: ## H4l
+      proc_str = "%s > %s QCD=2 QED=4 NP=2 NP^2==2" % (prod_str, decay_str)
+    else:
+      proc_str = "%s > %s QCD=2 QED=4 NP=2 NP^2==2" % (prod_str, decay_str)
+    ## launch reweighting
+    for rw_name in dict_para_rwt.keys():
+      # check if the reweighting parameter is included in the model or not
+      skipop=1
+      hasCPG=False ## special treatment for cpG (tree level)
+      for param in dict_para_rwt[rw_name]:
+        [param_name, value] = param
+        if param_name.upper() == "CPG": hasCPG=True
+        for bl in dim6_ops.keys():
+          if param_name in dim6_ops[bl]:
+            skipop=0
+            break
+      if skipop and (not hasCPG): continue 
+
+      # change process
+      ## cpG: tree*loop
+      if hasCPG:
+        # locally fixed model has cpG removed, need to use the original one
+        # only update the model if it's different than the nominal one
+        if curr_model_name != model_name_NLO_org: 
+
+          curr_model_name=model_name_NLO_org
+
+        rwcard.write("""
+change process %s [virt=QCD]
+""" % proc_str)
+
+      ## others: loop*loop
+      else:
+        # use locally fixed model with cpG removed
+        # only update the model if it's different than the nominal one
+        if curr_model_name != model_name_NLO: 
+          rwcard.write("""
+change model %s
+""" % model_name_NLO)
+          curr_model_name=model_name_NLO
+
+        rwcard.write("""
+change process %s %s
+""" % (proc_str, QCDmode))
+
+
+      # launch reweighting
+      rwcard.write("launch --rwgt_name=%s_lin\n" % rw_name)
+
+      params = dict() 
+      # initial setting, disable all EFT parameters
+      for bl in dim6_ops.keys():
+        params[bl] = dict()
+        for op in dim6_ops[bl]:
+          params[bl][op] = '1.000000e-20'
+    
+      params['dim6']['Lambda'] = '1000.'
+      # set the EFT operator under test
+      for param in dict_para_rwt[rw_name]:
+        [param_name, value] = param
+        for bl in dim6_ops.keys():
+          if param_name in dim6_ops[bl] or (bl=="dim6" and param_name.upper() == "CPG"):
+            params[bl][param_name] = value
+    
+      rwcard.write("set dim6 8 1.0\n")
+      rwcard.write("set dim62f 15 1.0e-20\n")
+      rwcard.write("set dim62f 19 1.0e-20")
+
+  ## quadratic
+  if doReweight==2 or doReweight==3:
+    dim6_ops=copy.deepcopy(dim6_ops_NLO)
+
+    if doZZ==1: ## on-shell ZZ
+      proc_str = "%s > z z QCD=2 QED=2 NP=2 NP^2==4" % (prod_str)
+    elif doZZ==2: ## HZZ
+      proc_str = "%s > h > z z QCD=2 QED=2 NP=2 NP^2==4" % (prod_str)
+    elif doZZ==3: ## no Higgs
+      proc_str = "%s > z z / h QCD=2 QED=2 NP=2 NP^2==4" % (prod_str)
+    elif doZZ==-1: ## H4l
+      proc_str = "%s > h > %s QCD=2 QED=4 NP=2 NP^2==4" % (prod_str, decay_str)
+    else:
+      proc_str = "%s > %s QCD=2 QED=4 NP=2 NP^2==4" % (prod_str, decay_str)
+    ## launch reweighting
+    for rw_name in dict_para_rwt.keys():
+      # check if the reweighting parameter is included in the model or not
+      skipop=1
+      hasCPG=False ## special treatment for cpG (tree level)
+      hasCPGonly=True  ## for cpG*cXX (tree*loop)
+      for param in dict_para_rwt[rw_name]:
+        [param_name, value] = param
+        if param_name.upper() == "CPG": hasCPG=True
+        else: hasCPGonly=False
+        for bl in dim6_ops.keys():
+          if param_name in dim6_ops[bl]:
+            skipop=0
+            break
+      if skipop and (not hasCPG): continue 
+
+      # change process
+      ## cpG: tree*tree
+      if hasCPG and hasCPGonly:
+        # TODO: this does not work, better to run it directly, not via reweighting
+        continue
+        # LO or NLO
+        ## NLO
+        if doNLO:
+          rwcard.write("""
+change mode NLO
+""")
+          # locally fixed model has cpG removed, need to use the original one
+          # only update the model if it's different than the nominal one
+          dim6_ops=copy.deepcopy(dim6_ops_NLO)
+          if curr_model_name != model_name_NLO_org: 
+            rwcard.write("""
+change model %s
+""" % model_name_NLO_org)
+            curr_model_name=model_name_NLO_org
+        
+          # change the process accordingly for cpG*cpG at NLO
+          if doZZ==1: ## on-shell ZZ
+            tmp_proc_str = "p p > z z QCD=2 QED=2 NP=2 NP^2==4 [QCD]"
+          elif doZZ==2: ## HZZ
+            tmp_proc_str = "p p > h > z z QCD=2 QED=2 NP=2 NP^2==4 [QCD]"
+          elif doZZ==3: ## no Higgs
+            tmp_proc_str = "p p > z z / h QCD=2 QED=2 NP=2 NP^2==4 [QCD]"
+          elif doZZ==-1: ## H4l
+            tmp_proc_str = "g g > h > %s QCD=2 QED=4 NP=2 NP^2==4 [QCD]" % (decay_str)
+          else:
+            tmp_proc_str = "g g > %s QCD=2 QED=4 NP=2 NP^2==4 [QCD]" % (decay_str)
+        
+          rwcard.write("""
+change process %s
+""" % tmp_proc_str)
+
+        ## LO
+        else:
+          dim6_ops=copy.deepcopy(dim6_ops_LO)
+          if curr_model_name != model_name_LO_org:
+            rwcard.write("""
+change model %s
+""" % model_name_LO_org)
+            curr_model_name=model_name_LO_org
+        
+          rwcard.write("""
+change process %s
+""" % proc_str)
+
+      ## cpG*cXX: tree*loop
+      elif hasCPG:
+        # NLO model has to be used here
+        # locally fixed model has cpG removed, need to use the original one
+        # only update the model if it's different than the nominal one
+        dim6_ops=copy.deepcopy(dim6_ops_NLO)
+        if curr_model_name != model_name_NLO_org: 
+          rwcard.write("""
+change model %s
+""" % model_name_NLO_org)
+          curr_model_name=model_name_NLO_org
+      
+        rwcard.write("""
+change process %s [virt=QCD]
+""" % (proc_str))
+
+      ## others: loop*loop
+      else:
+        # NLO model has to be used here
+        dim6_ops=copy.deepcopy(dim6_ops_NLO)
+        if curr_model_name != model_name_NLO: 
+          rwcard.write("""
+change model %s
+""" % model_name_NLO)
+          curr_model_name=model_name_NLO
+      
+        rwcard.write("""
+change process %s %s
+""" % (proc_str, QCDmode))
+
+      # launch reweighting
+      rwcard.write("launch --rwgt_name=%s_quad\n" % rw_name)
+    
+      params = dict() 
+      # initial setting, disable all EFT parameters
+      for bl in dim6_ops.keys():
+        params[bl] = dict()
+        for op in dim6_ops[bl]:
+          params[bl][op] = '1.000000e-20'
+    
+      params['dim6']['Lambda'] = '1000.'
+      # set the EFT operator under test
+      for param in dict_para_rwt[rw_name]:
+        [param_name, value] = param
+        for bl in dim6_ops.keys():
+          if param_name in dim6_ops[bl] or (bl=="dim6" and param_name.upper() == "CPG"):
+            params[bl][param_name] = value
+    
+      rwcard.write("set dim6 8 1.0\n") #cpG
+      rwcard.write("set dim62f 15 1.0e-20\n") #cpt 
+      rwcard.write("set dim62f 19 1.0") #ctp
+
+  rwcard.close()
+
+print_cards()
+
+#---------------------------------------------------------------------------
+# MG5 process (lhe) generation
+#---------------------------------------------------------------------------
+
+generate(process_dir=process_dir,runArgs=runArgs,grid_pack=gridpack_mode)
+
+arrange_output(process_dir=process_dir,runArgs=runArgs,lhe_version=3,saveProcDir=True)
+
+# Helper for resetting process number
+check_reset_proc_number(opts)
+
+
+#---------------------------------------------------------------------------
+# Shower
+#---------------------------------------------------------------------------
+
+include("Pythia8_i/Pythia8_A14_NNPDF23LO_EvtGen_Common.py")
+include("Pythia8_i/Pythia8_aMcAtNlo.py")
+
+## Z boson decay