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