diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..cf46103 --- /dev/null +++ b/.gitignore @@ -0,0 +1,36 @@ +# Prerequisites +*.d + +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll + +# Fortran module files +*.mod +*.smod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +# Other +*.pdf +*~ \ No newline at end of file diff --git a/MG5_aMC_v3_1_1/bin/mg5_aMC b/MG5_aMC_v3_1_1/bin/mg5_aMC new file mode 100755 index 0000000..634c0fc --- /dev/null +++ b/MG5_aMC_v3_1_1/bin/mg5_aMC @@ -0,0 +1,227 @@ +#! /usr/bin/env python3 +from __future__ import absolute_import +from __future__ import print_function +import time +start = time.time() + +################################################################################ +# +# Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors +# +# This file is a part of the MadGraph5_aMC@NLO project, an application which +# automatically generates Feynman diagrams and matrix elements for arbitrary +# high-energy processes in the Standard Model and beyond. +# +# It is subject to the MadGraph5_aMC@NLO license which should accompany this +# distribution. +# +# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch +# +################################################################################ +"""This is the main executable, a simple frontend to set up the PYTHONPATH +and call immediately the command line interface scripts""" + +import sys +if sys.version_info[1] < 7: + if sys.version_info[0] == 2: + sys.exit("MadGraph5_aMC@NLO works only with python 2.7 or python 3.7 (and later).\n"+\ + " You are currently using Python2.%s. Please use a more recent version of Python." % sys.version_info[1]) + if sys.version_info[0] == 3: + sys.exit("MadGraph5_aMC@NLO works only with python 2.7 or python 3.7 (and later).\n"+\ + " You are currently using Python 3.%i. So please upgrade your version of Python." % sys.version_info[1] +\ + " If you have python2.7 installed you need to run the code as\n"+\ + " python27 ./bin/mg5_aMC \n") + +try: + import six +except ImportError: + message = 'madgraph requires the six module. The easiest way to install it is to run "pip%s install six --user"\n' % (sys.version_info[0] if sys.version_info[0]==3 else '') + message += 'in case of problem with pip, you can download the file at https://pypi.org/project/six/ . It has a single python file that you just need to put inside a directory of your $PYTHONPATH environment variable.' + sys.exit(message) + +import os +import optparse + +# Get the parent directory (mg root) of the script real path (bin) +# and add it to the current PYTHONPATH + +root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0] +sys.path.insert(0, root_path) + + +# Write out nice usage message if called with -h or --help +usage = "usage: %prog [options] [FILE] " +parser = optparse.OptionParser(usage=usage) +parser.add_option("-l", "--logging", default='INFO', + help="logging level (DEBUG|INFO|WARNING|ERROR|CRITICAL) [%default]") +parser.add_option("-f", "--file", default='', + help="Use script file FILE") +parser.add_option("-d", "--mgme_dir", default='', dest = 'mgme_dir', + help="Use MG_ME directory MGME_DIR") +parser.add_option("","--web", action="store_true", default=False, dest='web', \ + help='force to be in secure mode') +parser.add_option("","--debug", action="store_true", default=False, dest='debug', \ + help='force to launch debug mode') +parser.add_option("-m", "--mode", dest="plugin", + help="Define some additional command provide by a PLUGIN") +parser.add_option("-s", "--nocaffeinate", action="store_false", default=True, dest='nosleep', + help='For mac user, forbids to use caffeinate when running with a script') +(options, args) = parser.parse_args() +if len(args) == 0: + args = '' + +import subprocess + +# Check if optimize mode is (and should be) activated +if __debug__ and not options.debug and \ + (not os.path.exists(os.path.join(root_path, 'bin','create_release.py')) or options.web): + os.system('%s -O -W ignore::DeprecationWarning %s' % (sys.executable, ' '.join(sys.argv))) + sys.exit() + +import logging +import logging.config +import madgraph.interface.coloring_logging + +if sys.version_info[0] ==2: + logging.warning("\033[91mpython2 support will be removed in last quarter 2021. If you use python2 due to issue with Python3, please report them on https://bugs.launchpad.net/mg5amcnlo\033[0m") + + +if ' ' in os.getcwd(): + logging.warning("\033[91mPath does contains spaces. We advise that you change your current path to avoid to have space in the path.\033[0m") + +try: + import readline +except ImportError: + try: + import pyreadline as readline + except: + print( "For tab completion and history, install module readline.") +else: + import rlcompleter + + if 'r261:67515' in sys.version and 'GCC 4.2.1 (Apple Inc. build 5646)' in sys.version: + readline.parse_and_bind("bind ^I rl_complete") + readline.__doc__ = 'libedit' + + elif hasattr(readline, '__doc__'): + if 'libedit' not in readline.__doc__: + readline.parse_and_bind("tab: complete") + else: + readline.parse_and_bind("bind ^I rl_complete") + else: + readline.__doc__ = 'GNU' + readline.parse_and_bind("tab: complete") + + # charge history file + try: + history_file = os.path.join(os.environ['HOME'], '.mg5', 'mg5history') + readline.read_history_file(history_file) + except: + pass + +try: + import psyco + psyco.full() +except: + pass + +if __debug__: + print( 'Running MG5 in debug mode') + +# Set logging level according to the logging level given by options +#logging.basicConfig(level=vars(logging)[options.logging]) +try: + if __debug__ and options.logging == 'INFO': + options.logging = 'DEBUG' + if options.logging.isdigit(): + level = int(options.logging) + else: + level = eval('logging.' + options.logging) + + logging.config.fileConfig(os.path.join(root_path, 'madgraph', 'interface', '.mg5_logging.conf')) + logging.root.setLevel(level) + logging.getLogger('madgraph').setLevel(level) + logging.getLogger('madevent').setLevel(level) +except: + pass + +import madgraph.interface.master_interface as interface + +if __debug__ and time.time() - start > 0.5: + print(( 'WARNING: loading of madgraph too slow!!!', time.time() - start)) + +try: + import madgraph.various.misc as misc + misc.EasterEgg('loading') +except: + pass + + +if options.plugin: + if not os.path.exists(os.path.join(root_path, 'PLUGIN', options.plugin)): + print( "ERROR: %s is not present in the PLUGIN directory. Please install it first") + __import__('PLUGIN.%s' % options.plugin) + plugin = sys.modules['PLUGIN.%s' % options.plugin] + if not plugin.new_interface: + logging.warning("Plugin: %s do not define dedicated interface and should be used without the --mode options" % options.plugin) + sys.exit() + import madgraph.various.misc as misc + if not misc.is_plugin_supported(plugin): + sys.exit() + cmd_line = plugin.new_interface(mgme_dir = options.mgme_dir) + cmd_line.plugin=options.plugin +elif options.web: + cmd_line = interface.MasterCmdWeb() +else: + cmd_line = interface.MasterCmd(mgme_dir = options.mgme_dir) + + + +# Call the cmd interface main loop +try: + if options.file or args: + if sys.platform == "darwin" and options.nosleep: + logging.getLogger('madgraph').warning("launching caffeinate to prevent idle sleep when MG5aMC is running. Run './bin/mg5_aMC -s' to prevent this.") + pid = os.getpid() + subprocess.Popen(['caffeinate', '-i', '-w', str(pid)]) + # They are an input file + if args: + input_file = os.path.realpath(args[0]) + else: + input_file = os.path.realpath(options.file) + if options.web: + cmd_line.debug_output = os.path.join(os.path.dirname(input_file),'generation.log') + cmd_line.use_rawinput = False + cmd_line.haspiping = False + cmd_line.run_cmd('import ' + input_file) + cmd_line.run_cmd('quit') + else: + cmd_line.use_rawinput = False + cmd_line.haspiping = False + cmd_line.run_cmd('import ' + input_file) + cmd_line.run_cmd('quit') + else: + # Interactive mode + if options.web: + if 'MADGRAPH_DATA' not in os.environ: + os.environ['MADGRAPH_DATA'] = root_path + os.environ['MADGRAPH_BASE'] = os.path.join(root_path,'input') + os.environ['REMOTE_USER'] = 'webmode' + cmd_line.cmdloop() + else: + cmd_line.cmdloop() +except KeyboardInterrupt: + print('writting history and quit on KeyboardInterrupt') + pass + +try: + cmd_line.exec_cmd('quit all', printcmd=False) + readline.set_history_length(100) + if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')): + os.mkdir(os.path.join(os.environ['HOME'], '.mg5')) + readline.write_history_file(history_file) +except Exception as error: + pass + + + diff --git a/README.md b/README.md index bf8b2ef..7bbbe09 100644 --- a/README.md +++ b/README.md @@ -1 +1,9 @@ -# aQGC \ No newline at end of file +# aQGC + +How to run: + +``` +lsetup "root 6.20.06-x86_64-centos7-gcc8-opt" +git clone -b cleanUp https://github.com/jstupak/aQGC.git +./aQGC/test.sh +``` diff --git a/analyzeDelphes.py b/analyzeDelphes.py index 87e4f32..f9848ae 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -1,11 +1,45 @@ maxEvents=9E9 -DEBUG=True +DEBUG=False + +jetType='VLCjetR15N2' + +particle={ + 11:"e", + 13:"mu", + 22:"gamma", #using this for all jets + 21:"j", + 23:"z", + 24:"w", + 12:"nu" +} + +#OBJECT SELECTION +pTMin={ + 11:1, + 13:1, + 21:5, #using this for all jets + 22:1 +} + +etaMax={11:2.5, + 13:2.5, + 21:10, + 22:1 + } +#--------------------------------------------------------- +#EVENT SELECTION + +#ADD EVENT SELECTION CRITERIA HERE + +#--------------------------------------------------------- import pdb +import operator from sys import argv #--------------------------------------------------------- from ROOT import * -from ROOT import TLorentzVector +TH1.SetDefaultSumw2() +import math gSystem.Load("libDelphes") gStyle.SetOptStat(0) @@ -14,7 +48,7 @@ def printHist(h): for i in range(h.GetNbinsX()+2): print h.GetBinContent(i), - + #--------------------------------------------------------- #return TLorentzVector corresponding to sum of inputs @@ -34,11 +68,11 @@ def getParents(p): result=[p] motherIndices=[] - if p.M1!=-1 and tree.Particle[p.M1].PID==p.PID: + if p.M1!=-1 and event.Particle[p.M1].PID==p.PID: motherIndices.append(p.M1) - if p.M2!=-1 and tree.Particle[p.M2].PID==p.PID: + if p.M2!=-1 and event.Particle[p.M2].PID==p.PID: motherIndices.append(p.M2) - result+=[getParents(tree.Particle[i]) for i in motherIndices] + result+=[getParents(event.Particle[i]) for i in motherIndices] return result @@ -48,297 +82,423 @@ def isBeamRemnant(p): return parents.Status==4 #--------------------------------------------------------- +def sortFunc(x): + return (x.PT)*math.cosh(x.Eta) +#--------------------------------------------------------- if __name__=='__main__': - - #f=TFile('/Users/jstupak/ou/Snowmass2021/workArea/aQGC/MG5_aMC_v2_9_2/PROC_sm_10/Events/run_01/tag_1_delphes_events.root') - #f=TFile('/raid01/users/kawale/MG5_aMC_v2_7_2/delphes/delphes.root') - #f=TFile('/raid01/users/azartash/delphes/mumumumuww_vbseft10.root ') - print str(argv[0]) - print str(argv[1]) - print str(argv[2]) - print argv - f=TFile(str(argv[1])) - #f=TFile('~/tmp/mumumumuww_schaneft10.root') - #path = str(argv[1]) - #prefix = "/raid01/users/cwaits/MG5_aMC_v2_7_2/"+str(argv[2]) - #suffix = "/Events/run_01/delphes.root" - #name = path[len(prefix):] - #name = name[:-len(suffix)] - - - tree=f.Get("Delphes") - name=str(argv[2]) - output=TFile(name+".delphes.root","RECREATE") - + + inputName=str(argv[1]) + if len(argv)>2: + tag=str(argv[2]) + else: + tag='' + outputName=inputName.replace('.root',(bool(tag)*('.%s'%tag))+'.hist.root') + + print "inputName:",inputName + print "tag:",tag + + f=TFile(inputName) + t=f.Delphes + output=TFile(outputName,"RECREATE") + + t.GetEntry(0) + sqrtS=0 + sqrtS+=t.Particle[0].E + sqrtS+=t.Particle[1].E #sets upper limit for bin range - bin_range = 15000 + bin_range = sqrtS/2 #declares truth-level pT, p, eta, and multiplicity histograms for e,mu,W,Z,gamma - T_e_pT = TH1F('T_e_Pt','Truth Electron pT;pT (GeV);Events', 200, 0, bin_range) - T_e_p = TH1F('T_e_p', 'Truth Electron Momentum;p (GeV);Events', 200, 0, bin_range) - T_e_eta = TH1F('T_e_eta', 'Truth Electron Eta;Eta;Events', 20, -4, 4) - T_e_multiplicity = TH1F('T_e_multiplicity', 'Truth Electron Multiplicity;Multiplicity;Events', 7, -.5, 6.5) - - T_mu_pT = TH1F('T_mu_Pt', 'Truth Muon pT;pT (GeV);Events', 200, 0, bin_range) - T_mu_p = TH1F('T_mu_p', 'Truth Muon Momentum;p (GeV);Events', 200, 0, bin_range) - T_mu_eta = TH1F('T_mu_eta', 'Truth Muon Eta;Eta;Events', 20, -4, 4) - T_mu_multiplicity = TH1F('T_mu_multiplicity', 'Truth Muon Multiplicity;Multiplicity;Events', 7, -.5, 6.5) - - T_W_pT = TH1F('T_W_Pt', 'Truth W pT;pT (GeV);Events', 200, 0, bin_range) - T_W_p = TH1F('T_W_p', 'Truth W Momentum;p (GeV);Events', 200, 0, bin_range) - T_W_eta = TH1F('T_W_eta', 'Truth W Eta;Eta;Events', 20, -4, 4) - T_W_multiplicity = TH1F('T_W_multiplicity', 'Truth W Multiplicity;Multiplicity;Events', 7, -.5, 6.5) - - T_z_pT = TH1F('T_z_Pt', 'Truth Z pT;pT (GeV);Events', 200, 0, bin_range) - T_z_p = TH1F('T_z_p', 'Truth Z Momentum;p (GeV);Events', 200, 0, bin_range) - T_z_eta = TH1F('T_z_eta', 'Truth Z Eta;Eta;Events', 20, -4, 4) - T_z_multiplicity = TH1F('T_z_multiplicity', 'Truth Z Multiplicity;Multiplicity;Events', 7, -.5, 6.5) - - T_gamma_pT = TH1F('T_gamma_Pt', 'Truth Photon pT;pT (GeV);Events', 200, 0, bin_range) - T_gamma_p = TH1F('T_gamma_p', 'Truth Photon Momentum;p (GeV);Events', 200, 0, bin_range) - T_gamma_eta = TH1F('T_gamma_eta', 'Truth Photon Eta;Eta;Events', 20, -4, 4) - T_gamma_multiplicity = TH1F('T_gamma_multiplicity', 'Truth Photon Mulitiplcity;Mulitplicity;Events', 7, -.5, 6.5) - + h_truth={} + for p in [11,13,21,22,23,24,12]: + h_truth[p]={} + h_truth[p]['pT'] =TH1F('T_%s_pT'%particle[p],';Truth %s pT [GeV];Events'%particle[p], 10, 0, bin_range) + h_truth[p]['p'] =TH1F('T_%s_p'%particle[p],';Truth %s p [GeV];Events'%particle[p], 10, 0, bin_range) + h_truth[p]['eta']=TH1F('T_%s_eta'%particle[p],';Truth %s #eta;Events'%particle[p], 20, -4, 4) + h_truth[p]['mult']=TH1F('T_%s_mult'%particle[p],';Truth %s multiplicity;Events'%particle[p], 7, -0.5, 6.5) + + #invariant mass of WW pair histogram + T_WW_mass = TH1F('WW_mass', ';Mass [GeV];Events', 200, 0, sqrtS) + #Eta vs m(WW) 2-D histogram + T_eta_vs_WWmass = TH2F('T_eta_vs_WWmass', ';Beam Remnant Eta; WW mass [GeV]', 10, 0, 10, 10, 500, sqrtS) #missing energy histograms - T_missingEt = TH1F('T_missingEt', ';Missing Transverse Energy (GeV);Events', 200, 0, bin_range) - T_missingE = TH1F('T_missingE', ';Missing Energy (GeV);Events', 200, 0, bin_range) - - #declares Reco-level inclusive, leading, 2nd leading, and 3rd leanding histograms for e,mu,gamma - R_e_pT_I = TH1F('R_e_Pt_I', ';pT (GeV);Events', 200, 0, bin_range) - R_e_p_I = TH1F('R_e_p_I', ';p (GeV);Events', 200, 0, bin_range) - R_e_eta_I = TH1F('R_e_eta_I', ';Eta;Events', 20, -4, 4) - R_e_pT_L = TH1F('R_e_Pt_L', ';pT (GeV);Events', 200, 0, bin_range) - R_e_eta_L = TH1F('R_e_eta_L', ';Eta;Events', 20, -4, 4) - R_e_pT_2 = TH1F('R_e_Pt_2', ';pT (GeV);Events', 200, 0, bin_range) - R_e_eta_2 = TH1F('R_e_eta_2', ';Eta;Events', 20, -4, 4) - R_e_pT_3 = TH1F('R_e_Pt_3', ';pT (GeV);Events', 200, 0, bin_range) - R_e_eta_3 = TH1F('R_e_eta_3', ';Eta;Events', 20, -4, 4) - R_e_multiplicity = TH1F('R_e_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) - - R_mu_pT_I = TH1F('R_mu_Pt_I', ';pT (GeV);Events', 200, 0, bin_range) - R_mu_p_I = TH1F('R_mu_p_I', ';p (GeV);Events', 200, 0, bin_range) - R_mu_eta_I = TH1F('R_mu_eta_I', ';Eta;Events', 20, -10, 10) - R_mu_pT_L = TH1F('R_mu_Pt_L', ';pT (GeV);Events', 200, 0, bin_range) - R_mu_eta_L = TH1F('R_mu_eta_L', ';Eta;Events', 20, -4, 4) - R_mu_pT_2 = TH1F('R_mu_Pt_2', ';pT (GeV);Events', 200, 0, bin_range) - R_mu_eta_2 = TH1F('R_mu_eta_2', ';Eta;Events', 20, -4, 4) - R_mu_pT_3 = TH1F('R_mu_Pt_3', ';pT (GeV);Events', 200, 0, bin_range) - R_mu_eta_3 = TH1F('R_mu_eta_3', ';Eta;Events', 20, -4, 4) - R_mu_multiplicity = TH1F('R_mu_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) - - R_gamma_pT_I = TH1F('R_gamma_Pt_I', ';pT (GeV);Events', 200, 0, bin_range) - R_gamma_p_I = TH1F('R_gamma_p_I', ';p (GeV);Events', 200, 0, bin_range) - R_gamma_E_I = TH1F('R_gamma_E_I', ';Energy (GeV);Events', 200, 0, bin_range) - R_gamma_eta_I = TH1F('R_gamma_eta_I', ';Eta;Events', 20, -4, 4) - R_gamma_pT_L = TH1F('R_gamma_Pt_L', ';pT (GeV);Events', 200, 0, bin_range) - R_gamma_E_L = TH1F('R_gamma__E_L', ';Energy (GeV);Events', 200, 0, bin_range) - R_gamma_eta_L = TH1F('R_gamma_eta_L', ';Eta;Events', 20, -4, 4) - R_gamma_pT_2 = TH1F('R_gamma_Pt_2', ';pT (GeV);Events', 200, 0, bin_range) - R_gamma_eta_2 = TH1F('R_gamma_eta_2', ';Eta;Events', 20, -4, 4) - R_gamma_pT_3 = TH1F('R_gamma_Pt_3', ';pT (GeV);Events', 200, 0, bin_range) - R_gamma_eta_3 = TH1F('R_gamma_eta_3', ';Eta;Events', 20, -4, 4) - R_gamma_multiplicity = TH1F('R_gamma_multiplicity', ';Multiplicity;Events', 5, -.5, 4.5) - + T_missingEt = TH1F('T_missingEt', ';Missing Transverse Energy [GeV];Events', 200, 0, bin_range) + T_missingE = TH1F('T_missingE', ';Missing Energy [GeV];Events', 200, 0, bin_range) + + h_reco={} + for p in [11,13,21,22]: + h_reco[p]={} + h_reco[p]['pT']={} + h_reco[p]['p']={} + h_reco[p]['eta']={} + h_reco[p]['mass']={} + h_reco[p]['mult']=TH1F('R_%s_mult'%particle[p],';%s multiplicity;Events'%particle[p], 7, -0.5, 6.5) + for i in ['I']+range(4): + h_reco[p]['pT'][i]=TH1F('R_%s_Pt_%s'%(particle[p],i),';%s pT [GeV];Events'%particle[p], 200, 0, bin_range) + h_reco[p]['p'][i]=TH1F('R_%s_P_%s'%(particle[p],i),';%s P [GeV];Events'%particle[p], 20, 0, bin_range) + h_reco[p]['eta'][i]=TH1F('R_%s_eta_%s'%(particle[p],i),';%s #eta;Events'%particle[p], 20, -4, 4) + if (p == 21): + h_reco[p]['mass'][i]=TH1F('R_%s_mass_%s'%(particle[p],i),';%s mass [GeV];Events'%particle[p], 20, 0, bin_range/2) + + R_reco_muon_leading_cos = TH1F('R_muon_leading_cos', ';cos(theta);Events', 10, -1, 1) + R_reco_muon_sub_cos = TH1F('R_muon_sub_cos', ';cos(theta);Events', 10, -1, 1) #histograms for OS pairs - R_ee_pT = TH1F('ee_pT', ';pT (GeV);Events', 200, 0, bin_range) + R_ee_pT = TH1F('ee_pT', ';pT [GeV];Events', 200, 0, bin_range) R_ee_eta = TH1F('ee_Eta', ';Eta;Events', 20, -4, 4) - R_ee_mass = TH1F('ee_mass', ';Mass (GeV);Events', 200, 0, 200) + R_ee_mass = TH1F('ee_mass', ';Mass [GeV];Events', 200, 0, 200) R_ee_deltaEta = TH1F('ee_Delta Eta', ';Delta Eta;Events', 20, -4, 4) R_ee_multiplicity = TH1F('multiplicity', ';Multiplicity;Events', 5, -.5, 4.5) - R_mumu_pT = TH1F('mumu_pT', ';pT (GeV);Events', 200, 0, bin_range) + R_mumu_pT = TH1F('mumu_pT', ';pT [GeV];Events', 200, 0, bin_range) R_mumu_eta = TH1F('mumu_Eta', ';Eta;Events', 20, -4, 4) - R_mumu_mass = TH1F('mumu_mass', ';Mass (GeV);Events', 200, 0, 200) + R_mumu_mass = TH1F('mumu_mass', ';Mass [GeV];Events', 20, 0, sqrtS) R_mumu_deltaEta = TH1F('mumu_Delta Eta', ';Delta Eta;Events', 20, -4, 4) R_mumu_multiplicity = TH1F('mumu_multiplicity', ';Multiplicity;Events', 5, -.5, 4.5) + R_mumu_p = TH1F('mumu_P', ';P [GeV];Events', 20, 0, sqrtS) + R_mumu_cos = TH1F('mumu_cos', ';cos(theta);Events', 10, -1, 1) - R_emu_pT = TH1F('emu_pT', ';pT (GeV);Events', 200, 0, bin_range) + R_emu_pT = TH1F('emu_pT', ';pT [GeV];Events', 200, 0, bin_range) R_emu_eta = TH1F('emu_Eta', ';Eta;Events', 20, -4, 4) - R_emu_mass = TH1F('emu_mass', ';Mass (GeV);Events', 200, 0, 200) + R_emu_mass = TH1F('emu_mass', ';Mass [GeV];Events', 200, 0, 200) R_emu_deltaEta = TH1F('emu_Delta Eta', ';Delta Eta;Events', 20, -4, 4) R_emu_multiplicity = TH1F('emu_multplicity', ';Multiplicity;Events', 5, -.5, 4.5) #histograms for missing energy - R_missingET = TH1F('R_MissingET', ';Missing Transverse Energy (GeV);Events', 200, 0, bin_range) - R_missingE = TH1F('R_MissingE', ';Missing Energy (GeV);Events', 200, 0, bin_range) - R_missingMass = TH1F('R_MissingMass', ';Missing Mass (GeV);Events' , 200, 0, bin_range) + R_missingET = TH1F('R_MissingET', ';Missing Transverse Energy [GeV];Events', 200, 0, bin_range) + R_missingE = TH1F('R_MissingE', ';Missing Energy [GeV];Events', 200, 0, bin_range) + R_missingMass = TH1F('R_MissingMass', ';Missing Mass [GeV];Events' , 200, 0, bin_range) + + #histograms for beam remnants + R_Wjets_pT = TH1F('R_Wjets_Pt', ';pT [GeV];Events', 20, 0, bin_range) + R_Wjets_p = TH1F('R_Wjets_p', ';p [GeV];Events', 200, 0, bin_range) + R_Wjets_eta = TH1F('R_Wjets_eta', ';Eta;Events', 40, -10, 10) + R_Wjets_multiplicity = TH1F('R_Wjets_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) + R_Wjets_mass = TH1F('R_Wjets_mass', ';mass [GeV];Events', 20, 0, sqrtS) + R_Wjets_p_leading = TH1F('R_Wjets_p_leading', ';p [GeV];Events', 20, 0, bin_range) + R_Wjets_p_sub = TH1F('R_Wjets_p_sub', ';p [GeV];Events', 20, 0, bin_range) + R_Wjets_cos_leading = TH1F('R_Wjets_cos_leading', ';cos(theta);Events', 10, -1, 1) + R_Wjets_cos_sub = TH1F('R_Wjets_cos_sub', ';cos(theta);Events', 10, -1, 1) + R_WJets_p_pair = TH1F('R_Wjets_p_pairs', ';p [GeV];Events', 20, 0, sqrtS) + R_Wjets_cos_pair = TH1F('R_Wjets_cos_pairs', ';cos(theta);Events',10, -1, 1) #histograms for beam remnants - T_beamRemnants_pT = TH1F('T_beamRemnants_Pt', ';pT (GeV);Events', 200, 0, bin_range) - T_beamRemnants_p = TH1F('T_beamRemnants_p', ';p (GeV);Events', 200, 0, bin_range) - T_beamRemnants_eta = TH1F('T_beamRemnants_eta', ';Eta;Events', 20, -10, 10) + T_beamRemnants_pT = TH1F('T_beamRemnants_Pt', ';pT [GeV];Events', 20, 0, bin_range) + T_beamRemnants_p = TH1F('T_beamRemnants_p', ';p [GeV];Events', 200, 0, bin_range) + T_beamRemnants_eta = TH1F('T_beamRemnants_eta', ';Eta;Events', 40, -10, 10) T_beamRemnants_multiplicity = TH1F('T_beamRemnants_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) - Test_beamRemnants_pT = TH1F('Test_beamRemnants_Pt', ';pT (GeV);Events', 30, 0, 3) - - for event in range(min(tree.GetEntries(),maxEvents)): - tree.GetEntry(event) - - #truth-level - T_electrons = 0 - T_muons = 0 - W = 0 - z = 0 - T_photons = 0 - neutrino_Px = 0 - neutrino_Py = 0 - neutrino_Pz = 0 - beamRemnants = 0 + T_beamRemnants_cosh = TH1F('T_beamRemnants_cosh', ';Cosh(Eta)/pT;Events', 20, 0, math.cosh(10)/(10*bin_range)) + T_beamRemnants_mass = TH1F('T_beamRemnants_mass', ';mass [GeV];Events', 300, 0, bin_range) + + #histograms for non-beam remnants + T_nonBeamRemnants_pT = TH1F('T_nonBeamRemnants_Pt', ';pT [GeV];Events', 20, 0, bin_range) + T_nonBeamRemnants_p = TH1F('T_nonBeamRemnants_p', ';p [GeV];Events', 200, 0, bin_range) + T_nonBeamRemnants_eta = TH1F('T_nonBeamRemnants_eta', ';Eta;Events', 40, -10, 10) + T_nonBeamRemnants_multiplicity = TH1F('T_nonBeamRemnants_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) + T_nonBeamRemnants_cosh = TH1F('T_nonBeamRemnants_cosh', ';Cosh(Eta)/pT;Events', 20, 0, math.cosh(10)/(10*bin_range)) + + #histograms for reco beam remnants + R_beamRemnants_pT = TH1F('R_beamRemnants_Pt', ';pT [GeV];Events', 20, 100, bin_range) + R_beamRemnants_p = TH1F('R_beamRemnants_p', ';p [GeV];Events', 200, 100, bin_range) + R_beamRemnants_eta = TH1F('R_beamRemnants_eta', ';Eta;Events', 40, -2.5, 2.5) + R_beamRemnants_multiplicity = TH1F('R_beamRemnants_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) + R_beamRemnants_cosh = TH1F('R_beamRemnants_cosh', ';Cosh(Eta)/pT;Events', 20, 0, math.cosh(10)/(10*bin_range)) + R_beamRemnants_mass = TH1F('R_beamRemnants_mass', ';mass [GeV];Events', 600, 0, sqrtS) + + #histograms for reco non-beam remnants + R_nonBeamRemnants_pT = TH1F('R_nonBeamRemnants_Pt', ';pT [GeV];Events', 20, 0, 100) + R_nonBeamRemnants_p = TH1F('R_nonBeamRemnants_p', ';p [GeV];Events', 200, 0, 100) + R_nonBeamRemnants_eta = TH1F('R_nonBeamRemnants_eta', ';Eta;Events', 40, -2.5, 2.5) + R_nonBeamRemnants_multiplicity = TH1F('R_nonBeamRemnants_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) + R_nonBeamRemnants_cosh = TH1F('R_nonBeamRemnants_cosh', ';Cosh(Eta)/pT;Events', 20, 0, math.cosh(10)/(10*bin_range)) + + CutFlow = TH1F('CutFlow', ';Cut;Events', 5, .5, 5.5) + + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + m=0 + for event in f.Delphes: + m+=1 + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + #truth level + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + truthElectrons=selector(event.Particle,'x.Status==1 and abs(x.PID)==11') + truthMuons =selector(event.Particle,'x.Status==1 and abs(x.PID)==13') + truthWs =selector(event.Particle,'abs(x.Status)==22 and abs(x.PID)==24') + truthZs =selector(event.Particle,'abs(x.Status) in range(21, 30) and abs(x.PID)==23') + truthPhotons =selector(event.Particle,'abs(x.Status)==1 and abs(x.PID)==22') + truthNeutrinos=selector(event.Particle,'x.Status==1 and (abs(x.PID)==12 or abs(x.PID)==14 or abs(x.PID)==16)') + + beamRemnantMuons =selector(truthMuons,'isBeamRemnant(x)') + nonBeamRemnantMuons=selector(truthMuons,'not isBeamRemnant(x)') + truthLeptons=truthElectrons + truthMuons + + electrons=selector(event.Electron, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[11],etaMax[11])) + muons =selector(event.Muon, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[13],etaMax[13])) + photons =selector(event.Photon, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[22],etaMax[22])) + jets =selector(event.__getattr__(jetType), 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[21],etaMax[21])) + leptons = electrons + muons + + # requires either one or two BRM and W-jets + reco_beamRemnants=[] + reco_nonBeamRemnants=[] + W_jets=[] + for i in range(len(muons)): + if (muons[i].P4().P() >= 100): + reco_beamRemnants.append(muons[i]) + for i in range(len(jets)): + if (50 <= jets[i].P4().M() <= 100): + W_jets.append(jets[i]) + + # sort object lists + for collection in [electrons, muons, photons, jets, W_jets, reco_beamRemnants]: collection.sort(key=sortFunc, reverse=True) + + # increments cutFlow histogram + # for zmass or zmass-veto channels + if len(muons)>=2: + muon_4vec=muons[0].P4() + muons[1].P4() + CutFlow.Fill(1) + if (len(muons) == 2) and (abs(muon_4vec.M()-91)<15): + CutFlow.Fill(2) + if len(jets) == 2: + CutFlow.Fill(3) + if len(W_jets) >= 1: + CutFlow.Fill(4) + if len(W_jets) == 2: + CutFlow.Fill(5) + + if len(muons) != 2: + continue + if len(W_jets) != 2: + continue + # zmass veto channel + #if (muon_4vec.M()-91)<15: + # continue + # zmass channel + if abs(muon_4vec.M()-91)>15: + continue + R_mumu_mass.Fill(muon_4vec.M()) + + h_truth[11]['mult'].Fill(len(truthElectrons)) + h_truth[13]['mult'].Fill(len(truthMuons)) + h_truth[24]['mult'].Fill(len(truthWs)) + h_truth[23]['mult'].Fill(len(truthZs)) + h_truth[22]['mult'].Fill(len(truthPhotons)) + h_truth[12]['mult'].Fill(len(truthNeutrinos)) + T_beamRemnants_multiplicity.Fill(len(beamRemnantMuons)) + T_nonBeamRemnants_multiplicity.Fill(len(nonBeamRemnantMuons)) + + if (len(truthWs) >= 2): + WW_4vec=truthWs[0].P4() + truthWs[1].P4() + T_WW_mass.Fill(WW_4vec.M()) + if (len(beamRemnantMuons) >= 2) and (len(truthWs) >= 2): + T_eta_vs_WWmass.Fill(abs(beamRemnantMuons[0].Eta), WW_4vec.M()) + T_eta_vs_WWmass.Fill(abs(beamRemnantMuons[1].Eta), WW_4vec.M()) + + for i in range(len(truthElectrons)): + h_truth[11]['pT'].Fill(truthElectrons[i].PT) + h_truth[11]['p'].Fill(truthElectrons[i].P4().P()) + h_truth[11]['eta'].Fill(truthElectrons[i].Eta) + + for i in range(len(truthMuons)): + h_truth[13]['pT'].Fill(truthMuons[i].PT) + h_truth[13]['p'].Fill(truthMuons[i].P4().P()) + h_truth[13]['eta'].Fill(truthMuons[i].Eta) + + for i in range(len(truthWs)): + h_truth[24]['pT'].Fill(truthWs[i].PT) + h_truth[24]['p'].Fill(truthWs[i].P4().P()) + h_truth[24]['eta'].Fill(truthWs[i].Eta) - for p in tree.Particle: - if (p.Status==1 and abs(p.PID)==11): - T_electrons = T_electrons + 1 - T_e_pT.Fill(p.PT) - T_e_p.Fill(p.P4().P()) - T_e_eta.Fill(p.Eta) - elif (p.Status==1 and abs(p.PID)==13): - T_muons = T_muons + 1 - T_mu_pT.Fill(p.PT) - T_mu_p.Fill(p.P4().P()) - T_mu_eta.Fill(p.Eta) - if isBeamRemnant(p): - beamRemnants = beamRemnants + 1 - T_beamRemnants_pT.Fill(p.PT) - T_beamRemnants_p.Fill(p.P4().P()) - T_beamRemnants_eta.Fill(p.Eta) - if (p.PT <= 2.0): - Test_beamRemnants_pT.Fill(p.PT) - elif (abs(p.Status)==22 and (abs(p.PID)==24)): - W = W + 1 - T_W_pT.Fill(p.PT) - T_W_p.Fill(p.P4().P()) - T_W_eta.Fill(p.Eta) - elif (abs(p.Status)==22 and abs(p.PID)==23): - z = z + 1 - T_z_pT.Fill(p.PT) - T_z_p.Fill(p.P4().P()) - T_z_eta.Fill(p.Eta) - elif (p.Status==1 and abs(p.PID==22)): - T_photons = T_photons + 1 - T_gamma_pT.Fill(p.PT) - T_gamma_p.Fill(p.P4().P()) - T_gamma_eta.Fill(p.Eta) - elif (p.Status ==1 and (abs(p.PID)==12 or abs(p.PID)==14 or abs(p.PID)==16)): - neutrino_Px += p.Px - neutrino_Py += p.Py - neutrino_Pz += p.Pz - - T_e_multiplicity.Fill(T_electrons) - T_mu_multiplicity.Fill(T_muons) - T_W_multiplicity.Fill(W) - T_z_multiplicity.Fill(z) - T_gamma_multiplicity.Fill(T_photons) - T_beamRemnants_multiplicity.Fill(beamRemnants) - T_missingEt.Fill((neutrino_Px**2 + neutrino_Py**2)**.5) - T_missingE.Fill((neutrino_Px**2 + neutrino_Py**2 + neutrino_Pz**2)**.5) - - #reco-level - R_missingET.Fill(tree.MissingET[0].MET) - R_missingE.Fill(tree.MissingET[0].P4().P()) + for i in range(len(truthZs)): + h_truth[23]['pT'].Fill(truthZs[i].PT) + h_truth[23]['p'].Fill(truthZs[i].P4().P()) + h_truth[23]['eta'].Fill(truthZs[i].Eta) + + for i in range(len(truthPhotons)): + h_truth[22]['pT'].Fill(truthPhotons[i].PT) + h_truth[22]['p'].Fill(truthPhotons[i].P4().P()) + h_truth[22]['eta'].Fill(truthPhotons[i].Eta) + + for i in range(len(truthNeutrinos)): + h_truth[12]['pT'].Fill(truthNeutrinos[i].PT) + h_truth[12]['p'].Fill(truthNeutrinos[i].P4().P()) + h_truth[12]['eta'].Fill(truthNeutrinos[i].Eta) + + if len(beamRemnantMuons)>=2: + T_BRM_4vec=beamRemnantMuons[0].P4()+beamRemnantMuons[1].P4() + T_beamRemnants_mass.Fill(T_BRM_4vec.M()) + for i in range(len(beamRemnantMuons)): + T_beamRemnants_pT.Fill(beamRemnantMuons[i].PT) + T_beamRemnants_p.Fill(beamRemnantMuons[i].P4().P()) + T_beamRemnants_eta.Fill(beamRemnantMuons[i].Eta) + T_beamRemnants_cosh.Fill(math.cosh(beamRemnantMuons[i].Eta)/beamRemnantMuons[i].PT) + + for i in range(len(nonBeamRemnantMuons)): + T_nonBeamRemnants_pT.Fill(nonBeamRemnantMuons[i].PT) + T_nonBeamRemnants_p.Fill(nonBeamRemnantMuons[i].P4().P()) + T_nonBeamRemnants_eta.Fill(nonBeamRemnantMuons[i].Eta) + T_nonBeamRemnants_cosh.Fill(math.cosh(nonBeamRemnantMuons[i].Eta)/nonBeamRemnantMuons[i].PT) + + truthMissingP=TLorentzVector() + for nu in truthNeutrinos: + truthMissingP+=nu.P4() + + T_missingEt.Fill(truthMissingP.Et()) + T_missingE.Fill(truthMissingP.E()) + + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + #reco level + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + R_missingET.Fill(event.MissingET[0].MET) + R_missingE.Fill(event.MissingET[0].P4().P()) + #final and initial-state 4-vectors - P4_i = TLorentzVector(0,0,0,10000) - P4_f = [] - #multiplicity counters - electrons = 0 - muons = 0 - photons = 0 - e_list1 = [] - e_list2 = [] - mu_list1 = [] - mu_list2 = [] - gamma_list1 = [] - gamma_list2 = [] - #for OS pairs - leptons = [] - #fills electron histograms - for p in tree.Electron: - electrons = electrons + 1 - R_e_pT_I.Fill(p.PT) - R_e_p_I.Fill(p.P4().P()) - R_e_eta_I.Fill(p.Eta) - #assigns each particle to a list to be sorted by pT after going through the event to get leading, 2nd leading, and 3rd leading - e_list1.append(p.PT) - e_list2.append(p.Eta) - leptons.append(p) - P4_f.append(p.P4()) - - for p in tree.Muon: - muons = muons + 1 - R_mu_pT_I.Fill(p.PT) - R_mu_p_I.Fill(p.P4().P()) - R_mu_eta_I.Fill(p.Eta) - mu_list1.append(p.PT) - mu_list2.append(p.Eta) - leptons.append(p) - P4_f.append(p.P4()) - - for p in tree.Photon: - photons = photons + 1 - R_gamma_pT_I.Fill(p.PT) - R_gamma_p_I.Fill(p.P4().P()) - R_gamma_E_I.Fill(p.E) - R_gamma_eta_I.Fill(p.Eta) - gamma_list1.append(p.PT) - gamma_list2.append(p.Eta) - P4_f.append(p.P4()) - - R_e_multiplicity.Fill(electrons) - R_mu_multiplicity.Fill(muons) - R_gamma_multiplicity.Fill(photons) - final_P4 = TLorentzVector(0,0,0,0) - for i in P4_f: - final_P4 = final_P4 + i - missingMass = P4_i - final_P4 - R_missingMass.Fill(missingMass.M()) - - #sorts particles by pT and then fills the leading, 2nd order, and 3rd order histograms for pT and eta - if (len(e_list1) != 0): - z = zip(e_list1, e_list2) - zed = sorted(z, key=lambda x:x[0]) - zed.reverse() - if (len(e_list1) == 1): - R_e_pT_L.Fill(zed[0][0]) - R_e_eta_L.Fill(zed[0][1]) - if (len(e_list1) == 2): - R_e_pT_2.Fill(zed[1][0]) - R_e_eta_2.Fill(zed[1][1]) - if (len(e_list1) == 3): - R_e_pT_3.Fill(zed[2][0]) - R_e_eta_3.Fill(zed[2][1]) - - if (len(mu_list1) != 0): - z = zip(mu_list1, mu_list2) - zed = sorted(z, key=lambda x:x[0]) - zed.reverse() - if (len(mu_list1) == 1): - R_mu_pT_L.Fill(zed[0][0]) - R_mu_eta_L.Fill(zed[0][1]) - if (len(mu_list1) == 2): - R_mu_pT_2.Fill(zed[1][0]) - R_mu_eta_2.Fill(zed[1][1]) - if (len(mu_list1) == 3): - R_mu_pT_3.Fill(zed[2][0]) - R_mu_eta_3.Fill(zed[2][1]) - - if (len(gamma_list1) != 0): - z = zip(gamma_list1, gamma_list2) - zed = sorted(z, key=lambda x:x[0]) - zed.reverse() - if (len(gamma_list1) == 1): - R_gamma_pT_L.Fill(zed[0][0]) - R_gamma_eta_L.Fill(zed[0][1]) - if (len(gamma_list1) == 2): - R_gamma_pT_2.Fill(zed[1][0]) - R_gamma_eta_2.Fill(zed[1][1]) - if (len(gamma_list1) == 3): - R_gamma_pT_3.Fill(zed[2][0]) - R_gamma_eta_3.Fill(zed[2][1]) - - #electrons=selector(tree.Electron,'x.PT>5 and abs(x.Eta)<2') - #muons=selector(tree.Muon,'x.PT>5 and abs(x.Eta)<2') - - #fMuons=selector(tree.Muon,'abs(x.Eta)>2') - #fMuons+=selector(tree.Electron,'abs(x.Eta)>2 and x.Particle.GetObject().PID==11 and x.Particle.GetObject().Status==1 and x.Particle.GetObject().M1<5') #this is a hack - not a typo + P4_i=TLorentzVector() + P4_i+=event.Particle[0].P4() + P4_i+=event.Particle[1].P4() + + P4_f = TLorentzVector() + + #electrons=selector(event.Electron,'x.PT>5 and abs(x.Eta)<2') + #muons=selector(event.Muon,'x.PT>5 and abs(x.Eta)<2') + + #fMuons=selector(event.Muon,'abs(x.Eta)>2') + #fMuons+=selector(event.Electron,'abs(x.Eta)>2 and x.Particle.GetObject().PID==11 and x.Particle.GetObject().Status==1 and x.Particle.GetObject().M1<5') #this is a hack - not a typo + + h_reco[11]['mult'].Fill(len(electrons)) + for i in range(len(electrons)): + P4_f+=electrons[i].P4() + h_reco[11]['pT']['I'].Fill(electrons[i].PT) + h_reco[11]['p']['I'].Fill(electrons[i].P4().P()) + h_reco[11]['eta']['I'].Fill(electrons[i].Eta) + if i<4: + h_reco[11]['pT'][i].Fill(electrons[i].PT) + h_reco[11]['p'][i].Fill(electrons[i].P4().P()) + h_reco[11]['eta'][i].Fill(electrons[i].Eta) + + h_reco[13]['mult'].Fill(len(muons)) + if len(muons)>=2: + R_muon_leading_cos.Fill(muons[0].P4().CosTheta()) + R_muon_sub_cos.Fill(muons[1].P4().CosTheta()) + if len(muons)==1: + R_muon_leading_cos.Fill(muons[0].P4().CosTheta()) + R_mumu_p.Fill(muon_4vec.P()) + R_mumu_cos.Fill(muon_4vec.CosTheta()) + for i in range(len(muons)): + P4_f+=muons[i].P4() + h_reco[13]['pT']['I'].Fill(muons[i].PT) + h_reco[13]['p']['I'].Fill(muons[i].P4().P()) + h_reco[13]['eta']['I'].Fill(muons[i].Eta) + if i<4: + h_reco[13]['pT'][i].Fill(muons[i].PT) + h_reco[13]['p'][i].Fill(muons[i].P4().P()) + h_reco[13]['eta'][i].Fill(muons[i].Eta) + + R_beamRemnants_multiplicity.Fill(len(reco_beamRemnants)) + if len(reco_beamRemnants)>=2: + BRM_4vec=reco_beamRemnants[0].P4() + reco_beamRemnants[1].P4() + R_beamRemnants_mass.Fill(BRM_4vec.M()) + for i in range(len(reco_beamRemnants)): + R_beamRemnants_pT.Fill(reco_beamRemnants[i].PT) + R_beamRemnants_p.Fill(reco_beamRemnants[i].P4().P()) + R_beamRemnants_eta.Fill(reco_beamRemnants[i].Eta) + R_beamRemnants_cosh.Fill(math.cosh(reco_beamRemnants[i].Eta)/reco_beamRemnants[i].PT) + + for i in range(len(leptons)): + if ((abs(leptons[i].Eta) <= 2) and (leptons[i].P4().P() < 100)): + reco_nonBeamRemnants.append(leptons[i]) + + R_nonBeamRemnants_multiplicity.Fill(len(reco_nonBeamRemnants)) + for i in range(len(reco_nonBeamRemnants)): + R_nonBeamRemnants_pT.Fill(reco_nonBeamRemnants[i].PT) + R_nonBeamRemnants_p.Fill(reco_nonBeamRemnants[i].P4().P()) + R_nonBeamRemnants_eta.Fill(reco_nonBeamRemnants[i].Eta) + R_nonBeamRemnants_cosh.Fill(math.cosh(reco_nonBeamRemnants[i].Eta)/reco_nonBeamRemnants[i].PT) + + h_reco[22]['mult'].Fill(len(photons)) + for i in range(len(photons)): + P4_f+=photons[i].P4() + h_reco[22]['pT']['I'].Fill(photons[i].PT) + h_reco[22]['p']['I'].Fill(photons[i].P4().P()) + h_reco[22]['eta']['I'].Fill(photons[i].Eta) + if i<4: + h_reco[22]['pT'][i].Fill(photons[i].PT) + h_reco[22]['p'][i].Fill(photons[i].P4().P()) + h_reco[22]['eta'][i].Fill(photons[i].Eta) + + h_reco[21]['mult'].Fill(len(jets)) + for i in range(len(jets)): + P4_f+=jets[i].P4() + h_reco[21]['pT']['I'].Fill(jets[i].PT) + h_reco[21]['p']['I'].Fill(jets[i].P4().P()) + h_reco[21]['eta']['I'].Fill(jets[i].Eta) + h_reco[21]['mass']['I'].Fill(jets[i].P4().M()) + if i<4: + h_reco[21]['pT'][i].Fill(jets[i].PT) + h_reco[21]['p'][i].Fill(jets[i].P4().P()) + h_reco[21]['eta'][i].Fill(jets[i].Eta) + h_reco[21]['mass'][i].Fill(jets[i].P4().M()) + + R_Wjets_multiplicity.Fill(len(W_jets)) + if len(W_jets)>=2: + Wjets_4vec=W_jets[0].P4() + W_jets[1].P4() + R_Wjets_mass.Fill(Wjets_4vec.M()) + R_Wjets_cos_leading.Fill(W_jets[0].P4().CosTheta()) + R_Wjets_cos_sub.Fill(W_jets[1].P4().CosTheta()) + R_Wjets_p_leading.Fill(W_jets[0].P4().P()) + R_Wjets_p_sub.Fill(W_jets[1].P4().P()) + R_Wjets_cos_pair.Fill(Wjets_4vec.P()) + R_Wjets_p_pairs.Fill(Wjets_4vec.CosTheta()) + if len(W_jets)==1: + R_Wjets_cos_leading.Fill(W_jets[0].P4().CosTheta()) + R_Wjets_p_leading.Fill(W_jets[0].P4().CosTheta()) + for i in range(len(W_jets)): + R_Wjets_pT.Fill(W_jets[i].PT) + R_Wjets_p.Fill(W_jets[i].P4().P()) + R_Wjets_eta.Fill(W_jets[i].Eta) + + R_missingMass.Fill((P4_f-P4_i).M()) + + #------------------------------------------------------------------------------------------------- + leptons=electrons+muons + + Zs={} + for i1 in range(len(leptons)-1): + l1=leptons[i1] + Zs[i1]={} + for i2 in range(i1+1,len(leptons)): + l2=leptons[i2] + + if (l1.Charge==l2.Charge) or (type(l1) != (type(l2))): Zs[i1][i2]=None + else: Zs[i1][i2]=parentConstructor(l1,l2) + + minimum=9E9 + theZ=None + for i1 in range(len(Zs.keys())): + Zcand=Zs[i1]={} + for i2 in range(len(Zs[i1].keys())): + Zcand=Zs[i1][i2] + if Zcand: + if abs(Zcand.Mass()-91.1876) [delphes card]" - echo "Example: ./generateEvents.sh $PWD/foo $PWD/delphes/cards/gen_card.tcl" +if [[ $# < 1 ]]; +then + echo "Usage: ./generateEvents.sh [-m madgraphScript] [-d delphesCard] [-a analysisScript]" + echo "Example: ./generateEvents.sh -m $PWD/test.mg -d $PWD/delphes/cardsdelphes_card_MuonColliderDet.tcl" exit else - mgScript=$1 - if [[ $# -gt 1 ]]; then - delphesCard=$2 - else - delphesCard="cards/gen_card.tcl" - fi + while getopts m:d:a: flag; do + echo "flag -$flag, Argument $OPTARG"; + case "$flag" in + m) mgScript=$OPTARG;; + d) delphesCard=$OPTARG;; + a) analysisScript=$OPTARG;; + esac + done fi -source setup.sh +#restore environment +. ~/.aQGCEnv -touch dummy +touch $aQGCWorkDir/.dummy ######################################################################### -cd MG5_aMC_v3_1_1 -python ./bin/mg5_aMC < "${mgScript}" -gzs=`find . -newer ../dummy -name "unweighted_events.lhe.gz" -exec echo $PWD/{} \;` -echo $gzs +python $madgraphDir/bin/mg5_aMC < $mgScript -cd .. +gzs=`find . -newer $aQGCWorkDir/.dummy -name "unweighted_events.lhe.gz" -exec echo {} \;` +echo $gzs #------------------------------------------------------------------------ -cd delphes for gz in $gzs; do lhe=${gz%%.gz} output=${lhe%%.lhe}.root #set delphes output path/name - + gunzip $gz n=`grep -c \ $lhe` echo $n - sed s%examples/Pythia8/events.lhe%$lhe% examples/Pythia8/configLHE.cmnd > configLHE.cmnd #this will create a new config pointing to your lhe - sed "s%Main:numberOfEvents = 10%Main:numberOfEvents = $n%" --in-place configLHE.cmnd - ./DelphesPythia8 $delphesCard configLHE.cmnd $output #this runs delphes using your new config + sed s%examples/Pythia8/events.lhe%$lhe% $delphesDir/examples/Pythia8/configLHE.cmnd > configLHE.cmnd #this will create a new config pointing to your lhe + sed -in-place "s%Main:numberOfEvents = 10%Main:numberOfEvents = $n%" configLHE.cmnd + $delphesDir/DelphesPythia8 $delphesCard configLHE.cmnd $output #this runs delphes using your new config + python $analysisScript $output done diff --git a/install.sh b/install.sh index 33e5439..2f2f415 100755 --- a/install.sh +++ b/install.sh @@ -1,7 +1,22 @@ -if [[ `basename $PWD` != "aQGC" ]]; then echo "Execute from aQGC dir"; exit; fi +export aQGCWorkDir=$PWD -baseDir=$PWD -source setup.sh +if [[ $HOSTNAME == "login.snowmass21.io" ]]; then + #important + module load python/2.7.15 + module load py-numpy/1.15.2-py2.7 + module load py-six/1.11.0-py2.7 + . /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.08/x86_64-centos7-gcc48-opt/bin/thisroot.sh + + #convenient + module load emacs +else + which root >> /dev/null + if [[ $? -ne 0 ]]; then + echo "Please ensure ROOT is available and retry" + fi +fi + +#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #MadGraph wget https://launchpad.net/mg5amcnlo/3.0/3.1.x/+download/MG5_aMC_v3.1.1.tar.gz @@ -9,32 +24,42 @@ tar -xzvf MG5_aMC_v3.1.1.tar.gz if [[ $? -ne 0 ]]; then echo "ERROR getting madgraph" exit +else + export madgraphDir=$aQGCWorkDir/MG5_aMC_v3_1_1 + rm MG5_aMC_v3.1.1.tar.gz fi -rm MG5_aMC_v3.1.1.tar.gz #Pythia wget https://pythia.org/download/pythia83/pythia8306.tgz tar -xzvf pythia8306.tgz rm pythia8306.tgz cd pythia8306 -./configure --prefix=$baseDir +./configure --prefix=$aQGCWorkDir make -j make install if [[ $? -ne 0 ]]; then echo "ERROR compiling pythia" exit +else + export pythiaDir=$aQGCWorkDir/pythia8306 + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${aQGCWorkDir}/lib + export PYTHIA8=$aQGCWorkDir + rm pythia8306.tgz + cd .. fi -cd .. #Delphes -export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${baseDir}/lib -export PYTHIA8=$baseDir git clone https://github.com/delphes/delphes.git cd delphes -make HAS_PYTHIA8=true -j -I${baseDir}/include/Pythia8 +make HAS_PYTHIA8=true -j -I${aQGCWorkDir}/include/Pythia8 if [[ $? -ne 0 ]]; then echo "ERROR compiling delphes" exit +else + export delphesDir=$aQGCWorkDir/delphes + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${delphesDir} + export ROOT_INCLUDE_PATH=ROOT_INCLUDE_PATH:${delphesDir}/external/ + cd .. fi -#echo "export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${baseDir}/lib" >> ~/.bash_profile +export -p > ~/.aQGCEnv diff --git a/overlay.py b/overlay.py new file mode 100644 index 0000000..f49f195 --- /dev/null +++ b/overlay.py @@ -0,0 +1,36 @@ +from ROOT import * +from sys import argv + +norm=False + +inputs=argv[1:] + +f0=TFile(inputs[0]) +keys=f0.GetListOfKeys() + +for k in keys: + c=TCanvas() + m=0 + for f in inputs: + + i=inputs.index(f) + f=TFile(f) + SetOwnership(f,False) + o=f.Get(k.GetName()) + + if type(o)!=type(TH1F()): continue + + try: + if norm: o.Scale(1./o.Integral(0,o.GetNbinsX()+2)) + except: continue + o.SetLineColor(i+1) + o.SetLineStyle(i+1) + if i==0: + first=o + o.Draw() + else: o.Draw("SAME") + + m=max(m,o.GetMaximum()) + first.SetMaximum(1.25*m) + c.Modified() + c.SaveAs(o.GetName()+'.pdf') diff --git a/runAnalysis.sh b/runAnalysis.sh new file mode 100755 index 0000000..8897875 --- /dev/null +++ b/runAnalysis.sh @@ -0,0 +1,19 @@ +#!/bin/bash -x + +if [ $# -lt 2 ]; then + echo "Usage: ./runAnalysis.sh [analysis script]" + echo "Example: ./runAnalysis.sh 'PROC_sm_*/Events/run_*/unweighted_events.root' foo" + exit +fi + +inputs=$1 +tag=$2 +if [ $# -gt 2 ]; then + analysisScript=$3 +else + analysisScript=$aQGCWorkDir/aQGC/analyzeDelphes.py +fi + +for input in `ls $inputs`; do + python $analysisScript $input $tag +done diff --git a/setup.sh b/setup.sh index af0dbf9..9126209 100644 --- a/setup.sh +++ b/setup.sh @@ -1,20 +1,3 @@ -if [[ `basename $PWD` != "aQGC" ]]; then echo "Execute from aQGC dir"; exit; fi +#!/bin/bash -export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${PWD}/lib -export prodBase=$PWD - -if [[ $HOSTNAME == "login.snowmass21.io" ]]; then - #important - module load python/2.7.15 - module load py-numpy/1.15.2-py2.7 - module load py-six/1.11.0-py2.7 - . /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.08/x86_64-centos7-gcc48-opt/bin/thisroot.sh - - #convenient - module load emacs -else - which root >> /dev/null - if [[ $? -ne 0 ]]; then - echo "Please ensure ROOT is available and retry" - fi -fi +. ~/.aQGCEnv diff --git a/test.mg b/test.mg new file mode 100644 index 0000000..3177008 --- /dev/null +++ b/test.mg @@ -0,0 +1,4 @@ +generate mu+ mu- > Z > mu+ mu- +add process mu+ mu- > Z > u u~ +output +launch \ No newline at end of file diff --git a/test.sh b/test.sh new file mode 100755 index 0000000..465a68f --- /dev/null +++ b/test.sh @@ -0,0 +1,8 @@ +#!/bin/bash -x + +git clone -b cleanUp https://github.com/jstupak/aQGC.git + +./aQGC/install.sh +source aQGC/setup.sh +./aQGC/generateEvents.sh -m $aQGCWorkDir/aQGC/test.mg -d $delphesDir/cards/delphes_card_MuonColliderDet.tcl -a $aQGCWorkDir/aQGC/analyzeDelphes.py +