diff --git a/src/ANL/ANL_CC1npip_Evt_1DcosmuStar_nu.cxx b/src/ANL/ANL_CC1npip_Evt_1DcosmuStar_nu.cxx
index 30ec684..7e7f5ab 100644
--- a/src/ANL/ANL_CC1npip_Evt_1DcosmuStar_nu.cxx
+++ b/src/ANL/ANL_CC1npip_Evt_1DcosmuStar_nu.cxx
@@ -1,97 +1,97 @@
// 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 "ANL_CC1npip_Evt_1DcosmuStar_nu.h"
//********************************************************************
ANL_CC1npip_Evt_1DcosmuStar_nu::ANL_CC1npip_Evt_1DcosmuStar_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "ANL CC1npip Event Rate 1DcosmuStar nu sample. \n" \
"Target: D2 \n" \
"Flux: ANL fhc numu \n" \
"Signal: CC1pi 3 Prong (SignalDef::isCC1pi3Prong) \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("ANL #nu_mu CC1n#pi^{+}");
- // fSettings.SetDescription(descrip);
- // fSettings.SetXTitle("cos(#theta*)");
- // fSettings.SetYTitle("Number of events");
- // fSettings.SetEnuRange(0.0, 1.5);
- // fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
- // fSettings.DefineAllowedTargets("D,H");
- // fSettings.DefineAllowedSpecies("numu");
- // fSettings.SetDataInput( FitPar::GetDataBase() + "/ANL/CC1pip_on_n/ANL_CC1npip_cosmuStar.csv" );
+ fSettings.SetDescription(descrip);
+ fSettings.SetXTitle("cos(#theta*)");
+ fSettings.SetYTitle("Number of events");
+ fSettings.SetEnuRange(0.0, 1.5);
+ fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
+ fSettings.DefineAllowedTargets("D,H");
+ fSettings.DefineAllowedSpecies("numu");
+ fSettings.SetDataInput( FitPar::GetDataBase() + "/ANL/CC1pip_on_n/ANL_CC1npip_cosmuStar.csv" );
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")/(fNEvents+0.)*2./1.;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetPoissonErrors();
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void ANL_CC1npip_Evt_1DcosmuStar_nu::FillEventVariables(FitEvent *event) {
if (event->NumISParticle(2112) == 0 || // Initial state particles
event->NumFSParticle(2112) == 0 || event->NumFSParticle(211) == 0 || event->NumFSParticle(13) == 0) { // Final state particles
return;
}
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pin = event->GetHMISParticle(2112)->fP;
TLorentzVector Pn = event->GetHMFSParticle(2112)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double hadMass = FitUtils::MpPi(Pn, Ppip);
double cosmuStar = -999;
// Now need to boost into center-of-mass frame
TLorentzVector CMS = Pnu + Pin;
// Boost the muon backwards
Pmu.Boost(-CMS.BoostVector());
// Boost the neutrino forwards
Pnu.Boost(CMS.BoostVector());
// ANL has a M(pi, p) < 1.4 GeV cut imposed
// Find angle in CMS frame
if (hadMass < 1400) cosmuStar = cos(FitUtils::th(Pmu, Pnu));
fXVar = cosmuStar;
return;
};
bool ANL_CC1npip_Evt_1DcosmuStar_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi3Prong(event, 14, 211, 2112, EnuMin, EnuMax);
}
diff --git a/src/ANL/ANL_CC1ppip_XSec_1DQ2_nu.cxx b/src/ANL/ANL_CC1ppip_XSec_1DQ2_nu.cxx
index 09cfc95..d73714b 100644
--- a/src/ANL/ANL_CC1ppip_XSec_1DQ2_nu.cxx
+++ b/src/ANL/ANL_CC1ppip_XSec_1DQ2_nu.cxx
@@ -1,92 +1,121 @@
// 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/>.
*******************************************************************************/
/**
* Radecky et al. Phys Rev D, 3rd series, volume 25, number 5, 1 March 1982, p 1161-1173
*/
#include "ANL_CC1ppip_XSec_1DQ2_nu.h"
//********************************************************************
ANL_CC1ppip_XSec_1DQ2_nu::ANL_CC1ppip_XSec_1DQ2_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "ANL CC1ppip XSec 1DQ2 nu sample. \n" \
"Target: D2 \n" \
"Flux: \n" \
"Signal: \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{CC#pi} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ^{2}_{CC#pi} (GeV^{2})");
fSettings.SetAllowedTypes("FIX/FREE,SHAPE/DIAG", "FIX/DIAG");
fSettings.SetEnuRange(0.5, 6.0);
fSettings.DefineAllowedTargets("D,H");
// plot information
fSettings.SetTitle("ANL #nu_mu CC1n#pi^{+}");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput( FitPar::GetDataBase()
+ "/ANL/CC1pip_on_p/ANL_CC1pip_on_p_dSigdQ2_W14_1982.txt" );
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
+ // this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/((fNEvents+0.)*TotalIntegratedFlux("width"))*16./8.;
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / ((fNEvents + 0.) * TotalIntegratedFlux("width")) * 2. / 1.;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void ANL_CC1ppip_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(211) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pp = event->GetHMFSParticle(2212)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double hadMass = FitUtils::MpPi(Pp, Ppip);
double q2CCpip = -1.0;
// I use the W < 1.4GeV cut ANL data to isolate single pion
// there is also a W < 1.6 GeV and an uncut spectrum ANL 1982
if (hadMass < 1400) q2CCpip = -1 * (Pnu - Pmu).Mag2() / 1.E6;
fXVar = q2CCpip;
return;
};
bool ANL_CC1ppip_XSec_1DQ2_nu::isSignal(FitEvent *event) {
+ // std::cout << "CC1ppip Enu " << EnuMin << " " << EnuMax << std::endl;
return SignalDef::isCC1pi3Prong(event, 14, 211, 2212, EnuMin, EnuMax);
}
+// void ANL_CC1ppip_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) {
+// 46
+// 47 if (event->NumFSParticle(2212) == 0 ||
+// 48 event->NumFSParticle(211) == 0 ||
+// 49 event->NumFSParticle(13) == 0)
+// 50 return;
+// 51
+// 52 TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+// 53 TLorentzVector Pp = event->GetHMFSParticle(2212)->fP;
+// 54 TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
+// 55 TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
+// 56
+// 57 double hadMass = FitUtils::MpPi(Pp, Ppip);
+// 58 double q2CCpip = -1.0;
+// 59
+// 60 // I use the W < 1.4GeV cut ANL data to isolate single pion
+// 61 // there is also a W < 1.6 GeV and an uncut spectrum ANL 1982
+// 62 if (hadMass < 1400) q2CCpip = -1*(Pnu-Pmu).Mag2()/1.E6;
+// 63
+// 64 fXVar = q2CCpip;
+// 65
+// 66 return;
+// 67 };
+// 68
+// 69 bool ANL_CC1ppip_XSec_1DQ2_nu::isSignal(FitEvent *event) {
+// 70 return SignalDef::isCC1pi3Prong(event, 14, 211, 2212, EnuMin, EnuMax);
+// 71 }
diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
index 35157dd..0ea14e1 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,967 +1,1035 @@
#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! " << endl;
if (fIterationTree && !name.compare(fIterationTree->GetName())) {
fIterationTree->Reset();
return;
}
fIterationTree = new TTree(name.c_str(), name.c_str());
// Setup Main Branches
fIterationTree->Branch("total_likelihood", &fLikelihood,
"total_likelihood/D");
fIterationTree->Branch("total_ndof", &fNDOF, "total_ndof/I");
// Setup Sample Arrays
int ninputs = fSamples.size() + fPulls.size();
fSampleLikes = new double[ninputs];
fSampleNDOF = new int[ninputs];
fNDOF = GetNDOF();
// Setup Sample Branches
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
std::string name = exp->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
std::string name = pull->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
// Add Dial Branches
std::vector<std::string> dials = rw->GetDialNames();
fNDials = dials.size();
fDialVals = new double[fNDials];
for (int i = 0; i < fNDials; i++) {
fIterationTree->Branch(dials[i].c_str(), &fDialVals[i],
(dials[i] + "/D").c_str());
}
}
//***************************************************
void JointFCN::DestroyIterationTree() {
//***************************************************
if (!fIterationTree) {
delete fIterationTree;
}
}
//***************************************************
void JointFCN::WriteIterationTree() {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Can't save empty iteration tree!" << endl;
throw;
}
fIterationTree->Write();
}
//***************************************************
void JointFCN::FillIterationTree(FitWeight* rw) {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << endl;
throw;
}
rw->GetAllDials(fDialVals, fNDials);
fIterationTree->Fill();
}
//***************************************************
double JointFCN::DoEval(const double* x) {
//***************************************************
// WEIGHT ENGINE
fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
FitBase::GetRW()->UpdateWeightEngine(x);
if (fDialChanged) {
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
}
if (LOG_LEVEL(REC)) {
FitBase::GetRW()->Print();
}
// SORT SAMPLES
ReconfigureSamples();
// GET TEST STAT
fLikelihood = GetLikelihood();
// PRINT PROGRESS
LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood
<< std::endl;
// UPDATE TREE
if (fIterationTree) FillIterationTree(FitBase::GetRW());
return fLikelihood;
}
//***************************************************
int JointFCN::GetNDOF() {
//***************************************************
int totaldof = 0;
int count = 0;
// Total number of Free bins in each MC prediction
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
int dof = exp->GetNDOF();
// Save Seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double dof = pull->GetLikelihood();
// Save seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Set Data Variable
fNDOF = totaldof;
return totaldof;
}
//***************************************************
double JointFCN::GetLikelihood() {
//***************************************************
LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..." << " : " << "-2logL" << endl;
// Loop and add up likelihoods in an uncorrelated way
double like = 0.0;
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
double newlike = exp->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : " << newlike << endl;
// Add Weight Scaling
// like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName());
// Add to total
like += newlike;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double newlike = pull->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
// Add to total
like += newlike;
count++;
}
// Set Data Variable
fLikelihood = like;
return like;
};
void JointFCN::LoadSamples(std::vector<nuiskey> samplekeys) {
LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl;
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys[i];
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype = key.GetS("type");
std::string fakeData = "";
LOG(MIN) << "Loading Sample : " << samplename << std::endl;
fOutputDir->cd();
MeasurementBase* NewLoadedSample
= SampleUtils::CreateSample(key);
if (!NewLoadedSample) {
ERR(FTL) << "Could not load sample provided: " << samplename << std::endl;
ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
<< endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
}
+void JointFCN::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();
+ 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"
// << endl;
// throw;
// } else {
// fSamples.push_back(NewLoadedSample);
// }
// }
// // Sample Inputs
// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") ||
// !samplevect[0].compare("throws")) {
// // Get all inputs
// std::string name = samplevect[1];
// std::string files = samplevect[2];
// std::string type = "DEFAULT";
// if (samplevect.size() > 3) {
// type = samplevect[3];
// } else if (!samplevect[0].compare("pull")) {
// type = "GAUSPULL";
// } else if (!samplevect[0].compare("throws")) {
// type = "GAUSTHROWS";
// }
// // Create Pull Class
// LOG(MIN) << "Loading up pull term: " << name << " << " << files << " ("
// << type << ")" << std::endl;
// std::string fakeData = "";
// fOutputDir->cd();
// fPulls.push_back(new ParamPull(name, files, type));
// }
// }
// card.close();
// };
//***************************************************
void JointFCN::ReconfigureSamples(bool fullconfig) {
//***************************************************
int starttime = time(NULL);
LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << endl;
-
+ // 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" << 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) {
+ std::cout << " Save Signal = " << savesignal << std::endl;
+ sleep(5);
+ if (savesignal) {
// Reset all of our event signal vectors
fSignalEventBoxes.clear();
fSignalEventFlags.clear();
fSampleSignalFlags.clear();
fSignalEventSplines.clear();
}
// Make sure we have a list of inputs
if (fInputList.empty()) {
fInputList = GetInputList();
fSubSampleList = GetSubSampleList();
}
// If all inputs are splines make sure the readers are told
// they need to be reconfigured.
std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell reader in each BaseEvent it needs a Reconfigure next weight calc.
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) {
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
- }
+ }
// MAIN INPUT LOOP ====================
int fillcount = 0;
int inputcount = 0;
inp_iter = fInputList.begin();
// Loop over each input in manager
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Get event information
FitEvent* curevent = curinput->FirstNuisanceEvent();
int i = 0;
int nevents = curinput->GetNEvents();
int countwidth = nevents / 20;
// Start event loop iterating until we get a NULL pointer.
while (curevent) {
// Get Event Weight
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
double rwweight = curevent->Weight;
// std::cout << "RWWeight = " << curevent->RWWeight << " " << curevent->InputWeight << std::endl;
-
+
// Logging
if (LOG_LEVEL(REC)) {
if (i % countwidth == 0) {
LOG(REC) << "Processed " << i << " events. [M, W] = [" << curevent->Mode << ", " << rwweight << "]" << std::endl;
}
}
// Setup flag for if signal found in at least one sample
bool foundsignal = false;
// Create a new signal bitset for this event
std::vector<bool> signalbitset(fSubSampleList.size());
// Create a new signal box vector for this event
std::vector<MeasurementVariableBox*> signalboxes;
// Start measurement iterator
size_t measitercount = 0;
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
// Loop over all subsamples (sub in JointMeas)
for (; meas_iter != fSubSampleList.end(); meas_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// Compare input pointers, to current input, skip if not.
// Pointer tells us if it matches without doing ID checks.
if (curinput != curmeas->GetInput()) {
if (savesignal) {
// Set bit to 0 as definitely not signal
signalbitset[measitercount] = 0;
}
// Count up what measurement we are on.
measitercount++;
// Skip sample as input not signal.
continue;
}
// Fill events for matching inputs.
MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);
bool signal = curmeas->isSignal(curevent);
curmeas->SetSignal(signal);
curmeas->FillHistograms(rwweight);
// If its Signal tally up fills
if (signal) {
fillcount++;
}
// If we are saving signal/splines fill the bitset
if (savesignal) {
signalbitset[measitercount] = signal;
}
// If signal save a clone of the event box for use later.
if (savesignal and signal) {
foundsignal = true;
signalboxes.push_back( box->CloneSignalBox() );
}
// Keep track of Measurement we are on.
measitercount++;
}
// Once we've filled the measurements, if saving signal
// push back if any sample flagged this event as signal
if (savesignal) {
fSignalEventFlags.push_back(foundsignal);
}
// Save the vector of signal boxes for this event
if (savesignal and foundsignal) {
fSignalEventBoxes.push_back(signalboxes);
fSampleSignalFlags.push_back(signalbitset);
}
// If all inputs are splines we can save the spline coefficients
// for fast in memory reconfigures later.
if (fIsAllSplines 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;
+ // std::cout << "Pushed Back Coeff " << splinecount << " : ";
+ // for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) {
+ // std::cout << " " << fSignalEventSplines[splinecount][l];
+ // }
+ // std::cout << std::endl;
// }
}
// Clean up vectors once done with this event
signalboxes.clear();
signalbitset.clear();
// Iterate to the next event.
curevent = curinput->NextNuisanceEvent();
i++;
}
// Keep track of what input we are on.
inputcount++;
}
// End of Event Loop ===============================
// Now event loop is finished loop over all Measurements
// Converting Binned events to XSec Distributions
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Print out statements on approximate memory usage for profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
if (savesignal) {
int mem = ( //sizeof(fSignalEventBoxes) +
// fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
sizeof(MeasurementVariableBox1D) * fillcount) * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " signal boxes for faster access. (~" << mem << " MB)" << std::endl;
if (fIsAllSplines and !fSignalEventSplines.empty()) {
int splmem = sizeof(float) * fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size() << " spline sets into memory. (~" << splmem << " MB)" << std::endl;
}
}
LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl;
// End of reconfigure
return;
};
//***************************************************
void JointFCN::ReconfigureFastUsingManager() {
//***************************************************
// Get Start time for profilling
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// Check for saved variables if not do a full reconfigure.
if (fSignalEventFlags.empty()) {
+ 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 / 5;
// 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);
}
}
- // Start of Fast Event Loop ============================
+ // Loop over all possible spline inputs
+ int coreeventcount = 0;
+ double* coreeventweights;
+ splinecount = 0;
+
+ if (fIsAllSplines) {
+ coreeventweights = new double[fSignalEventSplines.size()];
+ splinecount = 0;
+ }
- // Start input iterators
inp_iter = fInputList.begin();
+ inpsig_iter = fSignalEventFlags.begin();
+ spline_iter = fSignalEventSplines.begin();
- // Loop over all inputs
- for (; inp_iter != fInputList.end(); inp_iter++) {
- InputHandlerBase* curinput = (*inp_iter);
- // Get Only base event with RW info.
+ // 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 (int iinput = 0; iinput < fInputList.size(); iinput++) {
+
+ InputHandlerBase* curinput = fInputList[iinput];
BaseFitEvt* curevent = curinput->FirstBaseEvent();
- int i = 0;
- // Iterate over all events and signal flags.
- while (curevent != 0 and (inpsig_iter != fSignalEventFlags.end())) {
+ for (int i = 0; i < curinput->GetNEvents(); i++) {
- // Logging
- if (LOG_LEVEL(REC)) {
- if (i % countwidth == 0) {
- LOG(REC) << "Processed " << i << " signal events.";
+ 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];
}
- }
- // If event has not been flagged as signal in vector skip it.
- if (!(*inpsig_iter)) {
+ curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
+ curevent->Weight = curevent->RWWeight * curevent->InputWeight;
+ rwweight = curevent->Weight;
- if (LOG_LEVEL(REC)) {
- if (i % countwidth == 0) {
- std::cout << std::endl;
- }
+ // #pragma omp atomic
+ coreeventweights[splinecount] = rwweight;
+
+ if (splinecount % countwidth == 0) {
+ LOG(REC) << "Processed " << splinecount << " event weights." << std::endl;
}
- inpsig_iter++;
- i++;
- continue;
+ // #pragma omp atomic
+ splinecount++;
}
- // Get The Base Event
- if (fFillNuisanceEvent) curinput->GetNuisanceEvent(i);
- else curevent = curinput->GetBaseEvent(i);
+ // #pragma omp atomic
+ sigcount++;
- // End iterator if NULL pointer at this entry..
- if (!curevent) break;
+ }
+ }
+ LOG(SAM) << "Processed event weights." << std::endl;
- // Setup signal splines pointer for this event if required.
- if (fIsAllSplines) {
- curevent->fSplineCoeff = &(*spline_iter)[0];
- }
- // Get Event Weight
- curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
- curevent->Weight = curevent->RWWeight * curevent->InputWeight;
- double rwweight = curevent->Weight;
+ // #pragma omp barrier
- // Weight Logging
- if (LOG_LEVEL(REC)) {
- if (i % countwidth == 0) {
- std::cout << " W = " << rwweight << std::endl;
- }
- }
+ // Reset Iterators
+ inpsig_iter = fSignalEventFlags.begin();
+ spline_iter = fSignalEventSplines.begin();
+ box_iter = fSignalEventBoxes.begin();
+ samsig_iter = fSampleSignalFlags.begin();
+ int nsplineweights = splinecount;
+ splinecount = 0;
- // Get iterators for this event
- std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
- std::vector<bool>::iterator subsamsig_iter = (*samsig_iter).begin();
- std::vector<MeasurementVariableBox*>::iterator subbox_iter = (*box_iter).begin();
-
- // Loop over all sub measurements.
- for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
- MeasurementBase* curmeas = (*meas_iter);
- // If event flagged as signal for this sample fill from the box.
- if (*subsamsig_iter) {
- curmeas->SetSignal(true);
- curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);
+ // Start of Fast Event Loop ============================
- // Move onto next box if there is one.
- subbox_iter++;
- fillcount++;
- }
+ // 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 (splinecount % 1000 == 0) {
- // std::cout << "Read Back Coeff " << splinecount << " : ";
- // for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) {
- // std::cout << " " << fSignalEventSplines[splinecount][l];
- // }
- // std::cout << std::endl;
- // }
-
- // Iterate over the main signal event containers.
- samsig_iter++;
- box_iter++;
- spline_iter++;
- splinecount++;
- // iterate to next signal event
- inpsig_iter++;
- i++;
+ 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..." << endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->Write();
}
// Save Pull Terms
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Write();
}
if (FitPar::Config().GetParB("EventManager")) {
// Get list of inputs
std::map<int, InputHandlerBase*> fInputs = FitBase::EvtManager().GetInputs();
std::map<int, InputHandlerBase*>::const_iterator iterInp;
for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) {
InputHandlerBase* input = (iterInp->second);
input->GetFluxHistogram()->Write();
input->GetXSecHistogram()->Write();
input->GetEventHistogram()->Write();
}
}
};
//***************************************************
void JointFCN::SetFakeData(std::string fakeinput) {
//***************************************************
LOG(MIN) << "Setting fake data from " << fakeinput << endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->SetFakeDataValues(fakeinput);
}
return;
}
diff --git a/src/FCN/JointFCN.h b/src/FCN/JointFCN.h
index 93a7e74..031fae4 100755
--- a/src/FCN/JointFCN.h
+++ b/src/FCN/JointFCN.h
@@ -1,161 +1,162 @@
#ifndef _JOINT_FCN_H_
#define _JOINT_FCN_H_
/*!
* \addtogroup FCN
* @{
*/
#include <iostream>
#include <vector>
#include <fstream>
#include <list>
// ROOT headers
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include <TMatrixDSym.h>
#include "TGraphErrors.h"
#include <TVectorD.h>
#include <TMath.h>
#include "FitUtils.h"
// External fitter headers
#include "MeasurementBase.h"
#include "SampleList.h"
#define GetCurrentDir getcwd
#include "EventManager.h"
#include "Math/Functor.h"
#include "ParamPull.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#include "MeasurementVariableBox.h"
#include "MeasurementVariableBox1D.h"
using namespace FitUtils;
using namespace FitBase;
//! Main FCN Class which ROOT's joint function needs to evaulate the chi2 at each stage of the fit.
class JointFCN
{
public:
//! Constructor
//! cardfile = Path to input card file listing samples
JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile=NULL);
JointFCN(TFile* outfile=NULL); // Loads from global config
//! Destructor
~JointFCN();
//! Create sample list from cardfile
void LoadSamples(std::vector<nuiskey> samplekeys);
-
+ void LoadPulls(std::vector<nuiskey> pullkeys);
+
//! Main Likelihood evaluation FCN
double DoEval(const double *x);
//! Func Wrapper for ROOT
inline double operator() (const std::vector<double> & x) {
double* x_array = new double[x.size()];
return this->DoEval(x_array);
};
//! Func Wrapper for ROOT
inline double operator() (const double *x) {
return this->DoEval(x);
};
//! Create a TTree to save all dial value iterations for this FCN
void CreateIterationTree(std::string name, FitWeight* rw);
//! Fills dial values and sample likelihoods into tree
void FillIterationTree(FitWeight* rw);
//! Writes TTree to fOutput directory
void WriteIterationTree();
//! Deletes TTree
void DestroyIterationTree();
//! Get Degrees of Freedom for samples (NBins)
int GetNDOF();
//! Return NDOF wrapper
inline unsigned int NDim() {return this->GetNDOF();};
//! Reconfigure samples where we force all events to be looped over.
void ReconfigureAllEvents() ;
//! Call Reconfigure on samples.
//! Choice of reconfigure type depends on whether dials have changed
//! and the MC has been filled.
void ReconfigureSamples(bool fullconfig = false);
//! Call reconfigure on only signal events (defaults to all events if CurIter=0)
void ReconfigureSignal();
//! Gets likelihood for all samples in FCN (assuming uncorrelated)
double GetLikelihood();
//! Returns list of pointers to the samples
inline std::list<MeasurementBase*> GetSampleList(){ return fSamples; }
//! Return list of pointers to all the pulls
inline std::list<ParamPull*> GetPullList(){ return fPulls; };
//! Write all samples to output DIR
void Write();
//! Set Fake data from file/MC
void SetFakeData(std::string fakeinput);
//! Reconfigure looping over duplicate inputs
void ReconfigureUsingManager();
//! Reconfigure Fast looping over duplicate inputs
void ReconfigureFastUsingManager();
std::vector<MeasurementBase*> GetSubSampleList();
std::vector<InputHandlerBase*> GetInputList();
private:
//! Append the experiments to include in the fit to this list
std::list<MeasurementBase*> fSamples;
//! Append parameter pull terms to include penalties in the fit to this list
std::list<ParamPull*> fPulls;
TDirectory *fOutputDir; //!< Directory to save contents
std::string fCardFile; //!< Input Card text file
bool fDialChanged; //!< Flag for if RW dials changed
UInt_t fCurIter; //!< Counter for how many times reconfigure called
bool fMCFilled; //!< Check MC has at least been filled once
TTree* fIterationTree; //!< Tree to save RW values on each function call
int fNDials; //!< Number of RW Dials in FitWeight
double* fDialVals; //!< Current Values of RW Dials
double fLikelihood; //!< Current likelihood for joint sample likelihood
double fNDOF; //!< Total NDOF
double* fSampleLikes; //!< Likelihoods for each individual measurement in list
int * fSampleNDOF; //!< NDOF for each individual measurement in list
bool fUsingEventManager; //!< Flag for doing joint comparisons
std::vector< std::vector<float> > fSignalEventSplines;
std::vector< std::vector<MeasurementVariableBox*> > fSignalEventBoxes;
std::vector< bool > fSignalEventFlags;
std::vector< std::vector<bool> > fSampleSignalFlags;
std::vector<InputHandlerBase*> fInputList;
std::vector<MeasurementBase*> fSubSampleList;
bool fIsAllSplines;
};
/*! @} */
#endif // _JOINT_FCN_H_
diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx
index 946b423..47da54a 100644
--- a/src/FCN/SampleList.cxx
+++ b/src/FCN/SampleList.cxx
@@ -1,631 +1,633 @@
#include "SampleList.h"
//! Functions to make it easier for samples to be created and handled.
namespace SampleUtils {
//! Create a given sample given its name, file, type, fakdata(fkdt) file and the
//! current rw engine and push it back into the list fChain.
MeasurementBase* CreateSample(std::string name, std::string file,
std::string type, std::string fkdt, FitWeight* rw) {
nuiskey samplekey = Config::CreateKey("sample");
samplekey.AddS("name", name);
samplekey.AddS("input", file);
samplekey.AddS("type", type);
return CreateSample(samplekey);
}
MeasurementBase* CreateSample(nuiskey samplekey) {
FitWeight* rw = FitBase::GetRW();
std::string name = samplekey.GetS("name");
std::string file = samplekey.GetS("input");
std::string type = samplekey.GetS("type");
std::string fkdt = "";
/*
ANL CCQE Samples
*/
if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) {
return (new ANL_CCQE_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) {
return (new ANL_CCQE_Evt_1DQ2_nu(samplekey));
/*
ANL CC1ppip samples
*/
} else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) {
return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) {
return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) {
return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) {
return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) {
return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey));
/*
ANL CC1npip sample
*/
} else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) {
return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey));
/*
ANL CC1pi0 sample
*/
} else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey));
/*
ANL NC1npip sample
*/
} else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) {
return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey));
/*
ANL NC1ppim sample
*/
} else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) {
return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) {
return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey));
/*
ANL CC2pi sample
*/
} else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey));
/*
ArgoNeut Samples
*/
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) {
return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) {
return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) {
return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) {
return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey));
/*
BNL Samples
*/
} else if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) {
return (new BNL_CCQE_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) {
return (new BNL_CCQE_Evt_1DQ2_nu(samplekey));
/*
BNL CC1ppip samples
*/
} else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu") ||
- !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr")) {
+ !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr") ||
+ !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut") ||
+ !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr")) {
return (new BNL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu") ||
!name.compare("BNL_CC1ppip_Evt_1DQ2_nu_W14Cut")) {
return (new BNL_CC1ppip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) {
return (new BNL_CC1ppip_Evt_1DcosthAdler_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) {
return (new BNL_CC1ppip_Evt_1Dphi_nu(samplekey));
/*
BNL CC1npip samples
*/
} else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu") ||
!name.compare("BNL_CC1npip_XSec_1DEnu_nu_Uncorr")) {
return (new BNL_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) {
return (new BNL_CC1npip_Evt_1DQ2_nu(samplekey));
/*
BNL CC1pi0 samples
*/
} else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) {
return (new BNL_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) {
return (new BNL_CC1pi0_Evt_1DQ2_nu(samplekey));
/*
FNAL Samples
*/
} else if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) {
return (new FNAL_CCQE_Evt_1DQ2_nu(samplekey));
/*
FNAL CC1ppip
*/
} else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) {
return (new FNAL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) {
return (new FNAL_CC1ppip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) {
return (new FNAL_CC1ppip_Evt_1DQ2_nu(samplekey));
/*
FNAL CC1ppim
*/
} else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) {
return (new FNAL_CC1ppim_XSec_1DEnu_antinu(samplekey));
/*
BEBC Samples
*/
} else if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) {
return (new BEBC_CCQE_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1ppip samples
*/
} else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) {
return (new BEBC_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) {
return (new BEBC_CC1ppip_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1npip samples
*/
} else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) {
return (new BEBC_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) {
return (new BEBC_CC1npip_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1pi0 samples
*/
} else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) {
return (new BEBC_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) {
return (new BEBC_CC1pi0_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1npim samples
*/
} else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) {
return (new BEBC_CC1npim_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) {
return (new BEBC_CC1npim_XSec_1DQ2_antinu(samplekey));
/*
BEBC CC1ppim samples
*/
} else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) {
return (new BEBC_CC1ppim_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) {
return (new BEBC_CC1ppim_XSec_1DQ2_antinu(samplekey));
/*
GGM CC1ppip samples
*/
} else if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) {
return (new GGM_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) {
return (new GGM_CC1ppip_Evt_1DQ2_nu(samplekey));
/*
MiniBooNE Samples
*/
/*
CCQE
*/
} else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) {
return (new MiniBooNE_CCQE_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu") ||
!name.compare("MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu")) {
return (new MiniBooNE_CCQE_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) {
return (new MiniBooNE_CCQE_XSec_2DTcos_nu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) {
return (new MiniBooNE_CCQE_XSec_2DTcos_antinu(samplekey));
/*
MiniBooNE CC1pi+
*/
// 1D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DTpi_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DTu_nu(samplekey));
// 2D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(samplekey));
/*
MiniBooNE CC1pi0
*/
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DTu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_nu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_antinu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_rhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dppi0_antinu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_nu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_fhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dppi0_nu(samplekey));
/*
MiniBooNE NCEL
*/
} else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) {
ERR(FTL)
<< "MiniBooNE_NCEL_XSec_Treco_nu not implemented in current interface."
<< std::endl;
throw 5;
// return (new MiniBooNE_NCEL_XSec_Treco_nu(file, rw, type,
// fkdt));
/*
MINERvA Samples
*/
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) {
return (
new MINERvA_CCQE_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) {
return (
new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) {
return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) {
return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) {
return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) {
return (
new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) {
return (
new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey));
/*
CC1pi+
*/
// DONE
} else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_fluxcorr") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg_fluxcorr") ) {
return (new MINERvA_CC1pip_XSec_1DTpi_nu(samplekey));
// DONE
} else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_fluxcorr") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg_fluxcorr")) {
return (new MINERvA_CC1pip_XSec_1Dth_nu(samplekey));
/*
CCNpi+
*/
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2016") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_fluxcorr") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg_fluxcorr")) {
return (new MINERvA_CCNpip_XSec_1Dth_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2016") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_fluxcorr") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg_fluxcorr")) {
return (new MINERvA_CCNpip_XSec_1DTpi_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) {
return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) {
return (new MINERvA_CCNpip_XSec_1DEnu_nu(samplekey));
/*
CC1pi0
*/
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_fluxcorr") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015_fluxcorr") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016_fluxcorr")) {
return (new MINERvA_CC1pi0_XSec_1Dth_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu_fluxcorr")) {
return (new MINERvA_CC1pi0_XSec_1Dppi0_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi0_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DTpi0_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DQ2_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1Dthmu_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1Dpmu_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DEnu_antinu(samplekey));
/*
CCINC
*/
} else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) {
return (new MINERvA_CCinc_XSec_2DEavq3_nu(samplekey));
} else if (!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Pb208_CH")) {
return (
new MINERvA_CCinc_XSec_1Dx_ratio(samplekey));
} else if (!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Pb208_CH")) {
return (new MINERvA_CCinc_XSec_1DEnu_ratio(samplekey));
/*
T2K Samples
*/
} else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu") ||
!name.compare("T2K_CC0pi_XSec_2DPcos_nu_I") ||
!name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) {
return (new T2K_CC0pi_XSec_2DPcos_nu(samplekey));
/*
T2K CC1pi+ CH samples
*/
// Comment these out for now because we don't have the proper data
/*
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dpmu_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dpmu_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dppi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) {
return (new T2K_CC1pip_CH_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthq3pi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) {
return (new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt));
*/
/*
T2K CC1pi+ H2O samples
*/
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu")) {
return (new T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuMB_nu")) {
return (new T2K_CC1pip_H2O_XSec_1DEnuMB_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmu_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcosmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcospi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcospi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dpmu_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dppi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dppi_nu(samplekey));
/*
T2K CC0pi + np CH samples
*/
} else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) {
return (new T2K_CC0pinp_STV_XSec_1Ddpt_nu(samplekey));
/*
K2K Samples
*/
/*
NC1pi0
*/
} else if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) {
return (new K2K_NC1pi0_Evt_1Dppi0_nu(samplekey));
/*
Fake Studies
*/
} else if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
return (
new ExpMultDist_CCQE_XSec_1DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("ExpMultDist_CCQE_XSec_2D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
return (
new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("GenericFlux_") != std::string::npos) {
return (new GenericFlux_Tester(name, file, rw, type, fkdt));
} else if (!name.compare("T2K2017_FakeData")) {
return (new T2K2017_FakeData(samplekey));
} else if (!name.compare("ElectronFlux_FlatTree")) {
return (new ElectronFlux_FlatTree(name, file, rw, type, fkdt));
} else if (name.find("ElectronData_") != std::string::npos) {
return new ElectronScattering_DurhamData(samplekey);
//<<<<<<< HEAD
// } else if (name.find("MCStudy_KaonPreSelection") != std::string::npos) {
// return (new MCStudy_KaonPreSelection(name, file, rw, type, fkdt));
//=======
//} else if (name.find("MCStudy_KaonPreSelection") != std::string::npos) {
//fChain->push_back(new MCStudy_KaonPreSelection(name, file, rw, type, fkdt));
//>>>>>>> 96ed014ac03821c4f771d6c484740e8b25350aa1
} else if (name.find("MuonValidation_") != std::string::npos) {
return (new MCStudy_MuonValidation(name, file, rw, type, fkdt));
} else {
ERR(FTL) << "Error: No such sample: " << name << std::endl;
exit(-1);
return NULL;
}
// Return NULL if no sample loaded.
return NULL;
}
}
diff --git a/src/FitBase/InputHandler2.h b/src/FitBase/InputHandler2.h
index f5877b2..2de4759 100644
--- a/src/FitBase/InputHandler2.h
+++ b/src/FitBase/InputHandler2.h
@@ -1,183 +1,189 @@
#ifndef INPUTHANDLER2_H
#define INPUTHANDLER2_H
#include "TH1D.h"
#include "FitEvent.h"
#include "BaseFitEvt.h"
class InputHandlerBase {
public:
InputHandlerBase() {
fName = "";
fFluxHist = NULL;
fEventHist = NULL;
fNEvents = 0;
fNUISANCEEvent = NULL;
fBaseEvent = NULL;
};
~InputHandlerBase() {};
virtual FitEvent* GetNuisanceEvent(const UInt_t entry) = 0;
virtual BaseFitEvt* GetBaseEvent(const UInt_t entry) = 0;
virtual void Print(){};
inline int GetNEvents(void) { return fNEvents; }
inline TH1D* GetFluxHistogram(void) {return fFluxHist;};
inline TH1D* GetEventHistogram(void) {return fEventHist;};
inline std::string GetName(void) {return fName;};
inline int GetType(void) {return fEventType;};
inline TH1D* GetXSecHistogram(void) {
fXSecHist = (TH1D*)fFluxHist->Clone();
fXSecHist->Divide(fEventHist);
return fXSecHist;
};
//********************************************************************
inline double PredictedEventRate(double low, double high,
std::string intOpt) {
//********************************************************************
int minBin = fFluxHist->GetXaxis()->FindBin(low);
int maxBin = fFluxHist->GetXaxis()->FindBin(high);
return fEventHist->Integral(minBin, maxBin + 1, intOpt.c_str());
};
//********************************************************************
- inline double TotalIntegratedFlux(double low, double high,
- std::string intOpt) {
+ inline double TotalIntegratedFlux(double low = -9999.9, double high = -9999.9,
+ std::string intOpt="") {
//********************************************************************
+ std::cout << "Getting Total Integrated Flux Between : " << low << " - " << high << std::endl;
+
Int_t minBin = fFluxHist->GetXaxis()->FindFixBin(low);
Int_t maxBin = fFluxHist->GetXaxis()->FindFixBin(high);
if ((fFluxHist->IsBinOverflow(minBin) && (low != -9999.9))) {
minBin = 1;
}
if ((fFluxHist->IsBinOverflow(maxBin) && (high != -9999.9))) {
maxBin = fFluxHist->GetXaxis()->GetNbins() + 1;
}
-
+ std::cout << "Getting Between Bins " << minBin << " " << maxBin << std::endl;
// If we are within a single bin
if (minBin == maxBin) {
+ std::cout << "Getting minBin == maxBin " << std::endl;
// Get the contained fraction of the single bin's width
return ((high - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
}
double lowBinUpEdge = fFluxHist->GetXaxis()->GetBinUpEdge(minBin);
double highBinLowEdge = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin);
double lowBinfracIntegral =
((lowBinUpEdge - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
double highBinfracIntegral =
((high - highBinLowEdge) / fFluxHist->GetXaxis()->GetBinWidth(maxBin)) *
fFluxHist->Integral(maxBin, maxBin, intOpt.c_str());
// If they are neighbouring bins
if ((minBin + 1) == maxBin) {
+ std::cout << "Get lowfrac + highfrac" << std::endl;
// Get the contained fraction of the two bin's width
return lowBinfracIntegral + highBinfracIntegral;
}
+ std::cout << "Returning highBinFracIntegral and LowBinFracIntegral " << std::endl;
// If there are filled bins between them
return lowBinfracIntegral + highBinfracIntegral +
fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
+ // return fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
}
std::vector<TH1*> GetFluxList(void) { return std::vector<TH1*>(1, fFluxHist); };
std::vector<TH1*> GetEventList(void) { return std::vector<TH1*>(1, fEventHist); };
std::vector<TH1*> GetXSecList(void) { return std::vector<TH1*>(1, GetXSecHistogram()); };
virtual FitEvent* FirstNuisanceEvent() {
fCurrentIndex = 0;
return GetNuisanceEvent(fCurrentIndex);
};
virtual FitEvent* NextNuisanceEvent() {
fCurrentIndex++;
if (jointinput and fMaxEvents != -1) {
while ( fCurrentIndex < jointindexlow[jointindexswitch] ||
fCurrentIndex >= jointindexhigh[jointindexswitch] ) {
jointindexswitch++;
// Loop Around
if (jointindexswitch == jointindexlow.size()) {
jointindexswitch = 0;
}
}
if (fCurrentIndex > jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) {
fCurrentIndex = jointindexlow[jointindexswitch];
}
}
return GetNuisanceEvent(fCurrentIndex);
};
virtual BaseFitEvt* FirstBaseEvent() {
fCurrentIndex = 0;
return GetBaseEvent(fCurrentIndex);
};
virtual BaseFitEvt* NextBaseEvent() {
fCurrentIndex++;
if (jointinput and fMaxEvents != -1) {
while ( fCurrentIndex < jointindexlow[jointindexswitch] ||
fCurrentIndex >= jointindexhigh[jointindexswitch] ) {
jointindexswitch++;
// Loop Around
if (jointindexswitch == jointindexlow.size()) {
jointindexswitch = 0;
}
}
if (fCurrentIndex > jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) {
fCurrentIndex = jointindexlow[jointindexswitch];
}
}
return GetBaseEvent(fCurrentIndex);
};
std::vector<TH1D*> jointfluxinputs;
std::vector<TH1D*> jointeventinputs;
std::vector<int> jointindexlow;
std::vector<int> jointindexhigh;
std::vector<int> jointindexallowed;
size_t jointindexswitch;
bool jointinput;
std::vector<double> jointindexscale;
std::string fName;
TH1D* fFluxHist;
TH1D* fEventHist;
TH1D* fXSecHist;
int fNEvents;
int fMaxEvents;
FitEvent* fNUISANCEEvent;
BaseFitEvt* fBaseEvent;
int fEventType;
int fCurrentIndex;
};
#endif
\ No newline at end of file
diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx
index 60cfce2..3e2ab60 100644
--- a/src/FitBase/Measurement1D.cxx
+++ b/src/FitBase/Measurement1D.cxx
@@ -1,1778 +1,1782 @@
// 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";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = 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) 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!" << endl;
ERR(WRN) << "If you want to use it please add fNormError=VAL" << 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);
}
}
// 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);
}
}
// 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(),
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);
}
// 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("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;
std::cout << fSettings.fKeyValues.fNode << std::endl;
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);
// }
// 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() {
//********************************************************************
- return fDataHist->GetNbinsX() - fMaskHist->Integral();
+ 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;
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) {
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);
}
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 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) 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);
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/Measurement2D.cxx b/src/FitBase/Measurement2D.cxx
index e7f4840..27c9550 100644
--- a/src/FitBase/Measurement2D.cxx
+++ b/src/FitBase/Measurement2D.cxx
@@ -1,1935 +1,1937 @@
// 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 "Measurement2D.h"
#include "TDecompChol.h"
//********************************************************************
Measurement2D::Measurement2D(void) {
//********************************************************************
covar = NULL;
fDecomp = NULL;
fFullCovar = NULL;
fMCHist = NULL;
fMCFine = NULL;
fDataHist = NULL;
fMCHist_X = NULL;
fMCHist_Y = NULL;
fDataHist_X = NULL;
fDataHist_Y = NULL;
fMaskHist = NULL;
fMapHist = NULL;
fDataOrig = NULL;
fDataTrue = NULL;
fMCWeighted = NULL;
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/FITPROJX/"
"FITPROJY";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsDifXSec = false;
fIsEnu = false;
// 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";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsDifXSec = false;
fIsEnu1D = false;
// Inputs
fInput = NULL;
fRW = NULL;
// Extra Histograms
fMCHist_Modes = NULL;
}
//********************************************************************
Measurement2D::~Measurement2D(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 Measurement2D::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) 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!" << endl;
ERR(WRN) << "If you want to use it please add fNormError=VAL" << endl;
throw;
}
}
}
void Measurement2D::CreateDataHistogram(int dimx, double* binx, int dimy, double* biny) {
if (fDataHist) delete fDataHist;
LOG(SAM) << "Creating Data Histogram dim : " << dimx << " " << dimy << std::endl;
fDataHist = new TH2D( (fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(),
dimx - 1, binx, dimy - 1, biny );
}
void Measurement2D::SetDataFromTextFile(std::string datfile) {
// fDataHist = PlotUtils::GetTH2DFromTextFile(datfile,"");
}
void Measurement2D::SetDataFromRootFile(std::string datfile, std::string histname) {
fDataHist = PlotUtils::GetTH2DFromRootFile(datfile, histname);
}
void Measurement2D::SetDataValuesFromTextFile(std::string datfile, TH2D* hist) {
LOG(SAM) << "Setting data values from text file" << std::endl;
if (!hist) hist = fDataHist;
// Read TH2D From textfile
TH2D* valhist = (TH2D*) hist->Clone();
valhist->Reset();
PlotUtils::Set2DHistFromText(datfile, valhist, 1.0, true);
LOG(SAM) << " -> Filling values from read hist." << std::endl;
for (int i = 0; i < valhist->GetNbinsX(); i++) {
for (int j = 0; j < valhist->GetNbinsY(); j++) {
hist->SetBinContent(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1));
}
}
LOG(SAM) << " --> Done" << std::endl;
}
void Measurement2D::SetDataErrorsFromTextFile(std::string datfile, TH2D* hist) {
LOG(SAM) << "Setting data errors from text file" << std::endl;
if (!hist) hist = fDataHist;
// Read TH2D From textfile
TH2D* valhist = (TH2D*) hist->Clone();
valhist->Reset();
PlotUtils::Set2DHistFromText(datfile, valhist, 1.0);
// Fill Errors
LOG(SAM) << " -> Filling errors from read hist." << std::endl;
for (int i = 0; i < valhist->GetNbinsX(); i++) {
for (int j = 0; j < valhist->GetNbinsY(); j++) {
hist->SetBinError(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1));
}
}
LOG(SAM) << " --> Done" << std::endl;
}
void Measurement2D::SetMapValuesFromText(std::string dataFile) {
TH2D* hist = fDataHist;
std::vector<double> edgex;
std::vector<double> edgey;
for (int i = 0; i <= hist->GetNbinsX(); i++) edgex.push_back(hist->GetXaxis()->GetBinLowEdge(i+1));
for (int i = 0; i <= hist->GetNbinsY(); i++) edgey.push_back(hist->GetYaxis()->GetBinLowEdge(i+1));
fMapHist = new TH2I((fName + "_map").c_str(), (fName + fPlotTitles).c_str(),
edgex.size()-1, &edgex[0], edgey.size()-1, &edgey[0]);
LOG(SAM) << "Reading map from: " << dataFile << std::endl;
PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0);
}
//********************************************************************
void Measurement2D::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 Measurement2D::SetCovarFromDiagonal(TH2D* 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 Measurement2D::SetCovarFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = this->GetNDOF(false);
}
LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl;
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement2D::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 Measurement2D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = this->GetNDOF(false);
}
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 Measurement2D::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 Measurement2D::SetCorrelationFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) dim = this->GetNDOF(false);
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);
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement2D::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);
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement2D::SetCholDecompFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = this->GetNDOF(false);
}
LOG(SAM) << "Reading cholesky from text file: " << covfile << " " << dim << 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 Measurement2D::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 Measurement2D::ScaleData(double scale) {
//********************************************************************
fDataHist->Scale(scale);
}
//********************************************************************
void Measurement2D::ScaleDataErrors(double scale) {
//********************************************************************
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsY(); j++) {
fDataHist->SetBinError(i + 1, j + 1, fDataHist->GetBinError(i + 1, j + 1) * scale);
}
}
}
//********************************************************************
void Measurement2D::ScaleCovar(double scale) {
//********************************************************************
(*fFullCovar) *= scale;
(*covar) *= 1.0 / scale;
(*fDecomp) *= sqrt(scale);
}
//********************************************************************
void Measurement2D::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 nbinsx = fDataHist->GetNbinsX();
int nbinxy = fDataHist->GetNbinsY();
fMaskHist =
new TH2I((fSettings.GetName() + "_BINMASK").c_str(),
(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbinsx, 0, nbinsx, nbinxy, 0, nbinxy);
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) << "Measurement2D::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[2] > 0) val = 1;
fMaskHist->SetBinContent(entries[0], entries[1], val);
}
// Apply masking by setting masked data bins to zero
PlotUtils::MaskBins(fDataHist, fMaskHist);
return;
}
//********************************************************************
void Measurement2D::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 = (TH2D*)fDataHist->Clone();
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist->Reset();
// Setup fMCFine
fMCFine = new TH2D("mcfine", "mcfine", fDataHist->GetNbinsX() * 6,
fMCHist->GetXaxis()->GetBinLowEdge(1),
fMCHist->GetXaxis()->GetBinLowEdge(fDataHist->GetNbinsX() + 1),
fDataHist->GetNbinsY() * 6,
fMCHist->GetYaxis()->GetBinLowEdge(1),
fMCHist->GetYaxis()->GetBinLowEdge(fDataHist->GetNbinsY() + 1));
fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCFine->Reset();
// Setup MC Stat
fMCStat = (TH2D*)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);
}
// 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 = (TH2D*)fMCHist->Clone();
fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
(fName + "_MCWGHTS" + fPlotTitles).c_str());
fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
}
}
//********************************************************************
void Measurement2D::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("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;
// 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;
else
fIsChi2 = true;
// EXTRAS
if (opt.find("RAW") != std::string::npos) fIsRawEvents = true;
if (opt.find("DIF") != std::string::npos) fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) fIsEnu = true;
if (opt.find("NORM") != std::string::npos) fAddNormPen = true;
if (opt.find("MASK") != std::string::npos) fIsMask = true;
fIsProjFitX = (opt.find("FITPROJX") != std::string::npos);
fIsProjFitY = (opt.find("FITPROJY") != std::string::npos);
return;
};
/*
Reconfigure LOOP
*/
//********************************************************************
void Measurement2D::ResetAll() {
//********************************************************************
fMCHist->Reset();
fMCFine->Reset();
fMCStat->Reset();
return;
};
//********************************************************************
void Measurement2D::FillHistograms() {
//********************************************************************
if (Signal) {
fMCHist->Fill(fXVar, fYVar, Weight);
fMCFine->Fill(fXVar, fYVar, Weight);
fMCStat->Fill(fXVar, fYVar, 1.0);
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, fYVar, Weight);
}
return;
};
//********************************************************************
void Measurement2D::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);
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor);
// if (fMCHist_Modes) {
// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
// GetEventHistogram(), fScaleFactor,
// fNEvents);
// }
// 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 Measurement2D::ApplyNormScale(double norm) {
//********************************************************************
fCurrentNorm = norm;
fMCHist->Scale(1.0 / norm);
fMCFine->Scale(1.0 / norm);
return;
};
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int Measurement2D::GetNDOF(bool applymasking) {
//********************************************************************
// return fDataHist->GetNbinsX() - fMaskHist->Integral();
// Just incase it has gone...
if (!fDataHist) return 0;
int nDOF = 0;
// If datahist has no errors make sure we don't include those bins as they are
// not data points
for (int xBin = 0; xBin < fDataHist->GetNbinsX() + 1; ++xBin) {
for (int yBin = 0; yBin < fDataHist->GetNbinsY() + 1; ++yBin) {
if (fDataHist->GetBinContent(xBin, yBin) != 0)
++nDOF;
}
}
// Account for possible bin masking
int nMasked = 0;
if (fMaskHist and fIsMask)
if (fMaskHist->Integral() > 0)
for (int xBin = 0; xBin < fMaskHist->GetNbinsX() + 1; ++xBin)
for (int yBin = 0; yBin < fMaskHist->GetNbinsY() + 1; ++yBin)
if (fMaskHist->GetBinContent(xBin, yBin) > 0.5) ++nMasked;
// Take away those masked DOF
if (applymasking) {
nDOF -= nMasked;
}
return nDOF;
}
//********************************************************************
double Measurement2D::GetLikelihood() {
//********************************************************************
// If this is for a ratio, there is no data histogram to compare to!
if (fNoData || !fDataHist) return 0.;
// Fix weird masking bug
if (!fIsMask) {
if (fMaskHist) {
fMaskHist = NULL;
}
} else {
if (fMaskHist) {
PlotUtils::MaskBins(fMCHist, fMaskHist);
}
}
// if (fIsProjFitX or fIsProjFitY) return GetProjectedChi2();
// Scale up the results to match each other (Not using width might be
// inconsistent with Meas1D)
double scaleF = fDataHist->Integral() / fMCHist->Integral();
if (fIsShape) {
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF);
}
if (!fMapHist) {
fMapHist = StatUtils::GenerateMap(fDataHist);
}
// Get the chi2 from either covar or diagonals
double chi2;
if (fIsChi2) {
if (fIsDiag) {
chi2 =
StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMapHist, fMaskHist);
} else {
chi2 = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist,
fMaskHist);
}
}
// Add a normal penalty term
if (fAddNormPen) {
chi2 +=
(1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError);
LOG(REC) << "Norm penalty = "
<< (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) /
(fNormError * fNormError)
<< std::endl;
}
// Adjust the shape back to where it was.
if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
return chi2;
}
/*
Fake Data Functions
*/
//********************************************************************
void Measurement2D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Setup original/datatrue
TH2D* tempdata = (TH2D*) fDataHist->Clone();
if (!fIsFakeData) {
fIsFakeData = true;
// Make a copy of the original data histogram.
if (!fDataOrig) fDataOrig = (TH2D*)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 = (TH2D*)fMCHist->Clone((fName + "_MC").c_str());
// Fake File
} else {
if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
fDataHist = (TH2D*)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 = (TH2D*)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() * fDataHist->GetNbinsY();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
if (tempdata->GetBinContent(i + 1) && tempdata->GetBinContent(j + 1)) {
alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
} else {
alpha_i = 0.0;
alpha_j = 0.0;
}
(*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 Measurement2D::ResetFakeData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH2D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
}
}
//********************************************************************
void Measurement2D::ResetData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH2D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
}
fIsFakeData = false;
}
//********************************************************************
void Measurement2D::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
*/
//********************************************************************
TH2D* Measurement2D::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;
};
//********************************************************************
TH2D* Measurement2D::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 Measurement2D::Write(std::string drawOpt) {
//********************************************************************
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
// 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
// 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();
/// 2D VERSION
// If null pointer return
if (!fMCHist and !fDataHist) {
LOG(SAM) << fName << "Incomplete histogram set!" << std::endl;
return;
}
// FitPar::Config().out->cd();
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
bool drawData = (drawOpt.find("DATA") != std::string::npos);
bool drawNormal = (drawOpt.find("MC") != std::string::npos);
bool drawEvents = (drawOpt.find("EVT") != std::string::npos);
bool drawXSec = (drawOpt.find("XSEC") != std::string::npos);
bool drawFine = (drawOpt.find("FINE") != std::string::npos);
bool drawRatio = (drawOpt.find("RATIO") != std::string::npos);
bool drawModes = (drawOpt.find("MODES") != std::string::npos);
bool drawShape = (drawOpt.find("SHAPE") != std::string::npos);
bool residual = (drawOpt.find("RESIDUAL") != std::string::npos);
bool drawMatrix = (drawOpt.find("MATRIX") != std::string::npos);
bool drawFlux = (drawOpt.find("FLUX") != std::string::npos);
bool drawMask = (drawOpt.find("MASK") != std::string::npos);
bool drawMap = (drawOpt.find("MAP") != std::string::npos);
bool drawProj = (drawOpt.find("PROJ") != std::string::npos);
// bool drawCanvPDG = (drawOpt.find("CANVPDG") != std::string::npos);
bool drawCov = (drawOpt.find("COV") != std::string::npos);
bool drawSliceCanvYMC = (drawOpt.find("CANVYMC") != std::string::npos);
bool drawWeighted = (drawOpt.find("WGHT") != std::string::npos);
if (FitPar::Config().GetParB("EventManager")) {
drawFlux = false;
drawXSec = false;
drawEvents = false;
}
if (fMaskHist) fMaskHist->Write();
// Save standard plots
if (drawData) this->GetDataList().at(0)->Write();
if (drawNormal) this->GetMCList().at(0)->Write();
if (drawCov) {
TH2D(*fFullCovar).Write((fName + "_COV").c_str());
}
if (drawOpt.find("INVCOV") != std::string::npos) {
TH2D(*covar).Write((fName + "_INVCOV").c_str());
}
// Generate a simple map
if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist);
// Convert to 1D Lists
TH1D* data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist);
TH1D* mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist);
TH1I* mask_1D = StatUtils::MapToMask(fMaskHist, fMapHist);
data_1D->Write();
mc_1D->Write();
if (mask_1D) {
mask_1D->Write();
TMatrixDSym* calc_cov =
StatUtils::ApplyInvertedMatrixMasking(covar, mask_1D);
TH1D* calc_data = StatUtils::ApplyHistogramMasking(data_1D, mask_1D);
TH1D* calc_mc = StatUtils::ApplyHistogramMasking(mc_1D, mask_1D);
TH2D* bin_cov = new TH2D(*calc_cov);
bin_cov->Write();
calc_data->Write();
calc_mc->Write();
delete mask_1D;
delete calc_cov;
delete calc_data;
delete calc_mc;
delete bin_cov;
}
delete data_1D;
delete mc_1D;
// Save only mc and data if splines
if (fEventType == 4 or fEventType == 3) {
return;
}
// Draw Extra plots
if (drawFine) this->GetFineList().at(0)->Write();
if (drawFlux and GetFluxHistogram()) {
GetFluxHistogram()->Write();
}
if (drawEvents and GetEventHistogram()) {
GetEventHistogram()->Write();
}
if (fIsMask and drawMask) {
fMaskHist->Write((fName + "_MSK").c_str()); //< save mask
}
if (drawMap) fMapHist->Write((fName + "_MAP").c_str()); //< save map
// // Save neut stack
// if (drawModes) {
// THStack combo_fMCHist_PDG = PlotUtils::GetNeutModeStack(
// (fName + "_MC_PDG").c_str(), (TH1**)fMCHist_PDG, 0);
// combo_fMCHist_PDG.Write();
// }
// Save Matrix plots
if (drawMatrix and fFullCovar and covar and fDecomp) {
TH2D cov = TH2D((*fFullCovar));
cov.SetNameTitle((fName + "_cov").c_str(),
(fName + "_cov;Bins; Bins;").c_str());
cov.Write();
TH2D covinv = TH2D((*this->covar));
covinv.SetNameTitle((fName + "_covinv").c_str(),
(fName + "_cov;Bins; Bins;").c_str());
covinv.Write();
TH2D covdec = TH2D((*fDecomp));
covdec.SetNameTitle((fName + "_covdec").c_str(),
(fName + "_cov;Bins; Bins;").c_str());
covdec.Write();
}
// Save ratio plots if required
if (drawRatio) {
// Needed for error bars
for (int i = 0; i < fMCHist->GetNbinsX() * fMCHist->GetNbinsY(); i++)
fMCHist->SetBinError(i + 1, 0.0);
fDataHist->GetSumw2();
fMCHist->GetSumw2();
// Create Ratio Histograms
TH2D* dataRatio = (TH2D*)fDataHist->Clone((fName + "_data_RATIO").c_str());
TH2D* mcRatio = (TH2D*)fMCHist->Clone((fName + "_MC_RATIO").c_str());
mcRatio->Divide(fMCHist);
dataRatio->Divide(fMCHist);
// Cancel bin errors on MC
for (int i = 0; i < mcRatio->GetNbinsX() * mcRatio->GetNbinsY(); i++) {
mcRatio->SetBinError(
i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
}
mcRatio->SetMinimum(0);
mcRatio->SetMaximum(2);
dataRatio->SetMinimum(0);
dataRatio->SetMaximum(2);
mcRatio->Write();
dataRatio->Write();
delete mcRatio;
delete dataRatio;
}
// Save Shape Plots if required
if (drawShape) {
// Create Shape Histogram
TH2D* mcShape = (TH2D*)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);
mcShape->SetLineWidth(3);
mcShape->SetLineStyle(7); // dashes
mcShape->Write();
// Save shape ratios
if (drawRatio) {
// Needed for error bars
mcShape->GetSumw2();
// Create shape ratio histograms
TH2D* mcShapeRatio =
(TH2D*)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
TH2D* dataShapeRatio =
(TH2D*)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); // dashes
mcShapeRatio->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
delete mcShape;
}
// Save residual calculations of what contributed to the chi2 values.
if (residual) {
}
if (fIsProjFitX or fIsProjFitY or drawProj) {
// If not already made, make the projections
if (!fMCHist_X) {
PlotUtils::MatchEmptyBins(fDataHist, fMCHist);
fMCHist_X = PlotUtils::GetProjectionX(fMCHist, fMaskHist);
fMCHist_Y = PlotUtils::GetProjectionY(fMCHist, fMaskHist);
fDataHist_X = PlotUtils::GetProjectionX(fDataHist, fMaskHist);
fDataHist_Y = PlotUtils::GetProjectionY(fDataHist, fMaskHist);
double chi2X = StatUtils::GetChi2FromDiag(fDataHist_X, fMCHist_X);
double chi2Y = StatUtils::GetChi2FromDiag(fDataHist_Y, fMCHist_Y);
fMCHist_X->SetTitle(Form("%f", chi2X));
fMCHist_Y->SetTitle(Form("%f", chi2Y));
}
// Save the histograms
fDataHist_X->Write();
fMCHist_X->Write();
fDataHist_Y->Write();
fMCHist_Y->Write();
}
if (drawSliceCanvYMC or true) {
TCanvas* c1 = new TCanvas((fName + "_MC_CANV_Y").c_str(),
(fName + "_MC_CANV_Y").c_str(), 800, 600);
c1->Divide(int(sqrt(fDataHist->GetNbinsY() + 1)),
int(sqrt(fDataHist->GetNbinsY() + 1)));
TH2D* mcShape = (TH2D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale =
fDataHist->Integral("width") / fMCHist->Integral("width");
mcShape->Scale(shapeScale);
mcShape->SetLineStyle(7);
c1->cd(1);
TLegend* leg = new TLegend(0.6, 0.6, 0.9, 0.9);
leg->AddEntry(fDataHist, (fName + " Data").c_str(), "ep");
leg->AddEntry(fMCHist, (fName + " MC").c_str(), "l");
leg->AddEntry(mcShape, (fName + " Shape").c_str(), "l");
leg->Draw("SAME");
/*
// Make Y slices
for (int i = 0; i < fDataHist->GetNbinY(); i++){
c1->cd(i+2);
TH1D* fDataHist_SliceY = PlotUtils::GetSliceY(fDataHist, i);
fDataHist_SliceY->Draw("E1");
TH1D* fMCHist_SliceY = PlotUtils::GetSliceY(fMCHist, i);
fMCHist_SliceY->Draw("SAME HIST C");
TH1D* mcShape_SliceY = PlotUtils::GetSliceY(mcShape, i);
mcShape_SliceY->Draw("SAME HIST C");
}
*/
c1->Write();
}
if (drawWeighted) {
fMCWeighted->Write();
}
// Returning
LOG(SAM) << "Written Histograms: " << fName << std::endl;
return;
// Returning
LOG(SAM) << "Written Histograms: " << fName << std::endl;
return;
}
/*
Setup Functions
*/
//********************************************************************
void Measurement2D::SetupMeasurement(std::string inputfile, std::string type,
FitWeight* rw, std::string fkdt) {
//********************************************************************
// 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;
}
fIsEnu = false;
if ((fName.find("XSec") != std::string::npos) &&
(fName.find("Enu") != std::string::npos)) {
fIsEnu = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
"not flux averaged!"
<< std::endl;
if (FitPar::Config().GetParB("EventManager")) {
ERR(FTL) << "Enu Measurements do not yet work with the Event Manager!"
<< std::endl;
ERR(FTL) << "If you want decent flux unfolded results please run in "
"series mode (-q EventManager=0)"
<< std::endl;
sleep(2);
}
}
if (fIsEnu && 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);
}
// Reset everything to NULL
fRW = rw;
// Setting up 2D Inputs
this->SetupInputs(inputfile);
// Set Default Options
SetFitOptions(fDefaultTypes);
// Set Passed Options
SetFitOptions(type);
}
//********************************************************************
void Measurement2D::SetupDefaultHist() {
//********************************************************************
// Setup fMCHist
fMCHist = (TH2D*)fDataHist->Clone();
fMCHist->SetNameTitle((fName + "_MC").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
// Setup fMCFine
Int_t nBinsX = fMCHist->GetNbinsX();
Int_t nBinsY = fMCHist->GetNbinsY();
fMCFine = new TH2D((fName + "_MC_FINE").c_str(),
(fName + "_MC_FINE" + fPlotTitles).c_str(), nBinsX * 3,
fMCHist->GetXaxis()->GetBinLowEdge(1),
fMCHist->GetXaxis()->GetBinLowEdge(nBinsX + 1), nBinsY * 3,
fMCHist->GetYaxis()->GetBinLowEdge(1),
fMCHist->GetYaxis()->GetBinLowEdge(nBinsY + 1));
// Setup MC Stat
fMCStat = (TH2D*)fMCHist->Clone();
fMCStat->Reset();
// Setup the NEUT Mode Array
PlotUtils::CreateNeutModeArray(fMCHist, (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);
}
return;
}
//********************************************************************
void Measurement2D::SetDataValues(std::string dataFile, std::string TH2Dname) {
//********************************************************************
if (dataFile.find(".root") == std::string::npos) {
ERR(FTL) << "Error! " << dataFile << " is not a .root file" << std::endl;
ERR(FTL) << "Currently only .root file reading is supported (MiniBooNE "
"CC1pi+ 2D), but implementing .txt should be dirt easy"
<< std::endl;
ERR(FTL) << "See me at " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
} else {
TFile* inFile = new TFile(dataFile.c_str(), "READ");
fDataHist = (TH2D*)(inFile->Get(TH2Dname.c_str())->Clone());
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle((fName + "_data").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
delete inFile;
}
return;
}
//********************************************************************
void Measurement2D::SetDataValues(std::string dataFile, double dataNorm,
std::string errorFile, double errorNorm) {
//********************************************************************
// Make a counter to track the line number
int yBin = 0;
std::string line;
std::ifstream data(dataFile.c_str(), ifstream::in);
fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(),
fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);
if (data.is_open())
LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
while (std::getline(data >> std::ws, line, '\n')) {
int xBin = 0;
// Loop over entries and insert them into the histogram
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::iterator iter = entries.begin();
iter != entries.end(); iter++) {
fDataHist->SetBinContent(xBin + 1, yBin + 1, (*iter) * dataNorm);
xBin++;
}
yBin++;
}
yBin = 0;
std::ifstream error(errorFile.c_str(), ifstream::in);
if (error.is_open())
LOG(SAM) << "Reading errors from: " << errorFile.c_str() << std::endl;
while (std::getline(error >> std::ws, line, '\n')) {
int xBin = 0;
// Loop over entries and insert them into the histogram
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::iterator iter = entries.begin();
iter != entries.end(); iter++) {
fDataHist->SetBinError(xBin + 1, yBin + 1, (*iter) * errorNorm);
xBin++;
}
yBin++;
}
return;
};
//********************************************************************
void Measurement2D::SetDataValuesFromText(std::string dataFile,
double dataNorm) {
//********************************************************************
fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(),
fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);
LOG(SAM) << "Reading data from: " << dataFile << std::endl;
PlotUtils::Set2DHistFromText(dataFile, fDataHist, dataNorm, true);
return;
};
//********************************************************************
void Measurement2D::SetCovarMatrix(std::string covarFile) {
//********************************************************************
// Used to read a covariance matrix from a root file
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
// Make plots that we want
TH2D* covarPlot = new TH2D();
// TH2D* decmpPlot = new TH2D();
TH2D* covarInvPlot = new TH2D();
TH2D* fFullCovarPlot = new TH2D();
// Get covariance options for fake data studies
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("throw_covariance");
// Which matrix to get?
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());
// Throw either the sub matrix or the full matrix
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;
// Bin masking?
int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
int covdim = int(fDataHist->GetNbinsX());
// Make new covars
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
fDecomp = new TMatrixDSym(dim);
// Full covariance values
int row, column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++) {
// masking can be dodgy
// 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) {
for (Int_t i = 0; i < fDataHist->GetNbinsX(); i++) {
fDataHist->SetBinError(
i + 1, sqrt((covarPlot->GetBinContent(i + 1, i + 1))) * 1E-38);
}
}
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
tempFile->Close();
delete tempFile;
return;
};
//********************************************************************
void Measurement2D::SetCovarMatrixFromText(std::string covarFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream covar(covarFile.c_str(), ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covar.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
while (std::getline(covar >> 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) * fDataHist->GetBinError(row + 1) * 1E38 *
fDataHist->GetBinError(column + 1) * 1E38;
(*this->covar)(row, column) = val;
(*fFullCovar)(row, column) = val;
column++;
}
row++;
}
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
void Measurement2D::SetCovarMatrixFromChol(std::string covarFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream covarread(covarFile.c_str(), ifstream::in);
TMatrixD* newcov = new TMatrixD(dim, dim);
if (covarread.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
while (std::getline(covarread >> 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++) {
(*newcov)(row, column) = *iter;
column++;
}
row++;
}
covarread.close();
// Form full covariance
TMatrixD* trans = (TMatrixD*)(newcov)->Clone();
trans->T();
(*trans) *= (*newcov);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
delete newcov;
delete trans;
// Robust matrix inversion method
TDecompChol LU = TDecompChol(*this->fFullCovar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
// //********************************************************************
// void Measurement2D::SetMapValuesFromText(std::string dataFile) {
// //********************************************************************
// fMapHist = new TH2I((fName + "_map").c_str(), (fName + fPlotTitles).c_str(),
// fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);
// LOG(SAM) << "Reading map from: " << dataFile << std::endl;
// PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0);
// return;
// };
diff --git a/src/FitBase/MeasurementBase.cxx b/src/FitBase/MeasurementBase.cxx
index b9f49a4..3d6eec5 100644
--- a/src/FitBase/MeasurementBase.cxx
+++ b/src/FitBase/MeasurementBase.cxx
@@ -1,532 +1,537 @@
// 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 "MeasurementBase.h"
/*
Constructor/Destructors
*/
//********************************************************************
// 2nd Level Constructor (Inherits From MeasurementBase.h)
MeasurementBase::MeasurementBase(void) {
//********************************************************************
fScaleFactor = 1.0;
fMCFilled = false;
fNoData = false;
fInput = NULL;
// Set the default values
// After-wards this gets set in SetupMeasurement
EnuMin = 0.;
EnuMax = 1.E5;
fMeasurementSpeciesType = kSingleSpeciesMeasurement;
fEventVariables = NULL;
};
void MeasurementBase::FinaliseMeasurement() {
// Used to setup default data hists, covars, etc.
}
//********************************************************************
// 2nd Level Destructor (Inherits From MeasurementBase.h)
MeasurementBase::~MeasurementBase() {
//********************************************************************
};
//********************************************************************
double MeasurementBase::TotalIntegratedFlux(std::string intOpt, double low,
double high) {
//********************************************************************
// Set Energy Limits
if (low == -9999.9) low = this->EnuMin;
if (high == -9999.9) high = this->EnuMax;
return GetInput()->TotalIntegratedFlux(low, high, intOpt);
};
//********************************************************************
double MeasurementBase::PredictedEventRate(std::string intOpt, double low,
double high) {
//********************************************************************
// Set Energy Limits
if (low == -9999.9) low = this->EnuMin;
if (high == -9999.9) high = this->EnuMax;
return GetInput()->PredictedEventRate(low, high, intOpt) * 1E-38;
};
//********************************************************************
void MeasurementBase::SetupInputs(std::string inputfile) {
//********************************************************************
// Add this infile to the global manager
if (FitPar::Config().GetParB("EventManager")) {
fInput = FitBase::AddInput(fName, inputfile);
} else {
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfile, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inpType =
InputUtils::ParseInputType(file_descriptor[0]);
fInput = InputUtils::CreateInputHandler(fName, inpType, file_descriptor[1]);
}
fNEvents = fInput->GetNEvents();
// Expect INPUTTYPE:FileLocation(s)
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfile, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
fInputType = InputUtils::ParseInputType(file_descriptor[0]);
fInputFileName = file_descriptor[1];
if (EnuMin == 0 && EnuMax == 1.E5) {
EnuMin = fInput->GetFluxHistogram()->GetBinLowEdge(1);
EnuMax = fInput->GetFluxHistogram()->GetBinLowEdge(
fInput->GetFluxHistogram()->GetNbinsX() + 1);
}
fFluxHist = fInput->GetFluxHistogram();
fEventHist = fInput->GetEventHistogram();
}
//***********************************************
int MeasurementBase::GetInputID() {
//***********************************************
return FitBase::GetInputID(fInputFileName);
}
//***********************************************
SampleSettings MeasurementBase::LoadSampleSettings(nuiskey samplekey) {
//***********************************************
SampleSettings setting = SampleSettings(samplekey);
fName = setting.GetS("name");
// Used as an initial setup function incase we need to do anything here.
LOG(SAM) << "Loading Sample : " << setting.GetName() << std::endl;
SetupInputs( setting.GetS("input") );
return setting;
}
//***********************************************
SampleSettings MeasurementBase::LoadSampleSettings(std::string name, std::string input, std::string type) {
//***********************************************
nuiskey samplekey = Config::CreateKey("sample");
samplekey.SetS("name",name);
samplekey.SetS("input",input);
samplekey.SetS("type",type);
return LoadSampleSettings(samplekey);
}
void MeasurementBase::FinaliseSampleSettings() {
+ EnuMin = fSettings.GetD("enu_min");
+ EnuMax = fSettings.GetD("enu_max");
+
+ std::cout << "SetEnuMin = " << EnuMin << " "<< EnuMax << std::endl;
+
}
//***********************************************
void MeasurementBase::Reconfigure() {
//***********************************************
LOG(REC) << " Reconfiguring sample " << fName << std::endl;
// Reset Histograms
ResetExtraHistograms();
AutoResetExtraTH1();
this->ResetAll();
// FitEvent* cust_event = fInput->GetEventPointer();
int fNEvents = fInput->GetNEvents();
int countwidth = (fNEvents / 5);
// MAIN EVENT LOOP
FitEvent* cust_event = fInput->FirstNuisanceEvent();
int i = 0;
int npassed = 0;
while(cust_event){
cust_event->RWWeight = fRW->CalcWeight(cust_event);
cust_event->Weight = cust_event->RWWeight * cust_event->InputWeight;
Weight = cust_event->Weight;
// Initialize
fXVar = -999.9;
fYVar = -999.9;
fZVar = -999.9;
Signal = false;
Mode = cust_event->Mode;
// Extract Measurement Variables
this->FillEventVariables(cust_event);
Signal = this->isSignal(cust_event);
if (Signal) npassed++;
GetBox()->SetX(fXVar);
GetBox()->SetY(fYVar);
GetBox()->SetZ(fZVar);
GetBox()->SetMode(Mode);
// GetBox()->fSignal = Signal;
// Fill Histogram Values
GetBox()->FillBoxFromEvent(cust_event);
// this->FillExtraHistograms(GetBox(), Weight);
this->FillHistogramsFromBox(GetBox(), Weight);
// Print Out
if (LOG_LEVEL(REC) && countwidth > 0 && !(i % countwidth)) {
std::stringstream ss("");
ss.unsetf(ios_base::fixed);
ss << std::setw(7) << std::right << i << "/" << fNEvents << " events ("
<< std::setw(2) << double(i) / double(fNEvents) * 100. << std::left
<< std::setw(5) << "%) "
<< "[S,X,Y,Z,M,W] = [" << std::fixed << std::setprecision(2)
<< std::right << Signal << ", " << std::setw(5) << fXVar << ", "
<< std::setw(5) << fYVar << ", " << std::setw(5) << fYVar << ", "
<< std::setw(3) << (int)Mode << ", " << std::setw(5) << Weight << "] "
<< std::endl;
LOG(SAM) << ss.str();
}
// iterate
cust_event = fInput->NextNuisanceEvent();
i++;
}
LOG(SAM) << npassed << "/" << fNEvents << " passed selection " << std::endl;
if (npassed == 0) {
LOG(SAM) << "WARNING: NO EVENTS PASSED SELECTION!" << std::endl;
}
LOG(REC) << std::setw(10) << std::right << NSignal << "/"
<< fNEvents << " events passed selection + binning after reweight"
<< std::endl;
// Finalise Histograms
fMCFilled = true;
this->ConvertEventRates();
}
void MeasurementBase::FillHistogramsFromBox(MeasurementVariableBox* var, double weight) {
fXVar = var->GetX();
fYVar = var->GetY();
fZVar = var->GetZ();
// Signal = var->fSignal;
// Mode = var->fMode;
Weight = weight;
FillHistograms();
FillExtraHistograms(var, weight);
}
void MeasurementBase::FillHistograms(double weight){
Weight = weight;
FillHistograms();
FillExtraHistograms(GetBox(), Weight);
}
MeasurementVariableBox* MeasurementBase::FillVariableBox(FitEvent* event) {
GetBox()->Reset();
Mode = event->Mode;
this->FillEventVariables(event);
Signal = this->isSignal(event);
GetBox()->FillBoxFromEvent(event);
GetBox()->SetX(fXVar);
GetBox()->SetY(fYVar);
GetBox()->SetZ(fZVar);
GetBox()->SetMode(event->Mode);
// GetBox()->fSignal = Signal;
return GetBox();
}
MeasurementVariableBox* MeasurementBase::GetBox() {
if (!fEventVariables) fEventVariables = CreateBox();
return fEventVariables;
}
//***********************************************
void MeasurementBase::ReconfigureFast() {
//***********************************************
this->Reconfigure();
}
//***********************************************
void MeasurementBase::ConvertEventRates() {
//***********************************************
AutoScaleExtraTH1();
ScaleExtraHistograms(GetBox());
this->ScaleEvents();
double normval = fRW->GetSampleNorm(this->fName);
if (normval < 0.01 or normval > 10.0){
ERR(WRN) << "Norm Value inside MeasurementBase::ConvertEventRates() looks off!" << std::endl;
ERR(WRN) << "It could have become out of sync with the minimizer norm list." << std::endl;
ERR(WRN) << "Setting it to 1.0" << std::endl;
normval = 1.0;
}
AutoNormExtraTH1(normval);
NormExtraHistograms(GetBox(), normval);
this->ApplyNormScale(normval);
}
//***********************************************
InputHandlerBase* MeasurementBase::GetInput() {
//***********************************************
if (!fInput) {
ERR(FTL) << "MeasurementBase::fInput not set. Please submit your command "
"line options and input cardfile with a bug report to: "
"nuisance@projects.hepforge.org"
<< std::endl;
throw;
}
return fInput;
};
//***********************************************
void MeasurementBase::Renormalise() {
//***********************************************
// Called when the fitter has changed a measurements normalisation but not any
// reweight dials
// Means we don't have to call the time consuming reconfigure when this
// happens.
double norm = fRW->GetDialValue(this->fName + "_norm");
if ((this->fCurrentNorm == 0.0 and norm != 0.0) or not fMCFilled) {
this->ReconfigureFast();
return;
}
if (this->fCurrentNorm == norm) return;
this->ApplyNormScale(1.0 / this->fCurrentNorm);
this->ApplyNormScale(norm);
return;
};
//***********************************************
void MeasurementBase::SetSignal(bool sig) {
//***********************************************
Signal = sig;
}
//***********************************************
void MeasurementBase::SetSignal(FitEvent* evt) {
//***********************************************
Signal = this->isSignal(evt);
}
//***********************************************
void MeasurementBase::SetWeight(double wght) {
//***********************************************
Weight = wght;
}
//***********************************************
void MeasurementBase::SetMode(int md) {
//***********************************************
Mode = md;
}
//***********************************************
std::vector<TH1*> MeasurementBase::GetFluxList() {
//***********************************************
return GetInput()->GetFluxList();
}
//***********************************************
std::vector<TH1*> MeasurementBase::GetEventRateList() {
//***********************************************
return GetInput()->GetEventList();
}
//***********************************************
std::vector<TH1*> MeasurementBase::GetXSecList() {
//***********************************************
return GetInput()->GetXSecList();
}
void MeasurementBase::ProcessExtraHistograms(int cmd,
MeasurementVariableBox* vars,
double weight) {
// This should be overriden if we have extra histograms!!!
// Add a flag to tell user this...
return;
}
void MeasurementBase::FillExtraHistograms(MeasurementVariableBox* vars,
double weight) {
ProcessExtraHistograms(kCMD_Fill, vars, weight);
}
void MeasurementBase::ScaleExtraHistograms(MeasurementVariableBox* vars) {
ProcessExtraHistograms(kCMD_Scale, vars, 1.0);
}
void MeasurementBase::ResetExtraHistograms() {
ProcessExtraHistograms(kCMD_Reset, NULL, 1.0);
}
void MeasurementBase::NormExtraHistograms(MeasurementVariableBox* vars,
double norm) {
ProcessExtraHistograms(kCMD_Norm, vars, norm);
}
void MeasurementBase::WriteExtraHistograms() {
ProcessExtraHistograms(kCMD_Write, NULL, 1.00);
}
void MeasurementBase::SetAutoProcessTH1(TH1* hist, int c1, int c2, int c3, int c4, int c5) {
FakeStack* fake = new FakeStack(hist);
SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command!
}
void MeasurementBase::SetAutoProcessTH1(StackBase* hist, int c1, int c2, int c3, int c4, int c5) {
// Set Defaults
// int ncommands = kCMD_extraplotflags;
bool autoflags[5];
autoflags[0] = false;
autoflags[1] = false;
autoflags[2] = false;
autoflags[3] = false;
autoflags[4] = false;
int givenflags[5];
givenflags[0] = c1;
givenflags[1] = c2;
givenflags[2] = c3;
givenflags[3] = c4;
givenflags[4] = c5;
fExtraTH1s[hist] = std::vector<int>(5,0);
// Setup a default one.
if (c1 == -1 && c2 == -1 && c3 == -1 && c4 == -1 && c5 == -1){
fExtraTH1s[hist][kCMD_Reset] = 1;
fExtraTH1s[hist][kCMD_Scale] = 1;
fExtraTH1s[hist][kCMD_Norm] = 1;
fExtraTH1s[hist][kCMD_Write] = 1;
}
for (int i = 0; i < 5; i++) {
switch (givenflags[i]) {
// Skip over...
case -1:
break;
case kCMD_Reset:
case kCMD_Scale:
case kCMD_Norm:
case kCMD_Write:
fExtraTH1s[hist][givenflags[i]] = 1;
break;
case kCMD_Fill:
ERR(FTL) << "Can't auto fill yet!" << std::endl;
autoflags[givenflags[i]] = 1;
break;
default:
break;
}
}
// LOG(SAM) << "AutoProcessing " << hist->GetName() << std::endl;
};
void MeasurementBase::AutoFillExtraTH1() {
ERR(FTL) << "Can't auto fill yet! it's too inefficent!" << std::endl;
return;
}
void MeasurementBase::AutoResetExtraTH1() {
for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!((*iter).second)[kCMD_Reset]) continue;
(*iter).first->Reset();
}
};
void MeasurementBase::AutoScaleExtraTH1() {
for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!((*iter).second)[kCMD_Scale]) continue;
(*iter).first->Scale(fScaleFactor, "width");
}
};
void MeasurementBase::AutoNormExtraTH1(double norm) {
double sfactor = 0.0;
if (norm != 0.0) sfactor = 1.0 / norm;
for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!((*iter).second)[kCMD_Norm]) continue;
(*iter).first->Scale(sfactor);
}
};
void MeasurementBase::AutoWriteExtraTH1() {
for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!(((*iter).second)[kCMD_Write])) continue;
(*iter).first->Write();
}
};
diff --git a/src/FitBase/SampleSettings.cxx b/src/FitBase/SampleSettings.cxx
index 9e58c35..2c7b3af 100644
--- a/src/FitBase/SampleSettings.cxx
+++ b/src/FitBase/SampleSettings.cxx
@@ -1,154 +1,161 @@
#include "SampleSettings.h"
SampleSettings::SampleSettings() {
};
SampleSettings::SampleSettings(nuiskey key) {
fKeyValues = key;
if (!key.Has("type")) key.AddS("type", "DEFAULT");
}
std::string SampleSettings::GetName() {
return GetS("name");
}
void SampleSettings::SetS(std::string name, std::string val) {
fKeyValues.SetS(name, val);
};
bool SampleSettings::Found(std::string name, std::string substr) {
std::string compstring = fKeyValues.GetS(name);
return compstring.find(substr) != std::string::npos;
}
void SampleSettings::SetXTitle(std::string name) {
SetDefault("xtitle", name);
};
void SampleSettings::SetYTitle(std::string name) {
SetDefault("ytitle", name);
};
void SampleSettings::SetZTitle(std::string name) {
SetDefault("ztitle", name);
};
void SampleSettings::SetNormError(double norm) {
SetDefault("norm_error", GeneralUtils::DblToStr(norm));
};
std::string SampleSettings::GetCovarInput() {
return GetS("covar");
}
void SampleSettings::SetAllowedTypes(std::string allowed, std::string defaulttype) {
SetDefault("default_types", allowed);
SetDefault("allowed_types", defaulttype);
};
void SampleSettings::SetEnuRangeFromFlux(TH1D* fluxhist) {
+ std::cout << "Setting Enu frmo flux " << std::endl;
double enu_min = fluxhist->GetXaxis()->GetXmin();
double enu_max = fluxhist->GetXaxis()->GetXmax();
SetEnuRange(enu_min, enu_max);
};
void SampleSettings::SetEnuRange(double min, double max) {
+ std::cout << "Setting EnuRange = " << min << " " << max << std::endl;
SetDefault("enu_min", min);
SetDefault("enu_max", max);
};
void SampleSettings::DefineAllowedTargets(std::string targ) {
// fAllowedTargets = TargetUtils::ParseTargetsToIntVect(targ);
};
void SampleSettings::FoundFill(std::string name, std::string substr, bool& cont, bool def) {
std::string compstring = fKeyValues.GetS(name);
if (compstring.find(substr) != std::string::npos) {
cont = def;
} else {
cont = !def;
}
};
void SampleSettings::SetTitle(std::string val) {
SetDefault("title", val);
};
void SampleSettings::SetDataInput(std::string val) {
SetDefault("data", val);
};
std::string SampleSettings::GetMapInput() {
return GetS("map");
}
void SampleSettings::SetMapInput(std::string val) {
SetDefault("map", val);
}
void SampleSettings::SetErrorInput(std::string val) {
SetDefault("error", val);
};
void SampleSettings::SetCovarInput(std::string val) {
SetDefault("covar", val);
};
void SampleSettings::SetDefault(std::string name, std::string val) {
if (!fKeyValues.Has(name)) fKeyValues.AddS(name, val);
};
void SampleSettings::SetDefault(std::string name, double val ) {
+ std::cout << "SetDefault D = " << val << std::endl;
if (!fKeyValues.Has(name)) fKeyValues.AddD(name, val);
};
void SampleSettings::SetHasExtraHistograms(bool opt) {
fHasExtraHistograms = opt;
};
void SampleSettings::DefineAllowedSpecies(std::string species) {
fAllowedTargets = BeamUtils::ParseSpeciesToIntVect(species);
};
std::string SampleSettings::Title() {
return GetS("title");
}
std::string SampleSettings::GetDataInput() {
return GetS("data");
};
std::string SampleSettings::GetErrorInput() {
return GetS("error");
};
std::string SampleSettings::PlotTitles() {
return ";" + GetS("xtitle") + ";" + GetS("ytitle");
};
std::string SampleSettings::GetFullTitles() {
return Title() + PlotTitles();
}
int SampleSettings::GetI(std::string name) {
return fKeyValues.GetI(name);
}
+double SampleSettings::GetD(std::string name) {
+ return fKeyValues.GetD(name);
+}
+
std::string SampleSettings::GetS(std::string name) {
return fKeyValues.GetS(name);
}
void SampleSettings::SetSuggestedFlux(std::string str) {
SetS("suggested_flux", str);
}
void SampleSettings::SetDescription(std::string str) {
SetS("description", str);
}
diff --git a/src/FitBase/SampleSettings.h b/src/FitBase/SampleSettings.h
index f54af9a..607635d 100644
--- a/src/FitBase/SampleSettings.h
+++ b/src/FitBase/SampleSettings.h
@@ -1,61 +1,62 @@
#ifndef SAMPLESETTINGS_H
#define SAMPLESETTINGS_H
#include "NuisConfig.h"
#include "NuisKey.h"
#include "TH1D.h"
#include "BeamUtils.h"
class SampleSettings {
public:
SampleSettings();
SampleSettings(nuiskey key);
inline std::string Name(){return GetName();};
std::string GetName();
void SetS(std::string name, std::string val);
void SetXTitle(std::string name);
void SetYTitle(std::string name);
void SetZTitle(std::string name);
void SetNormError(double norm);
void SetAllowedTypes(std::string allowed, std::string defaulttype="FIX");
void SetEnuRangeFromFlux(TH1D* fluxhist);
void SetEnuRange(double min, double max);
std::string Title();
void DefineAllowedTargets(std::string targ);
void FoundFill(std::string name, std::string substr, bool& cont, bool def);
// void FoundFill(std::string name, std::string substr, double& val, double )
bool Found(std::string name, std::string substr);
void SetTitle(std::string val);
void SetDataInput(std::string val);
void SetErrorInput(std::string val);
std::string GetErrorInput();
std::string GetMapInput();
void SetMapInput(std::string val);
void SetCovarInput(std::string val);
void SetDefault(std::string name, std::string val);
void SetDefault(std::string name, double val);
void SetHasExtraHistograms(bool opt = true);
void DefineAllowedSpecies(std::string species);
void SetSuggestedFlux(std::string str);
void SetDescription(std::string str);
std::string GetFullTitles();
bool Has(std::string name){return fKeyValues.Has(name);};
std::string GetDataInput();
std::string PlotTitles();
std::string GetS(std::string name);
int GetI(std::string name);
+ double GetD(std::string name);
std::string GetCovarInput();
std::vector<int> fAllowedTargets;
std::vector<int> fAllowedSpecies;
nuiskey fKeyValues;
bool fHasExtraHistograms;
};
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
index adc28ed..7f4d58f 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,274 +1,278 @@
#include "NUISANCEWeightCalcs.h"
GaussianModeCorr::GaussianModeCorr() {
// Init
fApply_CCQE = false;
fGausVal_CCQE[kPosNorm] = 0.0;
fGausVal_CCQE[kPosTilt] = 0.0;
fGausVal_CCQE[kPosPq0] = 1.0;
fGausVal_CCQE[kPosWq0] = 1.0;
fGausVal_CCQE[kPosPq3] = 1.0;
fGausVal_CCQE[kPosWq3] = 1.0;
fApply_2p2h = false;
fGausVal_2p2h[kPosNorm] = 0.0;
fGausVal_2p2h[kPosTilt] = 0.0;
fGausVal_2p2h[kPosPq0] = 1.0;
fGausVal_2p2h[kPosWq0] = 1.0;
fGausVal_2p2h[kPosPq3] = 1.0;
fGausVal_2p2h[kPosWq3] = 1.0;
fApply_2p2h_PPandNN = false;
fGausVal_2p2h_PPandNN[kPosNorm] = 0.0;
fGausVal_2p2h_PPandNN[kPosTilt] = 0.0;
fGausVal_2p2h_PPandNN[kPosPq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosPq3] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq3] = 1.0;
fApply_2p2h_NP = false;
fGausVal_2p2h_NP[kPosNorm] = 0.0;
fGausVal_2p2h_NP[kPosTilt] = 0.0;
fGausVal_2p2h_NP[kPosPq0] = 1.0;
fGausVal_2p2h_NP[kPosWq0] = 1.0;
fGausVal_2p2h_NP[kPosPq3] = 1.0;
fGausVal_2p2h_NP[kPosWq3] = 1.0;
fApply_CC1pi = false;
fGausVal_CC1pi[kPosNorm] = 0.0;
fGausVal_CC1pi[kPosTilt] = 0.0;
fGausVal_CC1pi[kPosPq0] = 1.0;
fGausVal_CC1pi[kPosWq0] = 1.0;
fGausVal_CC1pi[kPosPq3] = 1.0;
fGausVal_CC1pi[kPosWq3] = 1.0;
- fAllowSuppression = false;
-
fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG");
}
double GaussianModeCorr::CalcWeight(BaseFitEvt* evt) {
FitEvent* fevt = static_cast<FitEvent*>(evt);
double rw_weight = 1.0;
// Get Neutrino
FitParticle* pnu = fevt->PartInfo(0);
int pdgnu = pnu->fPID;
FitParticle* plep = fevt->GetHMFSParticle(abs(pdgnu) - 1);
if (!plep) return 1.0;
TLorentzVector q = pnu->fP - plep->fP;
// Extra q0,q3
double q0 = fabs(q.E()) / 1.E3;
double q3 = fabs(q.Vect().Mag()) / 1.E3;
int initialstate = -1; // Undef
if (abs(fevt->Mode) == 2) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < fevt->Npart(); j++) {
if ((fevt->PartInfo(j))->fIsAlive) continue;
if (fevt->PartInfo(j)->fPID == 2212) npr++;
else if (fevt->PartInfo(j)->fPID == 2112) nne++;
}
// std::cout << "PN State = " << npr << " " << nne << std::endl;
if (fevt->Mode == 2 and npr == 1 and nne == 1) {
initialstate = 2;
} else if (fevt->Mode == 2 and ((npr == 0 and nne == 2) or (npr == 2 and nne == 0))) {
initialstate = 1;
}
}
// std::cout << "Got q0 q3 = " << q0 << " " << q3 << std::endl;
// Apply weighting
if (fApply_CCQE and abs(fevt->Mode) == 1) {
if (fDebugStatements) std::cout << "Getting CCQE Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_CCQE);
}
if (fApply_2p2h and abs(fevt->Mode) == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
}
if (fApply_2p2h_PPandNN and abs(fevt->Mode) == 2 and initialstate == 1) {
if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
}
if (fApply_2p2h_NP and abs(fevt->Mode) == 2 and initialstate == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
}
if (fApply_CC1pi and abs(fevt->Mode) >= 11 and abs(fevt->Mode) <= 13) {
if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
}
- if (!fAllowSuppression){
- if (rw_weight < 1.0) rw_weight = 1.0;
- }
// if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
return rw_weight;
}
double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) {
+ // // CCQE Without Suppression
+ // double Norm = 4.82788679036;
+ // double Tilt = 2.3501416116;
+ // double Pq0 = 0.363964889702;
+ // double Wq0 = 0.133976806938;
+ // double Pq3 = 0.431769740224;
+ // double Wq3 = 0.207666663434;
+
+ // // Also add for CCQE at the end
+ // return (w > 1.0) ? w : 1.0;
+
+ // // 2p2h with suppression
+ // double Norm = 15.967;
+ // double Tilt = -0.455655;
+ // double Pq0 = 0.214598;
+ // double Wq0 = 0.0291061;
+ // double Pq3 = 0.480194;
+ // double Wq3 = 0.134588;
+
+
double Norm = vals[kPosNorm];
double Tilt = vals[kPosTilt];
double Pq0 = vals[kPosPq0];
double Wq0 = vals[kPosWq0];
double Pq3 = vals[kPosPq3];
double Wq3 = vals[kPosWq3];
- // double w = Norm;
- // w *= TMath::Gaus(q0, Pq0, Wq0);
- // w *= TMath::Gaus(q3, Pq3, Wq3);
-
- // w = Norm * exp( -0.5 * ((x - [1])/[2])**2 - 0.5 *((y-[3])/[4])) );
-
double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);
double b = - sin(2 * Tilt) / (4 * Wq0 * Wq0);
b += sin(2 * Tilt) / (4 * Wq3 * Wq3);
double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);
double w = Norm;
w *= exp(-a * (q0 - Pq0) * (q0 - Pq0));
w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3));
w *= exp(-c * (q3 - Pq3) * (q3 - Pq3));
if (fDebugStatements) {
- // std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " << sin(Tilt) << std::endl;
- // std::cout << "abc = " << a << " " << b << " " << c << std::endl;
+ std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " << sin(Tilt) << std::endl;
+ std::cout << "abc = " << a << " " << b << " " << c << std::endl;
std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " " << Pq3 << " " << Wq3 << " " << w << std::endl;
}
return w;
}
void GaussianModeCorr::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void GaussianModeCorr::SetDialValue(int rwenum, double val) {
int curenum = rwenum % 1000;
// Check Handled
if (!IsHandled(curenum)) return;
// CCQE Setting
for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CCQE_norm;
fGausVal_CCQE[index] = val;
fApply_CCQE = true;
}
}
// 2p2h Setting
for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_norm;
fGausVal_2p2h[index] = val;
fApply_2p2h = true;
}
}
// 2p2h_PPandNN Setting
for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_PPandNN_norm;
fGausVal_2p2h_PPandNN[index] = val;
fApply_2p2h_PPandNN = true;
}
}
// 2p2h_NP Setting
for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_NP_norm;
fGausVal_2p2h_NP[index] = val;
fApply_2p2h_NP = true;
}
}
// CC1pi Setting
for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CC1pi_norm;
fGausVal_CC1pi[index] = val;
fApply_CC1pi = true;
}
}
- if (curenum == kGaussianCorr_AllowSuppression){
- fAllowSuppression = bool(val);
- }
-
}
bool GaussianModeCorr::IsHandled(int rwenum) {
int curenum = rwenum % 1000;
switch (curenum) {
case kGaussianCorr_CCQE_norm:
case kGaussianCorr_CCQE_tilt:
case kGaussianCorr_CCQE_Pq0:
case kGaussianCorr_CCQE_Wq0:
case kGaussianCorr_CCQE_Pq3:
case kGaussianCorr_CCQE_Wq3:
case kGaussianCorr_2p2h_norm:
case kGaussianCorr_2p2h_tilt:
case kGaussianCorr_2p2h_Pq0:
case kGaussianCorr_2p2h_Wq0:
case kGaussianCorr_2p2h_Pq3:
case kGaussianCorr_2p2h_Wq3:
case kGaussianCorr_2p2h_PPandNN_norm:
case kGaussianCorr_2p2h_PPandNN_tilt:
case kGaussianCorr_2p2h_PPandNN_Pq0:
case kGaussianCorr_2p2h_PPandNN_Wq0:
case kGaussianCorr_2p2h_PPandNN_Pq3:
case kGaussianCorr_2p2h_PPandNN_Wq3:
case kGaussianCorr_2p2h_NP_norm:
case kGaussianCorr_2p2h_NP_tilt:
case kGaussianCorr_2p2h_NP_Pq0:
case kGaussianCorr_2p2h_NP_Wq0:
case kGaussianCorr_2p2h_NP_Pq3:
case kGaussianCorr_2p2h_NP_Wq3:
case kGaussianCorr_CC1pi_norm:
case kGaussianCorr_CC1pi_tilt:
case kGaussianCorr_CC1pi_Pq0:
case kGaussianCorr_CC1pi_Wq0:
case kGaussianCorr_CC1pi_Pq3:
case kGaussianCorr_CC1pi_Wq3:
- case kGaussianCorr_AllowSuppression:
return true;
default:
return false;
}
}
diff --git a/src/Utils/GeneralUtils.cxx b/src/Utils/GeneralUtils.cxx
index 227d914..388c58a 100644
--- a/src/Utils/GeneralUtils.cxx
+++ b/src/Utils/GeneralUtils.cxx
@@ -1,183 +1,183 @@
// 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 "GeneralUtils.h"
std::string GeneralUtils::BoolToStr(bool val) {
std::ostringstream ss;
ss << val;
return ss.str();
}
std::string GeneralUtils::IntToStr(int val) {
std::ostringstream ss;
ss << val;
return ss.str();
};
-std::string GeneralUtils::DblToStr(int val) {
+std::string GeneralUtils::DblToStr(double val) {
std::ostringstream ss;
ss << val;
return ss.str();
};
std::vector<std::string> GeneralUtils::LoadCharToVectStr(int argc, char* argv[]){
std::vector<std::string> vect;
for (int i = 1; i < argc; i++){
vect.push_back( std::string(argv[i]) );
}
return vect;
}
std::vector<std::string> GeneralUtils::ParseToStr(std::string str, const char* del) {
std::istringstream stream(str);
std::string temp_string;
std::vector<std::string> vals;
while (std::getline(stream >> std::ws, temp_string, *del)) {
if (temp_string.empty()) continue;
vals.push_back(temp_string);
}
return vals;
}
std::vector<double> GeneralUtils::ParseToDbl(std::string str, const char* del) {
std::istringstream stream(str);
std::string temp_string;
std::vector<double> vals;
while (std::getline(stream >> std::ws, temp_string, *del)) {
if (temp_string.empty()) continue;
std::istringstream stream(temp_string);
double entry;
stream >> entry;
vals.push_back(entry);
}
return vals;
}
std::vector<int> GeneralUtils::ParseToInt(std::string str, const char* del) {
std::istringstream stream(str);
std::string temp_string;
std::vector<int> vals;
while (std::getline(stream >> std::ws, temp_string, *del)) {
if (temp_string.empty()) continue;
std::istringstream stream(temp_string);
int entry;
stream >> entry;
vals.push_back(entry);
}
return vals;
}
// To stop rooku's skin from crawling :p
double GeneralUtils::StrToDbl(std::string str) {
std::istringstream stream(str);
double val;
stream >> val;
return val;
}
int GeneralUtils::StrToInt(std::string str) {
std::istringstream stream(str);
int val;
stream >> val;
return val;
}
bool GeneralUtils::StrToBool(std::string str) {
// convert result to lower case
// for (int i = 0; i < str.size(); i++) str[i] = tolower(str[i]);
// Test for true/false
if (!str.compare("false")) return false;
else if (!str.compare("true") ) return true;
if (str.empty()) return false;
// Push into bool
std::istringstream stream(str);
bool val;
stream >> val;
return val;
}
std::vector<std::string> GeneralUtils::ParseFileToStr(std::string str, const char* del) {
std::vector<std::string> linevect;
std::string line;
std::ifstream read;
read.open(str.c_str());
if (!read.is_open()) {
ERR(FTL) << "Cannot open file " << str << " in ParseFileToStr" << std::endl;
throw;
}
while ( std::getline(read >> std::ws, line, *del) ) {
linevect.push_back(line);
}
read.close();
return linevect;
}
std::string GeneralUtils::GetTopLevelDir() {
static bool first = true;
static std::string topLevelVarVal;
if (first) {
char * const var = getenv("EXT_FIT");
if (!var) {
ERR(FTL) << "Cannot find top level directory! Set the EXT_FIT environmental variable" << std::endl;
exit(-1);
}
topLevelVarVal = std::string(var);
first = false;
}
return topLevelVarVal;
}
diff --git a/src/Utils/GeneralUtils.h b/src/Utils/GeneralUtils.h
index f6e030f..e34dadf 100644
--- a/src/Utils/GeneralUtils.h
+++ b/src/Utils/GeneralUtils.h
@@ -1,159 +1,159 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef GENERALUTILS_H_SEEN
#define GENERALUTILS_H_SEEN
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#include <cstring>
#include <fstream>
#include <iostream>
#include <iostream>
#include <numeric>
#include <limits>
#include <sstream>
#include <string>
#include <vector>
#include "FitLogger.h"
/*!
* \addtogroup Utils
* @{
*/
//! Functions which deal with basic string and file handling. They should have
//! no dependence on the other NUISANCE files!
namespace GeneralUtils {
/*!
String handling and file parsing functions
*/
//! Parse a string into a vector of doubles given a delimiter "del"
std::vector<double> ParseToDbl(std::string str, const char* del);
//! Parse a string into a vector of ints given a delimiter "del"
std::vector<int> ParseToInt(std::string str, const char* del);
//! Parse a string into a vector of strings given a delimiter "del"
std::vector<std::string> ParseToStr(std::string str, const char* del);
//! Parse text file into a vector of strings
std::vector<std::string> ParseFileToStr(std::string str, const char* del);
//! Convert a string to a double
double StrToDbl(std::string str);
//! Convert a string to an int
int StrToInt(std::string str);
//! Convert a string to an bool
bool StrToBool(std::string str);
std::string BoolToStr(bool val);
std::string IntToStr(int val);
-std::string DblToStr(int val);
+std::string DblToStr(double val);
//! Return the top level environmental variable for the fitter
std::string GetTopLevelDir();
// //! A utility function to return a std::vector from an array
// template <typename T, size_t N>
// std::vector<T> makeVector(const T (&data)[N]) {
// return std::vector<T>(data, data + N);
// }
std::vector<std::string> LoadCharToVectStr(int argc, char* argv[]);
template <typename T, size_t N>
size_t GetArraySize(const T (&data)[N]) {
return N;
}
template <typename T>
size_t GetHammingWeight(T const& d) {
T c = d;
size_t w = 0;
while (bool(c)) {
w += c & 1;
c = (c >> 1);
}
return w;
}
template <typename T>
size_t GetFirstOnBit(T const& d) {
T c = d;
size_t fob = 0;
while (bool(c)) {
if (c & 1) {
return fob;
} else {
c = (c >> 1);
}
fob++;
}
return fob;
}
template <typename T>
size_t IsSmallNum(T const& d) {
if (std::numeric_limits<T>::is_integer) {
return (d == 0);
}
return (((d > 0) && (d < std::numeric_limits<T>::epsilon())) ||
((d < 0) && (d > -std::numeric_limits<T>::epsilon())));
}
}
namespace PhysConst {
const double mass_proton = 0.93827203; // Proton mass in GeV
const double mass_neutron = 0.93956536; // Neutron mass in GeV
const double mass_delta = 1.232; // Delta mass in GeV
const double mass_muon = 0.10565837; // Muon mass in GeV
const int pdg_neutrinos[] = {12, -12, 14, -14 /*, 16, -16*/};
const int pdg_muons[] = {13, -13};
const int pdg_pions[] = {211, -211, 111};
const int pdg_charged_pions[] = {211, -211};
const int pdg_strangemesons[] = {
130, 310, 311, 321, 9000311, 9000321, 10311,
10321, 100311, 100321, 9010311, 9010321, 9020311, 9020321,
313, 323, 10313, 10323, 20313, 20323, 100313,
100323, 9000313, 9000323, 30313, 30323, 315, 325,
9000315, 9000325, 10315, 10325, 20315, 20325, 9010315,
9010325, 9020315, 9020325, 317, 327, 9010317, 9010327};
// Just *-1 to cover possibility
const int pdg_kplus = 321;
const int pdg_antistrangemesons[] = {
-130, -310, -311, -321, -9000311, -9000321, -10311,
-10321, -100311, -100321, -9010311, -9010321, -9020311, -9020321,
-313, -323, -10313, -10323, -20313, -20323, -100313,
-100323, -9000313, -9000323, -30313, -30323, -315, -325,
-9000315, -9000325, -10315, -10325, -20315, -20325, -9010315,
-9010325, -9020315, -9020325, -317, -327, -9010317, -9010327};
}
/*! @} */
#endif
diff --git a/src/Utils/StackBase.cxx b/src/Utils/StackBase.cxx
index dd2120c..1a2abc3 100644
--- a/src/Utils/StackBase.cxx
+++ b/src/Utils/StackBase.cxx
@@ -1,204 +1,204 @@
#include "StackBase.h"
void StackBase::AddMode(std::string name, std::string title,
int linecolor, int linewidth,
int fillstyle) {
// int ncur = fAllLabels.size();
fAllLabels.push_back(name);
fAllTitles.push_back(title);
std::vector<int> temp;
temp.push_back(linecolor);
temp.push_back(linewidth);
temp.push_back(fillstyle);
fAllStyles.push_back(temp);
}
void StackBase::AddMode(int index, std::string name, std::string title,
int linecolor, int linewidth, int fillstyle) {
while (fAllLabels.size() <= (UInt_t)index) {
fAllLabels.push_back("");
fAllTitles.push_back("");
fAllStyles.push_back(std::vector<int>(1, 1));
}
fAllLabels[index] = (name);
fAllTitles[index] = (title);
std::vector<int> temp;
temp.push_back(linecolor);
temp.push_back(linewidth);
temp.push_back(fillstyle);
fAllStyles[index] = temp;
}
bool StackBase::IncludeInStack(TH1* hist){
if (!FitPar::Config().GetParB("includeemptystackhists") and
hist->Integral() == 0.0) return false;
return true;
}
bool StackBase::IncludeInStack(int index) {
return true;
}
void StackBase::SetupStack(TH1* hist) {
fTemplate = (TH1*) hist->Clone(fName.c_str());
fTemplate->Reset();
// Determine template dim
fNDim = fTemplate->GetDimension();
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists.push_back( (TH1*) fTemplate->Clone( (fName + "_" + fAllLabels[i]).c_str() ) );
}
};
void StackBase::Scale(double sf, std::string opt) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
- std::cout << "Scaling Stack Hist " << i << " by " << sf << std::endl;
+ // std::cout << "Scaling Stack Hist " << i << " by " << sf << std::endl;
fAllHists[i]->Scale(sf, opt.c_str());
}
};
void StackBase::Reset() {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Reset();
}
};
void StackBase::FillStack(int index, double x, double y, double z, double weight) {
if (index < 0 or (UInt_t)index >= fAllLabels.size()) {
ERR(WRN) << "Returning Stack Fill Because Range = " << index << " " << fAllLabels.size() << std::endl;
return;
}
if (fNDim == 1) fAllHists[index]->Fill(x, y);
else if (fNDim == 2) {
// std::cout << "Filling 2D Stack " << index << " " << x << " " << y << " " << z << std::endl;
((TH2*)fAllHists[index])->Fill(x, y, z);
}
else if (fNDim == 3) ((TH3*)fAllHists[index])->Fill(x, y, z, weight);
}
void StackBase::Write() {
THStack* st = new THStack();
// Loop and add all histograms
bool saveseperate = FitPar::Config().GetParB("WriteSeperateStacks");
for (size_t i = 0; i < fAllLabels.size(); i++) {
if (!IncludeInStack(fAllHists[i])) continue;
if (!IncludeInStack(i)) continue;
fAllHists[i]->SetTitle( fAllTitles[i].c_str() );
fAllHists[i]->GetXaxis()->SetTitle( fXTitle.c_str() );
fAllHists[i]->GetYaxis()->SetTitle( fYTitle.c_str() );
fAllHists[i]->GetZaxis()->SetTitle( fZTitle.c_str() );
fAllHists[i]->SetLineColor( fAllStyles[i][0] );
fAllHists[i]->SetLineWidth( fAllStyles[i][1] );
fAllHists[i]->SetFillStyle( fAllStyles[i][2] );
fAllHists[i]->SetFillColor( fAllStyles[i][0] );
if (saveseperate) fAllHists[i]->Write();
st->Add(fAllHists[i]);
}
st->SetTitle(fTitle.c_str());
st->SetName(fName.c_str());
st->Write();
delete st;
};
void StackBase::Multiply(TH1* hist) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Multiply(hist);
}
}
void StackBase::Divide(TH1* hist) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Divide(hist);
}
}
void StackBase::Add(TH1* hist, double scale) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Add(hist, scale);
}
}
void StackBase::Add(StackBase* hist, double scale) {
if (hist->GetType() != fType) {
ERR(WRN) << "Trying to add two StackBases of different types!" << std::endl;
ERR(WRN) << fType << " + " << hist->GetType() << " = Undefined." << std::endl;
ERR(WRN) << "Doing nothing..." << std::endl;
return;
}
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Add(hist->GetHist(i));
}
}
TH1* StackBase::GetHist(int entry) {
return fAllHists[entry];
}
TH1* StackBase::GetHist(std::string label) {
TH1* hist = NULL;
std::vector<std::string> splitlabels = GeneralUtils::ParseToStr(label, "+");
for (size_t j = 0; j < splitlabels.size(); j++) {
std::string newlabel = splitlabels[j];
for (size_t i = 0; i < fAllLabels.size(); i++) {
if (newlabel == fAllLabels[i]) {
if (!hist) hist = (TH1*) fAllHists[i]->Clone();
else hist->Add(fAllHists[i]);
}
}
}
return hist;
}
THStack StackBase::GetStack() {
THStack st = THStack();
for (size_t i = 0; i < fAllLabels.size(); i++) {
st.Add(fAllHists[i]);
}
return st;
}
void StackBase::AddNewHist(std::string name, TH1* hist) {
AddMode( fAllLabels.size(), name, hist->GetTitle(), hist->GetLineColor());
fAllHists.push_back( (TH1*) hist->Clone() );
}
void StackBase::AddToCategory(std::string name, TH1* hist) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
if (name == fAllLabels[i]) {
fAllHists[i]->Add(hist);
break;
}
}
}
void StackBase::AddToCategory(int index, TH1* hist) {
fAllHists[index]->Add(hist);
}