ATLAS Offline Software
Loading...
Searching...
No Matches
APWeightSum.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5#define APWeightSum_cxx
6
7#include "ReweightUtils/APWeightSum.h"
8#include "ReweightUtils/APEvtWeight.h"
9#include "ReweightUtils/APWeightEntry.h"
10#include "ReweightUtils/APEvtWeight.h"
11#include "ReweightUtils/APReweightBase.h"
12#include "THnSparse.h"
13#include <cmath>
14#include <iostream>
15#include <limits>
16
17using namespace std;
18
31
32APWeightSum::~APWeightSum() {
33 for (vector<THnSparse*>::reverse_iterator it = m_linear_uncert.rbegin(); it != m_linear_uncert.rend(); ++it) delete *it;
34 m_linear_uncert.clear();
35}
36
37void APWeightSum::AddWeightToEvt(APWeightEntry* weight) {
38 m_current_evt_weights.push_back(weight);
39}
40
41double APWeightSum::GetSumW() const {
42 return m_k_evt_weight;
43}
44
45double APWeightSum::GetSumW2() const {
46 return m_k_evt_weight2;
47}
48
52
53double APWeightSum::GetStdDev() {
54 if ( !m_isComputed ) Compute();
55 return sqrt(m_variance);
56}
57
58double APWeightSum::GetVariance() {
59 if ( !m_isComputed ) Compute();
60 return m_variance;
61}
62
66
67double APWeightSum::GetVarianceFullCorr() {
68 if ( !m_isComputed ) Compute();
69 return m_variance_fullcorr;
70}
71
72double APWeightSum::GetSysUncert() const {
73 return sqrt(m_variance_sys);
74}
75
76unsigned long APWeightSum::GetKUnweighted() const {
77 return m_k_evt_orig;
78}
79
80void APWeightSum::FinishEvt(double ext_weight) {
81 ++m_k_evt_orig;
82 double evt_weight = 1.0;
83 double uncert = 0.0;
84 double uncert_sys = 0.0;
85 for (unsigned int i = 0, I = m_current_evt_weights.size(); i < I; ++i) {
86 double uncert_summand = 1.0;
87 evt_weight *= (1. - m_current_evt_weights[i]->GetExpectancy());
88 for (unsigned int k = 0; k < I; ++k) {
89 if (i != k) uncert_summand *= (1. - m_current_evt_weights[i]->GetExpectancy());
90 }
91 uncert += (uncert_summand * uncert_summand * m_current_evt_weights[i]->GetVariance());
92 uncert_sys += (uncert_summand * uncert_summand * m_current_evt_weights[i]->GetSysUncert2());
93 }
94 m_k_evt_weight += ext_weight * (1. - evt_weight);
95 m_k_evt_weight2 += m_k_evt_weight*m_k_evt_weight;
96 m_k_evt_weight_external += ext_weight;
97 m_variance += fabs(ext_weight) * uncert;
98 m_variance_sys += ext_weight * uncert_sys;
99 m_current_evt_weights.clear();
100}
101
102void APWeightSum::AddEvt(APEvtWeight* evt_weight, double ext_weight) {
103 ++m_k_evt_orig;
104 m_k_evt_weight += ext_weight * evt_weight->GetWeight();
105 m_k_evt_weight2 += m_k_evt_weight*m_k_evt_weight;
106 m_k_evt_weight_external += ext_weight;
107
108 vector<APWeightEntry*> temp_vec_mu = evt_weight->GetWeightObjects(APEvtWeight::kMuon);
109 vector<APWeightEntry*> temp_vec_dimu = evt_weight->GetWeightObjects(APEvtWeight::kDiMuon);
110 vector<APWeightEntry*> temp_vec_mumo = evt_weight->GetWeightObjects(APEvtWeight::kMuonMO);
111 vector<APWeightEntry*> temp_vec_tau = evt_weight->GetWeightObjects(APEvtWeight::kTau);
112 vector<APWeightEntry*> temp_vec_ditau = evt_weight->GetWeightObjects(APEvtWeight::kDiTau);
113 vector<APWeightEntry*> temp_vec_taumo = evt_weight->GetWeightObjects(APEvtWeight::kTauMO);
114 vector<APWeightEntry*> temp_vec_el = evt_weight->GetWeightObjects(APEvtWeight::kElectron);
115 vector<APWeightEntry*> temp_vec_diel = evt_weight->GetWeightObjects(APEvtWeight::kDiElectron);
116 vector<APWeightEntry*> temp_vec_elmo = evt_weight->GetWeightObjects(APEvtWeight::kElectronMO);
117 vector<APWeightEntry*> temp_vec_jet = evt_weight->GetWeightObjects(APEvtWeight::kJet);
118 vector<APWeightEntry*> temp_vec_dijet = evt_weight->GetWeightObjects(APEvtWeight::kDiJet);
119 vector<APWeightEntry*> temp_vec_jetmo = evt_weight->GetWeightObjects(APEvtWeight::kJetMO);
120
121 vector< vector<APWeightEntry*> > temp_vec_all{
122 temp_vec_mu, temp_vec_tau, temp_vec_el, temp_vec_jet,
123 std::move(temp_vec_mumo), std::move(temp_vec_taumo), std::move(temp_vec_elmo), std::move(temp_vec_jetmo),
124 temp_vec_dimu, temp_vec_ditau, temp_vec_diel, temp_vec_dijet
125 };
126
127
128 /* check if histogram for error propagation is already there; if not, create it */
129 for( unsigned int iAll = 0, IAll = temp_vec_all.size(); iAll < IAll; ++iAll ) {
130 for( unsigned int i = 0, I = temp_vec_all[iAll].size(); i < I; ++i ) {
131 unsigned int ID = temp_vec_all[iAll][i]->GetID();
132 if( m_linear_uncert.size() < ID+1 ) m_linear_uncert.resize(ID+1, 0);
133 if( m_linear_uncert[ID] == 0 ) {
134 vector<int> original_dimensions = temp_vec_all[iAll][i]->GetOriginalDimensions();
135 int *bins = new int[original_dimensions.size()];
136 double *xmin = new double[original_dimensions.size()];
137 double *xmax = new double[original_dimensions.size()];
138 for( unsigned int j = 0, J = original_dimensions.size(); j < J; ++j ) {
139 bins[j] = original_dimensions[j];
140 xmin[j] = 0.;
141 xmax[j] = 10.;
142 }
143 m_linear_uncert[ID] = new THnSparseD("","",original_dimensions.size(), bins, xmin, xmax);
144 }
145 }
146 }
147 /* calculate weight and derivatives */
148 /* these are the simple weights: kMuon, kTau, kElectron and kJet: single object trigger, OR of all elements */
149 if(evt_weight->GetType() <= APEvtWeight::kJet) {
150 vector<APWeightEntry*> temp_vec_rel;
151 if(evt_weight->GetType() == APEvtWeight::kMuon) temp_vec_rel = std::move(temp_vec_mu);
152 else if(evt_weight->GetType() == APEvtWeight::kTau) temp_vec_rel = std::move(temp_vec_tau);
153 else if(evt_weight->GetType() == APEvtWeight::kElectron) temp_vec_rel = std::move(temp_vec_el);
154 else if(evt_weight->GetType() == APEvtWeight::kJet) temp_vec_rel = std::move(temp_vec_jet);
155
156 for (unsigned int i = 0, I = temp_vec_rel.size(); i < I; ++i ) {
157 vector<int> coord = temp_vec_rel[i]->GetCoords();
158 double weight_uncert = sqrt(temp_vec_rel[i]->GetVariance());
159 for (unsigned int j = 0; j < I; ++j ) {
160 if (j == i) continue;
161 weight_uncert *= (1.0 - temp_vec_rel[j]->GetExpectancy());
162 }
163 m_linear_uncert[temp_vec_rel[i]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_rel[i]->GetID()]->GetBinContent(&coord.front())+weight_uncert);
164 m_variance_nocorr += weight_uncert*weight_uncert;
165 }
166 }
167
168 /* these are the DiMuon, DiTau, DiElectron and DiJet weights */
169 else if (evt_weight->GetType() >= APEvtWeight::kDiMuon && evt_weight->GetType() <= APEvtWeight::kDiJet) {
170 vector<APWeightEntry*> temp_vec_rel;
171 if(evt_weight->GetType() == APEvtWeight::kDiMuon) temp_vec_rel = std::move(temp_vec_dimu);
172 else if(evt_weight->GetType() == APEvtWeight::kDiTau) temp_vec_rel = std::move(temp_vec_ditau);
173 else if(evt_weight->GetType() == APEvtWeight::kDiElectron) temp_vec_rel = std::move(temp_vec_diel);
174 else if(evt_weight->GetType() == APEvtWeight::kDiJet) temp_vec_rel = std::move(temp_vec_dijet);
175
176 bool isAsymTrig = false;
177 vector<unsigned int> temp_vec_IDs;
178 temp_vec_IDs.push_back(temp_vec_rel[0]->GetID() );
179 for( unsigned int i = 1, I = temp_vec_rel.size(); i < I; ++i ) {
180 bool knownID = false;
181 for( unsigned int j = 0, J = temp_vec_IDs.size(); j < J; ++j ) {
182 if( temp_vec_rel[i]->GetID() == temp_vec_IDs[j] ) { knownID = true; break; }
183 }
184 if( !knownID ) temp_vec_IDs.push_back( temp_vec_rel[i]->GetID() );
185 }
186 if( temp_vec_IDs.size() != 1 ) isAsymTrig = true;
187
188 /* this is for symmetric dilepton triggers */
189 if( !isAsymTrig ) {
190 for (unsigned int i = 0, I = temp_vec_rel.size(); i < I; ++i ) {
191 vector<int> coord = temp_vec_rel[i]->GetCoords();
192 double weight_uncert = sqrt(temp_vec_rel[i]->GetVariance());
193 double weight_derivative = 0.;
194 for (unsigned int j = 0; j < I; ++j ) {
195 if (j == i) continue;
196 double weight_derivative_temp = temp_vec_rel[j]->GetExpectancy();
197 for (unsigned int k = 0; k < I; ++k ) {
198 if( k == j || k == i ) continue;
199 weight_derivative_temp *= (1.0 - temp_vec_rel[k]->GetExpectancy());
200 }
201 weight_derivative += weight_derivative_temp;
202 }
203 weight_uncert *= weight_derivative;
204 m_linear_uncert[temp_vec_rel[i]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_rel[i]->GetID()]->GetBinContent(&coord.front())+weight_uncert);
205 m_variance_nocorr += weight_uncert*weight_uncert;
206 }
207 }
208
209 /* this is for asymmetric triggers */
210 else {
211 /* at first the first leg of the trigger */
212 for( unsigned int k = 0, K = temp_vec_rel.size(); k < K; k += 4 ) {
213 vector<int> coord = temp_vec_rel[k]->GetCoords();
214 double variance_k = 0.;
215 double variance_temp = 1.;
216 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
217 if( i == k ) continue;
218 variance_temp *= (1. - temp_vec_rel[i]->GetExpectancy());
219 }
220 variance_k += variance_temp;
221 variance_temp = -1.;
222 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
223 variance_temp *= (1.0 - temp_vec_rel[i+2]->GetExpectancy());
224 if( i == k ) continue;
225 variance_temp *= (1.0 - temp_vec_rel[i]->GetExpectancy());
226 }
227 variance_k += variance_temp;
228 variance_temp = -temp_vec_rel[k+3]->GetExpectancy();
229 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
230 if( i == k ) continue;
231 variance_temp *= (1.0 - temp_vec_rel[i]->GetExpectancy())*(1.0 - temp_vec_rel[i+2]->GetExpectancy());
232 }
233 variance_k += variance_temp;
234 variance_temp = 0.;
235 for( unsigned int j = 0, J = temp_vec_rel.size(); j < J; j+= 4 ) {
236 if( j == k ) continue;
237 double variance_ijk_temp = temp_vec_rel[j]->GetExpectancy()*temp_vec_rel[j+3]->GetExpectancy();
238 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
239 if( i == j ) continue;
240 variance_ijk_temp *= (1.0 - temp_vec_rel[i+2]->GetExpectancy());
241 if( i == k ) continue;
242 variance_ijk_temp *= (1.0 - temp_vec_rel[i]->GetExpectancy());
243 }
244 variance_temp += variance_ijk_temp;
245 }
246 variance_k += variance_temp;
247 variance_k *= variance_k*temp_vec_rel[k]->GetVariance();
248 m_linear_uncert[temp_vec_rel[k]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_rel[k]->GetID()]->GetBinContent(&coord.front())+sqrt(variance_k));
249 m_variance_nocorr += variance_k;
250 }
251
252 /* second leg */
253 for( unsigned int k = 0, K = temp_vec_rel.size(); k < K; k += 4 ) {
254 vector<int> coord = temp_vec_rel[k+1]->GetCoords();
255 double variance_temp = 1.;
256 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
257 if( i == k ) continue;
258 variance_temp *= (1. - temp_vec_rel[i+1]->GetExpectancy());
259 }
260 variance_temp *= variance_temp*temp_vec_rel[k+1]->GetVariance();
261 m_linear_uncert[temp_vec_rel[k+1]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_rel[k+1]->GetID()]->GetBinContent(&coord.front())+sqrt(variance_temp));
262 m_variance_nocorr += variance_temp;
263 }
264
265 /* second leg | condition leg 1 failed */
266 for( unsigned int k = 0, K = temp_vec_rel.size(); k < K; k += 4 ) {
267 vector<int> coord = temp_vec_rel[k+2]->GetCoords();
268 double variance_k = 0.;
269 double variance_temp = 1.;
270 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
271 variance_temp *= (1. - temp_vec_rel[i]->GetExpectancy());
272 if( i == k ) continue;
273 variance_temp *= (1. - temp_vec_rel[i+2]->GetExpectancy());
274 }
275 variance_k += variance_temp;
276 variance_temp = 0.;
277 for( unsigned int j = 0, J = temp_vec_rel.size(); j < J; j+= 4 ) {
278 double variance_ijk_temp = temp_vec_rel[j]->GetExpectancy()*temp_vec_rel[j+3]->GetExpectancy();
279 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
280 if( j == i ) continue;
281 variance_ijk_temp *= (1.0 - temp_vec_rel[i]->GetExpectancy());
282 if( i == k ) continue;
283 variance_ijk_temp *= (1.0 - temp_vec_rel[i+2]->GetExpectancy());
284 }
285 variance_temp += variance_ijk_temp;
286 }
287 variance_k += variance_temp;
288 variance_k *= variance_k*temp_vec_rel[k+2]->GetVariance();
289 m_linear_uncert[temp_vec_rel[k+2]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_rel[k+2]->GetID()]->GetBinContent(&coord.front())+sqrt(variance_k));
290 m_variance_nocorr += variance_k;
291 }
292
293 /* second leg | condition leg 1 passed */
294 for( unsigned int k = 0, K = temp_vec_rel.size(); k < K; k += 4 ) {
295 vector<int> coord = temp_vec_rel[k+3]->GetCoords();
296 double variance_k = - temp_vec_rel[k]->GetExpectancy();
297 for( unsigned int i = 0, I = temp_vec_rel.size(); i < I; i += 4 ) {
298 if( i == k ) continue;
299 variance_k *= (1.0 - temp_vec_rel[i]->GetExpectancy())*(1.0 - temp_vec_rel[i+3]->GetExpectancy());
300 }
301 variance_k *= variance_k*temp_vec_rel[k+3]->GetVariance();
302 m_linear_uncert[temp_vec_rel[k+3]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_rel[k+3]->GetID()]->GetBinContent(&coord.front())+sqrt(variance_k));
303 m_variance_nocorr += variance_k;
304 }
305 }
306
307 }
308
309 /* these are the weights for ANDed triggers. Shouldn't be used but for examples. This will treat everything as ANDed. */
310 else if (evt_weight->GetType() == APEvtWeight::kANDed ) {
311 for (unsigned int iAll = 0, IAll = temp_vec_all.size(); iAll < IAll; ++iAll ) {
312 for (unsigned int i = 0, I = temp_vec_all[iAll].size(); i < I; ++i ) {
313 vector<int> coord = temp_vec_all[iAll][i]->GetCoords();
314 double weight_uncert = sqrt(temp_vec_all[iAll][i]->GetVariance());
315 weight_uncert *= (temp_vec_all[iAll][i]->GetExpectancy() <= numeric_limits<double>::epsilon() ) ? 0. : evt_weight->GetWeight()/temp_vec_all[iAll][i]->GetExpectancy();
316 m_linear_uncert[temp_vec_all[iAll][i]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_all[iAll][i]->GetID()]->GetBinContent(&coord.front())+weight_uncert);
317 m_variance_nocorr += weight_uncert*weight_uncert;
318 }
319 }
320 }
321
322 /* these are weights for ORed triggers, for instance dimuon OR single muon */
323 /* treat everything as ORed: calculate weight for each type, then combine weights */
324 /* objects of same type are ORed: kMuon, kTau, kElectron, kJet, kMuonMO, kTauMO, kElectronMO, kJetMo */
325 /* objects of different type are handled differently */
326 /* if type is kMOANDed, everything is ANDed.*/
327 /* if type is kMOORed, kMuonMO, kTauMO, kElectronMO, kJetMO are ANDed; all others are ORed with each other! */
328 else if (evt_weight->GetType() >= APEvtWeight::kORed ) {
329
330 vector<double> temp_weight_rel(12,1.);
331
332 /* calculate single object trigger weights: OR of all objects of same type: kMuon, kTau, kElectron, kJet, kMuonMO, kTauMO, kElectronMO, kJetMO */
333 for (unsigned int j = 0; j < 8; ++j) {
334 for (unsigned int i = 0, I = temp_vec_all[j].size(); i < I; ++i ) temp_weight_rel[j] *= (1.0 - temp_vec_all[j][i]->GetExpectancy());
335 // for triggers intended to be used in ORing, set weight to 1. if there is no object and an AND is require: (if no object is added, assume none is required)
336 if( j < 4 ) temp_weight_rel[j] = (temp_vec_all[j].size() > 0 || evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed) ? (1.0 - temp_weight_rel[j]) : 1.;
337 // for triggers inteded to be used in ANDing, set weight to 1. if there is no object added (if none is added assume none is required)
338 else temp_weight_rel[j] = (temp_vec_all[j].size() > 0) ? (1.0 - temp_weight_rel[j]) : 1.;
339 }
340
341
342 /* calculate diobject trigger weights: OR of all objects of same type but require at least 2: kDiMuon, kDiTau, kDiElectron, kDiJet */
343 for (unsigned int j = 8; j < 12; ++j) {
344 if( temp_vec_all[j].size() >= 2 ) {
345 for (unsigned int i = 0, I = temp_vec_all[j].size(); i < I; ++i) temp_weight_rel[j] *= (1.0 - temp_vec_all[j][i]->GetExpectancy());
346 for (unsigned int i = 0, I = temp_vec_all[j].size(); i < I; ++i) {
347 double temp_weight = temp_vec_all[j][i]->GetExpectancy();
348 for (unsigned int k = 0; k < I; ++k ) {
349 if( k == i ) continue;
350 temp_weight *= (1.0 - temp_vec_all[j][k]->GetExpectancy());
351 }
352 temp_weight_rel[j] += temp_weight;
353 }
354 }
355 temp_weight_rel[j] = (temp_vec_all[j].size() > 0 || evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed) ? (1.0 - temp_weight_rel[j]) : 1.;
356 }
357
358 /* calculat weight for multiobject trigger (AND of all MO objects) */
359 double temp_weight_MO = 1.;
360 int n_noObject_MO = 0;
361 for (unsigned int l = 4; l < 8; ++l ) {
362 temp_weight_MO *= temp_weight_rel[l];
363 if( temp_vec_all[l].size() == 0 ) n_noObject_MO += 1;
364 }
365 if( /*(evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed) &&*/ n_noObject_MO >= 2 ) temp_weight_MO = 0.;
366
367 /* all partial weights calculated; now calculate uncertainties */
368 /* single object triggers */
369 for (unsigned int j = 0; j < 8; ++j) {
370 for (unsigned int i = 0, I = temp_vec_all[j].size(); i < I; ++i ) {
371 vector<int> coord = temp_vec_all[j][i]->GetCoords();
372 double weight_uncert = sqrt(temp_vec_all[j][i]->GetVariance());
373 for (unsigned int k = 0; k < I; ++k ) {
374 if( k == i ) continue;
375 weight_uncert *= (1.0 - temp_vec_all[j][k]->GetExpectancy());
376 }
377 for( unsigned int l = 0; l < 4; ++l ) {
378 if( l == j ) continue;
379 if( evt_weight->GetType() == APEvtWeight::kMOANDed ) weight_uncert *= temp_weight_rel[l];
380 else if( evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed ) weight_uncert *= (1.0 - temp_weight_rel[l]);
381 else cout << "WARNING: handling for this weight type is unknown! uncertainties will be incorrect" << endl;
382 }
383 for( unsigned int l = 8; l < 12; ++l ) {
384 if( l == j ) continue;
385 if( evt_weight->GetType() == APEvtWeight::kMOANDed ) weight_uncert *= temp_weight_rel[l];
386 else if( evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed ) weight_uncert *= (1.0 - temp_weight_rel[l]);
387 else cout << "WARNING: handling for this weight type is unknown! uncertainties will be incorrect" << endl;
388 }
389
390 if( j < 4 || j > 7 ) {
391 if( evt_weight->GetType() == APEvtWeight::kMOANDed ) weight_uncert *= temp_weight_MO;
392 else if( evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed ) weight_uncert *= (1.0 - temp_weight_MO);
393 else cout << "WARNING: handling for this weight type is unknown! uncertainties will be incorrect" << endl;
394 }
395 else if( j >= 4 && j <= 7 && temp_weight_rel[j] > numeric_limits<double>::epsilon() ) weight_uncert = weight_uncert*temp_weight_MO/temp_weight_rel[j];
396 m_linear_uncert[temp_vec_all[j][i]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_all[j][i]->GetID()]->GetBinContent(&coord.front())+weight_uncert);
397 m_variance_nocorr += weight_uncert*weight_uncert;
398 }
399 }
400
401 /* diobject triggers */
402 for (unsigned int j = 8; j < 12; ++j) {
403 if( temp_vec_all[j].size() >= 2 ) {
404 for (unsigned int i = 0, I = temp_vec_all[j].size(); i < I; ++i ) {
405 vector<int> coord = temp_vec_all[j][i]->GetCoords();
406 double weight_uncert = sqrt(temp_vec_all[j][i]->GetVariance());
407 double weight_derivative = 0.;
408 for (unsigned int k = 0; k < I; ++k ) {
409 if( k == i ) continue;
410 double weight_derivative_temp = temp_vec_all[j][k]->GetExpectancy();
411 for (unsigned int l = 0; l < I; ++l ) {
412 if( l == k || l == i ) continue;
413 weight_derivative_temp *= (1.0 - temp_vec_all[j][l]->GetExpectancy());
414 }
415 weight_derivative += weight_derivative_temp;
416 }
417 weight_uncert *= weight_derivative;
418 for( unsigned int l = 0; l < 4; ++l ) {
419 // l cannot be equal to j at this point
420 if( evt_weight->GetType() == APEvtWeight::kMOANDed ) weight_uncert *= temp_weight_rel[l];
421 else if( evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed ) weight_uncert *= (1.0 - temp_weight_rel[l]);
422 else cout << "WARNING: handling for this weight type is unknown! uncertainties will be incorrect" << endl;
423 }
424 for( unsigned int l = 8; l < 12; ++l ) {
425 if( l == j ) continue;
426 if( evt_weight->GetType() == APEvtWeight::kMOANDed ) weight_uncert *= temp_weight_rel[l];
427 else if( evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed ) weight_uncert *= (1.0 - temp_weight_rel[l]);
428 else cout << "WARNING: handling for this weight type is unknown! uncertainties will be incorrect" << endl;
429 }
430 if( evt_weight->GetType() == APEvtWeight::kMOANDed ) weight_uncert *= temp_weight_MO;
431 else if( evt_weight->GetType() == APEvtWeight::kORed || evt_weight->GetType() == APEvtWeight::kMOORed ) weight_uncert *= (1.0 - temp_weight_MO);
432 else cout << "WARNING: handling for this weight type is unknown! uncertainties will be incorrect" << endl;
433
434 //j is equal or greater than 8 in this loop
435 m_linear_uncert[temp_vec_all[j][i]->GetID()]->SetBinContent(&coord.front(), m_linear_uncert[temp_vec_all[j][i]->GetID()]->GetBinContent(&coord.front())+weight_uncert);
436 m_variance_nocorr += weight_uncert*weight_uncert;
437 }
438 }
439 }
440 }
441
442 m_variance_sys += ext_weight * evt_weight->GetSysVariance();
443
444 m_isComputed = false;
445}
446
447void APWeightSum::Compute() {
448 m_variance = 0.;
449 m_variance_fullcorr = 0.;
450 for (unsigned int iLinearUncert = 0, ILinearUncert = m_linear_uncert.size(); iLinearUncert < ILinearUncert; ++iLinearUncert) {
451 if( m_linear_uncert[iLinearUncert] == 0 ) continue;
452 double temp_linear_uncerts = 0;
453 for (unsigned int i = 0, I = m_linear_uncert[iLinearUncert]->GetNbins(); i < I; ++i ) {
454 double bin_content = m_linear_uncert[iLinearUncert]->GetBinContent(i);
455 m_variance += bin_content*bin_content;
456 temp_linear_uncerts += bin_content;
457 }
458 m_variance_fullcorr += temp_linear_uncerts*temp_linear_uncerts;
459 }
460
461 m_isComputed = true;
462}
463
464THnSparse* APWeightSum::GetUncertHistogram( APReweightBase* weighter ) {
465 unsigned int temp_ID = weighter->GetID();
466 if( temp_ID > m_linear_uncert.size()-1 ) {
467 cout << "WARNING: ID unknown. Returning 0-pointer!" << endl;
468 return 0;
469 }
470 else return m_linear_uncert[temp_ID];
471}
472
473const vector<THnSparse*> & APWeightSum::GetAllUncertHistograms( ) {
474 return m_linear_uncert;
475}
coord
double coord
Type of coordination system.
Definition JetVoronoiDiagramHelpers.h:43
I
#define I(x, y, z)
Definition MD5.cxx:116
bins
static const std::vector< std::string > bins
Definition MuonTriggerSFFilesTest.cxx:22
size
size_t size() const
Number of registered mappings.
APEvtWeight
Class to calculate the sum of weights ("weighted counter").
Definition APEvtWeight.h:26
APEvtWeight::GetType
ObjType GetType()
Returns the type of the event weight (muon, electron, jet, ANDed, ORed).
Definition APEvtWeight.cxx:263
APEvtWeight::GetWeight
double GetWeight()
Returns the event weight.
Definition APEvtWeight.cxx:223
APEvtWeight::GetWeightObjects
std::vector< APWeightEntry * > GetWeightObjects(ObjType type)
Returns the vector of weight objects for a specific object type.
Definition APEvtWeight.cxx:251
APEvtWeight::GetSysVariance
double GetSysVariance()
Returns the systematic variance (from systematics assigned to weights).
Definition APEvtWeight.cxx:246
APReweightBase
Class to provide common methods for all Reweight classes.
Definition APReweightBase.h:23
APReweightBase::GetID
unsigned int GetID() const
Returns the unique ID for assignment of APWeightEntries to source.
Definition APReweightBase.cxx:23
APWeightEntry
Class to store a single weight entry (one bin).
Definition APWeightEntry.h:25
APWeightSum::m_k_evt_weight
double m_k_evt_weight
Holds the sum of weights.
Definition APWeightSum.h:61
APWeightSum::GetSumW2
double GetSumW2() const
Returns sum of (weights^2).
Definition APWeightSum.cxx:45
APWeightSum::m_k_evt_weight_external
double m_k_evt_weight_external
Holds the sum of external weights (no trigger weighting).
Definition APWeightSum.h:63
APWeightSum::AddEvt
void AddEvt(APEvtWeight *evt_weight, double ext_weight=1.0)
Adds an event with an externally calculated EvtWeight object.
Definition APWeightSum.cxx:102
APWeightSum::GetVarianceFullCorr
double GetVarianceFullCorr()
Returns the variance, assuming full correlation amongst objects.
Definition APWeightSum.cxx:67
APWeightSum::GetVariance
double GetVariance()
Returns the variance.
Definition APWeightSum.cxx:58
APWeightSum::GetSumWExternal
double GetSumWExternal() const
Returns the sum of weights without taking into account the trigger weighting (external weights only) ...
Definition APWeightSum.cxx:49
APWeightSum::GetVarianceNoCorr
double GetVarianceNoCorr()
Returns the variance, assuming no correlations.
Definition APWeightSum.cxx:63
APWeightSum::m_k_evt_orig
unsigned long int m_k_evt_orig
Holds the original amount of unweighted counts ("sum of 1's").
Definition APWeightSum.h:60
APWeightSum::GetUncertHistogram
THnSparse * GetUncertHistogram(APReweightBase *weighter)
Returns THnSparse holding the uncertainties for given APReweightBase instance.
Definition APWeightSum.cxx:464
APWeightSum::m_variance
double m_variance
Holds the variance.
Definition APWeightSum.h:64
APWeightSum::m_variance_sys
double m_variance_sys
Holds the systematic variance (from systematics assigned to weights).
Definition APWeightSum.h:67
APWeightSum::~APWeightSum
virtual ~APWeightSum()
Default destructor.
Definition APWeightSum.cxx:32
APWeightSum::m_isComputed
bool m_isComputed
Definition APWeightSum.h:68
APWeightSum::GetStdDev
double GetStdDev()
Returns the standard deviation.
Definition APWeightSum.cxx:53
APWeightSum::m_variance_fullcorr
double m_variance_fullcorr
Holds the variance, assuming full correlation amongst objects.
Definition APWeightSum.h:66
APWeightSum::m_k_evt_weight2
double m_k_evt_weight2
Holds the sum of squared weights.
Definition APWeightSum.h:62
APWeightSum::m_variance_nocorr
double m_variance_nocorr
Holds the variance, assuming no correlations.
Definition APWeightSum.h:65
APWeightSum::GetKUnweighted
unsigned long GetKUnweighted() const
Returns the unweighted sum of entries.
Definition APWeightSum.cxx:76
APWeightSum::AddWeightToEvt
void AddWeightToEvt(APWeightEntry *weight)
Adds a weight to the sum of weights.
Definition APWeightSum.cxx:37
APWeightSum::m_linear_uncert
std::vector< THnSparse * > m_linear_uncert
Holds all histograms for uncertainties.
Definition APWeightSum.h:59
APWeightSum::GetSysUncert
double GetSysUncert() const
Returns the systematic uncertainty (from systematics assigned to weights).
Definition APWeightSum.cxx:72
APWeightSum::GetSumW
double GetSumW() const
Returns the sum of weights.
Definition APWeightSum.cxx:41
APWeightSum::FinishEvt
void FinishEvt(double ext_weight=1.0)
Finishes the current event and calculates the event weight.
Definition APWeightSum.cxx:80
APWeightSum::APWeightSum
APWeightSum()
Default constructor.
Definition APWeightSum.cxx:19
APWeightSum::m_current_evt_weights
ClassDef(APWeightSum, 1) protected std::vector< APWeightEntry * > m_current_evt_weights
< Calculates the final uncertainties including correlations.
Definition APWeightSum.h:53
APWeightSum::GetAllUncertHistograms
const std::vector< THnSparse * > & GetAllUncertHistograms()
Returns vector of THnSparses holding the uncertainties for all APReweight IDs.
Definition APWeightSum.cxx:473
xmax
double xmax
Definition listroot.cxx:61
xmin
double xmin
Definition listroot.cxx:60
std
STL namespace.
Generated on for ATLAS Offline Software by doxygen 1.17.0