Page MenuHomeHEPForge
Files F10881680

No OneTemporary
Actions

View Options

diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 3983c0a..c260ba6 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -1,178 +1,192 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
################################################################################
set(TARGETS_TO_BUILD)
if(USE_MINIMIZER)
add_executable(nuismin nuismin.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin)
target_link_libraries(nuismin ${MODULETargets})
target_link_libraries(nuismin ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuismin ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuismin PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissplines nuissplines.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines)
target_link_libraries(nuissplines ${MODULETargets})
target_link_libraries(nuissplines ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissplines ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissplines PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/Routines)
include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler)
include_directories(${CMAKE_SOURCE_DIR}/src/Genie)
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/FCN)
include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies)
include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance)
include_directories(${EXP_INCLUDE_DIRECTORIES})
if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE)
add_executable(nuwro_nuisance nuwro_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance)
target_link_libraries(nuwro_nuisance ${MODULETargets})
target_link_libraries(nuwro_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuwro_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuwro_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (USE_NEUT)
add_executable(neut_nuisance neut_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};neut_nuisance)
target_link_libraries(neut_nuisance ${MODULETargets})
target_link_libraries(neut_nuisance ${CMAKE_DEPENDLIB_FLAGS})
target_link_libraries(neut_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(neut_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (BUILD_GEVGEN)
add_executable(gevgen_nuisance gEvGen_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance)
target_link_libraries(gevgen_nuisance ${MODULETargets})
target_link_libraries(gevgen_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(gevgen_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${GENIE_INCLUDES}/Apps)
include_directories(${GENIE_INCLUDES}/FluxDrivers)
include_directories(${GENIE_INCLUDES}/EVGDrivers)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(gevgen_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
+
+
+ add_executable(gevgen_nuisance_mixed gEvGen_NUISANCE_MIXED.cxx)
+ set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance_mixed)
+ target_link_libraries(gevgen_nuisance_mixed ${MODULETargets})
+ target_link_libraries(gevgen_nuisance_mixed ${CMAKE_DEPENDLIB_FLAGS})
+ # target_link_libraries(gevgen_nuisance_mixed ${ROOT_LIBS})
+ include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
+ include_directories(${GENIE_INCLUDES}/Apps)
+ include_directories(${GENIE_INCLUDES}/FluxDrivers)
+ include_directories(${GENIE_INCLUDES}/EVGDrivers)
+ if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
+ set_target_properties(gevgen_nuisance_mixed PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
+ endif()
endif()
if (USE_GiBUU)
add_executable(DumpGiBUUEvents DumpGiBUUEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};DumpGiBUUEvents)
target_link_libraries(DumpGiBUUEvents ${MODULETargets})
target_link_libraries(DumpGiBUUEvents ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(DumpGiBUUEvents ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(DumpGiBUUEvents PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
add_executable(nuiscomp nuiscomp.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuiscomp)
target_link_libraries(nuiscomp ${MODULETargets})
target_link_libraries(nuiscomp ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuiscomp ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuiscomp PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuisflat nuisflat.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisflat)
target_link_libraries(nuisflat ${MODULETargets})
target_link_libraries(nuisflat ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuisflat ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuisflat PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissmear nuissmear.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissmear)
target_link_libraries(nuissmear ${MODULETargets})
target_link_libraries(nuissmear ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissmear ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissmear PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissyst nuissyst.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissyst)
target_link_libraries(nuissyst ${MODULETargets})
target_link_libraries(nuissyst ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissyst ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissyst PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
if(USE_GENIE)
add_executable(PrepareGENIE PrepareGENIE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE)
target_link_libraries(PrepareGENIE ${MODULETargets})
target_link_libraries(PrepareGENIE ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareGENIE ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareGENIE PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if(USE_NEUT)
add_executable(PrepareNEUT PrepareNEUT.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT)
target_link_libraries(PrepareNEUT ${MODULETargets})
target_link_libraries(PrepareNEUT ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNEUT ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNEUT PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
# PREPARE NUWRO
# Commented out for the time being until it is finished..
if(USE_NuWro)
add_executable(PrepareNuwro PrepareNuwroEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro)
target_link_libraries(PrepareNuwro ${MODULETargets})
target_link_libraries(PrepareNuwro ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNuwro ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNuwro PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin)
diff --git a/app/gEvGen_NUISANCE_MIXED.cxx b/app/gEvGen_NUISANCE_MIXED.cxx
new file mode 100644
index 0000000..c0969c7
--- /dev/null
+++ b/app/gEvGen_NUISANCE_MIXED.cxx
@@ -0,0 +1,1001 @@
+#ifdef __GENIE_ENABLED__
+//____________________________________________________________________________
+/*!
+
+\program gevgen
+
+\brief A simple 'generic' GENIE v+A event generation driver (gevgen).
+
+ It handles:
+ a) event generation for a fixed init state (v+A) at fixed energy, or
+ b) event generation for simple fluxes (specified either via some
+ functional form, tabular text file or a ROOT histogram) and for
+ simple 'geometries' (a target mix with its corresponding weights)
+
+ See the GENIE manual for other apps handling experiment-specific
+ event generation cases using the outputs of detailed neutrino flux
+ simulations and realistic detector geometry descriptions.
+
+ Syntax :
+ gevgen [-h]
+ [-r run#]
+ -n nev
+ -e energy (or energy range)
+ -p neutrino_pdg
+ -t target_pdg
+ [-f flux_description]
+ [-w]
+ [--seed random_number_seed]
+ [--cross-sections xml_file]
+ [--event-generator-list list_name]
+ [--tune genie_tune]
+ [--message-thresholds xml_file]
+ [--unphysical-event-mask mask]
+ [--event-record-print-level level]
+ [--mc-job-status-refresh-rate rate]
+ [--cache-file root_file]
+
+ Options :
+ [] Denotes an optional argument.
+ -h
+ Prints-out help on using gevgen and exits.
+ -n
+ Specifies the number of events to generate.
+ -r
+ Specifies the MC run number.
+ -e
+ Specifies the neutrino energy.
+ If what follows the -e option is a comma separated pair of values
+ it will be interpreted as an energy range for the flux specified
+ via the -f option (see below).
+ -p
+ Specifies the neutrino PDG code.
+ -t
+ Specifies the target PDG code (pdg format: 10LZZZAAAI) _or_ a target
+ mix (pdg codes with corresponding weights) typed as a comma-separated
+ list of pdg codes with the corresponding weight fractions in brackets,
+ eg code1[fraction1],code2[fraction2],...
+ For example, to use a target mix of 95% O16 and 5% H type:
+ `-t 1000080160[0.95],1000010010[0.05]'.
+ -f
+ Specifies the neutrino flux spectrum.
+ It can be any of:
+ -- A function:
+ eg ` -f x*x+4*exp(-x)'
+ -- A vector file:
+ The vector file should contain 2 columns corresponding to
+ energy,flux (see $GENIE/data/flux/ for few examples).
+ -- A 1-D ROOT histogram (TH1D):
+ The general syntax is `-f /full/path/file.root,object_name'
+ -w
+ Forces generation of weighted events.
+ This option is relevant only if a neutrino flux is specified.
+ Note that 'weighted' refers to the selection of the primary
+ flux neutrino + target that were forced to interact. A weighting
+ scheme for the generated kinematics of individual processes can
+ still be in effect if enabled..
+ ** Only use that option if you understand what it means **
+ --seed
+ Random number seed.
+ --cross-sections
+ Name (incl. full path) of an XML file with pre-computed
+ cross-section values used for constructing splines.
+ --event-generator-list
+ List of event generators to load in event generation drivers.
+ [default: "Default"].
+ --tune
+ Specifies a GENIE comprehensive neutrino interaction model tune.
+ [default: "Default"].
+ --message-thresholds
+ Allows users to customize the message stream thresholds.
+ The thresholds are specified using an XML file.
+ See $GENIE/config/Messenger.xml for the XML schema.
+ --unphysical-event-mask
+ Allows users to specify a 16-bit mask to allow certain types of
+ unphysical events to be written in the output file.
+ [default: all unphysical events are rejected]
+ --event-record-print-level
+ Allows users to set the level of information shown when the event
+ record is printed in the screen. See GHepRecord::Print().
+ --mc-job-status-refresh-rate
+ Allows users to customize the refresh rate of the status file.
+ --cache-file
+ Allows users to specify a cache file so that the cache can be
+ re-used in subsequent MC jobs.
+
+ *** See the User Manual for more details and examples. ***
+
+\author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
+ University of Liverpool & STFC Rutherford Appleton Lab
+
+\created October 05, 2004
+
+\cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
+ For the full text of the license visit http://copyright.genie-mc.org
+ or see $GENIE/LICENSE
+*/
+//____________________________________________________________________________
+
+#include <cstdlib>
+#include <cassert>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <map>
+
+#if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
+#include <fenv.h> // for `feenableexcept`
+#endif
+
+#include <TFile.h>
+#include <TTree.h>
+#include <TSystem.h>
+#include <TVector3.h>
+#include <TH1.h>
+#include <TF1.h>
+
+#include "Conventions/XmlParserStatus.h"
+#include "Conventions/GBuild.h"
+#include "Conventions/Controls.h"
+#include "Conventions/Units.h"
+
+#include "EVGCore/EventRecord.h"
+#include "EVGDrivers/GFluxI.h"
+#include "EVGDrivers/GEVGDriver.h"
+//#include "EVGDrivers/GMCJDriver.h"
+#include "GNUISANCEMCJDriver.h"
+#include "EVGDrivers/GMCJMonitor.h"
+#include "Interaction/Interaction.h"
+#include "Messenger/Messenger.h"
+#include "Ntuple/NtpWriter.h"
+#include "Ntuple/NtpMCFormat.h"
+#include "Numerical/RandomGen.h"
+#include "Numerical/Spline.h"
+#include "PDG/PDGCodes.h"
+#include "PDG/PDGUtils.h"
+#include "Utils/AppInit.h"
+#include "Utils/RunOpt.h"
+#include "Utils/XSecSplineList.h"
+#include "Utils/StringUtils.h"
+#include "Utils/PrintUtils.h"
+#include "Utils/SystemUtils.h"
+#include "Utils/CmdLnArgParser.h"
+
+//#include "FitBase/GNUISANCEFlux.h"
+#include "GNUISANCEFlux.h"
+
+//#ifdef __GENIE_FLUX_DRIVERS_ENABLED__
+//#ifdef __GENIE_GEOM_DRIVERS_ENABLED__
+#define __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
+#include "GNUISANCEFlux.h"
+#include "FluxDrivers/GMonoEnergeticFlux.h"
+#include "Geo/PointGeomAnalyzer.h"
+//#endif
+//#endif
+
+using std::string;
+using std::vector;
+using std::map;
+using std::ostringstream;
+
+using namespace genie;
+using namespace genie::controls;
+using namespace genie::flux;
+
+void GetCommandLineArgs (int argc, char ** argv);
+void Initialize (void);
+void PrintSyntax (void);
+
+#ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
+void GenerateEventsUsingFluxOrTgtMix();
+GeomAnalyzerI * GeomDriver (void);
+GFluxI * FluxDriver (void);
+GFluxI * MonoEnergeticFluxDriver (void);
+GFluxI * TH1FluxDriver (void);
+string ConvertTargetIDs (string);
+string ConvertFluxIDs (string);
+void ListTargetIDs(void);
+void ListFluxIDs(void);
+#endif
+
+void GenerateEventsAtFixedInitState (void);
+
+//Default options (override them using the command line arguments):
+int kDefOptNevents = 0; // n-events to generate
+NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // ntuple format
+Long_t kDefOptRunNu = 0; // default run number
+
+//User-specified options:
+int gOptNevents; // n-events to generate
+double gOptNuEnergy; // neutrino E, or min neutrino energy in spectrum
+double gOptNuEnergyRange;// energy range in input spectrum
+int gOptNuPdgCode; // neutrino PDG code
+map<int,double> gOptTgtMix; // target mix (each with its relative weight)
+Long_t gOptRunNu; // run number
+string gOptFlux; //
+bool gOptWeighted; //
+bool gOptUsingFluxOrTgtMix = false;
+long int gOptRanSeed; // random number seed
+string gOptInpXSecFile; // cross-section splines
+
+
+vector<int> gOptNuPDGs; // Beam PDGS
+vector<TH1D*> gOptNuFluxs; // Beam Fluxes
+vector<int> gOptTgtPDGs; // Target PDGS
+int gOptNumberNucleons; // Total number of target nucleons
+
+//____________________________________________________________________________
+int main(int argc, char ** argv)
+{
+ GetCommandLineArgs(argc,argv);
+ Initialize();
+
+ // throw on NaNs and Infs...
+#if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+#endif
+ //
+ // Generate neutrino events
+ //
+
+ if(gOptUsingFluxOrTgtMix) {
+#ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
+ GenerateEventsUsingFluxOrTgtMix();
+#else
+ LOG("gevgen", pERROR)
+ << "\n To be able to generate neutrino events from a flux and/or a target mix"
+ << "\n you need to add the following config options at your GENIE installation:"
+ << "\n --enable-flux-drivers --enable-geom-drivers \n" ;
+#endif
+ } else {
+ GenerateEventsAtFixedInitState();
+ }
+ return 0;
+}
+//____________________________________________________________________________
+void Initialize()
+{
+ // Initialization of random number generators, cross-section table,
+ // messenger thresholds, cache file
+ utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
+ utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
+ utils::app_init::RandGen(gOptRanSeed);
+ utils::app_init::XSecTable(gOptInpXSecFile, false);
+
+ // Set GHEP print level
+ GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
+}
+//____________________________________________________________________________
+void GenerateEventsAtFixedInitState(void)
+{
+ int neutrino = gOptNuPdgCode;
+ int target = gOptTgtMix.begin()->first;
+ double Ev = gOptNuEnergy;
+ TLorentzVector nu_p4(0.,0.,Ev,Ev); // px,py,pz,E (GeV)
+
+ // Create init state
+ InitialState init_state(target, neutrino);
+
+ // Create/config event generation driver
+ GEVGDriver evg_driver;
+ evg_driver.SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
+ evg_driver.SetUnphysEventMask(*RunOpt::Instance()->UnphysEventMask());
+ evg_driver.Configure(init_state);
+
+ // Initialize an Ntuple Writer
+ NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu);
+ ntpw.Initialize();
+
+ // Create an MC Job Monitor
+ GMCJMonitor mcjmonitor(gOptRunNu);
+ mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
+
+ LOG("gevgen", pNOTICE)
+ << "\n ** Will generate " << gOptNevents << " events for \n"
+ << init_state << " at Ev = " << Ev << " GeV";
+
+ // Generate events / print the GHEP record / add it to the ntuple
+ int ievent = 0;
+ while (ievent < gOptNevents) {
+ LOG("gevgen", pNOTICE)
+ << " *** Generating event............ " << ievent;
+
+ // generate a single event
+ EventRecord * event = evg_driver.GenerateEvent(nu_p4);
+
+ if(!event) {
+ LOG("gevgen", pNOTICE)
+ << "Last attempt failed. Re-trying....";
+ continue;
+ }
+
+ LOG("gevgen", pNOTICE)
+ << "Generated Event GHEP Record: " << *event;
+
+ // add event at the output ntuple, refresh the mc job monitor & clean up
+ ntpw.AddEventRecord(ievent, event);
+ mcjmonitor.Update(ievent,event);
+ ievent++;
+ delete event;
+ }
+
+ // Save the generated MC events
+ ntpw.Save();
+}
+//____________________________________________________________________________
+
+#ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
+//............................................................................
+void GenerateEventsUsingFluxOrTgtMix(void)
+{
+ // Get flux and geom drivers
+ GFluxI * flux_driver = FluxDriver();
+ GeomAnalyzerI * geom_driver = GeomDriver();
+
+ // Create the monte carlo job driver
+ GNUISANCEMCJDriver * mcj_driver = new GNUISANCEMCJDriver;
+ mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
+ mcj_driver->SetUnphysEventMask(*RunOpt::Instance()->UnphysEventMask());
+ mcj_driver->UseFluxDriver(flux_driver);
+ mcj_driver->UseGeomAnalyzer(geom_driver);
+ mcj_driver->Configure();
+ mcj_driver->UseSplines();
+ if(!gOptWeighted)
+ mcj_driver->ForceSingleProbScale();
+
+ // Initialize an Ntuple Writer to save GHEP records into a TTree
+ NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu);
+ ntpw.Initialize();
+
+ // Now Create flux histograms
+ TH1D* totalflux = static_cast<GNUISANCEFlux*>(flux_driver)->GetTotalSpectrum();
+ vector<TH1D*> fluxinputs = static_cast<GNUISANCEFlux*>(flux_driver)->GetSpectrum();
+ totalflux->SetNameTitle("nuisance_flux","NUISANCE Total Flux");
+
+ TH1D* totalevent = (TH1D*) totalflux->Clone();
+ totalevent->SetNameTitle("nuisance_events", "NUISANCE Total Events");
+ // totalevent->Reset();
+
+ TH1D* totalxsec = (TH1D*)totalflux->Clone();
+ totalxsec->SetNameTitle("nuisance_xsec", "NUISANCE Total XSec");
+ totalxsec->Reset();
+
+ // Calculate fractional target nucleon integral
+ double tfracintegral = 0.0;
+ for (UInt_t j = 0; j < gOptTgtPDGs.size(); j++){
+ int tpdg = gOptTgtPDGs[j];
+ int tnucl = tpdg % 10000 / 10;
+ double mixfrac = gOptTgtMix[tpdg];
+ tfracintegral += mixfrac * tnucl;
+ }
+
+ // Save the cross-section histograms
+ for (UInt_t i = 0; i < gOptNuPDGs.size(); i++){
+ for (UInt_t j = 0; j < gOptTgtPDGs.size(); j++){
+
+ int npdg = gOptNuPDGs[i];
+ TH1D* fluxhist = gOptNuFluxs[i];
+
+ int tpdg = gOptTgtPDGs[j];
+
+ // Determine the total number of these targets by using the tgt ratios backwards
+ double mixfrac = gOptTgtMix[tpdg];
+
+ int tnucl = tpdg % 10000 / 10; // Total Single Target Nucleons
+ double enucl = double(mixfrac * double(gOptNumberNucleons)); // Expected nucleons from weight fractions
+
+ // double tfrac = double(enucl) / double(tnucl); // Scaling factor.
+ double tfrac = mixfrac;
+
+
+ // Get Initial state driver
+ InitialState init_state(tpdg, npdg);
+ GEVGPool* gpool = mcj_driver->GetConfigPool();
+ GEVGDriver * evgdriver = gpool->FindDriver(init_state);
+ const Spline * totxsecspl = evgdriver->XSecSumSpline();
+
+ // make new empty xsec hist
+ TH1D* xsechist = (TH1D*) totalflux->Clone();
+ xsechist->SetNameTitle((init_state.AsString()+"_xsec").c_str(),init_state.AsString().c_str());
+ xsechist->Reset();
+
+ TH1D* eventhist = (TH1D*) totalflux->Clone();
+ eventhist->SetNameTitle((init_state.AsString()+"_events").c_str(),init_state.AsString().c_str());
+ eventhist->Reset();
+
+ std::cout << "Total Frac for " << init_state.AsString() << " = " << tfrac << std::endl;
+ std::cout << "TNUCL = " << tnucl << std::endl;
+ std::cout << "ENUCL = " << enucl << std::endl;
+ std::cout << "MIXFRAC = " << mixfrac << std::endl;
+ std::cout << "TPDG = " << tpdg << std::endl;
+
+ // Fill xsec hist
+ for (int k = 0; k < totalflux->GetNbinsX(); k++){
+ double avgxsec = 0.0;
+
+
+ int res = 1000;
+ for (int l = 0; l < res; l++){
+ double E = fluxhist->GetXaxis()->GetBinLowEdge(k+1) + (double(l)*fluxhist->GetXaxis()->GetBinWidth(k+1) / double(res));
+ double xsec = totxsecspl->Evaluate( E ) / (1E-38 * genie::units::cm2);
+ avgxsec += xsec;
+ }
+ avgxsec /= double(res);
+
+ // double E = fluxhist->GetXaxis()->GetBinCenter(k+1);
+ // avgxsec = totxsecspl->Evaluate( E ) / (1E-38 * genie::units::cm2);
+
+ xsechist->SetBinContent(k+1, avgxsec);
+ eventhist->SetBinContent(k+1, avgxsec * fluxhist->GetBinContent(k+1));
+
+ }
+
+ // Scale totals by target fractions
+ xsechist->Scale(tfrac);
+ eventhist->Scale(tfrac);
+
+ // Add to totals
+ totalxsec->Add(xsechist);
+ //totalevent->Add(eventhist);
+
+ std::cout << "Total XSec Hist for " << init_state.AsString() << " = " << xsechist->Integral("width") << std::endl;
+
+ // Scale by nnucleons
+ // xsechist->Scale(1.0 / double(gOptNumberNucleons));
+ // eventhist->Scale(1.0 / double(gOptNumberNucleons));
+
+ // Write our new xsec histto file
+ fluxhist->SetNameTitle((init_state.AsString()+"_flux").c_str(),init_state.AsString().c_str());
+ fluxhist->Write();
+ xsechist->Write();
+ eventhist->Write();
+
+ delete xsechist;
+ delete eventhist;
+
+ // sleep(10);
+ }
+ }
+
+ // Determine NUISANCE Style Histograms
+ totalevent->SetNameTitle("nuisance_events", "NUISANCE Total Events");
+ totalevent->Multiply(totalxsec);
+ // totalevent->Scale(1.0 / double(gOptNumberNucleons));
+ totalevent->Scale(1.0 / tfracintegral);
+ std::cout << "GOPTNumberNucleons = " << gOptNumberNucleons << std::endl;
+ std::cout << "Inclusive XSec Per Nucleon = " << totalevent->Integral("width") * 1E-38 / totalflux->Integral("width") << std::endl;
+ std::cout << "XSec Hist Integral = " << totalxsec->Integral("width") << std::endl;
+ // sleep(20);
+
+ // totalxsec->Scale(1.0/double(gOptNumberNucleons));
+ // totalevent->Scale(1.0/double(gOptNumberNucleons));
+
+ // Create an MC Job Monitor
+ GMCJMonitor mcjmonitor(gOptRunNu);
+ mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
+
+ // Generate events / print the GHEP record / add it to the ntuple
+ int ievent = 0;
+ while ( ievent < gOptNevents) {
+
+ LOG("gevgen", pNOTICE) << " *** Generating event............ " << ievent;
+
+ // generate a single event for neutrinos coming from the specified flux
+ EventRecord * event = mcj_driver->GenerateEvent();
+
+ LOG("gevgen", pNOTICE) << "Generated Event GHEP Record: " << *event;
+
+ // add event at the output ntuple, refresh the mc job monitor & clean-up
+ ntpw.AddEventRecord(ievent, event);
+ mcjmonitor.Update(ievent,event);
+ ievent++;
+ delete event;
+ }
+
+ // Save the generated MC events
+ ntpw.Save();
+
+
+ delete flux_driver;
+ delete geom_driver;
+ delete mcj_driver;;
+}
+//____________________________________________________________________________
+GeomAnalyzerI * GeomDriver(void)
+{
+// create a trivial point geometry with the specified target or target mix
+
+ GeomAnalyzerI * geom_driver = new geometry::PointGeomAnalyzer(gOptTgtMix);
+ return geom_driver;
+}
+//____________________________________________________________________________
+GFluxI * FluxDriver(void)
+{
+// create & configure one of the generic flux drivers
+//
+ GFluxI * flux_driver = 0;
+
+ if(gOptNuEnergyRange<0) flux_driver = MonoEnergeticFluxDriver();
+ else flux_driver = TH1FluxDriver();
+
+ return flux_driver;
+}
+//____________________________________________________________________________
+GFluxI * MonoEnergeticFluxDriver(void)
+{
+//
+//
+ flux::GMonoEnergeticFlux * flux =
+ new flux::GMonoEnergeticFlux(gOptNuEnergy, gOptNuPdgCode);
+ GFluxI * flux_driver = dynamic_cast<GFluxI *>(flux);
+ return flux_driver;
+}
+//____________________________________________________________________________
+GFluxI * TH1FluxDriver(void)
+{
+
+ // GeVGEN NUISANCE Can only Read flux files from ROOT
+ flux::GNUISANCEFlux * flux = new flux::GNUISANCEFlux;
+
+ // Initial beam directions
+ TVector3 bdir (0,0,1);
+ TVector3 bspot(0,0,0);
+
+ flux->SetNuDirection (bdir);
+ flux->SetBeamSpot (bspot);
+ flux->SetTransverseRadius (-1);
+
+ // Flux inputs
+ int flux_entries = 100000;
+
+ double emin = gOptNuEnergy;
+ double emax = gOptNuEnergy+gOptNuEnergyRange;
+ double de = gOptNuEnergyRange;
+
+ // check whether the input flux is a file or a functional form
+ bool input_is_root_file = gOptFlux.find(".root") != string::npos &&
+ gOptFlux.find(",") != string::npos;
+
+ if (!input_is_root_file){
+ std::cout << "ONLY ROOT INPUTS ALLOWED" << std::endl;
+ throw;
+ }
+
+ //
+ // ** extract specified flux histogram from the input root file
+ //
+ vector<string> fv = utils::str::Split(gOptFlux,",");
+
+ assert(fv.size()>=2);
+ assert( !gSystem->AccessPathName(fv[0].c_str()) );
+
+ // First entry is the root file
+ LOG("gevgen", pNOTICE) << "Getting input flux from root file: " << fv[0];
+ TFile * flux_file = new TFile(fv[0].c_str(),"read");
+
+ // Later entries are flux inputs
+ for (UInt_t i = 1; i < fv.size(); i++){
+ string fluxid = fv[i];
+
+ // Parse out the beam pdg and flux name
+ string flux_with_pdg = fluxid;
+ string::size_type open_bracket = flux_with_pdg.find("[");
+ string::size_type close_bracket = flux_with_pdg.find("]");
+ string::size_type ibeg = 0;
+ string::size_type iend = open_bracket;
+ string::size_type jbeg = open_bracket+1;
+ string::size_type jend = close_bracket-1;
+
+ string fluxname = flux_with_pdg.substr(ibeg,iend);
+ int fluxpdg = atoi(flux_with_pdg.substr(jbeg,jend).c_str());
+
+ // Get histogram
+ LOG("gevgen", pNOTICE) << "Flux name: " << fluxname << " pdg = " << fluxpdg;
+ TH1D * hst = (TH1D *)flux_file->Get(fluxname.c_str());
+ assert(hst);
+
+ LOG("gevgen", pNOTICE) << hst->GetEntries();
+
+
+ // Copy in the flux histogram from the root file and remove bins outside the emin,emax range
+ TH1D* spectrum = (TH1D*)hst->Clone();
+ spectrum->SetNameTitle("spectrum","neutrino_flux");
+ spectrum->SetDirectory(0);
+ for(int ibin = 1; ibin <= hst->GetNbinsX(); ibin++) {
+ if(spectrum->GetBinCenter(ibin) > emax ||
+ spectrum->GetBinCenter(ibin) < emin){
+ spectrum->SetBinContent(ibin, 0);
+ }
+ }
+
+ LOG("gevgen", pNOTICE) << spectrum->GetEntries();
+
+
+ // Add energy spectrum to the flux driver.
+ flux->AddEnergySpectrum( fluxpdg, spectrum );
+
+ gOptNuPDGs.push_back(fluxpdg);
+ gOptNuFluxs.push_back(spectrum);
+
+ }
+ // flux->AddAllFluxes();
+
+ // Close inputs
+ flux_file->Close();
+ delete flux_file;
+
+ // Return flux driver.
+ GFluxI * flux_driver = dynamic_cast<GFluxI *>(flux);
+ return flux_driver;
+}
+
+///____________________________________________________________________________
+void ListTargetIDs(){
+
+ // Keep in sync with ConvertTargetIDs
+ LOG("gevgen", pNOTICE) << "Possible Target IDs: \n"
+ << "\n H : " << ConvertTargetIDs("H")
+ << "\n C : " << ConvertTargetIDs("C")
+ << "\n CH : " << ConvertTargetIDs("CH")
+ << "\n CH2 : " << ConvertTargetIDs("CH2")
+ << "\n H2O : " << ConvertTargetIDs("H2O")
+ << "\n Fe : " << ConvertTargetIDs("Fe")
+ << "\n Pb : " << ConvertTargetIDs("Pb")
+ << "\n D2 : " << ConvertTargetIDs("D2")
+ << "\n D2-free : " << ConvertTargetIDs("D2-free");
+}
+
+
+//____________________________________________________________________________
+string ConvertTargetIDs(string id){
+
+ if (!id.compare("H")) return "1000010010";
+ else if (!id.compare("C")) return "1000060120";
+ else if (!id.compare("CH")) return "13,1000060120[0.9231],1000010010[0.0769]";
+ else if (!id.compare("CH2")) return "14,1000060120[0.8571],1000010010[0.1429]";
+ else if (!id.compare("H2O")) return "18,1000080160[0.8888],1000010010[0.1111]";
+ else if (!id.compare("Fe")) return "1000260560";
+ else if (!id.compare("Pb")) return "1000822070";
+ else if (!id.compare("D2")) return "1000010020";
+ else if (!id.compare("D2-free")) return "2,1000010010[0.5],1000000010[0.5]";
+ else return "";
+
+};
+
+///____________________________________________________________________________
+void ListFluxIDs(){
+
+ // Keep in sync with ConvertTargetIDs
+ LOG("gevgen", pNOTICE) << "Possible Flux IDs: \n"
+ << "\n MINERvA_fhc_numu : " << ConvertFluxIDs("MINERvA_fhc_numu")
+ << "\n MINERvA_fhc_numunumubar : " << ConvertFluxIDs("MINERvA_fhc_numunumubar")
+ << "\n MINERvA_fhc_nue : " << ConvertFluxIDs("MINERvA_fhc_nue")
+ << "\n MINERvA_fhc_nuenuebar : " << ConvertFluxIDs("MINERvA_fhc_nuenuebar")
+ << "\n MINERvA_fhc_all : " << ConvertFluxIDs("MINERvA_fhc_all")
+
+ << "\n MINERvA_rhc_numubar : " << ConvertFluxIDs("MINERvA_rhc_numubar")
+ << "\n MINERvA_rhc_numubarnumu : " << ConvertFluxIDs("MINERvA_rhc_numubarnumu")
+ << "\n MINERvA_rhc_nuebar : " << ConvertFluxIDs("MINERvA_rhc_nuebar")
+ << "\n MINERvA_rhc_nuebarnue : " << ConvertFluxIDs("MINERvA_rhc_nuebarnue")
+ << "\n MINERvA_rhc_all : " << ConvertFluxIDs("MINERvA_rhc_all")
+
+ << "\n ANL_fhc_numu : " << ConvertFluxIDs("ANL_fhc_numu")
+ << "\n BNL_fhc_numu : " << ConvertFluxIDs("BNL_fhc_numu")
+ << "\n BNL_fhc_numu_ALT1986 : " << ConvertFluxIDs("BNL_fhc_numu_ALT1986")
+ << "\n BNL_fhc_numu_ALT1981 : " << ConvertFluxIDs("BNL_fhc_numu_ALT1981")
+ << "\n BEBC_fhc_numu : " << ConvertFluxIDs("BEBC_fhc_numu")
+ << "\n FNAL_fhc_numu : " << ConvertFluxIDs("FNAL_fhc_numu")
+ << "\n FNAL_rhc_numub : " << ConvertFluxIDs("FNAL_rhc_numub")
+ << "\n GGM_fhc_numu : " << ConvertFluxIDs("GGM_fhc_numu");
+
+}
+
+
+//____________________________________________________________________________
+string ConvertFluxIDs(string id){
+
+ char * const var = getenv("NUISANCE");
+ if (!var) {
+ std::cout << "Cannot find top level directory! Set the NUISANCE environmental variable" << std::endl;
+ exit(-1);
+ }
+ string topnuisancedir = string(var);
+ string fluxfolder = topnuisancedir + "/data/flux/";
+ string inputs = "";
+
+ if (!id.compare("MINERvA_fhc_numu")) inputs="minerva_flux.root,numu_fhc[14]";
+ else if (!id.compare("MINERvA_fhc_numunumubar")) inputs="minerva_flux.root,numu_fhc[14],numubar_fhc[-14]";
+ else if (!id.compare("MINERvA_fhc_numu")) inputs="minerva_flux.root,nue_fhc[12]";
+ else if (!id.compare("MINERvA_fhc_nuenuebar")) inputs="minerva_flux.root,nue_fhc[12],nuebar_fhc[-12]";
+ else if (!id.compare("MINERvA_fhc_all")) inputs="minerva_flux.root,numu_fhc[14],numubar_fhc[-14],nue_fhc[12],nuebar_fhc[-12]";
+
+ else if (!id.compare("MINERvA_rhc_numubar")) inputs="minerva_flux.root,numubar_rhc[-14]";
+ else if (!id.compare("MINERvA_rhc_numubarnumu")) inputs="minerva_flux.root,numubar_rhc[-14],numu_rhc[14]";
+ else if (!id.compare("MINERvA_rhc_nuebar")) inputs="minerva_flux.root,nuebar_rhc[-12]";
+ else if (!id.compare("MINERvA_rhc_nuebarnue")) inputs="minerva_flux.root,nuebar_rhc[-12],nue_rhc[12]";
+ else if (!id.compare("MINERvA_rhc_all")) inputs="minerva_flux.root,numu_rhc[14],numubar_rhc[-14],nue_rhc[12],nuebar_rhc[-12]";
+
+ else if (!id.compare("ANL_fhc_numu")) inputs="ANL_1977_2horn_rescan.root,numu_flux[14]";
+ else if (!id.compare("BNL_fhc_numu")) inputs="BNL_NuInt02_rescan.root,numu_flux[14]";
+ else if (!id.compare("BNL_fhc_numu_ALT1986")) inputs="BNL_1986_flux-ALTERNATIVE.root,numu_flux[14]";
+ else if (!id.compare("BNL_fhc_numu_ALT1981")) inputs="BNL_CCQE_1981_rescan-ALTERNATIVE.root,numu_flux[14]";
+
+ else if (!id.compare("BEBC_fhc_numu")) inputs="BEBC_Wachsmuth_numubar_table.root,numu_flux[14]";
+ else if (!id.compare("FNAL_fhc_numu")) inputs="FNAL_CCinc_1982_nu_MCadj.root,numu_flux[14]";
+ else if (!id.compare("FNAL_rhc_numub")) inputs="FNAL_coh_1993_anu.root,numu_flux[-14]";
+ else if (!id.compare("GGM_fhc_numu")) inputs="GGM_nu_flux_1979_rescan.root,numu_flux[14]";
+ else return "";
+
+ return fluxfolder + inputs;
+
+};
+
+
+
+//............................................................................
+#endif
+
+
+//____________________________________________________________________________
+void GetCommandLineArgs(int argc, char ** argv)
+{
+ LOG("gevgen", pINFO) << "Parsing command line arguments";
+
+ // Common run options. Set defaults and read.
+ RunOpt::Instance()->EnableBareXSecPreCalc(true);
+ RunOpt::Instance()->ReadFromCommandLine(argc,argv);
+
+ // Parse run options for this app
+
+ CmdLnArgParser parser(argc,argv);
+
+ // help?
+ bool help = parser.OptionExists('h');
+ if(help) {
+ PrintSyntax();
+ exit(0);
+ }
+
+ // number of events
+ if( parser.OptionExists('n') ) {
+ LOG("gevgen", pINFO) << "Reading number of events to generate";
+ gOptNevents = parser.ArgAsInt('n');
+ } else {
+ LOG("gevgen", pINFO)
+ << "Unspecified number of events to generate - Using default";
+ gOptNevents = kDefOptNevents;
+ }
+
+ // run number
+ if( parser.OptionExists('r') ) {
+ LOG("gevgen", pINFO) << "Reading MC run number";
+ gOptRunNu = parser.ArgAsLong('r');
+ } else {
+ LOG("gevgen", pINFO) << "Unspecified run number - Using default";
+ gOptRunNu = kDefOptRunNu;
+ }
+
+ // flux functional form
+ bool using_flux = false;
+ if( parser.OptionExists('f') ) {
+ LOG("gevgen", pINFO) << "Reading flux function";
+ gOptFlux = parser.ArgAsString('f');
+
+ // Convert for known strings.
+ string fluxid = ConvertFluxIDs(gOptFlux);
+ if (!fluxid.empty()) gOptFlux = fluxid;
+
+ using_flux = true;
+ }
+
+ if(parser.OptionExists('s')) {
+ LOG("gevgen", pWARN)
+ << "-s option no longer available. Please read the revised code documentation";
+ gAbortingInErr = true;
+ exit(1);
+ }
+
+
+ // generate weighted events option (only relevant if using a flux)
+ gOptWeighted = parser.OptionExists('w');
+
+ // neutrino energy
+ if( parser.OptionExists('e') ) {
+ LOG("gevgen", pINFO) << "Reading neutrino energy";
+ string nue = parser.ArgAsString('e');
+
+ // is it just a value or a range (comma separated set of values)
+ if(nue.find(",") != string::npos) {
+ // split the comma separated list
+ vector<string> nurange = utils::str::Split(nue, ",");
+ assert(nurange.size() == 2);
+ double emin = atof(nurange[0].c_str());
+ double emax = atof(nurange[1].c_str());
+ assert(emax>emin && emin>=0);
+ gOptNuEnergy = emin;
+ gOptNuEnergyRange = emax-emin;
+ if(!using_flux) {
+ LOG("gevgen", pWARN)
+ << "No flux was specified but an energy range was input!";
+ LOG("gevgen", pWARN)
+ << "Events will be generated at fixed E = " << gOptNuEnergy << " GeV";
+ gOptNuEnergyRange = -1;
+ }
+ } else {
+ gOptNuEnergy = atof(nue.c_str());
+ gOptNuEnergyRange = -1;
+ }
+ } else {
+ LOG("gevgen", pFATAL) << "Unspecified neutrino energy - Exiting";
+ PrintSyntax();
+ exit(1);
+ }
+
+ // neutrino PDG code
+ if( parser.OptionExists('p') && !parser.OptionExists('f')) {
+ LOG("gevgen", pINFO) << "Reading neutrino PDG code";
+ gOptNuPdgCode = parser.ArgAsInt('p');
+ } else if (!parser.OptionExists('p') && !parser.OptionExists('f')){
+ LOG("gevgen", pFATAL) << "Unspecified neutrino PDG code or Flux Inputs - Exiting";
+ PrintSyntax();
+ exit(1);
+ }
+
+ // target mix (their PDG codes with their corresponding weights)
+ bool using_tgtmix = false;
+ if( parser.OptionExists('t') ) {
+ LOG("gevgen", pINFO) << "Reading target mix";
+ string stgtmix = parser.ArgAsString('t');
+ gOptTgtMix.clear();
+
+ // Check for ID Strings
+ string tgtids = ConvertTargetIDs(stgtmix);
+ if (!tgtids.empty()){
+ stgtmix = tgtids;
+ }
+
+ // Parse Targets
+ vector<string> tgtmix = utils::str::Split(stgtmix,",");
+ if(tgtmix.size()==1) {
+ int pdg = atoi(tgtmix[0].c_str());
+ double wgt = 1.0;
+ gOptTgtPDGs.push_back(pdg);
+ gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
+
+ // For single target number of free nucleons
+ // automatic
+ gOptNumberNucleons = pdg % 10000/ 10;
+
+ } else {
+ using_tgtmix = true;
+ vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
+
+ // For multiple targets N nucleons must be specified first.
+ gOptNumberNucleons = atoi((*tgtmix_iter).c_str());
+ tgtmix_iter++;
+
+ // Now remainder is the full target list
+ for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
+ string tgt_with_wgt = *tgtmix_iter;
+ string::size_type open_bracket = tgt_with_wgt.find("[");
+ string::size_type close_bracket = tgt_with_wgt.find("]");
+ string::size_type ibeg = 0;
+ string::size_type iend = open_bracket;
+ string::size_type jbeg = open_bracket+1;
+ string::size_type jend = close_bracket-1;
+ int pdg = atoi(tgt_with_wgt.substr(ibeg,iend).c_str());
+ double wgt = atof(tgt_with_wgt.substr(jbeg,jend).c_str());
+ LOG("Main", pNOTICE)
+ << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
+ gOptTgtPDGs.push_back(pdg);
+ gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
+ }//tgtmix_iter
+ }//>1
+
+ } else {
+ LOG("gevgen", pFATAL) << "Unspecified target PDG code - Exiting";
+ PrintSyntax();
+ exit(1);
+ }
+
+ gOptUsingFluxOrTgtMix = using_flux || using_tgtmix;
+
+ // random number seed
+ if( parser.OptionExists("seed") ) {
+ LOG("gevgen", pINFO) << "Reading random number seed";
+ gOptRanSeed = parser.ArgAsLong("seed");
+ } else {
+ LOG("gevgen", pINFO) << "Unspecified random number seed - Using default";
+ gOptRanSeed = -1;
+ }
+
+ // input cross-section file
+ if( parser.OptionExists("cross-sections") ) {
+ LOG("gevgen", pINFO) << "Reading cross-section file";
+ gOptInpXSecFile = parser.ArgAsString("cross-sections");
+ } else {
+ LOG("gevgen", pINFO) << "Unspecified cross-section file";
+ gOptInpXSecFile = "";
+ }
+
+ //
+ // print-out the command line options
+ //
+ LOG("gevgen", pNOTICE)
+ << "\n"
+ << utils::print::PrintFramedMesg("gevgen job configuration");
+ LOG("gevgen", pNOTICE)
+ << "MC Run Number: " << gOptRunNu;
+ if(gOptRanSeed != -1) {
+ LOG("gevgen", pNOTICE)
+ << "Random number seed: " << gOptRanSeed;
+ } else {
+ LOG("gevgen", pNOTICE)
+ << "Random number seed was not set, using default";
+ }
+ LOG("gevgen", pNOTICE)
+ << "Number of events requested: " << gOptNevents;
+ if(gOptInpXSecFile.size() > 0) {
+ LOG("gevgen", pNOTICE)
+ << "Using cross-section splines read from: " << gOptInpXSecFile;
+ } else {
+ LOG("gevgen", pNOTICE)
+ << "No input cross-section spline file";
+ }
+ LOG("gevgen", pNOTICE)
+ << "Flux: " << gOptFlux;
+ LOG("gevgen", pNOTICE)
+ << "Generate weighted events? " << gOptWeighted;
+ if(gOptNuEnergyRange>0) {
+ LOG("gevgen", pNOTICE)
+ << "Neutrino energy: ["
+ << gOptNuEnergy << ", " << gOptNuEnergy+gOptNuEnergyRange << "]";
+ } else {
+ LOG("gevgen", pNOTICE)
+ << "Neutrino energy: " << gOptNuEnergy;
+ }
+ LOG("gevgen", pNOTICE)
+ << "Neutrino code (PDG): " << gOptNuPdgCode;
+ LOG("gevgen", pNOTICE)
+ << "Target code (PDG) & weight fraction (in case of multiple targets): ";
+ map<int,double>::const_iterator iter;
+ for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
+ int tgtpdgc = iter->first;
+ double wgt = iter->second;
+ LOG("gevgen", pNOTICE)
+ << " >> " << tgtpdgc << " (weight fraction = " << wgt << ")";
+ }
+ LOG("gevgen", pNOTICE) << "\n";
+
+ LOG("gevgen", pNOTICE) << *RunOpt::Instance();
+
+}
+//____________________________________________________________________________
+void PrintSyntax(void)
+{
+ LOG("gevgen", pNOTICE)
+ << "\n\n" << "Syntax:" << "\n"
+ << "\n gevgen [-h]"
+ << "\n [-r run#]"
+ << "\n -n nev"
+ << "\n -e energy (or energy range) "
+ << "\n -p neutrino_pdg"
+ << "\n -t target_pdg "
+ << "\n [-f flux_description]"
+ << "\n [-w]"
+ << "\n [--seed random_number_seed]"
+ << "\n [--cross-sections xml_file]"
+ << "\n [--event-generator-list list_name]"
+ << "\n [--message-thresholds xml_file]"
+ << "\n [--unphysical-event-mask mask]"
+ << "\n [--event-record-print-level level]"
+ << "\n [--mc-job-status-refresh-rate rate]"
+ << "\n [--cache-file root_file]"
+ << "\n\n" ;
+ ListTargetIDs();
+ ListFluxIDs();
+
+
+}
+//____________________________________________________________________________
+#endif
diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
index a4136e2..31023b2 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,1087 +1,1010 @@
#include "JointFCN.h"
#include <stdio.h>
#include "FitUtils.h"
-// //***************************************************
-// JointFCN::JointFCN(std::string cardfile, TFile* outfile) {
-// //***************************************************
-
-// fOutputDir = gDirectory;
-// FitPar::Config().out = outfile;
-
-// fCardFile = cardfile;
-
-// LoadSamples(fCardFile);
-
-// fCurIter = 0;
-// fMCFilled = false;
-
-// fOutputDir->cd();
-
-// fIterationTree = NULL;
-// fDialVals = NULL;
-// fNDials = 0;
-
-// fUsingEventManager = FitPar::Config().GetParB("EventManager");
-// fOutputDir->cd();
-// };
+//***************************************************
JointFCN::JointFCN(TFile* outfile) {
+//***************************************************
+
fOutputDir = gDirectory;
if (outfile) FitPar::Config().out = outfile;
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
LoadSamples(samplekeys);
std::vector<nuiskey> covarkeys = Config::QueryKeys("covar");
LoadPulls(covarkeys);
fCurIter = 0;
fMCFilled = false;
- fIterationTree = NULL;
+ fIterationTree = false;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
+//***************************************************
JointFCN::JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile) {
+//***************************************************
+
fOutputDir = gDirectory;
if (outfile) FitPar::Config().out = outfile;
LoadSamples(samplekeys);
fCurIter = 0;
fMCFilled = false;
fOutputDir->cd();
- fIterationTree = NULL;
+ fIterationTree = false;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
//***************************************************
JointFCN::~JointFCN() {
//***************************************************
// Delete Samples
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
delete exp;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
delete pull;
}
// Sort Tree
if (fIterationTree) DestroyIterationTree();
if (fDialVals) delete fDialVals;
if (fSampleLikes) delete fSampleLikes;
};
//***************************************************
void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) {
- //***************************************************
-
- LOG(FIT) << " Creating new iteration tree! " << std::endl;
- if (fIterationTree && !name.compare(fIterationTree->GetName())) {
- fIterationTree->Reset();
- return;
- }
-
- fIterationTree = new TTree(name.c_str(), name.c_str());
-
- // Setup Main Branches
- fIterationTree->Branch("total_likelihood", &fLikelihood,
- "total_likelihood/D");
- fIterationTree->Branch("total_ndof", &fNDOF, "total_ndof/I");
+//***************************************************
- // Setup Sample Arrays
- int ninputs = fSamples.size() + fPulls.size();
- fSampleLikes = new double[ninputs];
- fSampleNDOF = new int[ninputs];
- fNDOF = GetNDOF();
+ LOG(FIT) << " Creating new iteration container! " << std::endl;
+ DestroyIterationTree();
+ fIterationTreeName = name;
- // Setup Sample Branches
- int count = 0;
- for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
+ // Add sample likelihoods and ndof
+ for (MeasListConstIter iter = fSamples.begin();
+ iter != fSamples.end();
iter++) {
- MeasurementBase* exp = *iter;
+ MeasurementBase* exp = *iter;
std::string name = exp->GetName();
- std::string liketag = name + "_likelihood";
- std::string ndoftag = name + "_ndof";
- fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
- (liketag + "/D").c_str());
- fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
- (ndoftag + "/D").c_str());
+ std::string liketag = name + "_likelihood";
+ fNameValues.push_back(liketag);
+ fCurrentValues.push_back(0.0);
- count++;
+ std::string ndoftag = name + "_ndof";
+ fNameValues.push_back(ndoftag);
+ fCurrentValues.push_back(0.0);
}
- for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
- ParamPull* pull = *iter;
+ // Add Pull terms
+ for (PullListConstIter iter = fPulls.begin();
+ iter != fPulls.end(); iter++) {
+ ParamPull* pull = *iter;
std::string name = pull->GetName();
+
std::string liketag = name + "_likelihood";
+ fNameValues.push_back(liketag);
+ fCurrentValues.push_back(0.0);
+
std::string ndoftag = name + "_ndof";
+ fNameValues.push_back(ndoftag);
+ fCurrentValues.push_back(0.0);
+ }
- fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
- (liketag + "/D").c_str());
- fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
- (ndoftag + "/D").c_str());
+ // Add Likelihoods
+ fNameValues.push_back("total_likelihood");
+ fCurrentValues.push_back(0.0);
- count++;
- }
+ fNameValues.push_back("total_ndof");
+ fCurrentValues.push_back(0.0);
- // Add Dial Branches
- std::vector<std::string> dials = rw->GetDialNames();
- fNDials = dials.size();
- fDialVals = new double[fNDials];
+ // Setup Containers
+ fSampleN = fSamples.size() + fPulls.size();
+ fSampleLikes = new double[fSampleN];
+ fSampleNDOF = new int[fSampleN];
- for (int i = 0; i < fNDials; i++) {
- fIterationTree->Branch(dials[i].c_str(), &fDialVals[i],
- (dials[i] + "/D").c_str());
+ // Add Dials
+ std::vector<std::string> dials = rw->GetDialNames();
+ for (size_t i = 0; i < dials.size(); i++){
+ fNameValues.push_back( dials[i] );
+ fCurrentValues.push_back( 0.0 );
}
+ fNDials = dials.size();
+ fDialVals = new double[fNDials];
+
+ // Set IterationTree Flag
+ fIterationTree = true;
+
}
//***************************************************
void JointFCN::DestroyIterationTree() {
- //***************************************************
+//***************************************************
+
+ fIterationCount.clear();
+ fCurrentValues.clear();
+ fNameValues.clear();
+ fIterationValues.clear();
- if (!fIterationTree) {
- delete fIterationTree;
- }
}
//***************************************************
void JointFCN::WriteIterationTree() {
- //***************************************************
+//***************************************************
+ LOG(FIT) << "Writing iteration tree" << std::endl;
+
+ // Make a new TTree
+ TTree* itree = new TTree(fIterationTreeName.c_str(),
+ fIterationTreeName.c_str());
+
+ double* vals = new double[fNameValues.size()];
+ int count = 0;
+
+ itree->Branch("iteration",&count,"Iteration/I");
+ for (int i = 0; i < fNameValues.size(); i++) {
+ itree->Branch( fNameValues[i].c_str(),
+ &vals[i],
+ (fNameValues[i] + "/D").c_str() );
+ }
+
+ // Fill Iterations
+ for (size_t i = 0; i < fIterationValues.size(); i++){
+ std::vector<double> itervals = fIterationValues[i];
+
+ // Fill iteration state
+ count = fIterationCount[i];
+ for (size_t j = 0; j < itervals.size(); j++){
+ vals[j] = itervals[j];
+ }
- if (!fIterationTree) {
- ERR(FTL) << "Can't save empty iteration tree!" << std::endl;
- throw;
+ // Save to TTree
+ itree->Fill();
}
- fIterationTree->Write();
+
+ // Write to file
+ itree->Write();
}
//***************************************************
void JointFCN::FillIterationTree(FitWeight* rw) {
- //***************************************************
+//***************************************************
- if (!fIterationTree) {
- ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << std::endl;
- throw;
+ // Loop over samples count
+ int count = 0;
+ for (int i = 0; i < fSampleN; i++){
+ fCurrentValues[count++] = fSampleLikes[i];
+ fCurrentValues[count++] = double(fSampleNDOF[i]);
}
+ // Fill Totals
+ fCurrentValues[count++] = fLikelihood;
+ fCurrentValues[count++] = double(fNDOF);
+
+ // Loop Over Parameter Counts
rw->GetAllDials(fDialVals, fNDials);
- fIterationTree->Fill();
+ for (int i = 0; i < fNDials; i++){
+ fCurrentValues[count++] = double(fDialVals[i]);
+ }
+
+ // Push Back Into Container
+ fIterationCount.push_back( fCurIter );
+ fIterationValues.push_back(fCurrentValues);
}
//***************************************************
double JointFCN::DoEval(const double* x) {
//***************************************************
// WEIGHT ENGINE
fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
FitBase::GetRW()->UpdateWeightEngine(x);
if (fDialChanged) {
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
}
if (LOG_LEVEL(REC)) {
FitBase::GetRW()->Print();
}
// SORT SAMPLES
ReconfigureSamples();
// GET TEST STAT
fLikelihood = GetLikelihood();
-
+ fNDOF = GetNDOF();
+
// PRINT PROGRESS
LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood
<< std::endl;
// UPDATE TREE
if (fIterationTree) FillIterationTree(FitBase::GetRW());
return fLikelihood;
}
//***************************************************
int JointFCN::GetNDOF() {
//***************************************************
int totaldof = 0;
int count = 0;
// Total number of Free bins in each MC prediction
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
int dof = exp->GetNDOF();
// Save Seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double dof = pull->GetLikelihood();
// Save seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Set Data Variable
- fNDOF = totaldof;
+ fSampleNDOF[count] = totaldof;
return totaldof;
}
//***************************************************
double JointFCN::GetLikelihood() {
//***************************************************
LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..."
<< " : "
<< "-2logL" << std::endl;
// Loop and add up likelihoods in an uncorrelated way
double like = 0.0;
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
double newlike = exp->GetLikelihood();
int ndof = exp->GetNDOF();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : "
<< newlike << "/" << ndof << std::endl;
// Add Weight Scaling
// like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName());
// Add to total
like += newlike;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double newlike = pull->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
// Add to total
like += newlike;
count++;
}
// Set Data Variable
fLikelihood = like;
return like;
};
void JointFCN::LoadSamples(std::vector<nuiskey> samplekeys) {
LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl;
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys[i];
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype = key.GetS("type");
std::string fakeData = "";
LOG(MIN) << "Loading Sample : " << samplename << std::endl;
fOutputDir->cd();
MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(key);
if (!NewLoadedSample) {
ERR(FTL) << "Could not load sample provided: " << samplename << std::endl;
ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
<< std::endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
}
+//***************************************************
void JointFCN::LoadPulls(std::vector<nuiskey> pullkeys) {
+//***************************************************
for (size_t i = 0; i < pullkeys.size(); i++) {
nuiskey key = pullkeys[i];
std::string pullname = key.GetS("name");
std::string pullfile = key.GetS("input");
std::string pulltype = key.GetS("type");
fOutputDir->cd();
- std::cout << "Creating Pull Term : " << std::endl;
- sleep(1);
fPulls.push_back(new ParamPull(pullname, pullfile, pulltype));
}
-
- // // Sample Inputs
- // if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls")
- // ||
- // !samplevect[0].compare("throws")) {
- // // Get all inputs
- // std::string name = samplevect[1];
- // std::string files = samplevect[2];
- // std::string type = "DEFAULT";
-
- // if (samplevect.size() > 3) {
- // type = samplevect[3];
- // } else if (!samplevect[0].compare("pull")) {
- // type = "GAUSPULL";
- // } else if (!samplevect[0].compare("throws")) {
- // type = "GAUSTHROWS";
- // }
-
- // // Create Pull Class
- // LOG(MIN) << "Loading up pull term: " << name << " << " << files << "
- // ("
- // << type << ")" << std::endl;
- // std::string fakeData = "";
- // fOutputDir->cd();
- // fPulls.push_back(new ParamPull(name, files, type));
- // }
}
-// //***************************************************
-// void JointFCN::LoadSamples(std::string cardinput)
-// //***************************************************
-// {
-// LOG(MIN) << "Initializing Samples" << std::endl;
-
-// // Read the card file here and load objects
-// std::string line;
-// std::ifstream card(cardinput.c_str(), ifstream::in);
-
-// // Make sure they are created in correct working DIR
-// fOutputDir->cd();
-
-// while (std::getline(card >> std::ws, line, '\n')) {
-// // Skip Empties
-// if (line.c_str()[0] == '#') continue;
-// if (line.empty()) continue;
-
-// // Parse line
-// std::vector<std::string> samplevect = GeneralUtils::ParseToStr(line, "
-// ");
-
-// // Sample Inputs
-// if (!samplevect[0].compare("sample")) {
-// // Get all inputs
-// std::string name = samplevect[1];
-// std::string files = samplevect[2];
-// std::string type = "DEFAULT";
-// if (samplevect.size() > 3) type = samplevect[3];
-
-// // Create Sample Class
-// LOG(MIN) << "Loading up sample: " << name << " << " << files << " ("
-// << type << ")" << std::endl;
-// std::string fakeData = "";
-// fOutputDir->cd();
-// MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(name,
-// files, type,
-// fakeData, FitBase::GetRW());
-
-// if (!NewLoadedSample) {
-// ERR(FTL) << "Could not load sample provided: " << name << std::endl;
-// ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
-// << std::endl;
-// throw;
-// } else {
-// fSamples.push_back(NewLoadedSample);
-// }
-// }
-
-// // Sample Inputs
-// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") ||
-// !samplevect[0].compare("throws")) {
-// // Get all inputs
-// std::string name = samplevect[1];
-// std::string files = samplevect[2];
-// std::string type = "DEFAULT";
-
-// if (samplevect.size() > 3) {
-// type = samplevect[3];
-// } else if (!samplevect[0].compare("pull")) {
-// type = "GAUSPULL";
-// } else if (!samplevect[0].compare("throws")) {
-// type = "GAUSTHROWS";
-// }
-
-// // Create Pull Class
-// LOG(MIN) << "Loading up pull term: " << name << " << " << files << " ("
-// << type << ")" << std::endl;
-// std::string fakeData = "";
-// fOutputDir->cd();
-// fPulls.push_back(new ParamPull(name, files, type));
-// }
-// }
-// card.close();
-// };
-
//***************************************************
void JointFCN::ReconfigureSamples(bool fullconfig) {
- //***************************************************
+//***************************************************
int starttime = time(NULL);
LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl;
// std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled <<
// std::endl;
// Event Manager Reconf
if (fUsingEventManager) {
if (!fullconfig and fMCFilled)
ReconfigureFastUsingManager();
else
ReconfigureUsingManager();
} else {
// Loop over all Measurement Classes
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
// If RW Either do signal or full reconfigure.
if (fDialChanged or !fMCFilled or fullconfig) {
if (!fullconfig and fMCFilled)
exp->ReconfigureFast();
else
exp->Reconfigure();
// If RW Not needed just do normalisation
} else {
exp->Renormalise();
}
}
}
// Loop over pulls and update
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Reconfigure();
}
fMCFilled = true;
LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in "
<< time(NULL) - starttime << "s" << std::endl;
fCurIter++;
}
//***************************************************
void JointFCN::ReconfigureSignal() {
- //***************************************************
- this->ReconfigureSamples(false);
+//***************************************************
+ ReconfigureSamples(false);
}
//***************************************************
void JointFCN::ReconfigureAllEvents() {
//***************************************************
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
ReconfigureSamples(true);
}
std::vector<InputHandlerBase*> JointFCN::GetInputList() {
std::vector<InputHandlerBase*> InputList;
fIsAllSplines = true;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
InputHandlerBase* inp = subsamples[i]->GetInput();
if (std::find(InputList.begin(), InputList.end(), inp) ==
InputList.end()) {
if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER)
fIsAllSplines = false;
InputList.push_back(subsamples[i]->GetInput());
}
}
}
return InputList;
}
std::vector<MeasurementBase*> JointFCN::GetSubSampleList() {
std::vector<MeasurementBase*> SampleList;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
SampleList.push_back(subsamples[i]);
}
}
return SampleList;
}
//***************************************************
void JointFCN::ReconfigureUsingManager() {
- //***************************************************
+//***************************************************
// 'Slow' Event Manager Reconfigure
LOG(REC) << "Event Manager Reconfigure" << std::endl;
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// If we are siving signal, reset all containers.
bool savesignal = (FitPar::Config().GetParB("SignalReconfigures"));
if (savesignal) {
// Reset all of our event signal vectors
fSignalEventBoxes.clear();
fSignalEventFlags.clear();
fSampleSignalFlags.clear();
fSignalEventSplines.clear();
}
// Make sure we have a list of inputs
if (fInputList.empty()) {
fInputList = GetInputList();
fSubSampleList = GetSubSampleList();
}
// If all inputs are splines make sure the readers are told
// they need to be reconfigured.
std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell reader in each BaseEvent it needs a Reconfigure next weight calc.
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) {
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
}
// MAIN INPUT LOOP ====================
int fillcount = 0;
int inputcount = 0;
inp_iter = fInputList.begin();
// Loop over each input in manager
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Get event information
FitEvent* curevent = curinput->FirstNuisanceEvent();
curinput->CreateCache();
int i = 0;
int nevents = curinput->GetNEvents();
int countwidth = nevents / 5;
// Start event loop iterating until we get a NULL pointer.
while (curevent) {
// Get Event Weight
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
double rwweight = curevent->Weight;
// std::cout << "RWWeight = " << curevent->RWWeight << " " <<
// curevent->InputWeight << std::endl;
// Logging
// std::cout << CHECKLOG(1) << std::endl;
if (LOGGING(REC)) {
if (i % countwidth == 0) {
QLOG(REC, curinput->GetName()
- << " : Processed " << i << " events. [M, W] = ["
- << curevent->Mode << ", " << rwweight << "]");
+ << " : Processed " << i << " events. [M, W] = ["
+ << curevent->Mode << ", " << rwweight << "]");
}
}
// Setup flag for if signal found in at least one sample
bool foundsignal = false;
// Create a new signal bitset for this event
std::vector<bool> signalbitset(fSubSampleList.size());
// Create a new signal box vector for this event
std::vector<MeasurementVariableBox*> signalboxes;
// Start measurement iterator
size_t measitercount = 0;
std::vector<MeasurementBase*>::iterator meas_iter =
- fSubSampleList.begin();
+ fSubSampleList.begin();
// Loop over all subsamples (sub in JointMeas)
for (; meas_iter != fSubSampleList.end(); meas_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// Compare input pointers, to current input, skip if not.
// Pointer tells us if it matches without doing ID checks.
if (curinput != curmeas->GetInput()) {
if (savesignal) {
// Set bit to 0 as definitely not signal
signalbitset[measitercount] = 0;
}
// Count up what measurement we are on.
measitercount++;
// Skip sample as input not signal.
continue;
}
// Fill events for matching inputs.
MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);
bool signal = curmeas->isSignal(curevent);
curmeas->SetSignal(signal);
curmeas->FillHistograms(curevent->Weight);
// If its Signal tally up fills
if (signal) {
fillcount++;
}
// If we are saving signal/splines fill the bitset
if (savesignal) {
signalbitset[measitercount] = signal;
}
// If signal save a clone of the event box for use later.
if (savesignal and signal) {
foundsignal = true;
signalboxes.push_back(box->CloneSignalBox());
}
// Keep track of Measurement we are on.
measitercount++;
}
// Once we've filled the measurements, if saving signal
// push back if any sample flagged this event as signal
if (savesignal) {
fSignalEventFlags.push_back(foundsignal);
}
// Save the vector of signal boxes for this event
if (savesignal and foundsignal) {
fSignalEventBoxes.push_back(signalboxes);
fSampleSignalFlags.push_back(signalbitset);
}
// If all inputs are splines we can save the spline coefficients
// for fast in memory reconfigures later.
if (fIsAllSplines and savesignal and foundsignal) {
// Make temp vector to push back with
std::vector<float> coeff;
for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) {
coeff.push_back(curevent->fSplineCoeff[l]);
}
// Push back to signal event splines. Kept in sync with
// fSignalEventBoxes size.
// int splinecount = fSignalEventSplines.size();
fSignalEventSplines.push_back(coeff);
// if (splinecount % 1000 == 0) {
// std::cout << "Pushed Back Coeff " << splinecount << " : ";
// for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++)
// {
// std::cout << " " << fSignalEventSplines[splinecount][l];
// }
// std::cout << std::endl;
// }
}
// Clean up vectors once done with this event
signalboxes.clear();
signalbitset.clear();
// Iterate to the next event.
curevent = curinput->NextNuisanceEvent();
i++;
}
curinput->RemoveCache();
// Keep track of what input we are on.
inputcount++;
}
// End of Event Loop ===============================
// Now event loop is finished loop over all Measurements
// Converting Binned events to XSec Distributions
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Print out statements on approximate memory usage for profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
if (savesignal) {
int mem =
- ( // sizeof(fSignalEventBoxes) +
- // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
- sizeof(MeasurementVariableBox1D) * fillcount) *
- 1E-6;
+ ( // sizeof(fSignalEventBoxes) +
+ // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
+ sizeof(MeasurementVariableBox1D) * fillcount) *
+ 1E-6;
LOG(REC) << " -> Saved " << fillcount
<< " signal boxes for faster access. (~" << mem << " MB)"
<< std::endl;
if (fIsAllSplines and !fSignalEventSplines.empty()) {
int splmem = sizeof(float) * fSignalEventSplines.size() *
fSignalEventSplines[0].size() * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size()
<< " spline sets into memory. (~" << splmem << " MB)"
<< std::endl;
}
}
LOG(REC) << "Time taken ReconfigureUsingManager() : "
<< time(NULL) - timestart << std::endl;
// Check SignalReconfigures works for all samples
- if (savesignal){
+ if (savesignal) {
double likefull = GetLikelihood();
ReconfigureFastUsingManager();
double likefast = GetLikelihood();
-
+
if (fabs(likefull - likefast) > 0.0001)
- {
- ERROR(FTL,"Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast);
- ERROR(FTL,"This means some samples you are using are not setup to use SignalReconfigures=1");
- ERROR(FTL,"Please turn OFF signal reconfigures.");
- throw;
- } else {
+ {
+ ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast);
+ ERROR(FTL, "This means some samples you are using are not setup to use SignalReconfigures=1");
+ ERROR(FTL, "Please turn OFF signal reconfigures.");
+ throw;
+ } else {
LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl;
}
}
// End of reconfigure
return;
};
//***************************************************
void JointFCN::ReconfigureFastUsingManager() {
- //***************************************************
+//***************************************************
LOG(FIT) << " -> Doing FAST using manager" << std::endl;
// Get Start time for profilling
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// Check for saved variables if not do a full reconfigure.
if (fSignalEventFlags.empty()) {
ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl;
ReconfigureUsingManager();
return;
}
bool fFillNuisanceEvent =
- FitPar::Config().GetParB("FullEventOnSignalReconfigure");
+ FitPar::Config().GetParB("FullEventOnSignalReconfigure");
// Setup fast vector iterators.
std::vector<bool>::iterator inpsig_iter = fSignalEventFlags.begin();
std::vector<std::vector<MeasurementVariableBox*> >::iterator box_iter =
- fSignalEventBoxes.begin();
+ fSignalEventBoxes.begin();
std::vector<std::vector<float> >::iterator spline_iter =
- fSignalEventSplines.begin();
+ fSignalEventSplines.begin();
std::vector<std::vector<bool> >::iterator samsig_iter =
- fSampleSignalFlags.begin();
+ fSampleSignalFlags.begin();
int splinecount = 0;
// Setup stuff for logging
int fillcount = 0;
int nevents = fSignalEventFlags.size();
int countwidth = nevents / 20;
// If All Splines tell splines they need a reconfigure.
std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl;
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell each fSplineRead in BaseFitEvent to reconf next weight calc
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead)
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
// Loop over all possible spline inputs
double* coreeventweights = new double[fSignalEventBoxes.size()];
splinecount = 0;
inp_iter = fInputList.begin();
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
// Loop over all signal flags
// For each valid signal flag add one to splinecount
// Get Splines from that count and add to weight
// Add splinecount
int sigcount = 0;
splinecount = 0;
// #pragma omp parallel for shared(splinecount,sigcount)
for (uint iinput = 0; iinput < fInputList.size(); iinput++) {
InputHandlerBase* curinput = fInputList[iinput];
BaseFitEvt* curevent = curinput->FirstBaseEvent();
for (int i = 0; i < curinput->GetNEvents(); i++) {
double rwweight = 0.0;
if (fSignalEventFlags[sigcount]) {
// Get Event Info
if (!fIsAllSplines) {
if (fFillNuisanceEvent)
curinput->GetNuisanceEvent(i);
else
curevent = curinput->GetBaseEvent(i);
} else {
curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0];
}
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
rwweight = curevent->Weight;
coreeventweights[splinecount] = rwweight;
if (splinecount % countwidth == 0) {
LOG(REC) << "Processed " << splinecount
<< " event weights. W = " << rwweight << std::endl;
}
// #pragma omp atomic
splinecount++;
}
// #pragma omp atomic
sigcount++;
}
}
LOG(SAM) << "Processed event weights." << std::endl;
// #pragma omp barrier
// Reset Iterators
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
box_iter = fSignalEventBoxes.begin();
samsig_iter = fSampleSignalFlags.begin();
int nsplineweights = splinecount;
splinecount = 0;
// Start of Fast Event Loop ============================
// Start input iterators
// Loop over number of inputs
for (int ispline = 0; ispline < nsplineweights; ispline++) {
double rwweight = coreeventweights[ispline];
// Get iterators for this event
std::vector<bool>::iterator subsamsig_iter = (*samsig_iter).begin();
std::vector<MeasurementVariableBox*>::iterator subbox_iter =
- (*box_iter).begin();
+ (*box_iter).begin();
// Loop over all sub measurements.
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// If event flagged as signal for this sample fill from the box.
if (*subsamsig_iter) {
curmeas->SetSignal(true);
curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);
// Move onto next box if there is one.
subbox_iter++;
fillcount++;
}
}
if (ispline % countwidth == 0) {
LOG(REC) << "Filled " << ispline << " sample weights." << std::endl;
}
// Iterate over the main signal event containers.
samsig_iter++;
box_iter++;
spline_iter++;
splinecount++;
}
// End of Fast Event Loop ===================
LOG(SAM) << "Filled sample distributions." << std::endl;
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Cleanup coreeventweights
if (fIsAllSplines) {
delete coreeventweights;
}
// Print some reconfigure profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
LOG(REC) << "Time taken ReconfigureFastUsingManager() : "
<< time(NULL) - timestart << std::endl;
}
//***************************************************
void JointFCN::Write() {
- //***************************************************
+//***************************************************
// Save a likelihood/ndof plot
LOG(MIN) << "Writing likelihood plot.." << std::endl;
std::vector<double> likes;
std::vector<double> ndofs;
std::vector<std::string> names;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
- iter++){
+ iter++) {
MeasurementBase* exp = *iter;
double like = exp->GetLikelihood();
double ndof = exp->GetNDOF();
std::string name = exp->GetName();
likes.push_back(like);
ndofs.push_back(ndof);
names.push_back(name);
}
- TH1D likehist = TH1D("likelihood_hist","likelihood_hist",
- likes.size(), 0.0, double(likes.size()));
- TH1D ndofhist = TH1D("ndof_hist","ndof_hist",
- ndofs.size(), 0.0, double(ndofs.size()));
- TH1D divhist = TH1D("likedivndof_hist","likedivndof_hist",
- likes.size(), 0.0, double(likes.size()));
- for (size_t i = 0; i < likehist.GetNbinsX(); i++){
- likehist.SetBinContent(i+1, likes[i]);
- ndofhist.SetBinContent(i+1, ndofs[i]);
- if (ndofs[i] != 0.0){
- divhist.SetBinContent(i+1, likes[i]/ndofs[i]);
+ TH1D likehist = TH1D("likelihood_hist", "likelihood_hist",
+ likes.size(), 0.0, double(likes.size()));
+ TH1D ndofhist = TH1D("ndof_hist", "ndof_hist",
+ ndofs.size(), 0.0, double(ndofs.size()));
+ TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist",
+ likes.size(), 0.0, double(likes.size()));
+ for (size_t i = 0; i < likehist.GetNbinsX(); i++) {
+ likehist.SetBinContent(i + 1, likes[i]);
+ ndofhist.SetBinContent(i + 1, ndofs[i]);
+ if (ndofs[i] != 0.0) {
+ divhist.SetBinContent(i + 1, likes[i] / ndofs[i]);
}
- likehist.GetXaxis()->SetBinLabel(i+1, names[i].c_str());
- ndofhist.GetXaxis()->SetBinLabel(i+1, names[i].c_str());
- divhist.GetXaxis()->SetBinLabel(i+1, names[i].c_str());
+ likehist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
+ ndofhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
+ divhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
}
likehist.Write();
ndofhist.Write();
divhist.Write();
// Loop over individual experiments and call Write
LOG(MIN) << "Writing each of the data classes..." << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->Write();
}
// Save Pull Terms
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Write();
}
if (FitPar::Config().GetParB("EventManager")) {
// Get list of inputs
std::map<int, InputHandlerBase*> fInputs =
- FitBase::EvtManager().GetInputs();
+ FitBase::EvtManager().GetInputs();
std::map<int, InputHandlerBase*>::const_iterator iterInp;
for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) {
InputHandlerBase* input = (iterInp->second);
input->GetFluxHistogram()->Write();
input->GetXSecHistogram()->Write();
input->GetEventHistogram()->Write();
}
}
};
//***************************************************
void JointFCN::SetFakeData(std::string fakeinput) {
- //***************************************************
+//***************************************************
LOG(MIN) << "Setting fake data from " << fakeinput << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->SetFakeDataValues(fakeinput);
}
return;
}
//***************************************************
void JointFCN::ThrowDataToy() {
- //***************************************************
+//***************************************************
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->ThrowDataToy();
}
return;
}
diff --git a/src/FCN/JointFCN.h b/src/FCN/JointFCN.h
index 361a75a..bf87424 100755
--- a/src/FCN/JointFCN.h
+++ b/src/FCN/JointFCN.h
@@ -1,165 +1,176 @@
#ifndef _JOINT_FCN_H_
#define _JOINT_FCN_H_
/*!
* \addtogroup FCN
* @{
*/
#include <iostream>
#include <vector>
#include <fstream>
#include <list>
// ROOT headers
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include <TMatrixDSym.h>
#include "TGraphErrors.h"
#include <TVectorD.h>
#include <TMath.h>
#include "FitUtils.h"
// External fitter headers
#include "MeasurementBase.h"
#include "SampleList.h"
#define GetCurrentDir getcwd
#include "EventManager.h"
#include "Math/Functor.h"
#include "ParamPull.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#include "MeasurementVariableBox.h"
#include "MeasurementVariableBox1D.h"
using namespace FitUtils;
using namespace FitBase;
//! Main FCN Class which ROOT's joint function needs to evaulate the chi2 at each stage of the fit.
class JointFCN
{
public:
//! Constructor
//! cardfile = Path to input card file listing samples
JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile=NULL);
JointFCN(TFile* outfile=NULL); // Loads from global config
//! Destructor
~JointFCN();
//! Create sample list from cardfile
void LoadSamples(std::vector<nuiskey> samplekeys);
void LoadPulls(std::vector<nuiskey> pullkeys);
//! Main Likelihood evaluation FCN
double DoEval(const double *x);
//! Func Wrapper for ROOT
inline double operator() (const std::vector<double> & x) {
double* x_array = new double[x.size()];
return this->DoEval(x_array);
};
//! Func Wrapper for ROOT
inline double operator() (const double *x) {
return this->DoEval(x);
};
//! Create a TTree to save all dial value iterations for this FCN
void CreateIterationTree(std::string name, FitWeight* rw);
//! Fills dial values and sample likelihoods into tree
void FillIterationTree(FitWeight* rw);
//! Writes TTree to fOutput directory
void WriteIterationTree();
//! Deletes TTree
void DestroyIterationTree();
//! Get Degrees of Freedom for samples (NBins)
int GetNDOF();
//! Return NDOF wrapper
inline unsigned int NDim() {return this->GetNDOF();};
//! Reconfigure samples where we force all events to be looped over.
void ReconfigureAllEvents() ;
//! Call Reconfigure on samples.
//! Choice of reconfigure type depends on whether dials have changed
//! and the MC has been filled.
void ReconfigureSamples(bool fullconfig = false);
//! Call reconfigure on only signal events (defaults to all events if CurIter=0)
void ReconfigureSignal();
//! Gets likelihood for all samples in FCN (assuming uncorrelated)
double GetLikelihood();
//! Returns list of pointers to the samples
inline std::list<MeasurementBase*> GetSampleList(){ return fSamples; }
//! Return list of pointers to all the pulls
inline std::list<ParamPull*> GetPullList(){ return fPulls; };
//! Write all samples to output DIR
void Write();
//! Set Fake data from file/MC
void SetFakeData(std::string fakeinput);
//! Reconfigure looping over duplicate inputs
void ReconfigureUsingManager();
//! Reconfigure Fast looping over duplicate inputs
void ReconfigureFastUsingManager();
/// Throws data according to current stats
void ThrowDataToy();
std::vector<MeasurementBase*> GetSubSampleList();
std::vector<InputHandlerBase*> GetInputList();
private:
//! Append the experiments to include in the fit to this list
std::list<MeasurementBase*> fSamples;
//! Append parameter pull terms to include penalties in the fit to this list
std::list<ParamPull*> fPulls;
TDirectory *fOutputDir; //!< Directory to save contents
std::string fCardFile; //!< Input Card text file
bool fDialChanged; //!< Flag for if RW dials changed
UInt_t fCurIter; //!< Counter for how many times reconfigure called
bool fMCFilled; //!< Check MC has at least been filled once
- TTree* fIterationTree; //!< Tree to save RW values on each function call
+ bool fIterationTree; //!< Tree to save RW values on each function call
int fNDials; //!< Number of RW Dials in FitWeight
double* fDialVals; //!< Current Values of RW Dials
double fLikelihood; //!< Current likelihood for joint sample likelihood
double fNDOF; //!< Total NDOF
double* fSampleLikes; //!< Likelihoods for each individual measurement in list
int * fSampleNDOF; //!< NDOF for each individual measurement in list
bool fUsingEventManager; //!< Flag for doing joint comparisons
std::vector< std::vector<float> > fSignalEventSplines;
std::vector< std::vector<MeasurementVariableBox*> > fSignalEventBoxes;
std::vector< bool > fSignalEventFlags;
std::vector< std::vector<bool> > fSampleSignalFlags;
std::vector<InputHandlerBase*> fInputList;
std::vector<MeasurementBase*> fSubSampleList;
bool fIsAllSplines;
+
+ std::vector< int > fIterationCount;
+ std::vector< double > fCurrentValues;
+ std::vector< std::string > fNameValues;
+ std::vector< std::vector<double> > fIterationValues;
+ int fSampleN;
+ std::string fIterationTreeName;
+
+
+
+
};
/*! @} */
#endif // _JOINT_FCN_H_
diff --git a/src/FitBase/ParamPull.cxx b/src/FitBase/ParamPull.cxx
index 4a0c6bc..8f9cb6e 100644
--- a/src/FitBase/ParamPull.cxx
+++ b/src/FitBase/ParamPull.cxx
@@ -1,810 +1,844 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "ParamPull.h"
//*******************************************************************************
ParamPull::ParamPull(std::string name, std::string inputfile, std::string type, std::string dials) {
//*******************************************************************************
fMinHist = NULL;
fMaxHist = NULL;
fTypeHist = NULL;
fDialSelection = dials;
+ fLimitHist = NULL;
fName = name;
fInput = inputfile;
fType = type;
// Set the pull type
SetType(fType);
std::cout << fType << std::endl;
// Setup Histograms from input file
SetupHistograms(fInput);
};
//*******************************************************************************
void ParamPull::SetType(std::string type) {
//*******************************************************************************
fType = type;
// Assume Default if empty
if (type.empty() || type == "DEFAULT") {
ERR(WRN) << "No type specified for ParmPull class " << fName << std::endl;
ERR(WRN) << "Assuming GAUSTHROW/GAUSPULL" << std::endl;
type = "GAUSTHROW/GAUSPULL";
}
// Set Dial options
if (type.find("FRAC") != std::string::npos) {
fDialOptions = "FRAC";
fPlotTitles = ";; Fractional RW Value";
} else if (type.find("ABS") != std::string::npos) {
fDialOptions = "ABS";
fPlotTitles = ";; ABS RW Value";
} else {
fDialOptions = "";
fPlotTitles = ";; RW Value";
}
// Parse throw types
if (type.find("GAUSPULL") != std::string::npos) fCalcType = kGausPull;
else fCalcType = kNoPull;
if (type.find("GAUSTHROW") != std::string::npos) fThrowType = kGausThrow;
+ else if (type.find("FLATTHROW") != std::string::npos) fThrowType = kFlatThrow;
else fThrowType = kNoThrow;
// Extra check to see if throws or pulls are turned off
if (type.find("NOPULL") != std::string::npos) fCalcType = kNoPull;
if (type.find("NOTHROW") != std::string::npos) fThrowType = kNoThrow;
}
//*******************************************************************************
void ParamPull::SetupHistograms(std::string input) {
//*******************************************************************************
// Extract Types from Input
fFileType = "";
const int nfiletypes = 4;
const std::string filetypes[nfiletypes] = {"FIT", "ROOT", "TXT", "DIAL"};
for (int i = 0; i < nfiletypes; i++) {
std::string tempTypes = filetypes[i] + ":";
if (input.find(tempTypes) != std::string::npos) {
fFileType = filetypes[i];
input.replace(input.find(tempTypes), tempTypes.size(), "");
break;
}
}
// Read Files
if (!fFileType.compare("FIT")) ReadFitFile(input);
else if (!fFileType.compare("ROOT")) ReadRootFile(input);
else if (!fFileType.compare("VECT")) ReadVectFile(input);
else if (!fFileType.compare("DIAL")) ReadDialInput(input);
else {
ERR(FTL) << "Unknown ParamPull Type: " << input << std::endl;
throw;
}
// Check Dials are all good
CheckDialsValid();
// Setup MC Histogram
fMCHist = (TH1D*) fDataHist->Clone();
fMCHist->Reset();
fMCHist->SetNameTitle( (fName + "_MC").c_str(),
(fName + " MC" + fPlotTitles).c_str() );
// If no Covar input make an uncorrelated one
if (!fCovar) {
fCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist, 1.0);
}
// If no types or limits are provided give them a default option
if (!fMinHist) {
fMinHist = (TH1D*) fDataHist->Clone();
fMinHist->SetNameTitle( (fName + "_min").c_str(),
(fName + " min" + fPlotTitles).c_str() );
for (int i = 0; i < fMinHist->GetNbinsX(); i++) {
// TODO (P.Stowell) Change this to a NULL system where limits are actually free!
fMinHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
if (!fMaxHist) {
fMaxHist = (TH1D*) fDataHist->Clone();
fMaxHist->SetNameTitle( (fName + "_min").c_str(),
(fName + " min" + fPlotTitles).c_str() );
for (int i = 0; i < fMaxHist->GetNbinsX(); i++) {
fMaxHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
// Set types from state, or to unknown
if (!fTypeHist) {
int deftype = -1;
if (fType.find("T2K") != std::string::npos) { deftype = kT2K; }
else if (fType.find("NEUT") != std::string::npos) { deftype = kNEUT; }
else if (fType.find("NIWG") != std::string::npos) { deftype = kNIWG; }
else if (fType.find("GENIE") != std::string::npos) { deftype = kGENIE; }
else if (fType.find("NORM") != std::string::npos) { deftype = kNORM; }
else if (fType.find("NUWRO") != std::string::npos) { deftype = kNUWRO; }
fTypeHist = new TH1I( (fName + "_type").c_str(),
(fName + " type" + fPlotTitles).c_str(),
fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX() );
for (int i = 0; i < fTypeHist->GetNbinsX(); i++) {
fTypeHist->SetBinContent(i + 1, deftype );
}
}
// Sort Covariances
fInvCovar = StatUtils::GetInvert(fCovar);
fDecomp = StatUtils::GetDecomp(fCovar);
// Create DataTrue for Throws
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
fDataOrig = (TH1D*) fDataHist->Clone();
fDataOrig->SetNameTitle( (fName + "_origdata").c_str(),
(fName + " origdata" + fPlotTitles).c_str() );
// Select only dials we want
if (!fDialSelection.empty()) {
(*fDataHist) = RemoveBinsNotInString(*fDataHist, fDialSelection);
}
}
//*******************************************************************************
TH1D ParamPull::RemoveBinsNotInString(TH1D hist, std::string mystr) {
//*******************************************************************************
// Make list of allowed bins
std::vector<std::string> allowedbins;
for (int i = 0; i < hist.GetNbinsX(); i++) {
std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1D newhist = TH1D( hist.GetName(), hist.GetTitle(),
(Int_t)nbins, 0.0, (Double_t)nbins);
// Setup bins
for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
//*******************************************************************************
TH1I ParamPull::RemoveBinsNotInString(TH1I hist, std::string mystr) {
//*******************************************************************************
// Make list of allowed bins
std::vector<std::string> allowedbins;
for (int i = 0; i < hist.GetNbinsX(); i++) {
std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1I newhist = TH1I( hist.GetName(), hist.GetTitle(),
(Int_t)nbins, 0.0, (Int_t)nbins);
// Setup bins
for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
//*******************************************************************************
void ParamPull::ReadFitFile(std::string input) {
//*******************************************************************************
TFile* tempfile = new TFile(input.c_str(), "READ");
// Read Data
fDataHist = (TH1D*) tempfile->Get("fit_dials_free");
if (!fDataHist) {
ERR(FTL) << "Can't find TH1D hist fit_dials in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
// Read Covar
TH2D* tempcov = (TH2D*) tempfile->Get("covariance_free");
if (!tempcov) {
ERR(FTL) << "Can't find TH2D covariance_free in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
return;
}
//*******************************************************************************
void ParamPull::ReadRootFile(std::string input) {
//*******************************************************************************
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input, ";");
// Check all given
if (inputlist.size() < 2) {
ERR(FTL) << "Covar supplied in 'ROOT' format should have 3 semi-colon seperated entries!" << std::endl
<< "ROOT:filename;histname[;covarname]" << std::endl;
ERR(FTL) << "histname = TH1D, covarname = TH2D" << std::endl;
throw;
}
// Get Entries
std::string filename = inputlist[0];
LOG(DEB) << filename << std::endl;
std::string histname = inputlist[1];
LOG(DEB) << histname << std::endl;
LOG(DEB) << input << std::endl;
// Read File
TFile* tempfile = new TFile(filename.c_str(), "READ");
if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << filename << std::endl;
throw;
}
// Read Hist
fDataHist = (TH1D*) tempfile->Get(histname.c_str());
if (!fDataHist) {
ERR(FTL) << "Can't find TH1D hist " << histname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
LOG(DEB) << "READING COVAR" << std::endl;
// Read Covar
if (inputlist.size() > 2) {
std::string covarname = inputlist[2];
LOG(DEB) << "COVARNAME = " << covarname << std::endl;
TH2D* tempcov = (TH2D*) tempfile->Get(covarname.c_str());
if (!tempcov) {
ERR(FTL) << "Can't find TH2D covar " << covarname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
// Uncorrelated
} else {
LOG(SAM) << "No Covar provided so using diagonal errors for "
<< fName << std::endl;
fCovar = NULL;
}
}
//*******************************************************************************
void ParamPull::ReadVectFile(std::string input) {
//*******************************************************************************
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input, ";");
if (inputlist.size() < 4) {
ERR(FTL) << "Need 3 inputs for vector input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
// Open File
std::string rootname = inputlist[0];
TFile* tempfile = new TFile(rootname.c_str(), "READ");
if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << rootname << std::endl;
throw;
}
// Get Name
std::string tagname = inputlist[1];
// TVector<std::string> dialtags = tempfile->Get(tagname.c_str());
// if (!dialtags){
// ERR(FTL) << "Can't find list of dial names!" << std::endl;
// }
// Get Values
std::string valuename = inputlist[2];
TVectorD* dialvals = (TVectorD*)tempfile->Get(valuename.c_str());
if (!dialvals) {
ERR(FTL) << "Can't find dial values" << std::endl;
}
// Get Matrix
std::string matrixname = inputlist[3];
TMatrixD* matrixvals = (TMatrixD*)tempfile->Get(matrixname.c_str());
if (!matrixvals) {
ERR(FTL) << "Can't find matirx values" << std::endl;
}
// Get Types
if (inputlist.size() > 4) {
std::string typesname = inputlist[3];
}
// Get Minimum
if (inputlist.size() > 5) {
std::string minname = inputlist[4];
}
// Get Maximum
if (inputlist.size() > 6) {
std::string maxname = inputlist[5];
}
}
//*******************************************************************************
void ParamPull::ReadDialInput(std::string input) {
//*******************************************************************************
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(input, ";");
if (inputlist.size() < 3) {
ERR(FTL) << "Need 3 inputs for dial input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
std::vector<double> inputvals = GeneralUtils::ParseToDbl(input, ";");
std::string dialname = inputlist[0];
double val = inputvals[1];
double err = inputvals[2];
fDataHist = new TH1D( (fName + "_data").c_str(),
(fName + "_data" + fPlotTitles).c_str(), 1, 0, 1);
fDataHist->SetBinContent(1, val);
fDataHist->SetBinError(1, err);
fDataHist->GetXaxis()->SetBinLabel(1, dialname.c_str());
+
+ fLimitHist = new TH1D( (fName + "_limits").c_str(),
+ (fName + "_limits" + fPlotTitles).c_str(), 1, 0, 1);
+ fLimitHist->Reset();
+ if (inputvals.size() > 4) {
+ fLimitHist->SetBinContent(1, (inputvals[3] + inputvals[4]) / 2.0);
+ fLimitHist->SetBinError(1, (inputvals[4] - inputvals[3]) / 2.0);
+ }
+
fCovar = NULL;
}
//*******************************************************************************
std::map<std::string, int> ParamPull::GetAllDials() {
//*******************************************************************************
std::map<std::string, int> dialtypemap;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
std::string name = fDataHist->GetXaxis()->GetBinLabel(i + 1);
int type = fTypeHist->GetBinContent(i + 1);
dialtypemap[name] = type;
}
return dialtypemap;
}
//*******************************************************************************
bool ParamPull::CheckDialsValid() {
//*******************************************************************************
return true;
std::string helpstring = "";
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
std::string name = std::string(fDataHist->GetXaxis()->GetBinLabel(i + 1));
// If dial exists its all good
if (FitBase::GetRW()->DialIncluded(name)) continue;
// If it doesn't but its a sample norm also continue
if (name.find("_norm") != std::string::npos) {
ERR(WRN) << "Norm dial included in covar but not set in FitWeight." << std::endl;
ERR(WRN) << "Assuming its a sample norm and skipping..." << std::endl;
}
// Dial unknown so print a help statement
std::ostringstream tempstr;
tempstr << "unknown_parameter " << name << " "
<< fDataHist->GetBinContent(i + 1) << " "
<< fDataHist->GetBinContent(i + 1) - fDataHist->GetBinError(i + 1) << " "
<< fDataHist->GetBinContent(i + 1) + fDataHist->GetBinError(i + 1) << " "
<< fDataHist->GetBinError(i + 1) << " ";
if (!fType.empty()) tempstr << fType << std::endl;
else tempstr << "FREE" << std::endl;
helpstring += tempstr.str();
}
// Show statement before failing
if (!helpstring.empty()) {
ERR(WRN) << "Dial(s) included in covar but not set in FitWeight." << std::endl
<< "ParamPulls needs to know how you want it to be treated." << std::endl
<< "Include the following lines into your card:" << std::endl;
ERR(WRN) << helpstring << std::endl;
throw;
return false;
} else {
return true;
}
}
//*******************************************************************************
void ParamPull::Reconfigure() {
//*******************************************************************************
FitWeight* rw = FitBase::GetRW();
// Get Dial Names that are valid
std::vector<std::string> namevec = rw->GetDialNames();
std::vector<double> valuevec = rw->GetDialValues();
// Set Bin Values from RW
for (UInt_t i = 0; i < namevec.size(); i++) {
// Loop over bins and check name matches
std::string syst = namevec.at(i);
double systval = valuevec.at(i);
std::vector<std::string> allsyst = GeneralUtils::ParseToStr(syst, ",");
// Proper Reconf using RW
for (int j = 0; j < fMCHist->GetNbinsX(); j++) {
// Search for the name of this bin in the corrent dial
std::string binname = std::string(fMCHist->GetXaxis()->GetBinLabel(j + 1) );
// Check Full Name
if (!syst.compare(binname.c_str())) {
fMCHist->SetBinContent(j + 1, systval);
break;
}
std::vector<std::string> splitbinname = GeneralUtils::ParseToStr(binname, ",");
for (size_t l = 0; l < splitbinname.size(); l++) {
std::string singlebinname = splitbinname[l];
for (size_t k = 0; k < allsyst.size(); k++) {
if (!allsyst[k].compare(singlebinname.c_str())) {
fMCHist->SetBinContent(j + 1, systval);
}
}
}
}
}
return;
};
//*******************************************************************************
void ParamPull::ResetToy(void) {
//*******************************************************************************
if (fDataHist) delete fDataHist;
LOG(DEB) << "Resetting toy" << std::endl;
LOG(DEB) << fDataTrue << std::endl;
fDataHist = (TH1D*)fDataTrue->Clone();
LOG(DEB) << "Setting name" << std::endl;
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
}
//*******************************************************************************
void ParamPull::SetFakeData(std::string fakeinput) {
//*******************************************************************************
// Set from MC Setting
if (!fakeinput.compare("MC")) {
// Copy MC into data
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fMCHist->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " fakedata" + fPlotTitles).c_str() );
// Copy original data errors
for (int i = 0; i < fDataOrig->GetNbinsX(); i++) {
fDataHist->SetBinError(i + 1, fDataOrig->GetBinError(i + 1) );
}
// Make True Toy Central Value Hist
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
} else {
ERR(FTL) << "Trying to set fake data for ParamPulls not from MC!" << std::endl;
ERR(FTL) << "Not currently implemented.." << std::endl;
throw;
}
}
//*******************************************************************************
void ParamPull::RemoveFakeData() {
//*******************************************************************************
delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
}
//*******************************************************************************
double ParamPull::GetLikelihood() {
//*******************************************************************************
double like = 0.0;
switch (fCalcType) {
// Gaussian Calculation with correlations
case kGausPull:
like = StatUtils::GetChi2FromCov(fDataHist, fMCHist, fInvCovar, NULL);
like *= 1E-76;
break;
// Default says this has no pull
case kNoThrow:
default:
like = 0.0;
break;
}
LOG(DEB) << "Likelihood = " << like << " " << fCalcType << std::endl;
return like;
};
//*******************************************************************************
int ParamPull::GetNDOF() {
//*******************************************************************************
int ndof = 0;
if (fCalcType != kNoThrow) {
- ndof += fDataHist->GetNbinsX();
+ ndof = fDataHist->GetNbinsX();
}
return ndof;
};
//*******************************************************************************
void ParamPull::ThrowCovariance() {
//*******************************************************************************
// Reset toy for throw
ResetToy();
- LOG(DEB) << "Toy Reset " << std::endl;
+ LOG(FIT) << "Creating new toy dataset" << std::endl;
// Generate random Gaussian throws
std::vector<double> randthrows;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
double randtemp = 0.0;
switch (fThrowType) {
// Gaussian Throws
case kGausThrow:
randtemp = gRandom->Gaus(0.0, 1.0);
break;
+ // Uniform Throws
+ case kFlatThrow:
+ randtemp = gRandom->Uniform(0.0, 1.0);
+ if (fLimitHist) {
+ randtemp = fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) * ( randtemp * 2 - 1 );
+ }
+ break;
+
// No Throws (DEFAULT)
default:
break;
}
randthrows.push_back(randtemp);
}
// Create Bin Modifications
double totalres = 0.0;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
// Calc Bin Mod
double binmod = 0.0;
- for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
- // std::cout << "DECOMP " << j << " " << i << " " << randthrows.at(j) << std::endl;
- binmod += (*fDecomp)(j, i) * randthrows.at(j);
+
+ if (fThrowType == kGausThrow) {
+ for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
+ binmod += (*fDecomp)(j, i) * randthrows.at(j);
+ }
+ } else if (fThrowType == kFlatThrow) {
+ binmod = randthrows.at(i) - fDataHist->GetBinContent(i + 1);
}
+
// Add up fraction dif
totalres += binmod;
// Add to current data
fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) + binmod);
}
// Rename
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " toydata" + fPlotTitles).c_str() );
- // Print Status
- LOG(REC) << "Created new toy histogram. Total Dif = "
- << totalres << std::endl;
+ // Check Limits
+ if (fLimitHist) {
+ for (int i = 0; i < fLimitHist->GetNbinsX(); i++) {
+ if (fLimitHist->GetBinError(i + 1) == 0.0) continue;
+ if (fDataHist->GetBinContent(i + 1) > fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) ||
+ fDataHist->GetBinContent(i + 1) < fLimitHist->GetBinContent(i + 1) - fLimitHist->GetBinError(i + 1)) {
+ this->ThrowCovariance();
+ }
+ }
+ }
+
return;
};
//*******************************************************************************
TH2D ParamPull::GetCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fInvCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_INVCOV").c_str(),
(fName + " InvertedCovariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
TH2D ParamPull::GetFullCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_COV").c_str(),
(fName + " Covariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
TH2D ParamPull::GetDecompCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_DEC").c_str(),
(fName + " Decomposition;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
void ParamPull::Write(std::string writeoptt) {
//*******************************************************************************
fDataHist->Write();
fMCHist->Write();
-
+ if (fLimitHist) {
+ fLimitHist->Write();
+ }
GetCovar().Write();
GetFullCovar().Write();
GetDecompCovar().Write();
return;
};
diff --git a/src/FitBase/ParamPull.h b/src/FitBase/ParamPull.h
index 22e707d..7da9eac 100644
--- a/src/FitBase/ParamPull.h
+++ b/src/FitBase/ParamPull.h
@@ -1,188 +1,193 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef PARAM_PULL_H_SEEN
#define PARAM_PULL_H_SEEN
/*!
* \addtogroup FitBase
* @{
*/
#include <stdlib.h>
#include <numeric>
#include <math.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <list>
#include <vector>
// ROOT includes
#include <TROOT.h>
#include <TH1D.h>
#include <TTree.h>
#include <TFile.h>
#include <TMatrixDSym.h>
#include <TDecompSVD.h>
#include <TVectorD.h>
// Fit Includes
#include "PlotUtils.h"
#include "StatUtils.h"
#include "FitWeight.h"
#include "FitLogger.h"
#include "EventManager.h"
#include "TVector.h"
using namespace std;
//! Enums to allow Non Guassian Pulls in Future
enum FitPullTypes {
kUnknownPull = -1,
kNoPull = 0,
kGausPull = 1
};
enum FitThrowTypes {
kUnknownThrow = -1,
kNoThrow = 0,
- kGausThrow = 1
+ kGausThrow = 1,
+ kFlatThrow = 2
};
//! Used to produce gaussian penalty terms in the fit.
class ParamPull {
public:
//! Default Constructor
ParamPull(std::string name, std::string inputfile, std::string type, std::string dials="");
//! Default destructor
virtual ~ParamPull(void) {};
// Set dial types (DEFAULT,ABS,FRAC)
void SetType(std::string type);
// Setup Histogram inputs (from previous fit file, or ROOT file)
void SetupHistograms(std::string input);
// Read a previous NUISANCE file
void ReadFitFile(std::string input);
// Read a ROOT file with any histograms in
void ReadRootFile(std::string input);
// Read Text file
void ReadVectFile(std::string input);
// Read a single dial central value and error
void ReadDialInput(std::string input);
//! Reconfigure function reads in current RW engine dials and sets their value to MC
void Reconfigure(void);
//! Get likelihood given the current values
double GetLikelihood(void);
//! Get NDOF if used in likelihoods
int GetNDOF(void);
// Get Covariance Matrices as plots
TH2D GetCovar(void);
TH2D GetFullCovar(void);
TH2D GetDecompCovar(void);
// Get Covariance Matrices
inline TMatrixDSym GetCovarMatrix (void) const { return *fInvCovar; };
inline TMatrixDSym GetFullCovarMatrix (void) const { return *fCovar; };
inline TMatrixDSym GetDecompCovarMatrix (void) const { return *fDecomp; };
//! Save the histograms
void Write(std::string writeopt="");
//! Throw the dial values using the current covariance. Useful for parameter throws.
void ThrowCovariance(void);
//! Compare dials to RW
bool CheckDialsValid(void);
//! Reset toy data back to original data
void ResetToy(void);
//! Read fake data from MC
void SetFakeData(std::string fakeinput);
//! Reset fake data back to original data (before fake data or throws)
void RemoveFakeData();
// Get Functions
inline std::string GetName (void) const { return fName; };
inline std::string GetInput (void) const { return fInput; };
inline std::string GetType (void) const { return fType; };
inline std::string GetFileType (void) const { return fFileType; };
inline std::string GetDialOptions (void) const { return fDialOptions; };
std::map<std::string, int> GetAllDials();
inline TH1D GetDataHist (void) const { return *fDataHist; };
inline TH1D GetDataTrue (void) const { return *fDataTrue; };
inline TH1D GetDataOrig (void) const { return *fDataOrig; };
inline TH1D GetMCHist (void) const { return *fMCHist; };
inline TH1D GetMaxHist (void) const { return *fMaxHist; };
inline TH1D GetMinHist (void) const { return *fMinHist; };
inline TH1I GetDialTypes (void) const { return *fTypeHist; };
+ inline TH1D GetLimitHist (void) const { return *fLimitHist; };
+
private:
TH1D RemoveBinsNotInString(TH1D hist, std::string mystr);
TH1I RemoveBinsNotInString(TH1I hist, std::string mystr);
std::string fName; //!< Pull Name
std::string fInput; //!< Pull input string
std::string fType; //!< Pull options type
std::string fFileType; //!< Pull input file types
std::string fPlotTitles; //! Axis format
std::string fDialOptions; //!< Dial handling options
std::string fDialSelection; //!< Dial Selection
TH1D* fMCHist; //!< Current MC Histogram
TH1D* fDataHist; //!< Current data Histogram
TH1D* fDataTrue; //!< True Data (before histogram throws)
TH1D* fDataOrig; //!< Orig Data (without toys or fake data)
TH1D* fMaxHist; //!< Maximum limit on MC/Data
TH1D* fMinHist; //!< Maximum limit on MC/Data
TH1I* fTypeHist; //!< Dial Types
int fCalcType; //!< Method to calculate likelihood
int fThrowType; //!< Method to calculate throws
TMatrixDSym* fCovar; //!< Covariance
TMatrixDSym* fInvCovar; //!< Inverted Covariance
TMatrixDSym* fDecomp; //!< Decomposition
+
+ TH1D* fLimitHist;
};
// Class TypeDefs
typedef std::list<ParamPull*>::const_iterator PullListConstIter;
typedef std::list<ParamPull*>::iterator PullListIter;
typedef std::vector<ParamPull*>::const_iterator PullVectConstIter;
typedef std::vector<ParamPull*>::iterator PullVectIter;
/*! @} */
#endif
diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx
index 4c846a5..d33b35e 100755
--- a/src/Routines/SystematicRoutines.cxx
+++ b/src/Routines/SystematicRoutines.cxx
@@ -1,1516 +1,1517 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "SystematicRoutines.h"
void SystematicRoutines::Init(){
fInputFile = "";
fInputRootFile = NULL;
fOutputFile = "";
fOutputRootFile = NULL;
fCovar = fCovarFree = NULL;
fCorrel = fCorrelFree = NULL;
fDecomp = fDecompFree = NULL;
fStrategy = "ErrorBands";
fRoutines.clear();
fRoutines.push_back("ErrorBands");
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fAllowedRoutines = ("ErrorBands,PlotLimits");
};
SystematicRoutines::~SystematicRoutines(){
};
SystematicRoutines::SystematicRoutines(int argc, char* argv[]){
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
fNThrows = 250;
fStartThrows = 0;
fThrowString = "";
// Make easier to handle arguments.
std::vector<std::string> args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-s", fStartThrows, false, false);
ParserUtils::ParseArgument(args, "-t", fNThrows, false, false);
ParserUtils::ParseArgument(args, "-p", fThrowString, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
SETVERBOSITY(verbocount);
// Proper Setup
if (fStrategy.find("ErrorBands") != std::string::npos ||
fStrategy.find("MergeErrors") != std::string::npos){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
}
// fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupSystematicsFromXML();
SetupCovariance();
SetupRWEngine();
SetupFCN();
GetCovarFromFCN();
// Run();
return;
};
void SystematicRoutines::SetupSystematicsFromXML(){
LOG(FIT) << "Setting up nuismin" << std::endl;
// Setup Parameters ------------------------------------------
std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
// Override if state not given
if (!key.Has("state")){
key.SetS("state","FIX");
}
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Run Parameter Conversion if needed
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);;
fStartVals[parname] = parnom;
fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = parlow;
fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
}
// Setup Samples ----------------------------------------------
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read sample info " << i << " : "
<< samplename << std::endl
<< "\t\t input -> " << samplefile << std::endl
<< "\t\t state -> " << sampletype << std::endl
<< "\t\t norm -> " << samplenorm << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find(normname) != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fCurVals[normname] = samplenorm;
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
}
// Setup Fake Parameters -----------------------------
std::vector<nuiskey> fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nom")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nom");
// Push into vectors
fFakeVals[parname] = parnom;
}
}
void SystematicRoutines::ReadCard(std::string cardfile){
// Read cardlines into vector
std::vector<std::string> cardlines = GeneralUtils::ParseFileToStr(cardfile,"\n");
FitPar::Config().cardLines = cardlines;
// Read Samples first (norm params can be overridden)
int linecount = 0;
for (std::vector<std::string>::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Read Valid Samples
int samstatus = ReadSamples(line);
// Show line if bad to help user
if (samstatus == kErrorStatus) {
ERR(FTL) << "Bad Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
// Read Parameters second
linecount = 0;
for (std::vector<std::string>::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Try Parameter Reads
int parstatus = ReadParameters(line);
int fakstatus = ReadFakeDataPars(line);
// Show line if bad to help user
if (parstatus == kErrorStatus ||
fakstatus == kErrorStatus ){
ERR(FTL) << "Bad Parameter Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
return;
};
int SystematicRoutines::ReadParameters(std::string parstring){
std::string inputspec = "RW Dial Inputs Syntax \n"
"free input w/ limits: TYPE NAME START MIN MAX STEP [STATE] \n"
"fix input: TYPE NAME VALUE [STATE] \n"
"free input w/o limits: TYPE NAME START FREE,[STATE] \n"
"Allowed Types: \n"
"neut_parameter,niwg_parameter,t2k_parameter,"
"nuwro_parameter,gibuu_parameter";
// Check sample input
if (parstring.find("parameter") == std::string::npos) return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0].find("parameter") == std::string::npos){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Input rw dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string partype = strvct[0];
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
double minval = parval - 1.0;
double maxval = parval + 1.0;
double stepval = 1.0;
std::string state = "FIX"; //[DEFAULT]
// Check Type
if (FitBase::ConvDialType(partype) == kUNKNOWN){
ERR(FTL) << "Unknown parameter type! " << partype << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Check Parameter Name
if (FitBase::GetDialEnum(partype, parname) == -1){
ERR(FTL) << "Bad RW parameter name! " << partype << " " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (No Limits)
if (strvct.size() == 4){
state = strvct[3];
}
// Check for weirder inputs
if (strvct.size() > 4 && strvct.size() < 6){
ERR(FTL) << "Provided incomplete limits for " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (With limits and steps)
if (strvct.size() >= 6){
minval = GeneralUtils::StrToDbl(strvct[3]);
maxval = GeneralUtils::StrToDbl(strvct[4]);
stepval = GeneralUtils::StrToDbl(strvct[5]);
}
// Option Extra (dial state after limits)
if (strvct.size() == 7){
state = strvct[6];
}
// Run Parameter Conversion if needed
if (state.find("ABS") != std::string::npos){
parval = FitBase::RWAbsToSigma( partype, parname, parval );
minval = FitBase::RWAbsToSigma( partype, parname, minval );
maxval = FitBase::RWAbsToSigma( partype, parname, maxval );
stepval = FitBase::RWAbsToSigma( partype, parname, stepval );
} else if (state.find("FRAC") != std::string::npos){
parval = FitBase::RWFracToSigma( partype, parname, parval );
minval = FitBase::RWFracToSigma( partype, parname, minval );
maxval = FitBase::RWFracToSigma( partype, parname, maxval );
stepval = FitBase::RWFracToSigma( partype, parname, stepval );
}
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), parname) != fParams.end()){
ERR(FTL) << "Duplicate parameter names given for " << parname << std::endl;
throw;
}
// Setup Containers
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);
fStartVals[parname] = parval;
fCurVals[parname] = fStartVals[parname];
fErrorVals[parname] = 0.0;
fStateVals[parname] = state;
bool fixstate = state.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = minval;
fMaxVals[parname] = maxval;
fStepVals[parname] = stepval;
// Print the parameter
LOG(MIN) << "Read Parameter " << parname << " " << parval << " "
<< minval << " " << maxval << " "
<< stepval << " " << state << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//*******************************************
int SystematicRoutines::ReadFakeDataPars(std::string parstring){
//******************************************
std::string inputspec = "Fake Data Dial Inputs Syntax \n"
"fake value: fake_parameter NAME VALUE \n"
"Name should match dialnames given in actual dial specification.";
// Check sample input
if (parstring.find("fake_parameter") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] == "fake_parameter"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Fake dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Read Inputs
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
// Setup Container
fFakeVals[parname] = parval;
// Print the fake parameter
LOG(MIN) << "Read Fake Parameter " << parname << " " << parval << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//******************************************
int SystematicRoutines::ReadSamples(std::string samstring){
//******************************************
const static std::string inputspec =
"\tsample <sample_name> <input_type>:inputfile.root [OPTS] "
"[norm]\nsample_name: Name "
"of sample to include. e.g. MiniBooNE_CCQE_XSec_1DQ2_nu\ninput_type: The "
"input event format. e.g. NEUT, GENIE, EVSPLN, ...\nOPTS: Additional, "
"optional sample options.\nnorm: Additional, optional sample "
"normalisation factor.";
// Check sample input
if (samstring.find("sample") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(samstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] != "sample"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Sample need to provide at least 3 inputs." << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string samname = strvct[1];
std::string samfile = strvct[2];
if (samfile == "FIX") {
ERR(FTL) << "Input filename was \"FIX\", this line is probably malformed "
"in the input card file. Line:\'"
<< samstring << "\'" << std::endl;
ERR(FTL) << "Expect sample lines to look like:\n\t" << inputspec
<< std::endl;
throw;
}
std::string samtype = "DEFAULT";
double samnorm = 1.0;
// Optional Type
if (strvct.size() > 3) {
samtype = strvct[3];
// samname += "_"+samtype;
// Also get rid of the / and replace it with underscore because it might not be supported character
// while (samname.find("/") != std::string::npos) {
// samname.replace(samname.find("/"), 1, std::string("_"));
// }
}
// Optional Norm
if (strvct.size() > 4) samnorm = GeneralUtils::StrToDbl(strvct[4]);
// Add Sample Names as Norm Dials
std::string normname = samname + "_norm";
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), normname) != fParams.end()){
ERR(FTL) << "Duplicate samples given for " << samname << std::endl;
throw;
}
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStartVals[normname] = samnorm;
fCurVals[normname] = fStartVals[normname];
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = samtype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
// Print read in
LOG(MIN) << "Read sample " << samname << " "
<< samfile << " " << samtype << " "
<< samnorm << std::endl;
// Tell reader its all good
return kGoodStatus;
}
/*
Setup Functions
*/
//*************************************
void SystematicRoutines::SetupRWEngine(){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
}
UpdateRWEngine(fStartVals);
return;
}
//*************************************
void SystematicRoutines::SetupFCN(){
//*************************************
LOG(FIT)<<"Making the jointFCN"<<std::endl;
if (fSampleFCN) delete fSampleFCN;
fSampleFCN = new JointFCN(fOutputRootFile);
SetFakeData();
return;
}
//*************************************
void SystematicRoutines::SetFakeData(){
//*************************************
if (fFakeDataInput.empty()) return;
if (fFakeDataInput.compare("MC") == 0){
LOG(FIT)<<"Setting fake data from MC starting prediction." <<std::endl;
UpdateRWEngine(fFakeVals);
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->SetFakeData("MC");
UpdateRWEngine(fCurVals);
LOG(FIT)<<"Set all data to fake MC predictions."<<std::endl;
} else {
fSampleFCN->SetFakeData(fFakeDataInput);
}
return;
}
//*****************************************
void SystematicRoutines::GetCovarFromFCN(){
//*****************************************
LOG(FIT) << "Loading ParamPull objects from FCN to build covar" << std::endl;
// Make helperstring
std::ostringstream helperstr;
// Keep track of what is being thrown
std::map<std::string, std::string> dialthrowhandle;
// Get Covariance Objects from FCN
std::list<ParamPull*> inputpulls = fSampleFCN->GetPullList();
for (PullListConstIter iter = inputpulls.begin();
iter != inputpulls.end(); iter++){
ParamPull* pull = (*iter);
if (pull->GetType().find("THROW")){
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
LOG(FIT) << "Read ParamPull: " << pull->GetName() << " " << pull->GetType() << std::endl;
}
TH1D dialhist = pull->GetDataHist();
TH1D minhist = pull->GetMinHist();
TH1D maxhist = pull->GetMaxHist();
TH1I typehist = pull->GetDialTypes();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
dialthrowhandle[name] = pull->GetName();
if (fCurVals.find(name) == fCurVals.end()){
// Add to Containers
fParams.push_back(name);
fCurVals[name] = dialhist.GetBinContent(i+1);
fStartVals[name] = dialhist.GetBinContent(i+1);
fMinVals[name] = minhist.GetBinContent(i+1);
fMaxVals[name] = maxhist.GetBinContent(i+1);
fStepVals[name] = 1.0;
fFixVals[name] = false;
fStartFixVals[name] = false;
fTypeVals[name] = typehist.GetBinContent(i+1);
fStateVals[name] = "FREE" + pull->GetType();
// Maker Helper
helperstr << std::string(16, ' ' ) << FitBase::ConvDialType(fTypeVals[name]) << " "
<< name << " " << fMinVals[name] << " "
<< fMaxVals[name] << " " << fStepVals[name] << " " << fStateVals[name]
<< std::endl;
}
}
}
// Check if no throws given
if (fInputThrows.empty()){
ERR(WRN) << "No covariances given to nuissyst" << std::endl;
ERR(WRN) << "Pushing back an uncorrelated gaussian throw error for each free parameter using step size" << std::endl;
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (fFixVals[syst]) continue;
// Make Terms
std::string name = syst + "_pull";
std::ostringstream pullterm;
pullterm << "DIAL:" << syst << ";"
<< fStartVals[syst] << ";"
<< fStepVals[syst];
std::string type = "GAUSTHROW/NEUT";
// Push Back Pulls
ParamPull* pull = new ParamPull( name, pullterm.str(), type );
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
// Print Whats added
ERR(WRN) << "Added ParamPull : " << name << " " << pullterm.str() << " " << type << std::endl;
// Add helper string for future fits
helperstr << std::string(16, ' ' ) << "covar " << name << " " << pullterm.str() << " " << type << std::endl;
// Keep Track of Throws
dialthrowhandle[syst] = pull->GetName();
}
}
// Print Helper String
if (!helperstr.str().empty()){
LOG(FIT) << "To remove these statements in future studies, add the lines below to your card:" << std::endl;
// Can't use the logger properly because this can be multi-line. Use cout and added spaces to look better!
std::cout << helperstr.str();
sleep(2);
}
// Print Throw State
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (dialthrowhandle.find(syst) != dialthrowhandle.end()){
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = THROWING with " << dialthrowhandle[syst] << std::endl;
} else {
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = FIXED" << std::endl;
}
}
// Pause anyway
sleep(1);
return;
}
/*
Fitting Functions
*/
//*************************************
void SystematicRoutines::UpdateRWEngine(std::map<std::string,double>& updateVals){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name,updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void SystematicRoutines::PrintState(){
//*************************************
LOG(FIT)<<"------------"<<std::endl;
// Count max size
int maxcount = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
maxcount = max(int(fParams[i].size()), maxcount);
}
// Header
LOG(FIT) << " # " << left << setw(maxcount) << "Parameter "
<< " = "
<< setw(10) << "Value" << " +- "
<< setw(10) << "Error" << " "
<< setw(8) << "(Units)" << " "
<< setw(10) << "Conv. Val" << " +- "
<< setw(10) << "Conv. Err" << " "
<< setw(8) << "(Units)" << std::endl;
// Parameters
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams.at(i);
std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
std::string curunits = "(sig.)";
double curval = fCurVals[syst];
double curerr = fErrorVals[syst];
if (fStateVals[syst].find("ABS") != std::string::npos){
curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
curunits = "(Abs.)";
} else if (fStateVals[syst].find("FRAC") != std::string::npos){
curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
FitBase::RWSigmaToFrac(typestr, syst, 0.0));
curunits = "(Frac)";
}
std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
std::ostringstream curparstring;
curparstring << " " << setw(3) << left
<< i << ". "
<< setw(maxcount) << syst << " = "
<< setw(10) << curval << " +- "
<< setw(10) << curerr << " "
<< setw(8) << curunits << " "
<< setw(10) << convval << " +- "
<< setw(10) << converr << " "
<< setw(8) << convunits;
LOG(FIT) << curparstring.str() << std::endl;
}
LOG(FIT)<<"------------"<<std::endl;
double like = fSampleFCN->GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT)<<"------------"<<std::endl;
}
/*
Write Functions
*/
//*************************************
void SystematicRoutines::SaveResults(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
SaveCurrentState();
}
//*************************************
void SystematicRoutines::SaveCurrentState(std::string subdir){
//*************************************
LOG(FIT)<<"Saving current full FCN predictions" <<std::endl;
// Setup DIRS
TDirectory* curdir = gDirectory;
if (!subdir.empty()){
TDirectory* newdir =(TDirectory*) gDirectory->mkdir(subdir.c_str());
newdir->cd();
}
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void SystematicRoutines::SaveNominal(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Nominal Predictions (be cautious with this)" <<std::endl;
FitBase::GetRW()->Reconfigure();
SaveCurrentState("nominal");
};
//*************************************
void SystematicRoutines::SavePrefit(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Prefit Predictions"<<std::endl;
UpdateRWEngine(fStartVals);
SaveCurrentState("prefit");
UpdateRWEngine(fCurVals);
};
/*
MISC Functions
*/
//*************************************
int SystematicRoutines::GetStatus(){
//*************************************
return 0;
}
//*************************************
void SystematicRoutines::SetupCovariance(){
//*************************************
// Remove covares if they exist
if (fCovar) delete fCovar;
if (fCovarFree) delete fCovarFree;
if (fCorrel) delete fCorrel;
if (fCorrelFree) delete fCorrelFree;
if (fDecomp) delete fDecomp;
if (fDecompFree) delete fDecompFree;
int NFREE = 0;
int NDIM = 0;
// Get NFREE from min or from vals (for cases when doing throws)
NDIM = fParams.size();
for (UInt_t i = 0; i < fParams.size(); i++){
if (!fFixVals[fParams[i]]) NFREE++;
}
if (NDIM == 0) return;
fCovar = new TH2D("covariance","covariance",NDIM,0,NDIM,NDIM,0,NDIM);
if (NFREE > 0){
fCovarFree = new TH2D("covariance_free",
"covariance_free",
NFREE,0,NFREE,
NFREE,0,NFREE);
}
// Set Bin Labels
int countall = 0;
int countfree = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
fCovar->GetXaxis()->SetBinLabel(countall+1,fParams[i].c_str());
fCovar->GetYaxis()->SetBinLabel(countall+1,fParams[i].c_str());
countall++;
if (!fFixVals[fParams[i]] and NFREE > 0){
fCovarFree->GetXaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
fCovarFree->GetYaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
countfree++;
}
}
fCorrel = PlotUtils::GetCorrelationPlot(fCovar,"correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar,"decomposition");
if (NFREE > 0)fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free");
if (NFREE > 0)fDecompFree = PlotUtils::GetDecompPlot(fCovarFree,"decomposition_free");
return;
};
//*************************************
void SystematicRoutines::ThrowCovariance(bool uniformly){
//*************************************
// Set fThrownVals to all values in currentVals
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams.at(i);
fThrownVals[name] = fCurVals[name];
}
for (PullListConstIter iter = fInputThrows.begin();
iter != fInputThrows.end(); iter++){
ParamPull* pull = *iter;
pull->ThrowCovariance();
TH1D dialhist = pull->GetDataHist();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
if (fCurVals.find(name) != fCurVals.end()){
fThrownVals[name] = dialhist.GetBinContent(i+1);
}
}
// Reset throw incase pulls are calculated.
pull->ResetToy();
}
return;
};
//*************************************
void SystematicRoutines::PlotLimits(){
//*************************************
std::cout << "Plotting Limits" << std::endl;
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
TDirectory* limfolder = (TDirectory*) fOutputRootFile->mkdir("Limits");
limfolder->cd();
// Set all parameters at their starting values
for (UInt_t i = 0; i < fParams.size(); i++){
fCurVals[fParams[i]] = fStartVals[fParams[i]];
}
TDirectory* nomfolder = (TDirectory*) limfolder->mkdir("nominal");
nomfolder->cd();
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
limfolder->cd();
std::vector<std::string> allfolders;
// Loop through each parameter
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
std::cout << "Starting Param " << syst << std::endl;
if (fFixVals[syst]) continue;
// Loop Downwards
while (fCurVals[syst] > fMinVals[syst]){
fCurVals[syst] = fCurVals[syst] - fStepVals[syst];
// Check Limit
if (fCurVals[syst] < fMinVals[syst])
fCurVals[syst] = fMinVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder for variation
TDirectory* minfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
minfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to folder
fSampleFCN->Write();
}
// Reset before next loop
fCurVals[syst] = fStartVals[syst];
// Loop Upwards now
while (fCurVals[syst] < fMaxVals[syst]){
fCurVals[syst] = fCurVals[syst] + fStepVals[syst];
// Check Limit
if (fCurVals[syst] > fMaxVals[syst])
fCurVals[syst] = fMaxVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder
TDirectory* maxfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
maxfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to file
fSampleFCN->Write();
}
// Reset before leaving
fCurVals[syst] = fStartVals[syst];
UpdateRWEngine(fCurVals);
}
return;
}
//*************************************
void SystematicRoutines::Run(){
//*************************************
std::cout << "Running routines "<< std::endl;
fRoutines = GeneralUtils::ParseToStr(fStrategy,",");
for (UInt_t i = 0; i < fRoutines.size(); i++){
std::string routine = fRoutines.at(i);
int fitstate = kFitUnfinished;
LOG(FIT)<<"Running Routine: "<<routine<<std::endl;
if (routine.compare("PlotLimits") == 0) PlotLimits();
else if (routine.compare("ErrorBands") == 0) GenerateErrorBands();
else if (routine.compare("ThrowErrors") == 0) GenerateThrows();
else if (routine.compare("MergeErrors") == 0) MergeThrows();
else {
std::cout << "Unknown ROUTINE : " << routine << std::endl;
}
// If ending early break here
if (fitstate == kFitFinished || fitstate == kNoChange){
LOG(FIT) << "Ending fit routines loop." << std::endl;
break;
}
}
return;
}
void SystematicRoutines::GenerateErrorBands(){
GenerateThrows();
MergeThrows();
}
//*************************************
void SystematicRoutines::GenerateThrows(){
//*************************************
TFile* tempfile = new TFile((fOutputFile + ".throws.root").c_str(),"RECREATE");
tempfile->cd();
int nthrows = fNThrows;
int startthrows = fStartThrows;
int endthrows = startthrows + nthrows;
if (nthrows < 0) nthrows = endthrows;
if (startthrows < 0) startthrows = 0;
if (endthrows < 0) endthrows = startthrows + nthrows;
int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + endthrows) + time(NULL))/35;
gRandom->SetSeed(seed);
LOG(FIT) << "Using Seed : " << seed << std::endl;
LOG(FIT) << "nthrows = " << nthrows << std::endl;
LOG(FIT) << "startthrows = " << startthrows << std::endl;
LOG(FIT) << "endthrows = " << endthrows << std::endl;
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
if (startthrows == 0){
LOG(FIT) << "Making nominal " << std::endl;
TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal");
nominal->cd();
fSampleFCN->Write();
}
LOG(SAM) << "nthrows = " << nthrows << std::endl;
LOG(SAM) << "startthrows = " << startthrows << std::endl;
LOG(SAM) << "endthrows = " << endthrows << std::endl;
TTree* parameterTree = new TTree("throws","throws");
double chi2;
for (UInt_t i = 0; i < fParams.size(); i++)
parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str());
parameterTree->Branch("chi2",&chi2,"chi2/D");
fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW());
// Would anybody actually want to do uniform throws of any parameter??
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Run Throws and save
for (Int_t i = 0; i < endthrows+1; i++){
LOG(FIT) << "Loop " << i << std::endl;
ThrowCovariance(uniformly);
if (i < startthrows) continue;
if (i == 0) continue;
LOG(FIT) << "Throw " << i << " ================================" << std::endl;
// Generate Random Parameter Throw
// ThrowCovariance(uniformly);
TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i",i));
throwfolder->cd();
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Save the FCN
fSampleFCN->Write();
parameterTree->Fill();
}
+ tempfile->cd();
fSampleFCN->WriteIterationTree();
tempfile->Close();
}
void SystematicRoutines::MergeThrows(){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
// Make a container folder
TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands");
errorDIR->cd();
TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw");
outnominal->cd();
// Split Input Files
if (!fThrowString.empty()) fThrowList = GeneralUtils::ParseToStr(fThrowString,",");
// Add default if no throwlist given
if (fThrowList.size() < 1) fThrowList.push_back( fOutputFile + ".throws.root" );
/// Save location of file containing nominal
std::string nominalfile;
bool nominalfound;
// Loop over files and check they exist.
for (uint i = 0; i < fThrowList.size(); i++){
std::string file = fThrowList[i];
bool found = false;
// normal
std::string newfile = file;
TFile* throwfile = new TFile(file.c_str(),"READ");
if (throwfile and !throwfile->IsZombie()){
found = true;
}
// normal.throws.root
if (!found){
newfile = file + ".throws.root";
throwfile = new TFile((file + ".throws.root").c_str(),"READ");
if (throwfile and !throwfile->IsZombie()) {
found = true;
}
}
// If its found save to throwlist, else save empty.
// Also search for nominal
if (found){
fThrowList[i] = newfile;
LOG(FIT) << "Throws File :" << newfile << std::endl;
// Find input which contains nominal
if (throwfile->Get("nominal")){
nominalfound = true;
nominalfile = newfile;
}
throwfile->Close();
} else {
fThrowList[i] = "";
}
delete throwfile;
}
// Make sure we have a nominal file
if (!nominalfound or nominalfile.empty()){
ERR(FTL) << "No nominal found when mergining! Exiting!" << std::endl;
throw;
}
// Get the nominal throws file
TFile* tempfile = new TFile((nominalfile).c_str(),"READ");
tempfile->cd();
TDirectory* nominal = (TDirectory*)tempfile->Get("nominal");
// int nthrows = FitPar::Config().GetParI("error_throws");
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Check percentage of bad files is okay.
int badfilecount = 0;
for (uint i = 0; i < fThrowList.size(); i++){
if (!fThrowList[i].empty()){
LOG(FIT) << "Loading Throws From File " << i << " : "
<< fThrowList[i] << std::endl;
} else {
badfilecount++;
}
}
// Check we have at least one good file
if ((uint)badfilecount == fThrowList.size()){
ERR(FTL) << "Found no good throw files for MergeThrows" << std::endl;
throw;
} else if (badfilecount > fThrowList.size()*0.25){
ERR(WRN) << "Over 25% of your throw files are dodgy. Please check this is okay!" << std::endl;
ERR(WRN) << "Will continue for the time being..." << std::endl;
sleep(5);
}
// Now go through the keys in the temporary file and look for TH1D, and TH2D plots
TIter next(nominal->GetListOfKeys());
TKey *key;
while ((key = (TKey*)next())) {
TClass *cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
TH1* baseplot = (TH1D*)key->ReadObj();
std::string plotname = std::string(baseplot->GetName());
LOG(FIT) << "Creating error bands for " << plotname;
if (LOG_LEVEL(FIT)){
if (!uniformly) std::cout << " : Using COVARIANCE Throws! " << std::endl;
else std::cout << " : Using UNIFORM THROWS!!! " << std::endl;
}
int nbins = 0;
if (cl->InheritsFrom("TH1D")) nbins = ((TH1D*)baseplot)->GetNbinsX();
else nbins = ((TH1D*)baseplot)->GetNbinsX()* ((TH1D*)baseplot)->GetNbinsY();
// Setup TProfile with RMS option
TProfile* tprof = new TProfile((plotname + "_prof").c_str(),(plotname + "_prof").c_str(),nbins, 0, nbins, "S");
// Setup The TTREE
double* bincontents;
bincontents = new double[nbins];
double* binlowest;
binlowest = new double[nbins];
double* binhighest;
binhighest = new double[nbins];
errorDIR->cd();
TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
for (Int_t i = 0; i < nbins; i++){
bincontents[i] = 0.0;
binhighest[i] = 0.0;
binlowest[i] = 0.0;
bintree->Branch(Form("content_%i",i),&bincontents[i],Form("content_%i/D",i));
}
// Make new throw plot
TH1* newplot;
// Run Throw Merging.
for (UInt_t i = 0; i < fThrowList.size(); i++){
TFile* throwfile = new TFile(fThrowList[i].c_str(), "READ");
// Loop over all throws in a folder
TIter nextthrow(throwfile->GetListOfKeys());
TKey *throwkey;
while ((throwkey = (TKey*)nextthrow())) {
// Skip non throw folders
if (std::string(throwkey->GetName()).find("throw_") == std::string::npos) continue;
// Get Throw DIR
TDirectory* throwdir = (TDirectory*)throwkey->ReadObj();
// Get Plot From Throw
newplot = (TH1*)throwdir->Get(plotname.c_str());
if (!newplot) continue;
// Loop Over Plot
for (Int_t j = 0; j < nbins; j++){
tprof->Fill(j+0.5, newplot->GetBinContent(j+1));
bincontents[j] = newplot->GetBinContent(j+1);
if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j];
if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j];
}
errorDIR->cd();
bintree->Fill();
}
}
errorDIR->cd();
if (uniformly){
LOG(FIT) << "Uniformly Calculating Plot Errors!" << std::endl;
}
TH1* statplot = (TH1*) baseplot->Clone();
for (Int_t j = 0; j < nbins; j++){
if (!uniformly){
// if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) < 1.0) {
// baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) + pow(baseplot->GetBinError(j+1),2)));
// } else {
baseplot->SetBinContent(j+1,tprof->GetBinContent(j+1));
baseplot->SetBinError(j+1,tprof->GetBinError(j+1));
// }
} else {
baseplot->SetBinContent(j+1, 0.0);//(binlowest[j] + binhighest[j]) / 2.0);
baseplot->SetBinError(j+1, 0.0); //(binhighest[j] - binlowest[j])/2.0);
}
}
errorDIR->cd();
baseplot->Write();
tprof->Write();
bintree->Write();
outnominal->cd();
for (int i = 0; i < nbins; i++){
baseplot->SetBinError(i+1, sqrt(pow(statplot->GetBinError(i+1),2) + pow(baseplot->GetBinError(i+1),2)));
}
baseplot->Write();
delete statplot;
delete baseplot;
delete tprof;
delete bintree;
delete [] bincontents;
}
return;
};

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:41 AM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4982757
Default Alt Text
(168 KB)

Event Timeline