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; }