Functions | Variables
UnfoldXSec.C File Reference
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Ratio.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Unfold/UnfoldIterative.h"
#include "CAFAna/XSec/TrivialBkgdEstimator.h"
#include "CAFAna/Core/HistCache.h"
#include "NDAna/numucc_inc/NumuCCIncAnalysis.h"
#include "NDAna/numucc_inc/NumuCCIncBins.h"
#include "NDAna/numucc_inc/NumuCCIncVars.h"
#include "NDAna/numucc_inc/NumuCCIncCuts.h"
#include "NDAna/numucc_inc/NDXSecMuonPID.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TUnfold.h"
#include "TUnfoldDensity.h"
#include "TUnfoldBinning.h"
#include "RooUnfoldResponse.h"
#include "RooUnfoldBayes.h"
#include <iostream>
#include <vector>
#include <string>
#include <map>
#include <exception>
#include <chrono>

Go to the source code of this file.

Functions

TH1 * UnfoldMethodRooUnfold (TH1 *measBinsHist, TH1 *trueBinsHist, TH2 *unfoldingMatrix, TH1 *dataMinusBkgdHist, TMatrixD *&covMatrix)
 
bool Unfold3D (string systName, string varName, TDirectory *varDir, TDirectory *dataDir, TFile *outFile)
 
bool Unfold1D (string systName, string varName, TDirectory *varDir, TDirectory *dataDir, TFile *outFile)
 
void UnfoldAllVars (string systName, TDirectory *systDir, map< string, TDirectory * > &dataDirs, TFile *outFile)
 
void UnfoldStandardSysts (TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile)
 
void UnfoldBeamSysts (TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile)
 
void UnfoldGenie (TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile, const int firstGenie, const int lastGenie)
 
void UnfoldPPFX (TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile, const int firstPPFX, const int lastPPFX)
 
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)
 

Variables

const double pot = 8.09E20
 Global constant settings. More...
 
const double ntargs = 3.89137e+31
 Scale to the FHC protons-on-target. More...
 
const string recoLabel = "reco"
 Number of targets, calcuted by ??? More...
 
const string trueLabel = "true"
 
const string ppfxXsecWgtLabel = "ppfxxsecwgt"
 
const string ppfxWgtLabel = "ppfxwgt"
 
const string unwgtLabel = "unwgt"
 
const string xs_label = ppfxXsecWgtLabel
 
const string fluxDirName = "fluxes/CV"
 Directory of flux information in flux file. More...
 
map< string, vector< ana::Binning > > unfoldVars3D
 Unfolding variables. More...
 
const vector< stringunfoldVars1D
 1D Unfolding variables More...
 
map< string, stringvarToFakeData
 Map the variable names to their folder in the fake data (could preferably switch the labelling of the fake data instead) More...
 
map< string, stringvarToRealData
 
const vector< stringstandardSysts
 Standard systmatics to check. More...
 
const vector< stringbeamSysts
 Beam systematics (each has _Up and _Dw variations) More...
 
const unsigned int DEFAULT_ITERATIONS = 3
 
unsigned int nIterations = DEFAULT_ITERATIONS
 
const int NUM_GENIE_UNIV = 1000
 
const int NUM_PPFX_UNIV = 100
 

Function Documentation

bool Unfold1D ( string  systName,
string  varName,
TDirectory *  varDir,
TDirectory *  dataDir,
TFile *  outFile 
)

Set up all input directories Set up all input directories

Load the input directories as cafe objects

Multiply the fake or real data by the purity

Convert to ROOT variables

Definition at line 245 of file UnfoldXSec.C.

References ana::TrivialBkgdEstimator::Background(), om::cout, covMatrix, allTimeWatchdog::endl, ana::TrivialBkgdEstimator::LoadFrom(), ana::ReweightableSpectrum::LoadFrom(), ana::Spectrum::LoadFrom(), pot, ana::Ratio::ToTH1(), ana::Spectrum::ToTH1(), ana::ReweightableSpectrum::ToTH2(), UnfoldMethodRooUnfold(), and Write().

Referenced by UnfoldAllVars().

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 }
float covMatrix[xbins][xbins]
Definition: MakePlots.C:95
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
static std::unique_ptr< TrivialBkgdEstimator > LoadFrom(TDirectory *dir, const std::string &name)
const double pot
Global constant settings.
Definition: UnfoldXSec.C:53
Spectrum with the value of a second variable, allowing for reweighting
static std::unique_ptr< ReweightableSpectrum > LoadFrom(TDirectory *dir, const std::string &name)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
TFile * outFile
Definition: PlotXSec.C:135
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Spectrum Background() const override
TH1 * UnfoldMethodRooUnfold(TH1 *measBinsHist, TH1 *trueBinsHist, TH2 *unfoldingMatrix, TH1 *dataMinusBkgdHist, TMatrixD *&covMatrix)
Definition: UnfoldXSec.C:139
std::vector< float > Spectrum
Definition: Constants.h:728
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
Just return the MC prediction for the background.
TH2D * ToTH2(double pot) const
gm Write()
bool Unfold3D ( string  systName,
string  varName,
TDirectory *  varDir,
TDirectory *  dataDir,
TFile *  outFile 
)

Set up all input directories

Load the input directories as cafe objects

Multiply the fake or real data by the purity

Grab the 3D binning info (3D vars only)

Convert to ROOT variables

Definition at line 155 of file UnfoldXSec.C.

References ana::TrivialBkgdEstimator::Background(), om::cout, covMatrix, allTimeWatchdog::endl, ana::kBinContent, ana::TrivialBkgdEstimator::LoadFrom(), ana::ReweightableSpectrum::LoadFrom(), ana::Spectrum::LoadFrom(), pot, ana::Ratio::ToTH1(), ana::Spectrum::ToTH1(), ana::ReweightableSpectrum::ToTH2(), ana::ToTH3Helper(), UnfoldMethodRooUnfold(), unfoldVars3D, and Write().

Referenced by UnfoldAllVars().

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 }
float covMatrix[xbins][xbins]
Definition: MakePlots.C:95
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
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
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
static std::unique_ptr< TrivialBkgdEstimator > LoadFrom(TDirectory *dir, const std::string &name)
const double pot
Global constant settings.
Definition: UnfoldXSec.C:53
Spectrum with the value of a second variable, allowing for reweighting
static std::unique_ptr< ReweightableSpectrum > LoadFrom(TDirectory *dir, const std::string &name)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
TFile * outFile
Definition: PlotXSec.C:135
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Spectrum Background() const override
TH1 * UnfoldMethodRooUnfold(TH1 *measBinsHist, TH1 *trueBinsHist, TH2 *unfoldingMatrix, TH1 *dataMinusBkgdHist, TMatrixD *&covMatrix)
Definition: UnfoldXSec.C:139
std::vector< float > Spectrum
Definition: Constants.h:728
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
Just return the MC prediction for the background.
TH2D * ToTH2(double pot) const
gm Write()
void UnfoldAllVars ( string  systName,
TDirectory *  systDir,
map< string, TDirectory * > &  dataDirs,
TFile *  outFile 
)

Definition at line 325 of file UnfoldXSec.C.

References Unfold1D(), Unfold3D(), unfoldVars1D, and unfoldVars3D.

Referenced by UnfoldBeamSysts(), UnfoldGenie(), UnfoldPPFX(), and UnfoldStandardSysts().

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 }
bool Unfold3D(string systName, string varName, TDirectory *varDir, TDirectory *dataDir, TFile *outFile)
Definition: UnfoldXSec.C:155
TFile * outFile
Definition: PlotXSec.C:135
const vector< string > unfoldVars1D
1D Unfolding variables
Definition: UnfoldXSec.C:77
map< string, vector< ana::Binning > > unfoldVars3D
Unfolding variables.
Definition: UnfoldXSec.C:72
bool Unfold1D(string systName, string varName, TDirectory *varDir, TDirectory *dataDir, TFile *outFile)
Definition: UnfoldXSec.C:245
void UnfoldBeamSysts ( TFile *  inputFile,
map< string, TDirectory * > &  dataDirs,
TFile *  outFile 
)

Definition at line 371 of file UnfoldXSec.C.

References beamSysts, om::cerr, om::cout, allTimeWatchdog::endl, and UnfoldAllVars().

Referenced by UnfoldXSec().

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 }
TFile * inputFile
Definition: PlotXSec.C:134
void UnfoldAllVars(string systName, TDirectory *systDir, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:325
OStream cerr
Definition: OStream.cxx:7
TFile * outFile
Definition: PlotXSec.C:135
OStream cout
Definition: OStream.cxx:6
const vector< string > beamSysts
Beam systematics (each has _Up and _Dw variations)
Definition: UnfoldXSec.C:116
void UnfoldGenie ( TFile *  inputFile,
map< string, TDirectory * > &  dataDirs,
TFile *  outFile,
const int  firstGenie,
const int  lastGenie 
)

Definition at line 397 of file UnfoldXSec.C.

References om::cout, allTimeWatchdog::endl, NUM_GENIE_UNIV, art::to_string(), and UnfoldAllVars().

Referenced by UnfoldXSec().

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 }
TFile * inputFile
Definition: PlotXSec.C:134
void UnfoldAllVars(string systName, TDirectory *systDir, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:325
TFile * outFile
Definition: PlotXSec.C:135
OStream cout
Definition: OStream.cxx:6
const int NUM_GENIE_UNIV
Definition: UnfoldXSec.C:135
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
TH1* UnfoldMethodRooUnfold ( TH1 *  measBinsHist,
TH1 *  trueBinsHist,
TH2 *  unfoldingMatrix,
TH1 *  dataMinusBkgdHist,
TMatrixD *&  covMatrix 
)

Definition at line 139 of file UnfoldXSec.C.

References confusionMatrixTree::count, om::cout, allTimeWatchdog::endl, nIterations, recentWatchdog::now, and art::to_string().

Referenced by Unfold1D(), and Unfold3D().

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 }
unsigned int nIterations
Definition: UnfoldXSec.C:132
::xsd::cxx::tree::duration< char, simple_type > duration
Definition: Database.h:188
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
OStream cout
Definition: OStream.cxx:6
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void UnfoldPPFX ( TFile *  inputFile,
map< string, TDirectory * > &  dataDirs,
TFile *  outFile,
const int  firstPPFX,
const int  lastPPFX 
)

Definition at line 419 of file UnfoldXSec.C.

References om::cout, allTimeWatchdog::endl, NUM_PPFX_UNIV, art::to_string(), and UnfoldAllVars().

Referenced by UnfoldXSec().

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 }
TFile * inputFile
Definition: PlotXSec.C:134
void UnfoldAllVars(string systName, TDirectory *systDir, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:325
const int NUM_PPFX_UNIV
Definition: UnfoldXSec.C:136
TFile * outFile
Definition: PlotXSec.C:135
OStream cout
Definition: OStream.cxx:6
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void UnfoldStandardSysts ( TFile *  inputFile,
map< string, TDirectory * > &  dataDirs,
TFile *  outFile 
)

Definition at line 355 of file UnfoldXSec.C.

References om::cout, allTimeWatchdog::endl, standardSysts, and UnfoldAllVars().

Referenced by UnfoldXSec().

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 }
TFile * inputFile
Definition: PlotXSec.C:134
void UnfoldAllVars(string systName, TDirectory *systDir, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:325
TFile * outFile
Definition: PlotXSec.C:135
OStream cout
Definition: OStream.cxx:6
const vector< string > standardSysts
Standard systmatics to check.
Definition: UnfoldXSec.C:96
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 at line 440 of file UnfoldXSec.C.

References ana::assert(), stan::math::ceil(), om::cerr, om::cout, allTimeWatchdog::endl, cet::getenv(), std::min(), nIterations, outFile, launch_batch_jobs::process, recoLabel, art::to_string(), UnfoldBeamSysts(), UnfoldGenie(), UnfoldPPFX(), UnfoldStandardSysts(), unfoldVars1D, unfoldVars3D, varToFakeData, varToRealData, and xs_label.

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 }
unsigned int nIterations
Definition: UnfoldXSec.C:132
map< string, string > varToRealData
Definition: UnfoldXSec.C:89
OStream cerr
Definition: OStream.cxx:7
void UnfoldBeamSysts(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:371
void UnfoldPPFX(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile, const int firstPPFX, const int lastPPFX)
Definition: UnfoldXSec.C:419
TFile * outFile
Definition: PlotXSec.C:135
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
void UnfoldStandardSysts(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile)
Definition: UnfoldXSec.C:355
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
OStream cout
Definition: OStream.cxx:6
map< string, vector< ana::Binning > > unfoldVars3D
Unfolding variables.
Definition: UnfoldXSec.C:72
const string xs_label
Definition: UnfoldXSec.C:62
assert(nhit_max >=nhit_nbins)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T min(const caf::Proxy< T > &a, T b)
void UnfoldGenie(TFile *inputFile, map< string, TDirectory * > &dataDirs, TFile *outFile, const int firstGenie, const int lastGenie)
Definition: UnfoldXSec.C:397
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11

Variable Documentation

const vector<string> beamSysts
Initial value:
= {
"2kA",
"02mmBeamSpotSize",
"1mmBeamShiftX",
"1mmBeamShiftY",
"3mmHorn1X",
"3mmHorn1Y",
"3mmHorn2X",
"3mmHorn2Y",
"7mmTargetZ",
"MagneticFieldinDecayPipe",
"1mmHornWater"
}

Beam systematics (each has _Up and _Dw variations)

Definition at line 116 of file UnfoldXSec.C.

Referenced by UnfoldBeamSysts().

const unsigned int DEFAULT_ITERATIONS = 3

Definition at line 131 of file UnfoldXSec.C.

const string fluxDirName = "fluxes/CV"

Directory of flux information in flux file.

Definition at line 65 of file UnfoldXSec.C.

unsigned int nIterations = DEFAULT_ITERATIONS

Definition at line 132 of file UnfoldXSec.C.

Referenced by UnfoldMethodRooUnfold(), and UnfoldXSec().

const double ntargs = 3.89137e+31

Scale to the FHC protons-on-target.

Definition at line 54 of file UnfoldXSec.C.

const int NUM_GENIE_UNIV = 1000

Definition at line 135 of file UnfoldXSec.C.

Referenced by UnfoldGenie().

const int NUM_PPFX_UNIV = 100

Definition at line 136 of file UnfoldXSec.C.

Referenced by UnfoldPPFX().

const double pot = 8.09E20

Global constant settings.

Definition at line 53 of file UnfoldXSec.C.

Referenced by Unfold1D(), and Unfold3D().

const string ppfxWgtLabel = "ppfxwgt"

Definition at line 60 of file UnfoldXSec.C.

const string ppfxXsecWgtLabel = "ppfxxsecwgt"

Definition at line 59 of file UnfoldXSec.C.

const string recoLabel = "reco"

Number of targets, calcuted by ???

Set whether to use ppfx-weighted or unweighted fake-data

Definition at line 57 of file UnfoldXSec.C.

Referenced by UnfoldXSec().

const vector<string> standardSysts
Initial value:
=
{
"cv",
"lightup",
"lightdown",
"calibpos",
"calibneg",
"calibshift",
"muoneup",
"muonedw",
"ckv",
"angleup",
"angledw",
"neutron_Up",
"neutron_Dw",
"neutron_nom",
"nom_good_seventh"
}

Standard systmatics to check.

Definition at line 96 of file UnfoldXSec.C.

Referenced by UnfoldStandardSysts().

const string trueLabel = "true"

Definition at line 58 of file UnfoldXSec.C.

const vector<string> unfoldVars1D
Initial value:
{
"ENu",
"Q2"
}

1D Unfolding variables

Definition at line 77 of file UnfoldXSec.C.

Referenced by UnfoldAllVars(), and UnfoldXSec().

map<string, vector<ana::Binning > > unfoldVars3D
Initial value:
{
}
const Binning eavailbins
const Binning mukebins
Definition: NumuCCIncBins.h:90
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72

Unfolding variables.

3D Unfolding variables. Need to include a mapping to the x,y,z binnings

Definition at line 72 of file UnfoldXSec.C.

Referenced by Unfold3D(), UnfoldAllVars(), and UnfoldXSec().

const string unwgtLabel = "unwgt"

Definition at line 61 of file UnfoldXSec.C.

map<string, string> varToFakeData
Initial value:
{
{"EAvail", "fakedata"},
{"ENu", "fakedata_enu"},
{"Q2", "fakedata_q2"}
}

Map the variable names to their folder in the fake data (could preferably switch the labelling of the fake data instead)

Definition at line 83 of file UnfoldXSec.C.

Referenced by UnfoldXSec().

map<string, string> varToRealData
Initial value:
{
{"EAvail", "data_eavail"},
{"ENu", "data_enu"},
{"Q2", "data_q2"}
}

Definition at line 89 of file UnfoldXSec.C.

Referenced by UnfoldXSec().

const string xs_label = ppfxXsecWgtLabel

Definition at line 62 of file UnfoldXSec.C.

Referenced by UnfoldXSec().