diff --git a/parameters/config.xml b/parameters/config.xml index 4e103d0..41cad23 100644 --- a/parameters/config.xml +++ b/parameters/config.xml @@ -1,214 +1,202 @@ - - - - - - - - - - - - - - - + + + + - + - - - + + + - + - - - + + + - diff --git a/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx b/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx index fde3809..fc41380 100755 --- a/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx +++ b/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx @@ -1,153 +1,153 @@ // 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 "MINERvA_SignalDef.h" #include "MINERvA_CCinc_XSec_2DEavq3_nu.h" //******************************************************************** MINERvA_CCinc_XSec_2DEavq3_nu::MINERvA_CCinc_XSec_2DEavq3_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCinc_XSec_2DEavq3_nu sample. \n" \ "Target: CH \n" \ "Flux: MINERvA Medium Energy FHC numu \n" \ "Signal: CC-inclusive with theta < 20deg \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("q_{3} (GeV)"); fSettings.SetYTitle("E_{avail} (GeV)"); fSettings.SetZTitle("d^{2}#sigma/dq_{3}dE_{avail} (cm^{2}/GeV^{2})"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL"); fSettings.SetEnuRange(2.0, 6.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("MINERvA_CCinc_XSec_2DEavq3_nu"); fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/data_2D.txt" ); fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/covar_2D.txt" ); fSettings.SetMapInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/map_2D.txt" ); fSettings.DefineAllowedSpecies("numu"); - hadroncut = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu.hadron_cut"); - useq3true = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu.useq3true"); - splitMEC_PN_NN = FitPar::Config().GetParB("Modes.split_PN_NN"); + hadroncut = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu_hadron_cut"); + useq3true = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu_useq3true"); + splitMEC_PN_NN = FitPar::Config().GetParB("Modes_split_PN_NN"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-42 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); // Plot Setup ------------------------------------------------------- Double_t binx[7] = {0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8}; Double_t biny[17] = {0.0, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.25, 0.30, 0.35, 0.40, 0.50, 0.60, 0.80}; CreateDataHistogram(7, binx, 17, biny); SetDataValuesFromTextFile( fSettings.GetDataInput() ); ScaleData(1E-42); SetMapValuesFromText( fSettings.GetMapInput() ); SetCholDecompFromTextFile( fSettings.GetCovarInput(), 67); ScaleCovar(1E-16); StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCinc_XSec_2DEavq3_nu::FillEventVariables(FitEvent *event) { //******************************************************************** // Seperate MEC if (splitMEC_PN_NN) { int npr = 0; int nne = 0; for (UInt_t j = 0; j < event->Npart(); j++) { if ((event->PartInfo(j))->fIsAlive) continue; if (event->PartInfo(j)->fPID == 2212) npr++; else if (event->PartInfo(j)->fPID == 2112) nne++; } if (event->Mode == 2 and npr == 1 and nne == 1) { event->Mode = 2; Mode = 2; } else if (event->Mode == 2 and npr == 0 and nne == 2) { event->Mode = 3; Mode = 3; } } // Set Defaults double Eav = -999.9; double q3 = -999.9; // If muon found get kinematics FitParticle* muon = event->GetHMFSParticle(13); FitParticle* neutrino = event->GetNeutrinoIn(); if (muon && neutrino) { // Set Q from Muon TLorentzVector q = neutrino->fP - muon->fP; double q0 = (q.E()) / 1.E3; //double q3_true = (q.Vect().Mag())/1.E3; double thmu = muon->fP.Vect().Angle(neutrino->fP.Vect()); double pmu = muon->fP.Vect().Mag() / 1.E3; double emu = muon->fP.E() / 1.E3; double mmu = muon->fP.Mag() / 1.E3; // Get Enu Rec double enu_rec = emu + q0; // Set Q2 QE double q2qe = 2 * enu_rec * (emu - pmu * cos(thmu)) - mmu * mmu; // Calc Q3 from Q2QE and EnuTree q3 = sqrt(q2qe + q0 * q0); // Get Eav too Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event) / 1.E3; } // Set Hist Variables fXVar = q3; fYVar = Eav; return; } //******************************************************************** bool MINERvA_CCinc_XSec_2DEavq3_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCCincLowRecoil_MINERvA(event, EnuMin, EnuMax); } diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx index 8e0aab6..4586ae4 100644 --- a/src/Reweight/GENIEWeightEngine.cxx +++ b/src/Reweight/GENIEWeightEngine.cxx @@ -1,244 +1,244 @@ #include "GENIEWeightEngine.h" GENIEWeightEngine::GENIEWeightEngine(std::string name) { #ifdef __GENIE_ENABLED__ // Setup the NEUT Reweight engien fCalcName = name; LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl; // Create RW Engine suppressing cout StopTalking(); fGenieRW = new genie::rew::GReWeight(); // Get List of Vetos (Just for debugging) - std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fGenieRW_veto"); + std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fGenieRW_veto"); bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos; bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos; bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos; bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos; bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos; bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos; bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos; bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos; // Now actually add the RW Calcs if (xsec_ncel) fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL); if (xsec_ccqe){ fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE); // (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") ); } if (xsec_coh){ fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH()); // (dynamic_cast (fGenieRW->WghtCalc("xsec_coh"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") ); } if (xsec_nnres) fGenieRW->AdoptWghtCalc("xsec_nonresbkg", new genie::rew::GReWeightNonResonanceBkg); if (xsec_nudis) fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod); if (xsec_resdec) fGenieRW->AdoptWghtCalc("hadro_res_decay", new genie::rew::GReWeightResonanceDecay); if (xsec_fzone) fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); if (xsec_intra) fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); if (xsec_agky) fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY); if (xsec_qevec) fGenieRW->AdoptWghtCalc("xsec_ccqe_vec", new genie::rew::GReWeightNuXSecCCQEvec); #if __GENIE_VERSION__ >= 212 if (xsec_qeaxial) fGenieRW->AdoptWghtCalc("xsec_ccqe_axial", new genie::rew::GReWeightNuXSecCCQEaxial); #endif if (xsec_dis) fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS); if (xsec_nc) fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC); if (xsec_ccres){ #if __GENIE_VERSION__ < 213 fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES); #else fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES( FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default")); #endif } if (xsec_ncres) fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES); if (xsec_nucqe) fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM); GReWeightNuXSecCCQE * rwccqe = dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe")); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa); // Default to include shape and normalization changes for CCRES (can be changed downstream if desired) GReWeightNuXSecCCRES * rwccres = dynamic_cast (fGenieRW->WghtCalc("xsec_ccres")); std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode"); if (!marestype.compare("kModeNormAndMaMvShape")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); } else if (!marestype.compare("kModeMaMv")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); } else { THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype ); } // Default to include shape and normalization changes for NCRES (can be changed downstream if desired) GReWeightNuXSecNCRES * rwncres = dynamic_cast (fGenieRW->WghtCalc("xsec_ncres")); rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv); // Default to include shape and normalization changes for DIS (can be changed downstream if desired) GReWeightNuXSecDIS * rwdis = dynamic_cast (fGenieRW->WghtCalc("xsec_dis")); rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); // allow cout again StartTalking(); #else ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl; #endif }; void GENIEWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __GENIE_ENABLED__ // Get First enum int nuisenum = Reweight::ConvDial(name, kGENIE); // Setup Maps fEnumIndex[nuisenum];// = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fGENIESysts.push_back(rwsyst); // Initialize dial std::cout << "Registering " << singlename << " from " << name << std::endl; fGenieRW->Systematics().Init( fGENIESysts[index] ); // If Absolute if (fIsAbsTwk) { GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 ); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } #endif }; void GENIEWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::SetDialValue(std::string name, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::Reconfigure(bool silent) { #ifdef __GENIE_ENABLED__ // Hush now... if (silent) StopTalking(); // Reconf fGenieRW->Reconfigure(); fGenieRW->Print(); // Shout again if (silent) StartTalking(); #endif } double GENIEWeightEngine::CalcWeight(BaseFitEvt* evt) { double rw_weight = 1.0; #ifdef __GENIE_ENABLED__ // Skip Non GENIE if (evt->fType != kGENIE) return 1.0; // Make nom weight if (!evt) { THROW("evt not found : " << evt); } if (!(evt->genie_event)) { THROW("evt->genie_event not found!" << evt->genie_event); } if (!(evt->genie_event->event)) { THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event)); } if (!fGenieRW) { THROW("GENIE RW Not Found!" << fGenieRW); } rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event)); // std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl; #endif // Return rw_weight return rw_weight; } diff --git a/src/Reweight/NEUTWeightEngine.cxx b/src/Reweight/NEUTWeightEngine.cxx index 3bd0c2e..f9aaa78 100644 --- a/src/Reweight/NEUTWeightEngine.cxx +++ b/src/Reweight/NEUTWeightEngine.cxx @@ -1,186 +1,186 @@ #include "NEUTWeightEngine.h" NEUTWeightEngine::NEUTWeightEngine(std::string name) { #ifdef __NEUT_ENABLED__ // Setup the NEUT Reweight engien fCalcName = name; LOG(FIT) << "Setting up NEUT RW : " << fCalcName << std::endl; // Create RW Engine suppressing cout StopTalking(); fNeutRW = new neut::rew::NReWeight(); TDirectory* olddir = gDirectory; // get list of vetoed calc engines (just for debug really) - std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fNeutRW_veto"); + std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fNeutRW_veto"); bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_res = rw_engine_list.find("xsec_res") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool nucl_casc = rw_engine_list.find("nucl_casc") == std::string::npos; bool nucl_piless = rw_engine_list.find("nucl_piless") == std::string::npos; // Activate each calc engine if (xsec_ccqe) fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE); if (xsec_res) fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES); if (xsec_ccres) fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES); if (xsec_coh) fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH); if (xsec_dis) fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS); if (xsec_ncel) fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL); if (xsec_nc) fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC); if (xsec_ncres) fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES); if (nucl_casc) fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc); if (nucl_piless) fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless); fNeutRW->Reconfigure(); olddir->cd(); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); // allow cout again StartTalking(); #else THROW("NEUT RW NOT ENABLED!" ); #endif }; void NEUTWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __NEUT_ENABLED__ // Get First enum int nuisenum = Reweight::ConvDial(name, kNEUT); // Setup Maps fEnumIndex[nuisenum];// = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get Syst neut::rew::NSyst_t gensyst = NSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fNEUTSysts.push_back(gensyst); // Initialize dial LOG(FIT) << "Registering " << singlename << " dial." << std::endl; fNeutRW->Systematics().Init( fNEUTSysts[index] ); // If Absolute if (fIsAbsTwk) { NSystUncertainty::Instance()->SetUncertainty( fNEUTSysts[index], 1.0, 1.0 ); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } #endif } void NEUTWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __NEUT_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; std::cout << "Setting Dial Value = " << nuisenum << " " << i << " " << indices[i] << " " << fValues[indices[i]] << " Enum=" << fNEUTSysts[indices[i]] << std::endl; fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val); } #endif } void NEUTWeightEngine::SetDialValue(std::string name, double val) { #ifdef __NEUT_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; std::cout << "Setting Dial Value = " << name << " = " << i << " " << indices[i] << " " << fValues[indices[i]] << " Enum=" << fNEUTSysts[indices[i]] << std::endl; fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val); } #endif } void NEUTWeightEngine::Reconfigure(bool silent) { #ifdef __NEUT_ENABLED__ // Hush now... if (silent) StopTalking(); // Reconf fNeutRW->Reconfigure(); //if (LOG_LEVEL(DEB)){ fNeutRW->Print(); // } // Shout again if (silent) StartTalking(); #endif } double NEUTWeightEngine::CalcWeight(BaseFitEvt* evt) { double rw_weight = 1.0; #ifdef __NEUT_ENABLED__ // Skip Non NEUT if (evt->fType != kNEUT) return 1.0; // Hush now StopTalking(); // Fill NEUT Common blocks NEUTUtils::FillNeutCommons(evt->fNeutVect); // Call Weight calculation rw_weight = fNeutRW->CalcWeight(); // Speak Now StartTalking(); #endif // Return rw_weight return rw_weight; } diff --git a/src/Reweight/NIWGWeightEngine.cxx b/src/Reweight/NIWGWeightEngine.cxx index 94c53ec..14691b9 100644 --- a/src/Reweight/NIWGWeightEngine.cxx +++ b/src/Reweight/NIWGWeightEngine.cxx @@ -1,220 +1,220 @@ #include "NIWGWeightEngine.h" NIWGWeightEngine::NIWGWeightEngine(std::string name) { #ifdef __NIWG_ENABLED__ #ifdef __NEUT_ENABLED__ // Setup the NEUT Reweight engien fCalcName = name; LOG(FIT) << "Setting up NIWG RW : " << fCalcName << std::endl; // Create RW Engine suppressing cout StopTalking(); fNIWGRW = new niwg::rew::NIWGReWeight(); // Get List of Veto Calcs (For Debugging) std::string rw_engine_list = - FitPar::Config().GetParS("FitWeight.fNIWGRW_veto"); + FitPar::Config().GetParS("FitWeight_fNIWGRW_veto"); bool niwg_2012a = rw_engine_list.find("niwg_2012a") == std::string::npos; bool niwg_2014a = rw_engine_list.find("niwg_2014a") == std::string::npos; bool niwg_pimult = rw_engine_list.find("niwg_pimult") == std::string::npos; bool niwg_mec = rw_engine_list.find("niwg_mec") == std::string::npos; bool niwg_rpa = rw_engine_list.find("niwg_rpa") == std::string::npos; bool niwg_eff_rpa = rw_engine_list.find("niwg_eff_rpa") == std::string::npos; bool niwg_proton = rw_engine_list.find("niwg_protonFSIbug") == std::string::npos; bool niwg_hadron = rw_engine_list.find("niwg_HadronMultSwitch") == std::string::npos; // Add the RW Calcs if (niwg_2012a) fNIWGRW->AdoptWghtCalc("niwg_2012a", new niwg::rew::NIWGReWeight2012a); if (niwg_2014a) fNIWGRW->AdoptWghtCalc("niwg_2014a", new niwg::rew::NIWGReWeight2014a); if (niwg_pimult) fNIWGRW->AdoptWghtCalc("niwg_pimult", new niwg::rew::NIWGReWeightPiMult); if (niwg_mec) fNIWGRW->AdoptWghtCalc("niwg_mec", new niwg::rew::NIWGReWeightMEC); if (niwg_rpa) fNIWGRW->AdoptWghtCalc("niwg_rpa", new niwg::rew::NIWGReWeightRPA); if (niwg_eff_rpa) fNIWGRW->AdoptWghtCalc("niwg_eff_rpa", new niwg::rew::NIWGReWeightEffectiveRPA); if (niwg_proton) fNIWGRW->AdoptWghtCalc("niwg_protonFSIbug", new niwg::rew::NIWGReWeightProtonFSIbug); if (niwg_hadron) fNIWGRW->AdoptWghtCalc("niwg_HadronMultSwitch", new niwg::rew::NIWGReWeightHadronMultSwitch); fNIWGRW->Reconfigure(); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); // allow cout again StartTalking(); #else ERR(FTL) << "NIWG RW Enabled but NEUT RW is not!" << std::endl; #endif #else ERR(FTL) << "NIWG RW NOT ENABLED!" << std::endl; #endif }; void NIWGWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __NIWG_ENABLED__ #ifdef __NEUT_ENABLED__ // Get First enum int nuisenum = Reweight::ConvDial(name, kNIWG); // Setup Maps fEnumIndex[nuisenum];// = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW niwg::rew::NIWGSyst_t gensyst = niwg::rew::NIWGSyst::FromString(name); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fNIWGSysts.push_back(gensyst); // Initialize dial LOG(FIT) << "Registering " << singlename << " from " << name << std::endl; fNIWGRW->Systematics().Init( fNIWGSysts[index] ); // If Absolute if (fIsAbsTwk) { niwg::rew::NIWGSystUncertainty::Instance()->SetUncertainty( fNIWGSysts[index], 1.0, 1.0 ); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } #endif #endif } void NIWGWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __NIWG_ENABLED__ #ifdef __NEUT_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; std::cout << "Setting NIWG Dial Value " << i << " " << indices[i] << " " << fNIWGSysts[indices[i]] << std::endl; fNIWGRW->Systematics().Set( fNIWGSysts[indices[i]], val); } #endif #endif } void NIWGWeightEngine::SetDialValue(std::string name, double val) { #ifdef __NIWG_ENABLED__ #ifdef __NEUT_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fNIWGRW->Systematics().Set(fNIWGSysts[indices[i]], val); } #endif #endif } void NIWGWeightEngine::Reconfigure(bool silent) { #ifdef __NIWG_ENABLED__ #ifdef __NEUT_ENABLED__ // Hush now... if (silent) StopTalking(); // Reconf fNIWGRW->Reconfigure(); // Shout again if (silent) StartTalking(); #endif #endif } double NIWGWeightEngine::CalcWeight(BaseFitEvt* evt) { double rw_weight = 1.0; #ifdef __NEUT_ENABLED__ #ifdef __NIWG_ENABLED__ // Skip Non GENIE if (evt->fType != kNEUT) return 1.0; // Hush now StopTalking(); niwg::rew::NIWGEvent* niwg_event = NIWGWeightEngine::GetNIWGEventLocal(evt->fNeutVect); rw_weight *= fNIWGRW->CalcWeight(*niwg_event); delete niwg_event; // Speak Now StartTalking(); #endif #endif // Return rw_weight return rw_weight; } #ifdef __NEUT_ENABLED__ #ifdef __NIWG_ENABLED__ niwg::rew::NIWGEvent * NIWGWeightEngine::GetNIWGEventLocal(NeutVect* nvect) { niwg::rew::NIWGEvent * fDummyNIWGEvent = NULL; fDummyNIWGEvent = new niwg::rew::NIWGEvent(); fDummyNIWGEvent->detid = 1; // MiniBooNE (apply CCQE LowE variations) fDummyNIWGEvent->neutmode = nvect->Mode; fDummyNIWGEvent->targetA = nvect->TargetA; fDummyNIWGEvent->recenu_ccqe_sk = -1; if (nvect->Ibound == 0) fDummyNIWGEvent->targetA = 1; //RT: identifies as H, rather than O16 // Fill initial particle stack for (int ip = 0; ip < nvect->Npart(); ++ip) { niwg::rew::NIWGPartStack fDummyPartStack; fDummyPartStack.p = (nvect->PartInfo(ip)->fP) * 0.001; // Convert to GeV fDummyPartStack.pdg = nvect->PartInfo(ip)->fPID; fDummyPartStack.chase = nvect->PartInfo(ip)->fIsAlive; fDummyPartStack.parent = nvect->ParentIdx(ip) - 1; // WARNING: this needs to be tested with a NeutRoot file fDummyNIWGEvent->part_stack.push_back(fDummyPartStack); } fDummyNIWGEvent->CalcKinematics(); return fDummyNIWGEvent; } #endif #endif diff --git a/src/Reweight/NuWroWeightEngine.cxx b/src/Reweight/NuWroWeightEngine.cxx index 149e23a..c9dcc95 100644 --- a/src/Reweight/NuWroWeightEngine.cxx +++ b/src/Reweight/NuWroWeightEngine.cxx @@ -1,136 +1,136 @@ #include "NuWroWeightEngine.h" NuWroWeightEngine::NuWroWeightEngine(std::string name) { #ifdef __NUWRO_REWEIGHT_ENABLED__ // Setup the NEUT Reweight engien fCalcName = name; LOG(FIT) << "Setting up NuWro RW : " << fCalcName << std::endl; // Create RW Engine suppressing cout StopTalking(); // Create the engine fNuwroRW = new nuwro::rew::NuwroReWeight(); // Get List of Veto Calcs (For Debugging) std::string rw_engine_list = - FitPar::Config().GetParS("FitWeight.fNuwroRW_veto"); + FitPar::Config().GetParS("FitWeight_fNuwroRW_veto"); bool xsec_qel = rw_engine_list.find("nuwro_QEL") == std::string::npos; bool xsec_flag = rw_engine_list.find("nuwro_FlagNorm") == std::string::npos; bool xsec_res = rw_engine_list.find("nuwro_RES") == std::string::npos; // Add the RW Calcs if (xsec_qel) fNuwroRW->AdoptWghtCalc("nuwro_QEL", new nuwro::rew::NuwroReWeight_QEL); if (xsec_flag) fNuwroRW->AdoptWghtCalc("nuwro_FlagNorm", new nuwro::rew::NuwroReWeight_FlagNorm); if (xsec_res) fNuwroRW->AdoptWghtCalc("nuwro_RES", new nuwro::rew::NuwroReWeight_SPP); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); // allow cout again StartTalking(); #else ERR(FTL) << "NUWRO RW NOT ENABLED! " << std::endl; #endif }; void NuWroWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __NUWRO_REWEIGHT_ENABLED__ // Get RW Enum and name int nuisenum = Reweight::ConvDial(name, kNUWRO); // Initialise new vector fEnumIndex[nuisenum]; fNameIndex[name]; // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get Syst nuwro::rew::NuwroSyst_t gensyst = nuwro::rew::NuwroSyst::FromString(name); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fNUWROSysts.push_back(gensyst); // Initialise Dial LOG(FIT) << "Adding NuWro Syst " << fNUWROSysts[index] << std::endl; fNuwroRW->Systematics().Add(fNUWROSysts[index]); if (fIsAbsTwk) { nuwro::rew::NuwroSystUncertainty::Instance()->SetUncertainty( fNUWROSysts[index], 1.0, 1.0); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(name, startval); } #endif }; void NuWroWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __NUWRO_REWEIGHT_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fNuwroRW->Systematics().SetSystVal(fNUWROSysts[indices[i]], val); } #endif } void NuWroWeightEngine::SetDialValue(std::string name, double val) { #ifdef __NUWRO_REWEIGHT_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fNuwroRW->Systematics().SetSystVal(fNUWROSysts[indices[i]], val); } #endif } void NuWroWeightEngine::Reconfigure(bool silent) { #ifdef __NUWRO_REWEIGHT_ENABLED__ // Hush now... if (silent) StopTalking(); // Reconf fNuwroRW->Reconfigure(); // Shout again if (silent) StartTalking(); #endif } double NuWroWeightEngine::CalcWeight(BaseFitEvt* evt) { double rw_weight = 1.0; #ifdef __NUWRO_REWEIGHT_ENABLED__ // Skip Non GENIE if (evt->fType != kNUWRO) return 1.0; #ifdef __USE_NUWRO_SRW_EVENTS__ rw_weight = fNuwroRW->CalcWeight(evt->fNuwroSRWEvent, *evt->fNuwroParams); #else // Call Weight calculation rw_weight = fNuwroRW->CalcWeight(evt->fNuwroEvent); #endif #endif // Return rw_weight return rw_weight; } diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx index ee3e995..24cd30d 100755 --- a/src/Routines/MinimizerRoutines.cxx +++ b/src/Routines/MinimizerRoutines.cxx @@ -1,1506 +1,1506 @@ // 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 "MinimizerRoutines.h" #include "Simple_MH_Sampler.h" /* Constructor/Destructor */ //************************ void MinimizerRoutines::Init() { //************************ fInputFile = ""; fInputRootFile = NULL; fOutputFile = ""; fOutputRootFile = NULL; fCovar = NULL; fCovFree = NULL; fCorrel = NULL; fCorFree = NULL; fDecomp = NULL; fDecFree = NULL; fStrategy = "Migrad,FixAtLimBreak,Migrad"; fRoutines.clear(); fCardFile = ""; fFakeDataInput = ""; fSampleFCN = NULL; fMinimizer = NULL; fMinimizerFCN = NULL; fCallFunctor = NULL; fAllowedRoutines = ("Migrad,Simplex,Combined," "Brute,Fumili,ConjugateFR," "ConjugatePR,BFGS,BFGS2," "SteepDesc,GSLSimAn,FixAtLim,FixAtLimBreak," "Chi2Scan1D,Chi2Scan2D,Contours,ErrorBands," "DataToys,MCMC"); }; //************************************* MinimizerRoutines::~MinimizerRoutines(){ //************************************* }; /* Input Functions */ //************************************* MinimizerRoutines::MinimizerRoutines(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"); std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl; std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl; // FitPar::log_verb = verbocount; SETVERBOSITY(verbocount); // ERR_VERB(errorcount); // Minimizer Setup ======================================== fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); SetupMinimizerFromXML(); SetupCovariance(); SetupRWEngine(); SetupFCN(); return; }; //************************************* void MinimizerRoutines::SetupMinimizerFromXML() { //************************************* 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 state if none 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; fStartVals[normname] = samplenorm; 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 MinimizerRoutines::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 MinimizerRoutines::SetupFCN() { //************************************* LOG(FIT) << "Making the jointFCN" << std::endl; if (fSampleFCN) delete fSampleFCN; // fSampleFCN = new JointFCN(fCardFile, fOutputRootFile); fSampleFCN = new JointFCN(fOutputRootFile); SetFakeData(); fMinimizerFCN = new MinimizerFCN(fSampleFCN); fCallFunctor = new ROOT::Math::Functor(*fMinimizerFCN, fParams.size()); fSampleFCN->CreateIterationTree("fit_iterations", FitBase::GetRW()); return; } //****************************************** void MinimizerRoutines::SetupFitter(std::string routine) { //****************************************** // Make the fitter std::string fitclass = ""; std::string fittype = ""; bool UseMCMC = false; // Get correct types if (!routine.compare("Migrad")) { fitclass = "Minuit2"; fittype = "Migrad"; } else if (!routine.compare("Simplex")) { fitclass = "Minuit2"; fittype = "Simplex"; } else if (!routine.compare("Combined")) { fitclass = "Minuit2"; fittype = "Combined"; } else if (!routine.compare("Brute")) { fitclass = "Minuit2"; fittype = "Scan"; } else if (!routine.compare("Fumili")) { fitclass = "Minuit2"; fittype = "Fumili"; } else if (!routine.compare("ConjugateFR")) { fitclass = "GSLMultiMin"; fittype = "ConjugateFR"; } else if (!routine.compare("ConjugatePR")) { fitclass = "GSLMultiMin"; fittype = "ConjugatePR"; } else if (!routine.compare("BFGS")) { fitclass = "GSLMultiMin"; fittype = "BFGS"; } else if (!routine.compare("BFGS2")) { fitclass = "GSLMultiMin"; fittype = "BFGS2"; } else if (!routine.compare("SteepDesc")) { fitclass = "GSLMultiMin"; fittype = "SteepestDescent"; // } else if (!routine.compare("GSLMulti")) { fitclass = "GSLMultiFit"; // fittype = ""; // Doesn't work out of the box } else if (!routine.compare("GSLSimAn")) { fitclass = "GSLSimAn"; fittype = ""; } else if (!routine.compare("MCMC")) { UseMCMC = true; } // make minimizer if (fMinimizer) delete fMinimizer; if (UseMCMC) { fMinimizer = new Simple_MH_Sampler(); } else { fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype); } fMinimizer->SetMaxFunctionCalls( - FitPar::Config().GetParI("minimizer.maxcalls")); + FitPar::Config().GetParI("MAXCALLS")); if (!routine.compare("Brute")) { fMinimizer->SetMaxFunctionCalls(fParams.size() * fParams.size() * 4); fMinimizer->SetMaxIterations(fParams.size() * fParams.size() * 4); } fMinimizer->SetMaxIterations( - FitPar::Config().GetParI("minimizer.maxiterations")); - fMinimizer->SetTolerance(FitPar::Config().GetParD("minimizer.tolerance")); - fMinimizer->SetStrategy(FitPar::Config().GetParI("minimizer.strategy")); + FitPar::Config().GetParI("MAXITERATIONS")); + fMinimizer->SetTolerance(FitPar::Config().GetParD("TOLERANCE")); + fMinimizer->SetStrategy(FitPar::Config().GetParI("STRATEGY")); fMinimizer->SetFunction(*fCallFunctor); int ipar = 0; // Add Fit Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); bool fixed = true; double vstart, vstep, vlow, vhigh; vstart = vstep = vlow = vhigh = 0.0; if (fCurVals.find(syst) != fCurVals.end()) vstart = fCurVals.at(syst); if (fMinVals.find(syst) != fMinVals.end()) vlow = fMinVals.at(syst); if (fMaxVals.find(syst) != fMaxVals.end()) vhigh = fMaxVals.at(syst); if (fStepVals.find(syst) != fStepVals.end()) vstep = fStepVals.at(syst); if (fFixVals.find(syst) != fFixVals.end()) fixed = fFixVals.at(syst); // fix for errors if (vhigh == vlow) vhigh += 1.0; fMinimizer->SetVariable(ipar, syst, vstart, vstep); fMinimizer->SetVariableLimits(ipar, vlow, vhigh); if (fixed) { fMinimizer->FixVariable(ipar); LOG(FIT) << "Fixed Param: " << syst << std::endl; } else { LOG(FIT) << "Free Param: " << syst << " Start:" << vstart << " Range:" << vlow << " to " << vhigh << " Step:" << vstep << std::endl; } ipar++; } LOG(FIT) << "Setup Minimizer: " << fMinimizer->NDim() << "(NDim) " << fMinimizer->NFree() << "(NFree)" << std::endl; return; } //************************************* // Set fake data from user input void MinimizerRoutines::SetFakeData() { //************************************* // If the fake data input field (-d) isn't provided, return to caller if (fFakeDataInput.empty()) return; // If user specifies -d MC we set the data to the MC // User can also specify fake data parameters to reweight by doing // "fake_parameter" in input card file // "fake_parameter" gets read in ReadCard function (reads to fFakeVals) if (fFakeDataInput.compare("MC") == 0) { LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl; // fFakeVals get read in in ReadCard UpdateRWEngine(fFakeVals); // Reconfigure the reweight engine FitBase::GetRW()->Reconfigure(); // Reconfigure all the samples to the new reweight fSampleFCN->ReconfigureAllEvents(); // Feed on and set the fake-data in each measurement class fSampleFCN->SetFakeData("MC"); // Changed the reweight engine values back to the current values // So we start the fit at a different value than what we set the fake-data // to UpdateRWEngine(fCurVals); LOG(FIT) << "Set all data to fake MC predictions." << std::endl; } else { fSampleFCN->SetFakeData(fFakeDataInput); } return; } /* Fitting Functions */ //************************************* void MinimizerRoutines::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 MinimizerRoutines::Run() { //************************************* LOG(FIT) << "Running MinimizerRoutines : " << 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 MinimizerRoutines with no routines given!" << std::endl; throw; } for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); int fitstate = kFitUnfinished; LOG(FIT) << "Running Routine: " << routine << std::endl; // Try Routines if (routine.find("LowStat") != std::string::npos) LowStatRoutine(routine); else if (routine == "FixAtLim") FixAtLimit(); else if (routine == "FixAtLimBreak") fitstate = FixAtLimit(); else if (routine.find("ErrorBands") != std::string::npos) GenerateErrorBands(); else if (routine.find("DataToys") != std::string::npos) ThrowDataToys(); else if (!routine.compare("Chi2Scan1D")) Create1DScans(); else if (!routine.compare("Chi2Scan2D")) Chi2Scan2D(); else fitstate = RunFitRoutine(routine); // If ending early break here if (fitstate == kFitFinished || fitstate == kNoChange) { LOG(FIT) << "Ending fit routines loop." << std::endl; break; } } return; } //************************************* int MinimizerRoutines::RunFitRoutine(std::string routine) { //************************************* int endfits = kFitUnfinished; // set fitter at the current start values fOutputRootFile->cd(); SetupFitter(routine); // choose what to do with the minimizer depending on routine. if (!routine.compare("Migrad") or !routine.compare("Simplex") or !routine.compare("Combined") or !routine.compare("Brute") or !routine.compare("Fumili") or !routine.compare("ConjugateFR") or !routine.compare("ConjugatePR") or !routine.compare("BFGS") or !routine.compare("BFGS2") or !routine.compare("SteepDesc") or // !routine.compare("GSLMulti") or !routine.compare("GSLSimAn") or !routine.compare("MCMC")) { if (fMinimizer->NFree() > 0) { LOG(FIT) << fMinimizer->Minimize() << std::endl; GetMinimizerState(); } } // other otptions else if (!routine.compare("Contour")) { CreateContours(); } return endfits; } //************************************* void MinimizerRoutines::PrintState() { //************************************* LOG(FIT) << "------------" << std::endl; // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header LOG(FIT) << " # " << left << setw(maxcount) << "Parameter " << " = " << setw(10) << "Value" << " +- " << setw(10) << "Error" << " " << setw(8) << "(Units)" << " " << setw(10) << "Conv. Val" << " +- " << setw(10) << "Conv. Err" << " " << setw(8) << "(Units)" << std::endl; // Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); std::string typestr = FitBase::ConvDialType(fTypeVals[syst]); std::string curunits = "(sig.)"; double curval = fCurVals[syst]; double curerr = fErrorVals[syst]; if (fStateVals[syst].find("ABS") != std::string::npos) { curval = FitBase::RWSigmaToAbs(typestr, syst, curval); curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); curunits = "(Abs.)"; } else if (fStateVals[syst].find("FRAC") != std::string::npos) { curval = FitBase::RWSigmaToFrac(typestr, syst, curval); curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) - FitBase::RWSigmaToFrac(typestr, syst, 0.0)); curunits = "(Frac)"; } std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")"; double convval = FitBase::RWSigmaToAbs(typestr, syst, curval); double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); std::ostringstream curparstring; curparstring << " " << setw(3) << left << i << ". " << setw(maxcount) << syst << " = " << setw(10) << curval << " +- " << setw(10) << curerr << " " << setw(8) << curunits << " " << setw(10) << convval << " +- " << setw(10) << converr << " " << setw(8) << convunits; LOG(FIT) << curparstring.str() << std::endl; } LOG(FIT) << "------------" << std::endl; double like = fSampleFCN->GetLikelihood(); LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl; LOG(FIT) << "------------" << std::endl; } //************************************* void MinimizerRoutines::GetMinimizerState() { //************************************* LOG(FIT) << "Minimizer State: " << std::endl; // Get X and Err const double* values = fMinimizer->X(); const double* errors = fMinimizer->Errors(); // int ipar = 0; for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); fCurVals[syst] = values[i]; fErrorVals[syst] = errors[i]; } PrintState(); // Covar SetupCovariance(); if (fMinimizer->CovMatrixStatus() > 0) { // Fill Full Covar std::cout << "Filling covariance" << std::endl; for (int i = 0; i < fCovar->GetNbinsX(); i++) { for (int j = 0; j < fCovar->GetNbinsY(); j++) { fCovar->SetBinContent(i + 1, j + 1, fMinimizer->CovMatrix(i, j)); } } int freex = 0; int freey = 0; for (int i = 0; i < fCovar->GetNbinsX(); i++) { if (fMinimizer->IsFixedVariable(i)) continue; freey = 0; for (int j = 0; j < fCovar->GetNbinsY(); j++) { if (fMinimizer->IsFixedVariable(j)) continue; fCovFree->SetBinContent(freex + 1, freey + 1, fMinimizer->CovMatrix(i, j)); freey++; } freex++; } fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation"); fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition"); if (fMinimizer->NFree() > 0) { fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free"); fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free"); } } std::cout << "Got STATE" << std::endl; return; }; //************************************* void MinimizerRoutines::LowStatRoutine(std::string routine) { //************************************* LOG(FIT) << "Running Low Statistics Routine: " << routine << std::endl; - int lowstatsevents = FitPar::Config().GetParI("minimizer.lowstatevents"); - int maxevents = FitPar::Config().GetParI("input.maxevents"); + int lowstatsevents = FitPar::Config().GetParI("LOWSTATEVENTS"); + int maxevents = FitPar::Config().GetParI("MAXEVENTS"); int verbosity = FitPar::Config().GetParI("VERBOSITY"); std::string trueroutine = routine; std::string substring = "LowStat"; trueroutine.erase(trueroutine.find(substring), substring.length()); // Set MAX EVENTS=1000 - Config::SetPar("input.maxevents", lowstatsevents); + Config::SetPar("MAXEVENTS", lowstatsevents); Config::SetPar("VERBOSITY", 3); SetupFCN(); RunFitRoutine(trueroutine); - Config::SetPar("input.maxevents", maxevents); + Config::SetPar("MAXEVENTS", maxevents); SetupFCN(); Config::SetPar("VERBOSITY", verbosity); return; } //************************************* void MinimizerRoutines::Create1DScans() { //************************************* // 1D Scan Routine // Steps through all free parameters about nominal using the step size // Creates a graph for each free parameter // At the current point create a 1D Scan for all parametes (Uncorrelated) for (UInt_t i = 0; i < fParams.size(); i++) { if (fFixVals[fParams[i]]) continue; LOG(FIT) << "Running 1D Scan for " << fParams[i] << std::endl; fSampleFCN->CreateIterationTree(fParams[i] + "_scan1D_iterations", FitBase::GetRW()); double scanmiddlepoint = fCurVals[fParams[i]]; // Determine N points needed double limlow = fMinVals[fParams[i]]; double limhigh = fMaxVals[fParams[i]]; double step = fStepVals[fParams[i]]; int npoints = int(fabs(limhigh - limlow) / (step + 0.)); TH1D* contour = new TH1D(("Chi2Scan1D_" + fParams[i]).c_str(), ("Chi2Scan1D_" + fParams[i] + ";" + fParams[i]).c_str(), npoints, limlow, limhigh); // Fill bins for (int x = 0; x < contour->GetNbinsX(); x++) { // Set X Val fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1); // Run Eval double* vals = FitUtils::GetArrayFromMap(fParams, fCurVals); double chi2 = fSampleFCN->DoEval(vals); delete vals; // Fill Contour contour->SetBinContent(x + 1, chi2); } // Save contour contour->Write(); // Reset Parameter fCurVals[fParams[i]] = scanmiddlepoint; // Save TTree fSampleFCN->WriteIterationTree(); } return; } //************************************* void MinimizerRoutines::Chi2Scan2D() { //************************************* // Chi2 Scan 2D // Creates a 2D chi2 scan by stepping through all free parameters // Works for all pairwise combos of free parameters // Scan I for (UInt_t i = 0; i < fParams.size(); i++) { if (fFixVals[fParams[i]]) continue; // Scan J for (UInt_t j = 0; j < i; j++) { if (fFixVals[fParams[j]]) continue; fSampleFCN->CreateIterationTree( fParams[i] + "_" + fParams[j] + "_" + "scan2D_iterations", FitBase::GetRW()); double scanmid_i = fCurVals[fParams[i]]; double scanmid_j = fCurVals[fParams[j]]; double limlow_i = fMinVals[fParams[i]]; double limhigh_i = fMaxVals[fParams[i]]; double step_i = fStepVals[fParams[i]]; double limlow_j = fMinVals[fParams[j]]; double limhigh_j = fMaxVals[fParams[j]]; double step_j = fStepVals[fParams[j]]; int npoints_i = int(fabs(limhigh_i - limlow_i) / (step_i + 0.)) + 1; int npoints_j = int(fabs(limhigh_j - limlow_j) / (step_j + 0.)) + 1; TH2D* contour = new TH2D( ("Chi2Scan2D_" + fParams[i] + "_" + fParams[j]).c_str(), ("Chi2Scan2D_" + fParams[i] + "_" + fParams[j] + ";" + fParams[i] + ";" + fParams[j]) .c_str(), npoints_i, limlow_i, limhigh_i, npoints_j, limlow_j, limhigh_j); // Begin Scan LOG(FIT) << "Running scan for " << fParams[i] << " " << fParams[j] << std::endl; // Fill bins for (int x = 0; x < contour->GetNbinsX(); x++) { // Set X Val fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1); // Loop Y for (int y = 0; y < contour->GetNbinsY(); y++) { // Set Y Val fCurVals[fParams[j]] = contour->GetYaxis()->GetBinCenter(y + 1); // Run Eval double* vals = FitUtils::GetArrayFromMap(fParams, fCurVals); double chi2 = fSampleFCN->DoEval(vals); delete vals; // Fill Contour contour->SetBinContent(x + 1, y + 1, chi2); fCurVals[fParams[j]] = scanmid_j; } fCurVals[fParams[i]] = scanmid_i; fCurVals[fParams[j]] = scanmid_j; } // Save contour contour->Write(); // Save Iterations fSampleFCN->WriteIterationTree(); } } return; } //************************************* void MinimizerRoutines::CreateContours() { //************************************* // Use MINUIT for this if possible ERR(FTL) << " Contours not yet implemented as it is really slow!" << std::endl; throw; return; } //************************************* int MinimizerRoutines::FixAtLimit() { //************************************* bool fixedparam = false; for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); if (fFixVals[syst]) continue; double curVal = fCurVals.at(syst); double minVal = fMinVals.at(syst); double maxVal = fMinVals.at(syst); if (fabs(curVal - minVal) < 0.0001) { fCurVals[syst] = minVal; fFixVals[syst] = true; fixedparam = true; } if (fabs(maxVal - curVal) < 0.0001) { fCurVals[syst] = maxVal; fFixVals[syst] = true; fixedparam = true; } } if (!fixedparam) { LOG(FIT) << "No dials needed fixing!" << std::endl; return kNoChange; } else return kStateChange; } /* Write Functions */ //************************************* void MinimizerRoutines::SaveResults() { //************************************* fOutputRootFile->cd(); if (fMinimizer) { SetupCovariance(); SaveMinimizerState(); } SaveCurrentState(); } //************************************* void MinimizerRoutines::SaveMinimizerState() { //************************************* std::cout << "Saving Minimizer State" << std::endl; if (!fMinimizer) { ERR(FTL) << "Can't save minimizer state without min object" << std::endl; throw; } // Save main fit tree fSampleFCN->WriteIterationTree(); // Get Vals and Errors GetMinimizerState(); // Save tree with fit status std::vector nameVect; std::vector valVect; std::vector errVect; std::vector minVect; std::vector maxVect; std::vector startVect; std::vector endfixVect; std::vector startfixVect; // int NFREEPARS = fMinimizer->NFree(); int NPARS = fMinimizer->NDim(); int ipar = 0; // Dial Vals for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams.at(i); nameVect.push_back(name); valVect.push_back(fCurVals.at(name)); errVect.push_back(fErrorVals.at(name)); minVect.push_back(fMinVals.at(name)); maxVect.push_back(fMaxVals.at(name)); startVect.push_back(fStartVals.at(name)); endfixVect.push_back(fFixVals.at(name)); startfixVect.push_back(fStartFixVals.at(name)); ipar++; } int NFREE = fMinimizer->NFree(); int NDIM = fMinimizer->NDim(); double CHI2 = fSampleFCN->GetLikelihood(); int NBINS = fSampleFCN->GetNDOF(); int NDOF = NBINS - NFREE; // Write fit results TTree* fit_tree = new TTree("fit_result", "fit_result"); fit_tree->Branch("parameter_names", &nameVect); fit_tree->Branch("parameter_values", &valVect); fit_tree->Branch("parameter_errors", &errVect); fit_tree->Branch("parameter_min", &minVect); fit_tree->Branch("parameter_max", &maxVect); fit_tree->Branch("parameter_start", &startVect); fit_tree->Branch("parameter_fix", &endfixVect); fit_tree->Branch("parameter_startfix", &startfixVect); fit_tree->Branch("CHI2", &CHI2, "CHI2/D"); fit_tree->Branch("NDOF", &NDOF, "NDOF/I"); fit_tree->Branch("NBINS", &NBINS, "NBINS/I"); fit_tree->Branch("NDIM", &NDIM, "NDIM/I"); fit_tree->Branch("NFREE", &NFREE, "NFREE/I"); fit_tree->Fill(); fit_tree->Write(); // Make dial variables TH1D dialvar = TH1D("fit_dials", "fit_dials", NPARS, 0, NPARS); TH1D startvar = TH1D("start_dials", "start_dials", NPARS, 0, NPARS); TH1D minvar = TH1D("min_dials", "min_dials", NPARS, 0, NPARS); TH1D maxvar = TH1D("max_dials", "max_dials", NPARS, 0, NPARS); TH1D dialvarfree = TH1D("fit_dials_free", "fit_dials_free", NFREE, 0, NFREE); TH1D startvarfree = TH1D("start_dials_free", "start_dials_free", NFREE, 0, NFREE); TH1D minvarfree = TH1D("min_dials_free", "min_dials_free", NFREE, 0, NFREE); TH1D maxvarfree = TH1D("max_dials_free", "max_dials_free", NFREE, 0, NFREE); int freecount = 0; for (UInt_t i = 0; i < nameVect.size(); i++) { std::string name = nameVect.at(i); dialvar.SetBinContent(i + 1, valVect.at(i)); dialvar.SetBinError(i + 1, errVect.at(i)); dialvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); startvar.SetBinContent(i + 1, startVect.at(i)); startvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); minvar.SetBinContent(i + 1, minVect.at(i)); minvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); maxvar.SetBinContent(i + 1, maxVect.at(i)); maxvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); if (NFREE > 0) { if (!startfixVect.at(i)) { freecount++; dialvarfree.SetBinContent(freecount, valVect.at(i)); dialvarfree.SetBinError(freecount, errVect.at(i)); dialvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); startvarfree.SetBinContent(freecount, startVect.at(i)); startvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); minvarfree.SetBinContent(freecount, minVect.at(i)); minvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); maxvarfree.SetBinContent(freecount, maxVect.at(i)); maxvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); } } } // Save Dial Plots dialvar.Write(); startvar.Write(); minvar.Write(); maxvar.Write(); if (NFREE > 0) { dialvarfree.Write(); startvarfree.Write(); minvarfree.Write(); maxvarfree.Write(); } // Save fit_status plot TH1D statusplot = TH1D("fit_status", "fit_status", 8, 0, 8); std::string fit_labels[8] = {"status", "cov_status", "maxiter", "maxfunc", "iter", "func", "precision", "tolerance"}; double fit_vals[8]; fit_vals[0] = fMinimizer->Status() + 0.; fit_vals[1] = fMinimizer->CovMatrixStatus() + 0.; fit_vals[2] = fMinimizer->MaxIterations() + 0.; fit_vals[3] = fMinimizer->MaxFunctionCalls() + 0.; fit_vals[4] = fMinimizer->NIterations() + 0.; fit_vals[5] = fMinimizer->NCalls() + 0.; fit_vals[6] = fMinimizer->Precision() + 0.; fit_vals[7] = fMinimizer->Tolerance() + 0.; for (int i = 0; i < 8; i++) { statusplot.SetBinContent(i + 1, fit_vals[i]); statusplot.GetXaxis()->SetBinLabel(i + 1, fit_labels[i].c_str()); } statusplot.Write(); // Save Covars if (fCovar) fCovar->Write(); if (fCovFree) fCovFree->Write(); if (fCorrel) fCorrel->Write(); if (fCorFree) fCorFree->Write(); if (fDecomp) fDecomp->Write(); if (fDecFree) fDecFree->Write(); return; } //************************************* void MinimizerRoutines::SaveCurrentState(std::string subdir) { //************************************* LOG(FIT) << "Saving current full FCN predictions" << std::endl; // Setup DIRS TDirectory* curdir = gDirectory; if (!subdir.empty()) { TDirectory* newdir = (TDirectory*)gDirectory->mkdir(subdir.c_str()); newdir->cd(); } FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void MinimizerRoutines::SaveNominal() { //************************************* fOutputRootFile->cd(); LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl; FitBase::GetRW()->Reconfigure(); SaveCurrentState("nominal"); }; //************************************* void MinimizerRoutines::SavePrefit() { //************************************* fOutputRootFile->cd(); LOG(FIT) << "Saving Prefit Predictions" << std::endl; UpdateRWEngine(fStartVals); SaveCurrentState("prefit"); UpdateRWEngine(fCurVals); }; /* MISC Functions */ //************************************* int MinimizerRoutines::GetStatus() { //************************************* return 0; } //************************************* void MinimizerRoutines::SetupCovariance() { //************************************* // Remove covares if they exist if (fCovar) delete fCovar; if (fCovFree) delete fCovFree; if (fCorrel) delete fCorrel; if (fCorFree) delete fCorFree; if (fDecomp) delete fDecomp; if (fDecFree) delete fDecFree; LOG(FIT) << "Building covariance matrix.." << std::endl; int NFREE = 0; int NDIM = 0; // Get NFREE from min or from vals (for cases when doing throws) if (fMinimizer) { std::cout << "NFREE FROM MINIMIZER" << std::endl; NFREE = fMinimizer->NFree(); NDIM = fMinimizer->NDim(); } else { NDIM = fParams.size(); for (UInt_t i = 0; i < fParams.size(); i++) { std::cout << "Getting Param " << fParams[i] << std::endl; if (!fFixVals[fParams[i]]) NFREE++; } } if (NDIM == 0) return; LOG(FIT) << "NFREE == " << NFREE << std::endl; fCovar = new TH2D("covariance", "covariance", NDIM, 0, NDIM, NDIM, 0, NDIM); if (NFREE > 0) { fCovFree = new TH2D("covariance_free", "covariance_free", NFREE, 0, NFREE, NFREE, 0, NFREE); } else { fCovFree = NULL; } // Set Bin Labels int countall = 0; int countfree = 0; for (UInt_t i = 0; i < fParams.size(); i++) { std::cout << "Getting Param " << i << std::endl; std::cout << "ParamI = " << fParams[i] << std::endl; 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) { fCovFree->GetXaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); fCovFree->GetYaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); countfree++; } } std::cout << "Filling Matrices" << std::endl; fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation"); fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition"); if (NFREE > 0) { fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free"); fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free"); } else { fCorFree = NULL; fDecFree = NULL; } std::cout << " Set the covariance" << std::endl; return; }; //************************************* void MinimizerRoutines::ThrowCovariance(bool uniformly) { //************************************* std::vector rands; if (!fDecFree) { ERR(WRN) << "Trying to throw 0 free parameters" << std::endl; return; } // Generate Random Gaussians for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) { rands.push_back(gRandom->Gaus(0.0, 1.0)); } // Reset Thrown Values for (UInt_t i = 0; i < fParams.size(); i++) { fThrownVals[fParams[i]] = fCurVals[fParams[i]]; } // Loop and get decomp for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) { std::string parname = std::string(fDecFree->GetXaxis()->GetBinLabel(i + 1)); double mod = 0.0; if (!uniformly) { for (Int_t j = 0; j < fDecFree->GetNbinsY(); j++) { mod += rands[j] * fDecFree->GetBinContent(j + 1, i + 1); } } if (fCurVals.find(parname) != fCurVals.end()) { if (uniformly) fThrownVals[parname] = gRandom->Uniform(fMinVals[parname], fMaxVals[parname]); else { fThrownVals[parname] = fCurVals[parname] + mod; } } } // Check Limits for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams[i]; if (fFixVals[syst]) continue; if (fThrownVals[syst] < fMinVals[syst]) fThrownVals[syst] = fMinVals[syst]; if (fThrownVals[syst] > fMaxVals[syst]) fThrownVals[syst] = fMaxVals[syst]; } return; }; //************************************* void MinimizerRoutines::GenerateErrorBands() { //************************************* TDirectory* errorDIR = (TDirectory*)fOutputRootFile->mkdir("error_bands"); errorDIR->cd(); // Make a second file to store throws std::string tempFileName = fOutputFile; if (tempFileName.find(".root") != std::string::npos) tempFileName.erase(tempFileName.find(".root"), 5); tempFileName += ".throws.root"; TFile* tempfile = new TFile(tempFileName.c_str(), "RECREATE"); tempfile->cd(); int nthrows = FitPar::Config().GetParI("error_throws"); UpdateRWEngine(fCurVals); fSampleFCN->ReconfigureAllEvents(); TDirectory* nominal = (TDirectory*)tempfile->mkdir("nominal"); nominal->cd(); fSampleFCN->Write(); TDirectory* outnominal = (TDirectory*)fOutputRootFile->mkdir("nominal_throw"); outnominal->cd(); fSampleFCN->Write(); errorDIR->cd(); 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"); bool uniformly = FitPar::Config().GetParB("error_uniform"); // Run Throws and save for (Int_t i = 0; i < nthrows; i++) { TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i", i)); throwfolder->cd(); // Generate Random Parameter Throw ThrowCovariance(uniformly); // Run Eval double* vals = FitUtils::GetArrayFromMap(fParams, fThrownVals); chi2 = fSampleFCN->DoEval(vals); delete vals; // Save the FCN fSampleFCN->Write(); parameterTree->Fill(); } errorDIR->cd(); fDecFree->Write(); fCovFree->Write(); parameterTree->Write(); delete parameterTree; // 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; TH1D* baseplot = (TH1D*)key->ReadObj(); std::string plotname = std::string(baseplot->GetName()); int nbins = baseplot->GetNbinsX() * 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)); } for (Int_t i = 0; i < nthrows; i++) { TH1* newplot = (TH1*)tempfile->Get(Form(("throw_%i/" + plotname).c_str(), i)); 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(); delete newplot; } errorDIR->cd(); for (Int_t j = 0; j < nbins; j++) { if (!uniformly) { baseplot->SetBinError(j + 1, tprof->GetBinError(j + 1)); } else { baseplot->SetBinContent(j + 1, (binlowest[j] + binhighest[j]) / 2.0); baseplot->SetBinError(j + 1, (binhighest[j] - binlowest[j]) / 2.0); } } errorDIR->cd(); baseplot->Write(); tprof->Write(); bintree->Write(); delete baseplot; delete tprof; delete bintree; delete[] bincontents; } return; }; void MinimizerRoutines::ThrowDataToys() { LOG(FIT) << "Generating Toy Data Throws" << std::endl; int verb = Config::GetParI("VERBOSITY"); SETVERBOSITY(FIT); int nthrows = FitPar::Config().GetParI("NToyThrows"); double maxlike = -1.0; double minlike = -1.0; std::vector values; for (int i = 0; i < 1.E4; i++) { fSampleFCN->ThrowDataToy(); double like = fSampleFCN->GetLikelihood(); values.push_back(like); if (maxlike == -1.0 or like > maxlike) maxlike = like; if (minlike == -1.0 or like < minlike) minlike = like; } SETVERBOSITY(verb); // Fill Histogram TH1D* likes = new TH1D("toydatalikelihood", "toydatalikelihood", int(sqrt(nthrows)), minlike, maxlike); for (size_t i = 0; i < values.size(); i++) { likes->Fill(values[i]); } // Save to file LOG(FIT) << "Writing toy data throws" << std::endl; fOutputRootFile->cd(); likes->Write(); } diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx index 2648470..1af6e2e 100644 --- a/src/Statistical/StatUtils.cxx +++ b/src/Statistical/StatUtils.cxx @@ -1,1303 +1,1303 @@ // 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 "TH1D.h" #include "StatUtils.h" #include "NuisConfig.h" #include "GeneralUtils.h" //******************************************************************* Double_t StatUtils::GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask) { //******************************************************************* Double_t Chi2 = 0.0; TH1D* calc_data = (TH1D*)data->Clone(); TH1D* calc_mc = (TH1D*)mc->Clone(); // Add MC Error to data if required - if (FitPar::Config().GetParB("statutils.addmcerror")) { + if (FitPar::Config().GetParB("addmcerror")) { for (int i = 0; i < calc_data->GetNbinsX(); i++) { double dterr = calc_data->GetBinError(i + 1); double mcerr = calc_mc->GetBinError(i + 1); if (dterr > 0.0) { calc_data->SetBinError(i + 1, sqrt(dterr * dterr + mcerr * mcerr)); } } } // Apply masking if required if (mask) { calc_data = ApplyHistogramMasking(data, mask); calc_mc = ApplyHistogramMasking(mc, mask); } // Iterate over bins in X for (int i = 0; i < calc_data->GetNbinsX(); i++) { // Ignore bins with zero data or zero bin error if (calc_data->GetBinError(i + 1) <= 0.0 || calc_data->GetBinContent(i + 1) == 0.0) continue; // Take mc data difference double diff = calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1); double err = calc_data->GetBinError(i + 1); Chi2 += (diff * diff) / (err * err); } // cleanup delete calc_data; delete calc_mc; return Chi2; }; //******************************************************************* Double_t StatUtils::GetChi2FromDiag(TH2D* data, TH2D* mc, TH2I* map, TH2I* mask) { //******************************************************************* // Generate a simple map if (!map) map = GenerateMap(data); // Convert to 1D Histograms TH1D* data_1D = MapToTH1D(data, map); TH1D* mc_1D = MapToTH1D(mc, map); TH1I* mask_1D = MapToMask(mask, map); // Calculate 1D chi2 from 1D Plots Double_t Chi2 = StatUtils:: GetChi2FromDiag(data_1D, mc_1D, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; }; //******************************************************************* Double_t StatUtils::GetChi2FromCov(TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask, double data_scale, double covar_scale) { //******************************************************************* Double_t Chi2 = 0.0; TMatrixDSym* calc_cov = (TMatrixDSym*) invcov->Clone(); TH1D* calc_data = (TH1D*) data->Clone(); TH1D* calc_mc = (TH1D*) mc->Clone(); // If a mask if applied we need to apply it before the matrix is inverted if (mask) { calc_cov = ApplyInvertedMatrixMasking(invcov, mask); calc_data = ApplyHistogramMasking(data, mask); calc_mc = ApplyHistogramMasking(mc, mask); } // Add MC Error to data if required - if (FitPar::Config().GetParB("statutils.addmcerror")) { + if (FitPar::Config().GetParB("addmcerror")) { // Make temp cov TMatrixDSym* newcov = StatUtils::GetInvert(calc_cov); // Add MC err to diag for (int i = 0; i < calc_data->GetNbinsX(); i++) { double mcerr = calc_mc->GetBinError(i + 1) * sqrt(covar_scale); double oldval = (*newcov)(i, i); std::cout << "Adding cov stat " << mcerr*mcerr << " to " << (*newcov)(i,i) << std::endl; (*newcov)(i, i) = oldval + mcerr * mcerr; } // Reset the calc_cov to new invert delete calc_cov; calc_cov = GetInvert(newcov); // Delete the tempcov delete newcov; } calc_data->Scale(data_scale); calc_mc ->Scale(data_scale); (*calc_cov) *= covar_scale; // iterate over bins in X (i,j) for (int i = 0; i < calc_data->GetNbinsX(); i++) { for (int j = 0; j < calc_data->GetNbinsX(); j++) { if (calc_data->GetBinContent(i + 1) != 0 || calc_mc->GetBinContent(i + 1) != 0) { LOG(DEB) << "i j = " << i << " " << j << std::endl; LOG(DEB) << "Calc_data mc i = " << calc_data->GetBinContent(i + 1) << " " << calc_mc->GetBinContent(i + 1) << " Dif = " << ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) << std::endl; LOG(DEB) << "Calc_data mc i = " << calc_data->GetBinContent(j + 1) << " " << calc_mc->GetBinContent(j + 1) << " Dif = " << ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1) ) << std::endl; LOG(DEB) << "Covar = " << (*calc_cov)(i, j) << std::endl; LOG(DEB) << "Cont chi2 = " \ << ( ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) \ * (*calc_cov)(i, j) \ * ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1))) << " " << Chi2 << std::endl; Chi2 += ( ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) \ * (*calc_cov)(i, j) \ * ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1) ) ); } else { LOG(DEB) << "Error on bin (i,j) = (" << i << "," << j << ")" << std::endl; LOG(DEB) << "data->GetBinContent(i+1) = " << calc_data->GetBinContent(i + 1) << std::endl; LOG(DEB) << "mc->GetBinContent(i+1) = " << calc_mc->GetBinContent(i + 1) << std::endl; LOG(DEB) << "Adding zero to chi2 instead of dying horrifically " << std::endl; Chi2 += 0.; } } } // Cleanup delete calc_cov; delete calc_data; delete calc_mc; return Chi2; } //******************************************************************* Double_t StatUtils::GetChi2FromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map, TH2I* mask) { //******************************************************************* // Generate a simple map if (!map) { map = StatUtils::GenerateMap(data); } // Convert to 1D Histograms TH1D* data_1D = MapToTH1D(data, map); TH1D* mc_1D = MapToTH1D(mc, map); TH1I* mask_1D = MapToMask(mask, map); // Calculate 1D chi2 from 1D Plots Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; } //******************************************************************* Double_t StatUtils::GetChi2FromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask) { //******************************************************************* Double_t Chi2 = 0.0; TMatrixDSym* calc_cov = (TMatrixDSym*) cov->Clone(); TH1D* calc_data = (TH1D*) data->Clone(); TH1D* calc_mc = (TH1D*) mc->Clone(); // If a mask if applied we need to apply it before the matrix is inverted if (mask) { calc_cov = StatUtils::ApplyMatrixMasking(cov, mask); calc_data = StatUtils::ApplyHistogramMasking(data, mask); calc_mc = StatUtils::ApplyHistogramMasking(mc, mask); } // Decompose matrix TDecompSVD LU = TDecompSVD((*calc_cov)); LU.Decompose(); TMatrixDSym* cov_U = new TMatrixDSym(calc_data->GetNbinsX(), LU .GetU().GetMatrixArray(), ""); TVectorD* cov_S = new TVectorD( LU.GetSig() ); // Apply basis rotation before adding up chi2 Double_t rotated_difference = 0.0; for (int i = 0; i < calc_data->GetNbinsX(); i++) { rotated_difference = 0.0; // Rotate basis of Data - MC for (int j = 0; j < calc_data->GetNbinsY(); j++) rotated_difference += ( calc_data->GetBinContent(j + 1) - calc_mc ->GetBinContent(j + 1) ) * (*cov_U)(j, i) ; // Divide by rotated error cov_S Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i); } // Cleanup delete calc_cov; delete calc_data; delete calc_mc; delete cov_U; delete cov_S; return Chi2; } //******************************************************************* Double_t StatUtils::GetChi2FromSVD( TH2D* data, TH2D* mc, TMatrixDSym* cov, TH2I* map, TH2I* mask) { //******************************************************************* // Generate a simple map if (!map) map = StatUtils::GenerateMap(data); // Convert to 1D Histograms TH1D* data_1D = MapToTH1D(data, map); TH1D* mc_1D = MapToTH1D(mc, map); TH1I* mask_1D = MapToMask(mask, map); // Calculate from 1D Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; } //******************************************************************* double StatUtils::GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask) { //******************************************************************* // If just an event rate, for chi2 just use Poission Likelihood to calculate the chi2 component double chi2 = 0.0; TH1D* calc_data = (TH1D*)data->Clone(); TH1D* calc_mc = (TH1D*)mc->Clone(); // Apply masking if required if (mask) { calc_data = ApplyHistogramMasking(data, mask); calc_mc = ApplyHistogramMasking(mc, mask); } // Iterate over bins in X for (int i = 0; i < calc_data->GetNbinsX(); i++) { double dt = calc_data->GetBinContent(i + 1); double mc = calc_mc->GetBinContent(i + 1); if (mc <= 0) continue; if (dt <= 0) { // Only add difference chi2 += 2 * (mc - dt); } else { // Do the chi2 for Poisson distributions chi2 += 2 * (mc - dt + (dt * log(dt / mc))); } /* LOG(REC)<<"Evt Chi2 cont = "<Clone(); // If a mask is provided we need to apply it before getting NDOF if (mask) { calc_hist = StatUtils::ApplyHistogramMasking(hist, mask); } // NDOF is defined as total number of bins with non-zero errors Int_t NDOF = 0; for (int i = 0; i < calc_hist->GetNbinsX(); i++) { if (calc_hist->GetBinError(i + 1) > 0.0) NDOF++; } delete calc_hist; return NDOF; }; //******************************************************************* Int_t StatUtils::GetNDOF(TH2D* hist, TH2I* map, TH2I* mask) { //******************************************************************* Int_t NDOF = 0; if (!map) map = StatUtils::GenerateMap(hist); for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (mask->GetBinContent(i + 1, j + 1)) continue; if (map->GetBinContent(i + 1, j + 1) <= 0) continue; NDOF++; } } return NDOF; }; //******************************************************************* TH1D* StatUtils::ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag, TH1I* mask) { //******************************************************************* TH1D* calc_hist = (TH1D*) hist->Clone( (std::string(hist->GetName()) + "_THROW" ).c_str() ); TMatrixDSym* calc_cov = (TMatrixDSym*) cov->Clone(); Double_t correl_val = 0.0; // If a mask if applied we need to apply it before the matrix is decomposed if (mask) { calc_cov = ApplyMatrixMasking(cov, mask); calc_hist = ApplyHistogramMasking(calc_hist, mask); } // If a covariance is provided we need a preset random vector and a decomp std::vector rand_val; TMatrixDSym* decomp_cov; if (cov) { for (int i = 0; i < hist->GetNbinsX(); i++) { rand_val.push_back(gRandom->Gaus(0.0, 1.0)); } // Decomp the matrix decomp_cov = StatUtils::GetDecomp(calc_cov); } // iterate over bins for (int i = 0; i < hist->GetNbinsX(); i++) { // By Default the errors on the histogram are thrown uncorrelated to the other errors // if (throwdiag) { // calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + \ // gRandom->Gaus(0.0, 1.0) * calc_hist->GetBinError(i + 1)) ); // } // If a covariance is provided that is also thrown if (cov) { correl_val = 0.0; for (int j = 0; j < hist->GetNbinsX(); j++) { correl_val += rand_val[j] * (*decomp_cov)(j, i) ; } calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38)); } } delete calc_cov; delete decomp_cov; // return this new thrown data return calc_hist; }; //******************************************************************* TH2D* StatUtils::ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map, bool throwdiag, TH2I* mask) { //******************************************************************* // PLACEHOLDER!!!!!!!!! // Currently no support for throwing 2D Histograms from a covariance (void) hist; (void) cov; (void) map; (void) throwdiag; (void) mask; // /todo // Sort maps if required // Throw the covariance for a 1D plot // Unmap back to 2D Histogram return hist; } //******************************************************************* TH1D* StatUtils::ApplyHistogramMasking(TH1D* hist, TH1I* mask) { //******************************************************************* if (!mask) return ( (TH1D*)hist->Clone() ); // This masking is only sufficient for chi2 calculations, and will have dodgy bin edges. // Get New Bin Count Int_t NBins = 0; for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) continue; NBins++; } // Make new hist std::string newmaskname = std::string(hist->GetName()) + "_MSKD"; TH1D* calc_hist = new TH1D( newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins); // fill new hist int binindex = 0; for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) { LOG(REC) << "Applying mask to bin " << i + 1 << " " << hist->GetName() << std::endl; continue; } calc_hist->SetBinContent(binindex + 1, hist->GetBinContent(i + 1)); calc_hist->SetBinError(binindex + 1, hist->GetBinError(i + 1)); binindex++; } return calc_hist; }; //******************************************************************* TH2D* StatUtils::ApplyHistogramMasking(TH2D* hist, TH2I* mask) { //******************************************************************* TH2D* newhist = (TH2D*) hist->Clone(); if (!mask) return newhist; for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (mask->GetBinContent(i + 1, j + 1) > 0) { newhist->SetBinContent(i + 1, j + 1, 0.0); newhist->SetBinContent(i + 1, j + 1, 0.0); } } } return newhist; } //******************************************************************* TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH1I* mask) { //******************************************************************* if (!mask) return (TMatrixDSym*)(mat->Clone()); // Get New Bin Count Int_t NBins = 0; for (int i = 0; i < mask->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) continue; NBins++; } // make new matrix TMatrixDSym* calc_mat = new TMatrixDSym(NBins); int col, row; // Need to mask out bins in the current matrix row = 0; for (int i = 0; i < mask->GetNbinsX(); i++) { col = 0; // skip if masked if (mask->GetBinContent(i + 1) > 0.5) continue; for (int j = 0; j < mask->GetNbinsX(); j++) { // skip if masked if (mask->GetBinContent(j + 1) > 0.5) continue; (*calc_mat)(row, col) = (*mat)(i, j); col++; } row++; } return calc_mat; }; //******************************************************************* TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1I* mask_1D = StatUtils::MapToMask(mask, map); TMatrixDSym* newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D); delete mask_1D; return newmat; } //******************************************************************* TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH1I* mask) { //******************************************************************* TMatrixDSym* new_mat = GetInvert(mat); TMatrixDSym* masked_mat = ApplyMatrixMasking(new_mat, mask); TMatrixDSym* inverted_mat = GetInvert(masked_mat); delete masked_mat; delete new_mat; return inverted_mat; }; //******************************************************************* TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1I* mask_1D = StatUtils::MapToMask(mask, map); TMatrixDSym* newmat = ApplyInvertedMatrixMasking(mat, mask_1D); delete mask_1D; return newmat; } //******************************************************************* TMatrixDSym* StatUtils::GetInvert(TMatrixDSym* mat) { //******************************************************************* TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone(); // Check for diagonal bool non_diagonal = false; for (int i = 0; i < new_mat->GetNrows(); i++) { for (int j = 0; j < new_mat->GetNrows(); j++) { if (i == j) continue; if ((*new_mat)(i, j) != 0.0) { non_diagonal = true; break; } } } // If diag, just flip the diag if (!non_diagonal or new_mat->GetNrows() == 1) { for (int i = 0; i < new_mat->GetNrows(); i++) { if ((*new_mat)(i, i) != 0.0) (*new_mat)(i, i) = 1.0 / (*new_mat)(i, i); else (*new_mat)(i, i) = 0.0; } return new_mat; } // Invert full matrix TDecompSVD LU = TDecompSVD((*new_mat)); new_mat = new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), ""); return new_mat; } //******************************************************************* TMatrixDSym* StatUtils::GetDecomp(TMatrixDSym* mat) { //******************************************************************* TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone(); int nrows = new_mat->GetNrows(); // Check for diagonal bool diagonal = true; for (int i = 0; i < nrows; i++) { for (int j = 0; j < nrows; j++) { if (i == j) continue; if ((*new_mat)(i, j) != 0.0) { diagonal = false; break; } } } // If diag, just flip the diag if (diagonal or nrows == 1) { for (int i = 0; i < nrows; i++) { if ((*new_mat)(i, i) > 0.0) (*new_mat)(i, i) = sqrt((*new_mat)(i, i)); else (*new_mat)(i, i) = 0.0; } return new_mat; } TDecompChol LU = TDecompChol(*new_mat); LU.Decompose(); delete new_mat; TMatrixDSym* dec_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), ""); return dec_mat; } //******************************************************************* void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH1D* hist, double norm) { //******************************************************************* if (!mat) mat = MakeDiagonalCovarMatrix(hist); int nbins = mat->GetNrows(); TMatrixDSym* new_mat = new TMatrixDSym(nbins); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { double valx = hist->GetBinContent(i + 1) * 1E38; double valy = hist->GetBinContent(j + 1) * 1E38; (*new_mat)(i, j) = (*mat)(i, j) + norm * norm * valx * valy; } } // Swap the two delete mat; mat = new_mat; return; }; //******************************************************************* void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, TH2I* map ) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1D* data_1D = MapToTH1D(data, map); StatUtils::ForceNormIntoCovar(mat, data_1D, norm); delete data_1D; return; } //******************************************************************* TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH1D* data, double scaleF) { //******************************************************************* TMatrixDSym* newmat = new TMatrixDSym(data->GetNbinsX()); for (int i = 0; i < data->GetNbinsX(); i++) { (*newmat)(i, i) = data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF; } return newmat; } //******************************************************************* TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH2D* data, TH2I* map, double scaleF) { //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1D* data_1D = MapToTH1D(data, map); return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF); }; //******************************************************************* void StatUtils::SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, double scale) { //******************************************************************* // Check if (cov->GetNrows() != data->GetNbinsX()) { ERR(WRN) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << std::endl; } // Set bin errors form cov diag for (int i = 0; i < data->GetNbinsX(); i++) { data->SetBinError(i + 1, sqrt((*cov)(i, i)) * scale ); } return; } //******************************************************************* void StatUtils::SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map, double scale) { //******************************************************************* // Create map if required if (!map) map = StatUtils::GenerateMap(data); // Set Bin Errors from cov diag int count = 0; for (int i = 0; i < data->GetNbinsX(); i++) { for (int j = 0; j < data->GetNbinsY(); j++) { if (data->GetBinContent(i + 1, j + 1) == 0.0) continue; count = map->GetBinContent(i + 1, j + 1) - 1; data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale ); } } return; } TMatrixDSym* StatUtils::ExtractShapeOnlyCovar(TMatrixDSym* full_covar, TH1* data_hist, double data_scale){ int nbins = full_covar->GetNrows(); TMatrixDSym* shape_covar = new TMatrixDSym(nbins); // Check nobody is being silly if (data_hist->GetNbinsX() != nbins){ ERR(WRN) << "Inconsistent matrix and data histogram passed to StatUtils::ExtractShapeOnlyCovar!" << std::endl; ERR(WRN) << "data_hist has " << data_hist->GetNbinsX() << " matrix has " << nbins << std::endl; int err_bins = data_hist->GetNbinsX(); if (nbins > err_bins) err_bins = nbins; for (int i = 0; i < err_bins; ++i){ ERR(WRN) << "Matrix diag. = " << (*full_covar)(i, i) << " data = " << data_hist->GetBinContent(i+1) << std::endl; } return NULL; } double total_data = 0; double total_covar = 0; // Initial loop to calculate some constants for (int i = 0; i < nbins; ++i) { total_data += data_hist->GetBinContent(i+1)*data_scale; for (int j = 0; j < nbins; ++j) { total_covar += (*full_covar)(i,j); } } if (total_data == 0 || total_covar == 0){ ERR(WRN) << "Stupid matrix or data histogram passed to StatUtils::ExtractShapeOnlyCovar! Ignoring..." << std::endl; return NULL; } LOG(SAM) << "Norm error = " << sqrt(total_covar)/total_data << std::endl; // Now loop over and calculate the shape-only matrix for (int i = 0; i < nbins; ++i) { double data_i = data_hist->GetBinContent(i+1)*data_scale; for (int j = 0; j < nbins; ++j) { double data_j = data_hist->GetBinContent(j+1)*data_scale; double norm_term = data_i*data_j*total_covar/total_data/total_data; double mix_sum1 = 0; double mix_sum2 = 0; for (int k = 0; k < nbins; ++k){ mix_sum1 += (*full_covar)(k,j); mix_sum2 += (*full_covar)(i,k); } double mix_term1 = data_i*(mix_sum1/total_data - total_covar*data_j/total_data/total_data); double mix_term2 = data_j*(mix_sum2/total_data - total_covar*data_i/total_data/total_data); (*shape_covar)(i, j) = (*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term; } } return shape_covar; } //******************************************************************* TH2I* StatUtils::GenerateMap(TH2D* hist) { //******************************************************************* std::string maptitle = std::string(hist->GetName()) + "_MAP"; TH2I* map = new TH2I( maptitle.c_str(), maptitle.c_str(), hist->GetNbinsX(), 0, hist->GetNbinsX(), hist->GetNbinsY(), 0, hist->GetNbinsY()); Int_t index = 1; for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (hist->GetBinContent(i + 1, j + 1) > 0 && hist->GetBinError(i + 1, j + 1) > 0) { map->SetBinContent(i + 1, j + 1, index); index++; } else { map->SetBinContent(i + 1, j + 1, 0); } } } return map; } //******************************************************************* TH1D* StatUtils::MapToTH1D(TH2D* hist, TH2I* map) { //******************************************************************* if (!hist) return NULL; // Get N bins for 1D plot Int_t Nbins = map->GetMaximum(); std::string name1D = std::string(hist->GetName()) + "_1D"; // Make new 1D Hist TH1D* newhist = new TH1D(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins); // map bin contents for (int i = 0; i < map->GetNbinsX(); i++) { for (int j = 0; j < map->GetNbinsY(); j++) { if (map->GetBinContent(i + 1, j + 1) == 0) continue; newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1)); newhist->SetBinError(map->GetBinContent(i + 1, j + 1), hist->GetBinError(i + 1, j + 1)); } } // return return newhist; } //******************************************************************* TH1I* StatUtils::MapToMask(TH2I* hist, TH2I* map) { //******************************************************************* TH1I* newhist = NULL; if (!hist) return newhist; // Get N bins for 1D plot Int_t Nbins = map->GetMaximum(); std::string name1D = std::string(hist->GetName()) + "_1D"; // Make new 1D Hist newhist = new TH1I(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins); // map bin contents for (int i = 0; i < map->GetNbinsX(); i++) { for (int j = 0; j < map->GetNbinsY(); j++) { if (map->GetBinContent(i + 1, j + 1) == 0) continue; newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1)); } } // return return newhist; } TMatrixDSym* StatUtils::GetCovarFromCorrel(TMatrixDSym* correl, TH1D* data) { int nbins = correl->GetNrows(); TMatrixDSym* covar = new TMatrixDSym(nbins); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { (*covar)(i, j) = (*correl)(i, j) * data->GetBinError(i + 1) * data->GetBinError(j + 1); } } return covar; } //******************************************************************* TMatrixD* StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx, int dimy) { //******************************************************************* // Determine dim if (dimx == -1 and dimy == -1) { std::string line; std::ifstream covar(covfile.c_str(), std::ifstream::in); int row = 0; while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " "entries on this line: " << row << std::endl; } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { column++; if (column > dimx) dimx = column; } row++; if (row > dimy) dimy = row; } } // Or assume symmetric if (dimx != -1 and dimy == -1) { dimy = dimx; } assert(dimy != -1 && " matrix dimy not set."); // Make new matrix TMatrixD* mat = new TMatrixD(dimx, dimy); std::string line; std::ifstream covar(covfile.c_str(), std::ifstream::in); int row = 0; while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " "entries on this line: " << row << std::endl; } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { // Check Rows //assert(row > mat->GetNrows() && " covar rows doesn't match matrix rows."); //assert(column > mat->GetNcols() && " covar cols doesn't match matrix cols."); // Fill Matrix (*mat)(row, column) = (*iter); column++; } row++; } return mat; } //******************************************************************* TMatrixD* StatUtils::GetMatrixFromRootFile(std::string covfile, std::string histname) { //******************************************************************* std::string inputfile = covfile + ";" + histname; std::vector splitfile = GeneralUtils::ParseToStr(inputfile, ";"); if (splitfile.size() < 2) { ERR(FTL) << "No object name given!" << std::endl; throw; } // Get file TFile* tempfile = new TFile(splitfile[0].c_str(), "READ"); // Get Object TObject* obj = tempfile->Get(splitfile[1].c_str()); if (!obj) { ERR(FTL) << "Object " << splitfile[1] << " doesn't exist!" << std::endl; throw; } // Try casting TMatrixD* mat = dynamic_cast(obj); if (mat) { TMatrixD* newmat = (TMatrixD*)mat->Clone(); delete mat; tempfile->Close(); return newmat; } TMatrixDSym* matsym = dynamic_cast(obj); if (matsym) { TMatrixD* newmat = new TMatrixD(matsym->GetNrows(), matsym->GetNrows()); for (int i = 0; i < matsym->GetNrows(); i++) { for (int j = 0; j < matsym->GetNrows(); j++) { (*newmat)(i, j) = (*matsym)(i, j); } } delete matsym; tempfile->Close(); return newmat; } TH2D* mathist = dynamic_cast(obj); if (mathist) { TMatrixD* newmat = new TMatrixD(mathist->GetNbinsX(), mathist->GetNbinsX()); for (int i = 0; i < mathist->GetNbinsX(); i++) { for (int j = 0; j < mathist->GetNbinsX(); j++) { (*newmat)(i, j) = mathist->GetBinContent(i + 1, j + 1); } } delete mathist; tempfile->Close(); return newmat; } return NULL; } //******************************************************************* TMatrixDSym* StatUtils::GetCovarFromTextFile(std::string covfile, int dim){ //******************************************************************* // Delete TempMat TMatrixD* tempmat = GetMatrixFromTextFile(covfile, dim, dim); // Make a symmetric covariance TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows()); for (int i = 0; i < tempmat->GetNrows(); i++){ for (int j = 0; j < tempmat->GetNrows(); j++){ (*newmat)(i,j) = (*tempmat)(i,j); } } delete tempmat; return newmat; } //******************************************************************* TMatrixDSym* StatUtils::GetCovarFromRootFile(std::string covfile, std::string histname){ //******************************************************************* TMatrixD* tempmat = GetMatrixFromRootFile(covfile, histname); TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows()); for (int i = 0; i < tempmat->GetNrows(); i++){ for (int j = 0; j < tempmat->GetNrows(); j++){ (*newmat)(i,j) = (*tempmat)(i,j); } } delete tempmat; return newmat; }