EW stuff (untouched for >half a year)

WMass/python/plotter/ccFiles/functions.cc

@@ -174,6 +174,26 @@ float invMLepPhotons(const Vec_f& pt1, const Vec_f& eta1, const Vec_f& phi1, flo
   return psum.mass();
 }
 
+float pt2ijk(const Vec_f& pt_i, const Vec_f& eta_i, const Vec_f& phi_i, float m_i, const Vec_f& pt_j, const Vec_f& eta_j, const Vec_f& phi_j, float m_j, const Vec_f& pt_k, const Vec_f& eta_k, const Vec_f& phi_k, float m_k) {
+  // PtEtaPhiMVector p1(pt1[0], eta1[0], phi1[0], m1);
+  PtEtaPhiMVector p_i(0,0,0,0);
+  for (unsigned int i = 0; i < pt_i.size(); ++i) {
+    PtEtaPhiMVector ptemp(pt_i[i], eta_i[i], phi_i[i], m_i);
+    p_i += ptemp;
+  }
+  PtEtaPhiMVector p_j(0,0,0,0);
+  for (unsigned int i = 0; i < pt_j.size(); ++i) {
+    PtEtaPhiMVector ptemp(pt_j[i], eta_j[i], phi_j[i], m_j);
+    p_j += ptemp;
+  }
+  PtEtaPhiMVector p_k(0,0,0,0);
+  for (unsigned int i = 0; i < pt_k.size(); ++i) {
+    PtEtaPhiMVector ptemp(pt_k[i], eta_k[i], phi_k[i], m_k);
+    p_k += ptemp;
+  }
+  return 2*(p_i.Dot(p_j))*(p_j.Dot(p_k))/(p_i.Dot(p_j) + p_j.Dot(p_k) + p_i.Dot(p_k));
+}
+
 float rapidity(const Vec_f& pt, const Vec_f& eta, const Vec_f& phi, const Vec_f& m) {
   //typedef ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > PtEtaPhiMVector;
   PtEtaPhiMVector p1(pt[0], eta[0], phi[0], m[0]);;

WMass/python/plotter/ewBornMass.py

@@ -0,0 +1,179 @@
+#!/usr/bin/env python
+from math import *
+import re
+import os, os.path
+from array import array
+import logging
+from collections import OrderedDict
+import numbaDefines
+import numpy as np
+from glob import glob
+from tqdm import tqdm
+
+## safe batch mode
+import sys
+args = sys.argv[:]
+sys.argv = ['-b']
+import ROOT
+sys.argv = args
+ROOT.gROOT.SetBatch(True)
+ROOT.PyConfig.IgnoreCommandLineOptions = True
+ROOT.gInterpreter.ProcessLine(".O3")
+ROOT.EnableImplicitMT()
+
+from copy import *
+from cutsFile import *
+from fakeRate import *
+from mcCorrections import *
+
+import argparse
+parser = argparse.ArgumentParser()
+parser.add_argument("-i", "--input",  type=str, default='minnlo', help="File to read histos from")
+parser.add_argument("-b", "--boson", type=str, default='Z', help="boson")
+parser.add_argument("-m", "--maxFiles", type=int, default=50, help="Maximum number of files")
+parser.add_argument('--snapshot', action='store_true')
+parser.add_argument('--binning', action='store_true')
+parser.add_argument('-z', '--deleteZombies', action='store_true')
+args = parser.parse_args()
+
+def makeBinning(mass = 91.1535, width = 2.4932, initialStep = 0.1):
+    maxVal = ROOT.Math.breitwigner_pdf(mass, width, mass)
+    bins = [mass]
+    currentMass = mass
+    while currentMass - mass < 100:
+        binSize = maxVal / ROOT.Math.breitwigner_pdf(currentMass, width, mass) * initialStep
+        currentMass += binSize
+        bins.append(currentMass)
+        lowMass = 2*mass - currentMass
+        if lowMass - binSize > 0:
+            bins.insert(0, lowMass)
+    bins.insert(0, 0.)
+    return bins
+
+def normalizeAndDivideByBinWidth(hist):
+    hist.Scale(1./hist.Integral())
+    for i in range(1, hist.GetNbinsX()+1):
+        for j in range(1, hist.GetNbinsY()+1):
+            hist.SetBinContent(i, j, hist.GetBinContent(i,j) / hist.GetXaxis().GetBinWidth(i))
+            hist.SetBinError  (i, j, hist.GetBinError(i,j)   / hist.GetXaxis().GetBinWidth(i))
+    return hist
+
+if args.binning:
+    print(makeBinning())
+    sys.exit()
+
+if "/functions_cc.so" not in ROOT.gSystem.GetLibraries(): 
+    compileMacro("ccFiles/functions.cc")
+
+# ROOT.gStyle.SetOptStat(0)
+ROOT.gStyle.SetPalette(ROOT.kViridis)
+ROOT.TColor.InvertPalette()
+
+def makePlot(name):
+    print('making plot for ' + name)
+
+    chain = ROOT.TChain('Events')
+
+    path = paths[args.boson][name]
+    files = glob(path+'/*.root')
+    nFilesAdded = 0
+    print('Adding files')
+    for f in tqdm(files):
+        if args.maxFiles > 0 and nFilesAdded > args.maxFiles:
+            break
+        if '.root' in f:
+            markedZombie = False
+            if args.deleteZombies and not 'powhegMiNNLO' in path:
+                try:
+                    thisFile = ROOT.TFile(f)
+                    thisFile.Close()
+                except OSError:
+                    os.remove(f)
+                    markedZombie = True
+            if not markedZombie:
+                chain.Add(f)
+                nFilesAdded += 1
+    print('Added %i files' % nFilesAdded)
+
+    rdf = ROOT.RDataFrame(chain)
+    rdf = rdf.Define('genPrefsrLeps', 'Numba::prefsrLeptons(GenPart_status, GenPart_statusFlags, GenPart_pdgId, GenPart_genPartIdxMother, GenPart_pt)')
+    rdf = rdf.Define('ptVgen', 'transversemomentum(GenPart_pt[genPrefsrLeps],GenPart_phi[genPrefsrLeps])')
+    rdf = rdf.Define('rapVgen', 'abs(rapidity(GenPart_pt[genPrefsrLeps],GenPart_eta[genPrefsrLeps],GenPart_phi[genPrefsrLeps],GenPart_mass[genPrefsrLeps]))')
+    rdf = rdf.Define('massVgen', 'invariantmass(GenPart_pt[genPrefsrLeps],GenPart_eta[genPrefsrLeps],GenPart_phi[genPrefsrLeps],GenPart_mass[genPrefsrLeps])')
+
+    if args.snapshot:
+        rdf.Snapshot('Events', 'snapshot_%s.root' % name, ['ptVgen', 'ewSel', 'nGenPart', 'GenPart_pdgId', 'GenPart_status', 'GenPart_statusFlags', 'GenPart_genPartIdxMother', 'GenPart_pt', 'GenPart_eta', 'GenPart_phi', 'ptVgen', 'rapVgen', 'massVgen'])
+
+    if args.boson == 'Z':
+        massBins = makeBinning(mass = 91.1535, width = 2.4932, initialStep=0.10)
+    else:
+        massBins = makeBinning(mass = 80.3815, width = 2.0904, initialStep=0.10)
+    ptBins = np.logspace(-1,3,9)
+    ptModel = ROOT.RDF.TH2DModel(name+'_pt', ';Boson mass [GeV];Boson p_{T} [GeV]', len(massBins)-1, array('d', massBins), len(ptBins)-1, array('d', ptBins))
+    ptHist = rdf.Histo2D(ptModel, 'massVgen', 'ptVgen')
+
+    rapBins = np.linspace(0, 5, 11)
+    rapModel = ROOT.RDF.TH2DModel(name+'_rap', ';Boson mass [GeV];Boson rapidity', len(massBins)-1, array('d', massBins), len(rapBins)-1, array('d', rapBins))
+    rapHist = rdf.Histo2D(rapModel, 'massVgen', 'rapVgen')
+
+    c = ROOT.TCanvas('c','c',500,500)
+    c.cd()
+    c.SetRightMargin(0.12)
+    c.SetLeftMargin(0.15)
+    c.SetTopMargin(0.1)
+    c.SetTopMargin(0.05)
+    c.SetLogz()
+    c.cd()
+    if not os.path.isdir('plots/ewBornMass/%s/' % args.boson):
+        os.makedirs('plots/ewBornMass/%s/' % args.boson)
+
+    c.SetLogy(True)
+    ptHist.Draw('colz')
+    c.Update()
+
+    st = ptHist.FindObject('stats')
+    st.SetX1NDC(0.60)
+    st.SetX2NDC(0.87)
+    c.Modified()
+    c.Update()
+
+    c.Print('plots/ewBornMass/%s/mass_pt_%s.pdf'  % (args.boson, name))
+    c.Print('plots/ewBornMass/%s/mass_pt_%s.png'  % (args.boson, name))
+    c.Print('plots/ewBornMass/%s/mass_pt_%s.root' % (args.boson, name))
+
+    c.SetLogy(False)
+    rapHist.Draw('colz')
+    c.Update()
+
+    st = rapHist.FindObject('stats')
+    st.SetX1NDC(0.60)
+    st.SetX2NDC(0.87)
+    c.Modified()
+    c.Update()
+
+    c.Print('plots/ewBornMass/%s/mass_rap_%s.pdf'  % (args.boson, name))
+    c.Print('plots/ewBornMass/%s/mass_rap_%s.png'  % (args.boson, name))
+    c.Print('plots/ewBornMass/%s/mass_rap_%s.root' % (args.boson, name))
+
+baseNano = '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoAOD/'
+baseNanoGen = '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/'
+paths = {}
+paths['Z'] = {
+    'minnlo': baseNano+'DYJetsToMuMu_M-50_TuneCP5_13TeV-powhegMiNNLO-pythia8-photos/NanoV8MC*VFP/210319_*/*/',
+    'horace-photos': baseNanoGen+'ZToMuMu_TuneCP5_13TeV-horace-born-fsr-photos-isr-pythia/',
+}
+paths['Wplus'] = {
+    'minnlo': baseNano+'WplusJetsToMuNu_TuneCP5_13TeV-powhegMiNNLO-pythia8-photos/MC*VFPWeightFix/210701_234251/*/',
+    'horace-photos': baseNanoGen+'WplusToMuNu_TuneCP5_13TeV-horace-born-fsr-photos-isr-pythia',
+}
+paths['Wminus'] = {
+    'minnlo': baseNano+'WminusJetsToMuNu_TuneCP5_13TeV-powhegMiNNLO-pythia8-photos/MC*VFPWeightFix/210701_234131/*/',
+    'horace-photos': baseNanoGen+'WminusToMuNu_TuneCP5_13TeV-horace-born-fsr-photos-isr-pythia',
+}
+
+
+if args.input == 'all':
+    for pathname in paths[args.boson]:
+        makePlot(pathname)
+else:
+    makePlot(args.input)

WMass/python/plotter/ewBornMassRatio.py

@@ -0,0 +1,146 @@
+#!/usr/bin/env python
+from math import *
+import re
+import os, os.path
+from array import array
+import logging
+from collections import OrderedDict
+import numbaDefines
+from math import sqrt
+
+## safe batch mode
+import sys
+import ROOT
+ROOT.gROOT.SetBatch(True)
+ROOT.PyConfig.IgnoreCommandLineOptions = True
+ROOT.gInterpreter.ProcessLine(".O3")
+ROOT.EnableImplicitMT()
+
+import argparse
+parser = argparse.ArgumentParser()
+parser.add_argument("-i", "--input",  type=str, default='horace-photos/minnlo', help="File to read histos from")
+parser.add_argument("-b", "--boson", type=str, default='Z', help="boson")
+parser.add_argument("-w", "--withISR", type=bool, default=True, help="Consider ISR photons")
+args = parser.parse_args()
+
+ROOT.gStyle.SetOptStat(0)
+
+def makeRatio(boson, name1, name2, obs='rap'):
+    print('making plot for %s: %s/%s' % (boson, name1, name2))
+
+    legend = ROOT.TLegend(0.5,0.6,0.95,0.9)
+    legend.SetLineWidth(0)
+    legend.SetFillStyle(0)
+
+    tfile1 = ROOT.TFile('plots/ewBornMass/%s/mass_rap_%s.root' % (boson, name1))
+    canvas1 = tfile1.Get('c')
+    hist1 = canvas1.GetPrimitive(name1+'_'+obs)
+
+    hist1.Scale(1./hist1.Integral())
+
+    colors = [0, ROOT.kRed+1, ROOT.kOrange+1, ROOT.kGreen+1, ROOT.kCyan+1, ROOT.kAzure+1,  ROOT.kBlue+1, ROOT.kViolet+1, ROOT.kMagenta+1]
+
+    projs1 = []
+    ybins = [1, 3, 5, 7]
+    for ybin in ybins:
+        low = hist1.GetYaxis().GetBinLowEdge(ybin)
+        high = hist1.GetYaxis().GetBinUpEdge(ybin+1)
+        if low >= 4.:
+            break
+        proj = hist1.ProjectionX(name1+'_'+obs+str(ybin), ybin, ybin+1)
+        proj.SetLineColor(colors[ybin])
+        proj.GetYaxis().SetRangeUser(0.5, 1.3)
+        proj.GetYaxis().SetTitle('%s / %s' % (name1, name2))
+        proj.Scale(1./proj.Integral())
+        projs1.append(proj)
+        legend.AddEntry(proj, '%.2f < |Y| < %.2f' % (low, high), "l")
+
+    tfile2 = ROOT.TFile('plots/ewBornMass/%s/mass_rap_%s.root' % (boson, name2))
+    canvas2 = tfile2.Get('c')
+    hist2 = canvas2.GetPrimitive(name2+'_'+obs)
+
+    hist2.Scale(1./hist2.Integral())
+
+    projs2 = []
+    for ybin in ybins:
+        proj = hist2.ProjectionX(name2+'_'+obs+str(ybin), ybin, ybin)
+        proj.SetLineColor(1+ybin)
+        proj.Scale(1./proj.Integral())
+        projs2.append(proj)
+
+    print(projs1)
+
+    hist = hist1.Clone('ratio')
+    hist.Divide(hist1, hist2)
+
+    ROOT.gStyle.SetPalette(ROOT.kThermometer)
+    ROOT.gStyle.SetNumberContours(100)
+
+    c = ROOT.TCanvas('c','c',500,500)
+    c.cd()
+    c.SetRightMargin(0.12)
+    c.SetLeftMargin(0.15)
+    c.SetTopMargin(0.1)
+    c.SetTopMargin(0.05)
+    c.cd()
+
+    hist.GetZaxis().SetRangeUser(0.5, 1.5)
+    hist.SetContour(100)
+    hist.Draw('colz')
+
+    c.Print('plots/ewBornMass/%s/mass_rap_ratio_%s_%s.pdf'  % (boson, name1, name2))
+    c.Print('plots/ewBornMass/%s/mass_rap_ratio_%s_%s.png'  % (boson, name1, name2))
+    # c.Print('plots/ewBornMass/%s/mass_rap_ratio__%s_%s.root' % (args.boson, name1, name2))
+    rootFile = ROOT.TFile('plots/ewBornMass/%s/mass_rap_ratio_%s_%s.root' % (boson, name1, name2), 'RECREATE')
+    hist.Write()
+    rootFile.Close()
+
+    c.SetRightMargin(0.05)
+
+    for i in range(len(projs1)):
+        projs1[i].Divide(projs2[i])
+        projs1[i].GetXaxis().SetRangeUser(87, 95)
+        projs1[i].GetYaxis().SetRangeUser(0.95, 1.10)
+        if i == 0:
+            projs1[i].Draw()
+        else:
+            projs1[i].Draw('same')
+    legend.Draw()
+
+    c.Print('plots/ewBornMass/%s/mass_rap_ratio1d_%s_%s.pdf'  % (boson, name1, name2))
+    c.Print('plots/ewBornMass/%s/mass_rap_ratio1d_%s_%s.png'  % (boson, name1, name2))
+
+    legend = ROOT.TLegend(0.5,0.6,0.95,0.9)
+    legend.SetLineWidth(0)
+    legend.SetFillStyle(0)
+
+    hist1rap = hist1.ProjectionY(name1+'_rap')
+    hist2rap = hist2.ProjectionY(name2+'_rap')
+    hist1rap.SetLineColor(ROOT.kRed+1)
+    hist2rap.SetLineColor(ROOT.kBlue+1)
+    legend.AddEntry(hist1rap, name1, "l")
+    legend.AddEntry(hist2rap, name2, "l")
+
+    hist1rap.Draw()
+    hist2rap.Draw('same')
+    legend.Draw()
+
+    c.Print('plots/ewBornMass/%s/rap_%s_%s.pdf'  % (boson, name1, name2))
+    c.Print('plots/ewBornMass/%s/rap_%s_%s.png'  % (boson, name1, name2))
+
+
+if args.boson == 'all':
+    bosons = ['Z', 'Wplus', 'Wminus']
+else:
+    bosons = [args.boson]
+
+hists = {}
+
+for boson in bosons:
+    if args.input == 'all':
+        makeRatio(boson, 'minnlo', 'horace-photos')
+        makeRatio(boson, 'horace-exp', 'horace-photos')
+        makeRatio(boson, 'horace-exp-old', 'horace-photos')
+    else:
+        name1,name2 = tuple(args.input.split('/'))
+        makeRatio(boson, name1, name2)