diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
index f568571..eefdb7b 100755
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
@@ -1,464 +1,466 @@
//Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
/*
Authors: Adrian Orea (v1 2017)
Clarence Wret (v2 2018)
*/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CC0pi_XSec_2D_nu.h"
//********************************************************************
void MINERvA_CC0pi_XSec_2D_nu::SetupDataSettings() {
//********************************************************************
// Set Distribution
// See header file for enum and some descriptions
std::string name = fSettings.GetS("name");
// Have two distributions as of summer 2018
if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) fDist = kPtPz;
else if (!name.compare("MINERvA_CC0pi_XSec_2DptQ2_nu")) fDist = kPtQ2;
// Define what files to use from the dist
std::string datafile = "";
std::string corrfile = "";
std::string titles = "";
std::string distdescript = "";
std::string histname = "";
switch (fDist) {
case (kPtPz):
datafile = "MINERvA/CC0pi_2D/MINERvA_LE_CCQELike_pzmu.root";
corrfile = "MINERvA/CC0pi_2D/MINERvA_LE_CCQELike_pzmu.root";
titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} (GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)";
distdescript = "MINERvA_CC0pi_XSec_2Dptpz_nu sample";
histname = "h_pzmu_ptmu_data_nobck_unfold_effcor_cross_section_CV_WithErr";
break;
case (kPtQ2):
datafile = "MINERvA/CC0pi_2D/pt_q2_data.root";
corrfile = "MINERvA/CC0pi_2D/pt_q2_cov.root";
titles = "MINERvA CC0#pi #nu_{#mu} p_{t} Q^{2}_{QE};p_t{z} (GeV);Q^{2}_{QE} (GeV^{2});d^{2}#sigma/dP_{t}dQ^{2}_{QE} (cm^{2}/GeV^{3}/nucleon)";
distdescript = "MINERvA_CC0pi_XSec_2DptQ2_nu sample";
histname = "h_q2_ptmu_data_nobck_unfold_effcor_cross_section_CV_WithErr";
break;
default:
THROW("Unknown Analysis Distribution : " << fDist);
}
fSettings.SetTitle( GeneralUtils::ParseToStr(titles,";")[0] );
fSettings.SetXTitle( GeneralUtils::ParseToStr(titles,";")[1] );
fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[2] );
fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[3] );
// Sample overview ---------------------------------------------------
std::string descrip = distdescript + "\n"\
"Target: CH \n" \
"Flux: MINERvA Low Energy FHC numu \n" \
"Signal: CC-0pi \n";
fSettings.SetDescription(descrip);
// The input ROOT file
fSettings.SetDataInput( FitPar::GetDataBase() + datafile);
// The covariance matrix ROOT file
//if (fDist == kPtPz) {
- fSettings.SetCovarInput( FitPar::GetDataBase() + corrfile);
+ fSettings.SetCovarInput( FitPar::GetDataBase() + corrfile);
//} else {
//ERR(WRN) << " no covariance matrix available for PtQ2" << std::endl;
//ERR(WRN) << " ask Dan Ruterbories nicely and he may provide one" << std::endl;
//}
// Set the data file
SetDataValues(fSettings.GetDataInput(), histname);
}
//********************************************************************
MINERvA_CC0pi_XSec_2D_nu::MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey) {
//********************************************************************
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
SetupDataSettings();
FinaliseSampleSettings();
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux();
// Data is __NOT__ bin width normalised, so override the default
- fIsNoWidth = true;
+ //fIsNoWidth = true;
// Set the mapping values and the covariance matrix files
//SetMapValuesFromText( fSettings.GetMapInput() );
// Also have to make our own covariance matrix to exclude the under and overflow
//if (fDist == kPtPz) {
//TMatrixDSym* tempmat = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "TMatrixDBase");
TMatrixDSym* tempmat = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "CovMatrix");
fFullCovar = tempmat;
// Decomposition is stable for entries that aren't E-xxx
- double ScalingFactor = 1E38;
+ double ScalingFactor = 1E38*1E38;
(*fFullCovar) *= ScalingFactor;
//} else {
//SetCovarFromDiagonal();
//}
/*
// Now we cut out every first and last to exclude under and overflow
fFullCovar = new TMatrixDSym(fDataHist->GetXaxis()->GetNbins()*fDataHist->GetYaxis()->GetNbins());
// Count the real covariance matrix x and y
int xcounter = 0;
int xbins = fDataHist->GetXaxis()->GetNbins();
int ybins = fDataHist->GetYaxis()->GetNbins();
// Loop over the x bins (underflow adds one, overflow adds one)
for (int i = 0; i < (xbins+2)*(ybins+2); ++i) {
// Skip the under and overflow
if (i < (ybins+2) || i % (ybins+1) == 0 || ((i+1)%(ybins+1) == 0) || i > (ybins+2)*(xbins+1)) {
//std::cout << "Skipping ibin " << i << std::endl;
continue;
}
// The ycounter
int ycounter = 0;
// For one bin of pT we have pZ^2 bins
for (int j = 0; j < (xbins+2)*(ybins+2); ++j) {
// Skip the under and overflow
if (j < (ybins+2) || j % (ybins+1) == 0 || ((j+1)%(ybins+1) == 0) || j > (ybins+2)*(xbins+1)) {
//std::cout << "Skipping jbin " << j << std::endl;
continue;
}
(*fFullCovar)(xcounter, ycounter) = (*tempmat)(i, j);
//std::cout << xcounter << ", " << ycounter << " === " << i << ", " << j << std::endl;
ycounter++;
}
// Check dimensions
if (ycounter != xbins*ybins) {
std::cerr << "Counted " << ycounter << " y bins in cov matrix" << std::endl;
std::cerr << "Whereas there should be " << xbins*ybins << std::endl;
}
xcounter++;
}
// Check dimensions
if (xcounter != xbins*ybins) {
std::cerr << "Counted " << xcounter << " x bins in cov matrix" << std::endl;
std::cerr << "Whereas there should be " << xbins*ybins << std::endl;
}
// Delete the temporary
delete tempmat;
*/
// Just check that the data error and covariance are the same
for (int i = 0; i < fFullCovar->GetNrows(); ++i) {
for (int j = 0; j < fFullCovar->GetNcols(); ++j) {
// Get the global bin
int xbin1, ybin1, zbin1;
fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1);
double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1);
double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1+1);
double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1);
double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1+1);
if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) ||
yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue;
double data_error = fDataHist->GetBinError(xbin1, ybin1);
double cov_error = sqrt((*fFullCovar)(i,i)/ScalingFactor);
if (fabs(data_error - cov_error) > 1E-5) {
std::cerr << "Error on data is different to that of covariance" << std::endl;
ERR(FTL) << "Data error: " << data_error << std::endl;
ERR(FTL) << "Cov error: " << cov_error << std::endl;
ERR(FTL) << "Data/Cov: " << data_error/cov_error << std::endl;
ERR(FTL) << "Data-Cov: " << data_error-cov_error << std::endl;
ERR(FTL) << "For x: " << xlo1 << "-" << xhi1 << std::endl;
ERR(FTL) << "For y: " << ylo1 << "-" << yhi1 << std::endl;
throw;
}
}
}
// Now can make the inverted covariance
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
// Now scale back
(*fFullCovar) *= 1.0/ScalingFactor;
(*covar) *= ScalingFactor;
(*fDecomp) *= ScalingFactor;
// Use a TH2D version of the covariance to be able to use the global bin numbering scheme
+ /*
covar_th2d = new TH2D((fSettings.Title()+"_th2").c_str(), (fSettings.Title()+"_th2").c_str(), covar->GetNrows(), 0, covar->GetNrows(), covar->GetNcols(), 0, covar->GetNcols());
for (int i = 0; i < covar_th2d->GetXaxis()->GetNbins(); ++i) {
for (int j = 0; j < covar_th2d->GetYaxis()->GetNbins(); ++j) {
covar_th2d->SetBinContent(i+1, j+1, (*covar)(i,j));
}
}
+ */
//std::cout << "covar is " << covar_th2d->GetXaxis()->GetNbins() << " x " << covar_th2d->GetYaxis()->GetNbins() << " = " << covar_th2d->GetXaxis()->GetNbins()*covar_th2d->GetYaxis()->GetNbins() << std::endl;
//std::cout << "data is " << fDataHist->GetXaxis()->GetNbins() << " x " << fDataHist->GetYaxis()->GetNbins() << " = " << fDataHist->GetXaxis()->GetNbins()*fDataHist->GetYaxis()->GetNbins() << std::endl;
// Let's make our own mapping histogram
// The covariance matrix is dominant in Pt and sub-dominant in Pz and includes all bins with under and overflow
// So we have 13x12 data/MC bins, and including the under/overflow we have 15x14=210 covariance bins
// i.e. need to cut out the first and last bins of covariance matrix
// Mapping histogram will have same dimensions as the data
/*
fMapHist = (TH2I*)(fDataHist->Clone());
fMapHist->Reset();
std::string MapTitle = std::string(fDataHist->GetName())+"_MAP";
fMapHist->SetNameTitle(MapTitle.c_str(), MapTitle.c_str());
int counter = 1;
for (int i = 0; i <= fDataHist->GetXaxis()->GetNbins()+1; ++i) {
for (int j = 0; j <= fDataHist->GetYaxis()->GetNbins()+1; ++j) {
if (i == 0 || i == fDataHist->GetXaxis()->GetNbins()+1 || j == 0 || j == fDataHist->GetYaxis()->GetNbins()+1) {
fMapHist->SetBinContent(i+1, j+1, 0);
} else {
fMapHist->SetBinContent(i+1, j+1, counter);
counter++;
}
std::cout << fMapHist->GetBinContent(i+1, j+1) << " " << fDataHist->GetBinContent(i+1, j+1) << std::endl;
}
std::cout << std::endl;
}
*/
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC0pi_XSec_2D_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Checking to see if there is a Muon
if (event->NumFSParticle(13) == 0) return;
// Get the muon kinematics
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
Double_t px = Pmu.X()/1000;
Double_t py = Pmu.Y()/1000;
Double_t pt = sqrt(px*px+py*py);
// y-axis is transverse momentum for both measurements
fYVar = pt;
// Now we set the x-axis
switch (fDist) {
case kPtPz:
{
// Don't want to assume the event generators all have neutrino coming along z
// pz is muon momentum projected onto the neutrino direction
Double_t pz = Pmu.Vect().Dot(Pnu.Vect()*(1.0/Pnu.Vect().Mag()))/1000.;
// Set Hist Variables
fXVar = pz;
break;
}
case kPtQ2:
{
// 34 MeV binding energy and neutrino mode (true)
double q2qe = FitUtils::Q2QErec(Pmu, Pnu, 34.0, true);
fXVar = q2qe;
break;
}
default:
THROW("DIST NOT FOUND : " << fDist);
break;
}
};
//********************************************************************
bool MINERvA_CC0pi_XSec_2D_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax);
};
//********************************************************************
// Custom likelihood calculator because binning of covariance matrix
double MINERvA_CC0pi_XSec_2D_nu::GetLikelihood() {
//********************************************************************
// The calculated chi2
double chi2 = 0.0;
// STRICTLY DEBUGGING WITH DAN
/*
- std::vector Names;
+ std::vector Names;
//Names.push_back(std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie_MC.root");
//Names.push_back(std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie2p2hrpa_MC.root");
//Names.push_back(std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/MnvGENIE_MC.root");
Names.push_back(std::string(std::getenv("NUISANCE"))+"/data/MINERvA/CC0pi_2D/MINERvA_LE_CCQELike_pzmu.root");
// Hack hack hack hack
double scaleF = 0.0;
for (size_t a = 0; a < Names.size(); ++a) {
- std::cout << Names[a] << std::endl;
- //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie_MC.root").c_str(), "OPEN");
- //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie2p2hrpa_MC.root").c_str(), "OPEN");
- //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/MnvGENIE_MC.root").c_str(), "OPEN");
- TFile *file = new TFile(Names[a].c_str(), "OPEN");
- //TH2D *mc = (TH2D*)(file->Get("h_pzmu_ptmu_mc_nobck_unfold_effcor_cross_section_CV_WithStatErr")->Clone());
- TH2D *mc = (TH2D*)(file->Get("h_pzmu_ptmu_mc_nobck_unfold_effcor_cross_section_CV_WithErr")->Clone());
- fMCHist = mc;
- //fMCHist->Scale(1., "width");
- //fDataHist->Scale(1., "width");
- */
-
- // Support shape comparisons
- double scaleF = fDataHist->Integral() / fMCHist->Integral();
- if (fIsShape) {
- fMCHist->Scale(scaleF);
- fMCFine->Scale(scaleF);
- //PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF);
- }
+ std::cout << Names[a] << std::endl;
+ //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie_MC.root").c_str(), "OPEN");
+ //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie2p2hrpa_MC.root").c_str(), "OPEN");
+ //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/MnvGENIE_MC.root").c_str(), "OPEN");
+ TFile *file = new TFile(Names[a].c_str(), "OPEN");
+ //TH2D *mc = (TH2D*)(file->Get("h_pzmu_ptmu_mc_nobck_unfold_effcor_cross_section_CV_WithStatErr")->Clone());
+ TH2D *mc = (TH2D*)(file->Get("h_pzmu_ptmu_mc_nobck_unfold_effcor_cross_section_CV_WithErr")->Clone());
+ fMCHist = mc;
+ //fMCHist->Scale(1., "width");
+ //fDataHist->Scale(1., "width");
+ */
- // Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA used for their measurement
- // Can be prettified in due time but for now keep
- bool chi2_use_overflow_err = false;
- const int lowBin = chi2_use_overflow_err?0:1; // Either they both use underflow, or neither of them does
-
- int nbinsx=fMCHist->GetNbinsX();
- int nbinsy=fMCHist->GetNbinsY();
- const int highBinX = nbinsx;
- const int highBinY = nbinsy;
- Int_t Nbins = (highBinX-lowBin+1)* (highBinY-lowBin+1); // from low to high inclusive in each dimension
-
- TMatrixD covMatrixTmp = (*fFullCovar);
-
- TMatrixD covMatrix(Nbins, Nbins);
- const double scaling = 1e80;//ROOT can't seem to handle small entries with lots of zeros? Suggested scaling the histogram and then rescaling the inverted matrix
- for( int i = 0; i != Nbins; ++i ) {
- for( int j = 0; j != Nbins; ++j ) {
- // I think this is right... it checks out on a small sample
- int old_i_bin = (i/nbinsx + 1)* (nbinsx +2) +1 + i%nbinsx;
- int old_j_bin = (j/nbinsx + 1)* (nbinsx +2) +1 + j%nbinsx;
- covMatrix[i][j] = covMatrixTmp[old_i_bin][old_j_bin]*scaling;
- //std::cout << Nbins*Nbins << " (" << Nbins << "*" << Nbins << ")" << std::endl;
- //std::cout << i << ", " << j << " = " << old_i_bin << " " << old_j_bin << std::endl;
- }
- //throw;
- }
- TDecompSVD error(covMatrix);
- TMatrixD errorMatrix(covMatrix);
- if( ! error.Invert( errorMatrix ) ) {
- std::cout << "Cannot invert total covariance matrix. You could use statistical errors only for Chi2 calculation. But it isn't implemented yet" << std::endl;
- }
+ // Support shape comparisons
+ double scaleF = fDataHist->Integral() / fMCHist->Integral();
+ if (fIsShape) {
+ fMCHist->Scale(scaleF);
+ fMCFine->Scale(scaleF);
+ //PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF);
+ }
- covMatrix *= 1/scaling;
- errorMatrix *= scaling;
-
- for( int i = 0; i != Nbins; ++i ) {
- int hist_i_bin = chi2_use_overflow_err?i:((i/nbinsx + 1)* (nbinsx +2) +1 + i%nbinsx); // Translate to the histogram bin, if we aren't using the overflow errors, meaning the covariance matrix is smaller than the histogram
- const Double_t x_data_i = fDataHist->GetBinContent(hist_i_bin);
- const Double_t x_mc_i = fMCHist->GetBinContent(hist_i_bin);
- for( int j = 0; j != Nbins; ++j ) {
- // Each element of the inverted covariance matrix corresponds to a pair of data and MC
- int hist_j_bin = chi2_use_overflow_err?j:((j/nbinsx + 1)* (nbinsx +2) +1 + j%nbinsx);
- const Double_t x_data_j = fDataHist->GetBinContent(hist_j_bin);
- const Double_t x_mc_j = fMCHist->GetBinContent(hist_j_bin);
- const double chi2_ij = (x_data_i - x_mc_i) * errorMatrix[i][j] * (x_data_j - x_mc_j);
- std::cout << x_data_i << "\t" << x_mc_i << "\t" << i << "\t" << errorMatrix[i][j] << "\t" << j << "\t" << x_data_j << "\t" << x_mc_j << "\t" << chi2_ij << std::endl;
-
- chi2 += chi2_ij;
- }
- }
- //std::cout << Names[a] << " chi2 " << chi2 << std::endl;
-
- /*
- * CWRET CALC
- // Calculate the test-statistic
- for (int i = 0; i < covar_th2d->GetXaxis()->GetNbins()+1; ++i) {
- // Get the global bin for x
- int xbin1, ybin1, zbin1;
- fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1);
- double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1);
- double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1+1);
- double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1);
- double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1+1);
- if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
- ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
- xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) ||
- yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue;
-
- // Get the data
- double data1 = fDataHist->GetBinContent(i);
- // Get the MC
- double mc1 = fMCHist->GetBinContent(i);
-
- for (int j = 0; j < covar_th2d->GetYaxis()->GetNbins()+1; ++j) {
-
- int xbin2, ybin2, zbin2;
- fDataHist->GetBinXYZ(j, xbin2, ybin2, zbin2);
- double xlo2 = fDataHist->GetXaxis()->GetBinLowEdge(xbin2);
- double xhi2 = fDataHist->GetXaxis()->GetBinLowEdge(xbin2+1);
- double ylo2 = fDataHist->GetYaxis()->GetBinLowEdge(ybin2);
- double yhi2 = fDataHist->GetYaxis()->GetBinLowEdge(ybin2+1);
-
- if (xlo2 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
- ylo2 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
- xhi2 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) ||
- yhi2 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue;
-
-
- //std::cout << "Correlating: (" << xlo1 << "-" << xhi1 << "," << ylo1 << "-" << yhi1 << ") with (" << xlo2 << "-" << xhi2 << "," << ylo2 << "-" << yhi2 << ")" << std::endl;
-
- // Get the data
- double data2 = fDataHist->GetBinContent(j);
- // Get the MC
- double mc2 = fMCHist->GetBinContent(j);
- //std::cout << data1 << " " << mc1 << std::endl;
- //std::cout << data2 << " " << mc2 << std::endl;
- //std::cout << std::endl;
-
- // Get the inverse covariance matrix entry
- double coventry = covar_th2d->GetBinContent(i, j);
-
- //std::cout << fDataHist->GetXaxis()->GetBinLowEdge(i+1) << " - " << fDataHist->GetXaxis()->GetBinLowEdge(i+2) << ", " << fDataHist->GetYaxis()->GetBinLowEdge(j+1) << " - " << fDataHist->GetYaxis()->GetBinLowEdge(j+2) << " = " << coventry << " (global = " << global << ")" << std::endl;
-
- chi2 += (data1-mc1)*coventry*(data2-mc2);
+ // Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA used for their measurement
+ // Can be prettified in due time but for now keep
+ bool chi2_use_overflow_err = false;
+ const int lowBin = chi2_use_overflow_err?0:1; // Either they both use underflow, or neither of them does
+
+ int nbinsx=fMCHist->GetNbinsX();
+ int nbinsy=fMCHist->GetNbinsY();
+ const int highBinX = nbinsx;
+ const int highBinY = nbinsy;
+ Int_t Nbins = (highBinX-lowBin+1)* (highBinY-lowBin+1); // from low to high inclusive in each dimension
+
+ TMatrixD covMatrixTmp = (*fFullCovar);
+
+ TMatrixD covMatrix(Nbins, Nbins);
+ const double scaling = 1e80;//ROOT can't seem to handle small entries with lots of zeros? Suggested scaling the histogram and then rescaling the inverted matrix
+ for( int i = 0; i != Nbins; ++i ) {
+ for( int j = 0; j != Nbins; ++j ) {
+ // I think this is right... it checks out on a small sample
+ int old_i_bin = (i/nbinsx + 1)* (nbinsx +2) +1 + i%nbinsx;
+ int old_j_bin = (j/nbinsx + 1)* (nbinsx +2) +1 + j%nbinsx;
+ covMatrix[i][j] = covMatrixTmp[old_i_bin][old_j_bin]*scaling;
+ //std::cout << Nbins*Nbins << " (" << Nbins << "*" << Nbins << ")" << std::endl;
+ //std::cout << i << ", " << j << " = " << old_i_bin << " " << old_j_bin << std::endl;
}
+ //throw;
+ }
+ TDecompSVD error(covMatrix);
+ TMatrixD errorMatrix(covMatrix);
+ if( ! error.Invert( errorMatrix ) ) {
+ std::cout << "Cannot invert total covariance matrix. You could use statistical errors only for Chi2 calculation. But it isn't implemented yet" << std::endl;
+ }
+
+ covMatrix *= 1/scaling;
+ errorMatrix *= scaling;
+
+ for( int i = 0; i != Nbins; ++i ) {
+ int hist_i_bin = chi2_use_overflow_err?i:((i/nbinsx + 1)* (nbinsx +2) +1 + i%nbinsx); // Translate to the histogram bin, if we aren't using the overflow errors, meaning the covariance matrix is smaller than the histogram
+ const Double_t x_data_i = fDataHist->GetBinContent(hist_i_bin);
+ const Double_t x_mc_i = fMCHist->GetBinContent(hist_i_bin);
+ for( int j = 0; j != Nbins; ++j ) {
+ // Each element of the inverted covariance matrix corresponds to a pair of data and MC
+ int hist_j_bin = chi2_use_overflow_err?j:((j/nbinsx + 1)* (nbinsx +2) +1 + j%nbinsx);
+ const Double_t x_data_j = fDataHist->GetBinContent(hist_j_bin);
+ const Double_t x_mc_j = fMCHist->GetBinContent(hist_j_bin);
+ const double chi2_ij = (x_data_i - x_mc_i) * errorMatrix[i][j] * (x_data_j - x_mc_j);
+ //std::cout << x_data_i << "\t" << x_mc_i << "\t" << i << "\t" << errorMatrix[i][j] << "\t" << j << "\t" << x_data_j << "\t" << x_mc_j << "\t" << chi2_ij << std::endl;
+
+ chi2 += chi2_ij;
}
- */
+ }
+ //std::cout << Names[a] << " chi2 " << chi2 << std::endl;
+
+ /*
+ * CWRET CALC
+ // Calculate the test-statistic
+ for (int i = 0; i < covar_th2d->GetXaxis()->GetNbins()+1; ++i) {
+ // Get the global bin for x
+ int xbin1, ybin1, zbin1;
+ fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1);
+ double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1);
+ double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1+1);
+ double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1);
+ double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1+1);
+ if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
+ ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
+ xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) ||
+ yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue;
+
+ // Get the data
+ double data1 = fDataHist->GetBinContent(i);
+ // Get the MC
+ double mc1 = fMCHist->GetBinContent(i);
+
+ for (int j = 0; j < covar_th2d->GetYaxis()->GetNbins()+1; ++j) {
+
+ int xbin2, ybin2, zbin2;
+ fDataHist->GetBinXYZ(j, xbin2, ybin2, zbin2);
+ double xlo2 = fDataHist->GetXaxis()->GetBinLowEdge(xbin2);
+ double xhi2 = fDataHist->GetXaxis()->GetBinLowEdge(xbin2+1);
+ double ylo2 = fDataHist->GetYaxis()->GetBinLowEdge(ybin2);
+ double yhi2 = fDataHist->GetYaxis()->GetBinLowEdge(ybin2+1);
+
+ if (xlo2 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
+ ylo2 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
+ xhi2 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) ||
+ yhi2 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue;
+
+
+ //std::cout << "Correlating: (" << xlo1 << "-" << xhi1 << "," << ylo1 << "-" << yhi1 << ") with (" << xlo2 << "-" << xhi2 << "," << ylo2 << "-" << yhi2 << ")" << std::endl;
+
+ // Get the data
+ double data2 = fDataHist->GetBinContent(j);
+ // Get the MC
+ double mc2 = fMCHist->GetBinContent(j);
+ //std::cout << data1 << " " << mc1 << std::endl;
+ //std::cout << data2 << " " << mc2 << std::endl;
+ //std::cout << std::endl;
+
+ // Get the inverse covariance matrix entry
+ double coventry = covar_th2d->GetBinContent(i, j);
+
+ //std::cout << fDataHist->GetXaxis()->GetBinLowEdge(i+1) << " - " << fDataHist->GetXaxis()->GetBinLowEdge(i+2) << ", " << fDataHist->GetYaxis()->GetBinLowEdge(j+1) << " - " << fDataHist->GetYaxis()->GetBinLowEdge(j+2) << " = " << coventry << " (global = " << global << ")" << std::endl;
+
+ chi2 += (data1-mc1)*coventry*(data2-mc2);
+ }
+ }
+ */
//}
// Normalisation penalty term if included
if (fAddNormPen) {
chi2 +=
(1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError);
LOG(REC) << "Norm penalty = "
<< (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) /
(fNormError * fNormError)
<< std::endl;
}
// Adjust the shape back to where it was.
if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
fLikelihood = chi2;
return chi2;
};
diff --git a/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx b/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx
index fc41380..e95e633 100755
--- a/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx
+++ b/src/MINERvA/MINERvA_CCinc_XSec_2DEavq3_nu.cxx
@@ -1,153 +1,153 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCinc_XSec_2DEavq3_nu.h"
//********************************************************************
MINERvA_CCinc_XSec_2DEavq3_nu::MINERvA_CCinc_XSec_2DEavq3_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCinc_XSec_2DEavq3_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Medium Energy FHC numu \n" \
"Signal: CC-inclusive with theta < 20deg \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("q_{3} (GeV)");
fSettings.SetYTitle("E_{avail} (GeV)");
fSettings.SetZTitle("d^{2}#sigma/dq_{3}dE_{avail} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
fSettings.SetEnuRange(2.0, 6.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CCinc_XSec_2DEavq3_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/data_2D.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/covar_2D.txt" );
fSettings.SetMapInput( FitPar::GetDataBase() + "/MINERvA/CCEavq3/map_2D.txt" );
fSettings.DefineAllowedSpecies("numu");
hadroncut = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu_hadron_cut");
useq3true = FitPar::Config().GetParB("MINERvA_CCinc_XSec_2DEavq3_nu_useq3true");
splitMEC_PN_NN = FitPar::Config().GetParB("Modes_split_PN_NN");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-42 / (fNEvents + 0.)) / this->TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
Double_t binx[7] = {0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8};
Double_t biny[17] = {0.0, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.25, 0.30, 0.35, 0.40, 0.50, 0.60, 0.80};
CreateDataHistogram(7, binx, 17, biny);
SetDataValuesFromTextFile( fSettings.GetDataInput() );
ScaleData(1E-42);
SetMapValuesFromText( fSettings.GetMapInput() );
SetCholDecompFromTextFile( fSettings.GetCovarInput(), 67);
ScaleCovar(1E-16);
- StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38);
+ StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38, false);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CCinc_XSec_2DEavq3_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Seperate MEC
if (splitMEC_PN_NN) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < event->Npart(); j++) {
if ((event->PartInfo(j))->fIsAlive) continue;
if (event->PartInfo(j)->fPID == 2212) npr++;
else if (event->PartInfo(j)->fPID == 2112) nne++;
}
if (event->Mode == 2 and npr == 1 and nne == 1) {
event->Mode = 2;
Mode = 2;
} else if (event->Mode == 2 and npr == 0 and nne == 2) {
event->Mode = 3;
Mode = 3;
}
}
// Set Defaults
double Eav = -999.9;
double q3 = -999.9;
// If muon found get kinematics
FitParticle* muon = event->GetHMFSParticle(13);
FitParticle* neutrino = event->GetNeutrinoIn();
if (muon && neutrino) {
// Set Q from Muon
TLorentzVector q = neutrino->fP - muon->fP;
double q0 = (q.E()) / 1.E3;
//double q3_true = (q.Vect().Mag())/1.E3;
double thmu = muon->fP.Vect().Angle(neutrino->fP.Vect());
double pmu = muon->fP.Vect().Mag() / 1.E3;
double emu = muon->fP.E() / 1.E3;
double mmu = muon->fP.Mag() / 1.E3;
// Get Enu Rec
double enu_rec = emu + q0;
// Set Q2 QE
double q2qe = 2 * enu_rec * (emu - pmu * cos(thmu)) - mmu * mmu;
// Calc Q3 from Q2QE and EnuTree
q3 = sqrt(q2qe + q0 * q0);
// Get Eav too
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event) / 1.E3;
}
// Set Hist Variables
fXVar = q3;
fYVar = Eav;
return;
}
//********************************************************************
bool MINERvA_CCinc_XSec_2DEavq3_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCCincLowRecoil_MINERvA(event, EnuMin, EnuMax);
}
diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx
index a787cea..c4a173c 100644
--- a/src/Statistical/StatUtils.cxx
+++ b/src/Statistical/StatUtils.cxx
@@ -1,1358 +1,1362 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "StatUtils.h"
#include "GeneralUtils.h"
#include "NuisConfig.h"
#include "TH1D.h"
//*******************************************************************
Double_t StatUtils::GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask) {
//*******************************************************************
Double_t Chi2 = 0.0;
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// Add MC Error to data if required
if (FitPar::Config().GetParB("addmcerror")) {
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double dterr = calc_data->GetBinError(i + 1);
double mcerr = calc_mc->GetBinError(i + 1);
if (dterr > 0.0) {
calc_data->SetBinError(i + 1, sqrt(dterr * dterr + mcerr * mcerr));
}
}
}
// Apply masking if required
if (mask) {
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Iterate over bins in X
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
// Ignore bins with zero data or zero bin error
if (calc_data->GetBinError(i + 1) <= 0.0 ||
calc_data->GetBinContent(i + 1) == 0.0)
continue;
// Take mc data difference
double diff =
calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1);
double err = calc_data->GetBinError(i + 1);
Chi2 += (diff * diff) / (err * err);
}
// cleanup
delete calc_data;
delete calc_mc;
return Chi2;
};
//*******************************************************************
Double_t StatUtils::GetChi2FromDiag(TH2D* data, TH2D* mc, TH2I* map,
TH2I* mask) {
//*******************************************************************
// Generate a simple map
if (!map) map = GenerateMap(data);
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate 1D chi2 from 1D Plots
Double_t Chi2 = StatUtils::GetChi2FromDiag(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
};
//*******************************************************************
Double_t StatUtils::GetChi2FromCov(TH1D* data, TH1D* mc, TMatrixDSym* invcov,
TH1I* mask, double data_scale,
double covar_scale) {
//*******************************************************************
Double_t Chi2 = 0.0;
TMatrixDSym* calc_cov = (TMatrixDSym*)invcov->Clone();
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// If a mask if applied we need to apply it before the matrix is inverted
if (mask) {
calc_cov = ApplyInvertedMatrixMasking(invcov, mask);
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Add MC Error to data if required
if (FitPar::Config().GetParB("statutils.addmcerror")) {
// Make temp cov
TMatrixDSym* newcov = StatUtils::GetInvert(calc_cov);
// Add MC err to diag
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double mcerr = calc_mc->GetBinError(i + 1) * sqrt(covar_scale);
double oldval = (*newcov)(i, i);
LOG(FIT) << "Adding cov stat " << mcerr * mcerr << " to "
<< (*newcov)(i, i) << std::endl;
(*newcov)(i, i) = oldval + mcerr * mcerr;
}
// Reset the calc_cov to new invert
delete calc_cov;
calc_cov = GetInvert(newcov);
// Delete the tempcov
delete newcov;
}
calc_data->Scale(data_scale);
calc_mc->Scale(data_scale);
(*calc_cov) *= covar_scale;
// iterate over bins in X (i,j)
QLOG(DEB, "START Chi2 Calculation=================");
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
QLOG(DEB,
"[CHI2] i = " << i << " ["
<< calc_data->GetXaxis()->GetBinLowEdge(i + 1) << " -- "
<< calc_data->GetXaxis()->GetBinUpEdge(i + 1) << "].");
for (int j = 0; j < calc_data->GetNbinsX(); j++) {
QLOG(DEB, "[CHI2]\t j = "
<< i << " [" << calc_data->GetXaxis()->GetBinLowEdge(j + 1)
<< " -- " << calc_data->GetXaxis()->GetBinUpEdge(j + 1)
<< "].");
if ((calc_data->GetBinContent(i + 1) != 0 ||
calc_mc->GetBinContent(i + 1) != 0) &&
((*calc_cov)(i, j) != 0)) {
QLOG(DEB, "[CHI2]\t\t Chi2 contribution (i,j) = (" << i << "," << j
<< ")");
QLOG(DEB, "[CHI2]\t\t Data - MC(i) = "
<< calc_data->GetBinContent(i + 1) << " - "
<< calc_mc->GetBinContent(i + 1) << " = "
<< (calc_data->GetBinContent(i + 1) -
calc_mc->GetBinContent(i + 1)));
QLOG(DEB, "[CHI2]\t\t Data - MC(j) = "
<< calc_data->GetBinContent(j + 1) << " - "
<< calc_mc->GetBinContent(j + 1) << " = "
<< (calc_data->GetBinContent(j + 1) -
calc_mc->GetBinContent(j + 1)));
QLOG(DEB, "[CHI2]\t\t Covar = " << (*calc_cov)(i, j));
QLOG(DEB, "[CHI2]\t\t Cont chi2 = "
<< ((calc_data->GetBinContent(i + 1) -
calc_mc->GetBinContent(i + 1)) *
(*calc_cov)(i, j) * (calc_data->GetBinContent(j + 1) -
calc_mc->GetBinContent(j + 1)))
<< " " << Chi2);
Chi2 +=
((calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1)) *
(*calc_cov)(i, j) *
(calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)));
} else {
QLOG(DEB, "Skipping chi2 contribution (i,j) = ("
<< i << "," << j
<< "), Data = " << calc_data->GetBinContent(i + 1)
<< ", MC = " << calc_mc->GetBinContent(i + 1)
<< ", Cov = " << (*calc_cov)(i, j));
Chi2 += 0.;
}
}
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromCov(TH2D* data, TH2D* mc, TMatrixDSym* invcov,
TH2I* map, TH2I* mask) {
//*******************************************************************
// Generate a simple map
if (!map) {
map = StatUtils::GenerateMap(data);
}
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate 1D chi2 from 1D Plots
Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD(TH1D* data, TH1D* mc, TMatrixDSym* cov,
TH1I* mask) {
//*******************************************************************
Double_t Chi2 = 0.0;
TMatrixDSym* calc_cov = (TMatrixDSym*)cov->Clone();
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// If a mask if applied we need to apply it before the matrix is inverted
if (mask) {
calc_cov = StatUtils::ApplyMatrixMasking(cov, mask);
calc_data = StatUtils::ApplyHistogramMasking(data, mask);
calc_mc = StatUtils::ApplyHistogramMasking(mc, mask);
}
// Decompose matrix
TDecompSVD LU = TDecompSVD((*calc_cov));
LU.Decompose();
TMatrixDSym* cov_U =
new TMatrixDSym(calc_data->GetNbinsX(), LU.GetU().GetMatrixArray(), "");
TVectorD* cov_S = new TVectorD(LU.GetSig());
// Apply basis rotation before adding up chi2
Double_t rotated_difference = 0.0;
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
rotated_difference = 0.0;
// Rotate basis of Data - MC
for (int j = 0; j < calc_data->GetNbinsY(); j++)
rotated_difference +=
(calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)) *
(*cov_U)(j, i);
// Divide by rotated error cov_S
Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i);
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
delete cov_U;
delete cov_S;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD(TH2D* data, TH2D* mc, TMatrixDSym* cov,
TH2I* map, TH2I* mask) {
//*******************************************************************
// Generate a simple map
if (!map) map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
double StatUtils::GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask) {
//*******************************************************************
// If just an event rate, for chi2 just use Poission Likelihood to calculate
// the chi2 component
double chi2 = 0.0;
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// Apply masking if required
if (mask) {
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Iterate over bins in X
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double dt = calc_data->GetBinContent(i + 1);
double mc = calc_mc->GetBinContent(i + 1);
if (mc <= 0) continue;
if (dt <= 0) {
// Only add difference
chi2 += 2 * (mc - dt);
} else {
// Do the chi2 for Poisson distributions
chi2 += 2 * (mc - dt + (dt * log(dt / mc)));
}
/*
LOG(REC)<<"Evt Chi2 cont = "<Clone();
// If a mask is provided we need to apply it before getting NDOF
if (mask) {
calc_hist = StatUtils::ApplyHistogramMasking(hist, mask);
}
// NDOF is defined as total number of bins with non-zero errors
Int_t NDOF = 0;
for (int i = 0; i < calc_hist->GetNbinsX(); i++) {
if (calc_hist->GetBinError(i + 1) > 0.0) NDOF++;
}
delete calc_hist;
return NDOF;
};
//*******************************************************************
Int_t StatUtils::GetNDOF(TH2D* hist, TH2I* map, TH2I* mask) {
//*******************************************************************
Int_t NDOF = 0;
if (!map) map = StatUtils::GenerateMap(hist);
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1)) continue;
if (map->GetBinContent(i + 1, j + 1) <= 0) continue;
NDOF++;
}
}
return NDOF;
};
//*******************************************************************
TH1D* StatUtils::ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag,
TH1I* mask) {
//*******************************************************************
TH1D* calc_hist =
(TH1D*)hist->Clone((std::string(hist->GetName()) + "_THROW").c_str());
TMatrixDSym* calc_cov = (TMatrixDSym*)cov->Clone();
Double_t correl_val = 0.0;
// If a mask if applied we need to apply it before the matrix is decomposed
if (mask) {
calc_cov = ApplyMatrixMasking(cov, mask);
calc_hist = ApplyHistogramMasking(calc_hist, mask);
}
// If a covariance is provided we need a preset random vector and a decomp
std::vector rand_val;
TMatrixDSym* decomp_cov;
if (cov) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
rand_val.push_back(gRandom->Gaus(0.0, 1.0));
}
// Decomp the matrix
decomp_cov = StatUtils::GetDecomp(calc_cov);
}
// iterate over bins
for (int i = 0; i < hist->GetNbinsX(); i++) {
// By Default the errors on the histogram are thrown uncorrelated to the other errors
/*
if (throwdiag) {
calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + \
gRandom->Gaus(0.0, 1.0) * calc_hist->GetBinError(i + 1)) );
}
*/
// If a covariance is provided that is also thrown
if (cov) {
correl_val = 0.0;
for (int j = 0; j < hist->GetNbinsX(); j++) {
correl_val += rand_val[j] * (*decomp_cov)(j, i);
}
calc_hist->SetBinContent(
i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38));
}
}
delete calc_cov;
delete decomp_cov;
// return this new thrown data
return calc_hist;
};
//*******************************************************************
TH2D* StatUtils::ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map,
bool throwdiag, TH2I* mask) {
//*******************************************************************
// PLACEHOLDER!!!!!!!!!
// Currently no support for throwing 2D Histograms from a covariance
(void)hist;
(void)cov;
(void)map;
(void)throwdiag;
(void)mask;
// /todo
// Sort maps if required
// Throw the covariance for a 1D plot
// Unmap back to 2D Histogram
return hist;
}
//*******************************************************************
TH1D* StatUtils::ApplyHistogramMasking(TH1D* hist, TH1I* mask) {
//*******************************************************************
if (!mask) return ((TH1D*)hist->Clone());
// This masking is only sufficient for chi2 calculations, and will have dodgy
// bin edges.
// Get New Bin Count
Int_t NBins = 0;
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1)) continue;
NBins++;
}
// Make new hist
std::string newmaskname = std::string(hist->GetName()) + "_MSKD";
TH1D* calc_hist =
new TH1D(newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins);
// fill new hist
int binindex = 0;
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1)) {
LOG(REC) << "Applying mask to bin " << i + 1 << " " << hist->GetName()
<< std::endl;
continue;
}
calc_hist->SetBinContent(binindex + 1, hist->GetBinContent(i + 1));
calc_hist->SetBinError(binindex + 1, hist->GetBinError(i + 1));
binindex++;
}
return calc_hist;
};
//*******************************************************************
TH2D* StatUtils::ApplyHistogramMasking(TH2D* hist, TH2I* mask) {
//*******************************************************************
TH2D* newhist = (TH2D*)hist->Clone();
if (!mask) return newhist;
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1) > 0) {
newhist->SetBinContent(i + 1, j + 1, 0.0);
newhist->SetBinContent(i + 1, j + 1, 0.0);
}
}
}
return newhist;
}
//*******************************************************************
TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH1I* mask) {
//*******************************************************************
if (!mask) return (TMatrixDSym*)(mat->Clone());
// Get New Bin Count
Int_t NBins = 0;
for (int i = 0; i < mask->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1)) continue;
NBins++;
}
// make new matrix
TMatrixDSym* calc_mat = new TMatrixDSym(NBins);
int col, row;
// Need to mask out bins in the current matrix
row = 0;
for (int i = 0; i < mask->GetNbinsX(); i++) {
col = 0;
// skip if masked
if (mask->GetBinContent(i + 1) > 0.5) continue;
for (int j = 0; j < mask->GetNbinsX(); j++) {
// skip if masked
if (mask->GetBinContent(j + 1) > 0.5) continue;
(*calc_mat)(row, col) = (*mat)(i, j);
col++;
}
row++;
}
return calc_mat;
};
//*******************************************************************
TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data,
TH2I* mask, TH2I* map) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1I* mask_1D = StatUtils::MapToMask(mask, map);
TMatrixDSym* newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D);
delete mask_1D;
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat,
TH1I* mask) {
//*******************************************************************
TMatrixDSym* new_mat = GetInvert(mat);
TMatrixDSym* masked_mat = ApplyMatrixMasking(new_mat, mask);
TMatrixDSym* inverted_mat = GetInvert(masked_mat);
delete masked_mat;
delete new_mat;
return inverted_mat;
};
//*******************************************************************
TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data,
TH2I* mask, TH2I* map) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1I* mask_1D = StatUtils::MapToMask(mask, map);
TMatrixDSym* newmat = ApplyInvertedMatrixMasking(mat, mask_1D);
delete mask_1D;
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetInvert(TMatrixDSym* mat) {
//*******************************************************************
TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone();
// Check for diagonal
bool non_diagonal = false;
for (int i = 0; i < new_mat->GetNrows(); i++) {
for (int j = 0; j < new_mat->GetNrows(); j++) {
if (i == j) continue;
if ((*new_mat)(i, j) != 0.0) {
non_diagonal = true;
break;
}
}
}
// If diag, just flip the diag
if (!non_diagonal or new_mat->GetNrows() == 1) {
for (int i = 0; i < new_mat->GetNrows(); i++) {
if ((*new_mat)(i, i) != 0.0)
(*new_mat)(i, i) = 1.0 / (*new_mat)(i, i);
else
(*new_mat)(i, i) = 0.0;
}
return new_mat;
}
// Invert full matrix
TDecompSVD LU = TDecompSVD((*new_mat));
new_mat =
new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), "");
return new_mat;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetDecomp(TMatrixDSym* mat) {
//*******************************************************************
TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone();
int nrows = new_mat->GetNrows();
// Check for diagonal
bool diagonal = true;
for (int i = 0; i < nrows; i++) {
for (int j = 0; j < nrows; j++) {
if (i == j) continue;
if ((*new_mat)(i, j) != 0.0) {
diagonal = false;
break;
}
}
}
// If diag, just flip the diag
if (diagonal or nrows == 1) {
for (int i = 0; i < nrows; i++) {
if ((*new_mat)(i, i) > 0.0)
(*new_mat)(i, i) = sqrt((*new_mat)(i, i));
else
(*new_mat)(i, i) = 0.0;
}
return new_mat;
}
TDecompChol LU = TDecompChol(*new_mat);
LU.Decompose();
delete new_mat;
TMatrixDSym* dec_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), "");
return dec_mat;
}
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH1D* hist, double norm) {
//*******************************************************************
if (!mat) mat = MakeDiagonalCovarMatrix(hist);
int nbins = mat->GetNrows();
TMatrixDSym* new_mat = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
double valx = hist->GetBinContent(i + 1) * 1E38;
double valy = hist->GetBinContent(j + 1) * 1E38;
(*new_mat)(i, j) = (*mat)(i, j) + norm * norm * valx * valy;
}
}
// Swap the two
delete mat;
mat = new_mat;
return;
};
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm,
TH2I* map) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1D* data_1D = MapToTH1D(data, map);
StatUtils::ForceNormIntoCovar(mat, data_1D, norm);
delete data_1D;
return;
}
//*******************************************************************
TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH1D* data, double scaleF) {
//*******************************************************************
TMatrixDSym* newmat = new TMatrixDSym(data->GetNbinsX());
for (int i = 0; i < data->GetNbinsX(); i++) {
(*newmat)(i, i) =
data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF;
}
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH2D* data, TH2I* map,
double scaleF) {
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1D* data_1D = MapToTH1D(data, map);
return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF);
};
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov,
double scale, bool ErrorCheck) {
//*******************************************************************
// Check
- if (cov->GetNrows() != data->GetNbinsX()) {
- ERR(FTL) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << std::endl;
- ERR(FTL) << "Nrows = " << cov->GetNrows() << std::endl;
- ERR(FTL) << "Nbins = " << data->GetNbinsX() << std::endl;
- throw;
+ if (ErrorCheck) {
+ if (cov->GetNrows() != data->GetNbinsX()) {
+ ERR(FTL) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << std::endl;
+ ERR(FTL) << "Nrows = " << cov->GetNrows() << std::endl;
+ ERR(FTL) << "Nbins = " << data->GetNbinsX() << std::endl;
+ throw;
+ }
}
// Set bin errors form cov diag
// Check if the errors are set
bool ErrorsSet = false;
for (int i = 0; i < data->GetNbinsX(); i++) {
if (ErrorsSet == true) break;
if (data->GetBinError(i+1) != 0) ErrorsSet = true;
}
// Now loop over
if (ErrorsSet && ErrorCheck) {
for (int i = 0; i < data->GetNbinsX(); i++) {
double dataerr = data->GetBinError(i + 1);
double coverr = sqrt((*cov)(i, i))*scale;
// Check that the errors are within 1% of eachother
if (fabs(dataerr-coverr)/dataerr > 0.01) {
ERR(FTL) << "Data error does not match covariance error for bin " << i+1 << " (" << data->GetXaxis()->GetBinLowEdge(i+1) << "-" << data->GetXaxis()->GetBinLowEdge(i+2) << ")" << std::endl;
ERR(FTL) << "Data error: " << dataerr << std::endl;
ERR(FTL) << "Cov error: " << coverr << std::endl;
}
}
// Else blindly trust the covariance
} else {
for (int i = 0; i < data->GetNbinsX(); i++) {
data->SetBinError(i+1, sqrt((*cov)(i,i))*scale);
}
}
return;
}
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map, double scale, bool ErrorCheck) {
//*******************************************************************
// Check
- if (cov->GetNrows() != data->GetNbinsX()*data->GetNbinsY()) {
- ERR(FTL) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << std::endl;
- ERR(FTL) << "Nrows = " << cov->GetNrows() << std::endl;
- ERR(FTL) << "Nbins = " << data->GetNbinsX() << std::endl;
- throw;
+ if (ErrorCheck) {
+ if (cov->GetNrows() != data->GetNbinsX()*data->GetNbinsY()) {
+ ERR(FTL) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << std::endl;
+ ERR(FTL) << "Nrows = " << cov->GetNrows() << std::endl;
+ ERR(FTL) << "Nbins = " << data->GetNbinsX() << std::endl;
+ throw;
+ }
}
// Set bin errors form cov diag
// Check if the errors are set
bool ErrorsSet = false;
for (int i = 0; i < data->GetNbinsX(); i++) {
for (int j = 0; j < data->GetNbinsX(); j++) {
if (ErrorsSet == true) break;
if (data->GetBinError(i+1, j+1) != 0) ErrorsSet = true;
}
}
// Create map if required
if (!map) map = StatUtils::GenerateMap(data);
// Set Bin Errors from cov diag
int count = 0;
for (int i = 0; i < data->GetNbinsX(); i++) {
for (int j = 0; j < data->GetNbinsY(); j++) {
if (data->GetBinContent(i + 1, j + 1) == 0.0) continue;
// Get the entry in the cov matrix
count = map->GetBinContent(i + 1, j + 1) - 1;
double dataerr = data->GetBinError(i+1, j+1);
double coverr = sqrt((*cov)(count,count))*scale;
// Check that the errors are within 1% of eachother
if (ErrorsSet && ErrorCheck) {
if (fabs(dataerr-coverr)/dataerr > 0.01) {
ERR(FTL) << "Data error does not match covariance error for bin " << i+1 << " (" << data->GetXaxis()->GetBinLowEdge(i+1) << "-" << data->GetXaxis()->GetBinLowEdge(i+2) << ")" << std::endl;
ERR(FTL) << "Data error: " << dataerr << std::endl;
ERR(FTL) << "Cov error: " << coverr << std::endl;
}
} else {
data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale);
}
}
}
return;
}
TMatrixDSym* StatUtils::ExtractShapeOnlyCovar(TMatrixDSym* full_covar,
TH1* data_hist,
double data_scale) {
int nbins = full_covar->GetNrows();
TMatrixDSym* shape_covar = new TMatrixDSym(nbins);
// Check nobody is being silly
if (data_hist->GetNbinsX() != nbins) {
ERR(WRN) << "Inconsistent matrix and data histogram passed to "
"StatUtils::ExtractShapeOnlyCovar!"
<< std::endl;
ERR(WRN) << "data_hist has " << data_hist->GetNbinsX() << " matrix has "
<< nbins << std::endl;
int err_bins = data_hist->GetNbinsX();
if (nbins > err_bins) err_bins = nbins;
for (int i = 0; i < err_bins; ++i) {
ERR(WRN) << "Matrix diag. = " << (*full_covar)(i, i)
<< " data = " << data_hist->GetBinContent(i + 1) << std::endl;
}
return NULL;
}
double total_data = 0;
double total_covar = 0;
// Initial loop to calculate some constants
for (int i = 0; i < nbins; ++i) {
total_data += data_hist->GetBinContent(i + 1) * data_scale;
for (int j = 0; j < nbins; ++j) {
total_covar += (*full_covar)(i, j);
}
}
if (total_data == 0 || total_covar == 0) {
ERR(WRN) << "Stupid matrix or data histogram passed to "
"StatUtils::ExtractShapeOnlyCovar! Ignoring..."
<< std::endl;
return NULL;
}
LOG(SAM) << "Norm error = " << sqrt(total_covar) / total_data << std::endl;
// Now loop over and calculate the shape-only matrix
for (int i = 0; i < nbins; ++i) {
double data_i = data_hist->GetBinContent(i + 1) * data_scale;
for (int j = 0; j < nbins; ++j) {
double data_j = data_hist->GetBinContent(j + 1) * data_scale;
double norm_term =
data_i * data_j * total_covar / total_data / total_data;
double mix_sum1 = 0;
double mix_sum2 = 0;
for (int k = 0; k < nbins; ++k) {
mix_sum1 += (*full_covar)(k, j);
mix_sum2 += (*full_covar)(i, k);
}
double mix_term1 =
data_i * (mix_sum1 / total_data -
total_covar * data_j / total_data / total_data);
double mix_term2 =
data_j * (mix_sum2 / total_data -
total_covar * data_i / total_data / total_data);
(*shape_covar)(i, j) =
(*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term;
}
}
return shape_covar;
}
//*******************************************************************
TH2I* StatUtils::GenerateMap(TH2D* hist) {
//*******************************************************************
std::string maptitle = std::string(hist->GetName()) + "_MAP";
TH2I* map =
new TH2I(maptitle.c_str(), maptitle.c_str(), hist->GetNbinsX(), 0,
hist->GetNbinsX(), hist->GetNbinsY(), 0, hist->GetNbinsY());
Int_t index = 1;
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (hist->GetBinContent(i + 1, j + 1) > 0 &&
hist->GetBinError(i + 1, j + 1) > 0) {
map->SetBinContent(i + 1, j + 1, index);
index++;
} else {
map->SetBinContent(i + 1, j + 1, 0);
}
}
}
return map;
}
//*******************************************************************
TH1D* StatUtils::MapToTH1D(TH2D* hist, TH2I* map) {
//*******************************************************************
if (!hist) return NULL;
// Get N bins for 1D plot
Int_t Nbins = map->GetMaximum();
std::string name1D = std::string(hist->GetName()) + "_1D";
// Make new 1D Hist
TH1D* newhist = new TH1D(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);
// map bin contents
for (int i = 0; i < map->GetNbinsX(); i++) {
for (int j = 0; j < map->GetNbinsY(); j++) {
if (map->GetBinContent(i + 1, j + 1) == 0) continue;
newhist->SetBinContent(map->GetBinContent(i + 1, j + 1),
hist->GetBinContent(i + 1, j + 1));
newhist->SetBinError(map->GetBinContent(i + 1, j + 1),
hist->GetBinError(i + 1, j + 1));
}
}
// return
return newhist;
}
//*******************************************************************
TH1I* StatUtils::MapToMask(TH2I* hist, TH2I* map) {
//*******************************************************************
TH1I* newhist = NULL;
if (!hist) return newhist;
// Get N bins for 1D plot
Int_t Nbins = map->GetMaximum();
std::string name1D = std::string(hist->GetName()) + "_1D";
// Make new 1D Hist
newhist = new TH1I(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);
// map bin contents
for (int i = 0; i < map->GetNbinsX(); i++) {
for (int j = 0; j < map->GetNbinsY(); j++) {
if (map->GetBinContent(i + 1, j + 1) == 0) continue;
newhist->SetBinContent(map->GetBinContent(i + 1, j + 1),
hist->GetBinContent(i + 1, j + 1));
}
}
// return
return newhist;
}
TMatrixDSym* StatUtils::GetCovarFromCorrel(TMatrixDSym* correl, TH1D* data) {
int nbins = correl->GetNrows();
TMatrixDSym* covar = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*covar)(i, j) =
(*correl)(i, j) * data->GetBinError(i + 1) * data->GetBinError(j + 1);
}
}
return covar;
}
//*******************************************************************
TMatrixD* StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx,
int dimy) {
//*******************************************************************
// Determine dim
if (dimx == -1 and dimy == -1) {
std::string line;
std::ifstream covar(covfile.c_str(), std::ifstream::in);
int row = 0;
while (std::getline(covar >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
"entries on this line: "
<< row << std::endl;
}
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
column++;
if (column > dimx) dimx = column;
}
row++;
if (row > dimy) dimy = row;
}
}
// Or assume symmetric
if (dimx != -1 and dimy == -1) {
dimy = dimx;
}
assert(dimy != -1 && " matrix dimy not set.");
// Make new matrix
TMatrixD* mat = new TMatrixD(dimx, dimy);
std::string line;
std::ifstream covar(covfile.c_str(), std::ifstream::in);
int row = 0;
while (std::getline(covar >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
"entries on this line: "
<< row << std::endl;
}
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
// Check Rows
// assert(row > mat->GetNrows() && " covar rows doesn't match matrix
// rows.");
// assert(column > mat->GetNcols() && " covar cols doesn't match matrix
// cols.");
// Fill Matrix
(*mat)(row, column) = (*iter);
column++;
}
row++;
}
return mat;
}
//*******************************************************************
TMatrixD* StatUtils::GetMatrixFromRootFile(std::string covfile,
std::string histname) {
//*******************************************************************
std::string inputfile = covfile + ";" + histname;
std::vector splitfile = GeneralUtils::ParseToStr(inputfile, ";");
if (splitfile.size() < 2) {
ERR(FTL) << "No object name given!" << std::endl;
throw;
}
// Get file
TFile* tempfile = new TFile(splitfile[0].c_str(), "READ");
// Get Object
TObject* obj = tempfile->Get(splitfile[1].c_str());
if (!obj) {
ERR(FTL) << "Object " << splitfile[1] << " doesn't exist!" << std::endl;
throw;
}
// Try casting
TMatrixD* mat = dynamic_cast(obj);
if (mat) {
TMatrixD* newmat = (TMatrixD*)mat->Clone();
delete mat;
tempfile->Close();
return newmat;
}
TMatrixDSym* matsym = dynamic_cast(obj);
if (matsym) {
TMatrixD* newmat = new TMatrixD(matsym->GetNrows(), matsym->GetNrows());
for (int i = 0; i < matsym->GetNrows(); i++) {
for (int j = 0; j < matsym->GetNrows(); j++) {
(*newmat)(i, j) = (*matsym)(i, j);
}
}
delete matsym;
tempfile->Close();
return newmat;
}
TH2D* mathist = dynamic_cast(obj);
if (mathist) {
TMatrixD* newmat = new TMatrixD(mathist->GetNbinsX(), mathist->GetNbinsX());
for (int i = 0; i < mathist->GetNbinsX(); i++) {
for (int j = 0; j < mathist->GetNbinsX(); j++) {
(*newmat)(i, j) = mathist->GetBinContent(i + 1, j + 1);
}
}
delete mathist;
tempfile->Close();
return newmat;
}
return NULL;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetCovarFromTextFile(std::string covfile, int dim) {
//*******************************************************************
// Delete TempMat
TMatrixD* tempmat = GetMatrixFromTextFile(covfile, dim, dim);
// Make a symmetric covariance
TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows());
for (int i = 0; i < tempmat->GetNrows(); i++) {
for (int j = 0; j < tempmat->GetNrows(); j++) {
(*newmat)(i, j) = (*tempmat)(i, j);
}
}
delete tempmat;
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::GetCovarFromRootFile(std::string covfile,
std::string histname) {
//*******************************************************************
TMatrixD* tempmat = GetMatrixFromRootFile(covfile, histname);
TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows());
for (int i = 0; i < tempmat->GetNrows(); i++) {
for (int j = 0; j < tempmat->GetNrows(); j++) {
(*newmat)(i, j) = (*tempmat)(i, j);
}
}
delete tempmat;
return newmat;
}