Page MenuHomeHEPForge
Files F10664299

No OneTemporary
Actions

View Options

diff --git a/data/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root b/data/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root
index 1bae70d..5449bcb 100644
Binary files a/data/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root and b/data/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root differ
diff --git a/data/T2K/CC0pi/makedatafile_t2kcc0pi.py b/data/T2K/CC0pi/makedatafile_t2kcc0pi.py
index 19c09b6..a8710b7 100644
--- a/data/T2K/CC0pi/makedatafile_t2kcc0pi.py
+++ b/data/T2K/CC0pi/makedatafile_t2kcc0pi.py
@@ -1,148 +1,186 @@
from ROOT import *
from array import *
def GetMiddle(mystr):
lims = mystr.strip().split(" - ")
val = (float(lims[0]) + float(lims[1]))/2.0
return val
+def GetLowEdge(mystr):
+
+ lims = mystr.strip().split(" - ")
+ val = (float(lims[0]) + 0.00001)
+
+ return val
+
+def GetHighEdge(mystr):
+
+ lims = mystr.strip().split(" - ")
+ val = (float(lims[1]) - 0.00001)
+
+ return val
+
def GetIndex(mystr):
lims = mystr.split("-")
return int(lims[0]), int(lims[1])
outfile = TFile("T2K_CC0PI_2DPmuCosmu_Data.root","RECREATE")
# ANALYSIS I
#______________________________
+xedge = [0.0, 0.3, 0.4, 0.5, 0.65, 0.8, 0.95, 1.1, 1.25, 1.5, 2.0, 3.0, 5.0, 30.0]
+yedge = [-1.0, 0.0, 0.6, 0.7, 0.8, 0.85, 0.9, 0.94, 0.98, 1.0]
+datahist = TH2D("analysis1_data","analysis1_data",
+ len(xedge)-1, array('f',xedge),
+ len(yedge)-1, array('f',yedge))
+
+maphist = datahist.Clone("analysis1_map")
+maphist.SetTitle("analysis1_map")
+
+counthist = datahist.Clone("analysis1_entrycount")
+
with open("cross-section_analysisI.txt") as f:
count = 0
for line in f:
count += 1
if (count < 4): continue
data = line.strip().split("|")
if (len(data) < 1): continue
- ibin = int( data[0] )
- xval = GetMiddle( data[2] )
- yval = GetMiddle( data[1] )
- xsec = float( data[3] ) * 1E-38
+ ibin = int( data[0] ) + 1
- print ibin, xval, yval, xsec
-
+ xval = GetLowEdge( data[2] )
+ yval = GetLowEdge( data[1] )
+ xhig = GetHighEdge( data[2] )
+ yhig = GetHighEdge( data[1] )
+ xsec = float( data[3] ) * 1E-38
+
+ datahist.Fill(xval, yval, xsec)
+ counthist.Fill(xval, yval, 1.0)
+
+ for i in range(maphist.GetNbinsX()):
+ for j in range(maphist.GetNbinsY()):
+ xcent = maphist.GetXaxis().GetBinCenter(i+1)
+ ycent = maphist.GetYaxis().GetBinCenter(j+1)
+
+ if (xcent > xval and xcent < xhig and
+ ycent > yval and ycent < yhig):
+ maphist.SetBinContent(i+1,j+1, ibin)
+
+outfile.cd()
+datahist.Write()
+counthist.Write()
+maphist.Write()
+
# ANALYSIS II
#______________________________
-
xedge = [0.2, 0.35, 0.5, 0.65, 0.8, 0.95, 1.1, 1.25, 1.5, 2.0, 3.0, 5.0, 30.0]
yedge = [0.6, 0.7, 0.8, 0.85, 0.9, 0.925, 0.95, 0.975, 1.0]
datahist = TH2D("analysis2_data","analysis2_data",
len(xedge)-1, array('f',xedge),
len(yedge)-1, array('f',yedge))
maphist = datahist.Clone("analysis2_map")
maphist.SetTitle("analysis2_map")
counthist = datahist.Clone("analysis2_entrycount")
# Get Data Entries
entries = []
count = 0
with open("rps_crossSection_analysis2.txt") as f:
for line in f:
count += 1
if (count < 4): continue
data = line.strip().split("|")
if (len(data) < 1): continue
- ibin = int( data[0] )
+ ibin = int( data[0] ) + 1
xval = GetMiddle( data[2] )
yval = GetMiddle( data[1] )
xsec = float( data[3] ) * 1E-38
- entries.append( [ibin, xval, yval, xsec] )
-
-# Fill data hist
-for vals in entries:
- ibin, xval, yval, xsec = vals
-
- datahist.Fill(xval, yval, xsec)
- maphist.Fill(xval, yval, ibin)
-
- counthist.Fill(xval, yval, 1.0)
+ datahist.Fill(xval, yval, xsec)
+ maphist.Fill(xval, yval, ibin)
+
+ counthist.Fill(xval, yval, 1.0)
+ print ibin, "Map Value"
+
# Get N Bins
nbins = int(maphist.GetMaximum())
print "NBins I = ", nbins
# Get Covariances (keep in 1E-38 cm^2)
statcov = TH2D("analysis2_statcov","analysis2_statcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
systcov = TH2D("analysis2_systcov","analysis2_systcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
normcov = TH2D("analysis2_normcov","analysis2_normcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
totcov = TH2D("analysis2_totcov","analysis2_totcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
with open("rps_statsCov_analysis2.txt") as f:
count = 0
for line in f:
count += 1
if (count < 4): continue
data = line.strip().split("|")
if (len(data) < 1): continue
xi, yi = GetIndex(data[0])
cov = float(data[1])
statcov.SetBinContent(xi, yi, cov)
with open("rps_systCov_analysis2.txt") as f:
count = 0
for line in f:
count += 1
if (count < 4): continue
data = line.strip().split("|")
if (len(data) < 1): continue
xi, yi = GetIndex(data[0])
cov = float(data[1])
systcov.SetBinContent(xi + 1, yi + 1, cov)
with open("rps_fluxNormCov_analysis2.txt") as f:
count = 0
for line in f:
count += 1
if (count < 4): continue
data = line.strip().split("|")
if (len(data) < 1): continue
xi, yi = GetIndex(data[0])
cov = float(data[1])
normcov.SetBinContent(xi + 1, yi + 1, cov)
totcov.Add(systcov)
totcov.Add(statcov)
totcov.Add(normcov)
outfile.cd()
datahist.Write()
maphist.Write()
counthist.Write()
statcov.Write()
systcov.Write()
totcov.Write()
normcov.Write()
diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.cxx
index 6c70245..641c696 100644
--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.cxx
+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.cxx
@@ -1,255 +1,261 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "T2K_CC0pi_XSec_2DPcos_nu.h"
T2K_CC0pi_XSec_2DPcos_nu::T2K_CC0pi_XSec_2DPcos_nu(std::string name,
std::string inputfile,
FitWeight *rw,
std::string type){
fName = name;
if (fName == "T2K_CC0pi_XSec_2DPcos_nu_I") fAnalysis = 1;
else fAnalysis = 2;
forwardgoing = (type.find("REST") != std::string::npos);
EnuMin = 0;
EnuMax = 10.0;
fBeamDistance = 0.280;
fDefaultTypes = "FIX";
fAllowedTypes = "DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV";
fNDataPointsX = 12;
fNDataPointsY = 10;
Measurement2D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
fIsSystCov = type.find("SYSTCOV") != std::string::npos;
fIsStatCov = type.find("STATCOV") != std::string::npos;
fIsNormCov = type.find("NORMCOV") != std::string::npos;
fPlotTitles = "; P_{#mu} (GeV); cos#theta_{#mu}; d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)";
SetHistograms();
SetupDefaultHist();
// Diagonal covar setup
if (!fIsShape) fAddNormPen = true;
fNormError = 0.089; // Set from covar mat instead...
// Get Scaling
fScaleFactor = ((fEventHist->Integral("width")/(fNEvents+0.)) * 1E-38 /
(TotalIntegratedFlux()));
};
bool T2K_CC0pi_XSec_2DPcos_nu::isSignal(FitEvent *event){
return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, forwardgoing);
};
void T2K_CC0pi_XSec_2DPcos_nu::FillEventVariables(FitEvent* event){
double pmu = -999.9;
double CosThetaMu = -999.9;
// Loop over all particles
for (UInt_t j = 2; j < event->Npart(); ++j){
// Muon section
if ((event->PartInfo(j))->fPID == 13){
// Now find the kinematic values and fill the histogram
pmu = (event->PartInfo(j))->fP.Vect().Mag()/1000;
CosThetaMu = cos(((event->PartInfo(0))->fP.Vect().Angle((event->PartInfo(j))->fP.Vect())));
continue;
}
}
fXVar = pmu;
fYVar = CosThetaMu;
return;
};
// Modification is needed after the full reconfigure to move bins around
// Otherwise this would need to be replaced by a TH2Poly which is too awkward.
void T2K_CC0pi_XSec_2DPcos_nu::ConvertEventRates(){
// Do standard conversion.
Measurement2D::ConvertEventRates();
if (fAnalysis == 1){
// Following code handles weird ND280 Binning
int nbins = this->fMCHist->GetNbinsX() + 1;
double total = 0.0;
// Y = 1
total = 0.0;
for (int i = 3; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(1);
total += this->fMCHist->GetBinContent(i, 1) * width;
this->fMCHist->SetBinContent(i,1,0);
}
this->fMCHist->SetBinContent(3, 1, total / (1.0 * 29.6));
// Y = 2
total = 0.0;
for (int i = 5; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(2);
total += this->fMCHist->GetBinContent(i, 2)* width;
this->fMCHist->SetBinContent(i,2,0);
}
this->fMCHist->SetBinContent(5, 2, total / (0.6 *29.4));
// Y = 3
total = 0.0;
for (int i = 7; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(3);
total += this->fMCHist->GetBinContent(i, 3)* width;
this->fMCHist->SetBinContent(i, 3,0);
}
this->fMCHist->SetBinContent(7, 3, total/ (0.1 * 29.2));
// Y = 4
total = 0.0;
for (int i = 7; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(4);
total += this->fMCHist->GetBinContent(i, 4)* width;
this->fMCHist->SetBinContent(i, 4,0);
}
this->fMCHist->SetBinContent(7, 4, total / (0.1 * 29.2));
// Y = 5
total = 0.0;
for (int i = 8; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(5);
total += this->fMCHist->GetBinContent(i, 5)* width;
this->fMCHist->SetBinContent(i,5,0);
}
this->fMCHist->SetBinContent(8, 5, total / (0.05 * 29.0));
// Y = 6
total = 0.0;
for (int i = 9; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(6);
total += this->fMCHist->GetBinContent(i, 6)* width;
this->fMCHist->SetBinContent(i, 6,0);
}
this->fMCHist->SetBinContent(9, 6, total / (0.05 * 28.5));
// Y = 7
total = 0.0;
for (int i = 8; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(7);
total += this->fMCHist->GetBinContent(i, 7)* width;
this->fMCHist->SetBinContent(i, 7,0);
}
this->fMCHist->SetBinContent(8, 7, total/ (0.04 * 28.0));
// Y = 8
total = 0.0;
for (int i = 11; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(8);
total += this->fMCHist->GetBinContent(i, 8)* width;
this->fMCHist->SetBinContent(i, 8,0);
}
this->fMCHist->SetBinContent(11, 8, total / (0.4 * 27.0));
// Y = 9
total = 0.0;
for (int i = 9; i < nbins; i++){
double width = this->fMCHist->GetXaxis()->GetBinWidth(i) * this->fMCHist->GetYaxis()->GetBinWidth(9);
total += this->fMCHist->GetBinContent(i, 9)* width;
this->fMCHist->SetBinContent(i,9,0);
}
this->fMCHist->SetBinContent(9, 9, total / (0.02 * 25.0));
}
return;
}
void T2K_CC0pi_XSec_2DPcos_nu::SetHistograms(){
// Open file
std::string infile = FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root";
TFile* rootfile = new TFile(infile.c_str(), "READ");
TH2D* tempcov;
// ANALYSIS 2
if (fAnalysis == 2){
// Get Data
fDataHist = (TH2D*) rootfile->Get("analysis2_data");
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle((fName + "_data").c_str(),
(fName + "_data" + fPlotTitles).c_str());
// Get Map
fMapHist = (TH2I*) rootfile->Get("analysis2_map");
fMapHist->SetDirectory(0);
fMapHist->SetNameTitle((fName + "_map").c_str(),
(fName + "_map" + fPlotTitles).c_str());
// Get Syst/Stat Covar
TH2D* tempsyst = (TH2D*) rootfile->Get("analysis2_systcov");
TH2D* tempstat = (TH2D*) rootfile->Get("analysis2_statcov");
TH2D* tempnorm = (TH2D*) rootfile->Get("analysis2_normcov");
// Create covar [Default is both]
tempcov = (TH2D*) tempsyst->Clone();
tempcov->Reset();
if (fIsSystCov) tempcov->Add(tempsyst);
if (fIsStatCov) tempcov->Add(tempstat);
if (fIsNormCov) tempcov->Add(tempnorm);
if (!fIsSystCov && !fIsStatCov && !fIsNormCov){
tempcov->Add(tempsyst);
tempcov->Add(tempstat);
tempcov->Add(tempnorm);
}
// SARAS ANALYSIS
} else if (fAnalysis == 1){
- cout << "HELP AWKWARD BINNING";
+
+ //TODO (P.Stowell) Add a TH2Poly Measurement class
+ ERR(FTL) << " Analysis 1 is not yet available due to its awkward binning!" << endl;
+ ERR(FTL) << "If you want to use it, add a TH2Poly Class!" << endl;
+ throw;
+
}
// Setup Covar
int nbins = tempcov->GetNbinsX();
fFullCovar = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++){
for (int j = 0; j < nbins; j++){
-
+
+ cout << "Filling row " << i << " " << j << " " << nbins <<endl;
(*fFullCovar)(i,j) = tempcov->GetBinContent(i+1,j+1);
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(covar);
// Set Data Errors
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38);
// Remove root file
rootfile->Close();
return;
};
diff --git a/src/Utils/StatUtils.cxx b/src/Utils/StatUtils.cxx
index a2e3b34..4c818ea 100644
--- a/src/Utils/StatUtils.cxx
+++ b/src/Utils/StatUtils.cxx
@@ -1,1047 +1,1049 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "TH1D.h"
#include "StatUtils.h"
//*******************************************************************
Double_t StatUtils::GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask){
//*******************************************************************
Double_t Chi2 = 0.0;
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// Add MC Error to data if required
if (FitPar::Config().GetParB("statutils.addmcerror")){
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_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)*1E38;
double oldval = (*newcov)(i,i);
(*newcov)(i,i) = oldval + mcerr*mcerr;
}
// Reset the calc_cov to new invert
delete calc_cov;
calc_cov = GetInvert(newcov);
// Delete the tempcov
delete newcov;
}
// iterate over bins in X (i,j)
for (int i = 0; i < calc_data->GetNbinsX(); i++){
for (int j = 0; j < calc_data->GetNbinsX(); j++){
if (calc_data->GetBinContent(i+1) != 0 || calc_mc->GetBinContent(i+1) != 0) {
LOG(DEB)<< "i j = "<<i<<" "<<j<<std::endl;
LOG(DEB)<<"Calc_data mc i = "<<calc_data->GetBinContent(i+1)
<<" "<<calc_mc->GetBinContent(i+1)
<<" Dif = "
<<( calc_data->GetBinContent(i+1) - calc_mc->GetBinContent(i+1) )
<<std::endl;
LOG(DEB)<<"Calc_data mc i = "<<calc_data->GetBinContent(j+1)
<<" "<<calc_mc->GetBinContent(j+1)
<<" Dif = "
<<( calc_data->GetBinContent(j+1) - calc_mc->GetBinContent(j+1) )
<<std::endl;
LOG(DEB)<<"Covar = "<< (*calc_cov)(i, j) <<std::endl;
LOG(DEB)<<"Cont chi2 = " \
<<( ( calc_data->GetBinContent(i+1) - calc_mc->GetBinContent(i+1) ) \
* (*calc_cov)(i, j) \
* 1E76 \
* ( calc_data->GetBinContent(j+1) - calc_mc->GetBinContent(j+1) ) )
<<" "<<Chi2<<std::endl;
Chi2 += ( ( calc_data->GetBinContent(i+1) - calc_mc->GetBinContent(i+1) ) * \
(*calc_cov)(i, j) * 1E76 * \
( calc_data->GetBinContent(j+1) - calc_mc->GetBinContent(j+1) ) );
} else {
/*
std::cerr << "Error on bin (i,j) = (" << i << "," << j << ")" << std::endl;
std::cerr << "data->GetBinContent(i+1) = " << calc_data->GetBinContent(i+1) << std::endl;
std::cerr << "mc->GetBinContent(i+1) = " << calc_mc->GetBinContent(i+1) << std::endl;
std::cerr << "Adding zero to chi2 instead of dying horrifically " << std::endl;
*/
Chi2 += 0.;
}
}
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromCov( TH2D* data, TH2D* mc,
TMatrixDSym* invcov, TH2I* map, TH2I* mask){
//*******************************************************************
// Generate a simple map
if (!map) {
map = StatUtils::GenerateMap(data);
}
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate 1D chi2 from 1D Plots
Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD( TH1D* data, TH1D* mc,
TMatrixDSym* cov, TH1I* mask){
//*******************************************************************
Double_t Chi2 = 0.0;
TMatrixDSym* calc_cov = (TMatrixDSym*) cov->Clone();
TH1D* calc_data = (TH1D*) data->Clone();
TH1D* calc_mc = (TH1D*) mc->Clone();
// If a mask if applied we need to apply it before the matrix is inverted
if (mask) {
calc_cov = StatUtils::ApplyMatrixMasking(cov, mask);
calc_data = StatUtils::ApplyHistogramMasking(data, mask);
calc_mc = StatUtils::ApplyHistogramMasking(mc, mask);
}
// Decompose matrix
TDecompSVD LU = TDecompSVD((*calc_cov));
LU.Decompose();
TMatrixDSym* cov_U = new TMatrixDSym(calc_data->GetNbinsX(), LU .GetU().GetMatrixArray(), "");
TVectorD* cov_S = new TVectorD( LU.GetSig() );
// Apply basis rotation before adding up chi2
Double_t rotated_difference = 0.0;
for (int i = 0; i < calc_data->GetNbinsX(); i++){
rotated_difference = 0.0;
// Rotate basis of Data - MC
for (int j = 0; j < calc_data->GetNbinsY(); j++)
rotated_difference += ( calc_data->GetBinContent(j+1) - calc_mc ->GetBinContent(j+1) ) * (*cov_U)(j,i) ;
// Divide by rotated error cov_S
Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i);
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
delete cov_U;
delete cov_S;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD( TH2D* data, TH2D* mc,
TMatrixDSym* cov, TH2I* map, TH2I* mask){
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D* data_1D = MapToTH1D(data, map);
TH1D* mc_1D = MapToTH1D(mc, map);
TH1I* mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
double StatUtils::GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask) {
//*******************************************************************
// If just an event rate, for chi2 just use Poission Likelihood to calculate the chi2 component
double chi2 = 0.0;
TH1D* calc_data = (TH1D*)data->Clone();
TH1D* calc_mc = (TH1D*)mc->Clone();
// Apply masking if required
if (mask) {
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Iterate over bins in X
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double dt = calc_data->GetBinContent(i+1);
double mc = calc_mc->GetBinContent(i+1);
if (dt == 0 or mc == 0) {
LOG(REC) << "Found no MC in bin " << i << " for " << calc_data->GetName() << " (" << calc_data->GetBinLowEdge(i) << " - " << calc_data->GetBinLowEdge(i+1) << ")" << std::endl;
continue;
}
// Do the chi2 for Poisson distributions
chi2 += 2 * (mc - dt + (dt*log(dt/mc)));
/*
LOG(REC)<<"Evt Chi2 cont = "<<i<<" "
<<mc<<" "<<dt<<" "
<<2 * (mc - dt + (dt+0.) * log((dt+0.) / (mc+0.)))
<<" "<<Chi2<<std::endl;
*/
}
// cleanup
delete calc_data;
delete calc_mc;
return chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromEventRate(TH2D* data, TH2D* mc, 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::GetChi2FromEventRate(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromDiag(TH1D* data, TH1D* mc, TH1I* mask){
//*******************************************************************
// Currently just a placeholder!
(void) data;
(void) mc;
(void) mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromDiag(TH2D* data, TH2D* mc, 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 MLE = StatUtils::GetLikelihoodFromDiag(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask){
//*******************************************************************
// Currently just a placeholder !
(void) data;
(void) mc;
(void) invcov;
(void) mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromCov( 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 from 1D
Double_t MLE = StatUtils::GetLikelihoodFromCov(data_1D, mc_1D, invcov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask){
//*******************************************************************
// Currently just a placeholder!
(void) data;
(void) mc;
(void) cov;
(void) mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromSVD( 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 MLE = StatUtils::GetLikelihoodFromSVD(data_1D, mc_1D, cov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromEventRate(TH1D* data, TH1D* mc, TH1I* mask){
//*******************************************************************
// Currently just a placeholder!
(void) data;
(void) mc;
(void) mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromEventRate(TH2D* data, TH2D* mc, 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 MLE = StatUtils::GetChi2FromEventRate(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Int_t StatUtils::GetNDOF(TH1D* hist, TH1I* mask){
//*******************************************************************
TH1D* calc_hist = (TH1D*)hist->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<Double_t> 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 non_diagonal = false;
- for (int i = 0; i < new_mat->GetNrows(); i++){
- for (int j = 0; j < new_mat->GetNrows(); j++){
+ 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) {
- non_diagonal=true;
+ diagonal=false;
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 (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();
- new_mat = new TMatrixDSym(new_mat->GetNrows(), LU.GetU().GetMatrixArray(), "");
+ new_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), "");
sleep(2);
return new_mat;
}
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH1D* hist, double norm){
//*******************************************************************
if (!mat) mat = MakeDiagonalCovarMatrix(hist);
int nbins = mat->GetNrows();
TMatrixDSym* new_mat = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++){
for (int j = 0; j < nbins; j++){
double valx = hist->GetBinContent(i+1)*1E38;
double valy = hist->GetBinContent(j+1)*1E38;
(*new_mat)(i,j) = (*mat)(i,j) + norm*norm*valx*valy;
}
}
// Swap the two
delete mat;
mat = new_mat;
return;
};
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, TH2I* map ){
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1D* data_1D = MapToTH1D(data, map);
StatUtils::ForceNormIntoCovar(mat, data_1D, norm);
delete data_1D;
return;
}
//*******************************************************************
TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH1D* data, double scaleF){
//*******************************************************************
TMatrixDSym* newmat = new TMatrixDSym(data->GetNbinsX());
for (int i = 0; i < data->GetNbinsX(); i++){
(*newmat)(i,i) = data->GetBinError(i+1) * data->GetBinError(i+1) * scaleF * scaleF;
}
return newmat;
}
//*******************************************************************
TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH2D* data, TH2I* map, double scaleF){
//*******************************************************************
if (!map) map = StatUtils::GenerateMap(data);
TH1D* data_1D = MapToTH1D(data, map);
return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF);
};
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, double scale){
//*******************************************************************
+ // Check
+ if (cov->GetNrows() != data->GetNbinsX()){
+ ERR(WRN) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << endl;
+ }
+
// Set bin errors form cov diag
for (int i = 0; i < data->GetNbinsX(); i++){
-
- if (data->GetBinContent(i+1) == 0.0) continue;
-
data->SetBinError(i+1, sqrt((*cov)(i,i)) * scale );
-
}
return;
}
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map, double scale){
//*******************************************************************
// Create map if required
if (!map) map = StatUtils::GenerateMap(data);
// Set Bin Errors from cov diag
int count = 0;
for (int i = 0; i < data->GetNbinsX(); i++){
for (int j = 0; j < data->GetNbinsY(); j++){
if (data->GetBinContent(i+1,j+1) == 0.0) continue;
count = map->GetBinContent(i+1,j+1)-1;
data->SetBinError(i+1,j+1, sqrt((*cov)(count,count)) * scale );
}
}
return;
}
//*******************************************************************
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;
}

File Metadata

Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:36 AM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4874165
Default Alt Text
(45 KB)

Event Timeline