From ee1bcd103db580e13e01b675f65a3d2f2e85fd7e Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 11:02:17 -0600 Subject: [PATCH 01/22] changing envirnment handling --- install.sh | 43 ++++++++++++++++++++++++++++++++----------- setup.sh | 21 ++------------------- 2 files changed, 34 insertions(+), 30 deletions(-) diff --git a/install.sh b/install.sh index 33e5439..6e6d478 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,38 @@ 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/pythia830g + 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 + cd .. + export -p > ~/.prodEnv fi - -#echo "export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${baseDir}/lib" >> ~/.bash_profile diff --git a/setup.sh b/setup.sh index af0dbf9..07f29b0 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 +. ~/.prodEnv From 30d04ba68f86e72bf112735350ddf73974ed3c51 Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 11:16:23 -0600 Subject: [PATCH 02/22] bugfix --- MG5_aMC_v3_1_1/bin/mg5_aMC | 227 +++++++++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100755 MG5_aMC_v3_1_1/bin/mg5_aMC 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 + + + From 6b2cf77c43e4c9159e7a9fb28725d56bda006e04 Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 11:17:37 -0600 Subject: [PATCH 03/22] work on environment handling --- generateEvents.sh | 11 +++++------ install.sh | 2 +- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/generateEvents.sh b/generateEvents.sh index 77f5a11..87fc6cd 100755 --- a/generateEvents.sh +++ b/generateEvents.sh @@ -1,6 +1,5 @@ #!/bin/bash -x -if [[ `basename $PWD` != "MCProd" ]]; then echo "Execute from MCProd dir"; exit; fi if [[ $# < 1 ]]; then echo "Usage: ./generateEvents.sh [delphes card]" @@ -17,20 +16,20 @@ fi source setup.sh -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/{} \;` +python $madgraphDir/bin/mg5_aMC < "${mgScript}" + +gzs=`find . -newer $aQGCWorkDir/.dummy -name "unweighted_events.lhe.gz" -exec echo $PWD/{} \;` echo $gzs cd .. #------------------------------------------------------------------------ -cd delphes +cd $delphesDir for gz in $gzs; do lhe=${gz%%.gz} output=${lhe%%.lhe}.root #set delphes output path/name diff --git a/install.sh b/install.sh index 6e6d478..68ef24e 100755 --- a/install.sh +++ b/install.sh @@ -25,7 +25,7 @@ if [[ $? -ne 0 ]]; then echo "ERROR getting madgraph" exit else - export madgraphDir=$aQGCWorkDir/MG5_aMC_v3.1.1 + export madgraphDir=$aQGCWorkDir/MG5_aMC_v3_1_1 rm MG5_aMC_v3.1.1.tar.gz fi From 93b2fc376e8a11960d0170daefe51317d0a8061f Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 11:31:03 -0600 Subject: [PATCH 04/22] bugfix --- generateEvents.sh | 9 +++++---- install.sh | 4 +++- setup.sh | 2 +- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/generateEvents.sh b/generateEvents.sh index 87fc6cd..983da64 100755 --- a/generateEvents.sh +++ b/generateEvents.sh @@ -14,7 +14,8 @@ else fi fi -source setup.sh +#source setup.sh +. ~/.aQGCEnv touch $aQGCWorkDir/.dummy @@ -29,7 +30,7 @@ cd .. #------------------------------------------------------------------------ -cd $delphesDir +#cd $delphesDir for gz in $gzs; do lhe=${gz%%.gz} output=${lhe%%.lhe}.root #set delphes output path/name @@ -37,7 +38,7 @@ for gz in $gzs; do 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%examples/Pythia8/events.lhe%$lhe% $delphesDir/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 + $delphesDir/DelphesPythia8 $delphesCard configLHE.cmnd $output #this runs delphes using your new config done diff --git a/install.sh b/install.sh index 68ef24e..0e8a2bf 100755 --- a/install.sh +++ b/install.sh @@ -56,6 +56,8 @@ if [[ $? -ne 0 ]]; then echo "ERROR compiling delphes" exit else + export delphesDir=$aQGCWorkDir/delphes cd .. - export -p > ~/.prodEnv fi + +export -p > ~/.aQGCEnv diff --git a/setup.sh b/setup.sh index 07f29b0..9126209 100644 --- a/setup.sh +++ b/setup.sh @@ -1,3 +1,3 @@ #!/bin/bash -. ~/.prodEnv +. ~/.aQGCEnv From e2ca11dba7ecf7eda09cecdd278a5acaa4826c7c Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 11:47:21 -0600 Subject: [PATCH 05/22] bugfix --- generateEvents.sh | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/generateEvents.sh b/generateEvents.sh index 983da64..d1fc36d 100755 --- a/generateEvents.sh +++ b/generateEvents.sh @@ -10,7 +10,7 @@ else if [[ $# -gt 1 ]]; then delphesCard=$2 else - delphesCard="cards/gen_card.tcl" + delphesCard="$delphesDir/cards/gen_card.tcl" fi fi @@ -23,18 +23,15 @@ touch $aQGCWorkDir/.dummy python $madgraphDir/bin/mg5_aMC < "${mgScript}" -gzs=`find . -newer $aQGCWorkDir/.dummy -name "unweighted_events.lhe.gz" -exec echo $PWD/{} \;` +gzs=`find . -newer $aQGCWorkDir/.dummy -name "unweighted_events.lhe.gz" -exec echo {} \;` echo $gzs -cd .. - #------------------------------------------------------------------------ -#cd $delphesDir 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 From 0094e481638531c5ee7f4ebe6078dd031f5f05fc Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 11:49:02 -0600 Subject: [PATCH 06/22] adding exmple amdgraph script for testing --- test.mg | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 test.mg diff --git a/test.mg b/test.mg new file mode 100644 index 0000000..b0ea87e --- /dev/null +++ b/test.mg @@ -0,0 +1,3 @@ +generate mu+ mu- > Z +output +launch \ No newline at end of file From bed2c7b020f7fc9c563dba9d9d3bfd2b519f90ac Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 12:23:47 -0600 Subject: [PATCH 07/22] adding delphes script to generateEvents.sh --- generateEvents.sh | 29 ++++++++++++++++------------- install.sh | 1 + 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/generateEvents.sh b/generateEvents.sh index d1fc36d..5a5b32b 100755 --- a/generateEvents.sh +++ b/generateEvents.sh @@ -1,27 +1,29 @@ #!/bin/bash -x -if [[ $# < 1 ]]; -then - echo "Usage: ./generateEvents.sh [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/foo -d $PWD/delphes/cards/gen_card.tcl" exit else - mgScript=$1 - if [[ $# -gt 1 ]]; then - delphesCard=$2 - else - delphesCard="$delphesDir/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 $aQGCWorkDir/.dummy ######################################################################### -python $madgraphDir/bin/mg5_aMC < "${mgScript}" +python $madgraphDir/bin/mg5_aMC < $mgScript gzs=`find . -newer $aQGCWorkDir/.dummy -name "unweighted_events.lhe.gz" -exec echo {} \;` echo $gzs @@ -31,11 +33,12 @@ echo $gzs 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% $delphesDir/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 $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 0e8a2bf..50b8199 100755 --- a/install.sh +++ b/install.sh @@ -57,6 +57,7 @@ if [[ $? -ne 0 ]]; then exit else export delphesDir=$aQGCWorkDir/delphes + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${delphes} cd .. fi From 99703394ebb02118ef31bfbdf92880503c5b238f Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 17:44:27 -0600 Subject: [PATCH 08/22] cleaning up/debugging analysis script --- analyzeDelphes.py | 467 +++++++++++++++++++++------------------------- 1 file changed, 214 insertions(+), 253 deletions(-) diff --git a/analyzeDelphes.py b/analyzeDelphes.py index 87e4f32..fbe76e9 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -1,11 +1,39 @@ maxEvents=9E9 -DEBUG=True +DEBUG=False + +jetType='KTjet' + +particle={ + 11:"e", + 13:"mu", + 21:"gamma", #using this for all jets + 22:"j", + 23:"z", + 24:"w" +} + +#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 + +#--------------------------------------------------------- import pdb from sys import argv #--------------------------------------------------------- from ROOT import * -from ROOT import TLorentzVector gSystem.Load("libDelphes") gStyle.SetOptStat(0) @@ -14,7 +42,7 @@ def printHist(h): for i in range(h.GetNbinsX()+2): print h.GetBinContent(i), - + #--------------------------------------------------------- #return TLorentzVector corresponding to sum of inputs @@ -50,295 +78,228 @@ def isBeamRemnant(p): #--------------------------------------------------------- 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)] - - + + inputName=str(argv[1]) + if len(argv)>2: + tag=str(argv[2]) + else: + tag='' + outputName=(bool(tag)*(tag+'_'))+inputName.replace('.root','.hist.root') + + print "inputName:",inputName + print "tag:",tag + + f=TFile(inputName) tree=f.Get("Delphes") - name=str(argv[2]) - output=TFile(name+".delphes.root","RECREATE") + output=TFile(outputName,"RECREATE") #sets upper limit for bin range bin_range = 15000 #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]: + h_truth[p]={} + h_truth[p]['pT'] =TH1F('T_%s_pT'%particle[p],';Truth %s pT [GeV];Events'%particle[p], 200, 0, bin_range) + h_truth[p]['p'] =TH1F('T_%s_p'%particle[p],';Truth %s p [GeV];Events'%particle[p], 200, 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) + #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]['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 pT [GeV];Events'%particle[p], 200, 0, bin_range) + h_reco[p]['eta'][i]=TH1F('R_%s_eta_%s'%(particle[p],i),';%s pT [GeV];Events'%particle[p], 20, -4, 4) + h_reco[p]['mult']=TH1F('R_%s_mult'%particle[p],';%s multiplicity;Events'%particle[p], 7, -0.5, 6.5) #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', 200, 0, 200) 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_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 - 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_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_multiplicity = TH1F('T_beamRemnants_multiplicity', ';Multiplicity;Events', 7, -.5, 6.5) - Test_beamRemnants_pT = TH1F('Test_beamRemnants_Pt', ';pT (GeV);Events', 30, 0, 3) + 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 - - 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 + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + #truth level + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + truthElectrons=selector(tree.Particle,'x.Status==1 and abs(x.PID)==11') + truthMuons =selector(tree.Particle,'x.Status==1 and abs(x.PID)==13') + truthWs =selector(tree.Particle,'x.Status==22 and abs(x.PID)==24') + truthZs =selector(tree.Particle,'x.Status==22 and abs(x.PID)==23') + truthPhotons =selector(tree.Particle,'x.Status==1 and abs(x.PID)==22') + truthNeutrinos=selector(tree.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)') + + T_e_multiplicity.Fill(len(truthElectrons)) + T_mu_multiplicity.Fill(len(truthMuons)) + T_W_multiplicity.Fill(len(truthWs)) + T_z_multiplicity.Fill(len(truthZs)) + T_gamma_multiplicity.Fill(len(truthPhotons)) + T_beamRemnants_multiplicity.Fill(len(beamRemnantMuons)) + + truthMissingP=TLorentzVector() + for nu in truthNeutrinos: + truthMissingP+=nu.P4() + + T_missingEt.Fill(truthMissingP.Et()) + T_missingE.Fill(truthMissingP.E()) + + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + #reco level + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + R_missingET.Fill(tree.MissingET[0].MET) R_missingE.Fill(tree.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]) + P4_i=TLorentzVector() + P4_i+=Particle[0].P4() + P4_i+=Particle[1].P4() + + P4_f = TLorentzVector() + electrons=selector(tree.Electron, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[11],etaMax[11])) + muons =selector(tree.Muon, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[13],etaMax[13])) + photons =selector(tree.Photon, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[22],etaMax[22])) + jets =selector(tree.__getattribute__(jetType),'x.PT>%f and abs(x.Eta)<%f'%(pTMin[21],etaMax[21])) + #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 + + #sort object lists + for collection in [electrons, muons, photons, jets]: collection.sort(key=operator.attrgetter('PT'),reverse=True) + + 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].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].P) + h_reco[11]['eta'][i].Fill(electrons[i].Eta) + + h_reco[13]['mult'].Fill(len(muons)) + 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].PT) + h_reco[13]['eta']['I'].Fill(muons[i].PT) + if i<4: + h_reco[13]['pT'][i].Fill(muons[i].PT) + h_reco[13]['p'][i].Fill(muons[i].P) + h_reco[13]['eta'][i].Fill(muons[i].Eta) + + 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].PT) + h_reco[22]['eta']['I'].Fill(photons[i].PT) + if i<4: + h_reco[22]['pT'][i].Fill(photons[i].PT) + h_reco[22]['p'][i].Fill(photons[i].P) + h_reco[22]['eta'][i].Fill(photons[i].Eta) + + h_reco[11]['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].PT) + h_reco[21]['eta']['I'].Fill(jets[i].PT) + if i<4: + h_reco[21]['pT'][i].Fill(jets[i].PT) + h_reco[21]['p'][i].Fill(jets[i].P) + h_reco[21]['eta'][i].Fill(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())): + for i2 in range(len(Zs[i1].keys())): + Zcand=Zs[i1][i2] + if Zcand: + if abs(Zcand.Mass()-91.1876) Date: Mon, 8 Nov 2021 18:06:17 -0600 Subject: [PATCH 09/22] bugfix --- install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/install.sh b/install.sh index 50b8199..91af35b 100755 --- a/install.sh +++ b/install.sh @@ -57,7 +57,7 @@ if [[ $? -ne 0 ]]; then exit else export delphesDir=$aQGCWorkDir/delphes - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${delphes} + export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${delphesDir} cd .. fi From 4ab971afd96903eae377edf61c13ccb75a97daf9 Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 18:13:08 -0600 Subject: [PATCH 10/22] bugfix --- install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/install.sh b/install.sh index 91af35b..404c1e2 100755 --- a/install.sh +++ b/install.sh @@ -41,7 +41,7 @@ if [[ $? -ne 0 ]]; then echo "ERROR compiling pythia" exit else - export pythiaDir=$aQGCWorkDir/pythia830g + export pythiaDir=$aQGCWorkDir/pythia8306 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${aQGCWorkDir}/lib export PYTHIA8=$aQGCWorkDir rm pythia8306.tgz From 9dd21a244582d57ae9016e41223f168f1d36c556 Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 18:44:28 -0600 Subject: [PATCH 11/22] analyzeDelphes.py --- analyzeDelphes.py | 1 - 1 file changed, 1 deletion(-) diff --git a/analyzeDelphes.py b/analyzeDelphes.py index fbe76e9..ad8a907 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -122,7 +122,6 @@ def isBeamRemnant(p): 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 pT [GeV];Events'%particle[p], 200, 0, bin_range) h_reco[p]['eta'][i]=TH1F('R_%s_eta_%s'%(particle[p],i),';%s pT [GeV];Events'%particle[p], 20, -4, 4) - h_reco[p]['mult']=TH1F('R_%s_mult'%particle[p],';%s multiplicity;Events'%particle[p], 7, -0.5, 6.5) #histograms for OS pairs R_ee_pT = TH1F('ee_pT', ';pT [GeV];Events', 200, 0, bin_range) From 284b1051fbe4cf3b3201f9616ffa4018b5fcb345 Mon Sep 17 00:00:00 2001 From: John Stupak Date: Mon, 8 Nov 2021 19:24:19 -0600 Subject: [PATCH 12/22] adding anther environment var --- install.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/install.sh b/install.sh index 404c1e2..fb9c7ca 100755 --- a/install.sh +++ b/install.sh @@ -58,6 +58,7 @@ if [[ $? -ne 0 ]]; then else export delphesDir=$aQGCWorkDir/delphes export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${delphesDir} + export export ROOT_INCLUDE_PATH=ROOT_INCLUDE_PATH:${delphesDir}/external/ cd .. fi From 2c8fae60fe78aa692376ce79a7147d27246367c2 Mon Sep 17 00:00:00 2001 From: John Stupak Date: Tue, 9 Nov 2021 10:08:13 -0600 Subject: [PATCH 13/22] bug fixes --- analyzeDelphes.py | 93 +++++++++++++++++++++++++---------------------- generateEvents.sh | 4 +- overlay.py | 36 ++++++++++++++++++ test.mg | 3 +- test.sh | 8 ++++ 5 files changed, 98 insertions(+), 46 deletions(-) create mode 100644 overlay.py create mode 100755 test.sh diff --git a/analyzeDelphes.py b/analyzeDelphes.py index ad8a907..35eaf08 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -29,8 +29,11 @@ #--------------------------------------------------------- #EVENT SELECTION +#ADD EVENT SELECTION CRITERIA HERE + #--------------------------------------------------------- import pdb +import operator from sys import argv #--------------------------------------------------------- from ROOT import * @@ -62,11 +65,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 @@ -90,11 +93,15 @@ def isBeamRemnant(p): print "tag:",tag f=TFile(inputName) - tree=f.Get("Delphes") + 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/10 #declares truth-level pT, p, eta, and multiplicity histograms for e,mu,W,Z,gamma h_truth={} @@ -156,28 +163,27 @@ def isBeamRemnant(p): #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - for event in range(min(tree.GetEntries(),maxEvents)): - tree.GetEntry(event) + for event in f.Delphes: #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #truth level #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - truthElectrons=selector(tree.Particle,'x.Status==1 and abs(x.PID)==11') - truthMuons =selector(tree.Particle,'x.Status==1 and abs(x.PID)==13') - truthWs =selector(tree.Particle,'x.Status==22 and abs(x.PID)==24') - truthZs =selector(tree.Particle,'x.Status==22 and abs(x.PID)==23') - truthPhotons =selector(tree.Particle,'x.Status==1 and abs(x.PID)==22') - truthNeutrinos=selector(tree.Particle,'x.Status==1 and (abs(x.PID)==12 or abs(x.PID)==14 or abs(x.PID)==16)') + 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,'x.Status==22 and abs(x.PID)==24') + truthZs =selector(event.Particle,'x.Status==22 and abs(x.PID)==23') + truthPhotons =selector(event.Particle,'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)') - - T_e_multiplicity.Fill(len(truthElectrons)) - T_mu_multiplicity.Fill(len(truthMuons)) - T_W_multiplicity.Fill(len(truthWs)) - T_z_multiplicity.Fill(len(truthZs)) - T_gamma_multiplicity.Fill(len(truthPhotons)) + + 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)) T_beamRemnants_multiplicity.Fill(len(beamRemnantMuons)) truthMissingP=TLorentzVector() @@ -191,26 +197,26 @@ def isBeamRemnant(p): #reco level #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - R_missingET.Fill(tree.MissingET[0].MET) - R_missingE.Fill(tree.MissingET[0].P4().P()) + R_missingET.Fill(event.MissingET[0].MET) + R_missingE.Fill(event.MissingET[0].P4().P()) #final and initial-state 4-vectors P4_i=TLorentzVector() - P4_i+=Particle[0].P4() - P4_i+=Particle[1].P4() + P4_i+=event.Particle[0].P4() + P4_i+=event.Particle[1].P4() P4_f = TLorentzVector() - electrons=selector(tree.Electron, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[11],etaMax[11])) - muons =selector(tree.Muon, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[13],etaMax[13])) - photons =selector(tree.Photon, 'x.PT>%f and abs(x.Eta)<%f'%(pTMin[22],etaMax[22])) - jets =selector(tree.__getattribute__(jetType),'x.PT>%f and abs(x.Eta)<%f'%(pTMin[21],etaMax[21])) - - #electrons=selector(tree.Electron,'x.PT>5 and abs(x.Eta)<2') - #muons=selector(tree.Muon,'x.PT>5 and abs(x.Eta)<2') + 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])) + + #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(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 + #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 #sort object lists for collection in [electrons, muons, photons, jets]: collection.sort(key=operator.attrgetter('PT'),reverse=True) @@ -219,45 +225,45 @@ def isBeamRemnant(p): 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].P) + 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].P) + 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)) 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].PT) - h_reco[13]['eta']['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].P) + h_reco[13]['p'][i].Fill(muons[i].P4().P()) h_reco[13]['eta'][i].Fill(muons[i].Eta) 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].PT) - h_reco[22]['eta']['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].P) + h_reco[22]['p'][i].Fill(photons[i].P4().P()) h_reco[22]['eta'][i].Fill(photons[i].Eta) h_reco[11]['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].PT) - h_reco[21]['eta']['I'].Fill(jets[i].PT) + h_reco[21]['p']['I'].Fill(jets[i].P4().P()) + h_reco[21]['eta']['I'].Fill(jets[i].Eta) if i<4: h_reco[21]['pT'][i].Fill(jets[i].PT) - h_reco[21]['p'][i].Fill(jets[i].P) + h_reco[21]['p'][i].Fill(jets[i].P4().P()) h_reco[21]['eta'][i].Fill(jets[i].Eta) R_missingMass.Fill((P4_f-P4_i).M()) @@ -278,6 +284,7 @@ def isBeamRemnant(p): 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: diff --git a/generateEvents.sh b/generateEvents.sh index 5a5b32b..f05a6be 100755 --- a/generateEvents.sh +++ b/generateEvents.sh @@ -3,7 +3,7 @@ if [[ $# < 1 ]]; then echo "Usage: ./generateEvents.sh [-m madgraphScript] [-d delphesCard] [-a analysisScript]" - echo "Example: ./generateEvents.sh -m $PWD/foo -d $PWD/delphes/cards/gen_card.tcl" + echo "Example: ./generateEvents.sh -m $PWD/test.mg -d $PWD/delphes/cardsdelphes_card_MuonColliderDet.tcl" exit else while getopts m:d:a: flag; do @@ -38,7 +38,7 @@ for gz in $gzs; do n=`grep -c \ $lhe` echo $n 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 "s%Main:numberOfEvents = 10%Main:numberOfEvents = $n%" --in-place configLHE.cmnd + 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/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/test.mg b/test.mg index b0ea87e..3177008 100644 --- a/test.mg +++ b/test.mg @@ -1,3 +1,4 @@ -generate mu+ mu- > Z +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..d6593f5 --- /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/test.mg -d $delphesDir/cards/delphes_card_MuonColliderDet.tcl -a $aQGCWorkDir/aQGC/analyzeDelphes.py + From e0056b9ddebd8d4c90f6968003fb0bf6164bef59 Mon Sep 17 00:00:00 2001 From: John Stupak Date: Tue, 9 Nov 2021 10:11:52 -0600 Subject: [PATCH 14/22] adding gitignore --- .gitignore | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 .gitignore 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 From 714a332fbf344b72fa66dcd8bedc23f4ea6fa65b Mon Sep 17 00:00:00 2001 From: jstupak Date: Tue, 9 Nov 2021 12:00:04 -0600 Subject: [PATCH 15/22] bugfix --- install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/install.sh b/install.sh index fb9c7ca..2f2f415 100755 --- a/install.sh +++ b/install.sh @@ -58,7 +58,7 @@ if [[ $? -ne 0 ]]; then else export delphesDir=$aQGCWorkDir/delphes export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${delphesDir} - export export ROOT_INCLUDE_PATH=ROOT_INCLUDE_PATH:${delphesDir}/external/ + export ROOT_INCLUDE_PATH=ROOT_INCLUDE_PATH:${delphesDir}/external/ cd .. fi From 7f860f1cf00d85eefa8da1280514f3106b53fd9e Mon Sep 17 00:00:00 2001 From: jstupak Date: Tue, 9 Nov 2021 15:04:00 -0600 Subject: [PATCH 16/22] bugfix --- test.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test.sh b/test.sh index d6593f5..465a68f 100755 --- a/test.sh +++ b/test.sh @@ -4,5 +4,5 @@ git clone -b cleanUp https://github.com/jstupak/aQGC.git ./aQGC/install.sh source aQGC/setup.sh -./aQGC/generateEvents.sh -m $aQGCWorkDir/test.mg -d $delphesDir/cards/delphes_card_MuonColliderDet.tcl -a $aQGCWorkDir/aQGC/analyzeDelphes.py +./aQGC/generateEvents.sh -m $aQGCWorkDir/aQGC/test.mg -d $delphesDir/cards/delphes_card_MuonColliderDet.tcl -a $aQGCWorkDir/aQGC/analyzeDelphes.py From ce5d76f3a18d6c0eddbbfef7248fbe9be8d20c54 Mon Sep 17 00:00:00 2001 From: jstupak Date: Tue, 9 Nov 2021 15:30:24 -0600 Subject: [PATCH 17/22] Update README.md --- README.md | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) 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 +``` From 6f300437ffb84c1a6d4d43926e97eb979f11878d Mon Sep 17 00:00:00 2001 From: jstupak Date: Tue, 9 Nov 2021 17:00:38 -0600 Subject: [PATCH 18/22] adding script for re-running analysis --- analyzeDelphes.py | 37 ++++++++++++++++++++++++++++++++++++- runAnalysis.sh | 19 +++++++++++++++++++ 2 files changed, 55 insertions(+), 1 deletion(-) create mode 100755 runAnalysis.sh diff --git a/analyzeDelphes.py b/analyzeDelphes.py index 35eaf08..4b1cacb 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -87,7 +87,7 @@ def isBeamRemnant(p): tag=str(argv[2]) else: tag='' - outputName=(bool(tag)*(tag+'_'))+inputName.replace('.root','.hist.root') + outputName=inputName.replace('.root',(bool(tag)*('.%s'%tag))+'.hist.root') print "inputName:",inputName print "tag:",tag @@ -186,6 +186,41 @@ def isBeamRemnant(p): h_truth[22]['mult'].Fill(len(truthPhotons)) T_beamRemnants_multiplicity.Fill(len(beamRemnantMuons)) + 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[11]['pT'].Fill(truthMuons[i].PT) + h_truth[11]['p'].Fill(truthMuons[i].P4().P()) + h_truth[11]['eta'].Fill(truthMuons[i].Eta) + + for i in range(len(truthWs)): + h_truth[11]['pT'].Fill(truthWs[i].PT) + h_truth[11]['p'].Fill(truthWs[i].P4().P()) + h_truth[11]['eta'].Fill(truthWs[i].Eta) + + for i in range(len(truthZs)): + h_truth[11]['pT'].Fill(truthZs[i].PT) + h_truth[11]['p'].Fill(truthZs[i].P4().P()) + h_truth[11]['eta'].Fill(truthZs[i].Eta) + + for i in range(len(truthPhotons)): + h_truth[11]['pT'].Fill(truthPhotons[i].PT) + h_truth[11]['p'].Fill(truthPhotons[i].P4().P()) + h_truth[11]['eta'].Fill(truthPhotons[i].Eta) + + for i in range(len(truthNeutrinos)): + h_truth[11]['pT'].Fill(truthNeutrinos[i].PT) + h_truth[11]['p'].Fill(truthNeutrinos[i].P4().P()) + h_truth[11]['eta'].Fill(truthNeutrinos[i].Eta) + + 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) + truthMissingP=TLorentzVector() for nu in truthNeutrinos: truthMissingP+=nu.P4() 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 From 9edc104747254ac89433e71d1a93734b708b88b4 Mon Sep 17 00:00:00 2001 From: Connor Waits Date: Mon, 15 Nov 2021 14:47:26 -0600 Subject: [PATCH 19/22] fixed W status code bug --- analyzeDelphes.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analyzeDelphes.py b/analyzeDelphes.py index 4b1cacb..1c87dc6 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -171,8 +171,8 @@ def isBeamRemnant(p): 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,'x.Status==22 and abs(x.PID)==24') - truthZs =selector(event.Particle,'x.Status==22 and abs(x.PID)==23') + truthWs =selector(event.Particle,'abs(x.Status)==22 and abs(x.PID)==24') + truthZs =selector(event.Particle,'abs(x.Status)==22 and abs(x.PID)==23') truthPhotons =selector(event.Particle,'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)') From e4bf299c1965624fc2c1e530f5562270091abde7 Mon Sep 17 00:00:00 2001 From: Connor Waits Date: Sun, 21 Nov 2021 17:16:22 -0600 Subject: [PATCH 20/22] fixed bug so that hists are filled for correct particle species, added neutrinos to particle dictionary --- analyzeDelphes.py | 52 ++++++++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/analyzeDelphes.py b/analyzeDelphes.py index 1c87dc6..15effb9 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -9,7 +9,8 @@ 21:"gamma", #using this for all jets 22:"j", 23:"z", - 24:"w" + 24:"w", + 12:"nu" } #OBJECT SELECTION @@ -105,7 +106,7 @@ def isBeamRemnant(p): #declares truth-level pT, p, eta, and multiplicity histograms for e,mu,W,Z,gamma h_truth={} - for p in [11,13,21,22,23,24]: + 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], 200, 0, bin_range) h_truth[p]['p'] =TH1F('T_%s_p'%particle[p],';Truth %s p [GeV];Events'%particle[p], 200, 0, bin_range) @@ -162,9 +163,9 @@ def isBeamRemnant(p): Test_beamRemnants_pT = TH1F('Test_beamRemnants_Pt', ';pT [GeV];Events', 30, 0, 3) #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + m=0 for event in f.Delphes: - + m+=1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #truth level #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -172,8 +173,8 @@ def isBeamRemnant(p): 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)==22 and abs(x.PID)==23') - truthPhotons =selector(event.Particle,'x.Status==1 and abs(x.PID)==22') + 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)') @@ -184,37 +185,42 @@ def isBeamRemnant(p): 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)) - + + if (m<100): + print m + for i in range(len(truthWs)): + print truthWs[i].PT 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[11]['pT'].Fill(truthMuons[i].PT) - h_truth[11]['p'].Fill(truthMuons[i].P4().P()) - h_truth[11]['eta'].Fill(truthMuons[i].Eta) + 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[11]['pT'].Fill(truthWs[i].PT) - h_truth[11]['p'].Fill(truthWs[i].P4().P()) - h_truth[11]['eta'].Fill(truthWs[i].Eta) - + 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 i in range(len(truthZs)): - h_truth[11]['pT'].Fill(truthZs[i].PT) - h_truth[11]['p'].Fill(truthZs[i].P4().P()) - h_truth[11]['eta'].Fill(truthZs[i].Eta) + 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[11]['pT'].Fill(truthPhotons[i].PT) - h_truth[11]['p'].Fill(truthPhotons[i].P4().P()) - h_truth[11]['eta'].Fill(truthPhotons[i].Eta) + h_truth[21]['pT'].Fill(truthPhotons[i].PT) + h_truth[21]['p'].Fill(truthPhotons[i].P4().P()) + h_truth[21]['eta'].Fill(truthPhotons[i].Eta) for i in range(len(truthNeutrinos)): - h_truth[11]['pT'].Fill(truthNeutrinos[i].PT) - h_truth[11]['p'].Fill(truthNeutrinos[i].P4().P()) - h_truth[11]['eta'].Fill(truthNeutrinos[i].Eta) + 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) for i in range(len(beamRemnantMuons)): T_beamRemnants_pT.Fill(beamRemnantMuons[i].PT) From 8dcd7a61f810db9aabe281eaca5618a8acbb1001 Mon Sep 17 00:00:00 2001 From: Connor Waits Date: Sun, 21 Nov 2021 17:19:54 -0600 Subject: [PATCH 21/22] took out print statements from debugging --- analyzeDelphes.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/analyzeDelphes.py b/analyzeDelphes.py index 15effb9..416f3bd 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -163,9 +163,9 @@ def isBeamRemnant(p): Test_beamRemnants_pT = TH1F('Test_beamRemnants_Pt', ';pT [GeV];Events', 30, 0, 3) #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - m=0 + for event in f.Delphes: - m+=1 + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #truth level #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -187,11 +187,7 @@ def isBeamRemnant(p): h_truth[22]['mult'].Fill(len(truthPhotons)) h_truth[12]['mult'].Fill(len(truthNeutrinos)) T_beamRemnants_multiplicity.Fill(len(beamRemnantMuons)) - - if (m<100): - print m - for i in range(len(truthWs)): - print truthWs[i].PT + for i in range(len(truthElectrons)): h_truth[11]['pT'].Fill(truthElectrons[i].PT) h_truth[11]['p'].Fill(truthElectrons[i].P4().P()) From 10ca993a0fa50783a34a970050603252cab8bfb8 Mon Sep 17 00:00:00 2001 From: Connor Waits Date: Mon, 7 Mar 2022 10:48:32 -0600 Subject: [PATCH 22/22] added cutflow histogram --- analyzeDelphes.py | 235 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 199 insertions(+), 36 deletions(-) diff --git a/analyzeDelphes.py b/analyzeDelphes.py index 416f3bd..f9848ae 100644 --- a/analyzeDelphes.py +++ b/analyzeDelphes.py @@ -1,13 +1,13 @@ maxEvents=9E9 DEBUG=False -jetType='KTjet' +jetType='VLCjetR15N2' particle={ 11:"e", 13:"mu", - 21:"gamma", #using this for all jets - 22:"j", + 22:"gamma", #using this for all jets + 21:"j", 23:"z", 24:"w", 12:"nu" @@ -38,6 +38,8 @@ from sys import argv #--------------------------------------------------------- from ROOT import * +TH1.SetDefaultSumw2() +import math gSystem.Load("libDelphes") gStyle.SetOptStat(0) @@ -80,6 +82,9 @@ def isBeamRemnant(p): return parents.Status==4 #--------------------------------------------------------- +def sortFunc(x): + return (x.PT)*math.cosh(x.Eta) +#--------------------------------------------------------- if __name__=='__main__': @@ -102,18 +107,21 @@ def isBeamRemnant(p): sqrtS+=t.Particle[0].E sqrtS+=t.Particle[1].E #sets upper limit for bin range - bin_range = sqrtS/10 + bin_range = sqrtS/2 #declares truth-level pT, p, eta, and multiplicity histograms for e,mu,W,Z,gamma 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], 200, 0, bin_range) - h_truth[p]['p'] =TH1F('T_%s_p'%particle[p],';Truth %s p [GeV];Events'%particle[p], 200, 0, bin_range) + 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) @@ -124,13 +132,17 @@ def isBeamRemnant(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 pT [GeV];Events'%particle[p], 200, 0, bin_range) - h_reco[p]['eta'][i]=TH1F('R_%s_eta_%s'%(particle[p],i),';%s pT [GeV];Events'%particle[p], 20, -4, 4) - + 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_eta = TH1F('ee_Eta', ';Eta;Events', 20, -4, 4) @@ -140,9 +152,11 @@ def isBeamRemnant(p): 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_eta = TH1F('emu_Eta', ';Eta;Events', 20, -4, 4) @@ -154,18 +168,56 @@ def isBeamRemnant(p): 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_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', 20, -10, 10) + 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) + 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 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -179,7 +231,54 @@ def isBeamRemnant(p): 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)) @@ -187,12 +286,20 @@ def isBeamRemnant(p): 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()) @@ -209,20 +316,30 @@ def isBeamRemnant(p): h_truth[23]['eta'].Fill(truthZs[i].Eta) for i in range(len(truthPhotons)): - h_truth[21]['pT'].Fill(truthPhotons[i].PT) - h_truth[21]['p'].Fill(truthPhotons[i].P4().P()) - h_truth[21]['eta'].Fill(truthPhotons[i].Eta) + 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() @@ -244,19 +361,11 @@ def isBeamRemnant(p): P4_f = TLorentzVector() - 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])) - #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 - - #sort object lists - for collection in [electrons, muons, photons, jets]: collection.sort(key=operator.attrgetter('PT'),reverse=True) h_reco[11]['mult'].Fill(len(electrons)) for i in range(len(electrons)): @@ -264,13 +373,19 @@ def isBeamRemnant(p): 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)) + 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) @@ -281,6 +396,27 @@ def isBeamRemnant(p): 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() @@ -292,17 +428,37 @@ def isBeamRemnant(p): h_reco[22]['p'][i].Fill(photons[i].P4().P()) h_reco[22]['eta'][i].Fill(photons[i].Eta) - h_reco[11]['mult'].Fill(len(jets)) + 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()) #------------------------------------------------------------------------------------------------- @@ -337,7 +493,7 @@ def isBeamRemnant(p): R_ee_eta.Fill(theZ.Eta()) R_ee_deltaEta.Fill(abs(l1.Eta - l2.Eta)) else: - R_mumu_mass.Fill(theZ.Mass()) + #R_mumu_mass.Fill(theZ.Mass()) R_mumu_pT.Fill(theZ.Pt()) R_mumu_eta.Fill(theZ.Eta()) R_mumu_deltaEta.Fill(abs(l1.Eta - l2.Eta)) @@ -428,4 +584,11 @@ def isBeamRemnant(p): R_mumu_multiplicity.Fill(mumu) R_emu_multiplicity.Fill(emu) """ + T_eta_vs_WWmass.ProfileY('Y_axis_profile') + + print CutFlow.GetBinContent(1) + print CutFlow.GetBinContent(2) + print CutFlow.GetBinContent(3) + print CutFlow.GetBinContent(4) + print CutFlow.GetBinContent(5) output.Write()