diff --git a/WMass/python/plotter/ccFiles/functions.cc b/WMass/python/plotter/ccFiles/functions.cc index a55b3f77757f..6a17c346272e 100644 --- a/WMass/python/plotter/ccFiles/functions.cc +++ b/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 > PtEtaPhiMVector; PtEtaPhiMVector p1(pt[0], eta[0], phi[0], m[0]);; diff --git a/WMass/python/plotter/ewBornMass.py b/WMass/python/plotter/ewBornMass.py new file mode 100755 index 000000000000..34e4bfa16161 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/WMass/python/plotter/ewBornMassRatio.py b/WMass/python/plotter/ewBornMassRatio.py new file mode 100755 index 000000000000..627619a1083f --- /dev/null +++ b/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) \ No newline at end of file diff --git a/WMass/python/plotter/ewPhotonKinematics.py b/WMass/python/plotter/ewPhotonKinematics.py index e44b84818e05..a1d0d6b24f7f 100755 --- a/WMass/python/plotter/ewPhotonKinematics.py +++ b/WMass/python/plotter/ewPhotonKinematics.py @@ -6,6 +6,9 @@ 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 @@ -23,17 +26,46 @@ 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=500, help="Maximum number of files") +parser.add_argument("-w", "--withISR", type=bool, default=True, help="Consider ISR photons") +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") -if len(args) > 1: - name = args[1] -else: - name = 'test' -maxFiles = 500 - -withISR = True - # ROOT.gStyle.SetOptStat(0) ROOT.gStyle.SetPalette(ROOT.kViridis) ROOT.TColor.InvertPalette() @@ -43,37 +75,56 @@ def makePlot(name): chain = ROOT.TChain('Events') - if name == 'test': - chain.Add('/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoAOD/DYJetsToMuMu_M-50_TuneCP5_13TeV-powhegMiNNLO-pythia8-photos/NanoV8MCPreVFP/210319_193015/0000/NanoV8MCPreVFP_weightFix_313.root') - - else: - path = paths[name] - - nFilesAdded = 0 - for f in os.listdir(path): - if nFilesAdded > maxFiles: - break - if '.root' in f: - chain.Add(path+f) + 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) - if withISR: + if args.withISR: name = name + "_withISR" rdf = ROOT.RDataFrame(chain) # rdf = ROOT.RDataFrame('Events', '/afs/cern.ch/work/m/mseidel/WMass/MC/CMSSW_10_6_19_patch2/src/Configuration/WMassNanoGen/NanoGen.root') 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])') - # FIXME: handle pair production - rdf = rdf.Define('ewSel', 'Numba::ewPhotonKinematicsSel(GenPart_status, GenPart_statusFlags, GenPart_pdgId, GenPart_genPartIdxMother, GenPart_pt, GenPart_eta, GenPart_phi, %s)' % str(withISR).lower()) + rdf = rdf.Define('rapVgen', 'rapidity(GenPart_pt[genPrefsrLeps],GenPart_phi[genPrefsrLeps],GenPart_phi[genPrefsrLeps],GenPart_mass[genPrefsrLeps])') + rdf = rdf.Define('massVgen', 'invariantmass(GenPart_pt[genPrefsrLeps],GenPart_phi[genPrefsrLeps],GenPart_phi[genPrefsrLeps],GenPart_mass[genPrefsrLeps])') + rdf = rdf.Define('ewSel', 'Numba::ewPhotonKinematicsSel(GenPart_status, GenPart_statusFlags, GenPart_pdgId, GenPart_genPartIdxMother, GenPart_pt, GenPart_eta, GenPart_phi, %s)' % str(args.withISR).lower()) rdf = rdf.Define('sij', 'log10(invMLepPhotons(GenPart_pt[ewSel==1],GenPart_eta[ewSel==1],GenPart_phi[ewSel==1],0.105658,GenPart_pt[ewSel==3],GenPart_eta[ewSel==3],GenPart_phi[ewSel==3]))') - rdf = rdf.Define('sik', 'log10(invMLepPhotons(GenPart_pt[ewSel==1],GenPart_eta[ewSel==1],GenPart_phi[ewSel==1],0.105658,GenPart_pt[ewSel==2],GenPart_eta[ewSel==2],GenPart_phi[ewSel==2]))') + rdf = rdf.Define('sik', '(invMLepPhotons(GenPart_pt[ewSel==1|ewSel==2],GenPart_eta[ewSel==1|ewSel==2],GenPart_phi[ewSel==1|ewSel==2],0.105658, GenPart_pt[ewSel==99],GenPart_eta[ewSel==99],GenPart_phi[ewSel==99]))') rdf = rdf.Define('sjk', 'log10(invMLepPhotons(GenPart_pt[ewSel==2],GenPart_eta[ewSel==2],GenPart_phi[ewSel==2],0.105658,GenPart_pt[ewSel==3],GenPart_eta[ewSel==3],GenPart_phi[ewSel==3]))') - # rdf = rdf.Define('sijk', '(invMLepPhotons(GenPart_pt[ewSel>=2],GenPart_eta[ewSel>=2],GenPart_phi[ewSel>=2],0.105658,GenPart_pt[ewSel==3],GenPart_eta[ewSel==3],GenPart_phi[ewSel==3]))') + rdf = rdf.Define('sijk', '(invMLepPhotons(GenPart_pt[ewSel==1|ewSel==2],GenPart_eta[ewSel==1|ewSel==2],GenPart_phi[ewSel==1|ewSel==2],0.105658, GenPart_pt[ewSel==3],GenPart_eta[ewSel==3],GenPart_phi[ewSel==3]))') + rdf = rdf.Define('logMassDiff', 'log10(sijk-sik+1e-5)') + rdf = rdf.Define('pt2ijk', 'pt2ijk(GenPart_pt[ewSel==1],GenPart_eta[ewSel==1],GenPart_phi[ewSel==1],0.105658, GenPart_pt[ewSel==3],GenPart_eta[ewSel==3],GenPart_phi[ewSel==3],0., GenPart_pt[ewSel==2],GenPart_eta[ewSel==2],GenPart_phi[ewSel==2],0.105658)') + + 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', 'sij', 'sik', 'sjk', 'sijk', 'pt2ijk']) - # rdf.Snapshot('Events', 'snapshot_%s.root' % name, ['ptVgen', 'ewSel', 'nGenPart', 'GenPart_pdgId', 'GenPart_status', 'GenPart_statusFlags', 'GenPart_genPartIdxMother', 'GenPart_pt', 'GenPart_eta', 'GenPart_phi', 'sij', 'sik', 'sjk']) + # hist = rdf.Histo2D((name, ';log_{10} s(#mu- #sum #gamma);log_{10} s(#mu+ #sum #gamma)', 100, -1.5, 3.5, 100, -1.5, 3.5), 'sij', 'sjk') - hist = rdf.Histo2D((name, ';log_{10} s(#mu- #sum #gamma);log_{10} s(#mu+ #sum #gamma)', 100, -1.5, 3.5, 100, -1.5, 3.5), 'sij', 'sjk') + if args.boson == 'Z': + massBins = makeBinning(mass = 91.1535, width = 2.4932, initialStep=0.010) + else: + massBins = makeBinning(mass = 80.3815, width = 2.0904, initialStep=0.010) + massDiffBins = np.linspace(-5, 5, 101) + hmodel = ROOT.RDF.TH2DModel(name, ';M(#mu+ #mu-) [GeV];log_{10} #{}{ M(#mu+ #mu- #sum #gamma) - M(#mu+ #mu-) + 10^{-5} } [GeV]', len(massBins)-1, array('d', massBins), len(massDiffBins)-1, array('d', massDiffBins)) + hist = rdf.Histo2D(hmodel, 'sik', 'logMassDiff') c = ROOT.TCanvas('c','c',500,500) c.cd() @@ -85,24 +136,48 @@ def makePlot(name): c.cd() hist.Draw('colz') - - c.Print('plots/ewPhotonKinematics/ewPK_%s.pdf' % name) - c.Print('plots/ewPhotonKinematics/ewPK_%s.png' % name) - c.Print('plots/ewPhotonKinematics/ewPK_%s.root' % name) - - -paths = { - 'minnlo': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoAOD/DYJetsToMuMu_M-50_TuneCP5_13TeV-powhegMiNNLO-pythia8-photos/NanoV8MCPreVFP/210319_193015/0000/', - 'horace-photos': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-born-fsr-photos-isr-pythia/', - 'horace-pythia': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-born-fsr-pythia-isr-pythia/', - 'horace-exp': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-exp-fsr-off-isr-off/', - 'horace-exp-old': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-exp-old-fsr-off-isr-pythia/', - 'horace-alpha-old': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-alpha-old-fsr-off-isr-pythia/', - 'horace-photoslow': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-born-fsr-photoslow-isr-pythia/', + c.Update() + + st = hist.FindObject('stats') + st.SetX1NDC(0.60) + st.SetX2NDC(0.87) + c.Modified() + c.Update() + + c.Print('plots/ewPhotonKinematics/%s/ewPK_%s.pdf' % (args.boson, name)) + c.Print('plots/ewPhotonKinematics/%s/ewPK_%s.png' % (args.boson, name)) + c.Print('plots/ewPhotonKinematics/%s/ewPK_%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/', + 'horace-pythia': baseNanoGen+'ZToMuMu_TuneCP5_13TeV-horace-born-fsr-pythia-isr-pythia/', + 'horace-exp': baseNanoGen+'ZToMuMu_TuneCP5_13TeV-horace-exp-fsr-off-isr-off/', + 'horace-exp-old': baseNanoGen+'ZToMuMu_TuneCP5_13TeV-horace-exp-old-fsr-off-isr-pythia/', + 'horace-alpha-old': baseNanoGen+'ZToMuMu_TuneCP5_13TeV-horace-alpha-old-fsr-off-isr-pythia/', + 'horace-photoslow': baseNanoGen+'ZToMuMu_TuneCP5_13TeV-horace-born-fsr-photoslow-isr-pythia/', + 'horace-photosnopair': baseNanoGen+'ZToMuMu_TuneCP5_13TeV-horace-born-fsr-photosnopair-isr-pythia/', } +paths['Wplus'] = { + 'minnlo': baseNano+'WplusJetsToMuNu_TuneCP5_13TeV-powhegMiNNLO-pythia8-photos/MC*VFPWeightFix/210701_234251/*/', + # 'horace-photosnopair': baseNanoGen+'WplusToMuNu_TuneCP5_13TeV-horace-born-fsr-photosnopair-isr-pythia', + 'horace-photos': baseNanoGen+'WplusToMuNu_TuneCP5_13TeV-horace-born-fsr-photos-isr-pythia', + 'horace-exp': baseNanoGen+'WplusToMuNu_TuneCP5_13TeV-horace-exp-fsr-off-isr-off', + 'horace-exp-old': baseNanoGen+'WplusToMuNu_TuneCP5_13TeV-horace-exp-old-fsr-off-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', + 'horace-exp': baseNanoGen+'WminusToMuNu_TuneCP5_13TeV-horace-exp-fsr-off-isr-off', + 'horace-exp-old': baseNanoGen+'WminusToMuNu_TuneCP5_13TeV-horace-exp-old-fsr-off-isr-pythia', +} + -if name == 'all': - for pathname in paths: +if args.input == 'all': + for pathname in paths[args.boson]: makePlot(pathname) else: - makePlot(name) \ No newline at end of file + makePlot(args.input) \ No newline at end of file diff --git a/WMass/python/plotter/ewPhotonKinematicsRatio.py b/WMass/python/plotter/ewPhotonKinematicsRatio.py index e94605ff4e9f..8acfe0656729 100755 --- a/WMass/python/plotter/ewPhotonKinematicsRatio.py +++ b/WMass/python/plotter/ewPhotonKinematicsRatio.py @@ -10,36 +10,46 @@ ## 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() -withISR = True +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(name1, name2): - print('making plot for %s/%s' % (name1, name2)) +def makeRatio(boson, name1, name2): + print('making plot for %s: %s/%s' % (boson, name1, name2)) - hist = hists[name1].Clone() - hist.Divide(hists[name2]) - # hist.Smooth(1) + hist = hists[name1].Clone('weights') + hist.Divide(hists[name1], hists[name2]) + ROOT.gStyle.SetPalette(ROOT.kThermometer) + ROOT.gStyle.SetNumberContours(100) + + weightHistOrig = ROOT.TH1D('weightHistOrig', ';weight;bins', 60, 0, 3) + weightHistSmooth = ROOT.TH1D('weightHistSmooth', ';weight;bins', 60, 0, 3) + + # smoothing for i in range(hist.GetNbinsX()+2): for j in range(hist.GetNbinsY()+2): val = hist.GetBinContent(i, j) if val == 0.: hist.SetBinContent(i, j, 1.) continue - errWeight = 2. + errWeight = 2. # approach 1 when error is 50% relErr = min(1., errWeight * hist.GetBinError(i, j) / val) ipolVal = (1.-relErr) * val + relErr hist.SetBinContent(i, j, ipolVal) + weightHistOrig.Fill(val) + weightHistSmooth.Fill(ipolVal) c = ROOT.TCanvas('c','c',500,500) c.cd() @@ -49,64 +59,90 @@ def makeRatio(name1, name2): c.SetTopMargin(0.05) c.cd() - ROOT.gStyle.SetPalette(ROOT.kRedBlue) - ROOT.gStyle.SetNumberContours(100) - ROOT.TColor.InvertPalette() - - hist.GetZaxis().SetRangeUser(0,3) + hist.GetZaxis().SetRangeUser(0., 2.) hist.SetContour(100) hist.Draw('colz') - c.Print('plots/ewPhotonKinematics/ewPKRatio_%s_%s.pdf' % (name1, name2)) - c.Print('plots/ewPhotonKinematics/ewPKRatio_%s_%s.png' % (name1, name2)) - c.Print('plots/ewPhotonKinematics/ewPKRatio_%s_%s.root' % (name1, name2)) + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s.pdf' % (boson, name1, name2)) + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s.png' % (boson, name1, name2)) + # c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s.root' % (args.boson, name1, name2)) + rootFile = ROOT.TFile('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s.root' % (boson, name1, name2), 'RECREATE') + hist.Write() + rootFile.Close() + + # Plot error histogram + ROOT.gStyle.SetPalette(ROOT.kDarkBodyRadiator) + hist_err = hist.Clone('weights_err') + for i in range(0, hist_err.GetNbinsX()+2): + for j in range(0, hist_err.GetNbinsY()+2): + hist_err.SetBinContent(i, j, hist.GetBinError(i,j)) + c.SetLogz() + hist_err.GetZaxis().SetRangeUser(1e-3, 10) + hist_err.Draw('colz') + + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s_err.pdf' % (boson, name1, name2)) + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s_err.png' % (boson, name1, name2)) + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s_err.root' % (boson, name1, name2)) + + # Plot weight histogram + c.SetRightMargin(0.05) + weightHistOrig.SetLineColor(ROOT.kRed+1) + weightHistOrig.SetLineStyle(7) + weightHistSmooth.SetLineColor(ROOT.kBlue+1) + + weightHistOrig.Draw() + weightHistSmooth.Draw('same') + + legend = ROOT.TLegend(0.5,0.6,0.95,0.9) + legend.SetLineWidth(0) + legend.SetFillStyle(0) + legend.AddEntry(weightHistOrig, 'Original weights', "l") + legend.AddEntry(weightHistSmooth, 'Smoothed weights', "l") + legend.Draw() + + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s_1d.png' % (boson, name1, name2)) + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s_1d.root' % (boson, name1, name2)) + c.Print('plots/ewPhotonKinematics/%s/ewPKRatio_%s_%s_1d.pdf' % (boson, name1, name2)) -def getHist(name): - tfile = ROOT.TFile('plots/ewPhotonKinematics/ewPK_%s.root' % name) +def getHist(boson, name): + tfile = ROOT.TFile('plots/ewPhotonKinematics/%s/ewPK_%s.root' % (boson, name)) canvas = tfile.Get('c') hist = canvas.GetPrimitive(name) # merge no/low-emission bins - maxBinLE = hist.GetXaxis().FindFixBin(0.) + maxBinLE = hist.GetYaxis().FindFixBin(-2.) # integralLE = hist.Integral(11, maxBinLE, 11, maxBinLE) - integralLE = 0. - integralLEerr = 0. - edgeRadius = 9 - for i in range(11, maxBinLE): - for j in range(11, maxBinLE): - dx = max(0., i-maxBinLE+edgeRadius) - dy = max(0., j-maxBinLE+edgeRadius) - if sqrt(dx*dx + dy*dy) < edgeRadius-0.5: - integralLE += hist.GetBinContent(i, j) - integralLEerr += hist.GetBinError(i, j) - for i in range(11, maxBinLE): - for j in range(11, maxBinLE): - dx = max(0., i-maxBinLE+edgeRadius) - dy = max(0., j-maxBinLE+edgeRadius) - if sqrt(dx*dx + dy*dy) < edgeRadius-0.5: - hist.SetBinContent(i, j, integralLE) - hist.SetBinError(i, j, integralLEerr) + for i in range(0, hist.GetNbinsX()+2): + integralLE = 0. + integralLEerr = 0. + for j in range(0, maxBinLE): + integralLE += hist.GetBinContent(i, j) + integralLEerr += hist.GetBinError(i, j) + for j in range(0, maxBinLE): + hist.SetBinContent(i, j, integralLE) + hist.SetBinError(i, j, integralLEerr) hist.Scale(1./hist.Integral()) return hist -paths = { - 'minnlo': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoAOD/DYJetsToMuMu_M-50_TuneCP5_13TeV-powhegMiNNLO-pythia8-photos/NanoV8MCPreVFP/210319_193015/0000/', - 'horace-photos': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-born-fsr-photos-isr-pythia/', - 'horace-pythia': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-born-fsr-pythia-isr-pythia/', - 'horace-exp': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-exp-fsr-off-isr-off/', - 'horace-exp-old': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-exp-old-fsr-off-isr-pythia/', - 'horace-alpha-old': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-alpha-old-fsr-off-isr-pythia/', - 'horace-photoslow': '/eos/cms/store/cmst3/group/wmass/w-mass-13TeV/NanoGen/ZToMuMu_TuneCP5_13TeV-horace-born-fsr-photoslow-isr-pythia/', -} - -name1 = args[1] -name2 = args[2] +if args.boson == 'all': + bosons = ['Z', 'Wplus', 'Wminus'] +else: + bosons = [args.boson] hists = {} -# for name in paths: -for name in [name1, name2]: - hists[name] = getHist(name + '_withISR') -makeRatio(name1, name2) \ No newline at end of file +for boson in bosons: + if args.input == 'all': + allHists = ['minnlo', 'horace-photos', 'horace-exp', 'horace-exp-old'] + for name in allHists: + hists[name] = getHist(boson, name + '_withISR') + 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('/')) + for name in [name1, name2]: + hists[name] = getHist(boson, name + '_withISR') + makeRatio(boson, name1, name2) \ No newline at end of file diff --git a/WMass/python/plotter/numbaDefines.py b/WMass/python/plotter/numbaDefines.py index 0f8ff02d63d3..328b73f4381b 100644 --- a/WMass/python/plotter/numbaDefines.py +++ b/WMass/python/plotter/numbaDefines.py @@ -144,13 +144,15 @@ def prefsrLeptons(status, statusFlags, pdgId, motherIdx, pts): def ewPhotonKinematicsSel(status, statusFlags, pdgId, motherIdx, pts, etas, phis, withISR = False): isLepton = (np.abs(pdgId) >= 11) & (np.abs(pdgId) <= 14) & (motherIdx >= 0) isMuon = isLepton & (np.abs(pdgId) == 13) + isMuonNu = isLepton & (np.abs(pdgId) == 14) isPhoton = pdgId == 22 status1 = status == 1 isPrompt = ((statusFlags >> 0 ) & 1).astype(np.bool_) isHardProcess = ((statusFlags >> 8 ) & 1).astype(np.bool_) - leptons = isMuon & status1 & isPrompt & isHardProcess + leptons = (isMuon | isMuonNu) & status1 & isPrompt & isHardProcess photons = isPhoton & status1 & isPrompt + # TODO: handle pair production? if not withISR: motherV = (pdgId[motherIdx] == 23) | (np.abs(pdgId[motherIdx]) == 24) leptons = leptons & motherV