background_level_CSC_CSC_debug.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c839451f",
   "metadata": {},
   "source": [
    "## Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "58b6e643",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.24/06\n"
     ]
    }
   ],
   "source": [
    "import ROOT as rt\n",
    "import csv\n",
    "import re\n",
    "import sys\n",
    "import os\n",
    "import collections\n",
    "from collections import OrderedDict\n",
    "import uproot\n",
    "import pandas as pd\n",
    "from root_numpy import array2tree\n",
    "import scipy\n",
    "import awkward\n",
    "import numpy as np\n",
    "import time\n",
    "import math\n",
    "from scipy.stats import norm\n",
    "#import PrettyTable\n",
    "\n",
    "sys.path.append('/afs/desy.de/user/l/lbenato/cms-lpc-llp_repo/run3_muon_system_analysis/lib/')\n",
    "from histo_utilities import create_TH1D, create_TH2D, std_color_list, create_TGraph, make_ratio_plot, addOverflow\n",
    "from helper_functions import deltaR, deltaPhi, drawCMS\n",
    "#from helper import *\n",
    "\n",
    "rt.gStyle.SetOptStat(0)\n",
    "rt.gROOT.SetBatch(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d97fb5",
   "metadata": {},
   "source": [
    "## Load ntuples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0c47f19",
   "metadata": {},
   "source": [
    "## Load variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "id": "368e657b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1118331.3125610352\n"
     ]
    }
   ],
   "source": [
    "fpath =OrderedDict()\n",
    "tree = OrderedDict()\n",
    "\n",
    "start_t = time.time()\n",
    "data_year = 2022\n",
    "\n",
    "tag  = 'V1p19'\n",
    "vers = 'v6'\n",
    "#v1 includes events without selections\n",
    "#v3 includes events with at least 2 clusters\n",
    "#v4: no event level cuts, clusters made with DBSCAN, HLT time definition\n",
    "#v5: no event level cuts, clusters made with CA\n",
    "\n",
    "data_path = '/nfs/dust/cms/group/cms-llp/muon_system_run3/'+tag+'/Data'+str(data_year)+'/'+vers+'/'\n",
    "lumi = 23.02*1000\n",
    "xsec = 48.\n",
    "\n",
    "if data_year == 2022:\n",
    "    #fpath['data'] = data_path + 'DisplacedJet-EXOCSCCluster_Run2022EFG-PromptReco-v1_goodLumi.root'\n",
    "    fpath['signal'] = data_path + 'ggH_HToSSTobbbb_MH-125_MS-15_CTau1000_13p6TeV_1pb_weighted.root'\n",
    "\n",
    "NEvents = {}\n",
    "for k,v in fpath.items():\n",
    "    root_dir = uproot.open(v)\n",
    "    tree[k] = root_dir['MuonSystem']\n",
    "    NEvents[k] = root_dir['NEvents'][1]\n",
    "    a = tree[k][\"weight\"].array()\n",
    "    #print(a)\n",
    "    if k=='signal':\n",
    "        signal_gen_weight = tree[k][\"weight\"].array()\n",
    "\n",
    "signal_gen_yield = np.sum(signal_gen_weight)*lumi\n",
    "print(signal_gen_yield)\n",
    "\n",
    "for k, T in tree.items():\n",
    "    branch_names = T.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "id": "899a72f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "JET_PT_CUT = 10.0\n",
    "MUON_PT_CUT = 20.0\n",
    "N_RECHIT_CUT = 90\n",
    "jetPt_cut = 50\n",
    "tightid = False\n",
    "ring_cut = 50\n",
    "\n",
    "#cut_based = True\n",
    "#cut_based_version = 'v4'\n",
    "\n",
    "intime = True\n",
    "DPHI_CUT = 1\n",
    "\n",
    "\n",
    "gLLP_csc = {}\n",
    "\n",
    "\n",
    "selections_cluster = {}\n",
    "sel_cluster = {}\n",
    "sel_jetveto_csc = {}\n",
    "sel_muonveto_csc = {}\n",
    "sel_jetveto_dt = {}\n",
    "sel_muonveto_dt = {}\n",
    "met_trigger = {}\n",
    "met = {}\n",
    "gLLP_beta = {}\n",
    "jetPt = {}\n",
    "jetPhi = {}\n",
    "metPhi = {}\n",
    "angle ={}\n",
    "nLeptons = {}\n",
    "genJetPt = {}\n",
    "genJetPhi = {}\n",
    "genMet = {}\n",
    "genMetPhi = {}\n",
    "pileupWeight = {}\n",
    "gLLP_ctau = {}\n",
    "npv = {}\n",
    "nRechitClusters = {}\n",
    "nJets = {}\n",
    "nJets_50gev = {}\n",
    "nCscRings = {}\n",
    "nDtRings = {}\n",
    "\n",
    "weight = {}\n",
    "all_weight = {}\n",
    "jetMet_dPhiMin = {}\n",
    "dphiMet_cluster = {}\n",
    "nRechits_sr = {}\n",
    "jetMet_dPhiMin30_sr = {}\n",
    "cscClusterTimeSpread = {}\n",
    "bdt_score = {}\n",
    "nCscChambers = {}\n",
    "\n",
    "sel_ev = {}\n",
    "sel_ev_post = {}\n",
    "hlt_sel_ev = {}\n",
    "cluster_index = ''\n",
    "nRings = {}\n",
    "\n",
    "#CSC cluster variables\n",
    "cscClusterMuonVetoPt = {}\n",
    "cscClusterJetVetoPt = {}\n",
    "cscClusterPhi = {}\n",
    "cscClusterEta = {}\n",
    "cscClusterEta2 = {}\n",
    "cscClusterEta3 = {}\n",
    "cscClusterR = {}\n",
    "cscClusterZ = {}\n",
    "cscClusterSize = {}\n",
    "cscClusterSize2 = {}\n",
    "cscClusterSize3 = {}\n",
    "cscClusterTimeSpread = {}\n",
    "cscClusterTime = {}\n",
    "cscClusterTime2 = {}\n",
    "cscClusterTime3 = {}\n",
    "cscClusterMet_dPhi = {}\n",
    "\n",
    "cscClusterNStation = {}\n",
    "cscClusterDphi = {}\n",
    "nCscClusters = {}\n",
    "\n",
    "cscClusterNRechitMinus11 = {}\n",
    "cscClusterNRechitPlus11 = {}\n",
    "cscClusterNRechitMinus12 = {}\n",
    "cscClusterNRechitPlus12 = {}\n",
    "\n",
    "\n",
    "#DT cluster variables\n",
    "dtClusterMuonVetoPt = {}\n",
    "dtClusterJetVetoPt = {}\n",
    "dtClusterTime = {}\n",
    "dtClusterPhi = {}\n",
    "dtClusterEta = {}\n",
    "dtClusterSize = {}\n",
    "dtClusterDphi = {}\n",
    "nDtClusters = {}\n",
    "dtClusterNStation = {}\n",
    "dtClusterMaxStation = {}\n",
    "dtClusterNHitStation1 = {}\n",
    "dtCluster_match_MB1hits_0p4 = {}\n",
    "dtCluster_match_MB1hits_0p5 = {}\n",
    "dtClusterMaxStationRatio = {}\n",
    "dtClusterMet_dPhi = {}\n",
    "\n",
    "deltaEta = {}        \n",
    "deltaRCluster = {}\n",
    "\n",
    "metPhi = {}\n",
    "evtNum = {}\n",
    "runNum = {}\n",
    "lumiNum = {}\n",
    "\n",
    "nMe11 = {}\n",
    "\n",
    "tree_keys = []\n",
    "nClusterRatio = {}\n",
    "\n",
    "ME11_veto = {}\n",
    "ME12_veto = {}\n",
    "MB1_veto = {}\n",
    "\n",
    "#L1 plateau\n",
    "first_in_ME11 = {}\n",
    "first_in_ME12 = {}\n",
    "first_in_ME13 = {}\n",
    "first_in_ME21 = {}\n",
    "first_in_ME22 = {}\n",
    "first_in_ME31 = {}\n",
    "first_in_ME32 = {}\n",
    "first_in_ME41 = {}\n",
    "first_in_ME42 = {}\n",
    "\n",
    "\n",
    "first_in_plateau_ME11 = {}\n",
    "first_in_plateau_ME12 = {}\n",
    "first_in_plateau_ME13 = {}\n",
    "first_in_plateau_ME21 = {}\n",
    "first_in_plateau_ME22 = {}\n",
    "first_in_plateau_ME31 = {}\n",
    "first_in_plateau_ME32 = {}\n",
    "first_in_plateau_ME41 = {}\n",
    "first_in_plateau_ME42 = {}\n",
    "first_in_plateau = {}\n",
    "\n",
    "#Angular variables\n",
    "dPhi_csc_csc = {}\n",
    "dEta_csc_csc = {}\n",
    "dR_csc_csc = {}\n",
    "dt_csc_csc = {}\n",
    "dPhi_csc_dt = {}\n",
    "dEta_csc_dt = {}\n",
    "dR_csc_dt = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8091520",
   "metadata": {},
   "source": [
    "## Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "id": "6f8fbffd",
   "metadata": {},
   "outputs": [],
   "source": [
    "category = 0#-1#0\n",
    "category = 1\n",
    "#-1: at least 1csc\n",
    "#0: exactly 1 csc\n",
    "#1: 2csc,\n",
    "#2: 1csc+1dt, \n",
    "\n",
    "do_jet_veto = False\n",
    "do_inverted_jet_veto = False\n",
    "if (do_jet_veto and do_inverted_jet_veto):\n",
    "    print(\"Please check jet veto! Aborting...\")\n",
    "    exit()\n",
    "\n",
    "do_muon_veto = False\n",
    "do_inverted_muon_veto = False\n",
    "if (do_muon_veto and do_inverted_muon_veto):\n",
    "    print(\"Please check muon veto! Aborting...\")\n",
    "    exit()\n",
    "\n",
    "do_csc_time_spread_cut = False\n",
    "do_csc_in_time_cut = False\n",
    "do_ME_veto = False\n",
    "do_MB_veto = True\n",
    "do_HLT_eta_cut = False#True\n",
    "do_L1_plateau = True\n",
    "\n",
    "plot_string = \"\"\n",
    "#plot_string = \"_jet_veto\"\n",
    "#plot_string = \"_csc_time_spread_cut\"\n",
    "#plot_string = \"_csc_in_time_cut\"\n",
    "#plot_string = \"\"\n",
    "#plot_string = \"_ME_veto\"\n",
    "#plot_string = \"_ME_veto_HLT_eta_cut\"\n",
    "#plot_string = \"_ME_veto_HLT_eta_cut_at_least_1_cluster\"\n",
    "#plot_string = \"_ME_veto_HLT_eta_cut_in_time_cut_at_least_1_cluster\"\n",
    "#plot_string+= \"_exactly_1_csc_cluster\"\n",
    "\n",
    "#plot_string+= \"_exactly_2_csc_clusters_ME_veto\"\n",
    "#plot_string+= \"_exactly_2_csc_clusters_ME_veto_time_spread_cut_jet_muon_veto_2nd\"\n",
    "#plot_string+= \"_exactly_2_csc_clusters_time_spread_cut_jet_muon_veto_2nd\"\n",
    "#plot_string+= \"_2_csc_clusters_ME_veto\"\n",
    "#plot_string+= \"_at_least_1_csc_cluster_ME_veto\"\n",
    "#plot_string+= \"_csc_dt_cluster_ME_veto\"\n",
    "#plot_string+= \"_csc_dt_cluster_ME_MB_veto\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "bb3672eb",
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "index 0 is out of bounds for jagged min size 0",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-200-f00803e3368b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     36\u001b[0m     \u001b[0;31m#ME11 and ME12 veto\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m     ME11_veto[k] = np.logical_and(\n\u001b[0;32m---> 38\u001b[0;31m         \u001b[0mT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'cscRechitCluster'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mcluster_index\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'NRechitChamberMinus11'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m0\u001b[0m  \u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m         \u001b[0mT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'cscRechitCluster'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mcluster_index\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'NRechitChamberPlus11'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m     )  \n",
      "\u001b[0;32m/nfs/dust/cms/user/lbenato/anaconda3/envs/run3/lib/python3.7/site-packages/awkward/array/jagged.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, where)\u001b[0m\n\u001b[1;32m    624\u001b[0m                         \u001b[0mhead\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcounts\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mhead\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    625\u001b[0m                     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbitwise_and\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0mhead\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhead\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mcounts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 626\u001b[0;31m                         \u001b[0;32mraise\u001b[0m \u001b[0mIndexError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"index {0} is out of bounds for jagged min size {1}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moriginal_head\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcounts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    627\u001b[0m                     \u001b[0mnode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_content\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_starts\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mhead\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    628\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: index 0 is out of bounds for jagged min size 0"
     ]
    }
   ],
   "source": [
    "for k, T in tree.items():\n",
    "    tree_keys.append(k)\n",
    "    #if not k == 'data' and not k[-4:] == '1000':continue\n",
    "########### SELECTION: CLUSTERS ############\n",
    "    jet_veto_csc = T.array('cscRechitCluster' + cluster_index + 'JetVetoPt')<30\n",
    "    muon_veto_csc = np.logical_not(np.logical_and(T.array('cscRechitClusterMuonVetoPt') >= 30, T.array('cscRechitClusterMuonVetoGlobal')))\n",
    "    \n",
    "    sel_csccluster = T.array('cscRechitClusterSize') >= 0\n",
    "    \n",
    "    if do_jet_veto:\n",
    "        sel_csccluster = jet_veto_csc\n",
    "    if do_inverted_jet_veto:\n",
    "        sel_csccluster = np.logical_not(jet_veto_csc)       \n",
    "\n",
    "    if do_muon_veto:\n",
    "        sel_csccluster = np.logical_and(sel_csccluster, muon_veto_csc)\n",
    "    if do_inverted_muon_veto:\n",
    "        sel_csccluster = np.logical_and(sel_csccluster, np.logical_not(muon_veto_csc))    \n",
    "          \n",
    "    #Missing branch:\n",
    "    #sel_csccluster = np.logical_and(sel_csccluster, np.abs(T.array('cscRechitCluster' + cluster_index + 'MetEENoise_dPhi'))<1.2)\n",
    "\n",
    "    #Missing branch:\n",
    "    #sel_csccluster = np.logical_and(sel_csccluster, T.array('cscRechitCluster' + cluster_index + 'Me11Ratio')<1)\n",
    "\n",
    "    #Time spread cut\n",
    "    if do_csc_time_spread_cut:\n",
    "        sel_csccluster = np.logical_and(sel_csccluster, T.array('cscRechitCluster' + cluster_index + 'TimeSpreadWeightedAll')<20)\n",
    "    \n",
    "    #Time cut\n",
    "    if do_csc_in_time_cut:\n",
    "        sel_csccluster = np.logical_and(sel_csccluster, np.logical_and(\n",
    "            T.array('cscRechitCluster' + cluster_index + 'TimeWeighted')< 12.5, \n",
    "            T.array('cscRechitCluster' + cluster_index + 'TimeWeighted') > -5))\n",
    "\n",
    "    #ME11 and ME12 veto\n",
    "    ME11_veto[k] = np.logical_and(\n",
    "        T.array('cscRechitCluster' + cluster_index + 'NRechitChamberMinus11')[:,0]==0  ,\n",
    "        T.array('cscRechitCluster' + cluster_index + 'NRechitChamberPlus11')[:,0]==0\n",
    "    )  \n",
    "\n",
    "    ME12_veto[k] = np.logical_and(\n",
    "        T.array('cscRechitCluster' + cluster_index + 'NRechitChamberMinus12')[:,0]==0  ,\n",
    "        T.array('cscRechitCluster' + cluster_index + 'NRechitChamberPlus12')[:,0]==0\n",
    "    )      \n",
    "\n",
    "    print(ME11_veto[k].sum())\n",
    "    print( ME11_veto[k][ ME11_veto[k].sum()>=1][:,0] )\n",
    "    print( ME11_veto[k][ ME11_veto[k].sum()>=1][:] )\n",
    "       \n",
    "    \n",
    "    if do_ME_veto:\n",
    "        sel_csccluster = np.logical_and(sel_csccluster,\n",
    "                                    np.logical_and(ME11_veto[k],ME12_veto[k])\n",
    "                                   ) \n",
    "    if do_HLT_eta_cut:\n",
    "        sel_csccluster = np.logical_and(sel_csccluster,sel_trgCluster)\n",
    "        \n",
    "        \n",
    "    sel_dtcluster = T.array('dtRechitClusterSize') >= 0\n",
    "    #Missing branch:\n",
    "    #sel_dtcluster = np.abs(T.array('dtRechitClusterMetEENoise_dPhi')) < 1\n",
    "    jet_veto_dt = np.abs(T.array('dtRechitClusterJetVetoPt')) < 50\n",
    "    muon_veto_dt = np.logical_not(np.logical_and(T.array('dtRechitClusterMuonVetoPt') >= 10, T.array('dtRechitClusterMuonVetoLooseId')))\n",
    "\n",
    "    if do_muon_veto:\n",
    "        sel_dtcluster = np.logical_and(sel_dtcluster, muon_veto_dt)\n",
    "    if do_jet_veto:\n",
    "        sel_dtcluster = np.logical_and(sel_dtcluster, jet_veto_dt)\n",
    "\n",
    "    MB1_veto[k] = T.array('dtRechitClusterNHitStation1') == 0\n",
    "    if do_MB_veto:\n",
    "        sel_dtcluster = np.logical_and(sel_dtcluster, MB1_veto[k])\n",
    "        \n",
    "    \n",
    "########### SELECTION: JETS ############\n",
    "    \n",
    "    sel_jet = np.logical_and(T.array('jetPt') > jetPt_cut, np.abs(T.array('jetEta')) < 2.4 )\n",
    "\n",
    "########### SELECTION: EVENTS ############\n",
    "    hlt = T['HLTDecision'].array()    \n",
    "    HLT_CscCluster_Loose = hlt[:,566]\n",
    "    HLT_CscCluster_Medium = hlt[:,567]\n",
    "    HLT_CscCluster_Tight = hlt[:,568]\n",
    "    HLT_L1CSCShower_DTCluster50 = hlt[:,569]\n",
    "    HLT_L1CSCShower_DTCluster75 = hlt[:,570]\n",
    "\n",
    "    #Or of the triggers: starting with the CSC ones\n",
    "    #hlt_sel_ev[k] = np.logical_or(HLT_CscCluster_Loose,np.logical_or(HLT_CscCluster_Medium,HLT_CscCluster_Tight))\n",
    "    if category==2:\n",
    "        hlt_sel_ev[k] = HLT_L1CSCShower_DTCluster50\n",
    "    else:\n",
    "        hlt_sel_ev[k] = HLT_CscCluster_Loose\n",
    "     \n",
    "    \n",
    "    #sel_ev[k]  = np.logical_and(sel_ev[k], T.array('category') == 0)\n",
    "    #sel_ev[k] = np.logical_and(sel_ev[k], (T.array('nDtRings')+T.array('nCscRings'))<10)\n",
    "    sel_ev[k] = T.array('nLeptons') == 0\n",
    "    \n",
    "    sel_ev[k] = (sel_ev[k] & (T.array('evtNum')==9430 ))\n",
    "    \n",
    "    #If trigger:\n",
    "    #We don't have correct trigger emulation in signal at the moment\n",
    "    if k == 'data': sel_ev[k] = np.logical_and(sel_ev[k],hlt_sel_ev[k])\n",
    "    if k == 'data': sel_ev[k] = np.logical_and(sel_ev[k], T.array('runNum')>=360019)\n",
    "\n",
    "    #At least one pre-selected csc cluster, always required\n",
    "    sel_ev[k] = np.logical_and(sel_ev[k],sel_csccluster.sum()>=1)\n",
    "    print(\"Looking for survivors pre category\")\n",
    "    print(sel_ev[k][ sel_ev[k] ])\n",
    "    print(\"n clusters\")\n",
    "    print(sel_csccluster[ sel_ev[k] ].sum())\n",
    "    \n",
    "    \n",
    "    if category == 0:\n",
    "        sel_ev[k]  = np.logical_and(sel_ev[k],sel_csccluster.sum() == 1)\n",
    "        hlt_sel_ev[k] = np.logical_and(sel_ev[k],hlt_sel_ev[k])\n",
    "    elif category == 1:\n",
    "        sel_ev[k]  = np.logical_and(sel_ev[k],sel_csccluster.sum() >= 2)\n",
    "        hlt_sel_ev[k] = np.logical_and(sel_ev[k],hlt_sel_ev[k])\n",
    "    elif category == 2:      \n",
    "        sel_ev[k]  = np.logical_and(sel_ev[k],sel_csccluster.sum() == 1)\n",
    "        sel_ev[k]  = np.logical_and(sel_ev[k],sel_dtcluster.sum() == 1)\n",
    "        hlt_sel_ev[k] = np.logical_and(sel_ev[k],hlt_sel_ev[k])\n",
    "    elif category == -1:      \n",
    "        sel_ev[k]  = np.logical_and(sel_ev[k],sel_csccluster.sum() >= 1)\n",
    "        hlt_sel_ev[k] = np.logical_and(sel_ev[k],hlt_sel_ev[k])\n",
    "\n",
    "    print(\"Looking for survivors after category\")\n",
    "    print(sel_ev[k][ sel_ev[k] ])\n",
    "    print(\"n clusters\")\n",
    "    print(sel_csccluster[ sel_ev[k] ].sum())\n",
    "########### BRANCHES ############\n",
    "\n",
    "    \n",
    "    dtClusterMuonVetoPt[k] = T.array('dtRechitClusterMuonVetoPt')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterJetVetoPt[k] = T.array('dtRechitClusterJetVetoPt')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterPhi[k] = T.array('dtRechitClusterPhi')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterEta[k] = T.array('dtRechitClusterEta')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterSize[k] =  T.array('dtRechitClusterSize')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterTime[k] =  T.array('dtRechitCluster_match_RPCBx_dPhi0p5')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterNStation[k] =  T.array('dtRechitClusterNStation10')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterNHitStation1[k] = T.array('dtRechitClusterNHitStation1')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtCluster_match_MB1hits_0p4[k] = T.array('dtRechitCluster_match_MB1hits_0p4')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtCluster_match_MB1hits_0p5[k] = T.array('dtRechitCluster_match_MB1hits_0p5')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterMet_dPhi[k] =  T.array('dtRechitCluster' + cluster_index + 'Met_dPhi')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    dtClusterMaxStationRatio[k] = T.array('dtRechitCluster' + cluster_index + 'MaxStationRatio')[sel_dtcluster][sel_ev[k]][:,:]\n",
    "        \n",
    "    cscClusterMuonVetoPt[k] = T.array('cscRechitCluster' + cluster_index + 'MuonVetoPt')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterJetVetoPt[k] = T.array('cscRechitCluster' + cluster_index + 'JetVetoPt')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterPhi[k] = T.array('cscRechitCluster' + cluster_index + 'Phi')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterEta[k] = T.array('cscRechitCluster' + cluster_index + 'Eta')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterR[k] = np.sqrt(  T.array('cscRechitCluster' + cluster_index + 'X')[sel_csccluster][sel_ev[k]][:,:]**2 + T.array('cscRechitCluster' + cluster_index + 'Y')[sel_csccluster][sel_ev[k]][:,:]**2 )\n",
    "    cscClusterZ[k] = T.array('cscRechitCluster' + cluster_index + 'Z')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterSize[k] =  T.array('cscRechitCluster' + cluster_index + 'Size')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterTime[k] =  T.array('cscRechitCluster' + cluster_index + 'Time')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterTimeSpread[k] =  T.array('cscRechitCluster' + cluster_index + 'TimeSpread')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterNStation[k] =  T.array('cscRechitCluster' + cluster_index + 'NStation10')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterMet_dPhi[k] =  T.array('cscRechitCluster' + cluster_index + 'Met_dPhi')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    nCscClusters[k] = sel_csccluster.sum()[sel_ev[k]]\n",
    "    nDtClusters[k] = sel_dtcluster.sum()[sel_ev[k]]\n",
    "    \n",
    "    cscClusterNRechitMinus11[k] = T.array('cscRechitCluster' + cluster_index + 'NRechitChamberMinus11')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterNRechitPlus11[k] = T.array('cscRechitCluster' + cluster_index + 'NRechitChamberPlus11')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterNRechitMinus12[k] = T.array('cscRechitCluster' + cluster_index + 'NRechitChamberMinus12')[sel_csccluster][sel_ev[k]][:,:]\n",
    "    cscClusterNRechitPlus12[k] = T.array('cscRechitCluster' + cluster_index + 'NRechitChamberPlus12')[sel_csccluster][sel_ev[k]][:,:]\n",
    "\n",
    "      \n",
    "    sel_jetveto_csc[k] = jet_veto_csc[sel_csccluster][sel_ev[k]][:,:]\n",
    "    sel_muonveto_csc[k] = muon_veto_csc[sel_csccluster][sel_ev[k]][:,:]\n",
    "       \n",
    "    sel_jetveto_dt[k] = jet_veto_dt[sel_dtcluster][sel_ev[k]][:,:]\n",
    "    sel_muonveto_dt[k] = muon_veto_dt[sel_dtcluster][sel_ev[k]][:,:]        \n",
    "\n",
    "    dtClusterMaxStation[k] = T.array('dtRechitClusterMaxStation')[sel_dtcluster][sel_ev[k]]\n",
    "       \n",
    "    hlt_sel_ev[k] = hlt_sel_ev[k][ sel_ev[k] ]\n",
    "    weight[k] = T.array('weight')[ sel_ev[k] ]\n",
    "    evtNum[k] = T.array('evtNum')[ sel_ev[k] ]\n",
    "    sel_ev[k] = sel_ev[k][ sel_ev[k] ]\n",
    "    \n",
    "    all_weight[k] = T.array('weight')\n",
    "\n",
    "    nClusterRatio[k] = np.divide(nCscClusters[k].astype(float),nDtClusters[k].astype(float),out=np.zeros_like(nDtClusters[k].astype(float)),where=(nDtClusters[k].astype(float) != 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2afacb34",
   "metadata": {},
   "source": [
    "* Study delta phi with met\n",
    "* Apply at least one csc cluster requirement\n",
    "* Apply L1 plateau"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38016838",
   "metadata": {},
   "source": [
    "# L1 plateau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "id": "7c541ab6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ True]\n",
      "[2]\n"
     ]
    }
   ],
   "source": [
    "first_in_tr1 = {}\n",
    "first_in_tr2 = {}\n",
    "first_in_tr3 = {}\n",
    "\n",
    "sel_ev_pl = {}\n",
    "\n",
    "for k in tree_keys:\n",
    "    first_in_ME11[k] = (cscClusterR[k][:,0]>100)&(cscClusterR[k][:,0]<275) &(abs(cscClusterZ[k][:,0])>580)&(abs(cscClusterZ[k][:,0])<632) \n",
    "    first_in_ME12[k] = (cscClusterR[k][:,0]>275)&(cscClusterR[k][:,0]<465) &(abs(cscClusterZ[k][:,0])>668)&(abs(cscClusterZ[k][:,0])<724)\n",
    "    first_in_ME13[k] = (cscClusterR[k][:,0]>505)&(cscClusterR[k][:,0]<700) &(abs(cscClusterZ[k][:,0])>668)&(abs(cscClusterZ[k][:,0])<724)\n",
    "\n",
    "    first_in_ME21[k] = (cscClusterR[k][:,0]>139)&(cscClusterR[k][:,0]<345) &(abs(cscClusterZ[k][:,0])>789)&(abs(cscClusterZ[k][:,0])<850)\n",
    "    first_in_ME22[k] = (cscClusterR[k][:,0]>357)&(cscClusterR[k][:,0]<700) &(abs(cscClusterZ[k][:,0])>791)&(abs(cscClusterZ[k][:,0])<850)\n",
    "\n",
    "    first_in_ME31[k] = (cscClusterR[k][:,0]>160)&(cscClusterR[k][:,0]<345) &(abs(cscClusterZ[k][:,0])>915)&(abs(cscClusterZ[k][:,0])<970)\n",
    "    first_in_ME32[k] = (cscClusterR[k][:,0]>357)&(cscClusterR[k][:,0]<700) &(abs(cscClusterZ[k][:,0])>911)&(abs(cscClusterZ[k][:,0])<970)\n",
    "\n",
    "    first_in_ME41[k] = (cscClusterR[k][:,0]>178)&(cscClusterR[k][:,0]<345) &(abs(cscClusterZ[k][:,0])>1002)&(abs(cscClusterZ[k][:,0])<1063)\n",
    "    first_in_ME42[k] = (cscClusterR[k][:,0]>357)&(cscClusterR[k][:,0]<700) &(abs(cscClusterZ[k][:,0])>1002)&(abs(cscClusterZ[k][:,0])<1063)\n",
    "    \n",
    "    first_in_plateau_ME11[k] = first_in_ME11[k] & (cscClusterSize[k][:,0]>=500)\n",
    "    first_in_plateau_ME21[k] = first_in_ME21[k] & (cscClusterSize[k][:,0]>=500)\n",
    "    first_in_plateau_ME31[k] = first_in_ME31[k] & (cscClusterSize[k][:,0]>=500)\n",
    "    first_in_plateau_ME41[k] = first_in_ME41[k] & (cscClusterSize[k][:,0]>=500)\n",
    "\n",
    "    first_in_plateau_ME12[k] = first_in_ME12[k] & (cscClusterSize[k][:,0]>=200)\n",
    "    first_in_plateau_ME13[k] = first_in_ME13[k] & (cscClusterSize[k][:,0]>=200)\n",
    "    first_in_plateau_ME22[k] = first_in_ME22[k] & (cscClusterSize[k][:,0]>=200)\n",
    "    first_in_plateau_ME32[k] = first_in_ME32[k] & (cscClusterSize[k][:,0]>=200)\n",
    "    first_in_plateau_ME42[k] = first_in_ME42[k] & (cscClusterSize[k][:,0]>=200)    \n",
    "    \n",
    "    first_in_plateau[k] = first_in_plateau_ME11[k] | first_in_plateau_ME12[k] | first_in_plateau_ME13[k] | first_in_plateau_ME21[k] | first_in_plateau_ME22[k] | first_in_plateau_ME31[k] | first_in_plateau_ME32[k] | first_in_plateau_ME41[k] | first_in_plateau_ME42[k]\n",
    "    \n",
    "    print(first_in_plateau[k])\n",
    "\n",
    "    first_in_tr1[k] = np.logical_and( cscClusterSize[k][:,0] >= 100, np.logical_and(cscClusterNStation[k][:,0]>=2, np.abs(cscClusterEta[k][:,0])<1.9))\n",
    "    first_in_tr2[k] = np.logical_and( cscClusterSize[k][:,0] >= 200, np.logical_and(cscClusterNStation[k][:,0]==1, np.abs(cscClusterEta[k][:,0])<1.9))\n",
    "    first_in_tr3[k] = np.logical_and( cscClusterSize[k][:,0] >= 500, np.abs(cscClusterEta[k][:,0])>=1.9)\n",
    " \n",
    "    ME1_veto_generic = ((cscClusterNRechitMinus11[k]==0) & (cscClusterNRechitPlus11[k]==0) & (cscClusterNRechitMinus12[k]==0) & (cscClusterNRechitPlus12[k]==0))#.all()\n",
    "    ME1_veto_cl_0 = (cscClusterNRechitMinus11[k][:,0]==0) & (cscClusterNRechitPlus11[k][:,0]==0) & (cscClusterNRechitMinus12[k][:,0]==0) & (cscClusterNRechitPlus12[k][:,0]==0)\n",
    "    ME1_veto_cl_1 = (cscClusterNRechitMinus11[k][:,1]==0) & (cscClusterNRechitPlus11[k][:,1]==0) & (cscClusterNRechitMinus12[k][:,1]==0) & (cscClusterNRechitPlus12[k][:,1]==0)\n",
    "   \n",
    "    print(ME1_veto_generic.sum())\n",
    "    sel_ev_pl[k] = np.logical_and(sel_ev[k],first_in_plateau[k])# & ME1_veto_generic#ME1_veto_cl_0 & ME1_veto_cl_1\n",
    "    \n",
    "    #since I have removed n clusters==2 above (to avoid blocking events with 3 clusters), need to apply it now \n",
    "    if category==1:\n",
    "        dPhi_csc_csc[k] = ( deltaPhi(  cscClusterPhi[k][:,0].flatten(), cscClusterPhi[k][:,1].flatten() ) )\n",
    "        dEta_csc_csc[k] = (  cscClusterEta[k][:,0].flatten() - cscClusterEta[k][:,1].flatten() )\n",
    "        dR_csc_csc[k]   = deltaR( cscClusterEta[k][:,0].flatten(), cscClusterPhi[k][:,0].flatten(), cscClusterEta[k][:,1].flatten(), cscClusterPhi[k][:,1].flatten() )\n",
    "        dt_csc_csc[k]   = (  cscClusterTime[k][:,0].flatten() - cscClusterTime[k][:,1].flatten() )\n",
    "        min_dphi = 1.8#.65\n",
    "        sel_ev_post[k] = np.logical_and(sel_ev_pl[k],np.abs(dPhi_csc_csc[k])>min_dphi)\n",
    "    if category==2:\n",
    "        dPhi_csc_dt[k] = ( deltaPhi(  cscClusterPhi[k][:,0].flatten(), dtClusterPhi[k][:,0].flatten() ) )\n",
    "        dEta_csc_dt[k] = (  cscClusterEta[k][:,0].flatten() - dtClusterEta[k][:,0].flatten() )\n",
    "        dR_csc_dt[k]   = deltaR( cscClusterEta[k][:,0].flatten(), cscClusterPhi[k][:,0].flatten(), dtClusterEta[k][:,0].flatten(), dtClusterPhi[k][:,0].flatten() )\n",
    "        sel_ev_post[k] = sel_ev_pl[k]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a06c0be",
   "metadata": {},
   "source": [
    "## Debug yields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "dd790a38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Selection \t Yield \t\t Eff. vs no cuts(%)\n",
      "No cuts\t\t 48.58\t 0.004\n",
      "2 CSC + 0 lep\t 0.00\t 0.000\n",
      "Trigger plateau\t 0.00\t 0.000\n",
      "Min dphi\t 0.00\t 0.000\n",
      "+ME1[0] veto\t 0.00\t 0.000\n",
      "+ME1[1] veto\t 0.00\t 0.000\n",
      "+ME1[01] veto\t 0.00\t 0.000\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "lumi=1\n",
    "weight['signal'] = weight['signal'].astype(bool)\n",
    "\n",
    "print(\"Selection \\t Yield \\t\\t Eff. vs no cuts(%)\")\n",
    "\n",
    "table_label = \"No cuts\"\n",
    "n_val = np.sum(all_weight['signal'])*lumi\n",
    "print(\"%s\\t\\t %.2f\\t %.3f\" % (table_label, n_val, (100*n_val)/signal_gen_yield))\n",
    "\n",
    "table_label = \"2 CSC + 0 lep\"\n",
    "n_val = (np.sum(weight['signal'])*lumi)\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (table_label, n_val, (100*n_val)/signal_gen_yield))\n",
    "\n",
    "table_label = \"Trigger plateau\"\n",
    "n_val = (np.sum(weight['signal'][ sel_ev_pl['signal']  ])*lumi)\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (table_label, n_val, (100*n_val)/signal_gen_yield))\n",
    "\n",
    "table_label = \"Min dphi\"\n",
    "n_val = (np.sum(weight['signal'][ sel_ev_post['signal']  ])*lumi)\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (table_label, n_val, (100*n_val)/signal_gen_yield))\n",
    "\n",
    "ME1_veto_cl_0 = (cscClusterNRechitMinus11[k][:,0]==0) & (cscClusterNRechitPlus11[k][:,0]==0) & (cscClusterNRechitMinus12[k][:,0]==0) & (cscClusterNRechitPlus12[k][:,0]==0)\n",
    "ME1_veto_cl_1 = (cscClusterNRechitMinus11[k][:,1]==0) & (cscClusterNRechitPlus11[k][:,1]==0) & (cscClusterNRechitMinus12[k][:,1]==0) & (cscClusterNRechitPlus12[k][:,1]==0)\n",
    "\n",
    "\n",
    "table_label = \"+ME1[0] veto\"\n",
    "n_val = (np.sum(weight['signal'][ sel_ev_post['signal'] & ME1_veto_cl_0  ])*lumi)\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (table_label, n_val, (100*n_val)/signal_gen_yield))\n",
    "\n",
    "\n",
    "table_label = \"+ME1[1] veto\"\n",
    "n_val = (np.sum(weight['signal'][ sel_ev_post['signal'] & ME1_veto_cl_1  ])*lumi)\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (table_label, n_val, (100*n_val)/signal_gen_yield))\n",
    "\n",
    "table_label = \"+ME1[01] veto\"\n",
    "n_val = (np.sum(weight['signal'][ sel_ev_post['signal'] & ME1_veto_cl_0  & ME1_veto_cl_1  ])*lumi)\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (table_label, n_val, (100*n_val)/signal_gen_yield))\n",
    "\n",
    "print(evtNum['signal'][ sel_ev_post['signal'] & ME1_veto_cl_0  & ME1_veto_cl_1  ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c944d53",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e398b24d",
   "metadata": {},
   "source": [
    "* Look only at 10 k entries and find the issue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7c90aac",
   "metadata": {},
   "source": [
    "Old method"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b111d72e",
   "metadata": {},
   "source": [
    "Selection \t Yield \t\t Eff. vs no cuts(%)\n",
    "No cuts\t\t 48.58\t 0.004\n",
    "2 CSC + 0 lep\t 22.00\t 0.002\n",
    "Trigger plateau\t 14.00\t 0.001\n",
    "Min dphi\t 12.00\t 0.001\n",
    "+ME1[0] veto\t 12.00\t 0.001\n",
    "+ME1[1] veto\t 12.00\t 0.001\n",
    "+ME1[01] veto\t 12.00\t 0.001\n",
    "[4137 4215 5325 5362 5746 7166 7778 8496 8654 9081 9430 9574]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b9ddcf1",
   "metadata": {},
   "source": [
    "New"
   ]
  },
  {
   "cell_type": "raw",
   "id": "81e69d7e",
   "metadata": {},
   "source": [
    "Selection \t Yield \t\t Eff. vs no cuts(%)\n",
    "No cuts\t\t 48.58\t 0.004\n",
    "2 CSC + 0 lep\t 151.00\t 0.014\n",
    "Trigger plateau\t 77.00\t 0.007\n",
    "Min dphi\t 62.00\t 0.006\n",
    "+ME1[0] veto\t 34.00\t 0.003\n",
    "+ME1[1] veto\t 23.00\t 0.002\n",
    "+ME1[01] veto\t 11.00\t 0.001\n",
    "[4137 4215 5325 5362 5746 7166 7778 8496 8654 9081 9574]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5078501e",
   "metadata": {},
   "source": [
    "Event 9430 is discarded by the latter method --> inspect!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "96927dd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "cut = sel_ev_post['signal'] & ME1_veto_cl_0  & ME1_veto_cl_1 & (evtNum['signal'] == 9430)\n",
    "\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43b9bc4f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6acd5156",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3f7ffd55",
   "metadata": {},
   "source": [
    "Previous ME1 veto in the first loop"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4a548ab1",
   "metadata": {},
   "source": [
    "Selection \t     Yield \t\t Eff. vs no cuts(%)\n",
    "No cuts\t\t     1118331.31\t 100.000\n",
    "2 CSC + 0 lep\t 2545.64\t 0.228\n",
    "Trigger plateau\t 1324.77\t 0.118\n",
    "Min dphi\t     1156.23\t 0.103\n",
    "+ME1[0] veto\t 1156.23\t 0.103\n",
    "+ME1[1] veto\t 1156.23\t 0.103\n",
    "+ME1[01] veto\t 1156.23\t 0.103"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b3be526",
   "metadata": {},
   "source": [
    "ME1 veto done after the trigger plateau with surviving clusters is different for some reasons"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7e570152",
   "metadata": {},
   "source": [
    "Selection \t     Yield \t\t Eff. vs no cuts(%)\n",
    "No cuts\t\t     1118331.31\t 100.000\n",
    "2 CSC + 0 lep\t 15543.23\t 1.390\n",
    "Trigger plateau\t 6951.87\t 0.622\n",
    "Min dphi\t     5470.03\t 0.489\n",
    "+ME1[0] veto\t 2795.12\t 0.250\n",
    "+ME1[1] veto\t 2232.72\t 0.200\n",
    "+ME1[01] veto\t 1122.40\t 0.100\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61d84f6e",
   "metadata": {},
   "source": [
    "If I include ME veto in the trigger plateau, I get the same as the independent computation"
   ]
  },
  {
   "cell_type": "raw",
   "id": "cbd0bff6",
   "metadata": {},
   "source": [
    "Selection \t     Yield \t\t Eff. vs no cuts(%)\n",
    "No cuts\t\t     1118331.31\t 100.000\n",
    "2 CSC + 0 lep\t 15543.23\t 1.390\n",
    "Trigger plateau\t 1277.05\t 0.114\n",
    "Min dphi\t     1122.40\t 0.100\n",
    "+ME1[0] veto\t 1122.40\t 0.100\n",
    "+ME1[1] veto\t 1122.40\t 0.100\n",
    "+ME1[01] veto\t 1122.40\t 0.100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e392ee4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86aed1de",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a393ee7e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e818b6af",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94b25ad2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1a32c9c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f709ee36",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87ec707f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28e8c316",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "id": "c934be1c",
   "metadata": {},
   "source": [
    "#print(ME11_veto['signal'])\n",
    "print(\"Selection \\t Yield \\t\\t Eff.(%)\")\n",
    "print(\"%s\\t\\t %.2f\\t %.2f\" % (\"No cuts\",(100*np.sum(all_weight['signal'])*lumi) , (100*np.sum(all_weight['signal'])*lumi)/signal_gen_yield))\n",
    "#print( (100*np.sum(all_weight['signal'][ ME11_veto['signal']==True  ])*lumi)/signal_gen_yield)\n",
    "\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"Presel 2 CSC\", (100*np.sum(weight['signal'])*lumi) , (100*np.sum(weight['signal'])*lumi)/signal_gen_yield))\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"Trigger plateau\", (100*np.sum(weight['signal'][ sel_ev['signal'] ])*lumi), (100*np.sum(weight['signal'][ sel_ev['signal'] ])*lumi)/signal_gen_yield))\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"Min dphi>1.8\", (100*np.sum(weight['signal'][ sel_ev_post['signal'] ])*lumi), (100*np.sum(weight['signal'][ sel_ev_post['signal'] ])*lumi)/signal_gen_yield))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8492d52b",
   "metadata": {},
   "source": [
    "* Tune min dphi cut to remove poorly-clustered events"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18a616e6",
   "metadata": {},
   "source": [
    "# ABCD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c5a7ccde",
   "metadata": {},
   "outputs": [],
   "source": [
    "PHI_MIN = 1.8\n",
    "PHI_CUT = 2.85\n",
    "#PHI_MAX = 2.\n",
    "PHI_MAX = 3.15\n",
    "\n",
    "N_MIN = 50\n",
    "N_MAX = 500000000\n",
    "N_CUT = 240\n",
    "\n",
    "n_ev = 5000\n",
    "\n",
    "pre_cut = {}\n",
    "a = {}\n",
    "b = {}\n",
    "c = {}\n",
    "d = {}"
   ]
  },
  {
   "cell_type": "raw",
   "id": "0e76cf50",
   "metadata": {},
   "source": [
    "def extrapolation(vec,min_x,max_x,nbins):\n",
    "    h = create_TH1D(vec,title=\"\",axis_title = ['',''], binning=[nbins,min_x,max_x])#,bin_list=True)\n",
    "    #can = rt.TCanvas('c','c', 1000, 800)\n",
    "    #can.cd()\n",
    "    h.SetLineWidth(2)\n",
    "    h.SetMarkerStyle(20)\n",
    "    #f1 = rt.TF1(\"f1\",\"[0]*x + [1]\",min_x,max_x)\n",
    "    #f1.SetParameter(0,-2)\n",
    "    #f1.SetParameter(1,60)\n",
    "    f1 = rt.TF1(\"f1\",\"exp([0]*x + [1]/(x))/[2]\",min_x,max_x)\n",
    "    f1.SetParameter(0,-0.03)\n",
    "    f1.SetParameter(1,200)\n",
    "    f1.SetParameter(2,74)\n",
    "    h.Fit('f1', 'LRSQ+', '', min_x,max_x)\n",
    "    my_f1 = h.GetFunction(\"f1\")            \n",
    "    h.GetListOfFunctions().Remove(h.GetFunction(\"f1\"))\n",
    "    #h.Draw(\"PE\")\n",
    "    #h.GetXaxis().SetTitle(\"#Delta#varphi\")\n",
    "    ##h.SetTitle(\"High n hits region: \"+str(N_CUT)+\" < N hits < \" + str(N_MAX))\n",
    "    #my_f1.SetLineStyle(2)\n",
    "    #my_f1.Draw('L same')\n",
    "    #can.SetLogy()\n",
    "    #can.Draw()\n",
    "    #can.Print(\"test.pdf\")\n",
    "    #can.Print(\"test.png\")\n",
    "    return h, my_f1"
   ]
  },
  {
   "cell_type": "raw",
   "id": "da241094",
   "metadata": {},
   "source": [
    "def run_abcd(var1,var2,MIN1,CUT1,MAX1,MIN2,CUT2,MAX2,sel,phi_th=2.2,n_th=150,do_b_extrapolation=False):\n",
    "            \n",
    "    d = np.sum([sel & (var1>=CUT1) & (var1<MAX1) & ( var2>MIN2 ) & ( var2<CUT2 )])#high n, low phi\n",
    "    #risk of unblinding D when phi gets too large. Try to predict it instead.\n",
    "    if CUT2>=phi_th:\n",
    "        low_d = np.sum([sel & (var1>=CUT1) & (var1<MAX1) & ( var2>MIN2 ) & ( var2<phi_th )])\n",
    "        high_d_flat = (CUT2-phi_th)*low_d/( phi_th - MIN2 )\n",
    "        d_pred = low_d+high_d_flat\n",
    "        print(\"d\\t d low\\t d pred\")\n",
    "        print(\"%i\\t %i\\t  %.2f\" % (d,low_d, d_pred))\n",
    "        d = low_d+high_d_flat  \n",
    "    b = np.sum([sel & (var1>=MIN1) & (var1<CUT1) & ( var2>=CUT2 ) & ( var2<MAX2 )])#low n, high phi\n",
    "    #risk of unblinding B when n hits gets too large. Try to predict it instead.\n",
    "    if CUT1>=n_th and do_b_extrapolation:\n",
    "        low_b_vec = var1[sel & (var1>=MIN1) & (var1<n_th) & ( var2>=CUT2 ) & ( var2<MAX2 )]\n",
    "        low_b = np.sum([sel & (var1>=MIN1) & (var1<n_th) & ( var2>=CUT2 ) & ( var2<MAX2 )])\n",
    "        h,my_f1 = extrapolation(low_b_vec,MIN1,CUT1,20)\n",
    "        b_pred = my_f1.Integral(n_th,CUT1) + low_b\n",
    "        print(\"b\\t b pred\")\n",
    "        print(\"%i\\t%.2f\" % (b,b_pred))\n",
    "        #high_b_flat = (CUT2-phi_th)*low_b/( phi_th - MIN2 )\n",
    "        b = my_f1.Integral(n_th,CUT1) + low_b\n",
    "        \n",
    "    c = np.sum([sel & (var1>=MIN1) & (var1<CUT1) & ( var2>MIN2 ) & ( var2<CUT2 )])#low n, low phi \n",
    "    #if CUT2>=phi_th:\n",
    "    #    #low_c = np.sum([sel & (var1>=MIN1) & (var1<CUT1) & ( var2>MIN2 ) & ( var2<phi_th )])\n",
    "    #    #high_c_flat = (CUT2-phi_th)*low_c/( phi_th - MIN2 )\n",
    "    #    #c_pred = low_c+high_c_flat\n",
    "    #    #print(\"c\\t c low\\t c pred\")\n",
    "    #    #print(\"%i\\t %i\\t  %.2f\" % (c,low_c, c_pred))\n",
    "    #    #c = low_c+high_c_flat  \n",
    "    a = np.sum([sel & (var1>=CUT1) & (var1<MAX1) & ( var2>=CUT2 ) & ( var2<MAX2 )])#SR  \n",
    "    pred = (b*d)/c\n",
    "    unc_pred = (1./b if b>0 else 0 + 1./d if d>0 else 0 + 1./c if c>0 else 0)**0.5*pred\n",
    "    if math.isnan(unc_pred): \n",
    "        p_value = float(\"nan\")\n",
    "        z_value = float(\"nan\")\n",
    "    else:\n",
    "        mu = np.random.normal(pred, unc_pred, n_ev)\n",
    "        p_value = 0.0\n",
    "        for i in mu:\n",
    "            if i < 0:continue\n",
    "            n = np.random.poisson(i, n_ev)\n",
    "            p_value += np.count_nonzero(n>=a)\n",
    "        p_value = p_value/n_ev**2\n",
    "        z_value = (norm.ppf(p_value))\n",
    "    return a,b,c,d,pred,unc_pred, p_value, z_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9b60309a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def return_abcd_masks(var1,var2,MIN1,CUT1,MAX1,MIN2,CUT2,MAX2,sel):\n",
    "    d = sel & (var1>=CUT1) & (var1<MAX1) & ( var2>MIN2 ) & ( var2<CUT2 )\n",
    "    b = sel & (var1>=MIN1) & (var1<CUT1) & ( var2>=CUT2 ) & ( var2<MAX2 ) \n",
    "    c = sel & (var1>=MIN1) & (var1<CUT1) & ( var2>MIN2 ) & ( var2<CUT2 )  \n",
    "    a = sel & (var1>=CUT1) & (var1<MAX1) & ( var2>=CUT2 ) & ( var2<MAX2 )\n",
    "    return a,b,c,d"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b7f519e6",
   "metadata": {},
   "source": [
    "def make_datacard_2tag(outDataCardsDir,modelName,  signal_rate, normalization, bkg_rate, observation, bkg_unc, bkg_unc_name, sig_unc, sig_unc_name,signal_region, prefix):\n",
    "    a,b,c,d = bkg_rate[0], bkg_rate[1], bkg_rate[2], bkg_rate[3]\n",
    "    nSig = len(signal_rate.keys())\n",
    "    text_file = open(outDataCardsDir+modelName+\".txt\", \"w\")\n",
    "    text_file.write('# signal norm {0} \\n'.format(normalization))\n",
    "\n",
    "    text_file.write('imax {0} \\n'.format(4))\n",
    "    text_file.write('jmax {0} \\n'.format(nSig))\n",
    "    text_file.write('kmax * \\n')\n",
    "    text_file.write('shapes * * FAKE \\n')\n",
    "\n",
    "\n",
    "    text_file.write('--------------- \\n')\n",
    "    text_file.write('--------------- \\n')\n",
    "    text_file.write('bin \\t chA \\t chB \\t chC \\t chD \\n')\n",
    "    text_file.write('observation \\t {0:6.2f} \\t {1:6.2f} \\t {2:6.2f} \\t {3:6.2f} \\n'.format(observation[0],observation[1],observation[2],observation[3]))\n",
    "    text_file.write('------------------------------ \\n')\n",
    "    text_file.write('bin '+'\\t chA ' * (1+nSig) + '\\t chB ' * (1+nSig) +'\\t chC '*(1+nSig) +'\\t chD '*(1+nSig) +'\\n')\n",
    "    process_name = '\\t '+ (' \\t ').join(list(signal_rate.keys())) + '\\t bkg '\n",
    "    text_file.write('process ' + process_name * 4 + '\\n')\n",
    "    process_number = '\\t '+ (' \\t ').join(list((np.arange(nSig)*-1).astype(str))) + '\\t 1'\n",
    "    text_file.write('process ' + process_number * 4 + '\\n')\n",
    "    rate_string = 'rate'\n",
    "    for i in range(4):# 4 bins\n",
    "        for k,v in signal_rate.items():\n",
    "            rate_string +='\\t {0:e} '.format(v[i])\n",
    "        rate_string += '\\t 1 '\n",
    "    text_file.write(rate_string+'\\n')\n",
    "    text_file.write('------------------------------ \\n')\n",
    "\n",
    "    text_file.write(prefix+'A   rateParam       chA     bkg      (@0*@2/@1)                    '+prefix+'B,'+prefix+'C,'+prefix+'D \\n')\n",
    "    if b == 0: text_file.write(prefix+'B   rateParam       chB     bkg     {0:.2f}        [0,{1:.2f}] \\n'.format(b, c*7))\n",
    "    else: text_file.write(prefix+'B   rateParam       chB     bkg     {0:.2f}        [0,{1:.2f}] \\n'.format(b, b*7))\n",
    "    text_file.write(prefix+'C   rateParam       chC     bkg     {0:.2f}        [0,{1:.2f}] \\n'.format(c, c*7))\n",
    "    if d == 0:text_file.write(prefix+'D   rateParam       chD     bkg     {0:.2f}        [0,{1:.2f}] \\n'.format(d, c*7))\n",
    "    else: text_file.write(prefix+'D   rateParam       chD     bkg     {0:.2f}        [0,{1:.2f}] \\n'.format(d, d*7))\n",
    "\n",
    "\n",
    "    #for k,v in signal_rate.items():\n",
    "    #    text_file.write('norm rateParam * {0} 1  \\n'.format(k))\n",
    "        \n",
    "    '''\n",
    "\n",
    "\n",
    "  #### uncertainties ####\n",
    "    for k,v in sig_unc.items():assert(len(sig_unc_name)==len(v))\n",
    "    for i in range(len(sig_unc_name)):\n",
    "        if 'mc_stats' in sig_unc_name[i]:\n",
    "            for j, bin in enumerate(['A', 'B', 'C', 'D']):#bin\n",
    "                    for l, k in enumerate(sig_unc.keys()): #channels\n",
    "                        before = (len(sig_unc.keys())+1)*j+l\n",
    "                        after = (len(sig_unc.keys())+1)*4-before-1\n",
    "                        if sig_unc[k][i][j] > 0.0: text_file.write(sig_unc_name[i]+'_'+k+'_'+bin+' \\t gmN ' +str(int(sig_unc[k][i][j]))+ '  '+'\\t -  '*before + str(signal_rate[k][j]/int(sig_unc[k][i][j])) + '\\t - '*after +'\\n')\n",
    "\n",
    "        else:\n",
    "\n",
    "            unc_text = sig_unc_name[i]+' \\t lnN'\n",
    "            if len(sig_unc[list(sig_unc.keys())[0]][i])==4:#symmetric uncertainties\n",
    "                for j in range(4):#bin\n",
    "                    for k,v in sig_unc.items():\n",
    "                        if v[i][j] == 0.0:unc_text += ' \\t -'\n",
    "                        else: unc_text += ' \\t '+str(v[i][j]+1)\n",
    "                    unc_text += '\\t - '\n",
    "            else:#asymmetric\n",
    "                for j in range(4):#bin A, B, C, D\n",
    "                    for k,v in sig_unc.items():\n",
    "                        if  v[i][j] == 0.0 and v[i][j+4] == 0.0: unc_text += ' \\t -'\n",
    "                        else:unc_text += ' \\t {0}/{1}'.format(1-v[i][j],1+v[i][j+4])\n",
    "                    unc_text += '\\t -'\n",
    "            text_file.write(unc_text + ' \\n')\n",
    "    '''\n",
    "    for i in range(len(bkg_unc_name)):\n",
    "        bkg_unc_text = bkg_unc_name[i] + ' \\t lnN ' + '\\t - '*(4*nSig+3) + '\\t ' + str(1+bkg_unc[i]) + ' \\n'\n",
    "        text_file.write(bkg_unc_text)\n",
    "    \n",
    "\n",
    "    text_file.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80aeefea",
   "metadata": {},
   "source": [
    "## SR background levels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bea109b",
   "metadata": {},
   "source": [
    "* Scan minimum delta phi threshold\n",
    "* Scan optimal delta phi cut"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9213a01e",
   "metadata": {},
   "source": [
    "## Scan phi cut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dc4a9b02",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k\t phi\t B\t C\t D\t A\t pred\t\t zval\t perc\n",
      "signal 2.85 376.3480396568775 590.7999835163355 520.1214528828859 289.96317841112614 1777.232654467225\n"
     ]
    }
   ],
   "source": [
    "#0, 0.5, 0.8, 1, 1.1, 1.2\n",
    "dr_min_list = [0.5]#[0.9]#[1.9]\n",
    "phi_min_scan_list = [0,0.2,0.4,0.5,0.6,0.7,0.75,0.8,0.85,0.9,0.95,1.,1.05,1.1,1.15,1.2,1.25,1.3,1.35,1.4,1.45,1.5]\n",
    "\n",
    "phi_min_scan_list = [1.1]\n",
    "phi_cut_scan_list = np.linspace(1.5,3.1,64+1)#[1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,3.]\n",
    "phi_cut_scan_list = np.linspace(1.5,3.1,32+1)\n",
    "\n",
    "phi_cut_scan_list = [2.85]\n",
    "\n",
    "time_cut_string = \"IT\"\n",
    "plot_string_plus = plot_string+\"_\"+time_cut_string+\"_dR\"+str( dr_min_list[0] ).replace(\".\",\"p\")\n",
    "leghead = \"|#Delta#eta|<1.9, \"+time_cut_string+\", #DeltaR>\"+str(dr_min_list[0])+\",\"\n",
    "\n",
    "#one per min dphi cut\n",
    "bkg_rate = {}\n",
    "#bkg_unc = {}\n",
    "observation = {}\n",
    "signal_rate = {}\n",
    "\n",
    "normalization = 1\n",
    "\n",
    "for phi_cut_scan in phi_cut_scan_list:\n",
    "    signal_rate[phi_cut_scan] = {}\n",
    "\n",
    "print('k\\t phi\\t B\\t C\\t D\\t A\\t pred\\t\\t zval\\t perc')\n",
    "if category==1:\n",
    "    #Signal\n",
    "    #Data\n",
    "    for k in ['signal']:\n",
    "        tmp_dphi = np.abs( dPhi_csc_csc[k] )\n",
    "        tmp_deta = np.abs(  dEta_csc_csc[k] )\n",
    "        tmp_dR = np.abs(  dR_csc_csc[k] )\n",
    "        tmp_dt = np.abs( dt_csc_csc[k] )\n",
    "        var = cscClusterSize[k][:,1]\n",
    "        t_var = cscClusterTime[k][:,1]\n",
    "        t_spread_var = cscClusterTimeSpread[k][:,1]\n",
    "        dphi = np.abs(dPhi_csc_csc[k])\n",
    "        tmp_cut = sel_ev_post[k]  & (np.abs(tmp_deta)<1.9) & (np.abs(t_var)<15) & (np.abs(t_spread_var)<20) & (sel_jetveto_csc[k][:,1]==True) & (sel_muonveto_csc[k][:,1]==True)#& (t_var<-15) #& (t_var<-15) # #& (tmp_dR>r) #(tmp_deta>=1.9)  & ( t_var<-15  )\n",
    "        for phi_cut_scan in phi_cut_scan_list:\n",
    "            if phi_cut_scan<=PHI_MIN: continue\n",
    "            if 'data' in k:\n",
    "                a[k],b[k],c[k],d[k],pred,unc_pred, p_value, z_value = run_abcd(\n",
    "                    var,dphi,\n",
    "                    N_MIN,N_CUT,N_MAX,\n",
    "                    PHI_MIN,phi_cut_scan,PHI_MAX,\n",
    "                    tmp_cut\n",
    "                    )\n",
    "                print(\"%s\\t %.2f\\t %i\\t %i\\t %i\\t %i\\t %.2f +- %.2f\\t %.2f\\t %.2f\" % (k,phi_cut_scan,b[k],c[k],d[k],a[k],pred,unc_pred,z_value, 100*(pred-a[k])/a[k]))\n",
    "                observation[phi_cut_scan] = [ pred,b[k],c[k],d[k] ]\n",
    "                bkg_rate[phi_cut_scan] = [ pred,b[k],c[k],d[k] ]\n",
    "                #bkg_unc[phi_min_scan] = [ unc_pred,np.sqrt(b[k]),np.sqrt(c[k]),np.sqrt(d[k]) ]\n",
    "            else:\n",
    "                ma,mb,mc,md = return_abcd_masks(\n",
    "                    var,dphi,\n",
    "                    N_MIN,N_CUT,N_MAX,\n",
    "                    PHI_MIN,phi_cut_scan,PHI_MAX,\n",
    "                    tmp_cut\n",
    "                    )\n",
    "                a[k] = (np.sum(weight[k][ ma ]))*lumi\n",
    "                b[k] = (np.sum(weight[k][ mb ]))*lumi\n",
    "                c[k] = (np.sum(weight[k][ mc ]))*lumi\n",
    "                d[k] = (np.sum(weight[k][ md ]))*lumi\n",
    "                print(k,phi_cut_scan,a[k],b[k],c[k],d[k], a[k]+b[k]+c[k]+d[k])\n",
    "                signal_rate[phi_cut_scan]['ggH'] = [ a[k]*normalization,b[k]*normalization,c[k]*normalization,d[k]*normalization ]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "67748929",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Selection \t Yield \t\t Eff. vs no cuts(%)\n",
      "No cuts\t\t 1118331.31\t 100.000\n",
      "2 CSC + 0 lep\t 15543.23\t 1.390\n",
      "ME1 veto [0]\t 7406.14\t 0.662\n",
      "ME1 veto [1]\t 4898.57\t 0.438\n",
      "ME1 veto [0&1]\t 2461.67\t 0.220\n",
      "Trigger plateau\t 15543.23\t 1.390\n",
      "Trigger + ME1[0]\t 7406.14\t 0.662\n",
      "Trigger + ME1[0&1]\t 2461.67\t 0.220\n",
      "Min dphi>1.8\t 2795.12\t 0.250\n",
      "Delta eta<1.9\t 2474.96\t 0.221\n",
      "In time cut\t 2341.45\t 0.209\n",
      "Time spread\t 2230.30\t 0.199\n",
      "Jet muon vetoes\t 1777.23\t 0.159\n",
      "ABCD plane\t 1777.23\t 0.16\n",
      "C\t 520.12\t 0.05\n",
      "B\t 590.80\t 0.05\n",
      "D\t 289.96\t 0.03\n",
      "A\t 376.35\t 0.03\n"
     ]
    }
   ],
   "source": [
    "#print(ME11_veto['signal'])\n",
    "print(\"Selection \\t Yield \\t\\t Eff. vs no cuts(%)\")\n",
    "print(\"%s\\t\\t %.2f\\t %.3f\" % (\"No cuts\",(np.sum(all_weight['signal'])*lumi) , (100*np.sum(all_weight['signal'])*lumi)/signal_gen_yield))\n",
    "#print( (100*np.sum(all_weight['signal'][ ME11_veto['signal'][:,0]==True  ])*lumi)/signal_gen_yield)\n",
    "\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"2 CSC + 0 lep\", (np.sum(weight['signal'])*lumi) , (100*np.sum(weight['signal'])*lumi)/signal_gen_yield))\n",
    "\n",
    "ME1_veto_cl_0 = (cscClusterNRechitMinus11[k][:,0]==0) & (cscClusterNRechitPlus11[k][:,0]==0) & (cscClusterNRechitMinus12[k][:,0]==0) & (cscClusterNRechitPlus12[k][:,0]==0)\n",
    "\n",
    "ME1_veto_cl_1 = (cscClusterNRechitMinus11[k][:,1]==0) & (cscClusterNRechitPlus11[k][:,1]==0) & (cscClusterNRechitMinus12[k][:,1]==0) & (cscClusterNRechitPlus12[k][:,1]==0)\n",
    "\n",
    "\n",
    "tmp_cut = ME1_veto_cl_0\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"ME1 veto [0]\", (np.sum(weight['signal'][ tmp_cut ])*lumi), (100*np.sum(weight['signal'][ tmp_cut ])*lumi)/signal_gen_yield))\n",
    "\n",
    "tmp_cut = ME1_veto_cl_1\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"ME1 veto [1]\", (np.sum(weight['signal'][ tmp_cut ])*lumi), (100*np.sum(weight['signal'][ tmp_cut ])*lumi)/signal_gen_yield))\n",
    "\n",
    "tmp_cut = ME1_veto_cl_1 & ME1_veto_cl_0\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"ME1 veto [0&1]\", (np.sum(weight['signal'][ tmp_cut ])*lumi), (100*np.sum(weight['signal'][ tmp_cut ])*lumi)/signal_gen_yield))\n",
    "\n",
    "\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"Trigger plateau\", (np.sum(weight['signal'][ sel_ev['signal']  ])*lumi), (100*np.sum(weight['signal'][ sel_ev['signal'] ])*lumi)/signal_gen_yield))\n",
    "\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"Trigger + ME1[0]\", (np.sum(weight['signal'][ sel_ev['signal'] & ME1_veto_cl_0 ])*lumi), (100*np.sum(weight['signal'][ sel_ev['signal'] & ME1_veto_cl_0 ])*lumi)/signal_gen_yield))\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"Trigger + ME1[0&1]\", (np.sum(weight['signal'][ sel_ev['signal'] & ME1_veto_cl_0 & ME1_veto_cl_1 ])*lumi), (100*np.sum(weight['signal'][ sel_ev['signal'] & ME1_veto_cl_0 & ME1_veto_cl_1 ])*lumi)/signal_gen_yield))\n",
    "\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"Min dphi>1.8\", (np.sum(weight['signal'][ sel_ev_post['signal'] ])*lumi), (100*np.sum(weight['signal'][ sel_ev_post['signal'] ])*lumi)/signal_gen_yield))\n",
    "\n",
    "\n",
    "\n",
    "k='signal'\n",
    "tmp_dphi = np.abs( dPhi_csc_csc[k] )\n",
    "tmp_deta = np.abs(  dEta_csc_csc[k] )\n",
    "var = cscClusterSize[k][:,1]\n",
    "t_var = cscClusterTime[k][:,1]\n",
    "t_spread_var = cscClusterTimeSpread[k][:,1]\n",
    "\n",
    "\n",
    "\n",
    "dphi = np.abs(dPhi_csc_csc[k])\n",
    "tmp_cut = sel_ev_post[k]  & (np.abs(tmp_deta)<1.9)# & (np.abs(t_var)<15) & (np.abs(t_spread_var)<20) & (sel_jetveto_csc[k][:,1]==True) & (sel_muonveto_csc[k][:,1]==True)#& (t_var<-15) #& (t_var<-15) # #& (tmp_dR>r) #(tmp_deta>=1.9)  & ( t_var<-15  )\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"Delta eta<1.9\", (np.sum(weight['signal'][ tmp_cut ])*lumi), (100*np.sum(weight['signal'][ tmp_cut ])*lumi)/signal_gen_yield))\n",
    "tmp_cut = sel_ev_post[k]  & (np.abs(tmp_deta)<1.9) & (np.abs(t_var)<15)# & (np.abs(t_spread_var)<20) & (sel_jetveto_csc[k][:,1]==True) & (sel_muonveto_csc[k][:,1]==True)#& (t_var<-15) #& (t_var<-15) # #& (tmp_dR>r) #(tmp_deta>=1.9)  & ( t_var<-15  )\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"In time cut\", (np.sum(weight['signal'][ tmp_cut ])*lumi), (100*np.sum(weight['signal'][ tmp_cut ])*lumi)/signal_gen_yield))\n",
    "tmp_cut = sel_ev_post[k]  & (np.abs(tmp_deta)<1.9) & (np.abs(t_var)<15) & (np.abs(t_spread_var)<20)# & (sel_jetveto_csc[k][:,1]==True) & (sel_muonveto_csc[k][:,1]==True)#& (t_var<-15) #& (t_var<-15) # #& (tmp_dR>r) #(tmp_deta>=1.9)  & ( t_var<-15  )\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"Time spread\", (np.sum(weight['signal'][ tmp_cut ])*lumi), (100*np.sum(weight['signal'][ tmp_cut ])*lumi)/signal_gen_yield))\n",
    "tmp_cut = sel_ev_post[k]  & (np.abs(tmp_deta)<1.9) & (np.abs(t_var)<15) & (np.abs(t_spread_var)<20) & (sel_jetveto_csc[k][:,1]==True) & (sel_muonveto_csc[k][:,1]==True)#& (t_var<-15) #& (t_var<-15) # #& (tmp_dR>r) #(tmp_deta>=1.9)  & ( t_var<-15  )\n",
    "print(\"%s\\t %.2f\\t %.3f\" % (\"Jet muon vetoes\", (np.sum(weight['signal'][ tmp_cut ])*lumi), (100*np.sum(weight['signal'][ tmp_cut ])*lumi)/signal_gen_yield))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ma,mb,mc,md = return_abcd_masks(\n",
    "                        var,dphi,\n",
    "                        N_MIN,N_CUT,N_MAX,\n",
    "                        PHI_MIN,2.85,PHI_MAX,\n",
    "                        tmp_cut\n",
    "                        )\n",
    "#all_abcd = np.array(ma or mb or mc or md).flatten()\n",
    "#ma = np.array(ma).flatten()\n",
    "#mb = np.array(mb).flatten()\n",
    "#mc = np.array(mc).flatten()\n",
    "#md = np.array(md).flatten()\n",
    "\n",
    "MIN1 = N_MIN\n",
    "CUT1 = N_CUT\n",
    "MAX1 = N_MAX\n",
    "MIN2 = PHI_MIN\n",
    "CUT2 = 2.85\n",
    "MAX2 = PHI_MAX\n",
    "#print(MIN1,CUT1,MAX1)\n",
    "#print(MIN2,CUT2,MAX2)\n",
    "manual_md = tmp_cut & (var>=CUT1) & (var<MAX1) & ( dphi>MIN2 ) & ( dphi<CUT2 )#.astype(bool)\n",
    "manual_mb = tmp_cut & (var>=MIN1) & (var<CUT1) & ( dphi>=CUT2 ) & ( dphi<MAX2 )#.astype(bool)\n",
    "manual_mc = tmp_cut & (var>=MIN1) & (var<CUT1) & ( dphi>MIN2 ) & ( dphi<CUT2 )#.astype(bool)\n",
    "manual_ma = tmp_cut & (var>=CUT1) & (var<MAX1) & ( dphi>=CUT2 ) & ( dphi<MAX2 )#.astype(bool)\n",
    "\n",
    "all_abcd = tmp_cut & (manual_ma | manual_mb | manual_mc | manual_md)\n",
    "tmp_cut2 = tmp_cut & all_abcd\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"ABCD plane\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "tmp_cut2 = tmp_cut & manual_mc\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"C\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "tmp_cut2 = tmp_cut & manual_mb\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"B\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "tmp_cut2 = tmp_cut & manual_md\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"D\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "tmp_cut2 = tmp_cut & manual_ma\n",
    "print(\"%s\\t %.2f\\t %.2f\" % (\"A\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "\n",
    "#print(\"-----\")\n",
    "#all_abcd = tmp_cut & (ma | mb | mc | md)\n",
    "#tmp_cut2 = tmp_cut & all_abcd\n",
    "#print(\"%s\\t %.2f\\t %.2f\" % (\"ABCD plane\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "#tmp_cut2 = tmp_cut & mc\n",
    "#print(\"%s\\t %.2f\\t %.2f\" % (\"C\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "#tmp_cut2 = tmp_cut & mb\n",
    "#print(\"%s\\t %.2f\\t %.2f\" % (\"B\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "#tmp_cut2 = tmp_cut & md\n",
    "#print(\"%s\\t %.2f\\t %.2f\" % (\"D\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n",
    "#tmp_cut2 = tmp_cut & ma\n",
    "#print(\"%s\\t %.2f\\t %.2f\" % (\"A\", (np.sum(weight['signal'][ tmp_cut2 ])*lumi), (100*np.sum(weight['signal'][ tmp_cut2 ])*lumi)/signal_gen_yield))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5677aac8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "run3",
   "language": "python",
   "name": "run3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}