Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8310440
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
200 KB
Subscribers
None
View Options
diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
index 0d2d9b0..b6eecf2 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,1032 +1,1034 @@
#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);
std::vector<nuiskey> covarkeys = Config::QueryKeys("covar");
LoadPulls(covarkeys);
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! " << std::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!" << std::endl;
throw;
}
fIterationTree->Write();
}
//***************************************************
void JointFCN::FillIterationTree(FitWeight* rw) {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << std::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" << std::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 << std::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 : " << 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"
<< std::endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
}
void JointFCN::LoadPulls(std::vector<nuiskey> pullkeys) {
for (size_t i = 0; i < pullkeys.size(); i++) {
nuiskey key = pullkeys[i];
std::string pullname = key.GetS("name");
std::string pullfile = key.GetS("input");
std::string pulltype = key.GetS("type");
fOutputDir->cd();
+ std::cout << "Creating Pull Term : " << std::endl;
+ sleep(1);
fPulls.push_back(new ParamPull(pullname, pullfile, pulltype));
}
// // 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));
// }
}
// //***************************************************
// 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"
// << std::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 << std::endl;
// std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled << std::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" << std::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();
curinput->CreateCache();
int i = 0;
int nevents = curinput->GetNEvents();
int countwidth = nevents / 5;
// 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;
// std::cout << "RWWeight = " << curevent->RWWeight << " " << curevent->InputWeight << std::endl;
// Logging
// std::cout << CHECKLOG(1) << std::endl;
if (LOGGING(REC)) {
if (i % countwidth == 0) {
QLOG(REC, curinput->GetName() << " : Processed " << i
<< " events. [M, W] = ["
<< curevent->Mode << ", "
<< rwweight << "]");
}
}
// 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(curevent->Weight);
// 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 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++;
}
curinput->RemoveCache();
// 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 << " " << 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() {
//***************************************************
LOG(FIT) << " -> Doing FAST using manager" << std::endl;
// 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()) {
ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl;
ReconfigureUsingManager();
return;
}
bool fFillNuisanceEvent = FitPar::Config().GetParB("FullEventOnSignalReconfigure");
// 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 / 20;
// If All Splines tell splines they need a reconfigure.
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);
}
}
// Loop over all possible spline inputs
double* coreeventweights = new double[fSignalEventBoxes.size()];
splinecount = 0;
inp_iter = fInputList.begin();
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
// Loop over all signal flags
// For each valid signal flag add one to splinecount
// Get Splines from that count and add to weight
// Add splinecount
int sigcount = 0;
splinecount = 0;
// #pragma omp parallel for shared(splinecount,sigcount)
for (uint iinput = 0; iinput < fInputList.size(); iinput++) {
InputHandlerBase* curinput = fInputList[iinput];
BaseFitEvt* curevent = curinput->FirstBaseEvent();
for (int i = 0; i < curinput->GetNEvents(); i++) {
double rwweight = 0.0;
if (fSignalEventFlags[sigcount]) {
// Get Event Info
if (!fIsAllSplines) {
if (fFillNuisanceEvent) curinput->GetNuisanceEvent(i);
else curevent = curinput->GetBaseEvent(i);
} else {
curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0];
}
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
rwweight = curevent->Weight;
coreeventweights[splinecount] = rwweight;
if (splinecount % countwidth == 0) {
LOG(REC) << "Processed " << splinecount << " event weights. W = " << rwweight << std::endl;
}
// #pragma omp atomic
splinecount++;
}
// #pragma omp atomic
sigcount++;
}
}
LOG(SAM) << "Processed event weights." << std::endl;
// #pragma omp barrier
// Reset Iterators
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
box_iter = fSignalEventBoxes.begin();
samsig_iter = fSampleSignalFlags.begin();
int nsplineweights = splinecount;
splinecount = 0;
// Start of Fast Event Loop ============================
// Start input iterators
// Loop over number of inputs
for (int ispline = 0; ispline < nsplineweights; ispline++) {
double rwweight = coreeventweights[ispline];
// Get iterators for this event
std::vector<bool>::iterator subsamsig_iter = (*samsig_iter).begin();
std::vector<MeasurementVariableBox*>::iterator subbox_iter = (*box_iter).begin();
// Loop over all sub measurements.
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
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 (ispline % countwidth == 0) {
LOG(REC) << "Filled " << ispline << " sample weights." << std::endl;
}
// Iterate over the main signal event containers.
samsig_iter++;
box_iter++;
spline_iter++;
splinecount++;
}
// End of Fast Event Loop ===================
LOG(SAM) << "Filled sample distributions." << std::endl;
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Cleanup coreeventweights
if (fIsAllSplines) {
delete coreeventweights;
}
// 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..." << std::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 << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->SetFakeDataValues(fakeinput);
}
return;
}
diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx
index 651097e..25c5c54 100644
--- a/src/FitBase/Measurement1D.cxx
+++ b/src/FitBase/Measurement1D.cxx
@@ -1,1798 +1,1798 @@
// Copyright 2016 L. Pickering, P caltowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This ile 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 "Measurement1D.h"
//********************************************************************
Measurement1D::Measurement1D(void) {
//********************************************************************
// XSec Scalings
fScaleFactor = -1.0;
fCurrentNorm = 1.0;
// Histograms
fDataHist = NULL;
fDataTrue = NULL;
fMCHist = NULL;
fMCFine = NULL;
fMCWeighted = NULL;
fMaskHist = NULL;
// Covar
covar = NULL;
fFullCovar = NULL;
fCovar = NULL;
fInvert = NULL;
fDecomp = NULL;
// Fake Data
fFakeDataInput = "";
fFakeDataFile = NULL;
// Options
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsNoWidth = false;
fIsDifXSec = false;
fIsEnu1D = false;
// Inputs
fInput = NULL;
fRW = NULL;
// Extra Histograms
fMCHist_Modes = NULL;
}
//********************************************************************
Measurement1D::~Measurement1D(void) {
//********************************************************************
if (fDataHist) delete fDataHist;
if (fDataTrue) delete fDataTrue;
if (fMCHist) delete fMCHist;
if (fMCFine) delete fMCFine;
if (fMCWeighted) delete fMCWeighted;
if (fMaskHist) delete fMaskHist;
if (covar) delete covar;
if (fFullCovar) delete fFullCovar;
if (fCovar) delete fCovar;
if (fInvert) delete fInvert;
if (fDecomp) delete fDecomp;
}
//********************************************************************
void Measurement1D::FinaliseSampleSettings() {
//********************************************************************
MeasurementBase::FinaliseSampleSettings();
// Setup naming + renaming
fName = fSettings.GetName();
fSettings.SetS("originalname", fName);
if (fSettings.Has("rename")) {
fName = fSettings.GetS("rename");
fSettings.SetS("name", fName);
}
// Setup all other options
LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl;
if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but using poisson likelihoods."
<< std::endl;
}
if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
<< "not flux averaged!" << std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
if (!fRW) fRW = FitBase::GetRW();
if (!fInput and !fIsJoint) SetupInputs(fSettings.GetS("input"));
// Setup options
SetFitOptions(fDefaultTypes); // defaults
SetFitOptions(fSettings.GetS("type")); // user specified
EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));
if (fAddNormPen) {
if (fNormError <= 0.0) {
ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl;
ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl;
throw;
}
}
}
//********************************************************************
void Measurement1D::CreateDataHistogram(int dimx, double* binx) {
//********************************************************************
if (fDataHist) delete fDataHist;
fDataHist = new TH1D( (fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(),
dimx, binx) ;
}
//********************************************************************
void Measurement1D::SetDataFromTextFile(std::string datafile) {
//********************************************************************
LOG(SAM) << "Reading data from text file: " << datafile << std::endl;
fDataHist = PlotUtils::GetTH1DFromFile(datafile,
fSettings.GetName() + "_data",
fSettings.GetFullTitles());
}
//********************************************************************
void Measurement1D::SetDataFromRootFile(std::string datafile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname);
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
return;
};
//********************************************************************
void Measurement1D::SetEmptyData(){
//********************************************************************
fDataHist = new TH1D("EMPTY_DATA","EMPTY_DATA",1,0.0,1.0);
}
//********************************************************************
void Measurement1D::SetPoissonErrors() {
//********************************************************************
if (!fDataHist) {
ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl;
ERR(FTL) << "Setup Data First!" << std::endl;
throw;
}
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
}
}
//********************************************************************
void Measurement1D::SetCovarFromDiagonal(TH1D* data) {
//********************************************************************
if (!data and fDataHist) {
data = fDataHist;
}
if (data) {
LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl;
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
} else {
ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl;
}
// if (!fIsDiag) {
// ERR(FTL) << "SetCovarMatrixFromDiag called for measurement "
// << "that is not set as diagonal." << std::endl;
// throw;
// }
}
//********************************************************************
void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl;
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl;
fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl;
covar = StatUtils::GetCovarFromTextFile(covfile, dim);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl;
covar = StatUtils::GetCovarFromRootFile(covfile, histname);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) dim = fDataHist->GetNbinsX();
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void Measurement1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromRootFile(covfile, histname);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void Measurement1D::ScaleData(double scale) {
//********************************************************************
fDataHist->Scale(scale);
}
//********************************************************************
void Measurement1D::ScaleDataErrors(double scale) {
//********************************************************************
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
fDataHist->SetBinError(i + 1, fDataHist->GetBinError(i + 1) * scale);
}
}
//********************************************************************
void Measurement1D::ScaleCovar(double scale) {
//********************************************************************
(*fFullCovar) *= scale;
(*covar) *= 1.0 / scale;
(*fDecomp) *= sqrt(scale);
}
//********************************************************************
void Measurement1D::SetBinMask(std::string maskfile) {
//********************************************************************
if (!fIsMask) return;
LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl;
// Create a mask histogram with dim of data
int nbins = fDataHist->GetNbinsX();
fMaskHist =
new TH1I((fSettings.GetName() + "_BINMASK").c_str(),
(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbins, 0, nbins);
std::string line;
std::ifstream mask(maskfile.c_str(), ifstream::in);
if (!mask.is_open()) {
LOG(FTL) << " Cannot find mask file." << std::endl;
throw;
}
while (std::getline(mask >> std::ws, line, '\n')) {
std::vector<int> entries = GeneralUtils::ParseToInt(line, " ");
// Skip lines with poorly formatted lines
if (entries.size() < 2) {
LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << line
<< std::endl;
continue;
}
// The first index should be the bin number, the second should be the mask
// value.
int val = 0;
if (entries[1] > 0) val = 1;
fMaskHist->SetBinContent(entries[0], val);
}
// Apply masking by setting masked data bins to zero
PlotUtils::MaskBins(fDataHist, fMaskHist);
return;
}
//********************************************************************
void Measurement1D::FinaliseMeasurement() {
//********************************************************************
LOG(SAM) << "Finalising Measurement: " << fName << std::endl;
// Make sure data is setup
if (!fDataHist) {
ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl;
throw;
}
// Make sure covariances are setup
if (!fFullCovar) {
SetCovarFromDiagonal(fDataHist);
}
if (!covar) {
covar = StatUtils::GetInvert(fFullCovar);
}
if (!fDecomp) {
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
// Setup fMCHist from data
fMCHist = (TH1D*)fDataHist->Clone();
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist->Reset();
// Setup fMCFine
- fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX(),
+ fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX() * 6,
fMCHist->GetBinLowEdge(1),
fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1));
fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCFine->Reset();
// Setup MC Stat
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
// Search drawopts for possible types to include by default
std::string drawopts = FitPar::Config().GetParS("drawopts");
if (drawopts.find("MODES") != std::string::npos) {
fMCHist_Modes = new TrueModeStack( (fSettings.GetName() + "_MODES").c_str(),
("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);
}
// Setup bin masks using sample name
if (fIsMask) {
std::string curname = fName;
std::string origname = fSettings.GetS("originalname");
// Check rename.mask
std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");
// Check origname.mask
if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask");
// Check database
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
}
// Setup Bin Mask
SetBinMask(maskloc);
}
if (fScaleFactor < 0) {
ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl;
ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl;
ERR(FTL) << "EXITING" << std::endl;
throw;
}
// Create and fill Weighted Histogram
if (!fMCWeighted) {
fMCWeighted = (TH1D*)fMCHist->Clone();
fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
(fName + "_MCWGHTS" + fPlotTitles).c_str());
fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
}
}
//********************************************************************
void Measurement1D::SetFitOptions(std::string opt) {
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT") return;
// CHECK Conflicting Fit Options
std::vector<std::string> fit_option_allow =
GeneralUtils::ParseToStr(fAllowedTypes, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
std::vector<std::string> fit_option_section =
GeneralUtils::ParseToStr(fit_option_allow.at(i), ",");
bool found_option = false;
for (UInt_t j = 0; j < fit_option_section.size(); j++) {
std::string av_opt = fit_option_section.at(j);
if (!found_option and opt.find(av_opt) != std::string::npos) {
found_option = true;
} else if (found_option and opt.find(av_opt) != std::string::npos) {
ERR(FTL) << "ERROR: Conflicting fit options provided: "
<< opt << std::endl
<< "Conflicting group = " << fit_option_section.at(i) << std::endl
<< "You should only supply one of these options in card file." << std::endl;
throw;
}
}
}
// Check all options are allowed
std::vector<std::string> fit_options_input =
GeneralUtils::ParseToStr(opt, "/");
for (UInt_t i = 0; i < fit_options_input.size(); i++) {
if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i)
<< "' Provided is not allowed for this measurement."
<< std::endl;
ERR(FTL) << "Fit Options should be provided as a '/' seperated list "
"(e.g. FREE/DIAG/NORM)"
<< std::endl;
ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes
<< "'" << std::endl;
throw;
}
}
// Set TYPE
fFitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos) {
fIsFree = fIsShape = false;
fIsFix = true;
} else if (opt.find("SHAPE") != std::string::npos) {
fIsFree = fIsFix = false;
fIsShape = true;
} else if (opt.find("FREE") != std::string::npos) {
fIsFix = fIsShape = false;
fIsFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos) {
fIsDiag = true;
fIsFull = false;
} else if (opt.find("FULL") != std::string::npos) {
fIsDiag = false;
fIsFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos) {
fIsChi2 = false;
ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl;
ERR(FTL) << "Try to use a chi2!" << std::endl;
throw;
} else {
fIsChi2 = true;
}
// EXTRAS
if (opt.find("RAW") != std::string::npos) fIsRawEvents = true;
if (opt.find("NOWIDTH") != std::string::npos) fIsNoWidth = true;
if (opt.find("DIF") != std::string::npos) fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true;
if (opt.find("NORM") != std::string::npos) fAddNormPen = true;
if (opt.find("MASK") != std::string::npos) fIsMask = true;
return;
};
//********************************************************************
void Measurement1D::SetSmearingMatrix(std::string smearfile, int truedim,
int recodim) {
//********************************************************************
// The smearing matrix describes the migration from true bins (rows) to reco
// bins (columns)
// Counter over the true bins!
int row = 0;
std::string line;
std::ifstream smear(smearfile.c_str(), ifstream::in);
// Note that the smearing matrix may be rectangular.
fSmearMatrix = new TMatrixD(truedim, recodim);
if (smear.is_open())
LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl;
else
ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile
<< std::endl;
while (std::getline(smear >> std::ws, line, '\n')) {
int column = 0;
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*fSmearMatrix)(row, column) =
(*iter) / 100.; // Convert to fraction from
// percentage (this may not be
// general enough)
column++;
}
row++;
}
return;
}
//********************************************************************
void Measurement1D::ApplySmearingMatrix() {
//********************************************************************
if (!fSmearMatrix) {
ERR(WRN) << fName
<< ": attempted to apply smearing matrix, but none was set"
<< std::endl;
return;
}
TH1D* unsmeared = (TH1D*)fMCHist->Clone();
TH1D* smeared = (TH1D*)fMCHist->Clone();
smeared->Reset();
// Loop over reconstructed bins
// true = row; reco = column
for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) {
// Sum up the constributions from all true bins
double rBinVal = 0;
// Loop over true bins
for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) {
rBinVal +=
(*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1);
}
smeared->SetBinContent(rbin + 1, rBinVal);
}
fMCHist = (TH1D*)smeared->Clone();
return;
}
/*
Reconfigure LOOP
*/
//********************************************************************
void Measurement1D::ResetAll() {
//********************************************************************
fMCHist->Reset();
fMCFine->Reset();
fMCStat->Reset();
return;
};
//********************************************************************
void Measurement1D::FillHistograms() {
//********************************************************************
if (Signal) {
fMCHist->Fill(fXVar, Weight);
fMCFine->Fill(fXVar, Weight);
fMCStat->Fill(fXVar, 1.0);
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight);
}
return;
};
//********************************************************************
void Measurement1D::ScaleEvents() {
//********************************************************************
// Fill MCWeighted;
// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
// fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
// fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
// }
// Setup Stat ratios for MC and MC Fine
double* statratio = new double[fMCHist->GetNbinsX()];
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
if (fMCHist->GetBinContent(i + 1) != 0) {
statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
} else {
statratio[i] = 0.0;
}
}
double* statratiofine = new double[fMCFine->GetNbinsX()];
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
if (fMCFine->GetBinContent(i + 1) != 0) {
statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
} else {
statratiofine[i] = 0.0;
}
}
// Scaling for raw event rates
if (fIsRawEvents) {
double datamcratio = fDataHist->Integral() / fMCHist->Integral();
fMCHist->Scale(datamcratio);
fMCFine->Scale(datamcratio);
if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio);
// Scaling for XSec as function of Enu
} else if (fIsEnu1D) {
PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
// if (fMCHist_Modes) {
// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
// GetEventHistogram(), fScaleFactor,
// fNEvents);
// }
} else if (fIsNoWidth) {
fMCHist->Scale(fScaleFactor);
fMCFine->Scale(fScaleFactor);
if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor);
// Any other differential scaling
} else {
fMCHist->Scale(fScaleFactor, "width");
fMCFine->Scale(fScaleFactor, "width");
if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width");
}
// Proper error scaling - ROOT Freaks out with xsec weights sometimes
for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
}
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]);
}
// Clean up
delete statratio;
delete statratiofine;
return;
};
//********************************************************************
void Measurement1D::ApplyNormScale(double norm) {
//********************************************************************
fCurrentNorm = norm;
fMCHist->Scale(1.0 / norm);
fMCFine->Scale(1.0 / norm);
return;
};
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int Measurement1D::GetNDOF() {
//********************************************************************
int ndof = fDataHist->GetNbinsX();
if (fMaskHist) ndof -= fMaskHist->Integral();
return ndof;
}
//********************************************************************
double Measurement1D::GetLikelihood() {
//********************************************************************
// If this is for a ratio, there is no data histogram to compare to!
if (fNoData || !fDataHist) return 0.;
// Apply Masking to MC if Required.
if (fIsMask and fMaskHist) {
PlotUtils::MaskBins(fMCHist, fMaskHist);
}
// Sort Shape Scaling
double scaleF = 0.0;
// TODO Include !fIsRawEvents
if (fIsShape) {
if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) {
scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") /
fMCHist->Integral(1, fMCHist->GetNbinsX(), "width");
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
}
}
// Likelihood Calculation
double stat = 0.;
if (fIsChi2) {
if (fIsRawEvents) {
stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist);
} else if (fIsDiag) {
stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist);
} else if (!fIsDiag and !fIsRawEvents) {
std::cout << "Getting likelihood from covariance " << std::endl;
stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist);
}
}
// Sort Penalty Terms
if (fAddNormPen) {
double penalty =
(1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError);
stat += penalty;
}
// Return to normal scaling
if (fIsShape) { // and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
fLikelihood = stat;
return stat;
}
/*
Fake Data Functions
*/
//********************************************************************
void Measurement1D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Setup original/datatrue
TH1D* tempdata = (TH1D*) fDataHist->Clone();
if (!fIsFakeData) {
fIsFakeData = true;
// Make a copy of the original data histogram.
if (!fDataOrig) fDataOrig = (TH1D*)fDataHist->Clone((fName + "_data_original").c_str());
} else {
ResetFakeData();
}
// Setup Inputs
fFakeDataInput = fakeOption;
LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl;
// From MC
if (fFakeDataInput.compare("MC") == 0) {
fDataHist = (TH1D*)fMCHist->Clone((fName + "_MC").c_str());
// Fake File
} else {
if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
fDataHist = (TH1D*)fFakeDataFile->Get((fName + "_MC").c_str());
}
// Setup Data Hist
fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
(fName + fPlotTitles).c_str());
// Replace Data True
if (fDataTrue) delete fDataTrue;
fDataTrue = (TH1D*)fDataHist->Clone();
fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
(fName + fPlotTitles).c_str());
// Make a new covariance for fake data hist.
int nbins = fDataHist->GetNbinsX();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
(*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j);
}
}
// Setup Covariances
if (covar) delete covar;
covar = StatUtils::GetInvert(fFullCovar);
if (fDecomp) delete fDecomp;
fDecomp = StatUtils::GetInvert(fFullCovar);
delete tempdata;
return;
};
//********************************************************************
void Measurement1D::ResetFakeData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
}
}
//********************************************************************
void Measurement1D::ResetData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
}
fIsFakeData = false;
}
//********************************************************************
void Measurement1D::ThrowCovariance() {
//********************************************************************
// Take a fDecomposition and use it to throw the current dataset.
// Requires fDataTrue also be set incase used repeatedly.
if (fDataHist) delete fDataHist;
fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
return;
};
/*
Access Functions
*/
//********************************************************************
TH1D* Measurement1D::GetMCHistogram() {
//********************************************************************
if (!fMCHist) return fMCHist;
std::ostringstream chi2;
chi2 << std::setprecision(5) << this->GetLikelihood();
int linecolor = kRed;
int linestyle = 1;
int linewidth = 1;
int fillcolor = 0;
int fillstyle = 1001;
// if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor");
// if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle");
// if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth");
// if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor");
// if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle");
fMCHist->SetTitle(chi2.str().c_str());
fMCHist->SetLineColor(linecolor);
fMCHist->SetLineStyle(linestyle);
fMCHist->SetLineWidth(linewidth);
fMCHist->SetFillColor(fillcolor);
fMCHist->SetFillStyle(fillstyle);
return fMCHist;
};
//********************************************************************
TH1D* Measurement1D::GetDataHistogram() {
//********************************************************************
if (!fDataHist) return fDataHist;
int datacolor = kBlack;
int datastyle = 1;
int datawidth = 1;
// if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor");
// if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle");
// if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth");
fDataHist->SetLineColor(datacolor);
fDataHist->SetLineWidth(datawidth);
fDataHist->SetMarkerStyle(datastyle);
return fDataHist;
};
/*
Write Functions
*/
// Save all the histograms at once
//********************************************************************
void Measurement1D::Write(std::string drawOpt) {
//********************************************************************
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
// Write Settigns
if (drawOpt.find("SETTINGS") != std::string::npos){
fSettings.Set("#chi^{2}",fLikelihood);
fSettings.Set("NDOF", this->GetNDOF() );
fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF() );
fSettings.Write();
}
// Write Data/MC
GetDataList().at(0)->Write();
GetMCList().at(0)->Write();
// Write Fine Histogram
if (drawOpt.find("FINE") != std::string::npos)
GetFineList().at(0)->Write();
// Write Weighted Histogram
if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted)
fMCWeighted->Write();
// Save Flux/Evt if no event manager
if (!FitPar::Config().GetParB("EventManager")) {
if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram())
GetFluxHistogram()->Write();
if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
}
// Write Mask
if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) {
fMaskHist->Write();
}
// Write Covariances
if (drawOpt.find("COV") != std::string::npos && fFullCovar) {
PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName());
}
if (drawOpt.find("INVCOV") != std::string::npos && covar) {
PlotUtils::GetInvCovarPlot(covar, fSettings.GetName());
}
if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) {
PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName());
}
// // Likelihood residual plots
// if (drawOpt.find("RESIDUAL") != std::string::npos) {
// WriteResidualPlots();
// }
// Ratio and Shape Plots
if (drawOpt.find("RATIO") != std::string::npos) {
WriteRatioPlot();
}
if (drawOpt.find("SHAPE") != std::string::npos) {
WriteShapePlot();
if (drawOpt.find("RATIO") != std::string::npos)
WriteShapeRatioPlot();
}
// // RATIO
// if (drawOpt.find("CANVMC") != std::string::npos) {
// TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist);
// c1->Write();
// delete c1;
// }
// // PDG
// if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) {
// TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes);
// c2->Write();
// delete c2;
// }
// Write Extra Histograms
AutoWriteExtraTH1();
WriteExtraHistograms();
// Returning
LOG(SAM) << "Written Histograms: " << fName << std::endl;
return;
}
//********************************************************************
void Measurement1D::WriteRatioPlot() {
//********************************************************************
// Setup mc data ratios
TH1D* dataRatio = (TH1D*)fDataHist->Clone((fName + "_data_RATIO").c_str());
TH1D* mcRatio = (TH1D*)fMCHist->Clone((fName + "_MC_RATIO").c_str());
// Extra MC Data Ratios
for (int i = 0; i < mcRatio->GetNbinsX(); i++) {
dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1));
dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1));
mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
}
// Write ratios
mcRatio->Write();
dataRatio->Write();
delete mcRatio;
delete dataRatio;
}
//********************************************************************
void Measurement1D::WriteShapePlot() {
//********************************************************************
TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
std::stringstream ss;
ss << shapeScale;
mcShape->SetTitle(ss.str().c_str());
mcShape->SetLineWidth(3);
mcShape->SetLineStyle(7);
mcShape->Write();
delete mcShape;
}
//********************************************************************
void Measurement1D::WriteShapeRatioPlot() {
//********************************************************************
// Get a mcshape histogram
TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
// Create shape ratio histograms
TH1D* mcShapeRatio = (TH1D*)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
TH1D* dataShapeRatio = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());
// Divide the histograms
mcShapeRatio->Divide(mcShape);
dataShapeRatio->Divide(mcShape);
// Colour the shape ratio plots
mcShapeRatio->SetLineWidth(3);
mcShapeRatio->SetLineStyle(7);
mcShapeRatio->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
//// CRAP TO BE REMOVED
//********************************************************************
void Measurement1D::SetupMeasurement(std::string inputfile, std::string type,
FitWeight * rw, std::string fkdt) {
//********************************************************************
//nuiskey samplekey = Config::CreateKey("sample");
// samplekey.AddS("name", fName);
// samplekey.AddS("type",type);
// samplekey.AddS("input",inputfile);
// fSettings = LoadSampleSettings(samplekey);
// Reset everything to NULL
// Init();
// Check if name contains Evt, indicating that it is a raw number of events
// measurements and should thus be treated as once
fIsRawEvents = false;
if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!"
<< std::endl;
LOG(SAM) << "Overriding this and setting fIsRawEvents == true!"
<< std::endl;
}
fIsEnu1D = false;
if (fName.find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
"not flux averaged!"
<< std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
fRW = rw;
if (!fInput and !fIsJoint) SetupInputs(inputfile);
// Set Default Options
SetFitOptions(fDefaultTypes);
// Set Passed Options
SetFitOptions(type);
// Still adding support for flat flux inputs
// // Set Enu Flux Scaling
// if (isFlatFluxFolding) this->Input()->ApplyFluxFolding(
// this->defaultFluxHist );
// FinaliseMeasurement();
}
//********************************************************************
void Measurement1D::SetupDefaultHist() {
//********************************************************************
// Setup fMCHist
fMCHist = (TH1D*)fDataHist->Clone();
fMCHist->SetNameTitle((fName + "_MC").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
// Setup fMCFine
Int_t nBins = fMCHist->GetNbinsX();
fMCFine = new TH1D(
(fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(),
nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1));
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
fMCHist->Reset();
fMCFine->Reset();
// Setup the NEUT Mode Array
PlotUtils::CreateNeutModeArray((TH1D*)fMCHist, (TH1**)fMCHist_PDG);
PlotUtils::ResetNeutModeArray((TH1**)fMCHist_PDG);
// Setup bin masks using sample name
if (fIsMask) {
std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
}
SetBinMask(maskloc);
}
fMCHist_Modes = new TrueModeStack( (fName + "_MODES").c_str(), ("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);
return;
}
//********************************************************************
void Measurement1D::SetDataValues(std::string dataFile) {
//********************************************************************
// Override this function if the input file isn't in a suitable format
LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
fDataHist =
PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles);
fDataTrue = (TH1D*)fDataHist->Clone();
// Number of data points is number of bins
fNDataPointsX = fDataHist->GetXaxis()->GetNbins();
return;
};
//********************************************************************
void Measurement1D::SetDataFromDatabase(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(
(GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void Measurement1D::SetDataFromFile(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void Measurement1D::SetCovarMatrix(std::string covarFile) {
//********************************************************************
// Covariance function, only really used when reading in the MB Covariances.
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
TH2D* covarPlot = new TH2D();
// TH2D* decmpPlot = new TH2D();
TH2D* covarInvPlot = new TH2D();
TH2D* fFullCovarPlot = new TH2D();
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("thrown_covariance");
if (fIsShape || fIsFree) covName = "shp_";
if (fIsDiag)
covName += "diag";
else
covName += "full";
covarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").c_str());
if (!covOption.compare("SUB"))
fFullCovarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
else if (!covOption.compare("FULL"))
fFullCovarPlot = (TH2D*)tempFile->Get("fullcov");
else
ERR(WRN) << "Incorrect thrown_covariance option in parameters."
<< std::endl;
int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
int covdim = int(fDataHist->GetNbinsX());
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
fDecomp = new TMatrixDSym(dim);
int row, column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++) {
// if (this->masked->GetBinContent(i+1) > 0) continue;
for (Int_t j = 0; j < covdim; j++) {
// if (this->masked->GetBinContent(j+1) > 0) continue;
(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);
column++;
}
column = 0;
row++;
}
// Set bin errors on data
if (!fIsDiag) {
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar);
}
// Get Deteriminant and inverse matrix
// fCovDet = this->covar->Determinant();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// Sets the covariance matrix from a provided file in a text format
// scale is a multiplicative pre-factor to apply in the case where the
// covariance is given in some unit (e.g. 1E-38)
void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim,
double scale) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream covarread(covarFile.c_str(), ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covarread.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
else
ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile
<< std::endl;
// Loop over the lines in the file
while (std::getline(covarread >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
ERR(WRN) << "SetCovarMatrixFromText -> Covariance matrix only has <= 1 "
"entries on this line: "
<< row << std::endl;
}
for (std::vector<double>::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*covar)(row, column) = *iter;
(*fFullCovar)(row, column) = *iter;
column++;
}
row++;
}
covarread.close();
// Scale the actualy covariance matrix by some multiplicative factor
(*fFullCovar) *= scale;
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
// THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
// Now need to multiply by the scaling factor
// If the covariance
(*this->covar) *= 1. / (scale);
return;
};
//********************************************************************
void Measurement1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream corr(corrFile.c_str(), ifstream::in);
this->covar = new TMatrixDSym(dim);
this->fFullCovar = new TMatrixDSym(dim);
if (corr.is_open())
LOG(SAM) << "Reading and converting correlation matrix from file: "
<< corrFile << std::endl;
else {
ERR(FTL) << "Correlation matrix provided is incorrect: " << corrFile
<< std::endl;
exit(-1);
}
while (std::getline(corr >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::iterator iter = entries.begin();
iter != entries.end(); iter++) {
double val = (*iter) * this->fDataHist->GetBinError(row + 1) * 1E38 *
this->fDataHist->GetBinError(column + 1) * 1E38;
if (val == 0) {
ERR(FTL) << "Found a zero value in the covariance matrix, assuming "
"this is an error!"
<< std::endl;
exit(-1);
}
(*this->covar)(row, column) = val;
(*this->fFullCovar)(row, column) = val;
column++;
}
row++;
}
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// FullUnits refers to if we have "real" unscaled units in the covariance matrix, e.g. 1E-76.
// If this is the case we need to scale it so that the chi2 contribution is correct
// NUISANCE internally assumes the covariance matrix has units of 1E76
void Measurement1D::SetCovarFromDataFile(std::string covarFile,
std::string covName, bool FullUnits) {
//********************************************************************
LOG(SAM) << "Getting covariance from " << covarFile << "->" << covName
<< std::endl;
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
TH2D* covPlot = (TH2D*)tempFile->Get(covName.c_str());
covPlot->SetDirectory(0);
// Scale the covariance matrix if it comes in normal units
if (FullUnits) {
covPlot->Scale(1.E76);
}
int dim = covPlot->GetNbinsX();
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < dim; i++) {
for (int j = 0; j < dim; j++) {
(*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1);
}
}
this->covar = (TMatrixDSym*)fFullCovar->Clone();
fDecomp = (TMatrixDSym*)fFullCovar->Clone();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
TDecompChol LUChol = TDecompChol(*fDecomp);
LUChol.Decompose();
fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), "");
return;
};
// //********************************************************************
// void Measurement1D::SetBinMask(std::string maskFile) {
// //********************************************************************
// // Create a mask histogram.
// int nbins = fDataHist->GetNbinsX();
// fMaskHist =
// new TH1I((fName + "_fMaskHist").c_str(),
// (fName + "_fMaskHist; Bin; Mask?").c_str(), nbins, 0, nbins);
// std::string line;
// std::ifstream mask(maskFile.c_str(), ifstream::in);
// if (mask.is_open())
// LOG(SAM) << "Reading bin mask from file: " << maskFile << std::endl;
// else
// LOG(FTL) << " Cannot find mask file." << std::endl;
// while (std::getline(mask >> std::ws, line, '\n')) {
// std::vector<int> entries = GeneralUtils::ParseToInt(line, " ");
// // Skip lines with poorly formatted lines
// if (entries.size() < 2) {
// LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << line
// << std::endl;
// continue;
// }
// // The first index should be the bin number, the second should be the mask
// // value.
// fMaskHist->SetBinContent(entries[0], entries[1]);
// }
// // Set masked data bins to zero
// PlotUtils::MaskBins(fDataHist, fMaskHist);
// return;
// }
// //********************************************************************
// void Measurement1D::GetBinContents(std::vector<double>& cont,
// std::vector<double>& err) {
// //********************************************************************
// // Return a vector of the main bin contents
// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
// cont.push_back(fMCHist->GetBinContent(i + 1));
// err.push_back(fMCHist->GetBinError(i + 1));
// }
// return;
// };
/*
XSec Functions
*/
// //********************************************************************
// void Measurement1D::SetFluxHistogram(std::string fluxFile, int minE, int
// maxE,
// double fluxNorm) {
// //********************************************************************
// // Note this expects the flux bins to be given in terms of MeV
// LOG(SAM) << "Reading flux from file: " << fluxFile << std::endl;
// TGraph f(fluxFile.c_str(), "%lg %lg");
// fFluxHist =
// new TH1D((fName + "_flux").c_str(), (fName + "; E_{#nu} (GeV)").c_str(),
// f.GetN() - 1, minE, maxE);
// Double_t* yVal = f.GetY();
// for (int i = 0; i < fFluxHist->GetNbinsX(); ++i)
// fFluxHist->SetBinContent(i + 1, yVal[i] * fluxNorm);
// };
// //********************************************************************
// double Measurement1D::TotalIntegratedFlux(std::string intOpt, double low,
// double high) {
// //********************************************************************
// if (fInput->GetType() == kGiBUU) {
// return 1.0;
// }
// // The default case of low = -9999.9 and high = -9999.9
// if (low == -9999.9) low = this->EnuMin;
// if (high == -9999.9) high = this->EnuMax;
// int minBin = fFluxHist->GetXaxis()->FindBin(low);
// int maxBin = fFluxHist->GetXaxis()->FindBin(high);
// // Get integral over custom range
// double integral = fFluxHist->Integral(minBin, maxBin + 1, intOpt.c_str());
// return integral;
// };
diff --git a/src/FitBase/ParamPull.cxx b/src/FitBase/ParamPull.cxx
index f443f82..4143afe 100644
--- a/src/FitBase/ParamPull.cxx
+++ b/src/FitBase/ParamPull.cxx
@@ -1,792 +1,810 @@
// 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 "ParamPull.h"
//*******************************************************************************
-ParamPull::ParamPull(std::string name, std::string inputfile, std::string type, std::string dials){
+ParamPull::ParamPull(std::string name, std::string inputfile, std::string type, std::string dials) {
//*******************************************************************************
fMinHist = NULL;
fMaxHist = NULL;
fTypeHist = NULL;
fDialSelection = dials;
-
+
fName = name;
fInput = inputfile;
fType = type;
-
+
// Set the pull type
SetType(fType);
std::cout << fType << std::endl;
-
+
// Setup Histograms from input file
- SetupHistograms(fInput);
+ SetupHistograms(fInput);
};
//*******************************************************************************
-void ParamPull::SetType(std::string type){
+void ParamPull::SetType(std::string type) {
//*******************************************************************************
fType = type;
// Assume Default if empty
- if (type.empty() || type == "DEFAULT"){
+ if (type.empty() || type == "DEFAULT") {
ERR(WRN) << "No type specified for ParmPull class " << fName << std::endl;
ERR(WRN) << "Assuming GAUSTHROW/GAUSPULL" << std::endl;
type = "GAUSTHROW/GAUSPULL";
}
-
+
// Set Dial options
- if (type.find("FRAC") != std::string::npos){
+ if (type.find("FRAC") != std::string::npos) {
fDialOptions = "FRAC";
fPlotTitles = ";; Fractional RW Value";
- } else if (type.find("ABS") != std::string::npos){
+ } else if (type.find("ABS") != std::string::npos) {
fDialOptions = "ABS";
- fPlotTitles= ";; ABS RW Value";
+ fPlotTitles = ";; ABS RW Value";
} else {
fDialOptions = "";
- fPlotTitles= ";; RW Value";
+ fPlotTitles = ";; RW Value";
}
// Parse throw types
if (type.find("GAUSPULL") != std::string::npos) fCalcType = kGausPull;
else fCalcType = kNoPull;
if (type.find("GAUSTHROW") != std::string::npos) fThrowType = kGausThrow;
else fThrowType = kNoThrow;
// Extra check to see if throws or pulls are turned off
if (type.find("NOPULL") != std::string::npos) fCalcType = kNoPull;
if (type.find("NOTHROW") != std::string::npos) fThrowType = kNoThrow;
-
+
}
//*******************************************************************************
-void ParamPull::SetupHistograms(std::string input){
+void ParamPull::SetupHistograms(std::string input) {
//*******************************************************************************
// Extract Types from Input
fFileType = "";
const int nfiletypes = 4;
- const std::string filetypes[nfiletypes] = {"FIT","ROOT","TXT","DIAL"};
+ const std::string filetypes[nfiletypes] = {"FIT", "ROOT", "TXT", "DIAL"};
for (int i = 0; i < nfiletypes; i++) {
std::string tempTypes = filetypes[i] + ":";
if (input.find(tempTypes) != std::string::npos) {
fFileType = filetypes[i];
input.replace(input.find(tempTypes), tempTypes.size(), "");
break;
}
}
// Read Files
if (!fFileType.compare("FIT")) ReadFitFile(input);
else if (!fFileType.compare("ROOT")) ReadRootFile(input);
else if (!fFileType.compare("VECT")) ReadVectFile(input);
else if (!fFileType.compare("DIAL")) ReadDialInput(input);
else {
ERR(FTL) << "Unknown ParamPull Type: " << input << std::endl;
throw;
}
-
+
// Check Dials are all good
CheckDialsValid();
// Setup MC Histogram
fMCHist = (TH1D*) fDataHist->Clone();
fMCHist->Reset();
fMCHist->SetNameTitle( (fName + "_MC").c_str(),
- (fName + " MC" + fPlotTitles).c_str() );
+ (fName + " MC" + fPlotTitles).c_str() );
// If no Covar input make an uncorrelated one
- if (!fCovar){
+ if (!fCovar) {
fCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist, 1.0);
}
// If no types or limits are provided give them a default option
- if (!fMinHist){
+ if (!fMinHist) {
fMinHist = (TH1D*) fDataHist->Clone();
fMinHist->SetNameTitle( (fName + "_min").c_str(),
- (fName + " min" + fPlotTitles).c_str() );
- for (int i = 0; i < fMinHist->GetNbinsX(); i++){
- // TODO (P.Stowell) Change this to a NULL system where limits are actually free!
- fMinHist->SetBinContent(i+1, fDataHist->GetBinContent(i+1) - 1E6);
+ (fName + " min" + fPlotTitles).c_str() );
+ for (int i = 0; i < fMinHist->GetNbinsX(); i++) {
+ // TODO (P.Stowell) Change this to a NULL system where limits are actually free!
+ fMinHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
- if (!fMaxHist){
+ if (!fMaxHist) {
fMaxHist = (TH1D*) fDataHist->Clone();
fMaxHist->SetNameTitle( (fName + "_min").c_str(),
- (fName + " min" + fPlotTitles).c_str() );
- for(int i = 0; i < fMaxHist->GetNbinsX(); i++){
- fMaxHist->SetBinContent(i+1, fDataHist->GetBinContent(i+1) - 1E6);
+ (fName + " min" + fPlotTitles).c_str() );
+ for (int i = 0; i < fMaxHist->GetNbinsX(); i++) {
+ fMaxHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
// Set types from state, or to unknown
- if (!fTypeHist){
+ if (!fTypeHist) {
int deftype = -1;
- if (fType.find("T2K") != std::string::npos){ deftype = kT2K; }
- else if (fType.find("NEUT") != std::string::npos){ deftype = kNEUT; }
- else if (fType.find("NIWG") != std::string::npos){ deftype = kNIWG; }
- else if (fType.find("GENIE") != std::string::npos){ deftype = kGENIE; }
- else if (fType.find("NORM") != std::string::npos){ deftype = kNORM; }
- else if (fType.find("NUWRO") != std::string::npos){ deftype = kNUWRO; }
-
+ if (fType.find("T2K") != std::string::npos) { deftype = kT2K; }
+ else if (fType.find("NEUT") != std::string::npos) { deftype = kNEUT; }
+ else if (fType.find("NIWG") != std::string::npos) { deftype = kNIWG; }
+ else if (fType.find("GENIE") != std::string::npos) { deftype = kGENIE; }
+ else if (fType.find("NORM") != std::string::npos) { deftype = kNORM; }
+ else if (fType.find("NUWRO") != std::string::npos) { deftype = kNUWRO; }
+
fTypeHist = new TH1I( (fName + "_type").c_str(),
- (fName + " type" + fPlotTitles).c_str(),
- fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX() );
-
- for(int i = 0; i < fTypeHist->GetNbinsX(); i++){
- fTypeHist->SetBinContent(i+1, deftype );
+ (fName + " type" + fPlotTitles).c_str(),
+ fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX() );
+
+ for (int i = 0; i < fTypeHist->GetNbinsX(); i++) {
+ fTypeHist->SetBinContent(i + 1, deftype );
}
}
-
-
+
+
// Sort Covariances
fInvCovar = StatUtils::GetInvert(fCovar);
fDecomp = StatUtils::GetDecomp(fCovar);
// Create DataTrue for Throws
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
- (fName + " truedata" + fPlotTitles).c_str() );
+ (fName + " truedata" + fPlotTitles).c_str() );
fDataOrig = (TH1D*) fDataHist->Clone();
fDataOrig->SetNameTitle( (fName + "_origdata").c_str(),
- (fName + " origdata" + fPlotTitles).c_str() );
+ (fName + " origdata" + fPlotTitles).c_str() );
// Select only dials we want
- if (!fDialSelection.empty()){
+ if (!fDialSelection.empty()) {
(*fDataHist) = RemoveBinsNotInString(*fDataHist, fDialSelection);
-
-
+
+
}
-
+
}
-//*******************************************************************************
-TH1D ParamPull::RemoveBinsNotInString(TH1D hist, std::string mystr){
-//*******************************************************************************
+//*******************************************************************************
+TH1D ParamPull::RemoveBinsNotInString(TH1D hist, std::string mystr) {
+//*******************************************************************************
// Make list of allowed bins
std::vector<std::string> allowedbins;
- for (int i = 0; i < hist.GetNbinsX(); i++){
- std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i+1));
-
- if (mystr.find(syst) != std::string::npos){
+ for (int i = 0; i < hist.GetNbinsX(); i++) {
+ std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
+
+ if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1D newhist = TH1D( hist.GetName(), hist.GetTitle(),
- (Int_t)nbins, 0.0, (Double_t)nbins);
+ (Int_t)nbins, 0.0, (Double_t)nbins);
// Setup bins
- for (UInt_t i = 0; i < nbins; i++){
+ for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
- newhist.GetXaxis()->SetBinLabel(i+1, allowedbins[i].c_str());
+ newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
- for (Int_t j = 0; j < hist.GetNbinsX(); j++){
- if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j+1))){
- newhist.SetBinContent(i+1, hist.GetBinContent(j+1) );
- newhist.SetBinError(i+1, hist.GetBinError(j+1) );
+ for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
+ if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
+ newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
+ newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
//*******************************************************************************
-TH1I ParamPull::RemoveBinsNotInString(TH1I hist, std::string mystr){
+TH1I ParamPull::RemoveBinsNotInString(TH1I hist, std::string mystr) {
//*******************************************************************************
// Make list of allowed bins
std::vector<std::string> allowedbins;
- for (int i = 0; i < hist.GetNbinsX(); i++){
- std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i+1));
+ for (int i = 0; i < hist.GetNbinsX(); i++) {
+ std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
- if (mystr.find(syst) != std::string::npos){
+ if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1I newhist = TH1I( hist.GetName(), hist.GetTitle(),
- (Int_t)nbins, 0.0, (Int_t)nbins);
+ (Int_t)nbins, 0.0, (Int_t)nbins);
// Setup bins
- for (UInt_t i = 0; i < nbins; i++){
+ for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
- newhist.GetXaxis()->SetBinLabel(i+1, allowedbins[i].c_str());
+ newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
- for (Int_t j = 0; j < hist.GetNbinsX(); j++){
- if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j+1))){
- newhist.SetBinContent(i+1, hist.GetBinContent(j+1) );
- newhist.SetBinError(i+1, hist.GetBinError(j+1) );
+ for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
+ if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
+ newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
+ newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
-//*******************************************************************************
-void ParamPull::ReadFitFile(std::string input){
+//*******************************************************************************
+void ParamPull::ReadFitFile(std::string input) {
//*******************************************************************************
- TFile* tempfile = new TFile(input.c_str(),"READ");
+ TFile* tempfile = new TFile(input.c_str(), "READ");
// Read Data
fDataHist = (TH1D*) tempfile->Get("fit_dials_free");
- if (!fDataHist){
+ if (!fDataHist) {
ERR(FTL) << "Can't find TH1D hist fit_dials in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
- (fName + " data" + fPlotTitles).c_str() );
+ (fName + " data" + fPlotTitles).c_str() );
// Read Covar
TH2D* tempcov = (TH2D*) tempfile->Get("covariance_free");
- if (!tempcov){
+ if (!tempcov) {
ERR(FTL) << "Can't find TH2D covariance_free in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
- for (int i = 0; i < nbins; i++){
- for(int j = 0; j < nbins; j++){
- (*fCovar)(i,j) = tempcov->GetBinContent(i+1,j+1);
+ for (int i = 0; i < nbins; i++) {
+ for (int j = 0; j < nbins; j++) {
+ (*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
-
+
return;
}
//*******************************************************************************
-void ParamPull::ReadRootFile(std::string input){
-//*******************************************************************************
+void ParamPull::ReadRootFile(std::string input) {
+//*******************************************************************************
- std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input,";");
+ std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input, ";");
// Check all given
- if (inputlist.size() < 2){
+ if (inputlist.size() < 2) {
ERR(FTL) << "Covar supplied in 'ROOT' format should have 3 semi-colon seperated entries!" << std::endl
- << "ROOT:filename;histname[;covarname]" << std::endl;
+ << "ROOT:filename;histname[;covarname]" << std::endl;
ERR(FTL) << "histname = TH1D, covarname = TH2D" << std::endl;
throw;
}
// Get Entries
std::string filename = inputlist[0];
LOG(DEB) << filename << std::endl;
std::string histname = inputlist[1];
LOG(DEB) << histname << std::endl;
- LOG(DEB) << input<< std::endl;
+ LOG(DEB) << input << std::endl;
// Read File
- TFile* tempfile = new TFile(filename.c_str(),"READ");
- if (tempfile->IsZombie()){
+ TFile* tempfile = new TFile(filename.c_str(), "READ");
+ if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << filename << std::endl;
throw;
}
// Read Hist
fDataHist = (TH1D*) tempfile->Get(histname.c_str());
- if (!fDataHist){
+ if (!fDataHist) {
ERR(FTL) << "Can't find TH1D hist " << histname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
- (fName + " data" + fPlotTitles).c_str() );
+ (fName + " data" + fPlotTitles).c_str() );
LOG(DEB) << "READING COVAR" << std::endl;
// Read Covar
- if (inputlist.size() > 2){
+ if (inputlist.size() > 2) {
std::string covarname = inputlist[2];
- LOG(DEB) << "COVARNAME = " <<covarname << std::endl;
-
+ LOG(DEB) << "COVARNAME = " << covarname << std::endl;
+
TH2D* tempcov = (TH2D*) tempfile->Get(covarname.c_str());
- if (!tempcov){
+ if (!tempcov) {
ERR(FTL) << "Can't find TH2D covar " << covarname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
-
+
throw;
}
-
+
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
-
- for (int i = 0; i < nbins; i++){
- for (int j = 0; j < nbins; j++){
- (*fCovar)(i,j) = tempcov->GetBinContent(i+1,j+1);
+
+ for (int i = 0; i < nbins; i++) {
+ for (int j = 0; j < nbins; j++) {
+ (*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
-
- // Uncorrelated
+
+ // Uncorrelated
} else {
- LOG(SAM) <<"No Covar provided so using diagonal errors for "
- << fName << std::endl;
+ LOG(SAM) << "No Covar provided so using diagonal errors for "
+ << fName << std::endl;
fCovar = NULL;
}
}
-//*******************************************************************************
-void ParamPull::ReadVectFile(std::string input){
+//*******************************************************************************
+void ParamPull::ReadVectFile(std::string input) {
//*******************************************************************************
- std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input,";");
- if (inputlist.size() < 4){
+ std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input, ";");
+ if (inputlist.size() < 4) {
ERR(FTL) << "Need 3 inputs for vector input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
// Open File
std::string rootname = inputlist[0];
- TFile* tempfile = new TFile(rootname.c_str(),"READ");
- if (tempfile->IsZombie()){
+ TFile* tempfile = new TFile(rootname.c_str(), "READ");
+ if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << rootname << std::endl;
throw;
}
// Get Name
std::string tagname = inputlist[1];
// TVector<std::string> dialtags = tempfile->Get(tagname.c_str());
// if (!dialtags){
// ERR(FTL) << "Can't find list of dial names!" << std::endl;
// }
-
- // Get Values
+
+ // Get Values
std::string valuename = inputlist[2];
TVectorD* dialvals = (TVectorD*)tempfile->Get(valuename.c_str());
- if (!dialvals){
+ if (!dialvals) {
ERR(FTL) << "Can't find dial values" << std::endl;
}
- // Get Matrix
+ // Get Matrix
std::string matrixname = inputlist[3];
TMatrixD* matrixvals = (TMatrixD*)tempfile->Get(matrixname.c_str());
- if (!matrixvals){
+ if (!matrixvals) {
ERR(FTL) << "Can't find matirx values" << std::endl;
}
-
+
// Get Types
- if (inputlist.size() > 4){
+ if (inputlist.size() > 4) {
std::string typesname = inputlist[3];
}
// Get Minimum
- if (inputlist.size() > 5){
+ if (inputlist.size() > 5) {
std::string minname = inputlist[4];
}
// Get Maximum
- if (inputlist.size() > 6){
+ if (inputlist.size() > 6) {
std::string maxname = inputlist[5];
}
-
+
}
-//*******************************************************************************
-void ParamPull::ReadDialInput(std::string input){
-//*******************************************************************************
+//*******************************************************************************
+void ParamPull::ReadDialInput(std::string input) {
+//*******************************************************************************
- std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input,";");
- if (inputlist.size() < 3){
+ std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input, ";");
+ if (inputlist.size() < 3) {
ERR(FTL) << "Need 3 inputs for dial input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
- std::vector<double> inputvals = GeneralUtils::ParseToDbl(input,";");
+ std::vector<double> inputvals = GeneralUtils::ParseToDbl(input, ";");
std::string dialname = inputlist[0];
double val = inputvals[1];
double err = inputvals[2];
-
+
fDataHist = new TH1D( (fName + "_data").c_str(),
- (fName + "_data" + fPlotTitles).c_str(), 1, 0, 1);
+ (fName + "_data" + fPlotTitles).c_str(), 1, 0, 1);
fDataHist->SetBinContent(1, val);
fDataHist->SetBinError(1, err);
fDataHist->GetXaxis()->SetBinLabel(1, dialname.c_str());
fCovar = NULL;
}
-//*******************************************************************************
-std::map<std::string, int> ParamPull::GetAllDials(){
-//*******************************************************************************
+//*******************************************************************************
+std::map<std::string, int> ParamPull::GetAllDials() {
+//*******************************************************************************
std::map<std::string, int> dialtypemap;
-
- for (int i = 0; i < fDataHist->GetNbinsX(); i++){
- std::string name = fDataHist->GetXaxis()->GetBinLabel(i+1);
- int type = fTypeHist->GetBinContent(i+1);
-
+ for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
+
+ std::string name = fDataHist->GetXaxis()->GetBinLabel(i + 1);
+ int type = fTypeHist->GetBinContent(i + 1);
+
dialtypemap[name] = type;
-
+
}
return dialtypemap;
}
-//*******************************************************************************
-bool ParamPull::CheckDialsValid(){
+//*******************************************************************************
+bool ParamPull::CheckDialsValid() {
//*******************************************************************************
return true;
std::string helpstring = "";
-
- for (int i = 0; i < fDataHist->GetNbinsX(); i++){
- std::string name = std::string(fDataHist->GetXaxis()->GetBinLabel(i+1));
+
+ for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
+ std::string name = std::string(fDataHist->GetXaxis()->GetBinLabel(i + 1));
// If dial exists its all good
if (FitBase::GetRW()->DialIncluded(name)) continue;
// If it doesn't but its a sample norm also continue
- if (name.find("_norm") != std::string::npos){
+ if (name.find("_norm") != std::string::npos) {
ERR(WRN) << "Norm dial included in covar but not set in FitWeight." << std::endl;
ERR(WRN) << "Assuming its a sample norm and skipping..." << std::endl;
}
// Dial unknown so print a help statement
std::ostringstream tempstr;
tempstr << "unknown_parameter " << name << " "
- << fDataHist->GetBinContent(i+1) << " "
- << fDataHist->GetBinContent(i+1) - fDataHist->GetBinError(i+1) << " "
- << fDataHist->GetBinContent(i+1) + fDataHist->GetBinError(i+1) << " "
- << fDataHist->GetBinError(i+1) << " ";
+ << fDataHist->GetBinContent(i + 1) << " "
+ << fDataHist->GetBinContent(i + 1) - fDataHist->GetBinError(i + 1) << " "
+ << fDataHist->GetBinContent(i + 1) + fDataHist->GetBinError(i + 1) << " "
+ << fDataHist->GetBinError(i + 1) << " ";
if (!fType.empty()) tempstr << fType << std::endl;
else tempstr << "FREE" << std::endl;
- helpstring += tempstr.str();
+ helpstring += tempstr.str();
}
// Show statement before failing
- if (!helpstring.empty()){
-
- ERR(WRN) <<"Dial(s) included in covar but not set in FitWeight." << std::endl
+ if (!helpstring.empty()) {
+
+ ERR(WRN) << "Dial(s) included in covar but not set in FitWeight." << std::endl
<< "ParamPulls needs to know how you want it to be treated." << std::endl
- <<"Include the following lines into your card:" << std::endl;
+ << "Include the following lines into your card:" << std::endl;
ERR(WRN) << helpstring << std::endl;
throw;
return false;
} else {
return true;
}
-
+
}
-
+
//*******************************************************************************
-void ParamPull::Reconfigure(){
+void ParamPull::Reconfigure() {
//*******************************************************************************
-
+
FitWeight* rw = FitBase::GetRW();
-
+
// Get Dial Names that are valid
std::vector<std::string> namevec = rw->GetDialNames();
std::vector<double> valuevec = rw->GetDialValues();
// Set Bin Values from RW
- for (UInt_t i = 0; i < namevec.size(); i++){
+ for (UInt_t i = 0; i < namevec.size(); i++) {
// Loop over bins and check name matches
std::string syst = namevec.at(i);
+ double systval = valuevec.at(i);
+ std::vector<std::string> allsyst = GeneralUtils::ParseToStr(syst, ",");
// Proper Reconf using RW
- for (int j = 0; j < fMCHist->GetNbinsX(); j++){
+ for (int j = 0; j < fMCHist->GetNbinsX(); j++) {
- // If Match set value
- if (!syst.compare( fMCHist->GetXaxis()->GetBinLabel(j+1) )){
-
- double curval = rw->GetDialValue(syst);
- fMCHist->SetBinContent(j+1, curval);
+ // Search for the name of this bin in the corrent dial
+ std::string binname = std::string(fMCHist->GetXaxis()->GetBinLabel(j + 1) );
+
+ // Check Full Name
+ if (!syst.compare(binname.c_str())) {
+ fMCHist->SetBinContent(j + 1, systval);
+ break;
+ }
+
+
+ std::vector<std::string> splitbinname = GeneralUtils::ParseToStr(binname, ",");
+ for (int l = 0; l < splitbinname.size(); l++) {
+ std::string singlebinname = splitbinname[l];
+ for (int k = 0; k < allsyst.size(); k++) {
+ if (!allsyst[k].compare(singlebinname.c_str())) {
+ fMCHist->SetBinContent(j + 1, systval);
+ }
+ }
}
}
+
+
+
}
return;
};
-//*******************************************************************************
-void ParamPull::ResetToy(void){
-//*******************************************************************************
+//*******************************************************************************
+void ParamPull::ResetToy(void) {
+//*******************************************************************************
if (fDataHist) delete fDataHist;
LOG(DEB) << "Resetting toy" << std::endl;
LOG(DEB) << fDataTrue << std::endl;
fDataHist = (TH1D*)fDataTrue->Clone();
LOG(DEB) << "Setting name" << std::endl;
fDataHist->SetNameTitle( (fName + "_data").c_str(),
- (fName + " data" + fPlotTitles).c_str() );
-
+ (fName + " data" + fPlotTitles).c_str() );
+
}
-
-//*******************************************************************************
-void ParamPull::SetFakeData(std::string fakeinput){
-//*******************************************************************************
+
+//*******************************************************************************
+void ParamPull::SetFakeData(std::string fakeinput) {
+//*******************************************************************************
// Set from MC Setting
- if (!fakeinput.compare("MC")){
+ if (!fakeinput.compare("MC")) {
// Copy MC into data
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fMCHist->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
- (fName + " fakedata" + fPlotTitles).c_str() );
+ (fName + " fakedata" + fPlotTitles).c_str() );
// Copy original data errors
- for (int i = 0; i < fDataOrig->GetNbinsX(); i++){
- fDataHist->SetBinError(i+1, fDataOrig->GetBinError(i+1) );
+ for (int i = 0; i < fDataOrig->GetNbinsX(); i++) {
+ fDataHist->SetBinError(i + 1, fDataOrig->GetBinError(i + 1) );
}
// Make True Toy Central Value Hist
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
- (fName + " truedata" + fPlotTitles).c_str() );
-
+ (fName + " truedata" + fPlotTitles).c_str() );
+
} else {
ERR(FTL) << "Trying to set fake data for ParamPulls not from MC!" << std::endl;
ERR(FTL) << "Not currently implemented.." << std::endl;
throw;
-
+
}
}
-//*******************************************************************************
-void ParamPull::RemoveFakeData(){
+//*******************************************************************************
+void ParamPull::RemoveFakeData() {
//*******************************************************************************
delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
- (fName + " data" + fPlotTitles).c_str() );
-
+ (fName + " data" + fPlotTitles).c_str() );
+
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
- (fName + " truedata" + fPlotTitles).c_str() );
+ (fName + " truedata" + fPlotTitles).c_str() );
}
-
-//*******************************************************************************
-double ParamPull::GetLikelihood(){
+
+//*******************************************************************************
+double ParamPull::GetLikelihood() {
//*******************************************************************************
-
+
double like = 0.0;
- switch(fCalcType){
+ switch (fCalcType) {
// Gaussian Calculation with correlations
case kGausPull:
like = StatUtils::GetChi2FromCov(fDataHist, fMCHist, fInvCovar, NULL);
like *= 1E-76;
break;
// Default says this has no pull
case kNoThrow:
default:
like = 0.0;
break;
}
LOG(DEB) << "Likelihood = " << like << " " << fCalcType << std::endl;
return like;
-
+
};
//*******************************************************************************
-int ParamPull::GetNDOF(){
+int ParamPull::GetNDOF() {
//*******************************************************************************
int ndof = 0;
- if (fCalcType != kNoThrow){
+ if (fCalcType != kNoThrow) {
ndof += fDataHist->GetNbinsX();
}
return ndof;
};
-//*******************************************************************************
-void ParamPull::ThrowCovariance(){
+//*******************************************************************************
+void ParamPull::ThrowCovariance() {
//*******************************************************************************
// Reset toy for throw
ResetToy();
LOG(DEB) << "Toy Reset " << std::endl;
-
+
// Generate random Gaussian throws
std::vector<double> randthrows;
- for (int i = 0; i < fDataHist->GetNbinsX(); i++){
+ for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
double randtemp = 0.0;
- switch(fThrowType){
-
- // Gaussian Throws
+ switch (fThrowType) {
+
+ // Gaussian Throws
case kGausThrow:
- randtemp = gRandom->Gaus(0.0,1.0);
+ randtemp = gRandom->Gaus(0.0, 1.0);
break;
-
- // No Throws (DEFAULT)
+
+ // No Throws (DEFAULT)
default:
break;
}
-
+
randthrows.push_back(randtemp);
}
// Create Bin Modifications
double totalres = 0.0;
- for (int i = 0; i < fDataHist->GetNbinsX(); i++){
+ for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
// Calc Bin Mod
double binmod = 0.0;
- for (int j = 0; j < fDataHist->GetNbinsX(); j++){
+ for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
// std::cout << "DECOMP " << j << " " << i << " " << randthrows.at(j) << std::endl;
- binmod += (*fDecomp)(j,i) * randthrows.at(j);
+ binmod += (*fDecomp)(j, i) * randthrows.at(j);
}
// Add up fraction dif
totalres += binmod;
-
+
// Add to current data
- fDataHist->SetBinContent(i+1,fDataHist->GetBinContent(i+1) + binmod);
+ fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) + binmod);
}
// Rename
fDataHist->SetNameTitle( (fName + "_data").c_str(),
- (fName + " toydata" + fPlotTitles).c_str() );
-
+ (fName + " toydata" + fPlotTitles).c_str() );
+
// Print Status
LOG(REC) << "Created new toy histogram. Total Dif = "
- << totalres << std::endl;
+ << totalres << std::endl;
return;
};
//*******************************************************************************
-TH2D ParamPull::GetCovar(){
+TH2D ParamPull::GetCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fInvCovar);
- for (int i = 0; i < tempCov.GetNbinsX(); i++){
- tempCov.GetXaxis()->SetBinLabel(i+1,fDataHist->GetXaxis()->GetBinLabel(i+1));
- tempCov.GetYaxis()->SetBinLabel(i+1,fDataHist->GetXaxis()->GetBinLabel(i+1));
+ for (int i = 0; i < tempCov.GetNbinsX(); i++) {
+ tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
+ tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_INVCOV").c_str(),
- (fName + " InvertedCovariance;Dials;Dials").c_str());
+ (fName + " InvertedCovariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
-TH2D ParamPull::GetFullCovar(){
+TH2D ParamPull::GetFullCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fCovar);
- for (int i = 0; i < tempCov.GetNbinsX(); i++){
- tempCov.GetXaxis()->SetBinLabel(i+1,fDataHist->GetXaxis()->GetBinLabel(i+1));
- tempCov.GetYaxis()->SetBinLabel(i+1,fDataHist->GetXaxis()->GetBinLabel(i+1));
+ for (int i = 0; i < tempCov.GetNbinsX(); i++) {
+ tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
+ tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_COV").c_str(),
- (fName + " Covariance;Dials;Dials").c_str());
+ (fName + " Covariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
-TH2D ParamPull::GetDecompCovar(){
+TH2D ParamPull::GetDecompCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fCovar);
- for (int i = 0; i < tempCov.GetNbinsX(); i++){
- tempCov.GetXaxis()->SetBinLabel(i+1,fDataHist->GetXaxis()->GetBinLabel(i+1));
- tempCov.GetYaxis()->SetBinLabel(i+1,fDataHist->GetXaxis()->GetBinLabel(i+1));
+ for (int i = 0; i < tempCov.GetNbinsX(); i++) {
+ tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
+ tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_DEC").c_str(),
- (fName + " Decomposition;Dials;Dials").c_str());
+ (fName + " Decomposition;Dials;Dials").c_str());
return tempCov;
}
-//*******************************************************************************
-void ParamPull::Write(std::string writeoptt){
-//*******************************************************************************
+//*******************************************************************************
+void ParamPull::Write(std::string writeoptt) {
+//*******************************************************************************
fDataHist->Write();
fMCHist->Write();
-
+
GetCovar().Write();
GetFullCovar().Write();
GetDecompCovar().Write();
-
+
return;
};
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index 50714a6..c40a847 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,242 +1,246 @@
#include "FitWeight.h"
void FitWeight::AddRWEngine(int type) {
switch (type) {
case kNEUT:
fAllRW[type] = new NEUTWeightEngine("neutrw");
break;
case kNUWRO:
fAllRW[type] = new NuWroWeightEngine("nuwrorw");
break;
case kGENIE:
fAllRW[type] = new GENIEWeightEngine("genierw");
break;
case kNORM:
fAllRW[type] = new SampleNormEngine("normrw");
break;
case kLIKEWEIGHT:
fAllRW[type] = new LikelihoodWeightEngine("likerw");
break;
case kT2K:
fAllRW[type] = new T2KWeightEngine("t2krw");
break;
case kCUSTOM:
fAllRW[type] = new NUISANCEWeightEngine("nuisrw");
break;
case kSPLINEPARAMETER:
fAllRW[type] = new SplineWeightEngine("splinerw");
break;
case kNIWG:
fAllRW[type] = new NIWGWeightEngine("niwgrw");
break;
default:
THROW("CANNOT ADD RW Engine for unknown dial type: " << type);
break;
}
}
void FitWeight::IncludeDial(std::string name, std::string type, double val) {
// Should register the dial here.
int typeenum = Reweight::ConvDialType(type);
IncludeDial(name, typeenum, val);
}
void FitWeight::IncludeDial(std::string name, int dialtype, double val) {
// Get the dial enum
int nuisenum = Reweight::ConvDial(name, dialtype);
if (nuisenum == -1){
THROW("NUISENUM == " << nuisenum << " for " << name);
}
// Setup RW Engine Pointer
if (fAllRW.find(dialtype) == fAllRW.end()) {
AddRWEngine(dialtype);
}
WeightEngineBase* rw = fAllRW[dialtype];
// Include the dial
rw->IncludeDial(name, val);
// Set Dial Value
if (val != -9999.9) {
rw->SetDialValue(name, val);
}
// Sort Maps
fAllEnums[name] = nuisenum;
fAllValues[nuisenum] = val;
// Sort Lists
fNameList.push_back(name);
fEnumList.push_back(nuisenum);
fValueList.push_back(val);
}
void FitWeight::Reconfigure(bool silent) {
// Reconfigure all added RW engines
for (std::map<int, WeightEngineBase*>::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
+
+ // Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
SetDialValue(nuisenum, val);
}
// Allow for name aswell using GlobalList to determine sample name.
void FitWeight::SetDialValue(int nuisenum, double val) {
// Conv dial type
int dialtype = int(nuisenum - (nuisenum % 1000)) / 1000;
if (fAllRW.find(dialtype) == fAllRW.end()){
THROW("Cannot find RW Engine for dialtype = " << dialtype << " " << nuisenum << " " << (nuisenum - (nuisenum % 1000)) / 1000);
}
// Get RW Engine for this dial
fAllRW[dialtype]->SetDialValue(nuisenum, val);
fAllValues[nuisenum] = val;
// Update ValueList
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
fValueList[i] = val;
}
}
}
void FitWeight::SetAllDials(const double* x, int n) {
for (size_t i = 0; i < (UInt_t)n; i++) {
int rwenum = fEnumList[i];
SetDialValue(rwenum, x[i]);
}
Reconfigure();
}
double FitWeight::GetDialValue(std::string name) {
+
+ // Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
return GetDialValue(nuisenum);
}
double FitWeight::GetDialValue(int nuisenum) {
return fAllValues[nuisenum];
}
int FitWeight::GetDialPos(std::string name) {
int rwenum = fAllEnums[name];
return GetDialPos(rwenum);
}
int FitWeight::GetDialPos(int nuisenum) {
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
return i;
}
}
ERR(FTL) << "No Dial Found! " << std::endl;
throw;
return -1;
}
bool FitWeight::DialIncluded(std::string name) {
return (fAllEnums.find(name) != fAllEnums.end());
}
bool FitWeight::DialIncluded(int rwenum) {
return (fAllValues.find(rwenum) != fAllValues.end());
}
double FitWeight::CalcWeight(BaseFitEvt* evt) {
double rwweight = 1.0;
for (std::map<int, WeightEngineBase*>::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
// LOG(FIT) << "Iter " << (*iter).second->fCalcName << " = " << w << std::endl;
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double* x) {
size_t count = 0;
for (std::vector<int>::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
SetDialValue( (*iter), x[count] );
count++;
}
}
void FitWeight::GetAllDials(double* x, int n) {
for (int i = 0; i < n; i++) {
x[i] = GetDialValue( fEnumList[i] );
}
}
bool FitWeight::NeedsEventReWeight(const double* x) {
bool haschange = false;
size_t count = 0;
// Compare old to new and decide if RW needed.
for (std::vector<int>::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
int nuisenum = (*iter);
int type = (nuisenum / 1000) - (nuisenum % 1000);
// Compare old to new
double oldval = GetDialValue(nuisenum);
double newval = x[count];
if (oldval != newval) {
if (fAllRW[type]->NeedsEventReWeight()) {
haschange = true;
}
}
count++;
}
return haschange;
}
double FitWeight::GetSampleNorm(std::string name) {
if (name.empty()) return 1.0;
// Find norm dial
if (fAllEnums.find(name + "_norm") != fAllEnums.end()) {
if (fAllValues.find(fAllEnums[name+"_norm"]) != fAllValues.end()){
return fAllValues[ fAllEnums[name+"_norm"] ];
} else {
return 1.0;
}
} else {
return 1.0;
}
}
void FitWeight::Print() {
LOG(REC) << "Fit Weight State: " << std::endl;
for (size_t i = 0; i < fNameList.size(); i++) {
LOG(REC) << " -> Par " << i << ". " << fNameList[i] << " " << fValueList[i] << std::endl;
}
}
diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx
index 2496698..3eef0e4 100644
--- a/src/Reweight/GENIEWeightEngine.cxx
+++ b/src/Reweight/GENIEWeightEngine.cxx
@@ -1,203 +1,232 @@
#include "GENIEWeightEngine.h"
GENIEWeightEngine::GENIEWeightEngine(std::string name) {
#ifdef __GENIE_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
fGenieRW = new genie::rew::GReWeight();
// Get List of Vetos (Just for debugging)
std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fGenieRW_veto");
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos;
bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;
// Now actually add the RW Calcs
if (xsec_ncel)
fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
if (xsec_ccqe)
fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
if (xsec_coh)
fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH);
if (xsec_nnres)
fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
new genie::rew::GReWeightNonResonanceBkg);
if (xsec_nudis)
fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
if (xsec_resdec)
fGenieRW->AdoptWghtCalc("hadro_res_decay",
new genie::rew::GReWeightResonanceDecay);
if (xsec_fzone)
fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
if (xsec_intra)
fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
if (xsec_agky)
fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
if (xsec_qevec)
fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
new genie::rew::GReWeightNuXSecCCQEvec);
#if __GENIE_VERSION__ >= 212
if (xsec_qeaxial)
- fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
- new genie::rew::GReWeightNuXSecCCQEaxial);
+ fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
+ new genie::rew::GReWeightNuXSecCCQEaxial);
#endif
if (xsec_dis)
fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
if (xsec_nc)
fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
if (xsec_ccres)
fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
if (xsec_ncres)
fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
if (xsec_nucqe)
fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);
GReWeightNuXSecCCQE * rwccqe =
dynamic_cast<GReWeightNuXSecCCQE *> (fGenieRW->WghtCalc("xsec_ccqe"));
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
// Default to include shape and normalization changes for CCRES (can be changed downstream if desired)
GReWeightNuXSecCCRES * rwccres =
dynamic_cast<GReWeightNuXSecCCRES *> (fGenieRW->WghtCalc("xsec_ccres"));
rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv);
// Default to include shape and normalization changes for NCRES (can be changed downstream if desired)
GReWeightNuXSecNCRES * rwncres =
dynamic_cast<GReWeightNuXSecNCRES *> (fGenieRW->WghtCalc("xsec_ncres"));
rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
// Default to include shape and normalization changes for DIS (can be changed downstream if desired)
GReWeightNuXSecDIS * rwdis =
dynamic_cast<GReWeightNuXSecDIS *> (fGenieRW->WghtCalc("xsec_dis"));
rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);
// Final Reconfigure
fGenieRW->Reconfigure();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl;
#endif
};
-void GENIEWeightEngine::IncludeDial(int nuisenum, double startval) {
+void GENIEWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __GENIE_ENABLED__
- // Get RW Enum and name
- int rwenum = (nuisenum % 1000);
- genie::rew::GSyst_t rwsyst = static_cast<genie::rew::GSyst_t>(rwenum);
- std::string name = GSyst::AsString(rwsyst);
+ // Get First enum
+ int nuisenum = Reweight::ConvDial(name, kGENIE);
- // Fill Maps
- fGenieNameSysts[name] = rwsyst;
- fGenieEnumSysts[nuisenum] = rwsyst;
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
+ fNameIndex[name]; // = std::vector<size_t>(0);
- // Initialize dial
- fGenieRW->Systematics().Init( fGenieEnumSysts[nuisenum] );
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
- // Add line to check dial is actually handled.
- // if (fGenieRW->Systematics().IsHandled
+ // Get RW
+ genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);
+
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fGENIESysts.push_back(rwsyst);
+
+ // Initialize dial
+ std::cout << "Registering " << singlename << " from " << name << std::endl;
+ fGenieRW->Systematics().Init( fGENIESysts[index] );
+
+ // If Absolute
+ if (fIsAbsTwk) {
+ GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 );
+ }
+
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
- // If Absolute
- if (fIsAbsTwk) {
- GSystUncertainty::Instance()->SetUncertainty( fGenieEnumSysts[nuisenum], 1.0, 1.0 );
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
-
#endif
};
void GENIEWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __GENIE_ENABLED__
- // Set RW engine values
- int rwenum = (nuisenum % 1000);
- fGenieRW->Systematics().Set(static_cast<genie::rew::GSyst_t>(rwenum), val);
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ std::cout << "Setting GENIE Dial Value "
+ << i << " " << indices[i] << " "
+ << fGENIESysts[indices[i]] << std::endl;
+ fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val);
+ }
#endif
}
+void GENIEWeightEngine::SetDialValue(std::string name, double val) {
+#ifdef __GENIE_ENABLED__
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
+ }
+#endif
+}
void GENIEWeightEngine::Reconfigure(bool silent) {
#ifdef __GENIE_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fGenieRW->Reconfigure();
-
+ fGenieRW->Print();
+
// Shout again
if (silent) StartTalking();
#endif
}
double GENIEWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
-
+
#ifdef __GENIE_ENABLED__
// Skip Non GENIE
if (evt->fType != kGENIE) return 1.0;
// Make nom weight
- if (!evt){
- THROW("evt not found : " << evt);
+ if (!evt) {
+ THROW("evt not found : " << evt);
}
- if (!(evt->genie_event)){
- THROW("evt->genie_event not found!" << evt->genie_event);
+ if (!(evt->genie_event)) {
+ THROW("evt->genie_event not found!" << evt->genie_event);
}
- if (!(evt->genie_event->event)){
- THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event));
+ if (!(evt->genie_event->event)) {
+ THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event));
}
- if (!fGenieRW){
- THROW("GENIE RW Not Found!" << fGenieRW);
+ if (!fGenieRW) {
+ THROW("GENIE RW Not Found!" << fGenieRW);
}
rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
// std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl;
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/GENIEWeightEngine.h b/src/Reweight/GENIEWeightEngine.h
index 9d86c23..b0d3d03 100644
--- a/src/Reweight/GENIEWeightEngine.h
+++ b/src/Reweight/GENIEWeightEngine.h
@@ -1,63 +1,64 @@
#ifndef WEIGHT_ENGINE_GENIE_H
#define WEIGHT_ENGINE_GENIE_H
#include "FitLogger.h"
#ifdef __GENIE_ENABLED__
#include "EVGCore/EventRecord.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "GSyst.h"
#include "GSystUncertainty.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "ReWeight/GReWeight.h"
#include "ReWeight/GReWeightAGKY.h"
#include "ReWeight/GReWeightDISNuclMod.h"
#include "ReWeight/GReWeightFGM.h"
#include "ReWeight/GReWeightFZone.h"
#include "ReWeight/GReWeightINuke.h"
#include "ReWeight/GReWeightNonResonanceBkg.h"
#include "ReWeight/GReWeightNuXSecCCQE.h"
#include "ReWeight/GReWeightNuXSecCCQEvec.h"
#include "ReWeight/GReWeightNuXSecCCRES.h"
#include "ReWeight/GReWeightNuXSecCOH.h"
#include "ReWeight/GReWeightNuXSecDIS.h"
#include "ReWeight/GReWeightNuXSecNC.h"
#include "ReWeight/GReWeightNuXSecNCEL.h"
#include "ReWeight/GReWeightNuXSecNCRES.h"
#include "ReWeight/GReWeightResonanceDecay.h"
#if __GENIE_VERSION__ >= 212
#include "ReWeight/GReWeightNuXSecCCQEaxial.h"
#endif
using namespace genie;
using namespace genie::rew;
#endif
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
#include "FitWeight.h"
class GENIEWeightEngine : public WeightEngineBase {
public:
GENIEWeightEngine(std::string name);
~GENIEWeightEngine() {};
- void IncludeDial(int nuisenum, double startval);
+ void IncludeDial(std::string name, double startval);
void SetDialValue(int rwenum, double val);
+ void SetDialValue(std::string name, double val);
+
void Reconfigure(bool silent = false);
double CalcWeight(BaseFitEvt* evt);
inline bool NeedsEventReWeight() { return true; };
#ifdef __GENIE_ENABLED__
- std::map<std::string, genie::rew::GSyst_t> fGenieNameSysts;
- std::map<int, genie::rew::GSyst_t> fGenieEnumSysts;
+ std::vector<genie::rew::GSyst_t> fGENIESysts;
genie::rew::GReWeight* fGenieRW; //!< Genie RW Object
#endif
};
#endif
diff --git a/src/Reweight/LikelihoodWeightEngine.cxx b/src/Reweight/LikelihoodWeightEngine.cxx
index 4640f38..91cc285 100644
--- a/src/Reweight/LikelihoodWeightEngine.cxx
+++ b/src/Reweight/LikelihoodWeightEngine.cxx
@@ -1,58 +1,95 @@
#include "LikelihoodWeightEngine.h"
LikelihoodWeightEngine::LikelihoodWeightEngine(std::string name) {
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up Likelihood Weight RW : " << fCalcName << std::endl;
// Set Abs Twk Config
fIsAbsTwk = true;
};
void LikelihoodWeightEngine::IncludeDial(std::string name, double startval) {
- // Get NUISANCE Enum
- int nuisenum = Reweight::ConvDial(name, kNORM);
-
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(1.0);
-
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
-
- // Set Value if given
- if (startval != -999.9) {
- SetDialValue(name, startval);
- }
+ // Get NUISANCE Enum
+ int nuisenum = Reweight::ConvDial(name, kNORM);
+
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
+ fNameIndex[name]; // = std::vector<size_t>(0);
+
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++){
+ std::string singlename = allnames[i];
+
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(1.0);
+
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
+ }
+
+ // Set Value if given
+ if (startval != -999.9) {
+ SetDialValue(name, startval);
+ }
};
void LikelihoodWeightEngine::SetDialValue(int nuisenum, double val) {
- fValues[fEnumIndex[nuisenum]] = val;
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++){
+ fValues[indices[i]] = val;
+ }
}
void LikelihoodWeightEngine::SetDialValue(std::string name, double val){
- fValues[fNameIndex[name]] = val;
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++){
+ fValues[indices[i]] = val;
+ }
+}
+
+double LikelihoodWeightEngine::GetDialValue(std::string name){
+
+ // Check for exact dial names
+ if (fNameIndex.find(name) != fNameIndex.end()){
+ return fValues[ fNameIndex[name][0] ];
+
+ // If not iterate and check entry in one of the keys
+ } else {
+ for (std::map<std::string, std::vector<size_t> >::iterator iter = fNameIndex.begin();
+ iter != fNameIndex.end(); iter++){
+ std::string keyname = iter->first;
+
+ if ( keyname.find(name) != std::string::npos){
+ return fValues[ iter->second[0] ];
+ }
+ }
+ }
+
+ return -1.0;
}
void LikelihoodWeightEngine::Reconfigure(bool silent) {
- // Empty placeholder incase we want print statements...
+ // Empty placeholder incase we want print statements...
}
diff --git a/src/Reweight/LikelihoodWeightEngine.h b/src/Reweight/LikelihoodWeightEngine.h
index 8ae8917..69dea0a 100644
--- a/src/Reweight/LikelihoodWeightEngine.h
+++ b/src/Reweight/LikelihoodWeightEngine.h
@@ -1,25 +1,26 @@
#ifndef LikelihoodWeightEngine_H
#define LikelihoodWeightEngine_H
#include "FitLogger.h"
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
// Could probably just make a general weight engine for stuff like this...
class LikelihoodWeightEngine : public WeightEngineBase {
public:
LikelihoodWeightEngine(std::string name);
~LikelihoodWeightEngine(){};
void IncludeDial(std::string name, double startval);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
void Reconfigure(bool silent = false);
inline double CalcWeight(BaseFitEvt* evt) {return 1.0;};
inline bool NeedsEventReWeight(){ return false; };
+ double GetDialValue(std::string name);
};
#endif
diff --git a/src/Reweight/NEUTWeightEngine.cxx b/src/Reweight/NEUTWeightEngine.cxx
index 644a93b..3bd0c2e 100644
--- a/src/Reweight/NEUTWeightEngine.cxx
+++ b/src/Reweight/NEUTWeightEngine.cxx
@@ -1,155 +1,186 @@
#include "NEUTWeightEngine.h"
NEUTWeightEngine::NEUTWeightEngine(std::string name) {
#ifdef __NEUT_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up NEUT RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
fNeutRW = new neut::rew::NReWeight();
TDirectory* olddir = gDirectory;
// get list of vetoed calc engines (just for debug really)
std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fNeutRW_veto");
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_res = rw_engine_list.find("xsec_res") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool nucl_casc = rw_engine_list.find("nucl_casc") == std::string::npos;
bool nucl_piless = rw_engine_list.find("nucl_piless") == std::string::npos;
// Activate each calc engine
if (xsec_ccqe)
fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE);
if (xsec_res)
fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES);
if (xsec_ccres)
fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES);
if (xsec_coh)
fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH);
if (xsec_dis)
fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS);
if (xsec_ncel)
fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL);
if (xsec_nc)
fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC);
if (xsec_ncres)
fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES);
if (nucl_casc)
fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc);
if (nucl_piless)
fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless);
fNeutRW->Reconfigure();
olddir->cd();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
THROW("NEUT RW NOT ENABLED!" );
#endif
};
void NEUTWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __NEUT_ENABLED__
- // Get NEUT Syst.
- neut::rew::NSyst_t gensyst = NSyst::FromString(name);
+ // Get First enum
int nuisenum = Reweight::ConvDial(name, kNEUT);
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
- fNEUTSysts.push_back(gensyst);
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
+ fNameIndex[name]; // = std::vector<size_t>(0);
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
- // Initialize dial
- fNeutRW->Systematics().Init( fNEUTSysts[index] );
+ // Get Syst
+ neut::rew::NSyst_t gensyst = NSyst::FromString(singlename);
- // If Absolute
- if (fIsAbsTwk) {
- NSystUncertainty::Instance()->SetUncertainty( fNEUTSysts[index], 1.0, 1.0 );
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fNEUTSysts.push_back(gensyst);
+
+ // Initialize dial
+ LOG(FIT) << "Registering " << singlename << " dial." << std::endl;
+ fNeutRW->Systematics().Init( fNEUTSysts[index] );
+
+ // If Absolute
+ if (fIsAbsTwk) {
+ NSystUncertainty::Instance()->SetUncertainty( fNEUTSysts[index], 1.0, 1.0 );
+ }
+
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
+
#endif
}
void NEUTWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __NEUT_ENABLED__
- fValues[fEnumIndex[nuisenum]] = val;
- fNeutRW->Systematics().Set(fNEUTSysts[fEnumIndex[nuisenum]], val);
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ std::cout << "Setting Dial Value = " << nuisenum << " "
+ << i << " " << indices[i] << " " << fValues[indices[i]]
+ << " Enum=" << fNEUTSysts[indices[i]] << std::endl;
+ fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
+ }
#endif
}
void NEUTWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __NEUT_ENABLED__
- fValues[fNameIndex[name]] = val;
- fNeutRW->Systematics().Set(fNEUTSysts[fNameIndex[name]], val);
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ std::cout << "Setting Dial Value = " << name << " = "
+ << i << " " << indices[i] << " " << fValues[indices[i]]
+ << " Enum=" << fNEUTSysts[indices[i]] << std::endl;
+ fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
+ }
#endif
}
void NEUTWeightEngine::Reconfigure(bool silent) {
#ifdef __NEUT_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fNeutRW->Reconfigure();
+ //if (LOG_LEVEL(DEB)){
+ fNeutRW->Print();
+ // }
+
// Shout again
if (silent) StartTalking();
#endif
}
double NEUTWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __NEUT_ENABLED__
- // Skip Non GENIE
+ // Skip Non NEUT
if (evt->fType != kNEUT) return 1.0;
// Hush now
StopTalking();
// Fill NEUT Common blocks
NEUTUtils::FillNeutCommons(evt->fNeutVect);
// Call Weight calculation
rw_weight = fNeutRW->CalcWeight();
// Speak Now
StartTalking();
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/NIWGWeightEngine.cxx b/src/Reweight/NIWGWeightEngine.cxx
index 2c3a680..94c53ec 100644
--- a/src/Reweight/NIWGWeightEngine.cxx
+++ b/src/Reweight/NIWGWeightEngine.cxx
@@ -1,195 +1,220 @@
#include "NIWGWeightEngine.h"
NIWGWeightEngine::NIWGWeightEngine(std::string name) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up NIWG RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
fNIWGRW = new niwg::rew::NIWGReWeight();
// Get List of Veto Calcs (For Debugging)
std::string rw_engine_list =
FitPar::Config().GetParS("FitWeight.fNIWGRW_veto");
bool niwg_2012a = rw_engine_list.find("niwg_2012a") == std::string::npos;
bool niwg_2014a = rw_engine_list.find("niwg_2014a") == std::string::npos;
bool niwg_pimult = rw_engine_list.find("niwg_pimult") == std::string::npos;
bool niwg_mec = rw_engine_list.find("niwg_mec") == std::string::npos;
bool niwg_rpa = rw_engine_list.find("niwg_rpa") == std::string::npos;
bool niwg_eff_rpa = rw_engine_list.find("niwg_eff_rpa") == std::string::npos;
bool niwg_proton =
rw_engine_list.find("niwg_protonFSIbug") == std::string::npos;
bool niwg_hadron =
rw_engine_list.find("niwg_HadronMultSwitch") == std::string::npos;
// Add the RW Calcs
if (niwg_2012a)
fNIWGRW->AdoptWghtCalc("niwg_2012a", new niwg::rew::NIWGReWeight2012a);
if (niwg_2014a)
fNIWGRW->AdoptWghtCalc("niwg_2014a", new niwg::rew::NIWGReWeight2014a);
if (niwg_pimult)
fNIWGRW->AdoptWghtCalc("niwg_pimult", new niwg::rew::NIWGReWeightPiMult);
if (niwg_mec)
fNIWGRW->AdoptWghtCalc("niwg_mec", new niwg::rew::NIWGReWeightMEC);
if (niwg_rpa)
fNIWGRW->AdoptWghtCalc("niwg_rpa", new niwg::rew::NIWGReWeightRPA);
if (niwg_eff_rpa)
fNIWGRW->AdoptWghtCalc("niwg_eff_rpa",
new niwg::rew::NIWGReWeightEffectiveRPA);
if (niwg_proton)
fNIWGRW->AdoptWghtCalc("niwg_protonFSIbug",
new niwg::rew::NIWGReWeightProtonFSIbug);
if (niwg_hadron)
fNIWGRW->AdoptWghtCalc("niwg_HadronMultSwitch",
new niwg::rew::NIWGReWeightHadronMultSwitch);
fNIWGRW->Reconfigure();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
ERR(FTL) << "NIWG RW Enabled but NEUT RW is not!" << std::endl;
#endif
#else
ERR(FTL) << "NIWG RW NOT ENABLED!" << std::endl;
#endif
};
void NIWGWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
- // Get NEUT Syst.
- niwg::rew::NIWGSyst_t gensyst = niwg::rew::NIWGSyst::FromString(name);
+ // Get First enum
int nuisenum = Reweight::ConvDial(name, kNIWG);
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
- fNIWGSysts.push_back(gensyst);
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
+ fNameIndex[name]; // = std::vector<size_t>(0);
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
- // Initialize dial
- fNIWGRW->Systematics().Init( fNIWGSysts[index] );
+ // Get RW
+ niwg::rew::NIWGSyst_t gensyst = niwg::rew::NIWGSyst::FromString(name);
+
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fNIWGSysts.push_back(gensyst);
+
+
+ // Initialize dial
+ LOG(FIT) << "Registering " << singlename << " from " << name << std::endl;
+ fNIWGRW->Systematics().Init( fNIWGSysts[index] );
+
+ // If Absolute
+ if (fIsAbsTwk) {
+ niwg::rew::NIWGSystUncertainty::Instance()->SetUncertainty( fNIWGSysts[index], 1.0, 1.0 );
+ }
+
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
- // If Absolute
- if (fIsAbsTwk) {
- niwg::rew::NIWGSystUncertainty::Instance()->SetUncertainty( fNIWGSysts[index], 1.0, 1.0 );
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
#endif
}
void NIWGWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
- fValues[fEnumIndex[nuisenum]] = val;
- fNIWGRW->Systematics().Set(fNIWGSysts[fEnumIndex[nuisenum]], val);
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ std::cout << "Setting NIWG Dial Value "
+ << i << " " << indices[i] << " "
+ << fNIWGSysts[indices[i]] << std::endl;
+ fNIWGRW->Systematics().Set( fNIWGSysts[indices[i]], val);
+ }
#endif
#endif
}
void NIWGWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
- fValues[fNameIndex[name]] = val;
- fNIWGRW->Systematics().Set(fNIWGSysts[fNameIndex[name]], val);
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fNIWGRW->Systematics().Set(fNIWGSysts[indices[i]], val);
+ }
#endif
#endif
}
void NIWGWeightEngine::Reconfigure(bool silent) {
#ifdef __NIWG_ENABLED__
#ifdef __NEUT_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fNIWGRW->Reconfigure();
// Shout again
if (silent) StartTalking();
#endif
#endif
}
double NIWGWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __NEUT_ENABLED__
#ifdef __NIWG_ENABLED__
// Skip Non GENIE
if (evt->fType != kNEUT) return 1.0;
// Hush now
StopTalking();
niwg::rew::NIWGEvent* niwg_event = NIWGWeightEngine::GetNIWGEventLocal(evt->fNeutVect);
rw_weight *= fNIWGRW->CalcWeight(*niwg_event);
delete niwg_event;
// Speak Now
StartTalking();
#endif
#endif
// Return rw_weight
return rw_weight;
}
#ifdef __NEUT_ENABLED__
#ifdef __NIWG_ENABLED__
niwg::rew::NIWGEvent * NIWGWeightEngine::GetNIWGEventLocal(NeutVect* nvect)
{
- niwg::rew::NIWGEvent * fDummyNIWGEvent = NULL;
+ niwg::rew::NIWGEvent * fDummyNIWGEvent = NULL;
- fDummyNIWGEvent = new niwg::rew::NIWGEvent();
- fDummyNIWGEvent->detid = 1; // MiniBooNE (apply CCQE LowE variations)
- fDummyNIWGEvent->neutmode = nvect->Mode;
- fDummyNIWGEvent->targetA = nvect->TargetA;
- fDummyNIWGEvent->recenu_ccqe_sk = -1;
- if (nvect->Ibound==0) fDummyNIWGEvent->targetA = 1;//RT: identifies as H, rather than O16
+ fDummyNIWGEvent = new niwg::rew::NIWGEvent();
+ fDummyNIWGEvent->detid = 1; // MiniBooNE (apply CCQE LowE variations)
+ fDummyNIWGEvent->neutmode = nvect->Mode;
+ fDummyNIWGEvent->targetA = nvect->TargetA;
+ fDummyNIWGEvent->recenu_ccqe_sk = -1;
+ if (nvect->Ibound == 0) fDummyNIWGEvent->targetA = 1; //RT: identifies as H, rather than O16
- // Fill initial particle stack
- for (int ip=0; ip<nvect->Npart(); ++ip) {
+ // Fill initial particle stack
+ for (int ip = 0; ip < nvect->Npart(); ++ip) {
- niwg::rew::NIWGPartStack fDummyPartStack;
+ niwg::rew::NIWGPartStack fDummyPartStack;
- fDummyPartStack.p = (nvect->PartInfo(ip)->fP)*0.001; // Convert to GeV
+ fDummyPartStack.p = (nvect->PartInfo(ip)->fP) * 0.001; // Convert to GeV
- fDummyPartStack.pdg = nvect->PartInfo(ip)->fPID;
- fDummyPartStack.chase = nvect->PartInfo(ip)->fIsAlive;
- fDummyPartStack.parent = nvect->ParentIdx(ip)-1; // WARNING: this needs to be tested with a NeutRoot file
+ fDummyPartStack.pdg = nvect->PartInfo(ip)->fPID;
+ fDummyPartStack.chase = nvect->PartInfo(ip)->fIsAlive;
+ fDummyPartStack.parent = nvect->ParentIdx(ip) - 1; // WARNING: this needs to be tested with a NeutRoot file
- fDummyNIWGEvent->part_stack.push_back(fDummyPartStack);
- }
- fDummyNIWGEvent->CalcKinematics();
+ fDummyNIWGEvent->part_stack.push_back(fDummyPartStack);
+ }
+ fDummyNIWGEvent->CalcKinematics();
- return fDummyNIWGEvent;
+ return fDummyNIWGEvent;
}
#endif
#endif
diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx
index 1a45507..44b301a 100644
--- a/src/Reweight/NUISANCEWeightEngine.cxx
+++ b/src/Reweight/NUISANCEWeightEngine.cxx
@@ -1,106 +1,129 @@
#include "NUISANCEWeightEngine.h"
#include "NUISANCEWeightCalcs.h"
NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) {
- // Setup the NUISANCE Reweight engine
- fCalcName = name;
- LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl;
+ // Setup the NUISANCE Reweight engine
+ fCalcName = name;
+ LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl;
- // Load in all Weight Calculations
- fWeightCalculators.push_back( new GaussianModeCorr() );
+ // Load in all Weight Calculations
+ fWeightCalculators.push_back( new GaussianModeCorr() );
- // Set Abs Twk Config
- fIsAbsTwk = true;
+ // Set Abs Twk Config
+ fIsAbsTwk = true;
};
void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) {
- // Get NUISANCE Enum
- int nuisenum = Reweight::ConvDial(name, kCUSTOM);
+ // Get First enum
+ int nuisenum = Reweight::ConvDial(name, kCUSTOM);
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(1.0);
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
+ fNameIndex[name]; // = std::vector<size_t>(0);
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
- // Set Value if given
- if (startval != -999.9) {
- SetDialValue(name, startval);
- }
+ // Get RW
+ int singleenum = Reweight::ConvDial(singlename, kCUSTOM);
+
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fNUISANCEEnums.push_back(singleenum);
+
+ // Initialize dial
+ std::cout << "Registering " << singlename << " from " << name << std::endl;
+
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
+ }
+
+ // Set Value if given
+ if (startval != -999.9) {
+ SetDialValue(nuisenum, startval);
+ }
};
void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) {
- fValues[fEnumIndex[nuisenum]] = val;
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ }
}
void NUISANCEWeightEngine::SetDialValue(std::string name, double val) {
- fValues[fNameIndex[name]] = val;
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ }
}
void NUISANCEWeightEngine::Reconfigure(bool silent) {
- // Loop over all names
- for (std::map<int, size_t>::iterator enumiter = fEnumIndex.begin();
- enumiter != fEnumIndex.end(); enumiter++) {
+ for (size_t i = 0; i < fNUISANCEEnums.size(); i++) {
+
+ for (std::vector<NUISANCEWeightCalc*>::iterator calciter = fWeightCalculators.begin();
+ calciter != fWeightCalculators.end(); calciter++) {
- for (std::vector<NUISANCEWeightCalc*>::iterator calciter = fWeightCalculators.begin();
- calciter != fWeightCalculators.end(); calciter++) {
+ NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc*>(*calciter);
- NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc*>(*calciter);
- if (nuiscalc->IsHandled(enumiter->first)) {
- nuiscalc->SetDialValue(enumiter->first, fValues[enumiter->second]);
- }
- }
- }
+ if (nuiscalc->IsHandled(fNUISANCEEnums[i])) {
+ nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]);
+ }
+ }
+ }
}
-double NUISANCEWeightEngine::CalcWeight(BaseFitEvt* evt) {
- double rw_weight = 1.0;
+double NUISANCEWeightEngine::CalcWeight(BaseFitEvt * evt) {
+ double rw_weight = 1.0;
- // Cast as usable class
- for (std::vector<NUISANCEWeightCalc*>::iterator iter = fWeightCalculators.begin();
- iter != fWeightCalculators.end(); iter++) {
- NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc*>(*iter);
+ // Cast as usable class
+ for (std::vector<NUISANCEWeightCalc*>::iterator iter = fWeightCalculators.begin();
+ iter != fWeightCalculators.end(); iter++) {
+ NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc*>(*iter);
- rw_weight *= nuiscalc->CalcWeight(evt);
- }
+ rw_weight *= nuiscalc->CalcWeight(evt);
+ }
- // Return rw_weight
- return rw_weight;
+ // Return rw_weight
+ return rw_weight;
}
diff --git a/src/Reweight/NUISANCEWeightEngine.h b/src/Reweight/NUISANCEWeightEngine.h
index bd19b60..07e95b8 100644
--- a/src/Reweight/NUISANCEWeightEngine.h
+++ b/src/Reweight/NUISANCEWeightEngine.h
@@ -1,33 +1,35 @@
#ifndef WEIGHT_ENGINE_NUISANCE_H
#define WEIGHT_ENGINE_NUISANCE_H
#include "FitLogger.h"
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
#include "FitWeight.h"
#include "NUISANCESyst.h"
#include "NUISANCEWeightCalcs.h"
class NUISANCEWeightEngine : public WeightEngineBase {
public:
NUISANCEWeightEngine(std::string name);
~NUISANCEWeightEngine() {};
void IncludeDial(std::string name, double startval);
void SetDialValue(std::string name, double val);
void SetDialValue(int nuisenum, double val);
void Reconfigure(bool silent = false);
double CalcWeight(BaseFitEvt* evt);
inline bool NeedsEventReWeight() { return true; };
std::vector<NUISANCEWeightCalc*> fWeightCalculators;
+ std::vector<int> fNUISANCEEnums;
+
};
#endif
diff --git a/src/Reweight/NuWroWeightEngine.cxx b/src/Reweight/NuWroWeightEngine.cxx
index baab616..53694c9 100644
--- a/src/Reweight/NuWroWeightEngine.cxx
+++ b/src/Reweight/NuWroWeightEngine.cxx
@@ -1,126 +1,140 @@
#include "NuWroWeightEngine.h"
NuWroWeightEngine::NuWroWeightEngine(std::string name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up NuWro RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
// Create the engine
fNuwroRW = new nuwro::rew::NuwroReWeight();
// Get List of Veto Calcs (For Debugging)
std::string rw_engine_list =
FitPar::Config().GetParS("FitWeight.fNuwroRW_veto");
bool xsec_qel = rw_engine_list.find("nuwro_QEL") == std::string::npos;
bool xsec_flag = rw_engine_list.find("nuwro_FlagNorm") == std::string::npos;
bool xsec_res = rw_engine_list.find("nuwro_RES") == std::string::npos;
// Add the RW Calcs
if (xsec_qel)
fNuwroRW->AdoptWghtCalc("nuwro_QEL", new nuwro::rew::NuwroReWeight_QEL);
if (xsec_flag)
fNuwroRW->AdoptWghtCalc("nuwro_FlagNorm",
new nuwro::rew::NuwroReWeight_FlagNorm);
if (xsec_res) fNuwroRW->AdoptWghtCalc( "nuwro_RES", new
nuwro::rew::NuwroReWeight_SPP );
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
ERR(FTL) << "NUWRO RW NOT ENABLED! " << std::endl;
#endif
};
void NuWroWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Get RW Enum and name
- nuwro::rew::NuwroSyst_t gensyst = nuwro::rew::NuwroSyst::FromString(name);
int nuisenum = Reweight::ConvDial(name, kNUWRO);
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
- fNUWROSysts.push_back(gensyst);
+ // Initialise new vector
+ fEnumIndex[nuisenum];
+ fNameIndex[name];
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
- // Initialise Dial
- fNuwroRW->Systematics().Add( fNUWROSysts[index] );
+ // Get Syst
+ nuwro::rew::NuwroSyst_t gensyst = nuwro::rew::NuwroSyst::FromString(name);
- // If Absolute
- if (fIsAbsTwk) {
- nuwro::rew::NuwroSystUncertainty::Instance()->SetUncertainty( fNUWROSysts[index], 1.0, 1.0 );
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fNUWROSysts.push_back(gensyst);
+
+ // Initialise Dial
+ fNuwroRW->Systematics().Add( fNUWROSysts[index] );
+
+ if (fIsAbsTwk) {
+ nuwro::rew::NuwroSystUncertainty::Instance()->SetUncertainty( fNUWROSysts[index], 1.0, 1.0 );
+ }
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(name, startval);
}
#endif
};
void NuWroWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
- fValues[fEnumIndex[nuisenum]] = val;
- fNuwroRW->Systematics().SetSystVal(fNUWROSysts[fEnumIndex[nuisenum]], val);
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fNuwroRW->Systematics().SetSystVal(fNUWROSysts[indices[i]], val);
+ }
#endif
}
void NuWroWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
- fValues[fNameIndex[name]] = val;
- fNuwroRW->Systematics().SetSystVal(fNUWROSysts[fNameIndex[name]], val);
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fNuwroRW->Systematics().SetSystVal(fNUWROSysts[indices[i]], val);
+ }
#endif
}
void NuWroWeightEngine::Reconfigure(bool silent) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fNuwroRW->Reconfigure();
// Shout again
if (silent) StartTalking();
#endif
}
double NuWroWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __NUWRO_REWEIGHT_ENABLED__
// Skip Non GENIE
if (evt->fType != kNUWRO) return 1.0;
// Call Weight calculation
rw_weight = fNuwroRW->CalcWeight(evt->fNuwroEvent);
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/SampleNormEngine.cxx b/src/Reweight/SampleNormEngine.cxx
index 13d2acc..a9b4354 100644
--- a/src/Reweight/SampleNormEngine.cxx
+++ b/src/Reweight/SampleNormEngine.cxx
@@ -1,58 +1,94 @@
#include "SampleNormEngine.h"
SampleNormEngine::SampleNormEngine(std::string name) {
// Setup the NEUT Reweight engien
fCalcName = name;
- LOG(FIT) << "Setting up Sample Norm RW : " << fCalcName << std::endl;
+ LOG(FIT) << "Setting up Likelihood Weight RW : " << fCalcName << std::endl;
// Set Abs Twk Config
fIsAbsTwk = true;
};
void SampleNormEngine::IncludeDial(std::string name, double startval) {
- // Get NUISANCE Enum
- int nuisenum = Reweight::ConvDial(name, kNORM);
-
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(1.0);
-
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
-
- // Set Value if given
- if (startval != -999.9) {
- SetDialValue(name, startval);
- }
+ // Get NUISANCE Enum
+ int nuisenum = Reweight::ConvDial(name, kNORM);
+
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
+ fNameIndex[name]; // = std::vector<size_t>(0);
+
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++){
+ std::string singlename = allnames[i];
+
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(1.0);
+
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
+ }
+
+ // Set Value if given
+ if (startval != -999.9) {
+ SetDialValue(name, startval);
+ }
};
void SampleNormEngine::SetDialValue(int nuisenum, double val) {
- fValues[fEnumIndex[nuisenum]] = val;
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++){
+ fValues[indices[i]] = val;
+ }
}
void SampleNormEngine::SetDialValue(std::string name, double val){
- fValues[fNameIndex[name]] = val;
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++){
+ fValues[indices[i]] = val;
+ }
+}
+
+double SampleNormEngine::GetDialValue(std::string name){
+
+ // Check for exact dial names
+ if (fNameIndex.find(name) != fNameIndex.end()){
+ return fValues[ fNameIndex[name][0] ];
+
+ // If not iterate and check entry in one of the keys
+ } else {
+ for (std::map<std::string, std::vector<size_t> >::iterator iter = fNameIndex.begin();
+ iter != fNameIndex.end(); iter++){
+ std::string keyname = iter->first;
+ if (keyname.find(name) != std::string::npos){
+ return fValues[ iter->second[0] ];
+ }
+ }
+ }
+
+ return -1.0;
}
void SampleNormEngine::Reconfigure(bool silent) {
- // Empty placeholder incase we want print statements...
+ // Empty placeholder incase we want print statements...
}
diff --git a/src/Reweight/SampleNormEngine.h b/src/Reweight/SampleNormEngine.h
index c341a71..3d1bd48 100644
--- a/src/Reweight/SampleNormEngine.h
+++ b/src/Reweight/SampleNormEngine.h
@@ -1,22 +1,24 @@
#ifndef SampleNormEngine_H
#define SampleNormEngine_H
#include "FitLogger.h"
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
class SampleNormEngine : public WeightEngineBase {
public:
SampleNormEngine(std::string name);
~SampleNormEngine(){};
void IncludeDial(std::string name, double startval);
void SetDialValue(int rwenum, double val);
void SetDialValue(std::string name, double val);
void Reconfigure(bool silent = false);
inline double CalcWeight(BaseFitEvt* evt) {return 1.0;};
inline bool NeedsEventReWeight(){ return false; };
+
+ double GetDialValue(std::string name);
};
#endif
diff --git a/src/Reweight/SplineWeightEngine.cxx b/src/Reweight/SplineWeightEngine.cxx
index a2ec628..6e1df33 100644
--- a/src/Reweight/SplineWeightEngine.cxx
+++ b/src/Reweight/SplineWeightEngine.cxx
@@ -1,79 +1,94 @@
#include "SplineWeightEngine.h"
SplineWeightEngine::SplineWeightEngine(std::string name) {
// Setup the Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up Spline RW : " << fCalcName << std::endl;
// Set Abs Twk Config
fIsAbsTwk = true;
};
void SplineWeightEngine::IncludeDial(std::string name, double startval) {
- // Get NUISANCE Enum
+ // Get RW Enum and name
int nuisenum = Reweight::ConvDial(name, kSPLINEPARAMETER);
-
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
-
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
-
- // std::cout << "Inlcuded Spline Dial " << name << " " << nuisenum << " " << startval << " " << index << std::endl;
-
+
+ // Initialise new vector
+ fEnumIndex[nuisenum];
+ fNameIndex[name];
+
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
+
+ // Get Single Enum For this dial
+ int singleenum = Reweight::ConvDial(singlename, kCUSTOM);
+
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fSingleEnums.push_back(singleenum);
+ fSingleNames.push_back(singlename);
+
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
+ }
+
// Set Value if given
if (startval != -999.9) {
SetDialValue(name, startval);
}
}
void SplineWeightEngine::SetDialValue(int nuisenum, double val) {
- // LOG(FIT) << "Enum Val " << nuisenum << std::endl;
- // LOG(FIT) << fEnumIndex.size() << std::endl;
- // LOG(FIT) << "Set Dial Value to " << nuisenum << " " << fEnumIndex[nuisenum] << " " << val << std::endl;
- fValues[fEnumIndex[nuisenum]] = val;
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ }
}
void SplineWeightEngine::SetDialValue(std::string name, double val){
- fValues[fNameIndex[name]] = val;
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ }
}
void SplineWeightEngine::Reconfigure(bool silent) {
- for (std::map<std::string, size_t>::iterator iter = fNameIndex.begin();
- iter != fNameIndex.end(); iter++){
- // LOG(FIT) << "Reconfiguring Spline " << iter->first << " to be " << fValues[ iter->second ] << " Inside SPL RW" << std::endl;
- fSplineValueMap[ iter->first ] = fValues[ iter->second ];
+ for (size_t i = 0; i < fSingleNames.size(); i++){
+ fSplineValueMap[ fSingleNames[i] ] = fValues[i];
}
}
double SplineWeightEngine::CalcWeight(BaseFitEvt* evt) {
if (!evt->fSplineRead) return 1.0;
if (evt->fSplineRead->NeedsReconfigure()) {
evt->fSplineRead->Reconfigure(fSplineValueMap);
}
double rw_weight = evt->fSplineRead->CalcWeight( evt->fSplineCoeff );
if (rw_weight < 0.0) rw_weight = 0.0;
return rw_weight;
}
diff --git a/src/Reweight/SplineWeightEngine.h b/src/Reweight/SplineWeightEngine.h
index 974f9eb..9de1afa 100644
--- a/src/Reweight/SplineWeightEngine.h
+++ b/src/Reweight/SplineWeightEngine.h
@@ -1,24 +1,26 @@
#ifndef SplineWeightEngine_H
#define SplineWeightEngine_H
#include "FitLogger.h"
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
#include "SplineReader.h"
class SplineWeightEngine : public WeightEngineBase {
public:
SplineWeightEngine(std::string name);
~SplineWeightEngine(){};
void IncludeDial(std::string name, double startval);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
void Reconfigure(bool silent = false);
inline double CalcWeight(BaseFitEvt* evt);
inline bool NeedsEventReWeight(){ return true; };
std::map< std::string, double > fSplineValueMap;
+ std::vector<int> fSingleEnums;
+ std::vector<std::string> fSingleNames;
};
#endif
diff --git a/src/Reweight/T2KWeightEngine.cxx b/src/Reweight/T2KWeightEngine.cxx
index 8d1c761..ac0123b 100644
--- a/src/Reweight/T2KWeightEngine.cxx
+++ b/src/Reweight/T2KWeightEngine.cxx
@@ -1,126 +1,149 @@
#include "T2KWeightEngine.h"
T2KWeightEngine::T2KWeightEngine(std::string name) {
#ifdef __T2KREW_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up T2K RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
// Create Main RW Engine
fT2KRW = new t2krew::T2KReWeight();
// Setup Sub RW Engines (Only activated for neut and niwg)
fT2KNeutRW = new t2krew::T2KNeutReWeight();
fT2KNIWGRW = new t2krew::T2KNIWGReWeight();
fT2KRW->AdoptWghtEngine("fNeutRW", fT2KNeutRW);
fT2KRW->AdoptWghtEngine("fNIWGRW", fT2KNIWGRW);
fT2KRW->Reconfigure();
// allow cout again
StartTalking();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
-
-#else
+
+#else
ERR(FTL) << "T2K RW NOT ENABLED" << std::endl;
throw;
#endif
};
void T2KWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __T2KREW_ENABLED__
- // Get NEUT Syst.
- t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name);
+
+
+ // Get First enum
int nuisenum = Reweight::ConvDial(name, kT2K);
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
- fT2KSysts.push_back(gensyst);
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
+ fNameIndex[name]; // = std::vector<size_t>(0);
+
+ // Split by commas
+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
+
+ // Get RW
+ t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name);
- fEnumIndex[nuisenum] = index;
- fNameIndex[name] = index;
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fT2KSysts.push_back(gensyst);
- // Initialize dial
- fT2KRW->Systematics().Include(gensyst);
+ // Initialize dial
+ std::cout << "Registering " << singlename << " from " << name << std::endl;
+ fT2KRW->Systematics().Include(gensyst);
+
+ // If Absolute
+ if (fIsAbsTwk) {
+ fT2KRW->Systematics().SetAbsTwk(gensyst);
+ }
+
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
- // If Absolute
- if (fIsAbsTwk) {
- fT2KRW->Systematics().SetAbsTwk(gensyst);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
}
void T2KWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __T2KREW_ENABLED__
- fValues[fEnumIndex[nuisenum]] = val;
- fT2KRW->Systematics().SetTwkDial(fT2KSysts[fEnumIndex[nuisenum]], val);
+ std::vector<size_t> indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
+ }
#endif
}
void T2KWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __T2KREW_ENABLED__
- fValues[fNameIndex[name]] = val;
- fT2KRW->Systematics().SetTwkDial(fT2KSysts[fNameIndex[name]], val);
+ std::vector<size_t> indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
+ }
#endif
}
void T2KWeightEngine::Reconfigure(bool silent) {
#ifdef __T2KREW_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
StopTalking();
fT2KRW->Reconfigure();
StartTalking();
// Shout again
if (silent) StartTalking();
#endif
}
double T2KWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __T2KREW_ENABLED__
// Skip Non GENIE
if (evt->fType != kNEUT) return 1.0;
// Hush now
StopTalking();
// Get Weight For NEUT
rw_weight = fT2KRW->CalcWeight(evt->fNeutVect);
// Speak Now
StartTalking();
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/WeightEngineBase.cxx b/src/Reweight/WeightEngineBase.cxx
index 038a896..04bb6c0 100644
--- a/src/Reweight/WeightEngineBase.cxx
+++ b/src/Reweight/WeightEngineBase.cxx
@@ -1,25 +1,26 @@
#include "WeightEngineBase.h"
bool WeightEngineBase::IsDialIncluded(std::string name) {
return (fNameIndex.find(name) != fNameIndex.end());
}
bool WeightEngineBase::IsDialIncluded(int nuisenum) {
return (fEnumIndex.find(nuisenum) != fEnumIndex.end());
}
double WeightEngineBase::GetDialValue(std::string name) {
if (!IsDialIncluded(name)) {
ERR(FTL) << "Dial " << name
<< " not included in " << fCalcName << std::endl;
}
- return fValues[fNameIndex[name]];
+ return fValues[fNameIndex[name][0]];
}
double WeightEngineBase::GetDialValue(int nuisenum) {
if (!IsDialIncluded(nuisenum)) {
ERR(FTL) << "Dial Enum " << nuisenum
<< " not included in " << fCalcName << std::endl;
}
- return fValues[fEnumIndex[nuisenum]];
+ return fValues[fEnumIndex[nuisenum][0]];
}
+
diff --git a/src/Reweight/WeightEngineBase.h b/src/Reweight/WeightEngineBase.h
index c8adf71..a3a3c7e 100644
--- a/src/Reweight/WeightEngineBase.h
+++ b/src/Reweight/WeightEngineBase.h
@@ -1,59 +1,60 @@
#ifndef WEIGHTENGINE_BASE_H
#define WEIGHTENGINE_BASE_H
#include "BaseFitEvt.h"
#include "FitLogger.h"
#include "FitUtils.h"
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <deque>
#include <iomanip>
#include <iostream>
#include <map>
#include <numeric>
#include <sstream>
#include <string>
#include <vector>
#include "GeneratorUtils.h"
#include "TCanvas.h"
#include "TGraph2D.h"
#include "WeightUtils.h"
class WeightEngineBase {
public:
WeightEngineBase(){};
virtual ~WeightEngineBase(){};
// Functions requiring Override
virtual void IncludeDial(std::string name, double startval){};
virtual void SetDialValue(int nuisenum, double val){};
virtual void SetDialValue(std::string name, double val){};
virtual bool IsDialIncluded(std::string name);
virtual bool IsDialIncluded(int nuisenum);
virtual double GetDialValue(std::string name);
virtual double GetDialValue(int nuisenum);
virtual void Reconfigure(bool silent){};
virtual double CalcWeight(BaseFitEvt* evt){ return 1.0; };
virtual bool NeedsEventReWeight() = 0;
bool fHasChanged;
bool fIsAbsTwk;
- std::vector<double> fValues;
- std::map<int, size_t> fEnumIndex;
- std::map<std::string, size_t> fNameIndex;
+ std::vector< double > fValues;
+ std::map<int, std::vector<size_t> > fEnumIndex;
+ std::map<std::string, std::vector<size_t> > fNameIndex;
+
std::string fCalcName;
};
#endif
diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index 81d4fb3..758f91c 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,487 +1,489 @@
#include "WeightUtils.h"
//********************************************************************
TF1 FitBase::GetRWConvFunction(std::string type, std::string name) {
//********************************************************************
std::string dialfunc = "x";
std::string parType = type;
double low = -10000.0;
double high = 10000.0;
if (parType.find("parameter") == std::string::npos) parType += "_parameter";
string line;
ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card").c_str(),
ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 4) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
// inputlist[2] should be the units... ignore for now
dialfunc = inputlist[3];
// High and low are optional, check whether they exist
if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]);
if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]);
}
TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high);
return convfunc;
}
//********************************************************************
std::string FitBase::GetRWUnits(std::string type, std::string name) {
//********************************************************************
std::string unit = "sig.";
std::string parType = type;
if (parType.find("parameter") == std::string::npos) {
parType += "_parameter";
}
std::string line;
std::ifstream card((GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card").c_str(), ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 3) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
unit = inputlist[2];
break;
}
return unit;
}
//********************************************************************
double FitBase::RWAbsToSigma(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX(val);
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToAbs(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val);
return conv_val;
}
//********************************************************************
double FitBase::RWFracToSigma(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX((val * f1.Eval(0.0)));
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToFrac(std::string type, std::string name, double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val) / f1.Eval(0.0);
return conv_val;
}
int FitBase::ConvDialType(std::string type) {
if (!type.compare("neut_parameter")) return kNEUT;
else if (!type.compare("niwg_parameter")) return kNIWG;
else if (!type.compare("nuwro_parameter")) return kNUWRO;
else if (!type.compare("t2k_parameter")) return kT2K;
else if (!type.compare("genie_parameter")) return kGENIE;
else if (!type.compare("custom_parameter")) return kCUSTOM;
else if (!type.compare("norm_parameter")) return kNORM;
else if (!type.compare("modenorm_parameter")) return kMODENORM;
else if (!type.compare("likeweight_parameter")) return kLIKEWEIGHT;
else if (!type.compare("spline_parameter")) return kSPLINEPARAMETER;
else return kUNKNOWN;
}
std::string FitBase::ConvDialType(int type) {
switch (type) {
case kNEUT: { return "neut_parameter"; }
case kNIWG: { return "niwg_parameter"; }
case kNUWRO: { return "nuwro_parameter"; }
case kT2K: { return "t2k_parameter"; }
case kGENIE: { return "genie_parameter"; }
case kNORM: { return "norm_parameter"; }
case kCUSTOM: {return "custom_parameter"; }
case kMODENORM: { return "modenorm_parameter"; }
case kLIKEWEIGHT: { return "likeweight_parameter"; }
case kSPLINEPARAMETER: {return "spline_parameter";}
default: return "unknown_parameter";
}
}
int FitBase::GetDialEnum(std::string type, std::string name) {
return FitBase::GetDialEnum( FitBase::ConvDialType(type), name );
}
int FitBase::GetDialEnum(int type, std::string name) {
int offset = type * 1000;
int this_enum = -1; //Not Found
std::cout << "Getting dial enum " << type << " " << name << std::endl;
// Select Types
switch (type) {
// NEUT DIAL TYPE
case kNEUT: {
#ifdef __NEUT_ENABLED__
int neut_enum = (int)neut::rew::NSyst::FromString(name);
if (neut_enum != 0) { this_enum = neut_enum + offset; }
#else
this_enum = -2; //Not enabled
#endif
break;
}
// NIWG DIAL TYPE
case kNIWG: {
#ifdef __NIWG_ENABLED__
int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
if (niwg_enum != 0) { this_enum = niwg_enum + offset; }
#else
this_enum = -2;
#endif
break;
}
// NUWRO DIAL TYPE
case kNUWRO: {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
if (nuwro_enum > 0) { this_enum = nuwro_enum + offset; }
#else
this_enum = -2;
#endif
}
// GENIE DIAL TYPE
case kGENIE: {
#ifdef __GENIE_ENABLED__
int genie_enum = (int)genie::rew::GSyst::FromString(name);
if (genie_enum > 0) { this_enum = genie_enum + offset; }
#else
this_enum = -2;
#endif
break;
}
case kCUSTOM: {
int custom_enum = 0; // PLACEHOLDER
this_enum = custom_enum + offset;
break;
}
// T2K DIAL TYPE
case kT2K: {
#ifdef __T2KREW_ENABLED__
int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
if (t2k_enum > 0) { this_enum = t2k_enum + offset; }
#else
this_enum = -2;
#endif
break;
}
case kNORM: {
if (gNormEnums.find(name) == gNormEnums.end()) {
gNormEnums[name] = gNormEnums.size() + 1 + offset;
}
this_enum = gNormEnums[name];
break;
}
case kMODENORM: {
size_t us_pos = name.find_first_of('_');
std::string numstr = name.substr(us_pos + 1);
int mode_num = std::atoi(numstr.c_str());
LOG(FTL) << "Getting mode num " << mode_num << std::endl;
if (!mode_num) {
ERR(FTL) << "Attempting to parse dial name: \"" << name
<< "\" as a mode norm dial but failed." << std::endl;
throw;
}
this_enum = 60 + mode_num + offset;
break;
}
case kLIKEWEIGHT: {
if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
}
this_enum = gLikeWeightEnums[name];
break;
}
case kSPLINEPARAMETER: {
if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
}
this_enum = gSplineParameterEnums[name];
}
}
// If Not Enabled
if (this_enum == -2) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
}
// If Not Found
if (this_enum == -1) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
}
return this_enum;
}
int Reweight::ConvDialType(std::string type) {
if (!type.compare("neut_parameter")) return kNEUT;
else if (!type.compare("niwg_parameter")) return kNIWG;
else if (!type.compare("nuwro_parameter")) return kNUWRO;
else if (!type.compare("t2k_parameter")) return kT2K;
else if (!type.compare("genie_parameter")) return kGENIE;
else if (!type.compare("norm_parameter")) return kNORM;
else if (!type.compare("modenorm_parameter")) return kMODENORM;
else if (!type.compare("custom_parameter")) return kCUSTOM;
else if (!type.compare("likeweight_parameter")) return kLIKEWEIGHT;
else if (!type.compare("spline_parameter")) return kSPLINEPARAMETER;
else return kUNKNOWN;
}
std::string Reweight::ConvDialType(int type) {
switch (type) {
case kNEUT: { return "neut_parameter"; }
case kNIWG: { return "niwg_parameter"; }
case kNUWRO: { return "nuwro_parameter"; }
case kT2K: { return "t2k_parameter"; }
case kGENIE: { return "genie_parameter"; }
case kNORM: { return "norm_parameter"; }
case kCUSTOM: {return "custom_parameter"; }
case kMODENORM: { return "modenorm_parameter"; }
case kLIKEWEIGHT: { return "likeweight_parameter"; }
case kSPLINEPARAMETER: {return "spline_parameter";}
default: return "unknown_parameter";
}
}
int Reweight::NEUTEnumFromName(std::string name) {
#ifdef __NEUT_ENABLED__
int neutenum = (int)neut::rew::NSyst::FromString(name);
return (neutenum > 0) ? neutenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::NIWGEnumFromName(std::string name) {
#ifdef __NIWG_ENABLED__
int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name);
return (niwgenum != 0) ? niwgenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::NUWROEnumFromName(std::string name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name);
return (nuwroenum > 0) ? nuwroenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::GENIEEnumFromName(std::string name) {
#ifdef __GENIE_ENABLED__
int genieenum = (int)genie::rew::GSyst::FromString(name);
return (genieenum > 0) ? genieenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::T2KEnumFromName(std::string name) {
#ifdef __T2KREW_ENABLED__
int t2kenum = (int)t2krew::T2KSyst::FromString(name);
return (t2kenum > 0) ? t2kenum : kNoDialFound;
#else
return kGeneratorNotBuilt;
#endif
}
int Reweight::NUISANCEEnumFromName(std::string name, int type) {
int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type);
return nuisenum;
}
int Reweight::CustomEnumFromName(std::string name) {
int custenum = Reweight::ConvertNUISANCEDial(name);
return custenum;
}
int Reweight::ConvDial(std::string name, std::string type, bool exceptions) {
return Reweight::ConvDial( name, Reweight::ConvDialType(type), exceptions );
}
-int Reweight::ConvDial(std::string name, int type, bool exceptions) {
+int Reweight::ConvDial(std::string fullname, int type, bool exceptions) {
+
+ std::string name = GeneralUtils::ParseToStr(fullname,",")[0]; // Only use first dial given
// Produce offset seperating each type.
int offset = type * 1000;
int genenum = kNoDialFound;
switch (type) {
case kNEUT:
genenum = NEUTEnumFromName(name);
break;
case kNIWG:
genenum = NIWGEnumFromName(name);
break;
case kNUWRO:
genenum = NUWROEnumFromName(name);
break;
case kGENIE:
genenum = GENIEEnumFromName(name);
break;
case kT2K:
genenum = T2KEnumFromName(name);
break;
case kCUSTOM:
genenum = CustomEnumFromName(name);
break;
case kNORM:
case kMODENORM:
case kLIKEWEIGHT:
case kSPLINEPARAMETER:
case kNEWSPLINE:
genenum = NUISANCEEnumFromName(name, type);
break;
default:
genenum = kNoTypeFound;
break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == kGeneratorNotBuilt) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
throw;
}
// If no type enabled
if (genenum == kNoTypeFound) {
ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl;
throw;
}
// If Not Found
if (genenum == kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
throw;
}
}
// Add offset if no issue
int nuisenum = genenum;
if (genenum != kGeneratorNotBuilt and
genenum != kNoTypeFound and
genenum != kNoDialFound) {
nuisenum += offset;
}
// Now register dial
// std::cout << "Returning " << nuisenum << std::endl;
Reweight::DialList().RegisterDialEnum(name, type, nuisenum);
return nuisenum;
}
std::string Reweight::ConvDial(int nuisenum) {
//GlobalDialList* temp;
for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) {
if (Reweight::DialList().fAllDialEnums[i] == nuisenum) {
return Reweight::DialList().fAllDialNames[i];
}
}
LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl;
return "";
}
diff --git a/src/Utils/NuisKey.cxx b/src/Utils/NuisKey.cxx
index a5dd1ad..9a0a512 100644
--- a/src/Utils/NuisKey.cxx
+++ b/src/Utils/NuisKey.cxx
@@ -1,198 +1,198 @@
#include "NuisKey.h"
std::string nuiskey::GetS(std::string name){
return Config::Get().GetS(fNode,name);
};
int nuiskey::GetI(std::string name){
return Config::Get().GetI(fNode,name);
};
double nuiskey::GetD(std::string name){
return Config::Get().GetD(fNode,name);
};
bool nuiskey::GetB(std::string name){
return Config::Get().GetB(fNode,name);
};
std::vector<std::string> nuiskey::GetVS(std::string name, const char* del){
return Config::Get().GetVS(fNode,name,del);
};
std::vector<int> nuiskey::GetVI(std::string name, const char* del){
return Config::Get().GetVI(fNode,name,del);
};
std::vector<double> nuiskey::GetVD(std::string name, const char* del){
return Config::Get().GetVD(fNode,name,del);
};
std::vector<nuiskey> Config::QueryKeys(const std::string type, const std::string test1){
// Get Vector of nodes
std::vector<XMLNodePointer_t> nodelist = Config::Get().GetNodes(type);
// Convert List into a key list for easier access
std::vector<nuiskey> keylist;
for (std::vector<XMLNodePointer_t>::const_iterator iter = nodelist.begin();
iter != nodelist.end(); iter++){
// Create new key
nuiskey newkey = nuiskey(*iter);
// Add test1
if (!test1.empty()){
std::vector<std::string> testvect = GeneralUtils::ParseToStr(test1,"=");
if (testvect.size() < 2) continue;
if (newkey.GetS(testvect[0]) != testvect[1]) continue;
}
// Save node as nuiskey
keylist.push_back( newkey );
}
// Return list of keys
return keylist;
}
nuiskey Config::QueryLastKey(const std::string type, const std::string test1){
// Loop over all for now because I'm lazy...
std::vector<nuiskey> allkeys = Config::QueryKeys(type,test1);
if (allkeys.size() < 1) return nuiskey();
else return allkeys[allkeys.size()-1];
}
nuiskey Config::QueryFirstKey(const std::string type, const std::string test1){
// Loop over all for now because I'm lazy...
std::vector<nuiskey> allkeys = Config::QueryKeys(type,test1);
if (allkeys.size() < 1) return nuiskey();
else return allkeys[allkeys.size()-1];
}
nuiskey Config::CreateParameterKeyFromLine(const std::string line){
nuiskey parameterkey = Config::CreateKey("parameter");
// Parse
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
// Add to key
parameterkey.AddS("type", strvct[0]);
parameterkey.AddS("name", strvct[1]);
parameterkey.AddS("nominal", strvct[2]);
if (strvct.size() == 7){
parameterkey.AddS("low",strvct[3]);
parameterkey.AddS("high",strvct[4]);
parameterkey.AddS("step",strvct[5]);
parameterkey.AddS("state",strvct[6]);
} else if (strvct.size() == 3){
parameterkey.AddS("state","FIX");
}
return parameterkey;
}
bool nuiskey::Has(const std::string name){
return Config::Get().Has(fNode, name);
}
nuiskey Config::CreatePullKeyFromLine(const std::string line){
- nuiskey pullkey = Config::CreateKey("pull");
+ nuiskey pullkey = Config::CreateKey("covar");
// Parse
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
pullkey.AddS("name", strvct[1]);
pullkey.AddS("input", strvct[2]);
pullkey.AddS("type", strvct[3]);
return pullkey;
}
nuiskey Config::CreateOldConfigKeyFromLine(const std::string line){
nuiskey configkey = Config::CreateKey("config");
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
configkey.AddS( strvct[1], strvct[2] );
return configkey;
}
nuiskey Config::CreateSampleKeyFromLine(const std::string line){
// Create new key
nuiskey samplekey = Config::CreateKey("sample");
// Parse
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
// Get elements
samplekey.AddS("name" , strvct[1]);
samplekey.AddS("input" , strvct[2]);
if (strvct.size() > 3){
samplekey.AddS("type", strvct[3]);
}
if (strvct.size() > 4){
samplekey.AddS("norm", strvct[4]);
}
return samplekey;
}
nuiskey Config::CreateKey(const std::string name){
return nuiskey(Config::Get().CreateNode(name));
}
void nuiskey::AddS(std::string name, std::string newval){
Config::Get().AddS(fNode, name, newval);
}
void nuiskey::AddI(std::string name, int newval){
Config::Get().AddI(fNode, name, newval);
}
void nuiskey::AddD(std::string name, double newval){
Config::Get().AddD(fNode, name, newval);
}
void nuiskey::AddB(std::string name, bool newval){
Config::Get().AddB(fNode, name, newval);
}
void nuiskey::SetS(std::string name, std::string newval){
Config::Get().SetS(fNode, name, newval);
}
void nuiskey::SetI(std::string name, int newval){
Config::Get().SetI(fNode, name, newval);
}
void nuiskey::SetD(std::string name, double newval){
Config::Get().SetD(fNode, name, newval);
}
void nuiskey::SetB(std::string name, bool newval){
Config::Get().SetB(fNode, name, newval);
}
void nuiskey::ChangeS(std::string name, std::string newval){
Config::Get().ChangeS(fNode, name, newval);
}
void nuiskey::ChangeI(std::string name, int newval){
Config::Get().ChangeI(fNode, name, newval);
}
void nuiskey::ChangeD(std::string name, double newval){
Config::Get().ChangeD(fNode, name, newval);
}
void nuiskey::ChangeB(std::string name, bool newval){
Config::Get().ChangeB(fNode, name, newval);
}
std::vector<std::string> nuiskey::GetAllKeys(){
return Config::Get().GetAllKeysForNode(fNode);
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:30 PM (5 h, 35 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023813
Default Alt Text
(200 KB)
Attached To
rNUISANCEGIT nuisancegit
Event Timeline
Log In to Comment