interactiveCode.txt

### Interactive Plotting Code

### List of Contents
# ttSF plots for AN
# comparing limits from different files
# run expected limits on LPC
# make QCDSF plots for AN from betrees_mc.root
# compare pre-fit and post-fit bkg component in any region
# get MC yields for all bkg processes
# get signal MC yields for all masses
# get 1l electron eta comparisons
# plot GoF distribution from combine procedure
# compare bbVV and bbtautau H->WW reco
# print out limits from file
# HEM dilepton correction plots
# Bkg MC plots for AN
# Lepton efficiency plots in signal

## ----------------------------------------------------------------------------------------------------------------------------------------
## tt SF AN plots
## ----------------------------------------------------------------------------------------------------------------------------------------
TFile *fm = new TFile("HHbb1o2l_ttbarSF_mc_inputPlots.root");
TFile *fd = new TFile("HHbb1o2l_ttbarSF_data_inputPlots.root");

Plotter *p = new Plotter();
p->addStackHist((TH1*)fm->Get("ttbar_noHbb_hbbMass"),"t#bar{t} MC");
p->addStackHist((TH1*)fm->Get("other_noHbb_hbbMass"),"non-t#bar{t} MC");
p->addHist((TH1*)fd->Get("data_noHbb_hbbMass"),"data");
p->draw(false,"s");

TFile *f = TFile::Open("HHbb1o2l_ttbarSF.json.root");
TGraphErrors *g;
f->GetObject("RATIO",g);
g->GetXaxis()->SetNameTitle("mhh","m_{HH} [GeV]");
g->GetYaxis()->SetNameTitle("ratio","(N_{data}-N_{non-tt})/N_{tt}");
TCanvas *cg = new TCanvas("cg","cg");
cg->cd();
g->Draw();

## ----------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Compare limits from different file
## ----------------------------------------------------------------------------------------------------------------------------------------
double limit, mx;
float quant;
Plotter *p = new Plotter();
float limfactor = 1000;

std::vector<TString> fnames = {
	"~/Dropbox/Physics/HHbbWW/betrees/master/limits_b2g18008_radion.root",
	"~/Dropbox/Physics/HHbbWW/betrees/master/limits_b2g18008_blkgrav.root",
	"~/Dropbox/Physics/HHbbWW/betrees/master/expLimits_rad_030221.root",
	"~/Dropbox/Physics/HHbbWW/betrees/master/expLimits_blk_030221.root",
};

std::vector<TString> ids = {"B2G-18-008 spin-0","B2G-18-008 spin-2","Current spin-0","Current spin-2"};

for(unsigned iF=0; iF<fnames.size(); ++iF){
if(iF < 2) {
float num = 0.5824*2*(0.2137+.02619);
float den = 2*0.5824*(.2137+.002619);
limfactor *= (num/den);
}
TFile *f = TFile::Open(fnames[iF]);
TTree *t = (TTree*)f->Get("limit");
t->SetBranchAddress("quantileExpected",&quant);
t->SetBranchAddress("limit",&limit);
t->SetBranchAddress("mh",&mx);

std::vector<std::pair<double,double>> pts;
for(unsigned i=0; i<t->GetEntries(); ++i) {
t->GetEntry(i);
if(quant != 0.5) continue;
pts.emplace_back(mx,limfactor*limit);
}
delete f;
std::sort(pts.begin(),pts.end(),[](std::pair<double,double> a, std::pair<double,double> b) {return a.first < b.first;});
std::vector<double> lims, masses;
for(const auto& pt : pts) {
masses.push_back(pt.first);
lims.push_back(pt.second);
}
TGraph *g = new TGraph(pts.size(),&masses[0],&lims[0]);
p->addGraph(g,ids[iF],-1,1,4,20,0.8,true,true,false,"L");
limfactor = 1000;
}
p->setXTitle("#it{m}_{X} [GeV]");
p->setYTitle("#sigma#it{B}(X#rightarrowHH) [fb^{-1}]");
p->draw(false,"gs");

## ----------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Run some expected limits
## ----------------------------------------------------------------------------------------------------------------------------------------
for mx in 800 900 1000 1200 1400 1600 1800 2000 2500 3000 3500
do
  nohup combine -m ${mx} -M AsymptoticLimits --run expected -v 2 --rAbsAcc=0.00001 combinedCard.txt &
done
## ----------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## QCDSF plots from betrees_mc.root
## ----------------------------------------------------------------------------------------------------------------------------------------
treeMaker->cd()
Events->Draw("hbbMass>>wjets(30,30,210)","xsec*btag_N*pu_N*trig_N*lep_N*fjbtag_N*(process==4&&hbbMass>=30&&hbbMass<=210&&hhMass>=700)")
Events->Draw("hbbMass>>qcd(30,30,210)","xsec*btag_N*pu_N*trig_N*lep_N*fjbtag_N*(process==8&&hbbMass>=30&&hbbMass<=210&&hhMass>=700)")
Plotter *p = new Plotter();
p->addHist(wjets,"W+jets");
p->addHist(qcd,"QCD multi-jet");
p->normalize();
p->setYTitle("a.u.");
p->setXTitle("m_{bb} [GeV]")
p->draw(false,"s");
## ----------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Compare pre- and post-fit bkg components
## ----------------------------------------------------------------------------------------------------------------------------------------

std::vector<TString> bkgs = {"mw"};
std::vector<TString> bins = {"OF_L"};
TString reg = "SR";
TFile *fpreInt = 0, *fpre = 0;

TFile *fpst = TFile::Open("092120_comb_NonTopCR_rad_sepHbbScale/postFit.root");

for(const auto& bkg : bkgs){
if(reg.Contains("CR")) {
fpreInt = TFile::Open("bkgInputs"+reg+"/HHbb1o2l_"+reg+_"+bkg+"_distributions.root")
fpre    = TFile::Open("bkgInputs"+reg+"/HHbb1o2l_"+reg+"_"+bkg+"_2D_template_debug.root");
}else
fpreInt = TFile::Open("bkgInputs/HHbb1o2l_"+bkg+"_distributions.root")
fpre    = TFile::Open("bkgInputs/HHbb1o2l_"+bkg+"_2D_template_debug.root");
}

for(const auto& bin : bins){
TH2 *hpst = (TH2*)fpst->Get("postfit_"+bkg+"_"+bin+"_full");
TH2 *mc = (TH2*)fpreInt->Get(bkg+"_"+bin+"_full_hbbMass_hhMass");
TH2 *hpre = (TH2*)fpre->Get(bkg+"_"+bin+"_full");
hpre->Scale(mc->Integral()/hpre->Integral());

TH1 *hpstX = (TH1*)hpst->ProjectionX("pstX",1,hpst->GetNbinsY());
TH1 *hpstY = (TH1*)hpst->ProjectionY("pstY",1,hpst->GetNbinsX())
TH1 *hpreX = (TH1*)hpre->ProjectionX("preX",1,hpre->GetNbinsY());
TH1 *hpreY = (TH1*)hpre->ProjectionY("preY",1,hpre->GetNbinsX());

Plotter *pX = new Plotter();
Plotter *pY = new Plotter();

pX->addHistLine(hpreX,"prefit");
pX->addHistLine(hpstX,"postfit");
pY->addHistLine(hpreY,"prefit");
pY->addHistLine(hpstY,"postfit");

pX->drawSplitRatio(0,"mbb",false,false,"mbb");
pY->drawSplitRatio(0,"mhh",false,false,"mhh");
}
}

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Electron eta sel in 1l channel for B2G
## ----------------------------------------------------------------------------------------------------------------------------------------

TString var = "mhh";
TString smp = "qcd";
TFile *f = TFile::Open(smp+"_etaSel.root");
std::vector<TString> etas = {"2p5","2p1","1p479"};

Plotter *p = new Plotter();
for(const auto& eta : etas){
TH1 *h = (TH1*)f->Get(smp+"_eta"+eta+"_"+var);
TString lab = "|#eta_{e}| < "+eta;
lab.ReplaceAll("p",".");
p->addHist(h,lab);
}
p->rebin(10);
p->drawSplitRatio(0,smp,false,false,smp);
## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Print Bkg MC yields for all regions
## ----------------------------------------------------------------------------------------------------------------------------------------
TFile *f = TFile::Open("debugRun2.root");
std::vector<TString> procs = {"ttbar0","ttbar1","ttbar2","singlet","wjets","zjets","diboson","qcd","ttx","hx","qg","tw","mw","mt"};
std::vector<TString> regs= {"SR1","SR2","TopCR1","TopCR2","NonTopCR1","NonTopCR2"};

std::cout<<"{";
for(const auto& p : procs) std::cout<<p<<", ";
std::cout<<std::endl;

for(const TString r : regs){
std::cout<<r<<std::endl;

float bkgN = 0.0;
TH1 *hb = (TH1*)f->Get("bkg_"+r+"__mbb");
bkgN = hb->Integral(); std::cout<<std::endl; delete hb;

std::cout<<"bkg = "<<bkgN<<std::endl;

for(TString p : procs){
TString s = p+"_"+r+"__mbb";
TH1 *h = (TH1*)f->Get(p+"_"+r+"__mbb"); float procN = 0.0;
if(h) {
procN = h->Integral();
}
std::cout<<procN/bkgN<<std::endl;
}
std::cout<<std::endl;

}
## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Print signal yields for all masses and spins
## ----------------------------------------------------------------------------------------------------------------------------------------

std::vector<TString> sigs = {"radHH","blkHH"};
std::vector<TString> sels = {"emu_LMT_I","IF_LMT"};
std::vector<TString> masses = {"800","900","1000","1200","1400","1600","1800","2000","2500","3000","3500","4000","4500"};

for(const auto& sig : sigs) {
TFile *f = TFile::Open("signalInputs/HHbb1o2l_"+sig+"_exclM_distributions.root");
std::cout<<sig<<std::endl;
for(const auto& sel : sels) {
std::cout<<sel<<std::endl;
for(const auto& m : masses) {
TH2 *h = (TH2*)f->Get(sig+"_m"+m+"_"+sel+"_full_hbbMass_hhMass");
std::cout<<h->Integral()<<std::endl;
}
}

}
std::cout<<std::endl;

## ------------------------------------------------------------------------------------------------------------------------------------

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Plots for 1l electron eta selection studies
## ----------------------------------------------------------------------------------------------------------------------------------------

std::vector<TString> procs = {"qcd","radion","blkgrav"};
std::vector<TString> etaSs = {"2p5","2p1","1p479"};

for(const auto& proc : procs){
Plotter *p = new Plotter();
TFile *f = TFile::Open(proc+"_etaSel.root");
for(unsigned i=0;i<etaSs.size();i++){
TString eS = etaSs[i];
TH1 *h = (TH1*)f->Get(proc+"_eta"+eS+(proc=="qcd"?TString("_mhh"):TString("_mx")));
eS.ReplaceAll("p",".");
p->addHist(h,"|#eta| #leq "+eS);
}
p->setYTitle("Events");
p->setXTitle(proc=="qcd"?"M_{HH}":"M_{X}");
if(proc=="qcd") p->rebin(10);
p->setYTitleBot("Ratio");
p->drawSplitRatio(0,proc,false,false,proc);
}

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## GoF distribution w/ observed from the combine tool
## ----------------------------------------------------------------------------------------------------------------------------------------

TFile *fd = TFile::Open("higgsCombineTest.GoodnessOfFit.mH1000.root");
TFile *ft = TFile::Open("higgsCombineTest.GoodnessOfFit.mH1000.toys.root");

TTree *td = (TTree*)fd->Get("limit");
TTree *tt = (TTree*)ft->Get("limit");

double tsData, tsT;
double pval = 0.0;
td->SetBranchAddress("limit",&tsData);
tt->SetBranchAddress("limit",&tsT);
td->GetEntry(0);
std::vector<double> tsTs;

for(int i=0;i<tt->GetEntries();i++){
tt->GetEntry(i);
if(tsT > tsData) pval += 1.0;
tsTs.push_back(tsT);
}

std::sort(tsTs.begin(),tsTs.end(),[](double a, double b){return a < b;});

TH1D* h = new TH1D("h","h",40,tsTs.front()-100,tsTs.back()+100);
for(const auto& t : tsTs) h->Fill(t);

Plotter *p = new Plotter();
p->addStackHist(h,"toys");
p->setXTitle("GoF statistic");
p->setYTitle("N_{toys}");
TArrow *arr = new TArrow(tsData,h->GetMaximum(),tsData,0);
arr->SetLineWidth(2);
arr->SetLineColor(kRed);
//p->addLegendEntry(0,arr,"data");
TCanvas *c = p->draw(false,"s");
arr->Draw();


std::cout<<"data ts = "<<tsData<<std::endl;
std::cout<<"pval is = "<<pval/tt->GetEntries()<<std::endl;

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Compare H->WW reco. between bbVV and bbtautau signals
## ----------------------------------------------------------------------------------------------------------------------------------------

std::vector<TString> sigs = {"hww2l.root","hww2l_bbVV.root","hww2l_bbtt.root"};
std::vector<TString> sigNames = {"Inclusive","bbWW","bb#tau#tau"};

std::vector<TString> vars = {"mInv"};
TString hid = "base";
TString sig = "radion_m1000";

for(const auto& v : vars){
Plotter *p= new Plotter();

for(unsigned i=0;i<sigs.size();++i){
TFile *f = TFile::Open(sigs[i]);
TH1 *h = (TH1*)f->Get(sig+"_"+hid+"_"+v);
p->addHistLine(h,sigNames[i],-1,1,5);
}
p->rebin(5);
p->normalize();
p->setYTitle("a.u.");
p->draw(false,v);
}

sig = "radion";
Plotter *pp = new Plotter();
std::vector<TString> ms = {"m1000","m3000"};
std::vector<TString> massnames = {"1.0 TeV","3.0 TeV"};
for(unsigned i=1;i<=2;++i) {
TFile *f = TFile::Open(sigs[i]);
for(unsigned j=0;j<=1;++j) {
TH1 *h = (TH1*)f->Get(sig+"_"+ms[j]+"_"+hid+"_dTheta");
pp->addHistLine(h,sigNames[i]+" "+massnames[j]);
}
}
pp->normalize();
pp->draw(false,"decay comparison");

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Print out limits from a file
## ----------------------------------------------------------------------------------------------------------------------------------------

TFile *f = TFile::Open("higgsCombineTest.AsymptoticLimits.mH.root");
TTree *t = (TTree*)f->Get("limit");
double limit, mh;
float quantileExpected;
t->SetBranchAddress("limit",&limit);
t->SetBranchAddress("mh",&mh);
t->SetBranchAddress("quantileExpected",&quantileExpected);

std::vector<std::pair<double,double>> massLims;
for(int i=0; i<t->GetEntries(); ++i){
t->GetEntry(i);
if(quantileExpected != 0.5) continue;
massLims.push_back(std::make_pair(mh,limit));
}
std::sort(massLims.begin(),massLims.end(),[](std::pair<double,double> a, std::pair<double,double> b){return a.first < b.first;});

for(const auto& ml : massLims) {
std::cout<<"m"<<ml.first<<": "<<ml.second<<std::endl;
}

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## HEM dilepton plots
## ----------------------------------------------------------------------------------------------------------------------------------------

TFile *f = TFile::Open("debug2018.root");
std::vector<TString> sufs = {"eta","phi"};
for(const auto& suf : sufs){
Plotter *p = new Plotter();
TH1 *hd = (TH1*)f->Get("data_nosel2__e_"+suf);
TH1 *hd2 = (TH1*)f->Get("data_nosel2__e_noHEM_"+suf);
p->addHist(hd,"data w/ corr.",-1,1,4,20,1,true,true);
p->addHist(hd2,"data w/o corr.",-1,1,4,20,1,true,true);
p->rebin(2);
p->draw(false,suf);
}

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Bkg MC plots for AN
## ----------------------------------------------------------------------------------------------------------------------------------------
std::vector<TString> procs = {"qg","losttw","mw","mt"};
std::map<TString,TString> bmap { {"qg","q/g"}, {"mw","m_{W}"}, {"losttw","lost t/W"}, {"mt","m_{t}"} };
std::vector<TString> regs= {"emu_LMT_I_full","IF_LMT_full"};
std::vector<TString> vars= {"hbbMass","hhMass"};


for(const auto& r : regs) for(const auto& v : vars) {
std::cout<<r<<" "<<v<<std::endl;

double bw = 0;
Plotter *pp = new Plotter();

for(const auto& p : procs){
std::cout<<p<<std::endl;
TFile *f = TFile::Open("bkgInputs/HHbb1o2l_"+p+"_distributions.root");
TString sup = (p=="qg"&&r.Contains("emu")?"wQCD_noQCDSF_":"");
TH1 *h = (TH1*)f->Get(p+"_"+sup+r+"_"+v);

bw = h->GetBinWidth(1);
for(int i=2; i<=h->GetNbinsX(); ++i) {
double inCont = h->GetBinContent(i);
if(inCont <= 0) continue;
if(h->GetBinWidth(i) <= bw) continue;

double sf = bw / h->GetBinWidth(i);
h->SetBinContent(i,inCont*sf);

double inErr = h->GetBinError(i);
if(inErr == 0) continue;
h->SetBinError(i,inErr*sf);
}

pp->addStackHist(h,bmap[p]);
}

pp->setCMSLumi(0,10);
float rf = (v=="hhMass"?4:1);
pp->rebin(rf);
pp->setYTitle(TString::Format("Events / %d GeV",int(bw*rf)));

pp->draw(false,r+v);
}

## ------------------------------------------------------------------------------------------------------------------------------------

## ----------------------------------------------------------------------------------------------------------------------------------------
## Lepton efficiency plots in signal
## ----------------------------------------------------------------------------------------------------------------------------------------

std::vector<TString> chans{"1l","2l"};
std::vector<TString> flavs{"m","e"};
std::vector<TString> flavs2{"mm","ee","me","em"};

std::vector<TString> sels {"base","recoIDIso","recoID","reco"};
std::vector<TString> selSs {"den","reco. + ID + Iso","reco. + ID","reco."};

TString denID = "genpteta_";
TString sig = "radion";

TFile *f = TFile::Open("leptonEffs.root");
for(const auto& ch : chans) {
	for(const auto& fl : flavs) {
		Plotter *p = new Plotter();
		for(int i=0; i<sels.size();++i){
			TH1 *h = (TH1*)f->Get(sig+"_"+ch+"_"+denID+fl+"_"+sels[i]+"_mx");
			p->addHist(h,selSs[i],-1,1,4,20,1,true,false);
		}
		p->setYTitle("Efficiency");
		p->drawRatio(0,"stack",false,false,ch+" "+fl);
	}
	if(ch=="1l") continue;
	for(const auto& fl : flavs2) {
		Plotter *pp = new Plotter();
		for(int i=0; i<sels.size();++i){
			TH1 *h = (TH1*)f->Get(sig+"_"+ch+"_"+denID+fl+"_"+sels[i]+"_mx");
			pp->addHist(h,selSs[i],-1,1,4,20,1,true,false);
		}
		pp->setYTitle("Efficiency");
		pp->drawRatio(0,"stack",false,false,ch+" "+fl);
	}

}
## ------------------------------------------------------------------------------------------------------------------------------------