diff --git a/app/PrepareNuwroEvents.cxx b/app/PrepareNuwroEvents.cxx index c6daed4..2b44f8f 100644 --- a/app/PrepareNuwroEvents.cxx +++ b/app/PrepareNuwroEvents.cxx @@ -1,492 +1,492 @@ #include "event1.h" #include #include // Hopefully we don't need these as they're included above. // #include "params_all.h" // #include "params.h" #include "FitLogger.h" #include "PlotUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" void printInputCommands(char *argv[]) { std::cout << "[USAGE]: " << argv[0] << " [-h] [-f] [-F ,] [-o output.root] " "inputfile.root [file2.root ...]" << std::endl << "\t-h : Print this message." << std::endl << "\t-f : Pass -f argument to '$ hadd' invocation." << std::endl << "\t-F : Read input flux from input descriptor." << std::endl << "\t-o : Write full output to a new file." << std::endl << std::endl; }; void CreateRateHistograms(std::string inputs, bool force_out); void HaddNuwroFiles(std::vector &inputs, bool force_out); bool outputNewFile = false; std::string ofile = ""; bool haveFluxInputs = false; struct FluxInputBlob { FluxInputBlob(std::string _File, std::string _Hist, int _PDG, double _Fraction) : File(_File), Hist(_Hist), PDG(_PDG), Fraction(_Fraction) {} std::string File; std::string Hist; int PDG; double Fraction; }; std::vector FluxInputs; bool haddedFiles = false; TH1D *F2D(TH1F *f) { Double_t *bins = new Double_t[f->GetXaxis()->GetNbins() + 1]; for (Int_t bi_it = 0; bi_it < f->GetXaxis()->GetNbins(); ++bi_it) { bins[bi_it] = f->GetXaxis()->GetBinLowEdge(bi_it + 1); } bins[f->GetXaxis()->GetNbins()] = f->GetXaxis()->GetBinUpEdge(f->GetXaxis()->GetNbins()); TH1D *d = new TH1D((std::string(f->GetName()) + "_f").c_str(), f->GetTitle(), f->GetXaxis()->GetNbins(), bins); std::cout << "Converted TH1F with " << f->GetXaxis()->GetXbins() << " bins : " << std::endl; for (Int_t bi_it = 0; bi_it < f->GetXaxis()->GetNbins() + 2; ++bi_it) { d->SetBinContent(bi_it, f->GetBinContent(bi_it)); d->SetBinError(bi_it, f->GetBinError(bi_it)); std::cout << "\tF " << f->GetXaxis()->GetBinLowEdge(bi_it) << "--[" << bi_it << "]--" << f->GetXaxis()->GetBinUpEdge(bi_it) << ": " << f->GetBinContent(bi_it) << " D " << d->GetXaxis()->GetBinLowEdge(bi_it) << "--[" << bi_it << "]--" << d->GetXaxis()->GetBinUpEdge(bi_it) << ": " << d->GetBinContent(bi_it) << std::endl; } delete bins; return d; } TH1D *GetTH1DFromFile(std::string const &rootFile, std::string const &histName) { TFile *inpFile = new TFile(rootFile.c_str(), "READ"); if (!inpFile || !inpFile->IsOpen()) { NUIS_ABORT("Cannot open input root file: " << rootFile - << " to read input histo."); + << " to read input histo."); } TH1D *histD = dynamic_cast(inpFile->Get(histName.c_str())); if (!histD) { TH1F *histF = dynamic_cast(inpFile->Get(histName.c_str())); if (!histF) { - NUIS_ABORT("Cannot find TH1D/F: " << histName << " in root file: " << rootFile - << "."); + NUIS_ABORT("Cannot find TH1D/F: " << histName << " in root file: " + << rootFile << "."); } histD = F2D(histF); } else { histD = static_cast(histD->Clone()); } histD->SetDirectory(NULL); inpFile->Close(); return histD; } //******************************* int main(int argc, char *argv[]) { //******************************* // If No Arguments print commands if (argc == 1) { printInputCommands(argv); return 0; } std::vector inputfiles; bool force_output = false; // Get Inputs for (int i = 1; i < argc; ++i) { if (!std::strcmp(argv[i], "-h")) { printInputCommands(argv); return 0; } else if (!std::strcmp(argv[i], "-f")) { force_output = true; } else if (!std::strcmp(argv[i], "-o")) { outputNewFile = true; ofile = argv[++i]; } else if (!std::strcmp(argv[i], "-F")) { std::string inpLine = argv[++i]; std::vector fluxInputDescriptor = GeneralUtils::ParseToStr(inpLine, ","); if ((fluxInputDescriptor.size() != 2) && (fluxInputDescriptor.size() != 3) && (fluxInputDescriptor.size() != 4)) { NUIS_ABORT("Received -F argument with option: \"" - << inpLine - << "\", was expecting " - ",[,PDG[,speciesFraction]]."); + << inpLine + << "\", was expecting " + ",[,PDG[,speciesFraction]]."); } haveFluxInputs = true; FluxInputs.push_back( FluxInputBlob(fluxInputDescriptor[0], fluxInputDescriptor[1], (fluxInputDescriptor.size() > 2) ? GeneralUtils::StrToInt(fluxInputDescriptor[2]) : 14, (fluxInputDescriptor.size() > 3) ? GeneralUtils::StrToDbl(fluxInputDescriptor[3]) : 1)); if (!FluxInputs.back().File.length() || !FluxInputs.back().Hist.length()) { NUIS_ABORT("Received -F argument with option: \"" - << inpLine - << "\", was expecting " - ",[,PDG[,speciesFraction]]."); + << inpLine + << "\", was expecting " + ",[,PDG[,speciesFraction]]."); } } else { inputfiles.push_back(std::string(argv[i])); } } // If one input file just create flux histograms if (inputfiles.size() > (UInt_t)1) { HaddNuwroFiles(inputfiles, force_output); } else if (inputfiles.size() < (UInt_t)1) { printInputCommands(argv); } CreateRateHistograms(inputfiles[0], force_output); NUIS_LOG(FIT, "Finished NUWRO Prep."); }; //******************************* void CreateRateHistograms(std::string inputs, bool force_out) { //******************************* // Open root file TFile *outRootFile = 0; TTree *nuwrotree = 0; if (!haddedFiles && outputNewFile) { // we need to make the new file and clone the tree. TFile *inpFile = new TFile(inputs.c_str(), "READ"); if (!inpFile || !inpFile->IsOpen()) { NUIS_ABORT("Cannot open input root file: " << inputs); } TTree *inpTree = dynamic_cast(inpFile->Get("treeout")); if (!inpTree) { NUIS_ABORT("Cannot find TTree \"treeout\" in input root file: " - << inputs.c_str()); + << inputs.c_str()); } outRootFile = new TFile(ofile.c_str(), force_out ? "RECREATE" : "CREATE"); if (!outRootFile || !outRootFile->IsOpen()) { NUIS_ABORT("Couldn't open root file: " - << ofile << " for writing, does it already exist?"); + << ofile << " for writing, does it already exist?"); } nuwrotree = inpTree->CloneTree(-1, "fast"); nuwrotree->SetDirectory(outRootFile); nuwrotree->Write(nuwrotree->GetName()); } else { outRootFile = new TFile(inputs.c_str(), "UPDATE"); if (!outRootFile || !outRootFile->IsOpen()) { NUIS_ABORT("Cannot open input root file: " << inputs); } nuwrotree = dynamic_cast(outRootFile->Get("treeout")); if (!nuwrotree) { NUIS_ABORT("Cannot find TTree \"treeout\" in input root file: " - << inputs.c_str()); + << inputs.c_str()); } } // Get Flux Histogram event *evt = new event(); nuwrotree->SetBranchAddress("e", &evt); nuwrotree->GetEntry(0); int fluxtype = evt->par.beam_type; std::map fluxlist; std::map eventlist; std::vector allpdg; std::map nevtlist; std::map intxseclist; // Did the input file have a mono-energetic flux? bool isMono = false; nevtlist[0] = 0.0; intxseclist[0] = 0.0; allpdg.push_back(0); NUIS_LOG(FIT, "Nuwro fluxtype = " << fluxtype); if (haveFluxInputs) { double totalFraction = 0; for (size_t flux_it = 0; flux_it < FluxInputs.size(); ++flux_it) { FluxInputBlob &fb = FluxInputs[flux_it]; int pdg = fb.PDG; TH1D *fluxHist = GetTH1DFromFile(fb.File, fb.Hist); double pctg = fb.Fraction; totalFraction += pctg; double Elow = fluxHist->GetXaxis()->GetBinLowEdge(1); double Ehigh = fluxHist->GetXaxis()->GetBinLowEdge( fluxHist->GetXaxis()->GetNbins() + 1); NUIS_LOG(FIT, "Adding new nuwro flux " - << "pdg: " << pdg << " pctg: " << pctg << " Elow: " << Elow - << " Ehigh: " << Ehigh); + << "pdg: " << pdg << " pctg: " << pctg + << " Elow: " << Elow << " Ehigh: " << Ehigh); // Sort total flux plot if (!fluxlist[0]) { // Setup total flux fluxlist[0] = (TH1D *)fluxHist->Clone(); fluxlist[0]->SetNameTitle("FluxHist", "FluxHist"); // Prep empty total events eventlist[0] = (TH1D *)fluxHist->Clone(); eventlist[0]->SetNameTitle("EvtHist", "EvtHist"); eventlist[0]->Reset(); } else { // Add up each new plot fluxlist[0]->Add(fluxHist); } fluxHist->SetNameTitle(Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg), Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg)); TH1D *eventplot = (TH1D *)fluxHist->Clone(); eventplot->SetNameTitle(Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg), Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg)); eventplot->Reset(); fluxlist[pdg] = (TH1D *)fluxHist->Clone(); eventlist[pdg] = eventplot; nevtlist[pdg] = 0; intxseclist[pdg] = 0.0; allpdg.push_back(pdg); delete fluxHist; } if (fabs(totalFraction - 1) > 1E-5) { - NUIS_ABORT(FTL, "Total species fraction for input flux histos = " - << totalFraction << ", expected to sum to 1."); + NUIS_ABORT("Total species fraction for input flux histos = " + << totalFraction << ", expected to sum to 1."); } } else if (fluxtype == 0) { std::string fluxstring = evt->par.beam_energy; std::vector fluxvals = GeneralUtils::ParseToDbl(fluxstring, " "); int pdg = evt->par.beam_particle; double Elow = double(fluxvals[0]) / 1000.0; double Ehigh = double(fluxvals[1]) / 1000.0; TH1D *fluxplot = NULL; if (Elow > Ehigh) isMono = true; // For files produced with a flux distribution if (!isMono) { NUIS_LOG(FIT, "Adding new nuwro flux " - << "pdg: " << pdg << " Elow: " << Elow - << " Ehigh: " << Ehigh); + << "pdg: " << pdg << " Elow: " << Elow + << " Ehigh: " << Ehigh); fluxplot = new TH1D("fluxplot", "fluxplot", fluxvals.size() - 4, Elow, Ehigh); for (uint j = 2; j < fluxvals.size(); j++) { NUIS_LOG(DEB, j << " " << fluxvals[j]); fluxplot->SetBinContent(j - 1, fluxvals[j]); } } else { // For monoenergetic fluxes NUIS_LOG(FIT, "Adding mono-energetic nuwro flux " - << "pdg: " << pdg << " E: " << Elow); + << "pdg: " << pdg << " E: " << Elow); fluxplot = new TH1D("fluxplot", "fluxplot", 100, 0, Elow * 2); fluxplot->SetBinContent(fluxplot->FindBin(Elow), 1); } // Setup total flux fluxlist[0] = (TH1D *)fluxplot->Clone(); fluxlist[0]->SetNameTitle("FluxHist", "FluxHist"); // Prep empty total events eventlist[0] = (TH1D *)fluxplot->Clone(); eventlist[0]->SetNameTitle("EvtHist", "EvtHist"); eventlist[0]->Reset(); fluxplot->SetNameTitle(Form("nuwro_pdg%i_Flux", pdg), Form("nuwro_pdg%i_Flux", pdg)); TH1D *eventplot = (TH1D *)fluxplot->Clone(); eventplot->SetNameTitle(Form("nuwro_pdg%i_Evt", pdg), Form("nuwro_pdg%i_Evt", pdg)); eventplot->Reset(); fluxlist[pdg] = fluxplot; eventlist[pdg] = eventplot; nevtlist[pdg] = 0; intxseclist[pdg] = 0.0; allpdg.push_back(pdg); } else if (fluxtype == 1) { std::string fluxstring = evt->par.beam_content; std::vector fluxlines = GeneralUtils::ParseToStr(fluxstring, "\n"); for (uint i = 0; i < fluxlines.size(); i++) { std::vector fluxvals = GeneralUtils::ParseToDbl(fluxlines[i], " "); int pdg = int(fluxvals[0]); double pctg = double(fluxvals[1]) / 100.0; double Elow = double(fluxvals[2]) / 1000.0; double Ehigh = double(fluxvals[3]) / 1000.0; NUIS_LOG(FIT, "Adding new nuwro flux " - << "pdg: " << pdg << " pctg: " << pctg << " Elow: " << Elow - << " Ehigh: " << Ehigh); + << "pdg: " << pdg << " pctg: " << pctg + << " Elow: " << Elow << " Ehigh: " << Ehigh); TH1D *fluxplot = new TH1D("fluxplot", "fluxplot", fluxvals.size() - 4, Elow, Ehigh); for (uint j = 4; j < fluxvals.size(); j++) { fluxplot->SetBinContent(j + 1, fluxvals[j]); } // Sort total flux plot if (!fluxlist[0]) { // Setup total flux fluxlist[0] = (TH1D *)fluxplot->Clone(); fluxlist[0]->SetNameTitle("FluxHist", "FluxHist"); // Prep empty total events eventlist[0] = (TH1D *)fluxplot->Clone(); eventlist[0]->SetNameTitle("EvtHist", "EvtHist"); eventlist[0]->Reset(); } else { // Add up each new plot fluxlist[0]->Add(fluxplot); } fluxplot->SetNameTitle(Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg), Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg)); TH1D *eventplot = (TH1D *)fluxplot->Clone(); eventplot->SetNameTitle(Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg), Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg)); eventplot->Reset(); fluxlist[pdg] = fluxplot; eventlist[pdg] = eventplot; nevtlist[pdg] = 0; intxseclist[pdg] = 0.0; allpdg.push_back(pdg); } } // Start main event loop to fill plots int nevents = nuwrotree->GetEntries(); double Enu = 0.0; double TotXSec = 0.0; // double totaleventmode = 0.0; // double totalevents = 0.0; int pdg = 0; int countwidth = nevents / 50.0; countwidth = countwidth ? countwidth : 1; for (int i = 0; i < nevents; i++) { nuwrotree->GetEntry(i); // Get Variables Enu = evt->in[0].t / 1000.0; TotXSec = evt->weight; pdg = evt->in[0].pdg; eventlist[0]->Fill(Enu); eventlist[pdg]->Fill(Enu); nevtlist[0] += 1; nevtlist[pdg] += 1; intxseclist[0] += TotXSec; intxseclist[pdg] += TotXSec; if (i % countwidth == 0) { NUIS_LOG(FIT, "Processed " << i << " events " - << " (" << int(i * 100.0 / nevents) << "%)" - << " : E, W, PDG = " << Enu << ", " << TotXSec - << ", " << pdg) + << " (" << int(i * 100.0 / nevents) << "%)" + << " : E, W, PDG = " << Enu << ", " << TotXSec + << ", " << pdg) } } TH1D *zeroevents = (TH1D *)eventlist[0]->Clone(); outRootFile->cd(); // Loop over eventlist for (uint i = 0; i < allpdg.size(); i++) { int pdg = allpdg[i]; double AvgXSec = intxseclist[0] * 1E38 / double(nevtlist[0]); NUIS_LOG(FIT, pdg << " Avg XSec = " << AvgXSec); NUIS_LOG(FIT, pdg << " nevents = " << double(nevtlist[pdg])); if (!isMono) { // Convert events to PDF eventlist[pdg]->Scale(1.0 / zeroevents->Integral("width")); // Multiply by total predicted event rate eventlist[pdg]->Scale(fluxlist[0]->Integral("width") * AvgXSec); } else { // If a mono-energetic flux was used, width should not be used // The output is (now) forced to be flux = 1, evtrt = xsec (in 1E38 * nb // cm^2) eventlist[pdg]->Scale(1.0 / zeroevents->Integral()); eventlist[pdg]->Scale(fluxlist[0]->Integral() * AvgXSec); } // Save everything fluxlist[pdg]->Write("", TObject::kOverwrite); eventlist[pdg]->Write("", TObject::kOverwrite); } // Tidy up outRootFile->Close(); fluxlist.clear(); eventlist.clear(); // Exit Program return; } //******************************* void HaddNuwroFiles(std::vector &inputs, bool force_out) { //******************************* // Get output file name std::string outputname = inputs[0]; // Make command line string std::string cmd = "hadd "; if (outputNewFile) { cmd += ofile + " "; outputname = ofile; } else if (force_out) { cmd += "-f "; } for (UInt_t i = 0; i < inputs.size(); i++) { cmd += inputs[i] + " "; } NUIS_LOG(FIT, " Running HADD from PrepareNuwro: " << cmd); // Start HADD system(cmd.c_str()); // Return name of output file inputs.clear(); inputs.push_back(outputname); haddedFiles = true; return; } diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake index 43d8096..2b67e48 100644 --- a/cmake/GENIESetup.cmake +++ b/cmake/GENIESetup.cmake @@ -1,252 +1,251 @@ # 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 . ################################################################################ # TODO # check system for libxml2 # check whether we need the includes # check if we can use a subset of the GENIE libraries include(${CMAKE_SOURCE_DIR}/cmake/parseConfigApp.cmake) ################################################################################ # Check Dependencies ################################################################################ ################################# GENIE ###################################### if(GENIE STREQUAL "") cmessage(FATAL_ERROR "Variable GENIE is not defined. " "The location of a pre-built GENIE install must be defined either as" " $ cmake -DGENIE=/path/to/GENIE or as an environment variable" " $ export GENIE=/path/to/GENIE") endif() if (BUILD_GEVGEN) cmessage(STATUS "Building custom gevgen") LIST(APPEND EXTRA_CXX_FLAGS -D__GEVGEN_ENABLED__) endif() execute_process(COMMAND genie-config --version OUTPUT_VARIABLE GENIE_VER OUTPUT_STRIP_TRAILING_WHITESPACE) cmessage(STATUS "genie_ver: ${GENIE_VER}") if(GENIE_VER VERSION_GREATER 3.0.0) set(GENIE_POST_R3 1) string(REPLACE "." "" GENIE_VERSION ${GENIE_VER}) cmessage(STATUS "set genie_post_r3") endif() if(NOT GENIE_POST_R3) LIST(APPEND EXTRA_CXX_FLAGS -DGENIE_PRE_R3) cmessage(STATUS "setting genie_pre_r3 ${EXTRA_CXX_FLAGS}") if(GENIE_EMPMEC_REWEIGHT) cmessage(STATUS "Enable EMPMEC dials") LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED) endif() else() cmessage(STATUS "Enable EMPMEC dials") LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED) endif() execute_process (COMMAND genie-config --topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE) #Allows for external override in the case where genie-config lies. if(NOT DEFINED GENIE_LIB_DIR OR GENIE_LIB_DIR STREQUAL "") #This looks like it should call libdir, but it strips the argument with -L from the response of --libs GETLIBDIRS(genie-config --libs GENIE_LIB_DIR) endif() GETLIBS(genie-config --libs GENIE_LIBS) cmessage(STATUS "GENIE version : ${GENIE_VERSION}") cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}") cmessage(STATUS "GENIE libs : ${GENIE_LIBS}") string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS}) if(WASMATCHED AND GENIE_VERSION STREQUAL "210") set(GENIE_SEHGAL ${GENIE_LIBS}) STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS ${GENIE_SEHGAL}) cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS}") endif() if(NOT USE_REWEIGHT) SET(USING_GENIE_RW FALSE) elseif(NOT GENIE_POST_R3) LIST(FIND GENIE_LIBS GReWeight FOUND_GENIE_RW) if(FOUND_GENIE_RW EQUAL -1) cmessage(DEBUG "Did NOT find ReWeight library. Here are libs: ${GENIE_LIBS}") SET(USING_GENIE_RW FALSE) else() SET(USING_GENIE_RW TRUE) endif() elseif(DEFINED GENIE_REWEIGHT AND NOT GENIE_REWEIGHT STREQUAL "") LIST(FIND GENIE_LIBS GRwFwk FOUND_GENIE_RW) if(FOUND_GENIE_RW EQUAL -1) LIST(APPEND GENIE_LIBS GRwClc GRwFwk GRwIO) cmessage(DEBUG "Force added ReWeight library. Here are libs: ${GENIE_LIBS}") SET(USING_GENIE_RW TRUE) else() SET(USING_GENIE_RW FALSE) endif() endif() if(USING_GENIE_RW) cmessage(STATUS "Using GENIE ReWeight library.") else() cmessage(STATUS "Building without GENIE ReWeight support.") - LIST(APPEND EXTRA_CXX_FLAGS -D__NO_GENIE_REWEIGHT__) endif() LIST(APPEND GENIE_LIBS -Wl,--end-group ) LIST(REVERSE GENIE_LIBS) LIST(APPEND GENIE_LIBS -Wl,--start-group -Wl,--no-as-needed ) LIST(REVERSE GENIE_LIBS) cmessage(DEBUG "GENIE_LIBS: ${GENIE_LIBS}") ################################ LHAPDF ###################################### if(LHAPDF_LIB STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as an environment variable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries") endif() if(LHAPDF_INC STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as an environment variable $ export LHAPDF_INC=/path/to/LHAPDF_includes") endif() if(LHAPATH STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as an environment variable $ export LHAPATH=/path/to/LHAPATH") endif() ################################ LIBXML ###################################### if(LIBXML2_LIB STREQUAL "") GETLIBDIR(xml2-config --libs LIBXML2_LIB IGNORE_EMPTY_RESPONSE) if(LIBXML2_LIB STREQUAL "") message(WARNING "Variable LIBXML2_LIB is not defined, as xml2-config was found and didn't report a library include path, it is likely that libxml2.so can be found in the standard system location, lets hope so.") endif() endif() if(LIBXML2_INC STREQUAL "") GETINCDIR(xml2-config --cflags LIBXML2_INC) if(LIBXML2_INC STREQUAL "") message(FATAL_ERROR "Variable LIBXML2_INC is not defined and could not be found with xml2-config. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_includes or as an environment variable $ export LIBXML2_INC=/path/to/LIBXML2_includes") endif() endif() ############################### log4cpp ###################################### if(LOG4CPP_LIB STREQUAL "") GETLIBDIR(log4cpp-config --libs LOG4CPP_LIB) if(LOG4CPP_LIB STREQUAL "") message(FATAL_ERROR "Variable LOG4CPP_LIB is not defined and could not be found with log4cpp-config. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as an environment variable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries") endif() endif() if(LOG4CPP_INC STREQUAL "") GETINCDIR(log4cpp-config --cflags LOG4CPP_INC) if(LOG4CPP_INC STREQUAL "") message(FATAL_ERROR "Variable LOG4CPP_INC is not defined and could not be found with log4cpp-config. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_INC=/path/to/LOG4CPP_includes or as an environment variable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes") endif() endif() ################################################################################ # Set the compiler defines LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION}) LIST(APPEND EXTRA_LIBS ${GENIE_LIBS}) ############################### GSL ###################################### if(GENIE_POST_R3) if(GSL_LIB STREQUAL "") GETLIBDIR(gsl-config --libs GSL_LIB) if(GSL_LIB STREQUAL "") message(FATAL_ERROR "Variable GSL_LIB is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_LIB=/path/to/GSL_libraries or as an environment variable $ export GSL_LIB=/path/to/GSL_libraries") endif() endif() if(GSL_INC STREQUAL "") GETINCDIR(gsl-config --cflags GSL_INC) if(GSL_INC STREQUAL "") message(FATAL_ERROR "Variable GSL_INC is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_INC=/path/to/GSL_includes or as an environment variable $ export GSL_INC=/path/to/GSL_includes") endif() endif() GETLIBS(gsl-config --libs GSL_LIB_LIST) if(USING_GENIE_RW AND GENIE_REWEIGHT STREQUAL "") message(FATAL_ERROR "Variable GENIE_REWEIGHT is not defined. When using GENIE v3+, we require the reweight product to be built and accessible via the environment variable GENIE_REWEIGHT") endif() endif() ################################################################################ LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp) LIST(APPEND EXTRA_LINK_DIRS ${GENIE_LIB_DIR} ${LHAPDF_LIB} ${LOG4CPP_LIB}) # Append only if we have found GENIE ReWeight if(NOT GENIE_POST_R3) LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_INCLUDES_DIR} ${GENIE_INCLUDES_DIR}/GHEP ${GENIE_INCLUDES_DIR}/Ntuple) if(USING_GENIE_RW) LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_INCLUDES_DIR}/ReWeight) endif() LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_INCLUDES_DIR}/Apps ${GENIE_INCLUDES_DIR}/FluxDrivers ${GENIE_INCLUDES_DIR}/EVGDrivers ${LHAPDF_INC} ${LIBXML2_INC} ${LOG4CPP_INC}) else() LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_INCLUDES_DIR}) if(USING_GENIE_RW) LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_REWEIGHT}/src) endif() LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GSL_INC} ${LHAPDF_INC} ${LIBXML2_INC} ${LOG4CPP_INC}) if(USING_GENIE_RW) LIST(APPEND EXTRA_LINK_DIRS ${GENIE_REWEIGHT}/lib) endif() LIST(APPEND EXTRA_LINK_DIRS ${GSL_LIB} ) LIST(APPEND EXTRA_LIBS ${GSL_LIB_LIST}) endif() if(USE_PYTHIA8) set(NEED_PYTHIA8 TRUE) set(NEED_ROOTPYTHIA8 TRUE) else() set(NEED_PYTHIA6 TRUE) set(NEED_ROOTPYTHIA6 TRUE) endif() set(NEED_ROOTEVEGEN TRUE) SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE) diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake index 3cd85c6..eb4dd92 100644 --- a/cmake/NEUTSetup.cmake +++ b/cmake/NEUTSetup.cmake @@ -1,145 +1,209 @@ # 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 . ################################################################################ -if(NEUT_ROOT STREQUAL "") - cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.") -endif() +include(cmake/parseConfigApp.cmake) -if(CERN STREQUAL "") - cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.") -endif() +find_program(NEUTCONFIG NAMES neut-config) -if(CERN_LEVEL STREQUAL "") - cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).") -endif() +LIST(APPEND EXTRA_CXX_FLAGS -DNEED_FILL_NEUT_COMMONS) -if(${NEUT_VERSION} VERSION_LESS 5.4.0) - set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc) +SET(HAVENEUTCONFIG FALSE) +# We are dealing with shiny NEUT +if(NOT "${NEUTCONFIG}" STREQUAL "NEUTCONFIG-NOTFOUND") + SET(HAVENEUTCONFIG TRUE) + cmessage(STATUS "Found neut-config, using it to determine configuration.") else() - set(NEUT_LIB_DIR ${NEUT_ROOT}/lib) + cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.") endif() -set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass) - -LIST(APPEND EXTRA_CXX_FLAGS -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION}) - -LIST(APPEND EXTRA_CXX_FLAGS - -I${NEUT_ROOT}/include - -I${NEUT_ROOT}/src/neutclass - -I${NEUT_ROOT}/src/reweight) - -LIST(APPEND EXTRA_LINK_DIRS - ${NEUT_LIB_DIR} - ${CERN}/${CERN_LEVEL}/lib - ${NEUT_ROOT}/src/reweight) - -if(${NEUT_VERSION} VERSION_EQUAL 5.4.2) - LIST(APPEND EXTRA_LIBS - -Wl,--as-needed - NReWeight - -Wl,--start-group - neutcore_5.4.2 - nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear - nuceff_5.4.2 - partnuck_5.4.2 - skmcsvc_5.4.2 - tauola_5.4.2 - HT2p2h_5.4.0 - N1p1h_5.4.0 - -Wl,--end-group - jetset74 - pdflib804 - mathlib - packlib - pawlib) - - LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV) -elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0) - LIST(APPEND EXTRA_LIBS - -Wl,--as-needed - NReWeight - -Wl,--start-group - neutcore_5.4.0 - nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear - nuceff_5.4.0 - partnuck_5.4.0 - skmcsvc_5.4.0 - tauola_5.4.0 - HT2p2h_5.4.0 - N1p1h_5.4.0 - specfunc_5.4.0 - radcorr_5.4.0 - gfortran - -Wl,--end-group - jetset74 - pdflib804 - mathlib - packlib - pawlib) -else() - LIST(APPEND EXTRA_LIBS - -Wl,--as-needed - NReWeight - -Wl,--start-group - neutcore - nuccorrspl - nuceff - partnuck - skmcsvc - tauola - -Wl,--end-group - jetset74 - pdflib804 - mathlib - packlib - pawlib) -endif() +if(HAVENEUTCONFIG) + execute_process (COMMAND neut-config + --version OUTPUT_VARIABLE NEUT_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE) + + execute_process (COMMAND neut-config + --incdir OUTPUT_VARIABLE NEUT_INCLUDE_DIRS + OUTPUT_STRIP_TRAILING_WHITESPACE) + + execute_process (COMMAND neut-config + --libdir OUTPUT_VARIABLE NEUT_LINK_DIRS + OUTPUT_STRIP_TRAILING_WHITESPACE) + + GETLIBDIRS(neut-config --cernflags CERN_LIB_DIR) + LIST(APPEND NEUT_LINK_DIRS ${CERN_LIB_DIR}) + GETLIBS(neut-config --cernflags CERN_LIBS) + + if(USE_REWEIGHT) + execute_process (COMMAND neut-config + --rwlibflags OUTPUT_VARIABLE NEUT_RWLIBS + OUTPUT_STRIP_TRAILING_WHITESPACE) + GETLIBS(neut-config --rwlibflags NEUT_RWLIBS) + LIST(APPEND NEUT_LIBS ${NEUT_RWLIBS}) + else() + GETLIBS(neut-config --libflags NEUT_GENLIBS) + GETLIBS(neut-config --iolibflags NEUT_LIBS) + LIST(APPEND NEUT_LIBS ${NEUT_IOLIBS}) + LIST(APPEND NEUT_LIBS ${NEUT_GENLIBS}) + endif() + + LIST(APPEND NEUT_LIBS ${CERN_LIBS};gfortran) + LIST(APPEND EXTRA_LIBS ${NEUT_LIBS}) + + PrefixList(NEUT_INCLUDE_DIRS "-I" ${NEUT_INCLUDE_DIRS}) + LIST(APPEND EXTRA_CXX_FLAGS ${NEUT_INCLUDE_DIRS} -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION}) + + LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS}) + + cmessage(STATUS "NEUT") + cmessage(STATUS " Version : ${NEUT_VERSION}") + cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}") + cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}") + cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}") + cmessage(STATUS " Libs : ${NEUT_LIBS}") + + +else() # Everything better be set up already + + if(NEUT_ROOT STREQUAL "") + cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.") + endif() + + if(CERN STREQUAL "") + cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.") + endif() + + if(CERN_LEVEL STREQUAL "") + cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).") + endif() + + if(${NEUT_VERSION} VERSION_LESS 5.4.0) + set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc) + else() + set(NEUT_LIB_DIR ${NEUT_ROOT}/lib) + endif() + + set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass) + + LIST(APPEND EXTRA_CXX_FLAGS -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION}) + + LIST(APPEND EXTRA_CXX_FLAGS + -I${NEUT_ROOT}/include + -I${NEUT_ROOT}/src/neutclass + -I${NEUT_ROOT}/src/reweight) + + LIST(APPEND EXTRA_LINK_DIRS + ${NEUT_LIB_DIR} + ${CERN}/${CERN_LEVEL}/lib + ${NEUT_ROOT}/src/reweight) + + if(${NEUT_VERSION} VERSION_EQUAL 5.4.2) + LIST(APPEND EXTRA_LIBS + -Wl,--as-needed + NReWeight + -Wl,--start-group + neutcore_5.4.2 + nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear + nuceff_5.4.2 + partnuck_5.4.2 + skmcsvc_5.4.2 + tauola_5.4.2 + HT2p2h_5.4.0 + N1p1h_5.4.0 + -Wl,--end-group + jetset74 + pdflib804 + mathlib + packlib + pawlib) + + LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV) + elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0) + LIST(APPEND EXTRA_LIBS + -Wl,--as-needed + NReWeight + -Wl,--start-group + neutcore_5.4.0 + nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear + nuceff_5.4.0 + partnuck_5.4.0 + skmcsvc_5.4.0 + tauola_5.4.0 + HT2p2h_5.4.0 + N1p1h_5.4.0 + specfunc_5.4.0 + radcorr_5.4.0 + gfortran + -Wl,--end-group + jetset74 + pdflib804 + mathlib + packlib + pawlib) + else() + LIST(APPEND EXTRA_LIBS + -Wl,--as-needed + NReWeight + -Wl,--start-group + neutcore + nuccorrspl + nuceff + partnuck + skmcsvc + tauola + -Wl,--end-group + jetset74 + pdflib804 + mathlib + packlib + pawlib) + endif() + + set(NEUT_ROOT_LIBS) -set(NEUT_ROOT_LIBS) + LIST(APPEND NEUT_ROOT_LIBS + ${NEUT_CLASS}/neutctrl.so + ${NEUT_CLASS}/neutfsivert.so) -LIST(APPEND NEUT_ROOT_LIBS - ${NEUT_CLASS}/neutctrl.so - ${NEUT_CLASS}/neutfsivert.so) + # Check for new versions of NEUT with NUCLEON FSI + if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so") + set(NEUT_NUCFSI 1) + LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED) -# Check for new versions of NEUT with NUCLEON FSI -if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so") - set(NEUT_NUCFSI 1) - LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED) + LIST(APPEND NEUT_ROOT_LIBS + ${NEUT_CLASS}/neutnucfsistep.so + ${NEUT_CLASS}/neutnucfsivert.so + ) + endif() + + if(${NEUT_VERSION} VERSION_LESS 5.4.0) + LIST(APPEND NEUT_ROOT_LIBS + ${NEUT_CLASS}/neutrootTreeSingleton.so) + endif() LIST(APPEND NEUT_ROOT_LIBS - ${NEUT_CLASS}/neutnucfsistep.so - ${NEUT_CLASS}/neutnucfsivert.so + ${NEUT_CLASS}/neutvtx.so + ${NEUT_CLASS}/neutfsipart.so + ${NEUT_CLASS}/neutpart.so + ${NEUT_CLASS}/neutvect.so ) -endif() -if(${NEUT_VERSION} VERSION_LESS 5.4.0) - LIST(APPEND NEUT_ROOT_LIBS - ${NEUT_CLASS}/neutrootTreeSingleton.so) + foreach(OBJ ${NEUT_ROOT_LIBS}) + LIST(APPEND EXTRA_SHAREDOBJS ${OBJ}) + endforeach() endif() - -LIST(APPEND NEUT_ROOT_LIBS - ${NEUT_CLASS}/neutvtx.so - ${NEUT_CLASS}/neutfsipart.so - ${NEUT_CLASS}/neutpart.so - ${NEUT_CLASS}/neutvect.so - ) - -foreach(OBJ ${NEUT_ROOT_LIBS}) - LIST(APPEND EXTRA_SHAREDOBJS ${OBJ}) -endforeach() diff --git a/cmake/NuWroSetup.cmake b/cmake/NuWroSetup.cmake index f649bff..af9d455 100644 --- a/cmake/NuWroSetup.cmake +++ b/cmake/NuWroSetup.cmake @@ -1,112 +1,111 @@ # 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 . ################################################################################ if(NOT NUWRO_INPUT_FILE STREQUAL "") if(NOT EXISTS ${NUWRO_INPUT_FILE}) cmessage(FATAL_ERROR "Expected -DNUWRO_INPUT_FILE to point to a valid input file. Cannot find: '${NUWRO_INPUT_FILE}'") endif() if(CMAKE_BUILD_TYPE MATCHES DEBUG) BuildROOTProject(NuWro_event1 ${NUWRO_INPUT_FILE} "event,vec,vect,particle,flags,params,line" STATIC) SET(ROOTLIBNAME "libNuWro_event1.a") else(CMAKE_BUILD_TYPE MATCHES RELEASE) BuildROOTProject(NuWro_event1 ${NUWRO_INPUT_FILE} "event,vec,vect,particle,flags,params,line" SHARED) SET(ROOTLIBNAME "libNuWro_event1.so") endif() ADD_CUSTOM_TARGET(NuWro_event1HeaderLink ALL COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_BINARY_DIR}/NuWro_event1/event.h ${CMAKE_BINARY_DIR}/NuWro_event1/event1.h DEPENDS NuWro_event1) LIST(APPEND EXTRA_CXX_FLAGS -D__NUWRO_ENABLED__) LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${CMAKE_BINARY_DIR}/NuWro_event1) LIST(APPEND EXTRA_LINK_DIRS ${CMAKE_BINARY_DIR}) LIST(APPEND EXTRA_LIBS NuWro_event1) LIST(APPEND PROJECTWIDE_EXTRA_DEPENDENCIES NuWro_event1HeaderLink ) install(TARGETS NuWro_event1 DESTINATION lib) SET(NUWRO_BUILT_FROM_FILE TRUE) else() if(NUWRO STREQUAL "") cmessage(FATAL_ERROR "Variable NUWRO is not defined. " "This must be set to point to a prebuilt NuWro instance.") endif() # If you are using a version of NuWro without reweighting use this to compile. if(USE_NuWro_RW) if(NUWRO_INC STREQUAL "") cmessage(FATAL_ERROR "Variable NUWRO_INC is not defined. " "This must be set to point to an installed NuWro instance.") endif() LIST(APPEND EXTRA_CXX_FLAGS -D__NUWRO_ENABLED__ -D__NUWRO_REWEIGHT_ENABLED__) if(USE_NuWro_SRW_Event) LIST(APPEND EXTRA_CXX_FLAGS -D__USE_NUWRO_SRW_EVENTS__) endif() LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NUWRO}/src ${NUWRO}/src/reweight ${NUWRO_INC}/nuwro) LIST(APPEND EXTRA_LINK_DIRS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib) LIST(APPEND EXTRA_LIBS reweight event) else () LIST(APPEND EXTRA_CXX_FLAGS -D__NUWRO_ENABLED__) LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NUWRO}/src) if(NOT EXISTS ${NUWRO}/bin/event1.so) if(EXISTS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib) if(NUWRO_INC STREQUAL "") cmessage(FATAL_ERROR "Variable NUWRO_INC is not defined. " "This must be set to point to an installed NuWro instance.") endif() LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NUWRO_INC}/nuwro) LIST(APPEND EXTRA_LINK_DIRS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib) - LIST(APPEND EXTRA_LIBS event) + LIST(APPEND EXTRA_LIBS -Wl,--no-as-needed event -Wl,--as-needed) else() cmessage(FATAL_ERROR "Expected to find the NuWro event library in: ${NUWRO}/bin/event1.so, or if using NuWro with reweight support: ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib/libevent.a. Is NuWro built?") endif() else() - LIST(APPEND EXTRA_SHAREDOBJS ${NUWRO}/bin/event1.so) + LIST(APPEND EXTRA_SHAREDOBJS -Wl,--no-as-needed ${NUWRO}/bin/event1.so -Wl,--as-needed) endif() endif() set(NEED_PYTHIA6 TRUE) set(NEED_ROOTPYTHIA6 TRUE) endif() - diff --git a/cmake/ReweightEnginesSetup.cmake b/cmake/ReweightEnginesSetup.cmake index f72d2ac..31dc580 100644 --- a/cmake/ReweightEnginesSetup.cmake +++ b/cmake/ReweightEnginesSetup.cmake @@ -1,93 +1,97 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ +if(NOT USE_REWEIGHT) + LIST(APPEND EXTRA_CXX_FLAGS -D__NO_REWEIGHT__) +endif() + ################################## T2K ###################################### if(USE_T2K) include(${CMAKE_SOURCE_DIR}/cmake/T2KSetup.cmake) cmessage(STATUS "Using T2K Reweight engine.") set(USE_T2K TRUE CACHE BOOL "Whether to enable T2KReWeight support. Requires external libraries. " FORCE) endif() ################################## NIWG ###################################### if(USE_NIWG) include(${CMAKE_SOURCE_DIR}/cmake/NIWGSetup.cmake) cmessage(STATUS "Using NIWG Reweight engine.") set(USE_NIWG TRUE CACHE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. " FORCE) endif() ################################## MINERvA ###################################### if(USE_MINERvA_RW) include(${CMAKE_SOURCE_DIR}/cmake/MINERvASetup.cmake) cmessage(STATUS "Using MINERvA Reweight engine.") set(USE_MINERvA_RW TRUE CACHE BOOL "Whether to enable MINERvA ReWeight support. " FORCE) endif() ################################## NEUT ###################################### if(USE_NEUT) include(${CMAKE_SOURCE_DIR}/cmake/NEUTSetup.cmake) cmessage(STATUS "Using NEUT Reweight engine.") set(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. " FORCE) endif() ################################# NuWro ###################################### if(USE_NuWro) include(${CMAKE_SOURCE_DIR}/cmake/NuWroSetup.cmake) cmessage(STATUS "Using NuWro Reweight engine.") set(USE_NuWro TRUE CACHE BOOL "Whether to enable NuWro support. " FORCE) endif() ################################## GENIE ##################################### if(USE_GENIE) include(${CMAKE_SOURCE_DIR}/cmake/GENIESetup.cmake) cmessage(STATUS "Using GENIE.") set(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE support. Requires external libraries. " FORCE) endif() ################################################################################ ################################ NOvARwgt #################################### if(USE_NOvARwgt) include(${CMAKE_SOURCE_DIR}/cmake/NOvARwgtSetup.cmake) cmessage(STATUS "Using NOvARwgt Reweight engine.") set(USE_NOvARwgt TRUE CACHE BOOL "Whether to enable NOvARwgt (reweight) support. Requires external libraries. " FORCE) endif() ################################################################################ ################################ Prob3++ #################################### include(${CMAKE_SOURCE_DIR}/cmake/Prob3++Setup.cmake) ################################################################################ cmessage(STATUS "Reweight engine include directories: ${RWENGINE_INCLUDE_DIRECTORIES}") if(NEED_ROOTEVEGEN) cmessage(STATUS "Require ROOT eve generation libraries") LIST(REVERSE ROOT_LIBS) LIST(APPEND ROOT_LIBS Gui Ged Geom TreePlayer EG Eve) LIST(REVERSE ROOT_LIBS) endif() if(NEED_ROOTPYTHIA6) cmessage(STATUS "Require ROOT Pythia6 libraries") LIST(APPEND ROOT_LIBS EGPythia6) endif() LIST(APPEND EXTRA_LIBS ${ROOT_LIBS}) diff --git a/cmake/parseConfigApp.cmake b/cmake/parseConfigApp.cmake index 2f72ce5..76238c0 100644 --- a/cmake/parseConfigApp.cmake +++ b/cmake/parseConfigApp.cmake @@ -1,148 +1,185 @@ function(GETFIRSTMATCHINGDELMIMMEDDIR DELIM CONFIGAPP ARG DIR_OUT FAILURE_IS_NOT_ERROR) - cmessage(DEBUG "CONFIGAPP: ${CONFIGAPP}, ARG: ${ARG} DIR_OUT: ${DIR_OUT}, FAILURE_IS_NOT_ERROR: ${FAILURE_IS_NOT_ERROR}") - if(DELIM STREQUAL "") cmessage(FATAL_ERROR "GETFIRSTMATCHINGDELMIMMEDDIR Passed no delimiter. This is a build configuration bug in NUISANCE, please report to the developers.") endif() if(CONFIGAPP STREQUAL "") cmessage(FATAL_ERROR "GETFIRSTMATCHINGDELMIMMEDDIR Passed no configuration application. This is a build configuration bug in NUISANCE, please report to the developers.") endif() SET(CONFIGAPP_LOCATION "CONFIGAPP_LOCATION-NOTFOUND") find_program(CONFIGAPP_LOCATION ${CONFIGAPP}) if(NOT CONFIGAPP_LOCATION STREQUAL "CONFIGAPP_LOCATION-NOTFOUND") execute_process (COMMAND ${CONFIGAPP_LOCATION} ${ARG} OUTPUT_VARIABLE CONFIGAPP_RESPONSE_RAW OUTPUT_STRIP_TRAILING_WHITESPACE) cmessage(DEBUG "${CONFIGAPP_LOCATION} ${ARG} responded with: \"${CONFIGAPP_RESPONSE_RAW}\"") if(CONFIGAPP_RESPONSE_RAW STREQUAL "") if(FAILURE_IS_NOT_ERROR) cmessage(DEBUG "\"${CONFIGAPP_LOCATION} ${ARG}\" produced no output and was expected to.") set(${DIR_OUT} "" PARENT_SCOPE) else() cmessage(FATAL_ERROR "\"${CONFIGAPP_LOCATION} ${ARG}\" produced no output and was required to.") endif() else() string(REGEX MATCH "${DELIM}\([^ ]+\)" PARSE_CONFIGAPP_RESPONSE_MATCH ${CONFIGAPP_RESPONSE_RAW}) if(NOT PARSE_CONFIGAPP_RESPONSE_MATCH) if(FAILURE_IS_NOT_ERROR) cmessage(DEBUG "Couldn't find ${DELIM} flag, found: \"${CONFIGAPP_RESPONSE_RAW}\"") set(${CMAKE_MATCH_1} "") else() cmessage(FATAL_ERROR "Expected to be able to parse the result of ${CONFIGAPP} ${ARG} to a lib directory, but couldn't find a ${DELIM} flag, found: \"${CONFIGAPP_RESPONSE_RAW}\"") endif() endif() set(${DIR_OUT} ${CMAKE_MATCH_1} PARENT_SCOPE) endif() else() cmessage(FATAL_ERROR "[ERROR]: Failed to find dependency configuration application: \"${CONFIGAPP}\"") endif() endfunction() #Uselike GETLIBDIR(gsl-config --libs GSL_LIB_DIR) function(GETLIBDIR CONFIGAPP ARG LIBDIR_OUT) if(ARGN) set(FAILURE_IS_NOT_ERROR TRUE) else() set(FAILURE_IS_NOT_ERROR FALSE) endif() GETFIRSTMATCHINGDELMIMMEDDIR( "-L" ${CONFIGAPP} ${ARG} MATCHING_DIR ${FAILURE_IS_NOT_ERROR}) set(${LIBDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE) endfunction() #Uselike GETINCDIR(gsl-config --cflags GSL_INC_DIR) function(GETINCDIR CONFIGAPP ARG INCDIR_OUT) if(ARGN) set(FAILURE_IS_NOT_ERROR TRUE) else() set(FAILURE_IS_NOT_ERROR FALSE) endif() GETFIRSTMATCHINGDELMIMMEDDIR( "-I" ${CONFIGAPP} ${ARG} MATCHING_DIR ${FAILURE_IS_NOT_ERROR}) set(${INCDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE) endfunction() function(GETALLMATCHINGDELMIMMEDDIR DELIM CONFIGAPP ARG LIST_OUT) if(DELIM STREQUAL "") cmessage(FATAL_ERROR "GETALLMATCHINGDELMIMMEDDIR Passed no delimiter. This is a build configuration bug in NUISANCE, please report to the developers.") endif() if(CONFIGAPP STREQUAL "") cmessage(FATAL_ERROR "GETALLMATCHINGDELMIMMEDDIR Passed no configuration application. This is a build configuration bug in NUISANCE, please report to the developers.") endif() - #cmessage(WARNING "DELIM: ${DELIM}, CONFIGAPP: ${CONFIGAPP}, ARG: ${ARG}, LIST_OUT: ${LIST_OUT}") - SET(CONFIGAPP_LOCATION "CONFIGAPP_LOCATION-NOTFOUND") find_program(CONFIGAPP_LOCATION ${CONFIGAPP}) if(NOT CONFIGAPP_LOCATION STREQUAL "CONFIGAPP_LOCATION-NOTFOUND") execute_process (COMMAND ${CONFIGAPP_LOCATION} ${ARG} OUTPUT_VARIABLE CONFIGAPP_RESPONSE_RAW OUTPUT_STRIP_TRAILING_WHITESPACE) string(REPLACE " " ";" CONFIGAPP_RESPONSE_LIST "${CONFIGAPP_RESPONSE_RAW}") - # cmessage(WARNING "CONFIGAPP_RESPONSE_RAW: ${CONFIGAPP_RESPONSE_RAW}, CONFIGAPP_RESPONSE_LIST: ${CONFIGAPP_RESPONSE_LIST}") - set(LIST_BUILD) foreach(I ${CONFIGAPP_RESPONSE_LIST}) if(I) string(REGEX MATCH "^${DELIM}" WASMATCHED ${I}) if(WASMATCHED) string(REPLACE "${DELIM}" "" I_STRIPPED "${I}") LIST(APPEND LIST_BUILD ${I_STRIPPED}) endif() endif() endforeach() set(${LIST_OUT} ${LIST_BUILD} PARENT_SCOPE) else() cmessage(FATAL_ERROR "[ERROR]: Failed to find dependency configuration application: \"${CONFIGAPP}\"") endif() endfunction() #Uselike GETLIBDIRS(gsl-config --libs GSL_LIB_DIR) function(GETLIBDIRS CONFIGAPP ARG LIBDIR_OUT) GETALLMATCHINGDELMIMMEDDIR( "-L" ${CONFIGAPP} ${ARG} MATCHING_DIR) set(${LIBDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE) endfunction() #Uselike GETINCDIRS(gsl-config --cflags GSL_INC_DIR) function(GETINCDIRS CONFIGAPP ARG INCDIR_OUT) GETALLMATCHINGDELMIMMEDDIR( "-I" ${CONFIGAPP} ${ARG} MATCHING_DIR) set(${INCDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE) endfunction() #Uselike GETLIBS(gsl-config --libs GSL_LIB_DIR) function(GETLIBS CONFIGAPP ARG LIBLIST_OUT) - #cmessage(WARNING "LIBLIST_OUT: ${LIBLIST_OUT}") GETALLMATCHINGDELMIMMEDDIR( "-l" ${CONFIGAPP} ${ARG} MATCHING_ITEMS) set(${LIBLIST_OUT} ${MATCHING_ITEMS} PARENT_SCOPE) - #cmessage(WARNING "LIBLIST_OUT: ${LIBLIST_OUT}: ${${LIBLIST_OUT}}") +endfunction() + +function(BuildFlagString OUTPUTLIST DELIMITER) + LIST(APPEND INPUT_LIST ${ARGN}) + if(NOT DELIMITER STREQUAL "") + string(STRIP ${DELIMITER} DELIMITER) + endif() + if(NOT "${INPUT_LIST}" STREQUAL "") + string(REPLACE ";" " ${DELIMITER}" LIST_STR "${DELIMITER}${INPUT_LIST}") + string(STRIP "${LIST_STR}" LIST_STR) + set(${OUTPUTLIST} ${LIST_STR} PARENT_SCOPE) + else() + set(${OUTPUTLIST} PARENT_SCOPE) + endif() +endfunction() + +#Here we want to look for linker flags in a list of libraries and escape them +function(BuildLibraryFlagString OUTPUTLIST) + set(INPUT_LIST ${ARGN}) + BuildFlagString(NEED_SCRUB_LINKOPTS "-l" ${INPUT_LIST}) + string(REPLACE "-l-" "-" LIST_STR ${NEED_SCRUB_LINKOPTS}) + set(${OUTPUTLIST} ${LIST_STR} PARENT_SCOPE) +endfunction() + +function(CatStringsIfNotEmpty OUTPUT_STRING) + set(BUILD_STR) + foreach(I ${ARGN}) + if(NOT I STREQUAL "") + set(BUILD_STR "${BUILD_STR} ${I}") + endif() + endforeach() + if(NOT BUILD_STR STREQUAL "") + string(STRIP "${BUILD_STR}" BUILD_STR) + endif() + set(${OUTPUT_STRING} ${BUILD_STR} PARENT_SCOPE) +endfunction() + +function(PrefixList OUTPUT_LIST PREFIX) + set(${OUTPUT_LIST}) + foreach(I ${ARGN}) + if(NOT I STREQUAL "") + LIST(APPEND ${OUTPUT_LIST} "${PREFIX}${I}") + endif() + endforeach() + set(${OUTPUT_LIST} ${${OUTPUT_LIST}} PARENT_SCOPE) endfunction() diff --git a/src/FitBase/Measurement2D.cxx b/src/FitBase/Measurement2D.cxx index 37316bc..0cc89f3 100644 --- a/src/FitBase/Measurement2D.cxx +++ b/src/FitBase/Measurement2D.cxx @@ -1,1982 +1,2018 @@ // 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 "Measurement2D.h" #include "TDecompChol.h" //******************************************************************** Measurement2D::Measurement2D(void) { //******************************************************************** covar = NULL; fDecomp = NULL; fFullCovar = NULL; fMCHist = NULL; fMCFine = NULL; fDataHist = NULL; fMCHist_X = NULL; fMCHist_Y = NULL; fDataHist_X = NULL; fDataHist_Y = NULL; fMaskHist = NULL; fMapHist = NULL; fDataOrig = NULL; fDataTrue = NULL; fMCWeighted = NULL; + fResidualHist = NULL; + fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/FITPROJX/" "FITPROJY"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsDifXSec = false; fIsEnu = false; // XSec Scalings fScaleFactor = -1.0; fCurrentNorm = 1.0; // Histograms fDataHist = NULL; fDataTrue = NULL; fMCHist = NULL; fMCFine = NULL; fMCWeighted = NULL; fMaskHist = NULL; // Covar covar = NULL; fFullCovar = NULL; fCovar = NULL; fInvert = NULL; fDecomp = NULL; // Fake Data fFakeDataInput = ""; fFakeDataFile = NULL; // Options fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsDifXSec = false; fIsEnu1D = false; // Inputs fInput = NULL; fRW = NULL; // Extra Histograms fMCHist_Modes = NULL; } //******************************************************************** Measurement2D::~Measurement2D(void) { //******************************************************************** if (fDataHist) delete fDataHist; if (fDataTrue) delete fDataTrue; if (fMCHist) delete fMCHist; if (fMCFine) delete fMCFine; if (fMCWeighted) delete fMCWeighted; if (fMaskHist) delete fMaskHist; if (covar) delete covar; if (fFullCovar) delete fFullCovar; if (fCovar) delete fCovar; if (fInvert) delete fInvert; if (fDecomp) delete fDecomp; + + delete fResidualHist; } //******************************************************************** void Measurement2D::FinaliseSampleSettings() { //******************************************************************** MeasurementBase::FinaliseSampleSettings(); // Setup naming + renaming fName = fSettings.GetName(); fSettings.SetS("originalname", fName); if (fSettings.Has("rename")) { fName = fSettings.GetS("rename"); fSettings.SetS("name", fName); } // Setup all other options NUIS_LOG(SAM, "Finalising Sample Settings: " << fName); if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods."); } if (fSettings.GetS("originalname").find("Enu") != std::string::npos) { fIsEnu1D = true; NUIS_LOG(SAM, "::" << fName << "::"); NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " << "not flux averaged!"); } if (fIsEnu1D && fIsRawEvents) { NUIS_ERR(FTL, "Found 2D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__); } if (!fRW) fRW = FitBase::GetRW(); if (!fInput) SetupInputs(fSettings.GetS("input")); // Setup options SetFitOptions(fDefaultTypes); // defaults SetFitOptions(fSettings.GetS("type")); // user specified EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min")); EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max")); if (fAddNormPen) { fNormError = fSettings.GetNormError(); if (fNormError <= 0.0) { NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!"); NUIS_ABORT("If you want to use it please add fNormError=VAL"); } } } void Measurement2D::CreateDataHistogram(int dimx, double *binx, int dimy, double *biny) { if (fDataHist) delete fDataHist; NUIS_LOG(SAM, "Creating Data Histogram dim : " << dimx << " " << dimy); fDataHist = new TH2D((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(), dimx - 1, binx, dimy - 1, biny); } void Measurement2D::SetDataFromTextFile(std::string data, std::string binx, std::string biny) { // Get the data hist fDataHist = PlotUtils::GetTH2DFromTextFile(data, binx, biny); // Set the name properly fDataHist->SetName((fSettings.GetName() + "_data").c_str()); fDataHist->SetTitle(fSettings.GetFullTitles().c_str()); } void Measurement2D::SetDataFromRootFile(std::string datfile, std::string histname) { NUIS_LOG(SAM, "Reading data from root file: " << datfile << ";" << histname); fDataHist = PlotUtils::GetTH2DFromRootFile(datfile, histname); fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str()); } void Measurement2D::SetDataValuesFromTextFile(std::string datfile, TH2D *hist) { NUIS_LOG(SAM, "Setting data values from text file"); if (!hist) hist = fDataHist; // Read TH2D From textfile TH2D *valhist = (TH2D *)hist->Clone(); valhist->Reset(); PlotUtils::Set2DHistFromText(datfile, valhist, 1.0, true); NUIS_LOG(SAM, " -> Filling values from read hist."); for (int i = 0; i < valhist->GetNbinsX(); i++) { for (int j = 0; j < valhist->GetNbinsY(); j++) { hist->SetBinContent(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1)); } } NUIS_LOG(SAM, " --> Done"); } void Measurement2D::SetDataErrorsFromTextFile(std::string datfile, TH2D *hist) { NUIS_LOG(SAM, "Setting data errors from text file"); if (!hist) hist = fDataHist; // Read TH2D From textfile TH2D *valhist = (TH2D *)hist->Clone(); valhist->Reset(); PlotUtils::Set2DHistFromText(datfile, valhist, 1.0); // Fill Errors NUIS_LOG(SAM, " -> Filling errors from read hist."); for (int i = 0; i < valhist->GetNbinsX(); i++) { for (int j = 0; j < valhist->GetNbinsY(); j++) { hist->SetBinError(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1)); } } NUIS_LOG(SAM, " --> Done"); } void Measurement2D::SetMapValuesFromText(std::string dataFile) { TH2D *hist = fDataHist; std::vector edgex; std::vector edgey; for (int i = 0; i <= hist->GetNbinsX(); i++) edgex.push_back(hist->GetXaxis()->GetBinLowEdge(i + 1)); for (int i = 0; i <= hist->GetNbinsY(); i++) edgey.push_back(hist->GetYaxis()->GetBinLowEdge(i + 1)); fMapHist = new TH2I((fName + "_map").c_str(), (fName + fPlotTitles).c_str(), edgex.size() - 1, &edgex[0], edgey.size() - 1, &edgey[0]); NUIS_LOG(SAM, "Reading map from: " << dataFile); PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0); } //******************************************************************** void Measurement2D::SetPoissonErrors() { //******************************************************************** if (!fDataHist) { NUIS_ERR(FTL, "Need a data hist to setup possion errors! "); NUIS_ABORT("Setup Data First!"); } for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) { fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1))); } } //******************************************************************** void Measurement2D::SetCovarFromDiagonal(TH2D *data) { //******************************************************************** if (!data and fDataHist) { data = fDataHist; } if (data) { NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName()); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } else { NUIS_ABORT("No data input provided to set diagonal covar from!"); } // if (!fIsDiag) { // ERR(FTL) << "SetCovarMatrixFromDiag called for measurement " // << "that is not set as diagonal." << std::endl; // throw; // } } //******************************************************************** void Measurement2D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = this->GetNDOF(); } NUIS_LOG(SAM, "Reading covariance from text file: " << covfile << " " << dim); fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCovarFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading covariance from text file: " << covfile << ";" << histname); fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCovarInvertFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = this->GetNDOF(); } NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile); covar = StatUtils::GetCovarFromTextFile(covfile, dim); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";" << histname); covar = StatUtils::GetCovarFromRootFile(covfile, histname); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCorrelationFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) dim = this->GetNDOF(); NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << dim); TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim); if (!fDataHist) { NUIS_ABORT("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"); } // 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 Measurement2D::SetCorrelationFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << histname); TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname); if (!fDataHist) { NUIS_ABORT("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"); } // 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 Measurement2D::SetCholDecompFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = this->GetNDOF(); } NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile << " " << dim); TMatrixD *temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement2D::SetCholDecompFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";" << histname); TMatrixD *temp = StatUtils::GetMatrixFromRootFile(covfile, histname); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement2D::ScaleData(double scale) { //******************************************************************** fDataHist->Scale(scale); } //******************************************************************** void Measurement2D::ScaleDataErrors(double scale) { //******************************************************************** for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsY(); j++) { fDataHist->SetBinError(i + 1, j + 1, fDataHist->GetBinError(i + 1, j + 1) * scale); } } } //******************************************************************** void Measurement2D::ScaleCovar(double scale) { //******************************************************************** (*fFullCovar) *= scale; (*covar) *= 1.0 / scale; (*fDecomp) *= sqrt(scale); } //******************************************************************** void Measurement2D::SetBinMask(std::string maskfile) { //******************************************************************** if (!fIsMask) return; NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile); // Create a mask histogram with dim of data int nbinsx = fDataHist->GetNbinsX(); int nbinxy = fDataHist->GetNbinsY(); fMaskHist = new TH2I((fSettings.GetName() + "_BINMASK").c_str(), (fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbinsx, 0, nbinsx, nbinxy, 0, nbinxy); std::string line; std::ifstream mask(maskfile.c_str(), std::ifstream::in); if (!mask.is_open()) { NUIS_LOG(FTL, " Cannot find mask file."); 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) { NUIS_LOG(WRN, "Measurement2D::SetBinMask(), couldn't parse line: " << line); continue; } // The first index should be the bin number, the second should be the mask // value. int val = 0; if (entries[2] > 0) val = 1; fMaskHist->SetBinContent(entries[0], entries[1], val); } // Apply masking by setting masked data bins to zero PlotUtils::MaskBins(fDataHist, fMaskHist); return; } //******************************************************************** void Measurement2D::FinaliseMeasurement() { //******************************************************************** NUIS_LOG(SAM, "Finalising Measurement: " << fName); if (fSettings.GetB("onlymc")) { if (fDataHist) delete fDataHist; fDataHist = new TH2D("empty_data", "empty_data", 1, 0.0, 1.0, 1, 0.0, 1.0); } // Make sure data is setup if (!fDataHist) { NUIS_ABORT("No data has been setup inside " << fName << " constructor!"); } // Make sure covariances are setup if (!fFullCovar) { fIsDiag = true; SetCovarFromDiagonal(fDataHist); } if (!covar) { covar = StatUtils::GetInvert(fFullCovar); } if (!fDecomp) { fDecomp = StatUtils::GetDecomp(fFullCovar); } // Setup fMCHist from data fMCHist = (TH2D *)fDataHist->Clone(); fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(), (fSettings.GetFullTitles()).c_str()); fMCHist->Reset(); // Setup fMCFine fMCFine = new TH2D( "mcfine", "mcfine", fDataHist->GetNbinsX() * 6, fMCHist->GetXaxis()->GetBinLowEdge(1), fMCHist->GetXaxis()->GetBinLowEdge(fDataHist->GetNbinsX() + 1), fDataHist->GetNbinsY() * 6, fMCHist->GetYaxis()->GetBinLowEdge(1), fMCHist->GetYaxis()->GetBinLowEdge(fDataHist->GetNbinsY() + 1)); fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(), (fSettings.GetFullTitles()).c_str()); fMCFine->Reset(); // Setup MC Stat fMCStat = (TH2D *)fMCHist->Clone(); fMCStat->Reset(); // Search drawopts for possible types to include by default std::string drawopts = FitPar::Config().GetParS("drawopts"); if (drawopts.find("MODES") != std::string::npos) { fMCHist_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(), ("True Channels"), fMCHist); SetAutoProcessTH1(fMCHist_Modes); } // Setup bin masks using sample name if (fIsMask) { std::string curname = fName; std::string origname = fSettings.GetS("originalname"); // Check rename.mask std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask"); // Check origname.mask if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask"); // Check database if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask"; } // Setup Bin Mask SetBinMask(maskloc); } if (fScaleFactor < 0) { NUIS_ERR(FTL, "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__); NUIS_ERR(FTL, "fScaleFactor = " << fScaleFactor); NUIS_ABORT("EXITING"); } // Create and fill Weighted Histogram if (!fMCWeighted) { fMCWeighted = (TH2D *)fMCHist->Clone(); fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(), (fName + "_MCWGHTS" + fPlotTitles).c_str()); fMCWeighted->GetYaxis()->SetTitle("Weighted Events"); } if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist); } //******************************************************************** void Measurement2D::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) { NUIS_ABORT("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."); } } } // 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) { NUIS_ERR(FTL, "ERROR: Fit Option '" << fit_options_input.at(i) << "' Provided is not allowed for this measurement."); NUIS_ERR(FTL, "Fit Options should be provided as a '/' seperated list " "(e.g. FREE/DIAG/NORM)"); NUIS_ABORT("Available options for " << fName << " are '" << fAllowedTypes << "'"); } } // 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; NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!"); NUIS_ABORT("Try to use a chi2!"); } else { fIsChi2 = true; } // EXTRAS if (opt.find("RAW") != std::string::npos) fIsRawEvents = true; if (opt.find("DIF") != std::string::npos) fIsDifXSec = true; if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true; if (opt.find("NORM") != std::string::npos) fAddNormPen = true; if (opt.find("MASK") != std::string::npos) fIsMask = true; // Set TYPE fFitType = opt; // FIX,SHAPE,FREE if (opt.find("FIX") != std::string::npos) { fIsFree = fIsShape = false; fIsFix = true; } else if (opt.find("SHAPE") != std::string::npos) { fIsFree = fIsFix = false; fIsShape = true; } else if (opt.find("FREE") != std::string::npos) { fIsFix = fIsShape = false; fIsFree = true; } // DIAG,FULL (or default to full) if (opt.find("DIAG") != std::string::npos) { fIsDiag = true; fIsFull = false; } else if (opt.find("FULL") != std::string::npos) { fIsDiag = false; fIsFull = true; } // CHI2/LL (OTHERS?) if (opt.find("LOG") != std::string::npos) fIsChi2 = false; else fIsChi2 = true; // EXTRAS if (opt.find("RAW") != std::string::npos) fIsRawEvents = true; if (opt.find("DIF") != std::string::npos) fIsDifXSec = true; if (opt.find("ENU1D") != std::string::npos) fIsEnu = true; if (opt.find("NORM") != std::string::npos) fAddNormPen = true; if (opt.find("MASK") != std::string::npos) fIsMask = true; fIsProjFitX = (opt.find("FITPROJX") != std::string::npos); fIsProjFitY = (opt.find("FITPROJY") != std::string::npos); return; }; /* Reconfigure LOOP */ //******************************************************************** void Measurement2D::ResetAll() { //******************************************************************** fMCHist->Reset(); fMCFine->Reset(); fMCStat->Reset(); return; }; //******************************************************************** void Measurement2D::FillHistograms() { //******************************************************************** if (Signal) { fMCHist->Fill(fXVar, fYVar, Weight); fMCFine->Fill(fXVar, fYVar, Weight); fMCStat->Fill(fXVar, fYVar, 1.0); if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, fYVar, Weight); } return; }; //******************************************************************** void Measurement2D::ScaleEvents() { //******************************************************************** // Fill MCWeighted; // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1)); // fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1)); // } // Setup Stat ratios for MC and MC Fine double *statratio = new double[fMCHist->GetNbinsX()]; for (int i = 0; i < fMCHist->GetNbinsX(); i++) { if (fMCHist->GetBinContent(i + 1) != 0) { statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1); } else { statratio[i] = 0.0; } } double *statratiofine = new double[fMCFine->GetNbinsX()]; for (int i = 0; i < fMCFine->GetNbinsX(); i++) { if (fMCFine->GetBinContent(i + 1) != 0) { statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1); } else { statratiofine[i] = 0.0; } } // Scaling for raw event rates if (fIsRawEvents) { double datamcratio = fDataHist->Integral() / fMCHist->Integral(); fMCHist->Scale(datamcratio); fMCFine->Scale(datamcratio); if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio); // Scaling for XSec as function of Enu } else if (fIsEnu1D) { PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(), GetEventHistogram(), fScaleFactor); PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(), GetEventHistogram(), fScaleFactor); // if (fMCHist_Modes) { // PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(), // GetEventHistogram(), fScaleFactor, // fNEvents); // } // Any other differential scaling } else { fMCHist->Scale(fScaleFactor, "width"); fMCFine->Scale(fScaleFactor, "width"); // if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width"); } // Proper error scaling - ROOT Freaks out with xsec weights sometimes for (int i = 0; i < fMCStat->GetNbinsX(); i++) { fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]); } for (int i = 0; i < fMCFine->GetNbinsX(); i++) { fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]); } // Clean up delete statratio; delete statratiofine; return; }; //******************************************************************** void Measurement2D::ApplyNormScale(double norm) { //******************************************************************** fCurrentNorm = norm; fMCHist->Scale(1.0 / norm); fMCFine->Scale(1.0 / norm); return; }; /* Statistic Functions - Outsources to StatUtils */ //******************************************************************** int Measurement2D::GetNDOF() { //******************************************************************** // Just incase it has gone... if (!fDataHist) return -1; int nDOF = 0; // If datahist has no errors make sure we don't include those bins as they are // not data points for (int xBin = 0; xBin < fDataHist->GetNbinsX(); ++xBin) { for (int yBin = 0; yBin < fDataHist->GetNbinsY(); ++yBin) { if (fDataHist->GetBinError(xBin + 1, yBin + 1) != 0) ++nDOF; } } // Account for possible bin masking int nMasked = 0; if (fMaskHist and fIsMask) if (fMaskHist->Integral() > 0) for (int xBin = 0; xBin < fMaskHist->GetNbinsX() + 1; ++xBin) for (int yBin = 0; yBin < fMaskHist->GetNbinsY() + 1; ++yBin) if (fMaskHist->GetBinContent(xBin, yBin) > 0.5) ++nMasked; // Take away those masked DOF if (fIsMask) { nDOF -= nMasked; } return nDOF; } //******************************************************************** double Measurement2D::GetLikelihood() { //******************************************************************** // If this is for a ratio, there is no data histogram to compare to! if (fNoData || !fDataHist) return 0.; // Fix weird masking bug if (!fIsMask) { if (fMaskHist) { fMaskHist = NULL; } } else { if (fMaskHist) { PlotUtils::MaskBins(fMCHist, fMaskHist); } } // if (fIsProjFitX or fIsProjFitY) return GetProjectedChi2(); // Scale up the results to match each other (Not using width might be // inconsistent with Meas1D) double scaleF = fDataHist->Integral() / fMCHist->Integral(); if (fIsShape) { fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); // PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF); } if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist); // Get the chi2 from either covar or diagonals double chi2 = 0.0; if (fIsChi2) { if (fIsDiag) { chi2 = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMapHist, fMaskHist); } else { chi2 = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist, - fMaskHist); + fMaskHist, fResidualHist); + if (fChi2LessBinHist) { + TH2I *binmask = new TH2I("mask", "", fDataHist->GetNbinsX(), 0, + fDataHist->GetNbinsX(), fDataHist->GetNbinsY(), + 0, fDataHist->GetNbinsY()); + binmask->SetDirectory(NULL); + for (int xi = 0; xi < fDataHist->GetNbinsX(); ++xi) { + for (int yi = 0; yi < fDataHist->GetNbinsY(); ++yi) { + binmask->Reset(); + binmask->SetBinContent(xi + 1, yi + 1, 1); + fChi2LessBinHist->SetBinContent( + xi + 1, yi + 1, + StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist, + binmask)); + } + } + delete binmask; + } } } // Add a normal penalty term if (fAddNormPen) { chi2 += (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError); NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError)); } // Adjust the shape back to where it was. if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) { fMCHist->Scale(1. / scaleF); fMCFine->Scale(1. / scaleF); } fLikelihood = chi2; return chi2; } /* Fake Data Functions */ //******************************************************************** void Measurement2D::SetFakeDataValues(std::string fakeOption) { //******************************************************************** // Setup original/datatrue TH2D *tempdata = (TH2D *)fDataHist->Clone(); if (!fIsFakeData) { fIsFakeData = true; // Make a copy of the original data histogram. if (!fDataOrig) fDataOrig = (TH2D *)fDataHist->Clone((fName + "_data_original").c_str()); } else { ResetFakeData(); } // Setup Inputs fFakeDataInput = fakeOption; NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput); // From MC if (fFakeDataInput.compare("MC") == 0) { fDataHist = (TH2D *)fMCHist->Clone((fName + "_MC").c_str()); // Fake File } else { if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ"); fDataHist = (TH2D *)fFakeDataFile->Get((fName + "_MC").c_str()); } // Setup Data Hist fDataHist->SetNameTitle((fName + "_FAKE").c_str(), (fName + fPlotTitles).c_str()); // Replace Data True if (fDataTrue) delete fDataTrue; fDataTrue = (TH2D *)fDataHist->Clone(); fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(), (fName + fPlotTitles).c_str()); // Make a new covariance for fake data hist. int nbins = fDataHist->GetNbinsX() * fDataHist->GetNbinsY(); double alpha_i = 0.0; double alpha_j = 0.0; for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { if (tempdata->GetBinContent(i + 1) && tempdata->GetBinContent(j + 1)) { alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1); alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1); } else { alpha_i = 0.0; alpha_j = 0.0; } (*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j); } } // Setup Covariances if (covar) delete covar; covar = StatUtils::GetInvert(fFullCovar); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetInvert(fFullCovar); delete tempdata; return; }; //******************************************************************** void Measurement2D::ResetFakeData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH2D *)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str()); } } //******************************************************************** void Measurement2D::ResetData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH2D *)fDataOrig->Clone((fSettings.GetName() + "_data").c_str()); } fIsFakeData = false; } //******************************************************************** void Measurement2D::ThrowCovariance() { //******************************************************************** // Take a fDecomposition and use it to throw the current dataset. // Requires fDataTrue also be set incase used repeatedly. if (fDataHist) delete fDataHist; fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); return; }; //******************************************************************** void Measurement2D::ThrowDataToy() { //******************************************************************** if (!fDataTrue) fDataTrue = (TH2D *)fDataHist->Clone(); if (fMCHist) delete fMCHist; fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); } /* Access Functions */ //******************************************************************** TH2D *Measurement2D::GetMCHistogram() { //******************************************************************** if (!fMCHist) return fMCHist; std::ostringstream chi2; chi2 << std::setprecision(5) << this->GetLikelihood(); int linecolor = kRed; int linestyle = 1; int linewidth = 1; int fillcolor = 0; int fillstyle = 1001; if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor"); if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle"); if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth"); if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor"); if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle"); fMCHist->SetTitle(chi2.str().c_str()); fMCHist->SetLineColor(linecolor); fMCHist->SetLineStyle(linestyle); fMCHist->SetLineWidth(linewidth); fMCHist->SetFillColor(fillcolor); fMCHist->SetFillStyle(fillstyle); return fMCHist; }; //******************************************************************** TH2D *Measurement2D::GetDataHistogram() { //******************************************************************** if (!fDataHist) return fDataHist; int datacolor = kBlack; int datastyle = 1; int datawidth = 1; if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor"); if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle"); if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth"); fDataHist->SetLineColor(datacolor); fDataHist->SetLineWidth(datawidth); fDataHist->SetMarkerStyle(datastyle); return fDataHist; }; /* Write Functions */ // Save all the histograms at once //******************************************************************** void Measurement2D::Write(std::string drawOpt) { //******************************************************************** // Get Draw Options drawOpt = FitPar::Config().GetParS("drawopts"); // Write 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(); } + if (fIsChi2 && !fIsDiag) { + fResidualHist = (TH2D *)fMCHist->Clone((fName + "_RESIDUAL").c_str()); + fResidualHist->GetYaxis()->SetTitle("#Delta#chi^{2}"); + fResidualHist->Reset(); + + fChi2LessBinHist = (TH2D *)fMCHist->Clone((fName + "_Chi2NMinusOne").c_str()); + fChi2LessBinHist->GetYaxis()->SetTitle("Total #chi^{2} without bin_{i}"); + fChi2LessBinHist->Reset(); + + (void)GetLikelihood(); + + fResidualHist->Write(); + fChi2LessBinHist->Write(); + } + // // Likelihood residual plots // if (drawOpt.find("RESIDUAL") != std::string::npos) { // WriteResidualPlots(); //} // // RATIO // if (drawOpt.find("CANVMC") != std::string::npos) { // TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist); // c1->Write(); // delete c1; // } // // PDG // if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) { // TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes); // c2->Write(); // delete c2; // } /// 2D VERSION // If null pointer return if (!fMCHist and !fDataHist) { NUIS_LOG(SAM, fName << "Incomplete histogram set!"); return; } // Config::Get().out->cd(); // Get Draw Options drawOpt = FitPar::Config().GetParS("drawopts"); bool drawData = (drawOpt.find("DATA") != std::string::npos); bool drawNormal = (drawOpt.find("MC") != std::string::npos); bool drawEvents = (drawOpt.find("EVT") != std::string::npos); bool drawFine = (drawOpt.find("FINE") != std::string::npos); bool drawRatio = (drawOpt.find("RATIO") != std::string::npos); // bool drawModes = (drawOpt.find("MODES") != std::string::npos); bool drawShape = (drawOpt.find("SHAPE") != std::string::npos); bool residual = (drawOpt.find("RESIDUAL") != std::string::npos); bool drawMatrix = (drawOpt.find("MATRIX") != std::string::npos); bool drawFlux = (drawOpt.find("FLUX") != std::string::npos); bool drawMask = (drawOpt.find("MASK") != std::string::npos); bool drawMap = (drawOpt.find("MAP") != std::string::npos); bool drawProj = (drawOpt.find("PROJ") != std::string::npos); // bool drawCanvPDG = (drawOpt.find("CANVPDG") != std::string::npos); bool drawCov = (drawOpt.find("COV") != std::string::npos); bool drawSliceMC = (drawOpt.find("CANVSLICEMC") != std::string::npos); bool drawWeighted = (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted); if (FitPar::Config().GetParB("EventManager")) { drawFlux = false; drawEvents = false; } // Save standard plots if (drawData) { GetDataList().at(0)->Write(); // Generate a simple map if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist); // Convert to 1D Lists TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist); data_1D->Write(); delete data_1D; } if (drawNormal) { GetMCList().at(0)->Write(); if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist); TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist); mc_1D->SetLineColor(kRed); mc_1D->Write(); delete mc_1D; } // Write Weighted Histogram if (drawWeighted) fMCWeighted->Write(); if (drawCov) { TH2D(*fFullCovar).Write((fName + "_COV").c_str()); } if (drawOpt.find("INVCOV") != std::string::npos) { TH2D(*covar).Write((fName + "_INVCOV").c_str()); } // Save only mc and data if splines if (fEventType == 4 or fEventType == 3) { return; } // Draw Extra plots if (drawFine) this->GetFineList().at(0)->Write(); if (drawFlux and GetFluxHistogram()) { GetFluxHistogram()->Write(); } if (drawEvents and GetEventHistogram()) { GetEventHistogram()->Write(); } if (fIsMask and drawMask) { fMaskHist->Write((fName + "_MSK").c_str()); //< save mask TH1I *mask_1D = StatUtils::MapToMask(fMaskHist, fMapHist); if (mask_1D) { mask_1D->Write(); TMatrixDSym *calc_cov = StatUtils::ApplyInvertedMatrixMasking(covar, mask_1D); TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist); TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist); TH1D *calc_data = StatUtils::ApplyHistogramMasking(data_1D, mask_1D); TH1D *calc_mc = StatUtils::ApplyHistogramMasking(mc_1D, mask_1D); TH2D *bin_cov = new TH2D(*calc_cov); bin_cov->Write(); calc_data->Write(); calc_mc->Write(); delete mask_1D; delete calc_cov; delete calc_data; delete calc_mc; delete bin_cov; delete data_1D; delete mc_1D; } } if (drawMap) fMapHist->Write((fName + "_MAP").c_str()); //< save map // // Save neut stack // if (drawModes) { // THStack combo_fMCHist_PDG = PlotUtils::GetNeutModeStack( // (fName + "_MC_PDG").c_str(), // (TH1**)fMCHist_PDG, 0); // combo_fMCHist_PDG.Write(); // } // Save Matrix plots if (drawMatrix and fFullCovar and covar and fDecomp) { TH2D cov = TH2D((*fFullCovar)); cov.SetNameTitle((fName + "_cov").c_str(), (fName + "_cov;Bins; Bins;").c_str()); cov.Write(); TH2D covinv = TH2D((*this->covar)); covinv.SetNameTitle((fName + "_covinv").c_str(), (fName + "_cov;Bins; Bins;").c_str()); covinv.Write(); TH2D covdec = TH2D((*fDecomp)); covdec.SetNameTitle((fName + "_covdec").c_str(), (fName + "_cov;Bins; Bins;").c_str()); covdec.Write(); } // Save ratio plots if required if (drawRatio) { // Needed for error bars for (int i = 0; i < fMCHist->GetNbinsX() * fMCHist->GetNbinsY(); i++) fMCHist->SetBinError(i + 1, 0.0); fDataHist->GetSumw2(); fMCHist->GetSumw2(); // Create Ratio Histograms TH2D *dataRatio = (TH2D *)fDataHist->Clone((fName + "_data_RATIO").c_str()); TH2D *mcRatio = (TH2D *)fMCHist->Clone((fName + "_MC_RATIO").c_str()); mcRatio->Divide(fMCHist); dataRatio->Divide(fMCHist); // Cancel bin errors on MC for (int i = 0; i < mcRatio->GetNbinsX() * mcRatio->GetNbinsY(); i++) { mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); } mcRatio->SetMinimum(0); mcRatio->SetMaximum(2); dataRatio->SetMinimum(0); dataRatio->SetMaximum(2); mcRatio->Write(); dataRatio->Write(); delete mcRatio; delete dataRatio; } // Save Shape Plots if required if (drawShape) { // Create Shape Histogram TH2D *mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); mcShape->SetLineWidth(3); mcShape->SetLineStyle(7); // dashes mcShape->Write(); // Save shape ratios if (drawRatio) { // Needed for error bars mcShape->GetSumw2(); // Create shape ratio histograms TH2D *mcShapeRatio = (TH2D *)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str()); TH2D *dataShapeRatio = (TH2D *)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str()); // Divide the histograms mcShapeRatio->Divide(mcShape); dataShapeRatio->Divide(mcShape); // Colour the shape ratio plots mcShapeRatio->SetLineWidth(3); mcShapeRatio->SetLineStyle(7); // dashes mcShapeRatio->Write(); dataShapeRatio->Write(); delete mcShapeRatio; delete dataShapeRatio; } delete mcShape; } // Save residual calculations of what contributed to the chi2 values. if (residual) { } if (fIsProjFitX || fIsProjFitY || drawProj) { // If not already made, make the projections if (!fMCHist_X) { PlotUtils::MatchEmptyBins(fDataHist, fMCHist); fMCHist_X = PlotUtils::GetProjectionX(fMCHist, fMaskHist); fMCHist_Y = PlotUtils::GetProjectionY(fMCHist, fMaskHist); fDataHist_X = PlotUtils::GetProjectionX(fDataHist, fMaskHist); fDataHist_Y = PlotUtils::GetProjectionY(fDataHist, fMaskHist); // This is not the correct way of doing it // double chi2X = StatUtils::GetChi2FromDiag(fDataHist_X, fMCHist_X); // double chi2Y = StatUtils::GetChi2FromDiag(fDataHist_Y, fMCHist_Y); // fMCHist_X->SetTitle(Form("%f", chi2X)); // fMCHist_Y->SetTitle(Form("%f", chi2Y)); } // Save the histograms fDataHist_X->Write(); fMCHist_X->Write(); fDataHist_Y->Write(); fMCHist_Y->Write(); } if (drawSliceMC) { TCanvas *c1 = new TCanvas((fName + "_MC_CANV_Y").c_str(), (fName + "_MC_CANV_Y").c_str(), 1024, 1024); c1->Divide(2, int(fDataHist->GetNbinsY() / 3. + 1)); TH2D *mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); double shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); mcShape->Scale(shapeScale); mcShape->SetLineStyle(7); c1->cd(1); TLegend *leg = new TLegend(0.0, 0.0, 1.0, 1.0); leg->AddEntry(fDataHist, (fName + " Data").c_str(), "lep"); leg->AddEntry(fMCHist, (fName + " MC").c_str(), "l"); leg->AddEntry(mcShape, (fName + " Shape").c_str(), "l"); leg->SetLineColor(0); leg->SetLineStyle(0); leg->SetFillColor(0); leg->SetLineStyle(0); leg->Draw("SAME"); // Make Y slices for (int i = 1; i < fDataHist->GetNbinsY() + 1; i++) { c1->cd(i + 1); TH1D *fDataHist_SliceY = PlotUtils::GetSliceY(fDataHist, i); fDataHist_SliceY->Draw("E1"); TH1D *fMCHist_SliceY = PlotUtils::GetSliceY(fMCHist, i); fMCHist_SliceY->Draw("SAME"); TH1D *mcShape_SliceY = PlotUtils::GetSliceY(mcShape, i); mcShape_SliceY->Draw("SAME"); mcShape_SliceY->SetLineStyle(mcShape->GetLineStyle()); } c1->Write(); delete c1; delete leg; TCanvas *c2 = new TCanvas((fName + "_MC_CANV_X").c_str(), (fName + "_MC_CANV_X").c_str(), 1024, 1024); c2->Divide(2, int(fDataHist->GetNbinsX() / 3. + 1)); mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); mcShape->Scale(shapeScale); mcShape->SetLineStyle(7); c2->cd(1); TLegend *leg2 = new TLegend(0.0, 0.0, 1.0, 1.0); leg2->AddEntry(fDataHist, (fName + " Data").c_str(), "lep"); leg2->AddEntry(fMCHist, (fName + " MC").c_str(), "l"); leg2->AddEntry(mcShape, (fName + " Shape").c_str(), "l"); leg2->SetLineColor(0); leg2->SetLineStyle(0); leg2->SetFillColor(0); leg2->SetLineStyle(0); leg2->Draw("SAME"); // Make Y slices for (int i = 1; i < fDataHist->GetNbinsX() + 1; i++) { c2->cd(i + 1); TH1D *fDataHist_SliceX = PlotUtils::GetSliceX(fDataHist, i); fDataHist_SliceX->Draw("E1"); TH1D *fMCHist_SliceX = PlotUtils::GetSliceX(fMCHist, i); fMCHist_SliceX->Draw("SAME"); TH1D *mcShape_SliceX = PlotUtils::GetSliceX(mcShape, i); mcShape_SliceX->Draw("SAME"); mcShape_SliceX->SetLineStyle(mcShape->GetLineStyle()); } c2->Write(); delete c2; delete leg2; } // Write Extra Histograms AutoWriteExtraTH1(); WriteExtraHistograms(); // Returning NUIS_LOG(SAM, "Written Histograms: " << fName); return; } /* Setup Functions */ //******************************************************************** void Measurement2D::SetupMeasurement(std::string inputfile, std::string type, FitWeight *rw, std::string fkdt) { //******************************************************************** // Check if name contains Evt, indicating that it is a raw number of events // measurements and should thus be treated as once fIsRawEvents = false; if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but fIsRawEvents == false!"); NUIS_LOG(SAM, "Overriding this and setting fIsRawEvents == true!"); } fIsEnu = false; if ((fName.find("XSec") != std::string::npos) && (fName.find("Enu") != std::string::npos)) { fIsEnu = true; NUIS_LOG(SAM, "::" << fName << "::"); NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " "not flux averaged!"); if (FitPar::Config().GetParB("EventManager")) { NUIS_ERR(FTL, "Enu Measurements do not yet work with the Event Manager!"); NUIS_ERR(FTL, "If you want decent flux unfolded results please run in " "series mode (-q EventManager=0)"); sleep(2); throw; } } if (fIsEnu && fIsRawEvents) { NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__); } // Reset everything to NULL fRW = rw; // Setting up 2D Inputs this->SetupInputs(inputfile); // Set Default Options SetFitOptions(fDefaultTypes); // Set Passed Options SetFitOptions(type); } //******************************************************************** void Measurement2D::SetupDefaultHist() { //******************************************************************** // Setup fMCHist fMCHist = (TH2D *)fDataHist->Clone(); fMCHist->SetNameTitle((fName + "_MC").c_str(), (fName + "_MC" + fPlotTitles).c_str()); // Setup fMCFine Int_t nBinsX = fMCHist->GetNbinsX(); Int_t nBinsY = fMCHist->GetNbinsY(); fMCFine = new TH2D((fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(), nBinsX * 3, fMCHist->GetXaxis()->GetBinLowEdge(1), fMCHist->GetXaxis()->GetBinLowEdge(nBinsX + 1), nBinsY * 3, fMCHist->GetYaxis()->GetBinLowEdge(1), fMCHist->GetYaxis()->GetBinLowEdge(nBinsY + 1)); // Setup MC Stat fMCStat = (TH2D *)fMCHist->Clone(); fMCStat->Reset(); // Setup the NEUT Mode Array // PlotUtils::CreateNeutModeArray(fMCHist, (TH1**)fMCHist_PDG); // Setup bin masks using sample name if (fIsMask) { std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask"); if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask"; } SetBinMask(maskloc); } return; } //******************************************************************** void Measurement2D::SetDataValues(std::string dataFile, std::string TH2Dname) { //******************************************************************** if (dataFile.find(".root") == std::string::npos) { NUIS_ERR(FTL, "Error! " << dataFile << " is not a .root file"); NUIS_ERR(FTL, "Currently only .root file reading is supported (MiniBooNE " "CC1pi+ 2D), but implementing .txt should be dirt easy"); NUIS_ABORT("See me at " << __FILE__ << ":" << __LINE__); } else { TFile *inFile = new TFile(dataFile.c_str(), "READ"); fDataHist = (TH2D *)(inFile->Get(TH2Dname.c_str())->Clone()); fDataHist->SetDirectory(0); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_MC" + fPlotTitles).c_str()); delete inFile; } return; } //******************************************************************** void Measurement2D::SetDataValues(std::string dataFile, double dataNorm, std::string errorFile, double errorNorm) { //******************************************************************** // Make a counter to track the line number int yBin = 0; std::string line; std::ifstream data(dataFile.c_str(), std::ifstream::in); fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(), fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins); if (data.is_open()) { NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str()); } while (std::getline(data >> std::ws, line, '\n')) { int xBin = 0; // Loop over entries and insert them into the histogram std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { fDataHist->SetBinContent(xBin + 1, yBin + 1, (*iter) * dataNorm); xBin++; } yBin++; } yBin = 0; std::ifstream error(errorFile.c_str(), std::ifstream::in); if (error.is_open()) { NUIS_LOG(SAM, "Reading errors from: " << errorFile.c_str()); } while (std::getline(error >> std::ws, line, '\n')) { int xBin = 0; // Loop over entries and insert them into the histogram std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { fDataHist->SetBinError(xBin + 1, yBin + 1, (*iter) * errorNorm); xBin++; } yBin++; } return; }; //******************************************************************** void Measurement2D::SetDataValuesFromText(std::string dataFile, double dataNorm) { //******************************************************************** fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(), fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins); NUIS_LOG(SAM, "Reading data from: " << dataFile); PlotUtils::Set2DHistFromText(dataFile, fDataHist, dataNorm, true); return; }; //******************************************************************** void Measurement2D::SetCovarMatrix(std::string covarFile) { //******************************************************************** // Used to read a covariance matrix from a root file TFile *tempFile = new TFile(covarFile.c_str(), "READ"); // Make plots that we want TH2D *covarPlot = new TH2D(); TH2D *fFullCovarPlot = new TH2D(); // Get covariance options for fake data studies std::string covName = ""; std::string covOption = FitPar::Config().GetParS("throw_covariance"); // Which matrix to get? if (fIsShape || fIsFree) covName = "shp_"; if (fIsDiag) covName += "diag"; else covName += "full"; covarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); // Throw either the sub matrix or the full matrix if (!covOption.compare("SUB")) fFullCovarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); else if (!covOption.compare("FULL")) fFullCovarPlot = (TH2D *)tempFile->Get("fullcov"); else { NUIS_ERR(WRN, " Incorrect thrown_covariance option in parameters."); } // Bin masking? int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral()); int covdim = int(fDataHist->GetNbinsX()); // Make new covars this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); fDecomp = new TMatrixDSym(dim); // Full covariance values int row, column = 0; row = 0; column = 0; for (Int_t i = 0; i < covdim; i++) { // masking can be dodgy // if (this->masked->GetBinContent(i+1) > 0) continue; for (Int_t j = 0; j < covdim; j++) { // if (this->masked->GetBinContent(j+1) > 0) continue; (*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1); (*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1); column++; } column = 0; row++; } // Set bin errors on data if (!fIsDiag) { for (Int_t i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinError( i + 1, sqrt((covarPlot->GetBinContent(i + 1, i + 1))) * 1E-38); } } TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); tempFile->Close(); delete tempFile; return; }; //******************************************************************** void Measurement2D::SetCovarMatrixFromText(std::string covarFile, int dim) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream covar(covarFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); if (covar.is_open()) { NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile); } while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation // matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { double val = (*iter) * fDataHist->GetBinError(row + 1) * 1E38 * fDataHist->GetBinError(column + 1) * 1E38; (*this->covar)(row, column) = val; (*fFullCovar)(row, column) = val; column++; } row++; } // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** void Measurement2D::SetCovarMatrixFromChol(std::string covarFile, int dim) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream covarread(covarFile.c_str(), std::ifstream::in); TMatrixD *newcov = new TMatrixD(dim, dim); if (covarread.is_open()) { NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile); } while (std::getline(covarread >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation // matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*newcov)(row, column) = *iter; column++; } row++; } covarread.close(); // Form full covariance TMatrixD *trans = (TMatrixD *)(newcov)->Clone(); trans->T(); (*trans) *= (*newcov); fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), ""); delete newcov; delete trans; // Robust matrix inversion method TDecompChol LU = TDecompChol(*this->fFullCovar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; // //******************************************************************** // void Measurement2D::SetMapValuesFromText(std::string dataFile) { // //******************************************************************** // fMapHist = new TH2I((fName + "_map").c_str(), (fName + // fPlotTitles).c_str(), // fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins); // LOG(SAM) << "Reading map from: " << dataFile << std::endl; // PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0); // return; // }; diff --git a/src/FitBase/Measurement2D.h b/src/FitBase/Measurement2D.h index 2ab738c..00d32f8 100644 --- a/src/FitBase/Measurement2D.h +++ b/src/FitBase/Measurement2D.h @@ -1,641 +1,644 @@ // 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 MEASUREMENT_2D_HXX_SEEN #define MEASUREMENT_2D_HXX_SEEN /*! * \addtogroup FitBase * @{ */ #include #include #include #include #include #include #include #include #include // ROOT includes #include #include #include #include #include #include #include #include #include // External data fit includes #include "FitEvent.h" #include "FitUtils.h" #include "MeasurementBase.h" #include "PlotUtils.h" #include "SignalDef.h" #include "StatUtils.h" #include "MeasurementVariableBox2D.h" //******************************************************************** //! 2D Measurement base class. Histogram handling is done in this base layer. class Measurement2D : public MeasurementBase { //******************************************************************** public: /* Constructor/Deconstuctor */ //! Default Constructor Measurement2D(); //! Default Destructor virtual ~Measurement2D(); /* Setup Functions */ /// \brief Setup all configs once initialised /// /// Should be called after all configs have been setup inside fSettings container. /// Handles the processing of inputs and setting up of types. /// Replaces the old 'SetupMeasurement' function. void FinaliseSampleSettings(); /// \brief Creates the 2D data distribution given the binning provided. virtual void CreateDataHistogram(int dimx, double* binx, int dimy, double* biny); /// \brief Set Data Histogram from a list of contents in a text file /// /// Assumes the format: \n /// x_low_1 y_low_1 cont_11 err_11 \n /// x_low_1 y_low_2 cont_12 err_12 \n /// x_low_2 y_low_1 cont_21 err_21 \n /// x_low_2 y_low_2 cont_22 err_22 \n /// x_low_2 y_low_3 cont_23 err_23 \n /// x_low_3 y_low_2 cont_32 err_32 \n virtual void SetDataFromTextFile(std::string data, std::string binx, std::string biny); /// \brief Set Data Histogram from a TH2D in a file /// /// - datfile = Full path to data file /// - histname = Name of histogram /// /// If histname not given it assumes that datfile /// is in the format: \n /// 'file.root;histname' virtual void SetDataFromRootFile(std::string datfile, std::string histname=""); /// \brief Set data values from a 2D array in text file /// /// \warning requires DATA HISTOGRAM TO BE SET FIRST /// /// Assumes form: \n /// cont_11 cont_12 ... cont_1N \n /// cont_21 cont_22 ... cont_2N \n /// ... ... ... ... \n /// cont_N1 cont_N2 ... cont_NN \n virtual void SetDataValuesFromTextFile(std::string datfile, TH2D* hist = NULL); /// \brief Set data errors from a 2D array in text file /// /// \warning requires DATA HISTOGRAM TO BE SET FIRST /// /// Assumes form: \n /// errs_11 errs_12 ... errs_1N \n /// errs_21 errs_22 ... errs_2N \n /// ... ... ... ... \n /// errs_N1 errs_N2 ... errs_NN \n virtual void SetDataErrorsFromTextFile(std::string datfile, TH2D* hist = NULL); /// \brief Set data bin errors to sqrt(entries) /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Sets the data errors as the sqrt of the bin contents /// Should be use for counting experiments virtual void SetPoissonErrors(); /// \brief Make diagonal covariance from data /// /// \warning If no histogram passed, data must be setup first! /// Setup the covariance inputs by taking the data histogram /// errors and setting up a diagonal covariance matrix. /// /// If no data is supplied, fDataHist is used if already set. virtual void SetCovarFromDiagonal(TH2D* data = NULL); /// \brief Read the data covariance from a text file. /// /// Inputfile should have the format: \n /// covariance_11 covariance_12 covariance_13 ... \n /// covariance_21 covariance_22 covariance_23 ... \n /// ... ... ... ... \n /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarFromTextFile(std::string covfile, int dim = -1); /// \brief Read the data covariance from a ROOT file. /// /// - covfile specifies the full path to the file /// - histname specifies the name of the covariance object. Both TMatrixDSym and TH2D are supported. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the inverted data covariance from a text file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarInvertFromTextFile(std::string covfile, int dim = -1); /// \brief Read the inverted data covariance from a ROOT file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarInvertFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the data correlations from a text file. /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCorrelationFromTextFile(std::string covfile, int dim = -1); /// \brief Read the data correlations from a ROOT file. /// /// \warning REQUIRES DATA TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCorrelationFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the cholesky decomposed covariance from a text file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCholDecompFromTextFile(std::string covfile, int dim = -1); /// \brief Read the cholesky decomposed covariance from a ROOT file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCholDecompFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the map values from a text file /// /// \warning Requires DATA hist to be set beforehand. - /// Format should be a 2D array of mappings. + /// Format should be a 2D array of mappings. /// -1 indicates empty bins. \n /// e.g.: \n /// 1 2 3 4 5 \n /// -1 6 7 8 9 \n /// -1 -1 10 11 -1 \n virtual void SetMapValuesFromText(std::string dataFile); /// \brief Scale the data by some scale factor virtual void ScaleData(double scale); /// \brief Scale the data error bars by some scale factor virtual void ScaleDataErrors(double scale); /// \brief Scale the covariaince and its invert/decomp by some scale factor. virtual void ScaleCovar(double scale); /// \brief Setup a bin masking histogram and apply masking to data /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Reads in a list of bins in a text file to be masked. Format is: \n /// bin_index_x_1 bin_index_y_1 1 \n /// bin_index_x_2 bin_index_y_2 1 \n /// bin_index_x_3 bin_index_y_3 1 \n /// /// If 0 is given then a bin entry will NOT be masked. So for example: \n\n /// 1 1 1 \n /// 2 0 1 \n /// 3 4 0 \n /// 4 0 1 \n\n /// Will mask only the (1,1), (2,0), and (4,0) bins. /// /// Masking can be turned on by specifiying the MASK option when creating a sample. /// When this is passed NUISANCE will look in the following locations for the mask file: /// - FitPar::Config().GetParS(fName + ".mask") /// - "data/masks/" + fName + ".mask"; virtual void SetBinMask(std::string maskfile); /// \brief Set the current fit options from a string. /// /// This is called twice for each sample, once to set the default /// and once to set the current setting (if anything other than default given) /// /// For this to work properly it requires the default and allowed types to be /// set correctly. These should be specified as a string listing options. /// /// To split up options so that NUISANCE can automatically detect ones that /// are conflicting. Any options seperated with the '/' symbol are non conflicting /// and can be given together, whereas any seperated with the ',' symbol cannot /// be specified by the end user at the same time. /// /// Default Type Examples: /// - DIAG/FIX = Default option will be a diagonal covariance, with FIXED norm. /// - MASK/SHAPE = Default option will be a masked hist, with SHAPE always on. /// /// Allowed Type examples: /// - 'FULL/DIAG/NORM/MASK' = Any of these options can be specified. /// - 'FULL,FREE,SHAPE/MASK/NORM' = User can give either FULL, FREE, or SHAPE as on option. /// MASK and NORM can also be included as options. virtual void SetFitOptions(std::string opt); /// \brief Final constructor setup /// \warning Should be called right at the end of the constructor. /// /// Contains a series of checks to ensure the data and inputs have been setup. /// Also creates the MC histograms needed for fitting. void FinaliseMeasurement(); /* Reconfigure */ /// \brief Create a Measurement1D box /// /// Creates a new 1D variable box containing just fXVar. /// /// This box is the bare minimum required by the JointFCN when /// running fast reconfigures during a routine. /// If for some reason a sample needs extra variables to be saved then /// it should override this function creating its own MeasurementVariableBox /// that contains the extra variables. virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox2D();}; /// \brief Reset all MC histograms /// /// Resets all standard histograms and those registered to auto /// process to zero. /// /// If extra histograms are not included in auto processing, then they must be reset /// by overriding this function and doing it manually if required. virtual void ResetAll(void); /// \brief Fill MC Histograms from XVar, YVar /// /// Fill standard histograms using fXVar, fYVar, Weight read from the variable box. /// /// WARNING : Any extra MC histograms need to be filled by overriding this function, /// even if they have been set to auto process. virtual void FillHistograms(void); // \brief Convert event rates to final histogram /// /// Apply standard scaling procedure to standard mc histograms to convert from /// raw events to xsec prediction. /// /// If any distributions have been set to auto process /// that is done during this function call, and a differential xsec is assumed. /// If that is not the case this function must be overriden. virtual void ScaleEvents(void); /// \brief Scale MC by a factor=1/norm /// /// Apply a simple normalisation scaling if the option FREE or a norm_parameter /// has been specified in the NUISANCE routine. virtual void ApplyNormScale(double norm); /* Statistical Functions */ /// \brief Get Number of degrees of freedom /// /// Returns the number bins inside the data histogram accounting for /// any bin masking applied. virtual int GetNDOF(); /// \brief Return Data/MC Likelihood at current state /// /// Returns the likelihood of the data given the current MC prediction. /// Diferent likelihoods definitions are used depending on the FitOptions. virtual double GetLikelihood(void); /* Fake Data */ /// \brief Set the fake data values from either a file, or MC /// /// - Setting from a file "path": \n /// When reading from a file the full path must be given to a standard /// nuisance output. The standard MC histogram should have a name that matches /// this sample for it to be read in. /// \n\n /// - Setting from "MC": \n /// If the MC option is given the current MC prediction is used as fake data. virtual void SetFakeDataValues(std::string fakeOption); /// \brief Reset fake data back to starting fake data /// /// Reset the fake data back to original fake data (Reset back to before /// ThrowCovariance was first called) virtual void ResetFakeData(void); /// \brief Reset fake data back to original data /// /// Reset the data histogram back to the true original dataset for this sample /// before any fake data was defined. virtual void ResetData(void); /// \brief Generate fake data by throwing the covariance. /// /// Can be used on fake MC data or just the original dataset. /// Call ResetFakeData or ResetData to return to values before the throw. virtual void ThrowCovariance(void); - /// \brief Throw the data by its assigned errors and assign this to MC - /// - /// Used when creating data toys by assign the MC to this thrown data - /// so that the likelihood is calculated between data and thrown data + /// \brief Throw the data by its assigned errors and assign this to MC + /// + /// Used when creating data toys by assign the MC to this thrown data + /// so that the likelihood is calculated between data and thrown data virtual void ThrowDataToy(void); /* Access Functions */ /// \brief Returns nicely formatted MC Histogram /// /// Format options can also be given in the samplesettings: /// - linecolor /// - linestyle /// - linewidth /// - fillcolor /// - fillstyle /// /// So to have a sample line colored differently in the xml cardfile put: \n /// virtual TH2D* GetMCHistogram(void); /// \brief Returns nicely formatted data Histogram /// /// Format options can also be given in the samplesettings: /// - datacolor /// - datastyle /// - datawidth /// /// So to have a sample data colored differently in the xml cardfile put: \n /// virtual TH2D* GetDataHistogram(void); /// \brief Returns a list of all MC histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMCList(void) { return std::vector(1, GetMCHistogram()); } /// \brief Returns a list of all Data histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetDataList(void) { return std::vector(1, GetDataHistogram()); } /// \brief Returns a list of all Mask histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMaskList(void) { return std::vector(1, fMaskHist); }; /// \brief Returns a list of all Fine histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetFineList(void) { return std::vector(1, fMCFine); }; /* Write Functions */ /// \brief Save the current state to the current TFile directory \n /// /// Data/MC are both saved by default. /// A range of other histograms can be saved by setting the /// config option 'drawopts'. /// /// Possible options: \n /// - FINE = Write Fine Histogram \n /// - WEIGHTS = Write Weighted MC Histogram (before scaling) \n /// - FLUX = Write Flux Histogram from MC Input \n /// - EVT = Write Event Histogram from MC Input \n /// - XSEC = Write XSec Histogram from MC Input \n /// - MASK = Write Mask Histogram \n /// - COV = Write Covariance Histogram \n /// - INVCOV = Write Inverted Covariance Histogram \n /// - DECMOP = Write Decomp. Covariance Histogram \n /// - RESIDUAL= Write Resudial Histograms \n /// - RATIO = Write Data/MC Ratio Histograms \n /// - SHAPE = Write MC Shape Histograms norm. to Data \n /// - CANVMC = Build MC Canvas Showing Data, MC, Shape \n /// - MODES = Write PDG Stack \n /// - CANVPDG = Build MC Canvas Showing Data, PDGStack \n /// /// So to save a range of these in parameters/config.xml set: \n /// virtual void Write(std::string drawOpt); //////// OLD FUNCTIONS //////////// //! Intial setup of common measurement variables. Parse input files, types, //! etc. virtual void SetupMeasurement(std::string input, std::string type, FitWeight* rw, std::string fkdt); //! Setup the default mc Hist given a data histogram virtual void SetupDefaultHist(); //! Set the data values and errors from two files virtual void SetDataValues(std::string dataFile, double dataNorm, std::string errorFile, double errorNorm); virtual void SetDataValues(std::string dataFile, std::string TH2Dname); //! Set the data values only from a text file virtual void SetDataValuesFromText(std::string dataFile, double norm); //! Read a covariance matrix from a file (Default name "covar" in file) virtual void SetCovarMatrix(std::string covarFile); //! Set the covariance matrix from a text file virtual void SetCovarMatrixFromText(std::string covarFile, int dim); //! Set the covariance matrix from a text file containing the cholesky //! fDecomposition virtual void SetCovarMatrixFromChol(std::string covarFile, int dim); protected: // The data histograms TH2D* fDataHist; //!< default data histogram (use in chi2 calculations) TH2D* fDataOrig; //!< histogram to store original data before throws. TH2D* fDataTrue; //!< histogram to store true dataset TH1D* fDataHist_X; //!< Projections onto X of the fDataHist TH1D* fDataHist_Y; //!< Projections onto Y of the fDataHist // The MC histograms TH2D* fMCHist; //!< MC Histogram (used in chi2 calculations) TH2D* fMCFine; //!< Finely binned MC Histogram TH2D* fMCHist_PDG[61]; //!< MC Histograms for each interaction mode TH1D* fMCHist_X; //!< Projections onto X of the fMCHist TH1D* fMCHist_Y; //!< Projections onto Y of the fMCHist TH2D* fMCWeighted; //!< Raw Event Weights TH2D* fMCStat; TH2I* fMaskHist; //!< mask histogram for the data TH2I* fMapHist; //!< map histogram used to convert 2D to 1D distributions + TH2D *fResidualHist; + TH2D *fChi2LessBinHist; + bool fIsFakeData; //!< is current data actually fake std::string fakeDataFile; //!< MC fake data input file std::string fPlotTitles; //!< X and Y plot titles. std::string fFitType; std::string fDefaultTypes; //!< Default Fit Options std::string fAllowedTypes; //!< Any allowed Fit Options TMatrixDSym* covar; //!< inverted covariance matrix TMatrixDSym* fFullCovar; //!< covariance matrix TMatrixDSym* fDecomp; //!< fDecomposed covariance matrix TMatrixDSym* fCorrel; //!< correlation matrix double fCovDet; //!< covariance deteriminant double fNormError; //!< Normalisation on the error on the data double fLikelihood; //!< Likelihood value Double_t* fXBins; //!< X Bin Edges Double_t* fYBins; //!< Y Bin Edges Int_t fNDataPointsX; //!< Number of X data points Int_t fNDataPointsY; //!< NUmber of Y data points // Fit specific flags bool fIsShape; //!< Flag: Perform shape-only fit bool fIsFree; //!< Flag: Perform normalisation free fit bool fIsDiag; //!< Flag: Only use diagonal bin errors in stats bool fIsMask; //!< Flag: Apply bin masking bool fIsRawEvents; //!< Flag: Only event rates in histograms bool fIsEnu; //!< Needs Enu Unfolding bool fIsChi2SVD; //!< Flag: Chi2 SVD Method (DO NOT USE) bool fAddNormPen; //!< Flag: Add normalisation penalty to fi bool fIsProjFitX; //!< Flag: Use 1D projections onto X and Y to calculate the //!Chi2 Method. If flagged X will be used to set the rate. bool fIsProjFitY; //!< Flag: Use 1D projections onto X and Y to calculate the //!Chi2 Method. If flagged Y will be used to set the rate. bool fIsFix; //!< Flag: Fixed Histogram Norm bool fIsFull; //!< Flag; Use Full Covar bool fIsDifXSec; //!< Flag: Differential XSec bool fIsEnu1D; //!< Flag: Flux Unfolded XSec bool fIsChi2; //!< Flag; Use Chi2 over LL TrueModeStack* fMCHist_Modes; ///< Optional True Mode Stack TMatrixDSym* fCovar; ///< New FullCovar TMatrixDSym* fInvert; ///< New covar // Fake Data std::string fFakeDataInput; ///< Input fake data file path TFile* fFakeDataFile; ///< Input fake data file // Arrays for data entries Double_t* fDataValues; ///< REMOVE data bin contents Double_t* fDataErrors; ///< REMOVE data bin errors }; /*! @} */ #endif diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx index 7990f68..fe5eb18 100644 --- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx +++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx @@ -1,245 +1,284 @@ // 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 . *******************************************************************************/ /* Authors: Adrian Orea (v1 2017) Clarence Wret (v2 2018) */ #include "MINERvA_CC0pi_XSec_2D_nu.h" #include "MINERvA_SignalDef.h" //******************************************************************** void MINERvA_CC0pi_XSec_2D_nu::SetupDataSettings() { //******************************************************************** // Set Distribution // See header file for enum and some descriptions std::string name = fSettings.GetS("name"); // Have one distribution as of summer 2018 if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) fDist = kPtPz; // Define what files to use from the dist std::string datafile = ""; std::string corrfile = ""; std::string titles = ""; std::string distdescript = ""; std::string histname = ""; switch (fDist) { case (kPtPz): datafile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root"; corrfile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root"; titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} " "(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)"; distdescript = "MINERvA_CC0pi_XSec_2Dptpz_nu sample"; histname = "pt_pl_cross_section"; break; default: NUIS_ABORT("Unknown Analysis Distribution : " << name); } fSettings.SetTitle(GeneralUtils::ParseToStr(titles, ";")[0]); fSettings.SetXTitle(GeneralUtils::ParseToStr(titles, ";")[1]); fSettings.SetYTitle(GeneralUtils::ParseToStr(titles, ";")[2]); fSettings.SetZTitle(GeneralUtils::ParseToStr(titles, ";")[3]); // Sample overview --------------------------------------------------- std::string descrip = distdescript + "\n" "Target: CH \n" "Flux: MINERvA Low Energy FHC numu \n" "Signal: CC-0pi \n"; fSettings.SetDescription(descrip); // The input ROOT file fSettings.SetDataInput(FitPar::GetDataBase() + datafile); fSettings.SetCovarInput(FitPar::GetDataBase() + corrfile); // Set the data file SetDataValues(fSettings.GetDataInput(), histname); } //******************************************************************** MINERvA_CC0pi_XSec_2D_nu::MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey) { //******************************************************************** // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); SetupDataSettings(); FinaliseSampleSettings(); fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); TMatrixDSym *tempmat = StatUtils::GetCovarFromRootFile( fSettings.GetCovarInput(), "TotalCovariance"); fFullCovar = tempmat; // Decomposition is stable for entries that aren't E-xxx double ScalingFactor = 1E38 * 1E38; (*fFullCovar) *= ScalingFactor; // Just check that the data error and covariance are the same for (int i = 0; i < fFullCovar->GetNrows(); ++i) { for (int j = 0; j < fFullCovar->GetNcols(); ++j) { // Get the global bin int xbin1, ybin1, zbin1; fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1); double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1); double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1 + 1); double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1); double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1 + 1); if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) || ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) || xhi1 > fDataHist->GetXaxis()->GetBinLowEdge( fDataHist->GetXaxis()->GetNbins() + 1) || yhi1 > fDataHist->GetYaxis()->GetBinLowEdge( fDataHist->GetYaxis()->GetNbins() + 1)) continue; double data_error = fDataHist->GetBinError(xbin1, ybin1); double cov_error = sqrt((*fFullCovar)(i, i) / ScalingFactor); if (fabs(data_error - cov_error) > 1E-5) { std::cerr << "Error on data is different to that of covariance" << std::endl; NUIS_ERR(FTL, "Data error: " << data_error); NUIS_ERR(FTL, "Cov error: " << cov_error); NUIS_ERR(FTL, "Data/Cov: " << data_error / cov_error); NUIS_ERR(FTL, "Data-Cov: " << data_error - cov_error); NUIS_ERR(FTL, "For x: " << xlo1 << "-" << xhi1); NUIS_ABORT("For y: " << ylo1 << "-" << yhi1); } } } // Now can make the inverted covariance covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); // Now scale back (*fFullCovar) *= 1.0 / ScalingFactor; (*covar) *= ScalingFactor; (*fDecomp) *= ScalingFactor; // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CC0pi_XSec_2D_nu::FillEventVariables(FitEvent *event) { //******************************************************************** // Checking to see if there is a Muon if (event->NumFSParticle(13) == 0) return; // Get the muon kinematics TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; Double_t px = Pmu.X() / 1000; Double_t py = Pmu.Y() / 1000; Double_t pt = sqrt(px * px + py * py); // y-axis is transverse momentum for both measurements fYVar = pt; // Now we set the x-axis switch (fDist) { case kPtPz: { // Don't want to assume the event generators all have neutrino coming along // z pz is muon momentum projected onto the neutrino direction Double_t pz = Pmu.Vect().Dot(Pnu.Vect() * (1.0 / Pnu.Vect().Mag())) / 1000.; // Set Hist Variables fXVar = pz; break; } default: NUIS_ABORT("DIST NOT FOUND : " << fDist); break; } }; //******************************************************************** bool MINERvA_CC0pi_XSec_2D_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax); }; //******************************************************************** // Custom likelihood calculator because binning of covariance matrix double MINERvA_CC0pi_XSec_2D_nu::GetLikelihood() { //******************************************************************** // The calculated chi2 double chi2 = 0.0; // Support shape comparisons double scaleF = fDataHist->Integral() / fMCHist->Integral(); if (fIsShape) { fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); } // Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA // used for their measurement Can be prettified in due time but for now keep int nbinsx = fMCHist->GetNbinsX(); int nbinsy = fMCHist->GetNbinsY(); Int_t Nbins = nbinsx * nbinsy; // Loop over the covariance matrix bins for (int i = 0; i < Nbins; ++i) { - int xbin = i % nbinsx + 1; - int ybin = i / nbinsx + 1; + int xbin = (i % nbinsx) + 1; + int ybin = (i / nbinsx) + 1; double datax = fDataHist->GetBinContent(xbin, ybin); double mcx = fMCHist->GetBinContent(xbin, ybin); + double chi2_bin = 0; for (int j = 0; j < Nbins; ++j) { - int xbin2 = j % nbinsx + 1; - int ybin2 = j / nbinsx + 1; + int xbin2 = (j % nbinsx) + 1; + int ybin2 = (j / nbinsx) + 1; double datay = fDataHist->GetBinContent(xbin2, ybin2); double mcy = fMCHist->GetBinContent(xbin2, ybin2); double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); - chi2 += chi2_xy; + chi2_bin += chi2_xy; + } + if (fResidualHist) { + fResidualHist->SetBinContent(xbin, ybin, chi2_bin); + } + chi2 += chi2_bin; + } + + if (fChi2LessBinHist) { + for (int igbin = 0; igbin < Nbins; ++igbin) { + int igxbin = (igbin % nbinsx) + 1; + int igybin = (igbin / nbinsx) + 1; + double tchi2 = 0; + for (int i = 0; i < Nbins; ++i) { + int xbin = (i % nbinsx) + 1; + int ybin = (i / nbinsx) + 1; + if ((xbin == igxbin) && (ybin == igybin)) { + continue; + } + double datax = fDataHist->GetBinContent(xbin, ybin); + double mcx = fMCHist->GetBinContent(xbin, ybin); + double chi2_bin = 0; + for (int j = 0; j < Nbins; ++j) { + int xbin2 = (j % nbinsx) + 1; + int ybin2 = (j / nbinsx) + 1; + if ((xbin2 == igxbin) && (ybin2 == igybin)) { + continue; + } + double datay = fDataHist->GetBinContent(xbin2, ybin2); + double mcy = fMCHist->GetBinContent(xbin2, ybin2); + + double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); + chi2_bin += chi2_xy; + } + tchi2 += chi2_bin; + } + + fChi2LessBinHist->SetBinContent(igxbin, igybin, tchi2); } } // Normalisation penalty term if included if (fAddNormPen) { chi2 += (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError); - NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / - (fNormError * fNormError)); + NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * + (1 - (fCurrentNorm)) / + (fNormError * fNormError)); } // Adjust the shape back to where it was. if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) { fMCHist->Scale(1. / scaleF); fMCFine->Scale(1. / scaleF); } fLikelihood = chi2; return chi2; }; diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx index 97c060e..827a00f 100644 --- a/src/Reweight/GENIEWeightEngine.cxx +++ b/src/Reweight/GENIEWeightEngine.cxx @@ -1,478 +1,478 @@ #include "GENIEWeightEngine.h" #ifdef __GENIE_ENABLED__ #ifdef GENIE_PRE_R3 #include "Algorithm/AlgConfigPool.h" #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" #include "Ntuple/NtpMCEventRecord.h" -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ #include "ReWeight/GReWeightAGKY.h" #include "ReWeight/GReWeightDISNuclMod.h" #include "ReWeight/GReWeightFGM.h" #include "ReWeight/GReWeightFZone.h" #include "ReWeight/GReWeightINuke.h" #include "ReWeight/GReWeightNonResonanceBkg.h" #include "ReWeight/GReWeightNuXSecCCQE.h" #include "ReWeight/GReWeightNuXSecCCQEvec.h" #include "ReWeight/GReWeightNuXSecCCRES.h" #include "ReWeight/GReWeightNuXSecCOH.h" #include "ReWeight/GReWeightNuXSecDIS.h" #include "ReWeight/GReWeightNuXSecNC.h" #include "ReWeight/GReWeightNuXSecNCEL.h" #include "ReWeight/GReWeightNuXSecNCRES.h" #include "ReWeight/GReWeightResonanceDecay.h" #include "ReWeight/GSystUncertainty.h" #ifdef __GENIE_EMP_MECRW_ENABLED #include "ReWeight/GReWeightXSecEmpiricalMEC.h" #endif #endif #if __GENIE_VERSION__ >= 212 #include "ReWeight/GReWeightNuXSecCCQEaxial.h" #endif #else // GENIE v3 #include "Framework/Algorithm/AlgConfigPool.h" #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/Ntuple/NtpMCEventRecord.h" #include "Framework/Utils/AppInit.h" #include "Framework/Utils/RunOpt.h" using namespace genie; -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ #include "RwCalculators/GReWeightAGKY.h" #include "RwCalculators/GReWeightDISNuclMod.h" #include "RwCalculators/GReWeightFGM.h" #include "RwCalculators/GReWeightFZone.h" #include "RwCalculators/GReWeightINuke.h" #include "RwCalculators/GReWeightNonResonanceBkg.h" #include "RwCalculators/GReWeightNuXSecCCQE.h" #include "RwCalculators/GReWeightNuXSecCCQEaxial.h" #include "RwCalculators/GReWeightNuXSecCCQEvec.h" #include "RwCalculators/GReWeightNuXSecCCRES.h" #include "RwCalculators/GReWeightNuXSecCOH.h" #include "RwCalculators/GReWeightNuXSecDIS.h" #include "RwCalculators/GReWeightNuXSecNC.h" #include "RwCalculators/GReWeightNuXSecNCEL.h" #include "RwCalculators/GReWeightNuXSecNCRES.h" #include "RwCalculators/GReWeightResonanceDecay.h" #include "RwCalculators/GReWeightXSecEmpiricalMEC.h" #include "RwFramework/GSystUncertainty.h" using namespace genie::rew; #endif #endif #endif #include "FitLogger.h" GENIEWeightEngine::GENIEWeightEngine(std::string name) { #ifdef __GENIE_ENABLED__ -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ // Setup the NEUT Reweight engien fCalcName = name; NUIS_LOG(DEB, "Setting up GENIE RW : " << fCalcName); // Create RW Engine suppressing cout StopTalking(); fGenieRW = new genie::rew::GReWeight(); // Get List of Vetos (Just for debugging) std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fGenieRW_veto"); bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos; bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos; bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos; bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos; bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos; bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos; bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos; bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos; #ifdef __GENIE_EMP_MECRW_ENABLED bool xsec_empMEC = rw_engine_list.find("xsec_empMEC") == std::string::npos; #endif #ifndef GENIE_PRE_R3 genie::RunOpt *grunopt = genie::RunOpt::Instance(); grunopt->EnableBareXSecPreCalc(true); grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList")); if(!Config::HasPar("GENIETune")){ NUIS_ABORT("GENIE tune was not specified, this is required when reweighting GENIE V3+ events. Add a config parameter like: to your nuisance card."); } grunopt->SetTuneName(Config::GetParS("GENIETune")); grunopt->BuildTune(); std::string genv = std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml"; genie::utils::app_init::MesgThresholds(genv); #endif // Now actually add the RW Calcs if (xsec_ncel) fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL); if (xsec_ccqe) { fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE); // (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") ); } #ifdef __GENIE_EMP_MECRW_ENABLED if (xsec_empMEC) { fGenieRW->AdoptWghtCalc("xsec_empMEC", new genie::rew::GReWeightXSecEmpiricalMEC); } #endif if (xsec_coh) { fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH()); // (dynamic_cast (fGenieRW->WghtCalc("xsec_coh"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") ); } if (xsec_nnres) fGenieRW->AdoptWghtCalc("xsec_nonresbkg", new genie::rew::GReWeightNonResonanceBkg); if (xsec_nudis) fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod); if (xsec_resdec) fGenieRW->AdoptWghtCalc("hadro_res_decay", new genie::rew::GReWeightResonanceDecay); if (xsec_fzone) fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); if (xsec_intra) fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); if (xsec_agky) fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY); if (xsec_qevec) fGenieRW->AdoptWghtCalc("xsec_ccqe_vec", new genie::rew::GReWeightNuXSecCCQEvec); #if __GENIE_VERSION__ >= 212 if (xsec_qeaxial) fGenieRW->AdoptWghtCalc("xsec_ccqe_axial", new genie::rew::GReWeightNuXSecCCQEaxial); #endif if (xsec_dis) fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS); if (xsec_nc) fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC); if (xsec_ccres) { #if __GENIE_VERSION__ < 213 fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES); #else fGenieRW->AdoptWghtCalc( "xsec_ccres", new genie::rew::GReWeightNuXSecCCRES( FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default")); #endif } if (xsec_ncres) fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES); if (xsec_nucqe) fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM); // Set the CCQE reweighting style GReWeightNuXSecCCQE *rwccqe = dynamic_cast(fGenieRW->WghtCalc("xsec_ccqe")); // For MaCCQE reweighting std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode"); if (ccqetype == "kModeMa") { NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeMa"); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa); } else if (ccqetype == "kModeNormAndMaShape") { NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeNormAndMaShape"); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeNormAndMaShape); // For z-expansion reweighting } else if (ccqetype == "kModeZExp") { NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeZExp"); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeZExp); } else { NUIS_ERR(FTL, "Did not find specified GENIE ReWeight CCQE mode"); NUIS_ABORT("You provided: " << ccqetype << " in parameters/config.xml"); } #if (__GENIE_VERSION__ >= 212) and \ (__GENIE_VERSION__ <= \ 300) // This doesn't currently work as is for GENIE v3, but the reweighting // in v3 supposedly does similar checks anyway. // Check the UserPhysicsOptions too! AlgConfigPool *Pool = genie::AlgConfigPool::Instance(); Registry *full = Pool->GlobalParameterList(); std::string name_ax = full->GetAlg("AxialFormFactorModel").name; std::string config_ax = full->GetAlg("AxialFormFactorModel").config; if (name_ax == "genie::DipoleAxialFormFactorModel" && ccqetype == "kModeZExp") { NUIS_ERR( FTL, "Trying to run Z Expansion reweighting with Llewelyn-Smith model."); NUIS_ERR(FTL, "Please check your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } if (name_ax == "genie::ZExpAxialFormFactorModel" && ccqetype != "kModeZExp") { NUIS_ERR( FTL, "Trying to run Llewelyn-Smith reweighting with Z Expansion model."); NUIS_ERR(FTL, "Please change your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } std::string name_qelcc = full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").name; std::string config_qelcc = full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").config; if (config_qelcc == "Default" && ccqetype == "kModeZExp") { NUIS_ERR( FTL, "Trying to run Z Expansion reweighting with Llewelyn-Smith model."); NUIS_ERR(FTL, "Please change your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } if (config_qelcc == "ZExp" && ccqetype != "kModeZExp") { NUIS_ERR( FTL, "Trying to run Llewelyn-Smith reweighting with Z Expansion model."); NUIS_ERR(FTL, "Please change your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } #endif if (xsec_ccres) { // Default to include shape and normalization changes for CCRES (can be // changed downstream if desired) GReWeightNuXSecCCRES *rwccres = dynamic_cast(fGenieRW->WghtCalc("xsec_ccres")); std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode"); if (!marestype.compare("kModeNormAndMaMvShape")) { rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); } else if (!marestype.compare("kModeMaMv")) { rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); } else { NUIS_ABORT("Unkown MARES Mode in GENIE Weight Engine : " << marestype); } } if (xsec_ncres) { // Default to include shape and normalization changes for NCRES (can be // changed downstream if desired) GReWeightNuXSecNCRES *rwncres = dynamic_cast(fGenieRW->WghtCalc("xsec_ncres")); rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv); } if (xsec_dis) { // Default to include shape and normalization changes for DIS (can be // changed downstream if desired) GReWeightNuXSecDIS *rwdis = dynamic_cast(fGenieRW->WghtCalc("xsec_dis")); rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); } // allow cout again StartTalking(); #else NUIS_ERR(FTL, "GENIE ReWeight is __NOT ENABLED__ in GENIE and you're trying to " "run NUISANCE with it enabled"); NUIS_ERR(FTL, "Check your genie-config --libs for reweighting"); NUIS_ERR(FTL, "If not present you need to recompile GENIE"); NUIS_ABORT("If present you need to contact NUISANCE authors"); #endif #endif }; void GENIEWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __GENIE_ENABLED__ -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ // Get First enum int nuisenum = Reweight::ConvDial(name, kGENIE); // Check ZExp sillyness in GENIE // If ZExpansion parameters are used we need to set a different mode in GENIE // ReWeight... GENIE doesn't have a setter either... #if (__GENIE_VERSION__ >= 212) and (__GENIE_VERSION__ <= 300) std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode"); if (ccqetype != "kModeZExp" && (name == "ZExpA1CCQE" || name == "ZExpA2CCQE" || name == "ZExpA3CCQE" || name == "ZExpA4CCQE")) { NUIS_ERR(FTL, "Found a Z-expansion parameter in GENIE although the GENIE " "ReWeighting engine is set to use Llewelyn-Smith and MaQE!"); NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match requirements"); } if ((ccqetype != "kModeMa" && ccqetype != "kModeMaNormAndMaShape") && (name == "MaCCQE")) { NUIS_ERR(FTL, "Found MaCCQE parameter in GENIE although the GENIE " "ReWeighting engine is set to not use this!"); NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match requirements"); } #endif // Setup Maps fEnumIndex[nuisenum]; // = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fGENIESysts.push_back(rwsyst); // Initialize dial NUIS_LOG(DEB, "Registering " << singlename << " from " << name); fGenieRW->Systematics().Init(fGENIESysts[index]); // If Absolute if (fIsAbsTwk) { GSystUncertainty::Instance()->SetUncertainty(rwsyst, 1.0, 1.0); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } #endif #endif } void GENIEWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __GENIE_ENABLED__ -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif #endif } void GENIEWeightEngine::SetDialValue(std::string name, double val) { #ifdef __GENIE_ENABLED__ -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif #endif } void GENIEWeightEngine::Reconfigure(bool silent) { #ifdef __GENIE_ENABLED__ -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ // Hush now... if (silent) StopTalking(); // Reconf fGenieRW->Reconfigure(); fGenieRW->Print(); // Shout again if (silent) StartTalking(); #endif #endif } double GENIEWeightEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; #ifdef __GENIE_ENABLED__ -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ // Make nom weight if (!evt) { NUIS_ABORT("evt not found : " << evt); } // Skip Non GENIE if (evt->fType != kGENIE) return 1.0; if (!(evt->genie_event)) { NUIS_ABORT("evt->genie_event not found!" << evt->genie_event); } if (!(evt->genie_event->event)) { NUIS_ABORT("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event)); } if (!fGenieRW) { NUIS_ABORT("GENIE RW Not Found!" << fGenieRW); } rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event)); // std::cout << "Returning GENIE Weight for electron scattering = " << // rw_weight << std::endl; // if (rw_weight != 1.0 )std::cout << "mode=" << evt->Mode << " rw_weight = " // << rw_weight << std::endl; #endif #endif // Return rw_weight return rw_weight; } diff --git a/src/Reweight/GENIEWeightEngine.h b/src/Reweight/GENIEWeightEngine.h index 4c39f68..b2e6dc5 100644 --- a/src/Reweight/GENIEWeightEngine.h +++ b/src/Reweight/GENIEWeightEngine.h @@ -1,46 +1,46 @@ #ifndef WEIGHT_ENGINE_GENIE_H #define WEIGHT_ENGINE_GENIE_H #ifdef __GENIE_ENABLED__ #ifdef GENIE_PRE_R3 -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ #include "ReWeight/GSyst.h" #include "ReWeight/GReWeight.h" #endif #else -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ #include "RwFramework/GSyst.h" #include "RwFramework/GReWeight.h" using namespace genie; using namespace genie::rew; #endif #endif #endif #include "GeneratorUtils.h" #include "WeightEngineBase.h" #include "FitWeight.h" class GENIEWeightEngine : public WeightEngineBase { public: GENIEWeightEngine(std::string name); ~GENIEWeightEngine() {}; void IncludeDial(std::string name, double startval); void SetDialValue(int rwenum, double val); void SetDialValue(std::string name, double val); void Reconfigure(bool silent = false); double CalcWeight(BaseFitEvt* evt); inline bool NeedsEventReWeight() { return true; }; #ifdef __GENIE_ENABLED__ -#ifndef __NO_GENIE_REWEIGHT__ +#ifndef __NO_REWEIGHT__ std::vector fGENIESysts; genie::rew::GReWeight* fGenieRW; //!< Genie RW Object #endif #endif }; #endif diff --git a/src/Reweight/NEUTWeightEngine.cxx b/src/Reweight/NEUTWeightEngine.cxx index 88b515f..966a09b 100644 --- a/src/Reweight/NEUTWeightEngine.cxx +++ b/src/Reweight/NEUTWeightEngine.cxx @@ -1,179 +1,178 @@ #include "NEUTWeightEngine.h" NEUTWeightEngine::NEUTWeightEngine(std::string name) { -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__) // Setup the NEUT Reweight engien fCalcName = name; NUIS_LOG(FIT, "Setting up NEUT RW : " << fCalcName); // Create RW Engine suppressing cout StopTalking(); fNeutRW = new neut::rew::NReWeight(); TDirectory *olddir = gDirectory; // get list of vetoed calc engines (just for debug really) std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fNeutRW_veto"); bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_res = rw_engine_list.find("xsec_res") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool nucl_casc = rw_engine_list.find("nucl_casc") == std::string::npos; bool nucl_piless = rw_engine_list.find("nucl_piless") == std::string::npos; // Activate each calc engine if (xsec_ccqe) fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE); if (xsec_res) fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES); if (xsec_ccres) fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES); if (xsec_coh) fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH); if (xsec_dis) fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS); if (xsec_ncel) fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL); if (xsec_nc) fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC); if (xsec_ncres) fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES); if (nucl_casc) fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc); if (nucl_piless) fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless); fNeutRW->Reconfigure(); olddir->cd(); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); // allow cout again StartTalking(); #else NUIS_ABORT("NEUT RW NOT ENABLED!"); #endif }; void NEUTWeightEngine::IncludeDial(std::string name, double startval) { -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__) // Get First enum int nuisenum = Reweight::ConvDial(name, kNEUT); // Setup Maps fEnumIndex[nuisenum]; // = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get Syst neut::rew::NSyst_t gensyst = NSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fNEUTSysts.push_back(gensyst); // Initialize dial NUIS_LOG(FIT, "Registering " << singlename << " dial."); fNeutRW->Systematics().Init(fNEUTSysts[index]); // If Absolute if (fIsAbsTwk) { NSystUncertainty::Instance()->SetUncertainty(fNEUTSysts[index], 1.0, 1.0); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } - #endif } void NEUTWeightEngine::SetDialValue(int nuisenum, double val) { -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__) std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; - NUIS_LOG(FIT, "Setting Dial Value for " << nuisenum << " " << i << " " - << indices[i] << " " - << fValues[indices[i]] - << " Enum: " << fNEUTSysts[indices[i]]); + NUIS_LOG(FIT, "Setting Dial Value for " + << nuisenum << " " << i << " " << indices[i] << " " + << fValues[indices[i]] + << " Enum: " << fNEUTSysts[indices[i]]); fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val); } #endif } void NEUTWeightEngine::SetDialValue(std::string name, double val) { -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__) std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; NUIS_LOG(FIT, "Setting Dial Value for " - << name << " = " << i << " " << indices[i] << " " - << fValues[indices[i]] - << " Enum: " << fNEUTSysts[indices[i]]); + << name << " = " << i << " " << indices[i] << " " + << fValues[indices[i]] + << " Enum: " << fNEUTSysts[indices[i]]); fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val); } #endif } void NEUTWeightEngine::Reconfigure(bool silent) { -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__) // Hush now... if (silent) StopTalking(); // Reconf fNeutRW->Reconfigure(); // if (LOG_LEVEL(DEB)){ fNeutRW->Print(); // } // Shout again if (silent) StartTalking(); #endif } double NEUTWeightEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__) // Skip Non NEUT if (evt->fType != kNEUT) return 1.0; // Hush now StopTalking(); // Fill NEUT Common blocks NEUTUtils::FillNeutCommons(evt->fNeutVect); // Call Weight calculation rw_weight = fNeutRW->CalcWeight(); // Speak Now StartTalking(); #endif // Return rw_weight return rw_weight; } diff --git a/src/Reweight/NEUTWeightEngine.h b/src/Reweight/NEUTWeightEngine.h index a4ad648..7cc7bb0 100644 --- a/src/Reweight/NEUTWeightEngine.h +++ b/src/Reweight/NEUTWeightEngine.h @@ -1,53 +1,55 @@ #ifndef WEIGHT_ENGINE_NEUT_H #define WEIGHT_ENGINE_NEUT_H #include "FitLogger.h" #ifdef __NEUT_ENABLED__ +#ifndef __NO_REWEIGHT__ +#include "NEUTInputHandler.h" #include "NReWeight.h" #include "NReWeightCasc.h" #include "NReWeightNuXSecCCQE.h" #include "NReWeightNuXSecCCRES.h" #include "NReWeightNuXSecCOH.h" #include "NReWeightNuXSecDIS.h" #include "NReWeightNuXSecNC.h" #include "NReWeightNuXSecNCEL.h" #include "NReWeightNuXSecNCRES.h" #include "NReWeightNuXSecRES.h" #include "NReWeightNuclPiless.h" #include "NSyst.h" #include "NSystUncertainty.h" #include "neutpart.h" #include "neutvect.h" -#include "NEUTInputHandler.h" +#endif #endif - +#include "FitWeight.h" #include "GeneratorUtils.h" #include "WeightEngineBase.h" -#include "FitWeight.h" class NEUTWeightEngine : public WeightEngineBase { public: - NEUTWeightEngine(std::string name); - ~NEUTWeightEngine() {}; + NEUTWeightEngine(std::string name); + ~NEUTWeightEngine(){}; - void IncludeDial(std::string name, double startval); + void IncludeDial(std::string name, double startval); - void SetDialValue(std::string name, double val); - void SetDialValue(int nuisenum, double val); + void SetDialValue(std::string name, double val); + void SetDialValue(int nuisenum, double val); - void Reconfigure(bool silent = false); + void Reconfigure(bool silent = false); - double CalcWeight(BaseFitEvt* evt); - - inline bool NeedsEventReWeight() { return true; }; + double CalcWeight(BaseFitEvt *evt); + inline bool NeedsEventReWeight() { return true; }; #ifdef __NEUT_ENABLED__ - std::vector fNEUTSysts; - neut::rew::NReWeight* fNeutRW; +#ifndef __NO_REWEIGHT__ + std::vector fNEUTSysts; + neut::rew::NReWeight *fNeutRW; +#endif #endif }; #endif diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx index 8429a6a..b381b80 100644 --- a/src/Reweight/WeightUtils.cxx +++ b/src/Reweight/WeightUtils.cxx @@ -1,592 +1,595 @@ #include "WeightUtils.h" #include "FitLogger.h" +#ifndef __NO_REWEIGHT__ #ifdef __T2KREW_ENABLED__ #include "T2KGenieReWeight.h" #include "T2KNIWGReWeight.h" #include "T2KNIWGUtils.h" #include "T2KNeutReWeight.h" #include "T2KNeutUtils.h" #include "T2KReWeight.h" using namespace t2krew; #endif #ifdef __NIWG_ENABLED__ #include "NIWGReWeight.h" #include "NIWGSyst.h" #endif #ifdef __NEUT_ENABLED__ #include "NReWeight.h" #include "NSyst.h" #endif #ifdef __NUWRO_REWEIGHT_ENABLED__ #include "NuwroReWeight.h" #include "NuwroSyst.h" #endif #ifdef __GENIE_ENABLED__ #ifdef GENIE_PRE_R3 #ifndef __NO_GENIE_REWEIGHT__ #include "ReWeight/GReWeight.h" #include "ReWeight/GSyst.h" #endif #else using namespace genie; #ifndef __NO_GENIE_REWEIGHT__ #include "RwFramework/GReWeight.h" #include "RwFramework/GSyst.h" using namespace genie::rew; #endif #endif #endif #ifdef __NOVA_ENABLED__ #include "NOvARwgtEngine.h" #endif +#endif // end of no reweight + #include "GlobalDialList.h" #include "ModeNormEngine.h" #include "NUISANCESyst.h" #include "OscWeightEngine.h" //******************************************************************** TF1 FitBase::GetRWConvFunction(std::string const &type, std::string const &name) { //******************************************************************** std::string dialfunc = "x"; std::string parType = type; double low = -10000.0; double high = 10000.0; if (parType.find("parameter") == std::string::npos) parType += "_parameter"; std::string line; std::ifstream card( (GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card") .c_str(), std::ifstream::in); while (std::getline(card >> std::ws, line, '\n')) { std::vector inputlist = GeneralUtils::ParseToStr(line, " "); // Check the line length if (inputlist.size() < 4) continue; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name if (inputlist[1].compare(name) != 0) continue; // inputlist[2] should be the units... ignore for now dialfunc = inputlist[3]; // High and low are optional, check whether they exist if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]); if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]); } TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high); return convfunc; } //******************************************************************** std::string FitBase::GetRWUnits(std::string const &type, std::string const &name) { //******************************************************************** std::string unit = "sig."; std::string parType = type; if (parType.find("parameter") == std::string::npos) { parType += "_parameter"; } std::string line; std::ifstream card( (GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card") .c_str(), std::ifstream::in); while (std::getline(card >> std::ws, line, '\n')) { std::vector inputlist = GeneralUtils::ParseToStr(line, " "); // Check the line length if (inputlist.size() < 3) continue; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name if (inputlist[1].compare(name) != 0) continue; unit = inputlist[2]; break; } return unit; } //******************************************************************** double FitBase::RWAbsToSigma(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.GetX(val); if (fabs(conv_val) < 1E-10) conv_val = 0.0; NUIS_LOG(FIT, "AbsToSigma(" << name << ") = " << val << " -> " << conv_val); return conv_val; } //******************************************************************** double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.Eval(val); return conv_val; } //******************************************************************** double FitBase::RWFracToSigma(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.GetX((val * f1.Eval(0.0))); if (fabs(conv_val) < 1E-10) conv_val = 0.0; return conv_val; } //******************************************************************** double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.Eval(val) / f1.Eval(0.0); return conv_val; } int FitBase::ConvDialType(std::string const &type) { if (!type.compare("neut_parameter")) return kNEUT; else if (!type.compare("niwg_parameter")) return kNIWG; else if (!type.compare("nuwro_parameter")) return kNUWRO; else if (!type.compare("t2k_parameter")) return kT2K; else if (!type.compare("genie_parameter")) return kGENIE; else if (!type.compare("custom_parameter")) return kCUSTOM; else if (!type.compare("norm_parameter")) return kNORM; else if (!type.compare("likeweight_parameter")) return kLIKEWEIGHT; else if (!type.compare("spline_parameter")) return kSPLINEPARAMETER; else if (!type.compare("osc_parameter")) return kOSCILLATION; else if (!type.compare("modenorm_parameter")) return kMODENORM; else if (!type.compare("nova_parameter")) return kNOvARWGT; else return kUNKNOWN; } std::string FitBase::ConvDialType(int type) { switch (type) { case kNEUT: { return "neut_parameter"; } case kNIWG: { return "niwg_parameter"; } case kNUWRO: { return "nuwro_parameter"; } case kT2K: { return "t2k_parameter"; } case kGENIE: { return "genie_parameter"; } case kNORM: { return "norm_parameter"; } case kCUSTOM: { return "custom_parameter"; } case kLIKEWEIGHT: { return "likeweight_parameter"; } case kSPLINEPARAMETER: { return "spline_parameter"; } case kOSCILLATION: { return "osc_parameter"; } case kMODENORM: { return "modenorm_parameter"; } case kNOvARWGT: { return "nova_parameter"; } default: return "unknown_parameter"; } } int FitBase::GetDialEnum(std::string const &type, std::string const &name) { return FitBase::GetDialEnum(FitBase::ConvDialType(type), name); } int FitBase::GetDialEnum(int type, std::string const &name) { int offset = type * 1000; int this_enum = Reweight::kNoDialFound; // Not Found NUIS_LOG(DEB, "Getting dial enum " << type << " " << name); // Select Types switch (type) { // NEUT DIAL TYPE case kNEUT: { -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__) int neut_enum = (int)neut::rew::NSyst::FromString(name); if (neut_enum != 0) { this_enum = neut_enum + offset; } #else this_enum = Reweight::kNoTypeFound; // Not enabled #endif break; } // NIWG DIAL TYPE case kNIWG: { -#ifdef __NIWG_ENABLED__ +#if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__) int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name); if (niwg_enum != 0) { this_enum = niwg_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } // NUWRO DIAL TYPE case kNUWRO: { -#ifdef __NUWRO_REWEIGHT_ENABLED__ +#if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__) int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name); if (nuwro_enum > 0) { this_enum = nuwro_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif } // GENIE DIAL TYPE case kGENIE: { -#if defined(__GENIE_ENABLED__) && !defined(__NO_GENIE_REWEIGHT__) +#if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__) int genie_enum = (int)genie::rew::GSyst::FromString(name); if (genie_enum > 0) { this_enum = genie_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } case kCUSTOM: { int custom_enum = 0; // PLACEHOLDER this_enum = custom_enum + offset; break; } // T2K DIAL TYPE case kT2K: { -#ifdef __T2KREW_ENABLED__ +#if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__) int t2k_enum = (int)t2krew::T2KSyst::FromString(name); if (t2k_enum > 0) { this_enum = t2k_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } case kNORM: { if (gNormEnums.find(name) == gNormEnums.end()) { gNormEnums[name] = gNormEnums.size() + 1 + offset; } this_enum = gNormEnums[name]; break; } case kLIKEWEIGHT: { if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) { gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset; } this_enum = gLikeWeightEnums[name]; break; } case kSPLINEPARAMETER: { if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) { gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset; } this_enum = gSplineParameterEnums[name]; break; } case kOSCILLATION: { #ifdef __PROB3PP_ENABLED__ int oscEnum = OscWeightEngine::SystEnumFromString(name); if (oscEnum != 0) { this_enum = oscEnum + offset; } #else this_enum = Reweight::kNoTypeFound; // Not enabled #endif } case kMODENORM: { size_t us_pos = name.find_first_of('_'); std::string numstr = name.substr(us_pos + 1); int mode_num = std::atoi(numstr.c_str()); NUIS_LOG(FTL, "Getting mode num " << mode_num); if (!mode_num) { NUIS_ABORT("Attempting to parse dial name: \"" << name << "\" as a mode norm dial but failed."); } this_enum = 60 + mode_num + offset; break; } } // If Not Enabled if (this_enum == Reweight::kNoTypeFound) { NUIS_ERR(FTL, "RW Engine not supported for " << FitBase::ConvDialType(type)); NUIS_ABORT("Check dial " << name); } // If Not Found if (this_enum == Reweight::kNoDialFound) { NUIS_ABORT("Dial " << name << " not found."); } return this_enum; } int Reweight::ConvDialType(std::string const &type) { return FitBase::ConvDialType(type); } std::string Reweight::ConvDialType(int type) { return FitBase::ConvDialType(type); } int Reweight::GetDialType(int type) { int t = (type / 1000); return t > kNOvARWGT ? Reweight::kNoDialFound : t; } int Reweight::RemoveDialType(int type) { return (type % 1000); } int Reweight::NEUTEnumFromName(std::string const &name) { -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__) int neutenum = (int)neut::rew::NSyst::FromString(name); return (neutenum > 0) ? neutenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NIWGEnumFromName(std::string const &name) { -#ifdef __NIWG_ENABLED__ +#if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__) int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name); return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NUWROEnumFromName(std::string const &name) { -#ifdef __NUWRO_REWEIGHT_ENABLED__ +#if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__) int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name); return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::GENIEEnumFromName(std::string const &name) { -#if defined(__GENIE_ENABLED__) && !defined(__NO_GENIE_REWEIGHT__) +#if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__) int genieenum = (int)genie::rew::GSyst::FromString(name); return (genieenum > 0) ? genieenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::T2KEnumFromName(std::string const &name) { -#ifdef __T2KREW_ENABLED__ +#if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__) int t2kenum = (int)t2krew::T2KSyst::FromString(name); return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::OscillationEnumFromName(std::string const &name) { #ifdef __PROB3PP_ENABLED__ int oscEnum = OscWeightEngine::SystEnumFromString(name); return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NUISANCEEnumFromName(std::string const &name, int type) { int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type); return nuisenum; } int Reweight::CustomEnumFromName(std::string const &name) { int custenum = Reweight::ConvertNUISANCEDial(name); return custenum; } int Reweight::ConvDial(std::string const &name, std::string const &type, bool exceptions) { return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions); } int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) { std::string name = GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given // Produce offset seperating each type. int offset = type * 1000; int genenum = Reweight::kNoDialFound; switch (type) { case kNEUT: genenum = NEUTEnumFromName(name); break; case kNIWG: genenum = NIWGEnumFromName(name); break; case kNUWRO: genenum = NUWROEnumFromName(name); break; case kGENIE: genenum = GENIEEnumFromName(name); break; case kT2K: genenum = T2KEnumFromName(name); break; case kCUSTOM: genenum = CustomEnumFromName(name); break; case kNORM: case kLIKEWEIGHT: case kSPLINEPARAMETER: case kNEWSPLINE: genenum = NUISANCEEnumFromName(name, type); break; case kOSCILLATION: genenum = OscillationEnumFromName(name); break; case kMODENORM: genenum = ModeNormEngine::SystEnumFromString(name); break; #ifdef __NOVA_ENABLED__ case kNOvARWGT: genenum = NOvARwgtEngine::GetWeightGeneratorIndex(name); break; #endif default: genenum = Reweight::kNoTypeFound; break; } // Throw if required. if (exceptions) { // If Not Enabled if (genenum == Reweight::kGeneratorNotBuilt) { NUIS_ERR(FTL, "RW Engine not supported for " << FitBase::ConvDialType(type)); NUIS_ABORT("Check dial " << name); } // If no type enabled if (genenum == Reweight::kNoTypeFound) { NUIS_ABORT("Type mismatch inside ConvDialEnum"); } // If Not Found if (genenum == Reweight::kNoDialFound) { NUIS_ABORT("Dial " << name << " not found."); } } // Add offset if no issue int nuisenum = genenum; if ((genenum != Reweight::kGeneratorNotBuilt) && (genenum != Reweight::kNoTypeFound) && (genenum != Reweight::kNoDialFound)) { nuisenum += offset; } // Now register dial Reweight::DialList().RegisterDialEnum(name, type, nuisenum); return nuisenum; } std::string Reweight::ConvDial(int nuisenum) { for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) { if (Reweight::DialList().fAllDialEnums[i] == nuisenum) { return Reweight::DialList().fAllDialNames[i]; } } NUIS_LOG(FIT, "Cannot find dial with enum = " << nuisenum); return ""; } diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx index c4fa68b..586cbed 100644 --- a/src/Statistical/StatUtils.cxx +++ b/src/Statistical/StatUtils.cxx @@ -1,1406 +1,1448 @@ // 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 "StatUtils.h" #include "GeneralUtils.h" #include "NuisConfig.h" #include "TH1D.h" //******************************************************************* Double_t StatUtils::GetChi2FromDiag(TH1D *data, TH1D *mc, TH1I *mask) { //******************************************************************* Double_t Chi2 = 0.0; - TH1D *calc_data = (TH1D *)data->Clone(); - TH1D *calc_mc = (TH1D *)mc->Clone(); + TH1D *calc_data = (TH1D *)data->Clone("calc_data"); + calc_data->SetDirectory(NULL); + TH1D *calc_mc = (TH1D *)mc->Clone("calc_mc"); + calc_mc->SetDirectory(NULL); // Add MC Error to data if required if (FitPar::Config().GetParB("addmcerror")) { for (int i = 0; i < calc_data->GetNbinsX(); i++) { double dterr = calc_data->GetBinError(i + 1); double mcerr = calc_mc->GetBinError(i + 1); if (dterr > 0.0) { calc_data->SetBinError(i + 1, sqrt(dterr * dterr + mcerr * mcerr)); } } } // Apply masking if required if (mask) { calc_data = ApplyHistogramMasking(data, mask); + calc_data->SetDirectory(NULL); calc_mc = ApplyHistogramMasking(mc, mask); + calc_mc->SetDirectory(NULL); } // Iterate over bins in X for (int i = 0; i < calc_data->GetNbinsX(); i++) { // Ignore bins with zero data or zero bin error if (calc_data->GetBinError(i + 1) <= 0.0 || calc_data->GetBinContent(i + 1) == 0.0) continue; // Take mc data difference double diff = calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1); double err = calc_data->GetBinError(i + 1); Chi2 += (diff * diff) / (err * err); } // cleanup delete calc_data; delete calc_mc; return Chi2; }; //******************************************************************* Double_t StatUtils::GetChi2FromDiag(TH2D *data, TH2D *mc, TH2I *map, TH2I *mask) { //******************************************************************* // Generate a simple map if (!map) map = GenerateMap(data); // Convert to 1D Histograms TH1D *data_1D = MapToTH1D(data, map); TH1D *mc_1D = MapToTH1D(mc, map); TH1I *mask_1D = MapToMask(mask, map); // Calculate 1D chi2 from 1D Plots Double_t Chi2 = StatUtils::GetChi2FromDiag(data_1D, mc_1D, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; }; //******************************************************************* Double_t StatUtils::GetChi2FromCov(TH1D *data, TH1D *mc, TMatrixDSym *invcov, TH1I *mask, double data_scale, - double covar_scale) { + double covar_scale, TH1D *outchi2perbin) { //******************************************************************* Double_t Chi2 = 0.0; TMatrixDSym *calc_cov = (TMatrixDSym *)invcov->Clone(); TH1D *calc_data = (TH1D *)data->Clone(); TH1D *calc_mc = (TH1D *)mc->Clone(); // If a mask if applied we need to apply it before the matrix is inverted if (mask) { calc_cov = ApplyInvertedMatrixMasking(invcov, mask); calc_data = ApplyHistogramMasking(data, mask); calc_mc = ApplyHistogramMasking(mc, mask); } // Add MC Error to data if required if (FitPar::Config().GetParB("statutils.addmcerror")) { // Make temp cov TMatrixDSym *newcov = StatUtils::GetInvert(calc_cov); // Add MC err to diag for (int i = 0; i < calc_data->GetNbinsX(); i++) { double mcerr = calc_mc->GetBinError(i + 1) * sqrt(covar_scale); double oldval = (*newcov)(i, i); NUIS_LOG(FIT, "Adding cov stat " << mcerr * mcerr << " to " << (*newcov)(i, i)); (*newcov)(i, i) = oldval + mcerr * mcerr; } // Reset the calc_cov to new invert delete calc_cov; calc_cov = GetInvert(newcov); // Delete the tempcov delete newcov; } calc_data->Scale(data_scale); calc_mc->Scale(data_scale); (*calc_cov) *= covar_scale; // iterate over bins in X (i,j) NUIS_LOG(DEB, "START Chi2 Calculation================="); for (int i = 0; i < calc_data->GetNbinsX(); i++) { + double ibin_contrib = 0; NUIS_LOG(DEB, "[CHI2] i = " << i << " [" << calc_data->GetXaxis()->GetBinLowEdge(i + 1) << " -- " << calc_data->GetXaxis()->GetBinUpEdge(i + 1) << "]."); for (int j = 0; j < calc_data->GetNbinsX(); j++) { NUIS_LOG(DEB, "[CHI2]\t j = " << i << " [" << calc_data->GetXaxis()->GetBinLowEdge(j + 1) << " -- " << calc_data->GetXaxis()->GetBinUpEdge(j + 1) << "]."); if ((calc_data->GetBinContent(i + 1) != 0 || calc_mc->GetBinContent(i + 1) != 0) && ((*calc_cov)(i, j) != 0)) { NUIS_LOG(DEB, "[CHI2]\t\t Chi2 contribution (i,j) = (" << i << "," << j << ")"); NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(i) = " << calc_data->GetBinContent(i + 1) << " - " << calc_mc->GetBinContent(i + 1) << " = " << (calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1))); NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(j) = " << calc_data->GetBinContent(j + 1) << " - " << calc_mc->GetBinContent(j + 1) << " = " << (calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1))); NUIS_LOG(DEB, "[CHI2]\t\t Covar = " << (*calc_cov)(i, j)); NUIS_LOG(DEB, "[CHI2]\t\t Cont chi2 = " << ((calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1)) * (*calc_cov)(i, j) * (calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1))) << " " << Chi2); - - Chi2 += + double bin_cont = ((calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1)) * (*calc_cov)(i, j) * (calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1))); + Chi2 += bin_cont; + ibin_contrib += bin_cont; } else { NUIS_LOG(DEB, "Skipping chi2 contribution (i,j) = (" << i << "," << j << "), Data = " << calc_data->GetBinContent(i + 1) << ", MC = " << calc_mc->GetBinContent(i + 1) << ", Cov = " << (*calc_cov)(i, j)); Chi2 += 0.; } } + if (outchi2perbin) { + outchi2perbin->SetBinContent(i + 1, ibin_contrib); + } } // Cleanup delete calc_cov; delete calc_data; delete calc_mc; return Chi2; } //******************************************************************* Double_t StatUtils::GetChi2FromCov(TH2D *data, TH2D *mc, TMatrixDSym *invcov, - TH2I *map, TH2I *mask) { + TH2I *map, TH2I *mask, TH2D *outchi2perbin) { //******************************************************************* // Generate a simple map if (!map) { map = StatUtils::GenerateMap(data); } // Convert to 1D Histograms TH1D *data_1D = MapToTH1D(data, map); TH1D *mc_1D = MapToTH1D(mc, map); TH1I *mask_1D = MapToMask(mask, map); + TH1D *outchi2perbin_1D = NULL; + TH1D *outchi2perbin_map_1D = NULL; + + if (outchi2perbin) { + outchi2perbin_1D = MapToTH1D(outchi2perbin, map); + for (Int_t xbi_it = 0; xbi_it < outchi2perbin->GetXaxis()->GetNbins(); + ++xbi_it) { + for (Int_t ybi_it = 0; ybi_it < outchi2perbin->GetYaxis()->GetNbins(); + ++ybi_it) { + int gbin = outchi2perbin->GetBin(xbi_it + 1, ybi_it + 1); + // std::cout << " gbin " << gbin << " corresponds to " + // << " x: " << (xbi_it + 1) << ", y: " << (ybi_it + 1) + // << std::endl; + outchi2perbin->SetBinContent(xbi_it + 1, ybi_it + 1, gbin); + } + } + outchi2perbin_map_1D = MapToTH1D(outchi2perbin, map); + } // Calculate 1D chi2 from 1D Plots - Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D); + Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D, 1, + 1E76, outchi2perbin_1D); + if (outchi2perbin) { + for (int xbi_it = 0; xbi_it < outchi2perbin_1D->GetXaxis()->GetNbins(); + ++xbi_it) { + // std::cout << " adding chi2 " + // << outchi2perbin_1D->GetBinContent(xbi_it + 1) + // << " from 1d bin " << (xbi_it + 1) << " to gbin " + // << outchi2perbin_map_1D->GetBinContent(xbi_it + 1) << std::endl; + outchi2perbin->SetBinContent( + outchi2perbin_map_1D->GetBinContent(xbi_it + 1), + outchi2perbin_1D->GetBinContent(xbi_it + 1)); + } + } // CleanUp delete data_1D; delete mc_1D; delete mask_1D; + delete outchi2perbin_1D; + delete outchi2perbin_map_1D; return Chi2; } //******************************************************************* Double_t StatUtils::GetChi2FromSVD(TH1D *data, TH1D *mc, TMatrixDSym *cov, TH1I *mask) { //******************************************************************* Double_t Chi2 = 0.0; TMatrixDSym *calc_cov = (TMatrixDSym *)cov->Clone(); TH1D *calc_data = (TH1D *)data->Clone(); TH1D *calc_mc = (TH1D *)mc->Clone(); // If a mask if applied we need to apply it before the matrix is inverted if (mask) { calc_cov = StatUtils::ApplyMatrixMasking(cov, mask); calc_data = StatUtils::ApplyHistogramMasking(data, mask); calc_mc = StatUtils::ApplyHistogramMasking(mc, mask); } // Decompose matrix TDecompSVD LU = TDecompSVD((*calc_cov)); LU.Decompose(); TMatrixDSym *cov_U = new TMatrixDSym(calc_data->GetNbinsX(), LU.GetU().GetMatrixArray(), ""); TVectorD *cov_S = new TVectorD(LU.GetSig()); // Apply basis rotation before adding up chi2 Double_t rotated_difference = 0.0; for (int i = 0; i < calc_data->GetNbinsX(); i++) { rotated_difference = 0.0; // Rotate basis of Data - MC for (int j = 0; j < calc_data->GetNbinsY(); j++) rotated_difference += (calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)) * (*cov_U)(j, i); // Divide by rotated error cov_S Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i); } // Cleanup delete calc_cov; delete calc_data; delete calc_mc; delete cov_U; delete cov_S; return Chi2; } //******************************************************************* Double_t StatUtils::GetChi2FromSVD(TH2D *data, TH2D *mc, TMatrixDSym *cov, TH2I *map, TH2I *mask) { //******************************************************************* // Generate a simple map if (!map) map = StatUtils::GenerateMap(data); // Convert to 1D Histograms TH1D *data_1D = MapToTH1D(data, map); TH1D *mc_1D = MapToTH1D(mc, map); TH1I *mask_1D = MapToMask(mask, map); // Calculate from 1D Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; } //******************************************************************* double StatUtils::GetChi2FromEventRate(TH1D *data, TH1D *mc, TH1I *mask) { //******************************************************************* // If just an event rate, for chi2 just use Poission Likelihood to calculate // the chi2 component double chi2 = 0.0; TH1D *calc_data = (TH1D *)data->Clone(); TH1D *calc_mc = (TH1D *)mc->Clone(); // Apply masking if required if (mask) { calc_data = ApplyHistogramMasking(data, mask); calc_mc = ApplyHistogramMasking(mc, mask); } // Iterate over bins in X for (int i = 0; i < calc_data->GetNbinsX(); i++) { double dt = calc_data->GetBinContent(i + 1); double mc = calc_mc->GetBinContent(i + 1); if (mc <= 0) continue; if (dt <= 0) { // Only add difference chi2 += 2 * (mc - dt); } else { // Do the chi2 for Poisson distributions chi2 += 2 * (mc - dt + (dt * log(dt / mc))); } /* LOG(REC)<<"Evt Chi2 cont = "<Clone(); // If a mask is provided we need to apply it before getting NDOF if (mask) { calc_hist = StatUtils::ApplyHistogramMasking(hist, mask); } // NDOF is defined as total number of bins with non-zero errors Int_t NDOF = 0; for (int i = 0; i < calc_hist->GetNbinsX(); i++) { if (calc_hist->GetBinError(i + 1) > 0.0) NDOF++; } delete calc_hist; return NDOF; }; //******************************************************************* Int_t StatUtils::GetNDOF(TH2D *hist, TH2I *map, TH2I *mask) { //******************************************************************* Int_t NDOF = 0; if (!map) map = StatUtils::GenerateMap(hist); for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (mask->GetBinContent(i + 1, j + 1)) continue; if (map->GetBinContent(i + 1, j + 1) <= 0) continue; NDOF++; } } return NDOF; }; //******************************************************************* TH1D *StatUtils::ThrowHistogram(TH1D *hist, TMatrixDSym *cov, bool throwdiag, TH1I *mask) { //******************************************************************* TH1D *calc_hist = (TH1D *)hist->Clone((std::string(hist->GetName()) + "_THROW").c_str()); TMatrixDSym *calc_cov = (TMatrixDSym *)cov->Clone(); Double_t correl_val = 0.0; // If a mask if applied we need to apply it before the matrix is decomposed if (mask) { calc_cov = ApplyMatrixMasking(cov, mask); calc_hist = ApplyHistogramMasking(calc_hist, mask); } // If a covariance is provided we need a preset random vector and a decomp std::vector rand_val; TMatrixDSym *decomp_cov = NULL; if (cov) { for (int i = 0; i < hist->GetNbinsX(); i++) { rand_val.push_back(gRandom->Gaus(0.0, 1.0)); } // Decomp the matrix decomp_cov = StatUtils::GetDecomp(calc_cov); } // iterate over bins for (int i = 0; i < hist->GetNbinsX(); i++) { // By Default the errors on the histogram are thrown uncorrelated to the // other errors /* if (throwdiag) { calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + \ gRandom->Gaus(0.0, 1.0) * calc_hist->GetBinError(i + 1)) ); } */ // If a covariance is provided that is also thrown if (cov) { correl_val = 0.0; for (int j = 0; j < hist->GetNbinsX(); j++) { correl_val += rand_val[j] * (*decomp_cov)(j, i); } calc_hist->SetBinContent( i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38)); } } delete calc_cov; delete decomp_cov; // return this new thrown data return calc_hist; }; //******************************************************************* TH2D *StatUtils::ThrowHistogram(TH2D *hist, TMatrixDSym *cov, TH2I *map, bool throwdiag, TH2I *mask) { //******************************************************************* // PLACEHOLDER!!!!!!!!! // Currently no support for throwing 2D Histograms from a covariance (void)hist; (void)cov; (void)map; (void)throwdiag; (void)mask; // /todo // Sort maps if required // Throw the covariance for a 1D plot // Unmap back to 2D Histogram return hist; } //******************************************************************* TH1D *StatUtils::ApplyHistogramMasking(TH1D *hist, TH1I *mask) { //******************************************************************* if (!mask) return ((TH1D *)hist->Clone()); // This masking is only sufficient for chi2 calculations, and will have dodgy // bin edges. // Get New Bin Count Int_t NBins = 0; for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) continue; NBins++; } // Make new hist std::string newmaskname = std::string(hist->GetName()) + "_MSKD"; TH1D *calc_hist = new TH1D(newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins); // fill new hist int binindex = 0; for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) { NUIS_LOG(REC, "Applying mask to bin " << i + 1 << " " << hist->GetName()); continue; } calc_hist->SetBinContent(binindex + 1, hist->GetBinContent(i + 1)); calc_hist->SetBinError(binindex + 1, hist->GetBinError(i + 1)); binindex++; } return calc_hist; }; //******************************************************************* TH2D *StatUtils::ApplyHistogramMasking(TH2D *hist, TH2I *mask) { //******************************************************************* TH2D *newhist = (TH2D *)hist->Clone(); if (!mask) return newhist; for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (mask->GetBinContent(i + 1, j + 1) > 0) { newhist->SetBinContent(i + 1, j + 1, 0.0); newhist->SetBinContent(i + 1, j + 1, 0.0); } } } return newhist; } //******************************************************************* TMatrixDSym *StatUtils::ApplyMatrixMasking(TMatrixDSym *mat, TH1I *mask) { //******************************************************************* if (!mask) return (TMatrixDSym *)(mat->Clone()); // Get New Bin Count Int_t NBins = 0; for (int i = 0; i < mask->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) continue; NBins++; } // make new matrix TMatrixDSym *calc_mat = new TMatrixDSym(NBins); int col, row; // Need to mask out bins in the current matrix row = 0; for (int i = 0; i < mask->GetNbinsX(); i++) { col = 0; // skip if masked if (mask->GetBinContent(i + 1) > 0.5) continue; for (int j = 0; j < mask->GetNbinsX(); j++) { // skip if masked if (mask->GetBinContent(j + 1) > 0.5) continue; (*calc_mat)(row, col) = (*mat)(i, j); col++; } row++; } return calc_mat; }; //******************************************************************* TMatrixDSym *StatUtils::ApplyMatrixMasking(TMatrixDSym *mat, TH2D *data, TH2I *mask, TH2I *map) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1I *mask_1D = StatUtils::MapToMask(mask, map); TMatrixDSym *newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D); delete mask_1D; return newmat; } //******************************************************************* TMatrixDSym *StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym *mat, TH1I *mask) { //******************************************************************* TMatrixDSym *new_mat = GetInvert(mat); TMatrixDSym *masked_mat = ApplyMatrixMasking(new_mat, mask); TMatrixDSym *inverted_mat = GetInvert(masked_mat); delete masked_mat; delete new_mat; return inverted_mat; }; //******************************************************************* TMatrixDSym *StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym *mat, TH2D *data, TH2I *mask, TH2I *map) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1I *mask_1D = StatUtils::MapToMask(mask, map); TMatrixDSym *newmat = ApplyInvertedMatrixMasking(mat, mask_1D); delete mask_1D; return newmat; } //******************************************************************* TMatrixDSym *StatUtils::GetInvert(TMatrixDSym *mat) { //******************************************************************* TMatrixDSym *new_mat = (TMatrixDSym *)mat->Clone(); // Check for diagonal bool non_diagonal = false; for (int i = 0; i < new_mat->GetNrows(); i++) { for (int j = 0; j < new_mat->GetNrows(); j++) { if (i == j) continue; if ((*new_mat)(i, j) != 0.0) { non_diagonal = true; break; } } } // If diag, just flip the diag if (!non_diagonal or new_mat->GetNrows() == 1) { for (int i = 0; i < new_mat->GetNrows(); i++) { if ((*new_mat)(i, i) != 0.0) (*new_mat)(i, i) = 1.0 / (*new_mat)(i, i); else (*new_mat)(i, i) = 0.0; } return new_mat; } // Invert full matrix TDecompSVD LU = TDecompSVD((*new_mat)); new_mat = new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), ""); return new_mat; } //******************************************************************* TMatrixDSym *StatUtils::GetDecomp(TMatrixDSym *mat) { //******************************************************************* TMatrixDSym *new_mat = (TMatrixDSym *)mat->Clone(); int nrows = new_mat->GetNrows(); // Check for diagonal bool diagonal = true; for (int i = 0; i < nrows; i++) { for (int j = 0; j < nrows; j++) { if (i == j) continue; if ((*new_mat)(i, j) != 0.0) { diagonal = false; break; } } } // If diag, just flip the diag if (diagonal or nrows == 1) { for (int i = 0; i < nrows; i++) { if ((*new_mat)(i, i) > 0.0) (*new_mat)(i, i) = sqrt((*new_mat)(i, i)); else (*new_mat)(i, i) = 0.0; } return new_mat; } TDecompChol LU = TDecompChol(*new_mat); LU.Decompose(); delete new_mat; TMatrixDSym *dec_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), ""); return dec_mat; } //******************************************************************* void StatUtils::ForceNormIntoCovar(TMatrixDSym *&mat, TH1D *hist, double norm) { //******************************************************************* if (!mat) mat = MakeDiagonalCovarMatrix(hist); int nbins = mat->GetNrows(); TMatrixDSym *new_mat = new TMatrixDSym(nbins); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { double valx = hist->GetBinContent(i + 1) * 1E38; double valy = hist->GetBinContent(j + 1) * 1E38; (*new_mat)(i, j) = (*mat)(i, j) + norm * norm * valx * valy; } } // Swap the two delete mat; mat = new_mat; return; }; //******************************************************************* void StatUtils::ForceNormIntoCovar(TMatrixDSym *mat, TH2D *data, double norm, TH2I *map) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1D *data_1D = MapToTH1D(data, map); StatUtils::ForceNormIntoCovar(mat, data_1D, norm); delete data_1D; return; } //******************************************************************* TMatrixDSym *StatUtils::MakeDiagonalCovarMatrix(TH1D *data, double scaleF) { //******************************************************************* TMatrixDSym *newmat = new TMatrixDSym(data->GetNbinsX()); for (int i = 0; i < data->GetNbinsX(); i++) { (*newmat)(i, i) = data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF; } return newmat; } //******************************************************************* TMatrixDSym *StatUtils::MakeDiagonalCovarMatrix(TH2D *data, TH2I *map, double scaleF) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1D *data_1D = MapToTH1D(data, map); return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF); }; //******************************************************************* void StatUtils::SetDataErrorFromCov(TH1D *DataHist, TMatrixDSym *cov, double scale, bool ErrorCheck) { //******************************************************************* // Check if (ErrorCheck) { if (cov->GetNrows() != DataHist->GetNbinsX()) { NUIS_ERR( FTL, "Nrows in cov don't match nbins in DataHist for SetDataErrorFromCov"); NUIS_ERR(FTL, "Nrows = " << cov->GetNrows()); NUIS_ABORT("Nbins = " << DataHist->GetNbinsX()); } } // Set bin errors form cov diag // Check if the errors are set bool ErrorsSet = false; for (int i = 0; i < DataHist->GetNbinsX(); i++) { if (ErrorsSet == true) break; if (DataHist->GetBinError(i + 1) != 0 && DataHist->GetBinContent(i + 1) > 0) ErrorsSet = true; } // Now loop over if (ErrorsSet && ErrorCheck) { for (int i = 0; i < DataHist->GetNbinsX(); i++) { double DataHisterr = DataHist->GetBinError(i + 1); double coverr = sqrt((*cov)(i, i)) * scale; // Check that the errors are within 1% of eachother if (fabs(DataHisterr - coverr) / DataHisterr > 0.01) { NUIS_ERR(WRN, "Data error does not match covariance error for bin " << i + 1 << " (" << DataHist->GetXaxis()->GetBinLowEdge(i + 1) << "-" << DataHist->GetXaxis()->GetBinLowEdge(i + 2) << ")"); NUIS_ERR(WRN, "Data error: " << DataHisterr); NUIS_ERR(WRN, "Cov error: " << coverr); } } // Else blindly trust the covariance } else { for (int i = 0; i < DataHist->GetNbinsX(); i++) { DataHist->SetBinError(i + 1, sqrt((*cov)(i, i)) * scale); } } return; } //******************************************************************* void StatUtils::SetDataErrorFromCov(TH2D *data, TMatrixDSym *cov, TH2I *map, double scale, bool ErrorCheck) { //******************************************************************* // Check if (ErrorCheck) { if (cov->GetNrows() != data->GetNbinsX() * data->GetNbinsY()) { NUIS_ERR( FTL, "Nrows in cov don't match nbins in data for SetDataNUIS_ERRorFromCov"); NUIS_ERR(FTL, "Nrows = " << cov->GetNrows()); NUIS_ABORT("Nbins = " << data->GetNbinsX()); } } // Set bin errors form cov diag // Check if the errors are set bool ErrorsSet = false; for (int i = 0; i < data->GetNbinsX(); i++) { for (int j = 0; j < data->GetNbinsX(); j++) { if (ErrorsSet == true) break; if (data->GetBinError(i + 1, j + 1) != 0) ErrorsSet = true; } } // Create map if required if (!map) map = StatUtils::GenerateMap(data); // Set Bin Errors from cov diag int count = 0; for (int i = 0; i < data->GetNbinsX(); i++) { for (int j = 0; j < data->GetNbinsY(); j++) { if (data->GetBinContent(i + 1, j + 1) == 0.0) continue; // If we have errors on our histogram the map is good count = map->GetBinContent(i + 1, j + 1) - 1; double dataerr = data->GetBinError(i + 1, j + 1); double coverr = sqrt((*cov)(count, count)) * scale; // Check that the errors are within 1% of eachother if (ErrorsSet && ErrorCheck) { if (fabs(dataerr - coverr) / dataerr > 0.01) { NUIS_ERR(WRN, "Data error does not match covariance error for bin " << i + 1 << " (" << data->GetXaxis()->GetBinLowEdge(i + 1) << "-" << data->GetXaxis()->GetBinLowEdge(i + 2) << ")"); NUIS_ERR(WRN, "Data error: " << dataerr); NUIS_ERR(WRN, "Cov error: " << coverr); } } else { data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale); } } } // Delete the map now that we don't need it // Woops, it's needed elsewhere there! (Grrrrr) // map->Delete(); } TMatrixDSym *StatUtils::ExtractShapeOnlyCovar(TMatrixDSym *full_covar, TH1 *data_hist, double data_scale) { int nbins = full_covar->GetNrows(); TMatrixDSym *shape_covar = new TMatrixDSym(nbins); // Check nobody is being silly if (data_hist->GetNbinsX() != nbins) { NUIS_ERR(WRN, "Inconsistent matrix and data histogram passed to " "StatUtils::ExtractShapeOnlyCovar!"); NUIS_ERR(WRN, "data_hist has " << data_hist->GetNbinsX() << " matrix has " << nbins); int err_bins = data_hist->GetNbinsX(); if (nbins > err_bins) err_bins = nbins; for (int i = 0; i < err_bins; ++i) { NUIS_ERR(WRN, "Matrix diag. = " << (*full_covar)(i, i) << " data = " << data_hist->GetBinContent(i + 1)); } return NULL; } double total_data = 0; double total_covar = 0; // Initial loop to calculate some constants for (int i = 0; i < nbins; ++i) { total_data += data_hist->GetBinContent(i + 1) * data_scale; for (int j = 0; j < nbins; ++j) { total_covar += (*full_covar)(i, j); } } if (total_data == 0 || total_covar == 0) { NUIS_ERR(WRN, "Stupid matrix or data histogram passed to " "StatUtils::ExtractShapeOnlyCovar! Ignoring..."); return NULL; } NUIS_LOG(SAM, "Norm error = " << sqrt(total_covar) / total_data); // Now loop over and calculate the shape-only matrix for (int i = 0; i < nbins; ++i) { double data_i = data_hist->GetBinContent(i + 1) * data_scale; for (int j = 0; j < nbins; ++j) { double data_j = data_hist->GetBinContent(j + 1) * data_scale; double norm_term = data_i * data_j * total_covar / total_data / total_data; double mix_sum1 = 0; double mix_sum2 = 0; for (int k = 0; k < nbins; ++k) { mix_sum1 += (*full_covar)(k, j); mix_sum2 += (*full_covar)(i, k); } double mix_term1 = data_i * (mix_sum1 / total_data - total_covar * data_j / total_data / total_data); double mix_term2 = data_j * (mix_sum2 / total_data - total_covar * data_i / total_data / total_data); (*shape_covar)(i, j) = (*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term; } } return shape_covar; } //******************************************************************* TH2I *StatUtils::GenerateMap(TH2D *hist) { //******************************************************************* std::string maptitle = std::string(hist->GetName()) + "_MAP"; TH2I *map = new TH2I(maptitle.c_str(), maptitle.c_str(), hist->GetNbinsX(), 0, hist->GetNbinsX(), hist->GetNbinsY(), 0, hist->GetNbinsY()); Int_t index = 1; for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (hist->GetBinContent(i + 1, j + 1) > 0) { map->SetBinContent(i + 1, j + 1, index); index++; } else { map->SetBinContent(i + 1, j + 1, 0); } } } return map; } //******************************************************************* TH1D *StatUtils::MapToTH1D(TH2D *hist, TH2I *map) { //******************************************************************* if (!hist) return NULL; // Get N bins for 1D plot Int_t Nbins = map->GetMaximum(); std::string name1D = std::string(hist->GetName()) + "_1D"; // Make new 1D Hist TH1D *newhist = new TH1D(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins); // map bin contents for (int i = 0; i < map->GetNbinsX(); i++) { for (int j = 0; j < map->GetNbinsY(); j++) { if (map->GetBinContent(i + 1, j + 1) == 0) continue; newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1)); newhist->SetBinError(map->GetBinContent(i + 1, j + 1), hist->GetBinError(i + 1, j + 1)); } } // return return newhist; } //******************************************************************* TH1I *StatUtils::MapToMask(TH2I *hist, TH2I *map) { //******************************************************************* TH1I *newhist = NULL; if (!hist) return newhist; // Get N bins for 1D plot Int_t Nbins = map->GetMaximum(); std::string name1D = std::string(hist->GetName()) + "_1D"; // Make new 1D Hist newhist = new TH1I(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins); // map bin contents for (int i = 0; i < map->GetNbinsX(); i++) { for (int j = 0; j < map->GetNbinsY(); j++) { if (map->GetBinContent(i + 1, j + 1) == 0) continue; newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1)); } } // return return newhist; } TMatrixDSym *StatUtils::GetCovarFromCorrel(TMatrixDSym *correl, TH1D *data) { int nbins = correl->GetNrows(); TMatrixDSym *covar = new TMatrixDSym(nbins); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { (*covar)(i, j) = (*correl)(i, j) * data->GetBinError(i + 1) * data->GetBinError(j + 1); } } return covar; } //******************************************************************* TMatrixD *StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx, int dimy) { //******************************************************************* // Determine dim if (dimx == -1 and dimy == -1) { std::string line; std::ifstream covar(covfile.c_str(), std::ifstream::in); int row = 0; while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " "entries on this line: " << row); } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { column++; if (column > dimx) dimx = column; } row++; if (row > dimy) dimy = row; } } // Or assume symmetric if (dimx != -1 and dimy == -1) { dimy = dimx; } assert(dimy != -1 && " matrix dimy not set."); // Make new matrix TMatrixD *mat = new TMatrixD(dimx, dimy); std::string line; std::ifstream covar(covfile.c_str(), std::ifstream::in); int row = 0; while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " "entries on this line: " << row); } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { // Check Rows // assert(row > mat->GetNrows() && " covar rows doesn't match matrix // rows."); // assert(column > mat->GetNcols() && " covar cols doesn't match matrix // cols."); // Fill Matrix (*mat)(row, column) = (*iter); column++; } row++; } return mat; } //******************************************************************* TMatrixD *StatUtils::GetMatrixFromRootFile(std::string covfile, std::string histname) { //******************************************************************* std::string inputfile = covfile + ";" + histname; std::vector splitfile = GeneralUtils::ParseToStr(inputfile, ";"); if (splitfile.size() < 2) { NUIS_ABORT("No object name given!"); } // Get file TFile *tempfile = new TFile(splitfile[0].c_str(), "READ"); // Get Object TObject *obj = tempfile->Get(splitfile[1].c_str()); if (!obj) { NUIS_ABORT("Object " << splitfile[1] << " doesn't exist!"); } // Try casting TMatrixD *mat = dynamic_cast(obj); if (mat) { TMatrixD *newmat = (TMatrixD *)mat->Clone(); delete mat; tempfile->Close(); return newmat; } TMatrixDSym *matsym = dynamic_cast(obj); if (matsym) { TMatrixD *newmat = new TMatrixD(matsym->GetNrows(), matsym->GetNrows()); for (int i = 0; i < matsym->GetNrows(); i++) { for (int j = 0; j < matsym->GetNrows(); j++) { (*newmat)(i, j) = (*matsym)(i, j); } } delete matsym; tempfile->Close(); return newmat; } TH2D *mathist = dynamic_cast(obj); if (mathist) { TMatrixD *newmat = new TMatrixD(mathist->GetNbinsX(), mathist->GetNbinsX()); for (int i = 0; i < mathist->GetNbinsX(); i++) { for (int j = 0; j < mathist->GetNbinsX(); j++) { (*newmat)(i, j) = mathist->GetBinContent(i + 1, j + 1); } } delete mathist; tempfile->Close(); return newmat; } return NULL; } //******************************************************************* TMatrixDSym *StatUtils::GetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************* // Delete TempMat TMatrixD *tempmat = GetMatrixFromTextFile(covfile, dim, dim); // Make a symmetric covariance TMatrixDSym *newmat = new TMatrixDSym(tempmat->GetNrows()); for (int i = 0; i < tempmat->GetNrows(); i++) { for (int j = 0; j < tempmat->GetNrows(); j++) { (*newmat)(i, j) = (*tempmat)(i, j); } } delete tempmat; return newmat; } //******************************************************************* TMatrixDSym *StatUtils::GetCovarFromRootFile(std::string covfile, std::string histname) { //******************************************************************* TMatrixD *tempmat = GetMatrixFromRootFile(covfile, histname); TMatrixDSym *newmat = new TMatrixDSym(tempmat->GetNrows()); for (int i = 0; i < tempmat->GetNrows(); i++) { for (int j = 0; j < tempmat->GetNrows(); j++) { (*newmat)(i, j) = (*tempmat)(i, j); } } delete tempmat; return newmat; } diff --git a/src/Statistical/StatUtils.h b/src/Statistical/StatUtils.h index 435e6a7..d2f82c9 100644 --- a/src/Statistical/StatUtils.h +++ b/src/Statistical/StatUtils.h @@ -1,255 +1,255 @@ // 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 STATUTILS_H #define STATUTILS_H // C Includes #include #include #include #include #include #include #include #include #include "assert.h" // Root Includes #include "TH1D.h" #include "TH2I.h" #include "TH2D.h" #include "TFile.h" #include "TMatrixDSym.h" #include "TDecompSVD.h" #include "TMath.h" #include "TRandom3.h" #include "TDecompChol.h" #include "TGraphErrors.h" // Fit Includes #include "FitLogger.h" /*! * \addtogroup Utils * @{ */ //! Functions for handling statistics calculations namespace StatUtils{ /* Chi2 Functions */ //! Get Chi2 using diagonal bin errors from the histogram. Masking applied before calculation if mask provided. Double_t GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Get Chi2 using diagonal bin errors from the histogram. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromDiag(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); //! Get Chi2 using an inverted covariance for the data - Double_t GetChi2FromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL, double data_scale=1, double covar_scale=1E76); + Double_t GetChi2FromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL, double data_scale=1, double covar_scale=1E76, TH1D *outchi2perbin=NULL); //! Get Chi2 using an inverted covariance for the data //! Plots converted to 1D histograms before using 1D calculation. - Double_t GetChi2FromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL); + Double_t GetChi2FromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL, TH2D *outchi2perbin=NULL); //! Get Chi2 using an SVD method on the covariance before calculation. //! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work. Double_t GetChi2FromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask=NULL); //! Get Chi2 using an SVD method on the covariance before calculation. //! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromSVD( TH2D* data, TH2D* mc, TMatrixDSym* cov, TH2I* map=NULL, TH2I* mask=NULL); //! Get Chi2 using only the raw event rates given in each bin using a -2LL method. Double_t GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Get Chi2 using only the raw event rates given in each bin using a -2LL method. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromEventRate(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); // Likelihood Functions //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromDiag(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromDiag(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromSVD( TH2D* data, TH2D* mc, TMatrixDSym* cov, TH2I* map=NULL, TH2I* mask=NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromEventRate(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromEventRate(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); /* NDOF Functions */ //! Return 1D Histogram NDOF considering masking and empty bins Int_t GetNDOF(TH1D* hist, TH1I* mask=NULL); //! Return 2D Histogram NDOF considering masking and empty bins Int_t GetNDOF(TH2D* hist, TH2I* map=NULL, TH2I* mask=NULL); /* Fake Data Functions */ //! Given a full covariance for a 1D data set throw the decomposition to generate fake data. //! throwdiag determines whether diagonal statistical errors are thrown. //! If no covariance is provided only statistical errors are thrown. TH1D* ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag=true, TH1I* mask=NULL); //! Given a full covariance for a 2D data set throw the decomposition to generate fake data. //! Plots are converted to 1D histograms and the 1D ThrowHistogram is used, before being converted back to 2D histograms. TH2D* ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map=NULL, bool throwdiag=true, TH2I* mask=NULL); /* Masking Functions */ //! Given a mask histogram, mask out any bins in hist with non zero entries in mask. TH1D* ApplyHistogramMasking(TH1D* hist, TH1I* mask); //! Given a mask histogram, mask out any bins in hist with non zero entries in mask. TH2D* ApplyHistogramMasking(TH2D* hist, TH2I* mask); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance, before recalculating its inverse. TMatrixDSym* ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH1I* mask); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance, before recalculating its inverse. //! Converts to 1D data before using the 1D ApplyInvertedMatrixMasking function and converting back to 2D. TMatrixDSym* ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map=NULL); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance TMatrixDSym* ApplyMatrixMasking(TMatrixDSym* mat, TH1I* mask); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance //! Converts to 1D data before using the 1D ApplyInvertedMatrixMasking function and converting back to 2D. TMatrixDSym* ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map=NULL); /* Covariance Handling Functions */ //! Return inverted matrix of TMatrixDSym TMatrixDSym* GetInvert(TMatrixDSym* mat); //! Return Cholesky Decomposed matrix of TMatrixDSym TMatrixDSym* GetDecomp(TMatrixDSym* mat); //! Return full covariances TMatrixDSym* GetCovarFromCorrel(TMatrixDSym* correl, TH1D* data); //! Given a normalisation factor for a dataset add in a new normalisation term to the covariance. void ForceNormIntoCovar(TMatrixDSym*& mat, TH1D* data, double norm); //! Given a normalisation factor for a dataset add in a new normalisation term to the covariance. //! Convertes 2D to 1D, before using 1D ForceNormIntoCovar void ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, TH2I* map=NULL); //! Given a dataset generate an uncorrelated covariance matrix using the bin errors. TMatrixDSym* MakeDiagonalCovarMatrix(TH1D* data, double scaleF=1E38); //! Given a dataset generate an uncorrelated covariance matrix using the bin errors. TMatrixDSym* MakeDiagonalCovarMatrix(TH2D* data, TH2I* map=NULL, double scaleF=1E38); //! Given a covariance set the errors in each bin on the data from the covariance diagonals. void SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, double scale=1.0, bool ErrorCheck = true); //! Given a covariance set the errors in each bin on the data from the covariance diagonals. void SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map=NULL, double scale=1.0, bool ErrorCheck = true); //! Given a covariance, extracts the shape-only matrix using the method from the MiniBooNE TN TMatrixDSym* ExtractShapeOnlyCovar(TMatrixDSym* full_covar, TH1* data_hist, double data_scale=1E38); /* Mapping Functions */ //! If no map is provided for the 2D histogram, generate one by counting up through the bins along x and y. TH2I* GenerateMap(TH2D* hist); //! Apply a map to a 2D histogram converting it into a 1D histogram. TH1D* MapToTH1D(TH2D* hist, TH2I* map); //! Apply a map to a 2D mask convering it into a 1D mask. TH1I* MapToMask(TH2I* hist, TH2I* map); /// \brief Read TMatrixD from a text file /// /// - covfile = full path to text file /// - dimx = x dimensions of matrix /// - dimy = y dimensions of matrix /// /// Format of textfile should be: \n /// cov_11 cov_12 ... cov_1N \n /// cov_21 cov_22 ... cov_2N \n /// ... ... ... ... \n /// cov_N1 ... ... cov_NN \n /// /// If no dimensions are given, dimx and dimy are determined from rows/columns /// inside textfile. /// /// If only dimx is given a symmetric matrix is assumed. TMatrixD* GetMatrixFromTextFile(std::string covfile, int dimx=-1, int dimy=-1); /// \brief Read TMatrixD from a ROOT file /// /// - covfile = full path to root file (+';histogram') /// - histname = histogram name /// /// If no histogram name is given function assumes it has been appended /// covfile path as: \n - /// 'covfile.root;histname' + /// 'covfile.root;histname' /// /// histname can point to a TMatrixD object, a TMatrixDSym object, or /// a TH2D object. TMatrixD* GetMatrixFromRootFile(std::string covfile, std::string histname=""); /// \brief Calls GetMatrixFromTextFile and turns it into a TMatrixDSym TMatrixDSym* GetCovarFromTextFile(std::string covfile, int dim); /// \brief Calls GetMatrixFromRootFile and turns it into a TMatrixDSym TMatrixDSym* GetCovarFromRootFile(std::string covfile, std::string histname); }; /*! @} */ #endif diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx index 0cc4155..956401e 100644 --- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx +++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx @@ -1,87 +1,90 @@ // 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 . -*******************************************************************************/ + * 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 "T2K_SignalDef.h" #include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h" - -//******************************************************************** -T2K_CC0pinp_STV_XSec_1Ddat_nu::T2K_CC0pinp_STV_XSec_1Ddat_nu(nuiskey samplekey) { //******************************************************************** +T2K_CC0pinp_STV_XSec_1Ddat_nu::T2K_CC0pinp_STV_XSec_1Ddat_nu( + nuiskey samplekey) { + //******************************************************************** // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddat_nu sample. \n" \ - "Target: CH \n" \ - "Flux: T2K 2.5 degree off-axis (ND280) \n" \ - "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" \ - " p_mu > 250 MeV \n" \ - " cth_p > 0.6 \n" \ - " cth_mu > -0.4 \n"; + std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddat_nu sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" + " p_mu > 250 MeV \n" + " cth_p > 0.6 \n" + " cth_mu > -0.4 \n" + "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#delta#it{#alpha}_{T} (GeV c^{-1})"); - fSettings.SetYTitle("#frac{d#sigma}{d#delta#it{#alpha}_{T}} (cm^{2} nucleon^{-1} rads^{-1})"); + fSettings.SetYTitle( + "#frac{d#sigma}{d#delta#it{#alpha}_{T}} (cm^{2} nucleon^{-1} rads^{-1})"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 50.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddat_nu"); - // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.dat"); + // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + + // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.dat"); - fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/datResults.root;Result" ); - fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/datResults.root;Correlation_Matrix" ); + fSettings.SetDataInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/datResults.root;Result"); + fSettings.SetCovarInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/datResults.root;Correlation_Matrix"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / - (double(fNEvents) * TotalIntegratedFlux("width")); + (double(fNEvents) * TotalIntegratedFlux("width")); // Plot Setup ------------------------------------------------------- - //SetDataFromTextFile( fSettings.GetDataInput() ); - //ScaleData(1E-38); - //SetCovarFromDiagonal(); - - SetDataFromRootFile( fSettings.GetDataInput() ); - SetCorrelationFromRootFile(fSettings.GetCovarInput() ); + // SetDataFromTextFile( fSettings.GetDataInput() ); + // ScaleData(1E-38); + // SetCovarFromDiagonal(); + SetDataFromRootFile(fSettings.GetDataInput()); + SetCorrelationFromRootFile(fSettings.GetCovarInput()); // Final setup --------------------------------------------------- FinaliseMeasurement(); - }; void T2K_CC0pinp_STV_XSec_1Ddat_nu::FillEventVariables(FitEvent *event) { fXVar = FitUtils::Get_STV_dalphat(event, 14, true); return; }; //******************************************************************** bool T2K_CC0pinp_STV_XSec_1Ddat_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax); } diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx index e14e2b6..cdc60dd 100644 --- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx +++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx @@ -1,88 +1,90 @@ // 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 . -*******************************************************************************/ + * 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 "T2K_SignalDef.h" #include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h" - -//******************************************************************** -T2K_CC0pinp_STV_XSec_1Ddphit_nu::T2K_CC0pinp_STV_XSec_1Ddphit_nu(nuiskey samplekey) { //******************************************************************** +T2K_CC0pinp_STV_XSec_1Ddphit_nu::T2K_CC0pinp_STV_XSec_1Ddphit_nu( + nuiskey samplekey) { + //******************************************************************** // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddphit_nu sample. \n" \ - "Target: CH \n" \ - "Flux: T2K 2.5 degree off-axis (ND280) \n" \ - "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n"\ - " p_mu > 250 MeV \n"\ - " cth_p > 0.6 \n"\ - " cth_mu > -0.4 \n";\ + std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddphit_nu sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" + " p_mu > 250 MeV \n" + " cth_p > 0.6 \n" + " cth_mu > -0.4 \n" + "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#delta#it{#phi}_{T} (GeV c^{-1})"); - fSettings.SetYTitle("#frac{d#sigma}{d#delta#it{#phi}_{T}} (cm^{2} nucleon^{-1} rads^{-1})"); + fSettings.SetYTitle( + "#frac{d#sigma}{d#delta#it{#phi}_{T}} (cm^{2} nucleon^{-1} rads^{-1})"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 50.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddphit_nu"); - //fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.dat"); - - fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dphitResults.root;Result" ); - fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dphitResults.root;Correlation_Matrix" ); + // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + + // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.dat"); + + fSettings.SetDataInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/dphitResults.root;Result"); + fSettings.SetCovarInput( + FitPar::GetDataBase() + + "/T2K/CC0pi/STV/dphitResults.root;Correlation_Matrix"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / - (double(fNEvents) * TotalIntegratedFlux("width")); + (double(fNEvents) * TotalIntegratedFlux("width")); // Plot Setup ------------------------------------------------------- - //SetDataFromTextFile( fSettings.GetDataInput() ); - //ScaleData(1E-38); - //SetCovarFromDiagonal(); + // SetDataFromTextFile( fSettings.GetDataInput() ); + // ScaleData(1E-38); + // SetCovarFromDiagonal(); - SetDataFromRootFile( fSettings.GetDataInput() ); - SetCorrelationFromRootFile(fSettings.GetCovarInput() ); + SetDataFromRootFile(fSettings.GetDataInput()); + SetCorrelationFromRootFile(fSettings.GetCovarInput()); // Final setup --------------------------------------------------- FinaliseMeasurement(); - }; - - - void T2K_CC0pinp_STV_XSec_1Ddphit_nu::FillEventVariables(FitEvent *event) { fXVar = FitUtils::Get_STV_dphit(event, 14, true); return; }; //******************************************************************** bool T2K_CC0pinp_STV_XSec_1Ddphit_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax); } diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx index 0b862ed..809db8a 100644 --- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx +++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx @@ -1,89 +1,90 @@ // 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 . -*******************************************************************************/ + * 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 "T2K_SignalDef.h" #include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h" - -//******************************************************************** -T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu(nuiskey samplekey) { //******************************************************************** +T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu( + nuiskey samplekey) { + //******************************************************************** // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n" \ - "Target: CH \n" \ - "Flux: T2K 2.5 degree off-axis (ND280) \n" \ - "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n"\ - " p_mu > 250 MeV \n"\ - " cth_p > 0.6 \n"\ - " cth_mu > -0.4 \n";\ + std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" + " p_mu > 250 MeV \n" + " cth_p > 0.6 \n" + " cth_mu > -0.4 \n" + "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#delta#it{p}_{T} (GeV c^{-1})"); - fSettings.SetYTitle("#frac{d#sigma}{d#delta#it{p}_{T}} (cm^{2} nucleon^{-1} GeV^{-1} c)"); + fSettings.SetYTitle( + "#frac{d#sigma}{d#delta#it{p}_{T}} (cm^{2} nucleon^{-1} GeV^{-1} c)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 50.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddpt_nu"); - //fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat"); - - fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dptResults.root;Result" ); - fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dptResults.root;Correlation_Matrix" ); + // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + + // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat"); + fSettings.SetDataInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/dptResults.root;Result"); + fSettings.SetCovarInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/dptResults.root;Correlation_Matrix"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / - (double(fNEvents) * TotalIntegratedFlux("width")); + (double(fNEvents) * TotalIntegratedFlux("width")); // Plot Setup ------------------------------------------------------- - //SetDataFromTextFile( fSettings.GetDataInput() ); - //SetCovarFromDiagonal(); - //ScaleData(1E-38); - - - SetDataFromRootFile( fSettings.GetDataInput() ); - SetCorrelationFromRootFile(fSettings.GetCovarInput() ); - //SetCovarianceFromRootFile(fSettings.GetCovarInput() ); + // SetDataFromTextFile( fSettings.GetDataInput() ); + // SetCovarFromDiagonal(); + // ScaleData(1E-38); + SetDataFromRootFile(fSettings.GetDataInput()); + SetCorrelationFromRootFile(fSettings.GetCovarInput()); + // SetCovarianceFromRootFile(fSettings.GetCovarInput() ); // Final setup --------------------------------------------------- FinaliseMeasurement(); - }; void T2K_CC0pinp_STV_XSec_1Ddpt_nu::FillEventVariables(FitEvent *event) { fXVar = FitUtils::Get_STV_dpt(event, 14, true) / 1000.0; return; }; //******************************************************************** bool T2K_CC0pinp_STV_XSec_1Ddpt_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax); }