UnfoldXSec.C
Go to the documentation of this file.
1 ////////////////////////////////////////
2 /// \brief Remake of script to unfold and calculate cross-section
3 /// \author Connor Johnson
4 /// \date July 2019
5 ////////////////////////////////////////
6 
7 /// CAFAna includes
8 #include "CAFAna/Vars/Vars.h"
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Core/Ratio.h"
17 #include "CAFAna/Core/HistCache.h"
18 
19 /// NumuCCInc Analysis includes
25 
26 /// ROOT includes
27 #include "TFile.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TH3.h"
31 #include "TUnfold.h"
32 #include "TUnfoldDensity.h"
33 #include "TUnfoldBinning.h"
34 
35 // RooUnfold Method
36 #include "RooUnfoldResponse.h"
37 #include "RooUnfoldBayes.h"
38 // #include "RooUnfoldSvd.h"
39 
40 /// Standard includes
41 #include <iostream>
42 #include <vector>
43 #include <string>
44 #include <map>
45 #include <exception>
46 #include <chrono>
47 
48 using std::string; using std::map; using std::vector; using std::cout; using std::endl; //using std::chrono;
49 
50 ////////////////////////////////////////
51 /// Global constant settings
52 ////////////////////////////////////////
53 const double pot = 8.09E20; /// Scale to the FHC protons-on-target
54 const double ntargs = 3.89137e+31; /// Number of targets, calcuted by ???
55 
56 /// Set whether to use ppfx-weighted or unweighted fake-data
57 const string recoLabel = "reco";
58 const string trueLabel = "true";
59 const string ppfxXsecWgtLabel = "ppfxxsecwgt";
60 const string ppfxWgtLabel = "ppfxwgt";
61 const string unwgtLabel = "unwgt";
62 const string xs_label = ppfxXsecWgtLabel;
63 
64 /// Directory of flux information in flux file
65 const string fluxDirName = "fluxes/CV";
66 
67 ////////////////////////////////////////
68 /// Unfolding variables
69 ////////////////////////////////////////
70 
71 /// 3D Unfolding variables. Need to include a mapping to the x,y,z binnings
72 map<string, vector<ana::Binning > > unfoldVars3D{
74 };
75 
76 /// 1D Unfolding variables
77 const vector<string> unfoldVars1D{
78  "ENu",
79  "Q2"
80 };
81 
82 /// Map the variable names to their folder in the fake data (could preferably switch the labelling of the fake data instead)
83 map<string, string> varToFakeData{
84  {"EAvail", "fakedata"},
85  {"ENu", "fakedata_enu"},
86  {"Q2", "fakedata_q2"}
87 };
88 
89 map<string, string> varToRealData{
90  {"EAvail", "data_eavail"},
91  {"ENu", "data_enu"},
92  {"Q2", "data_q2"}
93 };
94 
95 /// Standard systmatics to check
96 const vector<string> standardSysts =
97 {
98  "cv",
99  "lightup",
100  "lightdown",
101  "calibpos",
102  "calibneg",
103  "calibshift",
104  "muoneup",
105  "muonedw",
106  "ckv",
107  "angleup",
108  "angledw",
109  "neutron_Up",
110  "neutron_Dw",
111  "neutron_nom",
112  "nom_good_seventh"
113 };
114 
115 /// Beam systematics (each has _Up and _Dw variations)
116 const vector<string> beamSysts = {
117  "2kA",
118  "02mmBeamSpotSize",
119  "1mmBeamShiftX",
120  "1mmBeamShiftY",
121  "3mmHorn1X",
122  "3mmHorn1Y",
123  "3mmHorn2X",
124  "3mmHorn2Y",
125  "7mmTargetZ",
126  "MagneticFieldinDecayPipe",
127  "1mmHornWater"
128 };
129 
130 // Number of iterations
131 const unsigned int DEFAULT_ITERATIONS = 3;
133 
134 // Number of GENIE / PPFX
135 const int NUM_GENIE_UNIV = 1000;
136 const int NUM_PPFX_UNIV = 100;
137 
138 // Unfold using RooUnfoldResponse
139 TH1* UnfoldMethodRooUnfold(TH1* measBinsHist, TH1* trueBinsHist, TH2* unfoldingMatrix, TH1* dataMinusBkgdHist, TMatrixD*& covMatrix)
140 {
141  // Set up response matrix
142  RooUnfoldResponse * unfoldResponse = new RooUnfoldResponse(measBinsHist, trueBinsHist, unfoldingMatrix);
143  RooUnfoldBayes * unfoldIterative = new RooUnfoldBayes(unfoldResponse, dataMinusBkgdHist, nIterations);
144  chrono::steady_clock::time_point time0 = chrono::steady_clock::now();
145  TH1* unfoldResult = unfoldIterative->Hreco(RooUnfold::ErrorTreatment::kCovariance);
146  chrono::steady_clock::time_point time1 = chrono::steady_clock::now();
147  float duration = chrono::duration_cast<chrono::seconds>(time1 - time0).count();
148  cout << "\tUnfolding at " << nIterations << " iter took: " << duration << " sec.\n"
149  << "\tAverage time per iteration = " << to_string(duration / float(nIterations)) << " sec" << endl;
150  covMatrix = new TMatrixD(unfoldIterative->Ereco(RooUnfold::ErrorTreatment::kCovariance));
151  return unfoldResult;
152 }
153 
154 // Unfolding wrapper in 3D
155 bool Unfold3D(string systName,
156  string varName,
157  TDirectory * varDir,
158  TDirectory * dataDir,
159  TFile * outFile)
160 {
161  cout<<"Unfolding Var : "<< varName << endl;
162  /// Set up all input directories
163  TDirectory * mcBkgdEstDir = varDir->GetDirectory("BkgdEst", true);
164  TDirectory * unfoldMatrixDir = varDir->GetDirectory("RecoTrue", true);
165  TDirectory * mcSelRecoDir = varDir->GetDirectory("MCSel_Reco", true);
166  TDirectory * mcSigRecoDir = varDir->GetDirectory("MCSig_Reco", true);
167  TDirectory * mcSelSigDir = varDir->GetDirectory("MCSig", true);
168  TDirectory * mcSigTrueDir = varDir->GetDirectory("MCSigNuTrue",true);
169 
170  /// Load the input directories as cafe objects
171  ana::TrivialBkgdEstimator * mcBkgdEst = ana::TrivialBkgdEstimator::LoadFrom(mcBkgdEstDir ).release();
172  ana::ReweightableSpectrum * unfoldMatrix = ana::ReweightableSpectrum::LoadFrom(unfoldMatrixDir).release();
173  ana::Spectrum * mcSelRecoSpec = ana::Spectrum::LoadFrom( mcSelRecoDir ).release();
174  ana::Spectrum * mcSigRecoSpec = ana::Spectrum::LoadFrom( mcSigRecoDir ).release();
175  ana::Spectrum * mcSelSigSpec = ana::Spectrum::LoadFrom( mcSelSigDir ).release();
176  ana::Spectrum * mcSigTrueSpec = ana::Spectrum::LoadFrom( mcSigTrueDir ).release();
177  ana::Spectrum * dataSpec = ana::Spectrum::LoadFrom( dataDir ).release();
178 
179  /// Multiply the fake or real data by the purity
180  ana::Spectrum * dataMinusBkgdSpec = new ana::Spectrum(*dataSpec - mcBkgdEst->Background());
181  ana::Ratio * purityRatio = new ana::Ratio(*mcSigRecoSpec, *mcSelRecoSpec, true); // Purity, errors calculate differently.
182  ana::Spectrum * dataTimesPurity = new ana::Spectrum((*dataSpec) * (*purityRatio));
183 
184  // Calculate efficiency
185  ana::Ratio * efficiencyRatio = new ana::Ratio(*mcSelSigSpec, *mcSigTrueSpec, true);
186 
187  /// Grab the 3D binning info (3D vars only)
188  vector<ana::Binning>& bins3D = unfoldVars3D[varName];
189  ana::Binning& binsX = bins3D[0];
190  ana::Binning& binsY = bins3D[1];
191  ana::Binning& binsZ = bins3D[2];
192 
193  /// Convert to ROOT variables
194  TH2 * unfoldingHist = unfoldMatrix->ToTH2(pot);
195  TH1 * dataMinusBkgdHist = dataMinusBkgdSpec->ToTH1(pot);// ana::ToTH3(*dataMinusBkgdSpec, pot, ana::kPOT, binsX, binsY, binsZ);
196  TH1 * dataTimesPurHist = dataTimesPurity->ToTH1(pot);
197  TH1 * purityHist = purityRatio->ToTH1();
198  TH1 * efficiencyHist = efficiencyRatio->ToTH1();
199  TH1 * mcSelSigHist = mcSelSigSpec->ToTH1(pot);
200  TH1 * mcSigTrueHist = mcSigTrueSpec->ToTH1(pot);
201  TH1 * mcSigRecoHist = mcSigRecoSpec->ToTH1(pot);
202  TH1 * mcSelRecoHist = mcSelRecoSpec->ToTH1(pot);
203  unfoldingHist-> SetName("UnfoldingMatrix");
204  dataMinusBkgdHist->SetName("DataMinusBkgd");
205  dataTimesPurHist-> SetName("DataTimesPur");
206  purityHist-> SetName("Purity");
207  efficiencyHist-> SetName("Efficiency");
208  mcSelSigHist-> SetName("SelSig_True");
209  mcSigTrueHist-> SetName("Signal_True");
210  mcSigRecoDir-> SetName("Signal_Reco");
211  mcSelRecoDir-> SetName("Select_Reco");
212 
213  // Unfolding
214  // TH3D* meanEmptyBins = new TH3D((systName + "_" + varName + "_" "meanBinning").c_str(), "meanBinning", binsX.NBins(), binsX.Edges().data(), binsY.NBins(), binsY.Edges().data(), binsZ.NBins(), binsZ.Edges().data());
215  // TH3D* trueEmptyBins = new TH3D((systName + "_" + varName + "_" "trueBinning").c_str(), "trueBinning", binsX.NBins(), binsX.Edges().data(), binsY.NBins(), binsY.Edges().data(), binsZ.NBins(), binsZ.Edges().data());
216  TMatrixD * covMatrix = NULL;
217  std::unique_ptr<TH1> unfoldedResult(std::move(UnfoldMethodRooUnfold(0, 0, unfoldingHist, dataTimesPurHist, covMatrix)));
218  TH3 * unfolded3D = ana::ToTH3Helper(std::move(unfoldedResult), binsX, binsY, binsZ, ana::kBinContent);
219  unfolded3D->SetName("UnfoldedData");
220 
221  // Outputting
222  TDirectory * systDir = outFile->GetDirectory(systName.c_str());
223  if (!systDir || systDir == NULL)
224  systDir = outFile->mkdir(systName.c_str());
225  TDirectory * varOutDir = systDir->mkdir(varName.c_str());
226  varOutDir->cd();
227 
228  unfoldingHist-> Write();
229  dataMinusBkgdHist->Write();
230  unfolded3D-> Write();
231  covMatrix-> Write("Covariance");
232  dataTimesPurHist-> Write();
233  purityHist-> Write();
234  efficiencyHist-> Write();
235  mcSelSigHist-> Write();
236  mcSigTrueHist-> Write();
237  mcSigRecoDir-> Write();
238  mcSelRecoDir-> Write();
239 
240  cout << "-----------------" << endl;
241  return false;
242 }
243 
244 // Unfolding wrapper in 1D
245 bool Unfold1D(string systName,
246  string varName,
247  TDirectory * varDir,
248  TDirectory * dataDir,
249  TFile * outFile)
250 {
251  cout<<"Unfolding Var : "<< varName << endl;
252  /// Set up all input directories
253  /// Set up all input directories
254  TDirectory * mcBkgdEstDir = varDir->GetDirectory("BkgdEst", true);
255  TDirectory * unfoldMatrixDir = varDir->GetDirectory("RecoTrue", true);
256  TDirectory * mcSelRecoDir = varDir->GetDirectory("MCSel_Reco", true);
257  TDirectory * mcSigRecoDir = varDir->GetDirectory("MCSig_Reco", true);
258  TDirectory * mcSelSigDir = varDir->GetDirectory("MCSig", true);
259  TDirectory * mcSigTrueDir = varDir->GetDirectory("MCSigNuTrue",true);
260 
261  /// Load the input directories as cafe objects
262  ana::TrivialBkgdEstimator * mcBkgdEst = ana::TrivialBkgdEstimator::LoadFrom(mcBkgdEstDir ).release();
263  ana::ReweightableSpectrum * unfoldMatrix = ana::ReweightableSpectrum::LoadFrom(unfoldMatrixDir).release();
264  ana::Spectrum * mcSelRecoSpec = ana::Spectrum::LoadFrom( mcSelRecoDir ).release();
265  ana::Spectrum * mcSigRecoSpec = ana::Spectrum::LoadFrom( mcSigRecoDir ).release();
266  ana::Spectrum * mcSelSigSpec = ana::Spectrum::LoadFrom( mcSelSigDir ).release();
267  ana::Spectrum * mcSigTrueSpec = ana::Spectrum::LoadFrom( mcSigTrueDir ).release();
268  ana::Spectrum * dataSpec = ana::Spectrum::LoadFrom( dataDir ).release();
269 
270  /// Multiply the fake or real data by the purity
271  ana::Spectrum * dataMinusBkgdSpec = new ana::Spectrum(*dataSpec - mcBkgdEst->Background());
272  ana::Ratio * purityRatio = new ana::Ratio(*mcSigRecoSpec, *mcSelRecoSpec, true); // Purity, errors calculate differently.
273  ana::Spectrum * dataTimesPurity = new ana::Spectrum((*dataSpec) * (*purityRatio));
274 
275  // Calculate efficiency
276  ana::Ratio * efficiencyRatio = new ana::Ratio(*mcSelSigSpec, *mcSigTrueSpec, true);
277 
278  /// Convert to ROOT variables
279  TH2 * unfoldingHist = unfoldMatrix->ToTH2(pot);
280  TH1 * dataMinusBkgdHist = dataMinusBkgdSpec->ToTH1(pot);// ana::ToTH3(*dataMinusBkgdSpec, pot, ana::kPOT, binsX, binsY, binsZ);
281  TH1 * dataTimesPurHist = dataTimesPurity->ToTH1(pot);
282  TH1 * purityHist = purityRatio->ToTH1();
283  TH1 * efficiencyHist = efficiencyRatio->ToTH1();
284  TH1 * mcSelSigHist = mcSelSigSpec->ToTH1(pot);
285  TH1 * mcSigTrueHist = mcSigTrueSpec->ToTH1(pot);
286  TH1 * mcSigRecoHist = mcSigRecoSpec->ToTH1(pot);
287  TH1 * mcSelRecoHist = mcSelRecoSpec->ToTH1(pot);
288  unfoldingHist-> SetName("UnfoldingMatrix");
289  dataMinusBkgdHist->SetName("DataMinusBkgd");
290  dataTimesPurHist-> SetName("DataTimesPur");
291  purityHist-> SetName("Purity");
292  efficiencyHist-> SetName("Efficiency");
293  mcSelSigHist-> SetName("SelSig_True");
294  mcSigTrueHist-> SetName("Signal_True");
295  mcSigRecoDir-> SetName("Signal_Reco");
296  mcSelRecoDir-> SetName("Select_Reco");
297 
298  // Unfolding
299  TMatrixD * covMatrix = NULL;
300  TH1 * unfoldedResult = (TH1*) UnfoldMethodRooUnfold(0, 0, unfoldingHist, dataMinusBkgdHist, covMatrix);
301  unfoldedResult->SetName("UnfoldedData");
302  TH2D* covHist = new TH2D(*covMatrix);
303  covHist->SetName("Covariance");
304 
305  // Outputting
306  TDirectory * systDir = outFile->GetDirectory(systName.c_str());
307  if (!systDir || systDir == NULL)
308  systDir = outFile->mkdir(systName.c_str());
309  TDirectory * varOutDir = systDir->mkdir(varName.c_str());
310  varOutDir->cd();
311 
312  unfoldingHist-> Write();
313  dataMinusBkgdHist->Write();
314  unfoldedResult-> Write();
315  covHist-> Write();
316  dataTimesPurHist-> Write();
317  purityHist-> Write();
318  efficiencyHist-> Write();
319 
320  cout << "-----------------" << endl;
321  return false;
322 }
323 
324 // Loop through all vars and unfold them.
325 void UnfoldAllVars( string systName,
326  TDirectory * systDir,
327  map<string, TDirectory*>& dataDirs,
328  TFile * outFile)
329 {
330  // Loop through 3D vars
331  for (const pair<string, vector<ana::Binning > >& varInfo : unfoldVars3D)
332  {
333  string varName = varInfo.first;
334  TDirectory * varDir = systDir->GetDirectory(varName.c_str());
335  if (!varDir || varDir == NULL)
336  continue;
337  TDirectory * dataDir = dataDirs[varName];
338  Unfold3D(systName, varName, varDir, dataDir, outFile);
339  ana::HistCache::ClearCache();
340  }
341 
342  // Loop through 1D vars
343  for (string varName : unfoldVars1D)
344  {
345  TDirectory * varDir = systDir->GetDirectory(varName.c_str());
346  if (!varDir || varDir == NULL)
347  continue;
348  TDirectory * dataDir = dataDirs[varName];
349  Unfold1D(systName, varName, varDir, dataDir, outFile);
350  ana::HistCache::ClearCache();
351  }
352 }
353 
354 // Looks through an input file for standard systematics to unfold
356  map<string, TDirectory*>& dataDirs,
357  TFile * outFile)
358 {
359  // Loop through the standard systs
360  for (string syst : standardSysts)
361  {
362  TDirectory * systDir = inputFile->GetDirectory(syst.c_str());
363  if (systDir && systDir != NULL){
364  cout << "\n~~~~~~~~~~~~~~~~~~~\n~~~~~~~~~~~~~~~~~~~\nUnfolding systematic: " << syst << endl;
365  UnfoldAllVars(syst, systDir, dataDirs, outFile);
366  }
367  }
368 }
369 
370 // Looks through an input file for beam systematics to unfold
371 void UnfoldBeamSysts( TFile * inputFile,
372  map<string, TDirectory*>& dataDirs,
373  TFile * outFile)
374 {
375  // Loop through the standard systs
376  for (string syst : beamSysts)
377  {
378  string systUp = syst + "_Up";
379  string systDw = syst + "_Dw";
380 
381  TDirectory * systUpDir = inputFile->GetDirectory(systUp.c_str());
382  if (systUpDir && systUpDir != NULL){
383  cout << "Unfolding systematic: " << syst << endl;
384  UnfoldAllVars(systUp, systUpDir, dataDirs, outFile);
385  }
386  else
387  continue;
388 
389  TDirectory * systDwDir = inputFile->GetDirectory(systDw.c_str());
390  if (systDwDir && systDwDir != NULL)
391  UnfoldAllVars(systDw, systDwDir, dataDirs, outFile);
392  else cerr << "Error: Missing 'Dw' information on systematic " << syst << ". Skipping for now, fix this." << endl;
393  }
394 }
395 
396 // Looks through an input file for genie universes to unfold
397 void UnfoldGenie( TFile * inputFile,
398  map<string, TDirectory*>& dataDirs,
399  TFile * outFile,
400  const int firstGenie,
401  const int lastGenie)
402 {
403  if (firstGenie < 0 || lastGenie < firstGenie)
404  return;
405 
406  for (int iuniv = firstGenie; iuniv <= lastGenie && iuniv < NUM_GENIE_UNIV; ++iuniv)
407  {
408  string syst = "genie" + to_string(iuniv);
409 
410  TDirectory * systDir = inputFile->GetDirectory(syst.c_str());
411  if (systDir && systDir != NULL){
412  cout << "Unfolding systematic: " << syst << endl;
413  UnfoldAllVars(syst, systDir, dataDirs, outFile);
414  }
415  }
416 }
417 
418 // Looks through an input file for genie universes to unfold
419 void UnfoldPPFX( TFile * inputFile,
420  map<string, TDirectory*>& dataDirs,
421  TFile * outFile,
422  const int firstPPFX,
423  const int lastPPFX)
424 {
425  if (firstPPFX < 0 || lastPPFX < firstPPFX)
426  return;
427 
428  for (int iuniv = firstPPFX; iuniv <= lastPPFX && iuniv < NUM_PPFX_UNIV; ++iuniv)
429  {
430  string syst = "ppfx" + to_string(iuniv);
431 
432  TDirectory * systDir = inputFile->GetDirectory(syst.c_str());
433  if (systDir && systDir != NULL){
434  cout << "Unfolding systematic: " << syst << endl;
435  UnfoldAllVars(syst, systDir, dataDirs, outFile);
436  }
437  }
438 }
439 
441  string dataFile,
442  string systFileName,
443  string outputFile,
444  int firstGenie = -1,
445  int lastGenie = -1,
446  int firstPPFX = -1,
447  int lastPPFX = -1,
448  const int overrideIter = -1,
449  const unsigned int nJobs = 1,
450  const bool useFakeData = false)
451 {
452  // NJobs > 1 implies grid.
453  int process = -1;
454  if (nJobs > 1){
455  if (!getenv("PROCESS")){
456  cerr << "Number of jobs given > 1, but no $PROCESS found. Aborting..." << endl;
457  return;
458  }
459  process = atoi(getenv("PROCESS"));
460  if (process >= (int) nJobs){
461  cerr << "Process found greater than expected number of jobs. Exiting." << endl;
462  return;
463  }
464 
465  // Compute first and last GENIE for job
466  int geniePerJob = std::ceil(float(lastGenie - firstGenie + 1.0) / float(nJobs));
467  if (firstGenie >= 0){
468  lastGenie = std::min(firstGenie + (process * geniePerJob) + geniePerJob - 1, lastGenie);
469  firstGenie = firstGenie + (process * geniePerJob);
470  cout << "With $PROCESS:" << process << " of nJobs:" << nJobs << ", will run " << to_string(geniePerJob) << " GENIE.\n\tReset to run univs " << firstGenie << "-" << lastGenie << endl;
471  }
472 
473  // Compute first and last PPFX for job
474  int ppfxPerJob = std::ceil(float(lastPPFX - firstPPFX + 1.0) / float(nJobs));
475  if (firstPPFX >= 0){
476  lastPPFX = std::min(firstPPFX + (process * ppfxPerJob) + ppfxPerJob - 1, lastPPFX);
477  firstPPFX = firstPPFX + (process * ppfxPerJob);
478  cout << "With $PROCESS:" << process << " of nJobs:" << nJobs << ", will run " << to_string(ppfxPerJob) << " PPFX.\n\tReset to run univs " << firstPPFX << "-" << lastPPFX << endl;
479  }
480  }
481 
482  // Set the global iteration count
483  if (overrideIter >= 0)
484  ::nIterations = overrideIter;
485 
486  // Open input and output files, ensure they opened correctly
487  // TFile * systFile = new TFile((inputDir + coreDefinitions["systFile"]).c_str(), "READ");
488  TFile * systFile = TFile::Open(systFileName.c_str(), "READ");
489  // TFile * beamFile = new TFile((inputDir + coreDefinitions["beamFile"]).c_str(), "READ");
490  TFile * realDataFluxFile = TFile::Open(dataFile.c_str(), "READ");
491  TFile * outFile = new TFile(outputFile.c_str(), "RECREATE");
492  assert(systFile && systFile != NULL && !systFile->IsZombie() && systFile->IsOpen());
493  assert(realDataFluxFile && realDataFluxFile != NULL && !realDataFluxFile->IsZombie() && realDataFluxFile->IsOpen());
494  assert(outFile && outFile != NULL && !outFile->IsZombie());
495  outFile->cd();
496 
497  // Persistent Directory, Spectrum, Ratio storage (1 member for each 3D or 1D variable)
498  map<string, TDirectory*> dataDirs; // Real or fake data
499  // map<string, ana::Ratio*> mcInvPurityRatios;
500 
501  // For each 3D variable, store the real or fake data
502  map<string, string>& varToData = useFakeData ? varToFakeData : varToRealData;
503  for(const pair<string, vector<ana::Binning > >& varInfo : unfoldVars3D)
504  {
505  string varName = varInfo.first;
506  dataDirs[varName] = realDataFluxFile->GetDirectory((varToData[varName] + "/" + recoLabel + "/" + xs_label).c_str(), true);
507  }
508 
509  // For each 1D variable, do the same
510  for(string varName : unfoldVars1D)
511  dataDirs[varName] = realDataFluxFile->GetDirectory((varToData[varName] + "/" + recoLabel + "/" + xs_label).c_str(), true);
512 
513  // Run unfolding on the various types of systematics
514  UnfoldStandardSysts(systFile, dataDirs, outFile);
515  UnfoldBeamSysts(systFile, dataDirs, outFile);
516  UnfoldGenie(systFile, dataDirs, outFile, firstGenie, lastGenie);
517  UnfoldPPFX(systFile, dataDirs, outFile, firstPPFX, lastPPFX);
518 }
TFile * inputFile
Definition: PlotXSec.C:134
void UnfoldAllVars(string systName, TDirectory *systDir, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:325
float covMatrix[xbins][xbins]
Definition: MakePlots.C:95
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
const double ntargs
Scale to the FHC protons-on-target.
Definition: UnfoldXSec.C:54
unsigned int nIterations
Definition: UnfoldXSec.C:132
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
map< string, string > varToRealData
Definition: UnfoldXSec.C:89
bool Unfold3D(string systName, string varName, TDirectory *varDir, TDirectory *dataDir, TFile *outFile)
Definition: UnfoldXSec.C:155
TH3 * ToTH3Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Helper for ana::ToTH3.
Definition: UtilsExt.cxx:186
const string fluxDirName
Directory of flux information in flux file.
Definition: UnfoldXSec.C:65
static std::unique_ptr< TrivialBkgdEstimator > LoadFrom(TDirectory *dir, const std::string &name)
const double pot
Global constant settings.
Definition: UnfoldXSec.C:53
void UnfoldXSec(string dataFile, string systFileName, string outputFile, int firstGenie=-1, int lastGenie=-1, int firstPPFX=-1, int lastPPFX=-1, const int overrideIter=-1, const unsigned int nJobs=1, const bool useFakeData=false)
Definition: UnfoldXSec.C:440
Spectrum with the value of a second variable, allowing for reweighting
OStream cerr
Definition: OStream.cxx:7
void UnfoldBeamSysts(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:371
const Binning eavailbins
const string trueLabel
Definition: UnfoldXSec.C:58
static std::unique_ptr< ReweightableSpectrum > LoadFrom(TDirectory *dir, const std::string &name)
const int NUM_PPFX_UNIV
Definition: UnfoldXSec.C:136
void UnfoldPPFX(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile, const int firstPPFX, const int lastPPFX)
Definition: UnfoldXSec.C:419
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
::xsd::cxx::tree::duration< char, simple_type > duration
Definition: Database.h:188
TFile * outFile
Definition: PlotXSec.C:135
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const vector< string > unfoldVars1D
1D Unfolding variables
Definition: UnfoldXSec.C:77
std::string getenv(std::string const &name)
const string recoLabel
Number of targets, calcuted by ???
Definition: UnfoldXSec.C:57
Spectrum Background() const override
TH1 * UnfoldMethodRooUnfold(TH1 *measBinsHist, TH1 *trueBinsHist, TH2 *unfoldingMatrix, TH1 *dataMinusBkgdHist, TMatrixD *&covMatrix)
Definition: UnfoldXSec.C:139
void UnfoldStandardSysts(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:355
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
map< string, string > varToFakeData
Map the variable names to their folder in the fake data (could preferably switch the labelling of the...
Definition: UnfoldXSec.C:83
std::vector< float > Spectrum
Definition: Constants.h:610
const string unwgtLabel
Definition: UnfoldXSec.C:61
Regular histogram.
Definition: UtilsExt.h:30
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
map< string, vector< ana::Binning > > unfoldVars3D
Unfolding variables.
Definition: UnfoldXSec.C:72
const Binning mukebins
Definition: NumuCCIncBins.h:90
bool Unfold1D(string systName, string varName, TDirectory *varDir, TDirectory *dataDir, TFile *outFile)
Definition: UnfoldXSec.C:245
const string xs_label
Definition: UnfoldXSec.C:62
assert(nhit_max >=nhit_nbins)
const vector< string > beamSysts
Beam systematics (each has _Up and _Dw variations)
Definition: UnfoldXSec.C:116
const string ppfxWgtLabel
Definition: UnfoldXSec.C:60
const int NUM_GENIE_UNIV
Definition: UnfoldXSec.C:135
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T min(const caf::Proxy< T > &a, T b)
const string ppfxXsecWgtLabel
Definition: UnfoldXSec.C:59
void UnfoldGenie(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile, const int firstGenie, const int lastGenie)
Definition: UnfoldXSec.C:397
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72
const unsigned int DEFAULT_ITERATIONS
Definition: UnfoldXSec.C:131
Just return the MC prediction for the background.
const vector< string > standardSysts
Standard systmatics to check.
Definition: UnfoldXSec.C:96
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11
TH2D * ToTH2(double pot) const
gm Write()
enum BeamMode string