diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
index 88d2530..d2fc662 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
@@ -1,475 +1,379 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "Smear_SVDUnfold_Propagation_Osc.h"
#include "SmearceptanceUtils.h"
#include "OscWeightEngine.h"
//********************************************************************
Smear_SVDUnfold_Propagation_Osc::Smear_SVDUnfold_Propagation_Osc(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip =
"Simple measurement class for doing fake data oscillation studies.\n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("Osc Studies");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("XXX");
fSettings.SetYTitle("Number of events");
fSettings.SetEnuRange(0.0, 50);
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.DefineAllowedTargets("*");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") / (fNEvents + 0.);
// Plot Setup -------------------------------------------------------
// Check that we have the relevant histograms specified.
if (!samplekey.Has("NDDataHist") || !samplekey.Has("FDDataHist") ||
!samplekey.Has("NDToSpectrumSmearingMatrix")) {
THROW(
"Expected to attributes named: NDDataHist, FDDataHist, "
"NDToSpectrumSmearingMatrix on the Smear_SVDUnfold_Propagation_Osc "
"sample "
"tag.");
}
NDDataHist = NULL;
FDDataHist = NULL;
NDToSpectrumSmearingMatrix = NULL;
fMCHist = NULL;
fDataHist = NULL;
std::vector NDHistDescriptor =
GeneralUtils::ParseToStr(samplekey.GetS("NDDataHist"), ",");
if (NDHistDescriptor.size() == 2) {
NDDataHist = PlotUtils::GetTH1DFromRootFile(NDHistDescriptor[0],
NDHistDescriptor[1]);
ND_True_Spectrum_Hist =
PlotUtils::GetTH1DFromRootFile(NDHistDescriptor[0], "ELep_rate");
}
if (!NDDataHist) {
THROW("Attempted to load NDDataHist from the descriptor: \""
<< samplekey.GetS("NDDataHist")
<< "\", but failed. Does it look correct? \",\"");
}
std::vector FDHistDescriptor =
GeneralUtils::ParseToStr(samplekey.GetS("FDDataHist"), ",");
if (FDHistDescriptor.size() == 2) {
FDDataHist = PlotUtils::GetTH1DFromRootFile(FDHistDescriptor[0],
FDHistDescriptor[1]);
FD_True_Spectrum_Hist =
PlotUtils::GetTH1DFromRootFile(FDHistDescriptor[0], "ELep_rate");
}
if (!FDDataHist) {
THROW("Attempted to load FDDataHist from the descriptor: \""
<< samplekey.GetS("FDDataHist")
<< "\", but failed. Does it look correct? \",\"");
}
std::vector NDToEvSmearingDescriptor = GeneralUtils::ParseToStr(
samplekey.GetS("NDToSpectrumSmearingMatrix"), ",");
if (NDToEvSmearingDescriptor.size() == 2) {
NDToSpectrumSmearingMatrix = PlotUtils::GetTH2DFromRootFile(
NDToEvSmearingDescriptor[0], NDToEvSmearingDescriptor[1]);
}
if (!NDToSpectrumSmearingMatrix) {
THROW("Attempted to load NDToSpectrumSmearingMatrix from the descriptor: \""
<< samplekey.GetS("NDToSpectrumSmearingMatrix")
<< "\", but failed. Does it look correct? \",\"");
}
TMatrixD NDToSpectrumResponseMatrix_l = SmearceptanceUtils::GetMatrix(
SmearceptanceUtils::SVDGetInverse(NDToSpectrumSmearingMatrix));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
TMatrixD SpectrumToFDResponseMatrix_l =
SmearceptanceUtils::GetMatrix(NDToSpectrumSmearingMatrix);
SpectrumToFDResponseMatrix.ResizeTo(SpectrumToFDResponseMatrix_l);
SpectrumToFDResponseMatrix = SpectrumToFDResponseMatrix_l;
fDataHist = static_cast(FDDataHist->Clone());
-
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
+
+ fMCHist = static_cast(FDDataHist->Clone());
+ fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
+ (fSettings.GetFullTitles()).c_str());
+
SetCovarFromDiagonal();
TruncateUpTo = 0;
if (samplekey.Has("TruncateUpTo")) {
TruncateUpTo = samplekey.GetI("TruncateUpTo");
QLOG(SAM, "Allowed to truncate unfolding matrix by up to "
<< TruncateUpTo
<< " singular values to limit negative ENu spectra.");
}
if (samplekey.Has("FitRegion_Min")) {
FitRegion_Min = samplekey.GetD("FitRegion_Min");
} else {
FitRegion_Min = 0xdeadbeef;
}
if (samplekey.Has("FitRegion_Max")) {
FitRegion_Max = samplekey.GetD("FitRegion_Max");
} else {
FitRegion_Max = 0xdeadbeef;
}
+
+ //----------------------------
+ // Mask data hist if needed
+ fDataHist->SetBinContent(0, 0);
+ fDataHist->SetBinError(0, 0);
+ fDataHist->SetBinContent(fDataHist->GetXaxis()->GetNbins() + 1, 0);
+ fDataHist->SetBinError(fDataHist->GetXaxis()->GetNbins() + 1, 0);
+
+ for (Int_t bi_it = 1; bi_it < fDataHist->GetXaxis()->GetNbins() + 1;
+ ++bi_it) {
+ if ((FitRegion_Min != 0xdeadbeef) &&
+ (fDataHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
+ fDataHist->SetBinContent(bi_it, 0);
+ fDataHist->SetBinError(bi_it, 0);
+ }
+ if ((FitRegion_Max != 0xdeadbeef) &&
+ (fDataHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
+ fDataHist->SetBinContent(bi_it, 0);
+ fDataHist->SetBinError(bi_it, 0);
+ }
+ }
+
TruncateStart = 0;
if (Config::Get().GetConfigNode("smear.SVD.truncation")) {
TruncateStart = Config::Get().ConfI("smear.SVD.truncation");
TMatrixD NDToSpectrumResponseMatrix_l =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
NDToSpectrumSmearingMatrix, TruncateStart));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
}
if (TruncateStart >= TruncateUpTo) {
TruncateUpTo = TruncateStart + 1;
}
+ NDFDRatio = 1;
+ if (Config::Get().GetConfigNode("Osc.NDFDRatio")) {
+ NDFDRatio = Config::Get().ConfD("Osc.NDFDRatio");
+ }
+
+ UnfoldToNDETrueSpectrum();
+
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void Smear_SVDUnfold_Propagation_Osc::FillEventVariables(FitEvent *event){};
bool Smear_SVDUnfold_Propagation_Osc::isSignal(FitEvent *event) {
return false;
}
void Smear_SVDUnfold_Propagation_Osc::UnfoldToNDETrueSpectrum(void) {
ND_Unfolded_Spectrum_Hist = NDToSpectrumSmearingMatrix->ProjectionY();
ND_Unfolded_Spectrum_Hist->Clear();
-
- TRandom3 rnjesus;
+ ND_Unfolded_Spectrum_Hist->SetName("ND_Unfolded_Spectrum_Hist");
bool HasNegValue = false;
int truncations = TruncateStart;
do {
if (truncations >= TruncateUpTo) {
THROW("Unfolded enu spectrum had negative values even after "
<< truncations << " SVD singular value truncations.");
}
// Unfold ND ERec -> Enu spectrum
// ------------------------------
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
NDDataHist, ND_Unfolded_Spectrum_Hist, NDToSpectrumResponseMatrix, 1000,
false);
HasNegValue = false;
for (Int_t bi_it = 1;
bi_it < ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
++bi_it) {
if ((FitRegion_Min != 0xdeadbeef) &&
(ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(bi_it) <=
FitRegion_Min)) {
continue;
}
if ((FitRegion_Max != 0xdeadbeef) &&
(ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(bi_it) >
FitRegion_Max)) {
continue;
}
if (ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) < 0) {
HasNegValue = true;
break;
}
}
if (HasNegValue) {
TMatrixD NDToSpectrumResponseMatrix_l =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
NDToSpectrumSmearingMatrix, truncations));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
}
truncations++;
} while (HasNegValue);
-}
-void Smear_SVDUnfold_Propagation_Osc::ConvertEventRates(void) {
+ NDFD_Corrected_Spectrum_Hist =
+ static_cast(ND_Unfolded_Spectrum_Hist->Clone());
+ NDFD_Corrected_Spectrum_Hist->Clear();
+ NDFD_Corrected_Spectrum_Hist->SetName("NDFD_Corrected_Spectrum_Hist");
+ FD_Propagated_Spectrum_Hist =
+ static_cast(ND_Unfolded_Spectrum_Hist->Clone());
+ FD_Propagated_Spectrum_Hist->Clear();
+ FD_Propagated_Spectrum_Hist->SetName("FD_Propagated_Spectrum_Hist");
+
// Apply FD/ND weights
// ------------------------------
- OscHist->Write("ND_FD_Corrected_ENuTrue", TObject::kOverwrite);
+ for (Int_t bi_it = 1;
+ bi_it < ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins() + 1; ++bi_it) {
+ NDFD_Corrected_Spectrum_Hist->SetBinContent(
+ bi_it, ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) * NDFDRatio);
+ NDFD_Corrected_Spectrum_Hist->SetBinError(
+ bi_it, ND_Unfolded_Spectrum_Hist->GetBinError(bi_it) * NDFDRatio);
+ }
+}
+void Smear_SVDUnfold_Propagation_Osc::ConvertEventRates(void) {
// Apply Oscillations
// ------------------------------
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
- for (Int_t bi_it = 1; bi_it < OscHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
- double oscWeight =
- oscWE->CalcWeight(OscHist->GetXaxis()->GetBinCenter(bi_it), 14);
- OscHist->SetBinContent(bi_it, OscHist->GetBinContent(bi_it) * oscWeight);
+ for (Int_t bi_it = 1;
+ bi_it < NDFD_Corrected_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
+ ++bi_it) {
+ double oscWeight = oscWE->CalcWeight(
+ NDFD_Corrected_Spectrum_Hist->GetXaxis()->GetBinCenter(bi_it), 14);
+ FD_Propagated_Spectrum_Hist->SetBinContent(
+ bi_it, NDFD_Corrected_Spectrum_Hist->GetBinContent(bi_it) * oscWeight);
+ FD_Propagated_Spectrum_Hist->SetBinError(
+ bi_it, NDFD_Corrected_Spectrum_Hist->GetBinError(bi_it) * oscWeight);
}
- OscHist->Write("ENuTrue_FD", TObject::kOverwrite);
-
- TGraph POsc;
-
- POsc.Set(1E4 - 1);
-
- double min = OscHist->GetXaxis()->GetBinLowEdge(1);
- double step =
- (OscHist->GetXaxis()->GetBinUpEdge(OscHist->GetXaxis()->GetNbins()) -
- OscHist->GetXaxis()->GetBinLowEdge(1)) /
- double(1E4);
-
- for (size_t i = 1; i < 1E4; ++i) {
- double enu = min + i * step;
- double ow = oscWE->CalcWeight(enu, 14);
- if (ow != ow) {
- std::cout << "Bad osc weight for ENu: " << enu << std::endl;
- }
- POsc.SetPoint(i - 1, enu, ow);
- }
-
- POsc.Write("POsc", TObject::kOverwrite);
// Forward fold Spectrum -> ERec FD
// ------------------------------
- if (fMCHist) {
- fMCHist->SetDirectory(NULL);
- delete fMCHist;
- }
-
- fMCHist = static_cast(FDDataHist->Clone());
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
- OscHist, fMCHist, SpectrumToFDResponseMatrix, 1000, true);
+ FD_Propagated_Spectrum_Hist, fMCHist, SpectrumToFDResponseMatrix, 1000,
+ true);
fMCHist->SetBinContent(0, 0);
fMCHist->SetBinError(0, 0);
- fDataHist->SetBinContent(0, 0);
- fDataHist->SetBinError(0, 0);
fMCHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
fMCHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
- fDataHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
- fDataHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
-
for (Int_t bi_it = 1; bi_it < fMCHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
if ((FitRegion_Min != 0xdeadbeef) &&
(fMCHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
fMCHist->SetBinContent(bi_it, 0);
fMCHist->SetBinError(bi_it, 0);
- fDataHist->SetBinContent(bi_it, 0);
- fDataHist->SetBinError(bi_it, 0);
}
if ((FitRegion_Max != 0xdeadbeef) &&
(fMCHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
fMCHist->SetBinContent(bi_it, 0);
fMCHist->SetBinError(bi_it, 0);
- fDataHist->SetBinContent(bi_it, 0);
- fDataHist->SetBinError(bi_it, 0);
}
}
}
-void Smear_SVDUnfold_Propagation_Osc::Write(void) {
+void Smear_SVDUnfold_Propagation_Osc::Write(std::string drawOpt) {
NDToSpectrumSmearingMatrix->Write("SmearingMatrix_ND", TObject::kOverwrite);
- TH1D *OscHist = NDToSpectrumSmearingMatrix->ProjectionY();
- OscHist->Clear();
-
- TRandom3 rnjesus;
-
- NDDataHist->Write("ERec_ND", TObject::kOverwrite);
-
- bool HasNegValue = false;
- int truncations = TruncateStart;
- do {
- if (truncations >= TruncateUpTo) {
- THROW("Unfolded enu spectrum had negative values even after "
- << truncations << " SVD singular value truncations.");
- }
-
- // Unfold ND ERec -> Enu spectrum
- // ------------------------------
- SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
- NDDataHist, OscHist, NDToSpectrumResponseMatrix, 1000, false);
-
- std::stringstream ss("");
- ss << "Unfolded_ENuTrue_ND_trunc_" << truncations;
- OscHist->Write(ss.str().c_str(), TObject::kOverwrite);
- ss.str("");
- ss << "Response_ND_trunc_" << truncations;
- SmearceptanceUtils::SVDGetInverse(NDToSpectrumSmearingMatrix, truncations)
- ->Write(ss.str().c_str(), TObject::kOverwrite);
+ NDDataHist->Write("Obs_ND", TObject::kOverwrite);
- HasNegValue = false;
-
- for (Int_t bi_it = 1; bi_it < OscHist->GetXaxis()->GetNbins() + 1;
- ++bi_it) {
- if ((FitRegion_Min != 0xdeadbeef) &&
- (OscHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
- continue;
- }
- if ((FitRegion_Max != 0xdeadbeef) &&
- (OscHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
- continue;
- }
+ ND_Unfolded_Spectrum_Hist->Write(ND_Unfolded_Spectrum_Hist->GetName(),
+ TObject::kOverwrite);
+ NDFD_Corrected_Spectrum_Hist->Write(NDFD_Corrected_Spectrum_Hist->GetName(),
+ TObject::kOverwrite);
+ FD_Propagated_Spectrum_Hist->Write(FD_Propagated_Spectrum_Hist->GetName(),
+ TObject::kOverwrite);
- if (OscHist->GetBinContent(bi_it) < 0) {
- HasNegValue = true;
- break;
- }
- }
-
- if (HasNegValue) {
- TMatrixD NDToSpectrumResponseMatrix_l =
- SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
- NDToSpectrumSmearingMatrix, truncations));
-
- NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
- NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
- }
-
- truncations++;
- } while (HasNegValue);
-
- // Apply FD/ND weights
- // ------------------------------
- OscHist->Write("ND_FD_Corrected_ENuTrue", TObject::kOverwrite);
+ fMCHist->Write("Pred_FD", TObject::kOverwrite);
+ fDataHist->Write("Obs_FD", TObject::kOverwrite);
// Apply Oscillations
// ------------------------------
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
- for (Int_t bi_it = 1; bi_it < OscHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
- double oscWeight =
- oscWE->CalcWeight(OscHist->GetXaxis()->GetBinCenter(bi_it), 14);
- OscHist->SetBinContent(bi_it, OscHist->GetBinContent(bi_it) * oscWeight);
- }
- OscHist->Write("ENuTrue_FD", TObject::kOverwrite);
-
TGraph POsc;
POsc.Set(1E4 - 1);
- double min = OscHist->GetXaxis()->GetBinLowEdge(1);
- double step =
- (OscHist->GetXaxis()->GetBinUpEdge(OscHist->GetXaxis()->GetNbins()) -
- OscHist->GetXaxis()->GetBinLowEdge(1)) /
- double(1E4);
+ double min = ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1);
+ double step = (ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(
+ ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins()) -
+ ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1)) /
+ double(1E4);
for (size_t i = 1; i < 1E4; ++i) {
double enu = min + i * step;
double ow = oscWE->CalcWeight(enu, 14);
if (ow != ow) {
std::cout << "Bad osc weight for ENu: " << enu << std::endl;
}
POsc.SetPoint(i - 1, enu, ow);
}
POsc.Write("POsc", TObject::kOverwrite);
-
- // Forward fold Spectrum -> ERec FD
- // ------------------------------
- if (fMCHist) {
- fMCHist->SetDirectory(NULL);
- delete fMCHist;
- }
-
- fMCHist = static_cast(FDDataHist->Clone());
-
- SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
- OscHist, fMCHist, SpectrumToFDResponseMatrix, 1000, true);
-
- NDToSpectrumSmearingMatrix->Write("SmearingMatrix_FD", TObject::kOverwrite);
-
- fMCHist->Write("ERec_FD_pred", TObject::kOverwrite);
- fDataHist->Write("ERec_FD_obs", TObject::kOverwrite);
-
- fMCHist->SetBinContent(0, 0);
- fMCHist->SetBinError(0, 0);
- fDataHist->SetBinContent(0, 0);
- fDataHist->SetBinError(0, 0);
-
- fMCHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
- fMCHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
- fDataHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
- fDataHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
-
- for (Int_t bi_it = 1; bi_it < fMCHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
- if ((FitRegion_Min != 0xdeadbeef) &&
- (fMCHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
- fMCHist->SetBinContent(bi_it, 0);
- fMCHist->SetBinError(bi_it, 0);
- fDataHist->SetBinContent(bi_it, 0);
- fDataHist->SetBinError(bi_it, 0);
- }
- if ((FitRegion_Max != 0xdeadbeef) &&
- (fMCHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
- fMCHist->SetBinContent(bi_it, 0);
- fMCHist->SetBinError(bi_it, 0);
- fDataHist->SetBinContent(bi_it, 0);
- fDataHist->SetBinError(bi_it, 0);
- }
- }
-
- fMCHist->Write("ERec_FD_pred_masked", TObject::kOverwrite);
- fDataHist->Write("ERec_FD_obs_masked", TObject::kOverwrite);
-
- if (ND_True_Spectrum_Hist) {
- ND_True_Spectrum_Hist->Write("ETrue_ND_Spectrum", TObject::kOverwrite);
- }
- if (FD_True_Spectrum_Hist) {
- FD_True_Spectrum_Hist->Write("ETrue_FD_Spectrum", TObject::kOverwrite);
- }
}
diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
index db8f646..5fd754b 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
@@ -1,57 +1,61 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef Smear_SVDUnfold_Propagation_Osc_H_SEEN
#define Smear_SVDUnfold_Propagation_Osc_H_SEEN
#include "Measurement1D.h"
class Smear_SVDUnfold_Propagation_Osc : public Measurement1D {
public:
Smear_SVDUnfold_Propagation_Osc(nuiskey samplekey);
virtual ~Smear_SVDUnfold_Propagation_Osc(){};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
+ void UnfoldToNDETrueSpectrum(void);
void ConvertEventRates(void);
+ void Write(std::string drawOpt);
TH1D *NDDataHist;
TH1D *FDDataHist;
TH1D *ND_Unfolded_Spectrum_Hist;
TH1D *NDFD_Corrected_Spectrum_Hist;
TH1D *FD_Propagated_Spectrum_Hist;
TH1D *ND_True_Spectrum_Hist;
TH1D *FD_True_Spectrum_Hist;
TH2D *NDToSpectrumSmearingMatrix;
TMatrixD NDToSpectrumResponseMatrix;
TMatrixD SpectrumToFDResponseMatrix;
Int_t TruncateStart;
Int_t TruncateUpTo;
double FitRegion_Min;
double FitRegion_Max;
+
+ double NDFDRatio;
};
#endif
diff --git a/src/Reweight/OscWeightEngine.cxx b/src/Reweight/OscWeightEngine.cxx
index a36fd9a..7fde603 100644
--- a/src/Reweight/OscWeightEngine.cxx
+++ b/src/Reweight/OscWeightEngine.cxx
@@ -1,291 +1,291 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
-#define DEBUG_OSC_WE
+//#define DEBUG_OSC_WE
#include "OscWeightEngine.h"
#include
enum nuTypes {
kNuebarType = -1,
kNumubarType = -2,
kNutaubarType = -3,
kNueType = 1,
kNumuType = 2,
kNutauType = 3,
};
nuTypes GetNuType(int pdg) {
switch (pdg) {
case 16:
return kNutauType;
case 14:
return kNumuType;
case 12:
return kNueType;
case -16:
return kNutaubarType;
case -14:
return kNumubarType;
case -12:
return kNuebarType;
default: { THROW("Attempting to convert \"neutrino pdg\": " << pdg); }
}
}
OscWeightEngine::OscWeightEngine()
:
#ifdef __PROB3PP_ENABLED__
bp(),
#endif
theta12(0.825),
theta13(0.10),
theta23(1.0),
dm12(7.9e-5),
dm23(2.5e-3),
dcp(0.0),
LengthParam(0xdeadbeef),
TargetNuType(0) {
Config();
}
void OscWeightEngine::Config() {
std::vector OscParam = Config::QueryKeys("OscParam");
if (OscParam.size() < 1) {
ERROR(WRN,
"Oscillation parameters specified but no OscParam element "
"configuring the experimental characteristics found.\nExpect at "
"least . Pausing for "
"10...");
sleep(10);
return;
}
if (OscParam[0].Has("baseline_km")) {
LengthParamIsZenith = false;
LengthParam = OscParam[0].GetD("baseline_km");
constant_density = OscParam[0].Has("matter_density")
? OscParam[0].GetD("matter_density")
: 0xdeadbeef;
} else if (OscParam[0].Has("detection_zenith_deg")) {
LengthParamIsZenith = true;
static const double deg2rad = asin(1) / 90.0;
LengthParam = cos(OscParam[0].GetD("detection_zenith_deg") * deg2rad);
} else {
ERROR(WRN,
"It appeared that you wanted to set up an oscillation weight "
"branch, but it was not correctly configured. You need to specify "
"either: detection_zenith_deg or baseline_km attributes on the "
"OscParam element, and if baseline_km is specified, you can "
"optionally also set matter_density for oscillations through a "
"constant matter density. Pausing for 10...");
sleep(10);
return;
}
dm23 = OscParam[0].Has("dm23") ? OscParam[0].GetD("dm23") : dm23;
theta23 = OscParam[0].Has("sinsq_theta23") ? OscParam[0].GetD("sinsq_theta23")
: theta23;
theta13 = OscParam[0].Has("sinsq_theta13") ? OscParam[0].GetD("sinsq_theta13")
: theta13;
dm12 = OscParam[0].Has("dm12") ? OscParam[0].GetD("dm12") : dm12;
theta12 = OscParam[0].Has("sinsq_theta12") ? OscParam[0].GetD("sinsq_theta12")
: theta12;
dcp = OscParam[0].Has("dcp") ? OscParam[0].GetD("dcp") : dcp;
TargetNuType = OscParam[0].Has("TargetNuPDG")
? GetNuType(OscParam[0].GetI("TargetNuPDG"))
: 0;
QLOG(FIT, "Configured oscillation weighter:");
if (LengthParamIsZenith) {
QLOG(FIT,
"Earth density profile with detection cos(zenith) = " << LengthParam);
} else {
if (constant_density != 0xdeadbeef) {
QLOG(FIT,
"Constant density with experimental baseline = " << LengthParam);
} else {
QLOG(FIT,
"Vacuum oscillations with experimental baseline = " << LengthParam);
}
}
params[0] = dm23;
params[1] = theta23;
params[2] = theta13;
params[3] = dm12;
params[4] = theta12;
params[5] = dcp;
QLOG(FIT, "\tdm23 : " << params[0]);
QLOG(FIT, "\tsinsq_theta23: " << params[1]);
QLOG(FIT, "\tsinsq_theta13: " << params[2]);
QLOG(FIT, "\tdm12 : " << params[3]);
QLOG(FIT, "\tsinsq_theta12: " << params[4]);
QLOG(FIT, "\tdcp : " << params[5]);
}
void OscWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef DEBUG_OSC_WE
std::cout << "IncludeDial: " << name << " at " << startval << std::endl;
#endif
int dial = SystEnumFromString(name);
if (!dial) {
THROW("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
params[dial - 1] = startval;
}
void OscWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef DEBUG_OSC_WE
std::cout << "SetDial: " << (nuisenum % 1000) << " at " << val << std::endl;
#endif
fHasChanged = (params[(nuisenum % 1000) - 1] - val) >
std::numeric_limits::epsilon();
params[(nuisenum % 1000) - 1] = val;
}
void OscWeightEngine::SetDialValue(std::string name, double val) {
#ifdef DEBUG_OSC_WE
std::cout << "SetDial: " << name << " at " << val << std::endl;
#endif
int dial = SystEnumFromString(name);
if (!dial) {
THROW("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
fHasChanged =
(params[dial - 1] - val) > std::numeric_limits::epsilon();
params[dial - 1] = val;
}
bool OscWeightEngine::IsDialIncluded(std::string name) {
return SystEnumFromString(name);
}
bool OscWeightEngine::IsDialIncluded(int nuisenum) {
return ((nuisenum % 1000) > 0) && ((nuisenum % 1000) < 6);
}
double OscWeightEngine::GetDialValue(std::string name) {
int dial = SystEnumFromString(name);
if (!dial) {
THROW("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
return params[dial - 1];
}
double OscWeightEngine::GetDialValue(int nuisenum) {
if (!(nuisenum % 1000) || (nuisenum % 1000) > 6) {
THROW("OscWeightEngine passed dial enum: "
<< (nuisenum % 1000)
<< " that it does not understand, expected [1,6].");
}
return params[(nuisenum % 1000) - 1];
}
void OscWeightEngine::Reconfigure(bool silent) { fHasChanged = false; };
bool OscWeightEngine::NeedsEventReWeight() {
if (fHasChanged) {
return true;
}
return false;
}
double OscWeightEngine::CalcWeight(BaseFitEvt* evt) {
static bool Warned = false;
if (evt->probe_E == 0xdeadbeef) {
if (!Warned) {
ERROR(WRN,
"Oscillation weights asked for but using 'litemode' or "
"unsupported generator input. Pasuing for 10...");
sleep(10);
Warned = true;
}
return 1;
}
return CalcWeight(evt->probe_E * 1E-3, evt->probe_pdg);
}
double OscWeightEngine::CalcWeight(double ENu, int PDGNu) {
if (LengthParam == 0xdeadbeef) { // not configured.
return 1;
}
#ifdef __PROB3PP_ENABLED__
int NuType = GetNuType(PDGNu);
bp.SetMNS(params[theta12_idx], params[theta13_idx], params[theta23_idx],
params[dm12_idx], params[dm23_idx], params[dcp_idx], ENu,
true, NuType);
int pmt = 0;
double prob_weight = 1;
if (LengthParamIsZenith) { // Use earth density
bp.DefinePath(LengthParam, 0);
bp.propagate(NuType);
pmt = 0;
prob_weight = bp.GetProb(NuType, TargetNuType ? TargetNuType : NuType);
} else {
if (constant_density != 0xdeadbeef) {
bp.propagateLinear(NuType, LengthParam, constant_density);
pmt = 1;
prob_weight = bp.GetProb(NuType, TargetNuType ? TargetNuType : NuType);
} else {
pmt = 2;
prob_weight =
bp.GetVacuumProb(NuType, TargetNuType ? TargetNuType : NuType,
ENu * 1E-3, LengthParam);
}
}
#ifdef DEBUG_OSC_WE
if (prob_weight != prob_weight) {
THROW("Calculated bad prob weight: "
<< prob_weight << "(Osc Type: " << pmt << " -- " << NuType << " -> "
<< (TargetNuType ? TargetNuType : NuType) << ")");
}
#endif
return prob_weight;
#else
return 1;
#endif
}
int OscWeightEngine::SystEnumFromString(std::string const& name) {
if (name == "dm23") {
return 1;
} else if (name == "sinsq_theta23") {
return 2;
} else if (name == "sinsq_theta13") {
return 3;
} else if (name == "dm12") {
return 4;
} else if (name == "sinsq_theta12") {
return 5;
} else if (name == "dcp") {
return 6;
} else {
return 0;
}
}
diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx
index 9d1dc26..e5ad2a5 100755
--- a/src/Routines/MinimizerRoutines.cxx
+++ b/src/Routines/MinimizerRoutines.cxx
@@ -1,1510 +1,1511 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MinimizerRoutines.h"
+#include "Simple_MH_Sampler.h"
+
/*
Constructor/Destructor
*/
//************************
void MinimizerRoutines::Init() {
-//************************
+ //************************
fInputFile = "";
fInputRootFile = NULL;
fOutputFile = "";
fOutputRootFile = NULL;
- fCovar = NULL;
- fCovFree = NULL;
- fCorrel = NULL;
+ fCovar = NULL;
+ fCovFree = NULL;
+ fCorrel = NULL;
fCorFree = NULL;
- fDecomp = NULL;
+ fDecomp = NULL;
fDecFree = NULL;
fStrategy = "Migrad,FixAtLimBreak,Migrad";
fRoutines.clear();
fCardFile = "";
fFakeDataInput = "";
- fSampleFCN = NULL;
+ fSampleFCN = NULL;
- fMinimizer = NULL;
+ fMinimizer = NULL;
fMinimizerFCN = NULL;
- fCallFunctor = NULL;
-
- fAllowedRoutines = ("Migrad,Simplex,Combined,"
- "Brute,Fumili,ConjugateFR,"
- "ConjugatePR,BFGS,BFGS2,"
- "SteepDesc,GSLSimAn,FixAtLim,FixAtLimBreak,"
- "Chi2Scan1D,Chi2Scan2D,Contours,ErrorBands,"
- "DataToys");
+ fCallFunctor = NULL;
+
+ fAllowedRoutines =
+ ("Migrad,Simplex,Combined,"
+ "Brute,Fumili,ConjugateFR,"
+ "ConjugatePR,BFGS,BFGS2,"
+ "SteepDesc,GSLSimAn,FixAtLim,FixAtLimBreak,"
+ "Chi2Scan1D,Chi2Scan2D,Contours,ErrorBands,"
+ "DataToys,MCMC");
};
//*************************************
-MinimizerRoutines::~MinimizerRoutines() {
-//*************************************
+MinimizerRoutines::~MinimizerRoutines(){
+ //*************************************
};
/*
Input Functions
*/
//*************************************
MinimizerRoutines::MinimizerRoutines(int argc, char* argv[]) {
-//*************************************
+ //*************************************
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector xmlcmds;
std::vector configargs;
// Make easier to handle arguments.
std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
- ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
+ ERR(WRN) << "No output supplied so saving it to: " << fOutputFile
+ << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
nuiskey fCompKey;
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
- if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
+ if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
- if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
+ if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
- configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
+ configuration.LoadConfig(fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
- if (maxevents.compare("-1")){
+ if (maxevents.compare("-1")) {
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
// FitPar::log_verb = verbocount;
SETVERBOSITY(verbocount);
// ERR_VERB(errorcount);
// Minimizer Setup ========================================
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupMinimizerFromXML();
SetupCovariance();
SetupRWEngine();
SetupFCN();
return;
};
//*************************************
void MinimizerRoutines::SetupMinimizerFromXML() {
-//*************************************
+ //*************************************
LOG(FIT) << "Setting up nuismin" << std::endl;
// Setup Parameters ------------------------------------------
std::vector parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
- double parnom = key.GetD("nominal");
- double parlow = parnom - 1;
+ double parnom = key.GetD("nominal");
+ double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
// Override state if none given
- if (!key.Has("state")){
- key.SetS("state","FIX");
+ if (!key.Has("state")) {
+ key.SetS("state", "FIX");
}
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("low")) {
- parlow = key.GetD("low");
+ parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
- LOG(FIT) << "Read " << partype << " : "
- << parname << " = "
- << parnom << " : "
- << parlow << " < p < " << parhigh
- << " : " << parstate << std::endl;
+ LOG(FIT) << "Read " << partype << " : " << parname << " = " << parnom
+ << " : " << parlow << " < p < " << parhigh << " : " << parstate
+ << std::endl;
} else {
- LOG(FIT) << "Read " << partype << " : "
- << parname << " = "
- << parnom << " : "
- << parstate << std::endl;
+ LOG(FIT) << "Read " << partype << " : " << parname << " = " << parnom
+ << " : " << parstate << std::endl;
}
// Run Parameter Conversion if needed
if (parstate.find("ABS") != std::string::npos) {
- parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
- parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
- parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
- parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
+ parnom = FitBase::RWAbsToSigma(partype, parname, parnom);
+ parlow = FitBase::RWAbsToSigma(partype, parname, parlow);
+ parhigh = FitBase::RWAbsToSigma(partype, parname, parhigh);
+ parstep = FitBase::RWAbsToSigma(partype, parname, parstep);
} else if (parstate.find("FRAC") != std::string::npos) {
- parnom = FitBase::RWFracToSigma( partype, parname, parnom );
- parlow = FitBase::RWFracToSigma( partype, parname, parlow );
- parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
- parstep = FitBase::RWFracToSigma( partype, parname, parstep );
+ parnom = FitBase::RWFracToSigma(partype, parname, parnom);
+ parlow = FitBase::RWFracToSigma(partype, parname, parlow);
+ parhigh = FitBase::RWFracToSigma(partype, parname, parhigh);
+ parstep = FitBase::RWFracToSigma(partype, parname, parstep);
}
// Push into vectors
fParams.push_back(parname);
- fTypeVals[parname] = FitBase::ConvDialType(partype);;
+ fTypeVals[parname] = FitBase::ConvDialType(partype);
+ ;
fStartVals[parname] = parnom;
- fCurVals[parname] = parnom;
+ fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
- fStateVals[parname] = parstate;
+ fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
- fFixVals[parname] = fixstate;
+ fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
- fMinVals[parname] = parlow;
- fMaxVals[parname] = parhigh;
+ fMinVals[parname] = parlow;
+ fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
-
}
// Setup Samples ----------------------------------------------
- std::vector samplekeys = Config::QueryKeys("sample");
+ std::vector samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
- std::string sampletype =
- key.Has("type") ? key.GetS("type") : "DEFAULT";
+ std::string sampletype = key.Has("type") ? key.GetS("type") : "DEFAULT";
- double samplenorm =
- key.Has("norm") ? key.GetD("norm") : 1.0;
+ double samplenorm = key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
- LOG(FIT) << "Read sample info " << i << " : "
- << samplename << std::endl
- << "\t\t input -> " << samplefile << std::endl
+ LOG(FIT) << "Read sample info " << i << " : " << samplename << std::endl
+ << "\t\t input -> " << samplefile << std::endl
<< "\t\t state -> " << sampletype << std::endl
<< "\t\t norm -> " << samplenorm << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find(normname) != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fStartVals[normname] = samplenorm;
fCurVals[normname] = samplenorm;
fErrorVals[normname] = 0.0;
- fMinVals[normname] = 0.1;
- fMaxVals[normname] = 10.0;
+ fMinVals[normname] = 0.1;
+ fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
- fFixVals[normname] = state;
+ fFixVals[normname] = state;
fStartFixVals[normname] = state;
-
-
-
}
// Setup Fake Parameters -----------------------------
std::vector fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nom")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
- double parnom = key.GetD("nom");
+ double parnom = key.GetD("nom");
// Push into vectors
fFakeVals[parname] = parnom;
}
-
-
}
-
/*
Setup Functions
*/
//*************************************
void MinimizerRoutines::SetupRWEngine() {
-//*************************************
+ //*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
- FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
+ FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name));
}
UpdateRWEngine(fStartVals);
return;
}
//*************************************
void MinimizerRoutines::SetupFCN() {
-//*************************************
+ //*************************************
LOG(FIT) << "Making the jointFCN" << std::endl;
if (fSampleFCN) delete fSampleFCN;
// fSampleFCN = new JointFCN(fCardFile, fOutputRootFile);
fSampleFCN = new JointFCN(fOutputRootFile);
-
+
SetFakeData();
- fMinimizerFCN = new MinimizerFCN( fSampleFCN );
- fCallFunctor = new ROOT::Math::Functor( *fMinimizerFCN, fParams.size() );
+ fMinimizerFCN = new MinimizerFCN(fSampleFCN);
+ fCallFunctor = new ROOT::Math::Functor(*fMinimizerFCN, fParams.size());
- fSampleFCN->CreateIterationTree( "fit_iterations", FitBase::GetRW() );
+ fSampleFCN->CreateIterationTree("fit_iterations", FitBase::GetRW());
return;
}
-
//******************************************
void MinimizerRoutines::SetupFitter(std::string routine) {
-//******************************************
+ //******************************************
// Make the fitter
std::string fitclass = "";
- std::string fittype = "";
+ std::string fittype = "";
+ bool UseMCMC = false;
// Get correct types
- if (!routine.compare("Migrad")) {
- fitclass = "Minuit2"; fittype = "Migrad";
- } else if (!routine.compare("Simplex")) {
- fitclass = "Minuit2"; fittype = "Simplex";
- } else if (!routine.compare("Combined")) {
- fitclass = "Minuit2"; fittype = "Combined";
- } else if (!routine.compare("Brute")) {
- fitclass = "Minuit2"; fittype = "Scan";
- } else if (!routine.compare("Fumili")) {
- fitclass = "Minuit2"; fittype = "Fumili";
+ if (!routine.compare("Migrad")) {
+ fitclass = "Minuit2";
+ fittype = "Migrad";
+ } else if (!routine.compare("Simplex")) {
+ fitclass = "Minuit2";
+ fittype = "Simplex";
+ } else if (!routine.compare("Combined")) {
+ fitclass = "Minuit2";
+ fittype = "Combined";
+ } else if (!routine.compare("Brute")) {
+ fitclass = "Minuit2";
+ fittype = "Scan";
+ } else if (!routine.compare("Fumili")) {
+ fitclass = "Minuit2";
+ fittype = "Fumili";
} else if (!routine.compare("ConjugateFR")) {
- fitclass = "GSLMultiMin"; fittype = "ConjugateFR";
+ fitclass = "GSLMultiMin";
+ fittype = "ConjugateFR";
} else if (!routine.compare("ConjugatePR")) {
- fitclass = "GSLMultiMin"; fittype = "ConjugatePR";
- } else if (!routine.compare("BFGS")) {
- fitclass = "GSLMultiMin"; fittype = "BFGS";
- } else if (!routine.compare("BFGS2")) {
- fitclass = "GSLMultiMin"; fittype = "BFGS2";
- } else if (!routine.compare("SteepDesc")) {
- fitclass = "GSLMultiMin"; fittype = "SteepestDescent";
- // } else if (!routine.compare("GSLMulti")) { fitclass = "GSLMultiFit"; fittype = ""; // Doesn't work out of the box
- } else if (!routine.compare("GSLSimAn")) { fitclass = "GSLSimAn"; fittype = ""; }
+ fitclass = "GSLMultiMin";
+ fittype = "ConjugatePR";
+ } else if (!routine.compare("BFGS")) {
+ fitclass = "GSLMultiMin";
+ fittype = "BFGS";
+ } else if (!routine.compare("BFGS2")) {
+ fitclass = "GSLMultiMin";
+ fittype = "BFGS2";
+ } else if (!routine.compare("SteepDesc")) {
+ fitclass = "GSLMultiMin";
+ fittype = "SteepestDescent";
+ // } else if (!routine.compare("GSLMulti")) { fitclass = "GSLMultiFit";
+ // fittype = ""; // Doesn't work out of the box
+ } else if (!routine.compare("GSLSimAn")) {
+ fitclass = "GSLSimAn";
+ fittype = "";
+ } else if (!routine.compare("MCMC")) {
+ UseMCMC = true;
+ }
// make minimizer
if (fMinimizer) delete fMinimizer;
- fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype);
- fMinimizer->SetMaxFunctionCalls(FitPar::Config().GetParI("minimizer.maxcalls"));
+ if (UseMCMC) {
+ fMinimizer = new Simple_MH_Sampler();
+ } else {
+ fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype);
+ }
+
+ fMinimizer->SetMaxFunctionCalls(
+ FitPar::Config().GetParI("minimizer.maxcalls"));
if (!routine.compare("Brute")) {
fMinimizer->SetMaxFunctionCalls(fParams.size() * fParams.size() * 4);
fMinimizer->SetMaxIterations(fParams.size() * fParams.size() * 4);
}
- fMinimizer->SetMaxIterations(FitPar::Config().GetParI("minimizer.maxiterations"));
+ fMinimizer->SetMaxIterations(
+ FitPar::Config().GetParI("minimizer.maxiterations"));
fMinimizer->SetTolerance(FitPar::Config().GetParD("minimizer.tolerance"));
fMinimizer->SetStrategy(FitPar::Config().GetParI("minimizer.strategy"));
fMinimizer->SetFunction(*fCallFunctor);
int ipar = 0;
- //Add Fit Parameters
+ // Add Fit Parameters
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
bool fixed = true;
double vstart, vstep, vlow, vhigh;
vstart = vstep = vlow = vhigh = 0.0;
- if (fCurVals.find(syst) != fCurVals.end() ) vstart = fCurVals.at(syst);
- if (fMinVals.find(syst) != fMinVals.end() ) vlow = fMinVals.at(syst);
- if (fMaxVals.find(syst) != fMaxVals.end() ) vhigh = fMaxVals.at(syst);
- if (fStepVals.find(syst) != fStepVals.end()) vstep = fStepVals.at(syst);
- if (fFixVals.find(syst) != fFixVals.end() ) fixed = fFixVals.at(syst);
+ if (fCurVals.find(syst) != fCurVals.end()) vstart = fCurVals.at(syst);
+ if (fMinVals.find(syst) != fMinVals.end()) vlow = fMinVals.at(syst);
+ if (fMaxVals.find(syst) != fMaxVals.end()) vhigh = fMaxVals.at(syst);
+ if (fStepVals.find(syst) != fStepVals.end()) vstep = fStepVals.at(syst);
+ if (fFixVals.find(syst) != fFixVals.end()) fixed = fFixVals.at(syst);
// fix for errors
if (vhigh == vlow) vhigh += 1.0;
fMinimizer->SetVariable(ipar, syst, vstart, vstep);
fMinimizer->SetVariableLimits(ipar, vlow, vhigh);
if (fixed) {
-
fMinimizer->FixVariable(ipar);
LOG(FIT) << "Fixed Param: " << syst << std::endl;
} else {
-
- LOG(FIT) << "Free Param: " << syst
- << " Start:" << vstart
- << " Range:" << vlow << " to " << vhigh
- << " Step:" << vstep << std::endl;
+ LOG(FIT) << "Free Param: " << syst << " Start:" << vstart
+ << " Range:" << vlow << " to " << vhigh << " Step:" << vstep
+ << std::endl;
}
ipar++;
}
- LOG(FIT) << "Setup Minimizer: " << fMinimizer->NDim() << "(NDim) " << fMinimizer->NFree() << "(NFree)" << std::endl;
+ LOG(FIT) << "Setup Minimizer: " << fMinimizer->NDim() << "(NDim) "
+ << fMinimizer->NFree() << "(NFree)" << std::endl;
return;
}
//*************************************
// Set fake data from user input
void MinimizerRoutines::SetFakeData() {
-//*************************************
+ //*************************************
// If the fake data input field (-d) isn't provided, return to caller
if (fFakeDataInput.empty()) return;
// If user specifies -d MC we set the data to the MC
- // User can also specify fake data parameters to reweight by doing "fake_parameter" in input card file
+ // User can also specify fake data parameters to reweight by doing
+ // "fake_parameter" in input card file
// "fake_parameter" gets read in ReadCard function (reads to fFakeVals)
if (fFakeDataInput.compare("MC") == 0) {
-
LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl;
// fFakeVals get read in in ReadCard
UpdateRWEngine(fFakeVals);
// Reconfigure the reweight engine
FitBase::GetRW()->Reconfigure();
// Reconfigure all the samples to the new reweight
fSampleFCN->ReconfigureAllEvents();
// Feed on and set the fake-data in each measurement class
fSampleFCN->SetFakeData("MC");
// Changed the reweight engine values back to the current values
- // So we start the fit at a different value than what we set the fake-data to
+ // So we start the fit at a different value than what we set the fake-data
+ // to
UpdateRWEngine(fCurVals);
LOG(FIT) << "Set all data to fake MC predictions." << std::endl;
} else {
fSampleFCN->SetFakeData(fFakeDataInput);
}
return;
}
/*
Fitting Functions
*/
//*************************************
-void MinimizerRoutines::UpdateRWEngine(std::map& updateVals) {
-//*************************************
+void MinimizerRoutines::UpdateRWEngine(
+ std::map& updateVals) {
+ //*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name, updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void MinimizerRoutines::Run() {
-//*************************************
+ //*************************************
LOG(FIT) << "Running MinimizerRoutines : " << fStrategy << std::endl;
if (FitPar::Config().GetParB("save_nominal")) {
SaveNominal();
}
// Parse given routines
- fRoutines = GeneralUtils::ParseToStr(fStrategy,",");
- if (fRoutines.empty()){
- ERR(FTL) << "Trying to run MinimizerRoutines with no routines given!" << std::endl;
+ fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
+ if (fRoutines.empty()) {
+ ERR(FTL) << "Trying to run MinimizerRoutines with no routines given!"
+ << std::endl;
throw;
}
for (UInt_t i = 0; i < fRoutines.size(); i++) {
-
std::string routine = fRoutines.at(i);
int fitstate = kFitUnfinished;
LOG(FIT) << "Running Routine: " << routine << std::endl;
// Try Routines
- if (routine.find("LowStat") != std::string::npos) LowStatRoutine(routine);
- else if (routine == "FixAtLim") FixAtLimit();
- else if (routine == "FixAtLimBreak") fitstate = FixAtLimit();
- else if (routine.find("ErrorBands") != std::string::npos) GenerateErrorBands();
- else if (routine.find("DataToys") != std::string::npos) ThrowDataToys();
- else if (!routine.compare("Chi2Scan1D")) Create1DScans();
- else if (!routine.compare("Chi2Scan2D")) Chi2Scan2D();
- else fitstate = RunFitRoutine(routine);
+ if (routine.find("LowStat") != std::string::npos)
+ LowStatRoutine(routine);
+ else if (routine == "FixAtLim")
+ FixAtLimit();
+ else if (routine == "FixAtLimBreak")
+ fitstate = FixAtLimit();
+ else if (routine.find("ErrorBands") != std::string::npos)
+ GenerateErrorBands();
+ else if (routine.find("DataToys") != std::string::npos)
+ ThrowDataToys();
+ else if (!routine.compare("Chi2Scan1D"))
+ Create1DScans();
+ else if (!routine.compare("Chi2Scan2D"))
+ Chi2Scan2D();
+ else
+ fitstate = RunFitRoutine(routine);
// If ending early break here
if (fitstate == kFitFinished || fitstate == kNoChange) {
LOG(FIT) << "Ending fit routines loop." << std::endl;
break;
}
}
return;
}
//*************************************
int MinimizerRoutines::RunFitRoutine(std::string routine) {
-//*************************************
+ //*************************************
int endfits = kFitUnfinished;
// set fitter at the current start values
fOutputRootFile->cd();
SetupFitter(routine);
// choose what to do with the minimizer depending on routine.
- if (!routine.compare("Migrad") or
- !routine.compare("Simplex") or
- !routine.compare("Combined") or
- !routine.compare("Brute") or
- !routine.compare("Fumili") or
- !routine.compare("ConjugateFR") or
- !routine.compare("ConjugatePR") or
- !routine.compare("BFGS") or
- !routine.compare("BFGS2") or
- !routine.compare("SteepDesc") or
- // !routine.compare("GSLMulti") or
- !routine.compare("GSLSimAn")) {
-
+ if (!routine.compare("Migrad") or !routine.compare("Simplex") or
+ !routine.compare("Combined") or !routine.compare("Brute") or
+ !routine.compare("Fumili") or !routine.compare("ConjugateFR") or
+ !routine.compare("ConjugatePR") or !routine.compare("BFGS") or
+ !routine.compare("BFGS2") or !routine.compare("SteepDesc") or
+ // !routine.compare("GSLMulti") or
+ !routine.compare("GSLSimAn") or !routine.compare("MCMC")) {
if (fMinimizer->NFree() > 0) {
LOG(FIT) << fMinimizer->Minimize() << std::endl;
GetMinimizerState();
}
}
// other otptions
else if (!routine.compare("Contour")) {
CreateContours();
}
return endfits;
}
//*************************************
void MinimizerRoutines::PrintState() {
-//*************************************
+ //*************************************
LOG(FIT) << "------------" << std::endl;
// Count max size
int maxcount = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
maxcount = max(int(fParams[i].size()), maxcount);
}
// Header
LOG(FIT) << " # " << left << setw(maxcount) << "Parameter "
- << " = "
- << setw(10) << "Value" << " +- "
- << setw(10) << "Error" << " "
- << setw(8) << "(Units)" << " "
- << setw(10) << "Conv. Val" << " +- "
- << setw(10) << "Conv. Err" << " "
- << setw(8) << "(Units)" << std::endl;
+ << " = " << setw(10) << "Value"
+ << " +- " << setw(10) << "Error"
+ << " " << setw(8) << "(Units)"
+ << " " << setw(10) << "Conv. Val"
+ << " +- " << setw(10) << "Conv. Err"
+ << " " << setw(8) << "(Units)" << std::endl;
// Parameters
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
- std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
+ std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
std::string curunits = "(sig.)";
- double curval = fCurVals[syst];
- double curerr = fErrorVals[syst];
+ double curval = fCurVals[syst];
+ double curerr = fErrorVals[syst];
if (fStateVals[syst].find("ABS") != std::string::npos) {
curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
curunits = "(Abs.)";
} else if (fStateVals[syst].find("FRAC") != std::string::npos) {
curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
FitBase::RWSigmaToFrac(typestr, syst, 0.0));
curunits = "(Frac)";
}
std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
- double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
- double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
- FitBase::RWSigmaToAbs(typestr, syst, 0.0));
+ double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
+ double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
+ FitBase::RWSigmaToAbs(typestr, syst, 0.0));
std::ostringstream curparstring;
- curparstring << " " << setw(3) << left
- << i << ". "
- << setw(maxcount) << syst << " = "
- << setw(10) << curval << " +- "
- << setw(10) << curerr << " "
- << setw(8) << curunits << " "
- << setw(10) << convval << " +- "
- << setw(10) << converr << " "
- << setw(8) << convunits;
-
+ curparstring << " " << setw(3) << left << i << ". " << setw(maxcount)
+ << syst << " = " << setw(10) << curval << " +- " << setw(10)
+ << curerr << " " << setw(8) << curunits << " " << setw(10)
+ << convval << " +- " << setw(10) << converr << " " << setw(8)
+ << convunits;
LOG(FIT) << curparstring.str() << std::endl;
}
LOG(FIT) << "------------" << std::endl;
double like = fSampleFCN->GetLikelihood();
- LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
+ LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like
+ << std::endl;
LOG(FIT) << "------------" << std::endl;
}
//*************************************
void MinimizerRoutines::GetMinimizerState() {
-//*************************************
+ //*************************************
LOG(FIT) << "Minimizer State: " << std::endl;
// Get X and Err
- const double *values = fMinimizer->X();
- const double *errors = fMinimizer->Errors();
+ const double* values = fMinimizer->X();
+ const double* errors = fMinimizer->Errors();
// int ipar = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
- fCurVals[syst] = values[i];
+ fCurVals[syst] = values[i];
fErrorVals[syst] = errors[i];
}
PrintState();
// Covar
SetupCovariance();
if (fMinimizer->CovMatrixStatus() > 0) {
-
// Fill Full Covar
std::cout << "Filling covariance" << std::endl;
for (int i = 0; i < fCovar->GetNbinsX(); i++) {
for (int j = 0; j < fCovar->GetNbinsY(); j++) {
fCovar->SetBinContent(i + 1, j + 1, fMinimizer->CovMatrix(i, j));
}
}
int freex = 0;
int freey = 0;
for (int i = 0; i < fCovar->GetNbinsX(); i++) {
if (fMinimizer->IsFixedVariable(i)) continue;
freey = 0;
for (int j = 0; j < fCovar->GetNbinsY(); j++) {
if (fMinimizer->IsFixedVariable(j)) continue;
- fCovFree->SetBinContent(freex + 1, freey + 1, fMinimizer->CovMatrix(i, j));
+ fCovFree->SetBinContent(freex + 1, freey + 1,
+ fMinimizer->CovMatrix(i, j));
freey++;
}
freex++;
}
- fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation");
- fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition");
+ fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation");
+ fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition");
if (fMinimizer->NFree() > 0) {
fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free");
fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free");
}
}
std::cout << "Got STATE" << std::endl;
return;
};
//*************************************
void MinimizerRoutines::LowStatRoutine(std::string routine) {
-//*************************************
+ //*************************************
LOG(FIT) << "Running Low Statistics Routine: " << routine << std::endl;
int lowstatsevents = FitPar::Config().GetParI("minimizer.lowstatevents");
- int maxevents = FitPar::Config().GetParI("input.maxevents");
- int verbosity = FitPar::Config().GetParI("VERBOSITY");
+ int maxevents = FitPar::Config().GetParI("input.maxevents");
+ int verbosity = FitPar::Config().GetParI("VERBOSITY");
std::string trueroutine = routine;
std::string substring = "LowStat";
- trueroutine.erase( trueroutine.find(substring),
- substring.length() );
+ trueroutine.erase(trueroutine.find(substring), substring.length());
// Set MAX EVENTS=1000
FitPar::Config().SetParI("input.maxevents", lowstatsevents);
FitPar::Config().SetParI("VERBOSITY", 3);
SetupFCN();
RunFitRoutine(trueroutine);
FitPar::Config().SetParI("input.maxevents", maxevents);
SetupFCN();
FitPar::Config().SetParI("VERBOSITY", verbosity);
return;
}
//*************************************
void MinimizerRoutines::Create1DScans() {
-//*************************************
+ //*************************************
// 1D Scan Routine
// Steps through all free parameters about nominal using the step size
// Creates a graph for each free parameter
// At the current point create a 1D Scan for all parametes (Uncorrelated)
for (UInt_t i = 0; i < fParams.size(); i++) {
-
if (fFixVals[fParams[i]]) continue;
LOG(FIT) << "Running 1D Scan for " << fParams[i] << std::endl;
- fSampleFCN->CreateIterationTree(fParams[i] +
- "_scan1D_iterations",
+ fSampleFCN->CreateIterationTree(fParams[i] + "_scan1D_iterations",
FitBase::GetRW());
double scanmiddlepoint = fCurVals[fParams[i]];
// Determine N points needed
- double limlow = fMinVals[fParams[i]];
+ double limlow = fMinVals[fParams[i]];
double limhigh = fMaxVals[fParams[i]];
- double step = fStepVals[fParams[i]];
+ double step = fStepVals[fParams[i]];
- int npoints = int( fabs(limhigh - limlow) / (step + 0.) );
+ int npoints = int(fabs(limhigh - limlow) / (step + 0.));
- TH1D* contour = new TH1D(("Chi2Scan1D_" + fParams[i]).c_str(),
- ("Chi2Scan1D_" + fParams[i] +
- ";" + fParams[i]).c_str(),
- npoints, limlow, limhigh);
+ TH1D* contour =
+ new TH1D(("Chi2Scan1D_" + fParams[i]).c_str(),
+ ("Chi2Scan1D_" + fParams[i] + ";" + fParams[i]).c_str(),
+ npoints, limlow, limhigh);
// Fill bins
for (int x = 0; x < contour->GetNbinsX(); x++) {
-
// Set X Val
fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1);
// Run Eval
- double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
- double chi2 = fSampleFCN->DoEval( vals );
+ double* vals = FitUtils::GetArrayFromMap(fParams, fCurVals);
+ double chi2 = fSampleFCN->DoEval(vals);
delete vals;
// Fill Contour
contour->SetBinContent(x + 1, chi2);
}
// Save contour
contour->Write();
// Reset Parameter
fCurVals[fParams[i]] = scanmiddlepoint;
// Save TTree
fSampleFCN->WriteIterationTree();
}
return;
}
//*************************************
void MinimizerRoutines::Chi2Scan2D() {
-//*************************************
+ //*************************************
// Chi2 Scan 2D
// Creates a 2D chi2 scan by stepping through all free parameters
// Works for all pairwise combos of free parameters
// Scan I
for (UInt_t i = 0; i < fParams.size(); i++) {
if (fFixVals[fParams[i]]) continue;
// Scan J
for (UInt_t j = 0; j < i; j++) {
if (fFixVals[fParams[j]]) continue;
- fSampleFCN->CreateIterationTree( fParams[i] + "_" +
- fParams[j] + "_" +
- "scan2D_iterations",
- FitBase::GetRW() );
+ fSampleFCN->CreateIterationTree(
+ fParams[i] + "_" + fParams[j] + "_" + "scan2D_iterations",
+ FitBase::GetRW());
double scanmid_i = fCurVals[fParams[i]];
double scanmid_j = fCurVals[fParams[j]];
- double limlow_i = fMinVals[fParams[i]];
+ double limlow_i = fMinVals[fParams[i]];
double limhigh_i = fMaxVals[fParams[i]];
- double step_i = fStepVals[fParams[i]];
+ double step_i = fStepVals[fParams[i]];
- double limlow_j = fMinVals[fParams[j]];
+ double limlow_j = fMinVals[fParams[j]];
double limhigh_j = fMaxVals[fParams[j]];
- double step_j = fStepVals[fParams[j]];
+ double step_j = fStepVals[fParams[j]];
- int npoints_i = int( fabs(limhigh_i - limlow_i) / (step_i + 0.) ) + 1;
- int npoints_j = int( fabs(limhigh_j - limlow_j) / (step_j + 0.) ) + 1;
+ int npoints_i = int(fabs(limhigh_i - limlow_i) / (step_i + 0.)) + 1;
+ int npoints_j = int(fabs(limhigh_j - limlow_j) / (step_j + 0.)) + 1;
- TH2D* contour = new TH2D(("Chi2Scan2D_" + fParams[i] + "_" + fParams[j]).c_str(),
- ("Chi2Scan2D_" + fParams[i] + "_" + fParams[j] +
- ";" + fParams[i] + ";" + fParams[j]).c_str(),
- npoints_i, limlow_i, limhigh_i,
- npoints_j, limlow_j, limhigh_j );
+ TH2D* contour = new TH2D(
+ ("Chi2Scan2D_" + fParams[i] + "_" + fParams[j]).c_str(),
+ ("Chi2Scan2D_" + fParams[i] + "_" + fParams[j] + ";" + fParams[i] +
+ ";" + fParams[j])
+ .c_str(),
+ npoints_i, limlow_i, limhigh_i, npoints_j, limlow_j, limhigh_j);
// Begin Scan
- LOG(FIT) << "Running scan for " << fParams[i] << " " << fParams[j] << std::endl;
+ LOG(FIT) << "Running scan for " << fParams[i] << " " << fParams[j]
+ << std::endl;
// Fill bins
for (int x = 0; x < contour->GetNbinsX(); x++) {
-
// Set X Val
fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1);
// Loop Y
for (int y = 0; y < contour->GetNbinsY(); y++) {
-
// Set Y Val
fCurVals[fParams[j]] = contour->GetYaxis()->GetBinCenter(y + 1);
// Run Eval
- double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
- double chi2 = fSampleFCN->DoEval( vals );
+ double* vals = FitUtils::GetArrayFromMap(fParams, fCurVals);
+ double chi2 = fSampleFCN->DoEval(vals);
delete vals;
// Fill Contour
contour->SetBinContent(x + 1, y + 1, chi2);
fCurVals[fParams[j]] = scanmid_j;
}
fCurVals[fParams[i]] = scanmid_i;
fCurVals[fParams[j]] = scanmid_j;
}
// Save contour
contour->Write();
// Save Iterations
fSampleFCN->WriteIterationTree();
-
}
}
return;
}
//*************************************
void MinimizerRoutines::CreateContours() {
-//*************************************
+ //*************************************
// Use MINUIT for this if possible
- ERR(FTL) << " Contours not yet implemented as it is really slow!" << std::endl;
+ ERR(FTL) << " Contours not yet implemented as it is really slow!"
+ << std::endl;
throw;
return;
}
//*************************************
int MinimizerRoutines::FixAtLimit() {
-//*************************************
+ //*************************************
bool fixedparam = false;
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
if (fFixVals[syst]) continue;
double curVal = fCurVals.at(syst);
double minVal = fMinVals.at(syst);
double maxVal = fMinVals.at(syst);
if (fabs(curVal - minVal) < 0.0001) {
fCurVals[syst] = minVal;
fFixVals[syst] = true;
fixedparam = true;
}
if (fabs(maxVal - curVal) < 0.0001) {
fCurVals[syst] = maxVal;
fFixVals[syst] = true;
fixedparam = true;
}
}
if (!fixedparam) {
LOG(FIT) << "No dials needed fixing!" << std::endl;
return kNoChange;
- } else return kStateChange;
+ } else
+ return kStateChange;
}
-
/*
Write Functions
*/
//*************************************
void MinimizerRoutines::SaveResults() {
-//*************************************
+ //*************************************
fOutputRootFile->cd();
if (fMinimizer) {
SetupCovariance();
SaveMinimizerState();
}
SaveCurrentState();
-
}
//*************************************
void MinimizerRoutines::SaveMinimizerState() {
-//*************************************
+ //*************************************
std::cout << "Saving Minimizer State" << std::endl;
if (!fMinimizer) {
ERR(FTL) << "Can't save minimizer state without min object" << std::endl;
throw;
}
// Save main fit tree
fSampleFCN->WriteIterationTree();
-
+
// Get Vals and Errors
GetMinimizerState();
// Save tree with fit status
std::vector nameVect;
- std::vector valVect;
- std::vector errVect;
- std::vector minVect;
- std::vector maxVect;
- std::vector startVect;
- std::vector endfixVect;
- std::vector startfixVect;
+ std::vector valVect;
+ std::vector errVect;
+ std::vector minVect;
+ std::vector maxVect;
+ std::vector startVect;
+ std::vector endfixVect;
+ std::vector startfixVect;
// int NFREEPARS = fMinimizer->NFree();
int NPARS = fMinimizer->NDim();
int ipar = 0;
// Dial Vals
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams.at(i);
- nameVect .push_back( name );
+ nameVect.push_back(name);
- valVect .push_back( fCurVals.at(name) );
+ valVect.push_back(fCurVals.at(name));
- errVect .push_back( fErrorVals.at(name) );
+ errVect.push_back(fErrorVals.at(name));
- minVect .push_back( fMinVals.at(name) );
+ minVect.push_back(fMinVals.at(name));
- maxVect .push_back( fMaxVals.at(name) );
+ maxVect.push_back(fMaxVals.at(name));
- startVect .push_back( fStartVals.at(name) );
+ startVect.push_back(fStartVals.at(name));
- endfixVect .push_back( fFixVals.at(name) );
+ endfixVect.push_back(fFixVals.at(name));
- startfixVect.push_back( fStartFixVals.at(name) );
+ startfixVect.push_back(fStartFixVals.at(name));
ipar++;
}
int NFREE = fMinimizer->NFree();
- int NDIM = fMinimizer->NDim();
+ int NDIM = fMinimizer->NDim();
double CHI2 = fSampleFCN->GetLikelihood();
int NBINS = fSampleFCN->GetNDOF();
int NDOF = NBINS - NFREE;
// Write fit results
TTree* fit_tree = new TTree("fit_result", "fit_result");
fit_tree->Branch("parameter_names", &nameVect);
fit_tree->Branch("parameter_values", &valVect);
fit_tree->Branch("parameter_errors", &errVect);
fit_tree->Branch("parameter_min", &minVect);
fit_tree->Branch("parameter_max", &maxVect);
fit_tree->Branch("parameter_start", &startVect);
fit_tree->Branch("parameter_fix", &endfixVect);
fit_tree->Branch("parameter_startfix", &startfixVect);
fit_tree->Branch("CHI2", &CHI2, "CHI2/D");
fit_tree->Branch("NDOF", &NDOF, "NDOF/I");
fit_tree->Branch("NBINS", &NBINS, "NBINS/I");
fit_tree->Branch("NDIM", &NDIM, "NDIM/I");
fit_tree->Branch("NFREE", &NFREE, "NFREE/I");
fit_tree->Fill();
fit_tree->Write();
// Make dial variables
- TH1D dialvar = TH1D("fit_dials", "fit_dials", NPARS, 0, NPARS);
+ TH1D dialvar = TH1D("fit_dials", "fit_dials", NPARS, 0, NPARS);
TH1D startvar = TH1D("start_dials", "start_dials", NPARS, 0, NPARS);
- TH1D minvar = TH1D("min_dials", "min_dials", NPARS, 0, NPARS);
- TH1D maxvar = TH1D("max_dials", "max_dials", NPARS, 0, NPARS);
+ TH1D minvar = TH1D("min_dials", "min_dials", NPARS, 0, NPARS);
+ TH1D maxvar = TH1D("max_dials", "max_dials", NPARS, 0, NPARS);
- TH1D dialvarfree = TH1D("fit_dials_free", "fit_dials_free", NFREE, 0, NFREE);
- TH1D startvarfree = TH1D("start_dials_free", "start_dials_free", NFREE, 0, NFREE);
- TH1D minvarfree = TH1D("min_dials_free", "min_dials_free", NFREE, 0, NFREE);
- TH1D maxvarfree = TH1D("max_dials_free", "max_dials_free", NFREE, 0, NFREE);
+ TH1D dialvarfree = TH1D("fit_dials_free", "fit_dials_free", NFREE, 0, NFREE);
+ TH1D startvarfree =
+ TH1D("start_dials_free", "start_dials_free", NFREE, 0, NFREE);
+ TH1D minvarfree = TH1D("min_dials_free", "min_dials_free", NFREE, 0, NFREE);
+ TH1D maxvarfree = TH1D("max_dials_free", "max_dials_free", NFREE, 0, NFREE);
int freecount = 0;
for (UInt_t i = 0; i < nameVect.size(); i++) {
std::string name = nameVect.at(i);
dialvar.SetBinContent(i + 1, valVect.at(i));
dialvar.SetBinError(i + 1, errVect.at(i));
dialvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
startvar.SetBinContent(i + 1, startVect.at(i));
startvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
- minvar.SetBinContent(i + 1, minVect.at(i));
+ minvar.SetBinContent(i + 1, minVect.at(i));
minvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
- maxvar.SetBinContent(i + 1, maxVect.at(i));
+ maxvar.SetBinContent(i + 1, maxVect.at(i));
maxvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
if (NFREE > 0) {
if (!startfixVect.at(i)) {
freecount++;
dialvarfree.SetBinContent(freecount, valVect.at(i));
dialvarfree.SetBinError(freecount, errVect.at(i));
dialvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
startvarfree.SetBinContent(freecount, startVect.at(i));
startvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
- minvarfree.SetBinContent(freecount, minVect.at(i));
+ minvarfree.SetBinContent(freecount, minVect.at(i));
minvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
- maxvarfree.SetBinContent(freecount, maxVect.at(i));
+ maxvarfree.SetBinContent(freecount, maxVect.at(i));
maxvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
-
}
}
}
// Save Dial Plots
dialvar.Write();
startvar.Write();
minvar.Write();
maxvar.Write();
if (NFREE > 0) {
dialvarfree.Write();
startvarfree.Write();
minvarfree.Write();
maxvarfree.Write();
}
// Save fit_status plot
TH1D statusplot = TH1D("fit_status", "fit_status", 8, 0, 8);
- std::string fit_labels[8] = {"status", "cov_status", \
- "maxiter", "maxfunc", \
- "iter", "func", \
- "precision", "tolerance"
- };
+ std::string fit_labels[8] = {"status", "cov_status", "maxiter",
+ "maxfunc", "iter", "func",
+ "precision", "tolerance"};
double fit_vals[8];
fit_vals[0] = fMinimizer->Status() + 0.;
fit_vals[1] = fMinimizer->CovMatrixStatus() + 0.;
fit_vals[2] = fMinimizer->MaxIterations() + 0.;
fit_vals[3] = fMinimizer->MaxFunctionCalls() + 0.;
fit_vals[4] = fMinimizer->NIterations() + 0.;
fit_vals[5] = fMinimizer->NCalls() + 0.;
fit_vals[6] = fMinimizer->Precision() + 0.;
fit_vals[7] = fMinimizer->Tolerance() + 0.;
for (int i = 0; i < 8; i++) {
statusplot.SetBinContent(i + 1, fit_vals[i]);
statusplot.GetXaxis()->SetBinLabel(i + 1, fit_labels[i].c_str());
}
statusplot.Write();
// Save Covars
if (fCovar) fCovar->Write();
if (fCovFree) fCovFree->Write();
if (fCorrel) fCorrel->Write();
if (fCorFree) fCorFree->Write();
if (fDecomp) fDecomp->Write();
if (fDecFree) fDecFree->Write();
return;
}
//*************************************
void MinimizerRoutines::SaveCurrentState(std::string subdir) {
-//*************************************
+ //*************************************
LOG(FIT) << "Saving current full FCN predictions" << std::endl;
// Setup DIRS
TDirectory* curdir = gDirectory;
if (!subdir.empty()) {
- TDirectory* newdir = (TDirectory*) gDirectory->mkdir(subdir.c_str());
+ TDirectory* newdir = (TDirectory*)gDirectory->mkdir(subdir.c_str());
newdir->cd();
}
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void MinimizerRoutines::SaveNominal() {
-//*************************************
+ //*************************************
fOutputRootFile->cd();
LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl;
FitBase::GetRW()->Reconfigure();
SaveCurrentState("nominal");
-
};
//*************************************
void MinimizerRoutines::SavePrefit() {
-//*************************************
+ //*************************************
fOutputRootFile->cd();
LOG(FIT) << "Saving Prefit Predictions" << std::endl;
UpdateRWEngine(fStartVals);
SaveCurrentState("prefit");
UpdateRWEngine(fCurVals);
-
};
-
/*
MISC Functions
*/
//*************************************
int MinimizerRoutines::GetStatus() {
-//*************************************
+ //*************************************
return 0;
}
//*************************************
void MinimizerRoutines::SetupCovariance() {
-//*************************************
+ //*************************************
// Remove covares if they exist
if (fCovar) delete fCovar;
if (fCovFree) delete fCovFree;
if (fCorrel) delete fCorrel;
if (fCorFree) delete fCorFree;
if (fDecomp) delete fDecomp;
if (fDecFree) delete fDecFree;
LOG(FIT) << "Building covariance matrix.." << std::endl;
int NFREE = 0;
int NDIM = 0;
-
// Get NFREE from min or from vals (for cases when doing throws)
if (fMinimizer) {
std::cout << "NFREE FROM MINIMIZER" << std::endl;
NFREE = fMinimizer->NFree();
- NDIM = fMinimizer->NDim();
+ NDIM = fMinimizer->NDim();
} else {
NDIM = fParams.size();
for (UInt_t i = 0; i < fParams.size(); i++) {
std::cout << "Getting Param " << fParams[i] << std::endl;
if (!fFixVals[fParams[i]]) NFREE++;
}
}
if (NDIM == 0) return;
LOG(FIT) << "NFREE == " << NFREE << std::endl;
fCovar = new TH2D("covariance", "covariance", NDIM, 0, NDIM, NDIM, 0, NDIM);
if (NFREE > 0) {
- fCovFree = new TH2D("covariance_free",
- "covariance_free",
- NFREE, 0, NFREE,
+ fCovFree = new TH2D("covariance_free", "covariance_free", NFREE, 0, NFREE,
NFREE, 0, NFREE);
} else {
fCovFree = NULL;
}
// Set Bin Labels
int countall = 0;
int countfree = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
-
std::cout << "Getting Param " << i << std::endl;
std::cout << "ParamI = " << fParams[i] << std::endl;
fCovar->GetXaxis()->SetBinLabel(countall + 1, fParams[i].c_str());
fCovar->GetYaxis()->SetBinLabel(countall + 1, fParams[i].c_str());
countall++;
if (!fFixVals[fParams[i]] and NFREE > 0) {
fCovFree->GetXaxis()->SetBinLabel(countfree + 1, fParams[i].c_str());
fCovFree->GetYaxis()->SetBinLabel(countfree + 1, fParams[i].c_str());
countfree++;
}
}
std::cout << "Filling Matrices" << std::endl;
fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition");
if (NFREE > 0) {
fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free");
fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free");
} else {
fCorFree = NULL;
fDecFree = NULL;
}
std::cout << " Set the covariance" << std::endl;
return;
};
//*************************************
void MinimizerRoutines::ThrowCovariance(bool uniformly) {
-//*************************************
+ //*************************************
std::vector rands;
if (!fDecFree) {
ERR(WRN) << "Trying to throw 0 free parameters" << std::endl;
return;
}
// Generate Random Gaussians
for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) {
rands.push_back(gRandom->Gaus(0.0, 1.0));
}
// Reset Thrown Values
for (UInt_t i = 0; i < fParams.size(); i++) {
fThrownVals[fParams[i]] = fCurVals[fParams[i]];
}
// Loop and get decomp
for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) {
-
std::string parname = std::string(fDecFree->GetXaxis()->GetBinLabel(i + 1));
double mod = 0.0;
if (!uniformly) {
for (Int_t j = 0; j < fDecFree->GetNbinsY(); j++) {
mod += rands[j] * fDecFree->GetBinContent(j + 1, i + 1);
}
}
if (fCurVals.find(parname) != fCurVals.end()) {
-
- if (uniformly) fThrownVals[parname] = gRandom->Uniform(fMinVals[parname], fMaxVals[parname]);
- else { fThrownVals[parname] = fCurVals[parname] + mod; }
-
+ if (uniformly)
+ fThrownVals[parname] =
+ gRandom->Uniform(fMinVals[parname], fMaxVals[parname]);
+ else {
+ fThrownVals[parname] = fCurVals[parname] + mod;
+ }
}
}
// Check Limits
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams[i];
if (fFixVals[syst]) continue;
if (fThrownVals[syst] < fMinVals[syst]) fThrownVals[syst] = fMinVals[syst];
if (fThrownVals[syst] > fMaxVals[syst]) fThrownVals[syst] = fMaxVals[syst];
}
return;
};
//*************************************
void MinimizerRoutines::GenerateErrorBands() {
-//*************************************
+ //*************************************
- TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands");
+ TDirectory* errorDIR = (TDirectory*)fOutputRootFile->mkdir("error_bands");
errorDIR->cd();
- // Make a second file to store throws
+ // Make a second file to store throws
std::string tempFileName = fOutputFile;
- if (tempFileName.find(".root") != std::string::npos) tempFileName.erase(tempFileName.find(".root"), 5);
+ if (tempFileName.find(".root") != std::string::npos)
+ tempFileName.erase(tempFileName.find(".root"), 5);
tempFileName += ".throws.root";
- TFile* tempfile = new TFile(tempFileName.c_str(),"RECREATE");
+ TFile* tempfile = new TFile(tempFileName.c_str(), "RECREATE");
tempfile->cd();
int nthrows = FitPar::Config().GetParI("error_throws");
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
- TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal");
+ TDirectory* nominal = (TDirectory*)tempfile->mkdir("nominal");
nominal->cd();
fSampleFCN->Write();
-
- TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw");
+ TDirectory* outnominal = (TDirectory*)fOutputRootFile->mkdir("nominal_throw");
outnominal->cd();
fSampleFCN->Write();
-
errorDIR->cd();
TTree* parameterTree = new TTree("throws", "throws");
double chi2;
for (UInt_t i = 0; i < fParams.size(); i++)
- parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str());
+ parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]],
+ (fParams[i] + "/D").c_str());
parameterTree->Branch("chi2", &chi2, "chi2/D");
-
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Run Throws and save
for (Int_t i = 0; i < nthrows; i++) {
-
TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i", i));
throwfolder->cd();
// Generate Random Parameter Throw
ThrowCovariance(uniformly);
// Run Eval
- double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
- chi2 = fSampleFCN->DoEval( vals );
+ double* vals = FitUtils::GetArrayFromMap(fParams, fThrownVals);
+ chi2 = fSampleFCN->DoEval(vals);
delete vals;
// Save the FCN
fSampleFCN->Write();
parameterTree->Fill();
}
errorDIR->cd();
fDecFree->Write();
fCovFree->Write();
parameterTree->Write();
delete parameterTree;
- // Now go through the keys in the temporary file and look for TH1D, and TH2D plots
+ // Now go through the keys in the temporary file and look for TH1D, and TH2D
+ // plots
TIter next(nominal->GetListOfKeys());
- TKey *key;
+ TKey* key;
while ((key = (TKey*)next())) {
- TClass *cl = gROOT->GetClass(key->GetClassName());
+ TClass* cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
- TH1D *baseplot = (TH1D*)key->ReadObj();
+ TH1D* baseplot = (TH1D*)key->ReadObj();
std::string plotname = std::string(baseplot->GetName());
int nbins = baseplot->GetNbinsX() * baseplot->GetNbinsY();
// Setup TProfile with RMS option
- TProfile* tprof = new TProfile((plotname + "_prof").c_str(), (plotname + "_prof").c_str(), nbins, 0, nbins, "S");
+ TProfile* tprof =
+ new TProfile((plotname + "_prof").c_str(), (plotname + "_prof").c_str(),
+ nbins, 0, nbins, "S");
// Setup The TTREE
double* bincontents;
bincontents = new double[nbins];
double* binlowest;
binlowest = new double[nbins];
double* binhighest;
binhighest = new double[nbins];
errorDIR->cd();
- TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
+ TTree* bintree =
+ new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
for (Int_t i = 0; i < nbins; i++) {
bincontents[i] = 0.0;
binhighest[i] = 0.0;
binlowest[i] = 0.0;
- bintree->Branch(Form("content_%i", i), &bincontents[i], Form("content_%i/D", i));
+ bintree->Branch(Form("content_%i", i), &bincontents[i],
+ Form("content_%i/D", i));
}
for (Int_t i = 0; i < nthrows; i++) {
- TH1* newplot = (TH1*)tempfile->Get(Form(("throw_%i/" + plotname).c_str(), i));
+ TH1* newplot =
+ (TH1*)tempfile->Get(Form(("throw_%i/" + plotname).c_str(), i));
for (Int_t j = 0; j < nbins; j++) {
tprof->Fill(j + 0.5, newplot->GetBinContent(j + 1));
bincontents[j] = newplot->GetBinContent(j + 1);
- if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j];
- if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j];
+ if (bincontents[j] < binlowest[j] or i == 0)
+ binlowest[j] = bincontents[j];
+ if (bincontents[j] > binhighest[j] or i == 0)
+ binhighest[j] = bincontents[j];
}
errorDIR->cd();
bintree->Fill();
delete newplot;
}
errorDIR->cd();
for (Int_t j = 0; j < nbins; j++) {
-
if (!uniformly) {
baseplot->SetBinError(j + 1, tprof->GetBinError(j + 1));
} else {
baseplot->SetBinContent(j + 1, (binlowest[j] + binhighest[j]) / 2.0);
baseplot->SetBinError(j + 1, (binhighest[j] - binlowest[j]) / 2.0);
}
}
errorDIR->cd();
baseplot->Write();
tprof->Write();
bintree->Write();
delete baseplot;
delete tprof;
delete bintree;
- delete [] bincontents;
+ delete[] bincontents;
}
return;
};
-
-void MinimizerRoutines::ThrowDataToys(){
+void MinimizerRoutines::ThrowDataToys() {
LOG(FIT) << "Generating Toy Data Throws" << std::endl;
int verb = Config::Get().GetParI("VERBOSITY");
SETVERBOSITY(FIT);
int nthrows = FitPar::Config().GetParI("NToyThrows");
double maxlike = -1.0;
double minlike = -1.0;
std::vector values;
- for (int i = 0; i < 1.E4; i++){
+ for (int i = 0; i < 1.E4; i++) {
fSampleFCN->ThrowDataToy();
double like = fSampleFCN->GetLikelihood();
values.push_back(like);
if (maxlike == -1.0 or like > maxlike) maxlike = like;
if (minlike == -1.0 or like < minlike) minlike = like;
}
SETVERBOSITY(verb);
-
// Fill Histogram
- TH1D* likes =new TH1D("toydatalikelihood","toydatalikelihood",int(sqrt(nthrows)),minlike,maxlike);
- for (size_t i = 0; i < values.size(); i++){
+ TH1D* likes = new TH1D("toydatalikelihood", "toydatalikelihood",
+ int(sqrt(nthrows)), minlike, maxlike);
+ for (size_t i = 0; i < values.size(); i++) {
likes->Fill(values[i]);
}
// Save to file
LOG(FIT) << "Writing toy data throws" << std::endl;
fOutputRootFile->cd();
likes->Write();
}
diff --git a/src/Routines/Simple_MH_Sampler.h b/src/Routines/Simple_MH_Sampler.h
new file mode 100644
index 0000000..48d43fa
--- /dev/null
+++ b/src/Routines/Simple_MH_Sampler.h
@@ -0,0 +1,334 @@
+#include "Math/Minimizer.h"
+
+#include "FitLogger.h"
+
+using ROOT::Math::Minimizer;
+
+class Simple_MH_Sampler : public Minimizer {
+ TRandom3 RNJesus;
+
+ size_t step_i;
+ int moved;
+
+ size_t thin;
+ size_t thin_ctr;
+
+ size_t discard;
+
+ struct Param {
+ Param()
+ : IsFixed(false),
+ name(""),
+ Val(0xdeadbeef),
+ StepWidth(0xdeadbeef),
+ LowLim(0xdeadbeef),
+ UpLim(0xdeadbeef) {}
+ Param(bool i, std::string n, double v, double s, double l, double u)
+ : IsFixed(i), name(n), Val(v), StepWidth(s), LowLim(l), UpLim(u) {}
+ bool IsFixed;
+ std::string name;
+ double Val, StepWidth, LowLim, UpLim;
+ };
+
+ std::vector start_params;
+
+ double curr_value;
+ std::vector curr_params;
+
+ double propose_value;
+ std::vector propose_params;
+
+ double min_value;
+ std::vector min_params;
+
+ TGraph trace;
+
+ void RestartParams() {
+ curr_params.resize(start_params.size());
+ for (size_t p_it = 0; p_it < start_params.size(); ++p_it) {
+ curr_params[p_it] = start_params[p_it].Val;
+ }
+ min_params = curr_params;
+ propose_params = curr_params;
+ }
+
+ TTree *StepTree;
+
+ void Write();
+
+ ROOT::Math::IMultiGenFunction const *FCN;
+
+ public:
+ Simple_MH_Sampler() : Minimizer(), RNJesus(), trace() {
+ thin = Config::Get().ConfI("MCMC.thin");
+ thin_ctr = 0;
+ discard = Config::Get().ConfI("MCMC.BurnInSteps");
+ }
+
+ void SetFunction(ROOT::Math::IMultiGenFunction const &func) { FCN = &func; }
+
+ bool SetVariable(unsigned int ivar, std::string const &name, double val,
+ double step) {
+ if (start_params.size() <= ivar) {
+ start_params.resize(ivar + 1);
+ }
+ start_params[ivar] = Param(false, name, val, step, 0xdeadbeef, 0xdeadbeef);
+ RestartParams();
+ return true;
+ }
+
+ bool SetLowerLimitedVariable(unsigned int ivar, std::string const &name,
+ double val, double step, double lower) {
+ SetVariable(ivar, name, val, step);
+ start_params[ivar].LowLim = lower;
+ return true;
+ }
+ bool SetUpperLimitedVariable(unsigned int ivar, std::string const &name,
+ double val, double step, double upper) {
+ SetVariable(ivar, name, val, step);
+ start_params[ivar].UpLim = upper;
+ return true;
+ }
+ bool SetLimitedVariable(unsigned int ivar, std::string const &name,
+ double val, double step, double lower, double upper) {
+ SetLowerLimitedVariable(ivar, name, val, step, lower);
+ start_params[ivar].UpLim = upper;
+ return true;
+ }
+ bool SetFixedVariable(unsigned int ivar, std::string const &name,
+ double val) {
+ SetVariable(ivar, name, val, 0xdeadbeef);
+ start_params[ivar].IsFixed = true;
+ return true;
+ }
+ bool SetVariableValue(unsigned int ivar, double value) {
+ if (start_params.size() <= ivar) {
+ ERROR(WRN, "Tried to set uninitialised variable.");
+ return false;
+ }
+ start_params[ivar].Val = value;
+ return true;
+ }
+ bool SetVariableStepSize(unsigned int ivar, double value) {
+ if (start_params.size() <= ivar) {
+ ERROR(WRN, "Tried to set uninitialised variable.");
+ return false;
+ }
+ start_params[ivar].StepWidth = value;
+ return true;
+ }
+ bool SetVariableLowerLimit(unsigned int ivar, double lower) {
+ if (start_params.size() <= ivar) {
+ ERROR(WRN, "Tried to set uninitialised variable.");
+ return false;
+ }
+ start_params[ivar].LowLim = lower;
+ return true;
+ }
+ bool SetVariableUpperLimit(unsigned int ivar, double upper) {
+ if (start_params.size() <= ivar) {
+ ERROR(WRN, "Tried to set uninitialised variable.");
+ return false;
+ }
+ start_params[ivar].UpLim = upper;
+ return true;
+ }
+ bool SetVariableLimits(unsigned int ivar, double lower, double upper) {
+ SetVariableLowerLimit(ivar, lower);
+ SetVariableUpperLimit(ivar, upper);
+ return true;
+ }
+ bool FixVariable(unsigned int ivar) {
+ if (start_params.size() <= ivar) {
+ ERROR(WRN, "Tried to fix uninitialised variable.");
+ return false;
+ }
+ start_params[ivar].IsFixed = true;
+ return true;
+ }
+ bool ReleaseVariable(unsigned int ivar) {
+ if (start_params.size() <= ivar) {
+ ERROR(WRN, "Tried to fix uninitialised variable.");
+ return false;
+ }
+ start_params[ivar].IsFixed = false;
+ return true;
+ }
+ bool IsFixedVariable(unsigned int ivar) {
+ if (start_params.size() <= ivar) {
+ ERROR(WRN, "Tried to fix uninitialised variable.");
+ return false;
+ }
+ return start_params[ivar].IsFixed;
+ }
+
+ double MinValue() const { return min_value; }
+ const double *X() const { return min_params.data(); }
+ const double *Errors() const { return min_params.data(); }
+
+ unsigned int NDim() const { return start_params.size(); }
+
+ unsigned int NFree() const {
+ unsigned int NFree = 0;
+
+ for (size_t p_it = 0; p_it < start_params.size(); ++p_it) {
+ NFree += !start_params[p_it].IsFixed;
+ }
+ return NFree;
+ }
+
+ void AddBranches() {
+ TFile *ogf = gFile;
+ if (FitPar::Config().out && FitPar::Config().out->IsOpen()) {
+ FitPar::Config().out->cd();
+ }
+
+ StepTree = new TTree("MCMChain", "");
+ StepTree->Branch("Step", &step_i, "Step/I");
+ StepTree->Branch("Value", &curr_value, "Value/D");
+ StepTree->Branch("Moved", &moved, "Moved/I");
+
+ std::stringstream ss("");
+ for (size_t p_it = 0; p_it < curr_params.size(); ++p_it) {
+ ss.str("");
+ ss << "param_" << p_it;
+ StepTree->Branch(ss.str().c_str(), &curr_params[p_it],
+ (ss.str() + "/D").c_str());
+ }
+
+ if (ogf && ogf->IsOpen()) {
+ ogf->cd();
+ }
+ }
+
+ void Fill() { StepTree->Fill(); }
+
+ void Propose() {
+ for (size_t p_it = 0; p_it < start_params.size(); ++p_it) {
+ double propose_param = curr_params[p_it];
+
+ if (!start_params[p_it].IsFixed) {
+ size_t attempts = 0;
+ do {
+ if (attempts > 1000) {
+ THROW("After 1000 attempts, failed to throw Gaus("
+ << start_params[p_it].Val << ", "
+ << start_params[p_it].StepWidth << ") inside limits: [ "
+ << start_params[p_it].LowLim << " -- "
+ << start_params[p_it].UpLim << " ]");
+ }
+
+ double thr =
+ RNJesus.Gaus(curr_params[p_it], start_params[p_it].StepWidth);
+
+ if ((start_params[p_it].LowLim != 0xdeadbeef) &&
+ (thr < start_params[p_it].LowLim)) {
+ attempts++;
+ continue;
+ }
+
+ if ((start_params[p_it].UpLim != 0xdeadbeef) &&
+ (thr > start_params[p_it].UpLim)) {
+ attempts++;
+ continue;
+ }
+
+ propose_param = thr;
+ break;
+
+ } while (true);
+ }
+
+ propose_params[p_it] = propose_param;
+ }
+ }
+
+ void Evaluate() {
+ propose_value = exp(-(*FCN)(propose_params.data()) / 10000.0);
+ if (propose_value < min_value) {
+ min_params = propose_params;
+ }
+ }
+
+ void PrintResults() {
+ QLOG(FIT, "Simple_MH_Sampler State: ");
+ for (size_t p_it = 0; p_it < start_params.size(); ++p_it) {
+ QLOG(FIT, "\t[" << p_it
+ << "]: " << (start_params[p_it].IsFixed ? " FIX" : "FREE")
+ << " " << curr_params[p_it]);
+ }
+ QLOG(FIT, "Curr LHood: " << curr_value << ", Min LHood: " << min_value);
+ }
+
+ void Step() {
+ moved = false;
+
+ if (propose_value != propose_value) {
+ curr_params = propose_params;
+ curr_value = propose_value;
+ PrintResults();
+ THROW("Proposed a NAN value.");
+ }
+
+ std::cout << "[" << step_i << "] proposed: " << propose_value
+ << " | current: " << curr_value << std::endl;
+ double a = propose_value / curr_value;
+ std::cout << "\ta = " << a << std::endl;
+ if (a >= 1.0) {
+ moved = true;
+ std::cout << "\tMoved." << std::endl;
+ } else {
+ double b = RNJesus.Uniform(1);
+ if (b < a) {
+ moved = true;
+ std::cout << "\tMoved (" << b << ")" << std::endl;
+ } else {
+ std::cout << "\tStayed. (" << b << ")" << std::endl;
+ }
+ }
+
+ if (moved) {
+ curr_params = propose_params;
+ curr_value = propose_value;
+ }
+ }
+
+ bool Minimize() {
+ if (!start_params.size()) {
+ ERROR(FTL, "No Parameters passed to Simple_MH_Sampler.");
+ return false;
+ }
+
+ RestartParams();
+
+ AddBranches();
+
+ size_t NSteps = Options().MaxIterations();
+ trace.Set(NSteps);
+ QLOG(FIT, "Running chain for " << NSteps << " steps.");
+ step_i = 0;
+ while (step_i < NSteps) {
+ Propose();
+
+ Evaluate();
+
+ Step();
+
+ trace.SetPoint(step_i, step_i, curr_value);
+
+ if (step_i >= discard) {
+ thin_ctr++;
+ if (thin_ctr == thin) {
+ Fill();
+ thin_ctr = 0;
+ }
+ }
+ step_i++;
+ }
+
+ StepTree->Write();
+ trace.Write("MCMCTrace");
+
+ return true;
+ };
+};