Commits (59)
......@@ -24,7 +24,8 @@ atlas_depends_on_subdirs( PUBLIC
LArCalorimeter/LArIdentifier
LArCalorimeter/LArRawEvent
TileCalorimeter/TileCalib/TileCalibBlobObjs
TileCalorimeter/TileIdentifier )
TileCalorimeter/TileIdentifier
Event/xAOD/xAODCaloEvent )
# External dependencies:
find_package( CLHEP )
......@@ -34,7 +35,7 @@ atlas_add_component( CaloJiveXML
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS}
LINK_LIBRARIES ${CLHEP_LIBRARIES} CaloIdentifier AthenaBaseComps GaudiKernel LArToolsLib TileConditionsLib TileEvent JiveXMLLib CaloDetDescrLib CaloEvent AthenaKernel Identifier LArIdentifier LArRawEvent TileCalibBlobObjs TileIdentifier )
LINK_LIBRARIES ${CLHEP_LIBRARIES} CaloIdentifier AthenaBaseComps GaudiKernel LArToolsLib TileConditionsLib TileEvent JiveXMLLib CaloDetDescrLib CaloEvent AthenaKernel Identifier LArIdentifier LArRawEvent TileCalibBlobObjs TileIdentifier xAODCaloEvent)
# Install files from the package:
atlas_install_headers( CaloJiveXML )
......
......@@ -13,6 +13,8 @@
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
#include "xAODCaloEvent/CaloClusterContainer.h"
class CaloClusterContainer;
namespace JiveXML{
......@@ -42,7 +44,7 @@ namespace JiveXML{
/// Retrieve all the data
virtual StatusCode retrieve(ToolHandle<IFormatTool> &FormatTool);
const DataMap getData(const CaloClusterContainer*);
const DataMap getData(const xAOD::CaloClusterContainer*);
/// Return the name of the data type
virtual std::string dataTypeName() const { return m_typeName; };
......
......@@ -16,7 +16,7 @@ theCaloClusterRetriever = JiveXML__CaloClusterRetriever (name = "CaloClusterRetr
#theCaloClusterRetriever.DoWriteHLT = True
## Default collection (most Electron have elementLink to this one):
theCaloClusterRetriever.FavouriteClusterCollection="egClusterCollection"
theCaloClusterRetriever.FavouriteClusterCollection="egammaTopoCluster"
## example how to set other collection. when commented out: all other, non-HLT
##
......
......@@ -3,10 +3,6 @@
*/
#include "CaloJiveXML/CaloClusterRetriever.h"
#include "CaloEvent/CaloClusterContainer.h"
#include "CaloEvent/CaloCell.h"
#include "AthenaKernel/Units.h"
using Athena::Units::GeV;
......@@ -43,8 +39,8 @@ namespace JiveXML {
if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "in retrieveAll()" << endreq;
const DataHandle<CaloClusterContainer> iterator, end;
const CaloClusterContainer* ccc;
const DataHandle<xAOD::CaloClusterContainer> iterator, end;
const xAOD::CaloClusterContainer* ccc;
//obtain the default collection first
if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Trying to retrieve " << dataTypeName() << " (" << m_sgKeyFavourite << ")" << endreq;
......@@ -91,7 +87,7 @@ namespace JiveXML {
//obtain all collections with keys provided by user: m_otherKeys
std::vector<std::string>::const_iterator keyIter,endIter;
for ( keyIter=m_otherKeys.begin(); keyIter!=m_otherKeys.end(); ++keyIter ){
if ( evtStore()->contains<CaloClusterContainer>(*keyIter) ){ // to avoid some SG dumps
if ( evtStore()->contains<xAOD::CaloClusterContainer>(*keyIter) ){ // to avoid some SG dumps
if ( !evtStore()->retrieve( ccc, (*keyIter) ).isFailure()) {
if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Trying to retrieve selected " << dataTypeName() << " (" << (*keyIter) << ")" << endreq;
DataMap data = getData(ccc);
......@@ -115,7 +111,7 @@ namespace JiveXML {
* back-navigation causes Athena crash).
* @param FormatTool the tool that will create formated output from the DataMap
*/
const DataMap CaloClusterRetriever::getData(const CaloClusterContainer* ccc) {
const DataMap CaloClusterRetriever::getData(const xAOD::CaloClusterContainer* ccc) {
if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "getData()" << endreq;
......@@ -128,11 +124,11 @@ namespace JiveXML {
DataVect numCellsVec; numCellsVec.reserve(ccc->size());
std::string tagCells;
CaloClusterContainer::const_iterator itr = ccc->begin();
xAOD::CaloClusterContainer::const_iterator itr = ccc->begin();
int noClu = ccc->size();
int noCells = 0;
for (; itr != ccc->end(); ++itr){
for(CaloCluster::cell_iterator it = (*itr)->cell_begin();
for(xAOD::CaloCluster::const_cell_iterator it = (*itr)->cell_begin();
it != (*itr)->cell_end() ; ++it){
++noCells;
}
......@@ -153,7 +149,7 @@ namespace JiveXML {
idVec.push_back(DataType( ++id ));
int numCells = 0;
CaloCluster::cell_iterator cell = (*itr)->cell_begin();
xAOD::CaloCluster::const_cell_iterator cell = (*itr)->cell_begin();
for(; cell != (*itr)->cell_end() ; ++cell){
cells.push_back(DataType((*cell)->ID().get_compact()));
++numCells;
......
# Running the code
Using a single 2017 run as an illustrative example:
```
source setup.sh
infile="/eos/user/m/miokeefe/public/FullRun2_Reproc/MergedOutputs/HighMu/data17_13TeV/tree_340030.root"
grl="/cvmfs/atlas.cern.ch/repo/sw/database/GroupData/GoodRunsLists/data17_13TeV/20190708/data17_13TeV.periodAllYear_DetStatus-v105-pro22-13_Unknown_PHYS_StandardGRL_All_Good_25ns_Triggerno17e33prim.xml"
outdir="/afs/cern.ch/user/m/miokeefe/public/Zcounting/CSVOutputs/HighMu/data17_13TeV"
python -u dqt_zlumi_pandas.py --dblivetime --useofficial --infile $infile --grl $grl --campaign mc16d --outdir $out_dir
```
# Making single run plots
```
# Time dependent efficiency and luminosity plots
python plotting/efficiency.py --infile ~/public/Zcounting/CSVOutputs/HighMu/data17_13TeV/run_340030.csv
python plotting/luminosity.py --infile ~/public/Zcounting/CSVOutputs/HighMu/data17_13TeV/run_340030.csv
# Pileup dependent efficiency and luminosity plots
python plotting/efficiency.py --infile ~/public/Zcounting/CSVOutputs/HighMu/data17_13TeV/run_340030.csv --usemu
python plotting/luminosity.py --infile ~/public/Zcounting/CSVOutputs/HighMu/data17_13TeV/run_340030.csv --usemu
```
# Running over the entire dataset
```
./run_code.sh 17 local # to run jobs locally in a loop
./run_code.sh 17 batch # to submit all jobs to the batch
```
# Making plots for the entire dataset
```
for year in 15 16 17 18
do
for channel in Zee Zmumu Zll
do
python plotting/yearwise_luminosity.py --year $year --channel $channel
done
done
```
# Making full Run-2 "super plots"
```
for channel in Zee Zmumu Zll
do
python plotting/yearwise_luminosity.py --year all --channel $channel
done
```
#!/usr/bin/env python
# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
import numpy as np
import csv
from array import array
import ROOT
import sys, os
import logging
import argparse
import math
from DQUtils import fetch_iovs
from DQUtils.iov_arrangement import inverse_lblb
import tools.zlumi_alleff as dq_eff
import tools.zlumi_mc_cf as dq_cf
# testing toy sampling
import tools.toys as toys
ROOT.gROOT.SetBatch(ROOT.kTRUE)
ROOT.gStyle.SetOptStat(0)
logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.INFO)
parser = argparse.ArgumentParser()
parser.add_argument('--infile', type=str, help='input HIST file')
parser.add_argument('--grl', type=str, help='Specify an input GRL')
parser.add_argument('--tag', type=str, help='Lumi tag', default='OflLumiAcct-001')
parser.add_argument('--useofficial', action='store_true', help='Use official lumi folder (otherwise, use OflLumiAcct')
parser.add_argument('--lumifolder', type=str, help='Lumi folder', default='/TRIGGER/OFLLUMI/OflPrefLumi')
parser.add_argument('--lumitag', type=str, help='Lumi tag', default='OflLumi-13TeV-010')
parser.add_argument('--plotdir', type=str, help='Directory to dump plots', default='plots')
parser.add_argument('--outdir', type=str, help='Directory to dump plots', default='plots')
parser.add_argument('--dblivetime', action='store_true', help='Look up livetime from DB')
parser.add_argument('--campaign', type=str, help='mc16a/d/e')
args = parser.parse_args()
campaign = args.campaign
BINWIDTH = 10
ZPURITYFACTOR = 0.9935
ZXSEC = 1.929
ZATIMESC = 0.29632
ntoys = 10000000
do_toys = False
fin = ROOT.TFile.Open(args.infile)
runname = None
for key in fin.GetListOfKeys():
if key.GetName().startswith('run_'):
runname = key.GetName()
break
if args.grl:
grlReader = ROOT.Root.TGoodRunsListReader(args.grl)
grlReader.Interpret()
grl = grlReader.GetMergedGRLCollection()
else:
grl = None
if not runname:
logging.critical("Can't find run_* directory in input file %s", args.infile)
sys.exit(1)
lb_length = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_lblength_lb' % runname)
lbmin, lbmax = lb_length.GetXaxis().GetXmin(), lb_length.GetXaxis().GetXmax()
logging.info('low, high LBs: %s, %s', lbmin, lbmax)
if args.dblivetime:
logging.info('Starting livetime lookup ... (remove when we have a proper in-file implementation ...)')
livetime = ROOT.TProfile('livetime', 'Livetime', int(lbmax-lbmin), lbmin, lbmax)
else:
livetime = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_livetime_lb' % runname)
official_lum_zero = ROOT.TProfile('official_lum_zero', 'official inst luminosity', int(lbmax-lbmin), lbmin, lbmax)
official_mu = ROOT.TProfile('official_mu', 'official mu', int(lbmax-lbmin), lbmin, lbmax)
lblb = fetch_iovs("LBLB", runs=int(runname[4:]))
lbtime = inverse_lblb(lblb)
iovs_acct = fetch_iovs('COOLOFL_TRIGGER::/TRIGGER/OFLLUMI/LumiAccounting', lbtime.first.since, lbtime.last.until, tag=args.tag)
if args.useofficial:
print("Using official lumitag", args.lumitag)
iovs_lum = fetch_iovs('COOLOFL_TRIGGER::%s' % args.lumifolder, lblb.first.since, lblb.last.until, tag=args.lumitag, channels=[0])
lb_start_end = {}
lb_lhcfill = {}
for iov in lblb:
lb_start_end[iov.since & 0xffffffff] = (iov.StartTime/1e9, iov.EndTime/1e9)
#print(iov.StartTime/1e9, iov.EndTime/1e9)
for iov in iovs_acct:
if not lbmin < iov.LumiBlock < lbmax:
continue
lb_lhcfill[iov.LumiBlock] = iov.FillNumber
if args.dblivetime:
livetime.Fill(iov.LumiBlock, iov.LiveFraction)
if not args.useofficial:
official_lum_zero.Fill(iov.LumiBlock, iov.InstLumi/1e3)
official_mu.Fill(iov.LumiBlock, iov.AvEvtsPerBX)
else:
offlumiov = [_ for _ in iovs_lum if _.since.lumi==iov.LumiBlock]
if len(offlumiov) != 1:
print('MAJOR PROBLEM, LUMI IOV MISMATCH', iov.LumiBlock)
continue
offlumiov = offlumiov[0]
official_lum_zero.Fill(iov.LumiBlock, offlumiov.LBAvInstLumi/1e3)
official_mu.Fill(iov.LumiBlock, offlumiov.LBAvEvtsPerBX)
lb_full = lb_length.Clone('lb_full').ProjectionX()
divisor = lb_length.Clone('divisor').ProjectionX()
px = livetime.ProjectionX()
divisor.Multiply(px)
# Get run-wise electron channel histos outside of loop
# Also do the container efficiency bkg fit here as it is done
# once per run, then normalised to the ratio of LB/run in the function
# container_efficiency
hto = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_ele_template_os' % (runname))
hts = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_ele_template_ss' % (runname))
hphotontotal = fin.Get("%s/GLOBAL/DQTGlobalWZFinder/m_elContainertp_nomatch" % (runname))
hphotontotal.GetXaxis().SetRangeUser(66000, 250000)
# ==== Set signal region == 0, then fit ====
h_fit = hphotontotal.Clone()
for xbin in range(1, h_fit.GetNbinsX()+1):
mass = h_fit.GetBinLowEdge(xbin)
if mass > 75000 and mass < 100000:
h_fit.SetBinContent(xbin, 0)
h_fit.SetBinError(xbin, 0)
#elif mass > 150000:
# h_fit.SetBinContent(xbin, 0)
# h_fit.SetBinError(xbin, 0)
h_fit.Fit("pol2", "q")
o_recoeff_fit = {}
o_recoerr_fit = {}
# Do electron channel reco eff calculation once,
# then apply iterative correction to bkg subtraction in the main loop
lb_minus_one_reco_eff = [1.0, 1.0, 1.0]
for ibin in xrange(1, int(lbmax-lbmin)+1):
this_lb = int(lb_full.GetBinCenter(ibin))
lb = "lb_" + str(this_lb)
pileup = round(official_mu[ibin])
hmo = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_os' % (runname, lb))
hms = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_ss' % (runname, lb))
hno = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_os' % (runname, lb))
hns = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_ss' % (runname, lb))
try:
eff, err = dq_eff.template_method(hmo, hms, hno, hns, hto, hts, False, lb_minus_one_reco_eff) # do_scale set to false here as we are calculating it, when applying set True (next loop)
lb_minus_one_reco_eff = [eff, err, this_lb]
except AttributeError:
eff = 0
err = 0
if err != 0 and eff != 0:
weight = 1/pow(err, 2)
if pileup not in o_recoeff_fit:
o_recoeff_fit[pileup] = weight * eff
o_recoerr_fit[pileup] = weight
else:
o_recoeff_fit[pileup] += weight * eff
o_recoerr_fit[pileup] += weight
arr_rec_eff = array('d')
arr_rec_err = array('d')
arr_mu = array('d')
for pileup in o_recoeff_fit:
arr_mu.append(float(pileup))
arr_rec_eff.append(o_recoeff_fit[pileup]/o_recoerr_fit[pileup])
arr_rec_err.append(1/pow(o_recoerr_fit[pileup], 0.5))
tg_fit = ROOT.TGraphErrors(len(arr_mu), arr_mu, arr_rec_eff, ROOT.nullptr, arr_rec_err)
if len(o_recoeff_fit) == 1:
fit_type = "pol0"
elif len(o_recoeff_fit) == 2:
fit_type = "pol1"
elif len(o_recoeff_fit) > 2:
fit_type = "pol2"
print("Fit type", fit_type, "pileup bins", len(o_recoeff_fit))
tg_fit.Fit(fit_type, "q")
if args.outdir:
out_dir = args.outdir
os.system("mkdir -p " + out_dir)
out_dir += "/" + runname + ".csv"
else:
out_dir = runname + ".csv"
csvfile = open(out_dir, 'w')
csvwriter = csv.writer(csvfile, delimiter=',')
csvwriter.writerow(['FillNum','RunNum','LBNum','LBStart','LBEnd','LBLive','LBFull','OffLumi','OffMu', 'PassGRL',
'ZeeRaw','ZeeRawErr','ZeeN1','ZeeN2','ZeeEffTrig','ZeeErrTrig','ZeeEffReco','ZeeErrReco','ZeeEffComb','ZeeErrComb','ZeeEffAComb','ZeeErrAComb','ZeeDefTrig','ZeeDefReco','ZeeLumi','ZeeLumiErr',
'ZmumuRaw','ZmumuRawErr','ZmumuN1','ZmumuN2','ZmumuEffTrig','ZmumuErrTrig','ZmumuEffReco','ZmumuErrReco','ZmumuEffComb','ZmumuErrComb','ZmumuEffAComb','ZmumuErrAComb',
'ZmumuDefTrig','ZmumuDefReco','ZmumuLumi','ZmumuLumiErr',
'ZllLumi', 'ZllLumiErr'])
lb_minus_one_reco_eff = {}
lb_minus_one_reco_eff["Zee"] = [1.0, 1.0, 1]
lb_minus_one_reco_eff["Zmumu"] = [1.0, 1.0, 1]
lb_minus_one_trig_eff = {}
lb_minus_one_trig_eff["Zee"] = [1.0, 1.0, 1]
lb_minus_one_trig_eff["Zmumu"] = [1.0, 1.0, 1]
for ibin in xrange(1, int(lbmax-lbmin)+1):
out_dict = {}
out_dict["Zee"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
out_dict["Zmumu"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
this_lb = int(lb_full.GetBinCenter(ibin))
try:
this_fill = lb_lhcfill[this_lb]
except KeyError:
print("Fill not set. NA.")
this_fill = "NA"
passgrl = 1
for channel in ["Zee", "Zmumu"]:
if grl and not grl.HasRunLumiBlock(int(runname[4:]), this_lb):
passgrl = 0
continue
lb = "lb_" + str(this_lb)
if channel == 'Zee':
# Nominal signal histogram
hname = runname + '/lb_'+str(int(ibin+lbmin-0.5))+'/GLOBAL/DQTGlobalWZFinder/m_Z_mass_opsele'
# Tag-and-probe histos
h = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_eltrigtp_matches_os' % (runname, lb))
hmo = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_os' % (runname, lb))
hms = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_ss' % (runname, lb))
hno = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_os' % (runname, lb))
hns = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_ss' % (runname, lb))
#hphoton = fin.Get("%s/%s/GLOBAL/DQTGlobalWZFinder/m_elContainertp_nomatch" % (runname, lb))
#hpass = fin.Get("%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_passkine" % (runname, lb))
#hphoton.GetXaxis().SetRangeUser(66000, 250000)
#hpass.GetXaxis().SetRangeUser(66000, 250000)
ACCEPTANCE = 0.2996
elif channel == 'Zmumu':
# Nominal signal histogram
hname = runname + '/lb_'+str(int(ibin+lbmin-0.5))+'/GLOBAL/DQTGlobalWZFinder/m_Z_mass_opsmu'
# Tag-and-probe histos
h = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_mutrigtp_matches' % (runname, lb))
hmo = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_os' % (runname, lb))
hms = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_ss' % (runname, lb))
hno = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_os' % (runname, lb))
hns = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_ss' % (runname, lb))
ACCEPTANCE = 0.3323224
try:
z_m = fin.Get(hname).Integral()
z_merr = math.sqrt(z_m)
except AttributeError:
logging.error('Something unfortunate has happened; LB %d missing from Z count' % (ibin + lbmin - 0.5))
continue
N1 = N2 = 0
try:
N1 = h.GetBinContent(2)
N2 = h.GetBinContent(3)
if do_toys:
eff_trig, err_trig, arr_trig, arr_NZ = toys.toy_trigeff(N1, N2, ntoys)
else:
eff_trig, err_trig = dq_eff.trig_tag_and_probe(h, lb_minus_one_trig_eff[channel])
except TypeError:
eff_trig, err_trig = 0.0, 0.0
except AttributeError:
eff_trig, err_trig = 0.0, 0.0
if channel == "Zmumu":
try:
bin1 = hmo.GetXaxis().FindBin(86000)
bin2 = hmo.GetXaxis().FindBin(95000)
if do_toys:
eff_reco, err_reco, arr_reco = toys.muon_toy_recoeff(hmo.Integral(bin1, bin2), hms.Integral(bin1, bin2), hno.Integral(bin1, bin2), hns.Integral(bin1, bin2), ntoys)
else:
eff_reco, err_reco = dq_eff.reco_tag_and_probe(hmo, hms, hno, hns, lb_minus_one_reco_eff[channel])
except TypeError:
eff_reco, err_reco = 0.0, 0.0
except AttributeError:
eff_reco, err_reco = 0.0, 0.0
elif channel == "Zee":
try:
if do_toys:
bin1 = hmo.GetXaxis().FindBin(75000)
bin2 = hmo.GetXaxis().FindBin(104000)
bin3 = hmo.GetXaxis().FindBin(120000)
bin4 = hmo.GetXaxis().FindBin(250000)
matchos_peak = hmo.Integral(bin1, bin2)
matchos_tail = hmo.Integral(bin3, bin4)
matchss_tail = hms.Integral(bin3, bin4)
nomatchos_peak = hno.Integral(bin1, bin2)
nomatchos_tail = hno.Integral(bin3, bin4)
templateos_peak = hto.Integral(bin1, bin2)
templateos_tail = hto.Integral(bin3, bin4)
templatess_tail = hts.Integral(bin3, bin4)
eff_reco, err_reco, arr_reco = toys.electron_toy_recoeff(matchos_peak, matchos_tail, matchss_tail, nomatchos_peak, nomatchos_tail, templateos_peak, templateos_tail, templatess_tail, ntoys)
else:
eff_reco, err_reco = dq_eff.template_method(hmo, hms, hno, hns, hto, hts, True, lb_minus_one_reco_eff[channel], tg_fit.GetFunction(fit_type).Eval(official_mu[ibin]))
except TypeError:
eff_reco, err_reco = 0.0, 0.0
except AttributeError:
eff_reco, err_reco = 0.0, 0.0
#container_efficiency(hphoton, hphotontotal, h_fit, hpass, hto, hts)
defaulted_reco_eff = 0
defaulted_trig_eff = 0
if eff_reco == lb_minus_one_reco_eff[channel][0]:
defaulted_reco_eff = 1
if eff_trig == lb_minus_one_trig_eff[channel][0]:
defaulted_trig_eff = 1
if eff_reco != 0.0:
lb_minus_one_reco_eff[channel] = [eff_reco, err_reco, this_lb]
if eff_trig != 0.0:
lb_minus_one_trig_eff[channel] = [eff_trig, err_trig, this_lb]
if do_toys:
arr_comb = (1.0 - (1.0 - arr_trig)**2) * (arr_reco)**2
nonan_arr_comb = arr_comb[~np.isnan(arr_comb)]
eff_comb = np.median(nonan_arr_comb)
err_comb = nonan_arr_comb.std()
else:
eff_comb = (1-(1-eff_trig)**2)*(eff_reco)**2
err_comb = ((eff_reco**2*2*(1-eff_trig)*err_trig)**2+(2*eff_reco*(1-(1-eff_trig)**2)*err_reco)**2)**.5
eff_Acomb = ACCEPTANCE * eff_comb
err_Acomb = ACCEPTANCE * err_comb
lblive = divisor[ibin]
if lblive < 10:
continue
loclivetime = lblive
run = int(runname.replace("run_", ""))
if do_toys:
effcy = arr_comb * dq_cf.correction(official_mu[ibin], channel, campaign)
else:
effcy = eff_comb * dq_cf.correction(official_mu[ibin], channel, campaign)
effcyerr = err_comb * dq_cf.correction(official_mu[ibin], channel, campaign)
if run == 302831 and this_lb < 11:
loclivetime = 0
elif run == 329835 and this_lb < 554:
loclivetime = 0
elif run == 310247 and (this_lb == 442 or this_lb == 462):
loclivetime = 0
elif run == 281385 and this_lb < 197:
loclivetime = lblive * 4.0/6.0 # ad-hoc scale factor
elif run == 281385 and this_lb < 375:
loclivetime = lblive * 5.0/6.0 # ad-hoc scale factor
elif run == 286367:
loclivetime = lblive * 5.0/6.0 # ad-hoc scale factor
else:
loclivetime = lblive
zlumi = zlumistat = 0.0
if do_toys and loclivetime != 0.0:
arr_zlumi = np.divide(arr_NZ, effcy) * (ZPURITYFACTOR/ACCEPTANCE/ZXSEC)/loclivetime
arr_zlumi = arr_zlumi[~np.isnan(arr_zlumi)]
zlumi = np.median(arr_zlumi)
zlumistat = arr_zlumi.std()
elif (loclivetime != 0.0 and effcy != 0.0):
zlumi = (z_m/effcy)*(ZPURITYFACTOR/ACCEPTANCE/ZXSEC)/loclivetime
zlumistat = math.sqrt(pow(z_merr/effcy, 2) + pow(z_m/effcy**2*effcyerr, 2))*ZPURITYFACTOR/ACCEPTANCE/ZXSEC/loclivetime
out_dict[channel] = [z_m, z_merr, N1, N2, eff_trig, err_trig, eff_reco, err_reco, eff_comb, err_comb, eff_Acomb, err_Acomb, defaulted_trig_eff, defaulted_reco_eff, zlumi, zlumistat]
lumi_index = len(out_dict['Zee'])-2
error_index = len(out_dict['Zee'])-1
zll_lumi = (out_dict['Zee'][lumi_index] + out_dict['Zmumu'][lumi_index])/2
zll_lumi_err = 0.5 * math.sqrt( pow(out_dict['Zee'][error_index], 2) + pow(out_dict['Zmumu'][error_index], 2) )
out_write = [this_fill, run, this_lb, lb_start_end[this_lb][0], lb_start_end[this_lb][1], loclivetime, lb_full[ibin], official_lum_zero[ibin], official_mu[ibin], passgrl] + out_dict["Zee"] + out_dict["Zmumu"] + [zll_lumi, zll_lumi_err]
csvwriter.writerow(out_write)
csvfile.close()
#include "AtlasLabels.h"
#include "TLatex.h"
#include "TLine.h"
#include "TPave.h"
#include "TPad.h"
#include "TMarker.h"
void ATLASLabel(Double_t x,Double_t y,const char* text,Color_t color)
{
TLatex l; //l.SetTextAlign(12); l.SetTextSize(tsize);
l.SetNDC();
l.SetTextFont(72);
l.SetTextColor(color);
double delx = 0.115*696*gPad->GetWh()/(472*gPad->GetWw());
l.DrawLatex(x,y,"ATLAS");
if (text) {
TLatex p;
p.SetNDC();
p.SetTextFont(42);
p.SetTextColor(color);
p.DrawLatex(x+delx,y,text);
// p.DrawLatex(x,y,"#sqrt{s}=900GeV");
}
}
void ATLASLabelOld(Double_t x,Double_t y,bool Preliminary,Color_t color)
{
TLatex l; //l.SetTextAlign(12); l.SetTextSize(tsize);
l.SetNDC();
l.SetTextFont(72);
l.SetTextColor(color);
l.DrawLatex(x,y,"ATLAS");
if (Preliminary) {
TLatex p;
p.SetNDC();
p.SetTextFont(42);
p.SetTextColor(color);
p.DrawLatex(x+0.115,y,"Preliminary");
}
}
void ATLASVersion(const char* version,Double_t x,Double_t y,Color_t color)
{
if (version) {
char versionString[100];
sprintf(versionString,"Version %s",version);
TLatex l;
l.SetTextAlign(22);
l.SetTextSize(0.04);
l.SetNDC();
l.SetTextFont(72);
l.SetTextColor(color);
l.DrawLatex(x,y,versionString);
}
}
//
// @file AtlasLabels.h
//
// @author M.Sutton
//
// Copyright (C) 2010 Atlas Collaboration
//
// $Id: AtlasLabels.h, v0.0 Thu 25 Mar 2010 10:34:20 CET $
#ifndef __ATLASLABELS_H
#define __ATLASLABELS_H
#include "Rtypes.h"
void ATLASLabel(Double_t x,Double_t y,const char* text=NULL,Color_t color=1);
void ATLASLabelOld(Double_t x,Double_t y,bool Preliminary=false,Color_t color=1);
void ATLASVersion(const char* version=NULL,Double_t x=0.88,Double_t y=0.975,Color_t color=1);
#endif // __ATLASLABELS_H
//
// ATLAS Style, based on a style file from BaBar
//
#include <iostream>
#include "AtlasStyle.h"
#include "TROOT.h"
void SetAtlasStyle ()
{
static TStyle* atlasStyle = 0;
std::cout << "\nApplying ATLAS style settings...\n" << std::endl ;
if ( atlasStyle==0 ) atlasStyle = AtlasStyle();
gROOT->SetStyle("ATLAS");
gROOT->ForceStyle();
}
TStyle* AtlasStyle()
{
TStyle *atlasStyle = new TStyle("ATLAS","Atlas style");
// use plain black on white colors
Int_t icol=0; // WHITE
atlasStyle->SetFrameBorderMode(icol);
atlasStyle->SetFrameFillColor(icol);
atlasStyle->SetCanvasBorderMode(icol);
atlasStyle->SetCanvasColor(icol);
atlasStyle->SetPadBorderMode(icol);
atlasStyle->SetPadColor(icol);
atlasStyle->SetStatColor(icol);
//atlasStyle->SetFillColor(icol); // don't use: white fill color for *all* objects
// set the paper & margin sizes
atlasStyle->SetPaperSize(20,26);
// set margin sizes
atlasStyle->SetPadTopMargin(0.05);
atlasStyle->SetPadRightMargin(0.05);
atlasStyle->SetPadBottomMargin(0.16);
atlasStyle->SetPadLeftMargin(0.16);
// set title offsets (for axis label)
atlasStyle->SetTitleXOffset(1.4);
atlasStyle->SetTitleYOffset(1.4);
// use large fonts
//Int_t font=72; // Helvetica italics
Int_t font=42; // Helvetica
Double_t tsize=0.05;
atlasStyle->SetTextFont(font);
atlasStyle->SetTextSize(tsize);
atlasStyle->SetLabelFont(font,"x");
atlasStyle->SetTitleFont(font,"x");
atlasStyle->SetLabelFont(font,"y");
atlasStyle->SetTitleFont(font,"y");
atlasStyle->SetLabelFont(font,"z");
atlasStyle->SetTitleFont(font,"z");
atlasStyle->SetLabelSize(tsize,"x");
atlasStyle->SetTitleSize(tsize,"x");
atlasStyle->SetLabelSize(tsize,"y");
atlasStyle->SetTitleSize(tsize,"y");
atlasStyle->SetLabelSize(tsize,"z");
atlasStyle->SetTitleSize(tsize,"z");
// use bold lines and markers
atlasStyle->SetMarkerStyle(20);
atlasStyle->SetMarkerSize(1.2);
atlasStyle->SetHistLineWidth(2.);
atlasStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes
// get rid of X error bars
//atlasStyle->SetErrorX(0.001);
// get rid of error bar caps
atlasStyle->SetEndErrorSize(0.);
// do not display any of the standard histogram decorations
atlasStyle->SetOptTitle(0);
//atlasStyle->SetOptStat(1111);
atlasStyle->SetOptStat(0);
//atlasStyle->SetOptFit(1111);
atlasStyle->SetOptFit(0);
// put tick marks on top and RHS of plots
atlasStyle->SetPadTickX(1);
atlasStyle->SetPadTickY(1);
return atlasStyle;
}
//
// @file AtlasStyle.h
//
// ATLAS Style, based on a style file from BaBar
//
//
// @author M.Sutton
//
// Copyright (C) 2010 Atlas Collaboration
//
// $Id: AtlasStyle.h, v0.0 Thu 25 Mar 2010 10:34:20 CET $
#ifndef __ATLASSTYLE_H
#define __ATLASSTYLE_H
#include "TStyle.h"
void SetAtlasStyle();
TStyle* AtlasStyle();
#endif // __ATLASSTYLE_H
from ROOT import *
ROOT.gROOT.LoadMacro("AtlasStyle.C")
#SetAtlasStyle()
#include <iostream>
#include <cmath>
#include "AtlasUtils.h"
#include "TLine.h"
#include "TLatex.h"
#include "TMarker.h"
#include "TPave.h"
#include "TH1.h"
void ATLAS_LABEL(Double_t x,Double_t y,Color_t color)
{
TLatex l; //l.SetTextAlign(12); l.SetTextSize(tsize);
l.SetNDC();
l.SetTextFont(72);
l.SetTextColor(color);
l.DrawLatex(x,y,"ATLAS");
}
TGraphErrors* myTGraphErrorsDivide(TGraphErrors* g1,TGraphErrors* g2) {
const Int_t debug=0;
if (!g1) printf("**myTGraphErrorsDivide: g1 does not exist ! \n");
if (!g2) printf("**myTGraphErrorsDivide: g2 does not exist ! \n");
Int_t n1=g1->GetN();
Int_t n2=g2->GetN();
if (n1!=n2) {
printf("**myTGraphErrorsDivide: vector do not have same number of entries ! \n");
}
TGraphErrors* g3= new TGraphErrors();
Double_t x1=0., y1=0., x2=0., y2=0.;
Double_t dx1=0.,dy1=0., dy2=0.;
Int_t iv=0;
for (Int_t i1=0; i1<n1; i1++) {
for (Int_t i2=0; i2<n2; i2++) {
//if (debug) printf("**myTGraphErrorsDivide: %d %d ! \n",i1,i2);
g1->GetPoint(i1,x1,y1);
g2->GetPoint(i2,x2,y2);
if (x1!=x2) {
//printf("**myTGraphErrorsDivide: %d x1!=x2 %f %f ! \n",iv,x1,x2);
}else{
//if (debug) printf("**myTGraphErrorsDivide: %d x1=x2 %f %f ! \n",iv,x1,x2);
dx1 = g1->GetErrorX(i1);
if (y1!=0) dy1 = g1->GetErrorY(i1)/y1;
if (y2!=0) dy2 = g2->GetErrorY(i2)/y2;
if (debug)
printf("**myTGraphErrorsDivide: %d x1=%f x2=%f y1=%f y2=%f \n",iv,x1,x2,y1,y2);
if (y2!=0.) g3->SetPoint(iv, x1,y1/y2);
else g3->SetPoint(iv, x1,y2);
Double_t e=0.;
if (y1!=0 && y2!=0) e=std::sqrt(dy1*dy1+dy2*dy2)*(y1/y2);
g3->SetPointError(iv,dx1,e);
if (debug) {
//Double_t g3y, g3x,g3e;
//g3->GetPoint(iv, g3y,g3x);
//g3e=g3->GetErrorY(iv);
//printf("%d g3y= %f g3e=%f \n",iv,g3y,g3e);
}
iv++;
}
// printf("**myTGraphErrorsDivide: ...next \n");
}
}
return g3;
}
TGraphAsymmErrors* myTGraphErrorsDivide(TGraphAsymmErrors* g1,TGraphAsymmErrors* g2) {
const Int_t debug=0;
TGraphAsymmErrors* g3= new TGraphAsymmErrors();
Int_t n1=g1->GetN();
Int_t n2=g2->GetN();
if (n1!=n2) {
printf(" vectors do not have same number of entries ! \n");
return g3;
}
Double_t x1=0., y1=0., x2=0., y2=0.;
Double_t dx1h=0., dx1l=0.;
Double_t dy1h=0., dy1l=0.;
Double_t dy2h=0., dy2l=0.;
Double_t* X1 = g1->GetX();
Double_t* Y1 = g1->GetY();
Double_t* EXhigh1 = g1->GetEXhigh();
Double_t* EXlow1 = g1->GetEXlow();
Double_t* EYhigh1 = g1->GetEYhigh();
Double_t* EYlow1 = g1->GetEYlow();
Double_t* X2 = g2->GetX();
Double_t* Y2 = g2->GetY();
Double_t* EXhigh2 = g2->GetEXhigh();
Double_t* EXlow2 = g2->GetEXlow();
Double_t* EYhigh2 = g2->GetEYhigh();
Double_t* EYlow2 = g2->GetEYlow();
for (Int_t i=0; i<g1->GetN(); i++) {
g1->GetPoint(i,x1,y1);
g2->GetPoint(i,x2,y2);
dx1h = EXhigh1[i];
dx1l = EXlow1[i];
if (y1!=0.) dy1h = EYhigh1[i]/y1;
else dy1h = 0.;
if (y2!=0.) dy2h = EYhigh2[i]/y2;
else dy2h = 0.;
if (y1!=0.) dy1l = EYlow1 [i]/y1;
else dy1l = 0.;
if (y2!=0.) dy2l = EYlow2 [i]/y2;
else dy2l = 0.;
//if (debug)
//printf("%d x1=%f x2=%f y1=%f y2=%f \n",i,x1,x2,y1,y2);
if (debug)
printf("%d dy1=%f %f dy2=%f %f sqrt= %f %f \n",i,dy1l,dy1h,dy2l,dy2h,
std::sqrt(dy1l*dy1l+dy2l*dy2l), std::sqrt(dy1h*dy1h+dy2h*dy2h));
if (y2!=0.) g3->SetPoint(i, x1,y1/y2);
else g3->SetPoint(i, x1,y2);
Double_t el=0.; Double_t eh=0.;
if (y1!=0. && y2!=0.) el=std::sqrt(dy1l*dy1l+dy2l*dy2l)*(y1/y2);
if (y1!=0. && y2!=0.) eh=std::sqrt(dy1h*dy1h+dy2h*dy2h)*(y1/y2);
if (debug) printf("dx1h=%f dx1l=%f el=%f eh=%f \n",dx1h,dx1l,el,eh);
g3->SetPointError(i,dx1h,dx1l,el,eh);
}
return g3;
}
TGraphAsymmErrors* myMakeBand(TGraphErrors* g0, TGraphErrors* g1,TGraphErrors* g2) {
// default is g0
//const Int_t debug=0;
TGraphAsymmErrors* g3= new TGraphAsymmErrors();
Double_t x1=0., y1=0., x2=0., y2=0., y0=0, x3=0.;
//Double_t dx1=0.;
Double_t dum;
for (Int_t i=0; i<g1->GetN(); i++) {
g0->GetPoint(i, x1,y0);
g1->GetPoint(i, x1,y1);
g2->GetPoint(i, x1,y2);
// if (y1==0) y1=1;
//if (y2==0) y2=1;
if (i==g1->GetN()-1) x2=x1;
else g2->GetPoint(i+1,x2,dum);
if (i==0) x3=x1;
else g2->GetPoint(i-1,x3,dum);
Double_t tmp=y2;
if (y1<y2) {y2=y1; y1=tmp;}
//Double_t y3=1.;
Double_t y3=y0;
g3->SetPoint(i,x1,y3);
Double_t binwl=(x1-x3)/2.;
Double_t binwh=(x2-x1)/2.;
if (binwl==0.) binwl= binwh;
if (binwh==0.) binwh= binwl;
g3->SetPointError(i,binwl,binwh,(y3-y2),(y1-y3));
}
return g3;
}
void myAddtoBand(TGraphErrors* g1, TGraphAsymmErrors* g2) {
Double_t x1=0., y1=0., y2=0., y0=0;
//Double_t dx1=0.;
//Double_t dum;
if (g1->GetN()!=g2->GetN())
std::cout << " graphs have not the same # of elements " << std::endl;
Double_t* EYhigh = g2-> GetEYhigh();
Double_t* EYlow = g2-> GetEYlow();
for (Int_t i=0; i<g1->GetN(); i++) {
g1->GetPoint(i, x1,y1);
g2->GetPoint(i, x1,y2);
if ( y1==0 || y2==0 ) {
std::cerr << "check these points very carefully : myAddtoBand() : point " << i << std::endl;
}
// if (y1==0) y1=1;
// if (y2==0) y2=1;
// if (i==g1->GetN()-1) x2=x1;
// else g2->GetPoint(i+1,x2,dum);
// if (i==0) x3=x1;
// else g2->GetPoint(i-1,x3,dum);
Double_t eyh=0., eyl=0.;
//if (y1<y2) {y2=y1; y1=tmp;}
//Double_t y3=1.;
//printf("%d: y1=%f y2=%f Eyhigh= %f Eylow= %f \n",i,y1,y2,EYhigh[i],EYlow[i]);
y0=y1-y2;
if (y0!=0) {
if (y0>0){
eyh=EYhigh[i];
eyh=std::sqrt(eyh*eyh+y0*y0);
//printf("high: %d: y0=%f eyh=%f \n",i,y0,eyh);
g2->SetPointEYhigh(i,eyh);
} else {
eyl=EYlow[i];
eyl=std::sqrt(eyl*eyl+y0*y0);
// printf("low: %d: y0=%f eyl=%f \n",i,y0,eyl);
g2->SetPointEYlow (i,eyl);
}
}
}
return;
}
TGraphErrors* TH1TOTGraph(TH1 *h1){
if (!h1) std::cout << "TH1TOTGraph: histogram not found !" << std::endl;
TGraphErrors* g1= new TGraphErrors();
Double_t x, y, ex, ey;
for (Int_t i=1 ; i<=h1->GetNbinsX(); i++) {
y=h1->GetBinContent(i);
ey=h1->GetBinError(i);
x=h1->GetBinCenter(i);
ex=h1->GetBinWidth(i);
// cout << " x,y = " << x << " " << y << " ex,ey = " << ex << " " << ey << endl;
g1->SetPoint(i-1,x,y);
g1->SetPointError(i-1,ex,ey);
}
//g1->Print();
return g1;
}
void myText(Double_t x,Double_t y,Color_t color, const char *text, bool set_text_size = false) {
//Double_t tsize=0.05;
TLatex l; //l.SetTextAlign(12); l.SetTextSize(tsize);
l.SetTextFont(43);
l.SetTextSize(27);
if(set_text_size != 0){
l.SetTextSize(39);
}
l.SetNDC();
l.SetTextColor(color);
l.DrawLatex(x,y,text);
}
void myBoxText(Double_t x, Double_t y,Double_t boxsize,Int_t mcolor,const char *text)
{
Double_t tsize=0.06;
TLatex l; l.SetTextAlign(12); //l.SetTextSize(tsize);
l.SetNDC();
l.DrawLatex(x,y,text);
Double_t y1=y-0.25*tsize;
Double_t y2=y+0.25*tsize;
Double_t x2=x-0.3*tsize;
Double_t x1=x2-boxsize;
printf("x1= %f x2= %f y1= %f y2= %f \n",x1,x2,y1,y2);