From b7aec0ec04aad41cbe6496eef29d69348ce99883 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Brochet?= Date: Tue, 17 Nov 2015 10:54:28 +0100 Subject: [PATCH] Add new weight for PU reweighting --- data/PUProfiles/pu_histogram_all_2015.root | Bin 0 -> 5261 bytes interface/EventProducer.h | 7 ++ interface/PUReweighter.h | 86 +++++++++++++++++++++ interface/Tools.h | 6 ++ python/EventProducer.py | 6 +- src/EventProducer.cc | 5 ++ src/PUReweighter.cc | 53 +++++++++++++ 7 files changed, 162 insertions(+), 1 deletion(-) create mode 100644 data/PUProfiles/pu_histogram_all_2015.root create mode 100644 interface/PUReweighter.h create mode 100644 src/PUReweighter.cc diff --git a/data/PUProfiles/pu_histogram_all_2015.root b/data/PUProfiles/pu_histogram_all_2015.root new file mode 100644 index 0000000000000000000000000000000000000000..b6fdda05a9f8d47b59cd406e85d44cd158a7c8d6 GIT binary patch literal 5261 zcmb7obx_pdxBk-I-CY7Ixim|IlpwKyEFif83oI>NOA9IuyGV;PC?VY?OG$`y2uMpv z2;SA-Z|44YzxRG-&b(*l{k+dPGoSgK=XoEvhldXUuqO@x02}}S+GI>)i&qNdI z=E$rcsBY(D_ZXA|5QhNZLI?=J`j{T1V7kjZxN`9NpwUr6 z??jfv&PBS2pQir7UI>H68(q~F^uxm-KYqW2rwhFIwu}p(x_%i*?XDYk*KVwm)r7|z z48=z~45sC5Oj&%}mO0fzDx*2JE%MJp^Us4&=L$#Lw{1t3l&QB^MA#P0sZjtoeBO0# zJ!Zg91P@nUcpT3kd@4Si*J-i5QO7{-pT%WFhFIDyT@f-0`8<7lY zUTY#evAFhdm;l-pt}4v!E{`M1*Ux1)^sRD=!}|odCD+Nq)*%wj{!2Hh#~WjXcXLjz zr)g<}kdS!?>j(}Hzx6*;N#I5yhoBU&FkpfPrLq|^D-LNief(6m>N#rh3_ zoSWL!LKPFl)%)r8LT|p>>Zim!AUUyhhq(H=#<`m1lHI!?f@J8ECy{r{ToR&V)oMj& zpQ{y_%v{>q`PyHPgWjY>mn-)(yglh-y}Hd=bAiKxZECGaUZCeg=YAxJy`ZY7DK!_= z>>Moa=lda8m|o|fSSF(Awbjcyqc*Pq0}Gg0guF_y>ZExUVo+GoOt?6PSlWk8JO!hOoZ?sa{LIW>nNlDHNhr4mHDKHtS0X$-Cz=t}1PAT#w;RLLFz z84E4yEc9bi3@SE3!7H%Hyc9Gde&T4;+N8^Ti!101l^44vUQ}zzxYYe4uM*c3%k)Np zCMpPHpNEcTYCA zVu(H-LC5U79f~N3QhV(+&~gE5L&0DnvXDC&=g_<`mIYzxLEJSzxn*tS`g8_y0py_Z zE+!;jg>KsBYAYvZUcN=i^ijRsMcO7)h}buZ8@h&+laY1%nLuL1h3b?2=h@{gE-_L4 zR{KisD;yDcn!8qV?@R3chOy8`HmHG)&nq3>>wA3rZ>sW_3)8Y|muK`%zlbo_BBb8H zfckRvmV#+zL``sux%a)BFM>$GCxht{qYYZ%4%Ihr8&j$1M-r|l2fPJzovW7~{OHm7 zek12j3lg;>H;iUI^(>Lt^X_q3{Kk@PVe&lX;Y&1W%8|FiFe2fV1`%2J--1WvRvtcWY%j z4lY~}-O%0D`8DU3vKH~I2(fof4Su(y{uMThin#H}5uuO3f%=crGC!$dAD*#dT}`3M zWS0iduOaa?;XciDk!MLE+J*|S*&YVY1;eBVca2c2gXX=O>)d%tv1$~OzXaQX#^Lv+ zN*z1g4k)aAR;3b(9`zdxsiOL&e){sWhGW_>#i`Sv8An>5jUG=n{9b5t^^d2jHc2%!QONWWJ-+-tUYc3Pxt2gco@rGU$rS|gMTAHqzGA);Wm-Lm)7Ag zC{q^YWnDJN;Srg&f{c5$eGc$OiC4VVCCD<}jv?U0vMmiQyebt6-@W7Jt7pc^`Z07m zh0opv%8Y3R0z~RmB+s2Ad&OpPXmF2AAz}$!;2S{|&NZ@r_ZoR3Cst#JChDZIyF!hg z)!`o7XF&~5A4y#35S9rM@Jhm0lXWh2+7X+XL)=i{9;E|1&MAQULhgzZpf_1}Aq9aP zzB=yw+6k~*#e)9R0T|-sas5Y}pufcVR~-Nx@J#eD+)et#2kvO+?g-a`IeAd4d-y(c z3*>Zz!dx94FpIYjhD&Wp0Iob3E}gMl@?f2zQuK1NW34pO{N zRXR3~zr%9hJ<1F6S7_LG@+YhM;ySjN{Z;;#f5&0VcAi(nVYx919_@5YyO2V>KE9~*x|782#y1KGfHrIKm^8Vc)_i znkpY^_%!6mc{esHkwn<}*tQw8DH#WQ+j%cQqr3G%GVm!aeFp3mnv}JH`im}Vi5|MK zDL`SjUM%$rZVPpLui}Gu8#m-3jur7cOK8d=^8I9r*qOkN@^iWccLm>w9Vt@VdSJ|c zPU>xdIL`nn-ME@_diMJJu_-T}X(H8W@sQBx_$B81a5wh#)*`3<0aOU%1LVLfD97#3 z+!bwQ?MtCE&KAY~7#XGQ>kHg`J>i#XO(rSdbPor6ozZ@(mac+Mqdoh4c7FV0m7Joj zDBzgDbY#0<#LZb$c9VQ}d*Rw6Y=5Iw^0GC}th%yF&oOcV(>c6@W2DDjJ6nzs=e0`3 zRqNLJ#iz=Mtbp(bxVHLh6z+h*TJH|GAc#+2&D-*e zbIN-1CaK@}OwCp_NW~n)%INWUYkI&rB&K6{&W5zPsJn)^*;0j~xcirkb8rBTe2zjs zxV8#KIm!MP0Fg#=pVCxnO4ni4par1@E#@j?oN>*a?&oY&@GgOD^w&dMXJ)eS76?p= z3ubXg=;M{r5X@$ePyHQh3zOx3Q+;UG1)KL;Vh3s1N}vsUVF%+~_voGPUknzNb>GAG zGwpg=61Z3tu$a^NQPi?CAY!0`7{wRRV9{1qOYeHVaPSFOV*y#p;Hy8x^V+|jxdZRx zO`X#f|4;@xn3unl+m5PjbW2|oFT0xhxwvPxkKeI7)O!aGE=v9bIp!m-h<{tby%XCI zsl;~a^_b$w*vo)5iw;d&$deq*IihJcK=f<)! z$VY9Ne}3T*%i5|K5!G%nC~cw;3w0t^&w%)54s|UhP`>^;u7&LXv@ySy5VXWTJ?9k7 zqG}sol*ze{OroTmZQV*Eh42*tRB^8ZEb)7p4=;p;m?4sZ-t*GBs&s@|{a#O;VrR(3 zqP%*cFqsMb<2_64npY96?p!<#vMKV3)t%H)@DeZ7S<2xOvYb*jv>)NQczBEztp=~V*|b{;V6~}e-x!jK_3U< z#M$0fPx!@uMq{+R#_M40%mTx+^kHMDc-wh?RzJ|9tr@Yy)GanGGruVt!!!+gJ$@F| z{XC$#;ESi;ms;g)_UDYV)I6f3v%QoEHHxz%Get(gKrpQ<&FN4Vl1JTHB&6oF@(A|j zph1^FYRR_+C;`8w_O*c7`zpgLlwQlxwarOZi#9(Xg7I^Y^nTcI2^&t<-|4A{6Px4l zmLL=cV7856QqN3Cz5Dl2il zOk2BMBt29q-lO6+>jd;2f`E{38GVC_&no-Ap{gEock~t8$RKdnh}E7QFN=Y%{iV{* zJPtTR&~Z9tt*J^Oktc4=98HgGHW3_dqPGmgT}!9WK-sHo-)O(n?i$MIbS!lezrdxY zD-OT@f-|}^hIl>sOGA@@=1I?j!Pd{>X`fnTYcz0PvvcWQNojJSHm9^DH)$-vrrgY= z!!vGXZxbkHvX_kFa`TFONobGZYPT}hC@N8OA1DV-s)OyBs||~q1 z($b*BnN-Ky8N4h#4l^OykQ>C!|CE*pAynD|0ik1xneO}VXD7TJRu7?dtw|oNaJ09m zORzxQcPry(Eb1$DyFUy~d3^NBPtpB3dq-ME!yklAt+1Z})uV|y_0V2l_3ayw`yr-k z?RuI87W-o+3ydQRa*u6nw9+ex5aLINwY*e1_Kp`Cv~^(P>BP7rMq#m?M;3YBJxp*Xar78e@Bc!+B>o^yh;{w(M3Z||3Fk4fIvHCS z4D)L~cIjTO-ALAq$#67S+#txPobRxiF!CCVDk^RF9m!J+42IoZ4!`-n_`95j0-3W# z>V&0LBU^&Qyo~=sE~s5VzxSbuwQ8sm+JfEt$FL~p@$4AGruK=Bux0u_Wz^{g*%mTZ zPMvnYTl+?^xwxA`jGg?&aZMHZ3)-$E@~dG;cE*Zvj^7vMdxuiM_0PYCUT!=>=|sT& z29)jC#f*fIC)-|Y<6aBj_Kdslb4wNrTTzrk8NIR|e66B*=gusibaJkat(O>Rrs?+f zX#Vx()MNlO#PMm@yQxfo0dIKx6k`dH*NSI}7eX@%)?{+->DH;Aw+4JXon=yfSPC|l zIn`IHOB81Ug&1Y{m97r!wzm19 zFDg!cd=Ow$<2tvIaJ|IsGn)BVWeH#9A&1qGDBv7>pHMDWM4l| z;C8pwIgI4|n51R{4_yy&J&zAi!D&!MtMK?6Gr;{N!yy5zPi?arMtelC!uxib^!X?hq=WK)Kjm)0Au@gTxKGLN6u z0pB`M*pRJF9CcYQU5qxyqpQHqtIpe#B|il+kYY>zJrrD-8U@omO+)Z0 z5*4F*>z(+vATU$vBSh9|m<4X1Y8tog^^J)x>9|H0&hHxhb=kU-Lon8mb=RaTEv<_` z@($v|{fd!)M1FrnW@LJ$-+jEgV1G!=mOevcMry#x0_7^0YejwgffoPV8E>TjcRmwl zf};_Zpc^Y+wW&LtaXfrD&5eBO35usJx#y literal 0 HcmV?d00001 diff --git a/interface/EventProducer.h b/interface/EventProducer.h index af3aabb..d56acc7 100644 --- a/interface/EventProducer.h +++ b/interface/EventProducer.h @@ -2,6 +2,7 @@ #define EVENT_PRODUCER #include +#include #include #include @@ -13,6 +14,9 @@ class EventProducer: public Framework::Producer { Producer(name, tree, config) { + if (config.getUntrackedParameter("compute_pu_weights", true)) + m_pu_reweighter = std::make_shared(config.getParameterSet("pu_reweighter"), Framework::PUProfile::Run2015_25ns); + } virtual ~EventProducer() {} @@ -38,6 +42,8 @@ class EventProducer: public Framework::Producer { float m_event_weight_sum = 0; + std::shared_ptr m_pu_reweighter; + public: // Tree members @@ -49,6 +55,7 @@ class EventProducer: public Framework::Producer { BRANCH(npu, int); BRANCH(true_interactions, float); + BRANCH(pu_weight, float); BRANCH(pt_hat, float); BRANCH(weight, float); diff --git a/interface/PUReweighter.h b/interface/PUReweighter.h new file mode 100644 index 0000000..cb51fca --- /dev/null +++ b/interface/PUReweighter.h @@ -0,0 +1,86 @@ +#pragma once + +#include +#include +#include + +#include + +#include + +namespace Framework { + enum PUProfile : uint8_t { + Run2015_50ns, + Run2015_25ns + }; + + class PUReweighter { + public: + PUReweighter(const edm::ParameterSet&, PUProfile mc_pu_profile); + + float getWeight(float n_interactions); + + private: + std::shared_ptr m_pu_weights; + + std::map> m_pu_profiles { + { + Run2015_25ns, + { + 4.8551E-07, + 1.74806E-06, + 3.30868E-06, + 1.62972E-05, + 4.95667E-05, + 0.000606966, + 0.003307249, + 0.010340741, + 0.022852296, + 0.041948781, + 0.058609363, + 0.067475755, + 0.072817826, + 0.075931405, + 0.076782504, + 0.076202319, + 0.074502547, + 0.072355135, + 0.069642102, + 0.064920999, + 0.05725576, + 0.047289348, + 0.036528446, + 0.026376131, + 0.017806872, + 0.011249422, + 0.006643385, + 0.003662904, + 0.001899681, + 0.00095614, + 0.00050028, + 0.000297353, + 0.000208717, + 0.000165856, + 0.000139974, + 0.000120481, + 0.000103826, + 8.88868E-05, + 7.53323E-05, + 6.30863E-05, + 5.21356E-05, + 4.24754E-05, + 3.40876E-05, + 2.69282E-05, + 2.09267E-05, + 1.5989E-05, + 4.8551E-06, + 2.42755E-06, + 4.8551E-07, + 2.42755E-07, + 1.21378E-07, + 4.8551E-08 + } + } + }; + }; +}; diff --git a/interface/Tools.h b/interface/Tools.h index 45c19ea..91dfce0 100644 --- a/interface/Tools.h +++ b/interface/Tools.h @@ -1,5 +1,7 @@ #pragma once +#include + namespace pat { class Jet; } @@ -23,4 +25,8 @@ namespace Tools { bool passTightId(const pat::Jet& jet); bool passTightLeptonVetoId(const pat::Jet& jet); }; + + inline bool caseInsensitiveEquals(const std::string& a, const std::string& b) { + return std::equal(a.begin(), a.end(), b.begin(), [](char a, char b) { return std::tolower(a) == std::tolower(b); }); + } }; diff --git a/python/EventProducer.py b/python/EventProducer.py index 506e72d..c7f71f7 100644 --- a/python/EventProducer.py +++ b/python/EventProducer.py @@ -5,6 +5,10 @@ prefix = cms.string('event_'), enable = cms.bool(True), parameters = cms.PSet( - pu_summary = cms.untracked.InputTag('slimmedAddPileupInfo') + pu_summary = cms.untracked.InputTag('slimmedAddPileupInfo'), + compute_pu_weights = cms.untracked.bool(True), + pu_reweighter = cms.PSet( + data_pu_profile = cms.untracked.FileInPath('cp3_llbb/Framework/data/PUProfiles/pu_histogram_all_2015.root') + ) ) ) diff --git a/src/EventProducer.cc b/src/EventProducer.cc index 6219686..5e6eea5 100644 --- a/src/EventProducer.cc +++ b/src/EventProducer.cc @@ -11,6 +11,8 @@ void EventProducer::produce(edm::Event& event_, const edm::EventSetup& eventSetu rho = *rho_handle; + pu_weight = 1; + edm::Handle> pu_infos; if (event_.getByToken(m_pu_info_token, pu_infos)) { for (const auto& pu_info: *pu_infos) { @@ -19,6 +21,9 @@ void EventProducer::produce(edm::Event& event_, const edm::EventSetup& eventSetu npu = pu_info.getPU_NumInteractions(); true_interactions = pu_info.getTrueNumInteractions(); + + if (m_pu_reweighter.get()) + pu_weight = m_pu_reweighter->getWeight(true_interactions); } } diff --git a/src/PUReweighter.cc b/src/PUReweighter.cc new file mode 100644 index 0000000..49014dc --- /dev/null +++ b/src/PUReweighter.cc @@ -0,0 +1,53 @@ +#include +#include + +#include +#include + +#include + +namespace Framework { + + PUReweighter::PUReweighter(const edm::ParameterSet& config, PUProfile mc_pu_profile) { + + std::string file = config.getUntrackedParameter("data_pu_profile").fullPath(); + + // Retrieve PU histogram from file + std::unique_ptr f(TFile::Open(file.c_str())); + TH1* data_pu_histogram = static_cast(f->Get("pileup")); + + size_t x_max = (size_t) data_pu_histogram->GetXaxis()->GetXmax(); + size_t n_bins = data_pu_histogram->GetNbinsX(); + + std::shared_ptr mc_pu_histogram = std::make_shared("pu_mc", "", n_bins, 0, x_max); + mc_pu_histogram->SetDirectory(nullptr); + + const std::vector& coeffs = m_pu_profiles[mc_pu_profile]; + for (size_t i = 1; i <= (size_t) data_pu_histogram->GetNbinsX(); i++) { + size_t coef_index = (size_t) data_pu_histogram->GetBinCenter(i); + float coef = (coef_index) < coeffs.size() ? coeffs[coef_index] : 0.; + + mc_pu_histogram->SetBinContent(i, coef); + } + + data_pu_histogram->Scale(1. / data_pu_histogram->Integral()); + mc_pu_histogram->Scale(1. / mc_pu_histogram->Integral()); + + std::shared_ptr pu_weights(static_cast(data_pu_histogram->Clone())); + pu_weights->SetDirectory(nullptr); + + pu_weights->Divide(mc_pu_histogram.get()); + + m_pu_weights = pu_weights; + } + + float PUReweighter::getWeight(float n_interactions) { + + TH1* pu_weights = m_pu_weights.get(); + + if (! pu_weights) + return 1; + + return pu_weights->GetBinContent(pu_weights->GetXaxis()->FindBin(n_interactions)); + } +};