| 1 |
|
#! /usr/bin/env python2 |
| 2 |
|
# -*- coding: utf-8 -*- |
| 3 |
|
from treeFunctions import * |
| 4 |
+ |
import Styles |
| 5 |
+ |
#Styles.tdrStyle() |
| 6 |
|
import numpy |
| 7 |
|
import ROOT |
| 8 |
|
import argparse |
| 35 |
|
v2.push_back(element) |
| 36 |
|
return v2 |
| 37 |
|
|
| 38 |
< |
def genMatching( photons, photonElectrons, genElectrons, hist ): |
| 38 |
> |
def recElectronMatching( genElectrons, photonElectrons ): |
| 39 |
> |
for recE in photonElectrons: |
| 40 |
> |
minDeltaR = 50 |
| 41 |
> |
minPt = 8000 |
| 42 |
> |
minIndex = -1 |
| 43 |
> |
for i in range( genElectrons.size()): |
| 44 |
> |
dr = deltaR( recE, genElectrons[i] ) |
| 45 |
> |
pt = ( recE.pt - genElectrons[i].pt ) / genElectrons[i].pt |
| 46 |
> |
if dr < minDeltaR: |
| 47 |
> |
minIndex = i |
| 48 |
> |
minDeltaR = dr |
| 49 |
> |
minPt = pt |
| 50 |
> |
|
| 51 |
> |
if minDeltaR < 0.02 and abs(minPt) < 0.5: |
| 52 |
> |
photonElectrons[minIndex].phi = 5 |
| 53 |
> |
genElectrons[minIndex].phi = 5 |
| 54 |
> |
|
| 55 |
> |
return genElectrons, photonElectrons |
| 56 |
> |
|
| 57 |
> |
def genElectronMatching( photons, photonElectrons, genElectrons, hist ): |
| 58 |
|
for genE in genElectrons: |
| 59 |
|
minDeltaR = 50 |
| 60 |
|
minPt = 8000 |
| 77 |
|
minPt = pt |
| 78 |
|
typeName = "e" |
| 79 |
|
|
| 80 |
< |
if minIndex != -1: |
| 80 |
> |
#if minDeltaR < 0.02 and abs(minPt) < 0.5: |
| 81 |
> |
if minDeltaR<50: |
| 82 |
|
hist.Fill( minPt, minDeltaR ) |
| 83 |
|
if typeName == "e": |
| 84 |
|
photonElectrons[minIndex].pixelseed = -2 |
| 86 |
|
photons[minIndex].pixelseed = -2 |
| 87 |
|
return photons, photonElectrons, hist |
| 88 |
|
|
| 89 |
< |
def process( inputFileName, nExpected ): |
| 89 |
> |
def weightTree( photonTree, photonJetTree ): |
| 90 |
> |
# a tree can't be modified, therefore a new tree is saved within the same file |
| 91 |
> |
|
| 92 |
> |
# this is the binning with witch the fo reweighting will be done and has to be studied |
| 93 |
> |
nBins = 10 |
| 94 |
> |
xMin = 80 |
| 95 |
> |
xMax = 800 |
| 96 |
> |
h_photonPt = createHistoFromTree( photonTree, "photon[0].pt", "", nBins, xMin, xMax ) |
| 97 |
> |
h_jetPt = createHistoFromTree( photonJetTree, "photon[0].pt", "", nBins, xMin, xMax ) |
| 98 |
> |
|
| 99 |
> |
# h_jetPt is now histogram with w^{-1} = h_jetPt / h_photonPt |
| 100 |
> |
h_jetPt.Divide( h_photonPt ) |
| 101 |
> |
h_jetPt.SetTitle(";p_{T,#gamma};w^{-1}") |
| 102 |
> |
#h_jetPt.Draw() |
| 103 |
> |
#raw_input() |
| 104 |
> |
|
| 105 |
> |
weightTree = ROOT.TTree("QCDWeightTree", "my TFriend for weights") |
| 106 |
> |
import numpy |
| 107 |
> |
|
| 108 |
> |
# a python float corresponds to a root double |
| 109 |
> |
weight = numpy.zeros(1, dtype=float) |
| 110 |
> |
weightError = numpy.zeros(1, dtype=float) |
| 111 |
> |
weightTree.Branch( "weight", weight, "weight/D" ) |
| 112 |
> |
weightTree.Branch( "weightError", weightError, "weightError/D" ) |
| 113 |
> |
|
| 114 |
> |
for event in photonJetTree: |
| 115 |
> |
try: |
| 116 |
> |
pt = event.photon.at(0).pt |
| 117 |
> |
except: |
| 118 |
> |
pt = 0 |
| 119 |
> |
bin = h_jetPt.FindBin( pt ) |
| 120 |
> |
weight[0] = h_jetPt.GetBinContent( bin ) |
| 121 |
> |
weightError[0] = h_jetPt.GetBinError( bin ) |
| 122 |
> |
weightTree.Fill() |
| 123 |
> |
|
| 124 |
> |
return weightTree |
| 125 |
> |
|
| 126 |
> |
def splitCandidates( inputFileName, nExpected, isQCD, isEWK): |
| 127 |
|
outputFileName = "slim"+inputFileName |
| 128 |
|
|
| 129 |
|
eventHisto = readHisto( inputFileName ) |
| 134 |
|
photonTree, photons = gammaSelectionClone( tree, "photonTree" ) |
| 135 |
|
photonJetTree, photonJets = gammaSelectionClone( tree, "photonJetTree" ) |
| 136 |
|
photonElectronTree, photonElectrons = gammaSelectionClone( tree, "photonElectronTree" ) |
| 137 |
+ |
#genElectronTree = gammaSelectionClone( tree, "genElectronTree" )[0] |
| 138 |
|
|
| 139 |
< |
pu_weight = numpy.zeros(1, dtype=float) |
| 140 |
< |
ROOT.gROOT.ProcessLine("float pu_weight;") |
| 141 |
< |
from ROOT import pu_weight |
| 139 |
> |
weight = numpy.zeros(1, dtype=float) |
| 140 |
> |
photonTree.SetBranchAddress("weight", weight) |
| 141 |
> |
photonElectronTree.SetBranchAddress("weight", weight) |
| 142 |
> |
photonJetTree.SetBranchAddress("weight", weight) |
| 143 |
> |
#genElectronTree.SetBranchAddress("weight", weight) |
| 144 |
|
|
| 145 |
< |
print type( pu_weight ) |
| 146 |
< |
photonTree.SetBranchAddress("pu_weight", pu_weight) |
| 145 |
> |
#genElectrons = ROOT.std.vector("tree::Particle")() |
| 146 |
> |
#photonElectronTree.SetBranchAddress("genElectron", genElectrons ) |
| 147 |
|
|
| 148 |
< |
hist2 = ROOT.TH2F("blub2", "", 100, -.2, 1, 100, 0, 0.06 ) |
| 148 |
> |
|
| 149 |
> |
hist2 = ROOT.TH2F( inputFileName, "", 100, -.30, .30, 100, 0, 0.06 ) |
| 150 |
|
hist2.SetTitle(";#frac{p_{T,rec}-p_{T,gen}}{p_{T,gen}};min #DeltaR") |
| 151 |
|
|
| 152 |
|
for event in tree: |
| 153 |
< |
pu_weight[0] = event.pu_weight * nExpected / nGenerated |
| 153 |
> |
if not event.GetReadEntry()%100000: |
| 154 |
> |
print 100*event.GetReadEntry()/event.GetEntries(), "%" |
| 155 |
> |
weight[0] = event.weight * nExpected / nGenerated |
| 156 |
|
photons.clear() |
| 157 |
|
photonElectrons.clear() |
| 158 |
|
photonJets.clear() |
| 177 |
|
photons.push_back(gamma) |
| 178 |
|
|
| 179 |
|
# QCD fake object definition |
| 180 |
< |
if gamma.pt > 80 \ |
| 180 |
> |
if gamma.ptJet > 75 \ |
| 181 |
|
and gamma.hadTowOverEm < 0.05 \ |
| 182 |
|
and gamma.sigmaIetaIeta < 0.014 \ |
| 183 |
|
and gamma.chargedIso < 15 \ |
| 187 |
|
and ( gamma.sigmaIetaIeta>=0.012 or gamma.chargedIso>=2.6): |
| 188 |
|
photonJets.push_back( gamma ) |
| 189 |
|
|
| 190 |
< |
photons, photonElectrons, hist2 = genMatching( photons, photonElectrons, event.genElectron, hist2 ) |
| 190 |
> |
photons, photonElectrons, hist2 = genElectronMatching( photons, photonElectrons, event.genElectron, hist2 ) |
| 191 |
|
|
| 192 |
+ |
#genElectrons = recElectronMatching( event.genElectron, photonElectrons ) |
| 193 |
|
|
| 194 |
|
if photons.size() > 0: |
| 195 |
|
photonTree.Fill() |
| 198 |
|
photonElectronTree.Fill() |
| 199 |
|
if photonJets.size() > 0: |
| 200 |
|
photonJetTree.Fill() |
| 201 |
+ |
#if (photons.size()>0 or photonElectrons.size()>0) and event.genElectron.size() > 0: |
| 202 |
+ |
# genElectronTree.Fill() |
| 203 |
|
|
| 204 |
< |
#hist2.Draw("colz") |
| 205 |
< |
#raw_input() |
| 204 |
> |
can = ROOT.TCanvas() |
| 205 |
> |
can.cd() |
| 206 |
> |
hist2.Draw("colz") |
| 207 |
> |
can.SaveAs("pt_r_%s_new.pdf"%inName[0:5]) |
| 208 |
> |
|
| 209 |
> |
if isQCD: |
| 210 |
> |
qcdWeights = weightTree( photonTree, photonJetTree ) |
| 211 |
|
|
| 139 |
– |
eventHisto = readHisto( inputFileName ) |
| 212 |
|
|
| 213 |
|
f = ROOT.TFile( outputFileName, "recreate" ) |
| 214 |
|
f.cd() |
| 215 |
|
photonTree.Write() |
| 216 |
|
photonJetTree.Write() |
| 217 |
|
photonElectronTree.Write() |
| 218 |
+ |
#genElectronTree.Write() |
| 219 |
+ |
if isQCD: |
| 220 |
+ |
qcdWeights.Write() |
| 221 |
|
f.Write() |
| 222 |
|
f.Close() |
| 223 |
|
|
| 230 |
|
arguments.add_argument("--input", default=["WJets_V01.03_tree.root"], nargs="+" ) |
| 231 |
|
opts = arguments.parse_args() |
| 232 |
|
|
| 233 |
+ |
ROOT.gROOT.SetBatch() |
| 234 |
|
|
| 235 |
|
|
| 236 |
|
datasetConf = ConfigParser.SafeConfigParser() |
| 239 |
|
integratedLumi = 19300 #pb |
| 240 |
|
|
| 241 |
|
for inName in opts.input: |
| 242 |
+ |
shortName = None |
| 243 |
|
for configName in datasetConf.sections(): |
| 244 |
|
if inName.count( configName ): |
| 245 |
+ |
shortName = configName |
| 246 |
|
crosssection = datasetConf.getfloat( configName, "crosssection" ) |
| 247 |
|
|
| 248 |
+ |
qcd = False |
| 249 |
+ |
ewk = False |
| 250 |
+ |
if shortName in ["QCD_250_500", "QCD_500_1000", "QCD_1000_inf","GJets"]: |
| 251 |
+ |
qcd = True |
| 252 |
+ |
if shortName in ["TTJets", "WJets"]: |
| 253 |
+ |
ewk = True |
| 254 |
+ |
|
| 255 |
|
# N = L * sigma |
| 256 |
|
nExpected = integratedLumi * crosssection |
| 257 |
|
|
| 258 |
< |
process( inName, nExpected ) |
| 258 |
> |
splitCandidates( inName, nExpected, isQCD=qcd, isEWK=ewk ) |
| 259 |
|
|
