diff --git a/app/nuissmear.cxx b/app/nuissmear.cxx
index e88184b..878f0d5 100644
--- a/app/nuissmear.cxx
+++ b/app/nuissmear.cxx
@@ -1,223 +1,222 @@
// 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 .
*******************************************************************************/
#include "ComparisonRoutines.h"
#include "InputUtils.h"
#include "Smearceptance_Tester.h"
// Global Arguments
std::string gOptInputFile = "";
std::string gOptOutputFile = "";
std::string gOptType = "DEFAULT";
std::string gOptNumberEvents = "NULL";
std::string gOptCardInput = "";
std::string gOptOptions = "";
//*******************************
void PrintSyntax() {
//*******************************
std::cout << "nuisflat -i input [-o outfile] [-n nevents] [-t "
"options] [-q con=val] \n";
std::cout
<< "\n Arguments : "
<< "\n\t -i input : Path to input vector of events to flatten"
<< "\n\t"
<< "\n\t This should be given in the same format a normal "
"input file"
<< "\n\t is given to NUISANCE. {e.g. NUWRO:eventsout.root}."
<< "\n\t"
<< "\n\t[-c crd.xml]: Input card file to override configs or define "
"smearcepters."
<< "\n\t "
<< "\n\t[-o outfile]: Optional output file path. "
<< "\n\t "
<< "\n\t If none given, input.smear.root is chosen."
<< "\n\t"
<< "\n\t[-n nevents]: Optional choice of Nevents to run over. Default is "
"all."
<< "\n\t"
<< "\n\t[-t options]: Pass OPTION to the smearception sample. "
<< "\n\t Similar to type field in comparison xml configs."
<< "\n\t"
<< "\n\t[-q con=val]: Configuration overrides." << std::endl;
exit(-1);
};
//____________________________________________________________________________
void GetCommandLineArgs(int argc, char** argv) {
// Check for -h flag.
for (int i = 0; i < argc; i++) {
if ((!std::string(argv[i]).compare("-h")) ||
(!std::string(argv[i]).compare("-?")) ||
(!std::string(argv[i]).compare("--help")))
PrintSyntax();
}
// Format is nuwro -r run_number -n n events
std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv);
// Parse input file
ParserUtils::ParseArgument(args, "-i", gOptInputFile, false);
if (gOptInputFile == "") {
THROW("Need to provide a valid input file to nuisflat using -i flag!");
} else {
LOG(FIT) << "Reading Input File = " << gOptInputFile << std::endl;
gOptInputFile = InputUtils::PrependGuessedInputTypeToName(gOptInputFile);
}
// Get Output File
ParserUtils::ParseArgument(args, "-o", gOptOutputFile, false);
if (gOptOutputFile == "") {
gOptOutputFile = gOptInputFile + ".smear.root";
LOG(FIT) << "No output file given so saving nuisflat output to:"
<< gOptOutputFile << std::endl;
} else {
LOG(FIT) << "Saving nuisflat output to " << gOptOutputFile << std::endl;
}
// Get N Events and Configs
nuisconfig configuration = Config::Get();
ParserUtils::ParseArgument(args, "-n", gOptNumberEvents, false);
if (gOptNumberEvents.compare("NULL")) {
configuration.OverrideConfig("MAXEVENTS=" + gOptNumberEvents);
}
std::vector configargs;
ParserUtils::ParseArgument(args, "-q", configargs);
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
ParserUtils::ParseArgument(args, "-c", gOptCardInput, false);
if (gOptCardInput != "") {
QLOG(FIT, "Reading cardfile: " << gOptCardInput);
configuration.LoadConfig(gOptCardInput, "");
}
ParserUtils::ParseArgument(args, "-t", gOptOptions, false);
if (gOptOptions != "") {
QLOG(FIT, "Read options: \"" << gOptOptions << "\'");
}
return;
}
void SetupRW() {
std::vector parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
std::vector Params;
std::map TypeVals;
std::map CurrVals;
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
ERR(FTL) << "type='PARAMETER_TYPE'" << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
ERR(FTL) << "name='SAMPLE_NAME'" << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
ERR(FTL) << "nominal='NOMINAL_VALUE'" << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
// Push into vectors
Params.push_back(parname);
TypeVals[parname] = FitBase::ConvDialType(partype);
CurrVals[parname] = parnom;
}
for (UInt_t i = 0; i < Params.size(); i++) {
FitBase::GetRW()->IncludeDial(Params[i], TypeVals[Params[i]]);
FitBase::GetRW()->SetDialValue(Params[i], CurrVals[Params[i]]);
}
FitBase::GetRW()->Reconfigure();
}
//*******************************
int main(int argc, char* argv[]) {
//*******************************
// Parse
GetCommandLineArgs(argc, argv);
// Make output file
TFile* f = new TFile(gOptOutputFile.c_str(), "RECREATE");
if (f->IsZombie()) {
THROW("Cannot create output file!");
}
f->cd();
FitPar::Config().out = f;
// Create a new measurementbase class depending on the Format
MeasurementBase* flattreecreator = NULL;
// Make a new sample key for the format of interest.
nuiskey samplekey = Config::CreateKey("sample");
- samplekey.AddS("name", "FlatTree");
- samplekey.AddS("input", gOptInputFile);
- samplekey.AddS("type", gOptType);
+ samplekey.SetS("name", std::string("FlatTree_") + gOptOptions);
+ samplekey.SetS("input", gOptInputFile);
+ samplekey.SetS("type", gOptType);
+
if (gOptOptions == "") {
THROW(
"Attempting to flatten with Smearceptor, but no Smearceptor given. "
"Please supply a -t option.");
}
if (gOptCardInput == "") {
THROW(
"Attempting to flatten with Smearceptor, but no card passed with "
"Smearceptors configured. Please supply a -c option.");
}
SetupRW();
- flattreecreator =
- new Smearceptance_Tester(std::string("FlatTree_") + gOptOptions,
- gOptInputFile, FitBase::GetRW(), gOptType, "");
+ flattreecreator = new Smearceptance_Tester(samplekey);
// Make the FlatTree reconfigure
flattreecreator->Reconfigure();
f->cd();
flattreecreator->Write();
f->Close();
// Show Final Status
LOG(FIT) << "-------------------------------------" << std::endl;
LOG(FIT) << "Flattree Generation Complete." << std::endl;
LOG(FIT) << "-------------------------------------" << std::endl;
return 0;
}
diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx
index 70a03a5..e1c9936 100644
--- a/src/FitBase/Measurement1D.cxx
+++ b/src/FitBase/Measurement1D.cxx
@@ -1,1896 +1,1903 @@
// Copyright 2016 L. Pickering, P. Stowell, 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 .
*******************************************************************************/
#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;
fShapeCovar = NULL;
-
+
fCovar = NULL;
fInvert = NULL;
fDecomp = NULL;
-
+
// Fake Data
fFakeDataInput = "";
fFakeDataFile = NULL;
// Options
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsNoWidth = false;
fIsDifXSec = false;
fIsEnu1D = false;
// Inputs
fInput = NULL;
fRW = NULL;
// Extra Histograms
fMCHist_Modes = NULL;
}
//********************************************************************
Measurement1D::~Measurement1D(void) {
//********************************************************************
if (fDataHist) delete fDataHist;
if (fDataTrue) delete fDataTrue;
if (fMCHist) delete fMCHist;
if (fMCFine) delete fMCFine;
if (fMCWeighted) delete fMCWeighted;
if (fMaskHist) delete fMaskHist;
if (covar) delete covar;
if (fFullCovar) delete fFullCovar;
if (fShapeCovar) delete fShapeCovar;
if (fCovar) delete fCovar;
if (fInvert) delete fInvert;
if (fDecomp) delete fDecomp;
}
//********************************************************************
void Measurement1D::FinaliseSampleSettings() {
//********************************************************************
MeasurementBase::FinaliseSampleSettings();
// Setup naming + renaming
fName = fSettings.GetName();
fSettings.SetS("originalname", fName);
if (fSettings.Has("rename")) {
fName = fSettings.GetS("rename");
fSettings.SetS("name", fName);
}
// Setup all other options
LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl;
if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but using poisson likelihoods."
<< std::endl;
}
if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
<< "not flux averaged!" << std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
if (!fRW) fRW = FitBase::GetRW();
if (!fInput and !fIsJoint) SetupInputs(fSettings.GetS("input"));
// Setup options
SetFitOptions(fDefaultTypes); // defaults
SetFitOptions(fSettings.GetS("type")); // user specified
EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));
if (fAddNormPen) {
if (fNormError <= 0.0) {
ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl;
ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl;
throw;
}
}
}
//********************************************************************
void Measurement1D::CreateDataHistogram(int dimx, double* binx) {
//********************************************************************
if (fDataHist) delete fDataHist;
fDataHist = new TH1D( (fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(),
dimx, binx) ;
}
//********************************************************************
void Measurement1D::SetDataFromTextFile(std::string datafile) {
//********************************************************************
LOG(SAM) << "Reading data from text file: " << datafile << std::endl;
fDataHist = PlotUtils::GetTH1DFromFile(datafile,
fSettings.GetName() + "_data",
fSettings.GetFullTitles());
}
//********************************************************************
void Measurement1D::SetDataFromRootFile(std::string datafile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname);
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
return;
};
//********************************************************************
void Measurement1D::SetEmptyData(){
//********************************************************************
fDataHist = new TH1D("EMPTY_DATA","EMPTY_DATA",1,0.0,1.0);
}
//********************************************************************
void Measurement1D::SetPoissonErrors() {
//********************************************************************
if (!fDataHist) {
ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl;
ERR(FTL) << "Setup Data First!" << std::endl;
throw;
}
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
}
}
//********************************************************************
void Measurement1D::SetCovarFromDiagonal(TH1D* data) {
//********************************************************************
if (!data and fDataHist) {
data = fDataHist;
}
if (data) {
LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl;
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
} else {
ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl;
}
// if (!fIsDiag) {
// ERR(FTL) << "SetCovarMatrixFromDiag called for measurement "
// << "that is not set as diagonal." << std::endl;
// throw;
// }
}
//********************************************************************
void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl;
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarFromMultipleTextFiles(std::string covfiles, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector covList = GeneralUtils::ParseToStr(covfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < covList.size(); ++i){
LOG(SAM) << "Reading covariance from text file: " << covList[i] << std::endl;
TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim);
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl;
fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl;
covar = StatUtils::GetCovarFromTextFile(covfile, dim);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl;
covar = StatUtils::GetCovarFromRootFile(covfile, histname);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) dim = fDataHist->GetNbinsX();
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim) {
-//********************************************************************
+//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector corrList = GeneralUtils::ParseToStr(corrfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < corrList.size(); ++i){
LOG(SAM) << "Reading covariance from text file: " << corrList[i] << std::endl;
TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim);
-
+
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void Measurement1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromRootFile(covfile, histname);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
void Measurement1D::SetShapeCovar(){
-
+
// Return if this is missing any pre-requisites
if (!fFullCovar) return;
if (!fDataHist) return;
// Also return if it's bloody stupid under the circumstances
if (fIsDiag) return;
-
+
fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist);
return;
}
//********************************************************************
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 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;
if (fSettings.GetB("onlymc")){
if (fDataHist) delete fDataHist;
fDataHist = new TH1D("empty_data","empty_data",1,0.0,1.0);
}
// 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) {
fIsDiag = true;
SetCovarFromDiagonal(fDataHist);
}
if (!covar) {
covar = StatUtils::GetInvert(fFullCovar);
}
if (!fDecomp) {
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
// Push the diagonals of fFullCovar onto the data histogram
// Comment this out until the covariance/data scaling is consistent!
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38);
// If shape only, set covar and fDecomp using the shape-only matrix (if set)
if (fIsShape && fShapeCovar){
if (covar) delete covar;
covar = StatUtils::GetInvert(fShapeCovar);
if (fDecomp) delete 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() * 8,
fMCHist->GetBinLowEdge(1),
fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1));
fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCFine->Reset();
// Setup MC Stat
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
// Search drawopts for possible types to include by default
std::string drawopts = FitPar::Config().GetParS("drawopts");
if (drawopts.find("MODES") != std::string::npos) {
fMCHist_Modes = new TrueModeStack( (fSettings.GetName() + "_MODES").c_str(),
("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);
}
// Setup bin masks using sample name
if (fIsMask) {
std::string curname = fName;
std::string origname = fSettings.GetS("originalname");
// Check rename.mask
std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");
// Check origname.mask
if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask");
// Check database
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
}
// Setup Bin Mask
SetBinMask(maskloc);
}
if (fScaleFactor < 0) {
ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl;
ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl;
ERR(FTL) << "EXITING" << std::endl;
throw;
}
// Create and fill Weighted Histogram
if (!fMCWeighted) {
fMCWeighted = (TH1D*)fMCHist->Clone();
fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
(fName + "_MCWGHTS" + fPlotTitles).c_str());
fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
}
}
//********************************************************************
void Measurement1D::SetFitOptions(std::string opt) {
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT") return;
// CHECK Conflicting Fit Options
std::vector fit_option_allow =
GeneralUtils::ParseToStr(fAllowedTypes, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
std::vector 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 fit_options_input =
GeneralUtils::ParseToStr(opt, "/");
for (UInt_t i = 0; i < fit_options_input.size(); i++) {
if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i)
<< "' Provided is not allowed for this measurement."
<< std::endl;
ERR(FTL) << "Fit Options should be provided as a '/' seperated list "
"(e.g. FREE/DIAG/NORM)"
<< std::endl;
ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes
<< "'" << std::endl;
throw;
}
}
// Set TYPE
fFitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos) {
fIsFree = fIsShape = false;
fIsFix = true;
} else if (opt.find("SHAPE") != std::string::npos) {
fIsFree = fIsFix = false;
fIsShape = true;
} else if (opt.find("FREE") != std::string::npos) {
fIsFix = fIsShape = false;
fIsFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos) {
fIsDiag = true;
fIsFull = false;
} else if (opt.find("FULL") != std::string::npos) {
fIsDiag = false;
fIsFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos) {
fIsChi2 = false;
ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl;
ERR(FTL) << "Try to use a chi2!" << std::endl;
throw;
} else {
fIsChi2 = true;
}
// EXTRAS
if (opt.find("RAW") != std::string::npos) fIsRawEvents = true;
if (opt.find("NOWIDTH") != std::string::npos) fIsNoWidth = true;
if (opt.find("DIF") != std::string::npos) fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true;
if (opt.find("NORM") != std::string::npos) fAddNormPen = true;
if (opt.find("MASK") != std::string::npos) fIsMask = true;
return;
};
//********************************************************************
void Measurement1D::SetSmearingMatrix(std::string smearfile, int truedim,
int recodim) {
//********************************************************************
// The smearing matrix describes the migration from true bins (rows) to reco
// bins (columns)
// Counter over the true bins!
int row = 0;
std::string line;
std::ifstream smear(smearfile.c_str(), ifstream::in);
// Note that the smearing matrix may be rectangular.
fSmearMatrix = new TMatrixD(truedim, recodim);
if (smear.is_open())
LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl;
else
ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile
<< std::endl;
while (std::getline(smear >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*fSmearMatrix)(row, column) =
(*iter) / 100.; // Convert to fraction from
// percentage (this may not be
// general enough)
column++;
}
row++;
}
return;
}
//********************************************************************
void Measurement1D::ApplySmearingMatrix() {
//********************************************************************
if (!fSmearMatrix) {
ERR(WRN) << fName
<< ": attempted to apply smearing matrix, but none was set"
<< std::endl;
return;
}
TH1D* unsmeared = (TH1D*)fMCHist->Clone();
TH1D* smeared = (TH1D*)fMCHist->Clone();
smeared->Reset();
// Loop over reconstructed bins
// true = row; reco = column
for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) {
// Sum up the constributions from all true bins
double rBinVal = 0;
// Loop over true bins
for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) {
rBinVal +=
(*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1);
}
smeared->SetBinContent(rbin + 1, rBinVal);
}
fMCHist = (TH1D*)smeared->Clone();
return;
}
/*
Reconfigure LOOP
*/
//********************************************************************
void Measurement1D::ResetAll() {
//********************************************************************
fMCHist->Reset();
fMCFine->Reset();
fMCStat->Reset();
return;
};
//********************************************************************
void Measurement1D::FillHistograms() {
//********************************************************************
if (Signal) {
fMCHist->Fill(fXVar, Weight);
fMCFine->Fill(fXVar, Weight);
fMCStat->Fill(fXVar, 1.0);
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight);
}
return;
};
//********************************************************************
void Measurement1D::ScaleEvents() {
//********************************************************************
// Fill MCWeighted;
// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
// fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
// fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
// }
// Setup Stat ratios for MC and MC Fine
double* statratio = new double[fMCHist->GetNbinsX()];
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
if (fMCHist->GetBinContent(i + 1) != 0) {
statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
} else {
statratio[i] = 0.0;
}
}
double* statratiofine = new double[fMCFine->GetNbinsX()];
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
if (fMCFine->GetBinContent(i + 1) != 0) {
statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
} else {
statratiofine[i] = 0.0;
}
}
// Scaling for raw event rates
if (fIsRawEvents) {
double datamcratio = fDataHist->Integral() / fMCHist->Integral();
fMCHist->Scale(datamcratio);
fMCFine->Scale(datamcratio);
if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio);
// Scaling for XSec as function of Enu
} else if (fIsEnu1D) {
PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
// if (fMCHist_Modes) {
// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
// GetEventHistogram(), fScaleFactor,
// fNEvents);
// }
} else if (fIsNoWidth) {
fMCHist->Scale(fScaleFactor);
fMCFine->Scale(fScaleFactor);
if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor);
// Any other differential scaling
} else {
fMCHist->Scale(fScaleFactor, "width");
fMCFine->Scale(fScaleFactor, "width");
if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width");
}
// Proper error scaling - ROOT Freaks out with xsec weights sometimes
for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
}
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]);
}
// Clean up
delete statratio;
delete statratiofine;
return;
};
//********************************************************************
void Measurement1D::ApplyNormScale(double norm) {
//********************************************************************
fCurrentNorm = norm;
fMCHist->Scale(1.0 / norm);
fMCFine->Scale(1.0 / norm);
return;
};
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int Measurement1D::GetNDOF() {
//********************************************************************
int ndof = fDataHist->GetNbinsX();
if (fMaskHist and fIsMask) ndof -= fMaskHist->Integral();
return ndof;
}
//********************************************************************
double Measurement1D::GetLikelihood() {
//********************************************************************
// If this is for a ratio, there is no data histogram to compare to!
if (fNoData || !fDataHist) return 0.;
// Apply Masking to MC if Required.
if (fIsMask and fMaskHist) {
PlotUtils::MaskBins(fMCHist, fMaskHist);
}
// Sort Shape Scaling
double scaleF = 0.0;
// TODO Include !fIsRawEvents
if (fIsShape) {
if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) {
scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") /
fMCHist->Integral(1, fMCHist->GetNbinsX(), "width");
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
}
}
// Likelihood Calculation
double stat = 0.;
if (fIsChi2) {
if (fIsRawEvents) {
stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist);
} else if (fIsDiag) {
stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist);
} else if (!fIsDiag and !fIsRawEvents) {
stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist);
}
}
// Sort Penalty Terms
if (fAddNormPen) {
double penalty =
(1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError);
stat += penalty;
}
// Return to normal scaling
if (fIsShape) { // and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
fLikelihood = stat;
return stat;
}
/*
Fake Data Functions
*/
//********************************************************************
void Measurement1D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Setup original/datatrue
TH1D* tempdata = (TH1D*) fDataHist->Clone();
if (!fIsFakeData) {
fIsFakeData = true;
// Make a copy of the original data histogram.
if (!fDataOrig) fDataOrig = (TH1D*)fDataHist->Clone((fName + "_data_original").c_str());
} else {
ResetFakeData();
}
// Setup Inputs
fFakeDataInput = fakeOption;
LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl;
// From MC
if (fFakeDataInput.compare("MC") == 0) {
fDataHist = (TH1D*)fMCHist->Clone((fName + "_MC").c_str());
// Fake File
} else {
if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
fDataHist = (TH1D*)fFakeDataFile->Get((fName + "_MC").c_str());
}
// Setup Data Hist
fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
(fName + fPlotTitles).c_str());
// Replace Data True
if (fDataTrue) delete fDataTrue;
fDataTrue = (TH1D*)fDataHist->Clone();
fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
(fName + fPlotTitles).c_str());
// Make a new covariance for fake data hist.
int nbins = fDataHist->GetNbinsX();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
(*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j);
}
}
// Setup Covariances
if (covar) delete covar;
covar = StatUtils::GetInvert(fFullCovar);
if (fDecomp) delete fDecomp;
fDecomp = StatUtils::GetInvert(fFullCovar);
delete tempdata;
return;
};
//********************************************************************
void Measurement1D::ResetFakeData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
}
}
//********************************************************************
void Measurement1D::ResetData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
}
fIsFakeData = false;
}
//********************************************************************
void Measurement1D::ThrowCovariance() {
//********************************************************************
// Take a fDecomposition and use it to throw the current dataset.
// Requires fDataTrue also be set incase used repeatedly.
-
+
if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone();
if (fDataHist) delete fDataHist;
fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
return;
};
//********************************************************************
void Measurement1D::ThrowDataToy(){
-//********************************************************************
+//********************************************************************
if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone();
if (fMCHist) delete fMCHist;
fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
}
/*
Access Functions
*/
//********************************************************************
TH1D* Measurement1D::GetMCHistogram() {
//********************************************************************
if (!fMCHist) return fMCHist;
std::ostringstream chi2;
chi2 << std::setprecision(5) << this->GetLikelihood();
int linecolor = kRed;
int linestyle = 1;
int linewidth = 1;
int fillcolor = 0;
int fillstyle = 1001;
// if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor");
// if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle");
// if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth");
// if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor");
// if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle");
fMCHist->SetTitle(chi2.str().c_str());
fMCHist->SetLineColor(linecolor);
fMCHist->SetLineStyle(linestyle);
fMCHist->SetLineWidth(linewidth);
fMCHist->SetFillColor(fillcolor);
fMCHist->SetFillStyle(fillstyle);
return fMCHist;
};
//********************************************************************
TH1D* Measurement1D::GetDataHistogram() {
//********************************************************************
if (!fDataHist) return fDataHist;
int datacolor = kBlack;
int datastyle = 1;
int datawidth = 1;
// if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor");
// if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle");
// if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth");
fDataHist->SetLineColor(datacolor);
fDataHist->SetLineWidth(datawidth);
fDataHist->SetMarkerStyle(datastyle);
return fDataHist;
};
/*
Write Functions
*/
// Save all the histograms at once
//********************************************************************
void Measurement1D::Write(std::string drawOpt) {
//********************************************************************
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
// Write Settigns
if (drawOpt.find("SETTINGS") != std::string::npos){
fSettings.Set("#chi^{2}",fLikelihood);
fSettings.Set("NDOF", this->GetNDOF() );
fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF() );
fSettings.Write();
}
// Write Data/MC
GetDataList().at(0)->Write();
GetMCList().at(0)->Write();
+ if(fEvtRateScaleFactor != 0xdeadbeef){
+ TH1D * PredictedEvtRate = static_cast(GetMCList().at(0)->Clone());
+ PredictedEvtRate->Scale(fEvtRateScaleFactor);
+ PredictedEvtRate->GetYaxis()->SetTitle("Predicted event rate");
+ PredictedEvtRate->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())
GetXSecHistogram()->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());
TH1D* dataShape = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE").c_str());
if (fShapeCovar) StatUtils::SetDataErrorFromCov(dataShape, fShapeCovar, 1E-38);
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();
dataShape->Write();
delete mcShape;
}
//********************************************************************
void Measurement1D::WriteShapeRatioPlot() {
//********************************************************************
// Get a mcshape histogram
TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
// Create shape ratio histograms
TH1D* mcShapeRatio = (TH1D*)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
TH1D* dataShapeRatio = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());
// Divide the histograms
mcShapeRatio->Divide(mcShape);
dataShapeRatio->Divide(mcShape);
// Colour the shape ratio plots
mcShapeRatio->SetLineWidth(3);
mcShapeRatio->SetLineStyle(7);
mcShapeRatio->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
//// CRAP TO BE REMOVED
//********************************************************************
void Measurement1D::SetupMeasurement(std::string inputfile, std::string type,
FitWeight * rw, std::string fkdt) {
//********************************************************************
nuiskey samplekey = Config::CreateKey("sample");
samplekey.AddS("name", fName);
samplekey.AddS("type",type);
samplekey.AddS("input",inputfile);
fSettings = LoadSampleSettings(samplekey);
-
+
// Reset everything to NULL
// Init();
// Check if name contains Evt, indicating that it is a raw number of events
// measurements and should thus be treated as once
fIsRawEvents = false;
if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!"
<< std::endl;
LOG(SAM) << "Overriding this and setting fIsRawEvents == true!"
<< std::endl;
}
fIsEnu1D = false;
if (fName.find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
"not flux averaged!"
<< std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
fRW = rw;
if (!fInput and !fIsJoint) SetupInputs(inputfile);
// Set Default Options
SetFitOptions(fDefaultTypes);
// Set Passed Options
SetFitOptions(type);
// Still adding support for flat flux inputs
// // Set Enu Flux Scaling
// if (isFlatFluxFolding) this->Input()->ApplyFluxFolding(
// this->defaultFluxHist );
// FinaliseMeasurement();
}
//********************************************************************
void Measurement1D::SetupDefaultHist() {
//********************************************************************
// Setup fMCHist
fMCHist = (TH1D*)fDataHist->Clone();
fMCHist->SetNameTitle((fName + "_MC").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
// Setup fMCFine
Int_t nBins = fMCHist->GetNbinsX();
fMCFine = new TH1D(
(fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(),
nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1));
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
fMCHist->Reset();
fMCFine->Reset();
// Setup the NEUT Mode Array
PlotUtils::CreateNeutModeArray((TH1D*)fMCHist, (TH1**)fMCHist_PDG);
PlotUtils::ResetNeutModeArray((TH1**)fMCHist_PDG);
// Setup bin masks using sample name
if (fIsMask) {
std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
}
SetBinMask(maskloc);
}
fMCHist_Modes = new TrueModeStack( (fName + "_MODES").c_str(), ("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);
return;
}
//********************************************************************
void Measurement1D::SetDataValues(std::string dataFile) {
//********************************************************************
// Override this function if the input file isn't in a suitable format
LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
fDataHist =
PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles);
fDataTrue = (TH1D*)fDataHist->Clone();
// Number of data points is number of bins
fNDataPointsX = fDataHist->GetXaxis()->GetNbins();
return;
};
//********************************************************************
void Measurement1D::SetDataFromDatabase(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(
(GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void Measurement1D::SetDataFromFile(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void Measurement1D::SetCovarMatrix(std::string covarFile) {
//********************************************************************
// Covariance function, only really used when reading in the MB Covariances.
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
TH2D* covarPlot = new TH2D();
// TH2D* decmpPlot = new TH2D();
TH2D* covarInvPlot = new TH2D();
TH2D* fFullCovarPlot = new TH2D();
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("thrown_covariance");
if (fIsShape || fIsFree) covName = "shp_";
if (fIsDiag)
covName += "diag";
else
covName += "full";
covarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").c_str());
if (!covOption.compare("SUB"))
fFullCovarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
else if (!covOption.compare("FULL"))
fFullCovarPlot = (TH2D*)tempFile->Get("fullcov");
else
ERR(WRN) << "Incorrect thrown_covariance option in parameters."
<< std::endl;
int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
int covdim = int(fDataHist->GetNbinsX());
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
fDecomp = new TMatrixDSym(dim);
int row, column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++) {
// if (this->masked->GetBinContent(i+1) > 0) continue;
for (Int_t j = 0; j < covdim; j++) {
// if (this->masked->GetBinContent(j+1) > 0) continue;
(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);
column++;
}
column = 0;
row++;
}
// Set bin errors on data
if (!fIsDiag) {
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar);
}
// Get Deteriminant and inverse matrix
// fCovDet = this->covar->Determinant();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// Sets the covariance matrix from a provided file in a text format
// scale is a multiplicative pre-factor to apply in the case where the
// covariance is given in some unit (e.g. 1E-38)
void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim,
double scale) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream covarread(covarFile.c_str(), ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covarread.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
else
ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile
<< std::endl;
// Loop over the lines in the file
while (std::getline(covarread >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
std::vector 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::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 entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::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 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& cont,
// std::vector& 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/MeasurementBase.cxx b/src/FitBase/MeasurementBase.cxx
index 1ac6a92..dbfcd0a 100644
--- a/src/FitBase/MeasurementBase.cxx
+++ b/src/FitBase/MeasurementBase.cxx
@@ -1,563 +1,596 @@
// 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 .
*******************************************************************************/
#include "MeasurementBase.h"
/*
Constructor/Destructors
*/
//********************************************************************
// 2nd Level Constructor (Inherits From MeasurementBase.h)
MeasurementBase::MeasurementBase(void) {
//********************************************************************
fScaleFactor = 1.0;
fMCFilled = false;
fNoData = false;
fInput = NULL;
NSignal = 0;
// Set the default values
// After-wards this gets set in SetupMeasurement
EnuMin = 0.;
EnuMax = 1.E5;
fMeasurementSpeciesType = kSingleSpeciesMeasurement;
fEventVariables = NULL;
fIsJoint = false;
+
+ fNPOT = 0xdeadbeef;
+ fFluxIntegralOverride = 0xdeadbeef;
+ fTargetVolume = 0xdeadbeef;
+ fTargetMaterialDensity = 0xdeadbeef;
+ fEvtRateScaleFactor = 0xdeadbeef;
};
void MeasurementBase::FinaliseMeasurement() {
-
// Used to setup default data hists, covars, etc.
-
-
-
}
//********************************************************************
// 2nd Level Destructor (Inherits From MeasurementBase.h)
-MeasurementBase::~MeasurementBase() {
- //********************************************************************
+MeasurementBase::~MeasurementBase(){
+ //********************************************************************
};
//********************************************************************
double MeasurementBase::TotalIntegratedFlux(std::string intOpt, double low,
- double high) {
-//********************************************************************
+ 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) {
+ 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 file_descriptor =
- GeneralUtils::ParseToStr(inputfile, ":");
+ 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]);
+ InputUtils::ParseInputType(file_descriptor[0]);
fInput = InputUtils::CreateInputHandler(fName, inpType, file_descriptor[1]);
}
-
fNEvents = fInput->GetNEvents();
// Expect INPUTTYPE:FileLocation(s)
std::vector file_descriptor =
- GeneralUtils::ParseToStr(inputfile, ":");
+ 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);
+ 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;
- if (!fIsJoint) SetupInputs( setting.GetS("input") );
+
+ fEvtRateScaleFactor = 0xdeadbeef;
+
+ if (!fIsJoint) {
+ SetupInputs(setting.GetS("input"));
+
+ fNPOT = samplekey.Has("NPOT") ? samplekey.GetD("NPOT") : 1;
+ fFluxIntegralOverride = samplekey.Has("FluxIntegralOverride")
+ ? samplekey.GetD("FluxIntegralOverride")
+ : 0xdeadbeef;
+ fTargetVolume = samplekey.Has("TargetVolume")
+ ? samplekey.GetD("TargetVolume")
+ : 0xdeadbeef;
+ fTargetMaterialDensity = samplekey.Has("TargetMaterialDensity")
+ ? samplekey.GetD("TargetMaterialDensity")
+ : 0xdeadbeef;
+ if ((fTargetVolume != 0xdeadbeef) &&
+ (fTargetMaterialDensity != 0xdeadbeef)) {
+ double TargetMass_kg = fTargetVolume * fTargetMaterialDensity;
+ double NNucleons = TargetMass_kg / PhysConst::mass_nucleon_kg;
+ double NNeutrinos =
+ ((fFluxIntegralOverride == 0xdeadbeef) ? TotalIntegratedFlux()
+ : fFluxIntegralOverride) *
+ fNPOT;
+ fEvtRateScaleFactor = NNeutrinos * NNucleons;
+
+ QLOG(SAM, "\tEvent rate prediction : ");
+ QLOG(SAM, "\t\tTarget volume : " << fTargetVolume << " m^3");
+ QLOG(SAM, "\t\tTarget density : " << fTargetMaterialDensity << " kg/m^3");
+ QLOG(SAM, "\t\tTarget mass : " << TargetMass_kg << " kg");
+ QLOG(SAM, "\t\tNTarget Nucleons : " << NNucleons);
+ QLOG(SAM, "\t\tNNeutrinos : " << NNeutrinos
+ << ((fNPOT == 1) ? "/cm^2" : "/POT /cm^2"));
+ QLOG(SAM, "\t\tXSec -> EvtRate scale factor : " << fEvtRateScaleFactor);
+ }
+ }
return setting;
}
//***********************************************
-SampleSettings MeasurementBase::LoadSampleSettings(std::string name, std::string input, std::string type) {
-//***********************************************
+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);
+ 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");
-
}
-
//***********************************************
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){
-
+ 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(std::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;
+ 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();
+void MeasurementBase::FillHistogramsFromBox(MeasurementVariableBox* var,
+ double weight) {
+ fXVar = var->GetX();
+ fYVar = var->GetY();
+ fZVar = var->GetZ();
// Signal = var->fSignal;
// Mode = var->fMode;
Weight = weight;
fEventVariables = var;
FillHistograms();
FillExtraHistograms(var, weight);
-
}
-void MeasurementBase::FillHistograms(double weight){
+void MeasurementBase::FillHistograms(double weight) {
Weight = weight * GetBox()->GetSampleWeight();
FillHistograms();
FillExtraHistograms(GetBox(), Weight);
}
-
MeasurementVariableBox* MeasurementBase::FillVariableBox(FitEvent* event) {
-
GetBox()->Reset();
Mode = event->Mode;
- Weight = 1.0; //event->Weight;
+ Weight = 1.0; // event->Weight;
this->FillEventVariables(event);
Signal = this->isSignal(event);
GetBox()->FillBoxFromEvent(event);
GetBox()->SetX(fXVar);
GetBox()->SetY(fYVar);
GetBox()->SetZ(fZVar);
GetBox()->SetMode(event->Mode);
GetBox()->SetSampleWeight(Weight);
// 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;
+ 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"
+ "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 MeasurementBase::GetFluxList() {
//***********************************************
return GetInput()->GetFluxList();
}
//***********************************************
std::vector MeasurementBase::GetEventRateList() {
//***********************************************
return GetInput()->GetEventList();
}
//***********************************************
std::vector MeasurementBase::GetXSecList() {
//***********************************************
return GetInput()->GetXSecList();
}
-
void MeasurementBase::ProcessExtraHistograms(int cmd,
- MeasurementVariableBox* vars,
- double weight) {
+ 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) {
+ 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) {
+ 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) {
+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!
+ SetAutoProcessTH1(fake, c1, c2, c3, c4,
+ c5); // Need to add a destroy command!
}
-void MeasurementBase::SetAutoProcess(TH1* hist, int c1, int c2, int c3, int c4, int c5) {
+void MeasurementBase::SetAutoProcess(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!
+ SetAutoProcessTH1(fake, c1, c2, c3, c4,
+ c5); // Need to add a destroy command!
}
-void MeasurementBase::SetAutoProcess(TGraph* g, int c1, int c2, int c3, int c4, int c5) {
+void MeasurementBase::SetAutoProcess(TGraph* g, int c1, int c2, int c3, int c4,
+ int c5) {
FakeStack* fake = new FakeStack(g);
- SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command!
+ SetAutoProcessTH1(fake, c1, c2, c3, c4,
+ c5); // Need to add a destroy command!
}
-void MeasurementBase::SetAutoProcess(TF1* f, int c1, int c2, int c3, int c4, int c5) {
+void MeasurementBase::SetAutoProcess(TF1* f, int c1, int c2, int c3, int c4,
+ int c5) {
FakeStack* fake = new FakeStack(f);
- SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command!
+ SetAutoProcessTH1(fake, c1, c2, c3, c4,
+ c5); // Need to add a destroy command!
}
-void MeasurementBase::SetAutoProcess(StackBase* hist, int c1, int c2, int c3, int c4, int c5){
+void MeasurementBase::SetAutoProcess(StackBase* hist, int c1, int c2, int c3,
+ int c4, int c5) {
SetAutoProcessTH1(hist, c1, c2, c3, c4, c5);
}
-void MeasurementBase::SetAutoProcessTH1(StackBase* hist, int c1, int c2, int c3, int c4, int c5) {
-
+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(5,0);
+ fExtraTH1s[hist] = std::vector(5, 0);
// Setup a default one.
- if (c1 == -1 && c2 == -1 && c3 == -1 && c4 == -1 && c5 == -1){
+ 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;
+ // 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 >::iterator iter = fExtraTH1s.begin();
+ for (std::map >::iterator iter =
+ fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
-
if (!((*iter).second)[kCMD_Reset]) continue;
(*iter).first->Reset();
}
};
void MeasurementBase::AutoScaleExtraTH1() {
- for (std::map >::iterator iter = fExtraTH1s.begin();
+ for (std::map >::iterator iter =
+ fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
-
if (!((*iter).second)[kCMD_Scale]) continue;
- if (fIsNoWidth){
+ if (fIsNoWidth) {
(*iter).first->Scale(fScaleFactor);
} else {
(*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 >::iterator iter = fExtraTH1s.begin();
+ for (std::map >::iterator iter =
+ fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
-
if (!((*iter).second)[kCMD_Norm]) continue;
(*iter).first->Scale(sfactor);
}
};
void MeasurementBase::AutoWriteExtraTH1() {
- for (std::map >::iterator iter = fExtraTH1s.begin();
+ for (std::map >::iterator iter =
+ fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
-
if (!(((*iter).second)[kCMD_Write])) continue;
(*iter).first->Write();
}
};
-
-
-
diff --git a/src/FitBase/MeasurementBase.h b/src/FitBase/MeasurementBase.h
index ea2e1e5..33a0ef6 100644
--- a/src/FitBase/MeasurementBase.h
+++ b/src/FitBase/MeasurementBase.h
@@ -1,357 +1,362 @@
// 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 .
*******************************************************************************/
#ifndef INPUTHANDLER_H_SEEN
#define INPUTHANDLER_H_SEEN
/*!
* \addtogroup FitBase
* @{
*/
// C Includes
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
// ROOT includes
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
// External data fit includes
#include "FitEvent.h"
#include "FitParameters.h"
#include "FitUtils.h"
#include "GeneralUtils.h"
#include "PlotUtils.h"
#include "StatUtils.h"
#include "InputFactory.h"
#include "FitWeight.h"
#include "TMultiDimFit.h"
#ifdef __GENIE_ENABLED__
#include "Conventions/Units.h"
#endif
#include "EventManager.h"
#include "TObject.h"
#include "InputHandler.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#include "SampleSettings.h"
#include "StackBase.h"
#include "StandardStacks.h"
/// Enumerations to help with extra plot functions
enum extraplotflags {
kCMD_Reset = 0,
kCMD_Fill,
kCMD_Scale,
kCMD_Norm,
kCMD_Write,
kCMD_Error,
kCMD_extraplotflags
};
enum MeasurementSpeciesClass {
kSingleSpeciesMeasurement = 0,
kNumuWithWrongSignMeasurement,
kNueWithWrongSignMeasurement,
kFourSpeciesMeasurement,
};
/// InputHandler Class
///
/// Inherits from Measurement base to handle whatever input is throwna t the
/// fitter automatically.
/// All functions here handle how files are read in, converted to custom formats
/// and reconfigures are called.
/// Used generally for the MC inputs.
//! 2nd level experiment class that handles converting MC into a common format
//! and calling reconfigure
class MeasurementBase {
public:
/*
Constructor/Destructors
*/
//! Default Constructor. Set everything to NULL
MeasurementBase();
//! Default virtual destructor
virtual ~MeasurementBase(void);
virtual void InitialSetup(void) {};
/*
Reconfigure Functions
*/
//! Function called if MC tuning dials haven't been changed and all we want to
//! do is update the normalisation.
virtual void Renormalise(void);
//! Call reconfigure only looping over signal events to save time.
virtual void ReconfigureFast(void);
virtual void FillHistograms(double weight);
//! Call reconfigure looping over all MC events including background
virtual void Reconfigure(void);
// virtual TH2D GetCovarMatrix(void) = 0;
virtual double GetLikelihood(void) { return 0.0; };
virtual int GetNDOF(void) { return 0; };
virtual void ThrowCovariance(void) = 0;
virtual void ThrowDataToy(void) = 0;
virtual void SetFakeDataValues(std::string fkdt) = 0;
//! Get the total integrated flux between this samples energy range
virtual double TotalIntegratedFlux(std::string intOpt = "width",
double low = -9999.9,
double high = -9999.9);
//! Get the predicted event rate for this sample
virtual double PredictedEventRate(std::string intOpt = "width",
double low = -9999.9,
double high = -9999.9);
virtual SampleSettings LoadSampleSettings(nuiskey samplekey);
virtual SampleSettings LoadSampleSettings(std::string name, std::string input, std::string type);
virtual void FinaliseSampleSettings();
virtual void FinaliseMeasurement();
virtual void ProcessExtraHistograms(int cmd, MeasurementVariableBox* vars,
double weight = 1.0);
virtual void FillExtraHistograms(MeasurementVariableBox* vars, double weight = 1.0);
virtual void ScaleExtraHistograms(MeasurementVariableBox* vars);
virtual void ResetExtraHistograms();
virtual void NormExtraHistograms(MeasurementVariableBox* vars, double norm = 1.0);
virtual void WriteExtraHistograms();
virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox();};
int GetPassed() {
int signalSize = 0;
return signalSize;
}
int GetTotal() { return fNEvents; }
/*
Reconfigure LOOP
*/
// All these should be virtual
///! Reset Histograms (Handled at Measurement Stage)
virtual void ResetAll(void) = 0;
///! Fill the event variables for this sample (Handled in each inherited
/// sample)
virtual void FillEventVariables(FitEvent* event) { (void)event; };
///! Check whether this event is signle (Handled in each inherited sample)
virtual bool isSignal(FitEvent* event) {
(void)event;
return false;
};
///! Fill the histogram for this event using fXVar and fYVar (Handled in each
/// inherited sample)
virtual void FillHistograms(void) {};
///! Convert event rates to whatever distributions you need.
virtual void ConvertEventRates(void);
///! Call scale events after the plots have been filled at the end of
/// reconfigure.
virtual void ScaleEvents(void) {};
///! Apply the scale factor at the end of reconfigure.
virtual void ApplyNormScale(double norm) { (void)norm; };
///! Save Histograms
virtual void Write(std::string drawOpt = "") = 0;
virtual MeasurementVariableBox* FillVariableBox(FitEvent* event);
virtual MeasurementVariableBox* GetBox();
void FillHistogramsFromBox(MeasurementVariableBox* var, double weight);
/*
Histogram Access Functions
*/
///! Virtual function to get data histogram
virtual std::vector GetDataList(void) = 0;
///! Virtual function to get MC histogram
virtual std::vector GetMCList(void) = 0;
virtual std::vector GetFineList(void) = 0;
virtual std::vector GetMaskList(void) = 0;
///! Return flux histograms in a vector
virtual std::vector GetFluxList(void);
virtual std::vector GetEventRateList(void);
virtual std::vector GetXSecList(void);
virtual TH1D* GetEventHistogram() { return fInput->GetEventHistogram(); };
virtual TH1D* GetXSecHistogram() { return fInput->GetXSecHistogram(); };
virtual TH1D* GetFluxHistogram() { return fInput->GetFluxHistogram(); };
///! Return input for this sample
InputHandlerBase* GetInput(void);
std::string GetName(void) { return fName; };
double GetScaleFactor(void) { return fScaleFactor; };
double GetXVar(void) { return fXVar; };
double GetYVar(void) { return fYVar; };
double GetZVar(void) { return fZVar; };
double GetMode(void) { return this->Mode; };
double GetEnu(void) { return this->Enu; };
void SetupInputs(std::string inputfile);
int GetInputID(void);
std::string GetInputFileName() { return fInputFileName; };
void SetSignal(bool sig);
void SetSignal(FitEvent* evt);
void SetWeight(double wght);
void SetMode(int md);
void SetNoData(bool isTrue = true) { fNoData = isTrue; };
inline void SetXVar(double xvar) { fXVar = xvar; };
inline void SetYVar(double yvar) { fYVar = yvar; };
inline void SetZVar(double zvar) { fZVar = zvar; };
virtual std::vector GetSubSamples() {
return std::vector(1, this);
}
void SetAutoProcessTH1(TH1* hist, int c1 = -1,
int c2 = -1, int c3 = -1,
int c4 = -1, int c5 = -1);
void SetAutoProcess(TH1* hist, int c1 = -1,
int c2 = -1, int c3 = -1,
int c4 = -1, int c5 = -1);
void SetAutoProcess(TGraph* g, int c1 = -1,
int c2 = -1, int c3 = -1,
int c4 = -1, int c5 = -1);
void SetAutoProcess(TF1* f, int c1 = -1,
int c2 = -1, int c3 = -1,
int c4 = -1, int c5 = -1);
void SetAutoProcess(StackBase* hist, int c1 = -1,
int c2 = -1, int c3 = -1,
int c4 = -1, int c5 = -1);
void SetAutoProcessTH1(StackBase* hist, int c1 = -1,
int c2 = -1, int c3 = -1,
int c4 = -1, int c5 = -1);
void AutoFillExtraTH1();
void AutoResetExtraTH1();
void AutoScaleExtraTH1();
void AutoNormExtraTH1(double norm);
void AutoWriteExtraTH1();
// functions that need to be added.
// - Initial Check
// - Check Target/Beam loop.
// - Check flux shape if suggested one given.
// - Return MeasurementList (returns )
protected:
// Minimum and maximum energies
double Enu; //!< Neutrino Energy
double EnuMin; //!< Minimum incoming particle energy of events to include
double EnuMax; //!< Maximum incoming particle energy of events to include
BaseFitEvt* signal_event;
FitEvent* cust_event;
FitWeight* fRW; //!< Pointer to the rw engine
InputHandlerBase* fInput; //!< Instance of the input handler
std::string fName; //!< Name of the sample
int fEventType;
double fBeamDistance; //!< Incoming Particle flight distance (for oscillation
//! analysis)
double fScaleFactor; //!< fScaleFactor applied to events to convert from
//! eventrate to final distribution
double
fCurrentNorm; //!< current normalisation factor applied if fit is "FREE"
bool fMCFilled; //!< flag whether MC plots have been filled (For
//! ApplyNormalisation)
bool fNoData; //!< flag whether data plots do not exist (for ratios)
bool fIsNoWidth; ///< Flag : Don't scale by bin width
// TEMP OBJECTS TO HANDLE MERGE
double fXVar, fYVar, fZVar, Mode, Weight;
bool Signal;
int ievt;
int fNEvents;
double Enu_rec, ThetaMu, CosThetaMu;
InputUtils::InputType fInputType;
std::string fInputFileName;
TH1D* fFluxHist;
TH1D* fEventHist;
MeasurementSpeciesClass fMeasurementSpeciesType;
SampleSettings fSettings;
MeasurementVariableBox* fEventVariables;
std::map > fExtraTH1s;
int NSignal;
// std::map fExtaStacks;
bool fIsJoint;
+
+
+ double fNPOT, fFluxIntegralOverride, fTargetVolume, fTargetMaterialDensity;
+ double fEvtRateScaleFactor;
+
};
// Class TypeDefs
typedef std::list::const_iterator MeasListConstIter;
typedef std::list::iterator MeasListIter;
typedef std::vector::const_iterator MeasVectConstIter;
typedef std::vector::iterator MeasVectIter;
/*! @} */
#endif
diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
index d2fc662..60818db 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
@@ -1,379 +1,380 @@
// 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 .
*******************************************************************************/
#include "Smear_SVDUnfold_Propagation_Osc.h"
#include "SmearceptanceUtils.h"
#include "OscWeightEngine.h"
//********************************************************************
Smear_SVDUnfold_Propagation_Osc::Smear_SVDUnfold_Propagation_Osc(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip =
"Simple measurement class for doing fake data oscillation studies.\n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("Osc Studies");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("XXX");
fSettings.SetYTitle("Number of events");
fSettings.SetEnuRange(0.0, 50);
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.DefineAllowedTargets("*");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") / (fNEvents + 0.);
// Plot Setup -------------------------------------------------------
// Check that we have the relevant histograms specified.
if (!samplekey.Has("NDDataHist") || !samplekey.Has("FDDataHist") ||
!samplekey.Has("NDToSpectrumSmearingMatrix")) {
THROW(
"Expected to attributes named: NDDataHist, FDDataHist, "
"NDToSpectrumSmearingMatrix on the Smear_SVDUnfold_Propagation_Osc "
"sample "
"tag.");
}
NDDataHist = NULL;
FDDataHist = NULL;
NDToSpectrumSmearingMatrix = NULL;
fMCHist = NULL;
fDataHist = NULL;
std::vector NDHistDescriptor =
GeneralUtils::ParseToStr(samplekey.GetS("NDDataHist"), ",");
if (NDHistDescriptor.size() == 2) {
NDDataHist = PlotUtils::GetTH1DFromRootFile(NDHistDescriptor[0],
NDHistDescriptor[1]);
ND_True_Spectrum_Hist =
PlotUtils::GetTH1DFromRootFile(NDHistDescriptor[0], "ELep_rate");
}
if (!NDDataHist) {
THROW("Attempted to load NDDataHist from the descriptor: \""
<< samplekey.GetS("NDDataHist")
<< "\", but failed. Does it look correct? \",\"");
}
std::vector FDHistDescriptor =
GeneralUtils::ParseToStr(samplekey.GetS("FDDataHist"), ",");
if (FDHistDescriptor.size() == 2) {
FDDataHist = PlotUtils::GetTH1DFromRootFile(FDHistDescriptor[0],
FDHistDescriptor[1]);
FD_True_Spectrum_Hist =
PlotUtils::GetTH1DFromRootFile(FDHistDescriptor[0], "ELep_rate");
}
if (!FDDataHist) {
THROW("Attempted to load FDDataHist from the descriptor: \""
<< samplekey.GetS("FDDataHist")
<< "\", but failed. Does it look correct? \",\"");
}
std::vector NDToEvSmearingDescriptor = GeneralUtils::ParseToStr(
samplekey.GetS("NDToSpectrumSmearingMatrix"), ",");
if (NDToEvSmearingDescriptor.size() == 2) {
NDToSpectrumSmearingMatrix = PlotUtils::GetTH2DFromRootFile(
NDToEvSmearingDescriptor[0], NDToEvSmearingDescriptor[1]);
}
if (!NDToSpectrumSmearingMatrix) {
THROW("Attempted to load NDToSpectrumSmearingMatrix from the descriptor: \""
<< samplekey.GetS("NDToSpectrumSmearingMatrix")
<< "\", but failed. Does it look correct? \",\"");
}
TMatrixD NDToSpectrumResponseMatrix_l = SmearceptanceUtils::GetMatrix(
SmearceptanceUtils::SVDGetInverse(NDToSpectrumSmearingMatrix));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
TMatrixD SpectrumToFDResponseMatrix_l =
SmearceptanceUtils::GetMatrix(NDToSpectrumSmearingMatrix);
SpectrumToFDResponseMatrix.ResizeTo(SpectrumToFDResponseMatrix_l);
SpectrumToFDResponseMatrix = SpectrumToFDResponseMatrix_l;
fDataHist = static_cast(FDDataHist->Clone());
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist = static_cast(FDDataHist->Clone());
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
SetCovarFromDiagonal();
TruncateUpTo = 0;
if (samplekey.Has("TruncateUpTo")) {
TruncateUpTo = samplekey.GetI("TruncateUpTo");
QLOG(SAM, "Allowed to truncate unfolding matrix by up to "
<< TruncateUpTo
<< " singular values to limit negative ENu spectra.");
}
if (samplekey.Has("FitRegion_Min")) {
FitRegion_Min = samplekey.GetD("FitRegion_Min");
} else {
FitRegion_Min = 0xdeadbeef;
}
if (samplekey.Has("FitRegion_Max")) {
FitRegion_Max = samplekey.GetD("FitRegion_Max");
} else {
FitRegion_Max = 0xdeadbeef;
}
//----------------------------
// Mask data hist if needed
fDataHist->SetBinContent(0, 0);
fDataHist->SetBinError(0, 0);
fDataHist->SetBinContent(fDataHist->GetXaxis()->GetNbins() + 1, 0);
fDataHist->SetBinError(fDataHist->GetXaxis()->GetNbins() + 1, 0);
for (Int_t bi_it = 1; bi_it < fDataHist->GetXaxis()->GetNbins() + 1;
++bi_it) {
if ((FitRegion_Min != 0xdeadbeef) &&
(fDataHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
fDataHist->SetBinContent(bi_it, 0);
fDataHist->SetBinError(bi_it, 0);
}
if ((FitRegion_Max != 0xdeadbeef) &&
(fDataHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
fDataHist->SetBinContent(bi_it, 0);
fDataHist->SetBinError(bi_it, 0);
}
}
TruncateStart = 0;
if (Config::Get().GetConfigNode("smear.SVD.truncation")) {
TruncateStart = Config::Get().ConfI("smear.SVD.truncation");
TMatrixD NDToSpectrumResponseMatrix_l =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
NDToSpectrumSmearingMatrix, TruncateStart));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
}
if (TruncateStart >= TruncateUpTo) {
TruncateUpTo = TruncateStart + 1;
}
NDFDRatio = 1;
- if (Config::Get().GetConfigNode("Osc.NDFDRatio")) {
- NDFDRatio = Config::Get().ConfD("Osc.NDFDRatio");
+ if (samplekey.Has("FDNDRatio")) {
+ NDFDRatio = samplekey.GetD("FDNDRatio");
}
UnfoldToNDETrueSpectrum();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void Smear_SVDUnfold_Propagation_Osc::FillEventVariables(FitEvent *event){};
bool Smear_SVDUnfold_Propagation_Osc::isSignal(FitEvent *event) {
return false;
}
void Smear_SVDUnfold_Propagation_Osc::UnfoldToNDETrueSpectrum(void) {
ND_Unfolded_Spectrum_Hist = NDToSpectrumSmearingMatrix->ProjectionY();
ND_Unfolded_Spectrum_Hist->Clear();
ND_Unfolded_Spectrum_Hist->SetName("ND_Unfolded_Spectrum_Hist");
bool HasNegValue = false;
int truncations = TruncateStart;
do {
if (truncations >= TruncateUpTo) {
THROW("Unfolded enu spectrum had negative values even after "
<< truncations << " SVD singular value truncations.");
}
// Unfold ND ERec -> Enu spectrum
// ------------------------------
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
NDDataHist, ND_Unfolded_Spectrum_Hist, NDToSpectrumResponseMatrix, 1000,
false);
HasNegValue = false;
for (Int_t bi_it = 1;
bi_it < ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
++bi_it) {
if ((FitRegion_Min != 0xdeadbeef) &&
(ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(bi_it) <=
FitRegion_Min)) {
continue;
}
if ((FitRegion_Max != 0xdeadbeef) &&
(ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(bi_it) >
FitRegion_Max)) {
continue;
}
if (ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) < 0) {
HasNegValue = true;
break;
}
}
if (HasNegValue) {
TMatrixD NDToSpectrumResponseMatrix_l =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
NDToSpectrumSmearingMatrix, truncations));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
}
truncations++;
} while (HasNegValue);
NDFD_Corrected_Spectrum_Hist =
static_cast(ND_Unfolded_Spectrum_Hist->Clone());
NDFD_Corrected_Spectrum_Hist->Clear();
NDFD_Corrected_Spectrum_Hist->SetName("NDFD_Corrected_Spectrum_Hist");
FD_Propagated_Spectrum_Hist =
static_cast(ND_Unfolded_Spectrum_Hist->Clone());
FD_Propagated_Spectrum_Hist->Clear();
FD_Propagated_Spectrum_Hist->SetName("FD_Propagated_Spectrum_Hist");
// Apply FD/ND weights
// ------------------------------
for (Int_t bi_it = 1;
bi_it < ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins() + 1; ++bi_it) {
NDFD_Corrected_Spectrum_Hist->SetBinContent(
bi_it, ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) * NDFDRatio);
NDFD_Corrected_Spectrum_Hist->SetBinError(
bi_it, ND_Unfolded_Spectrum_Hist->GetBinError(bi_it) * NDFDRatio);
}
}
void Smear_SVDUnfold_Propagation_Osc::ConvertEventRates(void) {
// Apply Oscillations
// ------------------------------
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
for (Int_t bi_it = 1;
bi_it < NDFD_Corrected_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
++bi_it) {
double oscWeight = oscWE->CalcWeight(
NDFD_Corrected_Spectrum_Hist->GetXaxis()->GetBinCenter(bi_it), 14);
FD_Propagated_Spectrum_Hist->SetBinContent(
bi_it, NDFD_Corrected_Spectrum_Hist->GetBinContent(bi_it) * oscWeight);
FD_Propagated_Spectrum_Hist->SetBinError(
bi_it, NDFD_Corrected_Spectrum_Hist->GetBinError(bi_it) * oscWeight);
}
// Forward fold Spectrum -> ERec FD
// ------------------------------
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
FD_Propagated_Spectrum_Hist, fMCHist, SpectrumToFDResponseMatrix, 1000,
true);
fMCHist->SetBinContent(0, 0);
fMCHist->SetBinError(0, 0);
fMCHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
fMCHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
for (Int_t bi_it = 1; bi_it < fMCHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
if ((FitRegion_Min != 0xdeadbeef) &&
(fMCHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
fMCHist->SetBinContent(bi_it, 0);
fMCHist->SetBinError(bi_it, 0);
}
if ((FitRegion_Max != 0xdeadbeef) &&
(fMCHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
fMCHist->SetBinContent(bi_it, 0);
fMCHist->SetBinError(bi_it, 0);
}
}
}
void Smear_SVDUnfold_Propagation_Osc::Write(std::string drawOpt) {
NDToSpectrumSmearingMatrix->Write("SmearingMatrix_ND", TObject::kOverwrite);
NDDataHist->Write("Obs_ND", TObject::kOverwrite);
ND_Unfolded_Spectrum_Hist->Write(ND_Unfolded_Spectrum_Hist->GetName(),
TObject::kOverwrite);
NDFD_Corrected_Spectrum_Hist->Write(NDFD_Corrected_Spectrum_Hist->GetName(),
TObject::kOverwrite);
FD_Propagated_Spectrum_Hist->Write(FD_Propagated_Spectrum_Hist->GetName(),
TObject::kOverwrite);
fMCHist->Write("Pred_FD", TObject::kOverwrite);
fDataHist->Write("Obs_FD", TObject::kOverwrite);
// Apply Oscillations
// ------------------------------
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
TGraph POsc;
POsc.Set(1E4 - 1);
double min = ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1);
double step = (ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(
ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins()) -
ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1)) /
double(1E4);
for (size_t i = 1; i < 1E4; ++i) {
double enu = min + i * step;
double ow = oscWE->CalcWeight(enu, 14);
if (ow != ow) {
std::cout << "Bad osc weight for ENu: " << enu << std::endl;
}
POsc.SetPoint(i - 1, enu, ow);
}
POsc.Write("POsc", TObject::kOverwrite);
+ Measurement1D::Write(drawOpt);
}
diff --git a/src/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx
index 28f498b..f93c17e 100644
--- a/src/MCStudies/Smearceptance_Tester.cxx
+++ b/src/MCStudies/Smearceptance_Tester.cxx
@@ -1,690 +1,711 @@
// 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 .
*******************************************************************************/
#include "Smearceptance_Tester.h"
#include "SmearceptanceUtils.h"
#include "Smearcepterton.h"
#define DEBUG_SMEARTESTER
//********************************************************************
/// @brief Class to perform smearceptance MC Studies on a custom measurement
-Smearceptance_Tester::Smearceptance_Tester(std::string name,
- std::string inputfile, FitWeight *rw,
- std::string type,
- std::string fakeDataFile) {
+Smearceptance_Tester::Smearceptance_Tester(nuiskey samplekey) {
//********************************************************************
- // Measurement Details
- std::vector splitName = GeneralUtils::ParseToStr(name, "_");
- size_t firstUS = name.find_first_of("_");
-
- std::string smearceptorName = name.substr(firstUS + 1);
-
- fName = name.substr(0, firstUS);
- eventVariables = NULL;
-
- // Define our energy range for flux calcs
- EnuMin = 0.;
- EnuMax = 100.; // Arbritrarily high energy limit
+ // Sample overview ---------------------------------------------------
+ std::string descrip =
+ "Simple measurement class for producing an event summary tree of smeared "
+ "events.\n";
- // Set default fitter flags
- fIsDiag = true;
- fIsShape = false;
- fIsRawEvents = false;
+ if (Config::Get().GetConfigNode("NPOT")) {
+ samplekey.SetS("NPOT", Config::Get().ConfS("NPOT"));
+ }
+ if (Config::Get().GetConfigNode("FluxIntegralOverride")) {
+ samplekey.SetS("FluxIntegralOverride",
+ Config::Get().ConfS("FluxIntegralOverride"));
+ }
+ if (Config::Get().GetConfigNode("TargetVolume")) {
+ samplekey.SetS("TargetVolume", Config::Get().ConfS("TargetVolume"));
+ }
+ if (Config::Get().GetConfigNode("TargetMaterialDensity")) {
+ samplekey.SetS("TargetMaterialDensity",
+ Config::Get().ConfS("TargetMaterialDensity"));
+ }
+ // Setup common settings
+ fSettings = LoadSampleSettings(samplekey);
+
+ fSettings.SetTitle("Smearceptance Studies");
+ fSettings.SetDescription(descrip);
+ fSettings.SetXTitle("XXX");
+ fSettings.SetYTitle("Number of events");
+ fSettings.SetEnuRange(0.0, 1E5);
+ fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
+ fSettings.DefineAllowedTargets("*");
+ fSettings.DefineAllowedSpecies("*");
+
+ FinaliseSampleSettings();
+
+ // Scaling Setup ---------------------------------------------------
+ // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
+ fScaleFactor =
+ (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
+ TotalIntegratedFlux();
- // This function will sort out the input files automatically and parse all the
- // inputs,flags,etc.
- // There may be complex cases where you have to do this by hand, but usually
- // this will do.
- Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
+ // Measurement Details
+ std::vector splitName = GeneralUtils::ParseToStr(fName, "_");
+ size_t firstUS = fName.find_first_of("_");
- eventVariables = NULL;
+ std::string smearceptorName = fName.substr(firstUS + 1);
- // Setup fDataHist as a placeholder
- this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
- this->SetupDefaultHist();
+ fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
+ SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
- this->fScaleFactor =
- (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
- this->TotalIntegratedFlux();
+ eventVariables = NULL;
LOG(SAM) << "Smearceptance Flux Scaling Factor = " << fScaleFactor
<< std::endl;
if (fScaleFactor <= 0.0) {
ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
sleep(20);
}
// Setup our TTrees
- this->AddEventVariablesToTree();
+ AddEventVariablesToTree();
smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName);
Int_t RecNBins = 20, TrueNBins = 20;
double RecBinL = 0, TrueBinL = 0, RecBinH = 10, TrueBinH = 10;
if (Config::Get().GetConfigNode("smear.reconstructed.binning")) {
std::vector args = GeneralUtils::ParseToStr(
Config::Get().ConfS("smear.reconstructed.binning"), ",");
RecNBins = GeneralUtils::StrToInt(args[0]);
RecBinL = GeneralUtils::StrToDbl(args[1]);
RecBinH = GeneralUtils::StrToDbl(args[2]);
TrueNBins = RecNBins;
TrueBinL = RecBinL;
TrueBinH = RecBinH;
}
if (Config::Get().GetConfigNode("smear.true.binning")) {
std::vector args = GeneralUtils::ParseToStr(
Config::Get().ConfS("smear.true.binning"), ",");
TrueNBins = GeneralUtils::StrToInt(args[0]);
TrueBinL = GeneralUtils::StrToDbl(args[1]);
TrueBinH = GeneralUtils::StrToDbl(args[2]);
}
SVDTruncation = 0;
if (Config::Get().GetConfigNode("smear.true.binning")) {
SVDTruncation = Config::Get().ConfI("smear.SVD.truncation");
QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation)
}
QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- "
<< TrueBinH << "], Rec: " << RecNBins
<< ", [" << RecBinL << " -- " << RecBinH
<< "]");
ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins,
TrueBinL, TrueBinH);
ERecDistrib = new TH1D("ELepRec_rate", ";Rec E_{#nu};Count", RecNBins,
RecBinL, RecBinH);
ETrueDistrib->Sumw2();
ERecDistrib->Sumw2();
RecoSmear =
new TH2D("ELepHadVis_Recon", ";True E_{#nu};Recon. E_{#nu}", RecNBins,
RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH);
RecoSmear->Sumw2();
+
+ // Final setup ---------------------------------------------------
+ FinaliseMeasurement();
}
void Smearceptance_Tester::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
FitPar::Config().out->cd();
- eventVariables = new TTree((this->fName + "_VARS").c_str(),
- (this->fName + "_VARS").c_str());
+ eventVariables =
+ new TTree((fName + "_VARS").c_str(), (fName + "_VARS").c_str());
}
LOG(SAM) << "Adding Event Variables" << std::endl;
eventVariables->Branch("Omega_true", &Omega_true, "Omega_true/F");
eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F");
eventVariables->Branch("Mode_true", &Mode_true, "Mode_true/I");
eventVariables->Branch("EISLep_true", &EISLep_true, "EISLep_true/F");
eventVariables->Branch("HMFS_mu_true", &HMFS_mu_true);
eventVariables->Branch("HMFS_pip_true", &HMFS_pip_true);
eventVariables->Branch("HMFS_pim_true", &HMFS_pim_true);
eventVariables->Branch("HMFS_cpi_true", &HMFS_cpi_true);
eventVariables->Branch("HMFS_p_true", &HMFS_p_true);
eventVariables->Branch("KEFSHad_cpip_true", &KEFSHad_cpip_true,
"KEFSHad_cpip_true/F");
eventVariables->Branch("KEFSHad_cpim_true", &KEFSHad_cpim_true,
"KEFSHad_cpim_true/F");
eventVariables->Branch("KEFSHad_cpi_true", &KEFSHad_cpi_true,
"KEFSHad_cpi_true/F");
eventVariables->Branch("TEFSHad_pi0_true", &TEFSHad_pi0_true,
"TEFSHad_pi0_true/F");
eventVariables->Branch("KEFSHad_p_true", &KEFSHad_p_true, "KEFSHad_p_true/F");
eventVariables->Branch("KEFSHad_n_true", &KEFSHad_n_true, "KEFSHad_n_true/F");
eventVariables->Branch("EFSHad_true", &EFSHad_true, "EFSHad_true/F");
eventVariables->Branch("EFSChargedEMHad_true", &EFSChargedEMHad_true,
"EFSChargedEMHad_true/F");
eventVariables->Branch("EFSLep_true", &EFSLep_true, "EFSLep_true/F");
eventVariables->Branch("EFSgamma_true", &EFSgamma_true, "EFSgamma_true/F");
eventVariables->Branch("PDGISLep_true", &PDGISLep_true, "PDGISLep_true/I");
eventVariables->Branch("PDGFSLep_true", &PDGFSLep_true, "PDGFSLep_true/I");
eventVariables->Branch("Nprotons_true", &Nprotons_true, "Nprotons_true/I");
eventVariables->Branch("Nneutrons_true", &Nneutrons_true, "Nneutrons_true/I");
eventVariables->Branch("Ncpiplus_true", &Ncpiplus_true, "Ncpiplus_true/I");
eventVariables->Branch("Ncpiminus_true", &Ncpiminus_true, "Ncpiminus_true/I");
eventVariables->Branch("Ncpi_true", &Ncpi_true, "Ncpi_true/I");
eventVariables->Branch("Npi0_true", &Npi0_true, "Npi0_true/I");
eventVariables->Branch("HMFS_mu_rec", &HMFS_mu_rec);
eventVariables->Branch("HMFS_pip_rec", &HMFS_pip_rec);
eventVariables->Branch("HMFS_pim_rec", &HMFS_pim_rec);
eventVariables->Branch("HMFS_cpi_rec", &HMFS_cpi_rec);
eventVariables->Branch("HMFS_p_rec", &HMFS_p_rec);
eventVariables->Branch("KEFSHad_cpip_rec", &KEFSHad_cpip_rec,
"KEFSHad_cpip_rec/F");
eventVariables->Branch("KEFSHad_cpim_rec", &KEFSHad_cpim_rec,
"KEFSHad_cpim_rec/F");
eventVariables->Branch("KEFSHad_cpi_rec", &KEFSHad_cpi_rec,
"KEFSHad_cpi_rec/F");
eventVariables->Branch("TEFSHad_pi0_rec", &TEFSHad_pi0_rec,
"TEFSHad_pi0_rec/F");
eventVariables->Branch("KEFSHad_p_rec", &KEFSHad_p_rec, "KEFSHad_p_rec/F");
eventVariables->Branch("KEFSHad_n_rec", &KEFSHad_n_rec, "KEFSHad_n_rec/F");
eventVariables->Branch("EFSHad_rec", &EFSHad_rec, "EFSHad_rec/F");
eventVariables->Branch("EFSLep_rec", &EFSLep_rec, "EFSLep_rec/F");
eventVariables->Branch("EFSVis_cpip", &EFSVis_cpip, "EFSVis_cpip/F");
eventVariables->Branch("EFSVis_cpim", &EFSVis_cpim, "EFSVis_cpim/F");
eventVariables->Branch("EFSVis_cpi", &EFSVis_cpi, "EFSVis_cpi/F");
eventVariables->Branch("EFSVis_pi0", &EFSVis_pi0, "EFSVis_pi0/F");
eventVariables->Branch("EFSVis_p", &EFSVis_p, "EFSVis_p/F");
eventVariables->Branch("EFSVis_n", &EFSVis_n, "EFSVis_n/F");
eventVariables->Branch("EFSVis_gamma", &EFSVis_gamma, "EFSVis_gamma/F");
eventVariables->Branch("EFSVis_other", &EFSVis_other, "EFSVis_other/F");
eventVariables->Branch("EFSVis", &EFSVis, "EFSVis/F");
eventVariables->Branch("FSCLep_seen", &FSCLep_seen, "FSCLep_seen/I");
eventVariables->Branch("Nprotons_seen", &Nprotons_seen, "Nprotons_seen/I");
eventVariables->Branch("Nneutrons_seen", &Nneutrons_seen, "Nneutrons_seen/I");
eventVariables->Branch("Ncpip_seen", &Ncpip_seen, "Ncpip_seen/I");
eventVariables->Branch("Ncpim_seen", &Ncpim_seen, "Ncpim_seen/I");
eventVariables->Branch("Ncpi_seen", &Ncpi_seen, "Ncpi_seen/I");
eventVariables->Branch("Npi0_seen", &Npi0_seen, "Npi0_seen/I");
eventVariables->Branch("Nothers_seen", &Nothers_seen, "Nothers_seen/I");
eventVariables->Branch("EISLep_QE_rec", &EISLep_QE_rec, "EISLep_QE_rec/F");
eventVariables->Branch("EISLep_LepHad_rec", &EISLep_LepHad_rec,
"EISLep_LepHad_rec/F");
eventVariables->Branch("EISLep_LepHadVis_rec", &EISLep_LepHadVis_rec,
"EISLep_LepHadVis_rec/F");
eventVariables->Branch("Nprotons_contributed", &Nprotons_contributed,
"Nprotons_contributed/I");
eventVariables->Branch("Nneutrons_contributed", &Nneutrons_contributed,
"Nneutrons_contributed/I");
eventVariables->Branch("Ncpip_contributed", &Ncpip_contributed,
"Ncpip_contributed/I");
eventVariables->Branch("Ncpim_contributed", &Ncpim_contributed,
"Ncpim_contributed/I");
eventVariables->Branch("Ncpi_contributed", &Ncpi_contributed,
"Ncpi_contributed/I");
eventVariables->Branch("Npi0_contributed", &Npi0_contributed,
"Npi0_contributed/I");
eventVariables->Branch("Ngamma_contributed", &Ngamma_contributed,
"Ngamma_contributed/I");
eventVariables->Branch("Nothers_contibuted", &Nothers_contibuted,
"Nothers_contibuted/I");
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F");
eventVariables->Branch("EffWeight", &EffWeight, "EffWeight/F");
+ xsecScaling = fScaleFactor;
eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F");
eventVariables->Branch("flagCCINC_true", &flagCCINC_true, "flagCCINC_true/O");
eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true, "flagCC0Pi_true/O");
eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O");
eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O");
+
+ if (fEvtRateScaleFactor != 0xdeadbeef) {
+ eventVariables->Branch("PredEvtRateWeight", &PredEvtRateWeight,
+ "PredEvtRateWeight/F");
+ PredEvtRateWeight = fScaleFactor * fEvtRateScaleFactor;
+ }
}
template
int CountNPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum +=
std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]);
}
return sum;
}
template
int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum +=
(std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it])
? 0
: 1);
}
return sum;
}
template
int CountNPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum += std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
pdgs[pdg_it]);
}
return sum;
}
template
int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum += (std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
pdgs[pdg_it])
? 0
: 1);
}
return sum;
}
TLorentzVector GetHMFSRecParticles(RecoInfo ri, int pdg) {
TLorentzVector mom(0, 0, 0, 0);
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if ((ri.RecObjClass[p_it] == pdg) &&
(mom.Mag() < ri.RecObjMom[p_it].Mag())) {
mom.SetXYZM(ri.RecObjMom[p_it].X(), ri.RecObjMom[p_it].Y(),
ri.RecObjMom[p_it].Z(),
PhysConst::GetMass(ri.RecObjClass[p_it]) * 1.0E3);
}
}
return mom;
}
template
double SumKE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass) - mass;
}
return sum;
}
template
double SumTE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass);
}
return sum;
}
template
double SumVisE_RecoInfo(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
template
double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we know about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
//********************************************************************
void Smearceptance_Tester::FillEventVariables(FitEvent *event) {
//********************************************************************
static int const cpipPDG[] = {211};
static int const cpimPDG[] = {-211};
static int const pi0PDG[] = {111};
static int const ProtonPDG[] = {2212};
static int const NeutronPDG[] = {2112};
static int const GammaPDG[] = {22};
static int const CLeptonPDGs[] = {11, 13, 15};
static int const ExplicitPDGs[] = {211, -211, 111, 2212, 2112, 22,
11, 13, 15, 12, 14, 16};
RecoInfo *ri = smearceptor->Smearcept(event);
HMFS_mu_true = TLorentzVector(0, 0, 0, 0);
HMFS_mu_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsMu = event->GetHMFSMuon();
if (fsMu) {
HMFS_mu_true = fsMu->P4();
HMFS_mu_rec = GetHMFSRecParticles(*ri, 13);
}
HMFS_pip_true = TLorentzVector(0, 0, 0, 0);
HMFS_pip_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPip = event->GetHMFSPiPlus();
if (fsPip) {
HMFS_pip_true = fsPip->P4();
HMFS_pip_rec = GetHMFSRecParticles(*ri, 211);
}
HMFS_pim_true = TLorentzVector(0, 0, 0, 0);
HMFS_pim_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPim = event->GetHMFSPiMinus();
if (fsPim) {
HMFS_pim_true = fsPim->P4();
HMFS_pim_rec = GetHMFSRecParticles(*ri, -211);
}
HMFS_cpi_true = TLorentzVector(0, 0, 0, 0);
HMFS_cpi_rec = TLorentzVector(0, 0, 0, 0);
if (fsPip || fsPim) {
if (!fsPip) {
HMFS_cpi_true = HMFS_pim_true;
HMFS_cpi_rec = HMFS_pim_rec;
} else if (!fsPim) {
HMFS_cpi_true = HMFS_pip_true;
HMFS_cpi_rec = HMFS_pip_rec;
} else {
HMFS_cpi_true =
(fsPip->p2() > fsPim->p2()) ? HMFS_pip_true : HMFS_pim_true;
HMFS_cpi_rec = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_rec : HMFS_pim_rec;
}
}
HMFS_p_true = TLorentzVector(0, 0, 0, 0);
HMFS_p_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsP = event->GetHMFSProton();
if (fsP) {
HMFS_p_true = fsP->P4();
HMFS_p_rec = GetHMFSRecParticles(*ri, 2212);
}
TLorentzVector FourMomentumTransfer =
(event->GetHMISAnyLeptons()->P4() - event->GetHMFSAnyLeptons()->P4());
Omega_true = FourMomentumTransfer.E();
Q2_true = -1 * FourMomentumTransfer.Mag2();
Mode_true = event->Mode;
EISLep_true = event->GetHMISAnyLeptons()->E();
KEFSHad_cpip_true = FitUtils::SumTE_PartVect(event->GetAllFSPiPlus());
KEFSHad_cpim_true = FitUtils::SumTE_PartVect(event->GetAllFSPiMinus());
KEFSHad_cpi_true = KEFSHad_cpip_true + KEFSHad_cpim_true;
TEFSHad_pi0_true = FitUtils::SumTE_PartVect(event->GetAllFSPiZero());
KEFSHad_p_true = FitUtils::SumKE_PartVect(event->GetAllFSProton());
KEFSHad_n_true = FitUtils::SumKE_PartVect(event->GetAllFSNeutron());
EFSHad_true =
KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_n_true;
EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true;
EFSLep_true = event->GetHMFSAnyLeptons()->E();
EFSgamma_true = FitUtils::SumTE_PartVect(event->GetAllFSPhoton());
PDGISLep_true = event->GetHMISAnyLeptons()->PDG();
PDGFSLep_true = event->GetHMFSAnyLeptons()->PDG();
Nprotons_true = event->GetAllFSProton().size();
Nneutrons_true = event->GetAllFSNeutron().size();
Ncpiplus_true = event->GetAllFSPiPlus().size();
Ncpiminus_true = event->GetAllFSPiMinus().size();
Ncpi_true = Ncpiplus_true + Ncpiminus_true;
Npi0_true = event->GetAllFSPiZero().size();
KEFSHad_cpip_rec =
SumKE_RecoInfo(*ri, cpipPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpim_rec =
SumKE_RecoInfo(*ri, cpimPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpi_rec = KEFSHad_cpip_rec + KEFSHad_cpim_rec;
TEFSHad_pi0_rec =
SumTE_RecoInfo(*ri, pi0PDG, PhysConst::mass_pi0 * PhysConst::mass_MeV);
KEFSHad_p_rec = SumKE_RecoInfo(*ri, ProtonPDG,
PhysConst::mass_proton * PhysConst::mass_MeV);
KEFSHad_n_rec = SumKE_RecoInfo(*ri, NeutronPDG,
PhysConst::mass_neutron * PhysConst::mass_MeV);
EFSHad_rec =
KEFSHad_cpi_rec + TEFSHad_pi0_rec + KEFSHad_p_rec + KEFSHad_n_rec;
TLorentzVector FSLepMom_rec(0, 0, 0, 0);
if (event->GetHMFSAnyLeptons()) {
FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG());
EFSLep_rec = FSLepMom_rec.E();
} else {
EFSLep_rec = 0;
}
EFSVis_cpip = SumVisE_RecoInfo(*ri, cpipPDG);
EFSVis_cpim = SumVisE_RecoInfo(*ri, cpimPDG);
EFSVis_cpi = EFSVis_cpip + EFSVis_cpim;
EFSVis_pi0 = SumVisE_RecoInfo(*ri, pi0PDG);
EFSVis_p = SumVisE_RecoInfo(*ri, ProtonPDG);
EFSVis_n = SumVisE_RecoInfo(*ri, NeutronPDG);
EFSVis_gamma = SumVisE_RecoInfo(*ri, GammaPDG);
EFSVis_other = SumVisE_RecoInfo_NotPdgs(*ri, ExplicitPDGs);
EFSVis = EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma;
FSCLep_seen = CountNPdgsSeen(*ri, CLeptonPDGs);
Nprotons_seen = CountNPdgsSeen(*ri, ProtonPDG);
Nneutrons_seen = CountNPdgsSeen(*ri, NeutronPDG);
Ncpip_seen = CountNPdgsSeen(*ri, cpipPDG);
Ncpim_seen = CountNPdgsSeen(*ri, cpimPDG);
Ncpi_seen = Ncpip_seen + Ncpim_seen;
Npi0_seen = CountNPdgsSeen(*ri, pi0PDG);
Nothers_seen = CountNNotPdgsSeen(*ri, ExplicitPDGs);
if (FSCLep_seen && (FSLepMom_rec.Mag() > 1E-8)) {
EISLep_QE_rec =
FitUtils::EnuQErec(FSLepMom_rec.Mag() / 1000.0, FSLepMom_rec.CosTheta(),
34, PDGFSLep_true > 0) *
1000.0;
} else {
EISLep_QE_rec = 0;
}
EISLep_LepHad_rec = EFSLep_rec + EFSHad_rec;
EISLep_LepHadVis_rec = EFSLep_rec + EFSHad_rec + EFSVis;
Nprotons_contributed = CountNPdgsContributed(*ri, ProtonPDG);
Nneutrons_contributed = CountNPdgsContributed(*ri, NeutronPDG);
Ncpip_contributed = CountNPdgsContributed(*ri, cpipPDG);
Ncpim_contributed = CountNPdgsContributed(*ri, cpimPDG);
Ncpi_contributed = Ncpip_contributed + Ncpim_contributed;
Npi0_contributed = CountNPdgsContributed(*ri, pi0PDG);
Ngamma_contributed = CountNPdgsContributed(*ri, GammaPDG);
Nothers_contibuted = CountNNotPdgsContributed(*ri, ExplicitPDGs);
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
FluxWeight = GetFluxHistogram()->GetBinContent(
GetFluxHistogram()->FindBin(EISLep_true)) /
GetFluxHistogram()->Integral();
EffWeight = ri->Weight;
- xsecScaling = fScaleFactor;
-
flagCCINC_true = PDGFSLep_true & 1;
flagCC0Pi_true = (Ncpi_true + Npi0_true) == 0;
flagCCINC_rec = FSCLep_seen && PDGFSLep_true & 1;
flagCC0Pi_rec = ((Ncpi_seen + Npi0_seen) == 0) && flagCCINC_rec;
// Fill the eventVariables Tree
eventVariables->Fill();
RecoSmear->Fill(EISLep_true / 1000.0,
flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, Weight);
ETrueDistrib->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight : 0);
ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0, flagCCINC_rec ? Weight : 0);
return;
};
//********************************************************************
void Smearceptance_Tester::Write(std::string drawOpt) {
//********************************************************************
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
TH2D *SmearMatrix_ev =
static_cast(RecoSmear->Clone("ELepHadVis_Smear_ev"));
for (Int_t trueAxis_it = 1;
trueAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++trueAxis_it) {
double NEISLep = ETrueDistrib->GetBinContent(trueAxis_it);
for (Int_t recoAxis_it = 1;
recoAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++recoAxis_it) {
if (NEISLep > std::numeric_limits::epsilon()) {
SmearMatrix_ev->SetBinContent(
trueAxis_it, recoAxis_it,
SmearMatrix_ev->GetBinContent(trueAxis_it, recoAxis_it) / NEISLep);
}
}
}
ETrueDistrib->Write();
ERecDistrib->Write();
RecoSmear->Write();
SmearMatrix_ev->Write();
TH2D *ResponseMatrix_ev =
SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation);
ResponseMatrix_ev->SetName("ResponseMatrix_ev");
ResponseMatrix_ev->Write();
#ifdef DEBUG_SMEARTESTER
TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(SmearMatrix_ev);
TH1D *SmearedEvt = static_cast(ERecDistrib->Clone());
SmearedEvt->SetNameTitle("SmearedEvt", ";Rec E_{#nu}; count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ETrueDistrib, SmearedEvt, SmearMatrix_ev_md, 5000, false);
SmearedEvt->Write();
SmearedEvt->Scale(1, "width");
SmearedEvt->SetName("SmearedEvt_bw");
SmearedEvt->Write();
#endif
TMatrixD ResponseMatrix_evt_md =
SmearceptanceUtils::GetMatrix(ResponseMatrix_ev);
TH1D *Unfolded_enu_obs = static_cast(ETrueDistrib->Clone());
Unfolded_enu_obs->SetNameTitle("UnfoldedENu_evt", ";True E_{#nu};count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ERecDistrib, Unfolded_enu_obs, ResponseMatrix_evt_md, 5000, false);
Unfolded_enu_obs->Write();
Unfolded_enu_obs->Scale(1, "width");
Unfolded_enu_obs->SetName("UnfoldedENu_evt_bw");
Unfolded_enu_obs->Write();
ETrueDistrib->Scale(1, "width");
ETrueDistrib->SetName("ELep_rate_bw");
ETrueDistrib->Write();
ERecDistrib->Scale(1, "width");
ERecDistrib->SetName("ELepRec_rate_bw");
ERecDistrib->Write();
}
// -------------------------------------------------------------------
// Purely MC Plot
// Following functions are just overrides to handle this
// -------------------------------------------------------------------
//********************************************************************
/// Everything is classed as signal...
bool Smearceptance_Tester::isSignal(FitEvent *event) {
//********************************************************************
(void)event;
return true;
};
//********************************************************************
void Smearceptance_Tester::ScaleEvents() {
//********************************************************************
// Saving everything to a TTree so no scaling required
return;
}
//********************************************************************
void Smearceptance_Tester::ApplyNormScale(float norm) {
//********************************************************************
// Saving everything to a TTree so no scaling required
- this->fCurrentNorm = norm;
+ fCurrentNorm = norm;
return;
}
//********************************************************************
void Smearceptance_Tester::FillHistograms() {
//********************************************************************
// No Histograms need filling........
return;
}
//********************************************************************
void Smearceptance_Tester::ResetAll() {
//********************************************************************
eventVariables->Reset();
return;
}
//********************************************************************
float Smearceptance_Tester::GetChi2() {
//********************************************************************
// No Likelihood to test, purely MC
return 0.0;
}
diff --git a/src/MCStudies/Smearceptance_Tester.h b/src/MCStudies/Smearceptance_Tester.h
index f0abbaa..3025580 100644
--- a/src/MCStudies/Smearceptance_Tester.h
+++ b/src/MCStudies/Smearceptance_Tester.h
@@ -1,176 +1,176 @@
// 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 .
*******************************************************************************/
#ifndef Smearceptance_Tester_H_SEEN
#define Smearceptance_Tester_H_SEEN
#include "Measurement1D.h"
#include "ISmearcepter.h"
#ifdef __PROB3PP_ENABLED__
#include "OscWeightEngine.h"
#endif
//********************************************************************
class Smearceptance_Tester : public Measurement1D {
//********************************************************************
public:
- Smearceptance_Tester(std::string name, std::string inputfile, FitWeight *rw,
- std::string type, std::string fakeDataFile);
+ Smearceptance_Tester(nuiskey samplekey);
virtual ~Smearceptance_Tester(){};
//! Grab info from event
void FillEventVariables(FitEvent *event);
//! Fill Custom Histograms
void FillHistograms();
//! ResetAll
void ResetAll();
//! Scale
void ScaleEvents();
//! Norm
void ApplyNormScale(float norm);
//! Define this samples signal
bool isSignal(FitEvent *nvect);
//! Write Files
void Write(std::string drawOpt);
//! Get Chi2
float GetChi2();
void AddEventVariablesToTree();
private:
ISmearcepter *smearceptor;
TTree *eventVariables;
float Omega_true;
float Q2_true;
int Mode_true;
float EISLep_true;
TLorentzVector HMFS_mu_true;
TLorentzVector HMFS_pip_true;
TLorentzVector HMFS_pim_true;
TLorentzVector HMFS_cpi_true;
TLorentzVector HMFS_p_true;
float KEFSHad_cpip_true;
float KEFSHad_cpim_true;
float KEFSHad_cpi_true;
float TEFSHad_pi0_true;
float KEFSHad_p_true;
float KEFSHad_n_true;
float EFSHad_true;
float EFSChargedEMHad_true;
float EFSLep_true;
float EFSgamma_true;
int PDGISLep_true;
int PDGFSLep_true;
int Nprotons_true;
int Nneutrons_true;
int Ncpiplus_true;
int Ncpiminus_true;
int Ncpi_true;
int Npi0_true;
TLorentzVector HMFS_mu_rec;
TLorentzVector HMFS_pip_rec;
TLorentzVector HMFS_pim_rec;
TLorentzVector HMFS_cpi_rec;
TLorentzVector HMFS_p_rec;
float KEFSHad_cpip_rec;
float KEFSHad_cpim_rec;
float KEFSHad_cpi_rec;
float TEFSHad_pi0_rec;
float KEFSHad_p_rec;
float KEFSHad_n_rec;
float EFSHad_rec;
float EFSLep_rec;
float EFSVis_cpip;
float EFSVis_cpim;
float EFSVis_cpi;
float EFSVis_pi0;
float EFSVis_p;
float EFSVis_n;
float EFSVis_gamma;
float EFSVis_other;
float EFSVis;
int FSCLep_seen;
int Nprotons_seen;
int Nneutrons_seen;
int Ncpip_seen;
int Ncpim_seen;
int Ncpi_seen;
int Npi0_seen;
int Nothers_seen;
float EISLep_QE_rec;
float EISLep_LepHad_rec;
float EISLep_LepHadVis_rec;
int Nprotons_contributed;
int Nneutrons_contributed;
int Ncpip_contributed;
int Ncpim_contributed;
int Ncpi_contributed;
int Npi0_contributed;
int Ngamma_contributed;
int Nothers_contibuted;
float Weight;
float RWWeight;
float InputWeight;
float FluxWeight;
float EffWeight;
+ float PredEvtRateWeight;
float xsecScaling;
bool flagCCINC_true;
bool flagCC0Pi_true;
bool flagCCINC_rec;
bool flagCC0Pi_rec;
int SVDTruncation;
TH2D *RecoSmear;
TH1D *ETrueDistrib;
TH1D *ERecDistrib;
};
#endif
diff --git a/src/Utils/PhysConst.h b/src/Utils/PhysConst.h
index f616cc6..ecd9772 100644
--- a/src/Utils/PhysConst.h
+++ b/src/Utils/PhysConst.h
@@ -1,102 +1,105 @@
// 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 .
*******************************************************************************/
#ifndef PHYSCONST_H_SEEN
#define PHYSCONST_H_SEEN
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "FitLogger.h"
/// namespace to contain all physical constants used by NUISANCE
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_nucleon = (mass_proton + mass_neutron)/2.;
+const double mass_proton_kg = 1.6727E-27; // Proton mass in kg
+const double mass_neutron_kg = 1.6750E-27; // Neutron mass in kg
+const double mass_nucleon_kg = (mass_proton_kg + mass_neutron_kg)/2.;
const double mass_delta = 1.232; // Delta mass in GeV
const double mass_muon = 0.10565837; // Muon mass in GeV
const double mass_electron = 0.000510998928; // Electron mass in GeV
const double mass_cpi = 0.13957; // charged pion mass in GeV
const double mass_pi0 = 0.13498; // neutral pion mass in GeV
inline double GetMass(int pdg) {
switch (pdg) {
case 11:
return mass_electron;
case 13:
return mass_muon;
case 111:
return mass_pi0;
case 211:
case -211:
return mass_cpi;
case 2112:
return mass_neutron;
case 2212:
return mass_proton;
default: {
ERROR(WRN, "Attempted to get mass for PDG: "
<< pdg << ", but it is not catered for. Please add it to "
"src/Utils/PhysConst.h");
return -1;
}
}
}
const double mass_MeV = 1000; // MeV/GeV
const int pdg_neutrinos[] = {12, -12, 14, -14 /*, 16, -16*/};
const int pdg_muons[] = {13, -13};
const int pdg_leptons[] = {11, -11, 13, -13, 15, -15};
const int pdg_all_leptons[] = {11, -11, 13, -13, 15, -15,
12, -12, 14, -14, 16, -16};
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