4 #include "TObjString.h" 20 "TrivialCrossSectionAnalysis object not yet ready to setup systematic spectra." );
22 fXSecSpectra.insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
61 std::vector<SystShifts> genie_shifts)
64 "CrossSectionAnalysis object not yet ready to setup genie multiverses." );
81 std::vector<Var> ppfx_weights,
82 std::vector<NuTruthVar> ppfx_weights_nt)
85 "CrossSectionAnalysis object not yet ready to setup PPFX multiverses.");
86 std::vector<Spectrum *> mv_flux;
87 for(
int i = 0;
i < (
int) ppfx_weights.size();
i++) {
111 for(
int i = 0;
i < (
int) ppfx_weights.size();
i++) {
123 num->Divide(num, denom, 1, 1,
"B");
136 std::vector<Cut> bkgd_cuts,
138 bool is_differential)
155 std::pair<TH1 *, TH1 *>
160 std::vector<TH1 *> shifts;
165 std::pair<TH1 *, TH1 *> flux_error =
FluxError();
166 std::pair<TH1 *, TH1 *> genie_error =
GenieError();
167 TH1 * flux_up = flux_error.first;
168 TH1 * flux_down = flux_error.second;
169 TH1 * genie_up = genie_error.first;
170 TH1 * genie_down = genie_error.second;
173 TH1 * errup = (TH1*) nominal->Clone();
174 TH1 * errdown = (TH1*) nominal->Clone();
175 double up, down,
delta;
177 for(
int i = 1;
i <= nominal->GetNbinsX();
i++) {
178 for(
int j = 1;
j <= nominal->GetNbinsY();
j++) {
179 for(
int k = 1; k <= nominal->GetNbinsZ(); k++) {
181 double totaldown = 0;
183 for(
int ishift = 0; ishift < (
int) shifts.size(); ishift++) {
185 shifts[ishift]->GetBinContent(
i,
j, k) - nominal->GetBinContent(
i,
j, k);
195 totaldown += down * down;
198 up = genie_up->GetBinContent(
i,
j, k) - nominal->GetBinContent(
i,
j, k);
199 down = genie_down->GetBinContent(
i,
j, k) - nominal->GetBinContent(
i,
j, k);
201 totaldown += down * down;
203 up = flux_up->GetBinContent(
i,
j, k) - nominal->GetBinContent(
i,
j, k);
204 down = flux_down->GetBinContent(
i,
j, k) - nominal->GetBinContent(
i,
j, k);
206 totaldown += down * down;
207 errup->SetBinContent(
i,
j, k, totalup);
208 errdown->SetBinContent(
i,
j, k, totaldown);
212 return std::pair<TH1 *, TH1 *> (errup, errdown);
221 for(
int i = 1;
i <= rel->GetNbinsX();
i++) {
222 for(
int j = 1;
j <= rel->GetNbinsY();
j++) {
223 for(
int k = 1; k <= rel->GetNbinsZ(); k++) {
225 shift =
fXSec[syst]->GetBinContent(
i,
j, k);
226 delta =
fabs(shift - cv);
227 rel->SetBinContent(
i,
j, k, delta / cv);
235 std::pair<TH1 *, TH1 *>
244 return std::pair<TH1 *, TH1 *>
250 std::pair<TH1 *, TH1 *>
260 return std::pair<TH1 *, TH1 *>
278 num->Divide(num, denom, 1, 1,
"B");
283 std::unique_ptr<TrivialCrossSectionAnalysis>
286 dir = dir->GetDirectory(name.c_str());
289 TObjString * ptag = (TObjString*) dir->Get(
"type");
292 const TString
tag = ptag->GetString();
293 assert(tag ==
"TrivialCrossSectionAnalysis" &&
"Type does not match TrivialCrossSectionAnalysis");
297 std::map<XSecSysts::Syst_t, CrossSectionSpectra *> xsec_spectra;
298 std::vector<XSecSysts::Syst_t>
systs 303 insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
308 TVector3 * fidmax = (TVector3*) dir->Get(
"fFidMax");
309 TVector3 * fidmin = (TVector3*) dir->Get(
"fFidMin");
316 bool isdifferential =
false;
317 if(dir->GetDirectory(
"fIsDifferential"))
318 isdifferential =
true;
319 if(dir->GetDirectory(
"fGenieMVEffNum")) {
323 if(dir->GetDirectory(
"fGenieMVEffDenom")) {
327 if(dir->GetDirectory(
"fPPFXMVFlux")) {
331 if(dir->GetDirectory(
"fPPFXMVEffNum")) {
335 if(dir->GetDirectory(
"fPPFXMVEffDenom")) {
340 std::unique_ptr<TrivialCrossSectionAnalysis>
ret 352 double POT = ret->fData->POT();
354 for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
xsec : xsec_spectra)
358 genie_mv_eff_num = NULL;
359 genie_mv_eff_denom = NULL;
361 ppfx_mv_eff_num = NULL;
362 ppfx_mv_eff_denom = NULL;
376 TDirectory *
tmp = gDirectory;
378 dir = dir->mkdir(name.c_str());
381 TObjString(
"TrivialCrossSectionAnalysis").Write(
"type");
383 TObjString(
"true").Write(
"fIsDifferential");
384 for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
xsec :
fXSecSpectra) {
392 ndims.Write(
"ndims");
426 [
this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
xsec)
429 .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
435 [
this](std::pair<XSecSysts::Syst_t, const Var *>
weight)
438 .insert(std::pair<XSecSysts::Syst_t, const Var *>
444 [
this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
447 .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
449 Copy(shift.second)));
467 this->
fPPFXXSec.push_back((TH1*) ppfx->Clone());
494 std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
495 [
this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
xsec)
498 .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
504 std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
505 [
this](std::pair<XSecSysts::Syst_t, const Var *>
weight)
508 .insert(std::pair<XSecSysts::Syst_t, const Var *>
514 std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
515 [
this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
518 .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
527 this->
fSelectionCut = rhs.fSelectionCut; rhs.fSelectionCut = NULL;
528 this->
fRecoAxis = rhs.fRecoAxis; rhs.fRecoAxis = NULL;
529 this->
fSignalCut = rhs.fSignalCut; rhs.fSignalCut = NULL;
530 this->
fTruthAxis = rhs.fTruthAxis; rhs.fTruthAxis = NULL;
533 this->
fNDims = rhs.fNDims;
535 if(rhs.fPPFXXSec.size() != 0) {
536 std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
543 if(rhs.fGenieXSec.size() != 0) {
544 std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
552 this->
fGenieMVEffNum = rhs.fGenieMVEffNum; rhs.fGenieMVEffNum = NULL;
554 this->
fPPFXMVFlux = rhs.fPPFXMVFlux; rhs.fPPFXMVFlux = NULL;
555 this->
fPPFXMVEffNum = rhs.fPPFXMVEffNum; rhs.fPPFXMVEffNum = NULL;
556 this->
fPPFXMVEffDenom = rhs.fPPFXMVEffDenom; rhs.fPPFXMVEffDenom = NULL;
557 this->
fData = rhs.fData; rhs.fData = NULL;
565 std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
566 [
this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
xsec)
569 .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
575 std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
576 [
this](std::pair<XSecSysts::Syst_t, const Var *>
weight)
579 .insert(std::pair<XSecSysts::Syst_t, const Var *>
585 std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
586 [
this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
589 .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
598 this->
fSelectionCut = rhs.fSelectionCut; rhs.fSelectionCut = NULL;
599 this->
fRecoAxis = rhs.fRecoAxis; rhs.fRecoAxis = NULL;
600 this->
fSignalCut = rhs.fSignalCut; rhs.fSignalCut = NULL;
601 this->
fTruthAxis = rhs.fTruthAxis; rhs.fTruthAxis = NULL;
604 this->
fNDims = rhs.fNDims;
606 if(rhs.fPPFXXSec.size() != 0) {
607 std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
614 if(rhs.fGenieXSec.size() != 0) {
615 std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
623 this->
fGenieMVEffNum = rhs.fGenieMVEffNum; rhs.fGenieMVEffNum = NULL;
625 this->
fPPFXMVFlux = rhs.fPPFXMVFlux; rhs.fPPFXMVFlux = NULL;
626 this->
fPPFXMVEffNum = rhs.fPPFXMVEffNum; rhs.fPPFXMVEffNum = NULL;
627 this->
fPPFXMVEffDenom = rhs.fPPFXMVEffDenom; rhs.fPPFXMVEffDenom = NULL;
628 this->
fData = rhs.fData; rhs.fData = NULL;
636 [
this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
xsec)
639 .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
645 [
this](std::pair<XSecSysts::Syst_t, const Var *>
weight)
648 .insert(std::pair<XSecSysts::Syst_t, const Var *>
654 [
this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
657 .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
659 Copy(shift.second)));
677 this->
fPPFXXSec.push_back((TH1*) ppfx->Clone());
698 std::map<XSecSysts::Syst_t, CrossSectionSpectra*> _fXSecSpectra,
710 std::for_each(_fXSecSpectra.begin(), _fXSecSpectra.end(),
711 [
this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
xsec)
714 .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
720 this->
fNDims = (*ndims)[0];
728 for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
xsec : _fXSecSpectra)
729 if(
xsec.second != NULL)
731 genie_mv_eff_num = NULL;
732 genie_mv_eff_denom = NULL;
734 ppfx_mv_eff_num = NULL;
735 ppfx_mv_eff_denom = NULL;
741 for(
int ipxsec = 0; ipxsec < (
int)
fPPFXXSec.size(); ipxsec++) {
744 for(
int igxsec = 0; igxsec < (
int)
fGenieXSec.size(); igxsec++) {
static TString LabelFromSystematic(Syst_t)
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
void SaveTo(TDirectory *dir, const std::string &name) const
THE MAIN GENIE PROJECT NAMESPACE
fvar< T > fabs(const fvar< T > &x)
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
TrivialCrossSectionAnalysis & operator=(const TrivialCrossSectionAnalysis &)
Simple record of shifts applied to systematic parameters.
Multiverse * fGenieMVEffDenom
Multiverse * fGenieMVEffNum
const Cut * fSelectionCut
TH1 * Flux(XSecSysts::Syst_t) override
ICrossSectionAnalysis * Clone() const override
TH1 * Efficiency(XSecSysts::Syst_t) override
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
_Cut< caf::SRNeutrinoProxy > NuTruthCut
Cut designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Representation of a spectrum in any variable, with associated POT.
const TH1 * UnfoldedSignal(XSecSysts::Syst_t)
TH1 * RelativeUncertainty(XSecSysts::Syst_t) override
Function that calculates relative uncertainty on the cross section result.
const XML_Char const XML_Char * data
std::pair< TH1 *, TH1 * > GenieError()
Function that calculates the up (first) and down (second) error on the genie universes.
static std::unique_ptr< CrossSectionSpectra > LoadFrom(TDirectory *, const std::string &)
void SetWeight(XSecSysts::Syst_t, const Var *)
const HistAxis * fRecoAxis
std::pair< TH1 *, TH1 * > TotalErrors() override
TrivialCrossSectionAnalysis()
default constructor. Do nothing
std::vector< Cut > fBkgdCuts
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
std::map< XSecSysts::Syst_t, CrossSectionSpectra * > fXSecSpectra
const std::string cv[Ncv]
void SaveTo(TDirectory *, const std::string &name) const override
Multiverse * fPPFXMVEffDenom
static std::vector< Syst_t > GetAllSystematics()
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
std::pair< TH1 *, TH1 * > FluxError() override
Function that calculates the up (first) and down (second) error on the flux estimation.
void SetSystShifts(XSecSysts::Syst_t, const SystShifts *)
static std::unique_ptr< TrivialCrossSectionAnalysis > LoadFrom(TDirectory *, const std::string &name)
std::map< XSecSysts::Syst_t, const NuTruthVar * > fNuTruthWeights
HistAxis HistAxisFromNuTruthHistAxis(NuTruthHistAxis ntha, double _default)
std::vector< TH1 * > fGenieXSec
void SaveTo(TDirectory *dir, const std::string &name) const
const Binning * fENuBinning
void InitPPFXUniverses(SpectrumLoaderBase *, std::vector< Var >, std::vector< NuTruthVar >)
std::vector< float > Spectrum
void SetNuTruthWeight(XSecSysts::Syst_t, const NuTruthVar *)
Multiverse * fPPFXMVEffNum
static std::unique_ptr< Multiverse > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< double > POT
Base class for the various types of spectrum loader.
std::vector< TH1 * > fPPFXXSec
unsigned int NDimensions() const
std::map< XSecSysts::Syst_t, const Var * > fWeights
TH1 * Purity(XSecSysts::Syst_t) override
Function that returns the purity of the specified sample.
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
void InitGenieUniverses(SpectrumLoaderBase *, std::vector< SystShifts >)
std::map< XSecSysts::Syst_t, const SystShifts * > fSystShifts
bool fHasGenieMultiverses
assert(nhit_max >=nhit_nbins)
const NuTruthHistAxis * fTruthAxis
std::map< XSecSysts::Syst_t, TH1 * > fXSec
TH1 * ToHist(const Spectrum &, double=-1)
const NuTruthCut * fSignalCut
bool CanFillSpectra() const
Count number of targets within the main part of the ND (vSuperBlockND)
void InitSystematicSpectra(XSecSysts::Syst_t, SpectrumLoaderBase *)
UnfoldMethod_t fUnfoldMethod
Cut CutFromNuTruthCut(const NuTruthCut &stc)
~TrivialCrossSectionAnalysis()
void Divide(Multiverse &, bool=false)