diff --git a/src/FitBase/ParamPull.cxx b/src/FitBase/ParamPull.cxx index cd101bc..449789a 100644 --- a/src/FitBase/ParamPull.cxx +++ b/src/FitBase/ParamPull.cxx @@ -1,849 +1,857 @@ // 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; + ERR(FTL) << "Need FIT, ROOT, VECT or DIAL" << std::endl; throw; } // Check Dials are all good - CheckDialsValid(); + if (!CheckDialsValid()) { + ERR(FTL) << "DIALS NOT VALID" << std::endl; + throw; + } // 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) { + LOG(FIT) << "No minimum histogram found for pull parameters, setting to be content - 1E6..." << std::endl; 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) { + LOG(FIT) << "No maximum histogram found for pull parameters, setting to be content - 1E6..." << std::endl; 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 + // Not really sure when or if this is ever used 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; - } - + CheckHist(fDataHist); fDataHist->SetDirectory(0); fDataHist->SetNameTitle( (fName + "_data").c_str(), (fName + " data" + fPlotTitles).c_str() ); + fMinHist = (TH1D*)tempfile->Get("min_dials_free"); + CheckHist(fMinHist); + fMinHist->SetDirectory(0); + fMinHist->SetNameTitle( (fName + "_min").c_str(), + (fName + " min" + fPlotTitles).c_str() ); + + fMaxHist = (TH1D*)tempfile->Get("max_dials_free"); + CheckHist(fMaxHist); + fMaxHist->SetDirectory(0); + fMaxHist->SetNameTitle( (fName + "_max").c_str(), + (fName + " max" + 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()){ LOG(FIT) << "Looking for ParamPull input inside database" << std::endl; filename = FitPar::GetDataBase() + "/" + filename; 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 (FitBase::GetRW()->DialIncluded(name)) { + LOG(DEB) << "Found dial " << name << " in covariance " << fInput << " and matched to reweight engine " << std::endl; + 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(FTL) << "Dial(s) included in covar but not set in FitWeight." << std::endl; + ERR(FTL) << "ParamPulls needs to know how you want it to be treated." << std::endl; + ERR(FTL) << "Include the following lines into your card to throw UNCORRELATED:" << std::endl - ERR(WRN) << helpstring << std::endl; + << helpstring << std::endl; throw; return false; } else { return true; } + return false; } //******************************************************************************* 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(); } return ndof; }; //******************************************************************************* void ParamPull::ThrowCovariance() { //******************************************************************************* // Reset toy for throw ResetToy(); 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; 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() ); // 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(); + LOG(FIT) << "Threw outside allowed region, rethrowing..." << std::endl; + 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); + TH2D tempCov = TH2D(*fDecomp); 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; }; +void ParamPull::CheckHist(TH1D* hist) { + if (!hist) { + ERR(FTL) << "Can't find TH1D hist fit_dials in " << fName << std::endl; + ERR(FTL) << "File Entries:" << std::endl; + TFile *temp = new TFile(fInput.c_str(), "open"); + temp->ls(); + throw; + } +} diff --git a/src/FitBase/ParamPull.h b/src/FitBase/ParamPull.h index aa99977..e62acb5 100644 --- a/src/FitBase/ParamPull.h +++ b/src/FitBase/ParamPull.h @@ -1,194 +1,195 @@ // 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 #include "TRandom3.h" // Fit Includes #include "PlotUtils.h" #include "StatUtils.h" #include "FitWeight.h" #include "FitLogger.h" #include "EventManager.h" #include "TVector.h" using namespace std; //! Enums to allow Non Guassian Pulls in Future enum FitPullTypes { kUnknownPull = -1, kNoPull = 0, kGausPull = 1 }; enum FitThrowTypes { kUnknownThrow = -1, kNoThrow = 0, kGausThrow = 1, 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: + void CheckHist(TH1D*); 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/ComparisonRoutines.cxx b/src/Routines/ComparisonRoutines.cxx index 48b8abb..011e70d 100755 --- a/src/Routines/ComparisonRoutines.cxx +++ b/src/Routines/ComparisonRoutines.cxx @@ -1,536 +1,536 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "ComparisonRoutines.h" /* Constructor/Destructor */ //************************ void ComparisonRoutines::Init() { //************************ fOutputFile = ""; fOutputRootFile = NULL; fStrategy = "Compare"; fRoutines.clear(); fCardFile = ""; fFakeDataInput = ""; fSampleFCN = NULL; fAllowedRoutines = ("Compare"); }; //************************************* ComparisonRoutines::~ComparisonRoutines() { //************************************* delete fOutputRootFile; }; /* Input Functions */ //************************************* ComparisonRoutines::ComparisonRoutines(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; // 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, "-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 nuiskey fCompKey; 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"); bool trace = Config::GetParB("TRACE"); std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl; std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl; SETVERBOSITY(verbocount); SETTRACE(trace); // Comparison Setup ======================================== // Proper Setup fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); SetupComparisonsFromXML(); SetupRWEngine(); SetupFCN(); return; }; //************************************* void ComparisonRoutines::SetupComparisonsFromXML() { //************************************* LOG(FIT) << "Setting up nuiscomp" << 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; ERR(FTL) << "type='PARAMETER_TYPE'" << std::endl; throw; } else if (!key.Has("name")) { ERR(FTL) << "No name given for parameter " << i << std::endl; ERR(FTL) << "name='SAMPLE_NAME'" << std::endl; throw; } else if (!key.Has("nominal")) { ERR(FTL) << "No nominal given for parameter " << i << std::endl; ERR(FTL) << "nominal='NOMINAL_VALUE'" << std::endl; throw; } // Get Inputs std::string partype = key.GetS("type"); std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); 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"); // Check for incomplete limtis int limdef = ((int)key.Has("low") + (int)key.Has("high") + (int)key.Has("step")); if (limdef > 0 and limdef < 3){ ERR(FTL) << "Incomplete limit set given for parameter : " << parname << std::endl; ERR(FTL) << "Requires: low='LOWER_LIMIT' high='UPPER_LIMIT' step='STEP_SIZE' " << std::endl; throw; } // 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; } // Convert if required 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);; + fTypeVals[parname] = FitBase::ConvDialType(partype); fCurVals[parname] = parnom; fStateVals[parname] = parstate; } // 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 " << i << ". : " << samplename << " (" << sampletype << ") [Norm=" << samplenorm<<"]"<< std::endl << " -> input='" << samplefile << "'" << std::endl; // If FREE add to parameters otherwise continue if (sampletype.find("FREE") == std::string::npos) { if (samplenorm != 1.0) { ERR(FTL) << "You provided a sample normalisation but did not specify that the sample is free" << std::endl; ERR(FTL) << "Change so sample contains type=\"FREE\" and re-run" << std::endl; throw; } 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; } // 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("nominal")) { ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl; throw; } // Get Inputs std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); // Push into vectors fFakeVals[parname] = parnom; } } //************************************* void ComparisonRoutines::SetupRWEngine() { //************************************* LOG(FIT) << "Setting up FitWeight Engine" << std::endl; for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name)); } return; } //************************************* void ComparisonRoutines::SetupFCN() { //************************************* LOG(FIT) << "Building the SampleFCN" << std::endl; if (fSampleFCN) delete fSampleFCN; Config::Get().out = fOutputRootFile; fOutputRootFile->cd(); fSampleFCN = new JointFCN(fOutputRootFile); SetFakeData(); return; } //************************************* void ComparisonRoutines::SetFakeData() { //************************************* if (fFakeDataInput.empty()) return; if (fFakeDataInput.compare("MC") == 0) { LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl; UpdateRWEngine(fFakeVals); FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->SetFakeData("MC"); UpdateRWEngine(fCurVals); LOG(FIT) << "Set all data to fake MC predictions." << std::endl; } else { LOG(FIT) << "Setting fake data from: " << fFakeDataInput << std::endl; fSampleFCN->SetFakeData(fFakeDataInput); } return; } /* Fitting Functions */ //************************************* void ComparisonRoutines::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 ComparisonRoutines::Run() { //************************************* LOG(FIT) << "Running ComparisonRoutines : " << fStrategy << std::endl; if (FitPar::Config().GetParB("save_nominal")) { SaveNominal(); } // Parse given routines fRoutines = GeneralUtils::ParseToStr(fStrategy, ","); if (fRoutines.empty()) { ERR(FTL) << "Trying to run ComparisonRoutines with no routines given!" << std::endl; throw; } for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); LOG(FIT) << "Routine: " << routine << std::endl; if (!routine.compare("Compare")) { UpdateRWEngine(fCurVals); GenerateComparison(); PrintState(); SaveCurrentState(); } } return; } //************************************* void ComparisonRoutines::GenerateComparison() { //************************************* LOG(FIT) << "Generating Comparison." << std::endl; // Main Event Loop from event Manager fSampleFCN->ReconfigureAllEvents(); return; } //************************************* void ComparisonRoutines::PrintState() { //************************************* LOG(FIT) << "------------" << std::endl; // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header LOG(FIT) << " # " << left << setw(maxcount) << "Parameter " << " = " << setw(10) << "Value" << " +- " << setw(10) << "Error" << " " << setw(8) << "(Units)" << " " << setw(10) << "Conv. Val" << " +- " << setw(10) << "Conv. Err" << " " << setw(8) << "(Units)" << std::endl; // Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); std::string typestr = FitBase::ConvDialType(fTypeVals[syst]); std::string curunits = "(sig.)"; double curval = fCurVals[syst]; double curerr = 0.0; if (fStateVals[syst].find("ABS") != std::string::npos) { curval = FitBase::RWSigmaToAbs(typestr, syst, curval); curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); curunits = "(Abs.)"; } else if (fStateVals[syst].find("FRAC") != std::string::npos) { curval = FitBase::RWSigmaToFrac(typestr, syst, curval); curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) - FitBase::RWSigmaToFrac(typestr, syst, 0.0)); curunits = "(Frac)"; } std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")"; double convval = FitBase::RWSigmaToAbs(typestr, syst, curval); double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); std::ostringstream curparstring; curparstring << " " << setw(3) << left << i << ". " << setw(maxcount) << syst << " = " << setw(10) << curval << " +- " << setw(10) << curerr << " " << setw(8) << curunits << " " << setw(10) << convval << " +- " << setw(10) << converr << " " << setw(8) << convunits; LOG(FIT) << curparstring.str() << std::endl; } LOG(FIT) << "------------" << std::endl; double like = fSampleFCN->GetLikelihood(); LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl; LOG(FIT) << "------------" << std::endl; } /* Write Functions */ //************************************* void ComparisonRoutines::SaveCurrentState(std::string subdir) { //************************************* LOG(FIT) << "Saving current full FCN predictions" << std::endl; // Setup DIRS TDirectory* curdir = gDirectory; if (!subdir.empty()) { TDirectory* newdir = (TDirectory*)gDirectory->mkdir(subdir.c_str()); newdir->cd(); } fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void ComparisonRoutines::SaveNominal() { //************************************* fOutputRootFile->cd(); LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl; FitBase::GetRW()->Reconfigure(); GenerateComparison(); SaveCurrentState("nominal"); }; diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx index 49988e6..04c98f6 100755 --- a/src/Routines/SystematicRoutines.cxx +++ b/src/Routines/SystematicRoutines.cxx @@ -1,1457 +1,1437 @@ // 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.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(){ 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; } + } /* 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; + LOG(FIT) << "Loading ParamPull objects from FCN to build covariance..." << 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++){ + for (PullListConstIter iter = inputpulls.begin(); iter != inputpulls.end(); iter++){ ParamPull* pull = (*iter); - - if (pull->GetType().find("THROW")){ + if (pull->GetType().find("THROW") != std::string::npos){ 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; + LOG(FIT) << "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. + // Reset throw in case 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; + // 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); 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; + LOG(FIT) << "Throw " << i << "/" << endthrows << " ================================" << 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(); +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; + 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 + // 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; + << 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(); } throwfile->Close(); delete throwfile; } 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)); - // } + // 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->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; }; -void SystematicRoutines::EigenErrors(){ +void SystematicRoutines::EigenErrors() { - fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); - fOutputRootFile->cd(); + 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++){ + iter != fInputThrows.end(); iter++){ ParamPull* pull = *iter; // Check pull is actualyl Gaussian std::string pulltype = pull->GetType(); if (pulltype.find("GAUSTHROW") == std::string::npos){ THROW("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); - } - - } - } - + + 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(); - */ + 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 LOG(FIT) << "Parameter Set " << count << std::endl; for (int j = 0; j < eigenVect.GetNrows(); j++){ std::string param = fParams[j]; LOG(FIT) << " " << j << ". " << param << " : " << fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i] << std::endl; fThrownVals[param] = fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i]; } // Run Eval double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals ); double chi2 = fSampleFCN->DoEval( vals ); delete vals; count++; fSampleFCN->Write(); - + throwfolder = (TDirectory*)throwsdir->mkdir(Form("throw_%i",count)); throwfolder->cd(); // Get New Parameter Vector LOG(FIT) << "Parameter Set " << count << std::endl; for (int j = 0; j < eigenVect.GetNrows(); j++){ std::string param = fParams[j]; LOG(FIT) << " " << j << ". " << param << " : " <DoEval( vals2 ); 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; - + LOG(FIT) << "Creating error bands for " << key->GetName() << std::endl; 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){ - ERR(WRN) << "Cannot find new plot : " << Form("throw_%i/%s",i,plotname.c_str()) << std::endl; - ERR(WRN) << "This plot will not have the correct errors!" << std::endl; - continue; + ERR(WRN) << "Cannot find new plot : " << Form("throw_%i/%s",i,plotname.c_str()) << std::endl; + ERR(WRN) << "This plot will not have the correct errors!" << std::endl; + continue; } newplot->SetDirectory(0); nbins = newplot->GetNbinsX(); - + for (int j = 0; j < nbins; j++){ - if (i % 2 == 0){ - // std::cout << plotname<< " : upper " << errorplot_upper->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl; - // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " <GetBinContent(j+1) << std::endl; - errorplot_upper->SetBinContent(j+1, errorplot_upper->GetBinContent(j+1) + - pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2)); - // newplot->Print(); - } else { - // std::cout << plotname << " : lower " << errorplot_lower->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl; - // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " << newplot->GetBinContent(j+1) << std::endl; - errorplot_lower->SetBinContent(j+1, errorplot_lower->GetBinContent(j+1) + - pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2)); - // newplot->Print(); - } - meanplot->SetBinContent(j+1, meanplot->GetBinContent(j+1) + baseplot->GetBinContent(j+1)); + if (i % 2 == 0){ + // std::cout << plotname<< " : upper " << errorplot_upper->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl; + // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " <GetBinContent(j+1) << std::endl; + errorplot_upper->SetBinContent(j+1, errorplot_upper->GetBinContent(j+1) + + pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2)); + // newplot->Print(); + } else { + // std::cout << plotname << " : lower " << errorplot_lower->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl; + // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " << newplot->GetBinContent(j+1) << std::endl; + errorplot_lower->SetBinContent(j+1, errorplot_lower->GetBinContent(j+1) + + pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2)); + // newplot->Print(); + } + 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; } }