-
Notifications
You must be signed in to change notification settings - Fork 0
/
PUEventWeightProcessor.cc
158 lines (132 loc) · 8.15 KB
/
PUEventWeightProcessor.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#include "PUEventWeightProcessor.h"
#include <iostream>
#include "TFile.h"
#include "TH1.h"
#include "ConfigFile.h"
#include "CalibData.h"
// -----------------------------------------------------------------
PUEventWeightProcessor::PUEventWeightProcessor(const std::string& configfile, Parameters* param)
: EventProcessor("PU weighting",configfile,param), weightEvents_(false) {
ConfigFile config(configfile.c_str());
std::string histFileName = config.read<string>(name()+" histogram","");
std::string PU_mixing_era = config.read<string>(name()+" era","Flat10");
if( histFileName != "" ) {
weightEvents_ = true;
std::cout << "Weighting events to specified PU scenario" << std::endl;
std::cout << " Reading PU scenario from '" << histFileName << "'" << std::endl;
TFile file(histFileName.c_str(),"READ");
TH1 *h = 0;
file.GetObject("pileup",h);
if( h ) {
h->SetDirectory(0);
} else {
std::cerr << "ERROR in PUEventWeightProcessor: Histogram 'pileup' does not exist in file '" << histFileName << "'\n.";
std::cerr << "See https://twiki.cern.ch/twiki/bin/view/CMS/HamburgWikiAnalysisCalibration#Pile_Up_Reweighting for available input distributions." << std::endl;
// exit(1);
}
file.Close();
std::cout << " Computing weights" << std::endl;
if(PU_mixing_era=="Fall11"){
weights_ = generate_fall11_weights(h);
std::cout << "using Fall11 distribution" << std::endl;
}
else if(PU_mixing_era=="Summer12"){
weights_ = generate_summer12_weights(h);
std::cout << "using Summer12 distribution" << std::endl;
}
else {
weights_ = generate_flat10_weights(h);
std::cout << "using Flat10 distribution" << std::endl;
}
delete h;
}
}
// -----------------------------------------------------------------
int PUEventWeightProcessor::preprocess(std::vector<Event*>& data,
std::vector<Event*>& control1,
std::vector<Event*>& control2) {
if( !weightEvents_ ) return data.size();
int nProcEvts = 0; // Number of processed events
std::vector<Event*>::iterator evt1 = control1.begin();
for(; evt1 != control1.end(); ++evt1, ++nProcEvts) {
if( (*evt1)->type() == ParLimit ) continue;
short int npu = (*evt1)->nPU();
if( npu < static_cast<short int>(weights_.size()) ) {
(*evt1)->setWeight( weights_.at(npu) * ((*evt1)->weight()) );
} else {
std::cerr << "WARNING in PUEventWeightProcessor::preprocess: Number of PU vertices = " << npu << " out of histogram binning." << std::endl;
}
}
std::vector<Event*>::iterator evt2 = control2.begin();
for(; evt2 != control2.end(); ++evt2, ++nProcEvts) {
if( (*evt2)->type() == ParLimit ) continue;
short int npu = (*evt2)->nPU();
if( npu < static_cast<short int>(weights_.size()) ) {
(*evt2)->setWeight( weights_.at(npu) * ((*evt2)->weight()) );
} else {
std::cerr << "WARNING in PUEventWeightProcessor::preprocess: Number of PU vertices = " << npu << " out of histogram binning." << std::endl;
}
}
std::cout << " Applied weights for " << nProcEvts << " events in control1 and control2. \n";
return nProcEvts;
}
// Generate weights for Flat10 PU scenario for given
// data PU distribution
// Code from: https://twiki.cern.ch/twiki/bin/viewauth/CMS/PileupReweighting
// --------------------------------------------------
std::vector<double> PUEventWeightProcessor::generate_flat10_weights(const TH1* data_npu_estimated) const {
// see SimGeneral/MixingModule/python/mix_E7TeV_FlatDist10_2011EarlyData_inTimeOnly_cfi.py; copy and paste from there:
const double npu_probs[25] = {0.0698146584, 0.0698146584, 0.0698146584,0.0698146584,0.0698146584,0.0698146584,0.0698146584,0.0698146584,0.0698146584,0.0698146584,0.0698146584 /* <-- 10*/,0.0630151648,0.0526654164,0.0402754482,0.0292988928,0.0194384503,0.0122016783,0.007207042,0.004003637,0.0020278322,0.0010739954,0.0004595759,0.0002229748,0.0001028162,4.58337152809607E-05 /* <-- 24 */};
std::vector<double> result(25);
double s = 0.0;
for(int npu=0; npu<25; ++npu) {
double npu_estimated = data_npu_estimated->GetBinContent(data_npu_estimated->GetXaxis()->FindBin(npu));
result[npu] = npu_estimated / npu_probs[npu];
s += npu_estimated;
}
// normalize weights such that the total sum of weights over the whole sample is 1.0, i.e., sum_i result[i] * npu_probs[i] should be 1.0 (!)
for(int npu=0; npu<25; ++npu) {
result[npu] /= s;
}
return result;
}
// Generate weights for Fall11 Distribution PU scenario for given
// data PU distribution
// Code from: https://twiki.cern.ch/twiki/bin/viewauth/CMS/PileupReweighting
// --------------------------------------------------
std::vector<double> PUEventWeightProcessor::generate_fall11_weights(const TH1* data_npu_estimated) const {
//Distribution extracted from Fall11_PYTHIA_S6 MC
Double_t npu_probs_Fall2011[50] = {
0.00867359,0.0188749,0.0311076,0.0424069,0.0527288,0.0586238,0.0623814,0.0618389,0.0605955,0.0569491,0.0529175,0.0493315,0.0454852,0.0423341,0.0389899,0.0362847,0.0334895,0.0309533,0.0282723,0.0258533,0.0233813,0.0209824,0.0185732,0.0163106,0.0142085,0.0122677,0.0105115,0.00880698,0.00739917,0.00610915,0.00499312,0.0040509,0.00325658,0.00258698,0.00202705,0.00158217,0.00123929,0.000932184,0.000730729,0.000534201,0.000398986,0.000301543,0.000216509,0.000162587,0.000110855,8.18409e-05,5.9305e-05,4.21522e-05,3.03824e-05,2.02549e-05};
// see https://twiki.cern.ch/twiki/bin/viewauth/CMS/PileupInformation and mix_E7TeV_Fall2011_Reprocess_50ns_PoissonOOTPU_cfi.py copy and paste from there:
// Double_t npu_probs_Fall2011[50] = { 0.003388501, 0.010357558, 0.024724258, 0.042348605, 0.058279812, 0.068851751, 0.072914824, 0.071579609, 0.066811668, 0.060672356, 0.054528356, 0.04919354, 0.044886042, 0.041341896, 0.0384679, 0.035871463, 0.03341952, 0.030915649, 0.028395374, 0.025798107, 0.023237445, 0.020602754, 0.0180688, 0.015559693, 0.013211063, 0.010964293, 0.008920993, 0.007080504, 0.005499239, 0.004187022, 0.003096474, 0.002237361, 0.001566428, 0.001074149, 0.000721755, 0.000470838, 0.00030268, 0.000184665, 0.000112883, 6.74043E-05, 3.82178E-05, 2.22847E-05, 1.20933E-05, 6.96173E-06, 3.4689E-06, 1.96172E-06, 8.49283E-07, 5.02393E-07, 2.15311E-07, 9.56938E-08};
std::vector<double> result(50);
double s = 0.0;
for(int npu=0; npu<50; ++npu) {
double npu_estimated = data_npu_estimated->GetBinContent(data_npu_estimated->GetXaxis()->FindBin(npu));
result[npu] = npu_estimated / npu_probs_Fall2011[npu];
s += npu_estimated;
}
// normalize weights such that the total sum of weights over the whole sample is 1.0, i.e., sum_i result[i] * npu_probs_Fall2011[i] should be 1.0 (!)
for(int npu=0; npu<50; ++npu) {
result[npu] /= s;
}
return result;
}
std::vector<double> PUEventWeightProcessor::generate_summer12_weights(const TH1* data_npu_estimated) const {
//Distribution extracted from QCD_Pt-15to3000_Tune23_Flat_7TeV_herwigpp_Fall11-PU_S6_START44_V9B-v1
Double_t npu_probs_Summer2012[60] =
{3.31414e-05,7.21901e-05,0.000163905,0.000310388,0.000518547,0.000830137,0.00121311,0.00177762,0.00253647,0.00354002,0.00489291,0.00660114,0.00860995,0.0110814,0.0138752,0.0168897,0.02018,0.0234687,0.0268086,0.0300392,0.033068,0.0357128,0.0380327,0.039777,0.0412705,0.0423629,0.0429243,0.0430274,0.0425976,0.0418433,0.0408154,0.0391155,0.0375611,0.0355296,0.0332208,0.030871,0.0282428,0.0256578,0.0231456,0.0206244,0.0181904,0.0159587,0.0137418,0.011751,0.0100118,0.00834922,0.00698552,0.00567559,0.00459524,0.00374077,0.0029594,0.00235304,0.00185121,0.00139865,0.00109377,0.000829436,0.000619373,0.00045647,0.000343128,0.000252716};
std::vector<double> result(60);
double s = 0.0;
for(int npu=0; npu<60; ++npu) {
double npu_estimated = data_npu_estimated->GetBinContent(data_npu_estimated->GetXaxis()->FindBin(npu));
result[npu] = npu_estimated / npu_probs_Summer2012[npu];
s += npu_estimated;
}
// normalize weights such that the total sum of weights over the whole sample is 1.0, i.e., sum_i result[i] * npu_probs_Summer2012[i] should be 1.0 (!)
for(int npu=0; npu<60; ++npu) {
result[npu] /= s;
}
return result;
}