diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx index cb205f7..a89fe2f 100755 --- a/src/Routines/SystematicRoutines.cxx +++ b/src/Routines/SystematicRoutines.cxx @@ -1,1490 +1,1492 @@ // Copyright 2016-2021 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()) { NUIS_ABORT("No input supplied!"); } if (fOutputFile.empty() and !fCardFile.empty()) { fOutputFile = fCardFile + ".root"; NUIS_ERR(WRN, "No output supplied so saving it to: " << fOutputFile); } else if (fOutputFile.empty()) { NUIS_ABORT("No output file or cardfile supplied!"); } // 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.Set("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile); if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy); // Load XML Cardfile configuration.LoadSettings(fCompKey.GetS("cardfile"), ""); // 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.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::GetParI("VERBOSITY"); errorcount += Config::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() { NUIS_LOG(FIT, "Setting up nuismin"); // Setup Parameters ------------------------------------------ std::vector parkeys = Config::QueryKeys("parameter"); if (!parkeys.empty()) { NUIS_LOG(FIT, "Number of parameters : " << parkeys.size()); } for (size_t i = 0; i < parkeys.size(); i++) { nuiskey key = parkeys.at(i); // Check for type,name,nom if (!key.Has("type")) { NUIS_ABORT("No type given for parameter " << i); } else if (!key.Has("name")) { NUIS_ABORT("No name given for parameter " << i); } else if (!key.Has("nominal")) { NUIS_ABORT("No nominal given for parameter " << i); } // 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"); NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom << " : " << parlow << " < p < " << parhigh << " : " << parstate); } else { NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom << " : " << parstate); } // 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()) { NUIS_LOG(FIT, "Number of samples : " << samplekeys.size()); } 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 NUIS_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); // 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()) { NUIS_LOG(FIT, "Number of fake parameters : " << fakekeys.size()); } for (size_t i = 0; i < fakekeys.size(); i++) { nuiskey key = fakekeys.at(i); // Check for type,name,nom if (!key.Has("name")) { NUIS_ABORT("No name given for fakeparameter " << i); } else if (!key.Has("nom")) { NUIS_ABORT("No nominal given for fakeparameter " << i); } // Get Inputs std::string parname = key.GetS("name"); double parnom = key.GetD("nom"); // Push into vectors fFakeVals[parname] = parnom; } } /* 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() { //************************************* NUIS_LOG(FIT, "Making the jointFCN"); if (fSampleFCN) delete fSampleFCN; fSampleFCN = new JointFCN(fOutputRootFile); SetFakeData(); return; } //************************************* void SystematicRoutines::SetFakeData() { //************************************* if (fFakeDataInput.empty()) return; if (fFakeDataInput.compare("MC") == 0) { NUIS_LOG(FIT, "Setting fake data from MC starting prediction."); UpdateRWEngine(fFakeVals); FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->SetFakeData("MC"); UpdateRWEngine(fCurVals); NUIS_LOG(FIT, "Set all data to fake MC predictions."); } else { fSampleFCN->SetFakeData(fFakeDataInput); } return; } //***************************************** // Setup the parameter covariances from the FCN void SystematicRoutines::GetCovarFromFCN() { //***************************************** NUIS_LOG(FIT, "Loading ParamPull objects from FCN to build covariance..."); // 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") != std::string::npos) { fInputThrows.push_back(pull); fInputCovar.push_back(pull->GetFullCovarMatrix()); fInputDials.push_back(pull->GetDataHist()); NUIS_LOG(FIT, "Read ParamPull: " << pull->GetName() << " " << pull->GetType()); } 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(); // Add to Containers // Set the starting values etc to the postfit fParams.push_back(name); fCurVals[name] = dialhist.GetBinContent(i + 1); // Set the starting values to be nominal MC fStartVals[name] = 0.0; 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(); // If we find the string if (fCurVals.find(name) == fCurVals.end()) { // 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()) { NUIS_ERR(WRN, "No covariances given to nuissyst"); NUIS_ERR(WRN, "Pushing back an uncorrelated gaussian throw error for each " "free parameter using step size"); 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 NUIS_LOG(FIT, "Added ParamPull : " << name << " " << pullterm.str() << " " << type); // 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()) { NUIS_LOG(FIT, "To remove these statements in future studies, add the lines " "below to your card:"); // 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()) { NUIS_LOG(FIT, "Dial " << i << ". " << setw(20) << syst << " = THROWing with " << dialthrowhandle[syst]); } else { NUIS_LOG(FIT, "Dial " << i << ". " << setw(20) << syst << " = is FIXED"); } } } /* 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(); } //************************************* void SystematicRoutines::PrintState() { //************************************* NUIS_LOG(FIT, "------------"); // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header NUIS_LOG(FIT, " # " << left << setw(maxcount) << "Parameter " << " = " << setw(10) << "Value" << " +- " << setw(10) << "Error" << " " << setw(8) << "(Units)" << " " << setw(10) << "Conv. Val" << " +- " << setw(10) << "Conv. Err" << " " << setw(8) << "(Units)"); // Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); std::string typestr = FitBase::ConvDialType(fTypeVals[syst]); std::string curunits = "(sig.)"; double curval = fCurVals[syst]; double curerr = fErrorVals[syst]; if (fStateVals[syst].find("ABS") != std::string::npos) { curval = FitBase::RWSigmaToAbs(typestr, syst, curval); curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); curunits = "(Abs.)"; } else if (fStateVals[syst].find("FRAC") != std::string::npos) { curval = FitBase::RWSigmaToFrac(typestr, syst, curval); curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) - FitBase::RWSigmaToFrac(typestr, syst, 0.0)); curunits = "(Frac)"; } std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")"; double convval = FitBase::RWSigmaToAbs(typestr, syst, curval); double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); std::ostringstream curparstring; curparstring << " " << setw(3) << left << i << ". " << setw(maxcount) << syst << " = " << setw(10) << curval << " +- " << setw(10) << curerr << " " << setw(8) << curunits << " " << setw(10) << convval << " +- " << setw(10) << converr << " " << setw(8) << convunits; NUIS_LOG(FIT, curparstring.str()); } NUIS_LOG(FIT, "------------"); double like = fSampleFCN->GetLikelihood(); NUIS_LOG(FIT, std::left << std::setw(46) << "Likelihood for JointFCN: " << like); NUIS_LOG(FIT, "------------"); } /* Write Functions */ //************************************* void SystematicRoutines::SaveResults() { //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); SaveCurrentState(); } //************************************* void SystematicRoutines::SaveCurrentState(std::string subdir) { //************************************* NUIS_LOG(FIT, "Saving current full FCN predictions"); // Setup DIRS TDirectory *curdir = gDirectory; if (!subdir.empty()) { TDirectory *newdir = (TDirectory *)gDirectory->mkdir(subdir.c_str()); newdir->cd(); } FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void SystematicRoutines::SaveNominal() { //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); NUIS_LOG(FIT, "Saving Nominal Predictions (be cautious with this)"); FitBase::GetRW()->Reconfigure(); SaveCurrentState("nominal"); }; //************************************* void SystematicRoutines::SavePrefit() { //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); NUIS_LOG(FIT, "Saving Prefit Predictions"); UpdateRWEngine(fStartVals); SaveCurrentState("prefit"); UpdateRWEngine(fCurVals); }; /* MISC Functions */ //************************************* int SystematicRoutines::GetStatus() { //************************************* return 0; } //************************************* void SystematicRoutines::SetupCovariance() { //************************************* // Remove covares if they exist if (fCovar) delete fCovar; if (fCovarFree) delete fCovarFree; if (fCorrel) delete fCorrel; if (fCorrelFree) delete fCorrelFree; if (fDecomp) delete fDecomp; if (fDecompFree) delete fDecompFree; int NFREE = 0; int NDIM = 0; // Get NFREE from min or from vals (for cases when doing throws) NDIM = fParams.size(); for (UInt_t i = 0; i < fParams.size(); i++) { if (!fFixVals[fParams[i]]) NFREE++; } if (NDIM == 0) return; fCovar = new TH2D("covariance", "covariance", NDIM, 0, NDIM, NDIM, 0, NDIM); if (NFREE > 0) { fCovarFree = new TH2D("covariance_free", "covariance_free", NFREE, 0, NFREE, NFREE, 0, NFREE); } // Set Bin Labels int countall = 0; int countfree = 0; for (UInt_t i = 0; i < fParams.size(); i++) { fCovar->GetXaxis()->SetBinLabel(countall + 1, fParams[i].c_str()); fCovar->GetYaxis()->SetBinLabel(countall + 1, fParams[i].c_str()); countall++; if (!fFixVals[fParams[i]] and NFREE > 0) { fCovarFree->GetXaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); fCovarFree->GetYaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); countfree++; } } fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation"); fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition"); if (NFREE > 0) fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free"); if (NFREE > 0) fDecompFree = PlotUtils::GetDecompPlot(fCovarFree, "decomposition_free"); return; }; //************************************* void SystematicRoutines::ThrowCovariance(bool uniformly) { //************************************* // Set fThrownVals to all values in currentVals for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams.at(i); fThrownVals[name] = fCurVals[name]; } for (PullListConstIter iter = fInputThrows.begin(); iter != fInputThrows.end(); iter++) { ParamPull *pull = *iter; pull->ThrowCovariance(); TH1D dialhist = pull->GetDataHist(); for (int i = 0; i < dialhist.GetNbinsX(); i++) { std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i + 1)); if (fCurVals.find(name) != fCurVals.end()) { fThrownVals[name] = dialhist.GetBinContent(i + 1); } } // Reset throw in case pulls are calculated. pull->ResetToy(); } }; //************************************* 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() { //************************************* fRoutines = GeneralUtils::ParseToStr(fStrategy, ","); for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); int fitstate = kFitUnfinished; NUIS_LOG(FIT, "Running Routine: " << routine); if (routine.compare("PlotLimits") == 0) PlotLimits(); else if (routine.compare("ErrorBands") == 0) GenerateErrorBands(); else if (routine.compare("ThrowErrors") == 0) GenerateThrows(); else if (routine.compare("MergeErrors") == 0) MergeThrows(); else if (routine.compare("EigenErrors") == 0) EigenErrors(); else { NUIS_ABORT("Unknown ROUTINE : " << routine); } // If ending early break here if (fitstate == kFitFinished || fitstate == kNoChange) { NUIS_LOG(FIT, "Ending fit routines loop."); break; } } } void SystematicRoutines::GenerateErrorBands() { GenerateThrows(); MergeThrows(); } //************************************* void SystematicRoutines::GenerateThrows() { //************************************* TFile *tempfile = new TFile((fOutputFile + ".throws.root").c_str(), "RECREATE"); tempfile->cd(); // For generating throws we check with the config int nthrows = Config::GetParI("error_throws"); int startthrows = fStartThrows; int endthrows = startthrows + nthrows; if (nthrows < 0) nthrows = endthrows; if (startthrows < 0) startthrows = 0; if (endthrows < 0) endthrows = startthrows + nthrows; // Setting Seed // Matteo Mazzanti's Fix struct timeval mytime; gettimeofday(&mytime, NULL); Double_t seed = time(NULL) + int(getpid()) + (mytime.tv_sec * 1000.) + (mytime.tv_usec / 1000.); gRandom->SetSeed(seed); // int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + // endthrows) + time(NULL) + int(getpid()) ); gRandom->SetSeed(seed); NUIS_LOG(FIT, "Using Seed : " << seed); NUIS_LOG(FIT, "nthrows = " << nthrows); NUIS_LOG(FIT, "startthrows = " << startthrows); NUIS_LOG(FIT, "endthrows = " << endthrows); UpdateRWEngine(fStartVals); fSampleFCN->ReconfigureAllEvents(); // Make the nominal if (startthrows == 0) { NUIS_LOG(FIT, "Making nominal "); TDirectory *nominal = (TDirectory *)tempfile->mkdir("nominal"); nominal->cd(); fSampleFCN->Write(); // Make the postfit reading from the pull NUIS_LOG(FIT, "Making postfit "); TDirectory *postfit = (TDirectory *)tempfile->mkdir("postfit"); postfit->cd(); UpdateRWEngine(fCurVals); fSampleFCN->ReconfigureSignal(); fSampleFCN->Write(); } 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++) { // Generate Random Parameter Throw ThrowCovariance(uniformly); if (i < startthrows) continue; if (i == 0) continue; NUIS_LOG(FIT, "Throw " << i << "/" << endthrows << " ================================"); TDirectory *throwfolder = (TDirectory *)tempfile->mkdir(Form("throw_%i", i)); throwfolder->cd(); // Run Eval double *vals = FitUtils::GetArrayFromMap(fParams, fThrownVals); fSampleFCN->DoEval(vals); delete vals; // Save the FCN fSampleFCN->Write(); } tempfile->cd(); fSampleFCN->WriteIterationTree(); tempfile->Close(); } // Merge throws together into one summary 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; NUIS_LOG(FIT, "Throws File :" << newfile); // 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()) { NUIS_ABORT("No nominal found when merging! Exiting!"); } // Get the nominal throws file TFile *tempfile = new TFile((nominalfile).c_str(), "READ"); tempfile->cd(); TDirectory *nominal = (TDirectory *)tempfile->Get("nominal"); 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()) { NUIS_LOG(FIT, "Loading Throws From File " << i << " : " << fThrowList[i]); } else { badfilecount++; } } // Check we have at least one good file if ((uint)badfilecount == fThrowList.size()) { NUIS_ABORT("Found no good throw files for MergeThrows"); throw; } else if (badfilecount > (fThrowList.size() * 0.25)) { NUIS_ERR( WRN, "Over 25% of your throw files are dodgy. Please check this is okay!"); NUIS_ERR(WRN, "Will continue for the time being..."); 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()); - NUIS_LOGN(FIT, "Creating error bands for " << plotname); + NUIS_LOG(FIT, "Creating error bands for " << plotname); if (LOG_LEVEL(FIT)) { if (!uniformly) { NUIS_LOG(FIT, " : Using COVARIANCE Throws! "); } else { NUIS_LOG(FIT, " : Using UNIFORM THROWS!!! "); } } 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(); } throwfile->Close(); delete throwfile; } errorDIR->cd(); if (uniformly) { NUIS_LOG(FIT, "Uniformly Calculating Plot Errors!"); } 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); } } baseplot->SetTitle("Profiled throws"); 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; } +fOutputRootFile->Write(); +fOutputRootFile->Close(); }; void SystematicRoutines::EigenErrors() { fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); // Make Covariance TMatrixDSym *fullcovar = new TMatrixDSym(fParams.size()); // Extract covariance from all loaded ParamPulls for (PullListConstIter iter = fInputThrows.begin(); iter != fInputThrows.end(); iter++) { ParamPull *pull = *iter; // Check pull is actualyl Gaussian std::string pulltype = pull->GetType(); if (pulltype.find("GAUSTHROW") == std::string::npos) { NUIS_ABORT("Can only calculate EigenErrors for Gaussian pulls!"); } // Get data and covariances TH1D dialhist = pull->GetDataHist(); TH2D covhist = pull->GetFullCovar(); // Loop over all dials and compare names for (size_t pari = 0; pari < fParams.size(); pari++) { for (size_t parj = 0; parj < fParams.size(); parj++) { std::string name_pari = fParams[pari]; std::string name_parj = fParams[parj]; // Compare names to those in the pull for (int pulli = 0; pulli < dialhist.GetNbinsX(); pulli++) { for (int pullj = 0; pullj < dialhist.GetNbinsX(); pullj++) { std::string name_pulli = dialhist.GetXaxis()->GetBinLabel(pulli + 1); std::string name_pullj = dialhist.GetXaxis()->GetBinLabel(pullj + 1); if (name_pulli == name_pari && name_pullj == name_parj) { (*fullcovar)[pari][parj] = covhist.GetBinContent(pulli + 1, pullj + 1); fCurVals[name_pari] = dialhist.GetBinContent(pulli + 1); fCurVals[name_parj] = dialhist.GetBinContent(pullj + 1); } } } } } } /* TFile* test = new TFile("testingcovar.root","RECREATE"); test->cd(); TH2D* joinedcov = new TH2D("COVAR","COVAR", fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows()), fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows())); for (int i = 0; i < fullcovar->GetNrows(); i++){ for (int j = 0; j < fullcovar->GetNcols(); j++){ joinedcov->SetBinContent(i+1, j+1, (*fullcovar)[i][j]); } } joinedcov->Write("COVAR"); test->Close(); */ // Calculator all EigenVectors and EigenValues TMatrixDSymEigen *eigen = new TMatrixDSymEigen(*fullcovar); const TVectorD eigenVals = eigen->GetEigenValues(); const TMatrixD eigenVect = eigen->GetEigenVectors(); eigenVals.Print(); eigenVect.Print(); TDirectory *outnominal = (TDirectory *)fOutputRootFile->mkdir("nominal"); outnominal->cd(); // double *valst = FitUtils::GetArrayFromMap(fParams, fCurVals); // double chi2 = fSampleFCN->DoEval( valst ); // delete valst; fSampleFCN->Write(); // Loop over all throws TDirectory *throwsdir = (TDirectory *)fOutputRootFile->mkdir("throws"); throwsdir->cd(); int count = 0; // Produce all error throws. for (int i = 0; i < eigenVect.GetNrows(); i++) { TDirectory *throwfolder = (TDirectory *)throwsdir->mkdir(Form("throw_%i", count)); throwfolder->cd(); // Get New Parameter Vector NUIS_LOG(FIT, "Parameter Set " << count); for (int j = 0; j < eigenVect.GetNrows(); j++) { std::string param = fParams[j]; NUIS_LOG(FIT, " " << j << ". " << param << " : " << fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i]); fThrownVals[param] = fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i]; } // Run Eval double *vals = FitUtils::GetArrayFromMap(fParams, fThrownVals); double chi2 = fSampleFCN->DoEval(vals); NUIS_LOG(DEB, "Chi2 = " << chi2); delete vals; count++; fSampleFCN->Write(); throwfolder = (TDirectory *)throwsdir->mkdir(Form("throw_%i", count)); throwfolder->cd(); // Get New Parameter Vector NUIS_LOG(FIT, "Parameter Set " << count); for (int j = 0; j < eigenVect.GetNrows(); j++) { std::string param = fParams[j]; NUIS_LOG(FIT, " " << j << ". " << param << " : " << fCurVals[param] - sqrt(eigenVals[i]) * eigenVect[j][i]); fThrownVals[param] = fCurVals[param] - sqrt(eigenVals[i]) * eigenVect[j][i]; } // Run Eval double *vals2 = FitUtils::GetArrayFromMap(fParams, fThrownVals); chi2 = fSampleFCN->DoEval(vals2); NUIS_LOG(DEB, "Chi2 = " << chi2); delete vals2; count++; // Save the FCN fSampleFCN->Write(); } fOutputRootFile->Close(); fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "UPDATE"); fOutputRootFile->cd(); throwsdir = (TDirectory *)fOutputRootFile->Get("throws"); outnominal = (TDirectory *)fOutputRootFile->Get("nominal"); // Loop through Error DIR TDirectory *outerr = (TDirectory *)fOutputRootFile->mkdir("errors"); outerr->cd(); TIter next(outnominal->GetListOfKeys()); TKey *key; while ((key = (TKey *)next())) { TClass *cl = gROOT->GetClass(key->GetClassName()); if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue; NUIS_LOG(FIT, "Creating error bands for " << key->GetName()); std::string plotname = std::string(key->GetName()); if (plotname.find("_EVT") != std::string::npos) continue; if (plotname.find("_FLUX") != std::string::npos) continue; if (plotname.find("_FLX") != std::string::npos) continue; TH1 *baseplot = (TH1D *)key->ReadObj()->Clone(Form("%s_ORIGINAL", key->GetName())); TH1 *errorplot_upper = (TH1D *)baseplot->Clone(Form("%s_ERROR_UPPER", key->GetName())); TH1 *errorplot_lower = (TH1D *)baseplot->Clone(Form("%s_ERROR_LOWER", key->GetName())); TH1 *meanplot = (TH1D *)baseplot->Clone(Form("%s_SET_MEAN", key->GetName())); TH1 *systplot = (TH1D *)baseplot->Clone(Form("%s_SYST", key->GetName())); TH1 *statplot = (TH1D *)baseplot->Clone(Form("%s_STAT", key->GetName())); TH1 *totlplot = (TH1D *)baseplot->Clone(Form("%s_TOTAL", key->GetName())); int nbins = 0; if (cl->InheritsFrom("TH1D")) nbins = ((TH1D *)baseplot)->GetNbinsX(); else nbins = ((TH1D *)baseplot)->GetNbinsX() * ((TH1D *)baseplot)->GetNbinsY(); meanplot->Reset(); errorplot_upper->Reset(); errorplot_lower->Reset(); for (int j = 0; j < nbins; j++) { errorplot_upper->SetBinError(j + 1, 0.0); errorplot_lower->SetBinError(j + 1, 0.0); } // Loop over throws and calculate mean and error for +- throws int addcount = 0; // Add baseplot first to slightly bias to central value meanplot->Add(baseplot); addcount++; for (int i = 0; i < count; i++) { TH1 *newplot = (TH1D *)throwsdir->Get(Form("throw_%i/%s", i, plotname.c_str())); if (!newplot) { NUIS_ERR(WRN, "Cannot find new plot : " << Form("throw_%i/%s", i, plotname.c_str())); NUIS_ERR(WRN, "This plot will not have the correct errors!"); continue; } newplot->SetDirectory(0); nbins = newplot->GetNbinsX(); for (int j = 0; j < nbins; j++) { if (i % 2 == 0) { errorplot_upper->SetBinContent( j + 1, errorplot_upper->GetBinContent(j + 1) + pow(baseplot->GetBinContent(j + 1) - newplot->GetBinContent(j + 1), 2)); } else { errorplot_lower->SetBinContent( j + 1, errorplot_lower->GetBinContent(j + 1) + pow(baseplot->GetBinContent(j + 1) - newplot->GetBinContent(j + 1), 2)); } meanplot->SetBinContent(j + 1, meanplot->GetBinContent(j + 1) + baseplot->GetBinContent(j + 1)); } delete newplot; addcount++; } // Get mean Average for (int j = 0; j < nbins; j++) { meanplot->SetBinContent(j + 1, meanplot->GetBinContent(j + 1) / double(addcount)); } for (int j = 0; j < nbins; j++) { errorplot_upper->SetBinContent( j + 1, sqrt(errorplot_upper->GetBinContent(j + 1))); errorplot_lower->SetBinContent( j + 1, sqrt(errorplot_lower->GetBinContent(j + 1))); statplot->SetBinError(j + 1, baseplot->GetBinError(j + 1)); systplot->SetBinError(j + 1, (errorplot_upper->GetBinContent(j + 1) + errorplot_lower->GetBinContent(j + 1)) / 2.0); totlplot->SetBinError(j + 1, sqrt(pow(statplot->GetBinError(j + 1), 2) + pow(systplot->GetBinError(j + 1), 2))); meanplot->SetBinError(j + 1, sqrt(pow(statplot->GetBinError(j + 1), 2) + pow(systplot->GetBinError(j + 1), 2))); } outerr->cd(); errorplot_upper->Write(); errorplot_lower->Write(); baseplot->Write(); meanplot->Write(); statplot->Write(); systplot->Write(); totlplot->Write(); delete errorplot_upper; delete errorplot_lower; delete baseplot; delete meanplot; delete statplot; delete systplot; delete totlplot; } }