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 . ################################################################################ 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 + 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 +#include +#include +#include +#include +#include + +#if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT) +#include // for `feenableexcept` +#endif + +#include +#include +#include +#include +#include +#include + +#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 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 gOptNuPDGs; // Beam PDGS +vector gOptNuFluxs; // Beam Fluxes +vector 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(flux_driver)->GetTotalSpectrum(); + vector fluxinputs = static_cast(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(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 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(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 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 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::value_type(pdg, wgt)); + + // For single target number of free nucleons + // automatic + gOptNumberNucleons = pdg % 10000/ 10; + + } else { + using_tgtmix = true; + vector::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::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::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 #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 samplekeys = Config::QueryKeys("sample"); LoadSamples(samplekeys); std::vector 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 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 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 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 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 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 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 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 JointFCN::GetInputList() { std::vector InputList; fIsAllSplines = true; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector 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 JointFCN::GetSubSampleList() { std::vector SampleList; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector 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::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 signalbitset(fSubSampleList.size()); // Create a new signal box vector for this event std::vector signalboxes; // Start measurement iterator size_t measitercount = 0; std::vector::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 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::iterator inpsig_iter = fSignalEventFlags.begin(); std::vector >::iterator box_iter = - fSignalEventBoxes.begin(); + fSignalEventBoxes.begin(); std::vector >::iterator spline_iter = - fSignalEventSplines.begin(); + fSignalEventSplines.begin(); std::vector >::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::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::iterator subsamsig_iter = (*samsig_iter).begin(); std::vector::iterator subbox_iter = - (*box_iter).begin(); + (*box_iter).begin(); // Loop over all sub measurements. std::vector::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 likes; std::vector ndofs; std::vector 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 fInputs = - FitBase::EvtManager().GetInputs(); + FitBase::EvtManager().GetInputs(); std::map::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 #include #include #include // ROOT headers #include "TTree.h" #include "TH1D.h" #include "TH2D.h" #include #include "TGraphErrors.h" #include #include #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 samplekeys, TFile* outfile=NULL); JointFCN(TFile* outfile=NULL); // Loads from global config //! Destructor ~JointFCN(); //! Create sample list from cardfile void LoadSamples(std::vector samplekeys); void LoadPulls(std::vector pullkeys); //! Main Likelihood evaluation FCN double DoEval(const double *x); //! Func Wrapper for ROOT inline double operator() (const std::vector & 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 GetSampleList(){ return fSamples; } //! Return list of pointers to all the pulls inline std::list 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 GetSubSampleList(); std::vector GetInputList(); private: //! Append the experiments to include in the fit to this list std::list fSamples; //! Append parameter pull terms to include penalties in the fit to this list std::list 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 > fSignalEventSplines; std::vector< std::vector > fSignalEventBoxes; std::vector< bool > fSignalEventFlags; std::vector< std::vector > fSampleSignalFlags; std::vector fInputList; std::vector fSubSampleList; bool fIsAllSplines; + + std::vector< int > fIterationCount; + std::vector< double > fCurrentValues; + std::vector< std::string > fNameValues; + std::vector< std::vector > 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 . *******************************************************************************/ #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 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 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 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 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 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 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 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 ParamPull::GetAllDials() { //******************************************************************************* std::map 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 namevec = rw->GetDialNames(); std::vector 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 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 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 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 . *******************************************************************************/ #ifndef PARAM_PULL_H_SEEN #define PARAM_PULL_H_SEEN /*! * \addtogroup FitBase * @{ */ #include #include #include #include #include #include #include #include // ROOT includes #include #include #include #include #include #include #include // 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 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::const_iterator PullListConstIter; typedef std::list::iterator PullListIter; typedef std::vector::const_iterator PullVectConstIter; typedef std::vector::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 . *******************************************************************************/ #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 xmlcmds; std::vector configargs; fNThrows = 250; fStartThrows = 0; fThrowString = ""; // Make easier to handle arguments. std::vector 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 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 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 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 cardlines = GeneralUtils::ParseFileToStr(cardfile,"\n"); FitPar::Config().cardLines = cardlines; // Read Samples first (norm params can be overridden) int linecount = 0; for (std::vector::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::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 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 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 :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 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"<Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->SetFakeData("MC"); UpdateRWEngine(fCurVals); LOG(FIT)<<"Set all data to fake MC predictions."<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 dialthrowhandle; // Get Covariance Objects from FCN std::list 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& 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)<<"------------"<GetLikelihood(); LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl; LOG(FIT)<<"------------"<cd(); SaveCurrentState(); } //************************************* void SystematicRoutines::SaveCurrentState(std::string subdir){ //************************************* LOG(FIT)<<"Saving current full FCN predictions" <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)" <Reconfigure(); SaveCurrentState("nominal"); }; //************************************* void SystematicRoutines::SavePrefit(){ //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); LOG(FIT)<<"Saving Prefit Predictions"< 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 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: "<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; };