diff --git a/src/Routines/SplineRoutines.cxx b/src/Routines/SplineRoutines.cxx
index 996122e..ea67d2d 100755
--- a/src/Routines/SplineRoutines.cxx
+++ b/src/Routines/SplineRoutines.cxx
@@ -1,1616 +1,1637 @@
// 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);
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("GenerateEventWeights")) { GenerateEventWeights(); }
// else if (!rout.compare("BuildEventSplines")) { BuildEventSplines(); }
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();
// Make an ugly list for N cores
int ncores = FitPar::Config().GetParI("NCORES");//omp_get_max_threads();
std::vector<SplineWriter*> splwriterlist;
for (int i = 0; i < ncores; i++){
SplineWriter* tmpwriter = new SplineWriter(fRW);
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
tmpwriter->AddSpline(splkey);
}
tmpwriter->SetupSplineSet();
splwriterlist.push_back(tmpwriter);
}
// 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);
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* weighttree = new TTree("weight_tree", "weight_tree");
splwrite->AddWeightsToTree(weighttree);
// Make container for all weights
- int nweights = splwrite->GetNWeights()+1;
+ int nweights = splwrite->GetNWeights();
double** weightcont = new double*[nevents];
for (int k = 0; k < nevents; k++) {
weightcont[k] = new double[nweights];
}
int npar = splwrite->GetNPars();
int lasttime = time(NULL);
// Could reorder this to save the weightconts in order instead of reconfiguring per event.
// 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);
splwrite->GetWeightsForEvent(nuisevent, weightcont[i]);
+ bool hasresponse = false;
+ for (int j = 0; j < nweights; j++){
+
+ if (weightcont[i][j] != 1.0){
+ // std::cout << "Non Zero Weight at " << i << " " << j << std::endl;
+ hasresponse = true;
+ } else {
+ // std::cout << "Empty Weight at " << i << " " << j << std::endl;
+ }
+ }
+ if (!hasresponse){
+ // std::cout << "Deleting flat response " << nuisevent->Mode << std::endl;
+ delete weightcont[i];
+ weightcont[i] = NULL;
+ }
// Save everything
eventtree->Fill();
weighttree->Fill();
// splinetree->Fill();
// nuisevent->Print();
// std::cout << "Done with event " << i << std::endl;
// Push weight sets into a the array
// 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 weights. " << timestring.str() << std::endl;
}
// Iterate
i++;
nuisevent = input->NextNuisanceEvent();
}
outputfile->cd();
eventtree->Write();
weighttree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
outputfile->Close();
outputfile = new TFile(outputfilename.c_str(), "UPDATE");
outputfile->cd();
weighttree = (TTree*) outputfile->Get("weight_tree");
splwrite->ReadWeightsFromTree(weighttree);
// Divide weights container into Ncores.
// Parrallelise this loop checking for what core we are on.
// for (int i = 0; i < nevents; i++){
// splwriterlist[int(i / (nevents/4))]->FitSplinesForEvent(coeff);
// }
// // Now loop over weights tree
// for (int i = 0; i < weighttree->GetEntries(); i++) {
// weighttree->GetEntry(i);
// splwrite->FitSplinesForEvent();
// splinetree->Fill();
// 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) << "Built " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
// }
// }
// Get Splines
float** allcoeff = new float*[nevents];
for (int k = 0; k < nevents; k++){
allcoeff[k] = new float[npar];
}
#pragma omp parallel for num_threads(ncores)
for (int i = 0; i < nevents; i++) {
//#pragma omp atomic
// printf("Using Thread %d to build event %d \n", int(omp_get_thread_num()), (int)i );
// std::cout<< " -> Writer = " << splwriterlist[ i / (nevents/ncores) ] << std::endl;
// #pragma omp atomic
- splwriterlist[ int(omp_get_thread_num()) ]->FitSplinesForEvent(weightcont[i], allcoeff[i]);
+ if (weightcont[i]){
+ splwriterlist[ int(omp_get_thread_num()) ]->FitSplinesForEvent(weightcont[i], allcoeff[i]);
+ } else {
+ for (int j = 0; j < npar; j++){
+ allcoeff[i][j] = float(0.0);
+ }
+ }
// splwrite->FitSplinesForEvent(weightcont[i], allcoeff[i]);
if (i % 500 == 0) {
if (LOG_LEVEL(REC)){
printf("Using Thread %d to build event %d \n", int(omp_get_thread_num()), (int)i );
}
}
/*
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.";
timestring << " Using Writer at " << i / (nevents/ncores) << " = " << splwriterlist[ i / (nevents/ncores) ] << std::endl;
}
LOG(REC) << "Built " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
}
*/
}
// Save Splines into TTree
float* coeff = new float[npar];
outputfile->cd();
TTree* splinetree = new TTree("spline_tree", "spline_tree");
splinetree->Branch("SplineCoeff", coeff, Form("SplineCoeff[%d]/F",npar));
std::cout << "Saving to the allcoeff" << std::endl;
for (int k = 0; k < nevents; k++) {
for (int l = 0; l < npar; l++) {
coeff[l] = allcoeff[k][l];
}
std::cout << "Coeff 0, 1, 2 = " << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
splinetree->Fill();
}
// Save flux and close file
outputfile->cd();
splinetree->Write();
// Delete the container.
for (int k = 0; k < nevents; k++) {
delete weightcont[k];
}
delete weightcont;
delete coeff;
// 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());
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;
}
// 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));
}
// 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->GetWeightsForEvent(nuisevent);
splwrite->FitSplinesForEvent(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/Splines/SplineWriter.cxx b/src/Splines/SplineWriter.cxx
index 5055b47..d4f8be8 100644
--- a/src/Splines/SplineWriter.cxx
+++ b/src/Splines/SplineWriter.cxx
@@ -1,765 +1,766 @@
#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::AddWeightsToTree(TTree* tr) {
// Add only the fitted spline coefficients to the ttree
std::cout << "Saving Spline Weights to TTree = " << Form("SplineWeights[%d]/F", (int)fValList.size()) << std::endl;
// sleep(1);
tr->Branch("SplineWeights", fWeightList, Form("SplineWeights[%d]/D", (int)fValList.size()));
}
void SplineWriter::ReadWeightsFromTree(TTree* tr) {
tr->SetBranchAddress("SplineWeights", fWeightList);
}
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(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;
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]);
for (size_t j = 0; j < vals.size(); 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);
}
}
fWeightList = new double[fValList.size()];
for (int i = 0; i < fValList.size(); i++) {
fWeightList[i] = 1.0;
}
// 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)) {
for (size_t j = 0; j < fParVect[i].size(); j++) {
std::cout << " " << fParVect[i][j];
}
std::cout << std::endl;
}
}
}
void SplineWriter::GetWeightsForEvent(FitEvent* event) {
// Get Starting Weight
fRW->SetAllDials(&fParVect[0][0], fParVect[0].size());
double nomweight = fRW->CalcWeight(event);
event->RWWeight = nomweight;
if (fDrawSplines) {
std::cout << "Nominal Spline Weight = " << nomweight << std::endl;
}
fWeightList[0] = nomweight;
// 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) {
// 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;
} else {
fWeightList[i] = 1.0;
}
}
}
void SplineWriter::GetWeightsForEvent(FitEvent* event, double* weights) {
// Get Starting Weight
fRW->SetAllDials(&fParVect[0][0], fParVect[0].size());
double nomweight = fRW->CalcWeight(event);
event->RWWeight = nomweight;
if (fDrawSplines) {
std::cout << "Nominal Spline Weight = " << nomweight << std::endl;
}
weights[0] = nomweight;
// 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) {
// Fill Weight Set
weights[i] = weight / nomweight;
if (fDrawSplines) std::cout << "Calculating values from weight set " << i << " " << fParVect[i][0] << " = " << weight << " " << weight / nomweight << std::endl;
} else {
weights[i] = 1.0;
}
+ fWeightList[i] = weights[i];
}
}
void SplineWriter::FitSplinesForEvent(double* inputweights, float* coeff) {
int n = fAllSplines.size();
int coeffcount = 0;
for (int i = 0; i < n; i++) {
// DialVals
std::vector<std::vector<double> > dialvals;
std::vector<double> weightvals;
bool hasresponse = false;
for (size_t j = 0; j < fSetIndex.size(); j++) {
if (fSetIndex[j] != i + 1) continue;
dialvals.push_back(fValList[j]);
double tempw = inputweights[j];
weightvals.push_back(tempw);
if (tempw != 1.0) hasresponse = true;
}
// Perform Fit
if (hasresponse) {
// std::cout << "Fitting Coeff" << std::endl;
FitCoeff(&fAllSplines[i], dialvals, weightvals, &coeff[coeffcount], fDrawSplines);
} else {
for (int j = 0; coeffcount + j < fNCoEff; j++) {
// std::cout << "Setting 0.0 response " << coeffcount + i << " " << fNCoEff << std::endl;
coeff[coeffcount + j] = 0.0;
}
}
// std::cout << "Adding coeffcount" << std::endl;
// Offset coeffcount
coeffcount += (fAllSplines[i]).GetNPar();
}
// std::cout << "FitEvent" << std::endl;
return;
}
void SplineWriter::FitSplinesForEvent(TCanvas * fitcanvas, bool saveplot) {
// 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();
}
// #pragma omp parallel for
for (int i = 0; i < n; i++) {
// Store X/Y Vals
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;
dialvals.push_back(fValList[j]);
double tempw = fWeightList[j];
weightvals.push_back(tempw);
if (tempw != 1.0) hasresponse = true;
}
// Make a new graph and fit coeff if response
if (hasresponse) {
FitCoeff(&fAllSplines[i], dialvals, weightvals, &fCoEffStorer[coeffcount], fDrawSplines);
} else {
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.size(), &x[0], &y[0], &w[0], coeff, draw);
break;
default:
break;
}
// fSplineFCNs[spl] = new SplineFCN(spl, v, w);
// fSplineFCNs[spl]->SaveAs("mysplinetest_" + spl->GetName() + ".pdf", coeff);
// sleep(1);
// delete fSplineFCNs[spl];
}
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()));
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);
}
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);
}
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]));
}
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;
}
// 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, int n, double * x, double * y, double * w, float * coeff, bool draw) {
#pragma omp critical
{
// if (gSpline != spl) gSpline = spl;
if (SplineUtils::gSpline != spl)
SplineUtils::gSpline = spl;
TF2* f2 = new TF2("f2", Func2DWrapper, -30.0, 30.0, -30.0, 30.0, spl->GetNPar());
TGraph2D* histmc = new TGraph2D(n, x, y, w);
f2->SetNpx(400);
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;
}
delete f2;
delete histmc;
}
}
void SplineWriter::FitCoeffNDGraph(Spline * spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float * coeff, bool draw) {
// 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);
}
}
// 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;
}