diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx index 16a8602..d962403 100644 --- a/src/FitBase/Measurement1D.cxx +++ b/src/FitBase/Measurement1D.cxx @@ -1,1901 +1,1903 @@ // 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; // 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; // 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; } //******************************************************************** 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 LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl; if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) { fIsRawEvents = true; LOG(SAM) << "Found event rate measurement but using poisson likelihoods." << std::endl; } if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; LOG(SAM) << "::" << fName << "::" << std::endl; LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, " << "not flux averaged!" << std::endl; } if (fIsEnu1D && fIsRawEvents) { LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!" << std::endl; LOG(SAM) << "Check experiment constructor for " << fName << " and correct this!" << std::endl; LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl; exit(-1); } 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) { ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl; ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl; 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) { //******************************************************************** LOG(SAM) << "Reading data from text file: " << datafile << std::endl; fDataHist = PlotUtils::GetTH1DFromFile(datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles()); } //******************************************************************** void Measurement1D::SetDataFromRootFile(std::string datafile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl; 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) { ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl; ERR(FTL) << "Setup Data First!" << std::endl; 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) { LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl; fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } else { ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl; } // if (!fIsDiag) { // ERR(FTL) << "SetCovarMatrixFromDiag called for measurement " // << "that is not set as diagonal." << std::endl; // throw; // } } //******************************************************************** void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl; 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){ LOG(SAM) << "Reading covariance from text file: " << covList[i] << std::endl; 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) { //******************************************************************** LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl; 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(); } LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl; covar = StatUtils::GetCovarFromTextFile(covfile, dim); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl; 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(); LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl; TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim); if (!fDataHist) { ERR(FTL) << "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" << std::endl; throw; } // 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){ LOG(SAM) << "Reading covariance from text file: " << corrList[i] << std::endl; 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) { //******************************************************************** LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl; TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname); if (!fDataHist) { ERR(FTL) << "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" << std::endl; throw; } // 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(); } LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl; 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) { //******************************************************************** LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl; 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; LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl; // 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()) { LOG(FTL) << " Cannot find mask file." << std::endl; 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) { 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. 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() { //******************************************************************** LOG(SAM) << "Finalising Measurement: " << fName << std::endl; 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) { ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl; throw; } // Make sure covariances are setup if (!fFullCovar) { fIsDiag = true; SetCovarFromDiagonal(fDataHist); } 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 and 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) { 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) { ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl; ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl; ERR(FTL) << "EXITING" << std::endl; 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) { ERR(FTL) << "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." << std::endl; throw; } } } // 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) { ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i) << "' Provided is not allowed for this measurement." << std::endl; ERR(FTL) << "Fit Options should be provided as a '/' seperated list " "(e.g. FREE/DIAG/NORM)" << std::endl; ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes << "'" << std::endl; throw; } } // 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; ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl; ERR(FTL) << "Try to use a chi2!" << std::endl; 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()) LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl; else ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile << std::endl; 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) { ERR(WRN) << fName << ": attempted to apply smearing matrix, but none was set" << std::endl; 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) { QLOG(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) { - // PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, 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); } } // 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; LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl; // 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()); } if (drawOpt.find("INVCOV") != std::string::npos && covar) { PlotUtils::GetInvCovarPlot(covar, fSettings.GetName()); } if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) { PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName()); } // // 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; // } // Write Extra Histograms AutoWriteExtraTH1(); WriteExtraHistograms(); // Returning LOG(SAM) << "Written Histograms: " << fName << std::endl; 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; LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!" << std::endl; LOG(SAM) << "Overriding this and setting fIsRawEvents == true!" << std::endl; } fIsEnu1D = false; if (fName.find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; LOG(SAM) << "::" << fName << "::" << std::endl; LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, " "not flux averaged!" << std::endl; } if (fIsEnu1D && fIsRawEvents) { LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!" << std::endl; LOG(SAM) << "Check experiment constructor for " << fName << " and correct this!" << std::endl; LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl; - exit(-1); + 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 LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl; 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) { //******************************************************************** LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname << std::endl; 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) { //******************************************************************** LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname << std::endl; 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* decmpPlot = new TH2D(); TH2D* covarInvPlot = 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()); covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").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 ERR(WRN) << "Incorrect thrown_covariance option in parameters." << std::endl; 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()) LOG(SAM) << "Reading covariance matrix from file: " << covarFile << std::endl; else ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile << std::endl; // 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) { ERR(WRN) << "SetCovarMatrixFromText -> Covariance matrix only has <= 1 " "entries on this line: " << row << std::endl; } 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()) LOG(SAM) << "Reading and converting correlation matrix from file: " << corrFile << std::endl; else { ERR(FTL) << "Correlation matrix provided is incorrect: " << corrFile << std::endl; exit(-1); } 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) { ERR(FTL) << "Found a zero value in the covariance matrix, assuming " "this is an error!" << std::endl; exit(-1); } (*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) { //******************************************************************** LOG(SAM) << "Getting covariance from " << covarFile << "->" << covName << std::endl; 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/StackBase.cxx b/src/FitBase/StackBase.cxx index 3abde11..d1a66e9 100644 --- a/src/FitBase/StackBase.cxx +++ b/src/FitBase/StackBase.cxx @@ -1,219 +1,216 @@ #include "StackBase.h" void StackBase::AddMode(std::string name, std::string title, - int linecolor, int linewidth, - int fillstyle) { + int linecolor, int linewidth, + int fillstyle) { - // int ncur = fAllLabels.size(); - fAllLabels.push_back(name); - fAllTitles.push_back(title); + // int ncur = fAllLabels.size(); + fAllLabels.push_back(name); + fAllTitles.push_back(title); - std::vector temp; - temp.push_back(linecolor); - temp.push_back(linewidth); - temp.push_back(fillstyle); + std::vector temp; + temp.push_back(linecolor); + temp.push_back(linewidth); + temp.push_back(fillstyle); - fAllStyles.push_back(temp); + fAllStyles.push_back(temp); } -void StackBase::FluxUnfold(TH1D* flux, TH1D* events, double scalefactor) { +void StackBase::FluxUnfold(TH1D* flux, TH1D* events, double scalefactor, int nevents) { + for (size_t i = 0; i < fAllLabels.size(); i++) { + if (fNDim == 1) { + PlotUtils::FluxUnfoldedScaling((TH1D*) fAllHists[i], flux, events, scalefactor, nevents); + } else if (fNDim == 2) { + PlotUtils::FluxUnfoldedScaling((TH2D*) fAllHists[i], flux, events, scalefactor); - - - for (size_t i = 0; i < fAllLabels.size(); i++) { - if (fNDim == 1) { - PlotUtils::FluxUnfoldedScaling((TH1D*) fAllHists[i], flux, events, scalefactor); - } else if (fNDim == 2) { - PlotUtils::FluxUnfoldedScaling((TH2D*) fAllHists[i], flux, events, scalefactor); - - } - } + } + } } void StackBase::AddMode(int index, std::string name, std::string title, - int linecolor, int linewidth, int fillstyle) { + int linecolor, int linewidth, int fillstyle) { - while (fAllLabels.size() <= (UInt_t)index) { - fAllLabels.push_back(""); - fAllTitles.push_back(""); - fAllStyles.push_back(std::vector(1, 1)); - } + while (fAllLabels.size() <= (UInt_t)index) { + fAllLabels.push_back(""); + fAllTitles.push_back(""); + fAllStyles.push_back(std::vector(1, 1)); + } - fAllLabels[index] = (name); - fAllTitles[index] = (title); + fAllLabels[index] = (name); + fAllTitles[index] = (title); - std::vector temp; - temp.push_back(linecolor); - temp.push_back(linewidth); - temp.push_back(fillstyle); + std::vector temp; + temp.push_back(linecolor); + temp.push_back(linewidth); + temp.push_back(fillstyle); - fAllStyles[index] = temp; + fAllStyles[index] = temp; } bool StackBase::IncludeInStack(TH1* hist) { - if (!FitPar::Config().GetParB("includeemptystackhists") and - hist->Integral() == 0.0) return false; - return true; + if (!FitPar::Config().GetParB("includeemptystackhists") and + hist->Integral() == 0.0) return false; + return true; } bool StackBase::IncludeInStack(int index) { - return true; + return true; } void StackBase::SetupStack(TH1* hist) { - fTemplate = (TH1*) hist->Clone(fName.c_str()); - fTemplate->Reset(); + fTemplate = (TH1*) hist->Clone(fName.c_str()); + fTemplate->Reset(); - // Determine template dim - fNDim = fTemplate->GetDimension(); + // Determine template dim + fNDim = fTemplate->GetDimension(); - for (size_t i = 0; i < fAllLabels.size(); i++) { - fAllHists.push_back( (TH1*) fTemplate->Clone( (fName + "_" + fAllLabels[i]).c_str() ) ); - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + fAllHists.push_back( (TH1*) fTemplate->Clone( (fName + "_" + fAllLabels[i]).c_str() ) ); + } }; void StackBase::Scale(double sf, std::string opt) { - for (size_t i = 0; i < fAllLabels.size(); i++) { - // std::cout << "Scaling Stack Hist " << i << " by " << sf << std::endl; - fAllHists[i]->Scale(sf, opt.c_str()); - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + // std::cout << "Scaling Stack Hist " << i << " by " << sf << std::endl; + fAllHists[i]->Scale(sf, opt.c_str()); + } }; void StackBase::Reset() { - for (size_t i = 0; i < fAllLabels.size(); i++) { - fAllHists[i]->Reset(); - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + fAllHists[i]->Reset(); + } }; void StackBase::FillStack(int index, double x, double y, double z, double weight) { - if (index < 0 or (UInt_t)index >= fAllLabels.size()) { - ERR(WRN) << "Returning Stack Fill Because Range = " << index << " " << fAllLabels.size() << std::endl; - return; - } + if (index < 0 or (UInt_t)index >= fAllLabels.size()) { + ERR(WRN) << "Returning Stack Fill Because Range = " << index << " " << fAllLabels.size() << std::endl; + return; + } - if (fNDim == 1) fAllHists[index]->Fill(x, y); - else if (fNDim == 2) { - // std::cout << "Filling 2D Stack " << index << " " << x << " " << y << " " << z << std::endl; - ((TH2*)fAllHists[index])->Fill(x, y, z); - } + if (fNDim == 1) fAllHists[index]->Fill(x, y); + else if (fNDim == 2) { + // std::cout << "Filling 2D Stack " << index << " " << x << " " << y << " " << z << std::endl; + ((TH2*)fAllHists[index])->Fill(x, y, z); + } - else if (fNDim == 3) ((TH3*)fAllHists[index])->Fill(x, y, z, weight); + else if (fNDim == 3) ((TH3*)fAllHists[index])->Fill(x, y, z, weight); } void StackBase::Write() { - THStack* st = new THStack(); - - // Loop and add all histograms - bool saveseperate = FitPar::Config().GetParB("WriteSeperateStacks"); - for (size_t i = 0; i < fAllLabels.size(); i++) { - - if (!IncludeInStack(fAllHists[i])) continue; - if (!IncludeInStack(i)) continue; - - fAllHists[i]->SetTitle( fAllTitles[i].c_str() ); - fAllHists[i]->GetXaxis()->SetTitle( fXTitle.c_str() ); - fAllHists[i]->GetYaxis()->SetTitle( fYTitle.c_str() ); - fAllHists[i]->GetZaxis()->SetTitle( fZTitle.c_str() ); - fAllHists[i]->SetLineColor( fAllStyles[i][0] ); - fAllHists[i]->SetLineWidth( fAllStyles[i][1] ); - fAllHists[i]->SetFillStyle( fAllStyles[i][2] ); - fAllHists[i]->SetFillColor( fAllStyles[i][0] ); - if (saveseperate) fAllHists[i]->Write(); - - st->Add(fAllHists[i]); - } - st->SetTitle(fTitle.c_str()); - st->SetName(fName.c_str()); - st->Write(); - delete st; + THStack* st = new THStack(); + + // Loop and add all histograms + bool saveseperate = FitPar::Config().GetParB("WriteSeperateStacks"); + for (size_t i = 0; i < fAllLabels.size(); i++) { + + if (!IncludeInStack(fAllHists[i])) continue; + if (!IncludeInStack(i)) continue; + + fAllHists[i]->SetTitle( fAllTitles[i].c_str() ); + fAllHists[i]->GetXaxis()->SetTitle( fXTitle.c_str() ); + fAllHists[i]->GetYaxis()->SetTitle( fYTitle.c_str() ); + fAllHists[i]->GetZaxis()->SetTitle( fZTitle.c_str() ); + fAllHists[i]->SetLineColor( fAllStyles[i][0] ); + fAllHists[i]->SetLineWidth( fAllStyles[i][1] ); + fAllHists[i]->SetFillStyle( fAllStyles[i][2] ); + fAllHists[i]->SetFillColor( fAllStyles[i][0] ); + if (saveseperate) fAllHists[i]->Write(); + + st->Add(fAllHists[i]); + } + st->SetTitle(fTitle.c_str()); + st->SetName(fName.c_str()); + st->Write(); + delete st; }; void StackBase::Multiply(TH1* hist) { - for (size_t i = 0; i < fAllLabels.size(); i++) { - fAllHists[i]->Multiply(hist); - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + fAllHists[i]->Multiply(hist); + } } void StackBase::Divide(TH1* hist) { - for (size_t i = 0; i < fAllLabels.size(); i++) { - fAllHists[i]->Divide(hist); - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + fAllHists[i]->Divide(hist); + } } void StackBase::Add(TH1* hist, double scale) { - for (size_t i = 0; i < fAllLabels.size(); i++) { - fAllHists[i]->Add(hist, scale); - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + fAllHists[i]->Add(hist, scale); + } } void StackBase::Add(StackBase* hist, double scale) { - if (hist->GetType() != fType) { - ERR(WRN) << "Trying to add two StackBases of different types!" << std::endl; - ERR(WRN) << fType << " + " << hist->GetType() << " = Undefined." << std::endl; - ERR(WRN) << "Doing nothing..." << std::endl; - return; - } + if (hist->GetType() != fType) { + ERR(WRN) << "Trying to add two StackBases of different types!" << std::endl; + ERR(WRN) << fType << " + " << hist->GetType() << " = Undefined." << std::endl; + ERR(WRN) << "Doing nothing..." << std::endl; + return; + } - for (size_t i = 0; i < fAllLabels.size(); i++) { - fAllHists[i]->Add(hist->GetHist(i)); - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + fAllHists[i]->Add(hist->GetHist(i)); + } } TH1* StackBase::GetHist(int entry) { - return fAllHists[entry]; + return fAllHists[entry]; } TH1* StackBase::GetHist(std::string label) { - TH1* hist = NULL; - std::vector splitlabels = GeneralUtils::ParseToStr(label, "+"); - for (size_t j = 0; j < splitlabels.size(); j++) { - std::string newlabel = splitlabels[j]; + TH1* hist = NULL; + std::vector splitlabels = GeneralUtils::ParseToStr(label, "+"); + for (size_t j = 0; j < splitlabels.size(); j++) { + std::string newlabel = splitlabels[j]; - for (size_t i = 0; i < fAllLabels.size(); i++) { - if (newlabel == fAllLabels[i]) { - if (!hist) hist = (TH1*) fAllHists[i]->Clone(); - else hist->Add(fAllHists[i]); - } - } - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + if (newlabel == fAllLabels[i]) { + if (!hist) hist = (TH1*) fAllHists[i]->Clone(); + else hist->Add(fAllHists[i]); + } + } + } - return hist; + return hist; } THStack StackBase::GetStack() { - THStack st = THStack(); - for (size_t i = 0; i < fAllLabels.size(); i++) { - st.Add(fAllHists[i]); - } - return st; + THStack st = THStack(); + for (size_t i = 0; i < fAllLabels.size(); i++) { + st.Add(fAllHists[i]); + } + return st; } void StackBase::AddNewHist(std::string name, TH1* hist) { - AddMode( fAllLabels.size(), name, hist->GetTitle(), hist->GetLineColor()); - fAllHists.push_back( (TH1*) hist->Clone() ); + AddMode( fAllLabels.size(), name, hist->GetTitle(), hist->GetLineColor()); + fAllHists.push_back( (TH1*) hist->Clone() ); } void StackBase::AddToCategory(std::string name, TH1* hist) { - for (size_t i = 0; i < fAllLabels.size(); i++) { - if (name == fAllLabels[i]) { - fAllHists[i]->Add(hist); - break; - } - } + for (size_t i = 0; i < fAllLabels.size(); i++) { + if (name == fAllLabels[i]) { + fAllHists[i]->Add(hist); + break; + } + } } void StackBase::AddToCategory(int index, TH1* hist) { - fAllHists[index]->Add(hist); + fAllHists[index]->Add(hist); } diff --git a/src/FitBase/StackBase.h b/src/FitBase/StackBase.h index 398afa6..c227813 100644 --- a/src/FitBase/StackBase.h +++ b/src/FitBase/StackBase.h @@ -1,126 +1,126 @@ #ifndef STACK_BASE_H #define STACK_BASE_H #include "MeasurementVariableBox.h" #include "FitLogger.h" #include "GeneralUtils.h" #include "TH1.h" #include "TH1D.h" #include "TH2D.h" #include "THStack.h" #include "TH2.h" #include "TH3.h" #include "PlotUtils.h" class StackBase { public: StackBase() {}; ~StackBase() {}; virtual void AddMode(std::string name, std::string title, int linecolor = 1, int linewidth = 1, int fillstyle = 1001); virtual void AddMode(int index, std::string name, std::string title, int linecolor = 1, int linewidth = 1, int fillstyle = 1001); virtual bool IncludeInStack(TH1* hist); virtual bool IncludeInStack(int index); virtual void SetupStack(TH1* hist); virtual void Scale(double sf, std::string opt = ""); - virtual void FluxUnfold(TH1D* flux, TH1D* events, double scalefactor); + virtual void FluxUnfold(TH1D* flux, TH1D* events, double scalefactor, int nevents); virtual void Reset(); virtual void FillStack(int index, double x, double y = 1.0, double z = 1.0, double weight = 1.0); virtual void Write(); virtual void Add(StackBase* hist, double scale); virtual void Add(TH1* hist, double scale); virtual void AddNewHist(std::string name, TH1* hist); virtual void AddToCategory(std::string name, TH1* hist); virtual void AddToCategory(int index, TH1* hist); virtual void Divide(TH1* hist); virtual void Multiply(TH1* hist); virtual TH1* GetHist(int entry); virtual TH1* GetHist(std::string label); virtual THStack GetStack(); std::string GetType(){return fType;}; std::string fName; std::string fTitle; std::string fXTitle; std::string fYTitle; std::string fZTitle; std::string fType; TH1* fTemplate; int fNDim; // Maps incase we want to be selective std::vector< std::vector > fAllStyles; std::vector fAllTitles; std::vector fAllLabels; std::vector fAllHists; }; /* class NuSpeciesStack : public StackBase { public: SetupStack(TH1* hist) { AddMode("numu", "numu", kBlue, 2, 3004); AddMode("numubar", "numubar", kRed, 2, 3004 ); AddMode("nue", "nue", kGreen, 2, 3004 ); StackBase::SetupStack(hist); }; void NuSpeciesStack::FillStack(int species, double x, double y = 1.0, double z = 1.0, double weight = 1.0) { Stackbase::FillStack(ConvertSpeciesToIndex(species), x, y, z, weight); } int ConvertSpeciesToIndex(int species) { switch (species) { case 14: return 0; case -14: return 1; case 11: return 2; default: return -1; } }; }; class TargetStack : public StackBase { public: SetupStack(TH1* hist) { AddMode("C", "C", kBlue, 2, 3004); AddMode("H", "H", kRed, 2, 3004 ); AddMode("O", "O", kGreen, 2, 3004 ); StackBase::SetupStack(hist); }; void NuSpeciesStack::FillStack(int species, double x, double y = 1.0, double z = 1.0, double weight = 1.0) { Stackbase::FillStack(ConvertTargetToIndex(species), x, y, z, weight); } int ConvertTargetToIndex(int target) { switch (species) { case 1000010010: return 0; case 1000: return 1; case 1000: return 2; default: return -1; } }; } */ #endif diff --git a/src/MCStudies/OfficialNIWGPlots.cxx b/src/MCStudies/OfficialNIWGPlots.cxx index 51349c8..3dc6239 100644 --- a/src/MCStudies/OfficialNIWGPlots.cxx +++ b/src/MCStudies/OfficialNIWGPlots.cxx @@ -1,679 +1,679 @@ // 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 "OfficialNIWGPlots.h" #include "T2K_SignalDef.h" #include "MINERvA_SignalDef.h" //******************************************************************** OfficialNIWGPlots::OfficialNIWGPlots(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "NIWG Official plots sample. \n" \ "Target: Any \n" \ "Flux: T2K Flux \n" \ "Signal: CC Inclusive \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 30.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); fSettings.SetTitle("NIWG Official Plots 2017"); fSettings.SetOnlyMC(1); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); fScaleFactorEnuXSec = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents); fScaleFactorDifXSec = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- fHist_NuMu_Enu = new TH1D("NuMu_Enu_MC", "NuMu_Enu_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_NuMu_Enu_Modes = new MCStudies::OfficialNIWGStack("NuMu_Enu_MC_MODES", "NuMu_Enu_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMu_Enu); fHist_NuMu_Enu_Pions = new MCStudies::OfficialPionStack("NuMu_Enu_MC_PIONS", "NuMu_Enu_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMu_Enu); fHist_NuMu_EnuRates = new TH1D("NuMu_EnuRates_MC", "NuMu_EnuRates_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_NuMu_EnuRates_Modes = new MCStudies::OfficialNIWGStack("NuMu_EnuRates_MC_MODES", "NuMu_EnuRates_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMu_Enu); fHist_NuMu_EnuRates_Pions = new MCStudies::OfficialPionStack("NuMu_EnuRates_MC_PIONS", "NuMu_EnuRates_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMu_Enu); fHist_NuMu_Q2 = new TH1D("NuMu_Q2_MC", "NuMu_Q2_MC;Q^{2} (GeV);d#sigma/dQ^{2} (cm^{2}/nucleon/GeV^{2})", 45, 0.0, 3.0); fHist_NuMu_Q2_Modes = new MCStudies::OfficialNIWGStack("NuMu_Q2_MC_MODES", "NuMu_Q2_MC_MODES;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_NuMu_Q2); fHist_NuMu_Q2_Pions = new MCStudies::OfficialPionStack("NuMu_Q2_MC_PIONS", "NuMu_Q2_MC_PIONS;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_NuMu_Q2); fHist_NuMu_Pmu = new TH1D("NuMu_Pmu_MC", "NuMu_Pmu_MC;P_{#mu} (GeV);#sigma (cm^{2}/nucleon/GeV)", 45, 0.0, 3.0); fHist_NuMu_Pmu_Modes = new MCStudies::OfficialNIWGStack("NuMu_Pmu_MC_MODES", "NuMu_Pmu_MC_MODES;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_NuMu_Pmu); fHist_NuMu_Pmu_Pions = new MCStudies::OfficialPionStack("NuMu_Pmu_MC_PIONS", "NuMu_Pmu_MC_PIONS;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_NuMu_Pmu); fHist_NuMu_Cosmu = new TH1D("NuMu_Cosmu_MC", "NuMu_Cosmu_MC;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", 45, -1.0, 1.0); fHist_NuMu_Cosmu_Modes = new MCStudies::OfficialNIWGStack("NuMu_Cosmu_MC_MODES", "NuMu_Cosmu_MC_MODES;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_NuMu_Cosmu); fHist_NuMu_Cosmu_Pions = new MCStudies::OfficialPionStack("NuMu_Cosmu_MC_PIONS", "NuMu_Cosmu_MC_PIONS;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_NuMu_Cosmu); fHist_NuMuBar_Enu = new TH1D("NuMuBar_Enu_MC", "NuMuBar_Enu_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_NuMuBar_Enu_Modes = new MCStudies::OfficialNIWGStack("NuMuBar_Enu_MC_MODES", "NuMuBar_Enu_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMuBar_Enu); fHist_NuMuBar_Enu_Pions = new MCStudies::OfficialPionStack("NuMuBar_Enu_MC_PIONS", "NuMuBar_Enu_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMuBar_Enu); fHist_NuMuBar_EnuRates = new TH1D("NuMuBar_EnuRates_MC", "NuMuBar_EnuRates_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_NuMuBar_EnuRates_Modes = new MCStudies::OfficialNIWGStack("NuMuBar_EnuRates_MC_MODES", "NuMuBar_EnuRates_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMuBar_Enu); fHist_NuMuBar_EnuRates_Pions = new MCStudies::OfficialPionStack("NuMuBar_EnuRates_MC_PIONS", "NuMuBar_EnuRates_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NuMuBar_Enu); fHist_NuMuBar_Q2 = new TH1D("NuMuBar_Q2_MC", "NuMuBar_Q2_MC;Q^{2} (GeV);d#sigma/dQ^{2} (cm^{2}/nucleon/GeV^{2})", 45, 0.0, 3.0); fHist_NuMuBar_Q2_Modes = new MCStudies::OfficialNIWGStack("NuMuBar_Q2_MC_MODES", "NuMuBar_Q2_MC_MODES;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_NuMuBar_Q2); fHist_NuMuBar_Q2_Pions = new MCStudies::OfficialPionStack("NuMuBar_Q2_MC_PIONS", "NuMuBar_Q2_MC_PIONS;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_NuMuBar_Q2); fHist_NuMuBar_Pmu = new TH1D("NuMuBar_Pmu_MC", "NuMuBar_Pmu_MC;P_{#mu} (GeV);#sigma (cm^{2}/nucleon/GeV)", 45, 0.0, 3.0); fHist_NuMuBar_Pmu_Modes = new MCStudies::OfficialNIWGStack("NuMuBar_Pmu_MC_MODES", "NuMuBar_Pmu_MC_MODES;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_NuMuBar_Pmu); fHist_NuMuBar_Pmu_Pions = new MCStudies::OfficialPionStack("NuMuBar_Pmu_MC_PIONS", "NuMuBar_Pmu_MC_PIONS;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_NuMuBar_Pmu); fHist_NuMuBar_Cosmu = new TH1D("NuMuBar_Cosmu_MC", "NuMuBar_Cosmu_MC;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", 45, -1.0, 1.0); fHist_NuMuBar_Cosmu_Modes = new MCStudies::OfficialNIWGStack("NuMuBar_Cosmu_MC_MODES", "NuMuBar_Cosmu_MC_MODES;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_NuMuBar_Cosmu); fHist_NuMuBar_Cosmu_Pions = new MCStudies::OfficialPionStack("NuMuBar_Cosmu_MC_PIONS", "NuMuBar_Cosmu_MC_PIONS;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_NuMuBar_Cosmu); fHist_Nue_Enu = new TH1D("Nue_Enu_MC", "Nue_Enu_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_Nue_Enu_Modes = new MCStudies::OfficialNIWGStack("Nue_Enu_MC_MODES", "Nue_Enu_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_Nue_Enu); fHist_Nue_Enu_Pions = new MCStudies::OfficialPionStack("Nue_Enu_MC_PIONS", "Nue_Enu_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_Nue_Enu); fHist_Nue_EnuRates = new TH1D("Nue_EnuRates_MC", "Nue_EnuRates_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_Nue_EnuRates_Modes = new MCStudies::OfficialNIWGStack("Nue_EnuRates_MC_MODES", "Nue_EnuRates_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_Nue_Enu); fHist_Nue_EnuRates_Pions = new MCStudies::OfficialPionStack("Nue_EnuRates_MC_PIONS", "Nue_EnuRates_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_Nue_Enu); fHist_Nue_Q2 = new TH1D("Nue_Q2_MC", "Nue_Q2_MC;Q^{2} (GeV);d#sigma/dQ^{2} (cm^{2}/nucleon/GeV^{2})", 45, 0.0, 3.0); fHist_Nue_Q2_Modes = new MCStudies::OfficialNIWGStack("Nue_Q2_MC_MODES", "Nue_Q2_MC_MODES;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_Nue_Q2); fHist_Nue_Q2_Pions = new MCStudies::OfficialPionStack("Nue_Q2_MC_PIONS", "Nue_Q2_MC_PIONS;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_Nue_Q2); fHist_Nue_Pmu = new TH1D("Nue_Pmu_MC", "Nue_Pmu_MC;P_{#mu} (GeV);#sigma (cm^{2}/nucleon/GeV)", 45, 0.0, 3.0); fHist_Nue_Pmu_Modes = new MCStudies::OfficialNIWGStack("Nue_Pmu_MC_MODES", "Nue_Pmu_MC_MODES;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_Nue_Pmu); fHist_Nue_Pmu_Pions = new MCStudies::OfficialPionStack("Nue_Pmu_MC_PIONS", "Nue_Pmu_MC_PIONS;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_Nue_Pmu); fHist_Nue_Cosmu = new TH1D("Nue_Cosmu_MC", "Nue_Cosmu_MC;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", 45, -1.0, 1.0); fHist_Nue_Cosmu_Modes = new MCStudies::OfficialNIWGStack("Nue_Cosmu_MC_MODES", "Nue_Cosmu_MC_MODES;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_Nue_Cosmu); fHist_Nue_Cosmu_Pions = new MCStudies::OfficialPionStack("Nue_Cosmu_MC_PIONS", "Nue_Cosmu_MC_PIONS;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_Nue_Cosmu); fHist_NueBar_Enu = new TH1D("NueBar_Enu_MC", "NueBar_Enu_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_NueBar_Enu_Modes = new MCStudies::OfficialNIWGStack("NueBar_Enu_MC_MODES", "NueBar_Enu_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NueBar_Enu); fHist_NueBar_Enu_Pions = new MCStudies::OfficialPionStack("NueBar_Enu_MC_PIONS", "NueBar_Enu_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NueBar_Enu); fHist_NueBar_EnuRates = new TH1D("NueBar_EnuRates_MC", "NueBar_EnuRates_MC;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", 45, 0.0, 3.0); fHist_NueBar_EnuRates_Modes = new MCStudies::OfficialNIWGStack("NueBar_EnuRates_MC_MODES", "NueBar_EnuRates_MC_MODES;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NueBar_Enu); fHist_NueBar_EnuRates_Pions = new MCStudies::OfficialPionStack("NueBar_EnuRates_MC_PIONS", "NueBar_EnuRates_MC_PIONS;E_{#nu}^{True} (GeV);#sigma (cm^{2}/nucleon)", fHist_NueBar_Enu); fHist_NueBar_Q2 = new TH1D("NueBar_Q2_MC", "NueBar_Q2_MC;Q^{2} (GeV);d#sigma/dQ^{2} (cm^{2}/nucleon/GeV^{2})", 45, 0.0, 3.0); fHist_NueBar_Q2_Modes = new MCStudies::OfficialNIWGStack("NueBar_Q2_MC_MODES", "NueBar_Q2_MC_MODES;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_NueBar_Q2); fHist_NueBar_Q2_Pions = new MCStudies::OfficialPionStack("NueBar_Q2_MC_PIONS", "NueBar_Q2_MC_PIONS;Q^{2} (GeV);#sigma (cm^{2}/nucleon/GeV^{2})", fHist_NueBar_Q2); fHist_NueBar_Pmu = new TH1D("NueBar_Pmu_MC", "NueBar_Pmu_MC;P_{#mu} (GeV);#sigma (cm^{2}/nucleon/GeV)", 45, 0.0, 3.0); fHist_NueBar_Pmu_Modes = new MCStudies::OfficialNIWGStack("NueBar_Pmu_MC_MODES", "NueBar_Pmu_MC_MODES;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_NueBar_Pmu); fHist_NueBar_Pmu_Pions = new MCStudies::OfficialPionStack("NueBar_Pmu_MC_PIONS", "NueBar_Pmu_MC_PIONS;P_{#mu} (GeV);d#sigma/dP_{#mu} (cm^{2}/nucleon/GeV)", fHist_NueBar_Pmu); fHist_NueBar_Cosmu = new TH1D("NueBar_Cosmu_MC", "NueBar_Cosmu_MC;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", 45, -1.0, 1.0); fHist_NueBar_Cosmu_Modes = new MCStudies::OfficialNIWGStack("NueBar_Cosmu_MC_MODES", "NueBar_Cosmu_MC_MODES;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_NueBar_Cosmu); fHist_NueBar_Cosmu_Pions = new MCStudies::OfficialPionStack("NueBar_Cosmu_MC_PIONS", "NueBar_Cosmu_MC_PIONS;cos#theta_{#mu};d#sigma/dcos#theta_{#mu} (cm^{2}/nucleon)", fHist_NueBar_Cosmu); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void OfficialNIWGPlots::FillEventVariables(FitEvent *event) { //******************************************************************** if (!event->GetNeutrinoIn() or abs(event->Mode) > 30) return; TLorentzVector vectnu = event->GetNeutrinoIn()->fP; int leptons[] = {12,-12,14,-14}; if (!event->GetHMFSParticle(leptons))return; TLorentzVector vectlep = event->GetHMFSParticle(leptons)->fP; double Q2 = fabs((vectlep - vectnu).Mag2()) / 1.E6; double Enu = vectnu.E() / 1.E3; double Pmu = vectlep.Vect().Mag() / 1.E3; double Cosmu = cos(vectlep.Vect().Angle(vectnu.Vect())); bool nue = (abs(event->GetNeutrinoIn()->fPID) == 12); bool nubar = (event->GetNeutrinoIn()->fPID > 0); if (OfficialNIWGPlots::isSignal(event)) { if (!nue and !nubar) { fHist_NuMu_Enu->Fill(Enu, Weight); fHist_NuMu_Enu_Modes->Fill(event, Enu, Weight); fHist_NuMu_Enu_Pions->Fill(event, Enu, Weight); fHist_NuMu_EnuRates->Fill(Enu, Weight); fHist_NuMu_EnuRates_Modes->Fill(event, Enu, Weight); fHist_NuMu_EnuRates_Pions->Fill(event, Enu, Weight); fHist_NuMu_Q2->Fill(Q2, Weight); fHist_NuMu_Q2_Modes->Fill(event, Q2, Weight); fHist_NuMu_Q2_Pions->Fill(event, Q2, Weight); fHist_NuMu_Pmu->Fill(Pmu, Weight); fHist_NuMu_Pmu_Modes->Fill(event, Pmu, Weight); fHist_NuMu_Pmu_Pions->Fill(event, Pmu, Weight); fHist_NuMu_Cosmu->Fill(Cosmu, Weight); fHist_NuMu_Cosmu_Modes->Fill(event, Cosmu, Weight); fHist_NuMu_Cosmu_Pions->Fill(event, Cosmu, Weight); } else if (!nue and nubar){ fHist_NuMuBar_Enu->Fill(Enu, Weight); fHist_NuMuBar_Enu_Modes->Fill(event, Enu, Weight); fHist_NuMuBar_Enu_Pions->Fill(event, Enu, Weight); fHist_NuMuBar_EnuRates->Fill(Enu, Weight); fHist_NuMuBar_EnuRates_Modes->Fill(event, Enu, Weight); fHist_NuMuBar_EnuRates_Pions->Fill(event, Enu, Weight); fHist_NuMuBar_Q2->Fill(Q2, Weight); fHist_NuMuBar_Q2_Modes->Fill(event, Q2, Weight); fHist_NuMuBar_Q2_Pions->Fill(event, Q2, Weight); fHist_NuMuBar_Pmu->Fill(Pmu, Weight); fHist_NuMuBar_Pmu_Modes->Fill(event, Pmu, Weight); fHist_NuMuBar_Pmu_Pions->Fill(event, Pmu, Weight); fHist_NuMuBar_Cosmu->Fill(Cosmu, Weight); fHist_NuMuBar_Cosmu_Modes->Fill(event, Cosmu, Weight); fHist_NuMuBar_Cosmu_Pions->Fill(event, Cosmu, Weight); } else if (nue and !nubar){ fHist_Nue_Enu->Fill(Enu, Weight); fHist_Nue_Enu_Modes->Fill(event, Enu, Weight); fHist_Nue_Enu_Pions->Fill(event, Enu, Weight); fHist_Nue_EnuRates->Fill(Enu, Weight); fHist_Nue_EnuRates_Modes->Fill(event, Enu, Weight); fHist_Nue_EnuRates_Pions->Fill(event, Enu, Weight); fHist_Nue_Q2->Fill(Q2, Weight); fHist_Nue_Q2_Modes->Fill(event, Q2, Weight); fHist_Nue_Q2_Pions->Fill(event, Q2, Weight); fHist_Nue_Pmu->Fill(Pmu, Weight); fHist_Nue_Pmu_Modes->Fill(event, Pmu, Weight); fHist_Nue_Pmu_Pions->Fill(event, Pmu, Weight); fHist_Nue_Cosmu->Fill(Cosmu, Weight); fHist_Nue_Cosmu_Modes->Fill(event, Cosmu, Weight); fHist_Nue_Cosmu_Pions->Fill(event, Cosmu, Weight); } else if (nue and nubar){ fHist_NueBar_Enu->Fill(Enu, Weight); fHist_NueBar_Enu_Modes->Fill(event, Enu, Weight); fHist_NueBar_Enu_Pions->Fill(event, Enu, Weight); fHist_NueBar_EnuRates->Fill(Enu, Weight); fHist_NueBar_EnuRates_Modes->Fill(event, Enu, Weight); fHist_NueBar_EnuRates_Pions->Fill(event, Enu, Weight); fHist_NueBar_Q2->Fill(Q2, Weight); fHist_NueBar_Q2_Modes->Fill(event, Q2, Weight); fHist_NueBar_Q2_Pions->Fill(event, Q2, Weight); fHist_NueBar_Pmu->Fill(Pmu, Weight); fHist_NueBar_Pmu_Modes->Fill(event, Pmu, Weight); fHist_NueBar_Pmu_Pions->Fill(event, Pmu, Weight); fHist_NueBar_Cosmu->Fill(Cosmu, Weight); fHist_NueBar_Cosmu_Modes->Fill(event, Cosmu, Weight); fHist_NueBar_Cosmu_Pions->Fill(event, Cosmu, Weight); } } return; }; //******************************************************************** void OfficialNIWGPlots::Write(std::string drawOpt) { //******************************************************************** LOG(FIT) << "Writing OfficialNIWGPlots " << std::endl; fHist_NuMu_Enu->Write(); fHist_NuMu_Enu_Modes->Write(); fHist_NuMu_Enu_Pions->Write(); fHist_NuMu_EnuRates->Write(); fHist_NuMu_EnuRates_Modes->Write(); fHist_NuMu_EnuRates_Pions->Write(); fHist_NuMu_Q2->Write(); fHist_NuMu_Q2_Modes->Write(); fHist_NuMu_Q2_Pions->Write(); fHist_NuMu_Pmu->Write(); fHist_NuMu_Pmu_Modes->Write(); fHist_NuMu_Pmu_Pions->Write(); fHist_NuMu_Cosmu->Write(); fHist_NuMu_Cosmu_Modes->Write(); fHist_NuMu_Cosmu_Pions->Write(); fHist_NuMuBar_Enu->Write(); fHist_NuMuBar_Enu_Modes->Write(); fHist_NuMuBar_Enu_Pions->Write(); fHist_NuMuBar_EnuRates->Write(); fHist_NuMuBar_EnuRates_Modes->Write(); fHist_NuMuBar_EnuRates_Pions->Write(); fHist_NuMuBar_Q2->Write(); fHist_NuMuBar_Q2_Modes->Write(); fHist_NuMuBar_Q2_Pions->Write(); fHist_NuMuBar_Pmu->Write(); fHist_NuMuBar_Pmu_Modes->Write(); fHist_NuMuBar_Pmu_Pions->Write(); fHist_NuMuBar_Cosmu->Write(); fHist_NuMuBar_Cosmu_Modes->Write(); fHist_NuMuBar_Cosmu_Pions->Write(); fHist_Nue_Enu->Write(); fHist_Nue_Enu_Modes->Write(); fHist_Nue_Enu_Pions->Write(); fHist_Nue_EnuRates->Write(); fHist_Nue_EnuRates_Modes->Write(); fHist_Nue_EnuRates_Pions->Write(); fHist_Nue_Q2->Write(); fHist_Nue_Q2_Modes->Write(); fHist_Nue_Q2_Pions->Write(); fHist_Nue_Pmu->Write(); fHist_Nue_Pmu_Modes->Write(); fHist_Nue_Pmu_Pions->Write(); fHist_Nue_Cosmu->Write(); fHist_Nue_Cosmu_Modes->Write(); fHist_Nue_Cosmu_Pions->Write(); fHist_NueBar_Enu->Write(); fHist_NueBar_Enu_Modes->Write(); fHist_NueBar_Enu_Pions->Write(); fHist_NueBar_EnuRates->Write(); fHist_NueBar_EnuRates_Modes->Write(); fHist_NueBar_EnuRates_Pions->Write(); fHist_NueBar_Q2->Write(); fHist_NueBar_Q2_Modes->Write(); fHist_NueBar_Q2_Pions->Write(); fHist_NueBar_Pmu->Write(); fHist_NueBar_Pmu_Modes->Write(); fHist_NueBar_Pmu_Pions->Write(); fHist_NueBar_Cosmu->Write(); fHist_NueBar_Cosmu_Modes->Write(); fHist_NueBar_Cosmu_Pions->Write(); return; } //******************************************************************** void OfficialNIWGPlots::ResetAll() { //******************************************************************** fHist_NuMu_Enu->Reset(); fHist_NuMu_Enu_Modes->Reset(); fHist_NuMu_Enu_Pions->Reset(); fHist_NuMu_EnuRates->Reset(); fHist_NuMu_EnuRates_Modes->Reset(); fHist_NuMu_EnuRates_Pions->Reset(); fHist_NuMu_Q2->Reset(); fHist_NuMu_Q2_Modes->Reset(); fHist_NuMu_Q2_Pions->Reset(); fHist_NuMu_Pmu->Reset(); fHist_NuMu_Pmu_Modes->Reset(); fHist_NuMu_Pmu_Pions->Reset(); fHist_NuMu_Cosmu->Reset(); fHist_NuMu_Cosmu_Modes->Reset(); fHist_NuMu_Cosmu_Pions->Reset(); fHist_NuMuBar_Enu->Reset(); fHist_NuMuBar_Enu_Modes->Reset(); fHist_NuMuBar_Enu_Pions->Reset(); fHist_NuMuBar_EnuRates->Reset(); fHist_NuMuBar_EnuRates_Modes->Reset(); fHist_NuMuBar_EnuRates_Pions->Reset(); fHist_NuMuBar_Q2->Reset(); fHist_NuMuBar_Q2_Modes->Reset(); fHist_NuMuBar_Q2_Pions->Reset(); fHist_NuMuBar_Pmu->Reset(); fHist_NuMuBar_Pmu_Modes->Reset(); fHist_NuMuBar_Pmu_Pions->Reset(); fHist_NuMuBar_Cosmu->Reset(); fHist_NuMuBar_Cosmu_Modes->Reset(); fHist_NuMuBar_Cosmu_Pions->Reset(); fHist_Nue_Enu->Reset(); fHist_Nue_Enu_Modes->Reset(); fHist_Nue_Enu_Pions->Reset(); fHist_Nue_EnuRates->Reset(); fHist_Nue_EnuRates_Modes->Reset(); fHist_Nue_EnuRates_Pions->Reset(); fHist_Nue_Q2->Reset(); fHist_Nue_Q2_Modes->Reset(); fHist_Nue_Q2_Pions->Reset(); fHist_Nue_Pmu->Reset(); fHist_Nue_Pmu_Modes->Reset(); fHist_Nue_Pmu_Pions->Reset(); fHist_Nue_Cosmu->Reset(); fHist_Nue_Cosmu_Modes->Reset(); fHist_Nue_Cosmu_Pions->Reset(); fHist_NueBar_Enu->Reset(); fHist_NueBar_Enu_Modes->Reset(); fHist_NueBar_Enu_Pions->Reset(); fHist_NueBar_EnuRates->Reset(); fHist_NueBar_EnuRates_Modes->Reset(); fHist_NueBar_EnuRates_Pions->Reset(); fHist_NueBar_Q2->Reset(); fHist_NueBar_Q2_Modes->Reset(); fHist_NueBar_Q2_Pions->Reset(); fHist_NueBar_Pmu->Reset(); fHist_NueBar_Pmu_Modes->Reset(); fHist_NueBar_Pmu_Pions->Reset(); fHist_NueBar_Cosmu->Reset(); fHist_NueBar_Cosmu_Modes->Reset(); fHist_NueBar_Cosmu_Pions->Reset(); return; } //******************************************************************** void OfficialNIWGPlots::ScaleEvents() { //******************************************************************** fHist_NuMu_EnuRates->Scale(fScaleFactorEnuXSec); fHist_NuMu_EnuRates_Modes->Scale(fScaleFactorEnuXSec); fHist_NuMu_EnuRates_Pions->Scale(fScaleFactorEnuXSec); PlotUtils::FluxUnfoldedScaling(fHist_NuMu_Enu, GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); - fHist_NuMu_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); - fHist_NuMu_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); + fHist_NuMu_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); + fHist_NuMu_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); fHist_NuMu_Q2->Scale(fScaleFactor, "width"); fHist_NuMu_Q2_Modes->Scale(fScaleFactor, "width"); fHist_NuMu_Q2_Pions->Scale(fScaleFactor, "width"); fHist_NuMu_Pmu->Scale(fScaleFactor, "width"); fHist_NuMu_Pmu_Modes->Scale(fScaleFactor, "width"); fHist_NuMu_Pmu_Pions->Scale(fScaleFactor, "width"); fHist_NuMu_Cosmu->Scale(fScaleFactor, "width"); fHist_NuMu_Cosmu_Modes->Scale(fScaleFactor, "width"); fHist_NuMu_Cosmu_Pions->Scale(fScaleFactor, "width"); fHist_NuMu_EnuRates->Scale(fScaleFactorEnuXSec); fHist_NuMu_EnuRates_Modes->Scale(fScaleFactorEnuXSec); fHist_NuMu_EnuRates_Pions->Scale(fScaleFactorEnuXSec); PlotUtils::FluxUnfoldedScaling(fHist_NuMuBar_Enu, GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); - fHist_NuMuBar_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); - fHist_NuMuBar_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); + fHist_NuMuBar_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); + fHist_NuMuBar_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); fHist_NuMuBar_Q2->Scale(fScaleFactor, "width"); fHist_NuMuBar_Q2_Modes->Scale(fScaleFactor, "width"); fHist_NuMuBar_Q2_Pions->Scale(fScaleFactor, "width"); fHist_NuMuBar_Pmu->Scale(fScaleFactor, "width"); fHist_NuMuBar_Pmu_Modes->Scale(fScaleFactor, "width"); fHist_NuMuBar_Pmu_Pions->Scale(fScaleFactor, "width"); fHist_NuMuBar_Cosmu->Scale(fScaleFactor, "width"); fHist_NuMuBar_Cosmu_Modes->Scale(fScaleFactor, "width"); fHist_NuMuBar_Cosmu_Pions->Scale(fScaleFactor, "width"); fHist_Nue_EnuRates->Scale(fScaleFactorEnuXSec); fHist_Nue_EnuRates_Modes->Scale(fScaleFactorEnuXSec); fHist_Nue_EnuRates_Pions->Scale(fScaleFactorEnuXSec); PlotUtils::FluxUnfoldedScaling(fHist_Nue_Enu, GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); - fHist_Nue_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); - fHist_Nue_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); + fHist_Nue_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); + fHist_Nue_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); fHist_Nue_Q2->Scale(fScaleFactor, "width"); fHist_Nue_Q2_Modes->Scale(fScaleFactor, "width"); fHist_Nue_Q2_Pions->Scale(fScaleFactor, "width"); fHist_Nue_Pmu->Scale(fScaleFactor, "width"); fHist_Nue_Pmu_Modes->Scale(fScaleFactor, "width"); fHist_Nue_Pmu_Pions->Scale(fScaleFactor, "width"); fHist_Nue_Cosmu->Scale(fScaleFactor, "width"); fHist_Nue_Cosmu_Modes->Scale(fScaleFactor, "width"); fHist_Nue_Cosmu_Pions->Scale(fScaleFactor, "width"); fHist_Nue_EnuRates->Scale(fScaleFactorEnuXSec); fHist_Nue_EnuRates_Modes->Scale(fScaleFactorEnuXSec); fHist_Nue_EnuRates_Pions->Scale(fScaleFactorEnuXSec); PlotUtils::FluxUnfoldedScaling(fHist_NueBar_Enu, GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); - fHist_NueBar_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); - fHist_NueBar_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec); + fHist_NueBar_Enu_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); + fHist_NueBar_Enu_Pions->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactorEnuXSec, fNEvents); fHist_NueBar_Q2->Scale(fScaleFactor, "width"); fHist_NueBar_Q2_Modes->Scale(fScaleFactor, "width"); fHist_NueBar_Q2_Pions->Scale(fScaleFactor, "width"); fHist_NueBar_Pmu->Scale(fScaleFactor, "width"); fHist_NueBar_Pmu_Modes->Scale(fScaleFactor, "width"); fHist_NueBar_Pmu_Pions->Scale(fScaleFactor, "width"); fHist_NueBar_Cosmu->Scale(fScaleFactor, "width"); fHist_NueBar_Cosmu_Modes->Scale(fScaleFactor, "width"); fHist_NueBar_Cosmu_Pions->Scale(fScaleFactor, "width"); return; } //******************************************************************** /// Select only events with final state Muons bool OfficialNIWGPlots::isSignal(FitEvent *event) { //******************************************************************** if (abs(event->Mode) > 30) return false; // Do we want any other signal? return true; }; /// Functions to deal with the SB mode stacks MCStudies::OfficialNIWGStack::OfficialNIWGStack(std::string name, std::string title, TH1* hist) { fName = name; fTitle = title; AddMode(0, "CC0PI", "CC-0#pi", kRed, 2, 1001); AddMode(1, "CC1PI", "CC-1#pi", kBlue, 2, 1001); AddMode(2, "CCOther", "CC-Other", kMagenta, 2, 1001); AddMode(3, "CCCOH", "CC-Coherent", kGreen, 2, 1001); AddMode(4, "OTHER", "Other", kYellow, 2, 1001); StackBase::SetupStack(hist); }; int MCStudies::OfficialNIWGStack::ConvertModeToIndex(FitEvent* event) { // Other if (((event->NumFSParticle(13) + event->NumFSParticle(-13)) != 1)) return 4; // CC 0 PI if ((event->NumFSMesons() == 0)) return 0; // CC Coherent if (abs(event->Mode) == 16) return 3; // CC 1 PI if ((event->NumFSParticle(PhysConst::pdg_charged_pions)) == 1) return 1; // CC Other return 2; }; void MCStudies::OfficialNIWGStack::Fill(FitEvent* evt, double x, double y, double z, double weight) { StackBase::FillStack(this->ConvertModeToIndex(evt), x, y, z, weight); }; /// Functions to deal with the SB mode stacks MCStudies::OfficialPionStack::OfficialPionStack(std::string name, std::string title, TH1* hist) { fName = name; fTitle = title; AddMode(0, "CC0PI", "CC-0#pi", kRed, 2, 1001); AddMode(1, "CC1PI", "CC-1#pi", kBlue, 2, 1001); AddMode(2, "CCNPI", "CC-N#pi", kGreen, 2, 1001); AddMode(3, "OTHER", "Other", kYellow, 2, 1001); StackBase::SetupStack(hist); }; int MCStudies::OfficialPionStack::ConvertModeToIndex(FitEvent* event) { // Other if (((event->NumFSParticle(13) + event->NumFSParticle(-13)) != 1)) return 3; // CC 0 PI int npi = (event->NumFSParticle(PhysConst::pdg_charged_pions)); if (npi == 0) return 0; if (npi == 1) return 1; if (npi > 1) return 2; return 3; }; void MCStudies::OfficialPionStack::Fill(FitEvent* evt, double x, double y, double z, double weight) { StackBase::FillStack(this->ConvertModeToIndex(evt), x, y, z, weight); }; diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h index a72e261..11e2d5f 100644 --- a/src/Utils/FitUtils.h +++ b/src/Utils/FitUtils.h @@ -1,188 +1,185 @@ // 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 FITUTILS_H_SEEN #define FITUTILS_H_SEEN #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "FitEvent.h" #include "TGraph.h" #include "TH2Poly.h" #include "FitEvent.h" #include "FitLogger.h" /*! * \addtogroup Utils * @{ */ /// Functions needed by individual samples for calculating kinematic quantities. namespace FitUtils { /// Return a vector of all values saved in map double *GetArrayFromMap(std::vector invals, std::map inmap); /// Returns kinetic energy of particle double T(TLorentzVector part); /// Returns momentum of particle double p(TLorentzVector part); double p(FitParticle* part); /// Returns angle between particles (_NOT_ cosine!) double th(TLorentzVector part, TLorentzVector part2); double th(FitParticle* part1, FitParticle* part2); /// Hadronic mass reconstruction double Wrec(TLorentzVector pnu, TLorentzVector pmu); /// Hadronic mass true from initial state particles and muon; useful if the full /// FSI vectors aren't not saved and we for some reasons need W_true double Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc); double SumKE_PartVect(std::vector const fps); double SumTE_PartVect(std::vector const fps); /// Return E Hadronic for all FS Particles in Hadronic System double GetErecoil_TRUE(FitEvent *event); /// Return E Hadronic for all Charged FS Particles in Hadronic System double GetErecoil_CHARGED(FitEvent *event); /* CCQE MiniBooNE/MINERvA */ /// Function to calculate the reconstructed Q^{2}_{QE} -double Q2QErec(TLorentzVector pmu, double costh, double binding, - bool neutrino = true); +double Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true); /// Function returns the reconstructed E_{nu} values -double EnuQErec(TLorentzVector pmu, double costh, double binding, - bool neutrino = true); +double EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true); //! Function to calculate the reconstructed Q^{2}_{QE} -double Q2QErec(double pl, double costh, double binding, - bool neutrino = true); +double Q2QErec(double pl, double costh, double binding, bool neutrino = true); //! Function to calculate the reconstructed Q^{2}_{QE} -double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino); +double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino = true); //! Function returns the reconstructed E_{nu} values double EnuQErec(double pl, double costh, double binding, bool neutrino = true); /* CCQE1p MINERvA */ /// Reconstruct Q2QE given just the maximum energy proton. double ProtonQ2QErec(double pE, double binding); /* E Recoil MINERvA */ double GetErecoil_MINERvA_LowRecoil(FitEvent *event); /* CC1pi0 MiniBooNE */ /// Reconstruct Enu from CCpi0 vectors and binding energy double EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0)); /// Reconstruct Q2 from CCpi0 vectors and binding energy double Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0)); /* CC1pi+ MiniBooNE */ /// returns reconstructed Enu a la MiniBooNE CCpi+ /// returns reconstructed Enu a la MiniBooNE CCpi+ // Also for when not having pion info (so when we have a Michel tag in T2K) double EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip, bool pionInfo = true); /// returns reconstructed Enu assumming resonance interaction where intermediate /// resonance was a Delta double EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu); /// returns reconstructed in a variety of flavours double Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip, int enuType = 0, bool pionInfo = true); /* T2K CC1pi+ on CH */ double thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi); double q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi); double WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip); double EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi); /* nucleon single pion */ double MpPi(TLorentzVector pp, TLorentzVector ppi); /// Gets delta p T as defined in Phys.Rev. C94 (2016) no.1, 015503 double Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi); /// Gets delta phi T as defined in Phys.Rev. C94 (2016) no.1, 015503 double Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi); /// Gets delta alpha T as defined in Phys.Rev. C94 (2016) no.1, 015503 double Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi); // As defined in PhysRevC.95.065501 // Using prescription from arXiv 1805.05486 double Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi); //For T2K inferred kinematics analyis - variables defined as on page 7 of T2K TN287v11 (and now arXiv 1802.05078) double ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino); TVector3 tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino); double cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino); double CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot); double PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot); } /*! @} */ #endif diff --git a/src/Utils/PlotUtils.cxx b/src/Utils/PlotUtils.cxx index 73d56fa..24f1c60 100644 --- a/src/Utils/PlotUtils.cxx +++ b/src/Utils/PlotUtils.cxx @@ -1,1033 +1,1033 @@ // 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)); // Now Get a flux PDF assuming X axis is Enu TH1D* pdfflux = (TH1D*)fMCHist->ProjectionX()->Clone(); // pdfflux->Write( (std::string(fMCHist->GetName()) + "_PROJX").c_str()); pdfflux->Reset(); // Awful MiniBooNE Check for the time being 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); } for (int i = 0; i < fMCHist->GetNbinsX(); i++) { for (int j = 0; j < fMCHist->GetNbinsY(); j++) { if (pdfflux->GetBinContent(i + 1) == 0.0) continue; double binWidth = fMCHist->GetYaxis()->GetBinLowEdge(j + 2) - fMCHist->GetYaxis()->GetBinLowEdge(j + 1); 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); } } delete eventhist; delete fFluxHist; return; }; 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; + //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; return; }; // 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 or 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()) { THROW(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) { ERR(FTL) << "Could not find distribution " << name << " in file " << file << std::endl; throw; } 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()) { THROW("Couldn't open root file: \"" << file << "\"."); } TH1* tempHist = dynamic_cast(rootHistFile->Get(name.c_str())->Clone()); if (!tempHist) { THROW("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()) { THROW("Couldn't open root file: \"" << file << "\"."); } TDirectory* newdir = gDirectory; TGraph* temp = dynamic_cast(rootHistFile->Get(name.c_str())->Clone()); if (!temp) { THROW("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()) { THROW("Couldn't find input file when attempting to parse : \"" << 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()) { THROW( "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()) { THROW("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) { THROW("Couldn't retrieve: \"" << histnames[i] << "\" from root file: \"" << fname[0] << "\"."); } tempHist->SetDirectory(0); hists.push_back(tempHist); } rootHistFile->Close(); } return hists; } TH2D* PlotUtils::GetTH2DFromTextFile(std::string file) { /// Contents should be /// Low Edfe return NULL; } 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); LOG(REC) << "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " to 0.0 +- 0.0" << std::endl; } 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); LOG(REC) << "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " " << j + 1 << " to 0.0 +- 0.0" << std::endl; } } 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); TH1D* hist_X = maskedhist->ProjectionX(); delete maskedhist; return hist_X; } //*************************************************** TH1D* PlotUtils::GetProjectionY(TH2D* hist, TH2I* mask) { //*************************************************** TH2D* maskedhist = StatUtils::ApplyHistogramMasking(hist, mask); TH1D* hist_Y = maskedhist->ProjectionY(); delete maskedhist; return hist_Y; }