Public Types | Public Member Functions | Static Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
ana::CutOptimization Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-02-24/CAFAna/XSec/CutOptimization.h"

Public Types

enum  OptimizeMethod_t { kBinByBin, kIntegratedUp, kIntegratedDown, kIntegratedInsideBounds }
 
enum  FOM_t { kdSigmaOverSigma }
 

Public Member Functions

 CutOptimization (const HistAxis *hist_axis, const Cut *signal_cut, const Cut *selection_cut)
 
void DefineNominal (SystematicDef *nominal)
 
void DefineSystematic (SystematicDef *syst)
 
void DefineMultiverseSystematic (SystematicDef *syst, std::vector< Var > mv_weights)
 
void DefineMultiverseSystematic (SystematicDef *syst, std::vector< SystShifts > mv_shifts)
 
void SetSpillCuts (const SpillCut &)
 
void Go ()
 
void SaveTo (TDirectory *dir, const std::string &name) const
 
 CutOptimization (const CutOptimization &rhs)
 
 CutOptimization (CutOptimization &&rhs)
 
CutOptimizationoperator= (const CutOptimization &rhs)
 
CutOptimizationoperator= (CutOptimization &&rhs)
 
 ~CutOptimization ()=default
 
void Optimize (OptimizeMethod_t method, FOM_t fom, double data_pot, TDirectory *save_dir, std::string plot_dump=".", bool debug=false)
 
 CutOptimization (Spectrum *nominal_signal, Spectrum *nominal_bkgd, Spectrum *nominal_selected_signal, Spectrum *nominal_selected_bkgd, Spectrum *nominal_selected, std::vector< SystematicDef * > syst_defs, std::map< SystematicDef *, UpDownPair< Spectrum > > syst_signal, std::map< SystematicDef *, UpDownPair< Spectrum > > syst_bkgd, std::map< SystematicDef *, UpDownPair< Spectrum > > syst_selected_signal, std::map< SystematicDef *, UpDownPair< Spectrum > > syst_selected_bkgd, std::map< SystematicDef *, UpDownPair< Spectrum > > syst_selected, std::vector< SystematicDef * > mv_defs, std::map< SystematicDef *, Multiverse * > mv_signal, std::map< SystematicDef *, Multiverse * > mv_bkgd, std::map< SystematicDef *, Multiverse * > mv_selected_signal, std::map< SystematicDef *, Multiverse * > mv_selected_bkgd, std::map< SystematicDef *, Multiverse * > mv_selected)
 
std::vector< SystematicDef * > GetSystDefs ()
 
std::vector< SystematicDef * > GetMVSystDefs ()
 
std::vector< SystematicDef * > GetAllSystDefs ()
 

Static Public Member Functions

static void PlotDebug (TH1 *hist, std::string draw_option, std::string title, std::string name)
 
static void PlotDebug (std::vector< TH1 * > hists, std::vector< std::string > labels, std::string draw_option, std::string title, std::string name, double ymax=-1)
 
static void PlotDebug (const std::vector< TH1 * > universes, Spectrum *nominal, std::string title, std::string name)
 
static void PlotDebug (const std::vector< TH1 * > universes, TH1 *h_nominal, std::string title, std::string name)
 
static std::unique_ptr< CutOptimizationLoadFrom (TDirectory *dir, const std::string &name)
 
static void Integrate (TH1 *&result, const TH1 *tmp, OptimizeMethod_t method)
 
static void Integrate (Multiverse *mv, OptimizeMethod_t method)
 

Private Member Functions

TH1 * AbsUncertainty (const UpDownPair< TH1 > syst, const TH1 *nominal) const
 
TH1 * AbsUncertaintySquared (const UpDownPair< TH1 > syst, const TH1 *nominal) const
 
UpDownPair< TH1 > ToUpDownHist (Multiverse *mv, const TH1 *h_nominal)
 
void OptimizedSigmaOverSigma (OptimizeMethod_t method, double data_pot, TDirectory *save_dir, std::string plot_dump=".", bool debug=false)
 
UpDownPair< TH1 > Efficiency (UpDownPair< TH1 > num, UpDownPair< TH1 > denom, OptimizeMethod_t method)
 
MultiverseEfficiency (Multiverse *num, Multiverse *denom, OptimizeMethod_t method)
 

Static Private Member Functions

static TH1 * ToHist (const Spectrum *spec, OptimizeMethod_t method, double POT=-1)
 
static TH1 * Efficiency (TH1 *num, TH1 *denom, OptimizeMethod_t method)
 

Private Attributes

SystematicDeffNominalSystDef
 
SpectrumfNominalSignal
 
SpectrumfNominalBkgd
 
SpectrumfNominalSelectedSignal
 
SpectrumfNominalSelectedBkgd
 
SpectrumfNominalSelected
 
std::vector< SystematicDef * > fSystDefs
 
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
 
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
 
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
 
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
 
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
 
std::vector< SystematicDef * > fMVSystDefs
 
std::map< SystematicDef *, Multiverse * > fMVSignal
 
std::map< SystematicDef *, Multiverse * > fMVBkgd
 
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
 
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
 
std::map< SystematicDef *, Multiverse * > fMVSelected
 
const CutfSignalCut = NULL
 
const CutfSelectionCut = NULL
 
const HistAxisfOptiHistAxis = NULL
 

Detailed Description

Definition at line 21 of file CutOptimization.h.

Member Enumeration Documentation

Enumerator
kdSigmaOverSigma 

Definition at line 55 of file CutOptimization.h.

Enumerator
kBinByBin 
kIntegratedUp 
kIntegratedDown 
kIntegratedInsideBounds 

Definition at line 49 of file CutOptimization.h.

49  {
50  kBinByBin, // bin by bin
51  kIntegratedUp, // integrate bin and above
52  kIntegratedDown, // integrate bin and below
53  kIntegratedInsideBounds, // integrate between upper and lower bounds
54  };

Constructor & Destructor Documentation

ana::CutOptimization::CutOptimization ( const HistAxis hist_axis,
const Cut signal_cut,
const Cut selection_cut 
)

Definition at line 134 of file CutOptimization.cxx.

137  : fSignalCut(signal_cut), fSelectionCut(selection_cut), fOptiHistAxis(opti)
138  {}
const HistAxis * fOptiHistAxis
const Cut signal_cut
ana::CutOptimization::CutOptimization ( const CutOptimization rhs)

Definition at line 460 of file CutOptimization.cxx.

References ana::CopyUpDownSpectrum(), fBkgd, fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fNominalBkgd, fNominalSelected, fNominalSelectedBkgd, fNominalSelectedSignal, fNominalSignal, fNominalSystDef, fOptiHistAxis, fSelected, fSelectedBkgd, fSelectedSignal, fSelectionCut, fSignal, fSignalCut, fSystDefs, and submit_hadd::u.

461  {
462  this->fSignalCut = new Cut(*rhs.fSignalCut);
463  this->fSelectionCut = new Cut(*rhs.fSignalCut);
464  this->fOptiHistAxis = new HistAxis(*rhs.fOptiHistAxis);
465 
466  this->fNominalSystDef = new SystematicDef(*rhs.fNominalSystDef);
467  this->fNominalSignal = new Spectrum(*rhs.fNominalSignal );
468  this->fNominalBkgd = new Spectrum(*rhs.fNominalBkgd );
469  this->fNominalSelectedSignal = new Spectrum(*rhs.fNominalSelectedSignal);
470  this->fNominalSelectedBkgd = new Spectrum(*rhs.fNominalSelectedBkgd );
471  this->fNominalSelected = new Spectrum(*rhs.fNominalSelected );
472  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
473  this->fSystDefs.push_back(new SystematicDef(*rhs.fSystDefs[isyst]));
474  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSignal .at(rhs.fSystDefs[isyst]));
475  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fBkgd .at(rhs.fSystDefs[isyst]));
476  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst]));
477  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst]));
478  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelected .at(rhs.fSystDefs[isyst]));
479  }
480  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
481  this->fMVSystDefs.push_back(new SystematicDef(*rhs.fMVSystDefs[imv]));
482  this->fMVSignal [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
483  this->fMVBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
484  this->fMVSelectedSignal[this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
485  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
486  this->fMVSelected [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
487  }
488  }
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
const HistAxis * fOptiHistAxis
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
GenericSystematicDef< caf::SRProxy > SystematicDef
UpDownPair< Spectrum > CopyUpDownSpectrum(const UpDownPair< Spectrum > &copy)
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
Spectrum * fNominalSelectedSignal
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
std::vector< float > Spectrum
Definition: Constants.h:728
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
SystematicDef * fNominalSystDef
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
ana::CutOptimization::CutOptimization ( CutOptimization &&  rhs)

Definition at line 492 of file CutOptimization.cxx.

References ana::CopyUpDownSpectrum(), fBkgd, fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fNominalBkgd, fNominalSelected, fNominalSelectedBkgd, fNominalSelectedSignal, fNominalSignal, fNominalSystDef, fOptiHistAxis, fSelected, fSelectedBkgd, fSelectedSignal, fSelectionCut, fSignal, fSignalCut, fSystDefs, and submit_hadd::u.

493  {
494  this->fSignalCut = std::move(rhs.fSignalCut);
495  this->fSelectionCut = std::move(rhs.fSignalCut);
496  this->fOptiHistAxis = std::move(rhs.fOptiHistAxis);
497 
498  this->fNominalSystDef = std::move(rhs.fNominalSystDef );
499  this->fNominalSignal = std::move(rhs.fNominalSignal );
500  this->fNominalBkgd = std::move(rhs.fNominalBkgd );
501  this->fNominalSelectedSignal = std::move(rhs.fNominalSelectedSignal);
502  this->fNominalSelectedBkgd = std::move(rhs.fNominalSelectedBkgd );
503  this->fNominalSelected = std::move(rhs.fNominalSelected );
504 
505  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
506  this->fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
507  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSignal .at(rhs.fSystDefs[isyst])));
508  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fBkgd .at(rhs.fSystDefs[isyst])));
509  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst])));
510  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst])));
511  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelected .at(rhs.fSystDefs[isyst])));
512  }
513  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
514  this->fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
515  this->fMVSignal [this->fMVSystDefs.back()] = std::move(rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
516  this->fMVBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
517  this->fMVSelectedSignal[this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
518  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
519  this->fMVSelected [this->fMVSystDefs.back()] = std::move(rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
520  }
521 
522  }
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
const HistAxis * fOptiHistAxis
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
UpDownPair< Spectrum > CopyUpDownSpectrum(const UpDownPair< Spectrum > &copy)
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
Spectrum * fNominalSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
SystematicDef * fNominalSystDef
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
ana::CutOptimization::~CutOptimization ( )
default
ana::CutOptimization::CutOptimization ( Spectrum nominal_signal,
Spectrum nominal_bkgd,
Spectrum nominal_selected_signal,
Spectrum nominal_selected_bkgd,
Spectrum nominal_selected,
std::vector< SystematicDef * >  syst_defs,
std::map< SystematicDef *, UpDownPair< Spectrum > >  syst_signal,
std::map< SystematicDef *, UpDownPair< Spectrum > >  syst_bkgd,
std::map< SystematicDef *, UpDownPair< Spectrum > >  syst_selected_signal,
std::map< SystematicDef *, UpDownPair< Spectrum > >  syst_selected_bkgd,
std::map< SystematicDef *, UpDownPair< Spectrum > >  syst_selected,
std::vector< SystematicDef * >  mv_defs,
std::map< SystematicDef *, Multiverse * >  mv_signal,
std::map< SystematicDef *, Multiverse * >  mv_bkgd,
std::map< SystematicDef *, Multiverse * >  mv_selected_signal,
std::map< SystematicDef *, Multiverse * >  mv_selected_bkgd,
std::map< SystematicDef *, Multiverse * >  mv_selected 
)

Definition at line 423 of file CutOptimization.cxx.

440  : fNominalSignal(nominal_signal),
441  fNominalBkgd(nominal_bkgd),
442  fNominalSelectedSignal(nominal_selected_signal),
443  fNominalSelectedBkgd(nominal_selected_bkgd),
444  fNominalSelected(nominal_selected),
445  fSystDefs(syst_defs),
446  fSignal(syst_signal),
447  fBkgd(syst_bkgd),
448  fSelectedSignal(syst_selected_signal),
449  fSelectedBkgd(syst_selected_bkgd),
450  fSelected(syst_selected),
451  fMVSystDefs(mv_defs),
452  fMVSignal(mv_signal),
453  fMVBkgd(mv_bkgd),
454  fMVSelectedSignal(mv_selected_signal),
455  fMVSelectedBkgd(mv_selected_bkgd),
456  fMVSelected(mv_selected)
457  {}
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
Spectrum * fNominalSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal

Member Function Documentation

TH1 * ana::CutOptimization::AbsUncertainty ( const UpDownPair< TH1 >  syst,
const TH1 *  nominal 
) const
private

Definition at line 1325 of file CutOptimization.cxx.

References std::abs(), ana::UpDownPair< W >::down, std::max(), runNovaSAM::ret, and ana::UpDownPair< W >::up.

Referenced by AbsUncertaintySquared(), and GetAllSystDefs().

1326  {
1327  TH1 * ret = (TH1*) syst.up->Clone();
1328  int NX = syst.up->GetNbinsX();
1329  int NY = syst.up->GetNbinsY();
1330  int NZ = syst.up->GetNbinsZ();
1331  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
1332  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
1333  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
1334  double nom = nominal->GetBinContent(ibinx, ibiny, ibinz);
1335 
1336  double up = syst.up->GetBinContent(ibinx, ibiny, ibinz);
1337  double down = nom;
1338  if(syst.down) down = syst.down->GetBinContent(ibinx, ibiny, ibinz);
1339 
1340  double error = std::max(std::abs(up - nom), std::abs(down - nom));
1341  ret->SetBinContent(ibinx, ibiny, ibinz, error);
1342  }
1343  }
1344  }
1345 
1346  return ret;
1347  }
T max(const caf::Proxy< T > &a, T b)
float abs(float number)
Definition: d0nt_math.hpp:39
TH1 * ana::CutOptimization::AbsUncertaintySquared ( const UpDownPair< TH1 >  syst,
const TH1 *  nominal 
) const
private

Definition at line 1351 of file CutOptimization.cxx.

References AbsUncertainty(), and runNovaSAM::ret.

Referenced by GetAllSystDefs().

1352  {
1353  TH1 * ret = AbsUncertainty(syst, nominal);
1354  ret->Multiply(ret);
1355  return ret;
1356  }
TH1 * AbsUncertainty(const UpDownPair< TH1 > syst, const TH1 *nominal) const
void ana::CutOptimization::DefineMultiverseSystematic ( SystematicDef syst,
std::vector< Var mv_weights 
)

Definition at line 194 of file CutOptimization.cxx.

References fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fOptiHistAxis, fSelectionCut, fSignalCut, ana::GenericSystematicDef< SRType >::Loaders(), ana::GenericSystematicDef< SRType >::Shifts(), and ana::GenericSystematicDef< SRType >::Weight().

Referenced by containment_optimization(), demo_cut_optimization(), fiducial_opt(), fiducial_optimization(), muonid_opt(), muonid_optimization(), and prongcvn_optimization().

196  {
197  fMVSystDefs.push_back(syst);
198  fMVSignal[syst] = new Multiverse(syst->Loaders(),
199  *fOptiHistAxis,
200  *fSignalCut,
201  *syst->Shifts(),
202  mv_weights,
203  *syst->Weight());
204 
205  fMVBkgd[syst] = new Multiverse(syst->Loaders(),
206  *fOptiHistAxis,
207  !*fSignalCut,
208  *syst->Shifts(),
209  mv_weights,
210  *syst->Weight());
211  fMVSelectedSignal[syst] = new Multiverse(syst->Loaders(),
212  *fOptiHistAxis,
214  *syst->Shifts(),
215  mv_weights,
216  *syst->Weight());
217  fMVSelectedBkgd[syst] = new Multiverse(syst->Loaders(),
218  *fOptiHistAxis,
220  *syst->Shifts(),
221  mv_weights,
222  *syst->Weight());
223  fMVSelected[syst] = new Multiverse(syst->Loaders(),
224  *fOptiHistAxis,
225  *fSelectionCut,
226  *syst->Shifts(),
227  mv_weights,
228  *syst->Weight());
229  }
const HistAxis * fOptiHistAxis
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
std::map< SystematicDef *, Multiverse * > fMVSignal
void ana::CutOptimization::DefineMultiverseSystematic ( SystematicDef syst,
std::vector< SystShifts mv_shifts 
)

Definition at line 233 of file CutOptimization.cxx.

References fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fOptiHistAxis, fSelectionCut, fSignalCut, ana::GenericSystematicDef< SRType >::Loaders(), and ana::GenericSystematicDef< SRType >::Weight().

234  {
235  fMVSystDefs.push_back(syst);
236  fMVSignal[syst] = new Multiverse(syst->Loaders(),
237  *fOptiHistAxis,
238  *fSignalCut,
239  mv_shifts,
240  *syst->Weight());
241  fMVBkgd[syst] = new Multiverse(syst->Loaders(),
242  *fOptiHistAxis,
243  !*fSignalCut,
244  mv_shifts,
245  *syst->Weight());
246  fMVSelectedSignal[syst] = new Multiverse(syst->Loaders(),
247  *fOptiHistAxis,
249  mv_shifts,
250  *syst->Weight());
251  fMVSelectedBkgd[syst] = new Multiverse(syst->Loaders(),
252  *fOptiHistAxis,
254  mv_shifts,
255  *syst->Weight());
256  fMVSelected[syst] = new Multiverse(syst->Loaders(),
257  *fOptiHistAxis,
258  *fSelectionCut,
259  mv_shifts,
260  *syst->Weight());
261  }
const HistAxis * fOptiHistAxis
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
std::map< SystematicDef *, Multiverse * > fMVSignal
void ana::CutOptimization::DefineNominal ( SystematicDef nominal)

Definition at line 141 of file CutOptimization.cxx.

References ana::GenericSystematicDef< SRType >::BuildSpectrum(), fNominalBkgd, fNominalSelected, fNominalSelectedBkgd, fNominalSelectedSignal, fNominalSignal, fNominalSystDef, fSelectionCut, and fSignalCut.

Referenced by containment_optimization(), demo_cut_optimization(), fiducial_opt(), fiducial_optimization(), muonid_opt(), muonid_optimization(), and prongcvn_optimization().

void ana::CutOptimization::DefineSystematic ( SystematicDef syst)

Definition at line 154 of file CutOptimization.cxx.

References ana::GenericSystematicDef< SRType >::BuildSpectrumDown(), ana::GenericSystematicDef< SRType >::BuildSpectrumUp(), fBkgd, fSelected, fSelectedBkgd, fSelectedSignal, fSelectionCut, fSignal, fSignalCut, fSystDefs, and ana::GenericSystematicDef< SRType >::IsTwoSided().

Referenced by containment_optimization(), demo_cut_optimization(), fiducial_opt(), fiducial_optimization(), muonid_opt(), muonid_optimization(), and prongcvn_optimization().

155  {
156  fSystDefs.push_back(syst);
157 
158  // for one sided systs, store in .up
159  fSignal [syst].up = syst->BuildSpectrumUp(*fSignalCut);
160  fBkgd [syst].up = syst->BuildSpectrumUp(!*fSignalCut);
161  fSelectedSignal[syst].up = syst->BuildSpectrumUp(*fSelectionCut && *fSignalCut);
162  fSelectedBkgd [syst].up = syst->BuildSpectrumUp(*fSelectionCut && !*fSignalCut);
163  fSelected [syst].up = syst->BuildSpectrumUp(*fSelectionCut);
164 
165  if(syst->IsTwoSided()) {
166  fSignal [syst].down = syst->BuildSpectrumDown(*fSignalCut);
167  fBkgd [syst].down = syst->BuildSpectrumDown(!*fSignalCut);
168  fSelectedSignal[syst].down = syst->BuildSpectrumDown(*fSelectionCut && *fSignalCut);
169  fSelectedBkgd [syst].down = syst->BuildSpectrumDown(*fSelectionCut && !*fSignalCut);
170  fSelected [syst].down = syst->BuildSpectrumDown(*fSelectionCut);
171  }
172  }
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
TH1 * ana::CutOptimization::Efficiency ( TH1 *  num,
TH1 *  denom,
OptimizeMethod_t  method 
)
staticprivate

Definition at line 789 of file CutOptimization.cxx.

References nd_projection_maker::eff, efficiency(), febshutoff_auto::integral, kIntegratedDown, kIntegratedInsideBounds, kIntegratedUp, std::sqrt(), and ana::UniqueName().

Referenced by Efficiency(), GetAllSystDefs(), and OptimizedSigmaOverSigma().

790  {
791  // numerator should already be integrated appropriately
792  TH1 * eff = (TH1*) num->Clone(UniqueName().c_str());
793 
794  // if using an integration method, denom should also be integrated appropriately
795  // and if so, scale numerator by total number of signal events, which is found
796  // in either the first or last bin of denominator
797  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
798  auto NX = denom->GetNbinsX();
799  auto NY = denom->GetNbinsY();
800  auto NZ = denom->GetNbinsZ();
801  double integral = 0;
802  if(method == kIntegratedUp || method == kIntegratedDown) {
803  integral = method == kIntegratedUp ? denom->GetBinContent(1) : denom->GetBinContent(NX, NY, NZ);
804  }
805  else if(method == kIntegratedInsideBounds) {
806  // in this case, the total number of signal events will be
807  // where the upper bound is at a max and lower bound is at a minimum
808  // y axis is upper bound,
809  // x axis is lower bound
810  integral = denom->GetBinContent(1, NY, 1);
811  }
812 
813  eff->Scale(1/integral);
814  // compute binomial errors by hand
815  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
816  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
817  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
818  double efficiency = eff->GetBinContent(ibinx, ibiny, ibinz);
819  double binom_error = 1/integral * std::sqrt( efficiency * (1 - efficiency / integral) );
820  eff->SetBinError(ibinx, ibiny, ibinz, binom_error);
821  }
822  }
823  }
824  }
825  else {
826  // else just divide numerator with denominator
827  eff->Divide(eff, denom, 1, 1, "B"); // divide with binomial error propagation
828  }
829 
830  eff->SetMaximum(1.2*eff->GetMaximum());
831 
832  return eff;
833  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
int num
Definition: f2_nu.C:119
void efficiency()
Definition: efficiency.C:58
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
UpDownPair< TH1 > ana::CutOptimization::Efficiency ( UpDownPair< TH1 >  num,
UpDownPair< TH1 >  denom,
OptimizeMethod_t  method 
)
private

Definition at line 837 of file CutOptimization.cxx.

References ana::UpDownPair< W >::down, Efficiency(), and ana::UpDownPair< W >::up.

838  {
839  TH1 * up = Efficiency(num.up, denom.up, method);
840  TH1 * down = num.down == NULL? NULL : Efficiency(num.down, denom.down, method);
841  return UpDownPair<TH1>{up, down};
842  }
static TH1 * Efficiency(TH1 *num, TH1 *denom, OptimizeMethod_t method)
int num
Definition: f2_nu.C:119
Multiverse * ana::CutOptimization::Efficiency ( Multiverse num,
Multiverse denom,
OptimizeMethod_t  method 
)
private

Definition at line 845 of file CutOptimization.cxx.

References ana::Multiverse::Divide(), nd_projection_maker::eff, Efficiency(), kIntegratedDown, kIntegratedInsideBounds, kIntegratedUp, and ana::Multiverse::Transform().

848  {
849  // calculate uncertainty on efficiency of a multiverse as the
850  // +- 1 sigma band of the efficiency multiverse
851  Multiverse * eff = new Multiverse(*num);
852  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
853  eff->Transform(*denom,
854  [method](TH1 * univ_num, TH1 * univ_denom)
855  {
856  return CutOptimization::Efficiency(univ_num, univ_denom, method);
857  });
858  }
859  else {
860  eff->Divide(*denom, true);
861  }
862  return eff;
863  }
static TH1 * Efficiency(TH1 *num, TH1 *denom, OptimizeMethod_t method)
int num
Definition: f2_nu.C:119
std::vector<SystematicDef*> ana::CutOptimization::GetAllSystDefs ( )
inline

Definition at line 108 of file CutOptimization.h.

References AbsUncertainty(), AbsUncertaintySquared(), gen_hdf5record::debug, Efficiency(), fMVSystDefs, fSystDefs, num, OptimizedSigmaOverSigma(), POT, runNovaSAM::ret, string, ToHist(), and ToUpDownHist().

108  {
109  std::vector<SystematicDef*> ret = fSystDefs;
110  ret.insert(fSystDefs.end(), fMVSystDefs.begin(), fMVSystDefs.end());
111  return ret;
112  }
std::vector< SystematicDef * > fSystDefs
std::vector< SystematicDef * > fMVSystDefs
std::vector<SystematicDef*> ana::CutOptimization::GetMVSystDefs ( )
inline

Definition at line 107 of file CutOptimization.h.

References fMVSystDefs.

107 { return fMVSystDefs; }
std::vector< SystematicDef * > fMVSystDefs
std::vector<SystematicDef*> ana::CutOptimization::GetSystDefs ( )
inline

Definition at line 106 of file CutOptimization.h.

References fSystDefs.

106 { return fSystDefs; }
std::vector< SystematicDef * > fSystDefs
void ana::CutOptimization::Go ( )

Definition at line 265 of file CutOptimization.cxx.

References fMVSystDefs, fNominalSystDef, fSystDefs, ana::GenericSystematicDef< SRType >::Go(), and submit_hadd::u.

Referenced by containment_optimization(), demo_cut_optimization(), fiducial_opt(), fiducial_optimization(), muonid_opt(), muonid_optimization(), and prongcvn_optimization().

266  {
267  fNominalSystDef->Go();
268  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++)
269  fSystDefs[isyst]->Go();
270  for(auto isyst = 0u; isyst < fMVSystDefs.size(); isyst++)
271  fMVSystDefs[isyst]->Go();
272  }
std::vector< SystematicDef * > fSystDefs
std::vector< SystematicDef * > fMVSystDefs
SystematicDef * fNominalSystDef
void ana::CutOptimization::Integrate ( TH1 *&  result,
const TH1 *  tmp,
OptimizeMethod_t  method 
)
static

Definition at line 595 of file CutOptimization.cxx.

References ana::assert(), genie::utils::style::Format(), febshutoff_auto::integral, kIntegratedDown, kIntegratedInsideBounds, fillBadChanDBTables::result, tmp, and ana::UniqueName().

Referenced by Integrate(), OptimizedSigmaOverSigma(), and ToHist().

596  {
597  auto ndims = tmp->GetDimension();
598  auto NX = tmp->GetNbinsX();
599  auto NY = tmp->GetNbinsY();
600  auto NZ = tmp->GetNbinsZ();
601 
602  // Useful for a 2-sided cut on a single variable like vertex position and fiducial
603  // constraints.
604  // Fills a 2D histogram where
605  // -- on the x axis is the value of the lower bound cut
606  // -- on the y axis is the value of the upper bound cut
607  // bin contents are the integral of the 1D input histogram between the bounds corresponding to
608  // upper bin edges
609  if(method == kIntegratedInsideBounds) {
610  assert(ndims == 1 && "I don't know how to handle this histogram for kIntegratedInsideBounds");
611 
612  // don't leak clone of 1D
613  delete result;
614 
615  // allow for variable bin width
616  auto array = tmp->GetXaxis()->GetXbins()->GetArray();
617  if(array) {
618  result = new TH2D(UniqueName().c_str(), "",
619  NX, tmp->GetXaxis()->GetXbins()->GetArray(),
620  NX, tmp->GetXaxis()->GetXbins()->GetArray());
621  }
622  else {
623  auto minx = tmp->GetXaxis()->GetBinLowEdge(1);
624  auto maxx = tmp->GetXaxis()->GetBinLowEdge(NX+1);
625  result = new TH2D(UniqueName().c_str(), "",
626  NX, minx, maxx,
627  NX, minx, maxx);
628 
629  }
630  result->GetYaxis()->SetTitle(TString::Format("%s %s", tmp->GetXaxis()->GetTitle(), "(Upper)"));
631  result->GetXaxis()->SetTitle(TString::Format("%s %s", tmp->GetXaxis()->GetTitle(), "(Lower)"));
632 
633  for(auto ibinlower = 1; ibinlower <= NX; ibinlower++) {
634  for(auto ibinupper = ibinlower; ibinupper <=NX; ibinupper++) {
635  double error = 0;
636  double integral = tmp->IntegralAndError(ibinlower, ibinupper, error);
637  result->SetBinContent(ibinlower, ibinupper, integral);
638  result->SetBinContent(ibinlower, ibinupper, integral);
639  }
640  }
641  return;
642  }
643  else {
644  switch(ndims) {
645  case 1:
646  {
647  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
648  auto x_a = method == kIntegratedDown? 0 : ibinx;
649  auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
650  double error = 0;
651  double integral = tmp->IntegralAndError(x_a, x_b, error);
652  result->SetBinContent(ibinx, integral);
653  result->SetBinError (ibinx, error );
654 
655  }
656  break;
657  }
658  case 2:
659  {
660  // this way, the casts are only performed twice
661  TH2 * result_2d = static_cast<TH2*> (result);
662  const TH2 * tmp_2d = static_cast<const TH2*>(tmp );
663  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
664  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
665  auto x_a = method == kIntegratedDown? 0 : ibinx;
666  auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
667  auto y_a = method == kIntegratedDown? 0 : ibiny;
668  auto y_b = method == kIntegratedDown? ibiny : NY+1 ;
669 
670  double error = 0;
671  double integral = tmp_2d->IntegralAndError(x_a, x_b,
672  y_a, y_b, error);
673  result_2d->SetBinContent(ibinx, ibiny, integral);
674  result_2d->SetBinError (ibinx, ibiny, error );
675  }
676  }
677  result = static_cast<TH1*>(result_2d);
678  break;
679  }
680  case 3:
681  {
682  // this way, the casts are only performed twice
683  TH3 * result_3d = static_cast<TH3*> (result);
684  const TH3 * tmp_3d = static_cast<const TH3*>(tmp );
685  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
686  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
687  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
688  auto x_a = method == kIntegratedDown? 0 : ibinx;
689  auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
690  auto y_a = method == kIntegratedDown? 0 : ibiny;
691  auto y_b = method == kIntegratedDown? ibiny : NY+1 ;
692  auto z_a = method == kIntegratedDown? 0 : ibinz;
693  auto z_b = method == kIntegratedDown? ibinz : NZ+1 ;
694 
695  double error = 0;
696  double integral = tmp_3d->IntegralAndError(x_a, x_b,
697  y_a, y_b,
698  z_a, z_b, error);
699  result_3d->SetBinContent(ibinx, ibiny, ibinz, integral);
700  result_3d->SetBinError (ibinx, ibiny, ibinz, error );
701  }
702  }
703  }
704  result = static_cast<TH1*>(result_3d);
705  break;
706  }
707  }
708  }
709 
710  }
Float_t tmp
Definition: plot.C:36
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void ana::CutOptimization::Integrate ( Multiverse mv,
OptimizeMethod_t  method 
)
static

Definition at line 867 of file CutOptimization.cxx.

References Integrate(), fillBadChanDBTables::result, and ana::Multiverse::Transform().

868  {
869  mv->Transform([method](TH1 * univ)
870  {
871  TH1 * result = (TH1*) univ->Clone();
872  CutOptimization::Integrate(result, univ, method);
873  return result;
874  });
875  }
static void Integrate(TH1 *&result, const TH1 *tmp, OptimizeMethod_t method)
std::unique_ptr< CutOptimization > ana::CutOptimization::LoadFrom ( TDirectory *  dir,
const std::string name 
)
static

Definition at line 343 of file CutOptimization.cxx.

References ana::assert(), dir, genie::utils::style::Format(), ana::Multiverse::LoadFrom(), ana::GenericSystematicDef< SRType >::LoadFrom(), ana::Spectrum::LoadFrom(), ana::LoadFromUpDownSpectra(), make_syst_table_plots::nsysts, and runNovaSAM::release.

Referenced by containment_optimization_plots(), demo_cut_optimization(), fiducial_optimization_plots(), fiducial_optimization_plots_2d(), and muonid_optimization_plots().

344  {
345  dir = dir->GetDirectory(name.c_str()); // switch to subdir
346  assert(dir);
347 
348  Spectrum * nominal_signal;
349  Spectrum * nominal_bkgd;
350  Spectrum * nominal_selected_signal;
351  Spectrum * nominal_selected_bkgd;
352  Spectrum * nominal_selected;
353 
354  std::vector<SystematicDef *> syst_defs;
355  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_signal;
356  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_bkgd;
357  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_signal;
358  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_bkgd;
359  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected;
360 
361  std::vector<SystematicDef *> mv_defs;
362  std::map<SystematicDef *, Multiverse * > mv_signal;
363  std::map<SystematicDef *, Multiverse * > mv_bkgd;
364  std::map<SystematicDef *, Multiverse * > mv_selected_signal;
365  std::map<SystematicDef *, Multiverse * > mv_selected_bkgd;
366  std::map<SystematicDef *, Multiverse * > mv_selected;
367 
368  TDirectory * nom_dir = dir->GetDirectory("Nominal");
369  nominal_signal = Spectrum::LoadFrom(nom_dir, "Signal" ).release();
370  nominal_bkgd = Spectrum::LoadFrom(nom_dir, "Bkgd" ).release();
371  nominal_selected_signal = Spectrum::LoadFrom(nom_dir, "SelectedSignal").release();
372  nominal_selected_bkgd = Spectrum::LoadFrom(nom_dir, "SelectedBkgd" ).release();
373  nominal_selected = Spectrum::LoadFrom(nom_dir, "Selected" ).release();
374 
375  TVectorD * vsysts = (TVectorD*) dir->Get("nsysts");
376  int nsysts = (*vsysts)[0];
377  for(int isyst = 0; isyst < nsysts; isyst++) {
378  syst_defs.push_back(SystematicDef::LoadFrom(dir, TString::Format("syst_defs/%d", isyst).Data()).release());
379 
380  TDirectory * syst_dir = dir->GetDirectory(TString::Format("syst_%s", syst_defs.back()->GetName().c_str()));
381  syst_signal [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSignal" ));
382  syst_bkgd [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fBkgd" ));
383  syst_selected_signal[syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelectedSignal" ));
384  syst_selected_bkgd [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelectedBkgd" ));
385  syst_selected [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelected" ));
386  }
387 
388  TVectorD * vmultiverse = (TVectorD*) dir->Get("nmultiverse");
389  int nmultiverse = (*vmultiverse)[0];
390  for(int imv = 0; imv < nmultiverse; imv++) {
391  mv_defs.push_back(SystematicDef::LoadFrom(dir, TString::Format("mv_defs/%d", imv).Data()).release());
392 
393  TDirectory * mv_dir = dir->GetDirectory(TString::Format("mv_%s", mv_defs.back()->GetName().c_str()));
394  mv_signal [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSignal" ).release();
395  mv_bkgd [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVBkgd" ).release();
396  mv_selected_signal[mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelectedSignal" ).release();
397  mv_selected_bkgd [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelectedBkgd" ).release();
398  mv_selected [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelected" ).release();
399  }
400 
401  delete dir;
402 
403  return std::make_unique<CutOptimization> (nominal_signal,
404  nominal_bkgd,
405  nominal_selected_signal,
406  nominal_selected_bkgd,
407  nominal_selected,
408  syst_defs,
409  syst_signal,
410  syst_bkgd,
411  syst_selected_signal,
412  syst_selected_bkgd,
413  syst_selected,
414  mv_defs,
415  mv_signal,
416  mv_bkgd,
417  mv_selected_signal,
418  mv_selected_bkgd,
419  mv_selected);
420  }
const XML_Char * name
Definition: expat.h:151
static std::unique_ptr< GenericSystematicDef< SRType > > LoadFrom(TDirectory *dir, const std::string &name)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
std::vector< float > Spectrum
Definition: Constants.h:728
static std::unique_ptr< Multiverse > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Multiverse.cxx:573
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
UpDownPair< Spectrum > LoadFromUpDownSpectra(TDirectory *dir)
CutOptimization & ana::CutOptimization::operator= ( const CutOptimization rhs)

Definition at line 526 of file CutOptimization.cxx.

References ana::CopyUpDownSpectrum(), fBkgd, fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fNominalBkgd, fNominalSelected, fNominalSelectedBkgd, fNominalSelectedSignal, fNominalSignal, fNominalSystDef, fOptiHistAxis, fSelected, fSelectedBkgd, fSelectedSignal, fSelectionCut, fSignal, fSignalCut, fSystDefs, and submit_hadd::u.

527  {
528  if(this == &rhs) return *this;
529 
530  this->fSignalCut = new Cut(*rhs.fSignalCut);
531  this->fSelectionCut = new Cut(*rhs.fSignalCut);
532  this->fOptiHistAxis = new HistAxis(*rhs.fOptiHistAxis);
533 
534  this->fNominalSystDef = new SystematicDef(*rhs.fNominalSystDef);
535  this->fNominalSignal = new Spectrum(*rhs.fNominalSignal );
536  this->fNominalBkgd = new Spectrum(*rhs.fNominalBkgd );
537  this->fNominalSelectedSignal = new Spectrum(*rhs.fNominalSelectedSignal);
538  this->fNominalSelectedBkgd = new Spectrum(*rhs.fNominalSelectedBkgd );
539  this->fNominalSelected = new Spectrum(*rhs.fNominalSelected );
540  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
541  this->fSystDefs.push_back(new SystematicDef(*rhs.fSystDefs[isyst]));
542  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSignal .at(rhs.fSystDefs[isyst]));
543  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fBkgd .at(rhs.fSystDefs[isyst]));
544  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst]));
545  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst]));
546  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelected .at(rhs.fSystDefs[isyst]));
547  }
548  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
549  this->fMVSystDefs.push_back(new SystematicDef(*rhs.fMVSystDefs[imv]));
550  this->fMVSignal [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
551  this->fMVBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
552  this->fMVSelectedSignal[this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
553  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
554  this->fMVSelected [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
555  }
556  return *this;
557  }
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
const HistAxis * fOptiHistAxis
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
GenericSystematicDef< caf::SRProxy > SystematicDef
UpDownPair< Spectrum > CopyUpDownSpectrum(const UpDownPair< Spectrum > &copy)
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
Spectrum * fNominalSelectedSignal
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
std::vector< float > Spectrum
Definition: Constants.h:728
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
SystematicDef * fNominalSystDef
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
CutOptimization & ana::CutOptimization::operator= ( CutOptimization &&  rhs)

Definition at line 561 of file CutOptimization.cxx.

References ana::CopyUpDownSpectrum(), fBkgd, fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fNominalBkgd, fNominalSelected, fNominalSelectedBkgd, fNominalSelectedSignal, fNominalSignal, fNominalSystDef, fOptiHistAxis, fSelected, fSelectedBkgd, fSelectedSignal, fSelectionCut, fSignal, fSignalCut, fSystDefs, and submit_hadd::u.

562  {
563  this->fSignalCut = std::move(rhs.fSignalCut);
564  this->fSelectionCut = std::move(rhs.fSignalCut);
565  this->fOptiHistAxis = std::move(rhs.fOptiHistAxis);
566 
567  this->fNominalSystDef = std::move(rhs.fNominalSystDef );
568  this->fNominalSignal = std::move(rhs.fNominalSignal );
569  this->fNominalBkgd = std::move(rhs.fNominalBkgd );
570  this->fNominalSelectedSignal = std::move(rhs.fNominalSelectedSignal);
571  this->fNominalSelectedBkgd = std::move(rhs.fNominalSelectedBkgd );
572  this->fNominalSelected = std::move(rhs.fNominalSelected );
573 
574  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
575  this->fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
576  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSignal .at(rhs.fSystDefs[isyst])));
577  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fBkgd .at(rhs.fSystDefs[isyst])));
578  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst])));
579  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst])));
580  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelected .at(rhs.fSystDefs[isyst])));
581  }
582  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
583  this->fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
584  this->fMVSignal [this->fMVSystDefs.back()] = std::move(rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
585  this->fMVBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
586  this->fMVSelectedSignal[this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
587  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
588  this->fMVSelected [this->fMVSystDefs.back()] = std::move(rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
589  }
590  return *this;
591  }
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
const HistAxis * fOptiHistAxis
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
UpDownPair< Spectrum > CopyUpDownSpectrum(const UpDownPair< Spectrum > &copy)
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
Spectrum * fNominalSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
SystematicDef * fNominalSystDef
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
void ana::CutOptimization::Optimize ( OptimizeMethod_t  method,
FOM_t  fom,
double  data_pot,
TDirectory *  save_dir,
std::string  plot_dump = ".",
bool  debug = false 
)

Definition at line 1361 of file CutOptimization.cxx.

References om::cout, allTimeWatchdog::endl, exit(), kdSigmaOverSigma, and OptimizedSigmaOverSigma().

1364  {
1365  switch(fom) {
1366  case kdSigmaOverSigma: OptimizedSigmaOverSigma(method, data_pot, save_dir, plot_dump, debug);
1367  break;
1368  default:
1369  std::cout << "Your provided FOM is not yet implemented." << std::endl;
1370  exit(0);
1371  }
1372  }
OStream cout
Definition: OStream.cxx:6
void OptimizedSigmaOverSigma(OptimizeMethod_t method, double data_pot, TDirectory *save_dir, std::string plot_dump=".", bool debug=false)
exit(0)
void ana::CutOptimization::OptimizedSigmaOverSigma ( OptimizeMethod_t  method,
double  data_pot,
TDirectory *  save_dir,
std::string  plot_dump = ".",
bool  debug = false 
)
private

Definition at line 879 of file CutOptimization.cxx.

References std::abs(), canvas(), plot_validation_datamc::Clone(), om::cout, Efficiency(), exit(), fBkgd, fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fNominalBkgd, fNominalSelected, fNominalSelectedBkgd, fNominalSelectedSignal, fNominalSignal, genie::utils::style::Format(), fSelected, fSelectedBkgd, fSelectedSignal, fSignal, fSystDefs, ana::Spectrum::GetLabels(), GetName(), Integrate(), kIntegratedDown, kIntegratedInsideBounds, kIntegratedUp, kRed, label, std::max(), maxy, mc_pot, ana::Spectrum::NDimensions(), ana::Spectrum::POT(), cet::pow(), std::sqrt(), string, systs, ToHist(), ToUpDownHist(), submit_hadd::u, and ana::UniqueName().

Referenced by GetAllSystDefs(), and Optimize().

882  {
883  double mc_pot = fNominalSignal->POT();
884  double s = data_pot / mc_pot;
885  // transform the data into histograms
886  // and integrate if method == kIntegratedUp || method == kIntegradedDown || method == kIntegratedInsideBounds
887  TH1 * hNominalSignal = ToHist(fNominalSignal , method, mc_pot ); // scaled to mc for eff
888  TH1 * hNominalSelectedSignal = ToHist(fNominalSelectedSignal, method, mc_pot ); // scaled to mc for eff
889  TH1 * hNominalBkgd = ToHist(fNominalBkgd , method, mc_pot ); // scaled to for mc statistics
890  TH1 * hNominalSelectedBkgd = ToHist(fNominalSelectedBkgd , method, mc_pot ); // scaled to for mc statistics
891  TH1 * hNominalSelected = ToHist(fNominalSelected , method, data_pot); // scaled to data
892  TH1 * hNominalEfficiency = Efficiency(hNominalSelectedSignal, hNominalSignal, method); // scaled to mc
893 
894 
895  // scale the multiverses and calculate efficiencies
896  std::map<SystematicDef*, Multiverse *> fMVEfficiency;
897  for(auto mv = fMVSystDefs.begin(); mv != fMVSystDefs.end(); mv++) {
898  fMVSignal [*mv]->Scale(mc_pot / mc_pot); // for posterity
899  fMVSelectedSignal[*mv]->Scale(mc_pot / mc_pot);
900  fMVBkgd [*mv]->Scale(mc_pot / mc_pot);
901  fMVSelectedBkgd [*mv]->Scale(mc_pot / mc_pot);
902  fMVSelected [*mv]->Scale(data_pot / mc_pot);
903 
904  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
905  Integrate(fMVSignal [*mv], method);
906  Integrate(fMVSelectedSignal[*mv], method);
907  Integrate(fMVBkgd [*mv], method);
908  Integrate(fMVSelectedBkgd [*mv], method);
909  Integrate(fMVSelected [*mv], method);
910  }
911 
912  fMVEfficiency[*mv] = Efficiency(fMVSelectedSignal[*mv], fMVSignal[*mv], method);
913  }
914 
915 
916  std::vector<UpDownPair<TH1 > > hSystSignal;
917  std::vector<UpDownPair<TH1 > > hSystBkgd;
918  std::vector<UpDownPair<TH1 > > hSystSelectedSignal;
919  std::vector<UpDownPair<TH1 > > hSystSelectedBkgd;
920  std::vector<UpDownPair<TH1 > > hSystSelected;
921  std::vector<UpDownPair<TH1 > > hSystEfficiency;
922 
923  std::vector<SystematicDef*> systs;
924  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
925  SystematicDef * syst = fSystDefs[isyst];
926  systs.push_back(syst);
927  hSystSignal .push_back({ToHist(fSignal .at(syst).up, method, mc_pot ), ToHist(fSignal .at(syst).down, method, mc_pot )});
928  hSystSelectedSignal.push_back({ToHist(fSelectedSignal .at(syst).up, method, mc_pot ), ToHist(fSelectedSignal .at(syst).down, method, mc_pot )});
929 
930  hSystBkgd .push_back({ToHist(fBkgd .at(syst).up, method, mc_pot ), ToHist(fBkgd .at(syst).down, method, mc_pot )});
931  hSystSelectedBkgd .push_back({ToHist(fSelectedBkgd .at(syst).up, method, mc_pot ), ToHist(fSelectedBkgd .at(syst).down, method, mc_pot )});
932  hSystSelected .push_back({ToHist(fSelected .at(syst).up, method, data_pot), ToHist(fSelected .at(syst).down, method, data_pot)});
933  hSystEfficiency .push_back(Efficiency(hSystSelectedSignal.back(), hSystSignal.back(), method));
934  }
935  for(auto imv = 0u; imv < fMVSystDefs.size(); imv++) {
936  SystematicDef * mv = fMVSystDefs[imv];
937  systs.push_back(mv);
938 
939  hSystSignal .push_back(ToUpDownHist(fMVSignal [mv], hNominalSignal ));
940  hSystSelectedSignal.push_back(ToUpDownHist(fMVSelectedSignal [mv], hNominalSelectedSignal));
941 
942  hSystBkgd .push_back(ToUpDownHist(fMVBkgd [mv], hNominalBkgd ));
943  hSystSelectedBkgd .push_back(ToUpDownHist(fMVSelectedBkgd [mv], hNominalSelectedBkgd ));
944  hSystSelected .push_back(ToUpDownHist(fMVSelected [mv], hNominalSelected ));
945  hSystEfficiency .push_back(ToUpDownHist(fMVEfficiency [mv], hNominalEfficiency ));
946  }
947 
948  TH1 * hFracUncertaintySelected = (TH1*) hNominalSelected ->Clone(UniqueName().c_str());
949  TH1 * hFracUncertaintyBkgdStat = (TH1*) hNominalSelectedBkgd ->Clone(UniqueName().c_str());
950  TH1 * hFracUncertaintyBkgdSyst = (TH1*) hNominalBkgd ->Clone(UniqueName().c_str());
951  TH1 * hFracUncertaintyEfficiency = (TH1*) hNominalSignal ->Clone(UniqueName().c_str());
952  TH1 * hFracUncertaintyXSec = (TH1*) hNominalSignal ->Clone(UniqueName().c_str());
953 
954  hFracUncertaintySelected ->Scale(0);
955  hFracUncertaintyBkgdStat ->Scale(0);
956  hFracUncertaintyBkgdSyst ->Scale(0);
957  hFracUncertaintyEfficiency ->Scale(0);
958  hFracUncertaintyXSec ->Scale(0);
959 
960  std::vector<TH1*> vabs_uncert_eff_syst (systs.size());
961  std::vector<TH1*> vabs_uncert_bkgd_syst (systs.size());
962  std::vector<TH1*> vfrac_uncert_eff_syst (systs.size());
963  std::vector<TH1*> vfrac_uncert_bkgd_syst(systs.size());
964  TH1* sel_minus_bkgd = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
965  sel_minus_bkgd->Scale(0);
966  for(auto isyst = 0u; isyst < systs.size(); isyst++) {
967  vabs_uncert_eff_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
968  vabs_uncert_bkgd_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
969  vfrac_uncert_eff_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
970  vfrac_uncert_bkgd_syst[isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
971 
972  if(hNominalSelected->GetDimension() == 1) {
973  vfrac_uncert_bkgd_syst[isyst] ->GetYaxis()->SetTitle("#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
974  vfrac_uncert_eff_syst [isyst] ->GetYaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}");
975  vabs_uncert_bkgd_syst [isyst] ->GetYaxis()->SetTitle("s#upoint#delta N^{syst}_{bkgd}");
976  vabs_uncert_eff_syst [isyst] ->GetYaxis()->SetTitle("#delta#varepsilon");
977  }
978 
979  else if(hNominalSelected->GetDimension() == 2) {
980  vfrac_uncert_bkgd_syst[isyst] ->GetZaxis()->SetTitle("#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
981  vfrac_uncert_eff_syst [isyst] ->GetZaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}");
982  vabs_uncert_bkgd_syst [isyst] ->GetZaxis()->SetTitle("s#upoint#delta N^{syst}_{bkgd}");
983  vabs_uncert_eff_syst [isyst] ->GetZaxis()->SetTitle("#delta#varepsilon");
984  }
985  vabs_uncert_eff_syst [isyst] ->Scale(0);
986  vabs_uncert_bkgd_syst [isyst] ->Scale(0);
987  vfrac_uncert_eff_syst [isyst] ->Scale(0);
988  vfrac_uncert_bkgd_syst[isyst] ->Scale(0);
989  }
990 
991  auto NX = hNominalSignal->GetNbinsX();
992  auto NY = hNominalSignal->GetNbinsY();
993  auto NZ = hNominalSignal->GetNbinsZ();
994  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
995  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
996  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
997  double abs_uncert_bkgd_syst = 0;
998  double abs_uncert_eff = 0;
999 
1000  double nom_sel = hNominalSelected ->GetBinContent(ibinx, ibiny, ibinz); // scaled to data
1001  double nom_bkgd = hNominalSelectedBkgd ->GetBinContent(ibinx, ibiny, ibinz); // scaled to mc
1002  double nom_eff = hNominalEfficiency ->GetBinContent(ibinx, ibiny, ibinz);
1003  double nom_signal = nom_sel - s * nom_bkgd;
1004 
1005  sel_minus_bkgd->SetBinContent(ibinx, ibiny, ibinz, nom_signal);
1006  for(auto isyst = 0u; isyst < systs.size(); isyst++) {
1007  double syst_bkgd = 0;
1008  double syst_eff = 0;
1009  // if two sided, take largest shift
1010  if(systs[isyst]->IsTwoSided()) {
1011 
1012  double up_bkgd = hSystSelectedBkgd[isyst].up ->GetBinContent(ibinx, ibiny, ibinz);
1013  double down_bkgd = hSystSelectedBkgd[isyst].down->GetBinContent(ibinx, ibiny, ibinz);
1014 
1015  double up_eff = hSystEfficiency [isyst].up ->GetBinContent(ibinx, ibiny, ibinz);
1016  double down_eff = hSystEfficiency [isyst].down->GetBinContent(ibinx, ibiny, ibinz);
1017 
1018  syst_bkgd = std::max(std::abs(up_bkgd - nom_bkgd), std::abs(down_bkgd - nom_bkgd));
1019  syst_eff = std::max(std::abs(up_eff - nom_eff ), std::abs(down_eff - nom_eff ));
1020  }
1021  else {
1022  syst_bkgd = std::abs(hSystSelectedBkgd[isyst].up ->GetBinContent(ibinx, ibiny, ibinz) - nom_bkgd);
1023  syst_eff = std::abs(hSystEfficiency [isyst].up ->GetBinContent(ibinx, ibiny, ibinz) - nom_eff );
1024  }
1025 
1026  // scale for MC statistics
1027  syst_bkgd = s * syst_bkgd;
1028 
1029  abs_uncert_bkgd_syst += std::pow(syst_bkgd, 2);
1030  abs_uncert_eff += std::pow(syst_eff , 2);
1031 
1032  vabs_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff );
1033  vabs_uncert_bkgd_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd);
1034  if(nom_signal)
1035  vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd / nom_signal);
1036  else
1037  vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);
1038  if(nom_eff)
1039  vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff / nom_eff);
1040  else
1041  vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);
1042  } // end loop over systematics
1043 
1044  // calculate fractional uncertainty on efficiency
1045  double frac_uncert_sel = 100;
1046  double frac_uncert_bkgd_stat = 100;
1047  double frac_uncert_bkgd_syst = 100;
1048  double frac_uncert_eff = 100;
1049  if(nom_signal) {
1050  frac_uncert_sel = std::sqrt(nom_sel) / nom_signal;
1051  frac_uncert_bkgd_stat = s * std::sqrt(nom_bkgd) / nom_signal;
1052  frac_uncert_bkgd_syst = std::sqrt(abs_uncert_bkgd_syst) / nom_signal;
1053  }
1054  if(nom_eff)
1055  frac_uncert_eff = std::sqrt(abs_uncert_eff) / nom_eff;
1056 
1057  // finally calculate FOM
1058  double frac_uncert_sigma = std::sqrt(std::pow(frac_uncert_sel , 2) +
1059  std::pow(frac_uncert_bkgd_stat, 2) +
1060  std::pow(frac_uncert_bkgd_syst, 2) +
1061  std::pow(frac_uncert_eff , 2) );
1062 
1063 
1064  // fill histograms
1065  hFracUncertaintyXSec ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sigma );
1066  hFracUncertaintySelected ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sel );
1067  hFracUncertaintyBkgdStat ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_stat);
1068  hFracUncertaintyBkgdSyst ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_syst);
1069  hFracUncertaintyEfficiency ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_eff );
1070  }// end loop over ibinz
1071  }// end loop over ibiny
1072  }// end loop over ibinx
1073 
1074 
1075  if(debug) {
1076  //////////////////////////////////// DEBUG //////////////////////////////////
1077  // save internal histograms to file for later plotting
1078  TFile * fdebug = new TFile((plot_dump + "/debug.root").c_str(), "recreate");
1079 
1080  hNominalSelectedBkgd->Write("hNominalSelectedBkgd");
1081  hNominalSelected->Write("hNominalSelected");
1082  hNominalSelectedSignal->Write("hNominalSelectedSignal");
1083  hNominalSignal->Write("hNominalSignal");
1084  hNominalBkgd->Write("hNominalBkgd");
1085  hNominalEfficiency->Write("hNominalEfficiency");
1086 
1087  hFracUncertaintyXSec ->Write("hFracUncertaintyXSec" );
1088  hFracUncertaintySelected ->Write("hFracUncertaintySelected" );
1089  hFracUncertaintyBkgdStat ->Write("hFracUncertaintyBkgdStat" );
1090  hFracUncertaintyBkgdSyst ->Write("hFracUncertaintyBkgdSyst" );
1091  hFracUncertaintyEfficiency ->Write("hFracUncertaintyEfficiency" );
1092 
1093 
1094  for(auto isyst = 0u; isyst < systs.size(); isyst++ ) {
1095  hSystSelectedBkgd [isyst].up->Write((systs[isyst]->GetName() + "_hSystSelectedBkgd" + "_UP").c_str());
1096  hSystSelected [isyst].up->Write((systs[isyst]->GetName() + "_hSystSelected" + "_UP").c_str());
1097  hSystSelectedSignal[isyst].up->Write((systs[isyst]->GetName() + "_hSystSelectedSignal" + "_UP").c_str());
1098  hSystBkgd [isyst].up->Write((systs[isyst]->GetName() + "_hSystBkgd" + "_UP").c_str());
1099  hSystSignal [isyst].up->Write((systs[isyst]->GetName() + "_hSystSignal" + "_UP").c_str());
1100  hSystEfficiency [isyst].up->Write((systs[isyst]->GetName() + "_hSystEfficiency" + "_UP").c_str());
1101 
1102  if(hSystSelectedBkgd[isyst].down != NULL) {
1103  hSystSelectedBkgd [isyst].down->Write((systs[isyst]->GetName() + "_hSystSelectedBkgd" + "_DOWN").c_str());
1104  hSystSelected [isyst].down->Write((systs[isyst]->GetName() + "_hSystSelected" + "_DOWN").c_str());
1105  hSystSelectedSignal[isyst].down->Write((systs[isyst]->GetName() + "_hSystSelectedSignal" + "_DOWN").c_str());
1106  hSystBkgd [isyst].down->Write((systs[isyst]->GetName() + "_hSystBkgd" + "_DOWN").c_str());
1107  hSystSignal [isyst].down->Write((systs[isyst]->GetName() + "_hSystSignal" + "_DOWN").c_str());
1108  hSystEfficiency [isyst].down->Write((systs[isyst]->GetName() + "_hSystEfficiency" + "_DOWN").c_str());
1109  }
1110 
1111  vabs_uncert_eff_syst [isyst]->Write((systs[isyst]->GetName() + "_abs_uncert_eff_syst" ).c_str());
1112  vabs_uncert_bkgd_syst [isyst]->Write((systs[isyst]->GetName() + "_abs_uncert_bkgd_syst" ).c_str());
1113  vfrac_uncert_eff_syst [isyst]->Write((systs[isyst]->GetName() + "_frac_uncert_eff_syst" ).c_str());
1114  vfrac_uncert_bkgd_syst[isyst]->Write((systs[isyst]->GetName() + "_frac_uncert_bkgd_syst").c_str());
1115 
1116 
1117 // debug_syst_sel_bkgd_up .push_back(hSystSelectedBkgd [isyst].up);
1118 // debug_syst_sel_bkgd_down .push_back(hSystSelectedBkgd [isyst].down);
1119 // debug_syst_sel_up .push_back(hSystSelected [isyst].up);
1120 // debug_syst_sel_down .push_back(hSystSelected [isyst].down);
1121 // debug_syst_sel_signal_up .push_back(hSystSelectedSignal[isyst].up);
1122 // debug_syst_sel_signal_down .push_back(hSystSelectedSignal[isyst].down);
1123 
1124 // debug_syst_bkgd_up .push_back(hSystBkgd[isyst].up);
1125 // debug_syst_bkgd_down .push_back(hSystBkgd[isyst].down);
1126 
1127  // debug_syst_signal_up .push_back(hSystSignal[isyst].up);
1128  // debug_syst_signal_down .push_back(hSystSignal[isyst].down);
1129 
1130  // debug_syst_eff_up .push_back(hSystEfficiency[isyst].up);
1131  // debug_syst_eff_down .push_back(hSystEfficiency[isyst].down);
1132 
1133  // debug_syst_bkgd_labels.push_back(systs[isyst]->GetName());
1134  // debug_syst_bkgd_labels_no_nominal.push_back(systs[isyst]->GetName());
1135  }
1136 
1137 
1138  /*
1139  PlotDebug(debug_syst_bkgd_up , debug_syst_bkgd_labels, "hist", "Background +1#sigma" , plot_dump + "/debug_syst_bkgd_up.pdf" );
1140  PlotDebug(debug_syst_bkgd_down , debug_syst_bkgd_labels, "hist", "Background -1#sigma" , plot_dump + "/debug_syst_bkgd_down.pdf" );
1141  PlotDebug(debug_syst_sel_up , debug_syst_bkgd_labels, "hist", "Selected +1#sigma" , plot_dump + "/debug_syst_sel_up.pdf" );
1142  PlotDebug(debug_syst_sel_down , debug_syst_bkgd_labels, "hist", "Selected -1#sigma" , plot_dump + "/debug_syst_sel_down.pdf" );
1143  PlotDebug(debug_syst_sel_bkgd_up , debug_syst_bkgd_labels, "hist", "Selected Background +1#sigma", plot_dump + "/debug_syst_sel_bkgd_up.pdf" );
1144  PlotDebug(debug_syst_sel_bkgd_down , debug_syst_bkgd_labels, "hist", "Selected Background -1#sigma", plot_dump + "/debug_syst_sel_bkgd_down.pdf" );
1145  PlotDebug(debug_syst_sel_signal_up , debug_syst_bkgd_labels, "hist", "Selected Signal +1#sigma" , plot_dump + "/debug_syst_sel_signal_up.pdf" );
1146  PlotDebug(debug_syst_sel_signal_down, debug_syst_bkgd_labels, "hist", "Selected Signal -1#sigma" , plot_dump + "/debug_syst_sel_signal_down.pdf");
1147  PlotDebug(debug_syst_signal_up , debug_syst_bkgd_labels, "hist", "Signal +1#sigma" , plot_dump + "/debug_syst_signal_up.pdf" );
1148  PlotDebug(debug_syst_signal_down , debug_syst_bkgd_labels, "hist", "Signal -1#sigma" , plot_dump + "/debug_syst_signal_down.pdf" );
1149  PlotDebug(debug_syst_eff_up , debug_syst_bkgd_labels, "hist", "Efficiency +1#sigma" , plot_dump + "/debug_syst_eff_up.pdf" );
1150  PlotDebug(debug_syst_eff_down , debug_syst_bkgd_labels, "hist", "Efficiency -1#sigma" , plot_dump + "/debug_syst_eff_down.pdf" );
1151 
1152  PlotDebug(vabs_uncert_eff_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Abs. Syst. Uncert Efficiency" , plot_dump + "/debug_syst_abs_uncert_eff.pdf" );
1153  PlotDebug(vabs_uncert_bkgd_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Abs. Syst. Uncert Background" , plot_dump + "/debug_syst_abs_uncert_bkgd.pdf" );
1154  PlotDebug(vfrac_uncert_eff_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Frac. Syst. Uncert Efficiency", plot_dump + "/debug_syst_frac_uncert_eff.pdf" , 2);
1155  PlotDebug(vfrac_uncert_bkgd_syst, debug_syst_bkgd_labels_no_nominal, "hist", "Frac. Syst. Uncert Background", plot_dump + "/debug_syst_frac_uncert_bkgd.pdf" , 2);
1156  */
1157 
1158  for(auto syst : fMVSystDefs) {
1159 
1160  fMVSignal [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Signal" );
1161  fMVBkgd [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Bkgd" );
1162  fMVSelected [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Selected" );
1163  fMVSelectedSignal[syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_SelectedSignal");
1164  fMVSelectedBkgd [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_SelectedBkgd" );
1165  fMVEfficiency [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Efficiency" );
1166 
1167 // const std::vector<TH1*> mv_signal_hists = fMVSignal [syst]->GetUniversesHist(kCyan);
1168 // const std::vector<TH1*> mv_bkgd_hists = fMVBkgd [syst]->GetUniversesHist(kCyan);
1169 // const std::vector<TH1*> mv_sel_hists = fMVSelected [syst]->GetUniversesHist(kCyan);
1170 // const std::vector<TH1*> mv_sel_signal_hists = fMVSelectedSignal[syst]->GetUniversesHist(kCyan);
1171 // const std::vector<TH1*> mv_sel_bkgd_hists = fMVSelectedBkgd [syst]->GetUniversesHist(kCyan);
1172 
1173 // std::string name = syst->GetName();
1174 
1175 // PlotDebug(mv_signal_hists , hNominalSignal , name + " Multiverse: Signal" , plot_dump + "/debug_" + name + "_signal.pdf" );
1176 // PlotDebug(mv_bkgd_hists , hNominalBkgd , name + " Multiverse: Bkgd" , plot_dump + "/debug_" + name + "_bkgd.pdf" );
1177 // PlotDebug(mv_sel_signal_hists, hNominalSelectedSignal, name + " Multiverse: SelectedSignal", plot_dump + "/debug_" + name + "_sel_signal.pdf");
1178 // PlotDebug(mv_sel_bkgd_hists , hNominalSelectedBkgd , name + " Multiverse: SelectedBkgd" , plot_dump + "/debug_" + name + "_sel_bkgd.pdf" );
1179 // PlotDebug(mv_sel_hists , hNominalSelected , name + " Multiverse: Selected" , plot_dump + "/debug_" + name + "_sel.pdf" );
1180  }
1181 
1182  /*
1183  PlotDebug(hNominalSignal , "hist", "hNominalSignal" , plot_dump + "/debug_nom_sig.pdf");
1184  PlotDebug(hNominalSelectedSignal , "hist", "hNominalSelectedSignal", plot_dump + "/debug_nom_sig_sel.pdf");
1185  PlotDebug(hNominalBkgd , "hist", "hNominalBkgd" , plot_dump + "/debug_nom_bkgd.pdf");
1186  PlotDebug(hNominalSelectedBkgd , "hist", "hNominalSelectedBkgd" , plot_dump + "/debug_nom_bkgd_sel.pdf");
1187  PlotDebug(hNominalSelected , "hist", "hNominalSelected" , plot_dump + "/debug_nom_sel.pdf");
1188  PlotDebug(hNominalEfficiency , "hist", "hNominalEfficiency" , plot_dump + "/debug_nom_eff.pdf");
1189 
1190  PlotDebug(sel_minus_bkgd , "hist", "Nominal Selected - Background", plot_dump + "/debug_nom_sel_minus_bkgd.pdf");
1191  */
1192  fdebug->Close();
1193  //////////////////////////////////// END DEBUG //////////////////////////////////
1194  }
1195 
1196  // plot histograms
1197  double maxy = 1;
1198  hFracUncertaintyXSec ->SetTitle("Fractional Uncertainty on Cross-section" );
1199  hFracUncertaintySelected ->SetTitle("Fractional Stat. Uncertainty on Selection" );
1200  hFracUncertaintyBkgdStat ->SetTitle("Fractional Stat. Uncertainty on Background");
1201  hFracUncertaintyBkgdSyst ->SetTitle("Fractional Syst. Uncertainty on Background");
1202  hFracUncertaintyEfficiency ->SetTitle("Fractional Syst. Uncertainty on Efficiency");
1203 
1204  if(fNominalSignal->NDimensions() == 1 && method != kIntegratedInsideBounds) {
1205 
1206  hFracUncertaintyXSec ->GetYaxis()->SetRangeUser(0, maxy);
1207  hFracUncertaintySelected ->GetYaxis()->SetRangeUser(0, maxy);
1208  hFracUncertaintyBkgdStat ->GetYaxis()->SetRangeUser(0, maxy);
1209  hFracUncertaintyBkgdSyst ->GetYaxis()->SetRangeUser(0, maxy);
1210  hFracUncertaintyEfficiency ->GetYaxis()->SetRangeUser(0, maxy);
1211 
1212  hFracUncertaintyXSec ->GetYaxis()->SetTitle("#frac{#delta#sigma}{#sigma}" );
1213  hFracUncertaintySelected ->GetYaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{N_{sel}-s#upoint N_{bkg}}" );
1214  hFracUncertaintyBkgdStat ->GetYaxis()->SetTitle("#frac{s#upoint#sqrt{N_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
1215  hFracUncertaintyBkgdSyst ->GetYaxis()->SetTitle("#frac{#sqrt{s#upoint#deltaN^{syst}_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
1216  hFracUncertaintyEfficiency ->GetYaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}" );
1217 
1218  hFracUncertaintyXSec ->GetYaxis()->CenterTitle();
1219  hFracUncertaintySelected ->GetYaxis()->CenterTitle();
1220  hFracUncertaintyBkgdStat ->GetYaxis()->CenterTitle();
1221  hFracUncertaintyBkgdSyst ->GetYaxis()->CenterTitle();
1222  hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();
1223 
1224  hFracUncertaintyXSec ->SetLineColor(kRed);
1225  hFracUncertaintySelected ->SetLineColor(kRed);
1226  hFracUncertaintyBkgdStat ->SetLineColor(kRed);
1227  hFracUncertaintyBkgdSyst ->SetLineColor(kRed);
1228  hFracUncertaintyEfficiency ->SetLineColor(kRed);
1229 
1230  // plot
1231  TCanvas * canvas = new TCanvas();
1232  canvas->SetLeftMargin(0.15);
1233  canvas->SetRightMargin(0.01);
1234  canvas->cd();
1236  hFracUncertaintyXSec ->Draw("hist");
1237  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_xsec.pdf",
1238  plot_dump.c_str(),
1239  label.c_str()));
1240  hFracUncertaintySelected ->Draw("hist");
1241  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_sel.pdf",
1242  plot_dump.c_str(),
1243  label.c_str()));
1244  hFracUncertaintyBkgdStat ->Draw("hist");
1245  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_bkgd_stat.pdf",
1246  plot_dump.c_str(),
1247  label.c_str()));
1248  hFracUncertaintyBkgdSyst ->Draw("hist");
1249  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_bkgd_syst.pdf",
1250  plot_dump.c_str(),
1251  label.c_str()));
1252  hFracUncertaintyEfficiency ->Draw("hist");
1253  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_eff.pdf",
1254  plot_dump.c_str(),
1255  label.c_str()));
1256 
1257 
1258  }
1259  else if(fNominalSignal->NDimensions() == 2 || method == kIntegratedInsideBounds) {
1260  hFracUncertaintyXSec ->GetZaxis()->SetRangeUser(0, maxy);
1261  hFracUncertaintySelected ->GetZaxis()->SetRangeUser(0, maxy);
1262  hFracUncertaintyBkgdStat ->GetZaxis()->SetRangeUser(0, maxy);
1263  hFracUncertaintyBkgdSyst ->GetZaxis()->SetRangeUser(0, maxy);
1264  hFracUncertaintyEfficiency ->GetZaxis()->SetRangeUser(0, maxy);
1265 
1266  hFracUncertaintyXSec ->GetZaxis()->SetTitle("#frac{#delta#sigma}{#sigma}" );
1267  hFracUncertaintySelected ->GetZaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}" );
1268  hFracUncertaintyBkgdStat ->GetZaxis()->SetTitle("#frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}" );
1269  hFracUncertaintyBkgdSyst ->GetZaxis()->SetTitle("#frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}");
1270  hFracUncertaintyEfficiency ->GetZaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}" );
1271 
1272  hFracUncertaintyXSec ->GetYaxis()->CenterTitle();
1273  hFracUncertaintySelected ->GetYaxis()->CenterTitle();
1274  hFracUncertaintyBkgdStat ->GetYaxis()->CenterTitle();
1275  hFracUncertaintyBkgdSyst ->GetYaxis()->CenterTitle();
1276  hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();
1277 
1278  // plot
1279  TCanvas * canvas = new TCanvas();
1280  canvas->cd();
1281  std::string labelx = fNominalSignal->GetLabels()[0];
1282  std::string labely = method == kIntegratedInsideBounds? labelx : fNominalSignal->GetLabels()[1];
1283  hFracUncertaintyXSec ->Draw("colz");
1284  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_xsec.pdf",
1285  plot_dump.c_str(),
1286  labelx.c_str(),
1287  labely.c_str()));
1288  hFracUncertaintySelected ->Draw("colz");
1289  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_sel.pdf",
1290  plot_dump.c_str(),
1291  labelx.c_str(),
1292  labely.c_str()));
1293  hFracUncertaintyBkgdStat ->Draw("colz");
1294  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_bkgd_stat.pdf",
1295  plot_dump.c_str(),
1296  labelx.c_str(),
1297  labely.c_str()));
1298  hFracUncertaintyBkgdSyst ->Draw("colz");
1299  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_bkgd_syst.pdf",
1300  plot_dump.c_str(),
1301  labelx.c_str(),
1302  labely.c_str()));
1303  hFracUncertaintyEfficiency ->Draw("colz");
1304  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_eff.pdf",
1305  plot_dump.c_str(),
1306  labelx.c_str(),
1307  labely.c_str()));
1308  }
1309  else {
1310  std::cout << "Don't know how to draw a " << fNominalSignal->NDimensions() << " dimensional histogram.";
1311  exit(1);
1312  }
1313 
1314  // save
1315  save_dir->cd();
1316  hFracUncertaintyXSec ->Write("frac_uncert_xsec" );
1317  hFracUncertaintySelected ->Write("frac_uncert_sel" );
1318  hFracUncertaintyBkgdStat ->Write("frac_uncert_bkgd_stat");
1319  hFracUncertaintyBkgdSyst ->Write("frac_uncert_bkgd_syst");
1320  hFracUncertaintyEfficiency ->Write("frac_uncert_eff" );
1321  }
T max(const caf::Proxy< T > &a, T b)
UpDownPair< TH1 > ToUpDownHist(Multiverse *mv, const TH1 *h_nominal)
enum BeamMode kRed
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
double maxy
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
const double mc_pot
std::map< SystematicDef *, Multiverse * > fMVBkgd
static TH1 * ToHist(const Spectrum *spec, OptimizeMethod_t method, double POT=-1)
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
TCanvas * canvas(const char *nm, const char *ti, const char *a)
Definition: pass1_plots.C:36
GenericSystematicDef< caf::SRProxy > SystematicDef
float abs(float number)
Definition: d0nt_math.hpp:39
const char * label
const XML_Char * s
Definition: expat.h:262
std::string GetName(int i)
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
static TH1 * Efficiency(TH1 *num, TH1 *denom, OptimizeMethod_t method)
Spectrum * fNominalSelectedSignal
unsigned int NDimensions() const
Definition: Spectrum.h:262
static void Integrate(TH1 *&result, const TH1 *tmp, OptimizeMethod_t method)
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
double POT() const
Definition: Spectrum.h:227
OStream cout
Definition: OStream.cxx:6
exit(0)
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:263
std::map< SystematicDef *, Multiverse * > fMVSelected
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
std::vector< SystematicDef * > fMVSystDefs
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
enum BeamMode string
void ana::CutOptimization::PlotDebug ( TH1 *  hist,
std::string  draw_option,
std::string  title,
std::string  name 
)
static

Definition at line 15 of file CutOptimization.cxx.

References plot_validation_datamc::c, and tmp.

Referenced by fiducial_optimization_plots(), muonid_optimization_plots(), and PlotDebug().

18  {
19  TDirectory * tmp = gDirectory;
20  TCanvas * c = new TCanvas();
21  c->SetLeftMargin(0.15);
22  c->SetRightMargin(0.01);
23 
24 
25  c->cd();
26 
27  hist->Draw(draw_option.c_str());
28  hist->SetTitle(title.c_str());
29  c->Print(name.c_str());
30 
31  tmp->cd();
32  }
const XML_Char * name
Definition: expat.h:151
Float_t tmp
Definition: plot.C:36
void ana::CutOptimization::PlotDebug ( std::vector< TH1 * >  hists,
std::vector< std::string labels,
std::string  draw_option,
std::string  title,
std::string  name,
double  ymax = -1 
)
static

Definition at line 35 of file CutOptimization.cxx.

References plot_validation_datamc::c, colors, make_mec_shifts_plots::GetMaximum(), kBlue, kGreen, kViolet, MECModelEnuComparisons::leg, std::max(), tmp, and submit_hadd::u.

40  {
41  TDirectory * tmp = gDirectory;
42  TCanvas * c = new TCanvas();
43  c->SetLeftMargin(0.15);
44  c->SetRightMargin(0.01);
45 
46  c->cd();
47  int colors[] = {kBlue, kViolet, kMagenta,
48  kPink, kGreen, kTeal,
49  kCyan, kAzure, kGray,
50  kGray+1, kGray+3};
51 
52  TLegend * leg = new TLegend(0.4,0.7,0.7,0.9);
53 
54  double max_val = -1e9;
55  for(auto ihist = 0u; ihist < hists.size(); ihist++) {
56  if(hists[ihist] == NULL) continue;
57  max_val = std::max(max_val, hists[ihist]->GetMaximum());
58  }
59 
60  int drawn = 0;
61  for(auto ihist = 0u; ihist < hists.size(); ihist++) {
62  // pointer could be null
63  if(hists[ihist] == NULL) continue;
64 
65  // if valid pointer, style hist
66  hists[ihist]->SetLineColor(colors[ihist]);
67  hists[ihist]->SetTitle(title.c_str());
68 
69  leg->AddEntry(hists[ihist], labels[ihist].c_str(), "l");
70  // use the axes of first valid histogram
71  if(drawn == 0) {
72  // adjust range so all hists will show up
73  if(ymax < 0)
74  hists[ihist]->GetYaxis()->SetRangeUser(0, max_val * 1.2);
75  else
76  hists[ihist]->GetYaxis()->SetRangeUser(0, ymax);
77  hists[ihist]->Draw(draw_option.c_str());
78  }
79  else {
80  hists[ihist]->Draw((draw_option + " same").c_str());
81  }
82  drawn++;
83  }
84 
85  leg->Draw("same");
86  c->Print(name.c_str());
87 
88  tmp->cd();
89  }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
Float_t tmp
Definition: plot.C:36
TString hists[nhists]
Definition: bdt_com.C:3
Double_t ymax
Definition: plot.C:25
int colors[6]
Definition: tools.h:1
enum BeamMode kViolet
enum BeamMode kGreen
enum BeamMode kBlue
void ana::CutOptimization::PlotDebug ( const std::vector< TH1 * >  universes,
Spectrum nominal,
std::string  title,
std::string  name 
)
static

Definition at line 92 of file CutOptimization.cxx.

References PlotDebug(), ana::Spectrum::POT(), and ana::Spectrum::ToTH1().

96  {
97  TH1 * h_nominal = nominal->ToTH1(nominal->POT());
98  PlotDebug(universes, h_nominal, title, name);
99  }
const XML_Char * name
Definition: expat.h:151
static void PlotDebug(TH1 *hist, std::string draw_option, std::string title, std::string name)
void ana::CutOptimization::PlotDebug ( const std::vector< TH1 * >  universes,
TH1 *  h_nominal,
std::string  title,
std::string  name 
)
static

Definition at line 102 of file CutOptimization.cxx.

References plot_validation_datamc::c, make_mec_shifts_plots::GetMaximum(), std::max(), tmp, and submit_hadd::u.

106  {
107  TDirectory * tmp = gDirectory;
108  TCanvas * c = new TCanvas();
109  c->SetLeftMargin(0.15);
110  c->SetRightMargin(0.01);
111  c->cd();
112 
113  double max_val = -1e9;
114  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
115  max_val = std::max(max_val, universes[ihist]->GetMaximum());
116  }
117  h_nominal->GetYaxis()->SetRangeUser(0, max_val * 1.2);
118  h_nominal->SetTitle(title.c_str());
119  h_nominal->Draw("hist");
120 
121  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
122  universes[ihist]->Draw("hist same");
123  }
124 
125  h_nominal->Draw("hist same");
126 
127 
128 
129  c->Print(name.c_str());
130 
131  tmp->cd();
132  }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
Float_t tmp
Definition: plot.C:36
void ana::CutOptimization::SaveTo ( TDirectory *  dir,
const std::string name 
) const

Definition at line 277 of file CutOptimization.cxx.

References dir, fBkgd, fMVBkgd, fMVSelected, fMVSelectedBkgd, fMVSelectedSignal, fMVSignal, fMVSystDefs, fNominalBkgd, fNominalSelected, fNominalSelectedBkgd, fNominalSelectedSignal, fNominalSignal, genie::utils::style::Format(), fSelected, fSelectedBkgd, fSelectedSignal, fSignal, fSystDefs, GetName(), make_syst_table_plots::nsysts, ana::Spectrum::SaveTo(), ana::SaveToUpDownSpectra(), tmp, and submit_hadd::u.

Referenced by containment_optimization(), demo_cut_optimization(), fiducial_opt(), fiducial_optimization(), muonid_opt(), muonid_optimization(), and prongcvn_optimization().

278  {
279  TDirectory * tmp = gDirectory;
280 
281  dir = dir->mkdir(name.c_str()); // switch to subdir
282  dir->cd();
283 
284  TObjString("CutOptimization").Write("type");
285 
286  TDirectory * nom_dir = dir->mkdir("Nominal");
287  // save nominal
288  fNominalSignal ->SaveTo(nom_dir, "Signal");
289  fNominalBkgd ->SaveTo(nom_dir, "Bkgd");
290  fNominalSelectedSignal->SaveTo(nom_dir, "SelectedSignal");
291  fNominalSelectedBkgd ->SaveTo(nom_dir, "SelectedBkgd");
292  fNominalSelected ->SaveTo(nom_dir, "Selected");
293 
294  dir->cd();
295 
296  TVectorD nsysts(1);
297  nsysts[0] = fSystDefs.size();
298  nsysts.Write("nsysts");
299 
300  // for each syst def, save up down pairs
301  // Iterating over the vector insures systematics will be saved in the same
302  // order everytime. Critical for hadd_cafana
303  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
304  // save systdefs
305  fSystDefs[isyst]->SaveTo(dir, TString::Format("syst_defs/%d", isyst).Data());
306 
307  TDirectory * syst_dir = dir->mkdir(TString::Format("syst_%s",fSystDefs[isyst]->GetName().c_str()).Data());
308  SaveToUpDownSpectra( fSignal .at(fSystDefs[isyst]), syst_dir->mkdir("fSignal" ));
309  SaveToUpDownSpectra( fBkgd .at(fSystDefs[isyst]), syst_dir->mkdir("fBkgd" ));
310  SaveToUpDownSpectra( fSelectedSignal.at(fSystDefs[isyst]), syst_dir->mkdir("fSelectedSignal"));
311  SaveToUpDownSpectra( fSelectedBkgd .at(fSystDefs[isyst]), syst_dir->mkdir("fSelectedBkgd" ));
312  SaveToUpDownSpectra( fSelected .at(fSystDefs[isyst]), syst_dir->mkdir("fSelected" ));
313  }
314 
315 
316  TVectorD nmultiverse(1);
317  nsysts[0] = fMVSystDefs.size();
318  nsysts.Write("nmultiverse");
319 
320  // for each mv syst def, save multiverse
321  // Iterating over the vector insures multiverses will be saved in the same
322  // order everytime. Critical for hadd_cafana
323  for(auto imv = 0u; imv < fMVSystDefs.size(); imv++) {
324  // save mv syst_defs
325  fMVSystDefs[imv]->SaveTo(dir, TString::Format("mv_defs/%d", imv).Data());
326 
327  TDirectory * mv_dir = dir->mkdir(TString::Format("mv_%s",fMVSystDefs[imv]->GetName().c_str()).Data());
328  fMVSignal .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSignal" );
329  fMVBkgd .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVBkgd" );
330  fMVSelectedSignal.at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelectedSignal");
331  fMVSelectedBkgd .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelectedBkgd" );
332  fMVSelected .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelected" );
333  }
334 
335  dir->Write();
336  delete dir;
337 
338  tmp->cd();
339 
340  }
const XML_Char * name
Definition: expat.h:151
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVBkgd
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
Float_t tmp
Definition: plot.C:36
std::string GetName(int i)
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
Spectrum * fNominalSelectedSignal
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
TDirectory * dir
Definition: macro.C:5
void SaveToUpDownSpectra(const UpDownPair< Spectrum > &save, TDirectory *dir)
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
std::vector< SystematicDef * > fSystDefs
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::map< SystematicDef *, Multiverse * > fMVSelected
std::vector< SystematicDef * > fMVSystDefs
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
void ana::CutOptimization::SetSpillCuts ( const SpillCut spillcut)

Definition at line 175 of file CutOptimization.cxx.

References fNominalSystDef, fSystDefs, ana::GenericSystematicDef< SRType >::Loaders(), ana::GenericSystematicDef< SRType >::LoadersUp(), and submit_hadd::u.

Referenced by containment_optimization(), demo_cut_optimization(), fiducial_optimization(), muonid_optimization(), and prongcvn_optimization().

176  {
177  for(auto iload = 0u; iload < fNominalSystDef->LoadersUp().size(); iload++)
178  fNominalSystDef->Loaders()[iload]->SetSpillCut(spillcut);
179 
180  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
181  for(auto iload = 0u; iload < fSystDefs[isyst]->LoadersUp().size(); iload++)
182  fSystDefs[isyst]->LoadersUp()[iload]->SetSpillCut(spillcut);
183  if(fSystDefs[isyst]->IsTwoSided()) {
184  for(auto iload = 0u; iload < fSystDefs[isyst]->LoadersDown().size(); iload++)
185  fSystDefs[isyst]->LoadersDown()[iload]->SetSpillCut(spillcut);
186  }
187  }
188  }
std::vector< SpectrumLoaderBase * > Loaders() const
std::vector< SystematicDef * > fSystDefs
std::vector< SpectrumLoaderBase * > LoadersUp() const
SystematicDef * fNominalSystDef
TH1 * ana::CutOptimization::ToHist ( const Spectrum spec,
OptimizeMethod_t  method,
double  POT = -1 
)
staticprivate

Definition at line 739 of file CutOptimization.cxx.

References om::cout, allTimeWatchdog::endl, exit(), Integrate(), kIntegratedDown, kIntegratedInsideBounds, kIntegratedUp, ana::Spectrum::NDimensions(), ana::Spectrum::POT(), runNovaSAM::ret, tmp, ana::Spectrum::ToTH1(), ana::Spectrum::ToTH2(), and ana::Spectrum::ToTH3().

Referenced by GetAllSystDefs(), and OptimizedSigmaOverSigma().

740  {
741  TH1 * tmp;
742  TH1 * ret;
743 
744  // this will happen for one sided systs
745  if(spec == NULL) return NULL;
746 
747  // if no POT provided, use the spectrum's POT
748  if(POT < 0) POT = spec->POT();
749 
750  // Call different function based on dimension
751  auto ndims = spec->NDimensions();
752  switch(ndims) {
753  case 1:
754  tmp = spec->ToTH1(POT);
755  break;
756  case 2:
757  tmp = spec->ToTH2(POT);
758  break;
759  case 3:
760  tmp = spec->ToTH3(POT);
761  break;
762  default:
763  std::cout << "Spectrum appears to be empty. Something went terribly wrong here" << std::endl;
764  exit(1);
765  }
766 
767  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
768  ret = (TH1*) tmp->Clone();
769 
770  // integrated result gets saved back to ret
771  Integrate(ret, tmp, method);
772  }
773  else {
774  // default is kBinByBin, just return the histogram then
775  ret = tmp;
776  }
777 
778  return ret;
779  }
Float_t tmp
Definition: plot.C:36
static void Integrate(TH1 *&result, const TH1 *tmp, OptimizeMethod_t method)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
exit(0)
UpDownPair< TH1 > ana::CutOptimization::ToUpDownHist ( Multiverse mv,
const TH1 *  h_nominal 
)
private

Definition at line 716 of file CutOptimization.cxx.

References std::abs(), ana::Multiverse::GetMinusOneSigmaShift(), ana::Multiverse::GetPlusOneSigmaShift(), std::max(), and runNovaSAM::ret.

Referenced by GetAllSystDefs(), and OptimizedSigmaOverSigma().

717  {
718  TH1 * h_up = mv->GetPlusOneSigmaShift (h_nominal);
719  TH1 * h_down = mv->GetMinusOneSigmaShift(h_nominal);
720  // transform into 1 sided shift as the most conservative shift
721  TH1 * ret = (TH1*) h_up->Clone();
722  for(auto ibinx = 1; ibinx <= ret->GetNbinsX(); ibinx++) {
723  for(auto ibiny = 1; ibiny <= ret->GetNbinsY(); ibiny++) {
724  for(auto ibinz = 1; ibinz <= ret->GetNbinsZ(); ibinz++) {
725  double up = h_up ->GetBinContent(ibinx, ibiny, ibinz);
726  double down = h_down ->GetBinContent(ibinx, ibiny, ibinz);
727  double nominal = h_nominal ->GetBinContent(ibinx, ibiny, ibinz);
728  double abs_shift = std::max(std::abs(up - nominal), std::abs(down - nominal));
729  ret->SetBinContent(ibinx, ibiny, ibinz, abs_shift + nominal);
730  }
731  }
732  }
733  return UpDownPair<TH1>{ret};
734  }
T max(const caf::Proxy< T > &a, T b)
float abs(float number)
Definition: d0nt_math.hpp:39

Member Data Documentation

std::map< SystematicDef*, UpDownPair<Spectrum> > ana::CutOptimization::fBkgd
private
std::map< SystematicDef*, Multiverse *> ana::CutOptimization::fMVBkgd
private
std::map< SystematicDef*, Multiverse *> ana::CutOptimization::fMVSelected
private
std::map< SystematicDef*, Multiverse *> ana::CutOptimization::fMVSelectedBkgd
private
std::map< SystematicDef*, Multiverse *> ana::CutOptimization::fMVSelectedSignal
private
std::map< SystematicDef*, Multiverse *> ana::CutOptimization::fMVSignal
private
std::vector<SystematicDef*> ana::CutOptimization::fMVSystDefs
private
Spectrum* ana::CutOptimization::fNominalBkgd
private
Spectrum* ana::CutOptimization::fNominalSelected
private
Spectrum* ana::CutOptimization::fNominalSelectedBkgd
private
Spectrum* ana::CutOptimization::fNominalSelectedSignal
private
Spectrum* ana::CutOptimization::fNominalSignal
private
SystematicDef* ana::CutOptimization::fNominalSystDef
private

Definition at line 129 of file CutOptimization.h.

Referenced by CutOptimization(), DefineNominal(), Go(), operator=(), and SetSpillCuts().

const HistAxis* ana::CutOptimization::fOptiHistAxis = NULL
private

Definition at line 153 of file CutOptimization.h.

Referenced by CutOptimization(), DefineMultiverseSystematic(), and operator=().

std::map< SystematicDef*, UpDownPair<Spectrum> > ana::CutOptimization::fSelected
private
std::map< SystematicDef*, UpDownPair<Spectrum> > ana::CutOptimization::fSelectedBkgd
private
std::map< SystematicDef*, UpDownPair<Spectrum> > ana::CutOptimization::fSelectedSignal
private
const Cut* ana::CutOptimization::fSelectionCut = NULL
private
std::map< SystematicDef*, UpDownPair<Spectrum> > ana::CutOptimization::fSignal
private
const Cut* ana::CutOptimization::fSignalCut = NULL
private
std::vector<SystematicDef*> ana::CutOptimization::fSystDefs
private

The documentation for this class was generated from the following files: