diff --git a/parameters/config.xml b/parameters/config.xml
index 4e103d0..41cad23 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,214 +1,202 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-
+
-
-
-
+
+
+
-
+
-
-
-
+
+
+
-
diff --git a/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx b/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx
index fde3809..fc41380 100755
--- a/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx
+++ b/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx
@@ -1,153 +1,153 @@
// 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 "MINERvA_SignalDef.h"
#include "MINERvA_CCinc_XSec_2DEavq3_nu.h"
//********************************************************************
MINERvA_CCinc_XSec_2DEavq3_nu::MINERvA_CCinc_XSec_2DEavq3_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCinc_XSec_2DEavq3_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Medium Energy FHC numu \n" \
"Signal: CC-inclusive with theta < 20deg \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("q_{3} (GeV)");
fSettings.SetYTitle("E_{avail} (GeV)");
fSettings.SetZTitle("d^{2}#sigma/dq_{3}dE_{avail} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
fSettings.SetEnuRange(2.0, 6.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CCinc_XSec_2DEavq3_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/data_2D.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/covar_2D.txt" );
fSettings.SetMapInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/map_2D.txt" );
fSettings.DefineAllowedSpecies("numu");
- hadroncut = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu.hadron_cut");
- useq3true = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu.useq3true");
- splitMEC_PN_NN = FitPar::Config().GetParB("Modes.split_PN_NN");
+ hadroncut = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu_hadron_cut");
+ useq3true = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu_useq3true");
+ splitMEC_PN_NN = FitPar::Config().GetParB("Modes_split_PN_NN");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-42 / (fNEvents + 0.)) / this->TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
Double_t binx[7] = {0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8};
Double_t biny[17] = {0.0, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.25, 0.30, 0.35, 0.40, 0.50, 0.60, 0.80};
CreateDataHistogram(7, binx, 17, biny);
SetDataValuesFromTextFile( fSettings.GetDataInput() );
ScaleData(1E-42);
SetMapValuesFromText( fSettings.GetMapInput() );
SetCholDecompFromTextFile( fSettings.GetCovarInput(), 67);
ScaleCovar(1E-16);
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CCinc_XSec_2DEavq3_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Seperate MEC
if (splitMEC_PN_NN) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < event->Npart(); j++) {
if ((event->PartInfo(j))->fIsAlive) continue;
if (event->PartInfo(j)->fPID == 2212) npr++;
else if (event->PartInfo(j)->fPID == 2112) nne++;
}
if (event->Mode == 2 and npr == 1 and nne == 1) {
event->Mode = 2;
Mode = 2;
} else if (event->Mode == 2 and npr == 0 and nne == 2) {
event->Mode = 3;
Mode = 3;
}
}
// Set Defaults
double Eav = -999.9;
double q3 = -999.9;
// If muon found get kinematics
FitParticle* muon = event->GetHMFSParticle(13);
FitParticle* neutrino = event->GetNeutrinoIn();
if (muon && neutrino) {
// Set Q from Muon
TLorentzVector q = neutrino->fP - muon->fP;
double q0 = (q.E()) / 1.E3;
//double q3_true = (q.Vect().Mag())/1.E3;
double thmu = muon->fP.Vect().Angle(neutrino->fP.Vect());
double pmu = muon->fP.Vect().Mag() / 1.E3;
double emu = muon->fP.E() / 1.E3;
double mmu = muon->fP.Mag() / 1.E3;
// Get Enu Rec
double enu_rec = emu + q0;
// Set Q2 QE
double q2qe = 2 * enu_rec * (emu - pmu * cos(thmu)) - mmu * mmu;
// Calc Q3 from Q2QE and EnuTree
q3 = sqrt(q2qe + q0 * q0);
// Get Eav too
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event) / 1.E3;
}
// Set Hist Variables
fXVar = q3;
fYVar = Eav;
return;
}
//********************************************************************
bool MINERvA_CCinc_XSec_2DEavq3_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCCincLowRecoil_MINERvA(event, EnuMin, EnuMax);
}
diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx
index 8e0aab6..4586ae4 100644
--- a/src/Reweight/GENIEWeightEngine.cxx
+++ b/src/Reweight/GENIEWeightEngine.cxx
@@ -1,244 +1,244 @@
#include "GENIEWeightEngine.h"
GENIEWeightEngine::GENIEWeightEngine(std::string name) {
#ifdef __GENIE_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
fGenieRW = new genie::rew::GReWeight();
// Get List of Vetos (Just for debugging)
- std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fGenieRW_veto");
+ std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fGenieRW_veto");
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos;
bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;
// Now actually add the RW Calcs
if (xsec_ncel)
fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
if (xsec_ccqe){
fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
// (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe")))
// ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") );
}
if (xsec_coh){
fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH());
// (dynamic_cast (fGenieRW->WghtCalc("xsec_coh")))
// ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") );
}
if (xsec_nnres)
fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
new genie::rew::GReWeightNonResonanceBkg);
if (xsec_nudis)
fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
if (xsec_resdec)
fGenieRW->AdoptWghtCalc("hadro_res_decay",
new genie::rew::GReWeightResonanceDecay);
if (xsec_fzone)
fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
if (xsec_intra)
fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
if (xsec_agky)
fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
if (xsec_qevec)
fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
new genie::rew::GReWeightNuXSecCCQEvec);
#if __GENIE_VERSION__ >= 212
if (xsec_qeaxial)
fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
new genie::rew::GReWeightNuXSecCCQEaxial);
#endif
if (xsec_dis)
fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
if (xsec_nc)
fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
if (xsec_ccres){
#if __GENIE_VERSION__ < 213
fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
#else
fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES( FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default"));
#endif
}
if (xsec_ncres)
fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
if (xsec_nucqe)
fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);
GReWeightNuXSecCCQE * rwccqe =
dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe"));
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
// Default to include shape and normalization changes for CCRES (can be changed downstream if desired)
GReWeightNuXSecCCRES * rwccres =
dynamic_cast (fGenieRW->WghtCalc("xsec_ccres"));
std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode");
if (!marestype.compare("kModeNormAndMaMvShape")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); }
else if (!marestype.compare("kModeMaMv")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); }
else {
THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype );
}
// Default to include shape and normalization changes for NCRES (can be changed downstream if desired)
GReWeightNuXSecNCRES * rwncres =
dynamic_cast (fGenieRW->WghtCalc("xsec_ncres"));
rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
// Default to include shape and normalization changes for DIS (can be changed downstream if desired)
GReWeightNuXSecDIS * rwdis =
dynamic_cast (fGenieRW->WghtCalc("xsec_dis"));
rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl;
#endif
};
void GENIEWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __GENIE_ENABLED__
// Get First enum
int nuisenum = Reweight::ConvDial(name, kGENIE);
// Setup Maps
fEnumIndex[nuisenum];// = std::vector(0);
fNameIndex[name]; // = std::vector(0);
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fGENIESysts.push_back(rwsyst);
// Initialize dial
std::cout << "Registering " << singlename << " from " << name << std::endl;
fGenieRW->Systematics().Init( fGENIESysts[index] );
// If Absolute
if (fIsAbsTwk) {
GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 );
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
};
void GENIEWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __GENIE_ENABLED__
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val);
}
#endif
}
void GENIEWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __GENIE_ENABLED__
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
}
#endif
}
void GENIEWeightEngine::Reconfigure(bool silent) {
#ifdef __GENIE_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fGenieRW->Reconfigure();
fGenieRW->Print();
// Shout again
if (silent) StartTalking();
#endif
}
double GENIEWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __GENIE_ENABLED__
// Skip Non GENIE
if (evt->fType != kGENIE) return 1.0;
// Make nom weight
if (!evt) {
THROW("evt not found : " << evt);
}
if (!(evt->genie_event)) {
THROW("evt->genie_event not found!" << evt->genie_event);
}
if (!(evt->genie_event->event)) {
THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event));
}
if (!fGenieRW) {
THROW("GENIE RW Not Found!" << fGenieRW);
}
rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
// std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl;
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/NEUTWeightEngine.cxx b/src/Reweight/NEUTWeightEngine.cxx
index 3bd0c2e..f9aaa78 100644
--- a/src/Reweight/NEUTWeightEngine.cxx
+++ b/src/Reweight/NEUTWeightEngine.cxx
@@ -1,186 +1,186 @@
#include "NEUTWeightEngine.h"
NEUTWeightEngine::NEUTWeightEngine(std::string name) {
#ifdef __NEUT_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up NEUT RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
fNeutRW = new neut::rew::NReWeight();
TDirectory* olddir = gDirectory;
// get list of vetoed calc engines (just for debug really)
- std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fNeutRW_veto");
+ std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fNeutRW_veto");
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_res = rw_engine_list.find("xsec_res") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool nucl_casc = rw_engine_list.find("nucl_casc") == std::string::npos;
bool nucl_piless = rw_engine_list.find("nucl_piless") == std::string::npos;
// Activate each calc engine
if (xsec_ccqe)
fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE);
if (xsec_res)
fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES);
if (xsec_ccres)
fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES);
if (xsec_coh)
fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH);
if (xsec_dis)
fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS);
if (xsec_ncel)
fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL);
if (xsec_nc)
fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC);
if (xsec_ncres)
fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES);
if (nucl_casc)
fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc);
if (nucl_piless)
fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless);
fNeutRW->Reconfigure();
olddir->cd();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
THROW("NEUT RW NOT ENABLED!" );
#endif
};
void NEUTWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __NEUT_ENABLED__
// Get First enum
int nuisenum = Reweight::ConvDial(name, kNEUT);
// Setup Maps
fEnumIndex[nuisenum];// = std::vector(0);
fNameIndex[name]; // = std::vector(0);
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get Syst
neut::rew::NSyst_t gensyst = NSyst::FromString(singlename);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNEUTSysts.push_back(gensyst);
// Initialize dial
LOG(FIT) << "Registering " << singlename << " dial." << std::endl;
fNeutRW->Systematics().Init( fNEUTSysts[index] );
// If Absolute
if (fIsAbsTwk) {
NSystUncertainty::Instance()->SetUncertainty( fNEUTSysts[index], 1.0, 1.0 );
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
}
void NEUTWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __NEUT_ENABLED__
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
std::cout << "Setting Dial Value = " << nuisenum << " "
<< i << " " << indices[i] << " " << fValues[indices[i]]
<< " Enum=" << fNEUTSysts[indices[i]] << std::endl;
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __NEUT_ENABLED__
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
std::cout << "Setting Dial Value = " << name << " = "
<< i << " " << indices[i] << " " << fValues[indices[i]]
<< " Enum=" << fNEUTSysts[indices[i]] << std::endl;
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::Reconfigure(bool silent) {
#ifdef __NEUT_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fNeutRW->Reconfigure();
//if (LOG_LEVEL(DEB)){
fNeutRW->Print();
// }
// Shout again
if (silent) StartTalking();
#endif
}
double NEUTWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __NEUT_ENABLED__
// Skip Non NEUT
if (evt->fType != kNEUT) return 1.0;
// Hush now
StopTalking();
// Fill NEUT Common blocks
NEUTUtils::FillNeutCommons(evt->fNeutVect);
// Call Weight calculation
rw_weight = fNeutRW->CalcWeight();
// Speak Now
StartTalking();
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/NIWGWeightEngine.cxx b/src/Reweight/NIWGWeightEngine.cxx
index 94c53ec..14691b9 100644
--- a/src/Reweight/NIWGWeightEngine.cxx
+++ b/src/Reweight/NIWGWeightEngine.cxx
@@ -1,220 +1,220 @@
#include "NIWGWeightEngine.h"
NIWGWeightEngine::NIWGWeightEngine(std::string name) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up NIWG RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
fNIWGRW = new niwg::rew::NIWGReWeight();
// Get List of Veto Calcs (For Debugging)
std::string rw_engine_list =
- FitPar::Config().GetParS("FitWeight.fNIWGRW_veto");
+ FitPar::Config().GetParS("FitWeight_fNIWGRW_veto");
bool niwg_2012a = rw_engine_list.find("niwg_2012a") == std::string::npos;
bool niwg_2014a = rw_engine_list.find("niwg_2014a") == std::string::npos;
bool niwg_pimult = rw_engine_list.find("niwg_pimult") == std::string::npos;
bool niwg_mec = rw_engine_list.find("niwg_mec") == std::string::npos;
bool niwg_rpa = rw_engine_list.find("niwg_rpa") == std::string::npos;
bool niwg_eff_rpa = rw_engine_list.find("niwg_eff_rpa") == std::string::npos;
bool niwg_proton =
rw_engine_list.find("niwg_protonFSIbug") == std::string::npos;
bool niwg_hadron =
rw_engine_list.find("niwg_HadronMultSwitch") == std::string::npos;
// Add the RW Calcs
if (niwg_2012a)
fNIWGRW->AdoptWghtCalc("niwg_2012a", new niwg::rew::NIWGReWeight2012a);
if (niwg_2014a)
fNIWGRW->AdoptWghtCalc("niwg_2014a", new niwg::rew::NIWGReWeight2014a);
if (niwg_pimult)
fNIWGRW->AdoptWghtCalc("niwg_pimult", new niwg::rew::NIWGReWeightPiMult);
if (niwg_mec)
fNIWGRW->AdoptWghtCalc("niwg_mec", new niwg::rew::NIWGReWeightMEC);
if (niwg_rpa)
fNIWGRW->AdoptWghtCalc("niwg_rpa", new niwg::rew::NIWGReWeightRPA);
if (niwg_eff_rpa)
fNIWGRW->AdoptWghtCalc("niwg_eff_rpa",
new niwg::rew::NIWGReWeightEffectiveRPA);
if (niwg_proton)
fNIWGRW->AdoptWghtCalc("niwg_protonFSIbug",
new niwg::rew::NIWGReWeightProtonFSIbug);
if (niwg_hadron)
fNIWGRW->AdoptWghtCalc("niwg_HadronMultSwitch",
new niwg::rew::NIWGReWeightHadronMultSwitch);
fNIWGRW->Reconfigure();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
ERR(FTL) << "NIWG RW Enabled but NEUT RW is not!" << std::endl;
#endif
#else
ERR(FTL) << "NIWG RW NOT ENABLED!" << std::endl;
#endif
};
void NIWGWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
// Get First enum
int nuisenum = Reweight::ConvDial(name, kNIWG);
// Setup Maps
fEnumIndex[nuisenum];// = std::vector(0);
fNameIndex[name]; // = std::vector(0);
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
niwg::rew::NIWGSyst_t gensyst = niwg::rew::NIWGSyst::FromString(name);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNIWGSysts.push_back(gensyst);
// Initialize dial
LOG(FIT) << "Registering " << singlename << " from " << name << std::endl;
fNIWGRW->Systematics().Init( fNIWGSysts[index] );
// If Absolute
if (fIsAbsTwk) {
niwg::rew::NIWGSystUncertainty::Instance()->SetUncertainty( fNIWGSysts[index], 1.0, 1.0 );
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
#endif
}
void NIWGWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
std::cout << "Setting NIWG Dial Value "
<< i << " " << indices[i] << " "
<< fNIWGSysts[indices[i]] << std::endl;
fNIWGRW->Systematics().Set( fNIWGSysts[indices[i]], val);
}
#endif
#endif
}
void NIWGWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fNIWGRW->Systematics().Set(fNIWGSysts[indices[i]], val);
}
#endif
#endif
}
void NIWGWeightEngine::Reconfigure(bool silent) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fNIWGRW->Reconfigure();
// Shout again
if (silent) StartTalking();
#endif
#endif
}
double NIWGWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __NEUT_ENABLED__
#ifdef __NIWG_ENABLED__
// Skip Non GENIE
if (evt->fType != kNEUT) return 1.0;
// Hush now
StopTalking();
niwg::rew::NIWGEvent* niwg_event = NIWGWeightEngine::GetNIWGEventLocal(evt->fNeutVect);
rw_weight *= fNIWGRW->CalcWeight(*niwg_event);
delete niwg_event;
// Speak Now
StartTalking();
#endif
#endif
// Return rw_weight
return rw_weight;
}
#ifdef __NEUT_ENABLED__
#ifdef __NIWG_ENABLED__
niwg::rew::NIWGEvent * NIWGWeightEngine::GetNIWGEventLocal(NeutVect* nvect)
{
niwg::rew::NIWGEvent * fDummyNIWGEvent = NULL;
fDummyNIWGEvent = new niwg::rew::NIWGEvent();
fDummyNIWGEvent->detid = 1; // MiniBooNE (apply CCQE LowE variations)
fDummyNIWGEvent->neutmode = nvect->Mode;
fDummyNIWGEvent->targetA = nvect->TargetA;
fDummyNIWGEvent->recenu_ccqe_sk = -1;
if (nvect->Ibound == 0) fDummyNIWGEvent->targetA = 1; //RT: identifies as H, rather than O16
// Fill initial particle stack
for (int ip = 0; ip < nvect->Npart(); ++ip) {
niwg::rew::NIWGPartStack fDummyPartStack;
fDummyPartStack.p = (nvect->PartInfo(ip)->fP) * 0.001; // Convert to GeV
fDummyPartStack.pdg = nvect->PartInfo(ip)->fPID;
fDummyPartStack.chase = nvect->PartInfo(ip)->fIsAlive;
fDummyPartStack.parent = nvect->ParentIdx(ip) - 1; // WARNING: this needs to be tested with a NeutRoot file
fDummyNIWGEvent->part_stack.push_back(fDummyPartStack);
}
fDummyNIWGEvent->CalcKinematics();
return fDummyNIWGEvent;
}
#endif
#endif
diff --git a/src/Reweight/NuWroWeightEngine.cxx b/src/Reweight/NuWroWeightEngine.cxx
index 149e23a..c9dcc95 100644
--- a/src/Reweight/NuWroWeightEngine.cxx
+++ b/src/Reweight/NuWroWeightEngine.cxx
@@ -1,136 +1,136 @@
#include "NuWroWeightEngine.h"
NuWroWeightEngine::NuWroWeightEngine(std::string name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up NuWro RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
// Create the engine
fNuwroRW = new nuwro::rew::NuwroReWeight();
// Get List of Veto Calcs (For Debugging)
std::string rw_engine_list =
- FitPar::Config().GetParS("FitWeight.fNuwroRW_veto");
+ FitPar::Config().GetParS("FitWeight_fNuwroRW_veto");
bool xsec_qel = rw_engine_list.find("nuwro_QEL") == std::string::npos;
bool xsec_flag = rw_engine_list.find("nuwro_FlagNorm") == std::string::npos;
bool xsec_res = rw_engine_list.find("nuwro_RES") == std::string::npos;
// Add the RW Calcs
if (xsec_qel)
fNuwroRW->AdoptWghtCalc("nuwro_QEL", new nuwro::rew::NuwroReWeight_QEL);
if (xsec_flag)
fNuwroRW->AdoptWghtCalc("nuwro_FlagNorm",
new nuwro::rew::NuwroReWeight_FlagNorm);
if (xsec_res)
fNuwroRW->AdoptWghtCalc("nuwro_RES", new nuwro::rew::NuwroReWeight_SPP);
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
ERR(FTL) << "NUWRO RW NOT ENABLED! " << std::endl;
#endif
};
void NuWroWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Get RW Enum and name
int nuisenum = Reweight::ConvDial(name, kNUWRO);
// Initialise new vector
fEnumIndex[nuisenum];
fNameIndex[name];
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get Syst
nuwro::rew::NuwroSyst_t gensyst = nuwro::rew::NuwroSyst::FromString(name);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNUWROSysts.push_back(gensyst);
// Initialise Dial
LOG(FIT) << "Adding NuWro Syst " << fNUWROSysts[index] << std::endl;
fNuwroRW->Systematics().Add(fNUWROSysts[index]);
if (fIsAbsTwk) {
nuwro::rew::NuwroSystUncertainty::Instance()->SetUncertainty(
fNUWROSysts[index], 1.0, 1.0);
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(name, startval);
}
#endif
};
void NuWroWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fNuwroRW->Systematics().SetSystVal(fNUWROSysts[indices[i]], val);
}
#endif
}
void NuWroWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fNuwroRW->Systematics().SetSystVal(fNUWROSysts[indices[i]], val);
}
#endif
}
void NuWroWeightEngine::Reconfigure(bool silent) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fNuwroRW->Reconfigure();
// Shout again
if (silent) StartTalking();
#endif
}
double NuWroWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Skip Non GENIE
if (evt->fType != kNUWRO) return 1.0;
#ifdef __USE_NUWRO_SRW_EVENTS__
rw_weight = fNuwroRW->CalcWeight(evt->fNuwroSRWEvent, *evt->fNuwroParams);
#else
// Call Weight calculation
rw_weight = fNuwroRW->CalcWeight(evt->fNuwroEvent);
#endif
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx
index ee3e995..24cd30d 100755
--- a/src/Routines/MinimizerRoutines.cxx
+++ b/src/Routines/MinimizerRoutines.cxx
@@ -1,1506 +1,1506 @@
// 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;
fCorFree = NULL;
fDecomp = NULL;
fDecFree = NULL;
fStrategy = "Migrad,FixAtLimBreak,Migrad";
fRoutines.clear();
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = 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,MCMC");
};
//*************************************
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;
} 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.Set("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadSettings( fCompKey.GetS("cardfile"), "");
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")) {
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::GetParI("VERBOSITY");
errorcount += Config::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 parhigh = parnom + 1;
double parstep = 1;
// Override state if none given
if (!key.Has("state")) {
key.SetS("state", "FIX");
}
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("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;
} else {
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);
} 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);
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);
;
fStartVals[parname] = parnom;
fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = parlow;
fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
}
// Setup Samples ----------------------------------------------
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";
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
<< "\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;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
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");
// 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));
}
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());
fSampleFCN->CreateIterationTree("fit_iterations", FitBase::GetRW());
return;
}
//******************************************
void MinimizerRoutines::SetupFitter(std::string routine) {
//******************************************
// Make the fitter
std::string fitclass = "";
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";
} else if (!routine.compare("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 = "";
} else if (!routine.compare("MCMC")) {
UseMCMC = true;
}
// make minimizer
if (fMinimizer) delete fMinimizer;
if (UseMCMC) {
fMinimizer = new Simple_MH_Sampler();
} else {
fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype);
}
fMinimizer->SetMaxFunctionCalls(
- FitPar::Config().GetParI("minimizer.maxcalls"));
+ FitPar::Config().GetParI("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->SetTolerance(FitPar::Config().GetParD("minimizer.tolerance"));
- fMinimizer->SetStrategy(FitPar::Config().GetParI("minimizer.strategy"));
+ FitPar::Config().GetParI("MAXITERATIONS"));
+ fMinimizer->SetTolerance(FitPar::Config().GetParD("TOLERANCE"));
+ fMinimizer->SetStrategy(FitPar::Config().GetParI("STRATEGY"));
fMinimizer->SetFunction(*fCallFunctor);
int ipar = 0;
// 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);
// 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;
}
ipar++;
}
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
// "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
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) {
//*************************************
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;
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 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") 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;
// 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 curunits = "(sig.)";
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));
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;
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::endl;
}
//*************************************
void MinimizerRoutines::GetMinimizerState() {
//*************************************
LOG(FIT) << "Minimizer State: " << std::endl;
// Get X and Err
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];
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));
freey++;
}
freex++;
}
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 lowstatsevents = FitPar::Config().GetParI("LOWSTATEVENTS");
+ int maxevents = FitPar::Config().GetParI("MAXEVENTS");
int verbosity = FitPar::Config().GetParI("VERBOSITY");
std::string trueroutine = routine;
std::string substring = "LowStat";
trueroutine.erase(trueroutine.find(substring), substring.length());
// Set MAX EVENTS=1000
- Config::SetPar("input.maxevents", lowstatsevents);
+ Config::SetPar("MAXEVENTS", lowstatsevents);
Config::SetPar("VERBOSITY", 3);
SetupFCN();
RunFitRoutine(trueroutine);
- Config::SetPar("input.maxevents", maxevents);
+ Config::SetPar("MAXEVENTS", maxevents);
SetupFCN();
Config::SetPar("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",
FitBase::GetRW());
double scanmiddlepoint = fCurVals[fParams[i]];
// Determine N points needed
double limlow = fMinVals[fParams[i]];
double limhigh = fMaxVals[fParams[i]];
double step = fStepVals[fParams[i]];
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);
// 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);
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());
double scanmid_i = fCurVals[fParams[i]];
double scanmid_j = fCurVals[fParams[j]];
double limlow_i = fMinVals[fParams[i]];
double limhigh_i = fMaxVals[fParams[i]];
double step_i = fStepVals[fParams[i]];
double limlow_j = fMinVals[fParams[j]];
double limhigh_j = fMaxVals[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;
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;
// 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);
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;
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;
}
/*
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;
// 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);
valVect.push_back(fCurVals.at(name));
errVect.push_back(fErrorVals.at(name));
minVect.push_back(fMinVals.at(name));
maxVect.push_back(fMaxVals.at(name));
startVect.push_back(fStartVals.at(name));
endfixVect.push_back(fFixVals.at(name));
startfixVect.push_back(fStartFixVals.at(name));
ipar++;
}
int NFREE = fMinimizer->NFree();
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 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 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.GetXaxis()->SetBinLabel(i + 1, name.c_str());
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.GetXaxis()->SetBinLabel(freecount, name.c_str());
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"};
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());
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();
} 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,
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;
}
}
}
// 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");
errorDIR->cd();
// Make a second file to store throws
std::string tempFileName = fOutputFile;
if (tempFileName.find(".root") != std::string::npos)
tempFileName.erase(tempFileName.find(".root"), 5);
tempFileName += ".throws.root";
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");
nominal->cd();
fSampleFCN->Write();
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("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);
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
TIter next(nominal->GetListOfKeys());
TKey* key;
while ((key = (TKey*)next())) {
TClass* cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
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");
// 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());
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));
}
for (Int_t i = 0; i < nthrows; 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];
}
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;
}
return;
};
void MinimizerRoutines::ThrowDataToys() {
LOG(FIT) << "Generating Toy Data Throws" << std::endl;
int verb = Config::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++) {
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++) {
likes->Fill(values[i]);
}
// Save to file
LOG(FIT) << "Writing toy data throws" << std::endl;
fOutputRootFile->cd();
likes->Write();
}
diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx
index 2648470..1af6e2e 100644
--- a/src/Statistical/StatUtils.cxx
+++ b/src/Statistical/StatUtils.cxx
@@ -1,1303 +1,1303 @@
// 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 "TH1D.h"
#include "StatUtils.h"
#include "NuisConfig.h"
#include "GeneralUtils.h"
//*******************************************************************
Double_t StatUtils::GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask) {
//*******************************************************************
Double_t Chi2 = 0.0;
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// Add MC Error to data if required
- if (FitPar::Config().GetParB("statutils.addmcerror")) {
+ if (FitPar::Config().GetParB("addmcerror")) {
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double dterr = calc_data->GetBinError(i + 1);
double mcerr = calc_mc->GetBinError(i + 1);
if (dterr > 0.0) {
calc_data->SetBinError(i + 1, sqrt(dterr * dterr + mcerr * mcerr));
}
}
}
// Apply masking if required
if (mask) {
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Iterate over bins in X
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
// Ignore bins with zero data or zero bin error
if (calc_data->GetBinError(i + 1) <= 0.0 ||
calc_data->GetBinContent(i + 1) == 0.0) continue;
// Take mc data difference
double diff = calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1);
double err = calc_data->GetBinError(i + 1);
Chi2 += (diff * diff) / (err * err);
}
// cleanup
delete calc_data;
delete calc_mc;
return Chi2;
};
//*******************************************************************
Double_t StatUtils::GetChi2FromDiag(TH2D* data, TH2D* mc,
TH2I* map, TH2I* mask) {
//*******************************************************************
// Generate a simple map
if (!map) map = GenerateMap(data);
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate 1D chi2 from 1D Plots
Double_t Chi2 = StatUtils:: GetChi2FromDiag(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
};
//*******************************************************************
Double_t StatUtils::GetChi2FromCov(TH1D* data, TH1D* mc,
TMatrixDSym* invcov, TH1I* mask,
double data_scale, double covar_scale) {
//*******************************************************************
Double_t Chi2 = 0.0;
TMatrixDSym* calc_cov = (TMatrixDSym*) invcov->Clone();
TH1D* calc_data = (TH1D*) data->Clone();
TH1D* calc_mc = (TH1D*) mc->Clone();
// If a mask if applied we need to apply it before the matrix is inverted
if (mask) {
calc_cov = ApplyInvertedMatrixMasking(invcov, mask);
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Add MC Error to data if required
- if (FitPar::Config().GetParB("statutils.addmcerror")) {
+ if (FitPar::Config().GetParB("addmcerror")) {
// Make temp cov
TMatrixDSym* newcov = StatUtils::GetInvert(calc_cov);
// Add MC err to diag
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double mcerr = calc_mc->GetBinError(i + 1) * sqrt(covar_scale);
double oldval = (*newcov)(i, i);
std::cout << "Adding cov stat " << mcerr*mcerr << " to " << (*newcov)(i,i) << std::endl;
(*newcov)(i, i) = oldval + mcerr * mcerr;
}
// Reset the calc_cov to new invert
delete calc_cov;
calc_cov = GetInvert(newcov);
// Delete the tempcov
delete newcov;
}
calc_data->Scale(data_scale);
calc_mc ->Scale(data_scale);
(*calc_cov) *= covar_scale;
// iterate over bins in X (i,j)
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
for (int j = 0; j < calc_data->GetNbinsX(); j++) {
if (calc_data->GetBinContent(i + 1) != 0 || calc_mc->GetBinContent(i + 1) != 0) {
LOG(DEB) << "i j = " << i << " " << j << std::endl;
LOG(DEB) << "Calc_data mc i = " << calc_data->GetBinContent(i + 1)
<< " " << calc_mc->GetBinContent(i + 1)
<< " Dif = "
<< ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) )
<< std::endl;
LOG(DEB) << "Calc_data mc i = " << calc_data->GetBinContent(j + 1)
<< " " << calc_mc->GetBinContent(j + 1)
<< " Dif = "
<< ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1) )
<< std::endl;
LOG(DEB) << "Covar = " << (*calc_cov)(i, j) << std::endl;
LOG(DEB) << "Cont chi2 = " \
<< ( ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) \
* (*calc_cov)(i, j) \
* ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)))
<< " " << Chi2 << std::endl;
Chi2 += ( ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) \
* (*calc_cov)(i, j) \
* ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1) ) );
} else {
LOG(DEB) << "Error on bin (i,j) = (" << i << "," << j << ")" << std::endl;
LOG(DEB) << "data->GetBinContent(i+1) = " << calc_data->GetBinContent(i + 1) << std::endl;
LOG(DEB) << "mc->GetBinContent(i+1) = " << calc_mc->GetBinContent(i + 1) << std::endl;
LOG(DEB) << "Adding zero to chi2 instead of dying horrifically " << std::endl;
Chi2 += 0.;
}
}
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromCov( TH2D* data, TH2D* mc,
TMatrixDSym* invcov, TH2I* map, TH2I* mask) {
//*******************************************************************
// Generate a simple map
if (!map) {
map = StatUtils::GenerateMap(data);
}
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate 1D chi2 from 1D Plots
Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD( TH1D* data, TH1D* mc,
TMatrixDSym* cov, TH1I* mask) {
//*******************************************************************
Double_t Chi2 = 0.0;
TMatrixDSym* calc_cov = (TMatrixDSym*) cov->Clone();
TH1D* calc_data = (TH1D*) data->Clone();
TH1D* calc_mc = (TH1D*) mc->Clone();
// If a mask if applied we need to apply it before the matrix is inverted
if (mask) {
calc_cov = StatUtils::ApplyMatrixMasking(cov, mask);
calc_data = StatUtils::ApplyHistogramMasking(data, mask);
calc_mc = StatUtils::ApplyHistogramMasking(mc, mask);
}
// Decompose matrix
TDecompSVD LU = TDecompSVD((*calc_cov));
LU.Decompose();
TMatrixDSym* cov_U = new TMatrixDSym(calc_data->GetNbinsX(), LU .GetU().GetMatrixArray(), "");
TVectorD* cov_S = new TVectorD( LU.GetSig() );
// Apply basis rotation before adding up chi2
Double_t rotated_difference = 0.0;
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
rotated_difference = 0.0;
// Rotate basis of Data - MC
for (int j = 0; j < calc_data->GetNbinsY(); j++)
rotated_difference += ( calc_data->GetBinContent(j + 1) - calc_mc ->GetBinContent(j + 1) ) * (*cov_U)(j, i) ;
// Divide by rotated error cov_S
Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i);
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
delete cov_U;
delete cov_S;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD( TH2D* data, TH2D* mc,
TMatrixDSym* cov, TH2I* map, TH2I* mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
double StatUtils::GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask) {
//*******************************************************************
// If just an event rate, for chi2 just use Poission Likelihood to calculate the chi2 component
double chi2 = 0.0;
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// Apply masking if required
if (mask) {
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Iterate over bins in X
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double dt = calc_data->GetBinContent(i + 1);
double mc = calc_mc->GetBinContent(i + 1);
if (mc <= 0) continue;
if (dt <= 0) {
// Only add difference
chi2 += 2 * (mc - dt);
} else {
// Do the chi2 for Poisson distributions
chi2 += 2 * (mc - dt + (dt * log(dt / mc)));
}
/*
LOG(REC)<<"Evt Chi2 cont = "<Clone();
// If a mask is provided we need to apply it before getting NDOF
if (mask) {
calc_hist = StatUtils::ApplyHistogramMasking(hist, mask);
}
// NDOF is defined as total number of bins with non-zero errors
Int_t NDOF = 0;
for (int i = 0; i < calc_hist->GetNbinsX(); i++) {
if (calc_hist->GetBinError(i + 1) > 0.0) NDOF++;
}
delete calc_hist;
return NDOF;
};
//*******************************************************************
Int_t StatUtils::GetNDOF(TH2D* hist, TH2I* map, TH2I* mask) {
//*******************************************************************
Int_t NDOF = 0;
if (!map) map = StatUtils::GenerateMap(hist);
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1)) continue;
if (map->GetBinContent(i + 1, j + 1) <= 0) continue;
NDOF++;
}
}
return NDOF;
};
//*******************************************************************
TH1D* StatUtils::ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag, TH1I* mask) {
//*******************************************************************
TH1D* calc_hist = (TH1D*) hist->Clone( (std::string(hist->GetName()) + "_THROW" ).c_str() );
TMatrixDSym* calc_cov = (TMatrixDSym*) cov->Clone();
Double_t correl_val = 0.0;
// If a mask if applied we need to apply it before the matrix is decomposed
if (mask) {
calc_cov = ApplyMatrixMasking(cov, mask);
calc_hist = ApplyHistogramMasking(calc_hist, mask);
}
// If a covariance is provided we need a preset random vector and a decomp
std::vector rand_val;
TMatrixDSym* decomp_cov;
if (cov) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
rand_val.push_back(gRandom->Gaus(0.0, 1.0));
}
// Decomp the matrix
decomp_cov = StatUtils::GetDecomp(calc_cov);
}
// iterate over bins
for (int i = 0; i < hist->GetNbinsX(); i++) {
// By Default the errors on the histogram are thrown uncorrelated to the other errors
// if (throwdiag) {
// calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + \
// gRandom->Gaus(0.0, 1.0) * calc_hist->GetBinError(i + 1)) );
// }
// If a covariance is provided that is also thrown
if (cov) {
correl_val = 0.0;
for (int j = 0; j < hist->GetNbinsX(); j++) {
correl_val += rand_val[j] * (*decomp_cov)(j, i) ;
}
calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38));
}
}
delete calc_cov;
delete decomp_cov;
// return this new thrown data
return calc_hist;
};
//*******************************************************************
TH2D* StatUtils::ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map, bool throwdiag, TH2I* mask) {
//*******************************************************************
// PLACEHOLDER!!!!!!!!!
// Currently no support for throwing 2D Histograms from a covariance
(void) hist;
(void) cov;
(void) map;
(void) throwdiag;
(void) mask;
// /todo
// Sort maps if required
// Throw the covariance for a 1D plot
// Unmap back to 2D Histogram
return hist;
}
//*******************************************************************
TH1D* StatUtils::ApplyHistogramMasking(TH1D* hist, TH1I* mask) {
//*******************************************************************
if (!mask) return ( (TH1D*)hist->Clone() );
// This masking is only sufficient for chi2 calculations, and will have dodgy bin edges.
// Get New Bin Count
Int_t NBins = 0;
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1)) continue;
NBins++;
}
// Make new hist
std::string newmaskname = std::string(hist->GetName()) + "_MSKD";
TH1D* calc_hist = new TH1D( newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins);
// fill new hist
int binindex = 0;
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1)) {
LOG(REC) << "Applying mask to bin " << i + 1 << " " << hist->GetName() << std::endl;
continue;
}
calc_hist->SetBinContent(binindex + 1, hist->GetBinContent(i + 1));
calc_hist->SetBinError(binindex + 1, hist->GetBinError(i + 1));
binindex++;
}
return calc_hist;
};
//*******************************************************************
TH2D* StatUtils::ApplyHistogramMasking(TH2D* hist, TH2I* mask) {
//*******************************************************************
TH2D* newhist = (TH2D*) hist->Clone();
if (!mask) return newhist;
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1) > 0) {
newhist->SetBinContent(i + 1, j + 1, 0.0);
newhist->SetBinContent(i + 1, j + 1, 0.0);
}
}
}
return newhist;
}
//*******************************************************************
TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH1I* mask) {
//*******************************************************************
if (!mask) return (TMatrixDSym*)(mat->Clone());
// Get New Bin Count
Int_t NBins = 0;
for (int i = 0; i < mask->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1)) continue;
NBins++;
}
// make new matrix
TMatrixDSym* calc_mat = new TMatrixDSym(NBins);
int col, row;
// Need to mask out bins in the current matrix
row = 0;
for (int i = 0; i < mask->GetNbinsX(); i++) {
col = 0;
// skip if masked
if (mask->GetBinContent(i + 1) > 0.5) continue;
for (int j = 0; j < mask->GetNbinsX(); j++) {
// skip if masked
if (mask->GetBinContent(j + 1) > 0.5) continue;
(*calc_mat)(row, col) = (*mat)(i, j);
col++;
}
row++;
}
return calc_mat;
};
//*******************************************************************
TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1I* mask_1D = StatUtils::MapToMask(mask, map);
TMatrixDSym* newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D);
delete mask_1D;
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH1I* mask) {
//*******************************************************************
TMatrixDSym* new_mat = GetInvert(mat);
TMatrixDSym* masked_mat = ApplyMatrixMasking(new_mat, mask);
TMatrixDSym* inverted_mat = GetInvert(masked_mat);
delete masked_mat;
delete new_mat;
return inverted_mat;
};
//*******************************************************************
TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1I* mask_1D = StatUtils::MapToMask(mask, map);
TMatrixDSym* newmat = ApplyInvertedMatrixMasking(mat, mask_1D);
delete mask_1D;
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetInvert(TMatrixDSym* mat) {
//*******************************************************************
TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone();
// Check for diagonal
bool non_diagonal = false;
for (int i = 0; i < new_mat->GetNrows(); i++) {
for (int j = 0; j < new_mat->GetNrows(); j++) {
if (i == j) continue;
if ((*new_mat)(i, j) != 0.0) {
non_diagonal = true;
break;
}
}
}
// If diag, just flip the diag
if (!non_diagonal or new_mat->GetNrows() == 1) {
for (int i = 0; i < new_mat->GetNrows(); i++) {
if ((*new_mat)(i, i) != 0.0)
(*new_mat)(i, i) = 1.0 / (*new_mat)(i, i);
else
(*new_mat)(i, i) = 0.0;
}
return new_mat;
}
// Invert full matrix
TDecompSVD LU = TDecompSVD((*new_mat));
new_mat = new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), "");
return new_mat;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetDecomp(TMatrixDSym* mat) {
//*******************************************************************
TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone();
int nrows = new_mat->GetNrows();
// Check for diagonal
bool diagonal = true;
for (int i = 0; i < nrows; i++) {
for (int j = 0; j < nrows; j++) {
if (i == j) continue;
if ((*new_mat)(i, j) != 0.0) {
diagonal = false;
break;
}
}
}
// If diag, just flip the diag
if (diagonal or nrows == 1) {
for (int i = 0; i < nrows; i++) {
if ((*new_mat)(i, i) > 0.0)
(*new_mat)(i, i) = sqrt((*new_mat)(i, i));
else
(*new_mat)(i, i) = 0.0;
}
return new_mat;
}
TDecompChol LU = TDecompChol(*new_mat);
LU.Decompose();
delete new_mat;
TMatrixDSym* dec_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), "");
return dec_mat;
}
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH1D* hist, double norm) {
//*******************************************************************
if (!mat) mat = MakeDiagonalCovarMatrix(hist);
int nbins = mat->GetNrows();
TMatrixDSym* new_mat = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
double valx = hist->GetBinContent(i + 1) * 1E38;
double valy = hist->GetBinContent(j + 1) * 1E38;
(*new_mat)(i, j) = (*mat)(i, j) + norm * norm * valx * valy;
}
}
// Swap the two
delete mat;
mat = new_mat;
return;
};
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, TH2I* map ) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1D* data_1D = MapToTH1D(data, map);
StatUtils::ForceNormIntoCovar(mat, data_1D, norm);
delete data_1D;
return;
}
//*******************************************************************
TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH1D* data, double scaleF) {
//*******************************************************************
TMatrixDSym* newmat = new TMatrixDSym(data->GetNbinsX());
for (int i = 0; i < data->GetNbinsX(); i++) {
(*newmat)(i, i) = data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF;
}
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH2D* data, TH2I* map, double scaleF) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1D* data_1D = MapToTH1D(data, map);
return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF);
};
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, double scale) {
//*******************************************************************
// Check
if (cov->GetNrows() != data->GetNbinsX()) {
ERR(WRN) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << std::endl;
}
// Set bin errors form cov diag
for (int i = 0; i < data->GetNbinsX(); i++) {
data->SetBinError(i + 1, sqrt((*cov)(i, i)) * scale );
}
return;
}
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map, double scale) {
//*******************************************************************
// Create map if required
if (!map) map = StatUtils::GenerateMap(data);
// Set Bin Errors from cov diag
int count = 0;
for (int i = 0; i < data->GetNbinsX(); i++) {
for (int j = 0; j < data->GetNbinsY(); j++) {
if (data->GetBinContent(i + 1, j + 1) == 0.0) continue;
count = map->GetBinContent(i + 1, j + 1) - 1;
data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale );
}
}
return;
}
TMatrixDSym* StatUtils::ExtractShapeOnlyCovar(TMatrixDSym* full_covar, TH1* data_hist, double data_scale){
int nbins = full_covar->GetNrows();
TMatrixDSym* shape_covar = new TMatrixDSym(nbins);
// Check nobody is being silly
if (data_hist->GetNbinsX() != nbins){
ERR(WRN) << "Inconsistent matrix and data histogram passed to StatUtils::ExtractShapeOnlyCovar!" << std::endl;
ERR(WRN) << "data_hist has " << data_hist->GetNbinsX() << " matrix has " << nbins << std::endl;
int err_bins = data_hist->GetNbinsX();
if (nbins > err_bins) err_bins = nbins;
for (int i = 0; i < err_bins; ++i){
ERR(WRN) << "Matrix diag. = " << (*full_covar)(i, i) << " data = " << data_hist->GetBinContent(i+1) << std::endl;
}
return NULL;
}
double total_data = 0;
double total_covar = 0;
// Initial loop to calculate some constants
for (int i = 0; i < nbins; ++i) {
total_data += data_hist->GetBinContent(i+1)*data_scale;
for (int j = 0; j < nbins; ++j) {
total_covar += (*full_covar)(i,j);
}
}
if (total_data == 0 || total_covar == 0){
ERR(WRN) << "Stupid matrix or data histogram passed to StatUtils::ExtractShapeOnlyCovar! Ignoring..." << std::endl;
return NULL;
}
LOG(SAM) << "Norm error = " << sqrt(total_covar)/total_data << std::endl;
// Now loop over and calculate the shape-only matrix
for (int i = 0; i < nbins; ++i) {
double data_i = data_hist->GetBinContent(i+1)*data_scale;
for (int j = 0; j < nbins; ++j) {
double data_j = data_hist->GetBinContent(j+1)*data_scale;
double norm_term = data_i*data_j*total_covar/total_data/total_data;
double mix_sum1 = 0;
double mix_sum2 = 0;
for (int k = 0; k < nbins; ++k){
mix_sum1 += (*full_covar)(k,j);
mix_sum2 += (*full_covar)(i,k);
}
double mix_term1 = data_i*(mix_sum1/total_data - total_covar*data_j/total_data/total_data);
double mix_term2 = data_j*(mix_sum2/total_data - total_covar*data_i/total_data/total_data);
(*shape_covar)(i, j) = (*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term;
}
}
return shape_covar;
}
//*******************************************************************
TH2I* StatUtils::GenerateMap(TH2D* hist) {
//*******************************************************************
std::string maptitle = std::string(hist->GetName()) + "_MAP";
TH2I* map = new TH2I( maptitle.c_str(), maptitle.c_str(),
hist->GetNbinsX(), 0, hist->GetNbinsX(),
hist->GetNbinsY(), 0, hist->GetNbinsY());
Int_t index = 1;
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (hist->GetBinContent(i + 1, j + 1) > 0 && hist->GetBinError(i + 1, j + 1) > 0) {
map->SetBinContent(i + 1, j + 1, index);
index++;
} else {
map->SetBinContent(i + 1, j + 1, 0);
}
}
}
return map;
}
//*******************************************************************
TH1D* StatUtils::MapToTH1D(TH2D* hist, TH2I* map) {
//*******************************************************************
if (!hist) return NULL;
// Get N bins for 1D plot
Int_t Nbins = map->GetMaximum();
std::string name1D = std::string(hist->GetName()) + "_1D";
// Make new 1D Hist
TH1D* newhist = new TH1D(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);
// map bin contents
for (int i = 0; i < map->GetNbinsX(); i++) {
for (int j = 0; j < map->GetNbinsY(); j++) {
if (map->GetBinContent(i + 1, j + 1) == 0) continue;
newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1));
newhist->SetBinError(map->GetBinContent(i + 1, j + 1), hist->GetBinError(i + 1, j + 1));
}
}
// return
return newhist;
}
//*******************************************************************
TH1I* StatUtils::MapToMask(TH2I* hist, TH2I* map) {
//*******************************************************************
TH1I* newhist = NULL;
if (!hist) return newhist;
// Get N bins for 1D plot
Int_t Nbins = map->GetMaximum();
std::string name1D = std::string(hist->GetName()) + "_1D";
// Make new 1D Hist
newhist = new TH1I(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);
// map bin contents
for (int i = 0; i < map->GetNbinsX(); i++) {
for (int j = 0; j < map->GetNbinsY(); j++) {
if (map->GetBinContent(i + 1, j + 1) == 0) continue;
newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1));
}
}
// return
return newhist;
}
TMatrixDSym* StatUtils::GetCovarFromCorrel(TMatrixDSym* correl, TH1D* data) {
int nbins = correl->GetNrows();
TMatrixDSym* covar = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*covar)(i, j) = (*correl)(i, j) * data->GetBinError(i + 1) * data->GetBinError(j + 1);
}
}
return covar;
}
//*******************************************************************
TMatrixD* StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx, int dimy) {
//*******************************************************************
// Determine dim
if (dimx == -1 and dimy == -1) {
std::string line;
std::ifstream covar(covfile.c_str(), std::ifstream::in);
int row = 0;
while (std::getline(covar >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
"entries on this line: " << row << std::endl;
}
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
column++;
if (column > dimx) dimx = column;
}
row++;
if (row > dimy) dimy = row;
}
}
// Or assume symmetric
if (dimx != -1 and dimy == -1) {
dimy = dimx;
}
assert(dimy != -1 && " matrix dimy not set.");
// Make new matrix
TMatrixD* mat = new TMatrixD(dimx, dimy);
std::string line;
std::ifstream covar(covfile.c_str(), std::ifstream::in);
int row = 0;
while (std::getline(covar >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
"entries on this line: " << row << std::endl;
}
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
// Check Rows
//assert(row > mat->GetNrows() && " covar rows doesn't match matrix rows.");
//assert(column > mat->GetNcols() && " covar cols doesn't match matrix cols.");
// Fill Matrix
(*mat)(row, column) = (*iter);
column++;
}
row++;
}
return mat;
}
//*******************************************************************
TMatrixD* StatUtils::GetMatrixFromRootFile(std::string covfile, std::string histname) {
//*******************************************************************
std::string inputfile = covfile + ";" + histname;
std::vector splitfile = GeneralUtils::ParseToStr(inputfile, ";");
if (splitfile.size() < 2) {
ERR(FTL) << "No object name given!" << std::endl;
throw;
}
// Get file
TFile* tempfile = new TFile(splitfile[0].c_str(), "READ");
// Get Object
TObject* obj = tempfile->Get(splitfile[1].c_str());
if (!obj) {
ERR(FTL) << "Object " << splitfile[1] << " doesn't exist!" << std::endl;
throw;
}
// Try casting
TMatrixD* mat = dynamic_cast(obj);
if (mat) {
TMatrixD* newmat = (TMatrixD*)mat->Clone();
delete mat;
tempfile->Close();
return newmat;
}
TMatrixDSym* matsym = dynamic_cast(obj);
if (matsym) {
TMatrixD* newmat = new TMatrixD(matsym->GetNrows(), matsym->GetNrows());
for (int i = 0; i < matsym->GetNrows(); i++) {
for (int j = 0; j < matsym->GetNrows(); j++) {
(*newmat)(i, j) = (*matsym)(i, j);
}
}
delete matsym;
tempfile->Close();
return newmat;
}
TH2D* mathist = dynamic_cast(obj);
if (mathist) {
TMatrixD* newmat = new TMatrixD(mathist->GetNbinsX(), mathist->GetNbinsX());
for (int i = 0; i < mathist->GetNbinsX(); i++) {
for (int j = 0; j < mathist->GetNbinsX(); j++) {
(*newmat)(i, j) = mathist->GetBinContent(i + 1, j + 1);
}
}
delete mathist;
tempfile->Close();
return newmat;
}
return NULL;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetCovarFromTextFile(std::string covfile, int dim){
//*******************************************************************
// Delete TempMat
TMatrixD* tempmat = GetMatrixFromTextFile(covfile, dim, dim);
// Make a symmetric covariance
TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows());
for (int i = 0; i < tempmat->GetNrows(); i++){
for (int j = 0; j < tempmat->GetNrows(); j++){
(*newmat)(i,j) = (*tempmat)(i,j);
}
}
delete tempmat;
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetCovarFromRootFile(std::string covfile, std::string histname){
//*******************************************************************
TMatrixD* tempmat = GetMatrixFromRootFile(covfile, histname);
TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows());
for (int i = 0; i < tempmat->GetNrows(); i++){
for (int j = 0; j < tempmat->GetNrows(); j++){
(*newmat)(i,j) = (*tempmat)(i,j);
}
}
delete tempmat;
return newmat;
}