diff --git a/parameters/config.xml b/parameters/config.xml
index 0a23a39..6c7ecfa 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,157 +1,158 @@
<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='3'/>
-<config VERBOSITY='3'/>
+<config verbosity='5'/>
+<config VERBOSITY='5'/>
<!-- # ERROR goes from -->
<!-- # 0 NONE -->
<!-- # 1 FATAL -->
<!-- # 2 WARN -->
<config ERROR='2'/>
<config cores='2' />
+<config spline_test_throws='50' />
<!-- # 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'/>
<!-- # 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 65f935c..3ba5a28 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,942 +1,962 @@
#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;
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;
+ LOG(MIN) << "Loading Samples : " << samplekeys.size() << 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 = "";
+ LOG(MIN) << "Loading Sample : " << samplename << std::endl;
+
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();
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() {
//***************************************************
// 'Slow' Event Manager Reconfigure
LOG(REC) << "Event Manager Reconfigure" << std::endl;
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// If we are siving signal, reset all containers.
bool savesignal = (FitPar::Config().GetParB("SignalReconfigures"));
if (savesignal) {
// Reset all of our event signal vectors
fSignalEventBoxes.clear();
fSignalEventFlags.clear();
fSampleSignalFlags.clear();
fSignalEventSplines.clear();
}
// Make sure we have a list of inputs
if (fInputList.empty()) {
fInputList = GetInputList();
fSubSampleList = GetSubSampleList();
}
// If all inputs are splines make sure the readers are told
// they need to be reconfigured.
std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell reader in each BaseEvent it needs a Reconfigure next weight calc.
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) {
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
-
- }
-
+ }
// MAIN INPUT LOOP ====================
int fillcount = 0;
int inputcount = 0;
inp_iter = fInputList.begin();
// Loop over each input in manager
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Get event information
FitEvent* curevent = curinput->FirstNuisanceEvent();
int i = 0;
int nevents = curinput->GetNEvents();
int countwidth = nevents / 20;
// Start event loop iterating until we get a NULL pointer.
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;
}
}
// 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;
// Start measurement iterator
size_t measitercount = 0;
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
// Loop over all subsamples (sub in JointMeas)
for (; meas_iter != fSubSampleList.end(); meas_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// Compare input pointers, to current input, skip if not.
// Pointer tells us if it matches without doing ID checks.
if (curinput != curmeas->GetInput()) {
if (savesignal) {
// Set bit to 0 as definitely not signal
signalbitset[measitercount] = 0;
}
// Count up what measurement we are on.
measitercount++;
// Skip sample as input not signal.
continue;
}
// Fill events for matching inputs.
MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);
bool signal = curmeas->isSignal(curevent);
curmeas->SetSignal(signal);
curmeas->FillHistograms(rwweight);
// If its Signal tally up fills
if (signal) {
fillcount++;
}
// If we are saving signal/splines fill the bitset
if (savesignal) {
signalbitset[measitercount] = signal;
}
// If signal save a clone of the event box for use later.
if (savesignal and signal) {
foundsignal = true;
signalboxes.push_back( box->CloneSignalBox() );
}
// Keep track of Measurement we are on.
measitercount++;
}
// Once we've filled the measurements, if saving signal
// push back if any sample flagged this event as signal
if (savesignal) {
fSignalEventFlags.push_back(foundsignal);
}
// Save the vector of signal boxes for this event
if (savesignal and foundsignal) {
fSignalEventBoxes.push_back(signalboxes);
fSampleSignalFlags.push_back(signalbitset);
}
// If all inputs are splines we can save the spline coefficients
// for fast in memory reconfigures later.
- if (fIsAllSplines) {
+ if (fIsAllSplines and savesignal and foundsignal) {
// Make temp vector to push back with
std::vector<float> coeff;
for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) {
coeff.push_back( curevent->fSplineCoeff[l] );
}
// Push back to signal event splines. Kept in sync with fSignalEventBoxes size.
+ int splinecount = fSignalEventSplines.size();
fSignalEventSplines.push_back(coeff);
+
+ // if (splinecount % 1000 == 0) {
+ // std::cout << "Pushed Back Coeff " << splinecount << " : ";
+ // for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) {
+ // std::cout << " " << fSignalEventSplines[splinecount][l];
+ // }
+ // std::cout << std::endl;
+ // }
+
}
// Clean up vectors once done with this event
signalboxes.clear();
signalbitset.clear();
// Iterate to the next event.
curevent = curinput->NextNuisanceEvent();
i++;
}
// Keep track of what input we are on.
inputcount++;
}
// End of Event Loop ===============================
// Now event loop is finished loop over all Measurements
// Converting Binned events to XSec Distributions
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Print out statements on approximate memory usage for profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
if (savesignal) {
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) << " -> Saved " << fillcount << " " << fSignalEventSplines.size() << " spline sets into memory. (~" << splmem << " MB)" << std::endl;
}
}
LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl;
// End of reconfigure
return;
};
//***************************************************
void JointFCN::ReconfigureFastUsingManager() {
//***************************************************
// Get Start time for profilling
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 not do a full reconfigure.
if (fSignalEventFlags.empty()) {
ReconfigureUsingManager();
return;
}
// Setup fast vector 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();
+ int splinecount = 0;
// Setup stuff for logging
int fillcount = 0;
int nevents = fSignalEventFlags.size();
int countwidth = nevents / 5;
// If All Splines tell splines they need a reconfigure.
- std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
- if (fIsAllSplines){
+ std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
+ if (fIsAllSplines) {
+ LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl;
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell each fSplineRead in BaseFitEvent to reconf next weight calc
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
// Start of Fast Event Loop ============================
// Start input iterators
inp_iter = fInputList.begin();
// Loop over all inputs
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Get Only base event with RW info.
BaseFitEvt* curevent = curinput->FirstBaseEvent();
int i = 0;
// Iterate over all events and signal flags.
while (curevent != 0 and (inpsig_iter != fSignalEventFlags.end())) {
// Logging
if (LOG_LEVEL(REC)) {
if (i % countwidth == 0) {
LOG(REC) << "Processed " << i << " signal events.";
}
}
// If event has not been flagged as signal in vector skip it.
if (!(*inpsig_iter)) {
- inpsig_iter++;
- i++;
-
if (LOG_LEVEL(REC)) {
if (i % countwidth == 0) {
std::cout << std::endl;
}
}
+ inpsig_iter++;
+ i++;
+
continue;
}
// Get The Base Event
curevent = curinput->GetBaseEvent(i);
// End iterator if NULL pointer at this entry..
if (!curevent) break;
- // Setup signal splines pointer for this event if required.
+ // Setup signal splines pointer for this event 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;
// Weight Logging
if (LOG_LEVEL(REC)) {
if (i % countwidth == 0) {
std::cout << " W = " << rwweight << std::endl;
}
}
// Get iterators 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();
// Loop over all sub measurements.
for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// If event flagged as signal for this sample fill from the box.
if (*subsamsig_iter) {
curmeas->SetSignal(true);
curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);
// Move onto next box if there is one.
subbox_iter++;
fillcount++;
}
}
+ // if (splinecount % 1000 == 0) {
+ // std::cout << "Read Back Coeff " << splinecount << " : ";
+ // for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) {
+ // std::cout << " " << fSignalEventSplines[splinecount][l];
+ // }
+ // std::cout << std::endl;
+ // }
+
// Iterate over the main signal event containers.
samsig_iter++;
box_iter++;
spline_iter++;
-
+ splinecount++;
// iterate to next signal event
inpsig_iter++;
i++;
}
}
// End of Fast Event Loop ===================
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Print some reconfigure profiling.
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/FitBase/EventManager.cxx b/src/FitBase/EventManager.cxx
index 57673b2..b159527 100644
--- a/src/FitBase/EventManager.cxx
+++ b/src/FitBase/EventManager.cxx
@@ -1,137 +1,142 @@
// 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 "EventManager.h"
EventManager* EventManager::m_evtmgrInstance = NULL;
EventManager& EventManager::Get(void) {
if (!m_evtmgrInstance) {
m_evtmgrInstance = new EventManager;
}
return *m_evtmgrInstance;
}
FitWeight* EventManager::GetRW() { return fRW; };
InputHandlerBase* EventManager::GetInput(int infile) { return finputs[infile]; };
FitEvent* EventManager::GetEvent(int infile, int i) {
/*
finputs[infile]->ReadEvent(i);
FitEvent* evtpt = finputs[infile]->GetEventPointer();
// If we don't need a full reweight
if (!frwneeded[infile][i]) {
evtpt->Weight = calc_rw[infile][i];
// if we do need a full reweight
} else {
evtpt->RWWeight = fRW->CalcWeight(evtpt);
evtpt->Weight = evtpt->RWWeight * evtpt->InputWeight;
calc_rw[infile][i] = evtpt->Weight;
frwneeded[infile][i] = false;
}
*/
return NULL;
}
double EventManager::GetEventWeight(int infile, int i) {
/*
if (!frwneeded[infile][i]) {
return calc_rw[infile][i];
}
finputs[infile]->GetTreeEntry(i);
FitEvent* evtpt = finputs[infile]->GetEventPointer();
double Weight = fRW->CalcWeight(evtpt) * evtpt->InputWeight;
calc_rw[infile][i] = Weight;
frwneeded[infile][i] = false;
*/
double Weight = 1.0;
return Weight;
}
std::map<int, InputHandlerBase*> EventManager::GetInputs() { return finputs; }
InputHandlerBase* EventManager::AddInput(std::string handle, std::string infile) {
// Expect INPUTTYPE:FileLocation(s)
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(infile, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << infile
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inpType =
InputUtils::ParseInputType(file_descriptor[0]);
int id = GetInputID(file_descriptor[1]);
if ((uint)id != fid.size()) {
LOG(SAM) << "Event manager already contains " << file_descriptor[1]
<< std::endl;
return finputs[id];
}
fid[file_descriptor[1]] = id;
finputs[id] = InputUtils::CreateInputHandler(handle, inpType, file_descriptor[1]);
frwneeded[id] = std::vector<bool>(finputs[id]->GetNEvents(), true);
calc_rw[id] = std::vector<double>(finputs[id]->GetNEvents(), 0.0);
LOG(SAM) << "Registered " << handle << " with EventManager." << std::endl;
return finputs[id];
}
// Reset the weight flags
// Should be called for every succesful event loop
void EventManager::ResetWeightFlags() {
/*
// Loop over the inpts
for (std::map<int, InputHandler*>::iterator iter = finputs.begin();
iter != finputs.end(); iter++) {
int id = iter->first;
frwneeded[id].clear();
// Reset so that all events need the reweight
frwneeded[id] = std::vector<bool>(finputs[id]->GetNEvents(), true);
}
*/
}
EventManager::EventManager() {
fRW = new FitWeight("FitWeight");
finputs.clear();
};
EventManager::~EventManager() {
delete fRW;
finputs.clear();
};
int EventManager::GetInputID(std::string infile) {
if (fid.find(infile) == fid.end()) {
return fid.size();
}
return fid[infile];
}
+void EventManager::SetRW(FitWeight* rw){
+ fRW = rw;
+}
+
+
diff --git a/src/FitBase/EventManager.h b/src/FitBase/EventManager.h
index 82ea12b..ff8c3ca 100644
--- a/src/FitBase/EventManager.h
+++ b/src/FitBase/EventManager.h
@@ -1,74 +1,78 @@
// 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 EVENTMANAGER_H
#define EVENTMANAGER_H
#include "InputHandler2.h"
#include "FitWeight.h"
#include "InputUtils.h"
// This class is menat to manage one input file for many distributions
class EventManager {
public:
static EventManager& Get(void);
FitWeight* GetRW();
+ void SetRW(FitWeight* rw);
+
InputHandlerBase* GetInput(int id);
FitEvent* GetEvent(int id, int i);
double GetEventWeight(int id, int i);
InputHandlerBase* AddInput(std::string handle, std::string infile);
void ResetWeightFlags();
int GetInputID(std::string infile);
std::map< int, InputHandlerBase* > GetInputs();
protected:
EventManager();
~EventManager();
static EventManager* m_evtmgrInstance;
FitWeight* fRW;
std::map< std::string, int > fid;
std::map< int, InputHandlerBase* > finputs;
std::map< int, std::vector< bool > > frwneeded;
std::map< int, std::vector< double > > calc_rw;
};
namespace FitBase {
inline
EventManager& EvtManager(){ return EventManager::Get(); };
inline
FitWeight* GetRW(){ return EvtManager().GetRW(); };
inline
+ void SetRW(FitWeight* rw){ EvtManager().SetRW(rw); };
+ inline
int GetInputID(std::string infile){ return EvtManager().GetInputID(infile); };
inline
InputHandlerBase* GetInput(int infile){ return EvtManager().GetInput(infile); };
inline
InputHandlerBase* AddInput(std::string handle, std::string infile){ return EvtManager().AddInput(handle, infile); };
}
#endif
diff --git a/src/FitBase/SplineInputHandler.cxx b/src/FitBase/SplineInputHandler.cxx
index 3581470..a4ac669 100644
--- a/src/FitBase/SplineInputHandler.cxx
+++ b/src/FitBase/SplineInputHandler.cxx
@@ -1,196 +1,195 @@
#include "SplineInputHandler.h"
SplineInputHandler::SplineInputHandler(std::string const& handle, std::string const& rawinputs) {
LOG(SAM) << "Creating SplineInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
jointinput = false;
jointindexswitch = 0;
// Setup the TChain
fFitEventTree = new TChain("nuisance_events");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// 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");
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Input File 0"
<< " : " << rawinputs << std::endl;
}
// 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);
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl
<< "\t\t|-> Event Integral : " << fEventHist->Integral("width") * 1.E-38 << " events/nucleon" << std::endl
<< "\t\t|-> Flux Integral : " << fFluxHist->Integral("width") << " /cm2" << std::endl
<< "\t\t|-> Event/Flux : "
<< fEventHist->Integral("width") * 1.E-38 / fFluxHist->Integral("width") << " cm2/nucleon" << std::endl;
}
- // Setup Reader
- fSplRead = new SplineReader();
- fSplRead->Read( (TTree*)inp_file->Get("spline_reader") );
- fNUISANCEEvent->fSplineRead = this->fSplRead;
-
// Load into memory
-
for (int j = 0; j < fNEvents; j++) {
std::vector<float> tempval;
fFitEventTree->GetEntry(j);
//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;
//}
fStartingWeights.push_back(fNUISANCEEvent->SavedRWWeight * fNUISANCEEvent->InputWeight);
}
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());
// Setup Max Events
if (fMaxEvents > 1 && fMaxEvents < fNEvents) {
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl;
}
fNEvents = fMaxEvents;
}
// Setup Spline Reader
LOG(SAM) << "Loading Spline Reader." << std::endl;
+
+ fSplRead = new SplineReader();
+ fSplRead->Read( (TTree*)inp_file->Get("spline_reader") );
+ fNUISANCEEvent->fSplineRead = this->fSplRead;
+
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);
};
FitEvent* SplineInputHandler::GetNuisanceEvent(const UInt_t entry) {
// Make sure events setup
if (entry >= (UInt_t)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 = fStartingWeights[entry] * GetInputWeight(entry);
} else {
fNUISANCEEvent->InputWeight = fStartingWeights[entry];
}
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 >= (UInt_t)fNEvents) return NULL;
fBaseEvent->eventid = entry;
// Set joint scaling if required
if (jointinput) {
fBaseEvent->InputWeight = fStartingWeights[entry] * GetInputWeight(entry);
} else {
fBaseEvent->InputWeight = fStartingWeights[entry];
}
return fBaseEvent;
}
void SplineInputHandler::Print() {}
diff --git a/src/Routines/SplineRoutines.cxx b/src/Routines/SplineRoutines.cxx
index ce08e95..df3e04e 100755
--- a/src/Routines/SplineRoutines.cxx
+++ b/src/Routines/SplineRoutines.cxx
@@ -1,877 +1,1468 @@
// 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();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
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::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
nuiskey fCompKey;
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
- if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
+ fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) 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++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")) {
std::cout << "[ NUISANCE ] : Overriding " << "MAXEVENTS=" + maxevents << std::endl;
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
FitPar::log_verb = verbocount;
LOG_VERB(verbocount);
ERR_VERB(errorcount);
// 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 (size_t 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();
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 (size_t 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();
else if (!rout.compare("TestSplines_1DEventScan")) TestSplines_1DEventScan();
+ else if (!rout.compare("TestSplines_NDEventThrow")) TestSplines_NDEventThrow();
+ else if (!rout.compare("SaveSplinePlots")) SaveSplinePlots();
+ else if (!rout.compare("TestSplines_1DLikelihoodScan")) TestSplines_1DLikelihoodScan();
+ else if (!rout.compare("TestSplines_NDLikelihoodThrow")) TestSplines_NDLikelihoodThrow();
}
}
//*************************************
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 (size_t 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 (size_t 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);
+ int countwidth = 1000; //(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.
SplineMerger* splmerge = new SplineMerger();
std::vector<nuiskey> splinekeys = Config::QueryKeys("splinemerge");
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
TFile* infile = new TFile(splkey.GetS("input").c_str(), "READ");
splmerge->AddSplineSetFromFile(infile);
}
splmerge->SetupSplineSet();
// Now get Event File
std::vector<nuiskey> eventkeys = Config::QueryKeys("eventmerge");
nuiskey key = eventkeys[0];
std::string inputfilename = key.GetS("input");
// 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]);
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();
// 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
splmerge->Write("spline_reader");
// Setup the spline TTree
TTree* splinetree = new TTree("spline_tree", "spline_tree");
splmerge->AddCoefficientsToTree(splinetree);
int lasttime = time(NULL);
int i = 0;
// Loop over all events and fill the TTree
while (nuisevent) {
// Calculate the weights for each parameter set
splmerge->FillMergedSplines(i);
// Save everything
eventtree->Fill();
splinetree->Fill();
// 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;
}
//*************************************
void SplineRoutines::TestSplines_1DEventScan() {
//*************************************
// Setup RW Engine
if (fRW) delete fRW;
SetupRWEngine();
// Make a spline RW Engine too.
FitWeight* splweight = new FitWeight("splinerwaweight");
// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
// Add Parameters
for (size_t 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");
- parhisttemplate->SetBinContent(i+1, nom);
- parhisttemplate->GetXaxis()->SetBinLabel(i+1, parname.c_str());
+ parhisttemplate->SetBinContent(i + 1, nom);
+ parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
splweight->IncludeDial( key.GetS("name"),
kSPLINEPARAMETER, nom);
splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
}
splweight->Reconfigure();
// Make a high resolution spline set.
std::vector<double> nomvals = fRW->GetDialValues();
int testres = FitPar::Config().GetParI("spline_test_resolution");
- std::vector< std::string > scanparset_names;
std::vector< std::vector<double> > scanparset_vals;
std::vector< TH1D* > scanparset_hists;
// Loop over all params
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
// Get Par Name
std::string name = key.GetS("name");
if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
continue;
}
// Push Back Scan
double low = key.GetD("low");
double high = key.GetD("high");
double cur = low;
double step = key.GetD("step");
while (cur <= high) {
// Make new set
std::vector<double> newvals = nomvals;
newvals[i] = cur;
// Add to vects
- scanparset_names.push_back(name);
scanparset_vals.push_back(newvals);
TH1D* parhist = (TH1D*)parhisttemplate->Clone();
for (size_t j = 0; j < newvals.size(); j++) {
parhist->SetBinContent(j + 1, newvals[j]);
}
scanparset_hists.push_back(parhist);
// Move to next one
cur += step;
}
}
// Print out the parameter set to test
- for (int i = 0; i < scanparset_names.size(); i++) {
+ for (int i = 0; i < scanparset_vals.size(); i++) {
std::cout << "Parset " << i;
for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
std::cout << " " << scanparset_vals[i][j];
}
std::cout << std::endl;
}
// Weight holders
double* rawweights = new double[scanparset_vals.size()];
double* splweights = new double[scanparset_vals.size()];
double* difweights = new double[scanparset_vals.size()];
int NParSets = scanparset_vals.size();
// Loop over all event I/O
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t 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 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]);
// Make handlers for input and output
InputHandlerBase* input = InputUtils::CreateInputHandler("rawevents", inptype, file_descriptor[1]);
InputHandlerBase* output = InputUtils::CreateInputHandler("splineevents", InputUtils::kEVSPLN_Input, outputfilename);
// Get Base Events for each case.
FitEvent* rawevent = input->FirstNuisanceEvent();
FitEvent* splevent = output->FirstNuisanceEvent();
// Setup outputfile
std::string outputtest = outputfilename + ".splinetest.1DEventScan.root";
TFile* outputtestfile = new TFile(outputtest.c_str(), "RECREATE");
outputtestfile->cd();
// Save Parameter Sets
for (size_t i = 0; i < scanparset_hists.size(); i++) {
- scanparset_hists[i]->Write(Form("Paramater_Set_%i",(int)i));
+ scanparset_hists[i]->Write(Form("Paramater_Set_%i", (int)i));
+ }
+
+ // Save a TTree of weights and differences.
+ TTree* weighttree = new TTree("weightscan", "weightscan");
+
+ // Make a branch for each weight set
+ for (size_t i = 0; i < scanparset_hists.size(); i++) {
+ weighttree->Branch(Form("RawWeights_Set_%i", (int)i), &rawweights[i], Form("RawWeights_Set_%i/D", (int)i) );
+ weighttree->Branch(Form("SplineWeights_Set_%i", (int)i), &splweights[i], Form("SplineWeights_Set_%i/D", (int)i) );
+ weighttree->Branch(Form("DifWeights_Set_%i", (int)i), &difweights[i], Form("DifWeights_Set_%i/D", (int)i) );
+
+ }
+
+ // Count
+ int i = 0;
+ int nevents = input->GetNEvents();
+ while (rawevent and splevent) {
+
+ // Loop over 1D parameter sets.
+ for (size_t j = 0; j < scanparset_vals.size(); j++) {
+
+ // Reconfigure
+ fRW->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
+ fRW->Reconfigure();
+
+ // Reconfigure spline RW
+ splweight->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
+ splweight->Reconfigure();
+
+ splevent->fSplineRead->SetNeedsReconfigure(true);
+
+ // Calc weight for both events
+ rawweights[j] = fRW->CalcWeight(rawevent);
+ splweights[j] = splweight->CalcWeight(splevent);
+ difweights[j] = splweights[j] - rawweights[j];
+ }
+
+
+ if (i % 1000 == 0) {
+ LOG(FIT) << "Processed " << i << "/" << nevents << std::endl;
+ }
+
+ // Fill Array
+ weighttree->Fill();
+
+ // Iterate to next event.
+ i++;
+ rawevent = input->NextNuisanceEvent();
+ splevent = output->NextNuisanceEvent();
+ }
+
+ outputtestfile->cd();
+ weighttree->Write();
+ outputtestfile->Close();
+ }
+}
+
+
+
+//*************************************
+void SplineRoutines::TestSplines_NDEventThrow() {
+//*************************************
+
+ // Setup RW Engine
+ if (fRW) delete fRW;
+ SetupRWEngine();
+
+ // Make a spline RW Engine too.
+ FitWeight* splweight = new FitWeight("splinerwaweight");
+ std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
+ TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
+
+ // Add Parameters
+ for (size_t 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");
+
+ parhisttemplate->SetBinContent(i + 1, nom);
+ parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
+
+ splweight->IncludeDial( key.GetS("name"),
+ kSPLINEPARAMETER, nom);
+ splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
+
+ }
+ splweight->Reconfigure();
+
+ // Make a high resolution spline set.
+ std::vector<double> nomvals = fRW->GetDialValues();
+ int testres = FitPar::Config().GetParI("spline_test_resolution");
+
+ std::vector< std::string > scanparset_names;
+ std::vector< std::vector<double> > scanparset_vals;
+ std::vector< TH1D* > scanparset_hists;
+
+ // Loop over all params
+ // Add Parameters
+ int nthrows = FitPar::Config().GetParI("spline_test_throws");
+ for (int i = 0; i < nthrows; i++) {
+
+ std::vector<double> newvals = nomvals;
+
+ for (size_t j = 0; j < parameterkeys.size(); j++) {
+ nuiskey key = parameterkeys[j];
+
+ if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
+ continue;
+ }
+
+ // Push Back Scan
+ double low = key.GetD("low");
+ double high = key.GetD("high");
+ newvals[j] = gRandom->Uniform(low, high);
+
+ }
+ // Add to vects
+ scanparset_vals.push_back(newvals);
+
+ TH1D* parhist = (TH1D*)parhisttemplate->Clone();
+ for (size_t j = 0; j < newvals.size(); j++) {
+ parhist->SetBinContent(j + 1, newvals[j]);
+ }
+ scanparset_hists.push_back(parhist);
+ }
+
+ // Print out the parameter set to test
+ for (int i = 0; i < scanparset_vals.size(); i++) {
+ std::cout << "Parset " << i;
+ for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
+ std::cout << " " << scanparset_vals[i][j];
+ }
+ std::cout << std::endl;
+ }
+
+
+ // Weight holders
+ double* rawweights = new double[scanparset_vals.size()];
+ double* splweights = new double[scanparset_vals.size()];
+ double* difweights = new double[scanparset_vals.size()];
+
+ int NParSets = scanparset_vals.size();
+
+ // Loop over all event I/O
+ std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
+ for (size_t 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 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]);
+
+
+ // Make handlers for input and output
+ InputHandlerBase* input = InputUtils::CreateInputHandler("rawevents", inptype, file_descriptor[1]);
+ InputHandlerBase* output = InputUtils::CreateInputHandler("splineevents", InputUtils::kEVSPLN_Input, outputfilename);
+
+ // Get Base Events for each case.
+ FitEvent* rawevent = input->FirstNuisanceEvent();
+ FitEvent* splevent = output->FirstNuisanceEvent();
+
+
+ // Setup outputfile
+ std::string outputtest = outputfilename + ".splinetest.NDEventThrow.root";
+ TFile* outputtestfile = new TFile(outputtest.c_str(), "RECREATE");
+ outputtestfile->cd();
+
+ // Save Parameter Sets
+ for (size_t i = 0; i < scanparset_hists.size(); i++) {
+ scanparset_hists[i]->Write(Form("Paramater_Set_%i", (int)i));
}
// Save a TTree of weights and differences.
TTree* weighttree = new TTree("weightscan", "weightscan");
// Make a branch for each weight set
for (size_t i = 0; i < scanparset_hists.size(); i++) {
weighttree->Branch(Form("RawWeights_Set_%i", (int)i), &rawweights[i], Form("RawWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("SplineWeights_Set_%i", (int)i), &splweights[i], Form("SplineWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("DifWeights_Set_%i", (int)i), &difweights[i], Form("DifWeights_Set_%i/D", (int)i) );
}
// Count
int i = 0;
int nevents = input->GetNEvents();
while (rawevent and splevent) {
// Loop over 1D parameter sets.
for (size_t j = 0; j < scanparset_vals.size(); j++) {
// Reconfigure
fRW->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
fRW->Reconfigure();
// Reconfigure spline RW
splweight->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
splweight->Reconfigure();
splevent->fSplineRead->SetNeedsReconfigure(true);
// Calc weight for both events
rawweights[j] = fRW->CalcWeight(rawevent);
splweights[j] = splweight->CalcWeight(splevent);
difweights[j] = splweights[j] - rawweights[j];
}
if (i % 1000 == 0) {
LOG(FIT) << "Processed " << i << "/" << nevents << std::endl;
}
// Fill Array
weighttree->Fill();
// Iterate to next event.
i++;
rawevent = input->NextNuisanceEvent();
splevent = output->NextNuisanceEvent();
}
outputtestfile->cd();
weighttree->Write();
outputtestfile->Close();
}
}
+void SplineRoutines::SaveSplinePlots() {
+
+ 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 (size_t 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
+ outputfilename += ".SplinePlots.root";
+ 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();
+
+ outputfile->cd();
+
+ int lasttime = time(NULL);
+ TCanvas* fitcanvas = 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, fitcanvas, true);
+
+ if (fitcanvas) {
+ outputfile->cd();
+ fitcanvas->Write(Form("Event_SplineCanvas_%i", (int)i));
+ }
+
+ // Logging
+ if (i % countwidth == 0) {
+ LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline plots. " << std::endl;
+ }
+
+ // Iterate
+ i++;
+ nuisevent = input->NextNuisanceEvent();
+ }
+ // Save flux and close file
+ outputfile->cd();
+
+ // Close Output
+ outputfile->Close();
+
+ // Delete Inputs
+ delete input;
+ }
+
+ // remove Keys
+ eventkeys.clear();
+
+}
+
+void SplineRoutines::TestSplines_NDLikelihoodThrow() {
+
+ // Setup RW Engine
+ if (fRW) delete fRW;
+ SetupRWEngine();
+
+ // Make a spline RW Engine too.
+ FitWeight* splweight = new FitWeight("splinerwaweight");
+ std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
+ TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
+
+ // Add Parameters
+ for (size_t 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");
+
+ parhisttemplate->SetBinContent(i + 1, nom);
+ parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
+
+ splweight->IncludeDial( key.GetS("name"),
+ kSPLINEPARAMETER, nom);
+ splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
+
+ }
+ splweight->Reconfigure();
+
+ // Make a high resolution spline set.
+ std::vector<double> nomvals = fRW->GetDialValues();
+ int testres = FitPar::Config().GetParI("spline_test_resolution");
+
+ std::vector< std::string > scanparset_names;
+ std::vector< std::vector<double> > scanparset_vals;
+ std::vector< TH1D* > scanparset_hists;
+
+ // Loop over all params
+ // Add Parameters
+ int nthrows = FitPar::Config().GetParI("spline_test_throws");
+ for (int i = 0; i < nthrows; i++) {
+
+ std::vector<double> newvals = nomvals;
+
+ for (size_t j = 0; j < parameterkeys.size(); j++) {
+ nuiskey key = parameterkeys[j];
+
+ if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
+ continue;
+ }
+
+ // Push Back Scan
+ double low = key.GetD("low");
+ double high = key.GetD("high");
+ newvals[j] = gRandom->Uniform(low, high);
+
+ }
+ // Add to vects
+ scanparset_vals.push_back(newvals);
+
+ TH1D* parhist = (TH1D*)parhisttemplate->Clone();
+ for (size_t j = 0; j < newvals.size(); j++) {
+ parhist->SetBinContent(j + 1, newvals[j]);
+ }
+ scanparset_hists.push_back(parhist);
+ }
+
+ // Print out the parameter set to test
+ for (int i = 0; i < scanparset_vals.size(); i++) {
+ std::cout << "Parset " << i;
+ for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
+ std::cout << " " << scanparset_vals[i][j];
+ }
+ std::cout << std::endl;
+ }
+
+ // Make a new set of Raw/Spline Sample Keys
+ std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
+ std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");
+
+ std::vector<nuiskey> rawkeys;
+ std::vector<nuiskey> splkeys;
+
+
+ for (std::vector<nuiskey>::iterator iter = testkeys.begin();
+ iter != testkeys.end(); iter++) {
+ nuiskey key = (*iter);
+
+ 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");
+ std::string splfile = eventskey.GetS("output");
+
+ nuiskey rawkeytemp = Config::CreateKey("sample");
+ rawkeytemp.SetS("name", samplename);
+ rawkeytemp.SetS("input", rawfile);
+
+ nuiskey splkeytemp = Config::CreateKey("sample");
+ splkeytemp.SetS("name", samplename);
+ splkeytemp.SetS("input", "EVSPLN:" + splfile);
+
+ rawkeys.push_back(rawkeytemp);
+ splkeys.push_back(splkeytemp);
+ }
+
+ if (fOutputRootFile) delete fOutputRootFile;
+ fOutputRootFile = new TFile(fOutputFile.c_str(), "RECREATE");
+
+ fOutputRootFile->ls();
+ // Make two new JointFCN
+ JointFCN* rawfcn = new JointFCN(rawkeys, fOutputRootFile);
+ JointFCN* splfcn = new JointFCN(splkeys, fOutputRootFile);
+
+ // Create iteration tree in output file
+ fOutputRootFile->cd();
+ rawfcn->CreateIterationTree("raw_iterations", fRW);
+ splfcn->CreateIterationTree("spl_iterations", splweight);
+
+ // Loop over parameter sets.
+ for (size_t j = 0; j < scanparset_vals.size(); j++) {
+
+ FitBase::SetRW(fRW);
+ double rawtotal = rawfcn->DoEval(&scanparset_vals[j][0]);
+
+ FitBase::SetRW(splweight);
+ double spltotal = splfcn->DoEval(&scanparset_vals[j][0]);
+
+ LOG(FIT) << "RAW SPLINE DIF = " << rawtotal << " " << spltotal << " " << spltotal - rawtotal << std::endl;
+ }
+
+ fOutputRootFile->cd();
+
+ rawfcn->WriteIterationTree();
+ splfcn->WriteIterationTree();
+
+}
+
+
+void SplineRoutines::TestSplines_1DLikelihoodScan() {
+
+ // Setup RW Engine.
+ if (fRW) delete fRW;
+ SetupRWEngine();
+
+ // Setup Parameter Set.
+ // Make a spline RW Engine too.
+ FitWeight* splweight = new FitWeight("splinerwaweight");
+ // std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
+ std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
+ TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
+
+ // Add Parameters
+ for (size_t 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");
+
+ parhisttemplate->SetBinContent(i + 1, nom);
+ parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
+
+ splweight->IncludeDial( key.GetS("name"),
+ kSPLINEPARAMETER, nom);
+ splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
+
+ }
+ splweight->Reconfigure();
+
+ // Make a high resolution spline set.
+ std::vector<double> nomvals = fRW->GetDialValues();
+ int testres = FitPar::Config().GetParI("spline_test_resolution");
+
+ std::vector< std::vector<double> > scanparset_vals;
+ std::vector< TH1D* > scanparset_hists;
+
+ // Loop over all params
+ // Add Parameters
+ for (size_t i = 0; i < parameterkeys.size(); i++) {
+ nuiskey key = parameterkeys[i];
+
+ // Get Par Name
+ std::string name = key.GetS("name");
+
+ if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
+ continue;
+ }
+
+ // Push Back Scan
+ double low = key.GetD("low");
+ double high = key.GetD("high");
+ double cur = low;
+ double step = key.GetD("step");
+
+ while (cur <= high) {
+
+ // Make new set
+ std::vector<double> newvals = nomvals;
+ newvals[i] = cur;
+
+ // Add to vects
+ scanparset_vals.push_back(newvals);
+
+ TH1D* parhist = (TH1D*)parhisttemplate->Clone();
+ for (size_t j = 0; j < newvals.size(); j++) {
+ parhist->SetBinContent(j + 1, newvals[j]);
+ }
+ scanparset_hists.push_back(parhist);
+
+
+ // Move to next one
+ cur += step;
+ }
+ }
+
+ // Print out the parameter set to test
+ for (int i = 0; i < scanparset_vals.size(); i++) {
+ std::cout << "Parset " << i;
+ for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
+ std::cout << " " << scanparset_vals[i][j];
+ }
+ std::cout << std::endl;
+ }
+
+ // Make a new set of Raw/Spline Sample Keys
+ std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
+ std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");
+
+ std::vector<nuiskey> rawkeys;
+ std::vector<nuiskey> splkeys;
+
+
+ for (std::vector<nuiskey>::iterator iter = testkeys.begin();
+ iter != testkeys.end(); iter++) {
+ nuiskey key = (*iter);
+
+ 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");
+ std::string splfile = eventskey.GetS("output");
+
+ nuiskey rawkeytemp = Config::CreateKey("sample");
+ rawkeytemp.SetS("name", samplename);
+ rawkeytemp.SetS("input", rawfile);
+
+ nuiskey splkeytemp = Config::CreateKey("sample");
+ splkeytemp.SetS("name", samplename);
+ splkeytemp.SetS("input", "EVSPLN:" + splfile);
+
+ rawkeys.push_back(rawkeytemp);
+ splkeys.push_back(splkeytemp);
+ }
+
+ if (fOutputRootFile) delete fOutputRootFile;
+ fOutputRootFile = new TFile(fOutputFile.c_str(), "RECREATE");
+
+ fOutputRootFile->ls();
+ // Make two new JointFCN
+ JointFCN* rawfcn = new JointFCN(rawkeys, fOutputRootFile);
+ JointFCN* splfcn = new JointFCN(splkeys, fOutputRootFile);
+
+ // Create iteration tree in output file
+ fOutputRootFile->cd();
+ rawfcn->CreateIterationTree("raw_iterations", fRW);
+ splfcn->CreateIterationTree("spl_iterations", splweight);
+
+ // Loop over parameter sets.
+ for (size_t j = 0; j < scanparset_vals.size(); j++) {
+
+ FitBase::SetRW(fRW);
+ double rawtotal = rawfcn->DoEval(&scanparset_vals[j][0]);
+
+ FitBase::SetRW(splweight);
+ double spltotal = splfcn->DoEval(&scanparset_vals[j][0]);
+
+ LOG(FIT) << "RAW SPLINE DIF = " << rawtotal << " " << spltotal << " " << spltotal - rawtotal << std::endl;
+ }
+
+ fOutputRootFile->cd();
+
+ rawfcn->WriteIterationTree();
+ splfcn->WriteIterationTree();
+
+}
+
+
/*
MISC Functions
*/
//*************************************
int SplineRoutines::GetStatus() {
//*************************************
return 0;
}
diff --git a/src/Routines/SplineRoutines.h b/src/Routines/SplineRoutines.h
index 58b05a0..54170a9 100755
--- a/src/Routines/SplineRoutines.h
+++ b/src/Routines/SplineRoutines.h
@@ -1,197 +1,194 @@
// 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 SPLINE_ROUTINES_H
#define SPLINE_ROUTINES_H
/*!
* \addtogroup Minimizer
* @{
*/
#include "TH1.h"
#include "TF1.h"
#include "TMatrixD.h"
#include "TVectorD.h"
#include "TSystem.h"
#include "TFile.h"
#include "TProfile.h"
#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <cstring>
#include "FitEvent.h"
#include "JointFCN.h"
#include "FitParameters.h"
#include "FitLogger.h"
#include "BaseFitEvt.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#include "SplineReader.h"
#include "SplineWriter.h"
#include "SplineMerger.h"
#include "ParserUtils.h"
enum minstate {
kErrorStatus = -1,
kGoodStatus,
kFitError,
kNoChange,
kFitFinished,
kFitUnfinished,
kStateChange,
};
//*************************************
//! Collects all possible fit routines into a single class to avoid repeated code
class SplineRoutines{
//*************************************
public:
/*
Constructor/Destructor
*/
//! Constructor reads in arguments given at the command line for the fit here.
SplineRoutines(int argc, char* argv[]);
//! Default destructor
~SplineRoutines();
//! Reset everything to default/NULL
void Init();
/*
Input Functions
*/
//! Splits the arguments ready for initial setup
void ParseArgs(int argc, char* argv[]);
/*
Setup Functions
*/
//! Setup the configuration given the arguments passed at the commandline and card file
void SetupConfig();
//! Setups up our custom RW engine with all the parameters passed in the card file
void SetupRWEngine();
void Run();
void SaveEvents();
void TestEvents();
void GenerateEventSplines();
/*
Testing Functions
*/
/// Scan parameter space in 1D at finer resolution than points. Compare Raw/Spline on an event by event basis.
void TestSplines_1DEventScan();
/// Scan parameter space in 1D at finer resolution than points. Compare likelihoods for all testsamples.
- void TestSplines_1DLikelihoodScan(){};
-
- /// Scan parameter space in 2D at 0.5 point resolution, compare average weight difference.
- void TestSplines_2DEventScan(){};
-
- /// Scan parameter space in 2D at 0.5 point resolution, compare likelihoods for all testsamples.
- void TestSplines_2DLikelihoodScan(){};
+ void TestSplines_1DLikelihoodScan();
/// Randomly throw in parameter space. For each throw, calc average weight difference.
- void TestSplines_NDEventThrow(){};
+ void TestSplines_NDEventThrow();
/// Randomly thow in parameter space. For each throw, calc likelihood difference for each sample.
- void TestSplines_NDLikelihoodThrow(){};
+ void TestSplines_NDLikelihoodThrow();
+
+ /// Generate a set of spline vs weight canvases and save to file
+ void SaveSplinePlots();
/*
Fitting Functions
*/
//! Given a new map change the values that the RW engine is currently set to
void UpdateRWEngine(std::map<std::string,double>& updateVals);
/*
MISC Functions
*/
//! Get previous fit status from a file
Int_t GetStatus();
void MergeSplines();
protected:
//! Our Custom ReWeight Object
FitWeight* rw;
FitWeight* fRW;
std::string fOutputFile;
std::string fInputFile;
TFile* fInputRootFile;
TFile* fOutputRootFile;
JointFCN* fSampleFCN;
std::list<MeasurementBase*> fSamples;
std::string fCardFile;
std::string fStrategy;
std::vector<std::string> fRoutines;
std::string fAllowedRoutines;
// Input Dial Vals
//! Vector of dial names
std::vector<std::string> fParams;
std::map<std::string, std::string> fStateVals;
std::map<std::string, double> fStartVals;
std::map<std::string, double> fCurVals;
std::map<std::string, double> fErrorVals;
std::map<std::string, double> fMinVals;
std::map<std::string, double> fMaxVals;
std::map<std::string, double> fStepVals;
std::map<std::string, int> fTypeVals;
std::map<std::string, bool> fFixVals;
std::map<std::string, bool> fStartFixVals;
std::vector<std::string> fGenericInputNames;
std::map<std::string, std::string> fGenericInputFiles;
std::map<std::string, std::string> fGenericOutputFiles;
std::map<std::string, std::string> fGenericOutputTypes;
std::map<std::string, InputHandlerBase*> fGenericInputs;
std::vector<std::string> fSplineNames;
std::map<std::string, std::string> fSplineTypes;
std::map<std::string, std::string> fSplinePoints;
nuiskey fCompKey;
};
/*! @} */
#endif
diff --git a/src/Splines/CMakeLists.txt b/src/Splines/CMakeLists.txt
index d62ef5a..aada294 100644
--- a/src/Splines/CMakeLists.txt
+++ b/src/Splines/CMakeLists.txt
@@ -1,62 +1,64 @@
# 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
+SplineUtils.cxx
Spline.cxx
)
set(HEADERFILES
FitSpline.h
FitSplineHead.h
SplineReader.h
SplineWriter.h
+SplineUtils.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 78692df..4367761 100644
--- a/src/Splines/Spline.cxx
+++ b/src/Splines/Spline.cxx
@@ -1,368 +1,569 @@
#include "Spline.h"
+using namespace SplineUtils;
+// Setup Functions
+// ----------------------------------------------
+Spline::Spline(std::string splname, std::string form, std::string points) {
-Spline::Spline(std::string splname, std::string form,
- std::vector<double> x) {
-
- // Initial setup
- fSplineNames = GeneralUtils::ParseToStr(splname, ":");
+ // Save Definition
fName = splname;
fForm = form;
+ fPoints = points;
+// functor = NULL;
+ minimizer = NULL;
- fXScan.clear();
- for (size_t i = 0; i < x.size(); i++){
- fXScan.push_back(x[i]);
- }
+ // Run Checks
+ // if ((UInt_t)fNDim != fSplineNames.size()) {
+ // ERR(FTL) << "Spline Dim:Names mismatch!" << std::endl;
+ // throw;
+ // }
- fX = 0.0;
+ // Setup Min Max for each Parameter
+ fSplitNames = GeneralUtils::ParseToStr(splname, ";");
+ std::vector< std::vector<double> > gridvals = SplineUtils::GetSplitDialPoints(fPoints);
- // 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];
- }
+ for (size_t i = 0; i < fSplitNames.size(); i++) {
- // Setup iters
- iter_low = fXScan.begin();
- iter_high = fXScan.begin();
- iter_high++;
- off = 0;
+ for (size_t j = 0; j < gridvals.size(); j++) {
+ if (i == 0) fXScan.push_back(gridvals[j][0]);
+ if (i == 1) fXScan.push_back(gridvals[j][1]);
+ }
+
+ double xmin = 9999.9;
+ double xmax = -9999.9;
+ for (size_t j = 0; j < gridvals.size(); j++) {
+ if (gridvals[j][i] < xmin) xmin = gridvals[j][i];
+ if (gridvals[j][i] > xmax) xmax = gridvals[j][i];
+ }
+
+ fVal.push_back(0.0);
+ fValMin.push_back(xmin);
+ fValMax.push_back(xmax);
+
+ // Define TSpline3 1D iterators here
+ if (i == 0) {
+ iter_low = fXScan.begin();
+ iter_high = fXScan.begin();
+ iter_high++;
+ off = 0;
+ }
+ }
// 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 if (!fForm.compare("2DPol6")) { Setup( k2DPol6, 2, 28 ); }
+ else if (!fForm.compare("2DGaus")) { Setup( k2DGaus, 2, 8 ); }
+ else if (!fForm.compare("2DTSpline3")) { Setup (k2DTSpline3, 2, fXScan.size() * fYScan.size() * 8); }
else {
ERR(FTL) << "Unknown spline form : " << fForm << std::endl;
throw;
}
- // Run Checks
- if ((UInt_t)fNDim != fSplineNames.size()) {
- ERR(FTL) << "Spline Dim:Names mismatch!" << std::endl;
- throw;
- }
-
+ std::cout << "Setup Spline " << fForm << " = " << fType << " " << fNPar << std::endl;
};
void Spline::Setup(int type, int ndim, int npar) {
fType = type;
fNDim = ndim;
fNPar = npar;
}
+
+// Reconfigure Functions
+// ----------------------------------------------
+void Spline::Reconfigure(float x, int index) {
+ // std::cout << "Reconfigured spline : " << fName << " : " << fForm << " to be " << x << " " << index << std::endl;
+ fVal[index] = x;
+ fOutsideLimits = false;
+
+ if (fVal[index] > fValMax[index]) fVal[index] = fValMax[index];
+ if (fVal[index] < fValMin[index]) fVal[index] = fValMin[index];
+ // std::cout << "Set at edge = " << fVal[index] << " " << index << std::endl;
+}
+
+void Spline::Reconfigure(std::string name, float x) {
+ for (size_t i = 0; i < fSplitNames.size(); i++) {
+ // std::cout << "-> Comparing in spline " << name << " to " << fSplitNames[i] << " = " << !fSplitNames[i].compare(name.c_str()) << std::endl;
+
+ if (!fSplitNames[i].compare(name.c_str())) {
+ // std::cout << "-> Reconfigured spline : " << fSplitNames[i] << " " << name << " to be " << x << " " << i << std::endl;
+ Reconfigure(x, i);
+ }
+ }
+
+}
+
+
+// Evaluation Functions
+// ----------------------------------------------
double Spline::operator()(const Double_t* x, const Double_t* par) const {
- Float_t tempx;
- tempx = x[0];
+ Float_t* tempx = new Float_t[fNDim];
+ for (size_t i = 0; i < (UInt_t)fNDim; i++) {
+ tempx[i] = x[i];
+ }
Float_t* tempp = new Float_t[fNPar];
- for (size_t i = 0; i < (UInt_t)fNPar; i++){
+ for (size_t i = 0; i < (UInt_t)fNPar; i++) {
tempp[i] = par[i];
}
- float val = DoEval(&tempx, tempp);
+ float val = DoEval(tempx, tempp);
delete tempp;
+ delete tempx;
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;
- fX = x;
- fOutsideLimits = false;
+// double Spline::operator()(const Double_t* x, const Double_t* y, const Double_t* par) const {
+// Float_t temp[2];
+// temp[0] = x[0];
+// temp[1] = y[0];
- if (fX > fXMax) fX = fXMax;
- if (fX < fXMin) fX = fXMin;
-}
-*/
+// float val = DoEval(temp, par);
+// if (val < 0.0) val = 0.0;
+// return val;
+// }
-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;
- 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;
- 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;
- if (fX < fXMin) fX = fXMin;
+ for (size_t i = 0; i < (UInt_t) fNDim; i++) {
+ fVal[i] = x[i];
+ if (fVal[i] > fValMax[i]) fVal[i] = fValMax[i];
+ if (fVal[i] < fValMin[i]) fVal[i] = fValMin[i];
+ // std::cout << "Set " << i << " fVal to " << x[i] << std::endl;
+ }
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 {
- if (!par) return 1.0;
+ if (!par)
+ return 1.0;
- // std::cout << "Spline::DoEval = " << par << std::endl;
+ // std::cout << "Spline::DoEval = " << par << std::endl;
// 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;
}
}
+ // std::cout << "TYpe = " << fType << " "<< fForm << std::endl;
// 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); }
+ case k2DPol6: { return Spline2DPol(par, 6); }
+ case k2DGaus: { return Spline2DGaus(par); }
+ case k2DTSpline3: { return Spline2DTSpline3(par); }
}
// Return nominal weight
return 1.0;
};
// Spline Functions
// ----------------------------------------------
+
+// 1D 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;
+ float xp = fVal[0];
return par[0] + par[1] * xp;
};
float Spline::Spline1DPol2(const Float_t* par) const {
- float xp = fX;
+ float xp = fVal[0];
return par[0] + par[1] * xp + par[2] * xp * xp;
};
float Spline::Spline1DPol3(const Float_t* par) const {
- float xp = fX;
+ float xp = fVal[0];
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;
+ float xp = fVal[0];
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;
+ float xp = fVal[0];
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 xp = fVal[0];
float w = 0.0;
// std::cout << "Pol Eval " << std::endl;
- for (int i = fNPar-1; i > 0; i--){
- w = xp * (par[0+i] + w);
+ for (int i = fNPar - 1; i > 0; i--) {
+ w = xp * (par[0 + i] + w);
}
- w += par[0];
+ w += par[0];
return w;
};
-float Spline::Spline1DTSpline3(const Float_t* par) const {
- // std::cout << "Doing Spline Eval " << fX << " " << par << std::endl;
+
+float Spline::GetMonomial(int p) const {
+
+ // // Copied From ROOT
+ // Int_t i = 0;
+ // Double_t p1 = 1;
+ // Double_t p2 = 0;
+ // Double_t p3 = 0;
+ // Double_t r = 0;
+
+ // switch (p) {
+ // case 1:
+ // r = 1;
+ // break;
+ // case 2:
+ // r = x;
+ // break;
+ // default:
+ // p2 = x;
+ // for (i = 3; i <= p; i++) {
+ // p3 = p2 * x;
+ // p1 = p2;
+ // p2 = p3;
+ // }
+ // r = p3;
+ // }
+
+ return 0.0;
+}
+
+float Spline::Spline2DTSpline3(const Float_t* par) const {
+
// Find matching point
- iter_low = fXScan.begin();
- iter_high = fXScan.begin();
- iter_high++;
+ std::vector<float>::iterator iter_low_x = fXScan.begin();
+ std::vector<float>::iterator iter_high_x = fXScan.begin();
+ std::vector<float>::iterator iter_low_y = fYScan.begin();
+ std::vector<float>::iterator iter_high_y = fYScan.begin();
+ iter_high_x++;
+ iter_high_y++;
+
off = 0;
-
- while ( iter_high != fXScan.end() and
- (fX < (*iter_low) or fX >= (*iter_high)) ) {
- off += 4;
- iter_low++;
- iter_high++;
+ fX = fVal[0];
+ fY = fVal[1];
+
+ while ( (iter_high_x != fXScan.end() and
+ iter_high_y != fYScan.end()) and
+ (fX < (*iter_low_y) or
+ fX >= (*iter_high) or
+ fY < (*iter_low_y) or
+ fY >= (*iter_low_y)) ) {
+ off += 9;
+ iter_low_x++;
+ iter_high_x++;
+
+ if (iter_high_x == fXScan.end()) {
+ iter_low_x = fXScan.begin();
+ iter_high_x = fXScan.begin();
+ iter_low_y++;
+ iter_high_y++;
+ std::cout << "Skipping to next tspline 3 rung " << *iter_low_y << std::endl;
+ }
}
- float dx = fX - (*iter_low);
+ std::cout << "Evaluting TSpline3 at " << fX << " " << (*iter_low_x) << " " << fY << " " << (*iter_low_y) << std::endl;
+ // sleep(1.0);
+ float dx = fX - (*iter_low_x);
+ float dy = fY - (*iter_low_y);
+
float weight = (par[off] + dx * (par[off + 1] + dx * (par[off + 2] + dx * par[off + 3])));
+ float weight2 = (par[off + 4] + dy * (par[off + 5] + dy * (par[off + 6] + dy * par[off + 7])));
- return weight;
+ return weight * weight2 * par[off + 8];
};
-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;
- }
-}
+float Spline::Spline2DPol(const Float_t* par, int n) const {
+ // float wx = 2.0 * (fVal[0] - fValMin[0] - (fValMax[0] - fValMin[0]) / 2.0) / (fValMax[0] - fValMin[0]);
+ // float wy = 2.0 * (fVal[1] - fValMin[1] - (fValMax[1] - fValMin[1]) / 2.0) / (fValMax[1] - fValMin[1]);
+ float wx = (fVal[0] - fValMin[0]) / (fValMax[0] - fValMin[0]);
+ float wy = (fVal[1] - fValMin[1]) / (fValMax[1] - fValMin[1]);
-// Fitting Functions
-void Spline::FitCoeff1DGraph(int n, double* x, double* y, float* coeff, bool draw) {
+ int count = 0;
+ bool speak = false;
- 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];
- }
+ // float w = 0.0;
+ // // std::cout << "Pol Eval " << std::endl;
+ // for (int i = 7 - 1; i > 0; i--) {
+ // w = xp * (par[0 + i] + w);
+ // }
+ // w += par[0];
- double xwidth = xmax - xmin;
- xmin = xmin - xwidth * 0.01;
- xmax = xmax + xwidth * 0.01;
+ // float w2 = 0.0;
+ // // std::cout << "Pol Eval " << std::endl;
+ // for (int i = 13 - 1; i > 0; i--) {
+ // w2 = xp * (par[7 + i] + w2);
+ // }
- // 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
+ // return w1 * w2;
- // Run the actual spline fit
- StopTalking();
-
- // If linear fit with two points
- if (n == 2 and fType == k1DPol1){
- float m = (y[1] - y[0]) / (x[1] - x[0]);
- float c = y[0] - (0.0 - x[0]) * m;
- func->SetParameter(0, c);
- func->SetParameter(1, m);
- } else if (fType == k1DPol1){
- gr->Fit(func, "WQ");
- } else {
- gr->Fit(func, "FMWQ");
- }
- StartTalking();
+ // for (int i = 0; i < n; i++){
+
+ // if (i == 0)
+
+ // }
+// float w = 0.0;
+// w += par[count++];
+// w += par[count++] * wx;
+// w += par[count++] * wy;
+
+// w += par[count++] * wx * wy; // / 10.0;
+// w += par[count++] * wx * wx; // / 10.0;
+// w += par[count++] * wy * wy; // / 10.0;
+
+// w += par[count++] * wx * wx * wy; // / 100.0;
+// w += par[count++] * wx * wy * wy; // / 100.0;
+// w += par[count++] * wx * wx * wx; // / 100.0;
+// w += par[count++] * wy * wy * wy; // / 100.0;
+
+
+// w += par[count++] * wx * wx * wx * wx; // / 100.0;
+// w += par[count++] * wx * wx * wx * wy; // / 100.0;
+// w += par[count++] * wx * wx * wy * wy; // / 100.0;
+// w += par[count++] * wx * wy * wy * wy; // / 100.0;
+// w += par[count++] * wy * wy * wy * wy; // / 100.0;
+
+// w += par[count++] * wx * wx * wx * wx * wx; // / 1000.0;
+// w += par[count++] * wx * wx * wx * wx * wy; // / 1000.0;
+// w += par[count++] * wx * wx * wx * wy * wy; // / 1000.0;
+// w += par[count++] * wx * wx * wy * wy * wy; // / 1000.0;
+// w += par[count++] * wx * wy * wy * wy * wy; // / 1000.0;
+// w += par[count++] * wy * wy * wy * wy * wy; // / 1000.0;
+
+// w += par[count++] * wx * wx * wx * wx * wx * wx; // / 10000.0;
+// w += par[count++] * wx * wx * wx * wx * wx * wy; // / 10000.0;
+// w += par[count++] * wx * wx * wx * wx * wy * wy; // / 10000.0;
+// w += par[count++] * wx * wx * wx * wy * wy * wy; // / 10000.0;
+// w += par[count++] * wx * wx * wy * wy * wy * wy; // / 10000.0;
+// w += par[count++] * wx * wy * wy * wy * wy * wy; // / 10000.0;
+// w += par[count++] * wy * wy * wy * wy * wy * wy; // / 10000.0;
+
+ // float w = 0.0;
+ // int power = 6;
+ // for (int i = 0; i < power; i++){
+ // int xpower = power;
+ // int ypower = 0;
+ // while (ypower <= power){
+ // w += par[count++] * pow(wx,xpower) * pow(wy, ypower);
+ // xpower--;
+ // ypower++;
+ // }
+ // }
- 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);
+ float w = 0.0;
+
+ w += par[count++] * wx;
+ w += par[count++] * wy;
+
+ w += par[count++] * wx * wy; // / 10.0;
+ w += par[count++] * wx * wx; // / 10.0;
+ w += par[count++] * wy * wy; // / 10.0;
+
+ w += par[count++] * wx * wx * wy; // / 100.0;
+ w += par[count++] * wx * wy * wy; // / 100.0;
+ w += par[count++] * wx * wx * wx; // / 100.0;
+ w += par[count++] * wy * wy * wy; // / 100.0;
+
+
+ w += par[count++] * wx * wx * wx * wx; // / 100.0;
+ w += par[count++] * wx * wx * wx * wy; // / 100.0;
+ w += par[count++] * wx * wx * wy * wy; // / 100.0;
+ w += par[count++] * wx * wy * wy * wy; // / 100.0;
+ w += par[count++] * wy * wy * wy * wy; // / 100.0;
+
+ w += par[count++] * wx * wx * wx * wx * wx; // / 1000.0;
+ w += par[count++] * wx * wx * wx * wx * wy; // / 1000.0;
+ w += par[count++] * wx * wx * wx * wy * wy; // / 1000.0;
+ w += par[count++] * wx * wx * wy * wy * wy; // / 1000.0;
+ w += par[count++] * wx * wy * wy * wy * wy; // / 1000.0;
+ w += par[count++] * wy * wy * wy * wy * wy; // / 1000.0;
+
+ w += par[count++] * wx * wx * wx * wx * wx * wx; // / 10000.0;
+ w += par[count++] * wx * wx * wx * wx * wx * wy; // / 10000.0;
+ w += par[count++] * wx * wx * wx * wx * wy * wy; // / 10000.0;
+ w += par[count++] * wx * wx * wx * wy * wy * wy; // / 10000.0;
+ w += par[count++] * wx * wx * wy * wy * wy * wy; // / 10000.0;
+ w += par[count++] * wx * wy * wy * wy * wy * wy; // / 10000.0;
+ w += par[count++] * wy * wy * wy * wy * wy * wy; // / 10000.0;
+
+
+
+ if (wx != 0.0 and speak) {
+ for (int i = 0; i < count; i++) {
+ std::cout << "Evaluated " << fName << " spline coeff " << i << " = " << par[i] << std::endl;
+ }
+ std::cout << "End Weight = " << w << " " << wx << " " << wy << std::endl;
+ // sleep(1);
}
+ // Add up all coefficients.
+ // for (int i = 0; i <= n; i++) {
+ // wx += pow(x, i) * par[count++];
+ // }
+ // for (int i = 0; i <= n; i++) {
+ // wy += pow(y, i) * par[count++];
+ // }
- delete func;
- delete gr;
+ return w;
}
+float Spline::Spline2DGaus(const Float_t* par) const {
-// Spline extraction Functions
-void Spline::GetCoeff1DTSpline3(int n, double* x, double* y, float* coeff, bool draw) {
+ double Norm = 5.0 + par[1] * 20.0;
+ double Tilt = par[2] * 10.0; //vals[kPosTilt];
+ double Pq0 = 1.0 + par[3] * 1.0; //vals[kPosPq0];
+ double Wq0 = 0.5 + par[4] * 1.0; //vals[kPosWq0];
+ double Pq3 = 1.0 + par[5] * 1.0; //vals[kPosPq3];
+ double Wq3 = 0.5 + par[6] * 1.0; //vals[kPosWq3];
+ double q0 = (fVal[0] - fValMin[0]) / (fValMax[0] - fValMin[0]);
+ double q3 = (fVal[1] - fValMin[1]) / (fValMax[1] - fValMin[1]);
- StopTalking();
- TSpline3 temp_spline = TSpline3("temp_spline", x, y, n);
- StartTalking();
+ // double w = Norm;
+ // w *= TMath::Gaus(q0, Pq0, Wq0);
+ // w *= TMath::Gaus(q3, Pq3*sin(Tilt), Wq3);
- for (size_t i = 0; i < (UInt_t)n; i++) {
+ double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
+ a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);
- double a, b, c, d, e;
- temp_spline.GetCoeff(i, a, b, c, d, e);
+ double b = - sin(2 * Tilt) / (4 * Wq0 * Wq0);
+ b += sin(2 * Tilt) / (4 * Wq3 * Wq3);
- coeff[i * 4] = y[i];
- coeff[i * 4 + 1] = c;
- coeff[i * 4 + 2] = d;
- coeff[i * 4 + 3] = e;
- }
+ double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
+ c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);
- if (draw) {
- TGraph* gr = new TGraph(n, x, y);
- temp_spline.Draw("CA");
- gr->Draw("PL SAME");
- gPad->Update();
- gPad->SaveAs(("plot_test_" + fName + ".pdf").c_str());
- // sleep(3);
- delete gr;
- }
+ double w = Norm;
+ w *= exp(-a * (q0 - Pq0) * (q0 - Pq0));
+ w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3));
+ w *= exp(-c * (q3 - Pq3) * (q3 - Pq3));
+
+ // Norm = par[7];
+ // Tilt = par[8]*50.0;//vals[kPosTilt];
+ // Pq0 = par[9]*10.0;//vals[kPosPq0];
+ // Wq0 = par[10]*10.0;//vals[kPosWq0];
+ // Pq3 = par[11]*10.0; //vals[kPosPq3];
+ // Wq3 = par[12]*10.0; //vals[kPosWq3];
+
+ // // double w = Norm;
+ // // w *= TMath::Gaus(q0, Pq0, Wq0);
+ // // w *= TMath::Gaus(q3, Pq3*sin(Tilt), Wq3);
- return;
+ // a = cos(Tilt)*cos(Tilt) / (2*Wq0*Wq0);
+ // a += sin(Tilt)*sin(Tilt) / (2*Wq3*Wq3);
+
+ // b = - sin(2*Tilt) / (4*Wq0*Wq0);
+ // b += sin(2*Tilt) / (4*Wq3*Wq3);
+
+ // c = sin(Tilt)*sin(Tilt) / (2*Wq0*Wq0);
+ // c += cos(Tilt)*cos(Tilt) / (2*Wq3*Wq3);
+
+ // double w2 = Norm;
+ // w2 *= exp(-a * (q0 - Pq0)*(q0 - Pq0));
+ // w2 *= exp(+2.0*b * (q0 - Pq0)*(q3 - Pq3));
+ // w2 *= exp(-c * (q3 - Pq3)*(q3 - Pq3));
+
+
+// std::cout << "GausWeight = " << w << std::endl;
+ return w;
}
+
+float Spline::Spline1DTSpline3(const Float_t* par) const {
+
+ // Find matching point
+ iter_low = fXScan.begin();
+ iter_high = fXScan.begin();
+ iter_high++;
+ off = 0;
+ fX = fVal[0];
+
+ while ( iter_high != fXScan.end() and
+ (fX < (*iter_low) or fX >= (*iter_high)) ) {
+ off += 4;
+ iter_low++;
+ iter_high++;
+ }
+
+ // std::cout << "Evaluting TSpline3 at " << fX << " " << (*iter_low) << std::endl;
+ // sleep(1.0);
+ float dx = fX - (*iter_low);
+ float weight = (par[off] + dx * (par[off + 1] + dx * (par[off + 2] + dx * par[off + 3])));
+
+ return weight;
+};
+
+
+// 2D Functions
+// ----------------------------------------------
+
+
+
+
+
+
+
diff --git a/src/Splines/Spline.h b/src/Splines/Spline.h
index 3936eec..82e4d90 100644
--- a/src/Splines/Spline.h
+++ b/src/Splines/Spline.h
@@ -1,101 +1,163 @@
#ifndef SPLINE_H
#define SPLINE_H
#include <vector>
#include "TObject.h"
#include "FitParameters.h"
// #include "PlotUtils.h"
// #include "FitUtils.h"
#include "stdint.h"
#include "stdlib.h"
#include "TCanvas.h"
#include <list>
#include "TF1.h"
#include "TSpline.h"
+#include "SplineUtils.h"
+#include "TGraph2D.h"
+#include "TF2.h"
+
+#include "Math/Minimizer.h"
+#include "Math/Factory.h"
+#include "Math/Functor.h"
+#include "TH1D.h"
// Spline Class
class Spline {
public:
- Spline(std::string splname, std::string form, std::vector<double> x);
+ Spline(std::string splname, std::string form, std::string points);
~Spline() {};
void Setup(int type, int ndim, int npar);
double operator()(const Double_t* x, const Double_t* par) const;
+ // double operator()(const Double_t* x, const Double_t* y, const Double_t* par) const;
+
// double DoEval(const Double_t* x, const Double_t* par) const;
// double DoEval(const Double_t* par, bool checkresponse=true) const;
float operator()(const Float_t* x, const Float_t* par) const;
+
float DoEval(const Float_t* x, const Float_t* par) const;
- float DoEval(const Float_t* par, bool checkresponse=true) const;
-
+ float DoEval(const Float_t* par, bool checkresponse = true) const;
+
// void FitCoeff(int n, double* x, double* y, double* par, bool draw);
- void FitCoeff(int n, double* x, double* y, float* par, bool draw);
+ void FitCoeff(std::vector< std::vector<double> > v, std::vector<double> w, float* coeff, bool draw);
inline std::string GetName(void) { return fName; };
inline int GetNDim(void) { return fNDim; };
inline int GetType(void) { return fType; };
inline int GetNPar(void) { return fNPar; };
inline std::string GetForm() {return fForm;};
//void Reconfigure(double x);
- void Reconfigure(float x);
+ void Reconfigure(float x, int index = 0);
+ void Reconfigure(std::string name, float x);
+
std::string fName;
- std::vector<std::string> fSplineNames;
int fType;
int fNDim;
int fNPar;
std::string fForm;
+ std::string fPoints;
bool fOutsideLimits;
+ std::vector<std::string> fSplitNames;
+ std::vector<std::string> fSplitPoints;
+
+ mutable std::vector<float> fVal;
+ mutable std::vector<float> fValMin;
+ mutable std::vector<float> fValMax;
+
+ mutable std::vector< std::vector<float> > fSplitScan;
+
+
+
+
mutable std::vector<float> fXScan;
mutable float fX;
mutable float fXMin;
mutable float fXMax;
+
+ mutable std::vector<float> fYScan;
+ mutable float fY;
+ mutable float fYMin;
+ mutable float fYMax;
+
int fSplineOffset;
+ // TSpline3 Objects.
mutable std::vector<float>::iterator iter_low;
mutable std::vector<float>::iterator iter_high;
mutable int off;
- // Spline List
- enum spline_types {
- k1DPol1 = 1,
- k1DPol2,
- k1DPol3,
- k1DPol4,
- k1DPol5,
- k1DPol6,
- k1DPol1C,
- k1DPol2C,
- k1DPol3C,
- k1DPol4C,
- k1DPol5C,
- k1DPol5C_LX,
- k1DPol10,
- k1DTSpline3,
- k1DPol25,
- k2DPol6
- };
+ // Create a new function for fitting.
+ ROOT::Math::Minimizer* minimizer;
// Available Spline Functions
float Spline1DPol1(const Float_t* par) const;
float Spline1DPol2(const Float_t* par) const;
float Spline1DPol3(const Float_t* par) const;
float Spline1DPol4(const Float_t* par) const;
float Spline1DPol5(const Float_t* par) const;
float Spline1DPol6(const Float_t* par) const;
+ float Spline2DPol(const Float_t* par, int n) const;
+ float Spline2DGaus(const Float_t* par) const;
+float GetMonomial(int p) const;
float Spline1DTSpline3(const Float_t* par) const;
+float Spline2DTSpline3(const Float_t* par) const;
- // Available Fitting Functions
- void FitCoeff1DGraph(int n, double* x, double* y, float* coeff, bool draw);
- void GetCoeff1DTSpline3(int n, double* x, double* y, float* coeff, bool draw);
+// void FitCoeff2DGraph(std::vector< std::vector<double> > v, std::vector<double> w, float* coeff, bool draw);
+};
+// class SplineFCN {
+// public:
+
+// SplineFCN(Spline* spl, std::vector<std::vector<double> > v, std::vector<double> w) { fSpl = spl; fVal = v; fWeight = w; };
+// ~SplineFCN() {};
+
+// double operator()(const double* x) const;
+// double DoEval(const double *x) const;
+// void SaveAs(std::string name, const double* x);
+
+// std::vector< std::vector<double> > fVal;
+// std::vector< double > fWeight;
+// Spline* fSpl;
+
+// };
+
+namespace SplineUtils {
+// // Available Fitting Functions
+// void FitCoeff1DGraph(Spline* spl, int n, double* x, double* y, float* coeff, bool draw);
+// void GetCoeff1DTSpline3(Spline* spl, int n, double* x, double* y, float* coeff, bool draw);
+// // void FitCoeff2DGraph(Spline* spl, int n, double* x, double* y, double* z, float* coeff, bool draw);
+
+// Spline List
+enum spline_types {
+ k1DPol1 = 1,
+ k1DPol2,
+ k1DPol3,
+ k1DPol4,
+ k1DPol5,
+ k1DPol6,
+ k1DPol1C,
+ k1DPol2C,
+ k1DPol3C,
+ k1DPol4C,
+ k1DPol5C,
+ k1DPol5C_LX,
+ k1DPol10,
+ k1DTSpline3,
+ k1DPol25,
+ k2DPol6,
+ k2DGaus,
+ k2DTSpline3
};
+}
+
#endif
diff --git a/src/Splines/SplineReader.cxx b/src/Splines/SplineReader.cxx
index 844677c..dca7ebf 100644
--- a/src/Splines/SplineReader.cxx
+++ b/src/Splines/SplineReader.cxx
@@ -1,135 +1,135 @@
#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;
// Add the spline to the list of all forms
- fAllSplines.push_back( Spline(splname, form, GeneralUtils::ParseToDbl(points, ",")) );
+ fAllSplines.push_back( Spline(splname, form, 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
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));
}
// Loop through and add splines from read type.
for (size_t i = 0; i < fSpline.size(); i++) {
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], ",")) );
+ fAllSplines.push_back( Spline(fSpline[i], fForm[i], fPoints[i]) );
}
}
void SplineReader::Reconfigure(std::map< std::string, double >& vals) {
+ // std::cout << "NEW SPLINE READER =========" << std::endl;
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;
+ // 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;
- fAllSplines[i].Reconfigure(iter->second);
- }
+ // std::cout << "Passing 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;
+ fAllSplines[i].Reconfigure(iter->first, 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;
// Form Offset list
std::vector<int> offsets;
for (size_t i = 0; i < fAllSplines.size(); i++) {
offsets.push_back(off);
off += fAllSplines[i].GetNPar();
}
// #pragma omp parrallel for
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[offsets[i]] );
rw_weight *= w;
// std::cout << "Spline RW Weight = " << rw_weight << " " << w << std::endl;
}
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 (size_t i = 0; i < fAllSplines.size(); i++){
n += fAllSplines[i].GetNPar();
}
return n;
}
diff --git a/src/Splines/SplineUtils.cxx b/src/Splines/SplineUtils.cxx
new file mode 100644
index 0000000..0c723ba
--- /dev/null
+++ b/src/Splines/SplineUtils.cxx
@@ -0,0 +1,111 @@
+#include "SplineUtils.h"
+
+
+// std::vector<int> SplineUtils::GetSplitDialPositions(FitWeight* rw, std::string names) {
+// std::vector<std::string> splitnames = GeneralUtils::ParseToStr(names, ";");
+// std::vector<int> splitpos;
+
+// for (size_t i = 0; i < splitnames.size(); i++) {
+// int pos = fRW->GetDialPos(splitnames[i]);
+// splitpos.push_back(pos);
+// }
+
+// return splitpos;
+// }
+
+std::vector< std::vector<double> > SplineUtils::GetSplitDialPoints(std::string points) {
+
+ // Determine Type
+ std::vector<std::string> splittype = GeneralUtils::ParseToStr(points, ":");
+ std::string deftype = "";
+ if (splittype.size() == 1) {
+ deftype = "GRID";
+ } else if (splittype.size() > 1) {
+ deftype = splittype[0];
+ splittype.erase(splittype.begin());
+ } else {
+ throw;
+ }
+
+ // Make Grid
+ std::vector<std::vector<double> > gridpoints;
+
+ // Possible Point Types
+ if (!deftype.compare("GRID")) {
+
+ std::vector<std::string> splitpoints = GeneralUtils::ParseToStr(splittype[0], ";");
+
+ // CBA writing an ND iterator to build this grids. 1D, 2D, and 3D are below..
+ if (splitpoints.size() == 1) {
+ std::vector<double> tempcont = std::vector<double>(1, 0.0);
+ std::vector<double> tempsplit = GeneralUtils::ParseToDbl(splitpoints[0], ",");
+
+ for (size_t i = 0; i < tempsplit.size(); i++) {
+ tempcont[0] = tempsplit[i];
+ gridpoints.push_back(tempcont);
+ }
+
+ // 2D
+ } else if (splitpoints.size() == 2) {
+
+ std::vector<double> tempcont = std::vector<double>(2, 0.0);
+ std::vector<double> tempsplit1 = GeneralUtils::ParseToDbl(splitpoints[0],",");
+ std::vector<double> tempsplit2 = GeneralUtils::ParseToDbl(splitpoints[1],",");
+
+ for (size_t i = 0; i < tempsplit1.size(); i++) {
+ for (size_t j = 0; j < tempsplit2.size(); j++) {
+ tempcont[0] = tempsplit1[i];
+ tempcont[1] = tempsplit2[j];
+ gridpoints.push_back(tempcont);
+ }
+ }
+ // 3D
+ } else if (splitpoints.size() == 3) {
+
+ std::vector<double> tempcont = std::vector<double>(3, 0.0);
+ std::vector<double> tempsplit1 = GeneralUtils::ParseToDbl(splitpoints[0],",");
+ std::vector<double> tempsplit2 = GeneralUtils::ParseToDbl(splitpoints[1],",");
+ std::vector<double> tempsplit3 = GeneralUtils::ParseToDbl(splitpoints[2],",");
+
+ for (size_t i = 0; i < tempsplit1.size(); i++) {
+ for (size_t j = 0; j < tempsplit2.size(); j++) {
+ for (size_t k = 0; k < tempsplit3.size(); k++) {
+ tempcont[0] = tempsplit1[i];
+ tempcont[1] = tempsplit2[j];
+ tempcont[2] = tempsplit3[j];
+ gridpoints.push_back(tempcont);
+ }
+ }
+ }
+ }
+ } else if (!deftype.compare("DIAG2D")) {
+ std::vector<std::string> splitpoints = GeneralUtils::ParseToStr(splittype[0], ";");
+
+ std::vector<double> tempcont = std::vector<double>(2, 0.0);
+ std::vector<double> tempsplit1 = GeneralUtils::ParseToDbl(splitpoints[0],",");
+
+ for (size_t i = 0; i < tempsplit1.size(); i++) {
+ tempcont[0] = tempsplit1[i];
+ tempcont[1] = 0.0;
+ gridpoints.push_back(tempcont);
+ }
+
+ for (size_t i = 0; i < tempsplit1.size(); i++) {
+ tempcont[0] = 0.0;
+ tempcont[1] = tempsplit1[i];
+ gridpoints.push_back(tempcont);
+ }
+ for (size_t i = 0; i < tempsplit1.size(); i++) {
+ tempcont[0] = tempsplit1[i];
+ tempcont[1] = tempsplit1[i];
+ gridpoints.push_back(tempcont);
+ }
+
+ for (size_t i = 0; i < tempsplit1.size(); i++) {
+ tempcont[0] = tempsplit1[i];
+ tempcont[1] = tempsplit1[ tempsplit1.size() - 1 - i];
+ gridpoints.push_back(tempcont);
+ }
+ }
+ return gridpoints;
+}
\ No newline at end of file
diff --git a/src/Splines/SplineUtils.h b/src/Splines/SplineUtils.h
new file mode 100644
index 0000000..859f92d
--- /dev/null
+++ b/src/Splines/SplineUtils.h
@@ -0,0 +1,12 @@
+#ifndef SPLINE_UTILS_H
+#define SPLINE_UTILS_H
+
+//#include "FitWeight.h"
+#include "GeneralUtils.h"
+
+namespace SplineUtils {
+ //std::vector<int> GetSplitDialPositions(FitWeight* rw, std::string names);
+ std::vector< std::vector<double> > GetSplitDialPoints(std::string points);
+};
+
+#endif
\ No newline at end of file
diff --git a/src/Splines/SplineWriter.cxx b/src/Splines/SplineWriter.cxx
index 0af9863..807d308 100644
--- a/src/Splines/SplineWriter.cxx
+++ b/src/Splines/SplineWriter.cxx
@@ -1,279 +1,677 @@
#include "SplineWriter.h"
-
+using namespace SplineUtils;
// 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 (size_t 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);
+ fValList.push_back(std::vector<double>(1, 0.0));
// Loop over all splines.
for (size_t 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]);
+ std::vector<int> pos; // = SplineUtils::GetSplitDialPositions(fRW, fSpline[i]);
+ std::vector<std::string> splitnames = GeneralUtils::ParseToStr(fSpline[i], ";");
+
+ for (size_t j = 0; j < splitnames.size(); j++) {
+ int temppos = fRW->GetDialPos(splitnames[j]);
+ pos.push_back(temppos);
+ }
+
+ std::vector< std::vector<double> > vals = SplineUtils::GetSplitDialPoints(fPoints[i]);
- // Split Points
- std::vector<double> vals = GeneralUtils::ParseToDbl(fPoints[i], ",");
for (size_t j = 0; j < vals.size(); j++) {
- newvals[pos] = vals[j];
+
+ for (size_t k = 0; k < pos.size(); k++) {
+ newvals[ pos[k] ] = vals[j][k];
+ }
fParVect.push_back(newvals);
fValList.push_back(vals[j]);
fWeightList.push_back(1.0);
- fSetIndex.push_back(i+1);
+ fSetIndex.push_back(i + 1);
}
}
// Print out the parameter set
LOG(FIT) << "Parset | Index | Pars --- " << std::endl;
for (size_t i = 0; i < fSetIndex.size(); i++) {
LOG(FIT) << " Set " << i << ". | " << fSetIndex[i] << " | ";
- if (LOG_LEVEL(FIT)){
+ if (LOG_LEVEL(FIT)) {
for (size_t j = 0; j < fParVect[i].size(); j++) {
- std::cout << " " << fParVect[i][j];
+ std::cout << " " << fParVect[i][j];
}
- std::cout << std::endl;
+ std::cout << std::endl;
}
}
}
-void SplineWriter::FitSplinesForEvent(FitEvent* event) {
+void SplineWriter::FitSplinesForEvent(FitEvent* event, TCanvas* fitcanvas, bool saveplot ) {
// Get Starting Weight
fRW->SetAllDials(&fParVect[0][0], fParVect[0].size());
double nomweight = fRW->CalcWeight(event);
event->RWWeight = nomweight;
- if (fDrawSplines){
+ if (fDrawSplines) {
std::cout << "Nominal Spline Weight = " << nomweight << std::endl;
}
// Loop over parameter sets
for (size_t 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){
+ if (weight >= 0.0 and weight < 200) {
// Fill Weight Set
- fWeightList[i] = weight/nomweight;
- if (fDrawSplines) std::cout << "Calculating values from weight set " << i << " " << fParVect[i][0] << " = " << weight << " " << weight/nomweight << std::endl;
+ fWeightList[i] = weight / nomweight;
+ 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;
int n = fAllSplines.size();
std::vector<int> splinecounter;
for (int i = 0; i < n; i++) {
splinecounter.push_back(coeffcount);
coeffcount += fAllSplines[i].GetNPar();
}
- // int parrallel = int(!fDrawSplines);
-
//#pragma omp parallel for if (parrallel)
for (int i = 0; i < n; i++) {
// Store X/Y Vals
- std::vector<double> xvals;
- std::vector<double> yvals;
+ std::vector<std::vector<double> > dialvals;
+ std::vector<double> weightvals;
bool hasresponse = false;
int npar = (fAllSplines[i]).GetNPar();
coeffcount = splinecounter[i];
for (size_t 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 (fSetIndex[j] != i + 1) continue;
+ dialvals.push_back(fValList[j]);
+ weightvals.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);
+ FitCoeff(&fAllSplines[i], dialvals, weightvals, &fCoEffStorer[coeffcount], fDrawSplines);
+ // for (int k = 0; k < npar; k++){
+ // std::cout << "Spline Coeff " << k << " saved as " << fCoEffStorer[coeffcount+k] << std::endl;
+ // }
+
+ // Move spline fitting to inside SplineWriter as thats where it is done.
+ // Make sure that one minimizer is created PER spline.
+ // Make sure that spline has the ability to calculate the difference between a list of weights and its current set.
+ // by fast calculation.
+
} else {
// std::cout << "Spline " << fSpline[i] << " has no response. " << std::endl;
for (int i = 0; i < npar; i++) {
fCoEffStorer[coeffcount + i] = 0.0;
}
}
+ }
+
+
+ // Check overrides
+ // if (saveplot) {
+ // coeffcount = 0;
+
+ // // Make new canvas to save stuff into
+ // if (fitcanvas) delete fitcanvas;
+ // fitcanvas = new TCanvas("c1", "c1", fAllSplines.size() * 400 , 600);
+ // fitcanvas->Divide(fAllSplines.size(), 1);
+
+ // // Clear out points
+ // for (size_t i = 0; i < fAllDrawnGraphs.size(); i++) {
+ // if (fAllDrawnGraphs[i]) delete fAllDrawnGraphs[i];
+ // if (fAllDrawnHists[i]) delete fAllDrawnHists[i];
+ // }
+ // fAllDrawnGraphs.clear();
+ // fAllDrawnHists.clear();
+
+
+ // for (size_t i = 0; i < fAllSplines.size(); i++) {
+
+ // fitcanvas->cd(i + 1);
+
+ // // Store X/Y Vals
+ // std::vector<std::vector<double> > dialvals;
+ // std::vector<double> weightvals;
+ // bool hasresponse = false;
+ // int npar = (fAllSplines[i]).GetNPar();
+
+ // for (size_t j = 0; j < fSetIndex.size(); j++) {
+ // if ((UInt_t)fSetIndex[j] != (UInt_t)i + 1) continue;
+ // dialvals.push_back(fValList[j]);
+ // weightvals.push_back(fWeightList[j] - 0.0);
+ // if (fWeightList[j] != 1.0) hasresponse = true;
+ // }
+
+ // if (hasresponse) {
+
+ // TGraph* gr = new TGraph(dialvals.size(), dialvals, weightvals);
+ // fAllDrawnGraphs.push_back(gr);
+
+ // // Get XMax Min
+ // int n = xvals.size();
+ // double xmin = 99999.9;
+ // double xmax = -99999.9;
+ // for (int j = 0; j < n; j++) {
+ // if (xvals[j] > xmax) xmax = xvals[j];
+ // if (xvals[j] < xmin) xmin = xvals[j];
+ // }
+
+ // TH1D* hist = new TH1D("temp", "temp", 100, xmin, xmax);
+ // fAllDrawnHists.push_back(hist);
+
+ // for (int k = 0; k < 100; k++) {
+ // double xtemp = hist->GetXaxis()->GetBinCenter(k + 1);
+ // fAllSplines[i].Reconfigure(xtemp);
+ // double ytemp = fAllSplines[i].DoEval(&fCoEffStorer[coeffcount]);
+ // hist->SetBinContent(k + 1, ytemp);
+ // }
+
+ // // 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(0.6, 0.6, 0.85, 0.85);
+ // gPad->Update();
+
+ // hist->SetTitle(fSpline[i].c_str());
+ // hist->GetXaxis()->SetTitle("Dial Variation");
+ // hist->GetYaxis()->SetTitle("Event Weight");
+ // fitcanvas->Update();
+ // }
+ // coeffcount += npar;
+ // }
+ // }
+}
+
+
+
+
+// Fitting
+// ----------------------------------------------
+void SplineWriter::FitCoeff(Spline* spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float* coeff, bool draw) {
+
+ std::vector<double> x;
+ std::vector<double> y;
+ std::vector<double> z;
+ for (size_t i = 0; i < v.size(); i++) {
+ x.push_back(v[i][0]);
+ if (v[i].size() > 1) y.push_back(v[i][1]);
+ if (v[i].size() > 2) z.push_back(v[i][2]);
+ }
+
+ switch (spl->GetType()) {
+
+ // Polynominal Graph Fits
+ case k1DPol1:
+ case k1DPol2:
+ case k1DPol3:
+ case k1DPol4:
+ case k1DPol5:
+ case k1DPol6:
+
+
+
+ FitCoeff1DGraph(spl, v.size(), &x[0], &w[0], coeff, draw);
+ break;
+
+ case k1DTSpline3:
+
+ GetCoeff1DTSpline3(spl, x.size(), &x[0], &w[0], coeff, draw);
+ break;
+
+ case k2DPol6:
+ case k2DGaus:
+ case k2DTSpline3:
+ FitCoeff2DGraph(spl, v, w, coeff, draw);
+ break;
+
+ default:
+ break;
+ }
+
+}
+
+void SplineWriter::FitCoeff1DGraph(Spline* spl, 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", spl, -30.0, 30.0, spl->GetNPar());
+ func->SetNpx(400);
+ func->FixParameter(0, 1.0); // Fix so 1.0 at nominal
+
+ // Run the actual spline fit
+ StopTalking();
+
+ // If linear fit with two points
+ if (n == 2 and spl->GetType() == k1DPol1) {
+
+ float m = (y[1] - y[0]) / (x[1] - x[0]);
+ float c = y[0] - (0.0 - x[0]) * m;
+
+ func->SetParameter(0, c);
+ func->SetParameter(1, m);
+
+ } else if (spl->GetType() == k1DPol1) {
+ gr->Fit(func, "WQ");
+ } else {
+ gr->Fit(func, "FMWQ");
+ }
+
+
+ StartTalking();
+
+ for (int i = 0; i < spl->GetNPar(); i++) {
+ coeff[i] = func->GetParameter(i);
+ }
+
+ if (draw) {
+ gr->Draw("APL");
+ gPad->Update();
+ gPad->SaveAs(("plot_test_" + spl->GetName() + ".pdf").c_str());
+ std::cout << "Saving Graph" << std::endl;
+ sleep(3);
+ }
+
+ delete func;
+ delete gr;
+}
+
+double SplineFCN::operator()(const double* x) const {
+ return DoEval(x);
+}
+
+double SplineFCN::DoEval(const double* x) const {
+
+ float* fx = new float[fSpl->GetNPar()];
+ for (int i = 0; i < fSpl->GetNPar(); i++) {
+ fx[i] = x[i];
+ }//
+
+ double tot = 0;
+ for (size_t i = 0; i < fVal.size(); i++) {
+
+ int nonzero = 0;
+ for (size_t j = 0; j < fVal[i].size(); j++) {
+ fSpl->Reconfigure(fVal[i][j], j);
+ // std::cout << "Reconfiguring " << fVal[i][j] << " " << j << std::endl;
+ if (fVal[i][j] != 0.0) nonzero++;
+ }
+ if (uncorrelated and nonzero > 1) continue;
+
+ double w = fSpl->DoEval(fx);
+ double wdif = (w - fWeight[i]);// / (fWeight[i] * 0.05);
+ tot += sqrt(wdif * wdif);
+ }
+
+ // delete fx;
+
+ return tot;
+}
+
+void SplineFCN::UpdateWeights(std::vector<double>& w) {
+ for (int i = 0; i < w.size(); i++) {
+ fWeight[i] = w[i];
+ }
+}
+
+void SplineFCN::SetCorrelated(bool state) {
+ uncorrelated = !state;
+}
+
+
+
+void SplineFCN::SaveAs(std::string name, const float* fx) {
+
+
+ if (fVal[0].size() != 2) {
+ TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);
+ c1->Divide(2, 1);
+ TH1D* histmc = new TH1D("hist", "hist", fVal.size(), 0.0, double(fVal.size()));
+ TH1D* histdt = new TH1D("histdt", "histdt", fVal.size(), 0.0, double(fVal.size()));
- // 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 j = 0; j < n; j++) {
- if (xvals[j] > xmax) xmax = xvals[j];
- if (xvals[j] < xmin) xmin = xvals[j];
+ for (size_t i = 0; i < fVal.size(); i++) {
+ for (size_t j = 0; j < fVal[i].size(); j++) {
+ fSpl->Reconfigure(fVal[i][j], 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);
- fAllSplines[i].Reconfigure(xtemp);
- double ytemp = fAllSplines[i].DoEval(&fCoEffStorer[coeffcount]);
- hist->SetBinContent(k + 1, ytemp);
- // std::cout << "Set Temp " << k << " " << ytemp << std::endl;
+ histmc->SetBinContent( i + 1, fSpl->DoEval(fx) );
+ histdt->SetBinContent( i + 1, fWeight[i] );
+ }
+
+ // histmc->Add(histdt,-1.0);
+
+ c1->cd(1);
+ histmc->SetTitle("Spline;Parameter Set;Weight");
+ histmc->Draw("HIST");
+ histdt->SetLineColor(kRed);
+ histdt->Draw("SAME HIST");
+
+ c1->cd(2);
+ TH1D* histdif = (TH1D*) histmc->Clone();
+ histdif->Add(histdt, -1.0);
+ histdif->Draw("HIST");
+
+ c1->Update();
+ c1->SaveAs(name.c_str());
+ delete c1;
+ } else {
+
+ TGraph2D* histmc = new TGraph2D();
+ TGraph2D* histdt = new TGraph2D();
+ TGraph2D* histdif = new TGraph2D();
+
+ for (size_t i = 0; i < fVal.size(); i++) {
+ for (size_t j = 0; j < fVal[i].size(); j++) {
+ fSpl->Reconfigure(fVal[i][j], j);
}
- // gr->Draw("APL");
- hist->SetLineColor(kRed);
- hist->Draw("HIST C");
- gr->SetMarkerStyle(20);
- gr->Draw("P SAME");
+ histmc->SetPoint(histmc->GetN(), fVal[i][0], fVal[i][1], fSpl->DoEval(fx));
+ histdt->SetPoint(histdt->GetN(), fVal[i][0], fVal[i][1], fWeight[i]);
+ histdif->SetPoint(histdif->GetN(), fVal[i][0], fVal[i][1], fabs(fSpl->DoEval(fx) - fWeight[i]));
- gPad->Update();
- std::cout << "Writing : " << ("F1eval_" + fSpline[i] + ".pdf") << std::endl;
- //gPad->SaveAs(("F1eval_" + fSpline[i] + ".pdf").c_str());
- delete gr;
- // delete f1;
}
+ TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);
+ c1->Divide(3, 1);
+ c1->cd(1);
+ histmc->SetTitle(("Spline;" + fSpl->GetName()).c_str());
+ histmc->Draw("COLZ");
+ gPad->Update();
+ c1->cd(2);
+ histdt->SetTitle(("Raw;" + fSpl->GetName()).c_str());
+ histdt->Draw("COLZ");
+ c1->cd(3);
+ histdif->SetTitle(("Dif;" + fSpl->GetName()).c_str());
+ histdif->Draw("COLZ");
+ gPad->SetRightMargin(0.15);
+ // gPad->SetLogz(1);
+
+ gPad->Update();
+ c1->SaveAs(name.c_str());
+ delete c1;
+ }
- // std::cout << "Finished Calcing splines for " << fAllSplines[i].GetName() << std::endl;
- // coeffcount += npar;
+
+
+ // gPad->SaveAs(name.c_str());
+}
+
+namespace SplineUtils {
+Spline* gSpline = NULL;
+}
+
+double SplineUtils::Func2DWrapper(double* x, double* p) {
+ // std::cout << "2D Func Wrapper " << x[0] << " " << x[1] << std::endl;
+ return (*gSpline)(x, p);
+}
+
+void SplineWriter::FitCoeff2DGraph(Spline* spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float* coeff, bool draw) {
+
+// Make a 2D Function
+ SplineUtils::gSpline = spl;
+ TF2* f2 = new TF2("f2", SplineUtils::Func2DWrapper, -30.0, 30.0, -30.0, 30.0, spl->GetNPar());
+ TGraph2D* histmc = new TGraph2D();
+ for (size_t i = 0; i < v.size(); i++) {
+ histmc->SetPoint(histmc->GetN(), v[i][0], v[i][1], w[i]);
}
+ StopTalking();
+ histmc->Fit(f2,"FMWQ");
+ StartTalking();
+ for (int i = 0; i < spl->GetNPar(); i++) {
+ coeff[i] = f2->GetParameter(i);
+ // std::cout << "Fit 2D Func " << i << " = " << coeff[i] << std::endl;
+ }
- // Check overrides
- if (fDrawSplines) {
- coeffcount = 0;
- for (size_t i = 0; i < fAllSplines.size(); i++) {
+ delete f2;
+ delete histmc;
+ // fSplineFCNs[spl] = new SplineFCN(spl, v, w);
+ // fSplineFCNs[spl]->SaveAs("mysplinetest_" + spl->GetName() + ".pdf", coeff);
+ // sleep(1);
+ // delete fSplineFCNs[spl];
- // Store X/Y Vals
- std::vector<double> xvals;
- std::vector<double> yvals;
- bool hasresponse = false;
- int npar = (fAllSplines[i]).GetNPar();
- for (size_t j = 0; j < fSetIndex.size(); j++) {
- if ((UInt_t)fSetIndex[j] != (UInt_t)i+1) continue;
- xvals.push_back(fValList[j]);
- yvals.push_back(fWeightList[j] - 0.0);
- if (fWeightList[j] != 1.0) hasresponse = true;
- }
+}
+
+
+void SplineWriter::FitCoeffNDGraph(Spline* spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float* coeff, bool draw) {
- 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 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);
- fAllSplines[i].Reconfigure(xtemp);
- double ytemp = fAllSplines[i].DoEval(&fCoEffStorer[coeffcount]);
- hist->SetBinContent(k + 1, ytemp);
-
- // 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;
+ // TGraph2D* gr = new TGraph2D(n, x, y);
+ // double xmin = 99999.9;
+ // double xmax = -99999.9;
+ // double ymin = 99999.9;
+ // double ymax = -99999.9;
+ // for (int i = 0; i < n; i++) {
+ // if (x[i] > xmax) xmax = x[i];
+ // if (x[i] < xmin) xmin = x[i];
+ // if (y[i] > ymax) ymax = y[i];
+ // if (y[i] < ymin) ymin = y[i];
+ // }
+ // double xwidth = xmax - xmin;
+ // xmin = xmin - xwidth * 0.01;
+ // xmax = xmax + xwidth * 0.01;
+ // double ywidth = ymax - ymin;
+ // ymin = ymin - ywidth * 0.01;
+ // ymax = ymax + ywidth * 0.01;
+
+
+ // StopTalking();
+
+ // Build if not already in writer
+ // if (fSplineFCNs.find(spl) != fSplineFCNs.end()) {
+ // fSplineFCNs.erase(spl);
+ // }
+ // fSplineFCNs[spl] = new SplineFCN(spl, v, w);
+
+
+ if (fSplineFunctors.find(spl) != fSplineFunctors.end()) {
+ delete fSplineFunctors[spl];
+ fSplineFunctors.erase(spl);
+ }
+
+ if (fSplineFCNs.find(spl) != fSplineFCNs.end()) {
+ delete fSplineFCNs[spl];
+ fSplineFCNs.erase(spl);
+ }
+
+ if (fSplineMinimizers.find(spl) != fSplineMinimizers.end()) {
+ delete fSplineMinimizers[spl];
+ fSplineMinimizers.erase(spl);
+ }
+
+
+
+ if (fSplineMinimizers.find(spl) == fSplineMinimizers.end()) {
+ std::cout << "Building new ND minimizer for " << spl << std::endl;
+ fSplineFCNs[spl] = new SplineFCN(spl, v, w);
+
+ // fSplineFCNs[spl] = new SplineFCN(spl, v, w);
+ fSplineFunctors[spl] = new ROOT::Math::Functor( *fSplineFCNs[spl], spl->GetNPar() );
+// fitclass = "GSLSimAn"; fittype = "";
+ fSplineMinimizers[spl] = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
+ fSplineMinimizers[spl]->SetMaxFunctionCalls(1E8);
+ fSplineMinimizers[spl]->SetMaxIterations(1E8);
+ fSplineMinimizers[spl]->SetTolerance(1.E-6);
+ fSplineMinimizers[spl]->SetStrategy(2);
+
+ for (int j = 0; j < spl->GetNPar(); j++) {
+ if (j != 0) {
+ fSplineMinimizers[spl]->SetVariable(j, Form("Par_%i", j), 0.1, 0.1);
+ // fSplineMinimizers[spl]->SetParameter(j, Form("Par_%i", j), 0.1, 0.1, -100.0, 100.0 );
+ } else {
+ fSplineMinimizers[spl]->SetVariable(j, Form("Par_%i", j), 1.0, 0.1);
+ // fSplineMinimizers[spl]->SetParameter(j, Form("Par_%i", j), 1.0, 0.1, -100.0, 100.0 );
+ // fSplineMinimizers[spl]->FixVariable(j);
}
- coeffcount += npar;
}
+ // fSplineMinimizers[spl]->SetFunction(*fSplineFunctors[spl]);
+ sleep(1);
}
+ // Create a new function for fitting.
+ // StopTalking();
+
+ // Update FCN
+ fSplineFCNs[spl]->UpdateWeights(w);
+ // fSplineMinimizers[spl]->SetDefaultOptions();
+ // fSplineMinimizers[spl]->Clear();
+ for (int j = 0; j < spl->GetNPar(); j++) {
+ if (j != 0) {
+ fSplineMinimizers[spl]->SetVariableValue(j, 0.1);
+ // fSplineMinimizers[spl]->SetParameter(j, Form("Par_%i", j), 0.1, 0.1, -100.0, 100.0 );
+
+ // fSplineMinimizers[spl]->SetVariableValue(j, 0.1);
+ }
+ }
+ // fSplineFCNs[spl]->SetCorrelated(false);
+ // fSplineFunctors[spl] = new ROOT::Math::Functor( *fSplineFCNs[spl], spl->GetNPar() );
+ // fSplineMinimizers[spl]->SetFunction(*fSplineFunctors[spl]);
+ // fSplineMinimizers[spl]->Minimize();
+
+ fSplineFCNs[spl]->SetCorrelated(true);
+ delete fSplineFunctors[spl];
+ fSplineFunctors[spl] = new ROOT::Math::Functor( *fSplineFCNs[spl], spl->GetNPar() );
+ fSplineMinimizers[spl]->SetFunction(*fSplineFunctors[spl]);
+ // ((TFitterMinuit*)fSplineMinimizers[spl])->CreateMinimizer(TFitterMinuit::kMigrad);
+ fSplineMinimizers[spl]->Minimize();
+ fSplineMinimizers[spl]->Minimize();
+
+ // ((TFitterMinuit*)fSplineMinimizers[spl])->CreateMinimizer(TFitterMinuit::kMigrad);
+ // fSplineMinimizers[spl]->Minimize();
+
+ // delete minimizer;
+ // StartTalking();
+
+ // Now Get Parameters
+ const double *values = fSplineMinimizers[spl]->X();
+ for (int i = 0; i < spl->GetNPar(); i++) {
+
+ std::cout << "Updated Coeff " << i << " " << values[i] << std::endl;
+ coeff[i] = values[i];
+ }
+
+ // Save a sample
+ fSplineFCNs[spl]->SaveAs("mysplinetest_" + spl->GetName() + ".pdf", coeff);
+ sleep(1);
+
+ // delete values;
+ // delete minimizer;
+}
+
+
+
+// Spline extraction Functions
+void SplineWriter::GetCoeff1DTSpline3(Spline* spl, int n, double* x, double* y, float* coeff, bool draw) {
+
+ StopTalking();
+ TSpline3 temp_spline = TSpline3("temp_spline", x, y, n);
+ StartTalking();
+
+ for (size_t i = 0; i < (UInt_t)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;
+
+ // std::cout << "Setting Spline Coeff Set " << i << " to " << y[i] << " " << c << " " << d << " " << e << std::endl;
+ }
+
+ if (draw) {
+ TGraph* gr = new TGraph(n, x, y);
+ temp_spline.Draw("CA");
+ gr->Draw("PL SAME");
+ gPad->Update();
+ gPad->SaveAs(("plot_test_" + spl->GetName() + ".pdf").c_str());
+
+ delete gr;
+ }
+
+ return;
}
+
diff --git a/src/Splines/SplineWriter.h b/src/Splines/SplineWriter.h
index ac2fe86..4cfcdaf 100644
--- a/src/Splines/SplineWriter.h
+++ b/src/Splines/SplineWriter.h
@@ -1,35 +1,82 @@
#ifndef SPLINEWRITER_H
#define SPLINEWRITER_H
#include "FitWeight.h"
#include "Spline.h"
#include "FitParameters.h"
+#include "SplineUtils.h"
+#include "TFitterMinuit.h"
+
+class SplineFCN {
+public:
+
+ SplineFCN(Spline* spl, std::vector<std::vector<double> > v, std::vector<double> w) { fSpl = spl; fVal = v; fWeight = w; };
+ ~SplineFCN() {};
+
+ double operator()(const double* x) const;
+ double DoEval(const double *x) const;
+ void SaveAs(std::string name, const float* fx);
+ void UpdateWeights(std::vector<double>& w);
+ void SetCorrelated(bool state = true);
+
+ bool uncorrelated;
+ std::vector< std::vector<double> > fVal;
+ std::vector< double > fWeight;
+ Spline* fSpl;
+
+};
+
class SplineWriter : public SplineReader {
- public:
- SplineWriter(FitWeight* fw){
+public:
+ SplineWriter(FitWeight* fw) {
fRW = fw;
fDrawSplines = FitPar::Config().GetParB("drawsplines");
};
- ~SplineWriter(){};
+ ~SplineWriter() {};
void SetupSplineSet();
void Write(std::string name);
void AddCoefficientsToTree(TTree* tree);
- void FitSplinesForEvent(FitEvent* event);
-
+ void FitSplinesForEvent(FitEvent* event, TCanvas* fitcanvas = NULL, bool saveplot = false);
int fNCoEff;
// double* fCoEffStorer;
float* fCoEffStorer;
std::vector< std::vector<double> > fParVect;
std::vector< int > fSetIndex;
std::vector< double > fWeightList;
- std::vector< double > fValList;
+ std::vector< std::vector<double> > fValList;
FitWeight* fRW;
bool fDrawSplines;
-
+
+ std::vector<TH1D*> fAllDrawnHists;
+ std::vector<TGraph*> fAllDrawnGraphs;
+
+ std::map<Spline*,SplineFCN*> fSplineFCNs;
+ std::map<Spline*,ROOT::Math::Functor*> fSplineFunctors;
+ std::map<Spline*,ROOT::Math::Minimizer*> fSplineMinimizers;
+
+ // Available Fitting Functions
+void FitCoeff(Spline* spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float* coeff, bool draw);
+void FitCoeff1DGraph(Spline* spl, int n, double* x, double* y, float* coeff, bool draw);
+void GetCoeff1DTSpline3(Spline* spl, int n, double* x, double* y, float* coeff, bool draw);
+void FitCoeff2DGraph(Spline* spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float* coeff, bool draw);
+void FitCoeffNDGraph(Spline* spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float* coeff, bool draw);
+
+
};
+
+namespace SplineUtils {
+
+double Func2DWrapper(double* x, double* p);
+extern Spline* gSpline;
+ // return 1.0;
+// }
+// void FitCoeff2DGraph(Spline* spl, int n, double* x, double* y, double* z, float* coeff, bool draw);
+
+}
+
#endif