ATLAS Offline Software
TPileupReweighting.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /******************************************************************************
6 Name: TPileupReweighting
7 
8 Author: Will Buttinger
9 Created: October 2011
10 
11 Description: Tool to get the calculated MC pileup weight.
12 ******************************************************************************/
13 
14 // Preprocessor magic for debugging
15 #define XXX std::cout << " I am here: " << __FILE__ << ":" << __LINE__ << std::endl;
16 
17 // This class' header
19 
20 
21 
22 
23 
24 // ROOT includes
25 #include <TROOT.h>
26 #include <TChain.h>
27 #include <TLeaf.h>
28 #include <TFile.h>
29 #include <TH1.h>
30 #include <TH1D.h>
31 #include <TH2D.h>
32 #include <TH3D.h>
33 #include <TAxis.h>
34 #include <TString.h>
35 #include <TRandom3.h>
36 #include <TObjString.h>
37 // include math
38 #include <cmath>
39 #include <algorithm>
40 
42 
43 
44 //=============================================================================
45 // Constructor
46 //=============================================================================
48  TNamed(name,"notitle"), m_parentTool(this),
49  m_SetWarnings(true),m_debugging(false),m_printInfo(true),
50  m_countingMode(true),m_unrepresentedDataAction(0),m_isInitialized(false),m_lumiVectorIsLoaded(false),
51  m_dataScaleFactorX(1.),m_dataScaleFactorY(1.),
52  m_mcScaleFactorX(1.),m_mcScaleFactorY(1.),
53  m_nextPeriodNumber(1),m_ignoreFilePeriods(false),m_metadatatree(0),m_unrepDataTolerance(0.05),m_doGlobalDataWeight(false),m_lumicalcRunNumberOffset(0), m_emptyHistogram(), m_random3(), m_ignoreBadChannels(false)
54 {
55  m_random3.reset( new TRandom3(0) );
56  m_random3->SetSeed(1);
57 
58  //create the global period
59  m_periodList.emplace_back( -1,0,9999999,0/* the global default channel*/ );
60  m_periods[ -1 ] = &( m_periodList.back() );
61 }
62 
63 void CP::TPileupReweighting::RemapPeriod(Int_t periodNumber1, Int_t periodNumber2) {
64  //check if periodNumber2 exists
65  if(m_periods.find(periodNumber2)==m_periods.end()) {
66  m_periodList.emplace_back( periodNumber2, 0, 0, GetDefaultChannel(-1) );
67  m_periods[periodNumber2] = &( m_periodList.back() );
68  }
69  m_periods[periodNumber1] = m_periods[periodNumber2];
70 }
71 
72 void CP::TPileupReweighting::SetDefaultChannel(Int_t channel, Int_t periodNumber) {
73  //check if periodNumber2 exists
74  if(m_periods.find(periodNumber)==m_periods.end()) {
75  m_periodList.emplace_back( periodNumber, 0, 0, channel );
76  m_periods[periodNumber] = &( m_periodList.back() );
77  } else {
78  m_periods[periodNumber]->SetDefaultChannel(channel);
79  }
80 }
81 
82 
83 Int_t CP::TPileupReweighting::SetBinning(Int_t nbinsx, Double_t* xbins, Int_t nbinsy, Double_t* ybins) {
84  if( nbinsy > 0 ) {
85  m_emptyHistogram.reset( new TH2D( "default", "default", nbinsx, xbins,
86  nbinsy, ybins ) );
87  } else {
88  m_emptyHistogram.reset( new TH1D( "default", "default", nbinsx, xbins ) );
89  }
90  m_emptyHistogram->SetDirectory(0);
91  return 0;
92 }
93 Int_t CP::TPileupReweighting::SetUniformBinning(Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy, Double_t ylow, Double_t yup) {
94  std::vector<double> xbins(nbinsx+1);
95  std::vector<double> ybins(nbinsy+1);
96  for(int i=0;i<nbinsx+1;i++) xbins[i] = xlow + i*(xup-xlow)/nbinsx;
97  if (nbinsy > 0)
98  for(int i=0;i<nbinsy+1;i++) ybins[i] = ylow + i*(yup-ylow)/nbinsy;
99  return SetBinning(nbinsx,&xbins[0],nbinsy,&ybins[0]);
100 }
102  if(!hist) return -1;
103  m_emptyHistogram.reset( dynamic_cast< TH1* >( hist->Clone( "default" ) ) );
104  if (m_emptyHistogram)
105  m_emptyHistogram->SetDirectory(0);
106  return 0;
107 }
108 
110  if(m_parentTool != this) return m_parentTool->GetDefaultChannel(mcRunNumber);
111  return m_periods[mcRunNumber]->defaultChannel;
112 }
113 
114 
115 Double_t CP::TPileupReweighting::GetIntegratedLumi(const TString& trigger) {
116  if(!m_isInitialized) {
117  if (m_printInfo) {
118  Info("GetIntegratedLumi","Initializing the subtool..");
119  }
120  Initialize();
121  }
122  if(!m_lumiVectorIsLoaded) {
123  Error("GetIntegratedLumi","No UNPRESCALED (Trigger=None) Lumicalc file loaded, so cannot get integrated lumi, returning 0");
124  return 0;
125  }
126  if(trigger=="" || trigger=="None") return GetSumOfEventWeights(-1)/1E6;
127 
128  if(m_triggerObjs.find(trigger)==m_triggerObjs.end()) {
129  CalculatePrescaledLuminosityHistograms(trigger);
130  }
131 
132  auto& h = m_triggerObjs[trigger]->triggerHists[-1][m_triggerObjs[trigger]->getBits(this)];
133 
134  return h->Integral(0,h->GetNbinsX()+1)/1E6;
135 }
136 
137 
138 
139 Double_t CP::TPileupReweighting::GetIntegratedLumi(Int_t periodNumber, UInt_t start, UInt_t end) {
140  if(!m_isInitialized) {
141  if (m_printInfo) {
142  Info("GetIntegratedLumi","Initializing the subtool..");
143  }
144  Initialize();
145  }
146  //look through dataPeriodRunTotals["pileup"][periodNumber] for runs inside the given period
147  double total = 0;
148  for(auto run : m_periods[periodNumber]->runNumbers) {
149  if(run >= start && run <= end) total += m_runs[run].inputHists["None"]->GetSumOfWeights();
150  }
151  return total*1E-6;
152 }
153 
155  if(!m_isInitialized) {
156  if (m_printInfo) {
157  Info("GetIntegratedLumi","Initializing the subtool..");
158  }
159  Initialize();
160  }
161  if(m_runs.find(runNumber)==m_runs.end()) return 0.;
162  auto& run = m_runs[runNumber];
163 
164  if(run.lumiByLbn.find(lb)==run.lumiByLbn.end()) return 0.;
165  return run.lumiByLbn[lb].first*1E-6;
166 }
167 
169  if(!m_isInitialized) {
170  if (m_printInfo) {
171  Info("GetLumiBlockMu","Initializing the subtool..");
172  }
173  Initialize();
174  }
175  if(m_runs.find(runNumber)==m_runs.end()) return -1.0;
176  auto& run = m_runs[runNumber];
177 
178  if(run.lumiByLbn.find(lb)==run.lumiByLbn.end()) return -1.0;
179  return run.lumiByLbn[lb].second;
180 }
181 
182 Double_t CP::TPileupReweighting::GetIntegratedLumiFraction(Int_t periodNumber, UInt_t start, UInt_t end) {
183  if(!m_isInitialized) {
184  if (m_printInfo) {
185  Info("GetIntegratedLumiFraction","Initializing the subtool..");
186  }
187  Initialize();
188  }
189 
190  if(!m_lumiVectorIsLoaded) {
191  Error("GetIntegratedLumiFraction","No UNPRESCALED (Trigger=None) Lumicalc file loaded, so no lumi fraction possible, returning 0");
192  return 0;
193  }
194  //total lumi in period is given by sumOfWeights[-1]
195  double total = m_periods[periodNumber]->sumOfWeights[-1];
196  //go through associated runs, adding up lumi from all within start and end (inclusive)
197  double numer(0);
198  for(auto run : m_periods[periodNumber]->runNumbers) {
199  if(run >= start && run <= end) numer += m_runs[run].inputHists["None"]->GetSumOfWeights();
200  }
201  return numer / total;
202 }
203 
204 Double_t CP::TPileupReweighting::GetIntegratedLumiFraction(Int_t periodNumber, Double_t mu, UInt_t start, UInt_t end) {
205  if(!m_isInitialized) {
206  if (m_printInfo) {
207  Info("GetIntegratedLumiFraction","Initializing the subtool..");
208  }
209  Initialize();
210  }
211  if(!m_lumiVectorIsLoaded) {
212  Error("GetIntegratedLumiFraction","No UNPRESCALED (Trigger=None) Lumicalc file loaded, so no lumi fraction possible, returning 0");
213  return 0;
214  }
215  //determine the mubin
216  if(!m_emptyHistogram) { Error("GetIntegratedLumiFraction", "Cannot do this without a lumicalc file!"); return 0; }
217  Int_t muBin = m_emptyHistogram->FindFixBin(mu);
218  //the total mu in this bin in this period is given by the triggerHist = "None" hist
219 
220 
221 
222  double total = m_triggerObjs["None"]->triggerHists[periodNumber][0]->GetBinContent(muBin);
223  //loop over assigned runs, if in range then include in numerator
224  double numer = 0;
225  for(auto run : m_periods[periodNumber]->runNumbers) {
226  if(run >= start && run <= end) numer += m_runs[run].muDist->GetBinContent(muBin);
227  }
228 
229  return numer / total;
230 }
231 
232 
233 
234 Int_t CP::TPileupReweighting::UsePeriodConfig(const TString& configName) {
235  if(configName=="MC11a") {
236  AddPeriod(180164, 177986,180481); //associates mc runnumber 180164 with data period 177986 to 180481 (period B-D)
237  AddPeriod(183003, 180614,184169); //period E-H
238  AddPeriod(185649, 185353,186934); //period I-K1. For I-K you would change the last number to 187815
239  AddPeriod(185761, 186935,191933); //everything else. Thanks Ellie!
240  SetUniformBinning(100,0,50);
241  if (m_printInfo) {
242  Info("UsePeriodConfig","Using MC11a Period configuration");
243  }
244  return 0;
245  } else if(configName=="MC11b" || configName=="MC11c") {
246  AddPeriod(180164, 177986, 180481);
247  AddPeriod(183003, 180614, 184169);
248  AddPeriod(186169, 185353, 187815);
249  AddPeriod(189751, 188902, 191933);
250  SetUniformBinning(100,0,50);
251  if (m_printInfo) {
252  Info("UsePeriodConfig","Using MC11b/c Period configuration");
253  }
254  return 0;
255  } else if(configName=="MC12a") {
256  AddPeriod(195847,200804,216432);
257  //mc12a binning is in integer values of mu
258 
259  if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
260  Error("UsePeriodConfig","Cannot use MC12a, an incompatible config has already been set up");
261  throw std::runtime_error("Throwing 13: Cannot use MC14, an incompatible config has already been set up");
262  return -1;
263  }
264 
265  SetUniformBinning(50,-0.5,49.5);
266  if (m_printInfo) {
267  Info("UsePeriodConfig","Using MC12a Period configuration");
268  }
269  return 0;
270  } else if(configName=="MC12b") {
271  AddPeriod(195848,200804,216432);
272  //mc12b binning is in integer values of mu
273 
274  if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
275  Error("UsePeriodConfig","Cannot use MC12b, an incompatible config has already been set up");
276  throw std::runtime_error("Throwing 13: Cannot use MC14, an incompatible config has already been set up");
277  return -1;
278  }
279 
280  SetUniformBinning(50,-0.5,49.5);
281  if (m_printInfo) {
282  Info("UsePeriodConfig","Using MC12b Period configuration");
283  }
284  return 0;
285  } else if(configName=="MC12ab") {
286  AddPeriod(195847,200804,216432);
287  AddPeriod(195848,200804,216432);
288 
289  if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
290  Error("UsePeriodConfig","Cannot use MC12ab, an incompatible config has already been set up");
291  throw std::runtime_error("Throwing 13: Cannot use MC14, an incompatible config has already been set up");
292  return -1;
293  }
294 
295  //mc12a/b binning is in integer values of mu
296  SetUniformBinning(50,-0.5,49.5);
297  if (m_printInfo) {
298  Info("UsePeriodConfig","Using MC12ab Period configuration");
299  }
300  return 0;
301  } else if(configName=="MC14_8TeV") {
302  AddPeriod(212272,200804,216432);
303  if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=50 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01 || fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(50)-49.5)>0.01) ) {
304  Error("UsePeriodConfig","Cannot use MC14_8TeV, an incompatible config has already been set up");
305  throw std::runtime_error("Throwing 13: Cannot use MC14_8TeV, an incompatible config has already been set up");
306  return -1;
307  }
308  SetUniformBinning(50,-0.5,49.5);
309  if (m_printInfo) {
310  Info("UsePeriodConfig","Using MC14_8TeV Period configuration");
311  }
312  return 0;
313  } else if(configName=="MC14_13TeV") {
314  AddPeriod(222222,222222,999999);
315  if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=100 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1)+0.5)>0.01|| fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(100)-99.5)>0.01) ) {
316  Error("UsePeriodConfig","Cannot use MC14_13TeV, an incompatible config has already been set up");
317  throw std::runtime_error("Throwing 13: Cannot use MC14_13TeV, an incompatible config has already been set up");
318  return -1;
319  }
320  SetUniformBinning(100,-0.5,99.5);
321  if (m_printInfo) {
322  Info("UsePeriodConfig","Using MC14_13TeV Period configuration");
323  }
324  return 0;
325  } else if(configName=="MC15") {
326  AddPeriod(222510,222222,999999);
327  AddPeriod(222525,222222,999999);
328  AddPeriod(222526,222222,999999);
329  AddPeriod(284500,222222,999999);
330  AddPeriod(295000,222222,999999); //mc15c for pPb data in 2016
331  if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=100 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1))>0.01 || fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(100)-100.)>0.01) ) {
332  Error("UsePeriodConfig","Cannot use MC15, an incompatible config has already been set up");
333  throw std::runtime_error("Throwing 13: Cannot use MC15, an incompatible config has already been set up");
334  return -1;
335  }
336  SetUniformBinning(100,0,100);
337  if (m_printInfo) {
338  Info("UsePeriodConfig","Using MC15 Period configuration");
339  }
340  return 0;
341  } else if(configName=="Run2") {
342  m_autoRunStart = 222222; m_autoRunEnd = 999999; //periods will be automatically added during Fill
343  if(m_emptyHistogram && (m_emptyHistogram->GetNbinsX()!=100 || fabs(m_emptyHistogram->GetXaxis()->GetBinLowEdge(1))>0.01 || fabs(m_emptyHistogram->GetXaxis()->GetBinUpEdge(100)-100.)>0.01) ) {
344  Error("UsePeriodConfig","Cannot use open Run2, an incompatible config has already been set up");
345  throw std::runtime_error("Throwing 13: Cannot use Run2 config, an incompatible config has already been set up");
346  return -1;
347  }
348  SetUniformBinning(100,0,100);
349  if (m_printInfo) {
350  Info("UsePeriodConfig","Using Run2 Period configuration, which assumes period assignment of 222222 to 999999");
351  }
352  return 0;
353  } else if(configName=="MC16") {
354  /* period configs are now assigned through the parent tool for MC16 */
355 
356  SetUniformBinning(100,0,100);
357  if (m_printInfo) {
358  Info("UsePeriodConfig","Using MC16 Period configuration");
359  }
360  return 0;
361  }
362  Error("UsePeriodConfig","Unrecognized period config");
363  return -1;
364 }
365 
366 
368 Int_t CP::TPileupReweighting::AddPeriod(Int_t periodNumber, UInt_t start, UInt_t end) {
369  //you can add more periods only in counting mode, but they must not be with subperiods!
370 
371  if(m_isInitialized && !m_countingMode) {
372  Error("AddPeriod","You cannot AddPeriod after initializing the tool, except when in config file generating mode. Reorder your code!");
373  throw std::runtime_error("Throwing 1: You cannot AddPeriod after initializing the tool, except when in config file generating mode. Reorder your code!");
374  }
375 
376  if(m_periods.find(periodNumber)==m_periods.end()) {
377  m_periodList.emplace_back( periodNumber, start, end,
378  GetDefaultChannel( -1 ) );
379  m_periods[periodNumber] = &( m_periodList.back() );
380  return periodNumber;
381  }
382 
383  //if period exists but both start and end are zero, then just assign that
384  Period* p = m_periods[periodNumber];
385 
386  if(p->subPeriods.size() == 0 && p->start==0 && p->end==0) { p->start=start; p->end=end; return periodNumber; }
387 
388  //check if period already contained (i.e. same period or one of the existing subperiods) ... if it is, do nothing
389  if(p->start==start && p->end==end) return periodNumber;
390  for(auto pp : p->subPeriods) {
391  if(pp->start==start && pp->end==end) return pp->id;
392  }
393 
394  //if get here, and in counting mode, this isn't allowed! .. i.e. we can't have subperiods
395  if(m_isInitialized && m_countingMode) {
396  Error("AddPeriod","You cannot have subperiods when in Config File Generating mode");
397  throw std::runtime_error("Throwing 44: In Config File Generating mode, but detected subperiods in the period definition. This is not allowed.");
398  }
399 
400  //if period is already defined, we create a new period with start and end, and assign our period to this one
401  //also create a new period for the existing period configuration of this period
402 
403  if(p->subPeriods.size()==0) {
404  while(m_periods.find(m_nextPeriodNumber) != m_periods.end()) m_nextPeriodNumber++;
405  m_periodList.emplace_back( m_nextPeriodNumber, p->start, p->end,
406  p->defaultChannel );
407  m_periods[m_nextPeriodNumber] = &( m_periodList.back() );
408  p->subPeriods.push_back(m_periods[m_nextPeriodNumber]);
409  p->start = 0; p->end=0;
410  }
411 
412  while(m_periods.find(m_nextPeriodNumber) != m_periods.end()) m_nextPeriodNumber++;
413  m_periodList.emplace_back( m_nextPeriodNumber, start, end,
414  p->defaultChannel );
415  m_periods[m_nextPeriodNumber] = &( m_periodList.back() );
416  p->subPeriods.push_back(m_periods[m_nextPeriodNumber]);
417 
418  return m_nextPeriodNumber;
419 }
420 
421 
423  //go through periods, first period that includes this runNumber wins!
424  //otherwise it returns the 'global period number': -1
425  for(auto period : m_periods) {
426  if(period.second->id==-1) continue;
427  if(period.second->contains(runNumber)) return period.second->id;
428  }
429  return -1;
430 
431 }
432 
433 
434 std::unique_ptr< TH1 >
436 
437  TString s("pileup");
438  if(channelNumber>=0) {
439  s += "_chan"; s += channelNumber;
440  } else {
441  s += "_data";
442  }
443  s+="_run"; s+= runNumber;
444  //get the empty histogram
445  if(!m_emptyHistogram) {
446  Error("CloneEmptyHistogram","There is no binning info - use SetBinning/SetUniformBinning or load a prw config file (This usually means you need to call AddConfigFile BEFORE AddLumiCalcFile)");
447  throw std::runtime_error("Throwing 47: There is no binning info - use SetBinning/SetUniformBinning or load a prw config file (This usually means you need to call AddConfigFile BEFORE AddLumiCalcFile)");
448  }
449 
450  std::unique_ptr< TH1 > out( dynamic_cast<TH1*>(m_emptyHistogram->Clone(s)) );
451  if (out) {
452  out->SetTitle(s);
453  out->SetDirectory(0); //take ownership
454  out->SetEntries(0);
455  }
456  return out;
457 
458 }
459 
460 Int_t CP::TPileupReweighting::GenerateMetaDataFile(const TString& fileName,const TString& channelBranchName) {
461 
462  TTree inTree("in","in");
463  inTree.ReadFile(fileName);
464  TTree outTree("ChannelMetaData","ChannelMetaData");
465  Int_t chanNum;
466  std::map<TString, Double_t> data;
467  TObjArray *leaves = inTree.GetListOfLeaves();
468  if(leaves==0) {Error("GenerateMetaDataFile","No leaves"); return -1; }
469  for(Int_t i=0;i<leaves->GetEntries();++i) {
470  TLeaf *leaf = (TLeaf *)leaves->At(i);
471  if(leaf) {
472  TBranch *branch = leaf->GetBranch();
473  if(strcmp(branch->GetName(),channelBranchName)==0) {
474  //this is the channel branch
475  if(strcmp(leaf->GetTypeName(),"Int_t")!=0) {
476  Error("GenerateMetaDataFile","Channel Branch must be type Int_t"); return -1;
477  }
478  branch->SetAddress(&chanNum);
479  outTree.Branch(channelBranchName,&chanNum);
480  } else if(strcmp(leaf->GetTypeName(),"Double_t")!=0) {
481  Warning("GenerateMetaDataFile","Cannot read non-double branch: %s",branch->GetName());
482  } else {
483  branch->SetAddress(&(data[branch->GetName()]));
484  outTree.Branch(branch->GetName(),&(data[branch->GetName()]));
485  }
486  }
487  }
488  //loop over tree entries and read in the values according to property type
489  for(Int_t i=0;i<inTree.GetEntries();++i) {
490  inTree.GetEntry(i);
491  outTree.Fill();
492  }
493 
494  //remove the file extension and replace with .root
495  TString outName = fileName(0,fileName.Last('.'));
496  outName += ".prw.root";
497  TFile f1(outName,"RECREATE");
498  outTree.Write();
499  f1.Close();
500  if (m_printInfo) {
501  Info("GenerateMetaDataFile","Succesfully Generated File %s",outName.Data());
502  }
503  return 0;
504 }
505 
506 Int_t CP::TPileupReweighting::AddMetaDataFile(const TString& fileName,const TString& channelBranchName) {
507  TDirectory* origDir = gDirectory;
508  TTree* tmp = 0;
509  TFile* rootFile = 0;
510  if(fileName.EndsWith(".root")) {
511  //looks for ChannelMetaData ttree
512  // Load the data ROOT file
513  rootFile = TFile::Open( fileName, "READ" );
514  if ( rootFile->IsZombie() ) {
515  Error("AddMetaDataFile","Could not open file: %s",fileName.Data());
516  throw std::runtime_error("Throwing 6");
517  } else {
518  //try to get the the known TTrees
519  tmp = (TTree*)rootFile->Get( "ChannelMetaData" );
520  if(!tmp) {
521  Error("AddMetaDataFile","%s is not a valid metadata file. Should have a ChannelMetaData TTree",fileName.Data());
522  throw std::runtime_error("Throwing 7");
523  }
524  rootFile->Close();
525  }
526  delete rootFile;
527  } else {
528  //try to make the TTree by reading the file
529  tmp = new TTree("ChannelMetaData","ChannelMetaData");
530  tmp->ReadFile(fileName);
531  }
532  Int_t chanNum;
533  std::map<TString, Double_t> data;
534  TObjArray *leaves = tmp->GetListOfLeaves();
535  if(leaves==0) {Error("AddMetaDataFile","No leaves"); return -1; }
536  for(Int_t i=0;i<leaves->GetEntries();++i) {
537  TLeaf *leaf = (TLeaf *)leaves->At(i);
538  if(leaf) {
539  TBranch *branch = leaf->GetBranch();
540  if(strcmp(branch->GetName(),channelBranchName)==0) {
541  //this is the channel branch
542  if(strcmp(leaf->GetTypeName(),"Int_t")!=0) {
543  Error("AddMetaDataFile","Channel Branch must be type Int_t");
544  throw std::runtime_error("Throwing 7");
545  }
546  branch->SetAddress(&chanNum);
547  } else if(strcmp(leaf->GetTypeName(),"Double_t")!=0) {
548  Warning("AddMetaDataFile","Cannot read non-double branch: %s",branch->GetName());
549  } else {
550  branch->SetAddress(&(data[branch->GetName()]));
551  }
552  }
553  }
554  //loop over tree entries and read in the values
555  for(Int_t i=0;i<tmp->GetEntries();++i) {
556  tmp->GetEntry(i);
557  for(std::map<TString, Double_t>::iterator it = data.begin();it!=data.end();++it) {
558  if(m_metadata.find(it->first)!=m_metadata.end()&&m_metadata[it->first].find(chanNum)!=m_metadata[it->first].end() && m_metadata[it->first][chanNum]!=it->second) {
559  Warning("AddMetaDataFile", "Overriding metadata [%s,%d]. %f becomes %f",(it->first).Data(),chanNum,m_metadata[it->first][chanNum],it->second);
560  }
561  m_metadata[it->first][chanNum] = it->second;
562  }
563  }
564 
565  delete tmp;
566  if(rootFile) { rootFile->Close();delete rootFile;}
567 
568  // Return to the original ROOT directory
570 
571  return 0;
572 }
573 
575  if(m_metadatatree) delete m_metadatatree;
576  m_metadatatree = new TTree(TString(this->GetName())+"MetaData",TString(this->GetName())+"MetaData");
577  m_metadatatree->SetDirectory(0);
578  Int_t channel=0;
579  m_metadatatree->Branch("mc_channel_number",&channel);
580  std::map<TString,Double_t> data;
581  std::map<Int_t,bool> channels;
582  //create a branch for each metadata item
583  for(std::map<TString,std::map<Int_t,Double_t> >::iterator it=m_metadata.begin();it!=m_metadata.end();++it) {
584  for(std::map<Int_t,Double_t>::iterator it2=(it->second).begin();it2!=(it->second).end();++it2) {
585  channels[it2->first]=true; //records which channels we have metadata for
586  }
587  const auto &[ptr, inserted ] = data.try_emplace(it->first, 0.);
588  if (inserted){
589  m_metadatatree->Branch(it->first,&(data[it->first]));
590  }
591  }
592  //also create branches for the NumberOfEvents and SumOfEventWeights
593  data["NumberOfEvents"]=0.;
594  m_metadatatree->Branch("NumberOfEvents",&(data["NumberOfEvents"]));
595  data["SumOfEventWeights"]=0.;
596  m_metadatatree->Branch("SumOfEventWeights",&(data["SumOfEventWeights"]));
597  //and add to channels list any channels that have events
598  auto& myMap = m_periods[-1]->sumOfWeights;
599  for(std::map<Int_t,Double_t>::iterator it=myMap.begin();it!=myMap.end();++it) {
600  if(it->first>=0 && it->second>0.) channels[it->first]=true;
601  }
602 
603  //now loop over the channels and fill the TTree
604  for(std::map<Int_t,bool>::iterator it=channels.begin();it!=channels.end();++it) {
605  channel=it->first;
606  for(std::map<TString,Double_t>::iterator it2=data.begin();it2!=data.end();++it2) {
607  if(it2->first=="NumberOfEvents") {
608  //check for this in globalNumberOfEntries
609  if(myMap.find(channel)==myMap.end()) {
610  data[it2->first]=0.;
611  Warning("GetMetaDataTree","Channel %d does not have MetaData %s",it->first,(it2->first).Data());
612  } else {
613  data[it2->first]=m_periods[-1]->numberOfEntries[channel];
614  }
615  } else if(it2->first=="SumOfEventWeights") {
616  //check for this in globalTotals
617  auto& myMap2 = m_periods[-1]->sumOfWeights;
618  if(myMap2.find(channel)==myMap2.end()) {
619  data[it2->first]=0.;
620  Warning("GetMetaDataTree","Channel %d does not have MetaData %s",it->first,(it2->first).Data());
621  } else {
622  data[it2->first]=m_periods[-1]->sumOfWeights[channel];
623  }
624  } else {
625  auto& myMap2 = m_metadata[it2->first];
626  if(myMap2.find(channel)==myMap2.end()) {
627  //this channel doesn't have this property
628  data[it2->first]=0.;
629  Warning("GetMetaDataTree","Channel %d does not have MetaData %s",it->first,(it2->first).Data());
630  } else {
631  data[it2->first]=myMap2[channel];
632  }
633  }
634  }
635  m_metadatatree->Fill();
636  }
637 
638  m_metadatatree->BuildIndex("mc_channel_number");
639  m_metadatatree->ResetBranchAddresses();
640 
641  return m_metadatatree;
642 }
643 
644 
645 
646 
647 
649 
650  bool isMC=true;
651  //should have the standard structure
652  Int_t channel = 0; UInt_t runNbr = 0;
653  auto pStarts = std::make_unique< std::vector< UInt_t > >();
654  auto pStartsPtr = pStarts.get();
655  auto pEnds = std::make_unique< std::vector< UInt_t > >();
656  auto pEndsPtr = pEnds.get();
657  Char_t histName[150];
658  Char_t customName[150];
659  Bool_t isDefaultForRunNumber(false); bool hasDefaultsBranch=false;
660  if(strcmp(tree->GetName(),"MCPileupReweighting")==0) {strcpy(customName,"pileup");isMC=true;}
661  else {
662  if(tree->SetBranchAddress("CustomName",&customName)!=0) {
663  Error("AddDistributionTree","Could not find CustomName branch in TTree");throw std::runtime_error("Throwing 18");
664  }
665  }
666  if(strcmp(tree->GetName(),"DataCustomReweighting")==0) {channel=-1;isMC=false;}
667  else {
668  if(tree->SetBranchAddress("Channel",&channel)!=0) {
669  Error("AddDistributionTree","Could not find Channel branch in TTree");throw std::runtime_error("Throwing 18");
670  }
671  }
672  if(tree->SetBranchAddress("RunNumber",&runNbr)!=0) {
673  Error("AddDistributionTree","Could not find RunNumber branch in TTree");throw std::runtime_error("Throwing 18");
674  }
675  if(isMC) {
676  if(tree->SetBranchAddress("PeriodStarts",&pStartsPtr)!=0) {
677  Error("AddDistributionTree","Could not find PeriodStarts branch in TTree");throw std::runtime_error("Throwing 18");
678  }
679  if(tree->SetBranchAddress("PeriodEnds",&pEndsPtr)!=0) {
680  Error("AddDistributionTree","Could not find PeriodEnds branch in TTree");throw std::runtime_error("Throwing 18");
681  }
682  if(tree->FindBranch("IsDefaultForRunNumber")!=0) {
683  tree->SetBranchAddress("IsDefaultForRunNumber",&isDefaultForRunNumber);
684  hasDefaultsBranch=true;
685  }
686  }
687 
688  if(tree->SetBranchAddress("HistName",&histName)!=0) {
689  Error("AddDistributionTree","Could not find HistName branch in TTree");throw std::runtime_error("Throwing 18");
690  }
691  long n = tree->GetEntries();
692  std::map<TString,bool> loadedHistos; //avoid double-loading from this file
693  for(long i=0;i<n;i++) {
694  tree->GetEntry(i);
695  TString sHistName(histName);
696  TString weightName(customName);
697  const auto &[ptr, inserted] = loadedHistos.try_emplace(sHistName, true);
698 
699  if(inserted) {
700  if(( (!m_ignoreFilePeriods) || m_periods.find(runNbr)==m_periods.end()) && isMC) {
701  //if ignoring file periods, will still add the period if it doesnt exist!
702  for(unsigned int j=0;j<pStarts->size();j++) {
703  unsigned int start = pStarts->at(j);
704  unsigned int end = pEnds->at(j);
705  AddPeriod(runNbr,start,end);
706  }
707  }
708  //load the histogram
709  TH1 *histo = (TH1*)file->Get( sHistName );
710  if(!histo) histo = (TH1*)file->Get( TString::Format("%sPileupReweighting/%s",m_prwFilesPathPrefix.c_str(),sHistName.Data()) ); //try from subdir
711  if(!histo) {
712  Error("AddDistributionTree","Unable to find the histogram %s in the File %s",sHistName.Data(),file->GetName());
713  throw std::runtime_error("Throwing 21");
714  }
715  //add it to the histograms
716  AddDistribution(histo, runNbr, channel);
717  //we also add it to the global count, if this isn't data
718  if(channel>=0) AddDistribution(histo,-1,channel);
719  }
720  if(hasDefaultsBranch) {
721  if(isDefaultForRunNumber && m_periods.find(runNbr)!=m_periods.end()) m_periods[runNbr]->defaultChannel=channel;
722  }
723  }
724 
725 }
726 
727 //data is signalled by a negative channel number
728 Int_t CP::TPileupReweighting::AddDistribution(TH1* hist,Int_t runNumber, Int_t channelNumber) {
729 
730  //NOTE: should probably just ignore histograms with 0 entries
731 
732  if(m_isInitialized) {
733  Error("AddDistribution","You cannot AddDistribution after initializing the tool. Reorder your code!");
734  throw std::runtime_error("Throwing 5: You cannot AddLumiCalcFile after initializing the tool. Reorder your code!");
735  }
736 
737  if(channelNumber>=0 && !m_periods[runNumber]) {
738  Error("AddDistribution","Unrecognised periodNumber: %d .. please use AddPeriod to define a period",runNumber);
739  throw std::runtime_error("Throwing 6: Unrecognised periodNumber. Please use AddPeriod to define a period");
740  }
741 
742  if(!m_useMultiPeriods && runNumber>=0) {
743  //not allowing multi periods in a single channel, so before we go any further, require no periods other than this incoming one to have a hist defined for this channel
744  //the global one can though
745  for(auto& period : m_periods) {
746  if(period.first<0) continue; //global sum period
747  if(period.first==runNumber) continue; //ok to add to this
748  if(period.second->inputHists.find(channelNumber)!=period.second->inputHists.end()) {
749  TString myMsg = TString::Format("Attempt to add distribution for channel %d to period %d, but this channels already has period %d defined.\nThis is indicative of use of incompatible PRW config files, please check your config files for multiple periods/runNumbers for the same channel.",channelNumber,runNumber,period.first);
750  Error("AddDistribution","%s",myMsg.Data());
751  throw std::runtime_error(myMsg.Data());
752  }
753  }
754  }
755 
756 
757  auto& inputHist = (channelNumber<0) ? m_runs[runNumber].inputHists["None"] : m_periods[runNumber]->inputHists[channelNumber];
758 
759  if(!inputHist) {
760  //if no emptyHistogram specified, we will use this histogram as the structure for the empty;
761  if(!m_emptyHistogram) {
762  //check if the input histogram is a TH1D or not. We like TH1D, not yucky TH1F
763  if(strcmp(hist->IsA()->GetName(),"TH1D")==0 || strcmp(hist->IsA()->GetName(),"TH2D")==0) {
764  m_emptyHistogram.reset( dynamic_cast<TH1*>(hist->Clone("default")) );
765  } else {
766  //check dimensionality
767  if(hist->GetDimension()==1) {
768  std::vector<Double_t> binsX;
769  for(int i=0;i<=hist->GetNbinsX();i++) binsX.push_back(hist->GetXaxis()->GetBinLowEdge(i+1));
770  TH1D tmpHist("tmpHist","tmpHist",binsX.size()-1,&binsX[0]);
771  m_emptyHistogram.reset( dynamic_cast<TH1*>(tmpHist.Clone("default")) );
772  } else if(hist->GetDimension()==2) {
773  std::vector<Double_t> binsX;std::vector<Double_t> binsY;
774  for(int i=0;i<=hist->GetNbinsX();i++) binsX.push_back(hist->GetXaxis()->GetBinLowEdge(i+1));
775  for(int i=0;i<=hist->GetNbinsY();i++) binsY.push_back(hist->GetYaxis()->GetBinLowEdge(i+1));
776  TH2D tmpHist("tmpHist","tmpHist",binsX.size()-1,&binsX[0],binsY.size()-1,&binsY[0]);
777  m_emptyHistogram.reset( dynamic_cast<TH1*>(tmpHist.Clone("default")) );
778  } else {
779  Error("AddDistribution","Unknown input histogram dimensionality: %d",hist->GetDimension());
780  throw std::runtime_error("Throwing 98");
781  }
782  }
783  if (m_emptyHistogram)
784  {
785  m_emptyHistogram->Reset();m_emptyHistogram->SetEntries(0);m_emptyHistogram->SetDirectory(0);//clear the histogram
786  }
787  }
788  inputHist = CloneEmptyHistogram(runNumber,channelNumber);
789 
790  }
791 
792 
793  //iterator over bins of histogram, filling the TH1 stored in the data map
794  Double_t numEntries = inputHist->GetEntries();
795 
796  std::unique_ptr<TH1> tmpHist;
797  if(channelNumber<0) {
798  m_runs[runNumber].nominalFromHists = true;
799  //for data we will do an interpolation to build the shape and then scale it to the integral ...
800  tmpHist.reset( static_cast<TH1*>(inputHist->Clone("tmpHist")) );
801  tmpHist->Reset();
802  Int_t bin,binx,biny;
803  for(biny=1; biny<=tmpHist->GetNbinsY(); biny++) {
804  for(binx=1; binx<=tmpHist->GetNbinsX(); binx++) {
805  bin = tmpHist->GetBin(binx,biny);
806  Double_t x = tmpHist->GetXaxis()->GetBinCenter(binx)/m_dataScaleFactorX;
807  Double_t y = tmpHist->GetYaxis()->GetBinCenter(biny)/m_dataScaleFactorY;
808  if(tmpHist->GetDimension()==1){
809  tmpHist->SetBinContent(bin, hist->Interpolate(x));
810  } else {
811  tmpHist->SetBinContent(bin, hist->Interpolate(x,y));
812  }
813  }
814  }
815  tmpHist->Scale( hist->Integral() / tmpHist->Integral() );
816  tmpHist->SetEntries( hist->GetEntries() );
817  hist = tmpHist.get();
818  }
819 
820  Int_t bin,binx,biny;
821  for(biny=1; biny<=hist->GetNbinsY(); biny++) {
822  for(binx=1; binx<=hist->GetNbinsX(); binx++) {
823  bin = hist->GetBin(binx,biny);
824  Double_t value = hist->GetBinContent(bin);
825  Double_t x = hist->GetXaxis()->GetBinCenter(binx);
826  Double_t y = hist->GetYaxis()->GetBinCenter(biny);
827  //shift x,y,z by the MCScaleFactors as appropriate
828  if(channelNumber>=0) {x *= m_mcScaleFactorX; y *= m_mcScaleFactorY;}
829  //data scale factor now dealt with in above interpolation
830  Int_t inBin = inputHist->FindFixBin(x,y);
831  Double_t inValue = inputHist->GetBinContent(inBin);
832  inputHist->SetBinContent(inBin,inValue+value);
833  }
834  }
835 
836  //also keep track of the number of entries
837  //SetBinContent screws with the entry count, so had to record it before the loops above
838  inputHist->SetEntries(numEntries+hist->GetEntries());
839 
840 
841  m_countingMode=false;
842  return 0;
843 
844 }
845 
846 
847 Int_t CP::TPileupReweighting::AddLumiCalcFile(const TString& fileName, const TString& trigger) {
848 
849  if(m_isInitialized) {
850  Error("AddLumiCalcFile","You cannot AddLumiCalcFile after initializing the tool. Reorder your code!");
851  throw std::runtime_error("Throwing 5: You cannot AddLumiCalcFile after initializing the tool. Reorder your code!");
852  }
853  TDirectory* origDir = gDirectory;
854  // Load the data ROOT file
855  TFile* rootFile = TFile::Open( fileName, "READ" );
856  if ( rootFile->IsZombie() ) {
857  Error("AddConfigFile","Could not open file: %s",fileName.Data());
858  std::string toThrow = "Throwing 6: Could not open file: "; toThrow += fileName.Data();
859  throw std::runtime_error(toThrow);
860  } else {
861  //try to get the the known TTrees
862  TTree *tmp = (TTree*)rootFile->Get( "LumiMetaData" );
863  if(tmp) {
864  m_lumicalcFiles[trigger].push_back(fileName);
865  if(trigger=="None") {
866  if (m_printInfo) {
867  Info("AddLumiCalcFile","Adding LumiMetaData (scale factor=%f)...",m_dataScaleFactorX);
868  }
869  //structure expected is as given by iLumiCalc:
870  // RunNbr, AvergeInteractionPerXing, IntLumi
871  UInt_t runNbr=0;Float_t intLumi=0;UInt_t lbn=0;TBranch *b_runNbr=0;TBranch *b_intLumi=0;TBranch *b_lbn=0;
872  Float_t mu=0.; TBranch *b_mu=0;
873  if(tmp->SetBranchAddress("RunNbr",&runNbr,&b_runNbr)!=0) {
874  Error("AddLumiCalcFile","Could not find RunNbr branch in Data TTree");throw std::runtime_error("Could not find RunNbr branch in Data TTree");
875  }
876  if(tmp->SetBranchAddress("AvergeInteractionPerXing",&mu,&b_mu)!=0) {
877  Error("AddLumiCalcFile","Could not find AvergeInteractionPerXing branch in Data TTree");throw std::runtime_error("Could not find AvergeInteractionPerXing branch in Data TTree");
878  }
879  if(tmp->SetBranchAddress("IntLumi",&intLumi,&b_intLumi)!=0) {
880  Error("AddLumiCalcFile","Could not find IntLumi branch in Data TTree");throw std::runtime_error("Could not find IntLumi branch in Data TTree");
881  }
882  if(tmp->SetBranchAddress("LBStart",&lbn,&b_lbn)!=0) {
883  Error("AddLumiCalcFile","Could not find LBStart branch in Data TTree");throw std::runtime_error("Could not find LBStart branch in Data TTree");
884  }
885  long nEntries = tmp->GetEntries();
886 
887  for(long i=0;i<nEntries;i++) {
888  b_runNbr->GetEntry(i);b_intLumi->GetEntry(i);b_mu->GetEntry(i);
889  runNbr += m_lumicalcRunNumberOffset;
890  //save the lumi by run and lbn
891  b_lbn->GetEntry(i);
892 
893  //before going any further, check the runnbr and lbn are ok
894  if(!m_parentTool->runLbnOK(runNbr,lbn)) continue;
895 
896  Run& r = m_runs[runNbr];
897  if(trigger=="None") {
898  r.lumiByLbn[lbn].first += intLumi;
899  r.lumiByLbn[lbn].second = mu;
900  if(r.nominalFromHists) continue; //don't fill runs that we already filled from hists ... only happens for the 'None' trigger
901  }
902 
903  //rescale the mu value ... do this *after* having stored the mu value in the lumiByLbn map
904  mu *= m_dataScaleFactorX;
905  //fill into input data histograms
906  //check if we need to create an empty histogram
907 
908  std::unique_ptr<TH1>& histptr = r.inputHists[trigger];
909  if (!histptr) {
910  histptr = CloneEmptyHistogram(runNbr,-1);
911  }
912  histptr->Fill(mu,intLumi);
913  }
914  m_countingMode=false;
915  m_lumiVectorIsLoaded=true;
916  } else {
917  if (m_printInfo) {
918  Info("AddLumiCalcFile","Adding LumiMetaData for DataWeight (trigger=%s) (scale factor=%f)...",trigger.Data(),m_dataScaleFactorX);
919  }
920  }
921 
922  } else {
923  Error("AddLumiCalcFile","No LumiMetaData found in file %s. not a LumiCalcFile?", fileName.Data());
924  throw std::runtime_error("No LumiMetaData found in file, not a LumiCalcFile?");
925  }
926  }
927 
928  delete rootFile;
929 
930  // Return to the original ROOT directory
932 
933  return 0;
934 }
935 
937 
938  //open the file and look for config TTrees
939  //recognised TTrees are: ChannelMetaData, MCPileupReweighting, MCCustomReweighting, DataCustomReweighting
940 
941  if(m_isInitialized) {
942  Error("AddConfigFile","You cannot AddConfigFile after initializing the tool. Reorder your code!");
943  throw std::runtime_error("Throwing 5: You cannot AddConfigFile after initializing the tool. Reorder your code!");
944  }
945 
946  TDirectory* origDir = gDirectory;
947  // Load the data ROOT file
948  TFile* rootFile = TFile::Open( fileName, "READ" );
949  if ( rootFile->IsZombie() ) {
950  Error("AddConfigFile","Could not open file: %s",fileName.Data());
951  std::string toThrow = "Throwing 6: Could not open file: "; toThrow += fileName.Data();
952  throw std::runtime_error(toThrow);
953  } else {
954  //try to get the the known TTrees
955  TTree *tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"MCPileupReweighting").c_str() );
956  if(tmp) {
957  //Info("AddConfigFile","Adding MCPileupReweighting...");
958  AddDistributionTree(tmp,rootFile);
959  m_countingMode=false;
960  }
961  tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"MCCustomReweighting").c_str() );
962  if(tmp) {
963  //Info("AddConfigFile","Adding MCCustomReweighting...");
964  AddDistributionTree(tmp,rootFile);
965  m_countingMode=false;
966  }
967  tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"DataCustomReweighting").c_str() );
968  if(tmp) {
969  //Info("AddConfigFile","Adding DataCustomReweighting...");
970  AddDistributionTree(tmp,rootFile);
971  m_countingMode=false;
972  }
973  tmp=0; tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"PileupReweighting/MCPileupReweighting").c_str() ); //allow trees in the PileupReweighting dir instead
974  if(tmp) {
975  //Info("AddConfigFile","Adding MCPileupReweighting...");
976  AddDistributionTree(tmp,rootFile);
977  m_countingMode=false;
978  }
979  tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"PileupReweighting/MCCustomReweighting").c_str() );
980  if(tmp) {
981  //Info("AddConfigFile","Adding MCCustomReweighting...");
982  AddDistributionTree(tmp,rootFile);
983  m_countingMode=false;
984  }
985  tmp = 0;tmp = (TTree*)rootFile->Get( std::string(m_prwFilesPathPrefix+"PileupReweighting/DataCustomReweighting").c_str() );
986  if(tmp) {
987  //Info("AddConfigFile","Adding DataCustomReweighting...");
988  AddDistributionTree(tmp,rootFile);
989  m_countingMode=false;
990  }
991  }
992 
993  delete rootFile;
994 
995  // Return to the original ROOT directory
997 
998  return 0;
999 }
1000 
1002  if(m_isInitialized) {
1003  Error("RemoveChannel","You cannot RemoveChannel after initializing the tool. Reorder your code!");
1004  throw std::runtime_error("Throwing 5: You cannot RemoveChannel after initializing the tool. Reorder your code!");
1005  }
1006 
1007  bool found(false);
1008  for(auto& period : m_periods) {
1009  auto itr = period.second->inputHists.begin();
1010  while(itr != period.second->inputHists.end()) {
1011  if(itr->first!=chanNum) ++itr;
1012  else {
1013  found=true;
1014  itr = period.second->inputHists.erase(itr);
1015  }
1016  }
1017  }
1018  return found;
1019 }
1020 
1022 
1023  //Need to calculate a period weight for each period
1024  //Need a reweighting histogram for each period
1025  //for merged mc run numbers, we must generate or modify the existing mc distributions
1026 
1027  //print out the period assignments if in debug mode
1028  if(m_debugging) {
1029  PrintPeriods();
1030  }
1031 
1032  if(m_countingMode) {
1033  if (m_printInfo) {
1034  Info("Initialize","In Config File Generating mode. Remember to call WriteToFile!");
1035  }
1036  //need to check no periods have subperiods, this is not allowed in counting mode
1037  for(auto period : m_periods) {
1038  if(period.second->subPeriods.size()!=0) {
1039  Error("Initialize","You cannot have subperiods when in Config File Generating mode");
1040  throw std::runtime_error("Throwing 44: In Config File Generating mode, but detected subperiods in the period definition. This is not allowed.");
1041  }
1042  }
1043  //histograms are instantiated in the event loop as the channels occur
1044  //can delay period definition here
1045  m_isInitialized=true;
1046  return 0;
1047  }
1048 
1049  //find all channels that have too much unrepresented data (more than tolerance). We will remove these channels entirely
1050 
1051  std::map<int,std::map<int,bool>> trackWarnings; //map used to avoid reprinting warnings as loop over histogram bins
1052 
1053  double totalData(0);
1054  //count up the unrepresented data, indexed by channel
1055  for(auto& run : m_runs) {
1056  if(run.second.inputHists.find("None") == run.second.inputHists.end()) continue;
1057  TH1* hist = run.second.inputHists["None"].get();
1058  Int_t bin,binx,biny;
1059  for(biny=1; biny<=hist->GetNbinsY(); biny++) {
1060  for(binx=1; binx<=hist->GetNbinsX(); binx++) {
1061  bin = hist->GetBin(binx,biny);
1062  Double_t value = hist->GetBinContent(bin);
1063  totalData += value;
1064  if(value==0.) continue;
1065  bool isUnrep(false);
1066  std::map<Int_t, bool> doneChannels;
1067  for(auto& period : m_periods) {
1068  if(period.first != period.second->id) continue; //skip remappings
1069  if(period.first == -1) continue; //don't look at the global period when calculating unrepresented data
1070  if(!period.second->contains(run.first)) continue;
1071  for(auto& inHist : m_periods[-1]->inputHists) { //use global period to iterate over channels
1072  if(inHist.first<0) continue; //skips any data hists (shouldn't happen really)
1073  if(doneChannels[inHist.first]) continue; //dont doublecount the data if the channel was to blame across multiple period
1074  if(period.second->inputHists.find(inHist.first)==period.second->inputHists.end()) {
1075  //all the data is missing ..
1076  isUnrep=true;doneChannels[inHist.first]=true;
1077  unrepDataByChannel[inHist.first] += value;
1078  if(!trackWarnings[inHist.first][period.first]) {
1079  trackWarnings[inHist.first][period.first]=true;
1080  Warning("Initialize","Channel %d has no configuration for period %d -- this is symptomatic of you trying to reweight an MC campaign to data not intended for that campaign (e.g. MC16a with 2017 data)",inHist.first,period.first);
1081  }
1082  } else if(period.second->inputHists[inHist.first]->GetBinContent(bin)==0) {
1083  //hist for this channel in this period exists, but this bin is empty
1084  isUnrep=true;doneChannels[inHist.first]=true;
1085  unrepDataByChannel[inHist.first] += value;
1086  }
1087  }
1088  }
1089  if(isUnrep) unrepDataByChannel[-1] += value; //store total unrep data
1090  }
1091  }
1092  }
1093 
1094  if(m_ignoreBadChannels && unrepDataByChannel[-1] && totalData && (unrepDataByChannel[-1]/totalData)>m_unrepDataTolerance) {
1095  Warning("Initialize","There is %f%% unrepresented data and 'IgnoreBadChannels' property is set to true. Will start ignoring channels until this is below the tolerance (%f%%)", 100.*(unrepDataByChannel[-1]/totalData),100.*m_unrepDataTolerance);
1096  //remove channels one-by-one until we are under tolerance
1097  while( (unrepDataByChannel[-1]/totalData)>m_unrepDataTolerance ) {
1098  int worstChannel = -1;
1099  double worstFraction = 0;
1100  for(auto& channel : unrepDataByChannel) {
1101  if(channel.first<0) continue; //ignore the data
1102  if(channel.second/totalData > worstFraction) { worstChannel=channel.first; worstFraction = channel.second/totalData; }
1103  }
1104  if(worstChannel==-1) {
1105  Error("Initialize","Run out of channels to remove. All your channels are bad ... oh dear oh dear oh dear. Please a better PRW config file");
1106  throw std::runtime_error("Throwing exception 333: Run out of channels to remove. All your channels are bad ... oh dear oh dear oh dear. Please a better PRW config file");
1107  }
1108 
1109  Warning("Initialize","Ignoring channel %d, which has %f%% unrepresented data (%f%% of the total unrep data)",worstChannel,100.*unrepDataByChannel[worstChannel]/totalData,100.*unrepDataByChannel[worstChannel]/unrepDataByChannel[-1]);
1110  //remove the bad channel
1111  RemoveChannel(worstChannel);
1112  unrepDataByChannel.erase(unrepDataByChannel.find(worstChannel));
1113 
1114  //recalculate the total unrep data
1115  unrepDataByChannel[-1] = 0;
1116  for(auto& run : m_runs) {
1117  if(run.second.inputHists.find("None") == run.second.inputHists.end()) continue;
1118  TH1* hist = run.second.inputHists["None"].get();
1119  Int_t bin,binx,biny;
1120  for(biny=1; biny<=hist->GetNbinsY(); biny++) {
1121  for(binx=1; binx<=hist->GetNbinsX(); binx++) {
1122  bin = hist->GetBin(binx,biny);
1123  Double_t value = hist->GetBinContent(bin);
1124  if(value==0.) continue;
1125  for(auto& period : m_periods) {
1126  if(period.first != period.second->id) continue; //skip remappings
1127  if(!period.second->contains(run.first)) continue;
1128  bool done(false);
1129  for(auto& inHist : period.second->inputHists) {
1130  if(inHist.first<0) continue; //skips any data hists (shouldn't happen really)
1131  if((inHist.second)->GetBinContent(bin)==0) {
1132  unrepDataByChannel[-1] += value; done=true; break;
1133  }
1134  }
1135  if(done) break;
1136  }
1137  }
1138  }
1139  }
1140  }
1141  }
1142 
1143 
1144  //need a CompositeTrigger to hold the unprescaled lumi hists ...
1145  m_triggerObjs["None"] = std::make_unique< CompositeTrigger >();
1146 
1147  //go through all periods, use inputHists to populate the primary and secondary dists, and fill sumOfWeights and numberOfEntries
1148  for(auto period : m_periods) {
1149  if(period.first != period.second->id) continue; //skips remappings
1150  for(auto& inputHist : period.second->inputHists) {
1151  int channel = inputHist.first;
1152  TH1* hist = inputHist.second.get();
1153  period.second->sumOfWeights[channel] += hist->GetSumOfWeights();
1154  period.second->numberOfEntries[channel] += hist->GetEntries();
1155  if(hist->GetDimension()==1) {
1156  period.second->primaryHists[channel] = CloneEmptyHistogram(period.first,channel);
1157  period.second->primaryHists[channel]->Add(hist);
1158  }
1159  }
1160  //also add all runs that belong to us ... only do the "None" trigger
1161  for(auto& run : m_runs) {
1162  if(!period.second->contains(run.first)) continue;
1163  if(run.second.inputHists.find("None") == run.second.inputHists.end()) continue;
1164  //add to list of runNumbers for this period
1165  period.second->runNumbers.push_back(run.first);
1166  TH1* hist = run.second.inputHists["None"].get();
1167 
1168  //sweep over bins, checking for unrepresented data and acting accordingly
1169  Int_t bin,binx,biny;
1170  for(biny=1; biny<=hist->GetNbinsY(); biny++) {
1171  for(binx=1; binx<=hist->GetNbinsX(); binx++) {
1172  bin = hist->GetBin(binx,biny);
1173  Double_t value = hist->GetBinContent(bin);
1174  if(value==0.) continue;
1175  //loop over channels, check for zero in this bin ... only need to do if we haven't checked ourselves
1176  if(period.second->inputBinRedirect.find(bin) == period.second->inputBinRedirect.end()) {
1177  for(auto& inputHist : period.second->inputHists) {
1178  if(inputHist.first<0) continue; //skips data
1179  Double_t mcValue = (inputHist.second)->GetBinContent(bin);
1180  if(mcValue==0.) {
1181  if(m_unrepresentedDataAction!=0 && m_debugging && m_printInfo) {
1182  Info("Initialize","Unrepresented data at coords [%f,%f] caused by periodNumber %d in channel %d",hist->GetXaxis()->GetBinCenter(binx),hist->GetYaxis()->GetBinCenter(biny),period.first,inputHist.first);
1183  }
1184  //if we are doing unrepaction=2, just set to 'not the bin number'. if we are doing action=3, find the nearest good bin
1185  if(m_unrepresentedDataAction==0) {period.second->inputBinRedirect[bin] = bin-1;Error("Initialize","Unrepresented data at coords [%f,%f] caused by periodNumber %d in channel %d",hist->GetXaxis()->GetBinCenter(binx),hist->GetYaxis()->GetBinCenter(biny),period.first,inputHist.first);}
1186  if(m_unrepresentedDataAction==1||m_unrepresentedDataAction==2) period.second->inputBinRedirect[bin] = bin-1;
1187  else if(m_unrepresentedDataAction==3) period.second->inputBinRedirect[bin] = GetNearestGoodBin(period.first, bin);
1188  break;
1189  }
1190  }
1191  if(period.second->inputBinRedirect.find(bin) == period.second->inputBinRedirect.end()) period.second->inputBinRedirect[bin] = bin; //a good bin
1192  }
1193  //if is a bad bin, act accordingly
1194  if(period.second->inputBinRedirect[bin] != bin) {
1195  run.second.badBins[bin]=value;
1196  if(m_unrepresentedDataAction==1) {
1197  //remove this data. we have to remember to do this for the trigger data too if this is the pileup weight
1198  hist->SetBinContent(bin,0);
1199  //also modify the m_inputTriggerHistograms ...
1200  for(auto& triggerDists : run.second.inputHists) triggerDists.second->SetBinContent(bin,0);
1201  } else if(m_unrepresentedDataAction==3) {
1202  //reassign the content to that bin
1203  hist->SetBinContent(period.second->inputBinRedirect[bin],hist->GetBinContent(period.second->inputBinRedirect[bin])+hist->GetBinContent(bin));
1204  hist->SetBinContent(bin,0);
1205  //also modify the m_inputTriggerHistograms ...
1206  for(auto& triggerDists : run.second.inputHists) {
1207  triggerDists.second->SetBinContent(period.second->inputBinRedirect[bin],triggerDists.second->GetBinContent(period.second->inputBinRedirect[bin])+triggerDists.second->GetBinContent(bin));
1208  triggerDists.second->SetBinContent(bin,0);
1209  }
1210  }
1211  }
1212  }
1213  }
1214 
1215  //ensure hist is normalized correctly, if was read in from the 'actual mu' config files (where the normalization can be slightly wrong)
1216  if(run.second.nominalFromHists) {
1217  //get the total lumi for the run ...
1218  double totLumi(0);
1219  for(auto lb : run.second.lumiByLbn) {
1220  totLumi += lb.second.first;
1221  }
1222  const double integral = hist->Integral();
1223  if (std::abs(integral) > 0.) {
1224  hist->Scale(totLumi / integral);
1225  }
1226 
1227  }
1228 
1229  //create the period's 'data' hist if necessary
1230  if( hist->GetDimension()==1 ) {
1231  if(!period.second->primaryHists[-1] ) {
1232  period.second->primaryHists[-1] = CloneEmptyHistogram(period.first,-1);
1233 
1234  //this will remain unnormalized, unlike the 'primaryHist'
1235  m_triggerObjs["None"]->triggerHists[period.second->id][0] =
1236  CloneEmptyHistogram( period. first, -1 );
1237  }
1238  }
1239  else if( hist->GetDimension()==2 ) {
1240  if(!period.second->secondaryHists[-1] ) {
1241  auto h = CloneEmptyHistogram(period.first,-1);
1242  period.second->secondaryHists[-1].reset( dynamic_cast< TH2* >( h.release() ) );
1243  period.second->primaryHists[-1].reset(
1244  period.second->secondaryHists[-1]->ProjectionX() );
1245  period.second->primaryHists[-1]->SetDirectory(0);
1246  period.second->primaryHists[-1]->Reset();
1247  m_triggerObjs["None"]->triggerHists[period.second->id][0].reset(
1248  static_cast<TH1*>(period.second->primaryHists[-1]->Clone("triggerHist")) );
1249  m_triggerObjs["None"]->triggerHists[period.second->id][0]->SetDirectory(0);
1250  }
1251 
1252  }
1253 
1254  }
1255 
1256  }
1257 
1258 
1259  //now that all periods and runs have been checked, with redirects set up where necessary, we actually accumulate the data across the runs
1260  //we should also check if there are any data runs that are not covered by the period assignments
1261 
1262  double ignoredData(0);
1263  for(auto& run : m_runs) {
1264  bool covered(false);
1265  TH1* hist = run.second.inputHists["None"].get();
1266  for(auto& period : m_periods) {
1267  if(!period.second->contains(run.first)) continue;
1268  if(period.first!=-1) covered=true; //don't count the global period
1269  if(hist->GetDimension()==1) {
1270  period.second->primaryHists[-1]->Add(hist);
1271  m_triggerObjs["None"]->triggerHists[period.second->id][0]->Add(hist);
1272  } else if(hist->GetDimension()==2) {
1273  period.second->secondaryHists[-1]->Add(hist);
1274  TH1* proj = ((TH2*)hist)->ProjectionX();
1275  period.second->primaryHists[-1]->Add(proj);
1276  m_triggerObjs["None"]->triggerHists[period.second->id][0]->Add(proj);
1277  delete proj;
1278  }
1279  period.second->sumOfWeights[-1] += hist->GetSumOfWeights();
1280  period.second->numberOfEntries[-1] += hist->GetEntries();
1281  }
1282  if(!covered && m_periods.size()>1) {
1283  Warning("Initialize","loaded data in run %d that is not covered by period assignments",run.first);
1284  ignoredData += hist->GetSumOfWeights();
1285  }
1286  }
1287  //double totalData = (m_unrepresentedDataAction==1) ? (m_periods[-1]->sumOfWeights[-1]+unrepDataByChannel[-1]) : m_periods[-1]->sumOfWeights[-1];
1288 
1289  if(ignoredData>0.) Warning("Initialize", "Period Assignments missed %f%% data",100.*ignoredData/totalData);
1290 
1291 
1292 
1293  if(unrepDataByChannel[-1]) {
1294  double frac = unrepDataByChannel[-1]/totalData;
1295  if( frac > m_unrepDataTolerance) {
1296  Error("Initialize", "%f%% unrepresented data, which suggests something is wrong with your prw config. Try EnableDebugging(true) to investigate",100.* (unrepDataByChannel[-1]/totalData));
1297  }
1298  if(m_unrepresentedDataAction==1) {
1299  Warning("Initialize","has %f%% unrepresented data. This was removed (UnrepresentedDataAction=1)",100.*frac);
1300  } else if(m_unrepresentedDataAction==2) {
1301  Warning("Initialize","has %f%% unrepresented data. This was kept in (UnrepresentedDataAction=2)",100.*frac);
1302  } else if(m_unrepresentedDataAction==3) {
1303  Warning("Initialize","has %f%% unrepresented data. This was reassigned (UnrepresentedDataAction=3)",100.*frac);
1304  } else if(m_unrepresentedDataAction==0) {
1305  Error("Initialize","has %f%% unrepresented data:",100.*frac);
1306  //print the report of which channels caused it
1307  for(auto& it : unrepDataByChannel) {
1308  if(it.first<0) continue; //ignore data
1309  Error("Initialize"," Channel %d caused %f%% of the unrepresented data (nEntries=%d)",it.first,100.*it.second/unrepDataByChannel[-1],m_periods[-1]->numberOfEntries[it.first]);
1310  }
1311  Error("Initialize","Exiting. You must decide how to proceed...");
1312  Error("Initialize","1) use AddPeriod or RemapPeriod to redefine the mc periods to include this data. You should not need to regenerate your mc config file");
1313  Error("Initialize","2) use SetUnrepresentedDataAction(1) to remove this data from the weight calculations. You should also veto such data events (using IsUnrepresentedData(..,..) method)");
1314  Error("Initialize","3) use SetUnrepresentedDataAction(2) to leave this data in the calculation. I hope you know what you're doing!!");
1315  Error("Initialize","4) use SetUnrepresentedDataAction(3) to reassign the data to the nearest representable bin");
1316  throw std::runtime_error("Throwing exception 22: Unrepresented data exists. Make a choice how to handle this. See err log for details");
1317  }
1318  //if get here, through the exception if frac is bad
1319  if( frac > m_unrepDataTolerance) {
1320  throw std::runtime_error("Throwing exception 222: Some channel had too much unrepresented data. You should fix your prw file");
1321  }
1322  }
1323 
1324 
1325  if(m_debugging && m_printInfo) Info("Initialize","Normalizing histograms and cleaning up...");
1326  //now that all the distributions are built. Normalize them all
1327  for(auto period : m_periods) {
1328  if(period.first != period.second->id) continue;
1329  for(auto& pHist : period.second->primaryHists) {
1330  normalizeHistogram(pHist.second.get());
1331  }
1332  for(auto& pHist : period.second->secondaryHists) {
1333  normalizeHistogram(pHist.second.get());
1334  }
1335  period.second->inputHists.clear();
1336  }
1337 
1338  //we keep the inputHists for data, because can use that to do random run numbers based on mu
1339 
1340  if(m_debugging && m_printInfo) Info("Initialize","...Done");
1341  //if no input histograms were added, we are in counting mode
1342 
1343 
1344  m_isInitialized=true;
1345 
1346  return 0;
1347 }
1348 
1349 //loop over the mc histograms, find bin closest to given bin that is "ok" (has all nonzero)
1350 //0 = ok bin, 1 = bad bin, 2 = out of range bin
1351 Int_t CP::TPileupReweighting::IsBadBin(Int_t periodNumber, Int_t bin) {
1352 
1353  if(bin<0) return 2; //out of range
1354 
1355  Period* p = m_periods[periodNumber];
1356  for(auto& inHist : p->inputHists) {
1357  int channelNumber = inHist.first;
1358  if(channelNumber<0) continue; //skip data hists
1359  if(bin > (inHist.second->GetNbinsX()+2)*(inHist.second->GetNbinsY()+2)*(inHist.second->GetNbinsZ()+2)) return 2; //definitely out of range by this point!
1360  if(inHist.second->GetBinContent(bin)==0) return 1;
1361  }
1362  return 0;
1363 
1364 }
1365 
1366 //loop over the mc histograms, find bin closest to given bin that is "ok" (has all nonzero)
1367 Int_t CP::TPileupReweighting::GetNearestGoodBin(Int_t thisMCRunNumber, Int_t bin) {
1368 
1369  int binDistance = 1;
1370  do {
1371  int res1 = IsBadBin(thisMCRunNumber,bin+binDistance);
1372  if(!res1) return bin+binDistance;
1373  int res2 = IsBadBin(thisMCRunNumber,bin-binDistance);
1374  if(!res2) return bin-binDistance;
1375  if(res1==2 && res2==2) { //run out of bins
1376  Error("GetNearestGoodBin", "None of the bins are good!!??");
1377  return -1;
1378  }
1379  binDistance++;
1380  } while(true); //scary!
1381 
1382  return -1;
1383 
1384 }
1385 
1387 
1388  //return unrepDataByChannel[channel]/GetSumOfEventWeights(-1); //the denominator is equivalent to asking for the total data!
1389 
1390 
1391  Period* p = m_periods[periodNumber];
1392  if(!p) {
1393  Error("GetUnrepresentedDataFraction","Unrecognised periodNumber: %d",periodNumber);
1394  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1395  }
1396 
1397  //look for channel through sumOfWeights - because if the channel has no unrepData, it will have no entry in the unrepData map
1398  if(p->sumOfWeights.find(channel) == p->sumOfWeights.end()) channel = GetDefaultChannel(periodNumber);
1399 
1400 
1401  if(p->sumOfWeights.find(channel) == p->sumOfWeights.end()) {
1402  Error("GetUnrepresentedDataFraction","Unrecognised channel: %d",channel);
1403  throw std::runtime_error("GetUnrepresentedDataFraction: Unrecognised channel");
1404  }
1405 
1406  return p->unrepData[channel]/p->sumOfWeights[-1];
1407 
1408 }
1409 
1410 
1411 
1413  if(!m_isInitialized) {
1414  if (m_printInfo)Info("GetRandomRunNumber","Initializing the subtool..");
1415  Initialize();
1416  }
1417  if(m_countingMode) { return 0; } //do nothing when in counting mode
1418 
1419  Period* p = m_periods[periodNumber];
1420  if(!p || !p->primaryHists[-1]) {
1421  Error("GetRandomRunNumber","Unrecognised periodNumber or no data loaded for: %d",periodNumber);
1422  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1423  }
1424 
1425  double lumi = p->sumOfWeights[-1] * m_random3->Rndm();
1426 
1427  //go through the runs assigned to this period ...
1428  double lumiSum=0;
1429  for(auto runNum : p->runNumbers) {
1430  lumiSum += m_runs[runNum].inputHists["None"]->GetSumOfWeights();
1431  if(lumiSum >= lumi-0.00001) return runNum; //small difference is to catch rounding errors
1432  }
1433  Error("GetRandomRunNumber","overran integrated luminosity for periodNumber=%d (%f vs %f)",periodNumber,lumiSum,lumi);
1434  throw std::runtime_error("Throwing 46.1: overran integrated luminosity for GetRandomRunNumber");
1435  return 0;
1436 }
1437 
1438 UInt_t CP::TPileupReweighting::GetRandomRunNumber(Int_t periodNumber, Double_t x) {
1439  if(!m_isInitialized) {
1440  if (m_printInfo) Info("GetRandomRunNumber","Initializing the subtool..");
1441  Initialize();
1442  }
1443  if(m_countingMode) { return 0; } //do nothing when in counting mode
1444 
1445  Period* p = m_periods[periodNumber];
1446  if(!p || !p->primaryHists[-1]) {
1447  Error("GetRandomRunNumber","Unrecognised periodNumber or no data loaded for: %d",periodNumber);
1448  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1449  }
1450 
1451  int bin = m_emptyHistogram->FindFixBin(x);
1452 
1453  double lumi = p->primaryHists[-1]->GetBinContent(bin) * p->sumOfWeights[-1] * m_random3->Rndm();
1454 
1455  if(!lumi) { return 0; } //no lumi with that mu
1456 
1457 
1458  //go through the runs assigned to this period ...
1459  double lumiSum=0;
1460  for(auto runNum : p->runNumbers) {
1461  lumiSum += m_runs[runNum].inputHists["None"]->GetBinContent(bin);
1462  if(lumiSum >= lumi-0.00001) return runNum; //small difference is to catch rounding errors
1463  }
1464  Error("GetRandomRunNumber","overran integrated luminosity for periodNumber=%d (%f vs %f)",periodNumber,lumiSum,lumi);
1465  throw std::runtime_error("Throwing 46.1: overran integrated luminosity for GetRandomRunNumber");
1466  return 0;
1467 }
1468 
1470  if(!m_isInitialized) {
1471  if (m_printInfo) {
1472  Info("GetRandomLumiBlockNumber","Initializing the subtool..");
1473  }
1474  Initialize();
1475  }
1476  if(m_countingMode) { return 0; } //do nothing when in counting mode
1477 
1478  double lumi = GetIntegratedLumi(runNumber,runNumber) * m_random3->Rndm() * 1E6 /* dont forget the lumi was divided by a million to get to pb */;
1479  double lumisum = 0.;
1480 
1481  for(auto& lbn : m_runs[runNumber].lumiByLbn) {
1482  if(m_unrepresentedDataAction==1 && IsUnrepresentedData(runNumber,lbn.second.second)) continue; //skip over bad lumi
1483  lumisum += lbn.second.first;
1484  if(lumisum >= lumi) return lbn.first;
1485  }
1486  Error("GetRandomLumiBlockNumber","overran integrated luminosity for RunNumber=%d (%f vs %f)",runNumber,lumi,lumisum);
1487  throw std::runtime_error("Throwing 46: overran integrated luminosity for runNumber");
1488  return 0;
1489 }
1490 
1491 //only considers periods assigned directly!
1493  if(!m_isInitialized) {
1494  if (m_printInfo) Info("GetRandomPeriodNumber","Initializing the subtool..");
1495  Initialize();
1496  }
1497  if(m_countingMode) { return 0; } //do nothing when in counting mode
1498 
1499  Period* p = m_periods[periodNumber];
1500  if(!p) {
1501  Error("GetRandomPeriodNumber","Unrecognised periodNumber: %d",periodNumber);
1502  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1503  }
1504 
1505  //if no subperiods, we just return ourselves
1506  if(p->subPeriods.size()==0) return p->id;
1507 
1508  //need the total lumi of this period
1509  double lumi = p->sumOfWeights[-1] * m_random3->Rndm();
1510 
1511  //loop over subPeriods ..
1512  double lumiSum=0;
1513  for(auto subp : p->subPeriods) {
1514  lumiSum += subp->sumOfWeights[-1];
1515  if(lumiSum >= lumi) return subp->id;
1516  }
1517 
1518  Error("GetRandomPeriodNumber","overran integrated luminosity for periodNumber=%d",periodNumber);
1519  throw std::runtime_error("Throwing 46.1: overran integrated luminosity for GetRandomPeriodNumber");
1520  return 0;
1521 
1522 }
1523 
1525  int bin = m_emptyHistogram->FindFixBin(x,y);
1526  return ( m_runs[runNumber].badBins[bin] );
1527 }
1528 
1529 //this method builds a file get can be friended to a TTree with a prw-hash branch
1530 Bool_t CP::TPileupReweighting::MakeWeightTree(TString channelNumbers, TString outFile, TString hashBranch, TString weightBranch) {
1531  if(!m_isInitialized ) {
1532  if (m_printInfo) Info("MakeWeightTree","Initializing the subtool..");
1533  Initialize();
1534  }
1535  TH1* hist = m_emptyHistogram.get();
1536  if(!hist) {
1537  Error("MakeWeightTree","Tool not configured properly ... please report this!");
1538  throw std::runtime_error("Throwing 47: Tool not configured properly ... please report this!");
1539  }
1540 
1541  //loop over given channels, and loop over all periods (except the global period) and get the pileup weight in each case for each bin
1542  TFile f1(outFile,"RECREATE");
1543 
1544  TTree* outTree = new TTree("prwTree","prwTree");
1545  ULong64_t prwHash(0);Float_t weight(0.);
1546  outTree->Branch(hashBranch,&prwHash);
1547  outTree->Branch(weightBranch,&weight);
1548 
1549  TObjArray *tx = channelNumbers.Tokenize(",");
1550  for (Int_t i = 0; i < tx->GetEntries(); i++) {
1551  int channelNumber = ((TObjString *)(tx->At(i)))->String().Atoi();
1552  if(channelNumber>999999 || channelNumber<0) {
1553  Error("MakeWeightTree","ChannelNumber can not be bigger than 999999 or less than 0 ... got %d",channelNumber);
1554  f1.Close();
1555  return 0;
1556  }
1557  for(auto& period : m_periods) {
1558  if(period.first==-1) continue;
1559  int periodNumber = period.first;
1560  for(int i=1;i<=hist->GetNbinsX();i++) {
1561  double x = hist->GetXaxis()->GetBinCenter(i);
1562  for(int j=1;j<=hist->GetNbinsY();j++) {
1563  double y = hist->GetYaxis()->GetBinCenter(j);
1564  weight = GetCombinedWeight(periodNumber,channelNumber,x,y);
1565  prwHash = GetPRWHash(periodNumber,channelNumber,x,y);
1566  outTree->Fill();
1567  }
1568  }
1569 
1570  }
1571  }
1572 
1573  outTree->BuildIndex(hashBranch.Data());
1574  outTree->Write();
1575  f1.Close();
1576 
1577  if (m_printInfo) {
1578  Info("MakeWeightTree","Successfully wrote prwTree to %s",outFile.Data());
1579  }
1580 
1581  return true;
1582 }
1583 
1584 ULong64_t CP::TPileupReweighting::GetPRWHash(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y) {
1585  if(!m_isInitialized) {
1586  if (m_printInfo) Info("GetPRWHash","Initializing the subtool..");
1587  Initialize();
1588  }
1589  TH1* hist = m_emptyHistogram.get();
1590  if(!hist) {
1591  Error("GetPRWHash","Tool not configured properly ... please report this!");
1592  throw std::runtime_error("Throwing 47: Tool not configured properly ... please report this!");
1593  }
1594 
1595  ULong64_t out = hist->FindFixBin(x,y);
1596  out += (unsigned long)(channelNumber)*10000000000;
1597  out += (unsigned long)(periodNumber)*10000;
1598 
1599  return out;
1600 
1601 }
1602 
1603 Float_t CP::TPileupReweighting::GetCombinedWeight(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y) {
1604  if(!m_isInitialized) {
1605  if (m_printInfo) Info("GetCombinedWeight","Initializing the subtool..");
1606  Initialize();
1607  }
1608  if(m_countingMode) return 0.;
1609 
1610  //decide how many dimensions this weight has - use the emptyHistogram to tell...
1611  TH1* hist = m_emptyHistogram.get();
1612  if(!hist) {
1613  Error("GetCombinedWeight","Tool not configured properly ... please report this!");
1614  throw std::runtime_error("Throwing 47: Tool not configured properly ... please report this!");
1615  }
1616 
1617  Float_t out = GetPeriodWeight(periodNumber,channelNumber)*GetPrimaryWeight(periodNumber,channelNumber,x);
1618  if(hist->GetDimension()>1) out *= GetSecondaryWeight(periodNumber,channelNumber,x,y);
1619 
1620  return out;
1621 }
1622 
1623 Float_t CP::TPileupReweighting::GetPeriodWeight(Int_t periodNumber, Int_t channelNumber) {
1624  //= L_A/L / N_A/N
1625  if(!m_isInitialized) {
1626  if (m_printInfo) Info("GetPeriodWeight","Initializing the subtool..");
1627  Initialize();
1628  }
1629  if(m_countingMode) return 0.;
1630 
1631  Period* p = m_periods[periodNumber];
1632  if(!p) {
1633  Error("GetPeriodWeight","Unrecognised periodNumber: %d",periodNumber);
1634  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1635  }
1636  if(p->sumOfWeights.find(channelNumber) == p->sumOfWeights.end()) {
1637  channelNumber = GetDefaultChannel(periodNumber);//p->defaultChannel;
1638  if(channelNumber!=0) {
1639  Warning("GetPeriodWeight","You're using a default config file ... you're gonna have a bad time!!");
1640  Warning("GetPeriodWeight","Please generate proper config files!");
1641  }
1642  }
1643 
1644  double n_a = p->sumOfWeights[channelNumber];
1645  double n = m_periods[-1]->sumOfWeights[channelNumber];
1646 
1647  double l_a = p->sumOfWeights[-1];
1648  double l = m_periods[-1]->sumOfWeights[-1];
1649 
1650  return (l_a/l) / (n_a/n);
1651 }
1652 
1653 Float_t CP::TPileupReweighting::GetPrimaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x) {
1654  //= L_i/L_A / N_i/N_A .. primaryHists have already been normalized
1655  if(!m_isInitialized ) {
1656  if (m_printInfo) Info("GetPrimaryWeight","Initializing the subtool..");
1657  Initialize();
1658  }
1659  if(m_countingMode) return 0.;
1660 
1661  Period* p = m_periods[periodNumber];
1662  if(!p) {
1663  Error("GetPrimaryWeight","Unrecognised periodNumber: %d",periodNumber);
1664  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1665  }
1666  int oChanNumber = channelNumber;
1667  if(p->sumOfWeights.find(channelNumber) == p->sumOfWeights.end()) channelNumber = GetDefaultChannel(periodNumber);//p->defaultChannel;
1668 
1669  if(!p->primaryHists[channelNumber]) {
1670  Error("GetPrimaryWeight","Unrecognised channelNumber %d for periodNumber %d",oChanNumber,periodNumber);
1671  throw std::runtime_error("Throwing 2: Unrecognised channelNumber");
1672  }
1673 
1674  int bin = p->primaryHists[channelNumber]->FindFixBin(x);
1675 
1676  double n = p->primaryHists[channelNumber]->GetBinContent(bin);
1677 
1678  // There is a special case that needs some additional protection
1679  // If you have the sum of weights ==0 for a given mu then you will have 0 in the PRW
1680  // profile leading to an exception or FPE later
1681  if (n==0){
1682  Warning("GetPrimaryWeight","No events expected in channelNumber %d for periodNumber %d with x=%g. Incorrect PRW profile or are you very unlucky???",channelNumber,periodNumber,x);
1683  return 1.;
1684  }
1685 
1686  if(!p->primaryHists[-1]) {
1687  Error("GetPrimaryWeight","No data loaded for period %d. Did you forget to load a lumicalc file or data config file?",periodNumber);
1688  throw std::runtime_error("Throwing 3: No data loaded. Did you forget to load a lumicalc file or data config file?");
1689  }
1690 
1691  double l = p->primaryHists[-1]->GetBinContent(bin);
1692 
1693  if (l==0 && n==0){
1694  Error("GetPrimaryWeight","No events expected with this mu. Incorrect PRW profile? Throwing exception ...");
1695  throw std::runtime_error(Form("Incorrect PRW profile detected. x=%g is not in the config file profile for channel %d, period %d", x, channelNumber, periodNumber ));
1696  }
1697 
1698  return l/n;
1699 }
1700 
1701 Float_t CP::TPileupReweighting::GetSecondaryWeight(Int_t periodNumber, Int_t channelNumber,Float_t x,Float_t y) {
1702  //= L_j/L_i / N_j/N_i .. secondary hists have already been normalized
1703  if(!m_isInitialized ) {
1704  if (m_printInfo) Info("GetSecondaryWeight","Initializing the subtool..");
1705  Initialize();
1706  }
1707  if(m_countingMode) return 0.;
1708 
1709  Period* p = m_periods[periodNumber];
1710  if(!p) {
1711  Error("GetSecondaryWeight","Unrecognised periodNumber: %d",periodNumber);
1712  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1713  }
1714  if(p->sumOfWeights.find(channelNumber) == p->sumOfWeights.end()) channelNumber = GetDefaultChannel(periodNumber);//p->defaultChannel;
1715  int bin = p->secondaryHists[channelNumber]->FindFixBin(x,y);
1716  double n = p->secondaryHists[channelNumber]->GetBinContent(bin);
1717  double l = p->secondaryHists[-1]->GetBinContent(bin);
1718 
1719  return l/n;
1720 }
1721 
1722 Double_t CP::TPileupReweighting::GetDataWeight(Int_t runNumber, const TString& trigger) {
1723  //special mu-independent version of GetDataWeight. Will just use the global luminosity
1724  m_doGlobalDataWeight=true;
1725  double out = GetDataWeight(runNumber, trigger, 0);
1726  m_doGlobalDataWeight=false;
1727  return out;
1728 }
1729 
1730 Double_t CP::TPileupReweighting::GetDataWeight(Int_t runNumber, const TString& trigger, Double_t x, bool runDependent) {
1731 
1732  if(!m_isInitialized ) {
1733  if (m_printInfo) Info("GetDataWeight","Initializing the subtool..");
1734  Initialize();
1735  }
1736  if(m_countingMode) return 0.;
1737 
1738  //determine which period this run number is in
1739  Int_t periodNumber = GetFirstFoundPeriodNumber(runNumber);
1740 
1741  Period* p = m_periods[periodNumber];
1742  if(!p) {
1743  Error("GetDataWeight","Unrecognised periodNumber: %d",periodNumber);
1744  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1745  }
1746 
1747  //see if we already made this trigger hist
1748  auto tobj = m_triggerObjs.find(trigger);
1749 
1750 
1751  if(tobj==m_triggerObjs.end()) {
1752  //try to do a subtrigger calculation we need to reopen all the lumicalc files and do the fiddily prescaled luminosity calculation
1753  //will throw errors if can't
1754  CalculatePrescaledLuminosityHistograms(trigger,0);
1755  tobj = m_triggerObjs.find(trigger);
1756  }
1757 
1758  //check
1759 
1760 
1761 
1762  if(tobj->second->triggerHists.find(p->id) == tobj->second->triggerHists.end()) {
1763  Error("GetDataWeight","Could not find trigger %s in period %d",trigger.Data(),p->id);
1764  throw std::runtime_error("GetDataWeight: Could not find trigger 1");
1765  }
1766 
1767 
1768  //now we need to evaluate the trigger bits, and if necessary calculate the new trigger hist for the new bits combination
1769  long tbits = tobj->second->getBits(this);
1770 
1771  if(tbits==0) return 1; //no trigger passed
1772 
1773  int idx = (runDependent) ? -runNumber : p->id;
1774 
1775 
1776  auto dItr = tobj->second->triggerHists[idx].find(tbits);
1777 
1778  TH1* denomHist = 0;
1779  if(dItr == tobj->second->triggerHists[idx].end()) {
1780  //possible that the tbits for this event have not been encountered before, so just redo the calculation ...
1781  calculateHistograms(tobj->second.get(),(runDependent) ? runNumber:0);
1782 
1783  denomHist = tobj->second->triggerHists[idx][tbits].get();
1784  if(denomHist==0) {
1785  Error("GetDataWeight","Could not find trigger %s in period %d with bits %ld",trigger.Data(),idx, tbits);
1786  throw std::runtime_error("GetDataWeight: Could not find trigger 2");
1787  }
1788  } else {
1789  denomHist = dItr->second.get();
1790  }
1791 
1792 
1793  TH1* numerHist = m_triggerObjs["None"]->triggerHists[idx][0].get();
1794  if(numerHist==0) {
1795  Error("GetDataWeight","Could not find unprescaled trigger in period %d",idx);
1796  throw std::runtime_error("GetDataWeight: Could not find unprescaled trigger. Please AddLumiCalc with a 'None' trigger");
1797  }
1798 
1799  if(m_doGlobalDataWeight) return numerHist->Integral(0,numerHist->GetNbinsX()+1)/denomHist->Integral(0,denomHist->GetNbinsX()+1);
1800 
1801  Int_t bin=(m_doPrescaleWeight) ? numerHist->FindFixBin(x) : //DO NOT SHIFT IF GETTING A PRESCALE WEIGHT - applied to MC, we only shift incoming mu if we are data
1802  numerHist->FindFixBin(x*m_dataScaleFactorX); //if getting a data weight (i.e. running on data) MUST SCALE BY THE DATA SCALE FACTOR! (assume incoming is raw unscaled)
1803 
1804  //we also need to redirect the binning if necessary (unrepData=3)
1805  if( (!m_doPrescaleWeight) && m_unrepresentedDataAction==3) bin = p->inputBinRedirect[bin];
1806 
1807  if(!denomHist->GetBinContent(bin)) {
1808  if(m_doPrescaleWeight) return -1; //happens if trigger was disabled/unavailable for that mu, even though that mu is in the dataset
1809  Error("GetDataWeight","Unrecognised mu value %f ... are you sure you included all lumicalc files",x);
1810  throw std::runtime_error("GetDataWeight: Unrecognised mu value. Please AddLumiCalc enough lumi with 'None' trigger");
1811  }
1812 
1813  return numerHist->GetBinContent(bin)/denomHist->GetBinContent(bin);
1814 }
1815 
1816 
1817 Double_t CP::TPileupReweighting::GetPrescaleWeight(Int_t runNumber, const TString& trigger) {
1818  //special mu-independent version of GetPrescaleWeight. Will just use the global luminosity
1819  m_doGlobalDataWeight=true;
1820  double out = GetPrescaleWeight(runNumber, trigger, 0);
1821  m_doGlobalDataWeight=false;
1822  return out;
1823 }
1824 
1825 Double_t CP::TPileupReweighting::GetPrescaleWeight(Int_t runNumber, const TString& trigger, Double_t x, bool runDependent) {
1826 
1827  m_doPrescaleWeight = true;
1828  double out = GetDataWeight(runNumber,trigger,x,runDependent);
1829  m_doPrescaleWeight = false;
1830  if(out<=0) return 0; //happens when triggers disabled/unavailable for a given mu ... therefore the prescale weight is 0
1831  return 1./out;
1832 
1833 }
1834 
1835 
1836 //fills the appropriate inputHistograms
1837 Int_t CP::TPileupReweighting::Fill(Int_t runNumber,Int_t channelNumber,Float_t w,Float_t x, Float_t y) {
1838  if(!m_isInitialized ) {
1839  if (m_printInfo) Info("Fill","Initializing the subtool..");
1840  Initialize();
1841  }
1842  //should only be given genuine mcRunNumbers if mc (channel>=0). We don't fill periodNumber distributions
1843 
1844  TH1* hist = 0;
1845 
1846  //get the period if mc, get the data run if data
1847  if(channelNumber>=0) {
1848  if(m_periods.find(runNumber)==m_periods.end()) {
1849  //if assume default binning is true, will add binning from 0 to 100 for this run number
1850  if(m_autoRunStart || m_autoRunEnd) {
1851  AddPeriod(runNumber,m_autoRunStart,m_autoRunEnd);
1852  }
1853  if(m_periods.find(runNumber)==m_periods.end()) {
1854  Error("Fill","Unrecognised runNumber: %d. Check your period configuration (AddPeriod or UsePeriodConfig)",runNumber);
1855  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1856  }
1857  }
1858  Period* p = m_periods[runNumber];
1859  if(!p) {
1860  Error("Fill","Unrecognised runNumber: %d. Check your period configuration (AddPeriod or UsePeriodConfig) ... but should never have got here so please report this!",runNumber);
1861  throw std::runtime_error("Throwing 1: Unrecognised periodNumber");
1862  }
1863  std::unique_ptr<TH1>& histptr = p->inputHists[channelNumber];
1864  if(!histptr) {
1865  //need to create my period histogram
1866  histptr = CloneEmptyHistogram(runNumber,channelNumber);
1867  }
1868  hist = histptr.get();
1869  } else {
1870  Run& r = m_runs[runNumber];
1871  if( ! r.inputHists["None"]) r.inputHists["None"]=CloneEmptyHistogram(runNumber,channelNumber);
1872  hist = r.inputHists["None"].get();
1873  }
1874 
1875  if(!hist) {
1876  Error("Fill","Unknown [run,channel] = [%d,%d]",runNumber,channelNumber);
1877  throw std::runtime_error("Throwing 45: Fill:: Unknown [run,channel] ");
1878  }
1879 
1880  if(hist->GetDimension()==1) {
1881  return hist->Fill(x,w);
1882  } else if(hist->GetDimension()==2) {
1883  return (static_cast<TH2*>(hist))->Fill(x,y,w);
1884  }
1885  return -1;
1886 }
1887 
1889 
1890  if(!m_countingMode) {Warning("WriteToFile","Not in counting mode, so no file will be written");return 0;}
1891 
1892 
1893  //build a TTree with the correct mc structure, and dump the mc histogram info in to it
1894  //also build aggregate histograms across all channels - shouldn't be used as input histograms unless for a single channel
1895  TString filename = fname;
1896  filename += (filename=="") ? TString(this->GetName()) + ".prw.root" : "";
1897 
1898  //loop over the weights. "pileup" gets it own MCPileupReweighting ttree. ... since we got rid of weightNames in this version, everything goes in MCPileupReweighting now!
1899  //data goes in to DataCustomReweighting. other goes in to MCCustomReweighting
1900  TFile* outFile = TFile::Open(filename,"RECREATE");
1901  Int_t r = WriteToFile(outFile);
1902  outFile->Close();
1903  return r;
1904 }
1905 
1906 
1908  if(!m_countingMode) {Warning("WriteToFile","Not in counting mode, so no file will be written");return 0;}
1909 
1910  TDirectory* origDir = gDirectory;
1911 
1912  outFile->cd();
1913 
1914  std::unique_ptr< TTree > outTreeMC;
1915  std::unique_ptr< TTree > outTreeData;
1916  Int_t channel = 0; UInt_t runNumber = 0;
1917  auto pStarts = std::make_unique< std::vector< UInt_t > >();
1918  auto pStartsPtr = pStarts.get();
1919  auto pEnds = std::make_unique< std::vector< UInt_t > >();
1920  auto pEndsPtr = pEnds.get();
1921  Char_t histName[150];
1922 
1923 
1924  //loop over periods ... periods only get entry in table if they have an input histogram
1925  for(auto period : m_periods) {
1926  if(period.first != period.second->id) continue; //skips redirects
1927  if(period.first<0) continue; //avoid the global run number
1928  runNumber = period.first;
1929  pStarts->clear();
1930  pEnds->clear();
1931  if(period.second->subPeriods.size()==0) {
1932  pStarts->push_back(period.second->start); pEnds->push_back(period.second->end);
1933  }
1934  else {
1935  for(auto subp : period.second->subPeriods) {
1936  pStarts->push_back(subp->start); pEnds->push_back(subp->end);
1937  }
1938  }
1939  for(auto& inHist : period.second->inputHists) {
1940  channel = inHist.first;
1941  TH1* hist = inHist.second.get();
1942  strncpy(histName,hist->GetName(),sizeof(histName)-1);
1943  hist->Write();
1944 
1945  if(!outTreeMC) {
1946  outTreeMC.reset( new TTree("MCPileupReweighting","MCPileupReweighting") );
1947  outTreeMC->Branch("Channel",&channel);
1948  outTreeMC->Branch("RunNumber",&runNumber);
1949  outTreeMC->Branch("PeriodStarts",&pStartsPtr);
1950  outTreeMC->Branch("PeriodEnds",&pEndsPtr);
1951  outTreeMC->Branch("HistName",&histName,"HistName/C");
1952  }
1953  outTreeMC->Fill();
1954  }
1955  }
1956 
1957  //loop over data
1958  for(auto& run : m_runs) {
1959  runNumber = run.first;
1960  if(run.second.inputHists.find("None")==run.second.inputHists.end()) continue;
1961 
1962  TH1* hist = run.second.inputHists["None"].get();
1963  strncpy(histName,hist->GetName(),sizeof(histName)-1);
1964  hist->Write();
1965  if(!outTreeData) {
1966  outTreeData.reset( new TTree("DataPileupReweighting","DataPileupReweighting") );
1967  outTreeData->Branch("RunNumber",&runNumber);
1968  outTreeData->Branch("HistName",&histName,"HistName/C");
1969  }
1970  outTreeData->Fill();
1971  }
1972 
1973 
1974  //write the non-zero ttrees
1975  if( outTreeMC ) {
1976  outTreeMC->Write();
1977  }
1978  if( outTreeData ) {
1979  outTreeData->Write();
1980  }
1981 
1982  if (m_printInfo) {
1983  Info("WriteToFile", "Successfully generated config file: %s",outFile->GetName());
1984  Info("WriteToFile", "Happy Reweighting :-)");
1985  }
1986 
1987  gDirectory = origDir;
1988 
1989  return 0;
1990 }
1991 
1992 
1993 //=============================================================================
1994 // Method to calculate the ratio histogram
1995 //=============================================================================
1997  // normalize the data histogram based on which sort of histogram it is
1998  if(hist){
1999  if(hist->InheritsFrom("TH3")) {
2000  Error("normalizeHistogram","3D reweighting not supported yet");
2001  throw std::runtime_error("Throwing 3: 3D reweighting not supported yet");
2002  }
2003  else if(hist->InheritsFrom("TH2")) {
2004  bool skipNorm=false;
2005  //normalize each bin according to the projection in x
2006  TH1D* proj = ((TH2*)hist)->ProjectionX();
2007  Int_t bin,binx,biny,binz;
2008  for(binz=1; binz<=hist->GetNbinsZ(); binz++) {
2009  for(biny=1; biny<=hist->GetNbinsY(); biny++) {
2010  for(binx=1; binx<=hist->GetNbinsX(); binx++) {
2011  bin = hist->GetBin(binx,biny,binz);
2012  Double_t value = hist->GetBinContent(bin);
2013  Double_t normalizer = proj->GetBinContent(binx);
2014  if(normalizer!=0.0) {
2015  hist->SetBinContent(bin,value/normalizer);
2016  } else {
2017  skipNorm=true;
2018  }
2019  }
2020  }
2021  }
2022  delete proj;
2023  if(skipNorm && m_debugging) Warning("normalizeHistogram","Skipped normalization in hist %s",hist->GetName());
2024  } else {
2025  //normalize to the sum of weights
2026  if(hist->GetSumOfWeights()!=0.0) {
2027  hist->Scale(1.0/hist->GetSumOfWeights());
2028  } else {
2029  if (m_debugging) Warning("normalizeHistogram","Skipping Normalizing histogram %s to ~zero: %f",hist->GetName(),hist->GetSumOfWeights());
2030  }
2031  }
2032  } else {
2033  Error("normalizeHistogram","Non existent histogram for normalizing");throw std::runtime_error("Throwing 56");
2034  }
2035 }
2036 
2037 
2038 //This allows PROOF to merge the generated histograms before the WriteToFile or WriteCustomToFile calls
2039 Int_t CP::TPileupReweighting::Merge(TCollection *coll) {
2040  if(!coll) return 0;
2041  if(coll->IsEmpty()) return 0;
2042 
2043  // Iterate over the elements of the collection:
2044  TIter next( coll );
2045  TObject* obj = 0;
2046  while( ( obj = next() ) ) {
2047 
2048  // Check that the element is of the right type:
2049  CP::TPileupReweighting* vobj = dynamic_cast< CP::TPileupReweighting* >( obj );
2050  if( ! vobj ) {
2051  Error( "Merge", "Unknown object type encountered: %s",obj->ClassName() );
2052  return 0;
2053  }
2054 
2055  //merge the inputHistograms ... all the periods should be identical
2056  for(auto period : vobj->m_periods) {
2057  if(period.first != period.second->id) continue;
2058  for(auto& iHist : period.second->inputHists) {
2059  if(GetInputHistogram(period.first,iHist.first)==0) {
2060  m_periods[period.first]->inputHists[iHist.first].reset(
2061  dynamic_cast<TH1*>(iHist.second->Clone(iHist.second->GetName())) );
2062  m_periods[period.first]->inputHists[iHist.first]->SetDirectory(0);
2063  } else {
2064  GetInputHistogram(period.first,iHist.first)->Add(iHist.second.get());
2065  }
2066  }
2067  }
2068  //also must remember to merge the runs too (where the data is held)
2069  for(auto& run : vobj->m_runs) {
2070  for(auto& iHist : run.second.inputHists) {
2071  if( ! m_runs[run.first].inputHists[iHist.first] ) {
2072  m_runs[run.first].inputHists[iHist.first].reset(
2073  dynamic_cast<TH1*>( iHist.second->Clone(iHist.second->GetName())) );
2074  m_runs[run.first].inputHists[iHist.first]->SetDirectory(0);
2075  } else {
2076  m_runs[run.first].inputHists[iHist.first]->Add(iHist.second.get());
2077  }
2078  }
2079  }
2080  }
2081 
2082  return 1;
2083 }
2084 
2085 
2086 void CP::TPileupReweighting::CalculatePrescaledLuminosityHistograms(const TString& trigger, int runDependentRun) {
2087  //check if this is a composite trigger (user has provided OR of triggers that are passed before prescale)
2088 
2089  TString triggerCopy = trigger; triggerCopy.ReplaceAll(" ",""); triggerCopy.ReplaceAll("&&","&");triggerCopy.ReplaceAll("||","|");
2090  auto t = makeTrigger(triggerCopy);
2091  if( ! t ) {
2092  Error("GetDataWeight","Error parsing composite trigger set: %s",trigger.Data());
2093  throw std::runtime_error("Could not parse composite trigger");
2094  }
2095 
2096  std::vector<TString> subTriggers;
2097 
2098  t->getTriggers(subTriggers); //fills the vector
2099 
2100  t->subTriggers = subTriggers; //cache the vector of triggers
2101 
2102  calculateHistograms(t.get(),runDependentRun);
2103 
2104  //save trigger object
2105  m_triggerObjs[trigger] = std::move( t );
2106 
2107 }
2108 
2110  TH1* unprescaledLumi = 0;
2111  if(runDependentRun) {
2112  unprescaledLumi = m_triggerObjs["None"]->triggerHists[-runDependentRun][0].get();
2113  }
2114  bool fillUnprescaledLumi = (unprescaledLumi==0); //only fill the unprescaled histogram once ... so if it already exists, it must have already been filled
2115 
2116  //now we need the trigger bits for this trigger for this event
2117  long tbits = t->getBits(this);
2118  //1. Open all the lumicalc files
2119  //2. For each entry in the 'None' lumicalc file, lookup the run and lb in the subTriggers, and get the prescales (if not found, assume trigger not in, prescale=infinity)
2120  //3. perform calculation of lumi as: Lumi*(1 - Prod_triggers(1-1/prescale)) .. sum this up into the output histogram
2121  // Note: Put each lumi into the period that contains the run ( use )
2122 
2123  TDirectory* origDir = gDirectory;
2124 
2125  std::map<TString, std::map<Int_t, std::map<Int_t, Float_t> > > prescaleByRunAndLbn;
2126 
2127  //Load the prescales for each lumicalc subTrigger
2128  for(std::vector<TString>::iterator it = t->subTriggers.begin();it!=t->subTriggers.end();++it) {
2129  //Info("aa","subtrigger %s Int lumi is %f", it->Data(), m_triggerPrimaryDistributions[*it][1]->Integral());
2130  for(auto& fileName : m_lumicalcFiles[*it]) {
2131  //if the file is also in the "None" list, then we will now the prescale must be 0
2132  //it seems that in these no-trigger files used for unprescaled lumi the prescales are all set to -1
2133  bool isUnprescaled = (std::find(m_lumicalcFiles["None"].begin(),m_lumicalcFiles["None"].end(),fileName)!=m_lumicalcFiles["None"].end());
2134 
2135  TFile* rootFile = TFile::Open( fileName, "READ" );
2136  if ( rootFile->IsZombie() ) {
2137  Error("CalculatePrescaledLuminosityHistograms","Could not open file: %s",fileName.Data());
2138  std::string toThrow = "Throwing 6: Could not open file: "; toThrow += fileName.Data();
2139  throw std::runtime_error(toThrow);
2140  } else {
2141  //try to get the the known TTrees
2142  TTree *tmp = (TTree*)rootFile->Get( "LumiMetaData" );
2143  if(tmp) {
2144  //structure expected is as given by iLumiCalc:
2145  // RunNbr, AvergeInteractionPerXing, IntLumi
2146  UInt_t runNbr=0;Float_t ps1=0;Float_t ps2=0; Float_t ps3=0;UInt_t lbn=0;TBranch *b_runNbr=0;TBranch *b_L1Presc=0;TBranch *b_L2Presc=0;TBranch *b_L3Presc=0;TBranch *b_lbn=0;
2147  if(tmp->SetBranchAddress("RunNbr",&runNbr,&b_runNbr)!=0) {
2148  Error("CalculatePrescaledLuminosityHistograms","Could not find RunNbr branch in Data TTree");throw std::runtime_error("Could not find RunNbr branch in Data TTree");
2149  }
2150  if(tmp->SetBranchAddress("L1Presc",&ps1,&b_L1Presc)!=0) {
2151  Error("CalculatePrescaledLuminosityHistograms","Could not find L1Presc branch in Data TTree");throw std::runtime_error("Could not find L1Presc branch in Data TTree");
2152  }
2153  if(tmp->SetBranchAddress("L2Presc",&ps2,&b_L2Presc)!=0) {
2154  Error("CalculatePrescaledLuminosityHistograms","Could not find L2Presc branch in Data TTree");throw std::runtime_error("Could not find L2Presc branch in Data TTree");
2155  }
2156  if(tmp->SetBranchAddress("L3Presc",&ps3,&b_L3Presc)!=0) {
2157  Error("CalculatePrescaledLuminosityHistograms","Could not find L3Presc branch in Data TTree");throw std::runtime_error("Could not find L3Presc branch in Data TTree");
2158  }
2159  if(tmp->SetBranchAddress("LBStart",&lbn,&b_lbn)!=0) {
2160  Error("CalculatePrescaledLuminosityHistograms","Could not find LBStart branch in Data TTree");throw std::runtime_error("Could not find LBStart branch in Data TTree");
2161  }
2162  long nEntries = tmp->GetEntries();
2163  for(long i=0;i<nEntries;i++) {
2164  b_runNbr->GetEntry(i);b_L1Presc->GetEntry(i);b_L2Presc->GetEntry(i);b_L3Presc->GetEntry(i);b_lbn->GetEntry(i);
2165  runNbr += m_lumicalcRunNumberOffset;
2166  if(runDependentRun && int(runNbr)!=runDependentRun) continue; //only use the given run with doing calculation run-dependently
2167  //save the prescale by run and lbn
2168  //if(runNbr==215643)
2169  if(m_debugging && m_printInfo) Info("...","prescale in [%d,%d] = %f %f %f", runNbr,lbn,ps1,ps2,ps3);
2170  if(ps1>0&&ps2>0&&ps3>0) prescaleByRunAndLbn[*it][runNbr][lbn] = ps1*ps2*ps3;
2171  else if(isUnprescaled) prescaleByRunAndLbn[*it][runNbr][lbn] = 1; //special case where the trigger is an unprescaled one and user is reusing the unprescaled lumicalc
2172  }
2173  }
2174  rootFile->Close();
2175  }
2176  delete rootFile;
2177  } //end of loop over lumicalc files for this trigger
2178  }
2179 
2180  for(auto& fileName : m_lumicalcFiles["None"]) {
2181  // Load the data ROOT file for the 'None' lumicalc file
2182  TFile* rootFile = TFile::Open( fileName, "READ" );
2183  if ( rootFile->IsZombie() ) {
2184  Error("CalculatePrescaledLuminosityHistograms","Could not open file: %s",fileName.Data());
2185  std::string toThrow = "Throwing 6: Could not open file: "; toThrow += fileName.Data();
2186  throw std::runtime_error(toThrow);
2187  } else {
2188  //try to get the the known TTrees
2189  TTree *tmp = (TTree*)rootFile->Get( "LumiMetaData" );
2190  if(tmp) {
2191  //structure expected is as given by iLumiCalc:
2192  // RunNbr, AvergeInteractionPerXing, IntLumi
2193  UInt_t runNbr=0;Float_t intLumi=0;UInt_t lbn=0;TBranch *b_runNbr=0;TBranch *b_intLumi=0;TBranch *b_lbn=0;
2194  Float_t mu=0.; TBranch *b_mu=0;
2195  if(tmp->SetBranchAddress("RunNbr",&runNbr,&b_runNbr)!=0) {
2196  Error("CalculatePrescaledLuminosityHistograms","Could not find RunNbr branch in Data TTree");throw std::runtime_error("Could not find RunNbr branch in Data TTree");
2197  }
2198  if(tmp->SetBranchAddress("AvergeInteractionPerXing",&mu,&b_mu)!=0) {
2199  Error("CalculatePrescaledLuminosityHistograms","Could not find AvergeInteractionPerXing branch in Data TTree");throw std::runtime_error("Could not find AvergeInteractionPerXing branch in Data TTree");
2200  }
2201  if(tmp->SetBranchAddress("IntLumi",&intLumi,&b_intLumi)!=0) {
2202  Error("CalculatePrescaledLuminosityHistograms","Could not find IntLumi branch in Data TTree");throw std::runtime_error("Could not find IntLumi branch in Data TTree");
2203  }
2204  if(tmp->SetBranchAddress("LBStart",&lbn,&b_lbn)!=0) {
2205  Error("CalculatePrescaledLuminosityHistograms","Could not find LBStart branch in Data TTree");throw std::runtime_error("Could not find LBStart branch in Data TTree");
2206  }
2207  long nEntries = tmp->GetEntries();
2208  //this loop ends up calling the triggerBeforePrescale method each time, but this should hopefully stop once the loading is done
2209  for(long i=0;i<nEntries;i++) {
2210  b_runNbr->GetEntry(i);b_intLumi->GetEntry(i);b_mu->GetEntry(i);
2211  b_lbn->GetEntry(i);
2212 
2213  runNbr += m_lumicalcRunNumberOffset;
2214  if(runDependentRun && int(runNbr)!=runDependentRun) continue; //only use the given run with doing calculation run-dependently
2215  //for each subtrigger, lookup prescale, and calculate pFactor
2216 
2217  double pFactor = t->eval(prescaleByRunAndLbn,runNbr,lbn,this);
2218 
2219  //Info("...","prescale in [%d,%d] = %f", runNbr,lbn,1./pFactor);
2220 
2221  int bin = m_emptyHistogram->FindFixBin(mu*m_dataScaleFactorX);
2222 
2223  //add into all periods that contain this runNbr
2224  bool firstFill=false;
2225  for(auto p : m_periods) {
2226  if(p.first != p.second->id) continue; //skips remappings
2227  if(!p.second->contains(runNbr)) continue;
2228  if(runDependentRun && firstFill) continue; //only fill runDependentRun once, since they aren't kept in the period they go in their own run-indexed histogram
2229 
2230  int idx = (runDependentRun) ? -runDependentRun : p.second->id; //use negative runNumber to avoid clash with periodID
2231 
2232 
2233  auto& triggerHists = t->triggerHists[idx];
2234  std::unique_ptr<TH1>& histptr = triggerHists[tbits];
2235  if(!histptr) {
2236  histptr = CloneEmptyHistogram(p.first,-1);
2237  if(m_debugging && m_printInfo) Info("CalculatePrescaledLuminosityHistograms","Created Data Weight Histogram for [%s,%d,%d,%ld]",t->val.Data(),p.first,idx,tbits);
2238  }
2239  //check if we were about to fill a bad bin ... if we are, we either skipp the fill (unrep action=1) or redirect (unrep action=3)
2240  if( (m_unrepresentedDataAction==1) && p.second->inputBinRedirect[bin]!=bin) { } //do nothing
2241  else if( m_unrepresentedDataAction==3 ) {histptr->Fill(triggerHists[tbits]->GetBinCenter(p.second->inputBinRedirect[bin]), intLumi*pFactor);}
2242  else triggerHists[tbits]->Fill(mu*m_dataScaleFactorX,intLumi*pFactor);
2243 
2244 
2245  if(runDependentRun && fillUnprescaledLumi) {
2246  //need to fill unprescaled lumi histogram for this particular run ...
2247  if(unprescaledLumi==0) {
2248  m_triggerObjs["None"]->triggerHists[-runDependentRun][0] = CloneEmptyHistogram(p.first,-1);
2249  unprescaledLumi = m_triggerObjs["None"]->triggerHists[-runDependentRun][0].get();
2250  }
2251  if( (m_unrepresentedDataAction==1) && p.second->inputBinRedirect[bin]!=bin) { } //do nothing
2252  else if( m_unrepresentedDataAction==3 ) {unprescaledLumi->Fill(unprescaledLumi->GetBinCenter(p.second->inputBinRedirect[bin]), intLumi);}
2253  else unprescaledLumi->Fill(mu*m_dataScaleFactorX,intLumi);
2254  }
2255 
2256  firstFill=true;
2257 
2258  }
2259  }
2260  rootFile->Close();
2261  }
2262  }
2263  delete rootFile;
2264  } //end loop over unprescaled lumicalc files
2265 
2266 
2267  //Info("aa","Int lumi is %f", m_triggerPrimaryDistributions[trigger][1]->Integral());
2268 
2269 
2270 
2271  // Return to the original ROOT directory
2272  gDirectory = origDir;
2273 
2274 
2275 }
2276 
2277 std::unique_ptr<CP::TPileupReweighting::CompositeTrigger>
2279 
2280  if( m_debugging && m_printInfo ) {
2281  Info( "makeTrigger", "Doing %s", s.Data() );
2282  }
2283 
2284  // Find the first operand
2285  TString oper1;
2286  std::unique_ptr< CompositeTrigger > cOper1;
2287  if( s.BeginsWith( "(" ) ) {
2288  // Find closing bracket
2289  int bCounter=1; int i=1;
2290  while( ( bCounter != 0 ) && ( i <= s.Length() ) ) {
2291  if( s( i ) == '(' ) {
2292  bCounter++;
2293  } else if( s( i ) == ')' ) {
2294  bCounter--;
2295  }
2296  i++;
2297  }
2298  if( bCounter != 0 ) {
2299  Error( "makeTrigger", "Missing closing bracket" );
2300  return std::unique_ptr< CompositeTrigger >();
2301  }
2302  oper1 = s( 1, i - 2 );
2303  if( i == s.Length() + 1 ) {
2304  return makeTrigger( oper1 ); //meainingless brackets, just return evaluation
2305  } else {
2306  cOper1 = makeTrigger( oper1 );
2307  if( ! cOper1 ) {
2308  return 0;
2309  }
2310  }
2311  }
2312 
2313  // Find the second operand. Do this by finding the first operation at
2314  // this level. If cannot find, then just set val
2315  int i=0;int bCounter=0; int op=0;
2316  while( i <= s.Length() ) {
2317  if( s( i ) == '(' ) {
2318  bCounter++;
2319  } else if( s( i ) == ')' ) {
2320  bCounter--;
2321  } else if( ( s( i ) == '&' ) && ( bCounter == 0 ) ) {
2322  op = 2;
2323  break;
2324  } else if( ( s( i ) == '|' ) && ( bCounter == 0 ) ) {
2325  op = 1;
2326  break;
2327  }
2328  i++;
2329  }
2330 
2331 
2332  if( op == 0 ) {
2333  if( m_lumicalcFiles.find( s ) == m_lumicalcFiles.end() ) {
2334  Error( "GetDataWeight", "Could not find subTrigger %s", s.Data() );
2335  return std::unique_ptr< CompositeTrigger >();
2336  }
2337  if( cOper1 ) {
2338  return cOper1;
2339  }
2340  auto out = std::make_unique< CompositeTrigger >();
2341  out->op = op;
2342  out->val = s;
2343  return out; //just a value trigger
2344  }
2345 
2346  auto out = std::make_unique< CompositeTrigger >();
2347  out->op = op;
2348 
2349  if( op == 1 ) { //an OR
2350  if( cOper1 ) {
2351  out->trig1 = std::move( cOper1 );
2352  } else {
2353  oper1 = s( 0, i );
2354  out->trig1 = makeTrigger( oper1 );
2355  if( ! out->trig1 ) {
2356  return std::unique_ptr< CompositeTrigger >();
2357  }
2358  }
2359  TString oper2 = s( i + 1, s.Length() );
2360  out->trig2 = makeTrigger( oper2 );
2361  if( ! out->trig2 ) {
2362  return std::unique_ptr< CompositeTrigger >();
2363  }
2364  return out;
2365  }
2366 
2367  // Got here, must be an AND keep going until we hit an OR at this level,
2368  // then make all that the first operand (delete cOper if necessary)
2369  int j = i;
2370  bCounter = 0;
2371  while( j <= s.Length() ) {
2372  if( s( j ) == '(' ) {
2373  bCounter++;
2374  } else if( s( j )== ')' ) {
2375  bCounter--;
2376  } else if( ( s( j ) == '|' ) && ( bCounter == 0 ) ) {
2377  op=1;
2378  break;
2379  }
2380  j++;
2381  }
2382  if( ( j == s.Length() + 1 ) && ( op == 2 ) ) { //no more OR found, set oper2 to remainder
2383  if( cOper1 ) {
2384  out->trig1 = std::move( cOper1 );
2385  } else {
2386  oper1 = s( 0, i );
2387  out->trig1 = makeTrigger( oper1 );
2388  }
2389  TString oper2 = s( i + 1, s.Length() );
2390  if( m_debugging && m_printInfo ) {
2391  Info( "makeTrigger", "Found & %s %s", oper1.Data(), oper2.Data() );
2392  }
2393  out->trig2 = makeTrigger( oper2 );
2394  if( ! out->trig2 ) {
2395  return std::unique_ptr< CompositeTrigger >();
2396  }
2397  return out;
2398  } else if( op ==1 ) { //found an OR, set oper1 to everything up to this
2399 
2400  oper1 = s( 0, j );
2401  TString oper2 = s( j + 1, s.Length() );
2402  if( m_debugging && m_printInfo ) {
2403  Info( "makeTrigger", "Found & then | %s %s", oper1.Data(),
2404  oper2.Data() );
2405  }
2406  out->op = op; //updates to an OR
2407  out->trig1 = makeTrigger( oper1 );
2408  if( ! out->trig1 ) {
2409  return std::unique_ptr< CompositeTrigger >();
2410  }
2411  out->trig2 = makeTrigger( oper2 );
2412  if( ! out->trig2 ) {
2413  return std::unique_ptr< CompositeTrigger >();
2414  }
2415  return out;
2416  }
2417  Error( "makeTrigger", "Should never have got here!, but did with %s",
2418  s.Data() );
2419  return std::unique_ptr< CompositeTrigger >();
2420 }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
python.AtlRunQueryAMI.period
period
Definition: AtlRunQueryAMI.py:225
beamspotman.r
def r
Definition: beamspotman.py:676
CP::TPileupReweighting::AddDistribution
Int_t AddDistribution(TH1 *hist, Int_t runNumber, Int_t channelNumber)
Definition: TPileupReweighting.cxx:728
python.selector.AtlRunQuerySelectorLhcOlc.runNumbers
string runNumbers
Definition: AtlRunQuerySelectorLhcOlc.py:586
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
CP::TPileupReweighting::Fill
Int_t Fill(Int_t runNumber, Int_t channelNumber, Float_t w, Float_t x, Float_t y=0.)
Definition: TPileupReweighting.cxx:1837
CP::TPileupReweighting::AddDistributionTree
void AddDistributionTree(TTree *tree, TFile *file)
Definition: TPileupReweighting.cxx:648
CP::TPileupReweighting::IsBadBin
Int_t IsBadBin(Int_t thisMCRunNumber, Int_t bin)
Definition: TPileupReweighting.cxx:1351
CP::TPileupReweighting::SetUniformBinning
Int_t SetUniformBinning(Int_t nbinsx, Double_t xlow, Double_t xup, Int_t nbinsy=0, Double_t ylow=0, Double_t yup=0)
Definition: TPileupReweighting.cxx:93
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
CP::TPileupReweighting::GetPrimaryWeight
Float_t GetPrimaryWeight(Int_t periodNumber, Int_t channelNumber, Float_t x)
Definition: TPileupReweighting.cxx:1653
CP::TPileupReweighting::AddLumiCalcFile
Int_t AddLumiCalcFile(const TString &fileName, const TString &trigger="None")
Definition: TPileupReweighting.cxx:847
CP::TPileupReweighting::GetUnrepresentedDataFraction
Double_t GetUnrepresentedDataFraction(Int_t periodNumber, Int_t channel)
return the unrepresented data fraction in a given channel .
Definition: TPileupReweighting.cxx:1386
AddEmptyComponent.histName
string histName
Definition: AddEmptyComponent.py:64
plotting.yearwise_efficiency.channel
channel
Definition: yearwise_efficiency.py:24
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
CP::TPileupReweighting::GetPeriodWeight
Float_t GetPeriodWeight(Int_t periodNumber, Int_t channelNumber)
Definition: TPileupReweighting.cxx:1623
integral
double integral(TH1 *h)
Definition: computils.cxx:58
CP::TPileupReweighting::IsUnrepresentedData
Bool_t IsUnrepresentedData(Int_t runNumber, Float_t x, Float_t y=0.)
Definition: TPileupReweighting.cxx:1524
CP::TPileupReweighting::AddMetaDataFile
Int_t AddMetaDataFile(const TString &fileName, const TString &channelBranchName="mc_channel_number")
Definition: TPileupReweighting.cxx:506
CP::TPileupReweighting
Definition: TPileupReweighting.h:50
CP::TPileupReweighting::SetBinning
Int_t SetBinning(Int_t nbinsx, Double_t *xbins, Int_t nbinsy=0, Double_t *ybins=0)
Add a histogram binning config.
Definition: TPileupReweighting.cxx:83
plotmaker.hist
hist
Definition: plotmaker.py:148
StateLessPT_NewConfig.Format
Format
Definition: StateLessPT_NewConfig.py:146
mergePhysValFiles.start
start
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:14
CP::TPileupReweighting::GenerateMetaDataFile
Int_t GenerateMetaDataFile(const TString &fileName, const TString &channelBranchName="mc_channel_number")
Definition: TPileupReweighting.cxx:460
tree
TChain * tree
Definition: tile_monitor.h:30
CP::TPileupReweighting::GetIntegratedLumi
Double_t GetIntegratedLumi(const TString &trigger="")
total luminosity loaded and accepted by the tool (in inverse pb)
Definition: TPileupReweighting.cxx:115
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
CP::TPileupReweighting::SetDefaultChannel
void SetDefaultChannel(Int_t channel, Int_t mcRunNumber=-1)
Set which channel should be used as a default when specific mc channel distributions cannot be found.
Definition: TPileupReweighting.cxx:72
skel.it
it
Definition: skel.GENtoEVGEN.py:396
CP::TPileupReweighting::Merge
Int_t Merge(TCollection *coll)
Definition: TPileupReweighting.cxx:2039
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
CP::TPileupReweighting::UsePeriodConfig
Int_t UsePeriodConfig(const TString &configName)
use a hardcoded period configuration
Definition: TPileupReweighting.cxx:234
PixelAthClusterMonAlgCfg.ybins
ybins
Definition: PixelAthClusterMonAlgCfg.py:169
run
int run(int argc, char *argv[])
Definition: ttree2hdf5.cxx:28
bin
Definition: BinsDiffFromStripMedian.h:43
CP::TPileupReweighting::MakeWeightTree
Bool_t MakeWeightTree(TString channelNumbers, TString outFile, TString hashBranch="PRWHash", TString weightBranch="PileupWeight")
Definition: TPileupReweighting.cxx:1530
athena.value
value
Definition: athena.py:124
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
RunTileMonitoring.WriteToFile
WriteToFile
Definition: RunTileMonitoring.py:428
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
dbg::ptr
void * ptr(T *p)
Definition: SGImplSvc.cxx:74
skel.runNum
runNum
Definition: skel.ABtoEVGEN.py:137
CP::TPileupReweighting::GetDataWeight
Double_t GetDataWeight(Int_t runNumber, const TString &trigger, Double_t x, bool run_dependent=false)
Method for weighting data to account for prescales and mu bias.
Definition: TPileupReweighting.cxx:1730
athena.ps1
ps1
rename ourselfs to athena, both the prompt and the process (for top & ps)
Definition: athena.py:135
CP::TPileupReweighting::RemoveChannel
Bool_t RemoveChannel(int chanNum)
Removes a channel from the inputs ...
Definition: TPileupReweighting.cxx:1001
x
#define x
dq_defect_copy_defect_database.channels
def channels
Definition: dq_defect_copy_defect_database.py:56
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:93
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
CP::TPileupReweighting::GetRandomRunNumber
UInt_t GetRandomRunNumber(Int_t mcRunNumber)
Gets a random data run number according to the integrated lumi distribution associated to this mcRunN...
Definition: TPileupReweighting.cxx:1412
CP::TPileupReweighting::GetRandomLumiBlockNumber
UInt_t GetRandomLumiBlockNumber(UInt_t runNumber)
Get a random lumi block from the run number given.
Definition: TPileupReweighting.cxx:1469
CP::TPileupReweighting::CompositeTrigger
Definition: TPileupReweighting.h:366
TPileupReweighting.h
CP::TPileupReweighting::AddPeriod
Int_t AddPeriod(Int_t periodNumber, UInt_t start, UInt_t end)
Assign an mc RunNumber to a data period.
Definition: TPileupReweighting.cxx:368
fitman.tmpHist
def tmpHist(what, wmin=-1e10, wmax=+1e10)
Definition: fitman.py:146
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
FortranAlgorithmOptions.fileName
fileName
Definition: FortranAlgorithmOptions.py:13
athena.ps2
ps2
Definition: athena.py:136
CP::TPileupReweighting::RemapPeriod
void RemapPeriod(Int_t periodNumber1, Int_t periodNumber2)
Combine two period numbers.
Definition: TPileupReweighting.cxx:63
python.BunchSpacingUtils.lb
lb
Definition: BunchSpacingUtils.py:88
checkCoolLatestUpdate.chanNum
chanNum
Definition: checkCoolLatestUpdate.py:27
fillPileUpNoiseLumi.next
next
Definition: fillPileUpNoiseLumi.py:52
CP::TPileupReweighting::Period
Definition: TPileupReweighting.h:408
ParseInputs.gDirectory
gDirectory
Definition: Final2012/ParseInputs.py:133
lumiFormat.i
int i
Definition: lumiFormat.py:85
CP::TPileupReweighting::CalculatePrescaledLuminosityHistograms
void CalculatePrescaledLuminosityHistograms(const TString &trigger, int run_dependent=0)
Definition: TPileupReweighting.cxx:2086
beamspotman.n
n
Definition: beamspotman.py:731
checkxAOD.frac
frac
Definition: Tools/PyUtils/bin/checkxAOD.py:259
ChangeHistoRange.binsY
list binsY
Definition: ChangeHistoRange.py:59
mergePhysValFiles.origDir
origDir
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:24
CP::TPileupReweighting::m_periodList
std::list< Period > m_periodList
List physically holding (owning) period objects.
Definition: TPileupReweighting.h:450
CP::TPileupReweighting::GetSecondaryWeight
Float_t GetSecondaryWeight(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y)
Definition: TPileupReweighting.cxx:1701
file
TFile * file
Definition: tile_monitor.h:29
CP::TPileupReweighting::GetCombinedWeight
Float_t GetCombinedWeight(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y=0.)
Definition: TPileupReweighting.cxx:1603
CP::TPileupReweighting::WriteToFile
Int_t WriteToFile(const TString &filename="")
Definition: TPileupReweighting.cxx:1888
CP::TPileupReweighting::m_periods
std::map< Int_t, Period * > m_periods
Definition: TPileupReweighting.h:451
CP::TPileupReweighting::GetPrescaleWeight
Double_t GetPrescaleWeight(Int_t runNumber, const TString &trigger, Double_t x, bool run_dependent=false)
Method for prescaling MC to account for prescales in data.
Definition: TPileupReweighting.cxx:1825
makeComparison.rootFile
rootFile
Definition: makeComparison.py:27
ClassImp
ClassImp(CP::TPileupReweighting) CP
Definition: TPileupReweighting.cxx:41
run
Definition: run.py:1
CP::TPileupReweighting::calculateHistograms
void calculateHistograms(CompositeTrigger *trigger, int run_dependent)
Definition: TPileupReweighting.cxx:2109
CP::TPileupReweighting::normalizeHistogram
void normalizeHistogram(TH1 *histo)
Normalize histograms.
Definition: TPileupReweighting.cxx:1996
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
CP::TPileupReweighting::GetNearestGoodBin
Int_t GetNearestGoodBin(Int_t thisMCRunNumber, Int_t bin)
Definition: TPileupReweighting.cxx:1367
make_coralServer_rep.proj
proj
Definition: make_coralServer_rep.py:48
CP::TPileupReweighting::m_runs
std::map< UInt_t, Run > m_runs
Definition: TPileupReweighting.h:452
DQPostProcessTest.outFile
outFile
Comment Out Those You do not wish to run.
Definition: DQPostProcessTest.py:37
TestSUSYToolsAlg.outName
string outName
Definition: TestSUSYToolsAlg.py:173
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:83
jobOptions.Initialize
Initialize
Definition: jobOptions.pA.py:28
CP::TPileupReweighting::makeTrigger
std::unique_ptr< CompositeTrigger > makeTrigger(const TString &s)
Definition: TPileupReweighting.cxx:2278
LArCellBinning.xbins
int xbins
Definition: LArCellBinning.py:163
CP::TPileupReweighting::AddConfigFile
Int_t AddConfigFile(const TString &fileName)
Definition: TPileupReweighting.cxx:936
python.AthDsoLogger.fname
string fname
Definition: AthDsoLogger.py:67
DeMoAtlasDataLoss.runNumber
string runNumber
Definition: DeMoAtlasDataLoss.py:64
CP::TPileupReweighting::GetIntegratedLumiFraction
Double_t GetIntegratedLumiFraction(Int_t periodNumber, UInt_t start, UInt_t end)
return fraction of lumi assigned to periodNumber (or mcRunNumber) that is between start and end data ...
Definition: TPileupReweighting.cxx:182
RTTAlgmain.branch
branch
Definition: RTTAlgmain.py:61
CP::TPileupReweighting::Initialize
Int_t Initialize()
Initialize this class once before the event loop starts If distribution information is provided,...
Definition: TPileupReweighting.cxx:1021
lumiFormat.lumi
lumi
Definition: lumiFormat.py:106
y
#define y
h
CP::TPileupReweighting::GetPRWHash
ULong64_t GetPRWHash(Int_t periodNumber, Int_t channelNumber, Float_t x, Float_t y=0.)
Definition: TPileupReweighting.cxx:1584
CondAlgsOpts.found
int found
Definition: CondAlgsOpts.py:101
CP::TPileupReweighting::GetMetaDataTree
TTree * GetMetaDataTree()
combines loaded metadata with channel sumsofweights and entry counts
Definition: TPileupReweighting.cxx:574
CP::TPileupReweighting::GetDefaultChannel
Int_t GetDefaultChannel(Int_t mcRunNumber=-1)
Definition: TPileupReweighting.cxx:109
DeMoScan.first
bool first
Definition: DeMoScan.py:536
CP::TPileupReweighting::CloneEmptyHistogram
std::unique_ptr< TH1 > CloneEmptyHistogram(Int_t runNumber, Int_t channelNumber)
Definition: TPileupReweighting.cxx:435
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:24
EventInfoRead.isMC
isMC
Definition: EventInfoRead.py:11
CP::TPileupReweighting::GetFirstFoundPeriodNumber
Int_t GetFirstFoundPeriodNumber(UInt_t runNumber)
Get the first period number with the data run number contained - assume all periods are disconnected ...
Definition: TPileupReweighting.cxx:422
ChangeHistoRange.binsX
list binsX
Definition: ChangeHistoRange.py:56
xAOD::JetConstituentVector::iterator
Definition: JetConstituentVector.h:121
CP::TPileupReweighting::Run
Definition: TPileupReweighting.h:439
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
CP::TPileupReweighting::GetLumiBlockMu
Float_t GetLumiBlockMu(Int_t runNumber, UInt_t lb)
get the lumiblock mu, useful for 'updating' the mu coming from data to account for new lumitags
Definition: TPileupReweighting.cxx:168
CP::TPileupReweighting::GetLumiBlockIntegratedLumi
Double_t GetLumiBlockIntegratedLumi(Int_t runNumber, UInt_t lb)
get integrated lumi for specific run and lumiblock number .
Definition: TPileupReweighting.cxx:154
plotBeamSpotCompare.histo
histo
Definition: plotBeamSpotCompare.py:415
CP::TPileupReweighting::GetRandomPeriodNumber
Int_t GetRandomPeriodNumber(Int_t mcRunNumber)
Get random period number from the sub-periods assigned to this run number.
Definition: TPileupReweighting.cxx:1492
python.PyAthena.obj
obj
Definition: PyAthena.py:132
dqBeamSpot.nEntries
int nEntries
Definition: dqBeamSpot.py:73
CaloNoise_fillDB.mu
mu
Definition: CaloNoise_fillDB.py:53
TileDCSDataPlotter.tx
tx
Definition: TileDCSDataPlotter.py:878
CP::TPileupReweighting::TPileupReweighting
TPileupReweighting(const char *name="TPileupReweighting")
Standard constructor.
read_hist_ntuple.f1
f1
Definition: read_hist_ntuple.py:4