diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx index d9e09cb..3bc668a 100644 --- a/src/FitBase/Measurement1D.cxx +++ b/src/FitBase/Measurement1D.cxx @@ -1,1940 +1,1966 @@ // Copyright 2016 L. Pickering, P. Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This ile 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 "Measurement1D.h" //******************************************************************** Measurement1D::Measurement1D(void) { //******************************************************************** // XSec Scalings fScaleFactor = -1.0; fCurrentNorm = 1.0; // Histograms fDataHist = NULL; fDataTrue = NULL; fMCHist = NULL; fMCFine = NULL; fMCWeighted = NULL; fMaskHist = NULL; // Covar covar = NULL; fFullCovar = NULL; fShapeCovar = NULL; fCovar = NULL; fInvert = NULL; fDecomp = NULL; fResidualHist = NULL; fChi2LessBinHist = NULL; // Fake Data fFakeDataInput = ""; fFakeDataFile = NULL; // Options fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsNoWidth = false; fIsDifXSec = false; fIsEnu1D = false; + fIsWriting = false; // Inputs fInput = NULL; fRW = NULL; // Extra Histograms fMCHist_Modes = NULL; } //******************************************************************** Measurement1D::~Measurement1D(void) { //******************************************************************** if (fDataHist) delete fDataHist; if (fDataTrue) delete fDataTrue; if (fMCHist) delete fMCHist; if (fMCFine) delete fMCFine; if (fMCWeighted) delete fMCWeighted; if (fMaskHist) delete fMaskHist; if (covar) delete covar; if (fFullCovar) delete fFullCovar; if (fShapeCovar) delete fShapeCovar; if (fCovar) delete fCovar; if (fInvert) delete fInvert; if (fDecomp) delete fDecomp; delete fResidualHist; delete fChi2LessBinHist; } //******************************************************************** void Measurement1D::FinaliseSampleSettings() { //******************************************************************** MeasurementBase::FinaliseSampleSettings(); // Setup naming + renaming fName = fSettings.GetName(); fSettings.SetS("originalname", fName); if (fSettings.Has("rename")) { fName = fSettings.GetS("rename"); fSettings.SetS("name", fName); } // Setup all other options NUIS_LOG(SAM, "Finalising Sample Settings: " << fName); if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods."); } if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; NUIS_LOG(SAM, "::" << fName << "::"); NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " << "not flux averaged!"); } if (fIsEnu1D && fIsRawEvents) { NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ERR(FTL, "I live in " << __FILE__ << ":" << __LINE__); throw; } if (!fRW) fRW = FitBase::GetRW(); if (!fInput and !fIsJoint) SetupInputs(fSettings.GetS("input")); // Setup options SetFitOptions(fDefaultTypes); // defaults SetFitOptions(fSettings.GetS("type")); // user specified EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min")); EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max")); if (fAddNormPen) { if (fNormError <= 0.0) { NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!"); NUIS_ERR(FTL, "If you want to use it please add fNormError=VAL"); throw; } } } //******************************************************************** void Measurement1D::CreateDataHistogram(int dimx, double *binx) { //******************************************************************** if (fDataHist) delete fDataHist; fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(), dimx, binx); } //******************************************************************** void Measurement1D::SetDataFromTextFile(std::string datafile) { //******************************************************************** NUIS_LOG(SAM, "Reading data from text file: " << datafile); fDataHist = PlotUtils::GetTH1DFromFile( datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles()); } //******************************************************************** void Measurement1D::SetDataFromRootFile(std::string datafile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading data from root file: " << datafile << ";" << histname); fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname); fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str()); return; }; //******************************************************************** void Measurement1D::SetEmptyData() { //******************************************************************** fDataHist = new TH1D("EMPTY_DATA", "EMPTY_DATA", 1, 0.0, 1.0); } //******************************************************************** void Measurement1D::SetPoissonErrors() { //******************************************************************** if (!fDataHist) { NUIS_ERR(FTL, "Need a data hist to setup possion errors! "); NUIS_ERR(FTL, "Setup Data First!"); throw; } for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) { fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1))); } } //******************************************************************** void Measurement1D::SetCovarFromDiagonal(TH1D *data) { //******************************************************************** if (!data and fDataHist) { data = fDataHist; } if (data) { NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName()); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } else { NUIS_ABORT("No data input provided to set diagonal covar from!"); } // if (!fIsDiag) { // ERR(FTL) << "SetCovarMatrixFromDiag called for measurement " // << "that is not set as diagonal." ); // throw; // } } //******************************************************************** void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } NUIS_LOG(SAM, "Reading covariance from text file: " << covfile); fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarFromMultipleTextFiles(std::string covfiles, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } std::vector covList = GeneralUtils::ParseToStr(covfiles, ";"); fFullCovar = new TMatrixDSym(dim); for (uint i = 0; i < covList.size(); ++i) { NUIS_LOG(SAM, "Reading covariance from text file: " << covList[i]); TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim); (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading covariance from text file: " << covfile << ";" << histname); fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile); covar = StatUtils::GetCovarFromTextFile(covfile, dim); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";" << histname); covar = StatUtils::GetCovarFromRootFile(covfile, histname); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) dim = fDataHist->GetNbinsX(); NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << dim); TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim); if (!fDataHist) { NUIS_ABORT("Trying to set correlations from text file but there is no " "data to build it from. \n" << "In constructor make sure data is set before " "SetCorrelationFromTextFile is called. \n"); } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(dim); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } std::vector corrList = GeneralUtils::ParseToStr(corrfiles, ";"); fFullCovar = new TMatrixDSym(dim); for (uint i = 0; i < corrList.size(); ++i) { NUIS_LOG(SAM, "Reading covariance from text file: " << corrList[i]); TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << histname); TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname); if (!fDataHist) { NUIS_ABORT("Trying to set correlations from text file but there is no " "data to build it from. \n" << "In constructor make sure data is set before " "SetCorrelationFromTextFile is called. \n"); } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX()); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile); TMatrixD *temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";" << histname); TMatrixD *temp = StatUtils::GetMatrixFromRootFile(covfile, histname); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } void Measurement1D::SetShapeCovar() { // Return if this is missing any pre-requisites if (!fFullCovar) return; if (!fDataHist) return; // Also return if it's bloody stupid under the circumstances if (fIsDiag) return; fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist); return; } //******************************************************************** void Measurement1D::ScaleData(double scale) { //******************************************************************** fDataHist->Scale(scale); } //******************************************************************** void Measurement1D::ScaleDataErrors(double scale) { //******************************************************************** for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinError(i + 1, fDataHist->GetBinError(i + 1) * scale); } } //******************************************************************** void Measurement1D::ScaleCovar(double scale) { //******************************************************************** (*fFullCovar) *= scale; (*covar) *= 1.0 / scale; (*fDecomp) *= sqrt(scale); } //******************************************************************** void Measurement1D::SetBinMask(std::string maskfile) { //******************************************************************** if (!fIsMask) return; NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile); // Create a mask histogram with dim of data int nbins = fDataHist->GetNbinsX(); fMaskHist = new TH1I((fSettings.GetName() + "_BINMASK").c_str(), (fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbins, 0, nbins); std::string line; std::ifstream mask(maskfile.c_str(), std::ifstream::in); if (!mask.is_open()) { NUIS_ABORT("Cannot find mask file."); } while (std::getline(mask >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToInt(line, " "); // Skip lines with poorly formatted lines if (entries.size() < 2) { NUIS_LOG(WRN, "Measurement1D::SetBinMask(), couldn't parse line: " << line); continue; } // The first index should be the bin number, the second should be the mask // value. int val = 0; if (entries[1] > 0) val = 1; fMaskHist->SetBinContent(entries[0], val); } // Apply masking by setting masked data bins to zero PlotUtils::MaskBins(fDataHist, fMaskHist); return; } //******************************************************************** void Measurement1D::FinaliseMeasurement() { //******************************************************************** NUIS_LOG(SAM, "Finalising Measurement: " << fName); if (fSettings.GetB("onlymc")) { if (fDataHist) delete fDataHist; fDataHist = new TH1D("empty_data", "empty_data", 1, 0.0, 1.0); } // Make sure data is setup if (!fDataHist) { NUIS_ABORT("No data has been setup inside " << fName << " constructor!"); } // Make sure covariances are setup if (!fFullCovar) { fIsDiag = true; SetCovarFromDiagonal(fDataHist); + } else if (fIsDiag) { // Have covariance but also set Diag + NUIS_LOG(SAM, "Have full covariance for sample " + << GetName() + << " but only using diagonal elements for likelihood"); + size_t nbins = fFullCovar->GetNcols(); + for (int i = 0; i < nbins; ++i) { + for (int j = 0; j < nbins; ++j) { + if (i != j) { + (*fFullCovar)[i][j] = 0; + } + } + } + delete covar; + covar = NULL; + delete fDecomp; + fDecomp = NULL; } if (!covar) { covar = StatUtils::GetInvert(fFullCovar); } if (!fDecomp) { fDecomp = StatUtils::GetDecomp(fFullCovar); } // Push the diagonals of fFullCovar onto the data histogram // Comment this out until the covariance/data scaling is consistent! StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38); // If shape only, set covar and fDecomp using the shape-only matrix (if set) if (fIsShape && fShapeCovar && FitPar::Config().GetParB("UseShapeCovar")) { if (covar) delete covar; covar = StatUtils::GetInvert(fShapeCovar); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetDecomp(fFullCovar); } // Setup fMCHist from data fMCHist = (TH1D *)fDataHist->Clone(); fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(), (fSettings.GetFullTitles()).c_str()); fMCHist->Reset(); // Setup fMCFine fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX() * 8, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1)); fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(), (fSettings.GetFullTitles()).c_str()); fMCFine->Reset(); // Setup MC Stat fMCStat = (TH1D *)fMCHist->Clone(); fMCStat->Reset(); // Search drawopts for possible types to include by default std::string drawopts = FitPar::Config().GetParS("drawopts"); if (drawopts.find("MODES") != std::string::npos) { fMCHist_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(), ("True Channels"), fMCHist); SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); } - // Setup bin masks using sample name - if (fIsMask) { + if (fSettings.Has("maskfile") && fSettings.Has("maskhist")) { + fMaskHist = dynamic_cast(PlotUtils::GetTH1FromRootFile( + fSettings.GetS("maskfile"), fSettings.GetS("maskhist"))); + fIsMask = bool(fMaskHist); + NUIS_LOG(SAM, "Loaded mask histogram: " << fSettings.GetS("maskhist") + << " from " + << fSettings.GetS("maskfile")); + } else if (fIsMask) { // Setup bin masks using sample name std::string curname = fName; std::string origname = fSettings.GetS("originalname"); // Check rename.mask std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask"); // Check origname.mask if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask"); // Check database if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask"; } // Setup Bin Mask SetBinMask(maskloc); } if (fScaleFactor < 0) { NUIS_ERR(FTL, "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__); NUIS_ERR(FTL, "fScaleFactor = " << fScaleFactor); NUIS_ERR(FTL, "EXITING"); throw; } // Create and fill Weighted Histogram if (!fMCWeighted) { fMCWeighted = (TH1D *)fMCHist->Clone(); fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(), (fName + "_MCWGHTS" + fPlotTitles).c_str()); fMCWeighted->GetYaxis()->SetTitle("Weighted Events"); } } //******************************************************************** void Measurement1D::SetFitOptions(std::string opt) { //******************************************************************** // Do nothing if default given if (opt == "DEFAULT") return; // CHECK Conflicting Fit Options std::vector fit_option_allow = GeneralUtils::ParseToStr(fAllowedTypes, "/"); for (UInt_t i = 0; i < fit_option_allow.size(); i++) { std::vector fit_option_section = GeneralUtils::ParseToStr(fit_option_allow.at(i), ","); bool found_option = false; for (UInt_t j = 0; j < fit_option_section.size(); j++) { std::string av_opt = fit_option_section.at(j); if (!found_option and opt.find(av_opt) != std::string::npos) { found_option = true; } else if (found_option and opt.find(av_opt) != std::string::npos) { NUIS_ABORT( "ERROR: Conflicting fit options provided: " << opt << std::endl << "Conflicting group = " << fit_option_section.at(i) << std::endl << "You should only supply one of these options in card file."); } } } // Check all options are allowed std::vector fit_options_input = GeneralUtils::ParseToStr(opt, "/"); for (UInt_t i = 0; i < fit_options_input.size(); i++) { if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) { NUIS_ERR(WRN, "ERROR: Fit Option '" << fit_options_input.at(i) << "' Provided is not allowed for this measurement."); NUIS_ERR(WRN, "Fit Options should be provided as a '/' seperated list " "(e.g. FREE/DIAG/NORM)"); NUIS_ABORT("Available options for " << fName << " are '" << fAllowedTypes << "'"); } } // Set TYPE fFitType = opt; // FIX,SHAPE,FREE if (opt.find("FIX") != std::string::npos) { fIsFree = fIsShape = false; fIsFix = true; } else if (opt.find("SHAPE") != std::string::npos) { fIsFree = fIsFix = false; fIsShape = true; } else if (opt.find("FREE") != std::string::npos) { fIsFix = fIsShape = false; fIsFree = true; } // DIAG,FULL (or default to full) if (opt.find("DIAG") != std::string::npos) { fIsDiag = true; fIsFull = false; } else if (opt.find("FULL") != std::string::npos) { fIsDiag = false; fIsFull = true; } // CHI2/LL (OTHERS?) if (opt.find("LOG") != std::string::npos) { fIsChi2 = false; NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!"); NUIS_ERR(FTL, "Try to use a chi2!"); throw; } else { fIsChi2 = true; } // EXTRAS if (opt.find("RAW") != std::string::npos) fIsRawEvents = true; if (opt.find("NOWIDTH") != std::string::npos) fIsNoWidth = true; if (opt.find("DIF") != std::string::npos) fIsDifXSec = true; if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true; if (opt.find("NORM") != std::string::npos) fAddNormPen = true; if (opt.find("MASK") != std::string::npos) fIsMask = true; return; }; //******************************************************************** void Measurement1D::SetSmearingMatrix(std::string smearfile, int truedim, int recodim) { //******************************************************************** // The smearing matrix describes the migration from true bins (rows) to reco // bins (columns) // Counter over the true bins! int row = 0; std::string line; std::ifstream smear(smearfile.c_str(), std::ifstream::in); // Note that the smearing matrix may be rectangular. fSmearMatrix = new TMatrixD(truedim, recodim); if (smear.is_open()) { NUIS_LOG(SAM, "Reading smearing matrix from file: " << smearfile); } else { NUIS_ABORT("Smearing matrix provided is incorrect: " << smearfile); } while (std::getline(smear >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*fSmearMatrix)(row, column) = (*iter) / 100.; // Convert to fraction from // percentage (this may not be // general enough) column++; } row++; } return; } //******************************************************************** void Measurement1D::ApplySmearingMatrix() { //******************************************************************** if (!fSmearMatrix) { NUIS_ERR(WRN, fName << ": attempted to apply smearing matrix, but none was set"); return; } TH1D *unsmeared = (TH1D *)fMCHist->Clone(); TH1D *smeared = (TH1D *)fMCHist->Clone(); smeared->Reset(); // Loop over reconstructed bins // true = row; reco = column for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) { // Sum up the constributions from all true bins double rBinVal = 0; // Loop over true bins for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) { rBinVal += (*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1); } smeared->SetBinContent(rbin + 1, rBinVal); } fMCHist = (TH1D *)smeared->Clone(); return; } /* Reconfigure LOOP */ //******************************************************************** void Measurement1D::ResetAll() { //******************************************************************** fMCHist->Reset(); fMCFine->Reset(); fMCStat->Reset(); return; }; //******************************************************************** void Measurement1D::FillHistograms() { //******************************************************************** if (Signal) { NUIS_LOG(DEB, "Fill MCHist: " << fXVar << ", " << Weight); fMCHist->Fill(fXVar, Weight); fMCFine->Fill(fXVar, Weight); fMCStat->Fill(fXVar, 1.0); if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight); } return; }; //******************************************************************** void Measurement1D::ScaleEvents() { //******************************************************************** // Fill MCWeighted; // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1)); // fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1)); // } // Setup Stat ratios for MC and MC Fine double *statratio = new double[fMCHist->GetNbinsX()]; for (int i = 0; i < fMCHist->GetNbinsX(); i++) { if (fMCHist->GetBinContent(i + 1) != 0) { statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1); } else { statratio[i] = 0.0; } } double *statratiofine = new double[fMCFine->GetNbinsX()]; for (int i = 0; i < fMCFine->GetNbinsX(); i++) { if (fMCFine->GetBinContent(i + 1) != 0) { statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1); } else { statratiofine[i] = 0.0; } } // Scaling for raw event rates if (fIsRawEvents) { double datamcratio = fDataHist->Integral() / fMCHist->Integral(); fMCHist->Scale(datamcratio); fMCFine->Scale(datamcratio); if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio); // Scaling for XSec as function of Enu } else if (fIsEnu1D) { PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); if (fMCHist_Modes) { // Loop over the modes fMCHist_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); // PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(), // GetEventHistogram(), fScaleFactor, // fNEvents); } } else if (fIsNoWidth) { fMCHist->Scale(fScaleFactor); fMCFine->Scale(fScaleFactor); if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor); // Any other differential scaling } else { fMCHist->Scale(fScaleFactor, "width"); fMCFine->Scale(fScaleFactor, "width"); if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width"); } // Proper error scaling - ROOT Freaks out with xsec weights sometimes for (int i = 0; i < fMCStat->GetNbinsX(); i++) { fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]); } for (int i = 0; i < fMCFine->GetNbinsX(); i++) { fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]); } // Clean up delete[] statratio; delete[] statratiofine; return; }; //******************************************************************** void Measurement1D::ApplyNormScale(double norm) { //******************************************************************** fCurrentNorm = norm; fMCHist->Scale(1.0 / norm); fMCFine->Scale(1.0 / norm); return; }; /* Statistic Functions - Outsources to StatUtils */ //******************************************************************** int Measurement1D::GetNDOF() { //******************************************************************** int ndof = fDataHist->GetNbinsX(); if (fMaskHist and fIsMask) ndof -= fMaskHist->Integral(); return ndof; } //******************************************************************** double Measurement1D::GetLikelihood() { //******************************************************************** // If this is for a ratio, there is no data histogram to compare to! if (fNoData || !fDataHist) return 0.; // Apply Masking to MC if Required. if (fIsMask and fMaskHist) { PlotUtils::MaskBins(fMCHist, fMaskHist); } // Sort Shape Scaling double scaleF = 0.0; // TODO Include !fIsRawEvents if (fIsShape) { if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) { scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") / fMCHist->Integral(1, fMCHist->GetNbinsX(), "width"); fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); } } // Likelihood Calculation double stat = 0.; if (fIsChi2) { if (fIsRawEvents) { stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist); } else if (fIsDiag) { stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist); } else if (!fIsDiag and !fIsRawEvents) { stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist); stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist, 1, - 1E76, fResidualHist); - if (fChi2LessBinHist) { - TH1I *binmask = new TH1I("mask", "", fDataHist->GetNbinsX(), 0, - fDataHist->GetNbinsX()); - binmask->SetDirectory(NULL); + 1E76, fIsWriting ? fResidualHist : NULL); + if (fChi2LessBinHist && fIsWriting) { for (int xi = 0; xi < fDataHist->GetNbinsX(); ++xi) { - binmask->Reset(); + TH1I *binmask = fMaskHist + ? static_cast(fMaskHist->Clone("mask")) + : new TH1I("mask", "", fDataHist->GetNbinsX(), 0, + fDataHist->GetNbinsX()); + binmask->SetDirectory(NULL); binmask->SetBinContent(xi + 1, 1); fChi2LessBinHist->SetBinContent( xi + 1, StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, binmask)); + delete binmask; } - delete binmask; } } } // Sort Penalty Terms if (fAddNormPen) { double penalty = (1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError); stat += penalty; } // Return to normal scaling if (fIsShape) { // and !FitPar::Config().GetParB("saveshapescaling")) { fMCHist->Scale(1. / scaleF); fMCFine->Scale(1. / scaleF); } fLikelihood = stat; return stat; } /* Fake Data Functions */ //******************************************************************** void Measurement1D::SetFakeDataValues(std::string fakeOption) { //******************************************************************** // Setup original/datatrue TH1D *tempdata = (TH1D *)fDataHist->Clone(); if (!fIsFakeData) { fIsFakeData = true; // Make a copy of the original data histogram. if (!fDataOrig) fDataOrig = (TH1D *)fDataHist->Clone((fName + "_data_original").c_str()); } else { ResetFakeData(); } // Setup Inputs fFakeDataInput = fakeOption; NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput); // From MC if (fFakeDataInput.compare("MC") == 0) { fDataHist = (TH1D *)fMCHist->Clone((fName + "_MC").c_str()); // Fake File } else { if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ"); fDataHist = (TH1D *)fFakeDataFile->Get((fName + "_MC").c_str()); } // Setup Data Hist fDataHist->SetNameTitle((fName + "_FAKE").c_str(), (fName + fPlotTitles).c_str()); // Replace Data True if (fDataTrue) delete fDataTrue; fDataTrue = (TH1D *)fDataHist->Clone(); fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(), (fName + fPlotTitles).c_str()); // Make a new covariance for fake data hist. int nbins = fDataHist->GetNbinsX(); double alpha_i = 0.0; double alpha_j = 0.0; for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1); alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1); (*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j); } } // Setup Covariances if (covar) delete covar; covar = StatUtils::GetInvert(fFullCovar); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetInvert(fFullCovar); delete tempdata; return; }; //******************************************************************** void Measurement1D::ResetFakeData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D *)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str()); } } //******************************************************************** void Measurement1D::ResetData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D *)fDataOrig->Clone((fSettings.GetName() + "_data").c_str()); } fIsFakeData = false; } //******************************************************************** void Measurement1D::ThrowCovariance() { //******************************************************************** // Take a fDecomposition and use it to throw the current dataset. // Requires fDataTrue also be set incase used repeatedly. if (!fDataTrue) fDataTrue = (TH1D *)fDataHist->Clone(); if (fDataHist) delete fDataHist; fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); return; }; //******************************************************************** void Measurement1D::ThrowDataToy() { //******************************************************************** if (!fDataTrue) fDataTrue = (TH1D *)fDataHist->Clone(); if (fMCHist) delete fMCHist; fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); } /* Access Functions */ //******************************************************************** TH1D *Measurement1D::GetMCHistogram() { //******************************************************************** if (!fMCHist) return fMCHist; std::ostringstream chi2; chi2 << std::setprecision(5) << this->GetLikelihood(); int linecolor = kRed; int linestyle = 1; int linewidth = 1; int fillcolor = 0; int fillstyle = 1001; // if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor"); // if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle"); // if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth"); // if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor"); // if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle"); fMCHist->SetTitle(chi2.str().c_str()); fMCHist->SetLineColor(linecolor); fMCHist->SetLineStyle(linestyle); fMCHist->SetLineWidth(linewidth); fMCHist->SetFillColor(fillcolor); fMCHist->SetFillStyle(fillstyle); return fMCHist; }; //******************************************************************** TH1D *Measurement1D::GetDataHistogram() { //******************************************************************** if (!fDataHist) return fDataHist; int datacolor = kBlack; int datastyle = 1; int datawidth = 1; // if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor"); // if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle"); // if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth"); fDataHist->SetLineColor(datacolor); fDataHist->SetLineWidth(datawidth); fDataHist->SetMarkerStyle(datastyle); return fDataHist; }; /* Write Functions */ // Save all the histograms at once //******************************************************************** void Measurement1D::Write(std::string drawOpt) { //******************************************************************** // Get Draw Options drawOpt = FitPar::Config().GetParS("drawopts"); // Write Settigns if (drawOpt.find("SETTINGS") != std::string::npos) { fSettings.Set("#chi^{2}", fLikelihood); fSettings.Set("NDOF", this->GetNDOF()); fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF()); fSettings.Write(); } // Write Data/MC if (drawOpt.find("DATA") != std::string::npos) GetDataList().at(0)->Write(); if (drawOpt.find("MC") != std::string::npos) { GetMCList().at(0)->Write(); if ((fEvtRateScaleFactor != 0xdeadbeef) && GetMCList().at(0)) { TH1D *PredictedEvtRate = static_cast(GetMCList().at(0)->Clone()); PredictedEvtRate->Scale(fEvtRateScaleFactor); PredictedEvtRate->GetYaxis()->SetTitle("Predicted event rate"); PredictedEvtRate->Write(); } } // Write Fine Histogram if (drawOpt.find("FINE") != std::string::npos) GetFineList().at(0)->Write(); // Write Weighted Histogram if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted) fMCWeighted->Write(); // Save Flux/Evt if no event manager if (!FitPar::Config().GetParB("EventManager")) { if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram()) GetFluxHistogram()->Write(); if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram()) GetEventHistogram()->Write(); if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram()) GetXSecHistogram()->Write(); } // Write Mask if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) { fMaskHist->Write(); } // Write Covariances if (drawOpt.find("COV") != std::string::npos && fFullCovar) { PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName())->Write(); } if (drawOpt.find("INVCOV") != std::string::npos && covar) { PlotUtils::GetInvCovarPlot(covar, fSettings.GetName())->Write(); } if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) { PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName())->Write(); } // // Likelihood residual plots // if (drawOpt.find("RESIDUAL") != std::string::npos) { // WriteResidualPlots(); // } // Ratio and Shape Plots if (drawOpt.find("RATIO") != std::string::npos) { WriteRatioPlot(); } if (drawOpt.find("SHAPE") != std::string::npos) { WriteShapePlot(); if (drawOpt.find("RATIO") != std::string::npos) WriteShapeRatioPlot(); } // // RATIO // if (drawOpt.find("CANVMC") != std::string::npos) { // TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist); // c1->Write(); // delete c1; // } // // PDG // if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) { // TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes); // c2->Write(); // delete c2; // } if (fIsChi2 && !fIsDiag) { fResidualHist = (TH1D *)fMCHist->Clone((fName + "_RESIDUAL").c_str()); fResidualHist->GetYaxis()->SetTitle("#Delta#chi^{2}"); fResidualHist->Reset(); fChi2LessBinHist = (TH1D *)fMCHist->Clone((fName + "_Chi2NMinusOne").c_str()); fChi2LessBinHist->GetYaxis()->SetTitle("Total #chi^{2} without bin_{i}"); fChi2LessBinHist->Reset(); + fIsWriting = true; (void)GetLikelihood(); + fIsWriting = false; - fResidualHist->Write(); - fChi2LessBinHist->Write(); + fResidualHist->Write((fName + "_RESIDUAL").c_str()); + fChi2LessBinHist->Write((fName + "_Chi2NMinusOne").c_str()); } // Write Extra Histograms AutoWriteExtraTH1(); WriteExtraHistograms(); // Returning NUIS_LOG(SAM, "Written Histograms: " << fName); return; } //******************************************************************** void Measurement1D::WriteRatioPlot() { //******************************************************************** // Setup mc data ratios TH1D *dataRatio = (TH1D *)fDataHist->Clone((fName + "_data_RATIO").c_str()); TH1D *mcRatio = (TH1D *)fMCHist->Clone((fName + "_MC_RATIO").c_str()); // Extra MC Data Ratios for (int i = 0; i < mcRatio->GetNbinsX(); i++) { dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1)); dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1)); mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); } // Write ratios mcRatio->Write(); dataRatio->Write(); delete mcRatio; delete dataRatio; } //******************************************************************** void Measurement1D::WriteShapePlot() { //******************************************************************** TH1D *mcShape = (TH1D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); TH1D *dataShape = (TH1D *)fDataHist->Clone((fName + "_data_SHAPE").c_str()); // Don't check error if (fShapeCovar) StatUtils::SetDataErrorFromCov(dataShape, fShapeCovar, 1E-38, false); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); std::stringstream ss; ss << shapeScale; mcShape->SetTitle(ss.str().c_str()); mcShape->SetLineWidth(3); mcShape->SetLineStyle(7); mcShape->Write(); dataShape->Write(); delete mcShape; } //******************************************************************** void Measurement1D::WriteShapeRatioPlot() { //******************************************************************** // Get a mcshape histogram TH1D *mcShape = (TH1D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); // Create shape ratio histograms TH1D *mcShapeRatio = (TH1D *)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str()); TH1D *dataShapeRatio = (TH1D *)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str()); // Divide the histograms mcShapeRatio->Divide(mcShape); dataShapeRatio->Divide(mcShape); // Colour the shape ratio plots mcShapeRatio->SetLineWidth(3); mcShapeRatio->SetLineStyle(7); mcShapeRatio->Write(); dataShapeRatio->Write(); delete mcShapeRatio; delete dataShapeRatio; } //// CRAP TO BE REMOVED //******************************************************************** void Measurement1D::SetupMeasurement(std::string inputfile, std::string type, FitWeight *rw, std::string fkdt) { //******************************************************************** nuiskey samplekey = Config::CreateKey("sample"); samplekey.Set("name", fName); samplekey.Set("type", type); samplekey.Set("input", inputfile); fSettings = LoadSampleSettings(samplekey); // Reset everything to NULL // Init(); // Check if name contains Evt, indicating that it is a raw number of events // measurements and should thus be treated as once fIsRawEvents = false; if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but fIsRawEvents == false!"); NUIS_LOG(SAM, "Overriding this and setting fIsRawEvents == true!"); } fIsEnu1D = false; if (fName.find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; NUIS_LOG(SAM, "::" << fName << "::"); NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " "not flux averaged!"); } if (fIsEnu1D && fIsRawEvents) { NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ERR(FTL, "I live in " << __FILE__ << ":" << __LINE__); throw; } fRW = rw; if (!fInput and !fIsJoint) SetupInputs(inputfile); // Set Default Options SetFitOptions(fDefaultTypes); // Set Passed Options SetFitOptions(type); // Still adding support for flat flux inputs // // Set Enu Flux Scaling // if (isFlatFluxFolding) this->Input()->ApplyFluxFolding( // this->defaultFluxHist ); // FinaliseMeasurement(); } //******************************************************************** void Measurement1D::SetupDefaultHist() { //******************************************************************** // Setup fMCHist fMCHist = (TH1D *)fDataHist->Clone(); fMCHist->SetNameTitle((fName + "_MC").c_str(), (fName + "_MC" + fPlotTitles).c_str()); // Setup fMCFine Int_t nBins = fMCHist->GetNbinsX(); fMCFine = new TH1D( (fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(), nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1)); fMCStat = (TH1D *)fMCHist->Clone(); fMCStat->Reset(); fMCHist->Reset(); fMCFine->Reset(); // Setup the NEUT Mode Array PlotUtils::CreateNeutModeArray((TH1D *)fMCHist, (TH1 **)fMCHist_PDG); PlotUtils::ResetNeutModeArray((TH1 **)fMCHist_PDG); // Setup bin masks using sample name if (fIsMask) { std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask"); if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask"; } SetBinMask(maskloc); } fMCHist_Modes = new TrueModeStack((fName + "_MODES").c_str(), ("True Channels"), fMCHist); SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); return; } //******************************************************************** void Measurement1D::SetDataValues(std::string dataFile) { //******************************************************************** // Override this function if the input file isn't in a suitable format NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str()); fDataHist = PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles); fDataTrue = (TH1D *)fDataHist->Clone(); // Number of data points is number of bins fNDataPointsX = fDataHist->GetXaxis()->GetNbins(); return; }; //******************************************************************** void Measurement1D::SetDataFromDatabase(std::string inhistfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname); fDataHist = PlotUtils::GetTH1DFromRootFile( (GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); return; }; //******************************************************************** void Measurement1D::SetDataFromFile(std::string inhistfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname); fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); return; }; //******************************************************************** void Measurement1D::SetCovarMatrix(std::string covarFile) { //******************************************************************** // Covariance function, only really used when reading in the MB Covariances. TFile *tempFile = new TFile(covarFile.c_str(), "READ"); TH2D *covarPlot = new TH2D(); TH2D *fFullCovarPlot = new TH2D(); std::string covName = ""; std::string covOption = FitPar::Config().GetParS("thrown_covariance"); if (fIsShape || fIsFree) covName = "shp_"; if (fIsDiag) covName += "diag"; else covName += "full"; covarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); if (!covOption.compare("SUB")) fFullCovarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); else if (!covOption.compare("FULL")) fFullCovarPlot = (TH2D *)tempFile->Get("fullcov"); else { NUIS_ERR(WRN, "Incorrect thrown_covariance option in parameters."); } int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral()); int covdim = int(fDataHist->GetNbinsX()); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); fDecomp = new TMatrixDSym(dim); int row, column = 0; row = 0; column = 0; for (Int_t i = 0; i < covdim; i++) { // if (this->masked->GetBinContent(i+1) > 0) continue; for (Int_t j = 0; j < covdim; j++) { // if (this->masked->GetBinContent(j+1) > 0) continue; (*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1); (*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1); column++; } column = 0; row++; } // Set bin errors on data if (!fIsDiag) { StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar); } // Get Deteriminant and inverse matrix // fCovDet = this->covar->Determinant(); TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** // Sets the covariance matrix from a provided file in a text format // scale is a multiplicative pre-factor to apply in the case where the // covariance is given in some unit (e.g. 1E-38) void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim, double scale) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream covarread(covarFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); if (covarread.is_open()) { NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile); } else { NUIS_ABORT("Covariance matrix provided is incorrect: " << covarFile); } // Loop over the lines in the file while (std::getline(covarread >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { NUIS_ERR(WRN, "SetCovarMatrixFromText -> Covariance matrix only has <= 1 " "entries on this line: " << row); } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*covar)(row, column) = *iter; (*fFullCovar)(row, column) = *iter; column++; } row++; } covarread.close(); // Scale the actualy covariance matrix by some multiplicative factor (*fFullCovar) *= scale; // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); // THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH delete this->covar; this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); // Now need to multiply by the scaling factor // If the covariance (*this->covar) *= 1. / (scale); return; }; //******************************************************************** void Measurement1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream corr(corrFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); this->fFullCovar = new TMatrixDSym(dim); if (corr.is_open()) { NUIS_LOG(SAM, "Reading and converting correlation matrix from file: " << corrFile); } else { NUIS_ABORT("Correlation matrix provided is incorrect: " << corrFile); } while (std::getline(corr >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation // matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { double val = (*iter) * this->fDataHist->GetBinError(row + 1) * 1E38 * this->fDataHist->GetBinError(column + 1) * 1E38; if (val == 0) { NUIS_ABORT("Found a zero value in the covariance matrix, assuming " "this is an error!"); } (*this->covar)(row, column) = val; (*this->fFullCovar)(row, column) = val; column++; } row++; } // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); delete this->covar; this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** // FullUnits refers to if we have "real" unscaled units in the covariance // matrix, e.g. 1E-76. If this is the case we need to scale it so that the chi2 // contribution is correct NUISANCE internally assumes the covariance matrix has // units of 1E76 void Measurement1D::SetCovarFromDataFile(std::string covarFile, std::string covName, bool FullUnits) { //******************************************************************** NUIS_LOG(SAM, "Getting covariance from " << covarFile << "->" << covName); TFile *tempFile = new TFile(covarFile.c_str(), "READ"); TH2D *covPlot = (TH2D *)tempFile->Get(covName.c_str()); covPlot->SetDirectory(0); // Scale the covariance matrix if it comes in normal units if (FullUnits) { covPlot->Scale(1.E76); } int dim = covPlot->GetNbinsX(); fFullCovar = new TMatrixDSym(dim); for (int i = 0; i < dim; i++) { for (int j = 0; j < dim; j++) { (*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1); } } this->covar = (TMatrixDSym *)fFullCovar->Clone(); fDecomp = (TMatrixDSym *)fFullCovar->Clone(); TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); TDecompChol LUChol = TDecompChol(*fDecomp); LUChol.Decompose(); fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), ""); return; }; // //******************************************************************** // void Measurement1D::SetBinMask(std::string maskFile) { // //******************************************************************** // // Create a mask histogram. // int nbins = fDataHist->GetNbinsX(); // fMaskHist = // new TH1I((fName + "_fMaskHist").c_str(), // (fName + "_fMaskHist; Bin; Mask?").c_str(), nbins, 0, nbins); // std::string line; // std::ifstream mask(maskFile.c_str(), std::ifstream::in); // if (mask.is_open()) // LOG(SAM) << "Reading bin mask from file: " << maskFile << std::endl; // else // LOG(FTL) << " Cannot find mask file." << std::endl; // while (std::getline(mask >> std::ws, line, '\n')) { // std::vector entries = GeneralUtils::ParseToInt(line, " "); // // Skip lines with poorly formatted lines // if (entries.size() < 2) { // LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << // line // << std::endl; // continue; // } // // The first index should be the bin number, the second should be the // mask // // value. // fMaskHist->SetBinContent(entries[0], entries[1]); // } // // Set masked data bins to zero // PlotUtils::MaskBins(fDataHist, fMaskHist); // return; // } // //******************************************************************** // void Measurement1D::GetBinContents(std::vector& cont, // std::vector& err) { // //******************************************************************** // // Return a vector of the main bin contents // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // cont.push_back(fMCHist->GetBinContent(i + 1)); // err.push_back(fMCHist->GetBinError(i + 1)); // } // return; // }; /* XSec Functions */ // //******************************************************************** // void Measurement1D::SetFluxHistogram(std::string fluxFile, int minE, int // maxE, // double fluxNorm) { // //******************************************************************** // // Note this expects the flux bins to be given in terms of MeV // LOG(SAM) << "Reading flux from file: " << fluxFile << std::endl; // TGraph f(fluxFile.c_str(), "%lg %lg"); // fFluxHist = // new TH1D((fName + "_flux").c_str(), (fName + "; E_{#nu} (GeV)").c_str(), // f.GetN() - 1, minE, maxE); // Double_t* yVal = f.GetY(); // for (int i = 0; i < fFluxHist->GetNbinsX(); ++i) // fFluxHist->SetBinContent(i + 1, yVal[i] * fluxNorm); // }; // //******************************************************************** // double Measurement1D::TotalIntegratedFlux(std::string intOpt, double low, // double high) { // //******************************************************************** // if (fInput->GetType() == kGiBUU) { // return 1.0; // } // // The default case of low = -9999.9 and high = -9999.9 // if (low == -9999.9) low = this->EnuMin; // if (high == -9999.9) high = this->EnuMax; // int minBin = fFluxHist->GetXaxis()->FindBin(low); // int maxBin = fFluxHist->GetXaxis()->FindBin(high); // // Get integral over custom range // double integral = fFluxHist->Integral(minBin, maxBin + 1, intOpt.c_str()); // return integral; // }; diff --git a/src/FitBase/Measurement1D.h b/src/FitBase/Measurement1D.h index 6c6b0c1..fada37f 100644 --- a/src/FitBase/Measurement1D.h +++ b/src/FitBase/Measurement1D.h @@ -1,659 +1,658 @@ // Copyright 2016 L. Pickering, P towell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MEASUREMENT_1D_H_SEEN #define MEASUREMENT_1D_H_SEEN /*! * \addtogroup FitBase * @{ */ #include #include #include #include #include #include #include #include // ROOT includes #include #include #include #include #include #include #include #include #include #include #include #include // External data fit includes #include "FitEvent.h" #include "FitUtils.h" #include "MeasurementBase.h" #include "PlotUtils.h" #include "StatUtils.h" #include "SignalDef.h" #include "MeasurementVariableBox.h" #include "MeasurementVariableBox1D.h" namespace NUISANCE { namespace FitBase { } } //******************************************************************** /// 1D Measurement base class. Histogram handling is done in this base layer. class Measurement1D : public MeasurementBase { //******************************************************************** public: /* Constructor/Deconstuctor */ Measurement1D(void); virtual ~Measurement1D(void); /* Setup Functions */ /// \brief Setup all configs once initialised /// /// Should be called after all configs have been setup inside fSettings container. /// Handles the processing of inputs and setting up of types. /// Replaces the old 'SetupMeasurement' function. void FinaliseSampleSettings(); /// \brief Creates the 1D data distribution given the binning provided. virtual void CreateDataHistogram(int dimx, double* binx); /// \brief Read 1D data inputs from a text file. /// /// Inputfile should have the format: \n /// low_binedge_1 bin_content_1 bin_error_1 \n /// low_binedge_2 bin_content_2 bin_error_2 \n /// .... .... .... \n /// high_bin_edge_N 0.0 0.0 virtual void SetDataFromTextFile(std::string datafile); /// \brief Read 1D data inputs from a root file. /// /// inhistfile specifies the path to the root file /// histname specifies the name of the histogram. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ';' so that: \n /// 'myhistfile.root;myhistname' \n /// will also work. virtual void SetDataFromRootFile(std::string inhistfile, std::string histname = ""); /// \brief Setup a default empty data histogram /// /// Only used for flattree creators. virtual void SetEmptyData(); /// \brief Set data bin errors to sqrt(entries) /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Sets the data errors as the sqrt of the bin contents /// Should be use for counting experiments virtual void SetPoissonErrors(); /// \brief Make diagonal covariance from data /// /// \warning If no histogram passed, data must be setup first! /// Setup the covariance inputs by taking the data histogram /// errors and setting up a diagonal covariance matrix. /// /// If no data is supplied, fDataHist is used if already set. virtual void SetCovarFromDiagonal(TH1D* data = NULL); /// \brief Read the data covariance from a text file. /// /// Inputfile should have the format: \n /// covariance_11 covariance_12 covariance_13 ... \n /// covariance_21 covariance_22 covariance_23 ... \n /// ... ... ... ... \n /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarFromTextFile(std::string covfile, int dim = -1); virtual void SetCovarFromMultipleTextFiles(std::string covfiles, int dim = -1); /// \brief Read the data covariance from a ROOT file. /// /// - covfile specifies the full path to the file /// - histname specifies the name of the covariance object. Both TMatrixDSym and TH2D are supported. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the inverted data covariance from a text file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarInvertFromTextFile(std::string covfile, int dim = -1); /// \brief Read the inverted data covariance from a ROOT file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarInvertFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the data correlations from a text file. /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCorrelationFromTextFile(std::string covfile, int dim = -1); /// \brief Read the data correlations from multiple text files. /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of the first corrfile. virtual void SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim = -1); /// \brief Read the data correlations from a ROOT file. /// /// \warning REQUIRES DATA TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCorrelationFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the cholesky decomposed covariance from a text file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCholDecompFromTextFile(std::string covfile, int dim = -1); /// \brief Read the cholesky decomposed covariance from a ROOT file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCholDecompFromRootFile(std::string covfile, std::string histname=""); /// \brief Try to extract a shape-only matrix from the existing covariance virtual void SetShapeCovar(); /// \brief Scale the data by some scale factor virtual void ScaleData(double scale); /// \brief Scale the data error bars by some scale factor virtual void ScaleDataErrors(double scale); /// \brief Scale the covariaince and its invert/decomp by some scale factor. virtual void ScaleCovar(double scale); /// \brief Setup a bin masking histogram and apply masking to data /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Reads in a list of bins in a text file to be masked. Format is: \n /// bin_index_1 1 \n /// bin_index_2 1 \n /// bin_index_3 1 \n /// /// If 0 is given then a bin entry will NOT be masked. So for example: \n\n /// 1 1 \n /// 2 1 \n /// 3 0 \n /// 4 1 \n\n /// Will mask only the 1st, 2nd, and 4th bins. /// /// Masking can be turned on by specifiying the MASK option when creating a sample. /// When this is passed NUISANCE will look in the following locations for the mask file: /// - FitPar::Config().GetParS(fName + ".mask") /// - "data/masks/" + fName + ".mask"; virtual void SetBinMask(std::string maskfile); /// \brief Set the current fit options from a string. /// /// This is called twice for each sample, once to set the default /// and once to set the current setting (if anything other than default given) /// /// For this to work properly it requires the default and allowed types to be /// set correctly. These should be specified as a string listing options. /// /// To split up options so that NUISANCE can automatically detect ones that /// are conflicting. Any options seperated with the '/' symbol are non conflicting /// and can be given together, whereas any seperated with the ',' symbol cannot /// be specified by the end user at the same time. /// /// Default Type Examples: /// - DIAG/FIX = Default option will be a diagonal covariance, with FIXED norm. /// - MASK/SHAPE = Default option will be a masked hist, with SHAPE always on. /// /// Allowed Type examples: /// - 'FULL/DIAG/NORM/MASK' = Any of these options can be specified. /// - 'FULL,FREE,SHAPE/MASK/NORM' = User can give either FULL, FREE, or SHAPE as on option. /// MASK and NORM can also be included as options. virtual void SetFitOptions(std::string opt); /// \brief Final constructor setup /// \warning Should be called right at the end of the constructor. /// /// Contains a series of checks to ensure the data and inputs have been setup. /// Also creates the MC histograms needed for fitting. void FinaliseMeasurement(); /* Smearing */ /// \brief Read in smearing matrix from file /// /// Set the smearing matrix from a text file given the size of the matrix virtual void SetSmearingMatrix(std::string smearfile, int truedim, int recodim); /// \brief Apply smearing to MC true to get MC reco /// /// Apply smearing matrix to fMCHist using fSmearingMatrix virtual void ApplySmearingMatrix(void); /* Reconfigure Functions */ /// \brief Create a Measurement1D box /// /// Creates a new 1D variable box containing just fXVar. /// /// This box is the bare minimum required by the JointFCN when /// running fast reconfigures during a routine. /// If for some reason a sample needs extra variables to be saved then /// it should override this function creating its own MeasurementVariableBox /// that contains the extra variables. virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox1D();}; /// \brief Reset all MC histograms /// /// Resets all standard histograms and those registered to auto /// process to zero. /// /// If extra histograms are not included in auto processing, then they must be reset /// by overriding this function and doing it manually if required. virtual void ResetAll(void); /// \brief Fill MC Histograms from XVar /// /// Fill standard histograms using fXVar, Weight read from the variable box. /// /// WARNING : Any extra MC histograms need to be filled by overriding this function, /// even if they have been set to auto process. virtual void FillHistograms(void); // \brief Convert event rates to final histogram /// /// Apply standard scaling procedure to standard mc histograms to convert from /// raw events to xsec prediction. /// /// If any distributions have been set to auto process /// that is done during this function call, and a differential xsec is assumed. /// If that is not the case this function must be overriden. virtual void ScaleEvents(void); /// \brief Scale MC by a factor=1/norm /// /// Apply a simple normalisation scaling if the option FREE or a norm_parameter /// has been specified in the NUISANCE routine. virtual void ApplyNormScale(double norm); /* Statistical Functions */ /// \brief Get Number of degrees of freedom /// /// Returns the number bins inside the data histogram accounting for /// any bin masking applied. virtual int GetNDOF(void); /// \brief Return Data/MC Likelihood at current state /// /// Returns the likelihood of the data given the current MC prediction. /// Diferent likelihoods definitions are used depending on the FitOptions. virtual double GetLikelihood(void); /* Fake Data */ /// \brief Set the fake data values from either a file, or MC /// /// - Setting from a file "path": \n /// When reading from a file the full path must be given to a standard /// nuisance output. The standard MC histogram should have a name that matches /// this sample for it to be read in. /// \n\n /// - Setting from "MC": \n /// If the MC option is given the current MC prediction is used as fake data. virtual void SetFakeDataValues(std::string fakeOption); /// \brief Reset fake data back to starting fake data /// /// Reset the fake data back to original fake data (Reset back to before /// ThrowCovariance was first called) virtual void ResetFakeData(void); /// \brief Reset fake data back to original data /// /// Reset the data histogram back to the true original dataset for this sample /// before any fake data was defined. virtual void ResetData(void); /// \brief Generate fake data by throwing the covariance. /// /// Can be used on fake MC data or just the original dataset. /// Call ResetFakeData or ResetData to return to values before the throw. virtual void ThrowCovariance(void); /// \brief Throw the data by its assigned errors and assign this to MC /// /// Used when creating data toys by assign the MC to this thrown data /// so that the likelihood is calculated between data and thrown data virtual void ThrowDataToy(void); /* Access Functions */ /// \brief Returns nicely formatted MC Histogram /// /// Format options can also be given in the samplesettings: /// - linecolor /// - linestyle /// - linewidth /// - fillcolor /// - fillstyle /// /// So to have a sample line colored differently in the xml cardfile put: \n /// virtual TH1D* GetMCHistogram(void); /// \brief Returns nicely formatted data Histogram /// /// Format options can also be given in the samplesettings: /// - datacolor /// - datastyle /// - datawidth /// /// So to have a sample data colored differently in the xml cardfile put: \n /// virtual TH1D* GetDataHistogram(void); /// \brief Returns a list of all MC histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMCList(void) { return std::vector(1, GetMCHistogram()); } /// \brief Returns a list of all Data histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetDataList(void) { return std::vector(1, GetDataHistogram()); } /// \brief Returns a list of all Mask histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMaskList(void) { return std::vector(1, fMaskHist); }; /// \brief Returns a list of all Fine histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetFineList(void) { return std::vector(1, fMCFine); }; /* Write Functions */ /// \brief Save the current state to the current TFile directory \n /// /// Data/MC are both saved by default. /// A range of other histograms can be saved by setting the /// config option 'drawopts'. /// /// Possible options: \n /// - FINE = Write Fine Histogram \n /// - WEIGHTS = Write Weighted MC Histogram (before scaling) \n /// - FLUX = Write Flux Histogram from MC Input \n /// - EVT = Write Event Histogram from MC Input \n /// - XSEC = Write XSec Histogram from MC Input \n /// - MASK = Write Mask Histogram \n /// - COV = Write Covariance Histogram \n /// - INVCOV = Write Inverted Covariance Histogram \n /// - DECMOP = Write Decomp. Covariance Histogram \n /// - RESIDUAL= Write Resudial Histograms \n /// - RATIO = Write Data/MC Ratio Histograms \n /// - SHAPE = Write MC Shape Histograms norm. to Data \n /// - CANVMC = Build MC Canvas Showing Data, MC, Shape \n /// - MODES = Write PDG Stack \n /// - CANVPDG = Build MC Canvas Showing Data, PDGStack \n /// /// So to save a range of these in parameters/config.xml set: \n /// virtual void Write(std::string drawOpt); virtual void WriteRatioPlot(); virtual void WriteShapePlot(); virtual void WriteShapeRatioPlot(); //* // OLD DEFUNCTIONS // /// OLD FUNCTION virtual void SetupMeasurement(std::string input, std::string type, FitWeight* rw, std::string fkdt); /// OLD FUNCTION virtual void SetupDefaultHist(void); /// OLD FUNCTION virtual void SetDataValues(std::string dataFile); /// OLD FUNCTION virtual void SetDataFromFile(std::string inhistfile, std::string histname); /// OLD FUNCTION virtual void SetDataFromDatabase(std::string inhistfile, std::string histname); /// OLD FUNCTION virtual void SetCovarMatrix(std::string covarFile); /// OLD FUNCTION virtual void SetCovarMatrixFromText(std::string covarFile, int dim, double scale = 1.0); /// OLD FUNCTION virtual void SetCovarMatrixFromCorrText(std::string covarFile, int dim); /// OLD FUNCTION virtual void SetCovarFromDataFile(std::string covarFile, std::string covName, bool FullUnits = false); /// OLD FUNCTION // virtual THStack GetModeStack(void); protected: // Data TH1D* fDataHist; ///< default data histogram TH1D* fDataOrig; ///< histogram to store original data before throws. TH1D* fDataTrue; ///< histogram to store true dataset std::string fPlotTitles; ///< Plot title x and y for the histograms // MC TH1D* fMCHist; ///< default MC Histogram used in the chi2 fits TH1D* fMCFine; ///< finely binned MC histogram TH1D* fMCStat; ///< histogram with unweighted events to properly calculate TH1D* fMCWeighted; ///< Weighted histogram before xsec scaling TH1I* fMaskHist; ///< Mask histogram for neglecting specific bins TMatrixD* fSmearMatrix; ///< Smearing matrix (note, this is not symmetric) TH1D *fResidualHist; TH1D *fChi2LessBinHist; TrueModeStack* fMCHist_Modes; ///< Optional True Mode Stack // Statistical TMatrixDSym* covar; ///< Inverted Covariance TMatrixDSym* fFullCovar; ///< Full Covariance TMatrixDSym* fDecomp; ///< Decomposed Covariance TMatrixDSym* fCorrel; ///< Correlation Matrix TMatrixDSym* fShapeCovar; ///< Shape-only covariance TMatrixDSym* fShapeDecomp; ///< Decomposed shape-only covariance TMatrixDSym* fShapeInvert; ///< Inverted shape-only covariance TMatrixDSym* fCovar; ///< New FullCovar TMatrixDSym* fInvert; ///< New covar double fNormError; ///< Sample norm error double fLikelihood; ///< Likelihood value // Fake Data bool fIsFakeData; ///< Flag: is the current data fake from MC std::string fFakeDataInput; ///< Input fake data file path TFile* fFakeDataFile; ///< Input fake data file // Fit specific flags std::string fFitType; ///< Current fit type std::string fAllowedTypes; ///< Fit Types Possible std::string fDefaultTypes; ///< Starting Default Fit Types bool fIsShape; ///< Flag : Perform Shape-only fit bool fIsFree; ///< Flag : Perform normalisation free fit bool fIsDiag; ///< Flag : only include uncorrelated diagonal errors bool fIsMask; ///< Flag : Apply bin masking bool fIsRawEvents; ///< Flag : Are bin contents just event rates bool fIsEnu1D; ///< Flag : Perform Flux Unfolded Scaling bool fIsChi2SVD; ///< Flag : Use alternative Chi2 SVD Method (Do not use) bool fAddNormPen; ///< Flag : Add a normalisation penalty term to the chi2. bool fIsFix; ///< Flag : keeping norm fixed bool fIsFull; ///< Flag : using full covariaince bool fIsDifXSec; ///< Flag : creating a dif xsec bool fIsChi2; ///< Flag : using Chi2 over LL methods bool fIsSmeared; ///< Flag : Apply smearing? - - + bool fIsWriting; /// OLD STUFF TO REMOVE TH1D* fMCHist_PDG[61]; ///< REMOVE OLD MC PDG Plot // Arrays for data entries Double_t* fXBins; ///< REMOVE xBin edges Double_t* fDataValues; ///< REMOVE data bin contents Double_t* fDataErrors; ///< REMOVE data bin errors Int_t fNDataPointsX; ///< REMOVE number of data points }; /*! @} */ #endif diff --git a/src/FitBase/Measurement2D.cxx b/src/FitBase/Measurement2D.cxx index d91a124..1a3186b 100644 --- a/src/FitBase/Measurement2D.cxx +++ b/src/FitBase/Measurement2D.cxx @@ -1,2032 +1,2057 @@ // 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 "Measurement2D.h" #include "TDecompChol.h" //******************************************************************** Measurement2D::Measurement2D(void) { //******************************************************************** covar = NULL; fDecomp = NULL; fFullCovar = NULL; fMCHist = NULL; fMCFine = NULL; fDataHist = NULL; fMCHist_X = NULL; fMCHist_Y = NULL; fDataHist_X = NULL; fDataHist_Y = NULL; fMaskHist = NULL; fMapHist = NULL; fDataOrig = NULL; fDataTrue = NULL; fMCWeighted = NULL; fResidualHist = NULL; fChi2LessBinHist = NULL; fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/FITPROJX/" "FITPROJY"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsDifXSec = false; fIsEnu = false; // XSec Scalings fScaleFactor = -1.0; fCurrentNorm = 1.0; // Histograms fDataHist = NULL; fDataTrue = NULL; fMCHist = NULL; fMCFine = NULL; fMCWeighted = NULL; fMaskHist = NULL; // Covar covar = NULL; fFullCovar = NULL; fCovar = NULL; fInvert = NULL; fDecomp = NULL; // Fake Data fFakeDataInput = ""; fFakeDataFile = NULL; // Options fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsDifXSec = false; fIsEnu1D = false; + fIsWriting = false; // Inputs fInput = NULL; fRW = NULL; // Extra Histograms fMCHist_Modes = NULL; } //******************************************************************** Measurement2D::~Measurement2D(void) { //******************************************************************** if (fDataHist) delete fDataHist; if (fDataTrue) delete fDataTrue; if (fMCHist) delete fMCHist; if (fMCFine) delete fMCFine; if (fMCWeighted) delete fMCWeighted; if (fMaskHist) delete fMaskHist; if (covar) delete covar; if (fFullCovar) delete fFullCovar; if (fCovar) delete fCovar; if (fInvert) delete fInvert; if (fDecomp) delete fDecomp; delete fResidualHist; delete fChi2LessBinHist; } //******************************************************************** void Measurement2D::FinaliseSampleSettings() { //******************************************************************** MeasurementBase::FinaliseSampleSettings(); // Setup naming + renaming fName = fSettings.GetName(); fSettings.SetS("originalname", fName); if (fSettings.Has("rename")) { fName = fSettings.GetS("rename"); fSettings.SetS("name", fName); } // Setup all other options NUIS_LOG(SAM, "Finalising Sample Settings: " << fName); if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods."); } if (fSettings.GetS("originalname").find("Enu") != std::string::npos) { fIsEnu1D = true; NUIS_LOG(SAM, "::" << fName << "::"); NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " << "not flux averaged!"); } if (fIsEnu1D && fIsRawEvents) { NUIS_ERR(FTL, "Found 2D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__); } if (!fRW) fRW = FitBase::GetRW(); if (!fInput) SetupInputs(fSettings.GetS("input")); // Setup options SetFitOptions(fDefaultTypes); // defaults SetFitOptions(fSettings.GetS("type")); // user specified EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min")); EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max")); if (fAddNormPen) { fNormError = fSettings.GetNormError(); if (fNormError <= 0.0) { NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!"); NUIS_ABORT("If you want to use it please add fNormError=VAL"); } } } void Measurement2D::CreateDataHistogram(int dimx, double *binx, int dimy, double *biny) { if (fDataHist) delete fDataHist; NUIS_LOG(SAM, "Creating Data Histogram dim : " << dimx << " " << dimy); fDataHist = new TH2D((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(), dimx - 1, binx, dimy - 1, biny); } void Measurement2D::SetDataFromTextFile(std::string data, std::string binx, std::string biny) { // Get the data hist fDataHist = PlotUtils::GetTH2DFromTextFile(data, binx, biny); // Set the name properly fDataHist->SetName((fSettings.GetName() + "_data").c_str()); fDataHist->SetTitle(fSettings.GetFullTitles().c_str()); } void Measurement2D::SetDataFromRootFile(std::string datfile, std::string histname) { NUIS_LOG(SAM, "Reading data from root file: " << datfile << ";" << histname); fDataHist = PlotUtils::GetTH2DFromRootFile(datfile, histname); fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str()); } void Measurement2D::SetDataValuesFromTextFile(std::string datfile, TH2D *hist) { NUIS_LOG(SAM, "Setting data values from text file"); if (!hist) hist = fDataHist; // Read TH2D From textfile TH2D *valhist = (TH2D *)hist->Clone(); valhist->Reset(); PlotUtils::Set2DHistFromText(datfile, valhist, 1.0, true); NUIS_LOG(SAM, " -> Filling values from read hist."); for (int i = 0; i < valhist->GetNbinsX(); i++) { for (int j = 0; j < valhist->GetNbinsY(); j++) { hist->SetBinContent(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1)); } } NUIS_LOG(SAM, " --> Done"); } void Measurement2D::SetDataErrorsFromTextFile(std::string datfile, TH2D *hist) { NUIS_LOG(SAM, "Setting data errors from text file"); if (!hist) hist = fDataHist; // Read TH2D From textfile TH2D *valhist = (TH2D *)hist->Clone(); valhist->Reset(); PlotUtils::Set2DHistFromText(datfile, valhist, 1.0); // Fill Errors NUIS_LOG(SAM, " -> Filling errors from read hist."); for (int i = 0; i < valhist->GetNbinsX(); i++) { for (int j = 0; j < valhist->GetNbinsY(); j++) { hist->SetBinError(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1)); } } NUIS_LOG(SAM, " --> Done"); } void Measurement2D::SetMapValuesFromText(std::string dataFile) { TH2D *hist = fDataHist; std::vector edgex; std::vector edgey; for (int i = 0; i <= hist->GetNbinsX(); i++) edgex.push_back(hist->GetXaxis()->GetBinLowEdge(i + 1)); for (int i = 0; i <= hist->GetNbinsY(); i++) edgey.push_back(hist->GetYaxis()->GetBinLowEdge(i + 1)); fMapHist = new TH2I((fName + "_map").c_str(), (fName + fPlotTitles).c_str(), edgex.size() - 1, &edgex[0], edgey.size() - 1, &edgey[0]); NUIS_LOG(SAM, "Reading map from: " << dataFile); PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0); } //******************************************************************** void Measurement2D::SetPoissonErrors() { //******************************************************************** if (!fDataHist) { NUIS_ERR(FTL, "Need a data hist to setup possion errors! "); NUIS_ABORT("Setup Data First!"); } for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) { fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1))); } } //******************************************************************** void Measurement2D::SetCovarFromDiagonal(TH2D *data) { //******************************************************************** if (!data and fDataHist) { data = fDataHist; } if (data) { NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName()); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } else { NUIS_ABORT("No data input provided to set diagonal covar from!"); } // if (!fIsDiag) { // ERR(FTL) << "SetCovarMatrixFromDiag called for measurement " // << "that is not set as diagonal." << std::endl; // throw; // } } //******************************************************************** void Measurement2D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = this->GetNDOF(); } NUIS_LOG(SAM, "Reading covariance from text file: " << covfile << " " << dim); fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCovarFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading covariance from text file: " << covfile << ";" << histname); fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCovarInvertFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = this->GetNDOF(); } NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile); covar = StatUtils::GetCovarFromTextFile(covfile, dim); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";" << histname); covar = StatUtils::GetCovarFromRootFile(covfile, histname); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement2D::SetCorrelationFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) dim = this->GetNDOF(); NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << dim); TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim); if (!fDataHist) { NUIS_ABORT("Trying to set correlations from text file but there is no " "data to build it from. \n" << "In constructor make sure data is set before " "SetCorrelationFromTextFile is called. \n"); } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(dim); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement2D::SetCorrelationFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << histname); TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname); if (!fDataHist) { NUIS_ABORT("Trying to set correlations from text file but there is no " "data to build it from. \n" << "In constructor make sure data is set before " "SetCorrelationFromTextFile is called. \n"); } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX()); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement2D::SetCholDecompFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = this->GetNDOF(); } NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile << " " << dim); TMatrixD *temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement2D::SetCholDecompFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";" << histname); TMatrixD *temp = StatUtils::GetMatrixFromRootFile(covfile, histname); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement2D::ScaleData(double scale) { //******************************************************************** fDataHist->Scale(scale); } //******************************************************************** void Measurement2D::ScaleDataErrors(double scale) { //******************************************************************** for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsY(); j++) { fDataHist->SetBinError(i + 1, j + 1, fDataHist->GetBinError(i + 1, j + 1) * scale); } } } //******************************************************************** void Measurement2D::ScaleCovar(double scale) { //******************************************************************** (*fFullCovar) *= scale; (*covar) *= 1.0 / scale; (*fDecomp) *= sqrt(scale); } //******************************************************************** void Measurement2D::SetBinMask(std::string maskfile) { //******************************************************************** if (!fIsMask) return; NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile); // Create a mask histogram with dim of data int nbinsx = fDataHist->GetNbinsX(); int nbinxy = fDataHist->GetNbinsY(); fMaskHist = new TH2I((fSettings.GetName() + "_BINMASK").c_str(), (fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbinsx, 0, nbinsx, nbinxy, 0, nbinxy); std::string line; std::ifstream mask(maskfile.c_str(), std::ifstream::in); if (!mask.is_open()) { NUIS_LOG(FTL, " Cannot find mask file."); throw; } while (std::getline(mask >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToInt(line, " "); // Skip lines with poorly formatted lines if (entries.size() < 2) { NUIS_LOG(WRN, "Measurement2D::SetBinMask(), couldn't parse line: " << line); continue; } // The first index should be the bin number, the second should be the mask // value. int val = 0; if (entries[2] > 0) val = 1; fMaskHist->SetBinContent(entries[0], entries[1], val); } // Apply masking by setting masked data bins to zero PlotUtils::MaskBins(fDataHist, fMaskHist); return; } //******************************************************************** void Measurement2D::FinaliseMeasurement() { //******************************************************************** NUIS_LOG(SAM, "Finalising Measurement: " << fName); if (fSettings.GetB("onlymc")) { if (fDataHist) delete fDataHist; fDataHist = new TH2D("empty_data", "empty_data", 1, 0.0, 1.0, 1, 0.0, 1.0); } // Make sure data is setup if (!fDataHist) { NUIS_ABORT("No data has been setup inside " << fName << " constructor!"); } // Make sure covariances are setup if (!fFullCovar) { fIsDiag = true; SetCovarFromDiagonal(fDataHist); + } else if (fIsDiag) { // Have covariance but also set Diag + NUIS_LOG(SAM, "Have full covariance for sample " + << GetName() + << " but only using diagonal elements for likelihood"); + size_t nbins = fFullCovar->GetNcols(); + for (int i = 0; i < nbins; ++i) { + for (int j = 0; j < nbins; ++j) { + if (i != j) { + (*fFullCovar)[i][j] = 0; + } + } + } + delete covar; + covar = NULL; + delete fDecomp; + fDecomp = NULL; } if (!covar) { covar = StatUtils::GetInvert(fFullCovar); } if (!fDecomp) { fDecomp = StatUtils::GetDecomp(fFullCovar); } // Setup fMCHist from data fMCHist = (TH2D *)fDataHist->Clone(); fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(), (fSettings.GetFullTitles()).c_str()); fMCHist->Reset(); // Setup fMCFine fMCFine = new TH2D( "mcfine", "mcfine", fDataHist->GetNbinsX() * 6, fMCHist->GetXaxis()->GetBinLowEdge(1), fMCHist->GetXaxis()->GetBinLowEdge(fDataHist->GetNbinsX() + 1), fDataHist->GetNbinsY() * 6, fMCHist->GetYaxis()->GetBinLowEdge(1), fMCHist->GetYaxis()->GetBinLowEdge(fDataHist->GetNbinsY() + 1)); fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(), (fSettings.GetFullTitles()).c_str()); fMCFine->Reset(); // Setup MC Stat fMCStat = (TH2D *)fMCHist->Clone(); fMCStat->Reset(); // Search drawopts for possible types to include by default std::string drawopts = FitPar::Config().GetParS("drawopts"); if (drawopts.find("MODES") != std::string::npos) { fMCHist_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(), ("True Channels"), fMCHist); SetAutoProcessTH1(fMCHist_Modes); } if (fSettings.Has("maskfile") && fSettings.Has("maskhist")) { fMaskHist = PlotUtils::GetTH2FromRootFile(fSettings.GetS("maskfile"), fSettings.GetS("maskhist")); - fIsMask = true; - } - - // Setup bin masks using sample name - if (fIsMask && !fMaskHist) { + fIsMask = bool(fMaskHist); + NUIS_LOG(SAM, "Loaded mask histogram: " << fSettings.GetS("maskhist") + << " from " + << fSettings.GetS("maskfile")); + } else if (fIsMask) { // Setup bin masks using sample name std::string curname = fName; std::string origname = fSettings.GetS("originalname"); // Check rename.mask std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask"); // Check origname.mask if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask"); // Check database if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask"; } // Setup Bin Mask SetBinMask(maskloc); } if (fScaleFactor < 0) { NUIS_ERR(FTL, "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__); NUIS_ERR(FTL, "fScaleFactor = " << fScaleFactor); NUIS_ABORT("EXITING"); } // Create and fill Weighted Histogram if (!fMCWeighted) { fMCWeighted = (TH2D *)fMCHist->Clone(); fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(), (fName + "_MCWGHTS" + fPlotTitles).c_str()); fMCWeighted->GetYaxis()->SetTitle("Weighted Events"); } if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist); } //******************************************************************** void Measurement2D::SetFitOptions(std::string opt) { //******************************************************************** // Do nothing if default given if (opt == "DEFAULT") return; // CHECK Conflicting Fit Options std::vector fit_option_allow = GeneralUtils::ParseToStr(fAllowedTypes, "/"); for (UInt_t i = 0; i < fit_option_allow.size(); i++) { std::vector fit_option_section = GeneralUtils::ParseToStr(fit_option_allow.at(i), ","); bool found_option = false; for (UInt_t j = 0; j < fit_option_section.size(); j++) { std::string av_opt = fit_option_section.at(j); if (!found_option and opt.find(av_opt) != std::string::npos) { found_option = true; } else if (found_option and opt.find(av_opt) != std::string::npos) { NUIS_ABORT( "ERROR: Conflicting fit options provided: " << opt << std::endl << "Conflicting group = " << fit_option_section.at(i) << std::endl << "You should only supply one of these options in card file."); } } } // Check all options are allowed std::vector fit_options_input = GeneralUtils::ParseToStr(opt, "/"); for (UInt_t i = 0; i < fit_options_input.size(); i++) { if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) { NUIS_ERR(FTL, "ERROR: Fit Option '" << fit_options_input.at(i) << "' Provided is not allowed for this measurement."); NUIS_ERR(FTL, "Fit Options should be provided as a '/' seperated list " "(e.g. FREE/DIAG/NORM)"); NUIS_ABORT("Available options for " << fName << " are '" << fAllowedTypes << "'"); } } // Set TYPE fFitType = opt; // FIX,SHAPE,FREE if (opt.find("FIX") != std::string::npos) { fIsFree = fIsShape = false; fIsFix = true; } else if (opt.find("SHAPE") != std::string::npos) { fIsFree = fIsFix = false; fIsShape = true; } else if (opt.find("FREE") != std::string::npos) { fIsFix = fIsShape = false; fIsFree = true; } // DIAG,FULL (or default to full) if (opt.find("DIAG") != std::string::npos) { fIsDiag = true; fIsFull = false; } else if (opt.find("FULL") != std::string::npos) { fIsDiag = false; fIsFull = true; } // CHI2/LL (OTHERS?) if (opt.find("LOG") != std::string::npos) { fIsChi2 = false; NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!"); NUIS_ABORT("Try to use a chi2!"); } else { fIsChi2 = true; } // EXTRAS if (opt.find("RAW") != std::string::npos) fIsRawEvents = true; if (opt.find("DIF") != std::string::npos) fIsDifXSec = true; if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true; if (opt.find("NORM") != std::string::npos) fAddNormPen = true; if (opt.find("MASK") != std::string::npos) fIsMask = true; // Set TYPE fFitType = opt; // FIX,SHAPE,FREE if (opt.find("FIX") != std::string::npos) { fIsFree = fIsShape = false; fIsFix = true; } else if (opt.find("SHAPE") != std::string::npos) { fIsFree = fIsFix = false; fIsShape = true; } else if (opt.find("FREE") != std::string::npos) { fIsFix = fIsShape = false; fIsFree = true; } // DIAG,FULL (or default to full) if (opt.find("DIAG") != std::string::npos) { fIsDiag = true; fIsFull = false; } else if (opt.find("FULL") != std::string::npos) { fIsDiag = false; fIsFull = true; } // CHI2/LL (OTHERS?) if (opt.find("LOG") != std::string::npos) fIsChi2 = false; else fIsChi2 = true; // EXTRAS if (opt.find("RAW") != std::string::npos) fIsRawEvents = true; if (opt.find("DIF") != std::string::npos) fIsDifXSec = true; if (opt.find("ENU1D") != std::string::npos) fIsEnu = true; if (opt.find("NORM") != std::string::npos) fAddNormPen = true; if (opt.find("MASK") != std::string::npos) fIsMask = true; fIsProjFitX = (opt.find("FITPROJX") != std::string::npos); fIsProjFitY = (opt.find("FITPROJY") != std::string::npos); return; }; /* Reconfigure LOOP */ //******************************************************************** void Measurement2D::ResetAll() { //******************************************************************** fMCHist->Reset(); fMCFine->Reset(); fMCStat->Reset(); return; }; //******************************************************************** void Measurement2D::FillHistograms() { //******************************************************************** if (Signal) { fMCHist->Fill(fXVar, fYVar, Weight); fMCFine->Fill(fXVar, fYVar, Weight); fMCStat->Fill(fXVar, fYVar, 1.0); if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, fYVar, Weight); } return; }; //******************************************************************** void Measurement2D::ScaleEvents() { //******************************************************************** // Fill MCWeighted; // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1)); // fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1)); // } // Setup Stat ratios for MC and MC Fine double *statratio = new double[fMCHist->GetNbinsX()]; for (int i = 0; i < fMCHist->GetNbinsX(); i++) { if (fMCHist->GetBinContent(i + 1) != 0) { statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1); } else { statratio[i] = 0.0; } } double *statratiofine = new double[fMCFine->GetNbinsX()]; for (int i = 0; i < fMCFine->GetNbinsX(); i++) { if (fMCFine->GetBinContent(i + 1) != 0) { statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1); } else { statratiofine[i] = 0.0; } } // Scaling for raw event rates if (fIsRawEvents) { double datamcratio = fDataHist->Integral() / fMCHist->Integral(); fMCHist->Scale(datamcratio); fMCFine->Scale(datamcratio); if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio); // Scaling for XSec as function of Enu } else if (fIsEnu1D) { PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(), GetEventHistogram(), fScaleFactor); PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(), GetEventHistogram(), fScaleFactor); // if (fMCHist_Modes) { // PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(), // GetEventHistogram(), fScaleFactor, // fNEvents); // } // Any other differential scaling } else { fMCHist->Scale(fScaleFactor, "width"); fMCFine->Scale(fScaleFactor, "width"); // if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width"); } // Proper error scaling - ROOT Freaks out with xsec weights sometimes for (int i = 0; i < fMCStat->GetNbinsX(); i++) { fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]); } for (int i = 0; i < fMCFine->GetNbinsX(); i++) { fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]); } // Clean up delete statratio; delete statratiofine; return; }; //******************************************************************** void Measurement2D::ApplyNormScale(double norm) { //******************************************************************** fCurrentNorm = norm; fMCHist->Scale(1.0 / norm); fMCFine->Scale(1.0 / norm); return; }; /* Statistic Functions - Outsources to StatUtils */ //******************************************************************** int Measurement2D::GetNDOF() { //******************************************************************** // Just incase it has gone... if (!fDataHist) return -1; int nDOF = 0; // If datahist has no errors make sure we don't include those bins as they are // not data points for (int xBin = 0; xBin < fDataHist->GetNbinsX(); ++xBin) { for (int yBin = 0; yBin < fDataHist->GetNbinsY(); ++yBin) { if (fDataHist->GetBinError(xBin + 1, yBin + 1) != 0) ++nDOF; } } // Account for possible bin masking int nMasked = 0; if (fMaskHist and fIsMask) if (fMaskHist->Integral() > 0) for (int xBin = 0; xBin < fMaskHist->GetNbinsX() + 1; ++xBin) for (int yBin = 0; yBin < fMaskHist->GetNbinsY() + 1; ++yBin) if (fMaskHist->GetBinContent(xBin, yBin) > 0.5) ++nMasked; // Take away those masked DOF if (fIsMask) { nDOF -= nMasked; } return nDOF; } //******************************************************************** double Measurement2D::GetLikelihood() { //******************************************************************** // If this is for a ratio, there is no data histogram to compare to! if (fNoData || !fDataHist) return 0.; // Fix weird masking bug if (!fIsMask) { if (fMaskHist) { fMaskHist = NULL; } } else { if (fMaskHist) { PlotUtils::MaskBins(fMCHist, fMaskHist); } } // if (fIsProjFitX or fIsProjFitY) return GetProjectedChi2(); // Scale up the results to match each other (Not using width might be // inconsistent with Meas1D) double scaleF = fDataHist->Integral() / fMCHist->Integral(); if (fIsShape) { fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); // PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF); } - if (!fMapHist) + if (!fMapHist) { fMapHist = StatUtils::GenerateMap(fDataHist); + } // Get the chi2 from either covar or diagonals double chi2 = 0.0; if (fIsChi2) { if (fIsDiag) { chi2 = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMapHist, fMaskHist); } else { chi2 = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist, - fMaskHist, fResidualHist); - if (fChi2LessBinHist) { - TH2I *binmask = new TH2I("mask", "", fDataHist->GetNbinsX(), 0, - fDataHist->GetNbinsX(), fDataHist->GetNbinsY(), - 0, fDataHist->GetNbinsY()); - binmask->SetDirectory(NULL); + fMaskHist, + fIsWriting ? fResidualHist : NULL); + if (fChi2LessBinHist && fIsWriting) { + NUIS_LOG(SAM, "Building n-1 chi2 contribution plot for " << GetName()); for (int xi = 0; xi < fDataHist->GetNbinsX(); ++xi) { for (int yi = 0; yi < fDataHist->GetNbinsY(); ++yi) { - binmask->Reset(); + TH2I *binmask = + fMaskHist + ? static_cast(fMaskHist->Clone("mask")) + : new TH2I("mask", "", fDataHist->GetNbinsX(), 0, + fDataHist->GetNbinsX(), fDataHist->GetNbinsY(), + 0, fDataHist->GetNbinsY()); + binmask->SetDirectory(NULL); + binmask->SetBinContent(xi + 1, yi + 1, 1); fChi2LessBinHist->SetBinContent( xi + 1, yi + 1, StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist, binmask)); + delete binmask; } } - delete binmask; } } } // Add a normal penalty term if (fAddNormPen) { chi2 += (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError); NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError)); } // Adjust the shape back to where it was. if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) { fMCHist->Scale(1. / scaleF); fMCFine->Scale(1. / scaleF); } fLikelihood = chi2; return chi2; } /* Fake Data Functions */ //******************************************************************** void Measurement2D::SetFakeDataValues(std::string fakeOption) { //******************************************************************** // Setup original/datatrue TH2D *tempdata = (TH2D *)fDataHist->Clone(); if (!fIsFakeData) { fIsFakeData = true; // Make a copy of the original data histogram. if (!fDataOrig) fDataOrig = (TH2D *)fDataHist->Clone((fName + "_data_original").c_str()); } else { ResetFakeData(); } // Setup Inputs fFakeDataInput = fakeOption; NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput); // From MC if (fFakeDataInput.compare("MC") == 0) { fDataHist = (TH2D *)fMCHist->Clone((fName + "_MC").c_str()); // Fake File } else { if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ"); fDataHist = (TH2D *)fFakeDataFile->Get((fName + "_MC").c_str()); } // Setup Data Hist fDataHist->SetNameTitle((fName + "_FAKE").c_str(), (fName + fPlotTitles).c_str()); // Replace Data True if (fDataTrue) delete fDataTrue; fDataTrue = (TH2D *)fDataHist->Clone(); fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(), (fName + fPlotTitles).c_str()); // Make a new covariance for fake data hist. int nbins = fDataHist->GetNbinsX() * fDataHist->GetNbinsY(); double alpha_i = 0.0; double alpha_j = 0.0; for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { if (tempdata->GetBinContent(i + 1) && tempdata->GetBinContent(j + 1)) { alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1); alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1); } else { alpha_i = 0.0; alpha_j = 0.0; } (*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j); } } // Setup Covariances if (covar) delete covar; covar = StatUtils::GetInvert(fFullCovar); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetInvert(fFullCovar); delete tempdata; return; }; //******************************************************************** void Measurement2D::ResetFakeData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH2D *)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str()); } } //******************************************************************** void Measurement2D::ResetData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH2D *)fDataOrig->Clone((fSettings.GetName() + "_data").c_str()); } fIsFakeData = false; } //******************************************************************** void Measurement2D::ThrowCovariance() { //******************************************************************** // Take a fDecomposition and use it to throw the current dataset. // Requires fDataTrue also be set incase used repeatedly. if (fDataHist) delete fDataHist; fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); return; }; //******************************************************************** void Measurement2D::ThrowDataToy() { //******************************************************************** if (!fDataTrue) fDataTrue = (TH2D *)fDataHist->Clone(); if (fMCHist) delete fMCHist; fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); } /* Access Functions */ //******************************************************************** TH2D *Measurement2D::GetMCHistogram() { //******************************************************************** if (!fMCHist) return fMCHist; std::ostringstream chi2; chi2 << std::setprecision(5) << this->GetLikelihood(); int linecolor = kRed; int linestyle = 1; int linewidth = 1; int fillcolor = 0; int fillstyle = 1001; if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor"); if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle"); if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth"); if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor"); if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle"); fMCHist->SetTitle(chi2.str().c_str()); fMCHist->SetLineColor(linecolor); fMCHist->SetLineStyle(linestyle); fMCHist->SetLineWidth(linewidth); fMCHist->SetFillColor(fillcolor); fMCHist->SetFillStyle(fillstyle); return fMCHist; }; //******************************************************************** TH2D *Measurement2D::GetDataHistogram() { //******************************************************************** if (!fDataHist) return fDataHist; int datacolor = kBlack; int datastyle = 1; int datawidth = 1; if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor"); if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle"); if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth"); fDataHist->SetLineColor(datacolor); fDataHist->SetLineWidth(datawidth); fDataHist->SetMarkerStyle(datastyle); return fDataHist; }; /* Write Functions */ // Save all the histograms at once //******************************************************************** void Measurement2D::Write(std::string drawOpt) { //******************************************************************** // Get Draw Options drawOpt = FitPar::Config().GetParS("drawopts"); // Write Settigns if (drawOpt.find("SETTINGS") != std::string::npos) { fSettings.Set("#chi^{2}", fLikelihood); fSettings.Set("NDOF", this->GetNDOF()); fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF()); fSettings.Write(); } if (fIsChi2 && !fIsDiag) { fResidualHist = (TH2D *)fMCHist->Clone((fName + "_RESIDUAL").c_str()); fResidualHist->GetYaxis()->SetTitle("#Delta#chi^{2}"); fResidualHist->Reset(); fChi2LessBinHist = (TH2D *)fMCHist->Clone((fName + "_Chi2NMinusOne").c_str()); fChi2LessBinHist->GetYaxis()->SetTitle("Total #chi^{2} without bin_{i}"); fChi2LessBinHist->Reset(); + fIsWriting = true; (void)GetLikelihood(); + fIsWriting = false; - fResidualHist->Write(); - fChi2LessBinHist->Write(); + fResidualHist->Write((fName + "_RESIDUAL").c_str()); + fChi2LessBinHist->Write((fName + "_Chi2NMinusOne").c_str()); } // // Likelihood residual plots // if (drawOpt.find("RESIDUAL") != std::string::npos) { // WriteResidualPlots(); //} // // RATIO // if (drawOpt.find("CANVMC") != std::string::npos) { // TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist); // c1->Write(); // delete c1; // } // // PDG // if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) { // TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes); // c2->Write(); // delete c2; // } /// 2D VERSION // If null pointer return if (!fMCHist and !fDataHist) { NUIS_LOG(SAM, fName << "Incomplete histogram set!"); return; } // Config::Get().out->cd(); // Get Draw Options drawOpt = FitPar::Config().GetParS("drawopts"); bool drawData = (drawOpt.find("DATA") != std::string::npos); bool drawNormal = (drawOpt.find("MC") != std::string::npos); bool drawEvents = (drawOpt.find("EVT") != std::string::npos); bool drawFine = (drawOpt.find("FINE") != std::string::npos); bool drawRatio = (drawOpt.find("RATIO") != std::string::npos); // bool drawModes = (drawOpt.find("MODES") != std::string::npos); bool drawShape = (drawOpt.find("SHAPE") != std::string::npos); bool residual = (drawOpt.find("RESIDUAL") != std::string::npos); bool drawMatrix = (drawOpt.find("MATRIX") != std::string::npos); bool drawFlux = (drawOpt.find("FLUX") != std::string::npos); bool drawMask = (drawOpt.find("MASK") != std::string::npos); bool drawMap = (drawOpt.find("MAP") != std::string::npos); bool drawProj = (drawOpt.find("PROJ") != std::string::npos); // bool drawCanvPDG = (drawOpt.find("CANVPDG") != std::string::npos); bool drawCov = (drawOpt.find("COV") != std::string::npos); bool drawSliceMC = (drawOpt.find("CANVSLICEMC") != std::string::npos); bool drawWeighted = (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted); if (FitPar::Config().GetParB("EventManager")) { drawFlux = false; drawEvents = false; } // Save standard plots if (drawData) { GetDataList().at(0)->Write(); // Generate a simple map if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist); // Convert to 1D Lists TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist); data_1D->Write(); delete data_1D; } if (drawNormal) { GetMCList().at(0)->Write(); if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist); TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist); mc_1D->SetLineColor(kRed); mc_1D->Write(); delete mc_1D; } // Write Weighted Histogram if (drawWeighted) fMCWeighted->Write(); if (drawCov) { TH2D(*fFullCovar).Write((fName + "_COV").c_str()); } if (drawOpt.find("INVCOV") != std::string::npos) { TH2D(*covar).Write((fName + "_INVCOV").c_str()); } // Save only mc and data if splines if (fEventType == 4 or fEventType == 3) { return; } // Draw Extra plots if (drawFine) this->GetFineList().at(0)->Write(); if (drawFlux and GetFluxHistogram()) { GetFluxHistogram()->Write(); } if (drawEvents and GetEventHistogram()) { GetEventHistogram()->Write(); } if (fIsMask and drawMask) { fMaskHist->Write((fName + "_MSK").c_str()); //< save mask TH1I *mask_1D = StatUtils::MapToMask(fMaskHist, fMapHist); if (mask_1D) { mask_1D->Write(); TMatrixDSym *calc_cov = StatUtils::ApplyInvertedMatrixMasking(covar, mask_1D); TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist); TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist); TH1D *calc_data = StatUtils::ApplyHistogramMasking(data_1D, mask_1D); TH1D *calc_mc = StatUtils::ApplyHistogramMasking(mc_1D, mask_1D); TH2D *bin_cov = new TH2D(*calc_cov); bin_cov->Write(); calc_data->Write(); calc_mc->Write(); delete mask_1D; delete calc_cov; delete calc_data; delete calc_mc; delete bin_cov; delete data_1D; delete mc_1D; } } if (drawMap) fMapHist->Write((fName + "_MAP").c_str()); //< save map // // Save neut stack // if (drawModes) { // THStack combo_fMCHist_PDG = PlotUtils::GetNeutModeStack( // (fName + "_MC_PDG").c_str(), // (TH1**)fMCHist_PDG, 0); // combo_fMCHist_PDG.Write(); // } // Save Matrix plots if (drawMatrix and fFullCovar and covar and fDecomp) { TH2D cov = TH2D((*fFullCovar)); cov.SetNameTitle((fName + "_cov").c_str(), (fName + "_cov;Bins; Bins;").c_str()); cov.Write(); TH2D covinv = TH2D((*this->covar)); covinv.SetNameTitle((fName + "_covinv").c_str(), (fName + "_cov;Bins; Bins;").c_str()); covinv.Write(); TH2D covdec = TH2D((*fDecomp)); covdec.SetNameTitle((fName + "_covdec").c_str(), (fName + "_cov;Bins; Bins;").c_str()); covdec.Write(); } // Save ratio plots if required if (drawRatio) { // Needed for error bars for (int i = 0; i < fMCHist->GetNbinsX() * fMCHist->GetNbinsY(); i++) fMCHist->SetBinError(i + 1, 0.0); fDataHist->GetSumw2(); fMCHist->GetSumw2(); // Create Ratio Histograms TH2D *dataRatio = (TH2D *)fDataHist->Clone((fName + "_data_RATIO").c_str()); TH2D *mcRatio = (TH2D *)fMCHist->Clone((fName + "_MC_RATIO").c_str()); mcRatio->Divide(fMCHist); dataRatio->Divide(fMCHist); // Cancel bin errors on MC for (int i = 0; i < mcRatio->GetNbinsX() * mcRatio->GetNbinsY(); i++) { mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); } mcRatio->SetMinimum(0); mcRatio->SetMaximum(2); dataRatio->SetMinimum(0); dataRatio->SetMaximum(2); mcRatio->Write(); dataRatio->Write(); delete mcRatio; delete dataRatio; } // Save Shape Plots if required if (drawShape) { // Create Shape Histogram TH2D *mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); mcShape->SetLineWidth(3); mcShape->SetLineStyle(7); // dashes mcShape->Write(); // Save shape ratios if (drawRatio) { // Needed for error bars mcShape->GetSumw2(); // Create shape ratio histograms TH2D *mcShapeRatio = (TH2D *)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str()); TH2D *dataShapeRatio = (TH2D *)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str()); // Divide the histograms mcShapeRatio->Divide(mcShape); dataShapeRatio->Divide(mcShape); // Colour the shape ratio plots mcShapeRatio->SetLineWidth(3); mcShapeRatio->SetLineStyle(7); // dashes mcShapeRatio->Write(); dataShapeRatio->Write(); delete mcShapeRatio; delete dataShapeRatio; } delete mcShape; } // Save residual calculations of what contributed to the chi2 values. if (residual) { } if (fIsProjFitX || fIsProjFitY || drawProj) { // If not already made, make the projections if (!fMCHist_X) { PlotUtils::MatchEmptyBins(fDataHist, fMCHist); fMCHist_X = PlotUtils::GetProjectionX(fMCHist, fMaskHist); fMCHist_Y = PlotUtils::GetProjectionY(fMCHist, fMaskHist); fDataHist_X = PlotUtils::GetProjectionX(fDataHist, fMaskHist); fDataHist_Y = PlotUtils::GetProjectionY(fDataHist, fMaskHist); // This is not the correct way of doing it // double chi2X = StatUtils::GetChi2FromDiag(fDataHist_X, fMCHist_X); // double chi2Y = StatUtils::GetChi2FromDiag(fDataHist_Y, fMCHist_Y); // fMCHist_X->SetTitle(Form("%f", chi2X)); // fMCHist_Y->SetTitle(Form("%f", chi2Y)); } // Save the histograms fDataHist_X->Write(); fMCHist_X->Write(); fDataHist_Y->Write(); fMCHist_Y->Write(); } if (drawSliceMC) { TCanvas *c1 = new TCanvas((fName + "_MC_CANV_Y").c_str(), (fName + "_MC_CANV_Y").c_str(), 1024, 1024); c1->Divide(2, int(fDataHist->GetNbinsY() / 3. + 1)); TH2D *mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); double shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); mcShape->Scale(shapeScale); mcShape->SetLineStyle(7); c1->cd(1); TLegend *leg = new TLegend(0.0, 0.0, 1.0, 1.0); leg->AddEntry(fDataHist, (fName + " Data").c_str(), "lep"); leg->AddEntry(fMCHist, (fName + " MC").c_str(), "l"); leg->AddEntry(mcShape, (fName + " Shape").c_str(), "l"); leg->SetLineColor(0); leg->SetLineStyle(0); leg->SetFillColor(0); leg->SetLineStyle(0); leg->Draw("SAME"); // Make Y slices for (int i = 1; i < fDataHist->GetNbinsY() + 1; i++) { c1->cd(i + 1); TH1D *fDataHist_SliceY = PlotUtils::GetSliceY(fDataHist, i); fDataHist_SliceY->Draw("E1"); TH1D *fMCHist_SliceY = PlotUtils::GetSliceY(fMCHist, i); fMCHist_SliceY->Draw("SAME"); TH1D *mcShape_SliceY = PlotUtils::GetSliceY(mcShape, i); mcShape_SliceY->Draw("SAME"); mcShape_SliceY->SetLineStyle(mcShape->GetLineStyle()); } c1->Write(); delete c1; delete leg; TCanvas *c2 = new TCanvas((fName + "_MC_CANV_X").c_str(), (fName + "_MC_CANV_X").c_str(), 1024, 1024); c2->Divide(2, int(fDataHist->GetNbinsX() / 3. + 1)); mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); mcShape->Scale(shapeScale); mcShape->SetLineStyle(7); c2->cd(1); TLegend *leg2 = new TLegend(0.0, 0.0, 1.0, 1.0); leg2->AddEntry(fDataHist, (fName + " Data").c_str(), "lep"); leg2->AddEntry(fMCHist, (fName + " MC").c_str(), "l"); leg2->AddEntry(mcShape, (fName + " Shape").c_str(), "l"); leg2->SetLineColor(0); leg2->SetLineStyle(0); leg2->SetFillColor(0); leg2->SetLineStyle(0); leg2->Draw("SAME"); // Make Y slices for (int i = 1; i < fDataHist->GetNbinsX() + 1; i++) { c2->cd(i + 1); TH1D *fDataHist_SliceX = PlotUtils::GetSliceX(fDataHist, i); fDataHist_SliceX->Draw("E1"); TH1D *fMCHist_SliceX = PlotUtils::GetSliceX(fMCHist, i); fMCHist_SliceX->Draw("SAME"); TH1D *mcShape_SliceX = PlotUtils::GetSliceX(mcShape, i); mcShape_SliceX->Draw("SAME"); mcShape_SliceX->SetLineStyle(mcShape->GetLineStyle()); } c2->Write(); delete c2; delete leg2; } // Write Extra Histograms AutoWriteExtraTH1(); WriteExtraHistograms(); // Returning NUIS_LOG(SAM, "Written Histograms: " << fName); return; } /* Setup Functions */ //******************************************************************** void Measurement2D::SetupMeasurement(std::string inputfile, std::string type, FitWeight *rw, std::string fkdt) { //******************************************************************** // Check if name contains Evt, indicating that it is a raw number of events // measurements and should thus be treated as once fIsRawEvents = false; if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but fIsRawEvents == false!"); NUIS_LOG(SAM, "Overriding this and setting fIsRawEvents == true!"); } fIsEnu = false; if ((fName.find("XSec") != std::string::npos) && (fName.find("Enu") != std::string::npos)) { fIsEnu = true; NUIS_LOG(SAM, "::" << fName << "::"); NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " "not flux averaged!"); if (FitPar::Config().GetParB("EventManager")) { NUIS_ERR(FTL, "Enu Measurements do not yet work with the Event Manager!"); NUIS_ERR(FTL, "If you want decent flux unfolded results please run in " "series mode (-q EventManager=0)"); sleep(2); throw; } } if (fIsEnu && fIsRawEvents) { NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__); } // Reset everything to NULL fRW = rw; // Setting up 2D Inputs this->SetupInputs(inputfile); // Set Default Options SetFitOptions(fDefaultTypes); // Set Passed Options SetFitOptions(type); } //******************************************************************** void Measurement2D::SetupDefaultHist() { //******************************************************************** // Setup fMCHist fMCHist = (TH2D *)fDataHist->Clone(); fMCHist->SetNameTitle((fName + "_MC").c_str(), (fName + "_MC" + fPlotTitles).c_str()); // Setup fMCFine Int_t nBinsX = fMCHist->GetNbinsX(); Int_t nBinsY = fMCHist->GetNbinsY(); fMCFine = new TH2D((fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(), nBinsX * 3, fMCHist->GetXaxis()->GetBinLowEdge(1), fMCHist->GetXaxis()->GetBinLowEdge(nBinsX + 1), nBinsY * 3, fMCHist->GetYaxis()->GetBinLowEdge(1), fMCHist->GetYaxis()->GetBinLowEdge(nBinsY + 1)); // Setup MC Stat fMCStat = (TH2D *)fMCHist->Clone(); fMCStat->Reset(); // Setup the NEUT Mode Array // PlotUtils::CreateNeutModeArray(fMCHist, (TH1**)fMCHist_PDG); // Setup bin masks using sample name if (fIsMask) { std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask"); if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask"; } SetBinMask(maskloc); } return; } //******************************************************************** void Measurement2D::SetDataValues(std::string dataFile, std::string TH2Dname) { //******************************************************************** if (dataFile.find(".root") == std::string::npos) { NUIS_ERR(FTL, "Error! " << dataFile << " is not a .root file"); NUIS_ERR(FTL, "Currently only .root file reading is supported (MiniBooNE " "CC1pi+ 2D), but implementing .txt should be dirt easy"); NUIS_ABORT("See me at " << __FILE__ << ":" << __LINE__); } else { TFile *inFile = new TFile(dataFile.c_str(), "READ"); fDataHist = (TH2D *)(inFile->Get(TH2Dname.c_str())->Clone()); fDataHist->SetDirectory(0); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_MC" + fPlotTitles).c_str()); delete inFile; } return; } //******************************************************************** void Measurement2D::SetDataValues(std::string dataFile, double dataNorm, std::string errorFile, double errorNorm) { //******************************************************************** // Make a counter to track the line number int yBin = 0; std::string line; std::ifstream data(dataFile.c_str(), std::ifstream::in); fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(), fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins); if (data.is_open()) { NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str()); } while (std::getline(data >> std::ws, line, '\n')) { int xBin = 0; // Loop over entries and insert them into the histogram std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { fDataHist->SetBinContent(xBin + 1, yBin + 1, (*iter) * dataNorm); xBin++; } yBin++; } yBin = 0; std::ifstream error(errorFile.c_str(), std::ifstream::in); if (error.is_open()) { NUIS_LOG(SAM, "Reading errors from: " << errorFile.c_str()); } while (std::getline(error >> std::ws, line, '\n')) { int xBin = 0; // Loop over entries and insert them into the histogram std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { fDataHist->SetBinError(xBin + 1, yBin + 1, (*iter) * errorNorm); xBin++; } yBin++; } return; }; //******************************************************************** void Measurement2D::SetDataValuesFromText(std::string dataFile, double dataNorm) { //******************************************************************** fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(), fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins); NUIS_LOG(SAM, "Reading data from: " << dataFile); PlotUtils::Set2DHistFromText(dataFile, fDataHist, dataNorm, true); return; }; //******************************************************************** void Measurement2D::SetCovarMatrix(std::string covarFile) { //******************************************************************** // Used to read a covariance matrix from a root file TFile *tempFile = new TFile(covarFile.c_str(), "READ"); // Make plots that we want TH2D *covarPlot = new TH2D(); TH2D *fFullCovarPlot = new TH2D(); // Get covariance options for fake data studies std::string covName = ""; std::string covOption = FitPar::Config().GetParS("throw_covariance"); // Which matrix to get? if (fIsShape || fIsFree) covName = "shp_"; if (fIsDiag) covName += "diag"; else covName += "full"; covarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); // Throw either the sub matrix or the full matrix if (!covOption.compare("SUB")) fFullCovarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); else if (!covOption.compare("FULL")) fFullCovarPlot = (TH2D *)tempFile->Get("fullcov"); else { NUIS_ERR(WRN, " Incorrect thrown_covariance option in parameters."); } // Bin masking? int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral()); int covdim = int(fDataHist->GetNbinsX()); // Make new covars this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); fDecomp = new TMatrixDSym(dim); // Full covariance values int row, column = 0; row = 0; column = 0; for (Int_t i = 0; i < covdim; i++) { // masking can be dodgy // if (this->masked->GetBinContent(i+1) > 0) continue; for (Int_t j = 0; j < covdim; j++) { // if (this->masked->GetBinContent(j+1) > 0) continue; (*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1); (*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1); column++; } column = 0; row++; } // Set bin errors on data if (!fIsDiag) { for (Int_t i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinError( i + 1, sqrt((covarPlot->GetBinContent(i + 1, i + 1))) * 1E-38); } } TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); tempFile->Close(); delete tempFile; return; }; //******************************************************************** void Measurement2D::SetCovarMatrixFromText(std::string covarFile, int dim) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream covar(covarFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); if (covar.is_open()) { NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile); } while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation // matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { double val = (*iter) * fDataHist->GetBinError(row + 1) * 1E38 * fDataHist->GetBinError(column + 1) * 1E38; (*this->covar)(row, column) = val; (*fFullCovar)(row, column) = val; column++; } row++; } // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** void Measurement2D::SetCovarMatrixFromChol(std::string covarFile, int dim) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream covarread(covarFile.c_str(), std::ifstream::in); TMatrixD *newcov = new TMatrixD(dim, dim); if (covarread.is_open()) { NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile); } while (std::getline(covarread >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation // matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*newcov)(row, column) = *iter; column++; } row++; } covarread.close(); // Form full covariance TMatrixD *trans = (TMatrixD *)(newcov)->Clone(); trans->T(); (*trans) *= (*newcov); fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), ""); delete newcov; delete trans; // Robust matrix inversion method TDecompChol LU = TDecompChol(*this->fFullCovar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; // //******************************************************************** // void Measurement2D::SetMapValuesFromText(std::string dataFile) { // //******************************************************************** // fMapHist = new TH2I((fName + "_map").c_str(), (fName + // fPlotTitles).c_str(), // fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins); // LOG(SAM) << "Reading map from: " << dataFile << std::endl; // PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0); // return; // }; diff --git a/src/FitBase/Measurement2D.h b/src/FitBase/Measurement2D.h index 00d32f8..e540454 100644 --- a/src/FitBase/Measurement2D.h +++ b/src/FitBase/Measurement2D.h @@ -1,644 +1,644 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MEASUREMENT_2D_HXX_SEEN #define MEASUREMENT_2D_HXX_SEEN /*! * \addtogroup FitBase * @{ */ #include #include #include #include #include #include #include #include #include // ROOT includes #include #include #include #include #include #include #include #include #include // External data fit includes #include "FitEvent.h" #include "FitUtils.h" #include "MeasurementBase.h" #include "PlotUtils.h" #include "SignalDef.h" #include "StatUtils.h" #include "MeasurementVariableBox2D.h" //******************************************************************** //! 2D Measurement base class. Histogram handling is done in this base layer. class Measurement2D : public MeasurementBase { //******************************************************************** public: /* Constructor/Deconstuctor */ //! Default Constructor Measurement2D(); //! Default Destructor virtual ~Measurement2D(); /* Setup Functions */ /// \brief Setup all configs once initialised /// /// Should be called after all configs have been setup inside fSettings container. /// Handles the processing of inputs and setting up of types. /// Replaces the old 'SetupMeasurement' function. void FinaliseSampleSettings(); /// \brief Creates the 2D data distribution given the binning provided. virtual void CreateDataHistogram(int dimx, double* binx, int dimy, double* biny); /// \brief Set Data Histogram from a list of contents in a text file /// /// Assumes the format: \n /// x_low_1 y_low_1 cont_11 err_11 \n /// x_low_1 y_low_2 cont_12 err_12 \n /// x_low_2 y_low_1 cont_21 err_21 \n /// x_low_2 y_low_2 cont_22 err_22 \n /// x_low_2 y_low_3 cont_23 err_23 \n /// x_low_3 y_low_2 cont_32 err_32 \n virtual void SetDataFromTextFile(std::string data, std::string binx, std::string biny); /// \brief Set Data Histogram from a TH2D in a file /// /// - datfile = Full path to data file /// - histname = Name of histogram /// /// If histname not given it assumes that datfile /// is in the format: \n /// 'file.root;histname' virtual void SetDataFromRootFile(std::string datfile, std::string histname=""); /// \brief Set data values from a 2D array in text file /// /// \warning requires DATA HISTOGRAM TO BE SET FIRST /// /// Assumes form: \n /// cont_11 cont_12 ... cont_1N \n /// cont_21 cont_22 ... cont_2N \n /// ... ... ... ... \n /// cont_N1 cont_N2 ... cont_NN \n virtual void SetDataValuesFromTextFile(std::string datfile, TH2D* hist = NULL); /// \brief Set data errors from a 2D array in text file /// /// \warning requires DATA HISTOGRAM TO BE SET FIRST /// /// Assumes form: \n /// errs_11 errs_12 ... errs_1N \n /// errs_21 errs_22 ... errs_2N \n /// ... ... ... ... \n /// errs_N1 errs_N2 ... errs_NN \n virtual void SetDataErrorsFromTextFile(std::string datfile, TH2D* hist = NULL); /// \brief Set data bin errors to sqrt(entries) /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Sets the data errors as the sqrt of the bin contents /// Should be use for counting experiments virtual void SetPoissonErrors(); /// \brief Make diagonal covariance from data /// /// \warning If no histogram passed, data must be setup first! /// Setup the covariance inputs by taking the data histogram /// errors and setting up a diagonal covariance matrix. /// /// If no data is supplied, fDataHist is used if already set. virtual void SetCovarFromDiagonal(TH2D* data = NULL); /// \brief Read the data covariance from a text file. /// /// Inputfile should have the format: \n /// covariance_11 covariance_12 covariance_13 ... \n /// covariance_21 covariance_22 covariance_23 ... \n /// ... ... ... ... \n /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarFromTextFile(std::string covfile, int dim = -1); /// \brief Read the data covariance from a ROOT file. /// /// - covfile specifies the full path to the file /// - histname specifies the name of the covariance object. Both TMatrixDSym and TH2D are supported. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the inverted data covariance from a text file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarInvertFromTextFile(std::string covfile, int dim = -1); /// \brief Read the inverted data covariance from a ROOT file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarInvertFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the data correlations from a text file. /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCorrelationFromTextFile(std::string covfile, int dim = -1); /// \brief Read the data correlations from a ROOT file. /// /// \warning REQUIRES DATA TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCorrelationFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the cholesky decomposed covariance from a text file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCholDecompFromTextFile(std::string covfile, int dim = -1); /// \brief Read the cholesky decomposed covariance from a ROOT file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCholDecompFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the map values from a text file /// /// \warning Requires DATA hist to be set beforehand. /// Format should be a 2D array of mappings. /// -1 indicates empty bins. \n /// e.g.: \n /// 1 2 3 4 5 \n /// -1 6 7 8 9 \n /// -1 -1 10 11 -1 \n virtual void SetMapValuesFromText(std::string dataFile); /// \brief Scale the data by some scale factor virtual void ScaleData(double scale); /// \brief Scale the data error bars by some scale factor virtual void ScaleDataErrors(double scale); /// \brief Scale the covariaince and its invert/decomp by some scale factor. virtual void ScaleCovar(double scale); /// \brief Setup a bin masking histogram and apply masking to data /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Reads in a list of bins in a text file to be masked. Format is: \n /// bin_index_x_1 bin_index_y_1 1 \n /// bin_index_x_2 bin_index_y_2 1 \n /// bin_index_x_3 bin_index_y_3 1 \n /// /// If 0 is given then a bin entry will NOT be masked. So for example: \n\n /// 1 1 1 \n /// 2 0 1 \n /// 3 4 0 \n /// 4 0 1 \n\n /// Will mask only the (1,1), (2,0), and (4,0) bins. /// /// Masking can be turned on by specifiying the MASK option when creating a sample. /// When this is passed NUISANCE will look in the following locations for the mask file: /// - FitPar::Config().GetParS(fName + ".mask") /// - "data/masks/" + fName + ".mask"; virtual void SetBinMask(std::string maskfile); /// \brief Set the current fit options from a string. /// /// This is called twice for each sample, once to set the default /// and once to set the current setting (if anything other than default given) /// /// For this to work properly it requires the default and allowed types to be /// set correctly. These should be specified as a string listing options. /// /// To split up options so that NUISANCE can automatically detect ones that /// are conflicting. Any options seperated with the '/' symbol are non conflicting /// and can be given together, whereas any seperated with the ',' symbol cannot /// be specified by the end user at the same time. /// /// Default Type Examples: /// - DIAG/FIX = Default option will be a diagonal covariance, with FIXED norm. /// - MASK/SHAPE = Default option will be a masked hist, with SHAPE always on. /// /// Allowed Type examples: /// - 'FULL/DIAG/NORM/MASK' = Any of these options can be specified. /// - 'FULL,FREE,SHAPE/MASK/NORM' = User can give either FULL, FREE, or SHAPE as on option. /// MASK and NORM can also be included as options. virtual void SetFitOptions(std::string opt); /// \brief Final constructor setup /// \warning Should be called right at the end of the constructor. /// /// Contains a series of checks to ensure the data and inputs have been setup. /// Also creates the MC histograms needed for fitting. void FinaliseMeasurement(); /* Reconfigure */ /// \brief Create a Measurement1D box /// /// Creates a new 1D variable box containing just fXVar. /// /// This box is the bare minimum required by the JointFCN when /// running fast reconfigures during a routine. /// If for some reason a sample needs extra variables to be saved then /// it should override this function creating its own MeasurementVariableBox /// that contains the extra variables. virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox2D();}; /// \brief Reset all MC histograms /// /// Resets all standard histograms and those registered to auto /// process to zero. /// /// If extra histograms are not included in auto processing, then they must be reset /// by overriding this function and doing it manually if required. virtual void ResetAll(void); /// \brief Fill MC Histograms from XVar, YVar /// /// Fill standard histograms using fXVar, fYVar, Weight read from the variable box. /// /// WARNING : Any extra MC histograms need to be filled by overriding this function, /// even if they have been set to auto process. virtual void FillHistograms(void); // \brief Convert event rates to final histogram /// /// Apply standard scaling procedure to standard mc histograms to convert from /// raw events to xsec prediction. /// /// If any distributions have been set to auto process /// that is done during this function call, and a differential xsec is assumed. /// If that is not the case this function must be overriden. virtual void ScaleEvents(void); /// \brief Scale MC by a factor=1/norm /// /// Apply a simple normalisation scaling if the option FREE or a norm_parameter /// has been specified in the NUISANCE routine. virtual void ApplyNormScale(double norm); /* Statistical Functions */ /// \brief Get Number of degrees of freedom /// /// Returns the number bins inside the data histogram accounting for /// any bin masking applied. virtual int GetNDOF(); /// \brief Return Data/MC Likelihood at current state /// /// Returns the likelihood of the data given the current MC prediction. /// Diferent likelihoods definitions are used depending on the FitOptions. virtual double GetLikelihood(void); /* Fake Data */ /// \brief Set the fake data values from either a file, or MC /// /// - Setting from a file "path": \n /// When reading from a file the full path must be given to a standard /// nuisance output. The standard MC histogram should have a name that matches /// this sample for it to be read in. /// \n\n /// - Setting from "MC": \n /// If the MC option is given the current MC prediction is used as fake data. virtual void SetFakeDataValues(std::string fakeOption); /// \brief Reset fake data back to starting fake data /// /// Reset the fake data back to original fake data (Reset back to before /// ThrowCovariance was first called) virtual void ResetFakeData(void); /// \brief Reset fake data back to original data /// /// Reset the data histogram back to the true original dataset for this sample /// before any fake data was defined. virtual void ResetData(void); /// \brief Generate fake data by throwing the covariance. /// /// Can be used on fake MC data or just the original dataset. /// Call ResetFakeData or ResetData to return to values before the throw. virtual void ThrowCovariance(void); /// \brief Throw the data by its assigned errors and assign this to MC /// /// Used when creating data toys by assign the MC to this thrown data /// so that the likelihood is calculated between data and thrown data virtual void ThrowDataToy(void); /* Access Functions */ /// \brief Returns nicely formatted MC Histogram /// /// Format options can also be given in the samplesettings: /// - linecolor /// - linestyle /// - linewidth /// - fillcolor /// - fillstyle /// /// So to have a sample line colored differently in the xml cardfile put: \n /// virtual TH2D* GetMCHistogram(void); /// \brief Returns nicely formatted data Histogram /// /// Format options can also be given in the samplesettings: /// - datacolor /// - datastyle /// - datawidth /// /// So to have a sample data colored differently in the xml cardfile put: \n /// virtual TH2D* GetDataHistogram(void); /// \brief Returns a list of all MC histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMCList(void) { return std::vector(1, GetMCHistogram()); } /// \brief Returns a list of all Data histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetDataList(void) { return std::vector(1, GetDataHistogram()); } /// \brief Returns a list of all Mask histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMaskList(void) { return std::vector(1, fMaskHist); }; /// \brief Returns a list of all Fine histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetFineList(void) { return std::vector(1, fMCFine); }; /* Write Functions */ /// \brief Save the current state to the current TFile directory \n /// /// Data/MC are both saved by default. /// A range of other histograms can be saved by setting the /// config option 'drawopts'. /// /// Possible options: \n /// - FINE = Write Fine Histogram \n /// - WEIGHTS = Write Weighted MC Histogram (before scaling) \n /// - FLUX = Write Flux Histogram from MC Input \n /// - EVT = Write Event Histogram from MC Input \n /// - XSEC = Write XSec Histogram from MC Input \n /// - MASK = Write Mask Histogram \n /// - COV = Write Covariance Histogram \n /// - INVCOV = Write Inverted Covariance Histogram \n /// - DECMOP = Write Decomp. Covariance Histogram \n /// - RESIDUAL= Write Resudial Histograms \n /// - RATIO = Write Data/MC Ratio Histograms \n /// - SHAPE = Write MC Shape Histograms norm. to Data \n /// - CANVMC = Build MC Canvas Showing Data, MC, Shape \n /// - MODES = Write PDG Stack \n /// - CANVPDG = Build MC Canvas Showing Data, PDGStack \n /// /// So to save a range of these in parameters/config.xml set: \n /// virtual void Write(std::string drawOpt); //////// OLD FUNCTIONS //////////// //! Intial setup of common measurement variables. Parse input files, types, //! etc. virtual void SetupMeasurement(std::string input, std::string type, FitWeight* rw, std::string fkdt); //! Setup the default mc Hist given a data histogram virtual void SetupDefaultHist(); //! Set the data values and errors from two files virtual void SetDataValues(std::string dataFile, double dataNorm, std::string errorFile, double errorNorm); virtual void SetDataValues(std::string dataFile, std::string TH2Dname); //! Set the data values only from a text file virtual void SetDataValuesFromText(std::string dataFile, double norm); //! Read a covariance matrix from a file (Default name "covar" in file) virtual void SetCovarMatrix(std::string covarFile); //! Set the covariance matrix from a text file virtual void SetCovarMatrixFromText(std::string covarFile, int dim); //! Set the covariance matrix from a text file containing the cholesky //! fDecomposition virtual void SetCovarMatrixFromChol(std::string covarFile, int dim); protected: // The data histograms TH2D* fDataHist; //!< default data histogram (use in chi2 calculations) TH2D* fDataOrig; //!< histogram to store original data before throws. TH2D* fDataTrue; //!< histogram to store true dataset TH1D* fDataHist_X; //!< Projections onto X of the fDataHist TH1D* fDataHist_Y; //!< Projections onto Y of the fDataHist // The MC histograms TH2D* fMCHist; //!< MC Histogram (used in chi2 calculations) TH2D* fMCFine; //!< Finely binned MC Histogram TH2D* fMCHist_PDG[61]; //!< MC Histograms for each interaction mode TH1D* fMCHist_X; //!< Projections onto X of the fMCHist TH1D* fMCHist_Y; //!< Projections onto Y of the fMCHist TH2D* fMCWeighted; //!< Raw Event Weights TH2D* fMCStat; TH2I* fMaskHist; //!< mask histogram for the data TH2I* fMapHist; //!< map histogram used to convert 2D to 1D distributions TH2D *fResidualHist; TH2D *fChi2LessBinHist; bool fIsFakeData; //!< is current data actually fake std::string fakeDataFile; //!< MC fake data input file std::string fPlotTitles; //!< X and Y plot titles. std::string fFitType; std::string fDefaultTypes; //!< Default Fit Options std::string fAllowedTypes; //!< Any allowed Fit Options TMatrixDSym* covar; //!< inverted covariance matrix TMatrixDSym* fFullCovar; //!< covariance matrix TMatrixDSym* fDecomp; //!< fDecomposed covariance matrix TMatrixDSym* fCorrel; //!< correlation matrix double fCovDet; //!< covariance deteriminant double fNormError; //!< Normalisation on the error on the data double fLikelihood; //!< Likelihood value Double_t* fXBins; //!< X Bin Edges Double_t* fYBins; //!< Y Bin Edges Int_t fNDataPointsX; //!< Number of X data points Int_t fNDataPointsY; //!< NUmber of Y data points // Fit specific flags bool fIsShape; //!< Flag: Perform shape-only fit bool fIsFree; //!< Flag: Perform normalisation free fit bool fIsDiag; //!< Flag: Only use diagonal bin errors in stats bool fIsMask; //!< Flag: Apply bin masking bool fIsRawEvents; //!< Flag: Only event rates in histograms bool fIsEnu; //!< Needs Enu Unfolding bool fIsChi2SVD; //!< Flag: Chi2 SVD Method (DO NOT USE) bool fAddNormPen; //!< Flag: Add normalisation penalty to fi bool fIsProjFitX; //!< Flag: Use 1D projections onto X and Y to calculate the //!Chi2 Method. If flagged X will be used to set the rate. bool fIsProjFitY; //!< Flag: Use 1D projections onto X and Y to calculate the //!Chi2 Method. If flagged Y will be used to set the rate. bool fIsFix; //!< Flag: Fixed Histogram Norm bool fIsFull; //!< Flag; Use Full Covar bool fIsDifXSec; //!< Flag: Differential XSec bool fIsEnu1D; //!< Flag: Flux Unfolded XSec bool fIsChi2; //!< Flag; Use Chi2 over LL - + bool fIsWriting; TrueModeStack* fMCHist_Modes; ///< Optional True Mode Stack TMatrixDSym* fCovar; ///< New FullCovar TMatrixDSym* fInvert; ///< New covar // Fake Data std::string fFakeDataInput; ///< Input fake data file path TFile* fFakeDataFile; ///< Input fake data file // Arrays for data entries Double_t* fDataValues; ///< REMOVE data bin contents Double_t* fDataErrors; ///< REMOVE data bin errors }; /*! @} */ #endif diff --git a/src/InputHandler/InputHandler.cxx b/src/InputHandler/InputHandler.cxx index 9101762..304b445 100644 --- a/src/InputHandler/InputHandler.cxx +++ b/src/InputHandler/InputHandler.cxx @@ -1,313 +1,311 @@ // 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 "InputHandler.h" #include "InputUtils.h" InputHandlerBase::InputHandlerBase() { fName = ""; fFluxHist = NULL; fEventHist = NULL; fNEvents = 0; fNUISANCEEvent = NULL; fBaseEvent = NULL; kRemoveUndefParticles = FitPar::Config().GetParB("RemoveUndefParticles"); kRemoveFSIParticles = FitPar::Config().GetParB("RemoveFSIParticles"); kRemoveNuclearParticles = FitPar::Config().GetParB("RemoveNuclearParticles"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); fTTreePerformance = NULL; fSkip = 0; if (FitPar::Config().HasConfig("NSKIPEVENTS")) { fSkip = FitPar::Config().GetParI("NSKIPEVENTS"); - std::cout << "Skipping " << fSkip << " events when reading input trees." - << std::endl; } }; InputHandlerBase::~InputHandlerBase() { if (fFluxHist) delete fFluxHist; if (fEventHist) delete fEventHist; // if (fXSecHist) delete fXSecHist; // if (fNUISANCEEvent) delete fNUISANCEEvent; jointfluxinputs.clear(); jointeventinputs.clear(); jointindexlow.clear(); jointindexhigh.clear(); jointindexallowed.clear(); jointindexscale.clear(); // if (fTTreePerformance) { // fTTreePerformance->SaveAs(("ttreeperfstats_" + fName + // ".root").c_str()); // } } void InputHandlerBase::Print(){}; TH1D *InputHandlerBase::GetXSecHistogram(void) { fXSecHist = (TH1D *)fFluxHist->Clone(); fXSecHist->Divide(fEventHist); return fXSecHist; }; double InputHandlerBase::PredictedEventRate(double low, double high, std::string intOpt) { Int_t minBin = fEventHist->GetXaxis()->FindFixBin(low); Int_t maxBin = fEventHist->GetXaxis()->FindFixBin(high); if ((fEventHist->IsBinOverflow(minBin) && (low != -9999.9))) { minBin = 1; } if ((fEventHist->IsBinOverflow(maxBin) && (high != -9999.9))) { maxBin = fEventHist->GetXaxis()->GetNbins() + 1; } // If we are within a single bin if (minBin == maxBin) { // Get the contained fraction of the single bin's width return ((high - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) * fEventHist->Integral(minBin, minBin, intOpt.c_str()); } double lowBinUpEdge = fEventHist->GetXaxis()->GetBinUpEdge(minBin); double highBinLowEdge = fEventHist->GetXaxis()->GetBinLowEdge(maxBin); double lowBinfracIntegral = ((lowBinUpEdge - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) * fEventHist->Integral(minBin, minBin, intOpt.c_str()); double highBinfracIntegral = ((high - highBinLowEdge) / fEventHist->GetXaxis()->GetBinWidth(maxBin)) * fEventHist->Integral(maxBin, maxBin, intOpt.c_str()); // If they are neighbouring bins if ((minBin + 1) == maxBin) { std::cout << "Get lowfrac + highfrac" << std::endl; // Get the contained fraction of the two bin's width return lowBinfracIntegral + highBinfracIntegral; } double ContainedIntegral = fEventHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str()); // If there are filled bins between them return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral; }; double InputHandlerBase::TotalIntegratedFlux(double low, double high, std::string intOpt) { Int_t minBin = fFluxHist->GetXaxis()->FindFixBin(low); Int_t maxBin = fFluxHist->GetXaxis()->FindFixBin(high); if ((fFluxHist->IsBinOverflow(minBin) && (low != -9999.9))) { minBin = 1; } if ((fFluxHist->IsBinOverflow(maxBin) && (high != -9999.9))) { maxBin = fFluxHist->GetXaxis()->GetNbins(); high = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin + 1); } // If we are within a single bin if (minBin == maxBin) { // Get the contained fraction of the single bin's width return ((high - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) * fFluxHist->Integral(minBin, minBin, intOpt.c_str()); } double lowBinUpEdge = fFluxHist->GetXaxis()->GetBinUpEdge(minBin); double highBinLowEdge = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin); double lowBinfracIntegral = ((lowBinUpEdge - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) * fFluxHist->Integral(minBin, minBin, intOpt.c_str()); double highBinfracIntegral = ((high - highBinLowEdge) / fFluxHist->GetXaxis()->GetBinWidth(maxBin)) * fFluxHist->Integral(maxBin, maxBin, intOpt.c_str()); // If they are neighbouring bins if ((minBin + 1) == maxBin) { std::cout << "Get lowfrac + highfrac" << std::endl; // Get the contained fraction of the two bin's width return lowBinfracIntegral + highBinfracIntegral; } double ContainedIntegral = fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str()); // If there are filled bins between them return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral; } std::vector InputHandlerBase::GetFluxList(void) { return std::vector(1, fFluxHist); }; std::vector InputHandlerBase::GetEventList(void) { return std::vector(1, fEventHist); }; std::vector InputHandlerBase::GetXSecList(void) { return std::vector(1, GetXSecHistogram()); }; FitEvent *InputHandlerBase::FirstNuisanceEvent() { fCurrentIndex = 0; return GetNuisanceEvent(fCurrentIndex); }; FitEvent *InputHandlerBase::NextNuisanceEvent() { fCurrentIndex++; if ((fMaxEvents != -1) && (fCurrentIndex > fMaxEvents)) { return NULL; } return GetNuisanceEvent(fCurrentIndex); }; BaseFitEvt *InputHandlerBase::FirstBaseEvent() { fCurrentIndex = 0; return GetBaseEvent(fCurrentIndex); }; BaseFitEvt *InputHandlerBase::NextBaseEvent() { fCurrentIndex++; if (jointinput and fMaxEvents != -1) { while (fCurrentIndex < jointindexlow[jointindexswitch] || fCurrentIndex >= jointindexhigh[jointindexswitch]) { jointindexswitch++; // Loop Around if (jointindexswitch == jointindexlow.size()) { jointindexswitch = 0; } } if (fCurrentIndex > jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) { fCurrentIndex = jointindexlow[jointindexswitch]; } } return GetBaseEvent(fCurrentIndex); }; void InputHandlerBase::RegisterJointInput(std::string input, int n, TH1D *f, TH1D *e) { if (jointfluxinputs.size() == 0) { jointindexswitch = 0; fNEvents = 0; } // Push into individual input vectors jointfluxinputs.push_back((TH1D *)f->Clone()); jointeventinputs.push_back((TH1D *)e->Clone()); jointindexlow.push_back(fNEvents); jointindexhigh.push_back(fNEvents + n); fNEvents += n; // Add to the total flux/event hist if (!fFluxHist) fFluxHist = (TH1D *)f->Clone(); else fFluxHist->Add(f); if (!fEventHist) fEventHist = (TH1D *)e->Clone(); else fEventHist->Add(e); } void InputHandlerBase::SetupJointInputs() { if (jointeventinputs.size() <= 1) { jointinput = false; } else if (jointeventinputs.size() > 1) { jointinput = true; jointindexswitch = 0; } fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); if (fMaxEvents != -1 and jointeventinputs.size() > 1) { NUIS_ABORT("Can only handle joint inputs when config MAXEVENTS = -1!"); } if (jointeventinputs.size() > 1) { NUIS_ERR( WRN, "GiBUU sample contains multiple inputs. This will only work for " "samples that expect multi-species inputs. If this sample does, you " "can ignore this warning."); } for (size_t i = 0; i < jointeventinputs.size(); i++) { double scale = double(fNEvents) / fEventHist->Integral("width"); scale *= jointeventinputs.at(i)->Integral("width"); scale /= double(jointindexhigh[i] - jointindexlow[i]); jointindexscale.push_back(scale); } fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str()); fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str()); // Setup Max Events if (fMaxEvents > 1 && fMaxEvents < fNEvents) { if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl; } fNEvents = fMaxEvents; } // Print out Status if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl << "\t\t|-> Event Integral : " << fEventHist->Integral("width") * 1.E-38 << " events/nucleon" << std::endl << "\t\t|-> Flux Integral : " << fFluxHist->Integral("width") << " /cm2" << std::endl << "\t\t|-> Event/Flux : " << fEventHist->Integral("width") * 1.E-38 / fFluxHist->Integral("width") << " cm2/nucleon" << std::endl; } } BaseFitEvt *InputHandlerBase::GetBaseEvent(const UInt_t entry) { // Do some light processing: don't calculate the kinematics return static_cast(GetNuisanceEvent(entry, true)); } double InputHandlerBase::GetInputWeight(int entry) { if (!jointinput) return 1.0; // Find Switch Scale while (entry < jointindexlow[jointindexswitch] || entry >= jointindexhigh[jointindexswitch]) { jointindexswitch++; // Loop Around if (jointindexswitch >= jointindexlow.size()) { jointindexswitch = 0; } } return jointindexscale[jointindexswitch]; }; diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx index 705b97e..519685d 100644 --- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx +++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx @@ -1,263 +1,288 @@ // 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 . *******************************************************************************/ /* Authors: Adrian Orea (v1 2017) Clarence Wret (v2 2018) */ #include "MINERvA_CC0pi_XSec_2D_nu.h" #include "MINERvA_SignalDef.h" //******************************************************************** void MINERvA_CC0pi_XSec_2D_nu::SetupDataSettings() { //******************************************************************** // Define what files to use from the dist std::string datafile = ""; std::string corrfile = ""; std::string titles = ""; std::string distdescript = ""; std::string histname = ""; datafile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root"; corrfile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root"; titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} " - "(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)"; + "(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)"; distdescript = "MINERvA_CC0pi_XSec_2Dptpz_nu sample"; histname = "pt_pl_cross_section"; fSettings.SetTitle(GeneralUtils::ParseToStr(titles, ";")[0]); fSettings.SetXTitle(GeneralUtils::ParseToStr(titles, ";")[1]); fSettings.SetYTitle(GeneralUtils::ParseToStr(titles, ";")[2]); fSettings.SetZTitle(GeneralUtils::ParseToStr(titles, ";")[3]); // Sample overview --------------------------------------------------- std::string descrip = distdescript + "\n" - "Target: CH \n" - "Flux: MINERvA Low Energy FHC numu \n" - "Signal: CC-0pi \n"; + "Target: CH \n" + "Flux: MINERvA Low Energy FHC numu \n" + "Signal: CC-0pi \n"; fSettings.SetDescription(descrip); // The input ROOT file fSettings.SetDataInput(FitPar::GetDataBase() + datafile); fSettings.SetCovarInput(FitPar::GetDataBase() + corrfile); // Set the data file SetDataValues(fSettings.GetDataInput(), histname); } //******************************************************************** MINERvA_CC0pi_XSec_2D_nu::MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey) { //******************************************************************** // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); SetupDataSettings(); FinaliseSampleSettings(); fScaleFactor = - (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / - this->TotalIntegratedFlux(); + (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / + this->TotalIntegratedFlux(); TMatrixDSym *tempmat = StatUtils::GetCovarFromRootFile( fSettings.GetCovarInput(), "TotalCovariance"); fFullCovar = tempmat; // Decomposition is stable for entries that aren't E-xxx double ScalingFactor = 1E38 * 1E38; (*fFullCovar) *= ScalingFactor; // Just check that the data error and covariance are the same for (int i = 0; i < fFullCovar->GetNrows(); ++i) { for (int j = 0; j < fFullCovar->GetNcols(); ++j) { // Get the global bin int xbin1, ybin1, zbin1; fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1); double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1); double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1 + 1); double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1); double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1 + 1); if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) || ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) || xhi1 > fDataHist->GetXaxis()->GetBinLowEdge( - fDataHist->GetXaxis()->GetNbins() + 1) || + fDataHist->GetXaxis()->GetNbins() + 1) || yhi1 > fDataHist->GetYaxis()->GetBinLowEdge( - fDataHist->GetYaxis()->GetNbins() + 1)) + fDataHist->GetYaxis()->GetNbins() + 1)) continue; double data_error = fDataHist->GetBinError(xbin1, ybin1); double cov_error = sqrt((*fFullCovar)(i, i) / ScalingFactor); if (fabs(data_error - cov_error) > 1E-5) { std::cerr << "Error on data is different to that of covariance" - << std::endl; + << std::endl; NUIS_ERR(FTL, "Data error: " << data_error); NUIS_ERR(FTL, "Cov error: " << cov_error); NUIS_ERR(FTL, "Data/Cov: " << data_error / cov_error); NUIS_ERR(FTL, "Data-Cov: " << data_error - cov_error); NUIS_ERR(FTL, "For x: " << xlo1 << "-" << xhi1); NUIS_ABORT("For y: " << ylo1 << "-" << yhi1); } } } // Now can make the inverted covariance covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); // Now scale back (*fFullCovar) *= 1.0 / ScalingFactor; - (*covar) *= ScalingFactor; (*fDecomp) *= ScalingFactor; + //Don't scale this back as GetLikelihood expects it to look like this + // (*covar) *= ScalingFactor; + + fMapHist = new TH2I("MINERvA_CC0pi_XSec_2D_nu_maphist", "", + fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX(), + fDataHist->GetNbinsY(), 0, fDataHist->GetNbinsY()); + + int nbinsx = fDataHist->GetNbinsX(); + int nbinsy = fDataHist->GetNbinsY(); + Int_t Nbins = nbinsx * nbinsy; + + // Loop over the covariance matrix bins + for (int i = 0; i < Nbins; ++i) { + int xbin = (i % nbinsx) + 1; + int ybin = (i / nbinsx) + 1; + + fMapHist->SetBinContent(xbin, ybin, i+1); + } // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CC0pi_XSec_2D_nu::FillEventVariables(FitEvent *event) { //******************************************************************** // Checking to see if there is a Muon if (event->NumFSParticle(13) == 0) return; // Get the muon kinematics TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; Double_t px = Pmu.X() / 1000; Double_t py = Pmu.Y() / 1000; Double_t pt = sqrt(px * px + py * py); // y-axis is transverse momentum for both measurements fYVar = pt; // Don't want to assume the event generators all have neutrino coming along // z pz is muon momentum projected onto the neutrino direction Double_t pz = Pmu.Vect().Dot(Pnu.Vect() * (1.0 / Pnu.Vect().Mag())) / 1000.; // Set Hist Variables fXVar = pz; }; //******************************************************************** bool MINERvA_CC0pi_XSec_2D_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax); }; -//******************************************************************** -// Custom likelihood calculator because binning of covariance matrix -double MINERvA_CC0pi_XSec_2D_nu::GetLikelihood() { - //******************************************************************** - - // The calculated chi2 - double chi2 = 0.0; - - // Support shape comparisons - double scaleF = fDataHist->Integral() / fMCHist->Integral(); - if (fIsShape) { - fMCHist->Scale(scaleF); - fMCFine->Scale(scaleF); - } - - // Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA - // used for their measurement Can be prettified in due time but for now keep - - int nbinsx = fMCHist->GetNbinsX(); - int nbinsy = fMCHist->GetNbinsY(); - Int_t Nbins = nbinsx * nbinsy; - - // Loop over the covariance matrix bins - for (int i = 0; i < Nbins; ++i) { - int xbin = (i % nbinsx) + 1; - int ybin = (i / nbinsx) + 1; - double datax = fDataHist->GetBinContent(xbin, ybin); - double mcx = fMCHist->GetBinContent(xbin, ybin); - double chi2_bin = 0; - for (int j = 0; j < Nbins; ++j) { - int xbin2 = (j % nbinsx) + 1; - int ybin2 = (j / nbinsx) + 1; - - double datay = fDataHist->GetBinContent(xbin2, ybin2); - double mcy = fMCHist->GetBinContent(xbin2, ybin2); - - double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); - chi2_bin += chi2_xy; - } - if (fResidualHist) { - fResidualHist->SetBinContent(xbin, ybin, chi2_bin); - } - chi2 += chi2_bin; - } - - if (fChi2LessBinHist) { - for (int igbin = 0; igbin < Nbins; ++igbin) { - int igxbin = (igbin % nbinsx) + 1; - int igybin = (igbin / nbinsx) + 1; - double tchi2 = 0; - for (int i = 0; i < Nbins; ++i) { - int xbin = (i % nbinsx) + 1; - int ybin = (i / nbinsx) + 1; - if ((xbin == igxbin) && (ybin == igybin)) { - continue; - } - double datax = fDataHist->GetBinContent(xbin, ybin); - double mcx = fMCHist->GetBinContent(xbin, ybin); - double chi2_bin = 0; - for (int j = 0; j < Nbins; ++j) { - int xbin2 = (j % nbinsx) + 1; - int ybin2 = (j / nbinsx) + 1; - if ((xbin2 == igxbin) && (ybin2 == igybin)) { - continue; - } - double datay = fDataHist->GetBinContent(xbin2, ybin2); - double mcy = fMCHist->GetBinContent(xbin2, ybin2); - - double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); - chi2_bin += chi2_xy; - } - tchi2 += chi2_bin; - } - - fChi2LessBinHist->SetBinContent(igxbin, igybin, tchi2); - } - } - - // Normalisation penalty term if included - if (fAddNormPen) { - chi2 += - (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError); - NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * - (1 - (fCurrentNorm)) / - (fNormError * fNormError)); - } - - // Adjust the shape back to where it was. - if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) { - fMCHist->Scale(1. / scaleF); - fMCFine->Scale(1. / scaleF); - } - - fLikelihood = chi2; - - return chi2; -}; +// //******************************************************************** +// // Custom likelihood calculator because binning of covariance matrix +// double MINERvA_CC0pi_XSec_2D_nu::GetLikelihood() { +// //******************************************************************** + +// return Measurement2D::GetLikelihood(); + +// // The calculated chi2 +// double chi2 = 0.0; + +// // Support shape comparisons +// double scaleF = fDataHist->Integral() / fMCHist->Integral(); +// if (fIsShape) { +// fMCHist->Scale(scaleF); +// fMCFine->Scale(scaleF); +// } + +// int nbinsx = fMCHist->GetNbinsX(); +// int nbinsy = fMCHist->GetNbinsY(); +// Int_t Nbins = nbinsx * nbinsy; + +// // Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA +// // used for their measurement Can be prettified in due time but for now keep + +// // Loop over the covariance matrix bins +// for (int i = 0; i < Nbins; ++i) { +// int xbin = (i % nbinsx) + 1; +// int ybin = (i / nbinsx) + 1; +// double datax = fDataHist->GetBinContent(xbin, ybin); +// double mcx = fMCHist->GetBinContent(xbin, ybin); + +// // std::cout << "CMap( " << xbin << ", " << ybin << ") = " << i +// // << ", mc = " << mcx << ", mapped(" +// // << fMapHist->GetBinContent(xbin, ybin) +// // << ") = " << Mapped_MC->GetBinContent(i) << std::endl; + +// double chi2_bin = 0; +// for (int j = 0; j < Nbins; ++j) { +// int xbin2 = (j % nbinsx) + 1; +// int ybin2 = (j / nbinsx) + 1; + +// double datay = fDataHist->GetBinContent(xbin2, ybin2); +// double mcy = fMCHist->GetBinContent(xbin2, ybin2); + +// double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); +// chi2_bin += chi2_xy; +// } +// if (fResidualHist) { +// fResidualHist->SetBinContent(xbin, ybin, chi2_bin); +// } +// chi2 += chi2_bin; +// } + +// if (fChi2LessBinHist) { +// for (int igbin = 0; igbin < Nbins; ++igbin) { +// int igxbin = (igbin % nbinsx) + 1; +// int igybin = (igbin / nbinsx) + 1; +// double tchi2 = 0; +// for (int i = 0; i < Nbins; ++i) { +// int xbin = (i % nbinsx) + 1; +// int ybin = (i / nbinsx) + 1; +// if ((xbin == igxbin) && (ybin == igybin)) { +// continue; +// } +// double datax = fDataHist->GetBinContent(xbin, ybin); +// double mcx = fMCHist->GetBinContent(xbin, ybin); +// double chi2_bin = 0; +// for (int j = 0; j < Nbins; ++j) { +// int xbin2 = (j % nbinsx) + 1; +// int ybin2 = (j / nbinsx) + 1; +// if ((xbin2 == igxbin) && (ybin2 == igybin)) { +// continue; +// } +// double datay = fDataHist->GetBinContent(xbin2, ybin2); +// double mcy = fMCHist->GetBinContent(xbin2, ybin2); + +// double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); +// chi2_bin += chi2_xy; +// } +// tchi2 += chi2_bin; +// } + +// fChi2LessBinHist->SetBinContent(igxbin, igybin, tchi2); +// } +// } + +// // Normalisation penalty term if included +// if (fAddNormPen) { +// chi2 += +// (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError); +// NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * +// (1 - (fCurrentNorm)) / +// (fNormError * fNormError)); +// } + +// // Adjust the shape back to where it was. +// if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) { +// fMCHist->Scale(1. / scaleF); +// fMCFine->Scale(1. / scaleF); +// } + +// fLikelihood = chi2; + +// return chi2; +// } diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h index cead984..2219083 100644 --- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h +++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h @@ -1,50 +1,50 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MINERVA_CC0PI_XSEC_2D_NU_H_SEEN #define MINERVA_CC0PI_XSEC_2D_NU_H_SEEN #include "Measurement2D.h" -//******************************************************************** +//******************************************************************** class MINERvA_CC0pi_XSec_2D_nu : public Measurement2D { -//******************************************************************** +//******************************************************************** public: // Constructor MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey); // Destructor virtual ~MINERvA_CC0pi_XSec_2D_nu() {}; // Required functions bool isSignal(FitEvent *nvect); void FillEventVariables(FitEvent *event); - + protected: // Converted covariance matrix to provide global binning method in GetLikelihood - double GetLikelihood(); + // double GetLikelihood(); // Set up settings based on distribution void SetupDataSettings(); }; - + #endif diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_3DptpzTp_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_3DptpzTp_nu.cxx index 43e4228..042a490 100644 --- a/src/MINERvA/MINERvA_CC0pi_XSec_3DptpzTp_nu.cxx +++ b/src/MINERvA/MINERvA_CC0pi_XSec_3DptpzTp_nu.cxx @@ -1,255 +1,255 @@ // 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 . *******************************************************************************/ /* Authors: Clarence Wret (v2 2018) */ #include "MINERvA_CC0pi_XSec_3DptpzTp_nu.h" #include "MINERvA_SignalDef.h" //******************************************************************** void MINERvA_CC0pi_XSec_3DptpzTp_nu::SetupDataSettings() { //******************************************************************** // Define what files to use from the dist std::string datafile = ""; std::string corrfile = ""; std::string titles = ""; std::string distdescript = ""; std::string histname = ""; datafile = "MINERvA/CC0pi_3D/cov_ptpl_3DptpzTp_qelike.root"; corrfile = "MINERvA/CC0pi_3D/cov_ptpl_3DptpzTp_qelike.root"; titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} " "(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)"; distdescript = "MINERvA_CC0pi_XSec_3DptpzTp_nu sample"; histname = "pt_pl_cross_section"; fSettings.SetTitle(GeneralUtils::ParseToStr(titles, ";")[0]); fSettings.SetXTitle(GeneralUtils::ParseToStr(titles, ";")[1]); fSettings.SetYTitle(GeneralUtils::ParseToStr(titles, ";")[2]); fSettings.SetZTitle(GeneralUtils::ParseToStr(titles, ";")[3]); // Sample overview --------------------------------------------------- std::string descrip = distdescript + "\n" "Target: CH \n" "Flux: MINERvA Low Energy FHC numu \n" "Signal: CC-0pi \n"; fSettings.SetDescription(descrip); // The input ROOT file //fSettings.SetDataInput(FitPar::GetDataBase() + datafile); //fSettings.SetCovarInput(FitPar::GetDataBase() + corrfile); nptbins = 7; ptbins = {0, 0.15, 0.25, 0.4, 0.7, 1, 2.5}; // GeV npzbins = 4; pzbins = {1.5, 3.5, 8, 20}; // GeV ntpbins = 15; sumTpbins = {0, 40, 80, 120, 160, 200, 240, 280, 320, 360, 400, 600, 800, 1000, 10000}; // MeV // Data is actually an array of 2D measurements // Have the 2D in pt pz and then Tp to be the left over axis //fDataHist = new TH2D("minerva_test", "minerva_test", nptbins, ptbins, npzbins, pzbins); for (int i = 0; i < ntpbins; ++i) { fDataHist_Slices.push_back(new TH2D(Form("minerva_data_test_%i", i), Form("minerva_data_test_%i", i), nptbins, ptbins, npzbins, pzbins); fMCHist_Slices.push_back(new TH2D(Form("minerva_mc_test_%i", i), Form("minerva_mc_test_%i", i), nptbins, ptbins, npzbins, pzbins); } //fDataHist->SetName((fSettings.GetName() + "_data").c_str()); //fDataHist->SetTitle(fSettings.GetFullTitles().c_str()); } //******************************************************************** MINERvA_CC0pi_XSec_3DptpzTp_nu::MINERvA_CC0pi_XSec_3DptpzTp_nu(nuiskey samplekey) { //******************************************************************** // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); SetupDataSettings(); FinaliseSampleSettings(); fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CC0pi_XSec_3DptpzTp_nu::FillEventVariables(FitEvent *event) { //******************************************************************** // Checking to see if there is a Muon if (event->NumFSParticle(13) == 0) return; // Get the muon kinematics TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; Double_t px = Pmu.X() / 1000; Double_t py = Pmu.Y() / 1000; Double_t pt = sqrt(px * px + py * py); // y-axis is transverse momentum for both measurements fYVar = pt; // Don't want to assume the event generators all have neutrino coming along // z pz is muon momentum projected onto the neutrino direction Double_t pz = Pmu.Vect().Dot(Pnu.Vect() * (1.0 / Pnu.Vect().Mag())) / 1000.; // Set Hist Variables fXVar = pz; // Sum up kinetic energy of protons double sum = 0.0; - for (std::vector::iterator it = event->GetAllFSProton().begin(); + for (std::vector::iterator it = event->GetAllFSProton().begin(); it != event->GetAllFSProton().end(); ++it) { sum += (*it)->KE(); } fZVar = sum; }; void MINERvA_CC0pi_XSec_3DptpzTp_nu::FillHistograms() { Measurement2D::FillHistograms(); if (Signal) { FillMCSlice(fXVar, fYVar, fZVar, Weight); } } void MINERvA_CC0pi_XSec_3DptpzTp_nu::FillMCSlice(double x, double y, double z, double w) { // Find the bin for (int i = 0; i < ntpbins; ++i) { if (z > sumTpbins[i] && z < sumTpbins[i+1]) fMCHist_Slices[i]->Fill(y, x, w); } } //******************************************************************** bool MINERvA_CC0pi_XSec_3DptpzTp_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax); // From Dan // if not (2212 or 2112 or 22 and E<10) BAD BAD BAD }; //******************************************************************** // Custom likelihood calculator because binning of covariance matrix double MINERvA_CC0pi_XSec_3DptpzTp_nu::GetLikelihood() { //******************************************************************** // The calculated chi2 double chi2 = 0.0; // Support shape comparisons double scaleF = fDataHist->Integral() / fMCHist->Integral(); if (fIsShape) { fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); } // Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA // used for their measurement Can be prettified in due time but for now keep int nbinsx = fMCHist->GetNbinsX(); int nbinsy = fMCHist->GetNbinsY(); Int_t Nbins = nbinsx * nbinsy; // Loop over the covariance matrix bins for (int i = 0; i < Nbins; ++i) { int xbin = (i % nbinsx) + 1; int ybin = (i / nbinsx) + 1; double datax = fDataHist->GetBinContent(xbin, ybin); double mcx = fMCHist->GetBinContent(xbin, ybin); double chi2_bin = 0; for (int j = 0; j < Nbins; ++j) { int xbin2 = (j % nbinsx) + 1; int ybin2 = (j / nbinsx) + 1; double datay = fDataHist->GetBinContent(xbin2, ybin2); double mcy = fMCHist->GetBinContent(xbin2, ybin2); double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); chi2_bin += chi2_xy; } if (fResidualHist) { fResidualHist->SetBinContent(xbin, ybin, chi2_bin); } chi2 += chi2_bin; } if (fChi2LessBinHist) { for (int igbin = 0; igbin < Nbins; ++igbin) { int igxbin = (igbin % nbinsx) + 1; int igybin = (igbin / nbinsx) + 1; double tchi2 = 0; for (int i = 0; i < Nbins; ++i) { int xbin = (i % nbinsx) + 1; int ybin = (i / nbinsx) + 1; if ((xbin == igxbin) && (ybin == igybin)) { continue; } double datax = fDataHist->GetBinContent(xbin, ybin); double mcx = fMCHist->GetBinContent(xbin, ybin); double chi2_bin = 0; for (int j = 0; j < Nbins; ++j) { int xbin2 = (j % nbinsx) + 1; int ybin2 = (j / nbinsx) + 1; if ((xbin2 == igxbin) && (ybin2 == igybin)) { continue; } double datay = fDataHist->GetBinContent(xbin2, ybin2); double mcy = fMCHist->GetBinContent(xbin2, ybin2); double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy); chi2_bin += chi2_xy; } tchi2 += chi2_bin; } fChi2LessBinHist->SetBinContent(igxbin, igybin, tchi2); } } // Normalisation penalty term if included if (fAddNormPen) { chi2 += (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError); NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError)); } // Adjust the shape back to where it was. if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) { fMCHist->Scale(1. / scaleF); fMCFine->Scale(1. / scaleF); } fLikelihood = chi2; return chi2; }; diff --git a/src/Reweight/ModeNormEngine.h b/src/Reweight/ModeNormEngine.h index 124b189..e28b455 100644 --- a/src/Reweight/ModeNormEngine.h +++ b/src/Reweight/ModeNormEngine.h @@ -1,86 +1,88 @@ #ifndef ModeNormEngine_H #define ModeNormEngine_H #include "FitLogger.h" #include "GeneratorUtils.h" #include "WeightEngineBase.h" class ModeNormEngine : public WeightEngineBase { public: ModeNormEngine(std::string name = "ModeNormEngine") : fName(name){}; ~ModeNormEngine(){}; void IncludeDial(std::string name, double startval) { int rwenum = Reweight::ConvDial(name, kMODENORM); int mode = Reweight::RemoveDialType(rwenum); if (fDialEnumIndex.count(mode)) { NUIS_ABORT("Mode dial: " << mode - << " already included. Cannot include twice."); + << " already included. Cannot include twice."); } fDialEnumIndex[mode] = fDialValues.size(); fDialValues.push_back(startval); - NUIS_LOG(FIT, "Added mode dial for mode: " << mode); + NUIS_LOG(FIT, "Added mode dial for mode: " << DialToMode(mode)); } void SetDialValue(int rwenum, double val) { int mode = Reweight::RemoveDialType(rwenum); if (!fDialEnumIndex.count(mode)) { NUIS_ABORT("Mode dial: " << mode - << " has not been included. Cannot set value."); + << " has not been included. Cannot set value."); } - NUIS_LOG(DEB, "[INFO]: ModeNormEngine ObsMode: " << mode << " weight " << val - << ", rwenum = " << rwenum); + NUIS_LOG(DEB, "[INFO]: ModeNormEngine ObsMode: " + << mode << " weight " << val << ", rwenum = " << rwenum); fDialValues[fDialEnumIndex[mode]] = val; } void SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kMODENORM), val); } void Reconfigure(bool silent = false) { (void)silent; } - static int ModeToDial(int mode) { return 60 + mode; } + static int ModeToDial(int mode) { return 100 + mode; } + static int DialToMode(int dial) { return dial - 100; } double CalcWeight(BaseFitEvt *evt) { int mode = ModeToDial(abs(evt->Mode)); if (!fDialEnumIndex.count(mode)) { return 1; } NUIS_LOG(DEB, "[INFO]: Ev mode " - << evt->Mode << ", ObsMode: " << mode - << ", weight = " << fDialValues[fDialEnumIndex[mode]]); + << evt->Mode << ", ObsMode: " << mode + << ", weight = " << fDialValues[fDialEnumIndex[mode]]); return fDialValues[fDialEnumIndex[mode]]; }; bool NeedsEventReWeight() { return false; }; double GetDialValue(std::string name) { int rwenum = Reweight::ConvDial(name, kMODENORM); int mode = Reweight::RemoveDialType(rwenum); if (fDialEnumIndex.count(mode)) { return fDialValues[fDialEnumIndex[mode]]; } else { return 0xdeadbeef; } } static int SystEnumFromString(std::string const &name) { std::vector splits = GeneralUtils::ParseToStr(name, "_"); if (splits.size() != 2) { - NUIS_ABORT("Attempting to parse dial name: \"" - << name - << "\" as a mode norm dial but failed. Expect e.g. \"mode_2\"."); + NUIS_ABORT( + "Attempting to parse dial name: \"" + << name + << "\" as a mode norm dial but failed. Expect e.g. \"mode_2\"."); } int mode_num = GeneralUtils::StrToInt(splits[1]); if (!mode_num) { NUIS_ABORT("Attempting to parse dial name: \"" - << name << "\" as a mode norm dial but failed."); + << name << "\" as a mode norm dial but failed."); } return ModeToDial(mode_num); } std::map fDialEnumIndex; std::vector fDialValues; std::string fName; }; #endif diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx index a12a28a..87a41fa 100644 --- a/src/Statistical/StatUtils.cxx +++ b/src/Statistical/StatUtils.cxx @@ -1,1451 +1,1437 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "StatUtils.h" #include "GeneralUtils.h" #include "NuisConfig.h" #include "TH1D.h" //******************************************************************* Double_t StatUtils::GetChi2FromDiag(TH1D *data, TH1D *mc, TH1I *mask) { //******************************************************************* Double_t Chi2 = 0.0; TH1D *calc_data = (TH1D *)data->Clone("calc_data"); calc_data->SetDirectory(NULL); TH1D *calc_mc = (TH1D *)mc->Clone("calc_mc"); calc_mc->SetDirectory(NULL); // Add MC Error to data if required 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_data->SetDirectory(NULL); calc_mc = ApplyHistogramMasking(mc, mask); calc_mc->SetDirectory(NULL); } // 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, TH1D *outchi2perbin) { //******************************************************************* Double_t Chi2 = 0.0; TMatrixDSym *calc_cov = (TMatrixDSym *)invcov->Clone("local_invcov"); TH1D *calc_data = (TH1D *)data->Clone("local_data"); TH1D *calc_mc = (TH1D *)mc->Clone("local_mc"); calc_data->SetDirectory(NULL); calc_mc->SetDirectory(NULL); // 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")) { // 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); NUIS_LOG(FIT, "Adding cov stat " << mcerr * mcerr << " to " << (*newcov)(i, i)); (*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) NUIS_LOG(DEB, "START Chi2 Calculation================="); for (int i = 0; i < calc_data->GetNbinsX(); i++) { double ibin_contrib = 0; NUIS_LOG(DEB, "[CHI2] i = " << i << " [" << calc_data->GetXaxis()->GetBinLowEdge(i + 1) << " -- " << calc_data->GetXaxis()->GetBinUpEdge(i + 1) << "]."); for (int j = 0; j < calc_data->GetNbinsX(); j++) { NUIS_LOG(DEB, "[CHI2]\t j = " << i << " [" << calc_data->GetXaxis()->GetBinLowEdge(j + 1) << " -- " << calc_data->GetXaxis()->GetBinUpEdge(j + 1) << "]."); if ((calc_data->GetBinContent(i + 1) != 0 || calc_mc->GetBinContent(i + 1) != 0) && ((*calc_cov)(i, j) != 0)) { NUIS_LOG(DEB, "[CHI2]\t\t Chi2 contribution (i,j) = (" << i << "," << j << ")"); NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(i) = " << calc_data->GetBinContent(i + 1) << " - " << calc_mc->GetBinContent(i + 1) << " = " << (calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1))); NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(j) = " << calc_data->GetBinContent(j + 1) << " - " << calc_mc->GetBinContent(j + 1) << " = " << (calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1))); NUIS_LOG(DEB, "[CHI2]\t\t Covar = " << (*calc_cov)(i, j)); NUIS_LOG(DEB, "[CHI2]\t\t 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); double bin_cont = ((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 += bin_cont; ibin_contrib += bin_cont; } else { NUIS_LOG(DEB, "Skipping chi2 contribution (i,j) = (" << i << "," << j << "), Data = " << calc_data->GetBinContent(i + 1) << ", MC = " << calc_mc->GetBinContent(i + 1) << ", Cov = " << (*calc_cov)(i, j)); Chi2 += 0.; } } if (outchi2perbin) { outchi2perbin->SetBinContent(i + 1, ibin_contrib); } } // 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, TH2D *outchi2perbin) { //******************************************************************* // 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); - TH1D *outchi2perbin_1D = NULL; - TH1D *outchi2perbin_map_1D = NULL; - - if (outchi2perbin) { - outchi2perbin_1D = MapToTH1D(outchi2perbin, map); - for (Int_t xbi_it = 0; xbi_it < outchi2perbin->GetXaxis()->GetNbins(); - ++xbi_it) { - for (Int_t ybi_it = 0; ybi_it < outchi2perbin->GetYaxis()->GetNbins(); - ++ybi_it) { - int gbin = outchi2perbin->GetBin(xbi_it + 1, ybi_it + 1); - // std::cout << " gbin " << gbin << " corresponds to " - // << " x: " << (xbi_it + 1) << ", y: " << (ybi_it + 1) - // << std::endl; - outchi2perbin->SetBinContent(xbi_it + 1, ybi_it + 1, gbin); - } - } - outchi2perbin_map_1D = MapToTH1D(outchi2perbin, map); - } + TH1D *outchi2perbin_1D = outchi2perbin ? MapToTH1D(outchi2perbin, map) : NULL; // Calculate 1D chi2 from 1D Plots Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D, 1, 1E76, outchi2perbin_1D); - if (outchi2perbin) { - for (int xbi_it = 0; xbi_it < outchi2perbin_1D->GetXaxis()->GetNbins(); - ++xbi_it) { - // std::cout << " adding chi2 " - // << outchi2perbin_1D->GetBinContent(xbi_it + 1) - // << " from 1d bin " << (xbi_it + 1) << " to gbin " - // << outchi2perbin_map_1D->GetBinContent(xbi_it + 1) << - // std::endl; - outchi2perbin->SetBinContent( - outchi2perbin_map_1D->GetBinContent(xbi_it + 1), - outchi2perbin_1D->GetBinContent(xbi_it + 1)); - } + if (outchi2perbin && outchi2perbin_1D) { + MapFromTH1D(outchi2perbin, outchi2perbin_1D, map); } // CleanUp delete data_1D; delete mc_1D; delete mask_1D; delete outchi2perbin_1D; - delete outchi2perbin_map_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 = NULL; 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)) { NUIS_LOG(DEB, "Applying mask to bin " << i + 1 << " " << hist->GetName()); 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 *DataHist, TMatrixDSym *cov, double scale, bool ErrorCheck) { //******************************************************************* // Check if (ErrorCheck) { if (cov->GetNrows() != DataHist->GetNbinsX()) { NUIS_ERR( FTL, "Nrows in cov don't match nbins in DataHist for SetDataErrorFromCov"); NUIS_ERR(FTL, "Nrows = " << cov->GetNrows()); NUIS_ABORT("Nbins = " << DataHist->GetNbinsX()); } } // Set bin errors form cov diag // Check if the errors are set bool ErrorsSet = false; for (int i = 0; i < DataHist->GetNbinsX(); i++) { if (ErrorsSet == true) break; if (DataHist->GetBinError(i + 1) != 0 && DataHist->GetBinContent(i + 1) > 0) ErrorsSet = true; } // Now loop over if (ErrorsSet && ErrorCheck) { for (int i = 0; i < DataHist->GetNbinsX(); i++) { double DataHisterr = DataHist->GetBinError(i + 1); double coverr = sqrt((*cov)(i, i)) * scale; // Check that the errors are within 1% of eachother if (fabs(DataHisterr - coverr) / DataHisterr > 0.01) { NUIS_ERR(WRN, "Data error does not match covariance error for bin " << i + 1 << " (" << DataHist->GetXaxis()->GetBinLowEdge(i + 1) << "-" << DataHist->GetXaxis()->GetBinLowEdge(i + 2) << ")"); NUIS_ERR(WRN, "Data error: " << DataHisterr); NUIS_ERR(WRN, "Cov error: " << coverr); } } // Else blindly trust the covariance } else { for (int i = 0; i < DataHist->GetNbinsX(); i++) { DataHist->SetBinError(i + 1, sqrt((*cov)(i, i)) * scale); } } return; } //******************************************************************* void StatUtils::SetDataErrorFromCov(TH2D *data, TMatrixDSym *cov, TH2I *map, double scale, bool ErrorCheck) { //******************************************************************* // Check if (ErrorCheck) { if (cov->GetNrows() != data->GetNbinsX() * data->GetNbinsY()) { NUIS_ERR(FTL, "Nrows in cov don't match nbins in data for " "SetDataNUIS_ERRorFromCov"); NUIS_ERR(FTL, "Nrows = " << cov->GetNrows()); NUIS_ABORT("Nbins = " << data->GetNbinsX()); } } // Set bin errors form cov diag // Check if the errors are set bool ErrorsSet = false; for (int i = 0; i < data->GetNbinsX(); i++) { for (int j = 0; j < data->GetNbinsX(); j++) { if (ErrorsSet == true) break; if (data->GetBinError(i + 1, j + 1) != 0) ErrorsSet = true; } } // 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; // If we have errors on our histogram the map is good count = map->GetBinContent(i + 1, j + 1) - 1; double dataerr = data->GetBinError(i + 1, j + 1); double coverr = sqrt((*cov)(count, count)) * scale; // Check that the errors are within 1% of eachother if (ErrorsSet && ErrorCheck) { if (fabs(dataerr - coverr) / dataerr > 0.01) { NUIS_ERR(WRN, "Data error does not match covariance error for bin " << i + 1 << " (" << data->GetXaxis()->GetBinLowEdge(i + 1) << "-" << data->GetXaxis()->GetBinLowEdge(i + 2) << ")"); NUIS_ERR(WRN, "Data error: " << dataerr); NUIS_ERR(WRN, "Cov error: " << coverr); } } else { data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale); } } } // Delete the map now that we don't need it // Woops, it's needed elsewhere there! (Grrrrr) // map->Delete(); } 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) { NUIS_ERR(WRN, "Inconsistent matrix and data histogram passed to " "StatUtils::ExtractShapeOnlyCovar!"); NUIS_ERR(WRN, "data_hist has " << data_hist->GetNbinsX() << " matrix has " << nbins); int err_bins = data_hist->GetNbinsX(); if (nbins > err_bins) err_bins = nbins; for (int i = 0; i < err_bins; ++i) { NUIS_ERR(WRN, "Matrix diag. = " << (*full_covar)(i, i) << " data = " << data_hist->GetBinContent(i + 1)); } 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) { NUIS_ERR(WRN, "Stupid matrix or data histogram passed to " "StatUtils::ExtractShapeOnlyCovar! Ignoring..."); return NULL; } NUIS_LOG(SAM, "Norm error = " << sqrt(total_covar) / total_data); // 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) { 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; } +void StatUtils::MapFromTH1D(TH2 *fillhist, TH1 *fromhist, TH2I *map) { + fillhist->Clear(); + + 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; + int gb = map->GetBinContent(i + 1, j + 1); + fillhist->SetBinContent(i + 1, j + 1, fromhist->GetBinContent(gb)); + fillhist->SetBinError(i + 1, j + 1, fromhist->GetBinError(gb)); + } + } +} + //******************************************************************* 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) { NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " "entries on this line: " << row); } 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) { NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " "entries on this line: " << row); } 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) { NUIS_ABORT("No object name given!"); } // Get file TFile *tempfile = new TFile(splitfile[0].c_str(), "READ"); // Get Object TObject *obj = tempfile->Get(splitfile[1].c_str()); if (!obj) { NUIS_ABORT("Object " << splitfile[1] << " doesn't exist!"); } // 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; } diff --git a/src/Statistical/StatUtils.h b/src/Statistical/StatUtils.h index d2f82c9..0e9e37b 100644 --- a/src/Statistical/StatUtils.h +++ b/src/Statistical/StatUtils.h @@ -1,255 +1,258 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef STATUTILS_H #define STATUTILS_H // C Includes #include #include #include #include #include #include #include #include #include "assert.h" // Root Includes #include "TH1D.h" #include "TH2I.h" #include "TH2D.h" #include "TFile.h" #include "TMatrixDSym.h" #include "TDecompSVD.h" #include "TMath.h" #include "TRandom3.h" #include "TDecompChol.h" #include "TGraphErrors.h" // Fit Includes #include "FitLogger.h" /*! * \addtogroup Utils * @{ */ //! Functions for handling statistics calculations namespace StatUtils{ /* Chi2 Functions */ //! Get Chi2 using diagonal bin errors from the histogram. Masking applied before calculation if mask provided. Double_t GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Get Chi2 using diagonal bin errors from the histogram. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromDiag(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); //! Get Chi2 using an inverted covariance for the data Double_t GetChi2FromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL, double data_scale=1, double covar_scale=1E76, TH1D *outchi2perbin=NULL); //! Get Chi2 using an inverted covariance for the data //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL, TH2D *outchi2perbin=NULL); //! Get Chi2 using an SVD method on the covariance before calculation. //! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work. Double_t GetChi2FromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask=NULL); //! Get Chi2 using an SVD method on the covariance before calculation. //! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromSVD( TH2D* data, TH2D* mc, TMatrixDSym* cov, TH2I* map=NULL, TH2I* mask=NULL); //! Get Chi2 using only the raw event rates given in each bin using a -2LL method. Double_t GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Get Chi2 using only the raw event rates given in each bin using a -2LL method. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromEventRate(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); // Likelihood Functions //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromDiag(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromDiag(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromSVD( TH2D* data, TH2D* mc, TMatrixDSym* cov, TH2I* map=NULL, TH2I* mask=NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromEventRate(TH1D* data, TH1D* mc, TH1I* mask=NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromEventRate(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL); /* NDOF Functions */ //! Return 1D Histogram NDOF considering masking and empty bins Int_t GetNDOF(TH1D* hist, TH1I* mask=NULL); //! Return 2D Histogram NDOF considering masking and empty bins Int_t GetNDOF(TH2D* hist, TH2I* map=NULL, TH2I* mask=NULL); /* Fake Data Functions */ //! Given a full covariance for a 1D data set throw the decomposition to generate fake data. //! throwdiag determines whether diagonal statistical errors are thrown. //! If no covariance is provided only statistical errors are thrown. TH1D* ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag=true, TH1I* mask=NULL); //! Given a full covariance for a 2D data set throw the decomposition to generate fake data. //! Plots are converted to 1D histograms and the 1D ThrowHistogram is used, before being converted back to 2D histograms. TH2D* ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map=NULL, bool throwdiag=true, TH2I* mask=NULL); /* Masking Functions */ //! Given a mask histogram, mask out any bins in hist with non zero entries in mask. TH1D* ApplyHistogramMasking(TH1D* hist, TH1I* mask); //! Given a mask histogram, mask out any bins in hist with non zero entries in mask. TH2D* ApplyHistogramMasking(TH2D* hist, TH2I* mask); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance, before recalculating its inverse. TMatrixDSym* ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH1I* mask); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance, before recalculating its inverse. //! Converts to 1D data before using the 1D ApplyInvertedMatrixMasking function and converting back to 2D. TMatrixDSym* ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map=NULL); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance TMatrixDSym* ApplyMatrixMasking(TMatrixDSym* mat, TH1I* mask); //! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance //! Converts to 1D data before using the 1D ApplyInvertedMatrixMasking function and converting back to 2D. TMatrixDSym* ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map=NULL); /* Covariance Handling Functions */ //! Return inverted matrix of TMatrixDSym TMatrixDSym* GetInvert(TMatrixDSym* mat); //! Return Cholesky Decomposed matrix of TMatrixDSym TMatrixDSym* GetDecomp(TMatrixDSym* mat); //! Return full covariances TMatrixDSym* GetCovarFromCorrel(TMatrixDSym* correl, TH1D* data); //! Given a normalisation factor for a dataset add in a new normalisation term to the covariance. void ForceNormIntoCovar(TMatrixDSym*& mat, TH1D* data, double norm); //! Given a normalisation factor for a dataset add in a new normalisation term to the covariance. //! Convertes 2D to 1D, before using 1D ForceNormIntoCovar void ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, TH2I* map=NULL); //! Given a dataset generate an uncorrelated covariance matrix using the bin errors. TMatrixDSym* MakeDiagonalCovarMatrix(TH1D* data, double scaleF=1E38); //! Given a dataset generate an uncorrelated covariance matrix using the bin errors. TMatrixDSym* MakeDiagonalCovarMatrix(TH2D* data, TH2I* map=NULL, double scaleF=1E38); //! Given a covariance set the errors in each bin on the data from the covariance diagonals. void SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, double scale=1.0, bool ErrorCheck = true); //! Given a covariance set the errors in each bin on the data from the covariance diagonals. void SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map=NULL, double scale=1.0, bool ErrorCheck = true); //! Given a covariance, extracts the shape-only matrix using the method from the MiniBooNE TN TMatrixDSym* ExtractShapeOnlyCovar(TMatrixDSym* full_covar, TH1* data_hist, double data_scale=1E38); /* Mapping Functions */ //! If no map is provided for the 2D histogram, generate one by counting up through the bins along x and y. TH2I* GenerateMap(TH2D* hist); //! Apply a map to a 2D histogram converting it into a 1D histogram. TH1D* MapToTH1D(TH2D* hist, TH2I* map); + //! Apply a map to fill a TH2D from a TH1D; + void MapFromTH1D(TH2* fillhist, TH1* fromhist, TH2I* map); + //! Apply a map to a 2D mask convering it into a 1D mask. TH1I* MapToMask(TH2I* hist, TH2I* map); /// \brief Read TMatrixD from a text file /// /// - covfile = full path to text file /// - dimx = x dimensions of matrix /// - dimy = y dimensions of matrix /// /// Format of textfile should be: \n /// cov_11 cov_12 ... cov_1N \n /// cov_21 cov_22 ... cov_2N \n /// ... ... ... ... \n /// cov_N1 ... ... cov_NN \n /// /// If no dimensions are given, dimx and dimy are determined from rows/columns /// inside textfile. /// /// If only dimx is given a symmetric matrix is assumed. TMatrixD* GetMatrixFromTextFile(std::string covfile, int dimx=-1, int dimy=-1); /// \brief Read TMatrixD from a ROOT file /// /// - covfile = full path to root file (+';histogram') /// - histname = histogram name /// /// If no histogram name is given function assumes it has been appended /// covfile path as: \n /// 'covfile.root;histname' /// /// histname can point to a TMatrixD object, a TMatrixDSym object, or /// a TH2D object. TMatrixD* GetMatrixFromRootFile(std::string covfile, std::string histname=""); /// \brief Calls GetMatrixFromTextFile and turns it into a TMatrixDSym TMatrixDSym* GetCovarFromTextFile(std::string covfile, int dim); /// \brief Calls GetMatrixFromRootFile and turns it into a TMatrixDSym TMatrixDSym* GetCovarFromRootFile(std::string covfile, std::string histname); }; /*! @} */ #endif diff --git a/src/Utils/PlotUtils.cxx b/src/Utils/PlotUtils.cxx index 89dc4f2..612c493 100644 --- a/src/Utils/PlotUtils.cxx +++ b/src/Utils/PlotUtils.cxx @@ -1,1197 +1,1206 @@ // 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 "PlotUtils.h" #include "FitEvent.h" #include "StatUtils.h" // MOVE TO MB UTILS! // This function is intended to be modified to enforce a consistent masking for // all models. TH2D *PlotUtils::SetMaskHist(std::string type, TH2D *data) { TH2D *fMaskHist = (TH2D *)data->Clone("fMaskHist"); for (int xBin = 0; xBin < fMaskHist->GetNbinsX(); ++xBin) { for (int yBin = 0; yBin < fMaskHist->GetNbinsY(); ++yBin) { if (data->GetBinContent(xBin + 1, yBin + 1) == 0) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0); else fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0.5); if (!type.compare("MB_numu_2D")) { if (yBin == 19 && xBin < 8) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0); } else { if (yBin == 19 && xBin < 11) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0); } if (yBin == 18 && xBin < 3) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0); if (yBin == 17 && xBin < 2) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0); if (yBin == 16 && xBin < 1) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0); } } return fMaskHist; }; // MOVE TO GENERAL UTILS? bool PlotUtils::CheckObjectWithName(TFile *inFile, std::string substring) { TIter nextkey(inFile->GetListOfKeys()); TKey *key; while ((key = (TKey *)nextkey())) { std::string test(key->GetName()); if (test.find(substring) != std::string::npos) return true; } return false; }; // MOVE TO GENERAL UTILS? std::string PlotUtils::GetObjectWithName(TFile *inFile, std::string substring) { TIter nextkey(inFile->GetListOfKeys()); TKey *key; std::string output = ""; while ((key = (TKey *)nextkey())) { std::string test(key->GetName()); if (test.find(substring) != std::string::npos) output = test; } return output; }; void PlotUtils::CreateNeutModeArray(TH1 *hist, TH1 *neutarray[]) { for (int i = 0; i < 60; i++) { neutarray[i] = (TH1 *)hist->Clone(Form("%s_NMODE_%i", hist->GetName(), i)); } return; }; void PlotUtils::DeleteNeutModeArray(TH1 *neutarray[]) { for (int i = 0; i < 60; i++) { delete neutarray[i]; } return; }; void PlotUtils::FillNeutModeArray(TH1D *hist[], int mode, double xval, double weight) { if (abs(mode) > 60) return; hist[abs(mode)]->Fill(xval, weight); return; }; void PlotUtils::FillNeutModeArray(TH2D *hist[], int mode, double xval, double yval, double weight) { if (abs(mode) > 60) return; hist[abs(mode)]->Fill(xval, yval, weight); return; }; THStack PlotUtils::GetNeutModeStack(std::string title, TH1 *ModeStack[], int option) { (void)option; THStack allmodes = THStack(title.c_str(), title.c_str()); for (int i = 0; i < 60; i++) { allmodes.Add(ModeStack[i]); } // Credit to Clarence for copying all this out. // CC ModeStack[1]->SetTitle("CCQE"); ModeStack[1]->SetFillColor(kBlue); // ModeStack[1]->SetFillStyle(3444); ModeStack[1]->SetLineColor(kBlue); ModeStack[2]->SetTitle("2p/2h Nieves"); ModeStack[2]->SetFillColor(kRed); // ModeStack[2]->SetFillStyle(3344); ModeStack[2]->SetLineColor(kRed); // ModeStack[11]->SetTitle("#it{#nu + p #rightarrow l^{-} + p + #pi^{+}}"); ModeStack[11]->SetTitle("CC1#pi^{+} on p"); ModeStack[11]->SetFillColor(kGreen); // ModeStack[11]->SetFillStyle(3004); ModeStack[11]->SetLineColor(kGreen); // ModeStack[12]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #pi^{0}}"); ModeStack[12]->SetTitle("CC1#pi^{0} on n"); ModeStack[12]->SetFillColor(kGreen + 3); // ModeStack[12]->SetFillStyle(3005); ModeStack[12]->SetLineColor(kGreen); // ModeStack[13]->SetTitle("#it{#nu + n #rightarrow l^{-} + n + #pi^{+}}"); ModeStack[13]->SetTitle("CC1#pi^{+} on n"); ModeStack[13]->SetFillColor(kGreen - 2); // ModeStack[13]->SetFillStyle(3004); ModeStack[13]->SetLineColor(kGreen); ModeStack[16]->SetTitle("CC coherent"); ModeStack[16]->SetFillColor(kBlue); // ModeStack[16]->SetFillStyle(3644); ModeStack[16]->SetLineColor(kBlue); // ModeStack[17]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #gamma}"); ModeStack[17]->SetTitle("CC1#gamma"); ModeStack[17]->SetFillColor(kMagenta); // ModeStack[17]->SetFillStyle(3001); ModeStack[17]->SetLineColor(kMagenta); ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)"); ModeStack[21]->SetFillColor(kYellow); // ModeStack[21]->SetFillStyle(3005); ModeStack[21]->SetLineColor(kYellow); // ModeStack[22]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #eta^{0}}"); ModeStack[22]->SetTitle("CC1#eta^{0} on n"); ModeStack[22]->SetFillColor(kYellow - 2); // ModeStack[22]->SetFillStyle(3013); ModeStack[22]->SetLineColor(kYellow - 2); // ModeStack[23]->SetTitle("#it{#nu + n #rightarrow l^{-} + #Lambda + // K^{+}}"); ModeStack[23]->SetTitle("CC1#Labda1K^{+}"); ModeStack[23]->SetFillColor(kYellow - 6); // ModeStack[23]->SetFillStyle(3013); ModeStack[23]->SetLineColor(kYellow - 6); ModeStack[26]->SetTitle("DIS (W > 2.0)"); ModeStack[26]->SetFillColor(kRed); // ModeStack[26]->SetFillStyle(3006); ModeStack[26]->SetLineColor(kRed); // NC // ModeStack[31]->SetTitle("#it{#nu + n #rightarrow #nu + n + #pi^{0}}"); ModeStack[31]->SetTitle("NC1#pi^{0} on n"); ModeStack[31]->SetFillColor(kBlue); // ModeStack[31]->SetFillStyle(3004); ModeStack[31]->SetLineColor(kBlue); // ModeStack[32]->SetTitle("#it{#nu + p #rightarrow #nu + p + #pi^{0}}"); ModeStack[32]->SetTitle("NC1#pi^{0} on p"); ModeStack[32]->SetFillColor(kBlue + 3); // ModeStack[32]->SetFillStyle(3004); ModeStack[32]->SetLineColor(kBlue + 3); // ModeStack[33]->SetTitle("#it{#nu + n #rightarrow #nu + p + #pi^{-}}"); ModeStack[33]->SetTitle("NC1#pi^{-} on n"); ModeStack[33]->SetFillColor(kBlue - 2); // ModeStack[33]->SetFillStyle(3005); ModeStack[33]->SetLineColor(kBlue - 2); // ModeStack[34]->SetTitle("#it{#nu + p #rightarrow #nu + n + #pi^{+}}"); ModeStack[34]->SetTitle("NC1#pi^{+} on p"); ModeStack[34]->SetFillColor(kBlue - 8); // ModeStack[34]->SetFillStyle(3005); ModeStack[34]->SetLineColor(kBlue - 8); ModeStack[36]->SetTitle("NC Coherent"); ModeStack[36]->SetFillColor(kBlue + 8); // ModeStack[36]->SetFillStyle(3644); ModeStack[36]->SetLineColor(kBlue + 8); // ModeStack[38]->SetTitle("#it{#nu + n #rightarrow #nu + n + #gamma}"); ModeStack[38]->SetTitle("NC1#gamma on n"); ModeStack[38]->SetFillColor(kMagenta); // ModeStack[38]->SetFillStyle(3001); ModeStack[38]->SetLineColor(kMagenta); // ModeStack[39]->SetTitle("#it{#nu + p #rightarrow #nu + p + #gamma}"); ModeStack[39]->SetTitle("NC1#gamma on p"); ModeStack[39]->SetFillColor(kMagenta - 10); // ModeStack[39]->SetFillStyle(3001); ModeStack[39]->SetLineColor(kMagenta - 10); ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)"); ModeStack[41]->SetFillColor(kBlue - 10); // ModeStack[41]->SetFillStyle(3005); ModeStack[41]->SetLineColor(kBlue - 10); // ModeStack[42]->SetTitle("#it{#nu + n #rightarrow #nu + n + #eta^{0}}"); ModeStack[42]->SetTitle("NC1#eta^{0} on n"); ModeStack[42]->SetFillColor(kYellow - 2); // ModeStack[42]->SetFillStyle(3013); ModeStack[42]->SetLineColor(kYellow - 2); // ModeStack[43]->SetTitle("#it{#nu + p #rightarrow #nu + p + #eta^{0}}"); ModeStack[43]->SetTitle("NC1#eta^{0} on p"); ModeStack[43]->SetFillColor(kYellow - 4); // ModeStack[43]->SetFillStyle(3013); ModeStack[43]->SetLineColor(kYellow - 4); // ModeStack[44]->SetTitle("#it{#nu + n #rightarrow #nu + #Lambda + K^{0}}"); ModeStack[44]->SetTitle("NC1#Lambda1K^{0} on n"); ModeStack[44]->SetFillColor(kYellow - 6); // ModeStack[44]->SetFillStyle(3014); ModeStack[44]->SetLineColor(kYellow - 6); // ModeStack[45]->SetTitle("#it{#nu + p #rightarrow #nu + #Lambda + K^{+}}"); ModeStack[45]->SetTitle("NC1#Lambda1K^{+}"); ModeStack[45]->SetFillColor(kYellow - 10); // ModeStack[45]->SetFillStyle(3014); ModeStack[45]->SetLineColor(kYellow - 10); ModeStack[46]->SetTitle("DIS (W > 2.0)"); ModeStack[46]->SetFillColor(kRed); // ModeStack[46]->SetFillStyle(3006); ModeStack[46]->SetLineColor(kRed); // ModeStack[51]->SetTitle("#it{#nu + p #rightarrow #nu + p}"); ModeStack[51]->SetTitle("NC on p"); ModeStack[51]->SetFillColor(kBlack); // ModeStack[51]->SetFillStyle(3444); ModeStack[51]->SetLineColor(kBlack); // ModeStack[52]->SetTitle("#it{#nu + n #rightarrow #nu + n}"); ModeStack[52]->SetTitle("NC on n"); ModeStack[52]->SetFillColor(kGray); // ModeStack[52]->SetFillStyle(3444); ModeStack[52]->SetLineColor(kGray); return allmodes; }; TLegend PlotUtils::GenerateStackLegend(THStack stack, int xlow, int ylow, int xhigh, int yhigh) { TLegend leg = TLegend(xlow, ylow, xhigh, yhigh); TObjArray *histarray = stack.GetStack(); int nhist = histarray->GetEntries(); for (int i = 0; i < nhist; i++) { TH1 *hist = (TH1 *)(histarray->At(i)); leg.AddEntry((hist), ((TH1 *)histarray->At(i))->GetTitle(), "fl"); } leg.SetName(Form("%s_LEG", stack.GetName())); return leg; }; void PlotUtils::ScaleNeutModeArray(TH1 *hist[], double factor, std::string option) { for (int i = 0; i < 60; i++) { if (hist[i]) hist[i]->Scale(factor, option.c_str()); } return; }; void PlotUtils::ResetNeutModeArray(TH1 *hist[]) { for (int i = 0; i < 60; i++) { if (hist[i]) hist[i]->Reset(); } return; }; //******************************************************************** // This assumes the Enu axis is the x axis, as is the case for MiniBooNE 2D // distributions void PlotUtils::FluxUnfoldedScaling(TH2D *fMCHist, TH1D *fhist, TH1D *ehist, double scalefactor) { //******************************************************************** // Make clones to avoid changing stuff TH1D *eventhist = (TH1D *)ehist->Clone(); TH1D *fFluxHist = (TH1D *)fhist->Clone(); // Undo width integral in SF fMCHist->Scale(scalefactor / eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width")); // Standardise The Flux eventhist->Scale(1.0 / fFluxHist->Integral()); fFluxHist->Scale(1.0 / fFluxHist->Integral()); // Do interpolation for 2D plots? // fFluxHist = PlotUtils::InterpolateFineHistogram(fFluxHist,100,"width"); // eventhist = PlotUtils::InterpolateFineHistogram(eventhist,100,"width"); // eventhist->Scale(1.0/fFluxHist->Integral()); // fFluxHist->Scale(1.0/fFluxHist->Integral()); // Scale fMCHist by eventhist integral fMCHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1)); // Find which axis is the Enu axis bool EnuOnXaxis = false; std::string xaxis = fMCHist->GetXaxis()->GetTitle(); if (xaxis.find("E") != std::string::npos && xaxis.find("nu") != std::string::npos) EnuOnXaxis = true; std::string yaxis = fMCHist->GetYaxis()->GetTitle(); if (yaxis.find("E") != std::string::npos && xaxis.find("nu") != std::string::npos) { // First check that xaxis didn't also find Enu if (EnuOnXaxis) { NUIS_ERR(FTL, fMCHist->GetTitle() << " error:"); NUIS_ERR(FTL, "Found Enu in xaxis title: " << xaxis); NUIS_ERR(FTL, "AND"); NUIS_ERR(FTL, "Found Enu in yaxis title: " << yaxis); NUIS_ABORT("Enu on x and Enu on y flux unfolded scaling isn't " - "implemented, please modify " - << __FILE__ << ":" << __LINE__); + "implemented, please modify " + << __FILE__ << ":" << __LINE__); } EnuOnXaxis = false; } // Now Get a flux PDF assuming X axis is Enu TH1D *pdfflux = NULL; // If xaxis is Enu if (EnuOnXaxis) pdfflux = (TH1D *)fMCHist->ProjectionX()->Clone(); // If yaxis is Enu else pdfflux = (TH1D *)fMCHist->ProjectionY()->Clone(); // pdfflux->Write( (std::string(fMCHist->GetName()) + "_PROJX").c_str()); pdfflux->Reset(); // Awful MiniBooNE Check for the time being // Needed because the flux is in GeV whereas the measurement is in MeV bool ismb = std::string(fMCHist->GetName()).find("MiniBooNE") != std::string::npos; for (int i = 0; i < pdfflux->GetNbinsX(); i++) { double Ml = pdfflux->GetXaxis()->GetBinLowEdge(i + 1); double Mh = pdfflux->GetXaxis()->GetBinLowEdge(i + 2); // double Mc = pdfflux->GetXaxis()->GetBinCenter(i+1); // double Mw = pdfflux->GetBinWidth(i+1); double fluxint = 0.0; // Scaling to match flux for MB if (ismb) { Ml /= 1.E3; Mh /= 1.E3; // Mc /= 1.E3; // Mw /= 1.E3; } for (int j = 0; j < fFluxHist->GetNbinsX(); j++) { // double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1); double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1); double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2); double Fe = fFluxHist->GetBinContent(j + 1); double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1); if (Fl >= Ml and Fh <= Mh) { fluxint += Fe; } else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) { fluxint += Fe * (Fh - Ml) / Fw; } else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) { fluxint += Fe * (Mh - Fl) / Fw; } else if (Ml >= Fl and Mh <= Fh) { fluxint += Fe * (Mh - Ml) / Fw; } else { continue; } } pdfflux->SetBinContent(i + 1, fluxint); } // Then finally divide by the bin-width in for (int i = 0; i < fMCHist->GetNbinsX(); i++) { for (int j = 0; j < fMCHist->GetNbinsY(); j++) { if (pdfflux->GetBinContent(i + 1) == 0.0) continue; // Different scaling depending on if Enu is on x or y axis double scaling = 1.0; // If Enu is on the x-axis, we want the ith entry of the flux // And to divide by the bin width of the jth bin if (EnuOnXaxis) { double binWidth = fMCHist->GetYaxis()->GetBinLowEdge(j + 2) - fMCHist->GetYaxis()->GetBinLowEdge(j + 1); scaling = pdfflux->GetBinContent(i + 1) * binWidth; } else { double binWidth = fMCHist->GetXaxis()->GetBinLowEdge(i + 2) - fMCHist->GetXaxis()->GetBinLowEdge(i + 1); scaling = pdfflux->GetBinContent(j + 1) * binWidth; } // fMCHist->SetBinContent(i + 1, j + 1, // fMCHist->GetBinContent(i + 1, j + 1) / // pdfflux->GetBinContent(i + 1) / binWidth); // fMCHist->SetBinError(i + 1, j + 1, fMCHist->GetBinError(i + 1, j + 1) / // pdfflux->GetBinContent(i + 1) / // binWidth); fMCHist->SetBinContent(i + 1, j + 1, fMCHist->GetBinContent(i + 1, j + 1) / scaling); fMCHist->SetBinError(i + 1, j + 1, fMCHist->GetBinError(i + 1, j + 1) / scaling); } } delete eventhist; delete fFluxHist; }; TH1D *PlotUtils::InterpolateFineHistogram(TH1D *hist, int res, std::string opt) { int nbins = hist->GetNbinsX(); double elow = hist->GetXaxis()->GetBinLowEdge(1); double ehigh = hist->GetXaxis()->GetBinLowEdge(nbins + 1); bool width = true; // opt.find("width") != std::string::npos; TH1D *fine = new TH1D("fine", "fine", nbins * res, elow, ehigh); TGraph *temp = new TGraph(); for (int i = 0; i < nbins; i++) { double E = hist->GetXaxis()->GetBinCenter(i + 1); double C = hist->GetBinContent(i + 1); double W = hist->GetXaxis()->GetBinWidth(i + 1); if (!width) W = 1.0; if (W != 0.0) temp->SetPoint(temp->GetN(), E, C / W); } for (int i = 0; i < fine->GetNbinsX(); i++) { double E = fine->GetXaxis()->GetBinCenter(i + 1); double W = fine->GetBinWidth(i + 1); if (!width) W = 1.0; fine->SetBinContent(i + 1, temp->Eval(E, 0, "S") * W); } fine->Scale(hist->Integral(1, hist->GetNbinsX() + 1) / fine->Integral(1, fine->GetNbinsX() + 1)); // std::cout << "Interpolation Difference = " //<< fine->Integral(1, fine->GetNbinsX() + 1) << "/" //<< hist->Integral(1, hist->GetNbinsX() + 1) << std::endl; return fine; } //******************************************************************** // This interpolates the flux by a TGraph instead of requiring the flux and MC // flux to have the same binning void PlotUtils::FluxUnfoldedScaling(TH1D *mcHist, TH1D *fhist, TH1D *ehist, double scalefactor, int nevents) { //******************************************************************** TH1D *eventhist = (TH1D *)ehist->Clone(); TH1D *fFluxHist = (TH1D *)fhist->Clone(); if (FitPar::Config().GetParB("save_flux_debug")) { std::string name = std::string(mcHist->GetName()); mcHist->Write((name + "_UNF_MC").c_str()); fFluxHist->Write((name + "_UNF_FLUX").c_str()); eventhist->Write((name + "_UNF_EVT").c_str()); TH1D *scalehist = new TH1D("scalehist", "scalehist", 1, 0.0, 1.0); scalehist->SetBinContent(1, scalefactor); scalehist->SetBinContent(2, nevents); scalehist->Write((name + "_UNF_SCALE").c_str()); } // Undo width integral in SF mcHist->Scale(scalefactor / eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width")); // Standardise The Flux eventhist->Scale(1.0 / fFluxHist->Integral()); fFluxHist->Scale(1.0 / fFluxHist->Integral()); // Scale mcHist by eventhist integral mcHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1)); // Now Get a flux PDF TH1D *pdfflux = (TH1D *)mcHist->Clone(); pdfflux->Reset(); for (int i = 0; i < mcHist->GetNbinsX(); i++) { double Ml = mcHist->GetXaxis()->GetBinLowEdge(i + 1); double Mh = mcHist->GetXaxis()->GetBinLowEdge(i + 2); // double Mc = mcHist->GetXaxis()->GetBinCenter(i+1); // double Me = mcHist->GetBinContent(i+1); // double Mw = mcHist->GetBinWidth(i+1); double fluxint = 0.0; for (int j = 0; j < fFluxHist->GetNbinsX(); j++) { // double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1); double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1); double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2); double Fe = fFluxHist->GetBinContent(j + 1); double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1); if (Fl >= Ml and Fh <= Mh) { fluxint += Fe; } else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) { fluxint += Fe * (Fh - Ml) / Fw; } else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) { fluxint += Fe * (Mh - Fl) / Fw; } else if (Ml >= Fl and Mh <= Fh) { fluxint += Fe * (Mh - Ml) / Fw; } else { continue; } } pdfflux->SetBinContent(i + 1, fluxint); } // Scale MC hist by pdfflux for (int i = 0; i < mcHist->GetNbinsX(); i++) { if (pdfflux->GetBinContent(i + 1) == 0.0) continue; mcHist->SetBinContent(i + 1, mcHist->GetBinContent(i + 1) / pdfflux->GetBinContent(i + 1)); mcHist->SetBinError(i + 1, mcHist->GetBinError(i + 1) / pdfflux->GetBinContent(i + 1)); } delete eventhist; delete fFluxHist; }; // MOVE TO GENERAL UTILS //******************************************************************** void PlotUtils::Set2DHistFromText(std::string dataFile, TH2 *hist, double norm, bool skipbins) { //******************************************************************** std::string line; std::ifstream data(dataFile.c_str(), std::ifstream::in); int yBin = 0; while (std::getline(data >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToDbl(line, " "); // Loop over entries and insert them into the histogram for (uint xBin = 0; xBin < entries.size(); xBin++) { if (!skipbins || entries[xBin] != -1.0) hist->SetBinContent(xBin + 1, yBin + 1, entries[xBin] * norm); } yBin++; } return; } // MOVE TO GENERAL UTILS TH1D *PlotUtils::GetTH1DFromFile(std::string dataFile, std::string title, std::string fPlotTitles, std::string alt_name) { TH1D *tempPlot; // If format is a root file if (dataFile.find(".root") != std::string::npos) { TFile *temp_infile = new TFile(dataFile.c_str(), "READ"); tempPlot = (TH1D *)temp_infile->Get(title.c_str()); tempPlot->SetDirectory(0); temp_infile->Close(); delete temp_infile; // Else its a space seperated txt file } else { // Make a TGraph Errors TGraphErrors *gr = new TGraphErrors(dataFile.c_str(), "%lg %lg %lg"); if (gr->IsZombie()) { - NUIS_ABORT(dataFile - << " is a zombie and could not be read. Are you sure it exists?" - << std::endl); + NUIS_ABORT( + dataFile + << " is a zombie and could not be read. Are you sure it exists?" + << std::endl); } double *bins = gr->GetX(); double *values = gr->GetY(); double *errors = gr->GetEY(); int npoints = gr->GetN(); // Fill the histogram from it tempPlot = new TH1D(title.c_str(), title.c_str(), npoints - 1, bins); for (int i = 0; i < npoints; ++i) { tempPlot->SetBinContent(i + 1, values[i]); // If only two columns are present in the input file, use the sqrt(values) // as the error equivalent to assuming that the error is statistical. Also // check that we're looking at an event rate rather than a cross section if (!errors[i] && values[i] > 1E-30) { tempPlot->SetBinError(i + 1, sqrt(values[i])); } else { tempPlot->SetBinError(i + 1, errors[i]); } } delete gr; } // Allow alternate naming for root files if (!alt_name.empty()) { tempPlot->SetNameTitle(alt_name.c_str(), alt_name.c_str()); } // Allow alternate axis titles if (!fPlotTitles.empty()) { tempPlot->SetNameTitle( tempPlot->GetName(), (std::string(tempPlot->GetTitle()) + fPlotTitles).c_str()); } return tempPlot; }; TH1D *PlotUtils::GetRatioPlot(TH1D *hist1, TH1D *hist2) { // make copy of first hist TH1D *new_hist = (TH1D *)hist1->Clone(); // Do bins and errors ourselves as scales can go awkward for (int i = 0; i < new_hist->GetNbinsX(); i++) { if (hist2->GetBinContent(i + 1) == 0.0) { new_hist->SetBinContent(i + 1, 0.0); } new_hist->SetBinContent(i + 1, hist1->GetBinContent(i + 1) / hist2->GetBinContent(i + 1)); new_hist->SetBinError(i + 1, hist1->GetBinError(i + 1) / hist2->GetBinContent(i + 1)); } return new_hist; }; TH1D *PlotUtils::GetRenormalisedPlot(TH1D *hist1, TH1D *hist2) { // make copy of first hist TH1D *new_hist = (TH1D *)hist1->Clone(); if (hist1->Integral("width") == 0 or hist2->Integral("width") == 0) { new_hist->Reset(); return new_hist; } Double_t scaleF = hist2->Integral("width") / hist1->Integral("width"); new_hist->Scale(scaleF); return new_hist; }; TH1D *PlotUtils::GetShapePlot(TH1D *hist1) { // make copy of first hist TH1D *new_hist = (TH1D *)hist1->Clone(); if (hist1->Integral("width") == 0) { new_hist->Reset(); return new_hist; } Double_t scaleF1 = 1.0 / hist1->Integral("width"); new_hist->Scale(scaleF1); return new_hist; }; TH1D *PlotUtils::GetShapeRatio(TH1D *hist1, TH1D *hist2) { TH1D *new_hist1 = GetShapePlot(hist1); TH1D *new_hist2 = GetShapePlot(hist2); // Do bins and errors ourselves as scales can go awkward for (int i = 0; i < new_hist1->GetNbinsX(); i++) { if (hist2->GetBinContent(i + 1) == 0) { new_hist1->SetBinContent(i + 1, 0.0); } new_hist1->SetBinContent(i + 1, new_hist1->GetBinContent(i + 1) / new_hist2->GetBinContent(i + 1)); new_hist1->SetBinError(i + 1, new_hist1->GetBinError(i + 1) / new_hist2->GetBinContent(i + 1)); } delete new_hist2; return new_hist1; }; TH2D *PlotUtils::GetCovarPlot(TMatrixDSym *cov, std::string name, std::string title) { TH2D *CovarPlot; if (cov) CovarPlot = new TH2D((*cov)); else CovarPlot = new TH2D(name.c_str(), title.c_str(), 1, 0, 1, 1, 0, 1); CovarPlot->SetName(name.c_str()); CovarPlot->SetTitle(title.c_str()); return CovarPlot; } TH2D *PlotUtils::GetFullCovarPlot(TMatrixDSym *cov, std::string name) { return PlotUtils::GetCovarPlot( cov, name + "_COV", name + "_COV;Bins;Bins;Covariance (#times10^{-76})"); } TH2D *PlotUtils::GetInvCovarPlot(TMatrixDSym *cov, std::string name) { return PlotUtils::GetCovarPlot( cov, name + "_INVCOV", name + "_INVCOV;Bins;Bins;Inv. Covariance (#times10^{-76})"); } TH2D *PlotUtils::GetDecompCovarPlot(TMatrixDSym *cov, std::string name) { return PlotUtils::GetCovarPlot( cov, name + "_DECCOV", name + "_DECCOV;Bins;Bins;Decomp Covariance (#times10^{-76})"); } TH1D *PlotUtils::GetTH1DFromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TFile *rootHistFile = new TFile(file.c_str(), "READ"); TH1D *tempHist = (TH1D *)rootHistFile->Get(name.c_str())->Clone(); if (tempHist == NULL) { NUIS_ABORT("Could not find distribution " << name << " in file " << file); } tempHist->SetDirectory(0); rootHistFile->Close(); return tempHist; } TH2D *PlotUtils::GetTH2DFromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TFile *rootHistFile = new TFile(file.c_str(), "READ"); TH2D *tempHist = (TH2D *)rootHistFile->Get(name.c_str())->Clone(); tempHist->SetDirectory(0); rootHistFile->Close(); delete rootHistFile; return tempHist; } TH1 *PlotUtils::GetTH1FromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TFile *rootHistFile = new TFile(file.c_str(), "READ"); if (!rootHistFile || rootHistFile->IsZombie()) { NUIS_ABORT("Couldn't open root file: \"" << file << "\"."); } TH1 *tempHist = dynamic_cast(rootHistFile->Get(name.c_str())->Clone()); if (!tempHist) { - NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \"" << file - << "\"."); + NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \"" + << file << "\"."); } tempHist->SetDirectory(0); rootHistFile->Close(); delete rootHistFile; return tempHist; } TGraph *PlotUtils::GetTGraphFromRootFile(std::string file, std::string name) { if (name.empty()) { std::vector tempfile = GeneralUtils::ParseToStr(file, ";"); file = tempfile[0]; name = tempfile[1]; } TDirectory *olddir = gDirectory; TFile *rootHistFile = new TFile(file.c_str(), "READ"); if (!rootHistFile || rootHistFile->IsZombie()) { NUIS_ABORT("Couldn't open root file: \"" << file << "\"."); } TDirectory *newdir = gDirectory; TGraph *temp = dynamic_cast(rootHistFile->Get(name.c_str())->Clone()); if (!temp) { - NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \"" << file - << "\"."); + NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \"" + << file << "\"."); } newdir->Remove(temp); olddir->Append(temp); rootHistFile->Close(); olddir->cd(); return temp; } /// Returns a vector of named TH1*s found in a single input file. /// /// Expects a descriptor like: file.root[hist1|hist2|...] std::vector PlotUtils::GetTH1sFromRootFile(std::string const &descriptor) { std::vector descriptors = GeneralUtils::ParseToStr(descriptor, ","); std::vector hists; for (size_t d_it = 0; d_it < descriptors.size(); ++d_it) { std::string &d = descriptors[d_it]; std::vector fname = GeneralUtils::ParseToStr(d, "["); if (!fname.size() || !fname[0].length()) { NUIS_ABORT("Couldn't find input file when attempting to parse : \"" - << d << "\". Expected input.root[hist1|hist2|...]."); + << d << "\". Expected input.root[hist1|hist2|...]."); } if (fname[1][fname[1].length() - 1] == ']') { fname[1] = fname[1].substr(0, fname[1].length() - 1); } std::vector histnames = GeneralUtils::ParseToStr(fname[1], "|"); if (!histnames.size()) { - NUIS_ABORT("Couldn't find any histogram name specifiers when attempting to " - "parse " - ": \"" - << fname[1] << "\". Expected hist1|hist2|..."); + NUIS_ABORT( + "Couldn't find any histogram name specifiers when attempting to " + "parse " + ": \"" + << fname[1] << "\". Expected hist1|hist2|..."); } TFile *rootHistFile = new TFile(fname[0].c_str(), "READ"); if (!rootHistFile || rootHistFile->IsZombie()) { NUIS_ABORT("Couldn't open root file: \"" << fname[0] << "\"."); } for (size_t i = 0; i < histnames.size(); ++i) { TH1 *tempHist = dynamic_cast(rootHistFile->Get(histnames[i].c_str())->Clone()); if (!tempHist) { NUIS_ABORT("Couldn't retrieve: \"" - << histnames[i] << "\" from root file: \"" << fname[0] << "\"."); + << histnames[i] << "\" from root file: \"" << fname[0] + << "\"."); } tempHist->SetDirectory(0); hists.push_back(tempHist); } rootHistFile->Close(); } return hists; } // Create an array from an input file std::vector PlotUtils::GetArrayFromTextFile(std::string DataFile) { std::string line; std::ifstream data(DataFile.c_str(), std::ifstream::in); // Get first line std::getline(data >> std::ws, line, '\n'); // Convert from a string into a vector of double std::vector entries = GeneralUtils::ParseToDbl(line, " "); return entries; } // Get a 2D array from a text file -std::vector > +std::vector> PlotUtils::Get2DArrayFromTextFile(std::string DataFile) { std::string line; - std::vector > DataArray; + std::vector> DataArray; std::ifstream data(DataFile.c_str(), std::ifstream::in); while (std::getline(data >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToDbl(line, " "); DataArray.push_back(entries); } return DataArray; } TH2D *PlotUtils::GetTH2DFromTextFile(std::string data, std::string binx, std::string biny) { // First read in the binning // Array of x binning std::vector xbins = GetArrayFromTextFile(binx); // Array of y binning std::vector ybins = GetArrayFromTextFile(biny); // Read in the data - std::vector > Data = Get2DArrayFromTextFile(data); + std::vector> Data = Get2DArrayFromTextFile(data); // And finally fill the data TH2D *DataPlot = new TH2D("TempHist", "TempHist", xbins.size() - 1, &xbins[0], ybins.size() - 1, &ybins[0]); int nBinsX = 0; int nBinsY = 0; - for (std::vector >::iterator it = Data.begin(); + for (std::vector>::iterator it = Data.begin(); it != Data.end(); ++it) { nBinsX++; // Get the inner vector std::vector temp = *it; // Save the previous number[of bins to make sure it's uniform binning int oldBinsY = nBinsY; // Reset the counter nBinsY = 0; for (std::vector::iterator jt = temp.begin(); jt != temp.end(); ++jt) { nBinsY++; DataPlot->SetBinContent(nBinsX, nBinsY, *jt); DataPlot->SetBinError(nBinsX, nBinsY, 0.0); } if (oldBinsY > 0 && oldBinsY != nBinsY) { NUIS_ERR(FTL, "Found non-uniform y-binning in " << data); NUIS_ERR(FTL, "Previous slice: " << oldBinsY); NUIS_ERR(FTL, "Current slice: " << nBinsY); NUIS_ABORT("Non-uniform binning is not supported in " - "PlotUtils::GetTH2DFromTextFile"); + "PlotUtils::GetTH2DFromTextFile"); } } // Check x bins if (size_t(nBinsX + 1) != xbins.size()) { - NUIS_ERR(FTL, "Number of x bins in data histogram does not match the binning " - "histogram!"); NUIS_ERR(FTL, - "Are they the wrong way around (i.e. xbinning should be ybinning)?"); + "Number of x bins in data histogram does not match the binning " + "histogram!"); + NUIS_ERR( + FTL, + "Are they the wrong way around (i.e. xbinning should be ybinning)?"); NUIS_ERR(FTL, "Data: " << nBinsX); NUIS_ABORT("From " << binx << " binning: " << xbins.size()); } // Check y bins if (size_t(nBinsY + 1) != ybins.size()) { - NUIS_ERR(FTL, "Number of y bins in data histogram does not match the binning " - "histogram!"); NUIS_ERR(FTL, - "Are they the wrong way around (i.e. xbinning should be ybinning)?"); + "Number of y bins in data histogram does not match the binning " + "histogram!"); + NUIS_ERR( + FTL, + "Are they the wrong way around (i.e. xbinning should be ybinning)?"); NUIS_ERR(FTL, "Data: " << nBinsY); NUIS_ABORT("From " << biny << " binning: " << ybins.size()); } return DataPlot; } TH1D *PlotUtils::GetSliceY(TH2D *Hist, int SliceNo) { TH1D *Slice = Hist->ProjectionX(Form("%s_SLICEY%i", Hist->GetName(), SliceNo), SliceNo, SliceNo, "e"); Slice->SetTitle(Form("%s, %.2f-%.2f", Hist->GetYaxis()->GetTitle(), Hist->GetYaxis()->GetBinLowEdge(SliceNo), Hist->GetYaxis()->GetBinLowEdge(SliceNo + 1))); Slice->GetYaxis()->SetTitle(Hist->GetZaxis()->GetTitle()); return Slice; } TH1D *PlotUtils::GetSliceX(TH2D *Hist, int SliceNo) { TH1D *Slice = Hist->ProjectionY(Form("%s_SLICEX%i", Hist->GetName(), SliceNo), SliceNo, SliceNo, "e"); Slice->SetTitle(Form("%s, %.2f-%.2f", Hist->GetXaxis()->GetTitle(), Hist->GetXaxis()->GetBinLowEdge(SliceNo), Hist->GetXaxis()->GetBinLowEdge(SliceNo + 1))); Slice->GetYaxis()->SetTitle(Hist->GetZaxis()->GetTitle()); return Slice; } void PlotUtils::AddNeutModeArray(TH1D *hist1[], TH1D *hist2[], double scaling) { for (int i = 0; i < 60; i++) { if (!hist2[i]) continue; if (!hist1[i]) continue; hist1[i]->Add(hist2[i], scaling); } return; } void PlotUtils::ScaleToData(TH1D *data, TH1D *mc, TH1I *mask) { double scaleF = GetDataMCRatio(data, mc, mask); mc->Scale(scaleF); return; } void PlotUtils::MaskBins(TH1D *hist, TH1I *mask) { for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1) <= 0.5) continue; hist->SetBinContent(i + 1, 0.0); hist->SetBinError(i + 1, 0.0); - NUIS_LOG(REC, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 - << " to 0.0 +- 0.0"); + NUIS_LOG(DEB, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 + << " to 0.0 +- 0.0"); } return; } void PlotUtils::MaskBins(TH2D *hist, TH2I *mask) { for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { if (mask->GetBinContent(i + 1, j + 1) <= 0.5) continue; hist->SetBinContent(i + 1, j + 1, 0.0); hist->SetBinError(i + 1, j + 1, 0.0); - NUIS_LOG(REC, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " " - << j + 1 << " to 0.0 +- 0.0"); + NUIS_LOG(DEB, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 + << " " << j + 1 << " to 0.0 +- 0.0"); } } return; } double PlotUtils::GetDataMCRatio(TH1D *data, TH1D *mc, TH1I *mask) { double rat = 1.0; TH1D *newmc = (TH1D *)mc->Clone(); TH1D *newdt = (TH1D *)data->Clone(); if (mask) { MaskBins(newmc, mask); MaskBins(newdt, mask); } rat = newdt->Integral() / newmc->Integral(); return rat; } TH2D *PlotUtils::GetCorrelationPlot(TH2D *cov, std::string name) { TH2D *cor = (TH2D *)cov->Clone(); cor->Reset(); for (int i = 0; i < cov->GetNbinsX(); i++) { for (int j = 0; j < cov->GetNbinsY(); j++) { if (cov->GetBinContent(i + 1, i + 1) != 0.0 and cov->GetBinContent(j + 1, j + 1) != 0.0) cor->SetBinContent(i + 1, j + 1, cov->GetBinContent(i + 1, j + 1) / (sqrt(cov->GetBinContent(i + 1, i + 1) * cov->GetBinContent(j + 1, j + 1)))); } } if (!name.empty()) { cor->SetNameTitle(name.c_str(), (name + ";;correlation").c_str()); } cor->SetMinimum(-1); cor->SetMaximum(1); return cor; } TH2D *PlotUtils::GetDecompPlot(TH2D *cov, std::string name) { TMatrixDSym *covarmat = new TMatrixDSym(cov->GetNbinsX()); for (int i = 0; i < cov->GetNbinsX(); i++) for (int j = 0; j < cov->GetNbinsY(); j++) (*covarmat)(i, j) = cov->GetBinContent(i + 1, j + 1); TMatrixDSym *decompmat = StatUtils::GetDecomp(covarmat); TH2D *dec = (TH2D *)cov->Clone(); for (int i = 0; i < cov->GetNbinsX(); i++) for (int j = 0; j < cov->GetNbinsY(); j++) dec->SetBinContent(i + 1, j + 1, (*decompmat)(i, j)); delete covarmat; delete decompmat; dec->SetNameTitle(name.c_str(), (name + ";;;decomposition").c_str()); return dec; } TH2D *PlotUtils::MergeIntoTH2D(TH1D *xhist, TH1D *yhist, std::string zname) { std::vector xedges, yedges; for (int i = 0; i < xhist->GetNbinsX() + 2; i++) { xedges.push_back(xhist->GetXaxis()->GetBinLowEdge(i + 1)); } for (int i = 0; i < yhist->GetNbinsX() + 2; i++) { yedges.push_back(yhist->GetXaxis()->GetBinLowEdge(i + 1)); } int nbinsx = xhist->GetNbinsX(); int nbinsy = yhist->GetNbinsX(); std::string name = std::string(xhist->GetName()) + "_vs_" + std::string(yhist->GetName()); std::string titles = ";" + std::string(xhist->GetXaxis()->GetTitle()) + ";" + std::string(yhist->GetXaxis()->GetTitle()) + ";" + zname; TH2D *newplot = new TH2D(name.c_str(), (name + titles).c_str(), nbinsx, &xedges[0], nbinsy, &yedges[0]); return newplot; } //*************************************************** void PlotUtils::MatchEmptyBins(TH1D *data, TH1D *mc) { //************************************************** for (int i = 0; i < data->GetNbinsX(); i++) { if (data->GetBinContent(i + 1) == 0.0 or data->GetBinError(i + 1) == 0.0) mc->SetBinContent(i + 1, 0.0); } return; } //*************************************************** void PlotUtils::MatchEmptyBins(TH2D *data, TH2D *mc) { //************************************************** 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 or data->GetBinError(i + 1, j + 1) == 0.0) mc->SetBinContent(i + 1, j + 1, 0.0); } } return; } //*************************************************** TH1D *PlotUtils::GetProjectionX(TH2D *hist, TH2I *mask) { //*************************************************** TH2D *maskedhist = StatUtils::ApplyHistogramMasking(hist, mask); // This includes the underflow/overflow - TH1D* hist_X = maskedhist->ProjectionX("_px", 1, maskedhist->GetXaxis()->GetNbins()); + TH1D *hist_X = + maskedhist->ProjectionX("_px", 1, maskedhist->GetXaxis()->GetNbins()); hist_X->SetTitle(Form("%s x no under/overflow", hist_X->GetTitle())); delete maskedhist; return hist_X; } //*************************************************** TH1D *PlotUtils::GetProjectionY(TH2D *hist, TH2I *mask) { //*************************************************** TH2D *maskedhist = StatUtils::ApplyHistogramMasking(hist, mask); // This includes the underflow/overflow - TH1D* hist_Y = maskedhist->ProjectionY("_py", 1, maskedhist->GetYaxis()->GetNbins()); + TH1D *hist_Y = + maskedhist->ProjectionY("_py", 1, maskedhist->GetYaxis()->GetNbins()); hist_Y->SetTitle(Form("%s y no under/overflow", hist_Y->GetTitle())); delete maskedhist; return hist_Y; }