Page MenuHomeHEPForge
Files F9501775

No OneTemporary
Actions

View Options

diff --git a/parameters/config.xml b/parameters/config.xml
index 11031cb..dc38035 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,157 +1,157 @@
<nuisance>
<!-- # ###################################################### -->
<!-- # NUISANCE CONFIGURATION OPTIONS -->
<!-- # This file is read in by default at runtime -->
<!-- # If you want to override on a case by case bases use -q at runtime -->
<!-- # ###################################################### -->
<!-- # MAIN Configs -->
<!-- # ###################################################### -->
<!-- # Logger goes from -->
<!-- # 1 Quiet -->
<!-- # 2 Fitter -->
<!-- # 3 Samples -->
<!-- # 4 Reconfigure Loops -->
<!-- # 5 Every Event print out (SHOUT) -->
<!-- # -1 DEBUGGING -->
<config verbosity='5'/>
<config VERBOSITY='5'/>
<!-- # ERROR goes from -->
<!-- # 0 NONE -->
<!-- # 1 FATAL -->
<!-- # 2 WARN -->
<config ERROR='2'/>
<config cores='2' />
<!-- # Input Configs -->
<!-- # ###################################################### -->
<!-- # Default Requirements file for the externalDataFitter Package -->
<!-- # MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions. -->
<config input.maxevents='-1'/>
<config MAXEVENTS='-1'/>
<config input.MAXEVENTS='-1'/>
<config includeemptystackhists='0'/>
<!-- # Turn on/off event manager -->
<!-- # EventManager enables us to only loop number of events once for multiple projections of the same measurements -->
<!-- # e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one -->
<config input.eventmanager='1'/>
<config EventManager='1'/>
<!-- # Event Directories -->
<!-- # Can setup default directories and use @EVENT_DIR/path to link to it -->
<config EVENT_DIR='/data2/stowell/NIWG/'/>
<config NEUT_EVENT_DIR='/data2/stowell/NIWG/neut/fit_samples_neut5.3.3/'/>
<config GENIE_EVENT_DIR='/data2/stowell/NIWG/genie/fit_samples_R.2.10.0/'/>
<config NUWRO_EVENT_DIR='/data2/stowell/NIWG/nuwro/fit_samples/'/>
<config GIBUU_EVENT_DIR='/data/GIBUU/DIR/'/>
<!-- # In PrepareGENIE the reconstructed splines can be saved into the file -->
<config save_genie_splines='1'/>
<!-- # In InputHandler the option to regenerate NuWro flux/xsec plots is available -->
<!-- # Going to move this to its own app soon -->
<config input.regen_nuwro_plots='0'/>
<!-- # DEVEL CONFIG OPTION, don't touch! -->
<config cachesize='0'/>
<!-- # ReWeighting Configuration Options -->
<!-- # ###################################################### -->
<!-- # Set absolute twkdial for parameters -->
<config params.setabstwk='0'/>
<!-- # Convert Dials in output statements using dial conversion card -->
<config convert_dials='0'/>
<!-- # Make RW Calculations be quiet -->
<config params.silentweighting='0'/>
<!-- # Vetos can be used to specify RW dials NOT to be loaded in -->
<!-- # Useful if one specific one has an issue -->
<config FitWeight.fNIWGRW_veto=''/>
<config FitWeight.fNuwroRW_veto=''/>
<config FitWeight.fNeutRW_veto=''/>
<config FitWeight.fGenieRW_veto=''/>
<!-- # Output Options -->
<!-- # ###################################################### -->
<!-- # Save Nominal prediction with all rw engines at default -->
<config savenominal='0'/>
<!-- # Save prefit with values at starting values -->
<config saveprefit='0'/>
<!-- # Here's the full list of drawing options -->
<!-- # See src/FitBase/Measurement1D::Write for more info -->
<!-- #config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP -->
<!-- #config drawopts DATA/MC -->
<config drawopts='DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC'/>
<!-- # Save the shape scaling applied with option SHAPE into the main MC hist -->
<config saveshapescaling='0'/>
<!-- # Set style of 1D output histograms -->
<config linecolour='1'/>
<config linestyle='1'/>
<config linewidth='1'/>
<!-- # For GenericFlux -->
<config isLiteMode='0'/>
<!-- # Statistical Options -->
<!-- # ###################################################### -->
<!-- # Add MC Statistical error to likelihoods -->
<config statutils.addmcerror='0'/>
<!-- # NUISMIN Configurations -->
<!-- # ###################################################### -->
<config minimizer.maxcalls='1000000'/>
<config minimizer.maxiterations='1000000'/>
<config minimizer.tolerance='0.001'/>
<!-- # Number of events required in low stats routines -->
<config minimizer.lowstatevents='25000'/>
<!-- # Error band generator configs -->
<!-- # ###################################################### -->
<!-- # For -f ErrorBands creates error bands for given measurements -->
<!-- # How many throws do we want (The higher the better precision) -->
-<config error_throws='250'/>
+<config error_throws='500'/>
<!-- # Are we throwing uniform or according to Gaussian? -->
<!-- # Only use uniform if wanting to study the limits of a dial. -->
<config error_uniform='0'/>
<!-- # Other Individual Case Configs -->
<!-- # ###################################################### -->
<!-- # Covariance throw options for fake data studies with MiniBooNE data. -->
<config thrown_covariance='FULL'/>
<config throw_mc_stat='0.0'/>
<config throw_diag_syst='0'/>
<config throw_corr_syst='0'/>
<config throw_mc_stat='0'/>
<!-- # Apply a shift to the muon momentum before calculation of Q2 -->
<config muon_momentum_shift='0.0'/>
<config muon_momentum_throw='0'/>
<!-- # MINERvA Specific Configs -->
<config MINERvA_XSec_CCinc_2DEavq3_nu.hadron_cut='0'/>
<config MINERvA_CCinc_XSec_2DEavq3_nu.useq3true='0'/>
<config Modes.split_PN_NN='0'/>
<config SignalReconfigures='0'/>
</nuisance>
diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
index 3e5dd8d..2dacd29 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,874 +1,913 @@
#include "JointFCN.h"
#include <stdio.h>
#include "FitUtils.h"
//***************************************************
JointFCN::JointFCN(std::string cardfile, TFile* outfile) {
//***************************************************
fOutputDir = gDirectory;
FitPar::Config().out = outfile;
fCardFile = cardfile;
LoadSamples(fCardFile);
fCurIter = 0;
fMCFilled = false;
fOutputDir->cd();
fIterationTree = NULL;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
};
JointFCN::JointFCN(TFile* outfile) {
fOutputDir = gDirectory;
if (outfile)
FitPar::Config().out = outfile;
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
LoadSamples(samplekeys);
fCurIter = 0;
fMCFilled = false;
// fOutputDir->cd();
fIterationTree = NULL;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
JointFCN::JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile) {
fOutputDir = gDirectory;
if (outfile)
FitPar::Config().out = outfile;
LoadSamples(samplekeys);
fCurIter = 0;
fMCFilled = false;
fOutputDir->cd();
fIterationTree = NULL;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
//***************************************************
JointFCN::~JointFCN() {
//***************************************************
// Delete Samples
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
delete exp;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
delete pull;
}
// Sort Tree
if (fIterationTree) DestroyIterationTree();
if (fDialVals) delete fDialVals;
if (fSampleLikes) delete fSampleLikes;
};
//***************************************************
void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) {
//***************************************************
LOG(FIT) << " Creating new iteration tree! " << endl;
if (fIterationTree && !name.compare(fIterationTree->GetName())) {
fIterationTree->Reset();
return;
}
fIterationTree = new TTree(name.c_str(), name.c_str());
// Setup Main Branches
fIterationTree->Branch("total_likelihood", &fLikelihood,
"total_likelihood/D");
fIterationTree->Branch("total_ndof", &fNDOF, "total_ndof/I");
// Setup Sample Arrays
int ninputs = fSamples.size() + fPulls.size();
fSampleLikes = new double[ninputs];
fSampleNDOF = new int[ninputs];
fNDOF = GetNDOF();
// Setup Sample Branches
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
std::string name = exp->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
std::string name = pull->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
// Add Dial Branches
std::vector<std::string> dials = rw->GetDialNames();
fNDials = dials.size();
fDialVals = new double[fNDials];
for (int i = 0; i < fNDials; i++) {
fIterationTree->Branch(dials[i].c_str(), &fDialVals[i],
(dials[i] + "/D").c_str());
}
}
//***************************************************
void JointFCN::DestroyIterationTree() {
//***************************************************
if (!fIterationTree) {
delete fIterationTree;
}
}
//***************************************************
void JointFCN::WriteIterationTree() {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Can't save empty iteration tree!" << endl;
throw;
}
fIterationTree->Write();
}
//***************************************************
void JointFCN::FillIterationTree(FitWeight* rw) {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << endl;
throw;
}
rw->GetAllDials(fDialVals, fNDials);
fIterationTree->Fill();
}
//***************************************************
double JointFCN::DoEval(const double* x) {
//***************************************************
// WEIGHT ENGINE
fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
FitBase::GetRW()->UpdateWeightEngine(x);
if (fDialChanged) {
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
}
if (LOG_LEVEL(REC)){
FitBase::GetRW()->Print();
}
// SORT SAMPLES
ReconfigureSamples();
// GET TEST STAT
fLikelihood = GetLikelihood();
// PRINT PROGRESS
LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood
<< std::endl;
// UPDATE TREE
if (fIterationTree) FillIterationTree(FitBase::GetRW());
return fLikelihood;
}
//***************************************************
int JointFCN::GetNDOF() {
//***************************************************
int totaldof = 0;
int count = 0;
// Total number of Free bins in each MC prediction
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
int dof = exp->GetNDOF();
// Save Seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double dof = pull->GetLikelihood();
// Save seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Set Data Variable
fNDOF = totaldof;
return totaldof;
}
//***************************************************
double JointFCN::GetLikelihood() {
//***************************************************
LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..." << " : " << "-2logL" << endl;
// Loop and add up likelihoods in an uncorrelated way
double like = 0.0;
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
double newlike = exp->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : " << newlike << endl;
// Add Weight Scaling
// like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName());
// Add to total
like += newlike;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double newlike = pull->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
// Add to total
like += newlike;
count++;
}
// Set Data Variable
fLikelihood = like;
return like;
};
void JointFCN::LoadSamples(std::vector<nuiskey> samplekeys) {
LOG(MIN) << "Loading Samples " << std::endl;
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys[i];
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype = key.GetS("type");
std::string fakeData = "";
fOutputDir->cd();
MeasurementBase* NewLoadedSample
= SampleUtils::CreateSample(key);
if (!NewLoadedSample) {
ERR(FTL) << "Could not load sample provided: " << samplename << std::endl;
ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
<< endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
}
//***************************************************
void JointFCN::LoadSamples(std::string cardinput)
//***************************************************
{
LOG(MIN) << "Initializing Samples" << std::endl;
// Read the card file here and load objects
std::string line;
std::ifstream card(cardinput.c_str(), ifstream::in);
// Make sure they are created in correct working DIR
fOutputDir->cd();
while (std::getline(card >> std::ws, line, '\n')) {
// Skip Empties
if (line.c_str()[0] == '#') continue;
if (line.empty()) continue;
// Parse line
std::vector<std::string> samplevect = GeneralUtils::ParseToStr(line, " ");
// Sample Inputs
if (!samplevect[0].compare("sample")) {
// Get all inputs
std::string name = samplevect[1];
std::string files = samplevect[2];
std::string type = "DEFAULT";
if (samplevect.size() > 3) type = samplevect[3];
// Create Sample Class
LOG(MIN) << "Loading up sample: " << name << " << " << files << " ("
<< type << ")" << std::endl;
std::string fakeData = "";
fOutputDir->cd();
MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(name, files, type,
fakeData, FitBase::GetRW());
if (!NewLoadedSample) {
ERR(FTL) << "Could not load sample provided: " << name << std::endl;
ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
<< endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
// Sample Inputs
if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") ||
!samplevect[0].compare("throws")) {
// Get all inputs
std::string name = samplevect[1];
std::string files = samplevect[2];
std::string type = "DEFAULT";
if (samplevect.size() > 3) {
type = samplevect[3];
} else if (!samplevect[0].compare("pull")) {
type = "GAUSPULL";
} else if (!samplevect[0].compare("throws")) {
type = "GAUSTHROWS";
}
// Create Pull Class
LOG(MIN) << "Loading up pull term: " << name << " << " << files << " ("
<< type << ")" << std::endl;
std::string fakeData = "";
fOutputDir->cd();
fPulls.push_back(new ParamPull(name, files, type));
}
}
card.close();
};
//***************************************************
void JointFCN::ReconfigureSamples(bool fullconfig) {
//***************************************************
int starttime = time(NULL);
LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << endl;
// Event Manager Reconf
if (fUsingEventManager) {
if (!fullconfig and fMCFilled)
ReconfigureFastUsingManager();
else
ReconfigureUsingManager();
} else {
// Loop over all Measurement Classes
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
// If RW Either do signal or full reconfigure.
if (fDialChanged or !fMCFilled or fullconfig) {
if (!fullconfig and fMCFilled)
exp->ReconfigureFast();
else
exp->Reconfigure();
// If RW Not needed just do normalisation
} else {
exp->Renormalise();
}
}
}
// Loop over pulls and update
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Reconfigure();
}
fMCFilled = true;
LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in "
<< time(NULL) - starttime << "s" << endl;
fCurIter++;
}
//***************************************************
void JointFCN::ReconfigureSignal() {
//***************************************************
this->ReconfigureSamples(false);
}
//***************************************************
void JointFCN::ReconfigureAllEvents() {
//***************************************************
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
this->ReconfigureSamples(true);
}
std::vector<InputHandlerBase*> JointFCN::GetInputList() {
std::vector<InputHandlerBase*> InputList;
+ fIsAllSplines = true;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
InputHandlerBase* inp = subsamples[i]->GetInput();
if (std::find(InputList.begin(), InputList.end(), inp) == InputList.end()) {
+
+ if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER) fIsAllSplines = false;
+
InputList.push_back(subsamples[i]->GetInput());
}
}
}
+
+
return InputList;
}
std::vector<MeasurementBase*> JointFCN::GetSubSampleList() {
std::vector<MeasurementBase*> SampleList;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
SampleList.push_back(subsamples[i]);
}
}
return SampleList;
}
//***************************************************
void JointFCN::ReconfigureUsingManager() {
//***************************************************
// New event manager plan
// ReconfigureEvtManager(){
// - Resets all our nice event vectors
// - Setflag: save variables
// - Gets list of inputs, and list of measurements, and list of input pointers for each measurement.
// - Loops over all inputs events
// - Create new SignalBoxVector
// - Compares event input pointer to current input pointer.
// - WeightCalc
// - if match: Run Event filling.
// - Do this using multithreading.
// - Call isSignal
// - Call FillEventVariables
// - GetPointerToBOX
// - Call FillHistogramsFromBox()
// - PushBack BOX->CloneSignalBox() into vector.
// - After looping over all measurements, if SignBoxVector not empty push it back into a vector and push bool back into another vector.
// - If pushing back signal, also push back double for event weight.
// - Call ScaleEvents, etc.
LOG(SAM) << "Reconfuguringusingeventmanager" << std::endl;
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// Check saving variables flag
bool savesignal = (FitPar::Config().GetParB("SignalReconfigures"));
if (savesignal) {
// Reset all of our event signal vectors
fSignalEventBoxes.clear();
fSignalEventFlags.clear();
fSampleSignalFlags.clear();
+ fSignalEventSplines.clear();
}
// Check out list of inputs
if (fInputList.empty()) {
fInputList = GetInputList();
fSubSampleList = GetSubSampleList();
}
+ std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
+ for (; inp_iter != fInputList.end(); inp_iter++) {
+ InputHandlerBase* curinput = (*inp_iter);
+ BaseFitEvt* curevent = curinput->FirstBaseEvent();
+ if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true);
+ }
+
int fillcount = 0;
// Loop over all inputs
int inputcount = 0;
- std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
+ inp_iter = fInputList.begin();
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
FitEvent* curevent = curinput->FirstNuisanceEvent();
int i = 0;
int nevents = curinput->GetNEvents();
int countwidth = nevents / 20;
// Start event loop
while (curevent) {
-
-
// Get Event Weight
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
double rwweight = curevent->Weight;
// Logging
if (LOG_LEVEL(REC)) {
if (i % countwidth == 0) {
LOG(REC) << "Processed " << i << " events. [M, W] = [" << curevent->Mode << ", " << rwweight << "]" << std::endl;
}
}
// std::cout << "Event " << i << " Weight = " << rwweight << std::endl;
// Setup flag for if signal found in at least one sample
bool foundsignal = false;
// Create a new signal bitset for this event
std::vector<bool> signalbitset(fSubSampleList.size());
// Create a new signal box vector for this event
std::vector<MeasurementVariableBox*> signalboxes;
// Loop over all samples
size_t measitercount = 0;
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
for (; meas_iter != fSubSampleList.end(); meas_iter++) {
// Get Measurement
MeasurementBase* curmeas = (*meas_iter);
// Compare input pointers, if != then skip
if (curinput != curmeas->GetInput()) {
if (savesignal) {
// Set bit to 0 as definitely not signal
signalbitset[measitercount] = 0;
}
measitercount++;
// Skip sample
continue;
}
// Fill events
MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);
bool signal = curmeas->isSignal(curevent);
// Get the event box after fill event variable
// std::cout << "Filling Meas Full = " << curmeas << std::endl;
curmeas->FillHistograms(rwweight);
if (signal) {
fillcount++;
}
if (savesignal) {
// Fill signal bitset
signalbitset[measitercount] = signal;
// If signal save a clone of the event box.
if (signal) {
foundsignal = true;
signalboxes.push_back( box->CloneSignalBox() );
}
}
measitercount++;
}
// push into vectors
if (savesignal) {
fSignalEventFlags.push_back(foundsignal);
if (foundsignal) {
fSignalEventBoxes.push_back(signalboxes);
fSampleSignalFlags.push_back(signalbitset);
- }
+
+ if (fIsAllSplines){
+ std::vector<float> coeff;
+ for (size_t l = 0; l < curevent->fSplineRead->GetNPar(); l++){
+ coeff.push_back( curevent->fSplineCoeff[l] );
+ }
+ fSignalEventSplines.push_back(coeff);
+ }
+ }
}
signalboxes.clear();
signalbitset.clear();
// iterate
curevent = curinput->NextNuisanceEvent();
i++;
}
inputcount++;
}
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
if (savesignal){
- double mem = ( //sizeof(fSignalEventBoxes) +
+ int mem = ( //sizeof(fSignalEventBoxes) +
// fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
sizeof(MeasurementVariableBox1D) * fillcount) * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " signal boxes for faster access. (~" << mem << " MB)" << std::endl;
+ if (fIsAllSplines and !fSignalEventSplines.empty()){
+ int splmem = sizeof(float)* fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6;
+ LOG(REC) << " -> Saved " << fillcount << " spline sets into memory. (~" << splmem << " MB)" << std::endl;
+ }
}
LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl;
};
//***************************************************
void JointFCN::ReconfigureFastUsingManager() {
//***************************************************
// ReconfigureFastEvtManager(){
// - Check for saved variables: Use normal if not.
// - Iterate over signalboolvector and look for true.
// - get tree entry and calculate new event weight.
// - Loop over measurements looking for ones that are actually signal.
// - MultiThread the FillHistogramsFromBox thing.
// - Create an array of what each box corresponds to.
// - If signal, get 'next' box and call FillHistogramsFromBox(weight)
// - Finish loop.
// - Call ScaleEvents, etc.
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// Check for saved variables
if (fSignalEventFlags.empty()) {
ReconfigureUsingManager();
return;
}
// Loop over all inputs
int curindex = 0;
// Get iterators
std::vector<bool>::iterator inpsig_iter = fSignalEventFlags.begin();
std::vector< std::vector<MeasurementVariableBox*> >::iterator box_iter = fSignalEventBoxes.begin();
+ std::vector< std::vector<float> >::iterator spline_iter = fSignalEventSplines.begin();
std::vector< std::vector<bool> >::iterator samsig_iter = fSampleSignalFlags.begin();
// Setup stuff for logging
int fillcount = 0;
int nevents = fSignalEventFlags.size();
int countwidth = nevents / 5;
- // Start input iterators
std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
+ BaseFitEvt* curevent = curinput->FirstBaseEvent();
+ if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true);
+ }
+
+
+ // Start input iterators
+ inp_iter = fInputList.begin();
+ for (; inp_iter != fInputList.end(); inp_iter++) {
+ InputHandlerBase* curinput = (*inp_iter);
// Get Events
BaseFitEvt* curevent = curinput->FirstBaseEvent();
int i = 0;
while (curevent != 0) {
// Logging
if (LOG_LEVEL(REC)){
if (i % countwidth == 0){
LOG(REC) << "Processed " << i << " signal events." << std::endl;
}
}
// If event is not signal skip it.
if (!(*inpsig_iter)) {
inpsig_iter++;
i++;
continue;
}
// Get Event
curevent = curinput->GetBaseEvent(i);
if (!curevent) break;
+ // Setup signal splines if required.
+ if (fIsAllSplines){
+ curevent->fSplineCoeff = &(*spline_iter)[0];
+ }
+
// Get Event Weight
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
double rwweight = curevent->Weight;
// Iterate over the measurements and get the corresponding signal boxes.
size_t measitercount = 0;
size_t boxitercount = 0;
// Get vectors for this event
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
std::vector<bool>::iterator subsamsig_iter = (*samsig_iter).begin();
std::vector<MeasurementVariableBox*>::iterator subbox_iter = (*box_iter).begin();
for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
// Get Measurement
MeasurementBase* curmeas = (*meas_iter);
// If not signal continue
if (*subsamsig_iter) {
// std::cout << "Filling histograms from box " << std::endl;
curmeas->SetSignal(true);
curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);
subbox_iter++;
fillcount++;
}
}
// Iterate over boxes
samsig_iter++;
box_iter++;
+ spline_iter++;
// iterate to next signal event
inpsig_iter++;
i++;
}
}
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
LOG(REC) << "Time taken ReconfigureFastUsingManager() : " << time(NULL) - timestart << std::endl;
}
//***************************************************
void JointFCN::Write() {
//***************************************************
// Loop over individual experiments and call Write
LOG(MIN) << "Writing each of the data classes..." << endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->Write();
}
// Save Pull Terms
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Write();
}
if (FitPar::Config().GetParB("EventManager")) {
// Get list of inputs
std::map<int, InputHandlerBase*> fInputs = FitBase::EvtManager().GetInputs();
std::map<int, InputHandlerBase*>::const_iterator iterInp;
for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) {
InputHandlerBase* input = (iterInp->second);
input->GetFluxHistogram()->Write();
input->GetXSecHistogram()->Write();
input->GetEventHistogram()->Write();
}
}
};
//***************************************************
void JointFCN::SetFakeData(std::string fakeinput) {
//***************************************************
LOG(MIN) << "Setting fake data from " << fakeinput << endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->SetFakeDataValues(fakeinput);
}
return;
}
diff --git a/src/FCN/JointFCN.h b/src/FCN/JointFCN.h
index 7c975d6..bf4aef1 100755
--- a/src/FCN/JointFCN.h
+++ b/src/FCN/JointFCN.h
@@ -1,169 +1,170 @@
#ifndef _JOINT_FCN_H_
#define _JOINT_FCN_H_
/*!
* \addtogroup FCN
* @{
*/
#include <iostream>
#include <vector>
#include <fstream>
#include <list>
// ROOT headers
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include <TMatrixDSym.h>
#include "TGraphErrors.h"
#include <TVectorD.h>
#include <TMath.h>
#include "FitUtils.h"
// External fitter headers
#include "MeasurementBase.h"
#include "SampleList.h"
#define GetCurrentDir getcwd
#include "EventManager.h"
#include "Math/Functor.h"
#include "ParamPull.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#include "MeasurementVariableBox.h"
#include "MeasurementVariableBox1D.h"
using namespace FitUtils;
using namespace FitBase;
//! Main FCN Class which ROOT's joint function needs to evaulate the chi2 at each stage of the fit.
class JointFCN
{
public:
//! Constructor
//! cardfile = Path to input card file listing samples
JointFCN(std::string cardfile, TFile *outfile=NULL);
JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile=NULL);
JointFCN(TFile* outfile=NULL); // Loads from global config
//! Destructor
~JointFCN();
//! Create sample list from cardfile
void LoadSamples(std::string cardFile);
void LoadSamples(std::vector<nuiskey> samplekeys);
//! Main Likelihood evaluation FCN
double DoEval(const double *x);
//! Func Wrapper for ROOT
inline double operator() (const std::vector<double> & x) {
double* x_array = new double[x.size()];
return this->DoEval(x_array);
};
//! Func Wrapper for ROOT
inline double operator() (const double *x) {
return this->DoEval(x);
};
//! Create a TTree to save all dial value iterations for this FCN
void CreateIterationTree(std::string name, FitWeight* rw);
//! Fills dial values and sample likelihoods into tree
void FillIterationTree(FitWeight* rw);
//! Writes TTree to fOutput directory
void WriteIterationTree();
//! Deletes TTree
void DestroyIterationTree();
//! Get Degrees of Freedom for samples (NBins)
int GetNDOF();
//! Return NDOF wrapper
inline unsigned int NDim() {return this->GetNDOF();};
//! Reconfigure samples where we force all events to be looped over.
void ReconfigureAllEvents() ;
//! Call Reconfigure on samples.
//! Choice of reconfigure type depends on whether dials have changed
//! and the MC has been filled.
void ReconfigureSamples(bool fullconfig = false);
//! Call reconfigure on only signal events (defaults to all events if CurIter=0)
void ReconfigureSignal();
//! Gets likelihood for all samples in FCN (assuming uncorrelated)
double GetLikelihood();
//! Returns list of pointers to the samples
inline std::list<MeasurementBase*> GetSampleList(){ return fSamples; }
//! Return list of pointers to all the pulls
inline std::list<ParamPull*> GetPullList(){ return fPulls; };
//! Write all samples to output DIR
void Write();
//! Set Fake data from file/MC
void SetFakeData(std::string fakeinput);
//! Reconfigure looping over duplicate inputs
void ReconfigureUsingManager();
//! Reconfigure Fast looping over duplicate inputs
void ReconfigureFastUsingManager();
std::vector<MeasurementBase*> GetSubSampleList();
std::vector<InputHandlerBase*> GetInputList();
private:
//! Append the experiments to include in the fit to this list
std::list<MeasurementBase*> fSamples;
//! Append parameter pull terms to include penalties in the fit to this list
std::list<ParamPull*> fPulls;
TDirectory *fOutputDir; //!< Directory to save contents
std::string fCardFile; //!< Input Card text file
bool fDialChanged; //!< Flag for if RW dials changed
UInt_t fCurIter; //!< Counter for how many times reconfigure called
bool fMCFilled; //!< Check MC has at least been filled once
TTree* fIterationTree; //!< Tree to save RW values on each function call
int fNDials; //!< Number of RW Dials in FitWeight
double* fDialVals; //!< Current Values of RW Dials
double fLikelihood; //!< Current likelihood for joint sample likelihood
double fNDOF; //!< Total NDOF
double* fSampleLikes; //!< Likelihoods for each individual measurement in list
int * fSampleNDOF; //!< NDOF for each individual measurement in list
bool fUsingEventManager; //!< Flag for doing joint comparisons
+ std::vector< std::vector<float> > fSignalEventSplines;
std::vector< std::vector<MeasurementVariableBox*> > fSignalEventBoxes;
std::vector< bool > fSignalEventFlags;
std::vector< std::vector<bool> > fSampleSignalFlags;
std::vector<InputHandlerBase*> fInputList;
std::vector<MeasurementBase*> fSubSampleList;
-
+ bool fIsAllSplines;
// iterator should go:
// Iterate over fSignalList
// If true get entry, and get 'next' in fSampleSignalList and fSampleSignalBoxes
// - Iterate over entry in fSampleSignalList
// - If true get 'next' in element in fSampleSignalBoxes.
// - - Fill input from MeasurementVariableBox*
};
/*! @} */
#endif // _JOINT_FCN_H_
diff --git a/src/FitBase/SplineInputHandler.cxx b/src/FitBase/SplineInputHandler.cxx
index 41686f6..befd0e2 100644
--- a/src/FitBase/SplineInputHandler.cxx
+++ b/src/FitBase/SplineInputHandler.cxx
@@ -1,148 +1,171 @@
#include "SplineInputHandler.h"
SplineInputHandler::SplineInputHandler(std::string const& handle, std::string const& rawinputs) {
// Run a joint input handling
fName = handle;
jointinput = false;
jointindexswitch = 0;
// Setup the TChain
fFitEventTree = new TChain("nuisance_events");
// Add to TChain
fFitEventTree->Add( rawinputs.c_str() );
// Open File for histogram access
int evtcounter = 0;
TFile* inp_file = new TFile(rawinputs.c_str(), "READ");
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "nuisance_fluxhist")).c_str());
TH1D* eventhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "nuisance_eventhist")).c_str());
// Get N events
TTree* neuttree = (TTree*)inp_file->Get("nuisance_events");
int nevents = neuttree->GetEntries();
// Push into individual input vectors
jointfluxinputs.push_back( (TH1D*) fluxhist->Clone() );
jointeventinputs.push_back( (TH1D*) eventhist->Clone() );
jointindexlow.push_back(evtcounter);
jointindexhigh.push_back(evtcounter + nevents);
evtcounter += nevents;
// Add to the total flux/event hist
if (!fFluxHist) fFluxHist = (TH1D*) fluxhist->Clone();
else fFluxHist->Add(fluxhist);
if (!fEventHist) fEventHist = (TH1D*) eventhist->Clone();
else fEventHist->Add(eventhist);
// Setup NEvents and the FitEvent
fNEvents = fFitEventTree->GetEntries();
fEventType = kNEWSPLINE;
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetBranchAddress(fFitEventTree);
// Setup Reader
fSplRead = new SplineReader();
fSplRead->Read( (TTree*)inp_file->Get("spline_reader") );
fNUISANCEEvent->fSplineRead = this->fSplRead;
+ // Load into memory
+ fNPar = fSplRead->GetNPar();
+
// Setup Friend Spline TTree
fSplTree = (TTree*)inp_file->Get("spline_tree");
fSplTree->SetBranchAddress( "SplineCoeff", fSplineCoeff );
fNUISANCEEvent->fSplineCoeff = this->fSplineCoeff;
fBaseEvent = static_cast<BaseFitEvt*>(fNUISANCEEvent);
+ // Load into memory
+ /*
+ for (int j = 0; j < fNEvents; j++){
+ std::vector<float> tempval;
+ fSplTree->GetEntry(j);
+ for (int i = 0; i < fNPar; i++){
+ tempval.push_back( fSplineCoeff[i] );
+ }
+ fAllSplineCoeff.push_back(tempval);
+ if (j % (fNEvents/10) == 0 and j != 0) {
+ LOG(SAM) << "Pushed " << int(j*100/fNEvents)+1 << "% of spline sets into memory for " << fName
+ << " (~" << int(sizeof(float)*tempval.size()*fAllSplineCoeff.size()/1.E6) <<"MB)" << std::endl;
+ }
+ }
+ */
+ fEventType = kSPLINEPARAMETER;
+
// Normalise event histograms for relative flux contributions.
for (size_t i = 0; i < jointeventinputs.size(); i++) {
TH1D* eventhist = (TH1D*) jointeventinputs.at(i)->Clone();
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= eventhist->Integral("width");
scale /= double(jointindexhigh[i] - jointindexlow[i]);
jointindexscale .push_back(scale);
}
fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str());
fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str());
};
FitEvent* SplineInputHandler::GetNuisanceEvent(const UInt_t entry) {
// Make sure events setup
if (entry >= fNEvents) return NULL;
// Reset all variables before tree read
fNUISANCEEvent->ResetEvent();
// Read NUISANCE Tree
fFitEventTree->GetEntry(entry);
fSplTree->GetEntry(entry);
+ fNUISANCEEvent->fSplineCoeff = fSplineCoeff;
fNUISANCEEvent->eventid = entry;
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent-> OrderStack();
// Setup Input scaling for joint inputs
if (jointinput) {
- fNUISANCEEvent->InputWeight *= GetInputWeight(entry);
+ fNUISANCEEvent->InputWeight *= fNUISANCEEvent->SavedRWWeight * GetInputWeight(entry);
} else {
- fNUISANCEEvent->InputWeight *= 1.0;
+ fNUISANCEEvent->InputWeight *= fNUISANCEEvent->SavedRWWeight;
}
return fNUISANCEEvent;
}
double SplineInputHandler::GetInputWeight(int entry) {
// Find Switch Scale
while ( entry < jointindexlow[jointindexswitch] ||
entry >= jointindexhigh[jointindexswitch] ) {
jointindexswitch++;
// Loop Around
if (jointindexswitch == jointindexlow.size()) {
jointindexswitch = 0;
}
}
return jointindexscale[jointindexswitch];
};
BaseFitEvt* SplineInputHandler::GetBaseEvent(const UInt_t entry) {
// Make sure events setup
if (entry >= fNEvents) return NULL;
// Read entry from TTree to fill NEUT Vect in BaseFitEvt;
// fFitEventTree->GetEntry(entry);
- fSplTree->GetEntry(entry);
+ //fSplTree->GetEntry(entry);
+
+ // fBaseEvent->fSplineCoeff = &fAllSplineCoeff[entry][0];
fBaseEvent->eventid = entry;
// Set joint scaling if required
if (jointinput) {
fBaseEvent->InputWeight *= GetInputWeight(entry);
} else {
fBaseEvent->InputWeight *= 1.0;
}
return fBaseEvent;
}
void SplineInputHandler::Print() {}
diff --git a/src/Reweight/GlobalDialList.cxx b/src/Reweight/GlobalDialList.cxx
index a4dea25..a3ef8f6 100644
--- a/src/Reweight/GlobalDialList.cxx
+++ b/src/Reweight/GlobalDialList.cxx
@@ -1,51 +1,56 @@
#include "GlobalDialList.h"
int GlobalDialList::EnumFromNameAndType(std::string name, int type){
// Setup Type Container
if (fTypeEnumCont.find(type) == fTypeEnumCont.end()){
std::map<std::string, int> temp;
temp[name] = 0;
fTypeEnumCont[type] = temp;
}
// Get Type Container
std::map<std::string, int> enumcont = fTypeEnumCont[type];
+ // LOG(FIT) << "Getting Enum From name and type " << name << " " << type << std::endl;
// Check name is in container, if its not add it
if (enumcont.find(name) == enumcont.end()){
+ // LOG(FIT) << " Name not found in Enum type " << type << std::endl;
int index = enumcont.size();
enumcont[name] = index;
+ // LOG(FIT) << "Returning new index of " << index << std::endl;
+ fTypeEnumCont[type][name] = index;
+ return index;
}
return enumcont[name];
}
void GlobalDialList::RegisterDialEnum(std::string name, int type, int nuisenum){
if (std::find(fAllDialNames.begin(), fAllDialNames.end(), name) != fAllDialNames.end()){
return;
}
- std::cout << "Registed Dial : " << name << " " << type << " " << nuisenum << std::endl;
+ LOG(FIT) << "Registed Dial Enum : " << name << " " << type << " " << nuisenum << std::endl;
fAllDialNames.push_back(name);
fAllDialTypes.push_back(type);
fAllDialEnums.push_back(nuisenum);
}
/// Singleton functions
GlobalDialList& Reweight::DialList(){
return GlobalDialList::Get();
}
GlobalDialList* GlobalDialList::m_diallistInstance = NULL;
GlobalDialList& GlobalDialList::Get(void){
if (!m_diallistInstance) m_diallistInstance = new GlobalDialList;
return *m_diallistInstance;
}
diff --git a/src/Reweight/SplineWeightEngine.cxx b/src/Reweight/SplineWeightEngine.cxx
index c0e9c9a..50505bd 100644
--- a/src/Reweight/SplineWeightEngine.cxx
+++ b/src/Reweight/SplineWeightEngine.cxx
@@ -1,80 +1,78 @@
#include "SplineWeightEngine.h"
SplineWeightEngine::SplineWeightEngine(std::string name) {
- // Setup the Reweight engien
- fCalcName = name;
- LOG(FIT) << "Setting up Spline RW : " << fCalcName << endl;
-
- // Set Abs Twk Config
- fIsAbsTwk = true;
-
+ // Setup the Reweight engien
+ fCalcName = name;
+ LOG(FIT) << "Setting up Spline RW : " << fCalcName << endl;
+
+ // Set Abs Twk Config
+ fIsAbsTwk = true;
+
};
void SplineWeightEngine::IncludeDial(std::string name, double startval) {
- // Get NUISANCE Enum
- int nuisenum = Reweight::ConvDial(name, kNEWSPLINE);
-
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
-
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
-
- std::cout << "Inlcuded Spline Dial " << name << " " << nuisenum << " " << startval << std::endl;
-
- // Set Value if given
- if (startval != -999.9) {
- SetDialValue(name, startval);
- }
-
+ // Get NUISANCE Enum
+ int nuisenum = Reweight::ConvDial(name, kSPLINEPARAMETER);
+
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+
+ fEnumIndex[nuisenum] = index;
+ fNameIndex[name] = index;
+
+ // std::cout << "Inlcuded Spline Dial " << name << " " << nuisenum << " " << startval << " " << index << std::endl;
+
+ // Set Value if given
+ if (startval != -999.9) {
+ SetDialValue(name, startval);
+ }
}
void SplineWeightEngine::SetDialValue(int nuisenum, double val) {
- LOG(FIT) << "Enum Val " << nuisenum << std::endl;
- LOG(FIT) << fEnumIndex.size() << std::endl;
- fValues[fEnumIndex[nuisenum]] = val;
+ // LOG(FIT) << "Enum Val " << nuisenum << std::endl;
+ // LOG(FIT) << fEnumIndex.size() << std::endl;
+ // LOG(FIT) << "Set Dial Value to " << nuisenum << " " << fEnumIndex[nuisenum] << " " << val << std::endl;
+ fValues[fEnumIndex[nuisenum]] = val;
}
void SplineWeightEngine::SetDialValue(std::string name, double val){
- fValues[fNameIndex[name]] = val;
+ fValues[fNameIndex[name]] = val;
}
void SplineWeightEngine::Reconfigure(bool silent) {
- for (std::map<std::string, size_t>::iterator iter = fNameIndex.begin();
- iter != fNameIndex.end(); iter++){
-
- LOG(FIT) << "Reconfiguring Spline " << iter->first << " to be " << fValues[ iter->second ] << " Inside SPL RW" << std::endl;
- fSplineValueMap[ iter->first ] = fValues[ iter->second ];
- }
+ for (std::map<std::string, size_t>::iterator iter = fNameIndex.begin();
+ iter != fNameIndex.end(); iter++){
+ // LOG(FIT) << "Reconfiguring Spline " << iter->first << " to be " << fValues[ iter->second ] << " Inside SPL RW" << std::endl;
+ fSplineValueMap[ iter->first ] = fValues[ iter->second ];
+ }
}
double SplineWeightEngine::CalcWeight(BaseFitEvt* evt) {
- if (!evt->fSplineRead) return 1.0;
-
- // if (!evt->fSplineRead->NeedsReconfigure()) {
- evt->fSplineRead->Reconfigure(fSplineValueMap);
- // return 1.0;
- // }
- double rw_weight = evt->fSplineRead->CalcWeight( evt->fSplineCoeff );
-
- return rw_weight;
-
+ if (!evt->fSplineRead) return 1.0;
+
+ if (evt->fSplineRead->NeedsReconfigure()) {
+ evt->fSplineRead->Reconfigure(fSplineValueMap);
+ }
+
+ double rw_weight = evt->fSplineRead->CalcWeight( evt->fSplineCoeff );
+
+ return rw_weight;
}
diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index 989f1e0..d3a1d0e 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,485 +1,486 @@
#include "WeightUtils.h"
//********************************************************************
TF1 FitBase::GetRWConvFunction(std::string type, std::string name) {
//********************************************************************
std::string dialfunc = "x";
std::string parType = type;
double low = -10000.0;
double high = 10000.0;
if (parType.find("parameter") == std::string::npos) parType += "_parameter";
string line;
ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card").c_str(),
ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 4) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
// inputlist[2] should be the units... ignore for now
dialfunc = inputlist[3];
// High and low are optional, check whether they exist
if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]);
if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]);
}
TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high);
return convfunc;
}
//********************************************************************
std::string FitBase::GetRWUnits(std::string type, std::string name) {
//********************************************************************
std::string unit = "sig.";
std::string parType = type;
if (parType.find("parameter") == std::string::npos) {
parType += "_parameter";
}
std::string line;
std::ifstream card((GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card").c_str(), ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 3) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
unit = inputlist[2];
break;
}
return unit;
}
//********************************************************************
double FitBase::RWAbsToSigma(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX(val);
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToAbs(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val);
return conv_val;
}
//********************************************************************
double FitBase::RWFracToSigma(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX((val * f1.Eval(0.0)));
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToFrac(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val) / f1.Eval(0.0);
return conv_val;
}
int FitBase::ConvDialType(std::string type) {
if (!type.compare("neut_parameter")) return kNEUT;
else if (!type.compare("niwg_parameter")) return kNIWG;
else if (!type.compare("nuwro_parameter")) return kNUWRO;
else if (!type.compare("t2k_parameter")) return kT2K;
else if (!type.compare("genie_parameter")) return kGENIE;
else if (!type.compare("custom_parameter")) return kCUSTOM;
else if (!type.compare("norm_parameter")) return kNORM;
else if (!type.compare("modenorm_parameter")) return kMODENORM;
else if (!type.compare("likeweight_parameter")) return kLIKEWEIGHT;
else if (!type.compare("spline_parameter")) return kSPLINEPARAMETER;
else return kUNKNOWN;
}
std::string FitBase::ConvDialType(int type) {
switch (type) {
case kNEUT: { return "neut_parameter"; }
case kNIWG: { return "niwg_parameter"; }
case kNUWRO: { return "nuwro_parameter"; }
case kT2K: { return "t2k_parameter"; }
case kGENIE: { return "genie_parameter"; }
case kNORM: { return "norm_parameter"; }
case kCUSTOM: {return "custom_parameter"; }
case kMODENORM: { return "modenorm_parameter"; }
case kLIKEWEIGHT: { return "likeweight_parameter"; }
case kSPLINEPARAMETER: {return "spline_parameter";}
default: return "unknown_parameter";
}
}
int FitBase::GetDialEnum(std::string type, std::string name) {
return FitBase::GetDialEnum( FitBase::ConvDialType(type), name );
}
int FitBase::GetDialEnum(int type, std::string name) {
int offset = type * 1000;
int this_enum = -1; //Not Found
std::cout << "Getting dial enum " << type << " " << name << std::endl;
// Select Types
switch (type) {
// NEUT DIAL TYPE
case kNEUT: {
#ifdef __NEUT_ENABLED__
int neut_enum = (int)neut::rew::NSyst::FromString(name);
if (neut_enum != 0) { this_enum = neut_enum + offset; }
#else
this_enum = -2; //Not enabled
#endif
break;
}
// NIWG DIAL TYPE
case kNIWG: {
#ifdef __NIWG_ENABLED__
int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
if (niwg_enum != 0) { this_enum = niwg_enum + offset; }
#else
this_enum = -2;
#endif
break;
}
// NUWRO DIAL TYPE
case kNUWRO: {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
if (nuwro_enum > 0) { this_enum = nuwro_enum + offset; }
#else
this_enum = -2;
#endif
}
// GENIE DIAL TYPE
case kGENIE: {
#ifdef __GENIE_ENABLED__
int genie_enum = (int)genie::rew::GSyst::FromString(name);
if (genie_enum > 0) { this_enum = genie_enum + offset; }
#else
this_enum = -2;
#endif
break;
}
case kCUSTOM: {
int custom_enum = 0; // PLACEHOLDER
this_enum = custom_enum + offset;
break;
}
// T2K DIAL TYPE
case kT2K: {
#ifdef __T2KREW_ENABLED__
int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
if (t2k_enum > 0) { this_enum = t2k_enum + offset; }
#else
this_enum = -2;
#endif
break;
}
case kNORM: {
if (gNormEnums.find(name) == gNormEnums.end()) {
gNormEnums[name] = gNormEnums.size() + 1 + offset;
}
this_enum = gNormEnums[name];
break;
}
case kMODENORM: {
size_t us_pos = name.find_first_of('_');
std::string numstr = name.substr(us_pos + 1);
int mode_num = std::atoi(numstr.c_str());
LOG(FTL) << "Getting mode num " << mode_num << endl;
if (!mode_num) {
ERR(FTL) << "Attempting to parse dial name: \"" << name
<< "\" as a mode norm dial but failed." << endl;
throw;
}
this_enum = 60 + mode_num + offset;
break;
}
case kLIKEWEIGHT: {
if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
}
this_enum = gLikeWeightEnums[name];
break;
}
case kSPLINEPARAMETER: {
if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
}
this_enum = gSplineParameterEnums[name];
}
}
// If Not Enabled
if (this_enum == -2) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << endl;
ERR(FTL) << "Check dial " << name << endl;
}
// If Not Found
if (this_enum == -1) {
ERR(FTL) << "Dial " << name << " not found." << endl;
}
return this_enum;
}
int Reweight::ConvDialType(std::string type) {
if (!type.compare("neut_parameter")) return kNEUT;
else if (!type.compare("niwg_parameter")) return kNIWG;
else if (!type.compare("nuwro_parameter")) return kNUWRO;
else if (!type.compare("t2k_parameter")) return kT2K;
else if (!type.compare("genie_parameter")) return kGENIE;
else if (!type.compare("norm_parameter")) return kNORM;
else if (!type.compare("modenorm_parameter")) return kMODENORM;
else if (!type.compare("custom_parameter")) return kCUSTOM;
else if (!type.compare("likeweight_parameter")) return kLIKEWEIGHT;
else if (!type.compare("spline_parameter")) return kSPLINEPARAMETER;
else return kUNKNOWN;
}
std::string Reweight::ConvDialType(int type) {
switch (type) {
case kNEUT: { return "neut_parameter"; }
case kNIWG: { return "niwg_parameter"; }
case kNUWRO: { return "nuwro_parameter"; }
case kT2K: { return "t2k_parameter"; }
case kGENIE: { return "genie_parameter"; }
case kNORM: { return "norm_parameter"; }
case kCUSTOM: {return "custom_parameter"; }
case kMODENORM: { return "modenorm_parameter"; }
case kLIKEWEIGHT: { return "likeweight_parameter"; }
case kSPLINEPARAMETER: {return "spline_parameter";}
default: return "unknown_parameter";
}
}
int Reweight::NEUTEnumFromName(std::string name) {
#ifdef __NEUT_ENABLED__
int neutenum = (int)neut::rew::NSyst::FromString(name);
return (neutenum > 0) ? neutenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::NIWGEnumFromName(std::string name) {
#ifdef __NIWG_ENABLED__
int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name);
return (niwgenum != 0) ? niwgenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::NUWROEnumFromName(std::string name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name);
return (nuwroenum > 0) ? nuwroenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::GENIEEnumFromName(std::string name) {
#ifdef __GENIE_ENABLED__
int genieenum = (int)genie::rew::GSyst::FromString(name);
return (genieenum > 0) ? genieenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::T2KEnumFromName(std::string name) {
#ifdef __T2KREW_ENABLED__
int t2kenum = (int)t2krew::T2KSyst::FromString(name);
return (t2kenum > 0) ? t2kenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::NUISANCEEnumFromName(std::string name, int type) {
int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type);
return nuisenum;
}
int Reweight::CustomEnumFromName(std::string name) {
int custenum = Reweight::ConvertNUISANCEDial(name);
return custenum;
}
int Reweight::ConvDial(std::string name, std::string type, bool exceptions) {
return Reweight::ConvDial( name, Reweight::ConvDialType(type), exceptions );
}
int Reweight::ConvDial(std::string name, int type, bool exceptions) {
// Produce offset seperating each type.
int offset = type * 1000;
int genenum = kNoDialFound;
switch (type) {
case kNEUT:
genenum = NEUTEnumFromName(name);
break;
case kNIWG:
genenum = NIWGEnumFromName(name);
break;
case kNUWRO:
genenum = NUWROEnumFromName(name);
break;
case kGENIE:
genenum = GENIEEnumFromName(name);
break;
case kT2K:
genenum = T2KEnumFromName(name);
break;
case kCUSTOM:
genenum = CustomEnumFromName(name);
break;
case kNORM:
case kMODENORM:
case kLIKEWEIGHT:
case kSPLINEPARAMETER:
- genenum = NUISANCEEnumFromName(name, kLIKEWEIGHT);
+ case kNEWSPLINE:
+ genenum = NUISANCEEnumFromName(name, type);
break;
default:
genenum = kNoTypeFound;
break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == kGeneratorNotBuilt) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << endl;
ERR(FTL) << "Check dial " << name << endl;
throw;
}
// If no type enabled
if (genenum == kNoTypeFound) {
ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl;
throw;
}
// If Not Found
if (genenum == kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << endl;
throw;
}
}
// Add offset if no issue
int nuisenum = genenum;
if (genenum != kGeneratorNotBuilt and
genenum != kNoTypeFound and
genenum != kNoDialFound) {
nuisenum += offset;
}
// Now register dial
// std::cout << "Returning " << nuisenum << std::endl;
Reweight::DialList().RegisterDialEnum(name, type, nuisenum);
return nuisenum;
}
std::string Reweight::ConvDial(int nuisenum) {
GlobalDialList* temp;
for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) {
if (Reweight::DialList().fAllDialEnums[i] == nuisenum) {
return Reweight::DialList().fAllDialNames[i];
}
}
LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl;
}
diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx
index 64d1231..19275e8 100755
--- a/src/Routines/MinimizerRoutines.cxx
+++ b/src/Routines/MinimizerRoutines.cxx
@@ -1,1436 +1,1437 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "StatusMessage.h"
#include "MinimizerRoutines.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");
};
//*************************************
MinimizerRoutines::~MinimizerRoutines() {
//*************************************
};
/*
Input Functions
*/
//*************************************
MinimizerRoutines::MinimizerRoutines(int argc, char* argv[]) {
//*************************************
// Set everything to defaults
Init();
nuisconfig configuration = Config::Get();
std::string cardfile = "";
int maxevents = -1;
- int errorcount = Config::Get().GetParI("ERROR");
- int verbocount = Config::Get().GetParI("VERBOSITY");
+ int errorcount = 0;//Config::Get().GetParI("ERROR");
+ int verbocount = 0;//Config::Get().GetParI("VERBOSITY");
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
// Make easier to handle arguments.
std::vector<std::string> 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::ParseSplitArgument(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 supplied!" << std::endl;
throw;
}
// Setup this configuration
fCompKey = Config::Get().CreateNode("nuiscomp");
fCompKey.AddS("cardfile", fCardFile);
fCompKey.AddS("outputfile", fOutputFile);
fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
// std::cout << "Adding XML Line " << xmlcmds[i] << std::endl;
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
// Add Error Verbo Lines
+ verbocount += FitPar::Config().GetParI("VERBOSITY");
FitPar::log_verb = verbocount;
LOG_VERB(verbocount);
ERR_VERB(errorcount);
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Proper 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<nuiskey> 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;
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<nuiskey> 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 + sampletype + "_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;
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<nuiskey> 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 = "";
// 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 = ""; }
// make minimizer
if (fMinimizer) delete fMinimizer;
fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype);
fMinimizer->SetMaxFunctionCalls(FitPar::Config().GetParI("minimizer.maxcalls"));
if (!routine.compare("Brute")) {
fMinimizer->SetMaxFunctionCalls(fParams.size() * fParams.size() * 4);
fMinimizer->SetMaxIterations(fParams.size() * fParams.size() * 4);
}
fMinimizer->SetMaxIterations(FitPar::Config().GetParI("minimizer.maxiterations"));
fMinimizer->SetTolerance(FitPar::Config().GetParD("minimizer.tolerance"));
fMinimizer->SetStrategy(FitPar::Config().GetParI("minimizer.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<std::string, double>& 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.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." << endl;
break;
}
}
return;
}
//*************************************
int MinimizerRoutines::RunFitRoutine(std::string routine) {
//*************************************
int endfits = kFitUnfinished;
// set fitter at the current start values
fOutputRootFile->cd();
SetupFitter(routine);
// choose what to do with the minimizer depending on routine.
if (!routine.compare("Migrad") or
!routine.compare("Simplex") or
!routine.compare("Combined") or
!routine.compare("Brute") or
!routine.compare("Fumili") or
!routine.compare("ConjugateFR") or
!routine.compare("ConjugatePR") or
!routine.compare("BFGS") or
!routine.compare("BFGS2") or
!routine.compare("SteepDesc") or
// !routine.compare("GSLMulti") or
!routine.compare("GSLSimAn")) {
if (fMinimizer->NFree() > 0) {
LOG(FIT) << StatusMessage(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() << endl;
}
LOG(FIT) << "------------" << std::endl;
double like = fSampleFCN->GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << 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
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");
}
}
return;
};
//*************************************
void MinimizerRoutines::LowStatRoutine(std::string routine) {
//*************************************
LOG(FIT) << "Running Low Statistics Routine: " << routine << std::endl;
int lowstatsevents = FitPar::Config().GetParI("minimizer.lowstatevents");
int maxevents = FitPar::Config().GetParI("input.maxevents");
int verbosity = FitPar::Config().GetParI("VERBOSITY");
std::string trueroutine = routine;
std::string substring = "LowStat";
trueroutine.erase( trueroutine.find(substring),
substring.length() );
// Set MAX EVENTS=1000
FitPar::Config().SetParI("input.maxevents", lowstatsevents);
FitPar::Config().SetParI("VERBOSITY", 3);
SetupFCN();
RunFitRoutine(trueroutine);
FitPar::Config().SetParI("input.maxevents", maxevents);
SetupFCN();
FitPar::Config().SetParI("VERBOSITY", verbosity);
return;
}
//*************************************
void MinimizerRoutines::Create1DScans() {
//*************************************
// 1D Scan Routine
// Steps through all free parameters about nominal using the step size
// Creates a graph for each free parameter
// At the current point create a 1D Scan for all parametes (Uncorrelated)
for (UInt_t i = 0; i < fParams.size(); i++) {
if (fFixVals[fParams[i]]) continue;
LOG(FIT) << "Running 1D Scan for " << fParams[i] << 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] << 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!" << 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!" << endl;
return kNoChange;
} else return kStateChange;
}
/*
Write Functions
*/
//*************************************
void MinimizerRoutines::SaveResults() {
//*************************************
fOutputRootFile->cd();
if (fMinimizer) {
SetupCovariance();
SaveMinimizerState();
}
SaveCurrentState();
}
//*************************************
void MinimizerRoutines::SaveMinimizerState() {
//*************************************
if (!fMinimizer) {
ERR(FTL) << "Can't save minimizer state without min object" << endl;
throw;
}
// Save main fit tree
fSampleFCN->WriteIterationTree();
// Get Vals and Errors
GetMinimizerState();
// Save tree with fit status
std::vector<std::string> nameVect;
std::vector<double> valVect;
std::vector<double> errVect;
std::vector<double> minVect;
std::vector<double> maxVect;
std::vector<double> startVect;
std::vector<int> endfixVect;
std::vector<int> 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.." << endl;
int NFREE = 0;
int NDIM = 0;
// Get NFREE from min or from vals (for cases when doing throws)
if (fMinimizer) {
NFREE = fMinimizer->NFree();
NDIM = fMinimizer->NDim();
} else {
NDIM = fParams.size();
for (UInt_t i = 0; i < fParams.size(); i++) {
if (!fFixVals[fParams[i]]) NFREE++;
}
}
if (NDIM == 0) return;
LOG(FIT) << "NFREE == " << NFREE << 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++) {
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++;
}
}
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;
}
return;
};
//*************************************
void MinimizerRoutines::ThrowCovariance(bool uniformly) {
//*************************************
std::vector<double> 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;
};
diff --git a/src/Routines/SplineRoutines.cxx b/src/Routines/SplineRoutines.cxx
index 983ee96..7bc2183 100755
--- a/src/Routines/SplineRoutines.cxx
+++ b/src/Routines/SplineRoutines.cxx
@@ -1,562 +1,592 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "StatusMessage.h"
#include "SplineRoutines.h"
/*
Constructor/Destructor
*/
//************************
void SplineRoutines::Init() {
//************************
fStrategy = "SaveEvents";
fRoutines.clear();
fCardFile = "";
fSampleFCN = NULL;
fRW = NULL;
fAllowedRoutines = ("SaveEvents,TestEvents,SaveSplineEvents");
};
//*************************************
SplineRoutines::~SplineRoutines() {
//*************************************
};
/*
Input Functions
*/
//*************************************
SplineRoutines::SplineRoutines(int argc, char* argv[]) {
//*************************************
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
std::string cardfile = "";
int maxevents = -1;
int errorcount = Config::Get().GetParI("ERROR");
int verbocount = Config::Get().GetParI("VERBOSITY");
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
// Make easier to handle arguments.
std::vector<std::string> 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::ParseSplitArgument(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 supplied!" << std::endl;
throw;
}
// Setup this configuration
fCompKey = Config::Get().CreateNode("nuiscomp");
fCompKey.AddS("cardfile", fCardFile);
fCompKey.AddS("outputfile", fOutputFile);
fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
// std::cout << "Adding XML Line " << xmlcmds[i] << std::endl;
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
// Add Error Verbo Lines
FitPar::log_verb = verbocount;
LOG_VERB(verbocount);
ERR_VERB(errorcount);
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Starting Setup
// ---------------------------
SetupRWEngine();
return;
};
/*
Setup Functions
*/
//*************************************
void SplineRoutines::SetupRWEngine() {
//*************************************
fRW = new FitWeight("splineweight");
// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
// Add Parameters
for (int i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
std::string parname = key.GetS("name");
std::string partype = key.GetS("type");
double nom = key.GetD("nominal");
fRW->IncludeDial( key.GetS("name"),
FitBase::ConvDialType(key.GetS("type")), nom);
fRW->SetDialValue( key.GetS("name"), key.GetD("nominal") );
}
fRW->Reconfigure();
// for (int i = 0; i < parameterkeys.size(); i++){
// nuiskey key = parameterkeys[i];
// fRW->IncludeDial( key.GetS("name"),
// FitBase::ConvDialType(key.GetS("type")));
// }
// LOG(FIT) << "Setting up FitWeight Engine" << std::endl;
// for (UInt_t i = 0; i < fParams.size(); i++) {
// std::string name = fParams[i];
// FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name));
// }
// UpdateRWEngine(fStartVals);
return;
}
/*
Fitting Functions
*/
//*************************************
void SplineRoutines::UpdateRWEngine(std::map<std::string, double>& updateVals) {
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
fRW->SetDialValue(name, updateVals.at(name));
}
fRW->Reconfigure();
return;
}
//*************************************
void SplineRoutines::Run() {
//*************************************
std::cout << "Running " << std::endl;
// Parse given routines
fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
if (fRoutines.empty()) {
ERR(FTL) << "Trying to run ComparisonRoutines with no routines given!" << std::endl;
throw;
}
for (int i = 0; i < fRoutines.size(); i++) {
LOG(FIT) << "Running Routine: " << fRoutines[i] << std::endl;
std::string rout = fRoutines[i];
if (!rout.compare("SaveEvents")) SaveEvents();
else if (!rout.compare("TestEvents")) TestEvents();
else if (!rout.compare("GenerateEventSplines")) GenerateEventSplines();
}
}
//*************************************
void SplineRoutines::SaveEvents() {
//*************************************
if (fRW) delete fRW;
SetupRWEngine();
fRW->Reconfigure();
fRW->Print();
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (int i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
// Get info from inputhandler
int nevents = input->GetNEvents();
int countwidth = (nevents / 10);
FitEvent* nuisevent = input->FirstNuisanceEvent();
// Setup a TTree to save the event
outputfile->cd();
TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
nuisevent->AddBranchesToTree(eventtree);
// Loop over all events and fill the TTree
int i = 0;
// int countwidth = nevents / 5;
while (nuisevent) {
// Get Event Weight
nuisevent->RWWeight = fRW->CalcWeight(nuisevent);
// if (nuisevent->RWWeight != 1.0){
// std::cout << "Weight = " << nuisevent->RWWeight << std::endl;
// }
// Save everything
eventtree->Fill();
// Logging
if (i % countwidth == 0) {
LOG(REC) << "Saved " << i << "/" << nevents
<< " nuisance events. [M, W] = ["
<< nuisevent->Mode << ", " << nuisevent->RWWeight << "]" << std::endl;
}
// iterate
nuisevent = input->NextNuisanceEvent();
i++;
}
// Save flux and close file
outputfile->cd();
eventtree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
// remove Keys
eventkeys.clear();
// Finished
LOG(FIT) << "Finished processing all nuisance events." << std::endl;
}
//*************************************
void SplineRoutines::TestEvents() {
//*************************************
LOG(FIT) << "Testing events." << std::endl;
// Create a new file for the test samples
if (!fOutputRootFile) {
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
}
// Loop over all tests
int count = 0;
std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");
for (std::vector<nuiskey>::iterator iter = testkeys.begin();
iter != testkeys.end(); iter++) {
nuiskey key = (*iter);
// 0. Create new measurement list
std::list<MeasurementBase*> samplelist;
// 1. Build Sample From Events
std::string samplename = key.GetS("name");
std::string eventsid = key.GetS("inputid");
nuiskey eventskey = Config::QueryLastKey("events", "id=" + eventsid);
std::string rawfile = eventskey.GetS("input");
LOG(FIT) << "Creating sample " << samplename << std::endl;
MeasurementBase* rawsample = SampleUtils::CreateSample(samplename, rawfile, "", "", FitBase::GetRW());
// 2. Build Sample From Nuisance Events
std::string eventsfile = eventskey.GetS("output");
LOG(FIT) << "Creating Fit Eevnt Sample " << samplename << " " << eventsfile << std::endl;
MeasurementBase* nuissample = SampleUtils::CreateSample(samplename, "FEVENT:" + eventsfile, "", "", FitBase::GetRW());
// 3. Make some folders to save stuff
TDirectory* sampledir = (TDirectory*) fOutputRootFile->mkdir(Form((samplename + "_test_%d").c_str(), count));
TDirectory* rawdir = (TDirectory*) sampledir->mkdir("raw");
TDirectory* nuisancedir = (TDirectory*) sampledir->mkdir("nuisance");
TDirectory* difdir = (TDirectory*) sampledir->mkdir("difference");
// 4. Reconfigure both
rawdir->cd();
rawsample->Reconfigure();
rawsample->Write();
nuisancedir->cd();
nuissample->Reconfigure();
nuissample->Write();
// 4. Compare Raw to Nuisance Events
// Loop over all keyse
TIter next(rawdir->GetListOfKeys());
TKey *dirkey;
while ((dirkey = (TKey*)next())) {
// If not a 1D/2D histogram skip
TClass *cl = gROOT->GetClass(dirkey->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
// Get TH1* from both dir
TH1 *rawplot = (TH1*)rawdir->Get(dirkey->GetName());
TH1 *nuisanceplot = (TH1*)nuisancedir->Get(dirkey->GetName());
// Take Difference
nuisanceplot->Add(rawplot, -1.0);
// Save to dif folder
difdir->cd();
nuisanceplot->Write();
}
// 5. Tidy Up
samplelist.clear();
// Iterator
count++;
}
}
//*************************************
void SplineRoutines::GenerateEventSplines() {
//*************************************
if (fRW) delete fRW;
SetupRWEngine();
// Setup the spline reader
SplineWriter* splwrite = new SplineWriter(fRW);
std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
// Add splines to splinewriter
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
splwrite->AddSpline(splkey);
}
splwrite->SetupSplineSet();
// Event Loop
// Loop over all events and calculate weights for each parameter set.
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (int i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
// Get info from inputhandler
int nevents = input->GetNEvents();
int countwidth = (nevents / 1000);
FitEvent* nuisevent = input->FirstNuisanceEvent();
// Setup a TTree to save the event
outputfile->cd();
TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
+
+ // Add a flag that allows just splines to be saved.
nuisevent->AddBranchesToTree(eventtree);
// Save the spline reader
splwrite->Write("spline_reader");
// Setup the spline TTree
TTree* splinetree = new TTree("spline_tree", "spline_tree");
splwrite->AddCoefficientsToTree(splinetree);
int lasttime = time(NULL);
// Loop over all events and fill the TTree
while (nuisevent) {
// std::cout << "Fitting event " << i << std::endl;
// Calculate the weights for each parameter set
splwrite->FitSplinesForEvent(nuisevent);
// Save everything
eventtree->Fill();
splinetree->Fill();
// nuisevent->Print();
// std::cout << "Done with event " << i << std::endl;
// sleep(4);
// Logging
if (i % countwidth == 0) {
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (i != 0 and timeelapsed){
lasttime = time(NULL);
int eventsleft = nevents - i;
float speed = float(countwidth) / float(timeelapsed);
float proj = (float(eventsleft)/float(speed))/60/60;
timestring << proj << " hours remaining.";
}
LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
}
// Iterate
i++;
nuisevent = input->NextNuisanceEvent();
}
// Save flux and close file
outputfile->cd();
eventtree->Write();
splinetree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
// remove Keys
eventkeys.clear();
}
+/*
+//*************************************
+void SplineRoutines::MergeSplines() {
+//*************************************
+ // Loop over all 'splinemerge' keys.
+ // Add them to the Merger.
+ // Call setup splines.
+
+ // Get the key with eventinput
+ // - remaining keys should have splineinput
+ // - Loop over number of entries.
+ // - FillEntry in merger.
+ // - Fill NUISANCEEvent into a new TTree.
+
+
+
+
+
+
+
+
+
+
+
+
+}
+*/
+
//*************************************
// void SplineRoutines::TestEventSplines() {
//*************************************
// Make a spline fit weight and a normal fit weight.
// Loop over test samples.
// Loop over splines
// Set spline parameters and rw parameters to same value
// Compare weight response by binning against sample set.
// }
/*
MISC Functions
*/
//*************************************
int SplineRoutines::GetStatus() {
//*************************************
return 0;
}
diff --git a/src/Splines/CMakeLists.txt b/src/Splines/CMakeLists.txt
index 870cb86..d62ef5a 100644
--- a/src/Splines/CMakeLists.txt
+++ b/src/Splines/CMakeLists.txt
@@ -1,60 +1,62 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
set(IMPLFILES
FitSpline.cxx
FitSplineHead.cxx
SplineReader.cxx
SplineWriter.cxx
+SplineMerger.cxx
Spline.cxx
)
set(HEADERFILES
FitSpline.h
FitSplineHead.h
SplineReader.h
SplineWriter.h
+SplineMerger.h
Spline.h
)
set(LIBNAME Splines)
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
add_library(${LIBNAME} STATIC ${IMPLFILES})
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
add_library(${LIBNAME} SHARED ${IMPLFILES})
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
set_target_properties(${LIBNAME} PROPERTIES VERSION
"${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}")
set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})
if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES)
add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES})
endif()
install(TARGETS ${LIBNAME} DESTINATION lib)
#Can uncomment this to install the headers... but is it really neccessary?
install(FILES ${HEADERFILES} DESTINATION include)
set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
diff --git a/src/Splines/Spline.cxx b/src/Splines/Spline.cxx
index 564f77e..f8f3717 100644
--- a/src/Splines/Spline.cxx
+++ b/src/Splines/Spline.cxx
@@ -1,326 +1,326 @@
#include "Spline.h"
Spline::Spline(std::string splname, std::string form,
std::vector<double> x) {
// Initial setup
fSplineNames = GeneralUtils::ParseToStr(splname, ":");
fName = splname;
fForm = form;
fXScan.clear();
for (size_t i = 0; i < x.size(); i++){
fXScan.push_back(x[i]);
}
fX = 0.0;
// Set X Min/Max from scan points
fXMax = -99999;
fXMin = 99999;
for (size_t i = 0; i < fXScan.size(); i++) {
if (fXScan[i] > fXMax) fXMax = fXScan[i];
if (fXScan[i] < fXMin) fXMin = fXScan[i];
}
// Set form from list
if (!fForm.compare("1DPol1")) { Setup( k1DPol1, 1, 2 ); }
else if (!fForm.compare("1DPol2")) { Setup( k1DPol2, 1, 3 ); }
else if (!fForm.compare("1DPol3")) { Setup( k1DPol3, 1, 4 ); }
else if (!fForm.compare("1DPol4")) { Setup( k1DPol4, 1, 5 ); }
else if (!fForm.compare("1DPol5")) { Setup( k1DPol5, 1, 6 ); }
else if (!fForm.compare("1DPol6")) { Setup( k1DPol6, 1, 7 ); }
else if (!fForm.compare("1DTSpline3")) { Setup( k1DTSpline3, 1, fXScan.size() * 4 ); }
else {
ERR(FTL) << "Unknown spline form : " << fForm << std::endl;
throw;
}
// Run Checks
if (fNDim != fSplineNames.size()) {
ERR(FTL) << "Spline Dim:Names mismatch!" << std::endl;
throw;
}
};
void Spline::Setup(int type, int ndim, int npar) {
fType = type;
fNDim = ndim;
fNPar = npar;
}
double Spline::operator()(const Double_t* x, const Double_t* par) const {
double val = DoEval(x, par);
if (val < 0.0) val = 0.0;
return val;
}
float Spline::operator()(const Float_t* x, const Float_t* par) const {
float val = DoEval(x, par);
if (val < 0.0) val = 0.0;
return val;
}
void Spline::Reconfigure(double x) {
- // std::cout << "Reconfigured spline : " << fName << " : " << fForm << " to be " << x << std::endl;
+ // std::cout << "Reconfigured spline : " << fName << " : " << fForm << " to be " << x << std::endl;
fX = x;
fOutsideLimits = false;
if (fX > fXMax) fX = fXMax;
else if (fX < fXMin) fX = fXMin;
}
void Spline::Reconfigure(float x) {
// std::cout << "Reconfigured spline : " << fName << " : " << fForm << " to be " << x << std::endl;
fX = x;
fOutsideLimits = false;
if (fX > fXMax) fX = fXMax;
else if (fX < fXMin) fX = fXMin;
}
double Spline::DoEval(const Double_t* x, const Double_t* par) const {
// Setup current fX to value
fX = x[0];
if (fX > fXMax) fX = fXMax;
else if (fX < fXMin) fX = fXMin;
double w = DoEval(&par[0], false);
if (w < 0.0) w = 0.0;
// Now evaluate spline how FitWeight will do it.
return w;
}
float Spline::DoEval(const Float_t* x, const Float_t* par) const {
// Setup current fX to value
fX = x[0];
if (fX > fXMax) fX = fXMax;
else if (fX < fXMin) fX = fXMin;
double w = DoEval(&par[0], false);
if (w < 0.0) w = 0.0;
// Now evaluate spline how FitWeight will do it.
return w;
}
double Spline::DoEval(const Double_t* par, bool checkresponse) const {
Float_t* temp = new Float_t[fNPar];
for (size_t i = 0; i < fNPar; i++){
temp[i] = par[i];
}
double val = DoEval(temp, checkresponse);
delete temp;
return val;
}
float Spline::DoEval(const Float_t* par, bool checkresponse) const {
// Check response
if (checkresponse) {
bool hasresponse = false;
for (int i = 0; i < fNPar; i++) {
if (par[i] != 0.0) {
hasresponse = true;
break;
}
}
if (!hasresponse) {
// std::cout << "No Response" << std::endl;
return 1.0;
}
}
// Now evaluate spline
switch (fType) {
case k1DPol1: { return Spline1DPol1(par); }
case k1DPol2: { return Spline1DPol2(par); }
case k1DPol3: { return Spline1DPol3(par); }
case k1DPol4: { return Spline1DPol4(par); }
case k1DPol5: { return Spline1DPol5(par); }
case k1DPol6: { return Spline1DPol6(par); }
case k1DTSpline3: { return Spline1DTSpline3(par); }
}
// Return nominal weight
return 1.0;
};
// Spline Functions
// ----------------------------------------------
float Spline::Spline1DPol1(const Float_t* par) const {
float xp = fX;
- // std::cout << "Eval 1DPol1 with " << par[0] << " " << par[1] << " " << xp << " " << par[0] + par[1] * xp << std::endl;
+ //std::cout << "Eval 1DPol1 with " << par[0] << " " << par[1] << " " << xp << " " << par[0] + par[1] * xp << std::endl;
return par[0] + par[1] * xp;
};
float Spline::Spline1DPol2(const Float_t* par) const {
float xp = fX;
return par[0] + par[1] * xp + par[2] * xp * xp;
};
float Spline::Spline1DPol3(const Float_t* par) const {
float xp = fX;
return par[0] + par[1] * xp + par[2] * xp * xp + par[3] * xp * xp * xp;
};
float Spline::Spline1DPol4(const Float_t* par) const {
float xp = fX;
return (par[0] + par[1] * xp + par[2] * xp * xp + par[3] * xp * xp * xp +
par[4] * xp * xp * xp * xp);
};
float Spline::Spline1DPol5(const Float_t* par) const {
float xp = fX;
return (par[0] + par[1] * xp + par[2] * xp * xp + par[3] * xp * xp * xp +
par[4] * xp * xp * xp * xp + par[5] * xp * xp * xp * xp * xp);
};
float Spline::Spline1DPol6(const Float_t* par) const {
float xp = fX;
float w = 0.0;
// std::cout << "Pol Eval " << std::endl;
for (int i = fNPar-1; i > 0; i--){
w = xp * (par[0+i] + w);
}
w += par[0];
return w;
};
float Spline::Spline1DTSpline3(const Float_t* par) const {
// iterator over knot values and find width
std::vector<float>::iterator iter_low = fXScan.begin();
std::vector<float>::iterator iter_high = fXScan.begin();
iter_high++;
int off = 0;
// Find matching point
while ( iter_high != fXScan.end() ) {
if (fX >= (*iter_low) and fX < (*iter_high)) {
break;
}
off += 4;
iter_low++;
iter_high++;
}
float dx = fX - (*iter_low);
float weight = (par[off] + dx * (par[off + 1] + dx * (par[off + 2] + dx * par[off + 3])));
return weight;
};
void Spline::FitCoeff(int n, double* x, double* y, float* coeff, bool draw) {
switch (fType) {
// Polynominal Graph Fits
case k1DPol1:
case k1DPol2:
case k1DPol3:
case k1DPol4:
case k1DPol5:
case k1DPol6:
FitCoeff1DGraph(n, x, y, coeff, draw);
break;
// Spline Fits use TSpline3 Class
case k1DTSpline3:
GetCoeff1DTSpline3(n, x, y, coeff, draw);
break;
}
}
// Fitting Functions
void Spline::FitCoeff1DGraph(int n, double* x, double* y, float* coeff, bool draw) {
TGraph* gr = new TGraph(n, x, y);
double xmin = 99999.9;
double xmax = -99999.9;
for (int i = 0; i < n; i++) {
if (x[i] > xmax) xmax = x[i];
if (x[i] < xmin) xmin = x[i];
}
double xwidth = xmax - xmin;
xmin = xmin - xwidth * 0.01;
xmax = xmax + xwidth * 0.01;
// Create a new function for fitting.
TF1* func = new TF1("f1", this, -30, 30, this->GetNPar());
func->SetNpx(400);
func->FixParameter(0, 1.0); // Fix so 1.0 at nominal
// Run the actual spline fit
// StopTalking();
- std::cout << "Fixing TGraph" << std::endl;
+ // std::cout << "Fixing TGraph" << std::endl;
gr->Fit(func, "FMWQ");
// StartTalking();
for (int i = 0; i < this->GetNPar(); i++) {
coeff[i] = func->GetParameter(i);
}
if (draw) {
gr->Draw("APL");
gPad->Update();
gPad->SaveAs(("plot_test_" + fName + ".pdf").c_str());
std::cout << "Saving Graph" << std::endl;
sleep(3);
}
delete func;
delete gr;
}
// Spline extraction Functions
void Spline::GetCoeff1DTSpline3(int n, double* x, double* y, float* coeff, bool draw) {
TSpline3 temp_spline = TSpline3("temp_spline", x, y, n);
for (size_t i = 0; i < n; i++) {
double a, b, c, d, e;
temp_spline.GetCoeff(i, a, b, c, d, e);
coeff[i * 4] = y[i];
coeff[i * 4 + 1] = c;
coeff[i * 4 + 2] = d;
coeff[i * 4 + 3] = e;
}
if (draw) {
TGraph* gr = new TGraph(n, x, y);
- temp_spline.Draw("C");
- gr->Draw("PL");
+ temp_spline.Draw("CA");
+ gr->Draw("PL SAME");
gPad->Update();
gPad->SaveAs(("plot_test_" + fName + ".pdf").c_str());
// sleep(3);
delete gr;
}
return;
}
diff --git a/src/Splines/SplineMerger.cxx b/src/Splines/SplineMerger.cxx
new file mode 100644
index 0000000..a1ab3ca
--- /dev/null
+++ b/src/Splines/SplineMerger.cxx
@@ -0,0 +1,90 @@
+#include "SplineMerger.h"
+
+void SplineMerger::AddSplineSetFromFile(TFile* file){
+
+ TTree* tr = (TTree*) file->Get("spline_reader");
+ std::vector<std::string>* tempspline = 0;
+ std::vector<std::string>* temptype = 0;
+ std::vector<std::string>* tempform = 0;
+ std::vector<std::string>* temppoints = 0;
+
+
+ // Loop over all splines and add a TString in the TTree for its inputs
+ tr->SetBranchAddress("Spline", &tempspline);
+ tr->SetBranchAddress("Type", &temptype);
+ tr->SetBranchAddress("Form", &tempform);
+ tr->SetBranchAddress("Points", &temppoints);
+ tr->GetEntry(0);
+
+ // Copy over
+ for (size_t i = 0; i < tempspline->size(); i++){
+ fSpline.push_back(tempspline->at(i));
+ fType.push_back(temptype->at(i));
+ fForm.push_back(tempform->at(i));
+ fPoints.push_back(temppoints->at(i));
+ }
+
+ delete tr;
+
+
+ // Now Get the coefficients setup.
+ size_t index = fSplineTreeList.size();
+ fSplineTreeList.push_back( (TTree*) file->Get("spline_tree") );
+ // float temp[1000];
+ // fSplineAddressList.push_back(temp);
+ fSplineTreeList[index]->SetBranchAddress("SplineCoeff", fSplineAddressList[index]);
+
+}
+
+void SplineMerger::SetupSplineSet(){
+
+ // Define NCoEff
+ fNCoEff = 0;
+ for (size_t i = 0; i < fSplineSizeList.size(); i++){
+ fNCoEff += fSplineSizeList[i];
+ }
+
+ // Define Storer
+ fCoEffStorer = new float[fNCoEff];
+
+}
+
+
+void SplineMerger::GetEntry(int entry){
+ for (size_t i = 0; i < fSplineTreeList.size(); i++){
+ fSplineTreeList[i]->GetEntry(entry);
+ }
+}
+
+void SplineMerger::Write(std::string name){
+
+ // Create a TTree with each form and scan points in it.
+ TTree* tr = new TTree(name.c_str(), name.c_str());
+
+ // Loop over all splines and add a TString in the TTree for its inputs
+ tr->Branch("Spline", &fSpline);
+ tr->Branch("Type", &fType);
+ tr->Branch("Form", &fForm);
+ tr->Branch("Points", &fPoints);
+ tr->Fill();
+ tr->Write();
+
+ delete tr;
+}
+
+void SplineMerger::AddCoefficientsToTree(TTree* tree){
+ tree->Branch("SplineCoeff", fCoEffStorer, Form("SplineCoeff[%d]/F", fNCoEff));
+}
+
+void SplineMerger::FillMergedSplines(int entry){
+ GetEntry(entry);
+
+ size_t count = 0;
+ for (size_t i = 0; i < fSplineSizeList.size(); i++){
+ for (size_t j = 0; j < fSplineSizeList[i]; j++){
+ fCoEffStorer[count] = fSplineAddressList[i][j];
+ count++;
+ }
+ }
+}
+
diff --git a/src/Splines/SplineMerger.h b/src/Splines/SplineMerger.h
index 4be14a8..34dff97 100644
--- a/src/Splines/SplineMerger.h
+++ b/src/Splines/SplineMerger.h
@@ -1,34 +1,38 @@
#ifndef SPLINEMerger_H
#define SPLINEMerger_H
-/*
+
#include "FitWeight.h"
#include "Spline.h"
#include "FitParameters.h"
#include "SplineReader.h"
-class SplineMerger {
+class SplineMerger : public SplineReader {
public:
- SplineMerger(FitWeight* fw){
- fDrawSplines = FitPar::Config().GetParB("drawsplines");
- };
+ SplineMerger(){};
~SplineMerger(){};
- void LoadSplinesForMerge(SplineReader* spl);
+ void AddSplineSetFromFile(TFile* file);
void SetupSplineSet();
+
void Write(std::string name);
void AddCoefficientsToTree(TTree* tree);
+ void GetEntry(int entry);
+
+ void FillMergedSplines(int entry);
+ float* fCoEffStorer;
int fNCoEff;
- double* fCoEffStorer;
+
+ std::vector< float[1000] > fSplineAddressList;
+ std::vector< int > fSplineSizeList;
+ std::vector< TTree* > fSplineTreeList;
std::vector< std::vector<double> > fParVect;
std::vector< int > fSetIndex;
std::vector< double > fWeightList;
std::vector< double > fValList;
- FitWeight* fRW;
- bool fDrawSplines;
};
-*/
+
#endif
diff --git a/src/Splines/SplineReader.cxx b/src/Splines/SplineReader.cxx
index cf35588..f0d0cbc 100644
--- a/src/Splines/SplineReader.cxx
+++ b/src/Splines/SplineReader.cxx
@@ -1,126 +1,125 @@
#include "SplineReader.h"
// Spline reader should have access to every spline.
// Should know when reconfigure is called what its limits are and adjust accordingly.
// Then should pass it the index required in the stack as &xvals[index], and &pars[parindex]
void SplineReader::AddSpline(nuiskey splinekey) {
// Read info from key
std::string splname = splinekey.GetS("name");
std::string type = splinekey.GetS("type");
std::string form = splinekey.GetS("form");
std::string points = splinekey.GetS("points");
std::string extrapolate = splinekey.GetS("extrapolate");
- std::cout << "AddSpline " << splname << " " << type << " " << form << " " << points << std::endl;
+ // std::cout << "AddSpline " << splname << " " << type << " " << form << " " << points << std::endl;
// Add the spline to the list of all forms
fAllSplines.push_back( Spline(splname, form, GeneralUtils::ParseToDbl(points, ",")) );
fSpline.push_back(splname);
fType.push_back(type);
fForm.push_back(form);
fPoints.push_back(points);
}
void SplineReader::Read(TTree* tr) {
std::vector<std::string>* tempspline = 0;
std::vector<std::string>* temptype = 0;
std::vector<std::string>* tempform = 0;
std::vector<std::string>* temppoints = 0;
// Loop over all splines and add a TString in the TTree for its inputs
tr->SetBranchAddress("Spline", &tempspline);
tr->SetBranchAddress("Type", &temptype);
tr->SetBranchAddress("Form", &tempform);
tr->SetBranchAddress("Points", &temppoints);
tr->GetEntry(0);
// Copy over
- std::cout << "TempSpline Size = " <<tempspline->size() << std::endl;
for (size_t i = 0; i < tempspline->size(); i++){
fSpline.push_back(tempspline->at(i));
fType.push_back(temptype->at(i));
fForm.push_back(tempform->at(i));
fPoints.push_back(temppoints->at(i));
- std::cout << " Found Spline In TTree " << tempspline->at(i) << std::endl;
}
- /*
- delete tempspline;
- delete tempform;
- delete temppoints;
- delete temptype;
-*/
+
// Loop through and add splines from read type.
for (size_t i = 0; i < fSpline.size(); i++) {
- std::cout << "Adding New Spline " << fSpline[i] << " " << fForm[i] << " " << fPoints[i] << std::endl;
+ LOG(SAM) << "Registering Input Spline " << fSpline[i] << " " << fForm[i] << " " << fPoints[i] << std::endl;
fAllSplines.push_back( Spline(fSpline[i], fForm[i],
GeneralUtils::ParseToDbl(fPoints[i], ",")) );
}
}
void SplineReader::Reconfigure(std::map< std::string, double >& vals) {
for (std::map<std::string, double>::iterator iter = vals.begin();
iter != vals.end(); iter++){
// std::cout << " Found " << iter->first << " in map handed to reader." << std::endl;
for (size_t i = 0; i < fSpline.size(); i++){
// std::cout << " Comparing it to : " << fSpline[i] << std::endl;
if (!fSpline[i].compare(iter->first.c_str())){
// std::cout << "Reconfiguring Value inside Reader to be " << fSpline[i] << " " << iter->second << std::endl;
// sleep(1);
fAllSplines[i].Reconfigure(iter->second);
}
}
}
fNeedsReconfigure = false;
}
bool SplineReader::NeedsReconfigure() {
return fNeedsReconfigure;
}
void SplineReader::SetNeedsReconfigure(bool val) {
fNeedsReconfigure = val;
}
double SplineReader::CalcWeight(float* coeffs) {
int off = 0;
double rw_weight = 1.0;
for (size_t i = 0; i < fAllSplines.size(); i++) {
// std::cout << "Evaluating spline " << fAllSplines[i].GetName() << " at coeff offset " << off << "(" << coeffs << ")" << std::endl;
// for (int j = 0; j < fAllSplines[i].GetNPar(); j++){
// std::cout << "Coeff " << j+off << " " << coeffs[off+j] << std::endl;
// }
double w = fAllSplines[i].DoEval( &coeffs[off] );
rw_weight *= w;
// std::cout << "Spline RW Weight = " << rw_weight << " " << w << std::endl;
off += fAllSplines[i].GetNPar();
}
if (rw_weight <= 0.0) rw_weight = 1.0;
// std::cout << "Returning spline rw_weight = " << rw_weight << std::endl;
return rw_weight;
}
-
+int SplineReader::GetNPar(){
+ int n = 0;
+ for (int i = 0; i < fAllSplines.size(); i++){
+ n += fAllSplines[i].GetNPar();
+ }
+ return n;
+}
diff --git a/src/Splines/SplineReader.h b/src/Splines/SplineReader.h
index 2a6e98d..fe1706c 100644
--- a/src/Splines/SplineReader.h
+++ b/src/Splines/SplineReader.h
@@ -1,37 +1,38 @@
#ifndef SPLINEREADER_H
#define SPLINEREADER_H
// #include "FitWeight.h"
#include "Spline.h"
#include "TTree.h"
// #include "GeneralUtils.h"
class SplineReader {
public:
SplineReader() {};
~SplineReader() {};
void AddSpline(nuiskey splkey);
void Read(TTree* tr);
void Reconfigure(std::map< std::string, double >& vals);
bool NeedsReconfigure();
void SetNeedsReconfigure(bool val = true);
+ int GetNPar();
double CalcWeight(float* coeffs);
std::vector<Spline> fAllSplines;
std::vector<std::string> fSpline;
std::vector<std::string> fType;
std::vector<std::string> fForm;
std::vector<std::string> fPoints;
std::vector<double> fDialValues;
std::vector<double> fParValues;
bool fNeedsReconfigure;
};
#endif
diff --git a/src/Splines/SplineWriter.cxx b/src/Splines/SplineWriter.cxx
index e61a23a..e423ebc 100644
--- a/src/Splines/SplineWriter.cxx
+++ b/src/Splines/SplineWriter.cxx
@@ -1,250 +1,262 @@
#include "SplineWriter.h"
// Spline reader should have access to every spline.
// Should know when reconfigure is called what its limits are and adjust accordingly.
// Then should pass it the index required in the stack as &xvals[index], and &pars[parindex]
void SplineWriter::Write(std::string name) {
// Create a TTree with each form and scan points in it.
TTree* tr = new TTree(name.c_str(), name.c_str());
// Loop over all splines and add a TString in the TTree for its inputs
tr->Branch("Spline", &fSpline);
tr->Branch("Type", &fType);
tr->Branch("Form", &fForm);
tr->Branch("Points", &fPoints);
tr->Fill();
tr->Write();
delete tr;
}
void SplineWriter::AddCoefficientsToTree(TTree* tr) {
// Add only the fitted spline coefficients to the ttree
std::cout << "Saving Spline Coeff to TTree = " << Form("SplineCoeff[%d]/F", fNCoEff) << std::endl;
// sleep(1);
tr->Branch("SplineCoeff", fCoEffStorer, Form("SplineCoeff[%d]/F", fNCoEff));
}
void SplineWriter::SetupSplineSet() {
std::cout << "Setting up spline set" << std::endl;
// Create the coefficients double*
fNCoEff = 0;
for (int i = 0; i < fAllSplines.size(); i++) {
fNCoEff += fAllSplines[i].GetNPar();
}
fCoEffStorer = new float[fNCoEff];
std::cout << "NCoeff = " << fNCoEff << std::endl;
// Calculate the grid of parsets
// Setup the list of parameter coefficients.
std::vector<double> nomvals = fRW->GetDialValues();
fParVect.clear();
fSetIndex.clear();
fParVect.push_back(nomvals);
fSetIndex.push_back(0);
fWeightList.push_back(1.0);
fValList.push_back(0.0);
// Loop over all splines.
for (int i = 0; i < fAllSplines.size(); i++) {
// For each dial loop over all points within a given position
std::vector<double> newvals = nomvals;
// Create a new set of nom vals for that dial set, and attribute it to the spline index.
int pos = fRW->GetDialPos(fSpline[i]);
// Split Points
std::vector<double> vals = GeneralUtils::ParseToDbl(fPoints[i], ",");
for (int j = 0; j < vals.size(); j++) {
newvals[pos] = vals[j];
fParVect.push_back(newvals);
fValList.push_back(vals[j]);
fWeightList.push_back(1.0);
fSetIndex.push_back(i+1);
}
}
// Print out the parameter set
std::cout << "Parset | Index | Pars --- " << std::endl;
for (int i = 0; i < fSetIndex.size(); i++) {
std::cout << " Set " << i << ". | " << fSetIndex[i] << " | ";
for (int j = 0; j < fParVect[i].size(); j++) {
std::cout << " " << fParVect[i][j];
}
std::cout << std::endl;
}
}
void SplineWriter::FitSplinesForEvent(FitEvent* event) {
fRW->SetAllDials(&fParVect[0][0], fParVect[0].size());
double nomweight = fRW->CalcWeight(event);
+ event->RWWeight = nomweight;
+ if (fDrawSplines){
+ std::cout << "Nominal Spline Weight = " << nomweight << std::endl;
+ }
+
// Loop over parameter sets
for (int i = 1; i < fParVect.size(); i++) {
// Update FRW
fRW->SetAllDials(&fParVect[i][0], fParVect[i].size());
// Calculate a weight for event
double weight = fRW->CalcWeight(event);
if (weight >= 0.0 and weight < 200){
// Fill Weight Set
fWeightList[i] = weight/nomweight;
- // std::cout << "Calculating values from weight set " << i << " " << fParVect[i][0] << " = " << weight << std::endl;
+ if (fDrawSplines) std::cout << "Calculating values from weight set " << i << " " << fParVect[i][0] << " = " << weight << " " << weight/nomweight << std::endl;
} else {
fWeightList[i] = 1.0;
}
}
// Loop over splines
int count = 0;
int coeffcount = 0;
for (int i = 0; i < fAllSplines.size(); i++) {
// Store X/Y Vals
std::vector<double> xvals;
std::vector<double> yvals;
bool hasresponse = false;
int npar = (fAllSplines[i]).GetNPar();
for (int j = 0; j < fSetIndex.size(); j++) {
if (fSetIndex[j] != i+1) continue;
xvals.push_back(fValList[j]);
yvals.push_back(fWeightList[j] - 0.0);
if (fWeightList[j] != 1.0) hasresponse = true;
}
// Make a new graph and fit coeff if response
if (hasresponse) {
(fAllSplines[i]).FitCoeff(int(xvals.size()), &xvals[0], &yvals[0], &fCoEffStorer[coeffcount], fDrawSplines);
} else {
// std::cout << "Spline " << fSpline[i] << " has no response. " << std::endl;
for (int i = 0; i < npar; i++) {
fCoEffStorer[coeffcount + i] = 0.0;
}
}
// Make a new plot
if (fDrawSplines and hasresponse) {
TGraph* gr = new TGraph(xvals.size(), &xvals[0], &yvals[0]);
// Get XMax Min
int n = xvals.size();
double xmin = 99999.9;
double xmax = -99999.9;
- for (int i = 0; i < n; i++) {
- if (xvals[i] > xmax) xmax = xvals[i];
- if (xvals[i] < xmin) xmin = xvals[i];
+ for (int j = 0; j < n; j++) {
+ if (xvals[j] > xmax) xmax = xvals[j];
+ if (xvals[j] < xmin) xmin = xvals[j];
}
double xwidth = xmax - xmin;
xmin = xmin - 0.01;
xmax = xmax + 0.01;
TH1D* hist = new TH1D("temp", "temp", 100, xmin, xmax);
for (int k = 0; k < 100; k++) {
double xtemp = hist->GetXaxis()->GetBinCenter(k + 1);
- double ytemp = fAllSplines[i].DoEval(&xtemp, &fCoEffStorer[coeffcount]);
+ fAllSplines[i].Reconfigure(xtemp);
+ double ytemp = fAllSplines[i].DoEval(&fCoEffStorer[coeffcount]);
hist->SetBinContent(k + 1, ytemp);
// std::cout << "Set Temp " << k << " " << ytemp << std::endl;
}
// gr->Draw("APL");
hist->SetLineColor(kRed);
hist->Draw("HIST C");
gr->SetMarkerStyle(20);
gr->Draw("P SAME");
gPad->Update();
gPad->SaveAs(("F1eval_" + fSpline[i] + ".pdf").c_str());
delete gr;
// delete f1;
}
coeffcount += npar;
}
// Check overrides
if (fDrawSplines) {
coeffcount = 0;
for (int i = 0; i < fAllSplines.size(); i++) {
// Store X/Y Vals
std::vector<double> xvals;
std::vector<double> yvals;
bool hasresponse = false;
int npar = (fAllSplines[i]).GetNPar();
for (int j = 0; j < fSetIndex.size(); j++) {
if (fSetIndex[j] != i+1) continue;
xvals.push_back(fValList[j]);
yvals.push_back(fWeightList[j] - 0.0);
if (fWeightList[j] != 1.0) hasresponse = true;
}
if (hasresponse) {
TGraph* gr = new TGraph(xvals.size(), &xvals[0], &yvals[0]);
// Get XMax Min
int n = xvals.size();
double xmin = 99999.9;
double xmax = -99999.9;
- for (int i = 0; i < n; i++) {
- if (xvals[i] > xmax) xmax = xvals[i];
- if (xvals[i] < xmin) xmin = xvals[i];
+ for (int j = 0; j < n; j++) {
+ if (xvals[j] > xmax) xmax = xvals[j];
+ if (xvals[j] < xmin) xmin = xvals[j];
}
// double xwidth = xmax - xmin;
// xmin = xmin - xwidth * 0.01;
//xmax = xmax + xwidth * 0.01;
TH1D* hist = new TH1D("temp", "temp", 100, xmin, xmax);
+ // for (int k = 0; k < fAllSplines[i].GetNPar(); k++){
+ // std::cout << fAllSplines[i].GetName() << " Coeff " << k << " = " << fCoEffStorer[coeffcount + k] << std::endl;
+ // }
+
for (int k = 0; k < 100; k++) {
double xtemp = hist->GetXaxis()->GetBinCenter(k + 1);
- double ytemp = fAllSplines[i].DoEval(&xtemp, &fCoEffStorer[coeffcount]);
+ fAllSplines[i].Reconfigure(xtemp);
+ double ytemp = fAllSplines[i].DoEval(&fCoEffStorer[coeffcount]);
hist->SetBinContent(k + 1, ytemp);
- // std::cout << "Set Temp " << k << " " << ytemp << std::endl;
+
+ // std::cout << "Set Temp " << k << " " << ytemp << " First Coeff = " << fCoEffStorer[coeffcount] << std::endl;
}
// gr->Draw("APL");
hist->SetLineColor(kRed);
hist->Draw("HIST C");
hist->SetTitle("Spline Response");
hist->GetYaxis()->SetRangeUser(0.0,3.0);
// gStyle->SetOptStat(0);
hist->SetStats(0);
gr->SetMarkerStyle(20);
gr->SetTitle("True Weight Points");
gr->Draw("P SAME");
gPad->BuildLegend();
gPad->Update();
hist->SetTitle(fSpline[i].c_str());
hist->GetXaxis()->SetTitle("Dial Variation");
hist->GetYaxis()->SetTitle("Event Weight");
gPad->Update();
gPad->SaveAs(("F2_eval_" + fSpline[i] + ".pdf").c_str());
delete gr;
std::cout << "Saved hist for " << fSpline[i] << std::endl;
sleep(5);
// delete f1;
}
coeffcount += npar;
}
}
}
diff --git a/src/Utils/FitLogger.cxx b/src/Utils/FitLogger.cxx
index a109f4e..25f7fd4 100644
--- a/src/Utils/FitLogger.cxx
+++ b/src/Utils/FitLogger.cxx
@@ -1,196 +1,218 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitLogger.h"
#include <fcntl.h>
#include <unistd.h>
namespace FitPar{
unsigned int log_verb = 0; //!< Current VERBOSITY
unsigned int err_verb = 0; //!< Current ERROR VERBOSITY
bool use_colors = true; //!< Use BASH Terminal Colors Flag
bool super_rainbow_mode = true; //!< For when fitting gets boring.
unsigned int super_rainbow_mode_colour = 0;
- bool showtrace = true;
+ bool showtrace = false;
// For redirecting various print outs
std::streambuf *default_cout = std::cout.rdbuf();
std::streambuf *default_cerr = std::cerr.rdbuf();
std::ofstream redirect_stream("/dev/null");
int silentfd = open("/dev/null",O_WRONLY);
int savedstdoutfd = dup(fileno(stdout));
int savedstderrfd = dup(fileno(stderr));
+
+ int nloggercalls = 0;
+ int timelastlog = 0;
+
}
std::ostream* logStream(&std::cout);
std::ostream* errStream(&std::cerr);
std::ofstream nullStream;
//******************************************
void LOG_VERB(std::string verb){
//******************************************
if (!verb.compare("DEB")) FitPar::log_verb=-1;
else if (!verb.compare("QUIET")) FitPar::log_verb=0;
else if (!verb.compare("FIT")) FitPar::log_verb=1;
else if (!verb.compare("MIN")) FitPar::log_verb=2;
else if (!verb.compare("SAM")) FitPar::log_verb=3;
else if (!verb.compare("REC")) FitPar::log_verb=4;
else if (!verb.compare("SIG")) FitPar::log_verb=5;
else if (!verb.compare("EVT")) FitPar::log_verb=6;
else FitPar::log_verb = GeneralUtils::StrToInt(verb);
std::cout << "Set logging verbosity to : " << FitPar::log_verb << std::endl;
return;
}
//******************************************
void ERR_VERB(std::string verb){
//******************************************
std::cout << "Setting ERROR VERB" << std::endl;
if (!verb.compare("ERRQUIET")) FitPar::err_verb=0;
else if (!verb.compare("FTL")) FitPar::err_verb=1;
else if (!verb.compare("WRN")) FitPar::err_verb=2;
else FitPar::err_verb = GeneralUtils::StrToInt(verb);
std::cout << "Set error verbosity to : " << FitPar::err_verb << std::endl;
return;
}
//******************************************
bool LOG_LEVEL(int level){
//******************************************
if (FitPar::log_verb == (unsigned int) DEB){
return true;
}
if (FitPar::log_verb < (unsigned int) level){
return false;
}
return true;
}
void SET_TRACE(bool val){
FitPar::showtrace = val;
}
// std::ostream& LOG(FIT){
// std::cout << "[ NUISANCE ]: ";
// return *logStream;
// }
//******************************************
std::ostream& _LOG(int level, const char* filename, const char* func, int line)
//******************************************
{
+ if (FitPar::nloggercalls > 1000){
+ std::cout << "[WARNING] : Large number of logger calls being piped to LOGGER." << std::endl;
+ std::cout << "This is super inefficient. : " << FitPar::nloggercalls << std::endl;
+ std::cout << "Occuring at " << filename << "::" << func << "[l. " << line << "]" << std::endl;
+ FitPar::nloggercalls = -1;
+ }
+
+ if (abs(time(NULL) - FitPar::timelastlog) > 10){
+ FitPar::timelastlog = time(NULL);
+ FitPar::nloggercalls = 0;
+ }
+
+ if (FitPar::nloggercalls != -1){
+ FitPar::nloggercalls++;
+ }
+
+
+
if (FitPar::log_verb < (unsigned int)level &&
FitPar::log_verb != (unsigned int)DEB){
return nullStream;
} else {
if (FitPar::showtrace){
- std::cout << filename << "::" << func << "[l" << line << "] : ";
+ std::cout << filename << "::" << func << "[l. " << line << "] : ";
}
if (FitPar::super_rainbow_mode and FitPar::use_colors){
switch(FitPar::super_rainbow_mode_colour){
case 1: std::cout<<RED;
case 2: std::cout<<GREEN;
case 3: std::cout<<YELLOW;
case 4: std::cout<<BLUE;
case 5: std::cout<<MAGENTA;
case 6: std::cout<<CYAN;
default: FitPar::super_rainbow_mode_colour = 0;
}
FitPar::super_rainbow_mode_colour++;
}
if (FitPar::use_colors){
switch(level){
case FIT: std::cout << BOLDGREEN; break;
case MIN: std::cout << BOLDBLUE; break;
case SAM: std::cout << MAGENTA; break;
case REC: std::cout << BLUE; break;
case SIG: std::cout << GREEN; break;
case DEB: std::cout << CYAN; break;
default: break;
}
}
switch(level){
case FIT: std::cout << "[LOG Fitter]: "; break;
case MIN: std::cout << "[LOG Minmzr]: "; break;
case SAM: std::cout << "[LOG Sample]: - "; break;
case REC: std::cout << "[LOG Reconf]: -- "; break;
case SIG: std::cout << "[LOG Signal]: --- "; break;
case EVT: std::cout << "[LOG Event ]: ---- "; break;
case DEB: std::cout << "[LOG DEBUG ]: "; break;
default: std::cout << "Log : "; break;
}
if (FitPar::use_colors or FitPar::super_rainbow_mode) std::cout << RESET;
return *logStream;
}
}
//******************************************
std::ostream& ERR(int level)
//******************************************
{
if (FitPar::use_colors) std::cerr << RED;
switch(level){
case FTL: std::cerr << "[ERR FATAL ]: "; break;
case WRN: std::cerr << "[ERR WARN ] : "; break;
}
if (FitPar::use_colors) std::cerr << RESET;
return *errStream;
}
void StopTalking(){
// Only redirect if we're not debugging
if (FitPar::log_verb == (unsigned int)DEB) return;
std::cout.rdbuf(FitPar::redirect_stream.rdbuf());
std::cerr.rdbuf(FitPar::redirect_stream.rdbuf());
shhnuisancepythiaitokay_();
fflush(stdout);
fflush(stderr);
dup2(FitPar::silentfd, fileno(stdout));
dup2(FitPar::silentfd, fileno(stderr));
}
void StartTalking(){
std::cout.rdbuf(FitPar::default_cout);
std::cerr.rdbuf(FitPar::default_cerr);
canihaznuisancepythia_();
fflush(stdout);
fflush(stderr);
dup2(FitPar::savedstdoutfd, fileno(stdout));
dup2(FitPar::savedstderrfd, fileno(stderr));
}
diff --git a/src/Utils/FitLogger.h b/src/Utils/FitLogger.h
index ffefb1b..9a5a6d9 100644
--- a/src/Utils/FitLogger.h
+++ b/src/Utils/FitLogger.h
@@ -1,127 +1,128 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef FITLOGGER_HPP
#define FITLOGGER_HPP
// We still need a forward declaration of 'ostream' in order to
// swallow templated manipulators such as 'endl'.
#include <iosfwd>
#include <iostream>
#include <fstream>
#include <sstream>
#include "GeneralUtils.h"
#define RESET "\033[0m"
#define BLACK "\033[30m" /* Black */
#define RED "\033[31m" /* Red */
#define GREEN "\033[32m" /* Green */
#define YELLOW "\033[33m" /* Yellow */
#define BLUE "\033[34m" /* Blue */
#define MAGENTA "\033[35m" /* Magenta */
#define CYAN "\033[36m" /* Cyan */
#define WHITE "\033[37m" /* White */
#define BOLDBLACK "\033[1m\033[30m" /* Bold Black */
#define BOLDRED "\033[1m\033[31m" /* Bold Red */
#define BOLDGREEN "\033[1m\033[32m" /* Bold Green */
#define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */
#define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */
#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */
#define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */
#define BOLDWHITE "\033[1m\033[37m" /* Bold White */
/*!
* \addtogroup FitBase
* @{
*/
namespace FitPar{
extern unsigned int log_verb; //!< Current VERBOSITY
extern unsigned int err_verb; //!< Current ERROR VERBOSITY
extern bool use_colors; //!< Use BASH Terminal Colors Flag
extern bool super_rainbow_mode; //!< For when fitting gets boring.
extern unsigned int super_rainbow_mode_colour;
extern bool showtrace; // Quick Tracing for debugging
-
+ extern int nloggercalls;
+ extern int timelastlog;
extern std::streambuf *default_cout; //!< Where the STDOUT stream is currently directed
extern std::streambuf *default_cerr; //!< Where the STDERR stream is currently directed
extern std::ofstream redirect_stream; //!< Where should unwanted messages be thrown
}
extern std::ostream* logStream;
extern std::ostream* errStream;
extern std::ofstream nullStream;
/// Fitter VERBOSITY Enumerations
/// These go through the different depths of the fitter.
///
/// 0 QUIET - Little output.
/// 1 FIT - Top Level Minimizer Status
/// 2 MIN - Output from the FCN Minimizer Functions
/// 3 SAM - Output from each of the samples during setup etc
/// 4 REC - Output during each reconfigure. Percentage progress etc.
/// 5 SIG - Output during every signal event that is found.
/// 6 EVT - Output during every event.
/// -1 DEB - Will print only debugging info wherever a LOG(DEB) statement was made
enum log_levels { DEB=-1, QUIET, FIT, MIN, SAM, REC, SIG, EVT };
/// Fitter ERROR VERBOSITY Enumerations
///
/// 0 QUIET - No Error Output
/// 1 FTL - Show errors only if fatal
/// 2 WRN - Show Warning messages
enum err_levels { ERRQUIET=0, FTL, WRN };
bool LOG_LEVEL(int level);
//! Set LOG VERBOSITY from a string
void LOG_VERB(std::string verb);
inline void LOG_VERB(int verb){ FitPar::log_verb = verb; };
void SET_TRACE(bool val);
//! Set ERROR VERBOSITY from a string
void ERR_VERB(std::string verb);
inline void ERR_VERB(int verb){ FitPar::err_verb = verb; };
#define __FILENAME__ (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
/// Logging Function. Use as a string stream. e.g. LOG(SAM) << "This sample is dope." << std::endl;
std::ostream& _LOG(int level, const char* filename, const char* funct, int line);
#define LOG(level) _LOG(level, __FILENAME__, __FUNCTION__, __LINE__)
//! Error Function. Use as a string stream. e.g. ERR(FTL) << "The fit is completely buggered." << std::endl;
std::ostream& ERR(int level);
void StopTalking();
void StartTalking();
extern "C" {
void shhnuisancepythiaitokay_(void);
void canihaznuisancepythia_(void);
}
/*! @} */
#endif // FILELOGGER_HPP

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 3:03 PM (38 m, 47 s)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486819
Default Alt Text
(164 KB)

Event Timeline